Category: ONA

Nasze drugie zadanie będzie polegało na tym, aby napisać program pozwalający  skonstruować interpolację splajnami dla danych o ekspresji genów.

Jako przykładu użyjemy danych z pracy: S. Hooper et al. 

Mamy tutaj geny pogrupowane według profilu ekspresji w pliku 1 oraz ich profile ekspresji w pliku Timecourse_normalized_mean.xls

W pracy tej wyróżniono 3 klasy genów, każda podzielona na podgrupy według ich profilu ekspresji podczas rozwoju zarodkowego:

Zadaniem Państwa jest napisanie programu, który na podstawie tych dwóch plików dokonuje następujących analiz:

1. Wczytuje nazwy genów w poszczególnych klasach i grupach (z pliku csv: 2 pkt, z pliku xls +1)
2. Wczytuje profile ekspresji wszystkich genów z pliku (csv: 2pkt, xls +1)
3. Znajduje średnie profile każdej z grup genów oraz dla każdej z klas (3 pkt)
4. Znajduje interpolacje splajnami (liniowymi lub kubicznymi) dla profilu każdego genu, grupy i klasy (4 pkt)
5. Tworzy po 3 wykresy zawierające gładkie (n=200) wykresy splajnów 1 i 2 rzędu (3 wykresy dla liniowych i 3 wykresy dla kubicznych) dla genów, grup i klas (4 pkt, +2 dla najładnieszych wykresów).

Jako rozwiązanie (do 8 V 2017) wysyłamy skrypt i wynikowe wykresy. Uruchomienie skryptu powinno spowodować powstanie wykresów w pdf. Jeśli korzystali Państwo z plikóœ pośrednich (csv) to proszę je załączyć. W temacie maila jak zwykle [ONA-2017-2].

Dzisiaj zajmiemy się Interpolacją funkcji przy pomocy wielomianów i funkcji sklejanych.

Teoria z wykładu jeśli chodzi o wielomiany znajduje się tu a jeśli chodzi o funkcje sklejane tu.

Większość interesujących nas dzisiaj funkcji znajdziemy w module scipy.interpolate, ale najprostsze funkcje polyfit  i poly1d znajdują się w module numpy.  W module interpolate interesują nas funkcje:

• InterpolatedUnivariateSpline
• UnivariateSpline
• lagrange
• barycentric_interpolate

Na laboratorium będziemy rozwiązywać następujące problemy:

1. Spróbuj zinterpolować funkcję sinus(x) na przedziale [0,math.pi] korzystając z równoodległych N węzłów (np. wygenerowanych używając np.linspace(0,math.pi,N)). Jak zachowuje się błąd średniokwadratowy tej interpolacji dla punktów np.linspace(0.math.pi,1000), gdy N rośnie?

