//共役勾配法(簡易版)
//////////////////////////////////////////////////////////////////////////////////////////
function [S]=aboutS(S_old,KB,A,A_old,WSet)
S=S_old+[[(A(WSet)-A_old(WSet)).*Y(WSet)]'*KB]';
endfunction
//////////////////////////////////////////////////////////////////////////////////////////
function [WSet,p_old,r_old,q2]=WorkingSet_CG(A,Y,S,Beta,p_old,ii)
i=1,up=1,down=N,WSet=[],top=[],bottom=[],Yp=[],q2=0,
r=(-ones(N,1)+Y.*S)
r_old=r;
p=r+Beta*r'*r*p_old//左右にどれだけ動くか
p_old=p;
Yp=(Y.*p)
[sort_Yp,ith]=sort(Yp)
while (i<=q/2)&(up<down),
while ( ( ( A(ith(up))<E0 )&( Y(ith(up))==1 ) )..
|( ( A(ith(up))>=C-E0 )&( Y(ith(up))==-1 ) ) )..
then,up=up+1
end
while ( ( ( A(ith(down))<E0 )&( Y(ith(down))==-1 ) )..
|( ( A(ith(down))>=C-E0)&( Y(ith(down))==1 ) ) )..
then,down=down-1
end
if (up<down) then,top(i)=ith(up),bottom(i)=ith(down),
i=i+1,up=up+1,down=down-1,end
end
WSet=[top;bottom],q2=size(WSet,1)
endfunction
//////////////////////////////////////////////////////////////////
function [AB,AN,YB,YN,KB,QBB,QBN]=separate_1(A,X,Y,WSet,q2)
Yi=[],Yj=[],KB=zeros(1,N)
left_right=X(WSet,:).*.ones(N,1)-ones(q2,1).*.X
normX2=sum( (left_right).^2,2 )
KB=matrix( exp( -(normX2)/(o^2) ),N,q2 )'
Yi=Y(WSet)*ones(1,N)
Yj=ones(q2,1)*Y'
QB=[Yi.*Yj.*KB]
QBB=QB(:,[WSet']), QB(:,[WSet'])=[], QBN=QB
AB=A([WSet]), A([WSet])=[], AN=A
YB=Y([WSet]), Y([WSet])=[], YN=Y
endfunction
//////////////////////////////////////////////////////////////////////////////////////////
function [AB,AN,YB,YN,KB,QBB,QBN]=separate(A,X,Y,WSet,q2,x2,Yj)
x1=X(WSet,:).*.ones(N,1)
if q2<>q then,
x2_=ones(q2,1).*.X
Yj_=ones(q2,1)*Y'
Kr=exp(-sum((x1-x2_).^2,2)/(o^2))
KB=matrix( Kr,N,q2 )'
Yi=Y(WSet)*ones(1,N)
QB=clean([Yi.*Yj_.*KB])
else,
Kr=exp(-sum((x1-x2).^2,2)/(o^2))
KB=matrix( Kr,N,q2 )'
Yi=Y(WSet)*ones(1,N)
QB=clean([Yi.*Yj.*KB])
end
QBB=QB(:,[WSet']), QB(:,[WSet'])=[], QBN=QB
AB=A([WSet]), A([WSet])=[], AN=A
YB=Y([WSet]), Y([WSet])=[], YN=Y
endfunction
//////////////////////////////////////////////////////////////////////////////////////////
function [fin,I_0,A0,AC,sa]=judgement(S,A,Y,Tau)
I_0=[],I_1=[],I_2=[],I_3=[],I_4=[],sa=0,
A0=[],AC=[],A_up=[],A_low=[],
I_0=find( (E0<=A)&(A<=C-E0) )
I_1=find( (Y==1)&(A<E0) )
I_2=find( (Y==-1)&(C-E0<A) )
I_3=find( (Y==1)&(C-E0<A) )
I_4=find( (Y==-1)&(A<E0) )
A0=[I_1,I_4]
AC=[I_2,I_3]
max_0=max(S(I_0)-Y(I_0))
min_0=min(S(I_0)-Y(I_0))
// A_up=find( ( (E0<=A)&(A<=C-E0) )|( (Y==1)&(A<E0) )|( (Y==-1)&(C-E0<A) ) )
// A_low=find( ( (E0<=A)&(A<=C-E0) )|( (Y==1)&(C-E0<A) )|( (Y==-1)&(A<E0) ) )
A_up=[I_1,I_2]
A_low=[I_3,I_4]
sa=max([max_0;S(A_low)-Y(A_low)])-min([min_0;S(A_up)-Y(A_up)])
// sa=max([max_0,max_3,max_4])-min([min_0,min_1,min_2])
//printf("sa=%f ",sa)
if ( sa<Tau ) then,fin=1,
else,fin=0,
end
endfunction
//////////////////////////////////////////////////////////////////////////////////////////
function [ers,Ps,ert,Pt]=error_rate(X,Y,x,y,A,i0C,bias)
x1=[],x3=[],Kr=0,K1=[],fs=0,cA=[],cB=[],ers=0,Ps=0
x1=X.*.ones(length(i0C),1);
x3=ones(N,1).*.X(i0C,:);
Kr=exp(-sum((x1-x3).^2,2)/(o^2));
K1=matrix(Kr,length(i0C),N);
fs=(A(i0C).*Y(i0C))'*K1+bias
cA=find(fs>0);
cB=find(fs<0);
ers=length(find(Y(cA)<>1))+length(find(Y(cB)<>-1));
Ps=ers/N*100;
x1=[],x3=[],Kr=0,K1=[],fs=0,cA=[],cB=[],ert=0,Pt=0
x1=x.*.ones(length(i0C),1);
x3=ones(Nt,1).*.X(i0C,:);
Kr=exp(-sum((x1-x3).^2,2)/(o^2));
K1=matrix(Kr,length(i0C),Nt);
fs=(A(i0C).*Y(i0C))'*K1+bias
cA=find(fs>0);
cB=find(fs<0);
ert=length(find(y(cA)<>1))+length(find(y(cB)<>-1));
Pt=ert/Nt*100;
endfunction
//////////////////////////////////////////////////////////////////////////////////////////
function [A,Ax,A0,AC,bias,ii,S]=pattern_CG(X,Y,C,N,q)
A=zeros(N,1),A_old=zeros(N,1),S=zeros(N,1),S_old=zeros(N,1),
Beta=0,p_old=zeros(N,1),//p_=zeros(q,1),
fin=0,mi=1,ii=0,sa=0,bias=0,
x2=ones(q,1).*.X,Yj=ones(q,1)*Y',
while (fin==0),//&(ii<15*N),
ii=ii+1
[WSet,p_old,r_old,q2]=WorkingSet_CG(A,Y,S,Beta,p_old,ii)
[AB,AN,YB,YN,KB,QBB,QBN]=separate(A,X,Y,WSet,q2,x2,Yj)
// [AB,AN,YB,YN,KB,QBB,QBN]=separate_1(A,X,Y,WSet,q2)
p_=-( ones(q2,1)-QBN*AN )
C2=YB'
b=-AN'*YN
ci=[zeros(q2,1)]
cs=C*[ones(q2,1)]
x0=AB
[AB]=quapro(QBB,p_,C2,b,ci,cs,mi,x0),
A(WSet)=clean(AB)
if (abs(A-A_old)<E0*ones(N,1))|(norm(r_old)<=E0*sqrt(N)) then,Beta=0,
else,Beta=1/(r_old'*r_old),
end
[S]=aboutS(S_old,KB,A,A_old,WSet)
A_old=A,S_old=S,
[fin,Ax,A0,AC,sa]=judgement(S,A,Y,Tau)
end
if length(Ax)<>0,
// bias=mean(Y-S);
bias=mean(Y(Ax)-S(Ax));
end,
endfunction
//////////////////////////////////////////////////////////////////////////////////////////
