Month: April 2017

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