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

function [old_parcorrf,pdag ,GG, maxScore]=MGraph_BGscore_hillClimb(d,G,pert_of_change)
%Input: d is the data, G is a directed graph or empty graph, 
%Output: pdag is best graph
global isordered

t0=clock;

%isordered =0 means input variables are not ordered, if isordered=1 then the input variables are ordered according the network
[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=100;
   
%check whether have edges go from high order to lower order and remove it
pdag=G;
[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
G=pdag;

for repeated=1:number_of_repeated
    repeated
%check this line
   
   while max_score<=current_score
       new_G=MGraph_add_best_edge_to_dag(d,G,sigma,alpha,L,T0,TL,isordered);
       current_score  = (gnt_scoring_uncompleteG(new_G,sigma,alpha,L,T0,TL));
       if current_score>max_score
           max_score=current_score
           G=new_G;
           disp('Add edge find new pattern');
       else
           break;
       end
   end
   
   current_score=max_score;
   while max_score<=current_score
       new_G=MGraph_remove_worst_edge_in_dag(d,G,sigma,alpha,L,T0,TL,isordered);
       current_score  = (gnt_scoring_uncompleteG(new_G,sigma,alpha,L,T0,TL));
       if current_score>max_score
           max_score=current_score
           G=new_G;
           disp('Remove edge find new pattern');
       else
           break;

       end
   end

   %end while of repeate 1     
    %repeated again
    all_maxScore(repeated)=max_score
    all_G{repeated}=G;
    G

    if repeated <number_of_repeated
            rand('seed',repeated);
            AA=[];BB=[];idx_n=[];
            [AA BB]=find(G);
            %pert_of_change=rand(1);
            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
            %
            for jj=1:number_of_changes
                temp_edge=all_edges(idx_n(jj),:);
                if ((abs(G(temp_edge(1),temp_edge(2)))==1 & G(temp_edge(2),temp_edge(1))==0) & mod(idx_n(jj),2)) 
                    tempG=G;
                    tempG(temp_edge(1),temp_edge(2))=0;
                    tempG(temp_edge(2),temp_edge(1))=-1;
                    [color,time,dd,ff,phi,back_edge]=dfs(tempG);
                    if isempty(back_edge)
                        G=tempG;
                    end
                elseif ((abs(G(temp_edge(1),temp_edge(2)))==0 & abs(G(temp_edge(2),temp_edge(1)))==1) & mod(idx_n(jj),2))
                     tempG=G;
                     tempG(temp_edge(1),temp_edge(2))=-1;
                     tempG(temp_edge(2),temp_edge(1))=0;
                     [color,time,dd,ff,phi,back_edge]=dfs(tempG);
                     if isempty(back_edge)
                         G=tempG;
                     end
                elseif ~mod(idx_n(jj),2) 
                     tempG=G;
                     tempG(temp_edge(1),temp_edge(2))=0;
                     tempG(temp_edge(2),temp_edge(1))=0;
                     [color,time,dd,ff,phi,back_edge]=dfs(tempG);
                     if isempty(back_edge)
                         G=tempG;
                     end     
                else
                    tempG=G;
                    tempG(temp_edge(1),temp_edge(2))=-1;
                    tempG(temp_edge(2),temp_edge(1))=0;
                    [color,time,dd,ff,phi,back_edge]=dfs(tempG);
                    if isempty(back_edge)
                         G=tempG;
                    end
                end
                %if (abs(G(temp_edge(1),temp_edge(2)))+abs(G(temp_edge(2),temp_edge(1)))>0)
                %       G(temp_edge(1),temp_edge(2))=0;
                %       G(temp_edge(2),temp_edge(1))=0;
                %end
            end
            
    end
            
     pdagG=G;
     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;
     G=pdagG;
     
   
    
end %end all
%end test

GG=[];
[maxScore idx_maxScore]=max(all_maxScore)
GG=all_G{idx_maxScore}
pdag=GG;
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)