function [old_parcorrf,pdag ,GG, maxScore]=MGraph_BGscore_average(d,G)
[nr nc]=size(d);
u0=median(d);
n=nc;
L=nr;
sigma=nc+2;
alpha=nc+2;
v=ones(1,n);
%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;
AA={};
gg={G}; %start model
gg_pscore(1)=current_score;
ln_of_gg=length(gg);
OL=-log10(5);
OR=5;
AA_pscore=[];
%up algorithm
while ~isempty(gg)
gg
gg_pscore
M=gg{1} %step 1
AA,AA_pscore
AA_pscore(end+1)=gg_pscore(1)
pause
current_score=gg_pscore(1);
AA(end+1)={M}; %accepted models; step 2
t_gg=gg(2:end);
gg={}; %consider models
gg=t_gg;
t_gg={};
t_gg_pscore=gg_pscore(2:end);
gg_pscore=[];
gg_pscore=t_gg_pscore;
t_gg_pscore=[];
A=[];B=[];
[A B]=find(abs(triu(setdiag(M,1)))==0); %potential add edges
idx=1;
tempG=[];
temp_graphs={};
for i=1:length(A); %step 3
%add edge in both direction
x=A(i);y=B(i);
tempG=M;
if (tempG(x,y)==0 & tempG(y,x)==0)
tempG(x,y)=-1;
tempG(y,x)=0;
[color,time,d,f,phi,back_edge]=dfs(tempG);
cond1= isempty(back_edge);
if cond1
current_score2=(gnt_scoring_uncompleteG(tempG,sigma,alpha,L,T0,TL));
%current_score, current_score2
BB=current_score-current_score2;
%OL,OR
if BB<OL %step 5
t_AA={};
t_AA_pscore=[];
t_AA_idx=1;
for j=1:length(AA)
if sum(sum(abs(M)-abs(AA{j})))~=0
t_AA{t_AA_idx}=AA{j};
t_AA_pscore(t_AA_idx)=AA_pscore(j);
t_AA_idx=t_AA_idx+1;
end
end
AA={};
AA_pscore=[];
AA=t_AA;
AA_pscore=t_AA_pscore;
%end if
isfound_M=0;
for j=1:length(gg)
if sum(sum(abs(tempG)-abs(gg{j})))==0
isfound_M=1;
break;
end
end
if (~isfound_M & current_score2>current_score)
gg(end+1)={tempG};
gg_pscore(end+1)=current_score2;
end
end %end step 5
if (BB<=OR & BB>=OL) %step 6
isfound_M=0;
for j=1:length(gg)
if sum(sum(abs(tempG)-abs(gg{j})))==0
isfound_M=1;
break;
end
end
% current_score2,current_score
if (~isfound_M & current_score2>current_score)
gg(end+1)={tempG}
gg_pscore(end+1)=current_score2
end
end %end step 6
end % end cond1
%add reversed edge
tempG(x,y)=0;
tempG(y,x)=-1;
[color,time,d,f,phi,back_edge]=dfs(tempG);
cond1= isempty(back_edge);
if cond1
current_score2=(gnt_scoring_uncompleteG(tempG,sigma,alpha,L,T0,TL));
BB=current_score-current_score2;
%OL,OR
if BB<OL %step 5
t_AA={};
t_AA_pscore=[];
t_AA_idx=1;
for j=1:length(AA)
if sum(sum(abs(M)-abs(AA{j})))~=0
t_AA{t_AA_idx}=AA{j};
t_AA_pscore(t_AA_idx)=AA_pscore(j);
t_AA_idx=t_AA_idx+1;
end
end
AA={};
AA_pscore=[];
AA=t_AA;
AA_pscore=t_AA_pscore;
%end if
isfound_M=0;
for j=1:length(gg)
if sum(sum(abs(tempG)-abs(gg{j})))==0
isfound_M=1;
break;
end
end
if (~isfound_M & current_score2>current_score)
gg(end+1)={tempG}
gg_pscore(end+1)=current_score2
end
end %end step 5
if (BB<=OR & BB>=OL) %step 6
isfound_M=0;
for j=1:length(gg)
if sum(sum(abs(tempG)-abs(gg{j})))==0
isfound_M=1;
break;
end
end
%current_score2,current_score
if (~isfound_M & current_score2>current_score)
gg(end+1)={tempG}
gg_pscore(end+1)=current_score2
end
end %end step 6
end % end cond1
end %end tempG
% i
%pause
end %end step 3
end %while
%down algorithm
%out put
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)'));
AA_pscore
AA{1}
[maxScore maxScoreI]=max(AA_pscore)
GG=AA{maxScoreI}
pdag=GG;
