%% Blatt 10 - Programmieraufgaben:
%% Lineares Ausgleichsproblem (4 Punkte)
% Vorname Nachname Mat.-Nr.
% Zu den äquidistanten Stützstellen 0 = x_0 < \ldots < x_N = 1 seien Datenpunkte y_0, \ldots, y_N \in \mathbb{R} gegeben. Um eine möglichst glatte Kurve durch diese Punkte zu legen, eignen sich kubische Splines. Es ist bekannt, dass durch den Ansatz
% y_k = s(x_k) = \sum_{j = -1}^{N+1} c_j B_3 \left (N(x_k - x_j)\right),
% mit passenden Randbedingungen, eine eindeutige Kurve definiert wird, welche die Punkte interpoliert. Allerdings interpoliert diese Kurve auch mögliche Messfehler, was wir vermeiden möchten. Deshalb reduzieren wir die Anzahl der Freiheitsgrade auf M + 3 \ll N, definieren h = 1/M, t_j = hj für j = -1, 0, 1, \ldots, M+1 und setzen s als
% s(x) = \sum_{j = -1}^{M+1} c_j B_3 \left( \frac{x - t_j}{h} \right)
% an. Die Bedingungen y_k = s(x_k), k = 0, \ldots, N, führen auf das lineare Ausgleichsproblem
% \pmatrix{B_3((x_0 - t_{-1})/h) & \ldots & B_x((x_0 - t_{M+1})/h) \cr B_3((x_1 - t_{-1})/h) & \ldots & B_3((x_1 - t_{M+1})/h) \cr \vdots & & \vdots \cr B_3((x_{N-1} - t_{-1})/h) & \ldots & B_3((x_{N-1} - t_{M+1})/h) \cr B_3((x_N - t_{-1})/h) & \ldots & B_3((x_N - t_{M+1})/h)} \pmatrix{c_{-1} \cr \vdots \cr c_{M+1}} = \pmatrix{y_0 \cr y_1 \cr \vdots \cr y_{N-1} \cr y_N}.
% Da dieses Gleichungssystem überbestimmt ist, kann es in der Regel nicht mehr exakt gelöst werden. Stattdessen betrachten wir das Least-Squares-Problem
% \sum_{k = 0}^N \left| y_k - s(x_k) \right|^2 \to \min.
% Aus der Vorlesung ist bekannt, dass der Koeffizientenvektor {\bf c} = (c_{-1}, \ldots, c_{M+1})^{\rm T} der Lösung des Gleichungssystems {\bf A}^{\rm T} {\bf A c} = {\bf A}^{\rm T} {\bf y} entspricht. Da die Matrix {\bf A}^{\rm T} {\bf A} schlecht konditioniert ist, möchten wir sie in der Praxis nicht explizit berechnen. Für die Lösung von Ausgleichsproblemen existieren Algorithmen, welche das umgehen, was praktischerweise bereits in den \-Solver von Matlab integriert wurde.
%
% Schreiben Sie hier am Ende dieses Live Scripts MATLAB-Funktionen
% cs = getCoefficients(xs, ys, M),
% ys = B_3(xs),
% ys = splineEval(xs, cs),
% wobei getCoefficients die Spline-Koeffizienten berechnen, B_3 den Spline B_3 an den Punkten xs = [x_1, \ldots, x_N] auswerten und splineEval die Funktion s auswerten soll. splineEval soll dabei auf B_3 zurückgreifen.
% Testen Sie Ihre Implementierung mit folgendem Code:
%

xs = linspace(0, 1, 200);
tol = 0.2;

ys = sin(8*pi*xs) + 2*tol*(rand(size(xs)) - 0.5);

for M = [6:13, 197]
    cs = getCoefficients(xs, ys, M);
    subplot(3, 3, min(M-5, 9))
    plot(xs, ys, 'x', xs, splineEval(xs, cs))
    title("M = "+num2str(M))
end

%
% Funktionen:

% TODO: implementieren