leaveOneOut_testGnet_pkc.m

close all
clear all
global isordered graph_form

file1='/u3/junbai/local/gnet_toolbox/demo_data/gaussInput_PKC_pathway_MissImput.dat'
[old_data,varnames,casenames] = tblread(file1,'\t') ;

file2='/u3/junbai/local/gnet_toolbox/demo_data/real_pkc_pathway.txt'
[net_data,net_varnames,net_casenames] = tblread(file2,'\t') ;
fid=fopen(['trick_for_PKC'],'w');

for loop=1:10
clear spareM old_parcorrf
%loop 10 times

data=old_data;
row_no=randperm(size(data,1));
data(row_no(1),:)=randn(1,size(data,2));

outdata.labels=casenames;
outdata.data=data;
outdata.N=size(outdata.data,2);
outdata.q=size(outdata.data,1);
cutoff_P=0.05;
outdata.cutoff_P=cutoff_P;
tdata=som_normalize(outdata.data','var');
old_var=cov(tdata);
N=outdata.N; % number of observations
q=outdata.q; % number of variables
cutoff=outdata.cutoff_P;
algormChoice=7
graph_form=1;
isordered =2

    if algormChoice==3
        pdag_wang = wang_learn_struct_pdag_pc('cond_indep_fisher_z', q, 5, old_var, N, cutoff);
    elseif algormChoice==6
        %IC forward fast
        marray_debuge('Run IG model with Forward fast algorizm')
        [old_parcorrf,delet_edge,record_edge,pdag_wang,G]=MGraph_IG_forward_fast2(old_var,N,q,cutoff);
    elseif algormChoice==7
        %IG with Bscore and IG % this one was not stable
        marray_debuge('Run IG plus BScore search algorizm')
         [old_parcorrf,delet_edge,record_edge,pdag_wang,G]=MGraph_IG_forwardFast_BGscore(tdata,old_var,N,q,cutoff);
     elseif algormChoice==9
        marray_debuge('Run BScore search algorizm with learning rate plus partial corr. coef.')
        tG= setdiag(triu(ones(q,q)),0).*(-1); %setdiag(triu(rand(q,q)>0.5).*(-1),0); %setdiag(triu(ones(q,q)),0).*(-1);
        [old_parcorrf,pdag_wang ,GG, maxScore]=MGraph_BGscore_hillClimb2(tdata,tG,0.15);
     end

   if algormChoice<3
        marray_debuge(['Estimated partial correlation coefficient= ']);
        %replace small valus as 0
        tempM=(abs(old_parcorrf)>0.0000000001);
        temp_old_parcorrf=tempM.*old_parcorrf;
        if exist('vpa')~=0
            marray_debuge([num2str(double((temp_old_parcorrf)))]);
        end
        marray_debuge(['Number of deleted edge= ', num2str(delet_edge)]);
        %marray_debuge(['Deleted edges= ']);
        %marray_debuge([ num2str(record_edge)]);
    end
    
    if algormChoice<3
        %draw figures
        spareM=(abs(old_parcorrf)>0.0000000001 & abs(old_parcorrf)~=1);
        %spareM=sparse(spareM);
        labels=strrep(cellstr( deblank(outdata.labels) ),'"','');
       % draw_graph(spareM,lower(labels));
    elseif (algormChoice==3 | algormChoice==4 | algormChoice==5 | algormChoice==6 | algormChoice==7| algormChoice==8| algormChoice==9)
        spareM=pdag_wang
        labels=strrep(cellstr( deblank(outdata.labels) ),'"','');
       % draw_graph(abs(spareM),lower(labels));
    end


 %start to compute the error rate
    %compare net_Data and abs(spareM)
    if (exist('file2')& file2~=0)
    AA=[]; BB=[];ispath_C=[];
    rel_A=[];rel_B=[];
    edge_existence_error_c=0;
    predicted_net=abs(spareM);
    [AA BB]=find(predicted_net);
    [rel_A rel_B]=find(net_data);
    len_rel_A=length(rel_A);
    %edge existence error of commission
    len_AA=length(AA);
    
    ispath_C = reachability_graph(abs(net_data));
    for ii=1:len_AA
        temp_edge=[AA(ii), BB(ii)];
        if ( net_data(temp_edge(1),temp_edge(2))==0 & (ispath_C(AA(ii),BB(ii))==0 & ispath_C(BB(ii),AA(ii))==0 ) )
            edge_existence_error_c=edge_existence_error_c+1;   
        end
    end
    ratio_of_edge_existence_error_c(loop)=edge_existence_error_c/len_rel_A
    
    %edge direction error of commission
    len_AA=length(AA);
    edge_direction_error_c=0;
    for ii=1:len_AA
        temp_edge=[AA(ii), BB(ii)];
        if ( abs(net_data(temp_edge(1),temp_edge(2))) + abs(net_data(temp_edge(2),temp_edge(1))) )>0
            if abs(net_data(temp_edge(1),temp_edge(2)))~=1
                edge_direction_error_c=edge_direction_error_c+1;
            end
        end
    end
    ratio_of_edge_direction_error_c(loop)=edge_direction_error_c/len_rel_A
    
  %edge existence error of ommission
    AA=[]; BB=[]; ispath_P=[];
    edge_existence_error_o=0;
    ispath_P = reachability_graph(abs(predicted_net));
    [AA BB]=find(setdiag(predicted_net,1)==0);
    len_AA=length(AA);
     for ii=1:len_AA
         temp_edge=[AA(ii),BB(ii)];
         %net_data(temp_edge)
         %pause
         if ( net_data(temp_edge(1),temp_edge(2))==-1 & (ispath_P(AA(ii),BB(ii))==0 & ispath_P(BB(ii),AA(ii))==0) )
              edge_existence_error_o= edge_existence_error_o+1;
         end
     end
     ratio_of_edge_existence_error_o(loop)=edge_existence_error_o/len_rel_A
 end      
 fprintf(fid,'%s %s\n','End loop: ', num2str(loop)); 
end
%end loop
fclose(fid);

mean_c=mean(ratio_of_edge_existence_error_c )
mean_dc=mean(ratio_of_edge_direction_error_c)
mean_o=mean(ratio_of_edge_existence_error_o)


save test_GNET_IG_GN_pkc2.mat
exit;