%TE-MAP均衡
%copyright:flyingxiang
function  [Lec,saidaequa]=Lblocoequa(y,var)
%得到均衡的輸出和外部信息


S=4;  %狀態數                                                   
P=zeros(4,4,length(y));
oute=[-1.629 -.815 0 0;0 0 .001 .815;-.815 -.001 0 0;0 0 .815 1.629];        %對應網格圖的輸出矩陣


A1=[0 1 0 0;0 0 0 1;0 1 0 0;0 0 0 1];
A0=[1 0 0 0;0 0 1 0;1 0 0 0;0 0 1 0];


pos=[1 1 0 0;0 0 1 1;1 1 0 0;0 0 1 1];
gamma=zeros(4,4,length(y));
B1=zeros(4,4,length(y));
B0=zeros(4,4,length(y));
num=zeros(1,length(y));
den=zeros(1,length(y));

%計算gamma
for k=1:length(y)   
   for i=1:4
      for j=1:4
         arg=-(y(k)-oute(i,j))^2;
       
          P(i,j,k)=exp(arg/(2*var))*pos(i,j);
      
           end
       end
   end
gamma=P;

N=length(y)+1;

%初始化前向遞規向量f
f=zeros(S,1,N);
f(:,1,1)=[1;0;0;0] ;
c=zeros(1,N);
c(1)=1;
%遞規計算前向向量f
for i=1:(N-1)
   f(:,1,i+1)=gamma(:,:,i)'*f(:,1,i);
   c(i+1)=sum(f(:,1,i+1));
   f(:,1,i+1)=f(:,1,i+1)./c(i+1);
end

%初始化后向遞規向量b
b=zeros(S,1,N);
b(:,1,N)=[1;1;1;1];
%遞規計算后向向量b
for i=(N-1):-1:1
   b(:,1,i)=gamma(:,:,i)*b(:,1,i+1);
   cb(i)=sum(b(:,1,i));
   b(:,1,i)=b(:,1,i)./cb(i);
end
b;


for n=1:length(y)
B1(:,:,n)=A1.*gamma(:,:,n);
B0(:,:,n)=A0.*gamma(:,:,n);
end

 
 
 L=zeros(1,N-1);
 saida=zeros(1,N-1);

for k=1:(N-1)
    
   num(k)=f(:,:,k)'*B1(:,:,k)*b(:,:,k+1);
    
   den(k)=f(:,:,k)'*B0(:,:,k)*b(:,:,k+1);
      Lec(k)=log((num(k))/(den(k)));            %得到外部信息（注意不是真正意義上的，應該說包含外部信息，用來相減）
end
 
for i=1:length(y)
   if Lec(i)>=0
      saidaequa(i)=1;
   elseif Lec(i)<0
      saidaequa(i)=-1;
   end
end