% Optimizacijski problem min_y ||Ay - b||_2^2 + alfa ||Ly||_2^2 lahko
% uporabimo tudi za interpolacijo tock, v katerih nam manjkajo nekatere
% vrednosti.

% Denimo, da imamo v tockah xs  = [xLevo xManjka xDesno]' izracunane vrednosti
% funkcije f, nato pa vrednosti v tockah xManjka izgubimo. Za ustrezno izbrane
% A, b in L lahko manjkajoce vrednosti dobimo preko regularizacije Tikhonova.

% Oglej si spodnje kandidate za matrike L. Denimo, da smo ze izracunali dela
% resitve yLevo in yDesno, ki pripadata tockam xLevo in xDesno.
% Razisci:
%    - Kaksne so lastnosti vrednosti yManjka pri spodnjih treh izberah za L?
%      Za idejo si oglej dobljeni graf in nato dokazi svojo domnevo.
%    - Ali ima parameter alfa kak vpliv na gornje lastnosti? Loci alfa > 0
%      in alfa = 0.
%    - Bi se dalo najti matrike L3, L4 ... ?
%
% Predpostavi, da so tocke v xs ekvidistantne.

n = 25;
% Izberemo, kateri podatki nam manjkajo
nL = 10;
nManjka = 10;
nD = n - (nL + nManjka);

% Spremeni spodnjo vrstico in izberi L1, da se prepricas o potrebnosti
% pogoja, da so tocke ekvidistantne.
ekvidistantne = true;
if ~ekvidistantne
    gosto = 15;
    xs = unique([linspace(0, pi/4, gosto) linspace(pi/4, pi, 1 + n - gosto)]');
else
    xs = linspace(0, pi, n)';
end

ys = sin(xs);  % neka funkcija

% Posebne matrike L
L0 = eye(n);       % Nas to morda spomni na priblizke za f?
L1 = diag(-1 * ones(n, 1)) + diag(ones(n - 1, 1), 1);
L2 = L1 + L1';
L1 = L1(1:n-1,:);  % Nas to morda spomni na priblizke za f'?
L2 = L2(2:n-1,:);  % Nas to morda spomni na priblizke za f''?

A = [eye(nL), zeros(nL, nManjka + nD);
     zeros(nD, nL +  nManjka), eye(nD)];
b = A * ys;  % izgubimo vrednosti v xManjka
L = L1;      % izberemo enega od L-jev

alfa = 0.1;  % izberemo neko vrednost za alfa
resitev = (A' * A + alfa * L' * L) \ (A' * b);  % naivno resimo sistem
plot(xs(1:nL), resitev(1:nL), 'bo;yLevo, yDesno;',
     xs(nL + 1: nL + nManjka), resitev(nL + 1: nL + nManjka), 'ro;yManjka;',
     xs(nL + nManjka + 1:end), resitev(nL + nManjka + 1:end), 'bo');
h = legend ('location', 'northeastoutside');
legend (h, 'location', 'northeast');