% login to bleecker and type
% export LD_LIBRARY_PATH=/local/mosek/mosek.4/4/tools/platform/linux32x86/bin/
% export MOSEKLM_LICENSE_FILE=/local/mosek/mosek.4/4/licenses/mosek.lic.orig
% export MATLABPATH=/local/mosek/mosek.4/4/toolbox/r14/
% matlab


addpath /homes/mosek/tools/platform/linuxintel/bin
addpath /homes/mosek/toolbox/solversr12
mosekopt
help quadprog
H=eye(3,3);f=[1;2;3];A=eye(3,3);b=[4;7;8];quadprog(H,f,A,b)

load LandmineData.mat
M = size(feature,2);
D = size(feature{1},2);
T = zeros(M,1);
for m=1:M
  T(m)=length(label{m});
end

maxF  = 10;
maxN  = 10;
maxC  = 10;
maxA  = 10;


aucTR = zeros(maxF,maxN,maxC,maxA,M);
aucCV = zeros(maxF,maxN,maxC,maxA,M);
aucTE = zeros(maxF,maxN,maxC,maxA,M);
timeX = zeros(maxF,maxN,maxC,maxA);
runsX = zeros(maxF,maxN,maxC,maxA);

for indF=1:maxF

  % After the first fold, randomly re-order the data
  ordering  = label;
  if (indF<2)
    for m=1:M
      ordering{m} = 1:T(m);
    end
  else
    for m=1:M
      [qs,qI]=sort(rand(T(m),1));
      ordering{m} = qI;
    end
  end

  labelTR   = label;
  labelTE   = label;
  kernelsTR = cell(M,D);
  kernelsTE = cell(M,D);

  for indN=1:maxN
    N = 20*indN;
    for m=1:M
      ord1 = ordering{m};
      lab1 = label{m};
      labelTR{m} = 2*(lab1(ord1(1:N))-0.5);
      labelTE{m} = 2*(lab1(ord1(N+1:end))-0.5);
      X   = feature{m};
      xTR = X(ord1(1:N),:);
      xTE = X(ord1(N+1:end),:); 
      for d=1:D
        kernelsTR{m,d} = xTR(:,d)*xTR(:,d)';
        kernelsTE{m,d} = xTE(:,d)*xTR(:,d)';
      end
    end

    for indC=1:maxC
      C=(sqrt(10))^indC;
      lambdas = labelTR;
      for m=1:M
        lambdas{m} = zeros(N,1);
      end
      for indA=1:maxA
        alpha=(sqrt(25))^(indA-1)-1;
        timenow = cputime;
        fprintf('Fold=%d N=%d C=%e alpha=%e time=%e\n',indF,N,C,alpha,timenow);
        [lambdas,y,S,runsnow] = multisvm(kernelsTR,labelTR,C,alpha,lambdas);
        timenow = cputime-timenow;
        timeX(indF,indN,indC,indA) = timenow;
        runsX(indF,indN,indC,indA) = runsnow;

        for indM=1:M
          aTR = auc(y{indM},((labelTR{indM}+1)/2));
          aucTR(indF,indN,indC,indA,indM) = aTR;
          lab1 = (labelTE{indM}+1)/2;
          half = floor(length(lab1)/2);
          y1  = zeros(length(lab1),1);
          for d=1:D
            y1  = y1+S(d)*(kernelsTE{indM,d}*(labelTR{indM}.*lambdas{indM}));
          end
          aCV = auc(y1(1:half),lab1(1:half));
          aucCV(indF,indN,indC,indA,indM) = aCV;
          aTE = auc(y1(half+1:end),lab1(half+1:end));
          aucTE(indF,indN,indC,indA,indM) = aTE;
          fprintf('Fold=%d N=%d C=%e alpha=%e M=%d TR=%f CV=%f TE=%f\n',indF,N,C,alpha,indM,aTR,aCV,aTE);
        end

        save aucTR.mat aucTR
        save aucCV.mat aucCV
        save aucTE.mat aucTE
        save timeX.mat timeX
        save runsX.mat runsX

      end
    end
  end
end

save aucTR.mat aucTR
save aucCV.mat aucCV
save aucTE.mat aucTE
save timeX.mat timeX
save runsX.mat runsX


scoreMETA = zeros(maxF,maxN);
scoreSVM  = zeros(maxF,maxN);

for indF=1:maxF
  for indN=1:maxN
    bestC = 1;
    bestA = 1;
    cross = 0;
    for indC=1:maxC
      for indA=1:maxA
        if (mean(aucCV(indF,indN,indC,indA,:))>cross)
          bestC=indC;
          bestA=indA;
          cross=mean(aucCV(indF,indN,indC,indA,:));
        end
      end
    end
    scoreMETA(indF,indN) = mean(aucTE(indF,indN,bestC,bestA,:));

    bestC = 1;    
    bestA = 1;
    cross = 0;
    for indC=1:maxC
      for indA=1:1
        if (mean(aucCV(indF,indN,indC,indA,:))>cross)
          bestC=indC;
          bestA=indA;
          cross=mean(aucCV(indF,indN,indC,indA,:));
        end
      end
    end
    scoreSVM(indF,indN) = mean(aucTE(indF,indN,bestC,bestA,:));
  end
end

muMETA = zeros(maxN,1);
cvMETA = zeros(maxN,1);
muSVM  = zeros(maxN,1);
cvSVM  = zeros(maxN,1);
axisN  = zeros(maxN,1);
for indN=1:maxN
  axisN(indN)  = 10*indN;
  muMETA(indN) = mean(scoreMETA(:,indN));
  cvMETA(indN) = sqrt(cov(scoreMETA(:,indN)));
  muSVM(indN)  = mean(scoreSVM(:,indN));
  cvSVM(indN)  = sqrt(cov(scoreSVM(:,indN)));
end

clf
errorbar(axisN,muMETA,cvMETA);
hold on
errorbar(axisN,muSVM,cvSVM,'r--');
h=gca;
set(h,'FontSize',16);
xlabel('Training Set Size per Task');
ylabel('Average AUC Across Tasks');
legend('Multi-Task Sparsity','Independent SVMs');
print -depsc average.eps



