% 07/19/2007: centering the data
%
%
isShift = 0;
useSparse = 0;
isImpaint = 0;
isAdaptive = 1;

% fromlevel is the smallest detail I will use
a_fromlevel=1;
% tolevel is the biggest approximation I will use
% tolevel is corresponding to the 'level' in wavedec2() fn
% previously I always decompose to level=5 which may not be
% good;
a_tolevel=2;
b_fromlevel=1;
b_tolevel=2;
wname='db1';
%noise level 0.1 is very high for level=7
relative_sigma=1;
%fixed random noises
state = 139137;
rand('state', state);
randn('state', state);

dimx=32;
dimy=32;

interval1=1/dimx;
domain1=interval1:interval1:1;
interval2=1/dimy;
domain2=interval2:interval2:1;
[t1,t2]=meshgrid(domain2,domain1);

Y=zeros(dimx,dimy);
A=2*ones(dimx,dimy);
[C1,S1] = wavedec2(A,b_tolevel,wname);
B=(5*t1-6*t2-5);
[C2,S2] = wavedec2(B,b_tolevel,wname);

%A: 3-d step function
% for i = 1:dimx/4
%   for j = 1:dimy/2+2
%     A(i,j)=1.0;
%   end
% end
% 
% for i = (dimx/4+1):dimx
%   for j = 1:dimy/2+2
%     A(i,j)=2.0;
%   end
% end
% 
% for i = 1:dimx/2 
%   for j = (dimy/2+2+1):dimy
%     A(i,j)=3.2;
%   end
% end
% 
% for i = dimx/2+1:dimx
%   for j = dimy/2+2+1:dimy
%     A(i,j)=4.0;
%   end
% end
% surf(t1,t2,A);
% [C1,S1] = wavedec2(A,a_tolevel,wname);
% A_nonezeros=sum(C1~=0);
% 
% %B: 3-d step function
% for i = 1:dimx/2
%   for j = 1:dimy
%     B(i,j)=2.0;
%   end
% end
% 
% for i = (dimx/2+1):dimx
%   for j = 1:dimy
%     B(i,j)=1.0;
%   end
% end
% 
% for i = 1:dimx/4 
%   for j = (dimy/2+1):dimy
%     B(i,j)=3.2;
%   end
% end
% 
% for i = dimx/4+1:dimx
%   for j = dimy/2+1:dimy
%     B(i,j)=1.5;
%   end
% end
% surf(t1,t2,B);
% [C2,S2] = wavedec2(B,b_tolevel,wname);
% B_nonezeros=sum(C2~=0);

%X=5*(4.*t1-2).*exp((-1)*((4.*t1-2).^2+(4.*t2-2).^2));

%X=5*exp((-1)*((4.*t1-2).^2+(4.*t2-2).^2))+20; %nonlinear and elevated;add big constant will
                                              % affect A's estimation
%X=5*exp((-1)*(4.*t1-2).^2)+5;
%X=5*t1.^2+6*t2.^2+3;
%X=5*cos(t1).*t2+2;
%X=3*cos(5*pi.*t1)+4*sin(2*pi.*t2).*t2+4;
%X=3*cos(2*pi.*t1)+4*sin(5*pi.*t2).*t2+4;
X=(3*cos(8*pi.*t1)+4*sin(8*pi.*t2).*t2+4);
X0=X;
%this is challenging since it is too flat, we have to cut
%from_level=2 and to_level=level-1 to make it work
%to do: maybe adaptive lasso can fix this;
%       or simply more lambda for level coefficients
%highresX=3*cos(2*pi.*t1)+4*sin(2*pi.*t2)+4;
%this is much better with high-freq and more variations everywhere
%highresX=3*cos(8*pi.*t1)+4*sin(8*pi.*t2)+4;

%lowresX=3*cos(2*pi.*k1)+4*sin(2*pi.*k2)+4;
%X=3*cos(5*pi.*t1+5*pi.*t2)+4;
%X=3*cos(5*pi.*t1)+4*sin(2*pi.*t2)+4*t1+6*t2+4;
%X=3*cos(10*pi.*t1)+4*sin(2*pi.*t2)+4*t1+6*t2;
%X=5*t1+6*t2+3; % linear will spoil the things.
%X=-5*t1+6*t2-5;
%X=5*t1-6*t2-5;
%X=-5*t1-6*t1+5;
%X=3*ones(dimx,dimy); %constant
% resizedX=imresize(lowresX,[dimx dimy],'bicubic',9);
% surf(k1,k2,lowresX)
% surf(t1,t2,resizedX)
surf(t1,t2,X)
%surf(t1,t2,highresX-resizedX);
Y=A+X.*B;
surf(t1,t2,Y);
%lowX0=resizedX;
% highX0=highresX;
% Y0=Y;
% Y=Y+randn([dimx dimy])*sqrt(0.01);
% highresX=highresX+randn([dimx dimy])*sqrt(0.0);
% [C,S] = wavedec2(highresX,2,'bior3.7');
% [thr,sorh,keepapp] = ddencmp('den','wv',highresX);
% [Xden,CXC,LXC,PERF0,PERFL2] = ...
% wdencmp('gbl',C,S,'bior3.7',2,thr,sorh,keepapp);
% surf(t1,t2,Xden);
% surf(t1,t2,Xden-highX0);

