test_baysein_network.m

clear all
d=[ -0.78 -1.55 0.11;...
        0.18 -3.04 -2.35;...
        1.87 1.04 0.48; ...
        -0.42 0.27 -0.68; ...
        1.23 1.52 0.31; ...
        0.51 -0.22 -0.6; ... 
        0.44 -0.18 0.13; ...
        0.57 -1.82 -2.76; ...
        0.64 0.47 0.74 ; ...
        1.05 0.15 0.2 ; ...
        0.43 2.13 0.63; ...
        0.16 -0.94 -1.96; ...
        1.64 1.25 1.03 ; ...
        -0.52 -2.18 -2.31; ...
        -0.37 -1.3 -0.7; ...
        1.35 0.87 0.23 ; ...
        1.44 -0.83 -1.61; ...
        -0.55 -1.33 -1.67; ...
        0.79 -0.62 -2.00 ; ...
        0.53 -0.93 -2.92]
[nr nc]=size(d);
u0=median(d)
n=nc
L=nr
sigma=nc+2
alpha=nc+2
v=ones(1,n)
%initial structure, 1,2,3
%initial structure, 1->2->3
temp_G=zeros(n,n);
G=[0 0 1; 0 0 1; 0 0 0]

b=gnt_graph_to_coeffiecent_b(G)

[current_score, T0, TL]=gnt_scoring_completeG(u0,v,b,sigma,alpha,n,L,d);
current_score=log10(current_score); %set a initial small value ??(current_score);
%[A B]=find(setdiag(G,1)==0);
%for i=1:length(A)
%    if A(i)<B(i)
%        G(A(i),B(i))=1;
%        break;
%    end
%end
%current_score = log10(gnt_scoring_uncompleteG(G,sigma,alpha,L,T0,TL));
max_score=current_score


while max_score<=current_score
        %add path to graph
        
        %dont test undirected graph !!
        %
        A=[];
        B=[];
        temp_edges={};   
        
        [A B]=find(setdiag(G,1)==0); %potential add edges    
        if isempty(A) break; end
     
        idx=1;
        pscore=[];
        for i=1:length(A);
            tempG=G;
            %if (tempG(B(i),A(i))==0 & A(i)<B(i))
                %force the order
                tempG(A(i),B(i))=1;
                temp_edges{idx}=[A(i) B(i)];
                pscore(idx) = log10(gnt_scoring_uncompleteG(tempG,sigma,alpha,L,T0,TL));
                idx=idx+1;
                %end
        end
        
        if ~isempty(temp_edges)
            [best_pscore, best_p] = max(pscore);
            best_edge = temp_edges{best_p};
            current_score=best_pscore;
            current_G=G;
            if current_G(best_edge(2),best_edge(1))==0
                current_G(best_edge(1),best_edge(2))=-1;
            else
                current_G(best_edge(1),best_edge(2))=1;
                current_G(best_edge(2),best_edge(1))=1;
            end
        end
        %[pdag, tG] = wang_learn_struct_undirect2pdag(1,'cond_indep_fisher_z', n, 5,current_G,cov(d), L, 0.05);
        %current_score=(gnt_scoring_uncompleteG(abs(pdag),sigma,alpha,L,T0,TL));
        if current_score>max_score
            max_score=current_score;
            G=current_G;
            best_edge
        else
            break;
        end
end

current_score=max_score;


while max_score<=current_score
    rA=[];
    rB=[];
    temp_r_edges={};
    rpscore=[];
    rev_score=[];
    temp_rev_edges={};
    idx=1;
     
    isundirected=MGraph_isundirected(G);
    if ~isundirected
    %star test
        [rA rB]=find(abs(G)==1); %potential remove edge/ another is revers edge ??
    
        for i=1:length(rA)
            temp_rG=G;
            %remove edge
            temp_rG(rA(i),rB(i))=0;
            temp_r_edges{i}=[rA(i) rB(i)];
            rpscore(i)=log10(gnt_scoring_uncompleteG(temp_rG,sigma,alpha,L,T0,TL));
        
            %reverse of edge
            temp_rev_G=G;
            temp_rev_edge=[rA(i),rB(i)];
            if ( temp_rev_G(temp_rev_edge(1),temp_rev_edge(2))==1 & temp_rev_G(temp_rev_edge(2),temp_rev_edge(1))==1)
                %try to reverse the edge 2 way
                temp_rev_edges{idx}=[temp_rev_edge(1),temp_rev_edge(2)];
                temp_rev_G1=temp_rev_G;
                temp_rev_G1(temp_rev_edge(1),temp_rev_edge(2))=-1;
                temp_rev_G1(temp_rev_edge(2),temp_rev_edge(1))=0;
                rev_score(idx)=log10(gnt_scoring_uncompleteG(temp_rev_G1,sigma,alpha,L,T0,TL));
        
                idx=idx+1;
                temp_rev_edges{idx}=[temp_rev_edge(2),temp_rev_edge(1)];
                temp_rev_G2=temp_rev_G;
                temp_rev_G2(temp_rev_edge(2),temp_rev_edge(1))=-1;
                temp_rev_G2(temp_rev_edge(1),temp_rev_edge(2))=0;
                rev_score(idx)=log10(gnt_scoring_uncompleteG(temp_rev_G2,sigma,alpha,L,T0,TL));
                idx=idx+1;
            else
                %reverse one way
                temp_rev_edges{idx}=[temp_rev_edge(2),temp_rev_edge(1)];
                temp_rev_G2=temp_rev_G;
                temp_rev_G2(temp_rev_edge(2),temp_rev_edge(1))=-1;
                temp_rev_G2(temp_rev_edge(1),temp_rev_edge(2))=0;
                rev_score(idx)=log10(gnt_scoring_uncompleteG(temp_rev_G2,sigma,alpha,L,T0,TL));
                idx=idx+1;
            end       
        end
      
    
        [best_rpscore best_rp]=max(rpscore);
        best_r_edge=temp_r_edges{best_rp};
    
        [best_revscore best_revp]=max(rev_score);
        best_rev_edge=temp_rev_edges{best_revp};
    
        both_score=[best_revscore,best_rpscore];
        both_p=[best_rev_edge;best_r_edge];
    
        [best_score best_p]=max(both_score);
        current_score=best_score; 
        if best_p==1
        %reverse
            current_G=G;
            need_r_edge=both_p(best_p,:);
            current_G(need_r_edge(1),need_r_edge(2))=-1;
            current_G(need_r_edge(2),need_r_edge(1))=0;
        else
        %remove
            current_G=G;
            need_r_edge=both_p(best_p,:);
            current_G(need_r_edge(1),need_r_edge(2))=0;
            current_G(need_r_edge(2),need_r_edge(1))=0;
        end
    
        %current_G=G;
        %if current_G(best_r_edge(2),best_r_edge(1))==1
        %    current_G(best_r_edge(1),best_r_edge(2))=0;
        %    current_G(best_r_edge(2),best_r_edge(1))=-1;
        %else
       %     current_G(best_r_edge(1),best_r_edge(2))=0;
       % end
        
        %[pdag, tG] = wang_learn_struct_undirect2pdag(1,'cond_indep_fisher_z', n, 5,current_G,cov(d), L, 0.05);
        %current_score=(gnt_scoring_uncompleteG(abs(pdag),sigma,alpha,L,T0,TL));
        %current_score=best_rpscore
    else
        %%change undirected graph to pattern
        [chA chB]=find(abs(G)==1);
        G(chA(end),chB(end))=0;
    end
    
    if current_score>max_score;
        max_score=current_score;
        G=current_G;
        %need_r_edge
        best_r_edge
    end
end