clear;
%load leleccum; 
%even always
n=100;
%s = leleccum(1:n); 
s=1:n;
l_s = length(s);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%wrap around is successful
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
wname='bior1.3';
[Lo_D,Hi_D,Lo_R,Hi_R] = wfilters(wname);
[C,L] = wavedec(s,2,wname);

% cA1 = appcoef(C,L,wname,1);
% % cD3 = detcoef(C,L,3); 
% cD2 = detcoef(C,L,2); 
% cD1 = detcoef(C,L,1);

k=2^1;
%a < k
a=mod(length(s)-length(Lo_D),k); %gap between the end of siginal and end of last base still in the range
%min b st. ext=b*k-a<len(base) <==> j*k<len(base)+a
if mod(length(Lo_D)+a,k)==0
    b=(length(Lo_D)+a)/k-1;
else
    b=floor((length(Lo_D)+a)/k);
end
ext=b*k-a;
s4=s;
s4(n+1:n+ext)=s(1:ext);
%since it biorthognal flip the signal first
s4=fliplr(s4);

d1=[];
for i=0:(length(s4)-length(Hi_D))/k;
  d1(i+1)=sum(Hi_D.*s4(k*i+1:k*i+length(Hi_D)));
end

a1=[];
for i=0:(length(s4)-length(Lo_D))/k;
  a1(i+1)=sum(Lo_D.*s4(k*i+1:k*i+length(Lo_D)));
end

s2=zeros(1,n);
%no need to flip if not using biorthogonal bases
a1=fliplr(a1);
d1=fliplr(d1);
%1-level reconstruction
base=zeros(1,length(s));
base(1:length(Hi_R))=Hi_R;
for i=0:length(d1)-1  % no need to go to end
 %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
 s2=base*d1(i+1)+s2;
base=circshift(base,[0 k]);
end

base=zeros(1,length(s));
base(1:length(Lo_R))=Lo_R;
for i=0:length(d1)-1  % no need to go to end
 %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
 s2=base*a1(i+1)+s2;
base=circshift(base,[0 k]);
end

%2D example
n=100;
interval=1/n;
domain=interval:interval:1;
[t1,t2]=meshgrid(domain);
%X=-5*t1+6*t2-5;
X=5*t1-6*t2-5;
%X=-5*t1-6*t1+5;
surf(t1,t2,X)

% [k1,k2]=meshgrid(1:length(scalebase));
% vv=kron(scalebase,scalebase');
% surf(k1,k2,vv);
% vw=kron(scalebase,waveletbase');%horizontal
% surf(k1,k2,vw);
% wv=kron(waveletbase,scalebase');%vertical
% surf(k1,k2,wv);
% ww=kron(waveletbase,waveletbase');%diag
% surf(k1,k2,ww);

[k1,k2]=meshgrid(1:length(Hi_R));
vv2=kron(Lo_R,Lo_R');
surf(k1,k2,vv2);
vw2=kron(Lo_R,Hi_R');%horizontal
surf(k1,k2,vw2);
wv2=kron(Hi_R,Lo_R');%vertical
surf(k1,k2,wv2);
ww2=kron(Hi_R,Hi_R');%diag
surf(k1,k2,ww2);

vv1=kron(Lo_D,Lo_D');
surf(k1,k2,vv1);
vw1=kron(Lo_D,Hi_D');%horizontal
surf(k1,k2,vw1);
wv1=kron(Hi_D,Lo_D');%vertical
surf(k1,k2,wv1);
ww1=kron(Hi_D,Hi_D');%diag
surf(k1,k2,ww1);
%1-level test
X0=X;
k=2^1;
%a < k
a=mod(length(s)-length(Lo_D),k); %gap between the end of siginal and end of last base still in the range
%min b st. ext=b*k-a<len(base) <==> j*k<len(base)+a
if mod(length(Lo_D)+a,k)==0
    b=(length(Lo_D)+a)/k-1;
else
    b=floor((length(Lo_D)+a)/k);
end
ext=b*k-a;
X(:,n+1:n+ext)=X(:,1:ext); %extend cols first
X(n+1:n+ext,:)=X(1:ext,:); %extend rows next

X=fliplr(X);
X=flipud(X);

h1=[];
for i=0:(size(X,2)-length(Hi_R))/k %col
    for j=0:(size(X,1)-length(Hi_R))/k %row
        h1(j+1,i+1)=sum(sum(vw1.*X(k*j+1:k*j+length(Hi_R),k*i+1:k*i+length(Hi_R))));
    end
end

v1=[];
for i=0:(size(X,2)-length(Hi_R))/k %col
    for j=0:(size(X,1)-length(Hi_R))/k %row
        v1(j+1,i+1)=sum(sum(wv1.*X(k*j+1:k*j+length(Hi_R),k*i+1:k*i+length(Hi_R))));
    end
end

diag1=[];
for i=0:(size(X,2)-length(Hi_R))/k %col
    for j=0:(size(X,1)-length(Hi_R))/k %row
        diag1(j+1,i+1)=sum(sum(ww1.*X(k*j+1:k*j+length(Hi_R),k*i+1:k*i+length(Hi_R))));
    end
end

approx1=[];
for i=0:(size(X,2)-length(Hi_R))/k %col
    for j=0:(size(X,1)-length(Hi_R))/k %row
        approx1(j+1,i+1)=sum(sum(vv1.*X(k*j+1:k*j+length(Hi_R),k*i+1:k*i+length(Hi_R))));
    end
end

%1-level recovery
h1=fliplr(h1);
h1=flipud(h1);
v1=fliplr(v1);
v1=flipud(v1);
diag1=fliplr(diag1);
diag1=flipud(diag1);
approx1=fliplr(approx1);
approx1=flipud(approx1);

base=zeros(size(X0,1),size(X0,2));
Xrecon=zeros(size(X0,1),size(X0,2));
base(1:length(Hi_R),1:length(Hi_R))=vw2;
for i=0:size(h1,2)-1 %col
    for j=0:size(h1,1)-1 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*h1(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
    base=circshift(base,[0 k]);
end

base=zeros(size(X0,1),size(X0,2));
base(1:length(Hi_R),1:length(Hi_R))=wv2;
for i=0:size(h1,2)-1 %col
    for j=0:size(h1,1)-1 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*v1(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
    base=circshift(base,[0 k]);
end

base=zeros(size(X0,1),size(X0,2));
base(1:length(Hi_R),1:length(Hi_R))=ww2;
for i=0:size(h1,2)-1 %col
    for j=0:size(h1,1)-1 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*diag1(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
    base=circshift(base,[0 k]);
end

base=zeros(size(X0,1),size(X0,2));
base(1:length(Hi_R),1:length(Hi_R))=vv2;
for i=0:size(h1,2)-1 %col
    for j=0:size(h1,1)-1 %row
     %  s5(k*i+1:k*i+length(Hi_R))=Hi_R*d1(i+1)+s5(k*i+1:k*i+length(Hi_R));
    Xrecon=base*approx1(j+1,i+1)+Xrecon;
    base=circshift(base,[k 0]);
    end
    base=circshift(base,[0 k]);
end
surf(Xrecon)
surf(Xrecon-X0)