%function [alpha,c,opt,b,flambda,mygckl,mygacv,dhat,trainMis, theoMis] = solve(mylambda, mysigma)
%function [mygckl,mygacv,rangacv,jxa1,jxa2,theoMis,trainMis] = solve(mylambda,mysigma); 
function [mygckl,mygacv,jxa1,jxa2,theoMis] = solve(mylambda,mysigma); 

global n;
global Nodes1;
global y;
global d;
global K;

%global noise;
%global perty;

%for  i =1:n
% for j=1:n
%      tmp = sum((Nodes1(i,:)-Nodes1(j,:)).^2);
%      K(i,j)= exp(-tmp/(2*mysigma*mysigma)); 
%  end
%end


scaling = 1e+4;
d1 = sqrt(8/pi - d^2)/2 - d/2;  
d2 = d1 + d;

%p1 = 1 - pi*d2^2/4;
%p2 = pi*(d2^2-d1^2)/4;
%p3 = pi*d1^2/4;
%vary = 1;
%varz = noise*noise;

alpha = zeros(n,1);
tol = 1e-10;


Y = diag(y);
H = Y*K*Y/(2*n*mylambda);
H = scaling*H;
p = -ones(n,1);
p = scaling*p;
A = [];
b1 = [];

Aeq = y';
beq = 0;
lb = zeros(n,1);
ub = ones(n,1);
x0 = zeros(n,1);    %starting point


%options = optimset('MaxIter',1000000,'LargeScale','off','DerivativeCheck','on','Diagnostics','on' );
options = optimset('MaxIter',1000000);
alpha = quadprog(H,p,A,b1,Aeq,beq,lb,ub,x0,options);
opt = 0.5*alpha'*H*alpha + p'*alpha;

c = Y*alpha/(2*n*mylambda);

Aph = diag(alpha);
b = ones(n,1)'*Aph*(eye(n)-Aph)*(y-K*c)/(sum(alpha.*(ones(n,1)-alpha)));

flambda = K*c+ones(n,1)*b;
%trainMis = sum(flambda.*y<tol)/n + sum(flambda.*y == 0)/(2*n);


%Compute mygckl 
mygckl = 0;
mygckl = gckl_est(flambda);


%Compute mygacv 
%compute dhat
dhat =0;
for i=1:n
   if y(i,1)*flambda(i,1) <= -1 + tol 
     dhat = dhat + alpha(i,1)/(n*mylambda);
   elseif y(i,1)*flambda(i,1) <= 1 + tol
     dhat = dhat + alpha(i,1)/(2*n*mylambda);
   end
end
dhat = dhat/n;
%compute obs, using a_+ = (a+abs(a))/2
xi = (ones(n,1) - y.*flambda + abs(ones(n,1) - y.*flambda))/2;
obs = sum(xi)/n;
%compute mygacv
mygacv = obs + dhat;


%Compute jxa1 and jxa2
jxa1 = sum(alpha/(2*n*mylambda) + xi - ones(n,1) > tol)/n;
jxa2 = sum(alpha/(n*mylambda) + xi - ones(n,1) > tol)/n;


%****** Compute rangacv
%ee = (y.*flambda <= -1+tol)+0.5*(y.*flambda>-1+tol & y.*flambda<=1+tol )+0*(y.*flambda>1+tol);
%newY = diag(perty);
%newH = newY*K*newY/(2*n*mylambda);
%newH = (newH + newH')/2;
%newH = scaling*newH;
%newp = -perty.*perty;
%newp = scaling*newp;
%newAeq = perty';

%options = optimset('MaxIter',10000,'LargeScale','off','DerivativeCheck','on','Diagnostics','on');
%options = optimset('MaxIter',10000);
%newalpha = quadprog(newH,newp,A,b1,newAeq,beq,lb,ub,alpha,options);
%newopt = 0.5*newalpha'*newH*newalpha + newp'*newalpha;

%newc = newY*newalpha/(2*n*mylambda);
%newAph = diag(newalpha);
%newb = ones(n,1)'*newAph*(eye(n)-newAph)*(perty-K*newc)/(sum(newalpha.*(ones(n,1)-newalpha)));

%newflambda = K*newc+ones(n,1)*newb;
%p1 = sum(ones(n,1)-y.*flambda + abs(ones(n,1)-y.*flambda))/(2*n);
%p2 = sum((perty-y).*ee.*(newflambda-flambda))*2/(n*varz);
%rangacv = p1+p2;


theoMis = 0;
for i = 1:n
  if  Nodes1(i,1)^2 + Nodes1(i,2)^2  >=  d2^2
     theoMis = theoMis + 0.05 * (flambda(i,1)>0) + 0.95 * (flambda(i,1)<0);
  elseif  Nodes1(i,1)^2 + Nodes1(i,2)^2  <=  d1^2
     theoMis = theoMis + 0.95 * (flambda(i,1)>0) + 0.05 * (flambda(i,1)<0);
  else  
      theoMis = theoMis + 0.5;
  end
end
theoMis = theoMis/n;


diary middle;
%[log2(mylambda) log2(mysigma) mygckl mygacv rangacv jxa1 jxa2 theoMis trainMis]
[log2(mylambda) log2(mysigma) log2(mygckl) log2(mygacv) log2(jxa1) log2(jxa2) log2(theoMis)]
diary off;



%draw_cont(mysigma,c,b,Nodes1,y);












