clear all
clc

tic

 load dataset.mat; %Partition Dataset for the 5 fold test


% rng('default');     % reset random generator.
opts.tFlag = 1;     % terminate after relative objective value does not changes much.
opts.tol = 10^-5;   % tolerance.  
opts.maxIter = 1000; % maximum iteration number of optimization.


for n=1:5
    test_VBM=Xt_VBM{1,n};
     test_CC=Xt_CC{1,n};
     test_Yt=Yt_SNP{1,n}(:,1); 

    task.DT{1}=X_VBM{1,n};
    task.DT{2}=X_CC{1,n};
    respons=Y_SNP{1,n}(:,1);
    task.target{1}=Y_SNP{1,n}(:,1); %1:rs429358
    task.target{2}=Y_SNP{1,n}(:,1); %2:rs429358
    task.lab{1}=Y{1,n};
    task.lab{2}=Y{1,n};
    gnd=task.lab{1};

    task.num=2;

 


 paraset=[0.0000001 0.0000003 0.000001 0.000003 0.00001];





    for j=1:length(paraset)
        opts.rho1=paraset(j);
        for k=1:length(paraset)
            opts.rho_L3=paraset(k);
            opts.init =2; % guess start point from data ZERO.
            kfold=5;
            kk=1;     
            % construct the index of cross_validation for each task.
            [tcv fcv]=f_myCV(gnd',kfold,kk); 
            %% begin to 5-fold.
            for cc=1:kfold 
                task.X = cell(task.num,1);
                task.Y = cell(task.num,1);
                for i=1:task.num                  
                    trLab=tcv{cc}';
                    % generate the task.              
                    task.X{i}=task.DT{i}(trLab,:);               
                    task.Y{i}=task.target{i}(trLab);
                    task.label{i}=task.lab{i}(trLab);
                end
                %----------Main Algorithm---------------        

                 
                  [S, D] = f_lapMatrix(task.Y{1});
                   [W, epsvalue] = f_MTM_APG(task.X,task.Y,opts,S, D);
                % find the selected features for each task. 
                trLab=tcv{cc}';
                teLab=fcv{cc}';
                pl=task.DT{1}(teLab,:)*W(:,1);
                pl2=task.DT{2}(teLab,:)*W(:,2);      

                et(cc)=sqrt(mean((pl-respons(teLab,1)).^2))+sqrt(mean((pl2-respons(teLab,1)).^2));                                    
                a=respons(teLab,1)-mean(respons(teLab,1));b=pl-mean(pl);
                a2=respons(teLab,1)-mean(respons(teLab,1));b2=pl2-mean(pl2);
              
                co(cc)=abs(sum(a.*b)/sqrt(sum(a.^2)*sum(b.^2)))+abs(sum(a2.*b2)/sqrt(sum(a2.^2)*sum(b2.^2)));%+abs(sum(a3.*b3)/sqrt(sum(a3.^2)*sum(b3.^2)))+abs(sum(a4.*b4)/sqrt(sum(a4.^2)*sum(b4.^2)));
            end       
            res_kfold_CO(kk)=mean(co);
            res_kfold_RMSE(kk)=mean(et);  
            res_CO(j,k)=mean(res_kfold_CO);
            res_RMSE(j,k)=mean(res_kfold_RMSE);
        end
    end
    ndim=size(res_RMSE);
    tempRMSE=10;
    tempCO=0;
    for ii=1:ndim(1)
        for jj=1:ndim(2)
               if  res_CO(ii,jj)>tempCO
                 tempCO=res_CO(ii,jj);
                 paraSet=[ii,jj];
               end
        end
    end
    paraSet

    test_opts.rho1=paraset(paraSet(1));
    test_opts.rho_L3=paraset(paraSet(2));
    Final_para{n}=paraSet;
    test_opts.tFlag = 1;     % terminate after relative objective value does not changes much.
    test_opts.tol = 10^-5;   % tolerance. 
    test_opts.maxIter = 1000; % maximum iteration number of optimization.
    test_opts.init = 2;  % guess start point from data ZERO.

        [newtask.S, newtask.D] = f_lapMatrix(task.DT{1});
        [newW, epsvalue] = f_MTM_APG(task.DT,task.target,opts,newtask.S, newtask.D);

    trainpl=task.DT{1}*newW(:,1);
    trainpl2=task.DT{2}*newW(:,2);

    trainRMSE(n)=sqrt(mean((trainpl-task.target{1}).^2));
    trainRMSE2(n)=sqrt(mean((trainpl2-task.target{2}).^2));

    aa=task.target{1}-mean(task.target{1});bb=trainpl-mean(trainpl);
    aa2=task.target{2}-mean(task.target{2});bb2=trainpl2-mean(trainpl2);

    trainCO(n)=sum(aa.*bb)/sqrt(sum(aa.^2)*sum(bb.^2));
    trainCO2(n)=sum(aa2.*bb2)/sqrt(sum(aa2.^2)*sum(bb2.^2));



    testpl=test_VBM*newW(:,1);  

     testpl2=test_CC*newW(:,2);

    testRMSE(n)=sqrt(mean((testpl-test_Yt).^2));
    testRMSE2(n)=sqrt(mean((testpl2-test_Yt).^2)); 


    Weight{n}=newW;

    aa=test_Yt-mean(test_Yt);bb=testpl-mean(testpl);
    aa2=test_Yt-mean(test_Yt);bb2=testpl2-mean(testpl2);

    testCO(n)=sum(aa.*bb)/sqrt(sum(aa.^2)*sum(bb.^2));
    testCO2(n)=sum(aa2.*bb2)/sqrt(sum(aa2.^2)*sum(bb2.^2));


    p{n}=testpl;
    p2{n}=testpl2;

end

toc
RMSE_VBM1=[mean( trainRMSE) std( trainRMSE)]
RMSE_CC2=[mean( trainRMSE2) std( trainRMSE2)]


CC_VBM1=[mean( trainCO) std( trainCO)]
CC_CC2=[mean( trainCO2) std( trainCO2)]


RMSE_VBM11=[mean(testRMSE) std(testRMSE)]
RMSE_CC22=[mean(testRMSE2) std(testRMSE2)]


CC_VBM11=[mean(testCO) std(testCO)]
CC_CC22=[mean(testCO2) std(testCO2)]