Category Archives: ONA

Obliczenia naukowe

ONA 12 – Całkowanie Monte Carlo i miejsca zerowe funkcji nieliniowych

Dziś przede wszystkim zajmowaliśmy się zerami funkcji nieliniowych, ale zrobiliśmy też dygresję (w zasadzie uzupełnienie do wykładu o całkowaniu numerycznym) n.t. metod Monte Carlo. Notatki do równań nieliniowych znajdą Państwo tu. Opis całkowania metodą MC znajdą Państwo np. tu

Zadania na dziś:

1. Napisz program, który znajduje  przybliżenie liczby Pi przy pomocy równania x**2+y**2=1

2. Napisz program, który korzysta z metody bisekcji oblicza wartość pierwiastka kwadratowego z 2. Jakiego równania zera szukamy? ile iteracji potrzeba?

3. Napisz program, który korzysta z metody stycznych, aby obliczyć tę samą wartość. O ile szybciej możemy znaleźć dokładną wartość?

4 (*). Zaimplementuj metodę iteracji prostej Banacha dla rozwiązania równania Keplera.  Możesz posiłkować się tym artykułem

ONA 11 – Obliczenia symboliczne

Dzisiejszy wykład porównujący cechy różnych pakietów do obliczeń jest tu: ONA11-Pakiety-oprogramowania

Zadania na dziś są dość proste i dotyczą przede wszystkim obliczeń symbolicznych i pakietu sympy (można to robić w platfoermie jupyter gdzie jest też zakładka zawierająca przykłady z sympy)

  1. Przy użyciu pakietu sympy, zbadaj podstawowe możliwości obliczeń symbolicznych: stwórz zmienną sumboliczną x (Symbol(‘x’)), stwórz wyrażenie wielomianowe (np. x**2-2) oraz wyrażenie trygonometryczne (np. sin(x)) i policz dla nich pochodne, całki (oznaczone i nieoznaczone) oraz granice (w 0, z prawej i lewej strony oraz , granice 1/x w nieskonczonosci: sympy.oo,).
  2. Zobacz jak wyniki tego typu zadań wyglądają w Wolfram Alpha
  3. Zarejestruj się w systemie SageMath i tam spróbuj rozwiązać te same problemy
  4. (*) Spróbuj rozwiązać symbolicznie układ równań różniczkowych zwyczajnych z poprzednich zajęć (model Lotki-Volterry albo Zombie Invasion) w jednym z wymienionych wcześniej pakietów

ONA – Zadanie zaliczeniowe 3

Naszym zadaniem jest napisanie prostej metody iteracyjnej, która znajdzie lokalnie optymalne przypisanie poziomów ekspresji dla różnych transkryptów genów.

Przypomnijmy – geny u eukariontów podlegają splicingowi – procesowi wycinania intronów zanim dojdzie do translacji. Jeśli transkrypt posiada wiele intronów, możemy otrzymać wiele różnych transkryptów – niekiedy, jak w przypadku genu dSCAM, możeich być naprawdę bardzo wiele.

Mamy więc 3 poziomy organizacji transkrypcji: Geny, transkrypty i eksony. Każdy ekson może należeć do wielu transkryptów, ale każdy transkrypt należy do jednego genu.

Kiedy sekwerncjonujemy mRNA, możemy policzyć ile odczytów przypadło na każdy z eksonów. Przykładowy wynik takiej operacji znajdą Państwo w pliku Reads_on_exons. Opisuje on dla każdego z kilkuset tysięcy exonów liczbę odczytów przypadających na niego w jednym z eksperymentów badania ekspresji genów w ludzkim mózgu prowadzonym w mojej grupie.

Opis przynależnośći eksonów do genów w ludzkim genomie znajdą Państwo w pliku hg38_full.gtf Dla każdego eksonu będzie tam linia przypisująca go do określonego genu i transkryptu.

Naszym zadaniem jest obliczenie najbardziej prawdopodobnych poziomów ekspresji transkryptów zakładając, że ekspresja eksonów jest w przybliżeniu równa sumie ekspresji transkryptów zawierających dany ekson. Oczywiście musimy wziąć pod uwagę pewien błąd, który jest nieuchronnie związany z losowością próbki odczytów oraz z niedoskonałością opisu transkryptów.

Punktacja wygląda następująco:

- zaprojektowanie sensownych struktur danych do przechowywania informacji o genach, transkryptach i exonach – 7 punktów
- wczytanie danych z plików do tychże struktur – 7 punktów
- iterując po genach rozwiązujemy zadanie minimalizacji błędu najmniejszych kwadratów, gdzie dana jest ekspresja na egzonach i przypisanie egzonów do transkryptów i genów, a szukanymi jest ekspresja transkryptów. – 7 punktów
- wypisujemy wyniki do plików tekstowych (4 pkt)
a) eskpresję transkryptów w formacie: nazwa genu, nazwa transkryptu, wartość ekspresji, lista egzonów, zliczenia dla egzonów
b) błąd średniokwadratowy dla egzonów do pliku tekstowego w formacie: egzon, wartosc ekspresji, suma ekspresji transkryptow, blad, lista transkryptow

