clear

clc


Data = [
3	60	350
4	85	400
5	125	478
6	140	520
%7	160	570
];

delta = 200
nn = size(Data);
n=nn(1,1);
m=nn(1,2);

W = zeros(n,n);
www=zeros(n,n);
% create the norm
for i = 1:n
   for j = i:n
if i==j
www(i,j)=0;
else
for k=1:1:m
temp= (Data(i,k)-Data(j,k))^2;
summ = summ+temp;
temp=0;
end
www(i,j)= sqrt(summ);
end
temp=0;
summ=0;
end
end
www

for i = 1:n
   for j = i:n
  W(i,j) = exp(-1*((www(i,j)/delta)^2));
end
end



% Weights in matrix are symmetric
for i = 1:n
    for j = 1:n
        if i>j
            W(i,j) = W(j,i);
        end
    end
end
W

% Compute the D matrix 
D = zeros(n,n);
for i = 1:n
    for j = 1:n
        D(i,i) = D(i,i) + W(i,j);
    end
end
D
% Compute P matrix 
P = inv(D)*W;
P
phi0 = eye(n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%CHANGING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%HERE
IminusL = (D^0.5)*P*(D^-0.5);
T = IminusL
gg=cjt(T)
% Sparse factorization provides the translations
% Columns of Q provide the orthogonal basis functions on the space V1
[Q,R] = qr(T);
Q
R
phi1 =Q
T2 = R*cjt(R)
T2=(R*Q)*(R*Q)
T2=  R*conj(R')

[Q1,R1] = qr(T2);
phi2 = Q1;
T4 = R1*cjt(R1);
[Q2,R2] = qr(T4);
phi3 = Q2;
T8 = R2*cjt(R2);

phi1
phi2
phi3
% Determining the wavelet function
[Qw,Rw] = qr(phi0 -( phi1*(phi1)));
wav1 = Qw
[Qw,Rw] = qr(phi0 - (phi2*(phi2)));
wav2 = Qw
[Qw,Rw] = qr(phi0 - (phi3*(phi3)));
wav3 = Qw


c=(Data'*phi1)
c2=(Data'*phi1)/2
d=(Data'*wav1)/2
f=phi1*c'
f2=phi1*c2' + wav1*d'

c3 = (Data'*phi2)
cc=(Data'*phi2)/2
dd=(Data'*wav2)/2
ff=phi2*c3'
f3=phi2*cc'+wav2*dd'
































































































































c4=(Data'*phi3)
ccc=(Data'*phi3)/2
ddd=(Data'*wav3)/2
fff=phi3*c4'
f4=phi3*ccc'+wav3*ddd'
% it is enough to use scaling and wavelet is not needed if the function is
% applied on static data. However, in a dynamic environment use of both scaling wavelet
% function is useful.


