# Rozwiązywanie układów równań

Będzie nam potrzebny cały kod z implementacją macierzy z poprzednich zajęć.
Jeżeli nie było Cię na poprzednich zajęciach, nie udało Ci się dokończyć wszystkich zadań
    lub masz wątpliwości, czy twoje rozwiązania są poprawne, możesz pobrać
    [stąd naszą implementację klasy `Matrix`](https://gist.github.com/Quazyrog/112c31fa763589ec732ff1be3b15cc95).

Następnie obok pliku `matrices.py` z implementacją macierzy z poprzednich labów utwórz plik
    `gauss.py` z poniższym kodem:

```py
from matrices import Matrix, read_matrix_file, write_matrix_file
import sys

def reduce_to_echelon(m):
    ... # TODO: funkcja schodkująca macierze

def compute_solution(echelon_m):
    ... # TODO: funkcja wyliczająca rozwiązanie z wyschodkowanej macierzy

if __name__ == "__main__" and len(sys.argv) >= 2:
    m = read_matrix_file(sys.argv[1])
    reduce_to_echelon(m)
    sol = compute_solution(m)
    print(f"Rozwiązanie: {sol}")
    if len(sys.argv) >= 3:
        write_matrix_file(m, sys.argv[2])
```

Celem na dziś będzie napisanie programu rozwiązującego układ równań:

$$
\begin{cases}
    m_{1\,1} x_1 + m_{1\,2} x_2 + \dots + m_{1\,n} x_n = m_{1\,n+1} \\
    m_{2\,1} x_1 + m_{2\,2} x_2 + \dots + m_{2\,n} x_n = m_{2\,n+1} \\
    \dots \\
    m_{n\,1} x_1 + m_{n\,2} x_2 + \dots + m_{n\,n} x_n = m_{n\,n+1} \\
\end{cases}
$$

Czyli układu opisanego macierzą:

$$
M = \begin{cases}
    m_{1\,1} & m_{1\,2} & \cdots & m_{1\,n} & m_{1\,n+1} \\
    m_{2\,1} & m_{2\,2} & \cdots & m_{2\,n} & m_{2\,n+1} \\
    \vdots   & \vdots   & \ddots & \vdots   & \vdots     \\
    m_{n\,1} & m_{n\,2} & \cdots & m_{n\,n} & m_{n\,n+1} \\
\end{cases}
$$

Napiszemy do tego dwie funkcje:
- `retuce_to_echelon`: która doprowadzi macierz do postaci schodkowej (nie zredukowanej).
- `compute_solution`: która dla macierzy w postaci schodkowej obliczy z niej rozwiązanie.

Wtedy przy pomocy powyższego programu będziemy mogli rozwiązywać układy równań:

```bash
$ python gauss.py plik_z_macierza.txt
```

:::{important}
Jeżeli w którymkolwiek momencie coś nie będzie Ci działać, możesz użyć funkcji
    `write_matrix_file(m)`, żeby wypisać macierz `m` i zobaczyć jej aktualną zawartość.
Najlepiej używaj z dodatkowym `print` przed wywołaniem, żeby odróżniać kolejne wypisywane
    wersje macierzy, jeżeli wypisujesz ją w pętli.
:::

## Schodkowanie macierzy

Zajmiemy się teraz implementacją metody Gaussa.

### Zadanie 1: redukcja kolumny

:::{topic} Zadanie 1
Napisz funkcję `reduce_column_down(m, k)`, która zmodyfikuje podaną macierz `m` przy pomocy
    operacji `Matrix.add_row` w taki sposób, żeby w kolumnie `k` począwszy od wiersza
    `k+1` były same zera.

```{tip}
:class: dropdown
Dla każdego wiersza o numerze $i\colon n \ge i > k$ dodawaj $k$-ty wiersz przemnożony przez
    $-\frac{m_{i\,k}}{m_{k\,k}}$, gdzie $n$ to liczba wierszy w $M$.
```

:::

Przetestuj swoją funkcję na jakiejś przykładowej macierzy.

### Zadanie 2: bardzo naiwny Gauss

:::{topic} Zadanie 2
Korzystając z funkcji `reduce_column_down` uzupełnij funkcję `reduce_to_echelon(m)` w następujący
    sposób tak, by schodkowała macierz.
:::

:::{topic} Zadanie 2.5
Co może pójść nie tak? Wymyśl macierz, dla której to podejście nie zadziała i zobacz, jaki wynik
    otrzymasz.

```{tip}
:class: dropdown
Kiedy ten program wykona dzielenie przez zero?
```

:::

### Zadanie 3: pół-pełny Gauss

Żeby uniknąć dzielenia przez zero (o ile to w ogóle możliwe — może się okazać, że rozwiązanie
    nie jest jednoznaczne i w którymś momencie dostaniemy same zera).
Przykładowo mając w drugim kroku (kolumna o indeksie `1`) macierz:

$$\begin{pmatrix}
    1 & 20 & 4 & 8 & 7 & 10 \\
    0 & 1 & 5 & 2 & 7 & -3 \\
    0 & -6 & 3 & 2 & 1 & 0 \\
    0 & 5 & 2 & 1 & 1 & -7
    \end{pmatrix}$$

Zanim wywoła `reduce_column_down(m, 1)` powinna wykonać `m.swap_rows(1, 2)`, ponieważ ze wszystkich
    wierszy jeszcze dostępnych trzeci ma największą wartość bezwzględną elementu w tej kolumnie.
(Pierwszy wiersz już nie może zostać użyty, chociaż ma 20).

:::{topic} Zadanie 3
Zmodyfikuj poprzednie rozwiązanie tak, żeby zamiast brać po prostu kolejny wiersz, najpierw
    zamieniało wiersze miejscami w taki sposób, żeby do wyzerowania użyć wiersza mającego
    w danej kolumnie największą co do wartości bezwzględnej wartość.

```{tip}
:class: dropdown
Dla wywołania `reduce_column_down(m, k)`:
1. Najpierw w pętli znajdź takie $i$, że $n \ge i \ge k$ oraz: 
    $$ \lvert m_{i\,k}\rvert = max\{\lvert m_{k\,k}\rvert, \lvert m_{k+1\,k}\rvert, \dots, \lvert m_{n\,k}\rvert\} $$
2. Następnie zamień wiersze $i$-ty z $k$-tym.
3. Wykonaj resztę operacji tak jak w poprzedniej wersji.
```

:::

## Obliczanie rozwiązania z wyschodkowanej macierzy

### Zadanie 4: obliczanie pojedynczej zmiennej

Zaczniemy od funkcji, która mając obliczone wartości zmiennych $x_{k+1}, \dots, x_n$ oraz
    pojedyncze równanie w *już wyschodkowanej macierzy* obliczy z niego kolejną zmienną.
Ta sytuacja jest przedstawiona na rysunku:

![Ilustracja obliczania pojedynczej zmiennej](static/gauss_solve_single_row.png)

Podświetlony wiersz odpowiada równaniu:

$$ 0x_1 + 0x_2 + \dots + m_{k\,k} x_k + m_{k\,k+1} x_{k+1} + m_{k\,k+2} x_{k+2} + \dots + m_{k\,n} x_n = m_{k\,n+1} $$

Na tym etapie rozwiązywania układu równań wartości zmiennych $x_{k+1}, \dots, x_n$ są już wyliczone.
Zatem kolejną zmienną $x_k$ wyliczamy jako:

$$ x_k = \frac{1}{m_{k\,k}}\left(m_{k\,n+1} - m_{k\,k+1} x_{k+1} - m_{k\,k+2} x_{k+2} - \dots - m_{k\,n} x_n\right) $$

W programie współczynniki będziemy brać z macierzy — $m_{k\,j} = \texttt{m.get(k - 1, j - 1)}$.
Wartości zmiennych będziemy pobierać i zapisywać do pojedynczej listy $x_j = \texttt{sol[j - 1]}$.

Później możemy kilka razy wywołać tę funkcję, aby wyliczyć całe rozwiązanie układu z wyschodkowanej
    macierzy:

```py
m = Matrix(3, 4)
m.set(0, 0, 1); m.set(0, 1, 5); m.set(0, 2, 4); m.set(0, 3, 9)
m.set(1, 1, 1); m.set(1, 2, -3); m.set(1, 3, 2)
m.set(2, 2, 1); m.set(2, 3, -1)
sol = [None] * 3  # Na początek nie znamy żadnej z trzech zmiennych
sol[2] = extend_solution(m, 2, sol)
sol[1] = extend_solution(m, 1, sol)
sol[0] = extend_solution(m, 0, sol)
if sol == [18, -1, -1]:
    print("extend_solution: OK")
else:
    print("extend_solution: ERROR")
```

:::{topic} Zadanie 4
Napisz funkcję `extend_solution(m, k, x)`, która na podstawie wyschodkowanej
    macierzy, numeru wiersza oraz częściowego rozwiązania obliczy kolejną zmienną w rozwiązaniu
    układu zadanego tą macierzą.
:::

### Zadanie 5: obliczanie wszystkich zmiennych

![Ilustracja obliczania pojedynczej zmiennej](static/gauss_solve_all.png)

:::{topic} Zadanie 5
Teraz korzystając z funkcji `extend_solution` uzupełnij funkcję `compute_solution`,
    która obliczy całe rozwiązanie (podobnie jak w przykładowym kodzie do pierwszego zadania).

```{tip}
:class: dropdown
1. Niech $n = \texttt{m.num_rows()}$.
2. Zacznij od zainicjowania listy na wynik jako `sol = [None] * n`.
3. Zaczynając od $n$-tego wiersza, kończąc na pierwszym przypisuj kolejne wartości rozwiązania
        korzystając z funkcji `extend_solution`. 
    (Pamiętaj, że Python indeksuje listy od zera).
```

:::

Poniżej masz kilka macierzy schodkowych, których możesz użyć do przetestowania:

1. Skopiuj i zapisz wartości macierzy do pliku tekstowego.
    *Nie kopiuj pierwszej linijki z rozwiązaniem*.
2. Wykonaj polecenie w terminalu: `python gauss.py <nazwa_pliku>`.
3. Porównaj otrzymane rozwiązanie.

```
# Rozwiązanie: [7]
1 7
```

```
# Rozwiązanie: [4, 2]
1  3  10
0  1   2
```

```
# Rozwiązanie: [0, 3, 2]
1  2  1   8
0  1  4  11
0  0  1   2
```

```
# Rozwiązanie: [2, -1, 3, 1, 3]
1  1  2  0  1   10
0  1  3  1  2   15
0  0  1  2  1    8
0  0  0  1  3   10
0  0  0  0  1    3
```

```
# Rozwiązanie: [-11868, 4832, -1446, 249, -55, 18, -1, -1]
1  2  -1   3   0   1  4  -2    5
0  1   3  -2   1   4  0   1   12
0  0   1   5  -3   2  1   3   -4
0  0   0   1   4  -1  2   2    7
0  0   0   0   1   3  0  -1    0
0  0   0   0   0   1  5   4    9
0  0   0   0   0   0  1  -3    2
0  0   0   0   0   0  0   1   -1
```

:::{topic} Zadanie 2.5
Przetestuj swoją funkcję `extend_solution` na jednej lub kilku z podanych wyżej macierzy.
:::

## Zadania domowe

:::{topic} Zadanie `g.err` **[1 punkt]**
Podaj przykład macierzy układu, który: a) nie ma rozwiązań; b) ma nieskończenie wiele rozwiązań.
Uruchom na tych macierzach zaimplementowaną na ćwiczeniach funkcję do schodkowania.
Jaki jest efekt?
:::

