MGraph_BGscore(d,G)

function MGraph_BGscore(d,G)
%Input: d is the data, G is the initial graph
%Output: best graph

[nr nc]=size(d);
u0=median(d);
n=nc;
L=nr;
sigma=nc+2;
alpha=nc+2;
v=ones(1,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);
max_score=current_score;
done=0;

while max_score<=current_score
        %add path to graph
        A=[];
        B=[];
        temp_edges={};
        idx=1;
        pscore=[];
        done=0;
        
        [A B]=find(abs(setdiag(G,1))==0); %potential add edges    
        if isempty(A) break; end
        
        for i=1:length(A);
            tempG=G;
            tempG(A(i),B(i))=1;
            if ~MGraph_isundirected(tempG); 
               %start add an undirected edge
               temp_edges{idx}=[A(i) B(i)];
               pscore(idx) = log10(gnt_scoring_uncompleteG(tempG,sigma,alpha,L,T0,TL));
               idx=idx+1;     
               %start to add an unshielded colliders at a if possibile ??     
               %start to add an unshielded colliders at b if possibile ??
            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
            done=1;
        else
            done=1;
        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
        elseif current_score==max_score & done
            break;
        end
end


current_score=max_score;
done=0;
while max_score<=current_score
    rA=[];
    rB=[];
    temp_r_edges={};
    rpscore=[];
    rev_score=[];
    temp_rev_edges={};
    both_p=[];
    idx1=1;
    idx2=1;
    
    need_r_edge=[];
    done=0;
    
    %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;
            %if ~MGraph_isundirected(temp_rG); 
                temp_r_edges{i}=[rA(i) rB(i)];
                rpscore(i)=log10(gnt_scoring_uncompleteG(temp_rG,sigma,alpha,L,T0,TL));
                %idx1=idx1+1;
                %end
            
            %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(idx2)=log10(gnt_scoring_uncompleteG(temp_rev_G1,sigma,alpha,L,T0,TL));
            %    idx2=idx2+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(idx2)=log10(gnt_scoring_uncompleteG(temp_rev_G2,sigma,alpha,L,T0,TL));
            %    idx2=idx2+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));
             %   idx2=idx2+1;
            %end       
      end
    
      
      if ~isempty(temp_r_edges)
            [best_rpscore best_rp]=max(rpscore);
            best_r_edge=temp_r_edges{best_rp};
            current_score=best_rpscore;
            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
            done=1;
        else
            done=1;
        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
        
      
    
       % [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,:);
       %     if ~isempty(need_r_edge)
       %         current_G(need_r_edge(1),need_r_edge(2))=-1;
       %         current_G(need_r_edge(2),need_r_edge(1))=0;
       %         done=1;
       %     else
      %          done=1;
      %      end
      %  else
      %  %remove
      %      current_G=G;
      %      need_r_edge=both_p(best_p,:);
      %      if ~isempty(need_r_edge)
      %          current_G(need_r_edge(1),need_r_edge(2))=0;
      %          current_G(need_r_edge(2),need_r_edge(1))=0;
      %          done=1;    
      %      else
     %           done=1;
     %       end
     %   end
    
        % 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
    elseif current_score==max_score & done
        break;
    end
end