Category Archives: ONA

Obliczenia naukowe

ONA 7 – Układy równań liniowych

Dziś zajmujemy się układami równań liniowych w reprezentacji macierzowej.

Wykład będzie prowadzony prz tablicy, więc slajdów nie ma, ale potrzebne materiały ( i dużo więcej niż nam potrzeba) są dostępne w materiałach do wykładu z Metod numerycznych . Nas w szczególności interesują wykłądy 5 (eliminacja Gaussa) i 7 (uwarunkowanie problemu). Można też zajrzeć do wykładu 8 (macierze rzadkie), ale jest on o dla nas zdecydowanie zbyt obszerny, a jest to dla nas tylko temat niejako poboczny.

Jeśli chodzi o operacje algebry liniowej w pythonie, przydadzą nam się następujące biblioteki:

numpy.linalg - podstawowe operacje na macierzach, m.in. cond (A), det(A)

scipy.linalg - nieco więcej operacji na macierzach, m.in. hilbert(N), solve(A,B)

scipy.sparse.linalg - operacje na macierach rzadkich m.in. spsolve(A,B)

scipy.io - wczytywanie różnych ciekawych formatów – np. macierzy rzadkich w formacie matlaba loadmat(f)

Zadania na laboratorium (ew. nierozwiązane jako praca domowa):

  1. Stwórz macierz Hilberta H[i,j]=1/i+j-1 dla sensownie dużego rozmiaru macierzy (100×100, 200×200) i wyświetl ją przy pomocy imshow()/
  2. Przypomnij sobie rozwiązywanie ukł. równań poprzez rozkład LU. Zaimplementuj ręcznię tę metodę (w wersji bez wyboru el. głównego) i spróbój rozwiązać tym sposobem układ postaci A= np.array([[1e-18,1.0],[1.0,2.0]]) B=np.array([1.0,4.0]). Porównaj wynik z wynikiem scipy.linalg.solve(A,B). Dlaczego w przypadku tej macierzy wybór elementu głównego ma takie znaczenie?
  3. Korzystając z metody scipy.rand(1000,1000) stwórz macierz losową. Użyj jej do przetestowania rozkładu macierzy  LU (scipy.linalg.lu) w wersji z permutacją i bez. Rozwiąż taki układ (dla b=scipy.rand(1000)) przy pomocy scipy.linalg.solve()
  4. Załóżmy, że mamy teraz jedną macierz A (my użyjemy losowej) i bardzo wiele różnych warunków brzegowych B (my wylosujemy 1000). Spróbuj rozwiązać wszystkie układy równań powstałe z przyrównania tej samej macierzy A do wielu różnych wektorów B. Czy można zamiast korzystać wielokrotnie z funkcji solve() coś przyspieszyć? Np. korzystając z funkcji lu()?
  5. Duże macierze równań liniowych często powstają w problemach inżynierskich. Pobierz jedną z macierzy z kolekcji układów powstałych przy projektowaniu elementów samolotów Boeing (np. nr 38 ) i spróbuj wczytać ją jako macierz rzadką do ipython’a. Następnie rozwiąż ją dla losowych warunków brzegowych b.

ONA – Zadanie zaliczeniowe 2

Naszym zadaniem będzienapisanie programu, który będzie przetwarzał strumienie danych genomicznych zapisane w plikach sgr . Są to bardzo proste pliki, które pozwalają opisywać zmienność parametru wzdłuż chromosomów. Każda linijka ma tylko 3 wartości:

chromosom pozycja wartość

Wartości są oddzielane znakiem tabluacji (“\t”), pozycja jest liczbą całkowitą dodatnią, a wartość jest liczbą zmiennoprzecinkową.

Zakładamy, że wartości pomiędzy punktami w pliku są odcinkami prostoliniowymi, co pozwala nam na podstawie pliku sgr, wyliczyć wartości dla wszystkich pozycji chromosomu.

