x=[-0.2163   -0.8328    0.0627    0.1438   -0.5732    0.5955    0.5946   -0.0188]+1;%u=0,sigma=0.5,gauss%
y=[0.5334    0.0296   -0.0478   -0.4162    0.1472   -0.6681    0.3572    0.8118]+1;%u=0,sigma=0.5,gauss%
a=zeros(8,8);                 % a(s,k)
b=zeros(8,8);   b(8,8)=1;   % b(s,k)
r=zeros(8,8,8);               % r(s',s,k)
%%%%%%%%%%%%%%%%%step 1 %%%%%%%%%%%%%%%%%%%
k=1;                 % s'=0,s=0,4  
s0=0.5*gau(x(k),y(k),0);
s1=0.5*gau(x(k),y(k),1);
r(8,8,k)=s0;
r(8,4,k)=s1;
a(8,k)=s0/(s0+s1);
a(4,k)=s1/(s0+s1);
%%%%%%%%%%%%%%%%%step 2 %%%%%%%%%%%%%%%%%%%
k=2;               %  s'=0,4  s=0,4,2,6
s0=0.5*gau(x(k),y(k),0);
s1=0.5*gau(x(k),y(k),1);
r(8,8,k)=s0;
r(8,4,k)=s1;
r(4,2,k)=s0;
r(4,6,k)=s1;
s=[8,4,2,6,1,5,3,7];
for j=1:8
a(j,k)=a(s,k-1)'*r(s,j,k)/(s0+s1);
end
%%%%%%%%%%%%%%%%%step 3 %%%%%%%%%%%%%%%%%%%
k=3;               %  s'=0,4,2,6  s=0,4,2,6,1,5,3,7
s0=0.5*gau(x(k),y(k),0);
s1=0.5*gau(x(k),y(k),1);
r(8,8,k)=s0;
r(8,4,k)=s1;
r(4,2,k)=s0;
r(4,6,k)=s1;
r(2,1,k)=s0;
r(2,5,k)=s1;
r(6,3,k)=s0;
r(6,7,k)=s1;
s=[8,4,2,6,1,5,3,7];
for j=1:8
a(j,k)=a(s,k-1)'*r(s,j,k)/(s0+s1);
end
%%%%%%%%%%%%%%%%%step 4 %%%%%%%%%%%%%%%%%%%
for k=4:8               % s'=0,4,2,6,1,5,3,7  s=0,4,2,6,1,5,3,7,
s0=0.5*gau(x(k),y(k),0);
s1=0.5*gau(x(k),y(k),1);
r(8,8,k)=s0;
r(8,4,k)=s1;
r(4,2,k)=s0;
r(4,6,k)=s1;
r(2,1,k)=s0;
r(2,5,k)=s1;
r(6,3,k)=s0;
r(6,7,k)=s1;
r(1,8,k)=s0;
r(1,4,k)=s1;
r(5,2,k)=s0;
r(5,6,k)=s1;
r(3,1,k)=s0;
r(3,5,k)=s1;
r(7,3,k)=s0;
r(7,7,k)=s1;
s=[8,4,2,6,1,5,3,7];
de=0;
for j=1:8
de=de+a(s,k-1)'*r(s,j,k);
end
for j=1:8
a(j,k)=a(s,k-1)'*r(s,j,k)/de;
end
end
%%%%%%%%%%%%%%%compute b(s',k-1)%%%%%%%%%%%%%%%%%%%%%%%%%%
for k=8:-1:2
    de=0;
    for j=1:8
        de=de+a(s,k-1)'*r(s,j,k);
    end
    for j=1:8
        b(j,k-1)=b(s)'*r(j,s,k)/de;
    end
end
ss1=[4 6 5 7 4 6 5 7];
ss0=[8 2 1 3 8 2 1 3];
l(k)=zeros(1,8);
for k=2:8
    nu=0;
    for j=1:8
        nu=nu+a(s(j),k-1)*r(s(j),ss1(j),k)*b(ss1(j),k);