2. korzystając z przykładu pokazanego na wykładzie, gdzie wartości y_i=[1,1,1,2], czy potrafisz dobrać pozycje x_i=[1,x2,x3,4] tak aby uzyskać dowolnie dużą amplitudę interpolowanego (polyfit(x_i,y_i,s=0) wielomianu na przedziale [1,4]? Czy błąd aproksymacji (polyfit, k=2) wielomianem stopnia 2 jest tak samo duży?

3. Spróbuj interpolować te same dane przy pomocy krzywych sklejanych. Czym różnią się splajny dla różnych k (k=1,2,3)?

4. Wygeneruj zaburzone obserwacje wg funkcji y=log_2(x) + scipy.stats.normal() w wielu (~1000) punktach na przedziale 1..100. Czy lepiej będzie interpolować, czy aproksymować aby znaleźć kształt funkcji?

5. Spróbuj użyć splajnów do “wygładzenia swojego rozwiązania zadania zaliczeniowego nr 1. Czy może przyda się funkcja RectBivariateSpline

Dzisiaj na wykładzie omówiliśmy zasadniczo tematy zawarte w wykładzie 12. z metod numerycznych. Jeśli ktoś chciałby doczytać to znajdzie materiały tutaj

Zadania na lab:

0. Rozważmy trzy punkty na płaszczyźnie: (0,6), (1,0) i (2,0). Jaka prosta przechodzi najbliżej nich? ułóż układ równań linioowych, który można rozwiązać metodą najmniejszych kwadratów. Wykorzystaj funkcje scipy.linalg.qr aby zobaczyć rozkład QR macierzy A. Użyj funkcji scipy. linalg.lstsq aby znaleźć rozwiązanie. jakie jest znaczenie wartośći zwróconych przez tę funkcje?

1. Rozważmy dane:

x	f(x)
0.00	4.00000000000000e+00
1.25	3.28650479686019e+00
2.50	3.08208499862390e+00
3.75	3.02351774585601e+00
5.00	3.00673794699909e+00
6.25	3.00193045413623e+00
7.50	0.00055308437015e+00
8.75	3.00015846132512e+00
10.00	3.00004539992976e+00

Jak będzie wyglądało dopasowanie met. najlepszych kwadratów funkcji f(x)=a+b*exp(-x) do tych danych?

2. Rozważmy dane o obwodzie pnia, trees-stripped ( do wczytywania przyda się funkcja scipy.loadtext). Kolejne kolumny oznaczają tu:

• obwód pnia
• wysokość drzewa
• objętość pozyskanego drewna.

Spróbuj dopasować (metodą najmniejszych kwadratów objętość drzewa jako funkcję:

• kombinację liniową obwodu pnia i wysokości drzewa
• iloczynu wysokości przez obwód
• kombinację liniową powyższych

Gdzie uzyskujemy najmniejszy błąd przybliżenia?

3. Przedstaw na wykresie wynik zadania 1 (wykres punktowy obserwacji i dużo gęściej próbkowany wykres liniowy znalezionej funkcji). Jaki jest problem naszego rozwiązania? czy można jakoś pomóc sobie używając ważonego problemu średnich kwadratów? jak to zrobić dla wagi jednej z obserwacji=0? a jak dla wagi 0.1?

 

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)

Dodatkowo:

Warto zajrzeć pod adres https://notch.mimuw.edu.pl:8888/ gdzie zainstalowałem dla Państwa serwer jupyter. Można tam korzystać z wyświetlania seaborn, matplotlib i bokeh. Każdy z Państwa ma swój “notebook” do wykonywania obliczeń.

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

1. Stwórz macierz Hilberta H[i,j]=1/i+j-1 i wyświetl ją przy pomocy serwera jupyter.
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.

 

http://www.cise.ufl.edu/research/sparse/matrices/Boeing/bcsstk38.html

Zgodnie z dzisiejszymi ustaleniami:

I projekt jest do oddania dziś do końca dnia.

II projekt (interpolacja splajnami) jest do oddania do 8 V

III projekt (iteracyjne dopasowywanie współczynników) będzie do oddania  do 9 VI

Kolokwium będzie 15 maja  (na wykładzie)

 

Zadanie polega na użyciu biblioteki freenect do napisania prostego programu, który w kreatywny sposób wykorzystuje informacje z czujnika głębokości kinect.

Przykładowo, w zeszłym roku powstały gry pozwalające użytkownikom odgadywać liczbę wylosowaną przez komputer poprzez ustawianie się przed komputerem w odległości i sprawdzanie czy to za blisko, czy za daleko.

W przyszłym tygodniu (3.IV) chciałbym, aby na zajęciach zgłosiły się do mnie 2 osobowe zespoły, które chcą ten projekt realizować.

Za 2 tygodnie (10 IV) chciałbym porozmawiać z Państwem o projektach, które chcą Państwo realizować.

Na zajęciach 5 VI będziemy prezentować projekty. Będzie można zdobyć maksymalnie 15 punktów.

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

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.

 

Slajdy  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.

Naszym zadaniem jest napisanie modułu (pliku .py) w języku python (wersja 3), który pozwoli na wykonanie normalizacji kwantylowej zadanych obserwacji.

Normalizacja kwantylowa, to operacja na macierzy, która pozwala zapewnić, że w każdej kolumnie macierzy znajdują się wartości o tym samym rozkładzie kwantyli, przy czym wartości są możliwie bliskie wejściowym.

Opis samej metody (dot. kolumn) można zobaczyć na przykładzie w wikipedii .

Nasz program powinien działać w trybie skryptu i pobierać opcje z linii komend (np. przy użyciu Argparse). Funkcjonalności wymagane to :

– wczytywanie wejścia i zapis wyniku w postaci macierzy (.npy) i tekstu (.csv) 3 pkt
– wczytywanie z excela (.xls) przy pomocy biblioteki pandas (bonus 1pkt)
– wykonanie procedury normalizacji kwantylowej wg kolumn (5 pkt)
– obsługa normalizacji wg wierszy (1 pkt)
– wyswietlenie macierzy wynikowej i pośrrednich, lub zapis do pliku obrazka w png lub pdf (1 pkt)
– dodatkowo można uwzględnić normalizację wg mniejszej liczby kwantyli (np. 100 kwantyli dla macierzy 1000 elementowej), przy założeniu liniowej zmiany wartości pomiędzy kwantylami (bonus 2 pkt)

Zadanie w postaci emaila z załączonym plikiem .py (nieskompresowanym) przesyłamy (z konta studenckiego) do 27 III  (godz. 00:00 w dowolnej strefie czasowej) na adres e-mail bartek@mimuw.edu.pl z dopiskiem [ONA-2017-1] w tytule (subject).

Po terminie można dostać maksymalnie 5 punktów.

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.