% written by Renqiang Min for generating new kernels using label information
% searching a good alpha using a validation set
addpath(spiderpath);
use_spider;

res=[]; res50=[]; 
rand('seed',1); randn('seed',1);
if restart==1
  load prof-kernel.txt;
  KK=prof_kernel;
  clear prof_kernel;
  k2=sqrt(diag(KK))+1e-12; k2=k2*k2'; KK=KK ./ k2; clear k2;
  Ys=load('allseq_membership-noheaders.table'); % load data splits
  load psiblast_homologs addpts;
  name='prof_hom_labelinfo2.mat'; % name of file to save
end;

% loop over data splits
for i=1:54
    global K; global X; global Y; 
    % allocate correct train and test set labels for this data split
    YY=Ys(:,i); 
    tmpY=YY;
    % test & test set sizes
    l=sum(YY==1)+sum(2==YY); lt=sum(YY==3)+sum(4==YY); 
    [trn]=find(YY==1 | YY==2); [tst]=find(YY==3 | YY==4);
    YY(trn(1:l))=sign((YY(trn(1:l))==1)-0.5); 
    YY(tst(1:lt))=sign((YY(tst(1:lt))==3)-0.5);

    %%--------------- pull in domain homologs -----------
    f=addpts{i};  % already cached these computations
    trn=[trn ;f']; trn=unique(trn); YY(f)=1;
    %% ------------------------------------------------
    

    % search a good alpha by training and validation
    % divide the original training set into a training set and a validation set to
    % decide a good value for alpha
    trpos = trn(find(YY(trn)==1));
    trneg = trn(find(YY(trn)==-1));
    
    numtrpos=size(trpos,1);
    if numtrpos>=5
       numfold=5;
    else
       numfold=2;
    end;
    perfoldpos = numtrpos/numfold;
    perfoldneg = size(trneg,1)/numfold;
    trpos=trpos(randperm(size(trpos,1)));
    trneg=trneg(randperm(size(trneg,1)));
    result50=[];
    alphalist =[1.5 2 2.5 3 3.5 7 9 100]; 
for cv=1:numfold
    fprintf(1,'cross validation round %d\n',cv);
    tmptrpos = trpos([1:floor((cv-1)*perfoldpos),floor(cv*perfoldpos)+1:end]);
    tmpttpos = trpos(floor((cv-1)*perfoldpos)+1:floor(cv*perfoldpos));
    tmptrneg = trneg([1:floor((cv-1)*perfoldneg),floor(cv*perfoldneg)+1:end]);;
    tmpttneg = trneg(floor((cv-1)*perfoldneg)+1:floor(cv*perfoldneg));
   

    trind = [tmptrpos;tmptrneg]; 
    ttind = [tmpttpos;tmpttneg];
    tmpKtr= KK(trind,trind);
    tmpKtt= KK(ttind,trind);

    tmpnumtr = size(tmpKtr,1);
    tmpnumtt = size(tmpKtt,1);
    tmpn = tmpnumtr + tmpnumtt;
    Y(1:tmpnumtr,1) = YY(trind,1);
    Y(tmpnumtr+1:tmpn,1) =YY(ttind,1);
    X = Y;           
    newtrn=[1:tmpnumtr];
    newtst=[tmpnumtr+1:tmpn];



    for tmpi = 1:size(alphalist,2)
        alpha=alphalist(tmpi); 
        tmpKtr2 = tmpKtr;

        
        % searching a good alpha
        
        for ii = 1:tmpnumtr
            for jj = 1:tmpnumtr
                if YY(trind(ii))~=YY(trind(jj))
                    tmpKtr2(ii,jj) = tmpKtr2(ii,jj)/alpha;
                end;
            end;
        end;
        
 
        Kextra=tmpKtt * inv(tmpKtr) * tmpKtr2;
        K = ones(tmpn, tmpn);
        K(1:tmpnumtr,1:tmpnumtr)=tmpKtr2;
        K(tmpnumtr+1:tmpn, 1:tmpnumtr)=Kextra;
        K(1:tmpnumtr,tmpnumtr+1:tmpn)=Kextra';
        
      
        d=data_global; 
        dtrn=get(d,newtrn); 
        dtst=get(d,newtst);  
        
        % now train svm
        a=svm; 
        a.optimizer='svmlight'; 
        a.child=kernel('custom_fast'); 
        a.verbosity=0; 
        a.use_signed_output=0;
        a.balanced_ridge=0.02*median(diag(K));
        [tr a]=train(a,dtrn); 
        r=test(a,dtst);
        
        % roc score      
        rr=loss(r,'roc50'); 
        result50(cv,tmpi)=rr.Y;
        K=[];
        clear a;
    end;
end;
    result50
    meanres50=mean(result50,1)
    [maxv tmpi]=max(meanres50);
    crossval(i,:)=meanres50;

    % Use one of the following two options
    % beginOption
  
    % the first option
    %tmpi=max(find(meanres50==maxv));

    % the second option
    tmpi=find(meanres50==maxv);
    indtmpi=size(tmpi,2);
    while((indtmpi>1) & (tmpi(indtmpi)-tmpi(indtmpi-1)>1))
       indtmpi=indtmpi-1;
    end;
    tmpi=tmpi(indtmpi);
    clear indtmpi;
    % endOption

    Alpha=alphalist(tmpi)
    Alphalist(i)=Alpha;
    clear a r rr dtst dtrn newtrn newtst d tmpnumtr tmpnumtt tmpn ...
            tmpKtr tmpKtt tmpKtr2 Kextra;
    X=[];
    Y=[]; 

    % traning using new kernel matrix based on the previously decided alpha

    Ktr = KK(trn, trn);
    Ktt = KK(tst, trn);
 
    Ktr2=Ktr;
    numtr = size(Ktr,1);
    numtt = size(Ktt,1);
    n = numtr+numtt;
    for ii = 1:numtr
        for jj = 1:numtr
            if YY(trn(ii))~=YY(trn(jj))
                Ktr2(ii,jj) = Ktr2(ii,jj)/Alpha;
            end;
        end;
    end;

    Kextra=Ktt * inv(Ktr) * Ktr2;
    K = ones(n, n);
    K(1:numtr,1:numtr)=Ktr2;
    K(numtr+1:n, 1:numtr)=Kextra;
    K(1:numtr,numtr+1:n)=Kextra';
    
    
    Y(1:numtr,1) = YY(trn');
    Y(numtr+1:n,1) =YY(tst');
    X = Y;
    d=data_global; 
    newtrn=[1:numtr];
    newtst=[numtr+1:n];
    dtrn=get(d,newtrn); 
    dtst=get(d,newtst);  
    
    % now train svm
    a=svm; 
    a.optimizer='svmlight'; 
    a.child=kernel('custom_fast'); 
    a.verbosity=2; 
    a.use_signed_output=0;
    a.balanced_ridge=0.02*median(diag(K));
    [tr a]=train(a,dtrn); 
    r=test(a,dtst);
    
    % roc score
    rr=loss(r,'roc'); 
    res=[res rr.Y];         
    rr=loss(r,'roc50'); 
    res50=[res50 rr.Y]
    clear dtst dtrn Ktr Kextra V1 D1 V2 D2 firstTerm W2tt W2;
end;
save(name,'res','res50','alphalist','Alphalist','crossval');


