global htxt2 ind1 init_fig ficonmessagesbg table table_fig plot_fig
global file1 outfname path1 data outdata graph_form
%bellow is object handle
global openBut1 htxt_cutoffP Brun htxt_RowOrCol RowOrCol_field Normalization_field deviceChoice_field isordered
set(0,'ShowHiddenHandles','on')
if ( ~exist('GUIFl') )
%close all hidden;
set(init_fig,'Name','MGraph - Graphical gaussian model');
clf(init_fig);
% if exist('ficonmessagesbg');
% delete(ficonmessagesbg);
% clear ficonmessagesbg;
%end
GUIFl=0;
plot_fig=11;
%Default input values
runcode=0;
cutoff_P=[];
RowOrCol=1;
deviceChoice=1;
algormChoice=1;
Normalization=1;
isordered=1;
graph_form=1;
%output
end;
if GUIFl==0
fmenu=uimenu('Label','MGraph');
uimenu(fmenu,'Label','MGraph Properties','Callback','MGraph_properties');
uimenu(fmenu,'Label','Graphical logline model','Callback',['clf(init_fig);MGraph_logline']);
uimenu(fmenu,'Label','Graphical gaussion model','Callback',[ 'clf(init_fig);MGraph_gauss']);
uimenu(fmenu,'Label','Quit','Callback',['clear all;close all hidden'],...
'Separator','on','Accelerator','Q');
dmenu=uimenu('Label','MGraph Debug');
uimenu(dmenu,'Label','MGraph Debug','Callback','marray_debug');
%Normalization=1;
%RowOrCol=1;
global t1 htxt
%if exist('l_openBut1')
% delete(l_openBut1);
% clear l_openBut1;
%end;
%if exist('l_htxt_cutoffP')
% delete(l_htxt_cutoffP);
% clear l_htxt_cutoffP;
%end;
%if exist('l_Brun')
% delete(l_Brun);
% clear l_Brun;
%end;
axHndl=init_fig;
set(init_fig,'Name','MGraph - Graphical gaussian model',...
'Position',[0.1 0.07 0.55 0.70]);
%text
htxt2 =uicontrol('Parent',axHndl,...% 'Units', 'normalized', ...
'Style', 'text', ...
'Position', [180 80 250 150], ... %[.3 .1 .55 0.7],...
'Backg', [1 1 1], ...
'Foreg', [0 0 0], ...
'String', [' Graphical gaussian model'], ...
'FontUnits', 'Normalized',...
'FontWeight', 'bold', ...
'HorizontalAlignment', 'center', ...
'Tag', 'genetext');
Brun=uicontrol('Parent',init_fig,'style','push','string','Run program',...
'position',[10 40 80 30 ], ...
'callback', ['runcode=1;clf(init_fig);MGraph_gauss;' ], ...
'Interruptible','on');
Bexport=uicontrol('Parent',init_fig,'style','push','string','Export to Pjart',...
'position',[10 10 80 30 ], ...
'callback', ['disp([''Export data '', path1,file1,''NET'','' ....'']);MGraph_export2Pjart([path1,file1,''NET''],outdata.labels,spareM,old_parcorrf);'], ...
'Interruptible','on');
%following is function
%Open Input file1
openBut1=uicontrol('Parent',init_fig,'style','push','string','Load input file',...
'position',[10 415 150 30], ...
'callback', [' [file1, path1]=uigetfile({''*.*'', ''All Files (*.*)''});clf(init_fig);MGraph_gauss;'],...
'Interruptible','on');
if exist('file1')
file_name1=uicontrol('Parent',init_fig,'Style','text','HorizontalAlignment','left',...
'Position',[10 375 150 30],...
'String', file1);
end
%Open Input real structure
openBut2=uicontrol('Parent',init_fig,'style','push','string','Load real pathway structure',...
'position',[210 415 150 30], ...
'callback', [' [file2, path2]=uigetfile({''*.*'', ''All Files (*.*)''});clf(init_fig);MGraph_gauss;'],...
'Interruptible','on');
if exist('file2')
file_name2=uicontrol('Parent',init_fig,'Style','text','HorizontalAlignment','left',...
'Position',[210 375 150 30],...
'String', file2);
end
%open input parameter
%%%%%%%
%cutoffP=uicontrol('Parent',properties_fig,'Style','popupmenu',...
% 'Position',[ 170 130 90 20],...
% 'String', 'Default(0) | Med(bk+2bk_std)| User select', ...
% 'Callback',['ch1_int_f2=get(ch1_intensity_file2,''Value'');']);
htxt_cutoffP=uicontrol('Parent',init_fig,... %'Units', 'normalized', ...
'Style', 'text', ...
'Position', [10 250 150 20], ... %[.02 .55 .24 0.05],...
'Backg', [1 1 1], ...
'Foreg', [0 0 0], ...
'String','Significance P=', ...
'FontUnits', 'Normalized',...
'FontWeight', 'bold', ...
'HorizontalAlignment', 'left', ...
'Tag', 'cutoffP','BackgroundColor','yellow');
cutoffP_field=uicontrol('Parent',init_fig,'Style','edit','Position', [10 230 150 20], ... %[.02 .48 .2 0.05],...
'BackgroundColor','white',...
'CallBack',['cutoff_P=str2num(get(cutoffP_field,''String''));']);
htxt_RowOrCol=uicontrol('Parent',init_fig,... % 'Units', 'normalized', ...
'Style', 'text', ...
'Position', [10 175 150 20], ... %[.02 .40 .24 0.05],...
'Backg', [1 1 1], ...
'Foreg', [0 0 0], ...
'String','Causal prediction', ...
'FontUnits', 'Normalized',...
'FontWeight', 'bold', ...
'HorizontalAlignment', 'left', ...
'Tag', 'cutoffP','BackgroundColor','yellow');
isordered_field=uicontrol('Parent',init_fig,'Style','popup',...
'Position',[10 320 100 30],...
'String', 'ordered variables|not ordered variables', ... %ch1=1, ch2=2
'Callback',['isordered=get(isordered_field,''Value'');']);
RowOrCol_field=uicontrol('Parent',init_fig,'Style','popup',...
'Position',[10 145 100 30],...
'String', 'Row variables|Column variables', ...
'Callback',['RowOrCol=get(RowOrCol_field,''Value'');']);
%ch1=1,ch2=2
Normalization_field=uicontrol('Parent',init_fig,'Style','popupmenu',...
'Position',[10 105 150 30],...
'String', 'Normalize row data to mean=0 variance=1|Normalize column data to mean=0 variance=1|Not normalize input data', ... %ch1=1, ch2=2
'Callback',['Normalization=get(Normalization_field,''Value'');']);
deviceChoice_field=uicontrol('Parent',init_fig,'Style','popupmenu',...
'Position',[10 65 150 30],...
'String', 'Device set as variance|junbai chi-square test', ...
'Callback',['deviceChoice=get(deviceChoice_field,''Value'');']);
algormChoice_field=uicontrol('Parent',init_fig,'Style','popupmenu',...
'Position',[10 285 150 30],...
'String', 'GGM forward algorizm|GGM backward algorizm|PC algorizm|IG forward algorizm|IG backward algorizm|IG forward fast algorizm|IG forward with Bscore|Gnet random restart|Gnet part corr coef', ...
'Callback',['algormChoice=get(algormChoice_field,''Value'');']);
if runcode==1 %start running
if exist('hwaite')
close(hwaite);
clear hwaite;
end
table=[];
marray_debuge(char(10));
marray_debuge('Start Running ...............');
%read input data
%for no lable input
%outdata=MGraph_loadDataforGauss([path1,file1]);
%labeled input
hwaite=waitbar(0,'Please wait....');
if exist('tblread')~=0
[data,varnames,casenames] = tblread([path1,file1],'\t') ;
net_data=[];net_varnames=[]; net_casenames=[];
if (exist('file2')& file2~=0)
net_data=[];net_varnames=[]; net_casenames=[];
[net_data,net_varnames,net_casenames] = tblread([path2,file2],'\t') ;
end
else
disp('No tblread available, use MGraph build in funcion read data...');
[data,varnames,casenames]=mgraph_tblread([path1,file1]);
net_data=[];net_varnames=[]; net_casenames=[];
if (exist('file2')& file2~=0)
[net_data,net_varnames,net_casenames] = tblread([path2,file2],'\t') ;
end
end
%check for causial prediction of row or column data
if ~exist('RowOrCol')
warndlg('Please select Row or Column variable then run the program again !','!! Warning !!')
end
if RowOrCol==1
outdata.labels=casenames;
outdata.data=data;
marray_debuge('Causal prediction of row variables');
elseif RowOrCol==2
outdata.labels=varnames;
outdata.data=data';
marray_debuge('Causal prediction of column variables');
else
error('Wrong in select Row or Column data');
end
set(RowOrCol_field,'Value',RowOrCol);
outdata.N=size(outdata.data,2);
outdata.q=size(outdata.data,1);
waitbar(1/4,hwaite);
%get cuttoff of P value
if ~isempty(cutoff_P)
set(cutoffP_field,'String',num2str(cutoff_P));
else
set(cutoffP_field,'String',num2str(0.05)); %default cut off P
cutoff_P=0.05;
marray_debuge('Use default P=0.05');
end
outdata.cutoff_P=cutoff_P;
%normalize raw data by mean 0 variance 1
%check normalization option
if ~exist('Normalization')
warndlg('Please select Normalization option then run the program again !','!! Warning !!')
end
if Normalization==1
tdata=som_normalize(outdata.data','var');
marray_debuge('Normalize row input data to mean=0, variance=1');
elseif Normalization==2
tdata=som_normalize(outdata.data,'var');
tdata=tdata';
marray_debuge('Normalize column input data to mean=0, variance=1');
else
tdata=outdata.data';
marray_debuge('Not normalize input data');
end
set(Normalization_field,'Value',Normalization);
marray_debuge(['Load file ',file1]);
marray_debuge(['Reading parameter' ]);
marray_debuge(['Number of data, N=',num2str(outdata.N)]);
marray_debuge(['Number of variables, q=',num2str(outdata.q)]);
marray_debuge(['Label of columns, ']);
marray_debuge([outdata.labels]);
marray_debuge(['Significance P=',num2str(outdata.cutoff_P)]);
old_var=cov(tdata);
N=outdata.N; % number of observations
q=outdata.q; % number of variables
cutoff=outdata.cutoff_P;
waitbar(2/4,hwaite);
%compute variance
%deviceChoice=0 use variance device
%deviceChoice=others use junbai's Chi-square test
%deviceChoice
if algormChoice==2
%Backward algorm
marray_debuge('Run Gaussian graphical model with Backward algorizm')
[old_parcorrf,delet_edge,record_edge]=MGraph_GaussModel(old_var,N,q,cutoff,deviceChoice);
elseif algormChoice==1
%forward algorm
marray_debuge('Run Gaussian graphical model with Forward algorizm')
[old_parcorrf,delet_edge,record_edge]=MGraph_GaussModel_forward(old_var,N,q,cutoff);
elseif algormChoice==3
pdag_wang = wang_learn_struct_pdag_pc('cond_indep_fisher_z', q, 5, old_var, N, cutoff);
elseif algormChoice==4
%IC forward
marray_debuge('Run IG model with Forward algorizm')
[old_parcorrf,delet_edge,record_edge,pdag_wang,G]=MGraph_IG_forward(old_var,N,q,cutoff);
elseif algormChoice==5
%IC backgward
marray_debuge('Run IG model with Backward algorizm')
[old_parcorrf,delet_edge,record_edge,pdag_wang,G]=MGraph_IG_Backward(old_var,N,q,cutoff,deviceChoice);
elseif algormChoice==6
%IC forward fast
marray_debuge('Run IG model with Forward fast algorizm')
[old_parcorrf,delet_edge,record_edge,pdag_wang,G]=MGraph_IG_forward_fast2(old_var,N,q,cutoff);
elseif algormChoice==7
%IG with Bscore and IG % this one was not stable
marray_debuge('Run IG plus BScore search algorizm')
[old_parcorrf,delet_edge,record_edge,pdag_wang,G]=MGraph_IG_forwardFast_BGscore(tdata,old_var,N,q,cutoff);
elseif algormChoice==8
% Bscore search
marray_debuge('Run BScore search algorizm with random restart')
%tG=zeros(q,q);
tG=setdiag(triu(ones(q,q)),0).*(-1);
%for ii=1:q-1
% for jj=ii+1:ii+1
% tG(ii,jj)=-1;
% end
%end
[old_parcorrf,pdag_wang ,GG, maxScore]=MGraph_BGscore_hillClimb2_rs(tdata,tG,0.15);
%test average 12.15
%[old_parcorrf,pdag_wang ,GG, maxScore]=MGraph_BGscore_average(tdata,tG);
%end test 12.15
elseif algormChoice==9
marray_debuge('Run BScore search algorizm with learning rate plus partial corr. coef.')
tG=setdiag(triu(ones(q,q)),0).*(-1);
[old_parcorrf,pdag_wang ,GG, maxScore]=MGraph_BGscore_hillClimb2(tdata,tG,0.15);
end
waitbar(3/4,hwaite);
set(deviceChoice_field,'Value',deviceChoice);
set(algormChoice_field,'Value',algormChoice);
set(isordered_field,'Value',isordered);
if algormChoice<3
marray_debuge(['Estimated partial correlation coefficient= ']);
%replace small valus as 0
tempM=(abs(old_parcorrf)>0.0000000001);
temp_old_parcorrf=tempM.*old_parcorrf;
if exist('vpa')~=0
marray_debuge([num2str(double((temp_old_parcorrf)))]);
end
marray_debuge(['Number of deleted edge= ', num2str(delet_edge)]);
%marray_debuge(['Deleted edges= ']);
%marray_debuge([ num2str(record_edge)]);
end
if algormChoice<3
%draw figures
spareM=(abs(old_parcorrf)>0.0000000001 & abs(old_parcorrf)~=1);
%spareM=sparse(spareM);
set(init_fig,'Visible','off');
figure(plot_fig);
clf(plot_fig);
labels=strrep(cellstr( deblank(outdata.labels) ),'"','');
draw_graph(spareM,lower(labels));
set(init_fig,'Visible','on');
elseif (algormChoice==3 | algormChoice==4 | algormChoice==5 | algormChoice==6 | algormChoice==7| algormChoice==8| algormChoice==9)
spareM=pdag_wang
set(init_fig,'Visible','off');
figure(plot_fig);
clf(plot_fig);
labels=strrep(cellstr( deblank(outdata.labels) ),'"','');
draw_graph(abs(spareM),lower(labels));
set(init_fig,'Visible','on');
end
%start to compute the error rate
%compare net_Data and abs(spareM)
if (exist('file2')& file2~=0)
AA=[]; BB=[];ispath_C=[];
rel_A=[];rel_B=[];
edge_existence_error_c=0;
predicted_net=abs(spareM);
[AA BB]=find(predicted_net);
[rel_A rel_B]=find(net_data);
len_rel_A=length(rel_A);
%edge existence error of commission
len_AA=length(AA);
ispath_C = reachability_graph(abs(net_data));
for ii=1:len_AA
temp_edge=[AA(ii), BB(ii)];
if ( net_data(temp_edge(1),temp_edge(2))==0 & (ispath_C(AA(ii),BB(ii))==0 & ispath_C(BB(ii),AA(ii))==0 ) )
edge_existence_error_c=edge_existence_error_c+1;
end
end
ratio_of_edge_existence_error_c=edge_existence_error_c/len_rel_A
%edge direction error of commission
len_AA=length(AA);
edge_direction_error_c=0;
for ii=1:len_AA
temp_edge=[AA(ii), BB(ii)];
if ( abs(net_data(temp_edge(1),temp_edge(2))) + abs(net_data(temp_edge(2),temp_edge(1))) )>0
if abs(net_data(temp_edge(1),temp_edge(2)))~=1
edge_direction_error_c=edge_direction_error_c+1;
end
end
end
ratio_of_edge_direction_error_c=edge_direction_error_c/len_rel_A
%edge existence error of ommission
AA=[]; BB=[]; ispath_P=[];
edge_existence_error_o=0;
ispath_P = reachability_graph(abs(predicted_net));
[AA BB]=find(setdiag(predicted_net,1)==0);
len_AA=length(AA);
for ii=1:len_AA
temp_edge=[AA(ii),BB(ii)];
%net_data(temp_edge)
%pause
if ( net_data(temp_edge(1),temp_edge(2))==-1 & (ispath_P(AA(ii),BB(ii))==0 & ispath_P(BB(ii),AA(ii))==0) )
edge_existence_error_o= edge_existence_error_o+1;
end
end
ratio_of_edge_existence_error_o=edge_existence_error_o/len_rel_A
end
%end of error rate
waitbar(1,hwaite);
close(hwaite);
clear hwaite;
runcode=0;
end %end run
end; %end all
