Author: bartek

Slajdy do wykładu można znaleźć tu: ONA1-Arytmetyka

Kiedy tylko LabKom uruchomi dla nas kursy w moodle’u, przeniesiemy się tam z materiałami do tego przedmiotu. Będę Państwa informował o tym e-mailem.

Pracujemy w narzędziu python3 (nie w octave) i wykonujemy zadania częściowo inspirowane zadaniami z labu 3 do metod numerycznych. Zachęcam do użycia interpretera ipython3, jeśli pracują Państwo na domowym komputerze. Jeśli ktoś ma kłopoty z interpreterem w domu, to może korzystać z serwera jupyter, używając hasła  podanego na wykładzie. Proszę wtedy pracować w osobnych podpisanych  notebookach. Nieco więcej technicznych/matematycznych szczegółów jest też dostępne w materiałach do wykładu z metod numerycznych.

Będziemy też korzystać z modułu numpy (pisząc “import numpy as np”)  a w szczególności z jego typów danych float64 i float128

Uwaga nt. liczb całkowitych. Obecnie python3 ukrywa zupełnie przed nami fakt istnienia ograniczeń wielkości rejestrów przy operacjach na liczbach całkowitych. Aby zobaczyć jak to działa, warto uruchomić interpreter python2.7 i wykonać operacje na dużych liczbach.

 

W szczególności:

