
%%======================= toy problem, 1d

c = 50;

Xtrue = [zeros(c/2,1); ones(c/2,1)]; 
G = speye(c) + sparse(2:c,1:c-1,1,c,c) +  sparse(1:c-1,2:c,1,c,c) ...
     +  sparse(3:c,1:c-2,1,c,c) +  sparse(1:c-2,3:c,1,c,c); 
G = G/5; 
r = c;

d = G * Xtrue;
plot(d)
  hold on
dn = d + 0.1*randn(c,1);
plot(dn,'r')

 X = (G + 0.01*eye(c))\dn;
 plot(X,'g')
 hold off
 
 
L = sparse(1:c-1,1:c-1,1,c-1,c) +  sparse(1:c-1,2:c,-1,c-1,c);

plot(L*Xtrue)
szL = size(L);
nrowsL = szL(1);
   


plot(dn,'r*');
  hold on
alpha = 2;           % alpha = 0.001 is still kinda noisy, and alpha=0.01 is nice but very smooth 

Gtemp=[G; alpha*L];
dtemp=[dn; zeros(nrowsL,1)];
n = length(dtemp);
m = size(G);

Niter = 50;
Lone = zeros(Niter,1);
sqR = speye(n);

for iter = 1:Niter,
   Gtil = sqR*Gtemp;
   dtil = sqR*dtemp;    
   X = Gtil\dtil;   
   
     absres = abs(dtemp - Gtemp*X);      % *not*  dtil, Gtil !!!
     Lone(iter) = sum(absres);
     restaper = max(sqrt(absres), 0.01);
     diagR = 1./restaper;
     %diagR(1:m) = 1;  
     %diagR((m+1):end) =   diagR((m+1):end)/mean(diagR((m+1):end));          % rescale the weights so they don't kill the G-fit
     sqR = sparse(1:n, 1:n, diagR);

 plot(X)
end
plot(X,'ok')
hold off