%here is a way: shrink a image, then enlarge it back
%so we will get rid of noise and still close to the
%true one.
% Xsmall=imresize(highresX,0.5,'bilinear');
% Xresize=imresize(Xsmall,2,'bilinear');
% surf(t1,t2,Xresize);
% surf(t1,t2,Xresize-highX0);
% surf(t1,t2,Xresize-Xden);
% k=2;
% mask=zeros(dimx,dimy);
% for i = dimx/2+k:dimx-3
%   for j = dimy/2+2+k:dimy-3
%     Y(i,j)=0;
%     resizedX(i,j)=0; %set both to zero can help
%     highresX(i,j)=0;
%     X1(i,j)=1; 
%     mask(i,j)=1;
%   end
%   k=k+1;
% end


% 2 missing parts, but this is very bad situation:
% k=2;
% for i = 3:dimx/2-k
%   for j = 3:dimy/2-k
%     Y(i,j)=0.0;
%     X(i,j)=0.0;
%     mask(i,j)=1;
%   end
%   k=k+1;
% end

% k=2;
% for i = 5+k:dimx/2-4;
%   for j = 4+k:dimy/2-4;
%     X(i,j)=1;
%     X2(i,j)=1;
%     mask(i,j)=1;
%   end
%   k=k+1;
% end

if isImpaint == 1
[paintedX, newmask]=simpleImpaint(X,mask,1);
surf(t1,t2,paintedX);

[paintedY, newmask]=simpleImpaint(Y,mask,1);
surf(t1,t2,paintedY);
X=paintedX;
Y=paintedY;
end

%This will stablise the calculation.
if(isShift==1)
  Xshift=mean(mean(X));
  X=X-Xshift;
else
    Xshift=0;
end

ybar=mean(mean(Y));
SD=sqrt((1/(dimx*dimy-1))*sum(sum((Y-ybar).^2)));
if(isShift==1)
  Yshift=mean(mean(Y));
  %This will stablise calculation
  Y=Y-Yshift;
else
  Yshift=0;
end

%Now data is ready, we have to make 
% the "A0" matrix 
% depends on the level for decomposition
% small detail -> large detail 
% d,v,h,...,d,v,h,a
% a is only for highest level
%let's delete A and B to free memory
% A=[];
% B=[];
%TODO: 1. basic_info can be put into function for all images
[A01 a_totalbase]=getHaarA([],dimx,dimy,a_fromlevel,a_tolevel,useSparse);
[A02 b_totalbase]=getHaarA(X,dimx,dimy,b_fromlevel,b_tolevel,useSparse);
A0=[A01 A02]; % A03 A04];

%TODO: padding and extention can be put into function again
Y0=Y;
% Y=padding(Y,basic_info,b_tolevel);
% Y=fliplr(Y);
% Y=flipud(Y);

y  = reshape(Y,size(Y,1)*size(Y,2),1);   % measurements with no noise
max_lambda = find_lambdamax_l1_ls(A0',y);
rel_tol = 1.0e-2;     % relative target duality gap
lambda = 1.0e-2;      % regularization parameter
[x1,status]=l1_ls(A0,y,lambda,rel_tol);

% if isAdaptive == 0
%     lambda=ones(size(A0,2),1)*1.0e-1;
%     lambda(1:length(lambda)*3/8) ...
%         =lambda(1:length(lambda)*3/8)*10;
%     lambda(length(lambda)/2+1:length(lambda)*7/8 )...
%         =lambda(length(lambda)/2+1:length(lambda)*7/8)*10;
% else
%     %next is adaptive lasso from Zou06
%     scale=1.0;
%     if size(A0,1) >= size(A0,2)
%         [lambda1,bint,r,rint,stats] = regress(y,A0);
%     else
%         lambda1 = myridge(y, A0, 1e-2);
%     end
%     lambda=ones(size(A0,2),1)./abs(lambda1)*scale;
% end
%[x1,status]=l1_ls_1(A0,y,lambda,rel_tol);

% lambda = 1.0e-1;      % regularization parameter
% [x1,status]=l1_ls(A0,y,lambda,rel_tol);

%[x1,bint,r,rint,stats] = regress(y,A0);
%let's do inverse dwt with x1
n=dimx;
c1=zeros(1,n*n);
len=size(x1,1);
c1(1:a_totalbase)=flipud(x1(1:a_totalbase,:))';
a=waverec2(c1,S1,'haar');
% a=recon(x1(1:length(x1)/2),basic_info,a_fromlevel,a_tolevel,dimx,dimy);
surf(t1,t2,a);
surf(t1,t2,A-a);

c2=zeros(1,n*n);
c2(1:b_totalbase)=flipud(x1(a_totalbase+1:a_totalbase+b_totalbase,:))';
b=waverec2(c2,S2,'haar');
% b=recon(x1(length(x1)/2:length(x1)),basic_info,b_fromlevel,b_tolevel,dimx,dimy);
surf(t1,t2,b);
surf(t1,t2,B-b);

yhat = a+b.*X0+Yshift;
surf(t1,t2,yhat);

%true Y
surf(t1,t2,Y0);
resid=(Y0-yhat)./(relative_sigma*SD);
%resid(mask==1)=0;
surf(t1,t2,resid);