1. Policz epsilon maszynowy (zarówno jako wartość dziesiętną ale też liczbe bitów cechy binarnej przy pomocy sprawdzania dla jakich n 1.0+2.0**(-n) >1.0. Co by się stało, gdybyśmy wykonali tę próbę w okolicy 0.0 zamiast 1.0. Czy taki epsilon maszynowy zależy od liczby bitów cechy, czy mantysy?

2. Oblicz wartość funkcji (pochodzącej od Siegfrieda Rump’a): f(x,y)= 333.75*y**6+x**2*(11*x**2*y**2-y**6-121*y**4-2)+5.5*y**8+x/(y*2) w punkcie x=77617, y=33096. Wykonaj obliczenia w pythonie2, pythonie3, np.float64, np.float32, np.float128 i module Fractions. Jakie błędy obliczeń tam zachodzą? Jeśli ktoś napisze program obliczający tę wartość prawidłowo, i prześle mi przez moodle program i prawidłowy wynik, to dostanie dodatkowy (bonusowy) punkt do zaliczenia.

3(*). Zastanów się nad zadaniem z szeregami zbieżnymi (lub rozbieżnymi) numerycznie ( w materiałach do labu 3). Zaimplementuj liczenie sum S(n) w pythonie, dla dowolnego n i dla kilku szeregów podanych w rozwiązaniu. Które z nich można dokładnie policzyć pry użyciu modułu Fractions?

4(*). Policz wartości e(x) wg rozwinięcia wzoru Taylora dla zadanego maksymalnego stopnia. Porównaj różnicę wyniku jaki otrzymujesz dla x=-90 i x=2.

5(*). Policz sumy częściowe szeregu harmonicznego przy użyciu liczb zmiennoprzecinkowych i przy użyciu modułu Fractions. Porównaj wyniki.

 

W naszym ostatnim i zarazem największym zadaniu zaliczeniowym odejdziemy od tematyki epidemiologicznej i zajmiemy się przetwarzaniem obrazów i ich kompresją. Konkretnie, naszym zadaniem będzie napisanie programu pozwalającego na “ukrywanie” informacji w obrazach, czyli steganografię. Wykorzystamy najprostszą metodę steganografii, czyli wykorzystanie modyfikacji najmniej znaczącego bitu do zakodowania “ukrytej informacji”.

Zasada działania:

Załóżmy, że mamy do przekazania komuś wektor bitów S, długości N  (może to być obraz, lub tekst – jest to dla nas w tej chwili bez znaczenia). Możemy wziąć jakiś obraz reprezentowany w standardowym formacie RGB, z 8 bitami na kanał i powiedzmy wymiarach 2000×1000 pikseli. Obraz ten, możemy wczytać do pamięci w postaci macierzy M o wymiarach 2000x1000x3 z typem danych int8uint8. Zauważmy, że każda z liczb w tej macierzy zawiera intensywność jednego z kolorów podstawowych w  postaci wartości od 0 do 255. Niewielkie zmiany tych wartości (np. o 1) nie są łatwo postrzegalne dla ludzkich oczu, ale program komputerowy może bardzo łatwo je wykryć.

Nasze zadanie będzie polegało na tym, że każdy kolejny bit z wektora S “wstawiać” jako ostatni bit kolejnego bajtu macierzy M. Jak go wstawiać? otóż wystarczy, że sprawdzimy czy reszta z dzielenia kolejnego bajtu macierzy M przez 2 jest równa wartości odpowiedniego bitu z wektora S.

M.reshape((2000*1000*3,))[i]%2 == S[i]

Jeśli tak, to nic nie musimy robić. Jeśli nie, to musimy odjąć lub dodać 1 do M.reshape((2000*1000*3,))[i], żeby nasze równanie było prawdziwe (uważając przy tym aby nie odjąć od zera ani nie dodawać do 255).

W ten sposób otrzymujemy nową macierz M, która ma tę własność, że ostatnie bity reprezentacji kolejnych bajtów M są równe kolejnym bitom S. Taki plik możemy zapisać jako obraz i będzie on bardzo subtelnie wizualnie zmieniony, natomiast odbiorca tego pliku, który wie czego szukać może bardzo łatwo odczytać wektor S

Aby rozwiązanie było praktyczne, musimy jeszcze być w stanie podać długość wektora S oraz typ danych do zapisu tego wektora. W naszym rozwiązaniu, pierwszy zakodowany bit (M[0]%2) będzie określał, czy mamy do czynienia z tekstem w alfabecie DNA (M[0]%2==0) czy z obrazem w skali szarości (M[0]%2==1). Kolejne 32 bity będą podawać nam długość ciągu DNA (w postaci jednej 32 bitowej całkowitej liczby dodatniej w kodzie dwójkowym) lub wymiary obrazka (w postaci dwóch 16bitowych całkowitych liczb dodatnich oznaczających odpowiednio szerokość x i wysokość y zakodowanego obrazka).

Zadanie: 

Napisz moduł pythonowy (z komentarzami!) o nazwie stegano.py,  który realizuje następujące funkcje:

1. read_image_greyscale (image_file) – zwracający wymiary obrazka oraz ciąg bitów wartości na podstawie pliku w formacie png (2 pkt)
2. write_image_greyscale (x,y,bit_vector,output_image_file) – zapisujący obrazek do formatu png w skali szarości (2pkt)
3. read_DNA(text_file) – zwracający ciąg bitów na podstawie pliku tekstowego zawierającego wyłącznie alfabet ACGT (kodowanie bitowe: A=00 ,C=01, G=10, T=11) (2 pkt)
4. write_DNA(bit_vector,output_text_file) – zapisuje ciąg liter DNA do pliku tekstowego na podstawie ciągu bitów wg powyższego kodowania (2 pkt)
5. encode_DNA(bit_vector,image_file, output_file) – wczytujący obraz RGB z pliku image file i zakodowujący wektor bit_vector reprezentujący DNA w ostatnich bitach obrazu (zgodnie z powyższym opisem, łącznie z 33 bitami nagłówka). Zmodyfikowany obraz powinien być zapisany do pliku output_file. Jeśli obraz jest za mały aby zakodować cały wektor, program powinien zwrócić błąd. (5 pkt)
6. encode_image(bit_vector,x,y,image_file,output_file) – wczytujący obraz RGB z pliku image file i zakodowujący wektor bit_vector reprezentujący obraz w ostatnich bitach obrazu (zgodnie z powyższym opisem, łącznie z 33 bitami nagłówka). Zmodyfikowany obraz powinien być zapisany do pliku output_file. Jeśli obraz jest za mały aby zakodować cały wektor, program powinien zwrócić błąd. (5 pkt)
7. decode(image_file, output_file) – powinien wczytywać obraz RGB z pliku, odczytywać z nagłówka (33 ostatnie bity z pierwszych 33 bajtów) informację o tym, czy zakodowane jest DNA, czy obraz, następnie odkodowywać informację i zapisać (korzystając z write_DNA lub write_image_grayscale) wynik do pliku output_file. (5 pkt)
8. main(), która powinna umożliwiać zakodowywanie i odkodowywanie informacji przy pomocy parametrów z linii komend obsługiwanych przy pomocy modułu Argparse. (2 pkt)

Terminy:

Rozwiązania (w postaci e-maila z normalnym załącznikami w postaci pliku stegano.py  oraz zakodowanych przykładowych danych i dopiskiem [ONA-3] w temacie) w ciągu 3 tygodni od dzisiaj, czyli do 18. VI 2020

Przykładowe dane:

Można popróbować np. zakodować taki obrazek:

lub kod DNA koronawirusa covid-19.fasta

W jakimś ulubionym pliku ze zdjęciem przedstawiającym swoją fizjonomię.

Wersja bonusowa (+5 pkt)

Można wyposażyć nasz program w dodatkową funkcjonalność polegającą na kompresji ciągu bitów przed zakodowaniem. Jeśli ktoś zaimplementuje kompresję (i dekompresję) metodą LZ77 to może otrzymać dodatkowe 5 punktów).

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 i 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

