PP1 laboratorium

tworzenie macierzy i indeksowanie elementów macierzy
operacje macierzowe i funkcje działające na macierzach
macierze logiczne i operacje na macierzach logicznych
optymalizacja kodu
wektoryzacja - zamiana pętli po elementach na mechanizmy indeksowania
prealokacja - dynamiczne modyfikowanie rozmiarów macierzy jest kosztowne

macierze można tworzyć za pomocą operatora indeksowania :
uwaga: ilość elementów w wierszach musi być jednakowa

A = [1 2 3; 4 5 6]
B = [5:8 ; 10:13]

A =

   1   2   3
   4   5   6

B =

    5    6    7    8
   10   11   12   13

rand(n,m) losowe wartość z przedziału $[0, 1]$
zeros(n,m) macierz zerowa
ones(n,m) macierz zawierająca jedynki
eye(n,m) macierz jednostkowa (jedynki na diagonali)

a = rand(3)
b = rand(3, 2)

a =

   0.306403   0.834324   0.986851
   0.031466   0.872381   0.695912
   0.942068   0.691957   0.619879

b =

   0.244435   0.975030
   0.421339   0.069198
   0.647612   0.217188

c = zeros(5)
d = ones(2,5)

c =

   0   0   0   0   0
   0   0   0   0   0
   0   0   0   0   0
   0   0   0   0   0
   0   0   0   0   0

d =

   1   1   1   1   1
   1   1   1   1   1

e = eye(4, 6)

e =

Diagonal Matrix

   1   0   0   0   0   0
   0   1   0   0   0   0
   0   0   1   0   0   0
   0   0   0   1   0   0

Ćwiczenie utwórz macierz o wymiarach 100×100, która posiada wartości 1 w pierwszej i czwartej ćwiartce oraz wartości 0 w pozostałych miejscach, zgodnie z poniższym przykładem macierzy 4×4

indeksowanie 2 wymiarów macierzy (N,M) od 1 do end
indeksowanie liniowe (N) od 1 do end po elementach kolejnych kolumn

X = [1, 2, 3; 4, 5, 6]

% indeksowanie 2-wymiarowe
a = X(1,2)           % 1-szy wiersz 2-ga kolumna
b = X(1:2, 1)
c = X(1:2, 1:2)

e = X(:, 1)          % cała kolumna
f = X(end, :)        % cały wiersz

% indeksowanie liniowe
a = X(1)
b = X(end)
c = X(1:end)    # wszystkie elementy w wektorze wierszowym
d = X(:)        # wszystkie elementy w wektorze kolumnowym

a =  1
b =  6
c =

   1   4   2   5   3   6

d =

   1
   4
   2
   5
   3
   6

length(w) ilość elementów wektora lub najdłuższy wymiar macierzy
size(m) zwraca ilość wierszy i kolumn
numel(m) ilość elementów

w = 1:100;

a = length(w)
b = size(w)
c = numel(w)

X = rand(13, 100);

a = length(X)
[w k] = size(X)
c = numel(X)

a =  100
w =  13
k =  100
c =  1300

W wielu przypadkach iteracja po elementach macierzy za pomocą pętli nie jest konieczna
Indeksowanie za pomocą :, operatory macierzowe (*, +, itd.) oraz funkcje dedykowane do obliczeń na macierzach (np. sum(X)) będą z reguły wydajniejsze od pętli

% Iterowanie po elementach macierzy
w = rand(1, 10);
for i=1:length(w)
    fprintf('Element %d wektora = %f\n', i, w(i))
end

Element 1 wektora = 0.468979
Element 2 wektora = 0.160596
Element 3 wektora = 0.254776
Element 4 wektora = 0.487126
Element 5 wektora = 0.721925
Element 6 wektora = 0.443289
Element 7 wektora = 0.345285
Element 8 wektora = 0.392169
Element 9 wektora = 0.011494
Element 10 wektora = 0.963213

% Iterowanie po elementach macierzy
m = rand(2 ,3);
[w k] = size(m);
for i=1:w
    for j=1:k
        fprintf('Element (%d,%d) macierzy = %f\n', i, j, m(i, j))
    end
end