Przykładowe pliki sgr dla dwóch sygnałów i dwóch filtrów są w pliku Zad2

Pliki te mogą być duże – np. większe niż dostępna pamięć, co oznacza, że nie możemy ich wczytać w całości do wektora.

Będziemy chcieli, żeby nasz program oferował następujące funkcje:

- Suma, różnica, iloczyn, iloraz – zwróć wynik operacji arytmetycznej na  dwóch sygnałach dla tych samych chromosomów. Zakładamy,  że na wejściu mamy dwa pliki, które opisują sygnał wzdłuż tych samych chromosomów, choć niekoniecznie w tych samych pozycjach. OPeracje są wykonywane “po pozycjach”. Zwracamy plik sgr, który opisuje funkcję wynikową (4 pkt)

- wygładzanie – przy pomocy średniej kroczącej o zadanej długości (w parach zasad), zwróć sygnał wygładzony (2 pkt)

- splot – Tym razem na wejściu bierzemy jeden plik sgr z sygnałem (wzdłuż chromosomów) i drugi plik sgr z wartościami funkcji filtra (np. filtr prostokątny, lub Gaussowski), wzdłuż sztucznego chromosomu o nazwie “filtr”. Wynikiem powinien być splot funkcji, czyli przefiltrowany sygnał 6  pkt)

- upraszczanie – dla zadanego pliku sgr wykryj, które punkty są niepotrzebne (tzn po ich usunięciu sygnał nie zmieni się) i zwróć sygnał pozbawiony tych linijek (3 pkt).

Rozwiązanie polega na napisani 7 funkcji, każda z nich powinna działać na 2 plikach (+ parametr szerokość dla wygładzania), wszystko w jednym pliku .py

Rozwiązania podpisane, z dopiskiem [ONA-2018-2] w temacie proszę wysyáć mailem do mnie, do 22. kwietnia 2018.

ONA 6 – Kompresja

Slajdy są tutaj ONA6-kompresja

Zadania na dziś:

1. Zapoznaj się z modułem gzip  zapisz plik tekstowy w pliku skompresowanym ze skryptu pythona, zdekompresuj go przy pomocy programu gunzip i odwrotnie

2. Weźmy takie pomocnicze funkcje w pythonie opisujące dyskretną transformatę kosinusową:

import urllib2,io
import Image
from scipy import fftpack

image_url='http://i.imgur.com/8vuLtqi.png'

def get_image_from_url(image_url='http://i.imgur.com/8vuLtqi.png', size=(128, 128)):
    file_descriptor = urllib2.urlopen(image_url)
    image_file = io.BytesIO(file_descriptor.read())
    image = Image.open(image_file)
    img_color = image.resize(size, 1)
    img_grey = img_color.convert('L')
    img = np.array(img_grey, dtype=np.float)
    return img

def get_2D_dct(img):
    """ Get 2D Cosine Transform of Image
    """
    return fftpack.dct(fftpack.dct(img.T, norm='ortho').T, norm='ortho')

def get_2d_idct(coefficients):
    """ Get 2D Inverse Cosine Transform of Image
    """
    return fftpack.idct(fftpack.idct(coefficients.T, norm='ortho').T, norm='ortho')

def get_reconstructed_image(raw):
    img = raw.clip(0, 255)
    img = img.astype('uint8')
    img = Image.fromarray(img)
    return img

wykonaj program, który dokonuje Ntego przybliżenia obrazu, przy pomocy zerowania wartośći powyżej Ntego wiersza i Ntej kolumny macierzy get_2D_dct(img) i wyświetl kilkanaście pierwszych przybliżeń.

3. Wyrysuj wykres Entropii dla kanału binarnego w zależności od P(1)

4. Napisz program, który dla zadanego pliku tekstowego tworzy tablicę częstotliwości znaków, kody Huffmana dla wszystkich znaków oraz wylicza entropię tego pliku i długość kodu Huffmana, który opisywałby cały plik

 

