[old_parcorrf,delet_edge,record_edges]=MGraph_GaussModel_forward(old_var,N,q,cutoff)

function [old_parcorrf,delet_edge,record_edges]=MGraph_GaussModel_forward(old_var,N,q,cutoff)
global graph_form
%P/181, Graphical models in  applied multivariate statistics
%old_var=[ 302.29 125.78 100.43 105.07 116.07; ...
%        125.78 170.88 84.19 93.6 97.89;...
%        100.43 84.19 111.6 110.84 120.49; ...
%        105.07 93.6 110.84 217.88 153.77; ...
%   116.07 97.89 120.49 153.77 294.37]
%N=88 % number of observations
%q=5 % number of variables
%cutoff=0.05;
%old_corrf=old_var./sqrt(diag(old_var)*diag(old_var)');
%inv_old_var=pinv(old_var); %bug at here sometime the inverse of variance equal inf Dec 22. 2002
%old_parcorrf=-inv_old_var./sqrt(abs(diag(inv_old_var)*diag(inv_old_var)'));
%devs=-N*log(1-old_parcorrf.^2) %device of the pair of edges
t0=clock;
%For forward algorithm
%inital model a+b+c+d...
initial_model=eye(q,q);
k_sets=1:q;
full_model=triu(ones(q,q));
[set_a1,set_a2]=find(initial_model);
set_a=[set_a1,set_a2];
ln_set_a=size(set_a,1);

%rest of edges
[set_b1,set_b2]=find(full_model);
set_b=[set_b1,set_b2];
ln_set_b=size(set_b,1);
need_test_sets=setdiff(set_b,set_a,'rows');

%compement set of set a
for i=1:ln_set_a
    comp_set_a{i}=setdiff(k_sets,set_a(i,:)); %complement of set a
end

%make the full model
%oldModel=MGraph_num2string((1:q)');
%temp_newModel=oldModel;


idx=1; %index for recorded edges
temp_G=zeros(q,q);
new_G=temp_G;

while 1
    %need update set_a, comp_set_a, need_test_sets in each loop
    
    fit_var=MGraph_GaussEstimate_variance(old_var,set_a,comp_set_a);
    %initial model device
    det_old=det(old_var);
    det_fit=det(fit_var);
    total_devs=abs(N*(log(det_old)-log(det_fit)));
    df=size(need_test_sets,1);
    total_p=Chisq(total_devs,df);
    

    %computer all pair of edges device
    new_model_set=set_a;
    index=1:df;
    devs=[];
    chtest_p=[];
    
    for i=1:df
        new_model_set(ln_set_a+1,:)=need_test_sets(i,:);
        for j=1:size(new_model_set,1)
            new_comp_set{j}=setdiff(k_sets,new_model_set(j,:));
        end
        new_fit_var=MGraph_GaussEstimate_variance(old_var,new_model_set,new_comp_set);
        det_old=det(fit_var);
        det_fit=det(new_fit_var);
        devs(i)=abs(N*(log(det_old)-log(det_fit)));
        chtest_p(i)=Chisq(devs(i),1) ;   
    end 
    %devs
    chtest_p(find(chtest_p<cutoff))
    %find maximu device and add this edge to initial model
    %[sortDev idx_sortDev]=sort(devs);
    %%flip from top to down
    %flip_sortDev=fliplr(sortDev);
    %flip_idx_sortDev=fliplr(idx_sortDev);
    %idx_maxDev=flip_idx_sortDev(1);
     %[maxDev idx_maxDev]=max(devs);
    %add sep.2003 for check decomposable modle
    %make new subgroup
    %if graph_form==1
    %    for ii=1:length(devs)
    %        idx_maxDev=flip_idx_sortDev(ii);
    %        delet_edge_location=need_test_sets(idx_maxDev,:);
    %   	    delet_edge_string=MGraph_num2string(delet_edge_location);
    %   	    new_model=MGraph_makeSubGraph(temp_newModel,delet_edge_string);
    %   	    clear numberVector
    %   	    for new_idx=1:length(new_model)
    %      		numberVector{new_idx}=num2str(MGraph_string2num(new_model(new_idx)));
    %   	    end
    %    	    input_k_cliques=(numberVector)';
    %        is_SDordering=MGraph_isSDordering(input_k_cliques)
    %        pause
    %        if is_SDordering
    %            break;
    %        end
    %     end
    %else
   %     idx_maxDev=flip_idx_sortDev(1);
   % end       
    %end add    
    
    %add test decomposable 
    sortedDev=[];
    sorted_need_test_sets=[];
    sorted_chtest_p=[];
    
    [sortedDev idx_sortDev]=sort(devs);   
    sorted_need_test_sets=need_test_sets(idx_sortDev,:);
    sorted_chtest_p=chtest_p(idx_sortDev);
    len_of_dev=length(sortedDev);
    for i=1:len_of_dev
         temp_G=new_G;
         temp_add_edge=sorted_need_test_sets(end-i+1,:);
         temp_G(temp_add_edge(1),temp_add_edge(2))=1;
         temp_G(temp_add_edge(2),temp_add_edge(1))=1;
         if graph_form==1
             cond1=(isDecomposableG(temp_G) & sorted_chtest_p(end-i+1)<cutoff);
         else
             cond1=sorted_chtest_p(end-i+1)<cutoff;
         end
         if cond1
             break;
         end
     end
     %end add
    
    if cond1
        %add it to the model
        %add_edge=need_test_sets(idx_maxDev,:);
        add_edge=temp_add_edge;
        new_G(temp_add_edge(1),temp_add_edge(2))=1;
        new_G(temp_add_edge(2),temp_add_edge(1))=1;
        temp_G=new_G;
        %end add in
        
        set_a(end+1,:)=add_edge;
        comp_set_a{end+1}=setdiff(k_sets,set_a(end,:));
        old_need_test_sets=need_test_sets;
        need_test_sets=setdiff(old_need_test_sets,add_edge,'rows');
        ln_set_a=size(set_a,1);
        record_old_var=fit_var;
        delet_edge= size(need_test_sets,1)+1;
        record_edges{idx}=num2str(add_edge);
        record_edges'
        idx=idx+1;
        %temp_newModel=new_model  ;
    else
        record_old_var=fit_var;
        delet_edge= size(need_test_sets,1);

        break;
    end
end %end while

old_corrf=record_old_var./sqrt(diag(record_old_var)*diag(record_old_var)');
inv_old_var=pinv(record_old_var); %bug at here sometime the inverse of variance equal inf Dec 22. 2002
old_parcorrf=-inv_old_var./sqrt(abs(diag(inv_old_var)*diag(inv_old_var)'));
etime(clock,t0)