Element (1,1) macierzy = 0.192715
Element (1,2) macierzy = 0.097834
Element (1,3) macierzy = 0.704280
Element (2,1) macierzy = 0.938989
Element (2,2) macierzy = 0.119784
Element (2,3) macierzy = 0.608151

Ćwiczene napisz funkcję suma, która wyznaczy sumę wszystkich elementów macierzy podanej w argumencie. Ćwiczenie wykonaj z użyciem instrukcji pętli, iterując po wszystkich elementach macierzy.

W jaki sposób osiągnąć to samo za pomocą funkcji matlabowej sum?

Ćwiczenie napisz skrypt/funkcję o nazwie suma2.m, która uruchamia operację sumowania z poprzedniego ćwiczenia wielokrotnie (np. n=10000) dla macierzy 100×100 o losowych wartościach i mierzy czas wykonania tej operacji. Tą samą operację powtórz za pomocą funkcji sum i porównaj czasy obliczeń z obu podejść.

Czas wykonywania operacji można wyznaczyć za pomocą poleceń tic i toc, przykład:

tic
sum(x)
toc

Elementy macierzy można również zsumować za pomocą iloczynu z macierzami jedynkowymi:

$$ \left[\begin{array}{ccc} 1 & 1 & 1 \\ \end{array}\right] \cdot\left[\begin{array}{ll} 3 & 1 \\ 2 & 1 \\ 1 & 0 \end{array}\right] \cdot\left[\begin{array}{ll} 1 \\ 1 \\ \end{array}\right] $$

Sprawdź również szybkość obliczeń sumy tą metodą.

reshape(A, M, N) lub reshape(A, [M N]) zmienia kształt macierzy na MxN (elementy układane kolumnami)
fliplr(m) , flipup(m), rot90() odbicie lewo-prawo, góra-dół i obrót 90 stopni

w = 1:12; 

a = reshape(w, 3, 4)
a1 = reshape(a, [2 6])

a =

    1    4    7   10
    2    5    8   11
    3    6    9   12

a1 =

    1    3    5    7    9   11
    2    4    6    8   10   12

b = fliplr(a)
c = flipud(a)
d = rot90(a)

b =

   10    7    4    1
   11    8    5    2
   12    9    6    3

c =

    3    6    9   12
    2    5    8   11
    1    4    7   10

d =

   10   11   12
    7    8    9
    4    5    6
    1    2    3

Ćwiczenie korzystając z funkcji randperm wygeneruj macierz 10×10 zawijającą liczby od 1 do 100 umieszczone w losowej kolejności.

m = [];    # pusta macierz (wektor)

a = size(m)
b = length(m)
c = numel(m)

m = [m 1]
size(m)

a =

   0   0

b = 0
c = 0
m =  1
ans =

   1   1

dla wektorów można usuwać dowolne elementy
dla macierzy można usuwać tylko całe wiersze lub całe kolumny

w = 1:12
m = reshape(w, 3, 4)

w(1) = []   % usuwanie elementu
w([3 6]) = []
w(6:end) = []

% m(1,1) = []
m(1, :) = []
m(:, 2) = []

w =

    1    2    3    4    5    6    7    8    9   10   11   12

m =

    1    4    7   10
    2    5    8   11
    3    6    9   12

w =

    2    3    4    5    6    7    8    9   10   11   12

w =

    2    3    5    6    8    9   10   11   12

w =

   2   3   5   6   8

m =

    2    5    8   11
    3    6    9   12

m =

    2    8   11
    3    9   12

operacje skalarne na elementach macierzy: +, -, *, /, .^, .\

A = [1 2 3; 4 5 6]

a = A + 3     # to samo co A .+ 3
b = A * 5     # to samo co A .* 5
c = A / 5     # to samo co A ./ 5

a =

   4   5   6
   7   8   9

b =

    5   10   15
   20   25   30

c =

   0.20000   0.40000   0.60000
   0.80000   1.00000   1.20000

d = A .^ 2     
e = A .\ 5    # to samo co 5 ./ A

d =

    1    4    9
   16   25   36

e =

   5.00000   2.50000   1.66667
   1.25000   1.00000   0.83333

operacje macierzowe: *, /, \
A * B iloczyn macierzy (iloczyn skalarny dla wektorów)
A / B dzielenie macierzy, wynikiem jest X takie, że X * B = A
A \ B dzielenie macierzy, wynikiem jest X takie, że A * X = B
potęgowanie A^n tylko dla macierzy A kwadratowych