ONA – 5 – Analiza obrazu

Dzisiejsze zajęcia poświęcimy na analizę obrazów 2d. Slajdy są tu ONA5-Obrazy

Przydać mogą się pakiety: scipy.ndimage i matplotlib

Zadania na dziś:

0. Pozyskaj przy pomocy kinect swoje zdjęcie w pracowni komputerowej zarówno w postaci macierzy głebokości, jak i obrazów rgb (przyklady tu i tu)

1. Przekształć obraz RGB do skali szarości poprzez uśrednienie składowych.

2. Wyświetl te obrazy na swoim komputerze przy pomocy biblioteki matplotlib

3. Wyświetl histogram obrazu i dokonaj jego “wyrównania” pisząc program w  języku python

f(x) = (x-min)/(max-min)*255 dla parametrów min i max będących maksymalną i minimalną jasnością punktu.

4. Zastosuj filtr Gaussowski rozmiaru k, aby “wygładzić” obraz i obejrzyj wyniki dla różnych k (3,5,7,9,11,…)

5. Napisz program wykrywający krawędzie przy pomocy filtru Sobel’a

6(*). Napisz program wykrywający cienie na obrazie z miernika odległości.

7(*). Użyj biblioteki pydicom do wczytania przekrojów przez głowę z projektu visible human. Napisz program, który zamieni te przekroje poziome na prekroje pionowe. Użyj wygładzania krawędzi.

 

ONA 4 – przetwarzanie sygnałów

Dziś zajmujemy się przetwarzaniem sygnałów, slajdy są dostępne tu: ONA4-DSP

Interesują nas przede wszystkim biblioteki scipy.signal i scipy.fftpack

1. Wygeneruj sygnał funkcji wielomianowej o zaburzeniu Gaussowskim, dla 10000 punktów. Wykonaj uśrednienia tego sygnału przy pomocy średniej kroczącej, albo splotu przy pomocy filtru kwadratowego, trójkątnego lub Gaussowskiego. WYniki zaprezentuj na wykresie

2. Pobierz sygnał o zajętości nukleosomowej z plików pochodzących z eksperymentu MNase-Seq ( dane tu w formacie bedgraph ) . Wczytaj go do numpy jako wektor. Przedstaw go na wykresie. Używając szybkiej trasformaty fouriera oblicz widmo Fourierowskie tego sygnału. Przedstaw je na innym wykresie. Wyzeruj część widma wysokich częstotliwości i dokonaj odwrotnej transformaty Fouriera, aby uzyskać wygładzony sygnał zajętości nukleosomami. Spróbuj tak dobrać parametry filtra, aby wykres po odwrotnej transformacie Fouriera miał okres zbliżony do 160-200 par zasad.

ONA – zadanie 1 – fraktale

Naszym zadaniem będzie napisanie programu, który oblicza kolejne przybliżenia zbiorów fraktalnych. Inspiracją będą dla nas zbiory Julii  i Mandelbrota, ale nasz program będzie bardziej ogólny.

Podstawową funkcją naszego programu, będzie funkcja:

M=approx(xrange,yrange,max_iter,max_value,f),

która zwraca macierz wyników iteracji funkcji f dla wszystkich par wartości z wektorów xrange  i yrange, interpretowanych jako liczby zespolone.

M[a,b] = liczba iteracji funkcji f dla wartości xrange[a]+j*yrange[b], dla których wartość funkcji nie wykracza poza max_value, przy czym przybliżamy wynik poprzez wykonanie maksymalnie max_iter iteracji. Tzn jeśli po wykonaniu max_iter iteracji nigdy nie wykroczyliśmy poza zakres max_value, to zwracamy max_iter, w przeciwnym wypadku zwracamy ten numer iteracji, dla którego nasza wartość f(z,n) przekroczyla max_value. argument f powinien być funkcją przyjmującą 2 argumenty: z i n, gdzie z to liczba zespolona a n to stopień iteracji.

