
function p4_1(mmtrain, nttrain, mmtest, nttest),

function testcoeffs = pcaproject(basisset, testset, NBASIS)
   % returns the PCA projection of testset on to the 
   % NBASIS principal components of basisset

   basisset = basisset - repmat(mean(basisset,2),[1 ...
	 size(basisset,2)]);
   testset = testset - repmat(mean(basisset,2),[1 ...
	 size(testset,2)]);
   basiscov = cov(basisset');
   [U,S,V] = svd(basiscov,'econ');

   basis = U(:,[1:NBASIS]);
   testcoeffs = (testset'*basis)';
end

function err = classify(svm, data, labels)

   [outlabels, dvs] = ...
    SVMClass(data,svm.AlphaY,svm.SVs,svm.Bias,svm.Parameters);

   err = (sum(abs(outlabels-labels)))/length(labels);
end;

function paramvec = makeParamVec(svmparams), 
paramvec = [svmparams.KERNEL_TYPE svmparams.DEGREE svmparams.GAMMA ...
      svmparams.COEFFICIENT svmparams.C svmparams.CACHE ...
      svmparams.EPSILON svmparams.SVM_TYPE];

end;

% first perform pca on the training set

NBASIS = 2;
svmparams.KERNEL_TYPE = 1;
svmparams.DEGREE = 3;
svmparams.GAMMA = 1;
svmparams.COEFFICIENT = 0;
svmparams.C = 1000;
svmparams.CACHE = 40000000;
svmparams.EPSILON = 0.001;
svmparams.SVM_TYPE = 0;

trainset = [mmtrain nttrain];
trainsetlabels = [ones(1,size(mmtrain,2)) -1*ones(1,size(nttrain,2))];
traincoeffs = pcaproject(trainset, trainset, NBASIS);

testset = [mmtest nttest];
testsetlabels = [ones(1,size(mmtest,2)) -1*ones(1,size(nttest,2))];
testcoeffs = pcaproject(trainset, testset, NBASIS);

% do C test 
Cs = [1 10 100 1000];
trainerrCs = [];
testerrCs = [];
trainerrRBF = [];
testerrRBF = [];
figure_i = 1;

for c=Cs,
  
  sprintf('POLY %d',c)
  svmparams.C = c;
  
  svmparams.KERNEL_TYPE = 1;
  svmparams.DEGREE = 3;
  % train the svm
  [svm.AlphaY, svm.SVs, svm.Bias, svm.Parameters, svm.nSV, svm.nLabel] = ...
      SVMTrain(traincoeffs,trainsetlabels,makeParamVec(svmparams), ...
      1.0);

  figure(figure_i);
  figure_i = figure_i+1;
  svm_plot(traincoeffs,trainsetlabels,svm.AlphaY,svm.SVs,svm.Bias,svm.Parameters,svm.nSV,svm.nLabel);
  title(sprintf('svm decision boundary for Kernel type %d C=%d',svmparams.KERNEL_TYPE,svmparams.C));
  
  % plot 
  
  trainerrCs = [trainerrCs classify(svm,traincoeffs,trainsetlabels)];
  testerrCs =  [testerrCs classify(svm,testcoeffs,testsetlabels)];
  
  

  sprintf('RBF %d',c)

  svmparams.KERNEL_TYPE = 3;
  [svm.AlphaY, svm.SVs, svm.Bias, svm.Parameters, nSV, svm.nLabel] = ...
      SVMTrain(traincoeffs,trainsetlabels,makeParamVec(svmparams), ...
      1.0);
  
  figure(figure_i);
  figure_i = figure_i+1;  
  svm_plot(traincoeffs,trainsetlabels,svm.AlphaY,svm.SVs,svm.Bias,svm.Parameters,svm.nSV,svm.nLabel);
  title(sprintf('svm decision boundary for Kernel type %d C=%d',svmparams.KERNEL_TYPE,svmparams.C));
  
  trainerrRBF = [trainerrRBF classify(svm,traincoeffs,trainsetlabels)];
  testerrRBF =  [testerrRBF classify(svm,testcoeffs,testsetlabels)];

end;

size(trainerrCs)

figure(figure_i);
figure_i = figure_i+1;
semilogx(Cs,[trainerrCs' testerrCs' trainerrRBF' testerrRBF']);
title('training and test error as a function of C');

%% gammas
Gammas = [0.001 1000];
trainerrGs = [];
testerrGs = [];
svmparams.C = 1; % fix C at 1 for this test
svmparams.KERNEL_TYPE=1;
svmparams.DEGREE=3;
for gamma = Gammas,
  
  svmparams.GAMMA = gamma;
  sprint('POLY GAMMA %d',gamma);
  % train the svm
  [svm.AlphaY, svm.SVs, svm.Bias, svm.Parameters, nSV, svm.nLabel] = ...
      SVMTrain(traincoeffs,trainsetlabels,makeParamVec(svmparams), ...
      1.0);

  trainerrGs = [trainerrGs classify(svm,traincoeffs,trainsetlabels)];
  testerrGs =  [testerrGs classify(svm,testcoeffs,testsetlabels)];
end;

figure(2);
semilogx(Gammas,[trainerrGs' testerrGs']);
title('training and test error as a function of Gamma');

end
