Category Archives: ONA

Obliczenia naukowe

ONA 13 – Metody Monte Carlo i rozwiązania równań 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 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

ONA 12 – Obliczenia symboliczne i pakiety do obliczeń

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

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

Tymczasem, za tydzień, na wykładzie mamy kolokwium. Zagadnienia do powtórki są tutaj: ONA_powtorka

Przykładowe kolokwium jest tu: kolokwium-2017 ( w tym roku dodatkowe są zagadnienia z kompresji).

Dodatkowe konsultacje przed kolokwium zrobimy w czwartek o 10:15 w sali 2041 (lab)

Egzamin będzie w poniedziałek 25. VI o godzinie 10tej.

Egzamin poprawkowy, 14. IX, też o 10tej

ONA – Zadanie 3

Naszym zadaniem będzie napisanie programu, który będzie wykonywał kompresję i testował jej zachowanie dla różnych dugości pliku i parametrów kompresji. Kolejne elementy podzadania to:

1) Zaimplementuj prosty algorytm kompresji oparty o metodę Lempela i Ziva (LZ77) omawianą na wykładzie. Powinna to być funkcja w pythonie lz77_compress(napis, rozmiar_słownika ), która pobiera napis w argumencie, a następnie zwraca jego wersję skompresowaną w postaci napisu. Zaimplementuj także funkcję odwrotną lz77_decompress(napis)>  (5 punktów)

2) Napisz dwie funkcje generujące napisy o zadanej długości. gen_losowy(dlugosc, zestaw_znakow), ktora generuje losowy ciąg znaków pochodzących z zadanego zestawu znaków. Zarówno zestaw dozwolonych znaków, jak i wynik jest w formie napisu. Drugi z generatorów powinien być mniej losowy. Najlepiej, gdyby losował kolejne znaki, ale za każdym razem wstawiał ciąg po 100 identycznych znaków. (5 punktów)

3) Napisz funkcję testuj(kompresor,generator, lista_długości), która używa podanego generatora do wygenerowania zestawu napisów o długościach zadanych w parametrze lista_dlugosci. Następnie używa podanego w pierwszym argumencie kompresora, aby dokonać kompresji wszystkich tych napisów. Ostatecznie, funkcja zwracałaby listę długości wszystkich skopmpresowanych napisów. (5 punktów)

4) Użyj napisanych funkcji, oraz swojej wiedzy o aproksymacji funkcji do tego aby napisać program pozwalający  znaleźć funkcje opisujące:

  • zależność rozmiaru skompresowanego od wejściowego jako funkcji liniowej lub kwadratowej dla obu generatorów
  • zależność współczynnika kompresji (wyestymowanego w poprzednim punkcie) od rozmiaru słownika

Obraną metodę i wyniki opisz nie tylko w komentarzach w pliku źródłowym, ale także w postaci krótkiego raportu (maksymalnie 2 strony), ilustrowanego sensowymi wykresami (10 pktów)

Rozwiązania ( w postaci 1 pliku .py i 1 pliku .pdf) przesyłamy do 3. czerwca, e-mailem do wykładowcy, proszę pamiętać o dopisku [ONA-2018-3] i o umieszczeniu plików jako standardowych załączników.

ONA 11 – Różniczkowanie i całkowanie numeryczne

Dziś będziemy zajmować się różniczkowaniem i całkowaniem numerycznym. Slajdy są tu: 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 właściwym zadaniem jest budowa podobnego modelu dla problemu inwazji Zombie (opisanego w artykule naukowym , kopia lokalna tu: Zombies). 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 10 – 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 9 – Interpolacja funkcji

Tym razem 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?

 

ONA 8 – Metoda 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ń liniowych, 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?

 

ONA 7 – Układy równań liniowych

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)

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

  1. Stwórz macierz Hilberta H[i,j]=1/i+j-1 dla sensownie dużego rozmiaru macierzy (100×100, 200×200) i wyświetl ją przy pomocy imshow()/
  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.

ONA – Zadanie zaliczeniowe 2

Naszym zadaniem będzienapisanie programu, który będzie przetwarzał strumienie danych genomicznych zapisane w plikach sgr . Są to bardzo proste pliki, które pozwalają opisywać zmienność parametru wzdłuż chromosomów. Każda linijka ma tylko 3 wartości:

chromosom pozycja wartość