Materiały do powtórki przed kolokwium: ONA_powtorka

Kwestie organizacyjne – Na ostatnim wykładzie za dwa tygodnie  (4. czerwca) będziemy mieli kolokwium (przeprowadzone online, podczas wykładu). Będzie to test, podobny do wcześniejszych: kolokwium-2017, z tym, że przeprowadzony przy pomocy systemu online. Podobnie zrobimy z egzaminem, w czasie sesji.

Ostatnie, trzecie, zadanie zaliczeniowe opublikuję dopiero za tydzień, ale będą Państwo mieli nadal 3 tygodnie na wykonanie go, czyli do 18. czerwca.

Dzisiejszy wykład porównujący cechy różnych pakietów do obliczeń jest tu: ONA12-Obliczenia-symboliczne-i-pakiety

Zadania na dziś są dość proste i dotyczą przede wszystkim obliczeń symbolicznych i pakietu sympy (można to robić w platformie jupyter)

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: ONA11-calkowanie ONA11-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 głównym zadaniem zwykle było zadanie o zombie (tutaj wersja zeszłoroczna). W związku z okolicznościami przyrody, może lepiej zajmiemy się modelem SEIR ewolucji epidemii. Jest to model opisany równaniami:

Gdzie:

• S to liczba osób niezarażonych,
• E to liczba osób zarażonych, ale jeszcze nie zarażających,
• I to liczba osób zarażonych i zarażających
• R to liczba osób, które już ozdrowiały i nie są podatne na zarażenie

Model ma następujące parametry:

• \beta to współczynnik “skuteczności zarażania”
• \sigma to współczynnik “inkubacji” im większy, tym krótsza inkubacja
• \gamma to współczynnik “zdrowienia” im większy, tym krótszy czas “zarażania”
• \nu to współczynnik “zyskiwania/utraty odporności”, jeśli jest dodatni, część populacji zyskuje odporność bez chorowania, jeśli jest ujemny, część populacji traci odporność po ozdrowieniu
• \mu to współczynnik normujący związany ze śmiertelnością populacji. Ten model zakłada stałe N w czasie, ale zawiera przepływy normujące \mu, które związane są z efektami śmiertelności na populacje

Interesują nas następujące pytania:

a) Jak wygląda ewolucja systemu w czasie (w zakresie 0..50),  jeśli na początku jest tylko 1000 ludzi niezarażonych, i jeden zarażony  i mamy następujące prametry:

• \beta=0.9
• \gamma=0.2
• \sigma=0.5
• \mu=0
• \nu=0

Rozwiązanie przedstaw na wykresie 4 zmiennych (S, E, I, R) od czasu

b) Jak na sytuację wpłyną zmiany współczynnika \nu? Co jeśli \nu jest ujemne, jak w przypadku grypy? Co jeśli \nu jest dodatnie, co oznaczałoby spontaniczne zyskiwanie odporności, bez fazy zarażania?

c) Zaprezentuj pole wektorowe populacji S i I względem parametrów \beta i \gamma?

d*) Zaprezentuj pole wektorowe 3d dla zmiennych SEI

Praca domowa (za 1 pkt)

Trzeba zaimplementować model z niezerową śmiertelnością (powiedzmy \mu=0.01), wykonać symulacje i opisać jednym zdaniem co się zmienia jeśli chodzi o liczbę chorych infekujących (I) w stosunku modelu z \mu=0.

Żeby dostać ten punkt, trzeba wysłać do mnie jako załącznik plik .py generujący wykresy (SEIR od czasu) dla dwóch modeli (\mu=0 i \mu=0.01) a w treści wiadomości napisać co zaobserwowaliśmy. Czas, jak zwykle, do następnych zajęć.