function results = graphinference(datapath,targetnetwork,genomicdata,algonames,resultpath,filesuffix,nfolds,nrepeats,Clocal)
%
% GRAPHINFERENCE : reconstruction of network
%
% Copyright 2007, Jean-Philippe Vert
%
%   Compare several methods for the reconstruction of a networks. A .mat
%   file of results with ROC and FDR performance is produced. Additionally
%   two figures (average ROC curves and averae FDR curves are produced.
%
%   This functions calls the following functions that must be on the path:
%       rocanalysis.m
%       metriclearning.m
%       simpleSVM package
%
%   USAGE:
%   results =
%   graphinference(datapath,targetnetwork,genomicdata,algonames,filesuffix,nfolds,nrepeats,Clocal)
%
%   datapath : path to the data directory (.mat files)
%
%   targetnetwork = ['metabolic'|'interaction'] : the target network (the
%   adjacency matrix must be in the file 'adj_<targetnetwork>.mat'
%
%   genomicdata = ['exp'|'phy'|'loc'|'y2h'|'int'] : the genomic data used.
%   The kernel matric must be in the fild
%   'K_<targetnetwork>_<genomicdata>.mat'
%
%   algoname = the list of algorithms to be tested, a sublist of
%   {'naive','ml','kcca','em','local'}
%
%   resultpath = the path to the directory where resuls are written 
%
%   filesuffix = the suffix to be added to the files created
%
%   nfolds = number of folds in cross-validation
%
%   nrepeats = number of times the cross-validation is repeated
%
%   Clocal = the C parameter for the 'local' method
%
%   Example:
%   r = graphinference('~/graphinference/data/matlab','metabolic','int', {'naive','ml','kcca','em','local'} , '~/graphinference/results/' , 'f5_r1_C1000', 5,1,1000)
%

if (isdeployed)
	Clocal = str2num(Clocal);
        nfolds = str2num(nfolds);
        nrepeats = str2num(nrepeats);
        eval(strcat(['algonames = ',algonames]));
end



% Suffix for the result files
filesuffix = strcat('_',targetnetwork,'_',genomicdata,'_',filesuffix);

% Some parameters

% Ratio of negative/positive examples for large-scale SVMs
negposratio = 3;
% exponent for the diffusion kernel
betadiffusion=1;
% length of the ROC curves (sampled uniformly over [0,1])
roclength = 1000;
% parameters for KCCA
lambdaKCCA = 0.01;
nfeaturesKCCA = 20;


%%% Load the data %%%

% Load network adjacency matrix in the variable 'adj'
eval(strcat(['load ',datapath,'adj_',targetnetwork,'; adj = adj_',targetnetwork,'; clear adj_',targetnetwork]))
nvertices = size(adj,1);
sizefold = round(nvertices/nfolds);


% Load data kernel in the variable 'mykernel'
kfile = strcat('K_',targetnetwork,'_',genomicdata);
eval(strcat(['load  ',datapath,kfile,'; mykernel = ',kfile,'; clear ',kfile]))

%%% Run the experiments %%%


nexp = 1;


for repeat = 1:nrepeats % Main loop (repeat it nrepeats times)

    % Randomize the data in the folds
    randomizer = randperm(nvertices);

    for fold = 1:nfolds % Cross-validation loop (repeat it nfolds times)

        disp(['Now starting iteration ',int2str(nexp),' out of ',int2str(nfolds*nrepeats)])

        % Create training and test sets
        if fold < nfolds
            TrainSet = [randomizer(1:(fold-1)*sizefold),randomizer(1+(fold*sizefold):end)];
            TestSet = randomizer((1+(fold-1)*sizefold):(fold*sizefold));
        else
            TrainSet = randomizer(1:(fold-1)*sizefold);
            TestSet = randomizer((1+(fold-1)*sizefold):end);
        end
        TrainSetSize = size(TrainSet,2);
        TestSetSize = size(TestSet,2);

        % Here we focus on the prediction of edges between one node in the
        % training set and one node in the test set.
        target_adj = adj(TrainSet,TestSet);


        % Some methods need the Laplacian or the diffusion kernel on the
        % training set only:

        % Laplace matrix of the training graph
        train_lap = adj(TrainSet,TrainSet);
        train_lap = diag(sum(train_lap)) - train_lap;
        % Normalize it to diagonal 1
        train_lap = train_lap ./ (sqrt(diag(train_lap))*sqrt(diag(train_lap))' + 0.00000001);
        % take the diffusion kernel and normalize it to diagonal 1
        train_diffusion_kernel = expm( - betadiffusion * train_lap );
        train_diffusion_kernel = train_diffusion_kernel ./ (sqrt(diag(train_diffusion_kernel))*sqrt(diag(train_diffusion_kernel))' + 0.00000001);


        for algo = algonames
            disp(['Now testing the algorithm ',algo])

            % Method 1 : Naive unsupervised method
            % Just predict an edge when the kernel is large
            if strcmp(algo,'naive')

                pred = mykernel(TrainSet,TestSet);
                [roc , fdr , auc , auf] = rocanalysis(target_adj(:),pred(:),roclength);
                results.naive.auc(nexp) = auc;
                results.naive.auf(nexp) = auf;
                results.naive.roc(nexp,:) = roc;
                results.naive.fdr(nexp,:) = fdr;
            end




            % Method 2 : Kernel CCA method (Yamanishi and Vert, 2003)
            if strcmp(algo,'kcca')
                % Perform kCCA
                A = [zeros(TrainSetSize,TrainSetSize) , mykernel(TrainSet,TrainSet)*train_diffusion_kernel ; train_diffusion_kernel * mykernel(TrainSet,TrainSet) , zeros(TrainSetSize,TrainSetSize)];
                B = [ (mykernel(TrainSet,TrainSet) + lambdaKCCA*eye(TrainSetSize))^2 , zeros(TrainSetSize,TrainSetSize) ; zeros(TrainSetSize,TrainSetSize) , ( train_diffusion_kernel + lambdaKCCA*eye(TrainSetSize))^2];
                [alpha,rho] = eig(A,B);
                % extract the dual features
                [m,i] = sort(diag(rho));
                alpha = alpha(1:TrainSetSize,i(1:nfeaturesKCCA));
                % Compute the feature for all points
                f = mykernel(:,TrainSet) * alpha;
                % Prediction by inner products
                pred = f(TrainSet,:) * f(TestSet,:)';
                % Estimate accuracy
                [roc , fdr , auc , auf] = rocanalysis(target_adj(:),pred(:),roclength);
                results.kcca.auc(nexp) = auc;
                results.kcca.auf(nexp) = auf;
                results.kcca.roc(nexp,:) = roc;
                results.kcca.fdr(nexp,:) = fdr;
            end


            % Method 3: Distance metric learning method (Vert and Yamanishi 2004)
            % Note :  here we take default parameters as suggested in the
            % paper, but they should be optimized by CV on the training set
            % (todo...)
            if strcmp(algo,'ml')
                newdistances = metriclearning(adj(TrainSet,TrainSet),mykernel([TrainSet,TestSet],TrainSet),2,10);
                pred = -newdistances(1:TrainSetSize,(TrainSetSize+1):end);
                [roc , fdr , auc , auf] = rocanalysis(target_adj(:),pred(:),roclength);
                results.ml.auc(nexp) = auc;
                results.ml.auf(nexp) = auf;
                results.ml.roc(nexp,:) = roc;
                results.ml.fdr(nexp,:) = fdr;
            end


            % Method 4: EM method of Kato et al. (Bioinformatics
            % 21(10):2488:2495, 2005).
            % Make prediction with em projection
            if strcmp(algo,'em')
                pred = train_diffusion_kernel * inv(mykernel(TrainSet,TrainSet)+eye(TrainSetSize)) * mykernel(TrainSet,TestSet);
                [roc , fdr , auc , auf] = rocanalysis(target_adj(:),pred(:),roclength);
                results.em.auc(nexp) = auc;
                results.em.auf(nexp) = auf;
                results.em.roc(nexp,:) = roc;
                results.em.fdr(nexp,:) = fdr;
            end



            % Method 5 : local method (the new one!)

            if strcmp(algo,'local')
                d = kernel('custom',mykernel);
                pred = zeros(TrainSetSize,TestSetSize);
                j=1;

                for i=TrainSet %train successively on all training points
                    trainindices = TrainSet([1:(j-1),(j+1):TrainSetSize]);
                    trainlabels = sign(adj(i,trainindices)-0.5)';
                    npos = size(find(trainlabels == 1),1);
                    nneg = size(find(trainlabels == -1),1);
                    Cbal = Clocal*ones(1,size(trainlabels,1));
                    Cbal(find(trainlabels == 1)) = Cbal(find(trainlabels == 1)) * nneg / (npos+nneg);
                    Cbal(find(trainlabels == -1)) = Cbal(find(trainlabels == -1)) * npos / (npos+nneg);
                    pa = param(Cbal',50,'binary','chol');
                    mydata = data(trainindices,trainlabels,TestSet,sign(adj(i,TestSet)-0.5)');
                    clear global svModel;
                    global svModel;
                    trainSVM(mydata,d,pa);
                    pred(j,:) = testSVM';
                    j=j+1;
                end
                [roc , fdr , auc , auf] = rocanalysis(target_adj(:),pred(:),roclength);
                results.local.auc(nexp) = auc;
                results.local.auf(nexp) = auf;
                results.local.roc(nexp,:) = roc;
                results.local.fdr(nexp,:) = fdr;
            end

            % Method 6: tensor product pairwise kernel (Ben-Hur and Noble)
            % This method requires a large-scale SVM implementation.


        end
        nexp = nexp+1;
    end
end

% save results
eval(char(strcat(['save ',resultpath,'results_',filesuffix,'.mat results'])))


% Print results in the file results.txt
eval(char(strcat('fid = fopen(''',resultpath,'results',filesuffix,'.txt'',''w'')')));
fprintf(fid,'AUC and AUF results on %d folds cross-validation repeated %d times\n\n',nfolds,nrepeats);
for nn=algonames
    eval(char(strcat('fprintf(fid,''', nn , '\t%2.2f \\pm %2.2f\t%2.2f \\pm %2.2f\n'',100*mean(results.' , nn , '.auc),100*std(results.' , nn ,'.auc),100*mean(results.' , nn , '.auf),100*std(results.' , nn ,'.auf)' , ')')));
end
fclose(fid);

% Plot ROC curves (of the last fold)

colors = {'b','g','r','c','k','m','y'};
clf
hold on
i=1;
for nn=algonames
    eval(char(strcat('plot([0:1/(roclength-1):1],mean(results.',nn,'.roc),''',colors(i),''')')));
    i=i+1;
end
legend(algonames,'Location','Southeast')
title('ROC curves')
xlabel('Ratio of false positives')
ylabel('Ratio of true positives')
hold off
eval(char(strcat(['print -depsc ',resultpath,'roc',filesuffix]))) % output roc.eps
eval(char(strcat(['print -dpng ',resultpath,'roc',filesuffix]))) % output roc.pdf

% Plot FDR curves
clf
hold on
i=1;
for nn=algonames
    eval(char(strcat('plot([0:1/(roclength-1):1],mean(results.',nn,'.fdr),''',colors(i),''')')));
    i=i+1;
end
legend(algonames,'Location','Southeast')
title('FDR curves')
xlabel('Ratio of true positives')
ylabel('False discovery rate')
hold off
eval(char(strcat(['print -depsc ',resultpath,'fdr',filesuffix]))) % output roc.eps
eval(char(strcat(['print -dpng ',resultpath,'fdr',filesuffix]))) % output roc.pdf


% Exit MATLAB (when the script is launched automatically from a shell)
%exit
