function alpha=svm_softassign
clear
figure(1)
clf
hold off;
figure(2)
clf
hold off;
figure(3)
clf
hold off;
n=100;
%Initializing the patterns (x) and outputs (y)
x=rand(n,2);
frac1=0.8;
xsum=sum(x')';
xminusone=find(xsum < 1);
x(xminusone)=x(xminusone).*frac1;
xplusone=find(xsum >= 1);
x(xplusone)=x(xplusone)./frac1;
y=ones(n,1);
y(xminusone)=-1*ones(length(xminusone),1);
%Initializing the free parameters
beta=0.5;
gamma=0.5;
anneal=0.98;
C=100;
deltaW=1;
deltaW_thr=0.5;
constraint_threshold=0.01;
alpha_thr=3;
alpha_init=1.0;
t1_min=10;
%Computing H
H=(x*x').*(y*y');
%Initializing alpha
alpha=alpha_init*ones(n,1);
sqrt_ratio=sqrt(sum(alpha(xplusone))/sum(alpha(xminusone)));
alpha(xminusone)=alpha(xminusone).*sqrt_ratio;
alpha(xplusone)=alpha(xplusone)./sqrt_ratio;
slack=C-alpha;
t1=0;
t2=0;
sumexpalpha_prod=ones(n,1);
while (deltaW >= deltaW_thr | t1 < t1_min),
   gamma=gamma*anneal;
   W_before=0.5*alpha'*H*alpha-ones(n,1)'*alpha;
   expalpha=alpha.*exp(-beta*(H*alpha-1))./sumexpalpha_prod; 
   expslack=slack./sumexpalpha_prod;
   constraint_err=1.0;
   inner_loop=0;
   while constraint_err >= constraint_threshold
	 sumexpalpha=(expalpha+expslack)/C;
	 sumexpalpha_prod=sumexpalpha_prod.*sumexpalpha;
	 expalpha=expalpha./sumexpalpha;
	 expslack=expslack./sumexpalpha;
	 sqrt_ratio=sqrt(sum(expalpha(xplusone))/sum(expalpha(xminusone)));
	 expalpha(xminusone)=expalpha(xminusone).*sqrt_ratio;
	 expalpha(xplusone)=expalpha(xplusone)./sqrt_ratio;
	 constraint_err=sum((expalpha+expslack-C).^2);
	 figure(3)
	 plot(t2,constraint_err,'x',t2,constraint_threshold,'-');
	 hold on;
	 t2=t2+1;
     inner_loop = inner_loop +1;
   end 
   alpha=alpha*(1-gamma)+gamma*expalpha;
   slack=slack*(1-gamma)+gamma*expslack;
   W=0.5*alpha'*H*alpha-ones(n,1)'*alpha;
   deltaW=W_before-W;
   disp(sprintf('W = %g deltaW = %g gamma = %g inner_loop = %d outer_loop = %d\n',W,deltaW, gamma, inner_loop, t1));
   figure(1)
   plot(t1,W);
   hold on;
   figure(2)
   plot(x(xminusone,1),x(xminusone,2),'x')
   hold on
   plot(x(xplusone,1),x(xplusone,2),'o')
   alphasupport=find(alpha > alpha_thr);
   plot(x(alphasupport,1),x(alphasupport,2),'r*');
   hold off
   t1=t1+1;
end
