function selected_model=MGraph_loglineModel(testdata,numoftype_ineachGp,orderof_Gp,columnStart,cutoff_P)
%For current version there is something wrong
%1) if the import data have zero entries.
%2) some bug may appeared in device or df calculations which confuse me the last 2 label were always wrong
%3) it only works for undicrect the graph.
%4) The device calculation have to be improved, currently each edge only appear in one subgroup will
%be considered to remove, so how to deal with ,if one edge appear in severl group and how to
%compute this device and degree freedom ?? 12.05.2002
%5) Therefore, it is fine for decomposiable model but have one potentional bug for unrestrited model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
%bellow is the main loop
%%Author junbai Wang
temp_model{1}=[]; %added sep 2006
temp_model{2}=[];
lnof_Gp=length(orderof_Gp);
for i=1:lnof_Gp
labelof_Gp(i)=char(64+orderof_Gp(i));
end
%
out=MGraph_MainGausslogline(testdata,numoftype_ineachGp,orderof_Gp,columnStart);
out.cutoff_P=cutoff_P;
out.labelof_Gp=labelof_Gp;
%added 12.05
out.new_model=labelof_Gp;
%end add
%I need add decomposable test in output too 12.06
[removeEG,needtestEG]=MGraph_loglineOutput(out) ;%calculate P value and select which edge going to remove
OldneedtestEG=needtestEG; %define the edge P value bigger than the cutoff P
%if in coherance step, the edge did not exist at this oldneedtestEG will not be removed in
%the later test, even if it's P are bigg than cutoff_P in the later model.
%if in noncoherence step, all the edge bigger than the cut off P will be removed, whenever it is found.
%
%needtestEG are edges need further test
%removeEG
if ~isempty(removeEG.name)
marray_debuge(['Remove Edge: ', strvcat(removeEG.name)]);
marray_debuge(['Degree freedom= ',num2str(removeEG.df)]);
marray_debuge(['Device= ',num2str(removeEG.device)]);
marray_debuge(['P= ',num2str(removeEG.P)]);
marray_debuge(['Location: ', num2str(removeEG.location)]);
marray_debuge(['Selected model: ']);
marray_debuge(strvcat(removeEG.new_model));
end
%marray_debuge(['Model: ', strvcat(removeEG.model)]);
marray_debuge(' ');
if ~isempty(removeEG.name)
total.device=removeEG.device;
total.df=removeEG.df;
total.P=Chisq(total.device,total.df);
marray_debuge(['Total device= ', num2str(total.device)]);
marray_debuge(['Total degree freedom= ', num2str(total.df)]);
marray_debuge(['Total P= ', num2str(total.P)]);
else
disp(['Model selection finished ']);
%disp(removeEG.model);
marray_debuge('Selected model=');
marray_debuge(removeEG.model);
selected_model={removeEG.model};
return;
end
if length(removeEG.model) > 2
while 1
%delete removed edge from needtestEG
[needtestEG.allname ti]=setdiff(needtestEG.allname,removeEG.name);
needtestEG.new_model=removeEG.new_model;
%if empty in needtestEG then stop
if isempty(needtestEG.allname)
disp(['Need test edge list is empty']);
disp(['Cut off P= ',num2str(cutoff_P)]);
disp(['Model selection finished, selected model: ']);
disp(removeEG.new_model);
%total
marray_debuge(['Total device= ', num2str(total.device)]);
marray_debuge(['Total degree freedom= ', num2str(total.df)]);
marray_debuge(['Total P= ', num2str(total.P)]);
marray_debuge('Selected model= ');
marray_debuge(strvcat(removeEG.new_model));
selected_model=removeEG.new_model;
break;
end %end if
lnof_newmodel=length(removeEG.new_model);
clear new_result new_out
for i=1:lnof_newmodel
new_out{i}=MGraph_loglineMergeEdge(i,removeEG,testdata,out); %make new group of data
new_result{i}=MGraph_MainGausslogline(new_out{i}.testdata,new_out{i}.numoftype_ineachGp,...
new_out{i}.orderof_Gp,new_out{i}.startnumof_Col); %calculate GUSS based on each new group
new_result{i}.labelof_Gp=new_out{i}.labelof_Gp;
new_result{i}.cutoff_P=cutoff_P;
end
%test P>0.05
clear test_edge;
loop=1;
loop2=1;
for i=1:lnof_newmodel
%added 12.05
if ~isempty(removeEG.new_model)
for il=1:length(removeEG.new_model) %added sep 2006
if ~isempty(removeEG.new_model{il})
newremoveEG.new_model{loop}=removeEG.new_model{il};
loop=loop+1;
%else
% newremoveEG.new_model{il}='';
end
end
new_result{i}.new_model=newremoveEG.new_model;
%end add
tp_test_edge=MGraph_loglineOutput2(new_result{i},OldneedtestEG) ;%1 as P>0.05
% if (~isempty(tp_test_edge.edge_name) & ~isempty(tp_test_edge.edge_df) )
% tp_test_edge
test_edge{i}=tp_test_edge;
% loop2=loop2+1;
%end
end
end
idx=1; %test the edge only exist in one group
clear temp_commedge temp_commedge_idx temp_commedge_jdx
for i=1:lnof_newmodel
ln_edge=length(test_edge{i}.edge_df);
test_edge{i}.edge_duplicat=zeros(1,ln_edge);
end
for i=1:lnof_newmodel-1
for j=i+1:lnof_newmodel
[temp_commedge{idx} temp_commedge_idx{idx} temp_commedge_jdx{idx} ]=...
intersect(test_edge{i}.edge_name,test_edge{j}.edge_name);
if ~isempty(temp_commedge_idx{idx})
test_edge{i}.edge_duplicat(temp_commedge_idx{idx})=1; %1 mark as duplicated
test_edge{j}.edge_duplicat(temp_commedge_jdx{idx})=1;
end
idx=idx+1;
end
end
%select need test edge from each subgropup
clear needtest_edge
idx=1;
for i=1:lnof_newmodel
%test_edge{i}
needtest_idx=test_edge{i}.edge_duplicat==0 ; %select edges only available in one group
if ~isempty(test_edge{i}.edge_Pmark)
needtest_idx2=test_edge{i}.edge_Pmark==1 ; %select edges's P >cutoff P and test whether did it
%exist in the old needtest edge list and only available in one group, if not add them in.
needtest_idx20=test_edge{i}.edge_Pmark==0; %select edges's P <cutoff P and test whether did it
onlyPgtCutoff_idx=find(needtest_idx20==1);
lnof_onlyPgtCutoff=length(onlyPgtCutoff_idx);
for ki=1:lnof_onlyPgtCutoff
%delete P<cutoff_P from need test edge list
if ~isempty(test_edge{i}.edge_name)
[needtestEG.allname ti]=setdiff(needtestEG.allname,test_edge{i}.edge_name(onlyPgtCutoff_idx(ki)));
%needtestEG.allLocation=needtestEG.allLocation(ti,:);
end %end if
end
PgreatCutoff_idx=find(needtest_idx2==1& needtest_idx==1);
%PgreatCutoff_idx=find(needtest_idx==1);
lnof_newPgreatCutoff=length(PgreatCutoff_idx);
needtest_idx22=[];
for ki=1:lnof_newPgreatCutoff
if ~ismember(test_edge{i}.edge_name(PgreatCutoff_idx(ki)),needtestEG.allname);
needtest_idx(PgreatCutoff_idx(ki))=2; %assign 2 as full fill all the requirements
iid=PgreatCutoff_idx(ki); %added sep 2006
needtest_idx22(ki)=iid;
needtestEG.allname(end+1)=test_edge{i}.edge_name(PgreatCutoff_idx(ki));%add in if it did not in the list
else
iid=PgreatCutoff_idx(ki); %added sep 2006
needtest_idx(iid)=2;
needtest_idx22(ki)=iid;
end
end
%needtest_idx22=find(needtest_idx==2); %added sep 2006
needtest_edge{i}.model=test_edge{i}.model;
needtest_edge{i}.disp_string=test_edge{i}.disp_string(needtest_idx22);
needtest_edge{i}.name=test_edge{i}.edge_name(needtest_idx22);
needtest_edge{i}.P=test_edge{i}.edge_P(needtest_idx22);
needtest_edge{i}.df=test_edge{i}.edge_df(needtest_idx22);
needtest_edge{i}.device=test_edge{i}.edge_device(needtest_idx22);
%added 12.05
needtest_edge{i}.edge_isSDordering=test_edge{i}.edge_isSDordering(needtest_idx22);
end %~isempty
end
%added 11.27
lnof_newmodel=length(needtest_edge);
%delete P<cutoff from needtestEG
for i=1:lnof_newmodel
%needtest_edge{i}
%disp('test 1')
if ~isempty(needtest_edge{i})
testP=needtest_edge{i}.P<cutoff_P;
end
if ~isempty(testP) & testP>0
[needtestEG.allname ti]=setdiff(needtestEG.allname,{needtest_edge{i}.name{find(testP==1)}});
end
end
%added 12.03
%if it is conherance step, delete edge if not exist in the first round of needtestEG
%default set as the conherance learning step
%in conherence step, delete edge which was not existed in the first round of needtestEG
%OldneedtestEG
%for i=1:length(OldneedtestEG.allname)
% OldneedtestEG.allname(i)=cellstr(deblank(OldneedtestEG.allname(i)));
%end
%needtestEG
% idx_needtestEG=1;
% clear temp_needtestEG
% for i=1:length(needtestEG.allname)
% needtestEG.allname(i)=cellstr(deblank(needtestEG.allname(i)));
% tempEdge=strvcat(needtestEG.allname(i));
% if (ismember(cellstr(tempEdge),OldneedtestEG.allname) | ismember(cellstr(tempEdge([1,3,2,4])),OldneedtestEG.allname))
% temp_needtestEG.allname(idx_needtestEG)=cellstr(tempEdge);
% idx_needtestEG=idx_needtestEG+1;
% end
% end
% needtestEG.allname=temp_needtestEG.allname;
% needtestEG
%for i=1:length(needtest_edge)
% needtest_edge{i}
%end
%%end add 12.03
if isempty(needtestEG.allname)
disp(['Need test edge list is empty']);
disp(['Cut off P= ',num2str(cutoff_P)]);
disp(['Model selection finished, selected model: ']);
%disp(removeEG.new_model);
total
if ~isempty(removeEG.name)
marray_debuge(['Selected model= ']);
marray_debuge(strvcat(removeEG.new_model));
end
%marray_debuge(removeEG.new_model);
marray_debuge(['Total device= ', num2str(total.device)]);
marray_debuge(['Total degree freedom= ', num2str(total.df)]);
marray_debuge(['Total P= ', num2str(total.P)]);
selected_model=removeEG.new_model;
break;
end
%print list of results and select one edge to remove and make new subgroup
print_string='';
for i=1:lnof_newmodel
if ~isempty(needtest_edge{i})
print_string=strvcat(strvcat(needtest_edge{i}.disp_string),print_string);
end
end
%select the edge with max P to remove
max_P=0;
clear maxEG;
maxEG.P=0;
%maxEG.df=0; %added sep 2006
%maxEG.device=0;
%maxEG.name=[];
%maxEG.modelidx=0;
mi=0;
mj=0;
for i=1:lnof_newmodel
% needtest_edge{i}
% disp('test 2')
if (~isempty(needtest_edge{i}) & ~isempty(needtest_edge{i}.name))
%added 12.05
for j=1:length(needtest_edge{i}.name)
if needtest_edge{i}.edge_isSDordering(j)==1
if needtest_edge{i}.P(j)>max_P
max_P=needtest_edge{i}.P(j);
mi=i;
mj=j;
end
end
end
end
end
%pause
if max_P~=maxEG.P
maxEG.name=needtest_edge{mi}.name(mj);
maxEG.P=needtest_edge{mi}.P(mj);
maxEG.df=needtest_edge{mi}.df(mj);
maxEG.device=needtest_edge{mi}.device(mj);
maxEG.model=needtest_edge{mi}.model;
maxEG.modelidx=mi;
else %if (max_P==0 & maxEG.P==0) %nothint was found
maxEG.modelidx=0; %added sep 2006
end
%bug at here may 09 2002!!!!!!!!!!!!!
isfind=0;
%initial value
%temp_removeEG
%needtest_edge
temp_removeEG.new_model={};
% lnof_newmodel,isfind, needtest_edge{1},maxEG.model
if maxEG.modelidx~=0 %added sep 2006
for i=1:lnof_newmodel
if i~=maxEG.modelidx
%here is to find out the rest untested modelbug at here nov 12.2002
if ~isempty(needtest_edge{i})
if isfind==0
temp_removeEG.new_model{i}=needtest_edge{i}.model;
else
temp_removeEG.new_model{i-1}=needtest_edge{i}.model;
end
end
else
if strcmp(needtest_edge{i}.model,maxEG.model)
tempname=strvcat(maxEG.name);
need_removed_edgeStart=findstr('[',tempname);
need_removed_edgeEnd=findstr(']',tempname);
need_removed_edge1=tempname(need_removed_edgeStart+1);
%bug at here if more than 2 nov.11.2002
temp_model{1}=setdiff(maxEG.model,need_removed_edge1);
need_removed_edge2=tempname(need_removed_edgeEnd-1);
temp_model{2}=setdiff(maxEG.model,need_removed_edge2);
% temp_removeEG.new_model={};
% temp_removeEG.new_model{i}=temp_model{1};
% temp_removeEG.new_model{i+1}=temp_model{2};
isfind=1;
else
disp('Wrong please check');
end
end
end %end for
%%if lnof_newmodel<removeEG.new_model then try to detect the lost edges
%in temp_removeEG.newmodel
tempR=zeros(size(removeEG.new_model));
if lnof_newmodel<length(removeEG.new_model)
for r_jdx=1:length(needtest_edge)
tempR=tempR+ismember(removeEG.new_model,needtest_edge{r_jdx}.model);
end
temp_modelStr=removeEG.new_model(find(tempR==0));
for r_idx=1:length(temp_modelStr)
temp_removeEG.new_model{end+r_idx}=temp_modelStr{r_idx};
end
end
%test the new subgroup whether included in one of the old group
idx=1;
clear record_model;
for ik=1:2
temp_idx=0;
for jk=1:length(temp_removeEG.new_model)
% temp_model{ik},temp_removeEG.new_model{jk}
is_include=intersect(temp_model{ik},temp_removeEG.new_model{jk});
%pause
if ~strcmp(is_include,temp_model{ik})
%record_model{idx}=temp_model{ik}
temp_idx=temp_idx+1;
end
end
if temp_idx==length(temp_removeEG.new_model)
record_model{idx}=temp_model{ik};
idx=idx+1;
end
end
%record_model
%temp_removeEG.new_model
%test record_model belong to which model
%for r_idx=1:length(test_edge)
%ismember(temp_removeEG.new_model,test_edge{r_idx}.model)
%add new modle from record_model to temp_removeEG.new_model
for lk=1:length(record_model)
%disp('record model')
%record_model{lk}
%lk
temp_removeEG.new_model{end+1}=record_model{lk};
end
%update removeEG
%marray_debuge(print_string);
removeEG.new_model=temp_removeEG.new_model;
removeEG.P=maxEG.P;
removeEG.df=maxEG.df;
removeEG.device=maxEG.device;
removeEG.name=maxEG.name;
total.device=removeEG.device+total.device;
total.df=removeEG.df+total.df;
total.P=Chisq(total.device,total.df);
end %end if maxEG
%display results
if ~isempty(removeEG.name)
marray_debuge(['Remove Edge: ', strvcat(removeEG.name)]);
marray_debuge(['Degree freedom= ',num2str(removeEG.df)]);
marray_debuge(['Device= ',num2str(removeEG.device)]);
marray_debuge(['P= ',num2str(removeEG.P)]);
marray_debuge(['Location: ', num2str(removeEG.location)]);
marray_debuge(['Selected model: ']);
marray_debuge(strvcat(removeEG.new_model));
end
%marray_debuge(['Model: ', strvcat(removeEG.model)]);
marray_debuge(['Total device= ', num2str(total.device)]);
marray_debuge(['Total degree freedom= ', num2str(total.df)]);
marray_debuge(['Total P= ', num2str(total.P)]);
marray_debuge(' ');
%pause
if total.P<cutoff_P
%total
marray_debuge(['Total device= ', num2str(total.device)]);
marray_debuge(['Total degree freedom= ', num2str(total.df)]);
marray_debuge(['Total P= ', num2str(total.P)]);
marray_debuge(['Selected model= ']);
marray_debuge(strvcat(removeEG.new_model));
selected_model=removeEG.new_model;
break;
elseif (maxEG.modelidx==0 &isfind==0) %nothing could be found
marray_debuge(['Total device= ', num2str(total.device)]);
marray_debuge(['Total degree freedom= ', num2str(total.df)]);
marray_debuge(['Total P= ', num2str(total.P)]);
marray_debuge(['Selected model= ']);
marray_debuge(strvcat(removeEG.new_model));
selected_model=removeEG.new_model;
break;
end
%pause
end %end while
end %end if
selected_model=removeEG.new_model;