Wartości są oddzielane znakiem tabluacji (“\t”), pozycja jest liczbą całkowitą dodatnią, a wartość jest liczbą zmiennoprzecinkową.

Zakładamy, że wartości pomiędzy punktami w pliku są odcinkami prostoliniowymi, co pozwala nam na podstawie pliku sgr, wyliczyć wartości dla wszystkich pozycji chromosomu.

Przykładowe pliki sgr dla dwóch sygnałów i dwóch filtrów są w pliku Zad2

Pliki te mogą być duże – np. większe niż dostępna pamięć, co oznacza, że nie możemy ich wczytać w całości do wektora.

Będziemy chcieli, żeby nasz program oferował następujące funkcje:

- Suma, różnica, iloczyn, iloraz – zwróć wynik operacji arytmetycznej na  dwóch sygnałach dla tych samych chromosomów. Zakładamy,  że na wejściu mamy dwa pliki, które opisują sygnał wzdłuż tych samych chromosomów, choć niekoniecznie w tych samych pozycjach. OPeracje są wykonywane “po pozycjach”. Zwracamy plik sgr, który opisuje funkcję wynikową (4 pkt)

- wygładzanie – przy pomocy średniej kroczącej o zadanej długości (w parach zasad), zwróć sygnał wygładzony (2 pkt)

- splot – Tym razem na wejściu bierzemy jeden plik sgr z sygnałem (wzdłuż chromosomów) i drugi plik sgr z wartościami funkcji filtra (np. filtr prostokątny, lub Gaussowski), wzdłuż sztucznego chromosomu o nazwie “filtr”. Wynikiem powinien być splot funkcji, czyli przefiltrowany sygnał 6  pkt)

- upraszczanie – dla zadanego pliku sgr wykryj, które punkty są niepotrzebne (tzn po ich usunięciu sygnał nie zmieni się) i zwróć sygnał pozbawiony tych linijek (3 pkt).

Rozwiązanie polega na napisani 7 funkcji, każda z nich powinna działać na 2 plikach (+ parametr szerokość dla wygładzania), wszystko w jednym pliku .py

Rozwiązania podpisane, z dopiskiem [ONA-2018-2] w temacie proszę wysyáć mailem do mnie, do 22. kwietnia 2018.

ONA 6 – Kompresja

Slajdy są tutaj ONA6-kompresja

Zadania na dziś:

1. Zapoznaj się z modułem gzip  zapisz plik tekstowy w pliku skompresowanym ze skryptu pythona, zdekompresuj go przy pomocy programu gunzip i odwrotnie

2. Weźmy takie pomocnicze funkcje w pythonie opisujące dyskretną transformatę kosinusową:

import urllib2,io
import Image
from scipy import fftpack

image_url='http://i.imgur.com/8vuLtqi.png'

def get_image_from_url(image_url='http://i.imgur.com/8vuLtqi.png', size=(128, 128)):
    file_descriptor = urllib2.urlopen(image_url)
    image_file = io.BytesIO(file_descriptor.read())
    image = Image.open(image_file)
    img_color = image.resize(size, 1)
    img_grey = img_color.convert('L')
    img = np.array(img_grey, dtype=np.float)
    return img

def get_2D_dct(img):
    """ Get 2D Cosine Transform of Image
    """
    return fftpack.dct(fftpack.dct(img.T, norm='ortho').T, norm='ortho')

def get_2d_idct(coefficients):
    """ Get 2D Inverse Cosine Transform of Image
    """
    return fftpack.idct(fftpack.idct(coefficients.T, norm='ortho').T, norm='ortho')

def get_reconstructed_image(raw):
    img = raw.clip(0, 255)
    img = img.astype('uint8')
    img = Image.fromarray(img)
    return img

wykonaj program, który dokonuje Ntego przybliżenia obrazu, przy pomocy zerowania wartośći powyżej Ntego wiersza i Ntej kolumny macierzy get_2D_dct(img) i wyświetl kilkanaście pierwszych przybliżeń.

3. Wyrysuj wykres Entropii dla kanału binarnego w zależności od P(1)

4. Napisz program, który dla zadanego pliku tekstowego tworzy tablicę częstotliwości znaków, kody Huffmana dla wszystkich znaków oraz wylicza entropię tego pliku i długość kodu Huffmana, który opisywałby cały plik