Oprócz funkcji approx(..), powinnismy zaimplementować funkcje Mandelbrot(z,n) i Julia(z,n), implementujące iterację dla zbioru Julii i Mandelbrota.

Nasz program powinien umożliwiać:

  • Wywołanie z linii komend i podawanie parametrów: zakresy x,y, n, typ funkcji, typ wyjścia, dodatkowe parametry (np. c dla zbioru Julii)
  • Wyświetlanie wyniku przy pomocy  funkcji imshow() wg wybranej skali kolorów
  • zapisywanie obrazka w pliku png, pdf lub macierzy w pliku  (przy pomocy operacji np.save())

Punktacja:

  • funkcja approx – 3 pkt
  • funkcje mandelbrot, julia – 2 pkt
  • obsługa z linii komend – 2 pkt
  • wyswietlanie – 2 pkt
  • zapis do pliku – 1 pkt

Rozwiązania wysyłamy do 25 III na adres e-mail wykładowcy z dopiskiem [ONA-1-2018]- w tytule 

ONA 3 – Wykresy

Dzisiaj zajmujemy się wykresami.

Slajdy są tu ONA3-wykresy

Zadania na dziś:

1. Narysuj wykres podobny do tego ze slajdu nr. 14  z wykładu. Uwzględnij opisy osi, zakresy wartości, kolory, znaczniki wartości, legendę i tytuły osi oraz wykresów.

2. Narysuj wykres 10*sin(x)+normal(0,1) dla 1000 punktów w zakresie 0,10. Przyda się funkcja normal().

3. Napisz program, który wykona 1000 prób rzutu 10 kostkami. Narysuj na jednym wykresie 3 przenikające się histogramy (alpha=0.25): suma oczek na kościach parzystych (0,2,4,6,8), suma oczek na kościach nieparzystych (1,3,5,7,9), suma oczek na wszystkich kościach. Może przydać się funkcja randint().

4. Przedstaw dane z wykresu 3. w postaci 3 wykresów pudełkowych (boxplot)

5. Narysuj powierzchnię 3d (x,y,z) in linspace(0,1,100) gdzie z=y*cos(x)+(1-y)*sin(x), Może być przydatny ten przykład

5. Narysuj mapę ciepła dla tej samej powierzchni. Użyj różnych map kolorów.

 

ONA 2 – Macierze i NumPy

Slajdy do wykładu – ONA2-Macierze

Jako interpretera pythona mogą Państwo użyć tego (Trochę porad jak używać tu )

Ćwiczenia:

  1. Stwórz macierze A i B o wymiarach  1000×1000 zawierające wartości A[i,j]==i*3-j*5 i B[i,j]==np.sqrt(A[i,j]) (Zauważ, że to wymaga wartości urojonych, a więc macierzy typu complex64 i wartość sqrt(-1) ->.j)
  2. Napisz funkcję poz(indeksy, shape), która dla zadanej krotki indeksów zwraca liniową pozycję elementu o zadanych indeksach w wielowymiarowaj macierzy o kształcie shape.
  3. Korzystając z broadcastingu i mnożenia macierzy zmień znak elementów w macierzy B, których suma indeksów jest nieparzysta
  4. korzystając z metody sort() , posortuj macierz 3-wymiarową wg. drugiej współrzędnej i obejrzyj wynik.
  5. Korzystając z funkcji frompyfunc() napisz funkcje wektorowe, które liczą a) sumę kwadratów dwóch macierzy  (dwie macierze do jednej)
    b) iloraz i resztę z dzielenia całkowitoliczbowego przez 17 (jedna macierz wejściowa i dwie macierze wyjściowe)
  6. Korzystając z modułu time, porównaj prędkość funkcji liczących N potęg dwójki – na liście i w wektorze. Dla jakich wartości N warto używać macierzy?