Category: ONA

Slajdy do wykładu są tutaj: ONA8-DSP

Interesują nas przede wszystkim biblioteki scipy.signal i scipy.fftpack

1. Wygeneruj sygnał funkcji wielomianowej o zaburzeniu Gaussowskim, dla 10000 punktów. Wykonaj uśrednienia tego sygnału przy pomocy średniej kroczącej, albo splotu przy pomocy filtru kwadratowego, trójkątnego lub Gaussowskiego. WYniki zaprezentuj na wykresie

2. Pobierz sygnał o zajętości nukleosomowej z plików pochodzących z eksperymentu MNase-Seq (tutaj). Wczytaj go do numpy jako wektor. Przedstaw go na wykresie. Używając szybkiej trasformaty fouriera oblicz widmo Fourierowskie tego sygnału. Przedstaw je na innym wykresie. Wyzeruj część widma wysokich częstotliwości i dokonaj odwrotnej transformaty Fouriera, aby uzyskać wygładzony sygnał zajętości nukleosomami. Spróbuj tak dobrać parametry filtra, aby wykres po odwrotnej transformacie Fouriera miał okres zbliżony do 160-200 par zasad.

 

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?

 

Dzisiaj zajmiemy się Interpolacją funkcji przy pomocy wielomianów i funkcji sklejanych.

Teoria z wykłądu 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?

 

 

 

Naszym zadaniem będzie napisanie programu generującego losowe krajobrazy fraktalne.

Losowy krajobraz jest płaszczyzną dwuwymiarową w przestrzeni trójwymiarowej. Zwkle rozważa się płąszczyzny trójkątne lub kwadratowe, my zajmiemy się tym drugim przypadkiem. Płaszczyznę będziemy reprezentować jako macierz wymiaru 2^N+1*2^N+1, gdzie M[i,j] oznacza wysokość (współrzędną Y) punktu o x=i i z=j.

Aby utworzyć losowy krajobraz tego rodzaju, będziemy postępować wg następującego rekurencyjnego schematu:

Zacznijmy od zainicjowania wartości narożnych macierzy (krok 0).

Następnie wykonujemy N zejść rekurencyjnych (k przebiega od N do 1):

• (diamond step na schemacie) ustaw wartość każdego punktu środkowego (żółty) na wartość średnią pomiędzy czterema narożnikami (czarnymi)zaburzoną o wartość 2*k*sigma* norm()
• (square step na schemacie) Ustaw wartości czterech wartości na środkach (żółtych) brzegów kwadratu na wartości średnie pomiędzy ich czterema sąsiadami (czarnymi). zaburzoną o ((2 *k)-1)*sigma*norm()

Naszym zadaniem jest nie tylko generowanie takiego krajobrazu, ale także jego wizualizacja w postaci mapy fizycznej (imshow) i płaszczyzny 3D.

Nasz program w wersji podstawowej (10pkt) powinien dać się uruchomić jako skrypt i pozwalać na ustawienie następujących parametrów:

• N – rząd macierzy (2^N+1). Domyślnie 10, powinno działać co najmniej do N=3014
• sigma – stopień “górzystości” terenu
• map_file – nazwa pliku (pdf lub png) z mapą. Jeśli podane, to mapa jest zapisywana do pliku, jeśli nie, otwieramy okienko z wykresem
• surf_file – nazwa pliku (pdf lub png), z wykresem płaszczyzny. Jeśli podane, to powierzchnia jest zapisywana do pliku, jeśli nie, otwieramy okienko z wykresem
• colormap – mapa kolorów z matplotlib do użycia
• matrix_file – zapis macierzy do pliku

Wersja nieco bardziej zaawansowana (dodatkowe 2 pkt) może obsługiwać także następujące opcje:

• wczytywanie macierzy już częśćiowo wypełnionej i losowanie tylko od kroku K
• Okresowe warunki brzegowe, tak, aby dało się “wykafelkować” większą powierzchnię

