Author: bartek

W związku z tym, że na ostatnich ćwiczeniach  wynegocjowali Państwo ze mną przesunięcie terminu prezentacji projektów bonusowych o tydzień, to dziś nie ma normalnych zajęć.

Jeśli ktoś ma ochotę skonsultować się w sprawie zadań (tych jeszcze nie wykonanych) albo tematów egzaminacyjnych, zapraszam w godzinach zajęć (10-14) do mnie do 5770.

Prezerntacje projektów, które były planowane na dziś, odbędą się 21. 6 w godzinach 12-15,

Slajdy do dzisiejszego wykładu są tu:ONA12_powtorka,  Znajdą tam Państwo tematy  do opanowania na egzamin pisemny.

Ustaliliśmy, że egzamin odbędzie się w postaci testowej, bez możliwości korzystania z materiałów, dnia 27. czerwca, o godzinie 10:00.

 

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:

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

 

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

 

Naszym zadaniem jest stworzenie pakietu umożliwiającego obróbkę macierzy kontaktów chromosomowych. Będzie on obejmował 3 zasadnicze zadania:

1. Wczytywanie listy kontaktów dla wybranego chromosomu do macierzy o zadanej rozdzielczości (5pkt). Mamy dany plik wejściowy, który zawiera 4 kolumny:

chromosom1 pozycja1 chromosom2 pozycja2

Każdy taki wpis (a będą ich miliony) opisuje kontakt pojedyńczej nici DNA pochodzącej z chromosomu1 na pozycji1 z nicią pochodzącą z chromosomu2 na pozycji2. Mając do dyspozycji taki plik oraz podany numer chromosomu C i długość segmentu S powinniśmy skonstruować kwadratową macierz symetryczną M kontaktów (rozmiaru N*N, gdzie N= L/S, gdzie L oznacza długość chromosomu C) pomiędzy segmentami zadanego chromosomu. Np. jeśli chromosom ma długość L=1000000 par zasad, a S=1000, otrzymamy macierz M rozmiaru 1000*1000, a wartość M[i,j] będzie oznaczać ile wierszy opisywało kontakt pomiędzy fragmentem [1000*i:1000*(i+1)] a fragmentem [1000*j:1000*(j+1)].

2. Normalizacja macierzy (10pkt). Mając daną macierz kwadratową M o zadanym rozmiarze N*N, często potrzebujemy dokonać jej normalizacji. W tym przypadku interesuje nas rozwiązanie równania M[i,j]= B[i]*[B[j]*T[i,j]+epsilon, gdzie wektor B (rozmiaru N) opisuje multiplikatywne obciążenia poszczególnych segmentów, zaś od Macierzy T (rozmiaru N*N) oczekujemy, że będzie symetryczna i znormalizowana (sum(T[:,i])=1.0 dla każdego i). Daje nam to układ równać, który wymaga rozwiązania iteracyjnego, minimalizującego błąd (epsilon). Wynikiem tej procedury powinna być macierz T, wektor B i wartość epsilon, a dozwolonym parametrem powinna być maksymalna liczba iteracji i maksymalny Epsilon.

3. Wyświetlanie  macierzy kontaktów (5 pkt). Kiedy już mamy macierz kontaktów (M lub T) to chcielibyśmy móc obejrzeć ją na ekranie lub zapisać do pliku. Potrzebny nam będzie program, który na podstawie macierzy i dodatkowych parametrów wyświetla zadaną macierz metodą imshow(). Program powinien pozwalać na następujące opcje:

• określanie podmacierzy do wyświetlenia (tylko prostokąt o współrzędnych x_1:x_2,y_1:y_2)
• Wygładzanie mapy przy pomocy filtra Gaussowskiego o określonym promieniu
• logarytmiczna lub liniowa mapa kolorów
• wybór schematu kolorów
• Zmniejszanie rozdzielczości dużych map do wizualizacji o stały czynnik (jeśli mapa ma np. 10000×10000 można zadać zmniejszenie jej przed wizualizacją – np. zmniejszenie o czynnik 10 powinno uśrednić dane z kwadratów rozmiaru 10×10 i wyświetlić mapę rozmiaru 1000×1000)

Wynik wysyłamy w postaci pojedyńczego skryptu (jako niespakowany załącznik) na adres prowadzącego zajęcia z dopiskiem ONA-3-2016 w tytule do dnia 20. VI 2016.

UWAGA: Przykłądowy plik z danymi (dla 1 chromosomu o długości 1mln par zasad) można znaleźć tutaj – uwaga, plik jest duży…)

UWAGA 2: Kontakty chromosomowe są dość często wzbogacone w kontakty blisko “przekątnej”, czyli kontakty niedużego zasięgu. W związku z tym warto rozważyć wersję programu, która ignoruje wartości w macierzy na samej przekątnej.

Wskazówki do rozwiązania po pierwszych pytaniach:

1) Czy dla dużych macierzy M (np. 1000×1000) możemy myśleć o rozwiązaniu dokłądnym tego układu równań?

2) Czy normalizacja kolumnowo/wierszowa macierzy T nasuwa nam jakieś oszacowanie na przybliżone wartości B[i]?

3) Czy jeśli ktoś dałby nam wartość przybliżoną B, to umielibyśmy znaleźć krok iteracji poprawiający to oszacowanie przy pomocy macierzy M?

4) Epsilon, o którym wspominam w zadaniu jest pewnego rodzaju skrótem myślowym. Kiedy występuje w równaniu, mamy na myśli różnicę między prawą, a lewą stroną tego konkretnego równania (powiedzmy, mały epsilon). Kiedy pojawia się jako kryterium zakończenia iteracji (duży Epsilon), jest to już różnica między wszystkimi lewymi i prawymi stronami równań. W zagadnieniu najmniejszych kwadratów, typowo ten duży Epsilon to suma kwadratów małych epsilonów. W naszym przypadku, gdzie liczba równań, a więc i składników sumy zależy istotnie od parametrów zadanych przez użytkownika, warto rozważyć duży Epsilon będący maksimum z wartości bezwzględnych małych epsilonów.