:::{topic} Zadanie `g.solve` **[1 punkt]**
Napisz funkcję `solve(m)`, która zwróci rozwiązanie układu równań opisanego macierzą `m`.
Funkcja może działać tylko w przypadkach, kiedy rozwiązanie jest dokładnie jedno.
:::

:::{topic} Zadanie `g.better_gauss` **[2 punkty]**
Popraw implementację funkcji `reduce_to_echelon` z ćwiczeń tak, żeby wykrywała macierze
    opisujące układy, które nie mają dokładnie jednego rozwiązania.
Jeżeli układ równań ma więcej rozwiązań, lub nie ma ich wcale, funkcja powinna zgłosić błąd
    instrukcją `raise ValueError("nieobsługiwany układ")`.
Weź pod uwagę błędy zaokrągleń: wartości mniejsze od $10^{-6}$ traktuj jako zera.

**Uwaga:** możesz dostać punkty tylko za jedno z zadań `g.better_gauss` i `g.better_gauss2`.
:::

:::{topic} Zadanie `g.coord` **[2 punkty]**
Napisz funkcję `compute_coordinates(vector, base)`, która zwróci współrzędne podanego wektora
    w podanej bazie.

- Argument `vector` to lista współrzędnych wektora w bazie standardowej.
- Argument `base` to lista wektorów bazowych, z których każdy jest opisany listą współrzędnych
    w bazie standardowej.