A = [1 2 3; 4 5 6];
B = [1 2 3; 4 5 6];

a = A * B'  # ilość wierszy A == ilość kolumn B
b = A / B   # wymiary A i B takie same
c = A \ B   # wymiary A i B takie same

a =

   14   32
   32   77

b =

   1.0000e+00   4.7673e-16
   2.2979e-16   1.0000e+00

c =

   0.83333   0.33333  -0.16667
   0.33333   0.33333   0.33333
  -0.16667   0.33333   0.83333

d = A(1:2,1:2)^2

opearcje macierzowe element po elemencie: +, -, .*, ./, .^, .\
wymiary macierzy muszą się zgadzać

A = [ 1 2 3; 4 5 6 ]; 
B = [ 3 4 6; 3 4 3 ];

a = A + B    # to samo co A .+ B
b = A .* B
c = A ./ B

a =

   4   6   9
   7   9   9

b =

    3    8   18
   12   20   18

c =

   0.33333   0.50000   0.50000
   1.33333   1.25000   2.00000

d = A .^ B
e = A .\ B  # to samo co B ./ A

d =

     1    16   729
    64   625   216

e =

   3.00000   2.00000   2.00000
   0.75000   0.80000   0.50000

operacje macierzy i wektora element po elemencie powielane wzdłuż pasującego wymiaru

A = ones(2, 3) + 1
w1 = [1 2 3]
w2 = [1; 2]

a = A + w1    # sumowanie elementów w wierszach macierzy
b = A + w2    # sumowanie elementów w kolumnach macierzy

a = A .* w1    
b = A .* w2

Wiele funkcji Matlaba jest przystosowana do działania na macierzach o dowolnych kształtach

w = 1:6
X = reshape(w, 2, 3)

a = sin(w)
b = cos(X)

w =

   1   2   3   4   5   6

X =

   1   3   5
   2   4   6

a =

   0.84147   0.90930   0.14112  -0.75680  -0.95892  -0.27942

b =

   0.54030  -0.98999   0.28366
  -0.41615  -0.65364   0.96017

Ćwiczenie popraw program pole_kola.m tak aby funkcja zwracała poprawny wynik dla argumentu macierzowego. Wynikiem funkcji będzie macierz zawierająca w każdym z elementów $P_{ij}$ wartość pola dla promienia $R_{ij}$. Postaraj się nie korzystać z instrukcji pętli.

Ćwiczenie zmierz czas wykonywania poniższych instrukcji.

porównaj szybkość działania przy zastosowaniu prealokowanego wektora x
przyśpiesz działanie stosując wektoryzację (pozbądź się instrukcji pętli)

n=100000
tic
for i=1:n
    x(i) = sqrt(i);
end
toc
sum(x)

n =  100000
Elapsed time is 0.43374 seconds.
ans =  21082008.97392

macierze zawierające wartości logiczne (typ logical), wartość 0 to fałsz, wartość 1 (lub niezerowa) to prawda
indeksowanie macierzy za pomocą macierzy logicznych
funkcje false() i true() pozwalają tworzyć macierze logiczne
funkcja logical() zamienia macierz liczbową na logiczną

X = [1 , -5, 7 ; -1, 42, -12];
ind = X < 0
whos ind

ind =

  0  1  0
  1  0  1

Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  ===== 
        ind         2x3                          6  logical

Total is 6 elements using 6 bytes

# indeksowanie macierzą logiczną
a = X(ind)    
b = X( X > 0)
X( X > 0) = 111

a =

   -1
   -5
  -12

b =

    1
   42
    7

X =

   111    -5   111
    -1   111   -12

A = true(2)
B = false(2, 3)
whos A B

A =

  1  1
  1  1

B =

  0  0  0
  0  0  0

Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  ===== 
        A           2x2                          4  logical
        B           2x3                          6  logical

Total is 10 elements using 10 bytes

# wartości logiczne mogą być traktowane jak liczby całkowite
C = A + 1
whos A C

C =

   2   2
   2   2

Variables in the current scope:

   Attr Name        Size                     Bytes  Class
   ==== ====        ====                     =====  ===== 
        A           2x2                          4  logical
        C           2x2                         32  double