Ze względu na opóźnienie pojawienia się zadania – proponuję termin oddania projektu na 15 VI.

ONA – Różniczkowanie i całkowanie numeryczne

Dziś będziemy zajmować się różniczkowaniem i całkowaniem numerycznym. Slajdy są tu: ONA10-calkowanie

Przydatny będzie nam przede wszystkim pakiet scipy.integrate ale będziemy też korzystać z nowych funkcji pakietu matplotlib, takich jak wykresy strzałkowe

0. Na początek proponuje obejrzeć sobie program opisujący model Lotki-Volterry ( lotka_volterra ) i spróbować wykonać go po kawałku ze zrozumieniem. Jeśli są pytania, to proszę pytać

1. Naszym właściwym zadaniem jest budowa podobnego modelu dla problemu inwazji Zombie (opisanego w artykule naukowym ). Mamy tu 3 zmienne: ludzi niezarażonych (S), zombie (Z) oraz ofiary konfliktu (R). System jest dynamiczny, a jego zmiany opisują równania różniczkowe:

  • dS/dt = P – B*S*Z -D*S,
  • dZ/dt = B*S*Z +G*R – A*S*Z,
  • dR/dt = D*S + A*S*Z – G*R,

gdzie

  • P oznacza przyrost naturalny ludzi,
  • B – współczynnik z którym zombie spotykając niezarażonego może go zarazić,
  • D – współczynnik śmiertelności ludzi z przyczyn naturalnych,
  • G – prawdopodobieństwo, że zombie powstanie z martwej ofiary,
  • A – szansa na to, że w pojedynku człowieka z zombie, to zombie zostanie pokonane

Interesują nas następujące pytania:

a) Jak wygląda ewolucja systemu w czasie (w zakresie 0..5000),  jeśli na początku jest tylko 1000 ludzi niezarażonych (nie ma zombie, ani ciał) i mamy następujące prametry:

  • P=0.0001
  • D=0.0001
  • B=0.0095
  • G=0.0001
  • A=0.0001

Rozwiązanie przedstaw na wykresie 3 zmiennych (S,Z i R) od czasu

b) Czy zagładę populacji można zatrzymać poprzez zwiększanie współczynnika urodzin P? Rozwiązanie uzasadnij. Który ze współczynników musiałby wynosić 0, żeby populacja mogła poradzić sobie z zombie?

c) Zaprezentuj pole wektorowe populacji zombie i ludzi względem parametrów A i B?

d*) Zaprezentuj pole wektorowe 3d dla zmiennych S,Z,R

ONA 2017 – Wartości i wektory własne

Dziś poznaliśmy numeryczne rozwiązania problemu znajdowania dominujacych wartości własnych metodą potęgową i wartości własnych bliskich zadanej wartości metodą odwrotną potęgową. Dużo więcej materiału można znaleźć na tej stronie

Jeśli chodzi o funkcje przydatne w pythonie, to przede wszystkim interesują nas funkcje:

scipy.linalg.eig(A)

oraz

scipy.linalg.eigvals(A)

Zadania na dziś:

1. Zaimplementuj metodę potęgową z wykładu i zastosuj ją dla rozważanej na wykładzie macierzy: M=array([[1,2.],[4,3]]). Czy Twoje rozwiązanie różni się od rozwiązania z funkcji eig? Jakie równanie spełniają wektory własne zwrócone przez tę funkcję?

2. Zaimplementuj metodę Rayleigha RQI (z wykładu ). zastosuje ją do znalezienia najmniejszej wartości własnej macierzy M z zad. 1.

3. Rozważ macierz A=array([[-1,2,2],[2,2,-1.],[2,-1,2]]). Spróbuj znaleźć jej wartości własne zaimplementowanymi przez siebie metodami. Skąd biorą się problemy? Czy podobne problemy spotkają nas dla macierzy rand(3,3)?

4. Rozważ macierz M=array([[a,1.],[0,b]]). Gdzie a i b są dodatnimi liczbami całkowitymi. Jak zmienią się wartości własne, gdy zamiast zera wstawimy do M[1,0] wartośći 1/k dla k rosnącego wykładniczo?

ONA 2017 – Zadanie zaliczeniowe 2 – Interpolacja splajnami

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

ONA 2017 – 8 Interpolacja funkcji

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:

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

ONA 2017 7 – Zagadnienie najmniejszych kwadratów

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?