Rozwiązania wysyłamy do 3. maja 2016 (włącznie) w postaci e-mail’a do wykładowcy z [ONA-1-2016] w tytule.

Rozwiązania po terminie otrzymują połowę punktów, na które zasługiwałyby w terminie i można je przysyłać do 6 VI 2016).

Dziś zajmujemy się ukłądami 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. 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?
2. 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()
3. 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()?
4. 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

Dzisiaj zajmujemy się wykresami.

Slajdy są tu ONA3-wykresy

Zadania na dziś:

1. Narysuj wykres podobny do tego ze slajdu nr. 14  z wykładu. Uwzględnij opisy osi, zakresy wartości, kolory, znaczniki wartości, legendę i tytuły osi oraz wykresów.

2. Narysuj wykres 10*sin(x)+normal(0,1) dla 1000 punktów w zakresie 0,10. Przyda się funkcja normal().

3. Napisz program, który wykona 1000 prób rzutu 10 kostkami. Narysuj na jednym wykresie 3 przenikające się histogramy (alpha=0.25): suma oczek na kościach parzystych (0,2,4,6,8), suma oczek na kościach nieparzystych (1,3,5,7,9), suma oczek na wszystkich kościach. Może przydać się funkcja randint().

4. Przedstaw dane z wykresu 3. w postaci 3 wykresów pudełkowych (boxplot)

5. Narysuj powierzchnię 3d (x,y,z) in linspace(0,1,100) gdzie z=y*cos(x)+(1-y)*sin(x), Może być przydatny ten przykład

5. Narysuj mapę ciepła dla tej samej powierzchni. Użyj różnych map kolorów.

 

Slajdy do wykładu są tu: ONA2-Macierze

Ćwiczenia:

1. Stwórz macierze A i B o wymiarach  1000×1000 zawierające wartości A[i,j]==i*3-j*5 i [i,j]==np.sqrt(A[i,j]) (Zauważ, że to wymaga wartości urojonych, a więc macierzy typu complex64 i wartość sqrt(-1) ->.j)
2. Napisz funkcję poz(indeksy, shape), która dla zadanej krotki indeksów zwraca liniową pozycję elementu o zadanych indeksach w wielowymiarowaj macierzy o kształcie shape.
3. Korzystając z broadcastingu i mnożenia macierzy zmień znak elementów w macierzy B, których suma indeksów jest nieparzysta
4. korzystając z sortowania , posortuj macierz 3-wymiarową wg. drugiej współrzędnej i obejrzyj wynik.
5. Korzystając z funkcji frompyfunc() napisz funkcje wektorowe, które liczą a) sumę kwadratów dwóch macierzy  (dwie macierze do jednej)
   b) iloraz i resztę z dzielenia całkowitoliczbowego przez 17 (jedna macierz wejściowa i dwie macierze wyjściowe)
6. Korzystając z modułu time, porównaj prędkość funkcji liczących N potęg dwójki – na liście i w wektorze. Dla jakich wartości N warto używać macierzy?

Tak jak wspominałem w zeszłym tygodniu – w związku z moim  wyjazdem do Berkeley w dzisiejszym terminie zajęcia się nie odbędą. Myślę, że pod koniec semestru zorganizujemy dodatkowe zajęcia, kiedy będziecie Państwo już realizować projekty.

 

Dzisiaj zajmujemy się arytmetyką komputerów.  Slajdy z naszego wykładu są tutaj (pdf).

Pracujemy w narzędziu ipython3 (nie w octave) i wykonujemy zadania z labu 3 metod numerycznych. Nieco więcej technicznych szczegółów jest też dostępne w materiałach do wykładu z metod numerycznych.

W szczególności:

Zastanów się nad zadaniem z szeregami zbieżnymi (lub rozbieżnymi) numerycznie. 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?

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.

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

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

Napisz funkcję znajdującą miejsca zerowe dla dowolnej funkcji w przedziale [a,b] ( def zero(f,a,b):) metodą bisekcji.  Sprawdź jak ona działa dla różnych równań typu 1/x-c.