[old_parcorrf,pdag ,GG, maxScore]=MGraph_BGscore_greedyLocal(d,G,pert_of_change)

function [old_parcorrf,pdag ,GG, maxScore]=MGraph_BGscore_greedyLocal(d,G,pert_of_change)
%Input: d is the data, G is the initial pattern graph, it should use pattern at here
%Output: pdag is best graph

t0=clock;
[nr nc]=size(d);
u0=median(d);
n=nc;
L=nr;
sigma=nc+2;
alpha=nc+2;
v=ones(1,n);

%test oct 12
all_edges=wang_kSubset(1:nc,2); 
total_number_of_edges=(nc^2-nc)/2;
number_of_repeated=10;
   
[pdag, G_patern] = wang_learn_struct_undirect2pdag(1,'cond_indep_fisher_z', nc, 2,G, cov(d), L, 0.05);
%check whether have edges go from high order to lower order and remove it
[pa pb]=find(pdag);
isreverse=pa>pb;
idx_remove=find(isreverse==1);
for i=1:length(idx_remove)
        pdag(pa(idx_remove(i)),pb(idx_remove(i))) =0;
        pdag(pb(idx_remove(i)),pa(idx_remove(i))) =-1;
end

b=gnt_graph_to_coeffiecent_b(pdag);
[current_score, T0, TL]=gnt_scoring_completeG(u0,v,b,sigma,alpha,n,L,d);
max_score=current_score

%test here
test_G=G;
for repeated=1:number_of_repeated
    repeated
    pdag=[];
    GG=[];
    pdag=[];G_patern=[];pa=[];pb=[];
    
     
    
    
%check this line
   
   while max_score<=current_score
       new_G=[];
       pdagG=[];
       GG=[];
       new_G=MGraph_add_best_edge_to_pat(d,test_G,sigma,alpha,L,T0,TL);
       [pdagG, GG] = wang_learn_struct_undirect2pdag(1,'cond_indep_fisher_z', nc, 2,new_G, cov(d), L, 0.05);
       current_score  = (gnt_scoring_uncompleteG(pdagG,sigma,alpha,L,T0,TL));
       if current_score>max_score
           max_score=current_score
           G=new_G;
           test_G=G;
           disp('Add edge find new pattern');
       else
           break;
       end
   end
   
   current_score=max_score;
   while max_score<=current_score
       new_G=[];
       pdagG=[];
       GG=[];
       new_G=MGraph_remove_worst_edge_in_pat(d,test_G,sigma,alpha,L,T0,TL);
       [pdagG, GG] = wang_learn_struct_undirect2pdag(1,'cond_indep_fisher_z', nc, 2,new_G, cov(d), L, 0.05);
       current_score  = (gnt_scoring_uncompleteG(pdagG,sigma,alpha,L,T0,TL));
       if current_score>max_score
           max_score=current_score
           G=new_G;
           test_G=G;
           disp('Remove edge find new pattern');
       else
           break;

       end
   end

   %end while of repeate 1     
    %repeated again
    if repeated>1
        if max_score>all_maxScore(repeated-1)
            all_maxScore(repeated)=max_score
            all_G{repeated}=G;
        else
            all_maxScore(repeated)=all_maxScore(repeated-1)
            all_G{repeated}=all_G{repeated-1};
            G=all_G{repeated-1};
        end
    else
         all_maxScore(repeated)=max_score
         all_G{repeated}=G;
    end
    
    G
    test_G=G;
    %[pdagG, GG] = wang_learn_struct_undirect2pdag(1,'cond_indep_fisher_z', nc, 2,G, cov(d), L, 0.05);
    if repeated <number_of_repeated
            %max_score=current_score;
            rand('seed',repeated);
            AA=[];BB=[];idx_n=[];
            [AA BB]=find(triu(test_G));
            if ~isempty(AA)
                all_edges=[AA BB]; 
                total_number_of_edges=length(AA);
                idx_n=randperm(total_number_of_edges);
                number_of_changes=ceil(total_number_of_edges*pert_of_change);
            else
                idx_n=randperm(total_number_of_edges);
                number_of_changes=ceil(total_number_of_edges*(1-pert_of_change));
            end
            %if start with zero G or very sparse graph, the number of changes should 
            %increase otherwise it should enoutgh
            %
            %all_edges=[AA BB];
            for jj=1:number_of_changes
                temp_edge=all_edges(idx_n(jj),:);
                if mod(idx_n(jj),2)
                   if (test_G(temp_edge(1),temp_edge(2))==1 & test_G(temp_edge(2),temp_edge(1))==1)
                       test_G(temp_edge(1),temp_edge(2))=-1;
                       test_G(temp_edge(2),temp_edge(1))=0;
                   else
                       test_G(temp_edge(1),temp_edge(2))=1;
                       test_G(temp_edge(2),temp_edge(1))=1;
                   end
               else
                   test_G(temp_edge(1),temp_edge(2))=0;
                   test_G(temp_edge(2),temp_edge(1))=0;
               end
                    %if ( G(temp_edge(1),temp_edge(2))==0 )
                   %     G(temp_edge(1),temp_edge(2))=1;
                   %     G(temp_edge(2),temp_edge(1))=1;
                   %elseif abs(G(temp_edge(1),temp_edge(2)))==1
                   %      G(temp_edge(1),temp_edge(2))=0;
                   %     G(temp_edge(2),temp_edge(1))=0;
                   % end
            end
            
    end
    %current_score=max_score;
     [pdagG, GG] = wang_learn_struct_undirect2pdag(1,'cond_indep_fisher_z', nc, 2,test_G, cov(d), L, 0.05);
     b={};
     isreverese=[];
     idx_remove=[];
     [pa pb]=find(pdagG);
     isreverse=pa>pb;
     idx_remove=find(isreverse==1);
     for i=1:length(idx_remove)
        pdagG(pa(idx_remove(i)),pb(idx_remove(i))) =0;
        pdagG(pb(idx_remove(i)),pa(idx_remove(i))) =-1;
    end
     b=gnt_graph_to_coeffiecent_b(pdagG);
    [current_score, T0, TL]=gnt_scoring_completeG(u0,v,b,sigma,alpha,n,L,d);
     max_score=current_score;
  
   
    
end %end all
%end test

GG=[];
[maxScore idx_maxScore]=max(all_maxScore)
GG=all_G{idx_maxScore}
[pdag, tG] = wang_learn_struct_undirect2pdag(1,'cond_indep_fisher_z', nc, 5,GG,cov(d), L, 0.05);

record_old_var=cov(d);
old_corrf=record_old_var./sqrt(diag(record_old_var)*diag(record_old_var)');
inv_old_var=pinv(record_old_var);
old_parcorrf=-inv_old_var./sqrt(abs(diag(inv_old_var)*diag(inv_old_var)'));


etime(clock,t0)