Total is 8 elements using 36 bytes

Ćwiczenie dlaczego poniższe wyrażenie nie działa? Spróbuj to naprawić.

w = 1:5;
w([1 0 1 0 1]) = -1

error: w(0): subscripts must be either integers 1 to (2^63)-1 or logicals

any() zwraca 1 gdy wektor posiada niezerowy element,
dla macierzy zwraca wektor z wartością 1 gdy w kolumnie jest niezerowa wartość

all() zwraca 1 dla każdej kolumny której wszystkie elementy są niezerowe
find() zwraca indeksy (indeksowanie liniowe) niezerowych elementów
isequal() sprawdza czy macierze mają jednakowe elementy

w1 = [ 0 0 1 2 ];
w2 = zeros(1,10);

a = all(w1)
b = any(w1)

c = all(w2)
d = any(w2)

a = 0
b = 1
c = 0
d = 0

E = [ 1 2 ; 0 0 ; 3 * eye(2)]'
all(E)
any(E)

E =

   1   0   3   0
   2   0   0   3

ans =

  1  0  0  0

ans =

  1  0  1  1

# indeksy liniowe
k = find(E)
a = E(k)

# indeksy wierszy i kolumn
[w k] = find(E)

# szukanie zerowych elementów
k2 = find(~E)

# szukanie elementów spełniających warunki
k3 = find(E > 2)

X = E
isequal(E, X)

% to samo za pomocą funkcji all()
all(all(X == E))   

X(1,1) = -X(1,1)
isequal(E, X)
all(all(X == E))

X =

   1   0   3   0
   2   0   0   3

ans = 1
ans = 1
X =

  -1   0   3   0
   2   0   0   3

ans = 0
ans = 0

operatory logiczne działające element po elemencie: & and, | or, ~ not

X = magic(5)
X(X > 5 & X < 10)

X =

   17   24    1    8   15
   23    5    7   14   16
    4    6   13   20   22
   10   12   19   21    3
   11   18   25    2    9

ans =

   6
   7
   8
   9

Ćwiczenie utwórz macierz X o wymiarach 10×10 posiadającą losowe wartości z zakresu od -1 do 1. Następnie zastąp wszystkie elementy o wartości ujemnej w macierzy X liczbą 0.

Napisz funkcję o nazwie mc_pi, która wyznaczy przybliżenie wartości liczby $\pi$ za pomocą całkowania metodą Monte Carlo. Zadanie sprowadza się do wyznaczenia pola koła o promieniu 1.

$$ P_{\circ} = \pi r^2 = \pi $$

Pole ćwierci koła można wyznaczyć licząc stosunek liczby $K$ losowych punktów wewnątrz ćwierci koła o promieniu 1 do liczby $N$ ilości wszystkich wylosowanych punktów w kwadracie o polu 1.

$$4 \cdot \frac{K}{N} \approx \pi $$

Algorytm:

Wylosuj $N$ losowych punktów $(x, y)$, gdzie $x, y \in [0, 1]$
Wyznacz liczbę punktów $K$, która trafiła w tarczę koła, tzn. spełniające warunek
$$\sqrt{x^2 + y^2} < 1$$
Wyznacz przybliżenie liczby $\pi$ równe $4 \cdot \frac{K}{N}$

Spróbuj wykonać to zadanie bez wykorzystywania instrukcji pętli (wektoryzacja kodu) oraz zadbaj o wcześniejsze zaalokowanie wykorzystywanych na potrzeby obliczeń macierzy lub wektorów.

Funkcja mc_pi() przyjmuje jeden argument N, który jest liczbą całkowitą określającą ilość losowań.
Gdy argument nie jest podany to używana jest domyślna wartość N=1000

Macierze c.d.

Definiowanie macierzy

Funkcje tworzące macierze

Indeksowanie macierzy

Rozmiary wektorów i macierzy

Iteracja po elementach

Zmiana rozmiaru macierzy

Puste macierze

Usuwanie elementów macierzy

Operacje na macierzach

Macierze i funkcje

Wektory i macierze logiczne

Przydatne funkcje operujące na wartościach logicznych

Operatory logiczne

Zadanie 5 Monte Carlo