%% Smoothing Spline
% NB: Dieses Programm verwendet die Matlab-Funktionen titanium, spmak und fnval, welche nur in der "Curve Fitting Toolbox" von Matlab verfügbar sind.
% Wir betrachten einen Datensatz (x_j,f_j), j=1,\ldots,m, mit 0=x_1 < x_2 < \ldots <x_m=1. Als Bespiel  wird ein in Matlab verfügbarer Test-Datensatz titanium  mit m=49 gewählt (umskaliert, so dass x_1=0, x_m=1 gilt). Dieser Datensatz wird gezeigt.

[x,Daten]=titanium;
m=length(x)
xmax=x(m); xmin=x(1);
x=-xmin/(xmax-xmin)+1/(xmax-xmin).*x;
x(1)=0; x(m)=1; %Korrektur für Rundungsfehler
plot(x,Daten,'+');
hold on;

% Wir betrachten die Smoothing-Spline-Methode mit kubischen Splines (siehe Buch). Hierbei wird eine Minimierungsaufgabe vom Typ
%        \min_{S \in \mathbb{P}_{4,\tau} }\sum_{j=1}^m  {\big(}S(x_j)-f_j\big)^2 + \theta^2 \|S''\|_2^2,                          (1)
% gelöst, wobei der Glättungsterm \|S''\|_2^2 durch eine einfacher berechenbare Approximation ersetzt wird (siehe Buch). Der Regularisierungsparameter \theta soll "geeignet" gewählt werden. Der Einfachheit halber werden für den Splineraum \mathbb{P}_{4,\tau} äquidistante Stützstellen
%      \tau = \{\tau_0,\ldots, \tau_{\ell+1}\}, \tau_j=jh,  j=0,\ldots,\ell+1,~~h:=1/(\ell+1)
% gewählt. Wir wählen \ell, so dass \ell+4 \leq m gilt.  Beachte: {\rm dim}(\mathbb{P}_{4,\tau})=\ell +4.

theta=10^-4
ell=15
Knoten=linspace(0,1,ell+2);

% Die Aufgabe (1) führt auf ein lineares Ausgleichsproblem. Dieses Ausgleichsproblem wird aufgestellt und gelöst. Die berechnete Splinefunktion wird gezeigt.

c=Ausgleichsproblem(x,Daten,Knoten,theta);
makeplotS(Knoten,c');
legend('Daten','Smoothing Spline','Knoten','Location','northwest');
hold off;

%% Aktivität:
% Sie können den Regularisierungsparameterwert \theta variieren oder ein anderen \ell-Wert setzen. Beobachten Sie die Effekte. Man könnte auch einen anderen Datensatz  selber definiieren.
%
%

function c=Ausgleichsproblem(x,Daten,Knoten,theta)
    m=length(x); 
    ell=length(Knoten)-2; h=1/(ell+1);
    A=zeros(m,4+ell);   % Matrix A_T aus Buch
    for j=1:(4+ell)
        Koeffizienten=zeros(1,4+ell); 
        Koeffizienten(j)=1; % j-ter Einheitsvektor
        A(:,j)= SplineAuswertung(x,Knoten,Koeffizienten);
    end
    C=zeros(2+ell,4+ell);   % Matrix C_T aus Buch
    for i=1:(2+ell)
        C(i,i)=1/3*sqrt(6)*h^(-1.5);
        C(i,i+2)=C(i,i);
        C(i,i+1)=-2*C(i,i);
    end
    Atot=[A; theta*C]; 
    rhs=zeros(m+ell+2,1);
    rhs(1:m)=Daten';
    c=Atot\rhs;
end

%

function y=SplineAuswertung(x,Knoten,Koeffizienten)
    Knotenerweiterung=[0 0 0 Knoten 1 1 1]; 
    sp=spmak(Knotenerweiterung,Koeffizienten);
    y=fnval(sp,x);
end

%

function makeplotS(Knoten,Koeffizienten)
    x=linspace(0,1,101);
    y=SplineAuswertung(x,Knoten,Koeffizienten);
    plot(x,y);
    xlim([0,1]);  ylim([0.4 2.2]);
    hold on;
    V=axis;
    plot(Knoten,V(3),'*','MarkerSize',12);
    hold off
end