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

function [old_parcorrf,pdag ,GG, maxScore]=MGraph_BGscore_hillClimb2(d,G,pert_of_change)
%Input: d is the data, G is a directed graph or empty graph, 
%Output: pdag is best graph
%This is the good one 12.15.2003
%this is learning rate + partial corr coef version Feb. 2005
global isordered

t0=clock;
%compute partial correlation
old_var=cov(d);
inv_old_var=pinv(old_var) ;%bug at here sometime the inverse of variance equal inf Dec 22. 2002
parcorrf=-inv_old_var./sqrt(abs(diag(inv_old_var)*diag(inv_old_var)'));

%

%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=20;
alpha0=0.8;
train_len=number_of_repeated;
learn_rate = alpha0 * (0.2/alpha0).^([0:(train_len-1)]/train_len);
%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,pdag);
max_score=current_score
oldest_current_score=current_score;
G=pdag;

for repeated=1:number_of_repeated
    repeated
%check this line
   ischanged=0;
   while max_score<=current_score
       [new_G,current_score]=MGraph_add_best_edge_to_dag(d,G,sigma,alpha,L,T0,TL,isordered);
       if current_score>max_score
           max_score=current_score
           G=new_G;
           ischanged=1;
           disp('Add edge find new pattern');
       else
           break;
       end
   end
   
   current_score=max_score;
   while max_score<=current_score
       [new_G,current_score]=MGraph_remove_worst_edge_in_dag(d,G,sigma,alpha,L,T0,TL,isordered);
       max_score
       if current_score>max_score
           max_score=current_score
           G=new_G;
           ischanged=1;
           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;
    tG=G;
    G

    if repeated <number_of_repeated
            rand('seed',repeated);
            AA=[];BB=[];idx_n=[];parCorf_of_edges=[];temp_edge=[];sorted_parCorf_of_edges=[];sorted_idx=[];
            
            [AA BB]=find(G);
            if (~isempty(AA) & ischanged)
            %if structure changes random remove
                all_edges=[];
                sorted_all_edges=[];
                all_edges=[AA BB]; 
                total_number_of_edges=length(AA);
                idx_n=1:total_number_of_edges;
                
                number_of_changes=ceil(total_number_of_edges*learn_rate(repeated));
                %sort all edges based on their partial corrf
                for i=1:total_number_of_edges
                    temp_edge=all_edges(i,:);
                    parCorf_of_edges(i)=abs(parcorrf(temp_edge(1),temp_edge(2)));
                end
                %sort from minmum to maximum
                [sorted_parCorf_of_edges,sorted_idx]=sort(parCorf_of_edges);
                sorted_all_edges=all_edges(sorted_idx,:);
            elseif (~isempty(AA) & ~ischanged)
             %if structure do not changes make new random input matrix
                all_edges=[];
                sorted_all_edges=[];
                number_of_changes=0;
               % all_edges=wang_kSubset(1:nc,2);
               % total_number_of_edges=length(all_edges);
               % idx_n=randperm(total_number_of_edges);
               % G=zeros(size(G));
               % for ii=1:ceil(length(all_edges)/ceil(rand(1)*10))
               %     G(idx_n(ii))=-1;
               % end
            else
                idx_n=randperm(total_number_of_edges);
                number_of_changes=ceil(total_number_of_edges*learn_rate(repeated));
            end
            %if start with zero G or very sparse graph, the number of changes should 
            %increase otherwise it should enoutgh
            %
            %delete less significant edges
         for jj=1:number_of_changes
                temp_edge=sorted_all_edges(idx_n(jj),:);
                     tempG=G;
%added jbw 10 2004
                     	tempG(temp_edge(1),temp_edge(2))=0;
                     	tempG(temp_edge(2),temp_edge(1))=0;
%added 
			[color,time,dd,ff,phi,back_edge]=dfs(tempG);
                     if isempty(back_edge)
                         G=tempG;
                     end     
            end
            
    end
      
    %removed Feb16.2004
     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
     %add jbw if have loop remive this edge
    [row,col,ss]=find(pdagG);
    tempG=zeros(size(pdagG));
    pdagG=tempG; 
    for i=1:length(row)
       tempG=pdagG;
       tempG(row(i),col(i))=-1;
       [color,time,dd,ff,phi,back_edge]=dfs(tempG);
       if isempty(back_edge)
	  pdagG=tempG;
	end 
    end

     b=gnt_graph_to_coeffiecent_b(pdagG);

     [current_score, T0, TL]=gnt_scoring_completeG(u0,v,b,sigma,alpha,n,L,d,pdagG);
     max_score=current_score;
     G=pdagG;
    %end removed
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)