- Wynik to lista współrzędnych, gdzie $i$-ty element odpowiada $i$-temu wektorowi z listy `base`.

Załóż, że podany argument jest poprawną bazą przestrzeni liniowej.

Na przykład `compute_coordinates([3, 6], [[1, 0], [0, 2]])` powinno zwrócić listę `[3, 3]`.
:::

:::{topic} Zadanie `g.better_gauss2` **[3 punkty]**
Rozwiąż zadanie `g.better_gauss`, ale zamiast traktować jako zera wartości mniejsze od $10^{-6}$
    wyznacz tę granicę na podstawie oryginalnej macierzy:

- Przed rozpoczęciem schodkowania znajdź największy co do wartości bezwzględnej element macierzy.
- Traktuj jako zera te wartości, które są mniejsze od $10^{-6}m$, gdzie $m$ to wartość bezwzględna
    znalezionego wcześniej największego elementu.

**Uwaga:** możesz dostać punkty tylko za jedno z zadań `g.better_gauss` i `g.better_gauss2`.
:::

:::{topic} Zadanie `g.raycast` **[3 punkty]**
Uzupełnij implementację poniższej funkcji, która oblicza przecięcie płaszczyzny z prostą:

```py
class Point3D:
    # Punkt w przestrzeni trójwymiarowej.
    def __init__(self, x, y, z):
        self.x = x
        self.y = y
        self.z = z


class Plane3D:
    # Płaszczyzna w przestrzeni trójwymiarowej dana za pomocą trzech punktów na niej leżących.
    def __init__(self, p1, p2, p3):
        self.p1 = p1
        self.p2 = p2
        self.p3 = p3

class Line3D:
    # Prosta w przestrzeni trójwymiarowej przechodząca przez podane punkty.

def raycast(plane, line):
    ... # TODO...
    return intersection  # zwraca Point3D, w którym prosta przecina linię.
```

Załóż, że prosta nie jest równoległa do tej płasczyzny.
Argumenty — płaszczyzna i prosta — są podawane to funkcji `raycast` jako obiekty klas zdefiniowanych
    w powyższym kodzie.

**Uwaga:** częścią rozwiązania jest również matematyczne uzasadnienie działania tej funkcji.
Możesz je dołączyć jako zdjęcie.
:::