function Connectivity_simulator_01
% Tool for simulation of connectivity in fMRI data. This tool enables to
% compare two methods for estimation of effective connectivity - Dynamic 
% causal modelling (DCM) and Granger causal modelling (GCM). This function
% initializes graphical interface of the tool. For more informations press
% "Help" button in the interface.
%
% -------------------------------------------------------------------------
% Created in MATLAB 7.8.0.347 (R2009a) with SPM8 toolbox and GCCA toolbox (Anil Seth code)
% -------------------------------------------------------------------------
% Written by Martin Gajdoš
% -------------------------------------------------------------------------

spm_defaults

ScreenSize = get(0,'ScreenSize');
h.gray_1=[0.9412 0.9412 0.9412]; % Define type of gray color (240/255)

% GUI window initialization
h.width = ScreenSize(3)-200;
h.height = ScreenSize(4)-200;
h.fpos = round((ScreenSize(3:4)-[h.width h.height])/2);
h.f = figure('Name','Connectivity_simulator','NumberTitle','off','Units','pixels','Position',[h.fpos h.width/3 h.height],'Resize','off','WindowStyle','normal', 'Color',h.gray_1,'MenuBar','none');

% Default settings of variables
h.M(1).A=[-0.5    0    0; ...   % Connectivity matrix
          0.9  -0.5    0; ...   
          0.7    0.8  -0.5];   

h.M(1).d=[0   0    0;      % Delay matrix [ms]
          50  0    0;     
          50   50   0]; 


h.ROI=3;                  % Number of ROIs
h.SNR=5;                  % Signal to noise ratio
h.TR =2;                  % Repeat time [s]
h.dt =5;                  % Time step [ms]
h.M(1).dt=h.dt/1000;      % Time step [s]
h.N  =200;                % Number of time points (BOLD signal after downsampling with TR)
h.T  =h.TR*h.N/h.M(1).dt; % Number of elements in fine sampling
h.M(1).Ud_dt = 0.1;       % Time step for DCM estimation 
h.inputs = 1;             % Number of inputs - stimulation vectors
h.no_of_events = 15;      % Number of stimuli for stimulation vector(default)
h.equal_length = 5;       % If selected, all stimulations in vector will last this time [s]
%---------------------------------------------------------------------
% input control (matrix)
h.M(1).C = zeros(h.ROI,h.ROI);
h.M(1).C(1,1)=1;          % input to the first node
%----------------------------------------------------------------------
% Stimulation vector
h.event_duration = 500;   % Shortest duration of stimuli[ms] 
% Positions of every event - scans or secs according to flag in CreateU
h.Pos = [2 11 34 48.6 60 90 160 200 240 280 300 320 340 360 380];
% Durations of every event - scans or secs according to flag in CreateU
h.Dur = [2 10 3 3 20 15 15 15 15 5 5 5 5 5 5 ];
% Amplitude of stimulation vector
h.Ampl = [1];
% Creating stimulation vector
h.xU = CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).dt,'secs', h.event_duration);
% Specifying the input
U = repmat(h.xU',h.ROI,1);
h.M(1).U = U;
% Creating stimulation vector for DCM estimation
h.xUd = CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).Ud_dt,'seconds', h.event_duration);
h.M(1).Ud = repmat(h.xUd,1,h.ROI);
% --------- GUI panel 1 - SIMULATION PARAMETERS - initialization ----------
h.panel_1 = uipanel('Title','Conditions for simulated data','FontSize',14,'FontWeight','bold','BackgroundColor',h.gray_1,'Position',[0.01 0.01 0.98 0.97]);

h.insert_vector = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.02 0.92 0.37 0.05],'Style','PushButton','String','Specify stimulation','FontSize',12,...
    'CallBack',@insert_vector,'BackGroundColor',h.gray_1);
h.defautlt_settings = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.4 0.92 0.31 0.05],'Style','PushButton','String','Default settings','FontSize',12,...
    'CallBack',@default_settings,'BackGroundColor',h.gray_1);
h.load_results2 = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.72 0.92 0.26 0.05],'Style','PushButton','String','Load results','FontSize',12,...
    'CallBack',@load_results2,'BackGroundColor',h.gray_1);
h.method = uibuttongroup('Parent',h.panel_1,'Title','Generating method','FontSize',12,'FontWeight','bold','BackgroundColor',h.gray_1,'Units','normalized','Position',[0.02 0.71 0.48 0.19],'SelectionChangeFcn',@method_change);
h.method_DCM = uicontrol('Style','Radio','String','DCM','FontSize',12,'parent',h.method,'Units','normalized','pos',[0.2 0.5 0.5 0.4],'BackgroundColor',h.gray_1);
h.method_GCM = uicontrol('Style','Radio','String','GCM','FontSize',12,'parent',h.method,'Units','normalized','pos',[0.2 0.1 0.5 0.4],'BackgroundColor',h.gray_1);
h.No_inputs_ctrl = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.83 0.85 0.1 0.04],'Style','Text','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',12,'HorizontalAlignment','center','string',num2str(h.inputs));
h.No_inputs_text = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.53 0.84 0.3 0.05],'Style','text',...
    'String','No. of inputs: ','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','right');
h.ROI_text = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.53 0.79 0.3 0.05],'Style','text',...
    'String','No. of ROI''s: ','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','right');
h.ROI_ctrl = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.83 0.805 0.1 0.04],'Style','Edit','FontWeight','normal','FontSize',12,'HorizontalAlignment','center','CallBack',@ROI_control,'BackGroundColor',[1 1 1],'string',num2str(h.ROI));
h.scan_ctrl = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.83 0.755 0.1 0.04],'Style','Edit','FontWeight','normal','FontSize',12,'HorizontalAlignment','center','CallBack',@scan_control,'BackGroundColor',[1 1 1],'string',num2str(h.N));
h.scan_text = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.53 0.74 0.3 0.05],'Style','text',...
    'String','No. of scans: ','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','right');
h.SNR_text = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.53 0.69 0.3 0.05],'Style','text',...
    'String','SNR: ','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','right');
h.SNR_ctrl = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.83 0.705 0.1 0.04],'Style','Edit','FontWeight','normal','FontSize',12,'HorizontalAlignment','center','CallBack',@SNR_control,'BackGroundColor',[1 1 1],'string',num2str(h.SNR));
h.TR_ctrl = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.83 0.655 0.1 0.04],'Style','Edit','FontWeight','normal','FontSize',12,'HorizontalAlignment','center','CallBack',@TR_control,'BackGroundColor',[1 1 1],'string',num2str(h.TR));
h.TR_text = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.53 0.64 0.3 0.05],'Style','text',...
    'String','TR: ','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','right');
h.TR_s = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.93 0.64 0.06 0.05],'Style','text',...
    'String','[s] ','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','left');
h.time_step_ctrl = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.3 0.655 0.1 0.04],'Style','Edit','FontWeight','normal','FontSize',12,'HorizontalAlignment','center','CallBack',@time_step_control,'BackGroundColor',[1 1 1],'string',num2str(h.dt));
h.time_step_ms = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.4 0.64 0.1 0.05],'Style','text',...
    'String','[ms] ','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','left');
h.time_step_text = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0 0.64 0.3 0.05],'Style','text',...
    'String','Time step: ','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','right');
h.panel_1_1_c_matrix = uipanel('Parent',h.panel_1,'FontSize',12,'FontWeight','bold','BackgroundColor',h.gray_1,'Position',[0.02 0.38 0.96 0.25]);
h.connectivity_matrix_string=uicontrol('Parent',h.panel_1_1_c_matrix,'Units','normalized','Position',[0 0.85 0.5 0.15],'Style','text',...
    'String','Connectivity matrix','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','center');
h.connectivity_matrix =  uitable('Parent',h.panel_1_1_c_matrix,'Data',h.M(1).A,'Units','normalized','Position',[0.05 0.1 0.43 0.73],'BackGroundColor',[1 1 1],'FontWeight','normal','FontSize',8,'ColumnWidth',{35},'ColumnFormat',repmat(cellstr('bank'),1,h.ROI),'CellEditCallback',@connectivity_m_fcn,'ColumnEditable', logical(ones(1,h.ROI)));
h.from_1= uicontrol('Parent',h.panel_1_1_c_matrix,'Units','normalized','Position',[0.05 0 0.43 0.11],'Style','text',...
    'String','From','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',8,'FontWeight','bold','HorizontalAlignment','center');
h.to_1= uicontrol('Parent',h.panel_1_1_c_matrix,'Units','normalized','Position',[0 0.4 0.05 0.11],'Style','text',...
    'String','To','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',8,'FontWeight','bold','HorizontalAlignment','left');
h.delay_matrix_string=uicontrol('Parent',h.panel_1_1_c_matrix,'Units','normalized','Position',[0.53 0.85 0.47 0.15],'Style','text',...
    'String','Delay matrix [ms]','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','center');
h.delay_matrix =  uitable('Parent',h.panel_1_1_c_matrix,'Data',h.M(1).d,'Units','normalized','Position',[0.55 0.1 0.43 0.73],'BackGroundColor',[1 1 1],'FontWeight','normal','FontSize',8,'ColumnWidth',{35},'ColumnFormat',repmat(cellstr('short'),1,h.ROI),'CellEditCallback',@delay_m_fcn,'ColumnEditable', logical(ones(1,h.ROI)));
h.from_2= uicontrol('Parent',h.panel_1_1_c_matrix,'Units','normalized','Position',[0.55 0 0.43 0.11],'Style','text',...
    'String','From','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',8,'FontWeight','bold','HorizontalAlignment','center');
h.to_1= uicontrol('Parent',h.panel_1_1_c_matrix,'Units','normalized','Position',[0.5 0.4 0.05 0.11],'Style','text',...
    'String','To','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',8,'FontWeight','bold','HorizontalAlignment','left');

h.panel_1_2_model_of_interactions = uipanel('Parent',h.panel_1,'Title','Model of interactions','FontSize',12,'FontWeight','bold','BackgroundColor',h.gray_1,'Position',[0.02 0.07 0.96 0.29]);
h.ax1 = axes('Units','normalized','Parent',h.panel_1_2_model_of_interactions,'Position',[0 0 1 1]);
subplot(h.ax1)
cca_plotcausality(h.M(1).A,[],1);

h.run = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.02 0.01 0.3 0.05],'Style','PushButton','String','Run','FontSize',12,...
    'CallBack',@simulation,'BackGroundColor',h.gray_1);
h.panel_1_3_state = uipanel('Parent',h.panel_1,'BackgroundColor','green','Position',[0.35 0.01 0.3 0.05]);
h.state=uicontrol('Parent',h.panel_1_3_state,'Units','normalized','Position',[0 0 1 0.8],'Style','Text','String','Ready','FontSize',10,'HorizontalAlignment','center','backgroundcolor','green','Fontweight','bold');
h.help = uicontrol('Parent',h.panel_1,'Units','normalized','Position',[0.68 0.01 0.3 0.05],'Style','PushButton','String','Help','FontSize',12,...
    'CallBack',@help,'BackGroundColor',h.gray_1);

% ---------------- GUI panel 2 - RESULTS -initialization ------------------

h.panel_2 = uipanel('Title','Results','FontSize',14,'FontWeight','bold','BackgroundColor',h.gray_1,'Position',[0.32 0.01 0.67 0.97],'Visible','off');
h.ax2 = axes('Units','normalized','Parent',h.panel_2,'Position',[0.05 0.71 0.9 0.29]);
h.x_scans_text =uicontrol('Parent',h.panel_2,'Units','normalized','Position',[0.465 0.66 0.05 0.02],'Style','text','String','Scans','FontSize',10,'HorizontalAlignment','center');
h.x_BOLD_text =uicontrol('Parent',h.panel_2,'Units','normalized','Position',[0 0.7 0.02 0.2],'Style','text','String','BOLD','FontSize',10);
h.load_results = uicontrol('Parent',h.panel_2,'Units','normalized','Position',[0.67 0.25 0.32 0.05],'Style','PushButton','String','Load results','FontSize',12,...
    'CallBack',@load_results,'BackGroundColor',h.gray_1);
h.save_results = uicontrol('Parent',h.panel_2,'Units','normalized','Position',[0.67 0.19 0.32 0.05],'Style','PushButton','String','Save results','FontSize',12,...
    'CallBack',@save_results,'BackGroundColor',h.gray_1);
h.version_text = uicontrol('Parent',h.panel_2,'Units','normalized','Position',[0.69 0.02 0.3 0.14],'Style','Text','String',{['Version:'];['Connectivity_simulator_DP'];[''];['Last update:'];['15.5.2012 11:05'];},'FontSize',10,...
   'BackGroundColor',h.gray_1,'HorizontalAlignment','Left');
h.panel_2_1_DCM = uipanel('Parent',h.panel_2,'Title','DCM analysis','FontSize',12,'FontWeight','bold','BackgroundColor',h.gray_1,'Position',[0.005 0.34 0.66 0.33]);

% ------------------- DCM results panel -----------------------

h.DCM_A_matrix =  uitable('Parent',h.panel_2_1_DCM,'Units','normalized','Position',[0.04 0.25 0.275 0.6],'BackGroundColor',[1 1 1],'FontWeight','normal','FontSize',8,'ColumnWidth',{35},'ColumnFormat',repmat(cellstr('bank'),1,h.ROI),'ColumnEditable', logical(ones(1,h.ROI)));
h.DCM_A_matrix_string=uicontrol('Parent',h.panel_2_1_DCM,'Units','normalized','Position',[0.02 0.9 0.3 0.1],'Style','text',...
    'String','Connectivity matrix','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','center');
h.from_3= uicontrol('Parent',h.panel_2_1_DCM,'Units','normalized','Position',[0.04 0.18 0.275 0.07],'Style','text',...
    'String','From','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',8,'FontWeight','bold','HorizontalAlignment','center');
h.to_3= uicontrol('Parent',h.panel_2_1_DCM,'Units','normalized','Position',[0.01 0.425 0.03 0.15],'Style','text',...
    'String','To','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',8,'FontWeight','bold','HorizontalAlignment','left');

h.DCM_Ap_matrix =  uitable('Parent',h.panel_2_1_DCM,'Units','normalized','Position',[0.35 0.25 0.275 0.6],'BackGroundColor',[1 1 1],'FontWeight','normal','FontSize',8,'ColumnWidth',{35},'ColumnFormat',repmat(cellstr('bank'),1,h.ROI),'ColumnEditable', logical(ones(1,h.ROI)));
h.DCM_Ap_matrix_string=uicontrol('Parent',h.panel_2_1_DCM,'Units','normalized','Position',[0.35 0.9 0.275 0.1],'Style','text',...
    'String','Probability matrix','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','center');
h.from_4= uicontrol('Parent',h.panel_2_1_DCM,'Units','normalized','Position',[0.35 0.18 0.275 0.07],'Style','text',...
    'String','From','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',8,'FontWeight','bold','HorizontalAlignment','center');
h.to_4= uicontrol('Parent',h.panel_2_1_DCM,'Units','normalized','Position',[0.32 0.425 0.03 0.15],'Style','text',...
    'String','To','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',8,'FontWeight','bold','HorizontalAlignment','left');
h.DCM_A_MOI=uicontrol('Parent',h.panel_2_1_DCM,'Units','normalized','Position',[0.625 0.9 0.375 0.1],'Style','text',...
        'String','Model of interactions','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','center');
h.ax3 = axes('Units','normalized','Parent',h.panel_2_1_DCM,'Position',[0.625 0.05 0.375 0.8]);

h.more_results_DCM= uicontrol('Parent',h.panel_2_1_DCM,'Units','normalized','Position',[0.04 0.05 0.275 0.125],'Style','PushButton','String','More results','FontSize',12,...
    'CallBack',@more_results_DCM,'BackGroundColor',h.gray_1);
h.slider_probability_DCM= uicontrol('Parent',h.panel_2_1_DCM,'Units','normalized','Position',[0.58 0.05 0.04 0.126],'Style','Slider',...
    'CallBack',@threshold_DCM,'BackGroundColor',h.gray_1);
h.slider_edit_probability_DCM= uicontrol('Parent',h.panel_2_1_DCM,'Units','normalized','Position',[0.52 0.05 0.05 0.125],'Style','edit',...
    'CallBack',@threshold_DCM,'BackGroundColor','white','String',num2str(0.9));
h.slider_text_DCM= uicontrol('Parent',h.panel_2_1_DCM,'Units','normalized','Position',[0.38 0.03 0.13 0.125],'Style','Text','String','Threshold:','FontSize',10,'HorizontalAlignment','right',...
    'BackGroundColor',h.gray_1);

% ------------------ GCM results panel --------------------

h.panel_2_2_GCM = uipanel('Parent',h.panel_2,'Title','GCM analysis','FontSize',12,'FontWeight','bold','BackgroundColor',h.gray_1,'Position',[0.005 0.01 0.66 0.33]);

h.GCM_A_matrix =  uitable('Parent',h.panel_2_2_GCM,'Units','normalized','Position',[0.04 0.25 0.275 0.6],'BackGroundColor',[1 1 1],'FontWeight','normal','FontSize',8,'ColumnWidth',{35},'ColumnFormat',repmat(cellstr('bank'),1,h.ROI),'ColumnEditable', logical(ones(1,h.ROI)));
h.GCM_A_matrix_string=uicontrol('Parent',h.panel_2_2_GCM,'Units','normalized','Position',[0.02 0.9 0.3 0.1],'Style','text',...
    'String','Connectivity matrix','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','center');
h.from_5= uicontrol('Parent',h.panel_2_2_GCM,'Units','normalized','Position',[0.04 0.18 0.275 0.07],'Style','text',...
    'String','From','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',8,'FontWeight','bold','HorizontalAlignment','center');
h.to_5= uicontrol('Parent',h.panel_2_2_GCM,'Units','normalized','Position',[0.01 0.425 0.03 0.15],'Style','text',...
    'String','To','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',8,'FontWeight','bold','HorizontalAlignment','left');

h.GCM_Ap_matrix =  uitable('Parent',h.panel_2_2_GCM,'Units','normalized','Position',[0.35 0.25 0.275 0.6],'BackGroundColor',[1 1 1],'FontWeight','normal','FontSize',8,'ColumnWidth',{35},'ColumnFormat',repmat(cellstr('short'),1,h.ROI),'ColumnEditable', logical(ones(1,h.ROI)));
h.GCM_Ap_matrix_string=uicontrol('Parent',h.panel_2_2_GCM,'Units','normalized','Position',[0.35 0.9 0.275 0.1],'Style','text',...
    'String','Significance matrix','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','center');
h.from_6 = uicontrol('Parent',h.panel_2_2_GCM,'Units','normalized','Position',[0.35 0.18 0.275 0.07],'Style','text',...
    'String','From','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',8,'FontWeight','bold','HorizontalAlignment','center');
h.to_6 = uicontrol('Parent',h.panel_2_2_GCM,'Units','normalized','Position',[0.32 0.425 0.03 0.15],'Style','text',...
    'String','To','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',8,'FontWeight','bold','HorizontalAlignment','left');
h.GCM_A_MOI = uicontrol('Parent',h.panel_2_2_GCM,'Units','normalized','Position',[0.625 0.9 0.375 0.1],'Style','text',...
        'String','Model of interactions','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','center');
h.ax4 = axes('Units','normalized','Parent',h.panel_2_2_GCM,'Position',[0.625 0.05 0.375 0.8]);

h.more_results_GCM = uicontrol('Parent',h.panel_2_2_GCM,'Units','normalized','Position',[0.04 0.05 0.275 0.125],'Style','PushButton','String','More results','FontSize',12,...
    'CallBack',@more_results_GCM,'BackGroundColor',h.gray_1,'Min',0,'Max',1,'SliderStep',[0.01 0.1]);
h.slider_probability_GCM = uicontrol('Parent',h.panel_2_2_GCM,'Units','normalized','Position',[0.58 0.05 0.04 0.126],'Style','Slider',...
    'CallBack',@significance_GCM,'BackGroundColor',h.gray_1);
h.slider_edit_probability_GCM = uicontrol('Parent',h.panel_2_2_GCM,'Units','normalized','Position',[0.52 0.05 0.05 0.125],'Style','edit',...
    'CallBack',@significance_GCM,'BackGroundColor','white','String',num2str(0.01));
h.slider_text_GCM = uicontrol('Parent',h.panel_2_2_GCM,'Units','normalized','Position',[0.32 0.03 0.19 0.125],'Style','Text','String','Significant at P <','FontSize',10,'HorizontalAlignment','right',...
    'BackGroundColor',h.gray_1);

h.panel_2_3_comparison      = uipanel('Parent',h.panel_2,'Title','Comparison','FontSize',12,'FontWeight','bold','BackgroundColor',h.gray_1,'Position',[0.67 0.34 0.32 0.33]);
h.used_method_text          = uicontrol('Parent',h.panel_2_3_comparison,'Units','normalized','Position',[0.1 0.35 0.8 0.1],'Style','text','String','Generated by DCM','FontSize',10,'HorizontalAlignment','center','backgroundcolor',h.gray_1,'Fontweight','bold','visible','off');
h.structural_similarity_text= uicontrol('Parent',h.panel_2_3_comparison,'Units','normalized','Position',[0.1 0.9 0.8 0.1],'Style','text','String','Structural similarity','FontSize',10,'HorizontalAlignment','center','backgroundcolor',h.gray_1,'Fontweight','bold','visible','on');
h.wrong_connections_text    = uicontrol('Parent',h.panel_2_3_comparison,'Units','normalized','Position',[0.05 0.8 0.9 0.1],'Style','text','String','No. of wrong connections estimated by:','FontSize',10,'HorizontalAlignment','center','backgroundcolor',h.gray_1,'Fontweight','normal','visible','on');
h.DCM_text   = uicontrol('Parent',h.panel_2_3_comparison,'Units','normalized','Position',[0.2 0.65 0.2 0.1],'Style','text','String','DCM :','FontSize',10,'HorizontalAlignment','center','backgroundcolor',h.gray_1,'Fontweight','normal','visible','on');
h.DCM_No     =  uicontrol('Parent',h.panel_2_3_comparison,'Units','normalized','Position',[0.4 0.65 0.1 0.1],'Style','edit','BackGroundColor','white','enable','inactive');
h.GCM_text   = uicontrol('Parent',h.panel_2_3_comparison,'Units','normalized','Position',[0.2 0.5 0.2 0.1],'Style','text','String','GCM :','FontSize',10,'HorizontalAlignment','center','backgroundcolor',h.gray_1,'Fontweight','normal','visible','on');
h.GCM_No     =  uicontrol('Parent',h.panel_2_3_comparison,'Units','normalized','Position',[0.4 0.5 0.1 0.1],'Style','edit','BackGroundColor','white','enable','inactive');
h.sum_of_squares = uicontrol('Parent',h.panel_2_3_comparison,'Units','normalized','Position',[0.05 0.25 0.9 0.1],'Style','text','String','Sum of squares of residuals:','FontSize',10,'HorizontalAlignment','center','backgroundcolor',h.gray_1,'Fontweight','normal','visible','on');
h.sum_of_squares_text   = uicontrol('Parent',h.panel_2_3_comparison,'Units','normalized','Position',[0.05 0.1 0.65 0.1],'Style','text','String','DCM model to estimation:','FontSize',10,'HorizontalAlignment','center','backgroundcolor',h.gray_1,'Fontweight','normal','visible','on');
h.sum_of_squares_No     =  uicontrol('Parent',h.panel_2_3_comparison,'Units','normalized','Position',[0.7 0.1 0.2 0.1],'Style','edit','BackGroundColor','white','enable','inactive');

% ---------- GUI panel - 3 STIMULATION VECTOR - Initialization ------------

h.panel_3 = uipanel('Title','Stimulation vector','FontSize',14,'FontWeight','bold','BackgroundColor',h.gray_1,'Position',[0.32 0.01 0.67 0.97],'Visible','off');
h.list_of_vectors = uicontrol('Parent',h.panel_3,'String',cellstr('Vector 1'),'Units','normalized','Position',[0.02 0.8 0.2 0.1],'Style','Listbox','BackGroundColor',[1 1 1],'Max',2,'Min',0,'Value',[]);
h.list_of_vectors_text = uicontrol('Parent',h.panel_3,'Units','normalized','Position',[0.02 0.9 0.2 0.05],'Style','text','String','List of stimulation vectors','FontSize',10,'HorizontalAlignment','center','BackGroundColor',h.gray_1,'FontWeight','bold');
h.no_of_events_ctrl = uicontrol('Parent',h.panel_3,'Units','normalized','Position',[0.17 0.708 0.05 0.05],'Style','Edit','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',12,'HorizontalAlignment','center','CallBack',@no_of_events_control,'BackGroundColor',[1 1 1],'string',num2str(h.no_of_events));
h.no_of_events_text = uicontrol('Parent',h.panel_3,'Units','normalized','Position',[0.02 0.7 0.14 0.05],'Style','text',...
    'String','No. of events: ','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','right');
h.equal_lengths_button = uicontrol('Parent',h.panel_3,'Units','normalized','Position',[0.02 0.59 0.14 0.05],'Style','Pushbutton','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',12,'HorizontalAlignment','center','CallBack',@equal_lengths_button,'BackGroundColor',h.gray_1, 'String','Equal lenghts');
h.equal_lengths_ctrl = uicontrol('Parent',h.panel_3,'Units','normalized','Position',[0.17 0.59 0.05 0.05],'Style','Edit','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',12,'HorizontalAlignment','center','CallBack',@equal_lengths,'BackGroundColor',[1 1 1],'string',num2str(h.equal_length),'enable','off');
h.amplitude_ctrl = uicontrol('Parent',h.panel_3,'Units','normalized','Position',[0.17 0.648 0.05 0.05],'Style','Edit','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',12,'HorizontalAlignment','center','CallBack',@amplitude,'BackGroundColor',[1 1 1],'string',num2str(h.Ampl));
h.amplitude_text = uicontrol('Parent',h.panel_3,'Units','normalized','Position',[0.02 0.64 0.14 0.05],'Style','text',...
    'String','Amplitude: ','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','right');
h.uniform_events_button = uicontrol('Parent',h.panel_3,'Units','normalized','Position',[0.02 0.52 0.2 0.05],'Style','Pushbutton','BackGroundColor',h.gray_1,'FontWeight','normal','FontSize',12,'HorizontalAlignment','center','CallBack',@uniform_distribution,'BackGroundColor',h.gray_1, 'String','Uniform events');
h.select_input_text = uicontrol('Parent',h.panel_3,'Units','normalized','Position',[0.02 0.45 0.2 0.05],'Style','text',...
    'String','Area of input: ROI 1','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','center');
h.pos_table_text = uicontrol('Parent',h.panel_3,'Units','normalized','Position',[0.31 0.92 0.6 0.05],'Style','text',...
    'String','Onset of every stimulus [s]','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','center');
h.pos_table = uitable('Parent',h.panel_3,'data',h.Pos,'Units','normalized','Position',[0.31 0.8 0.6 0.12],'BackGroundColor',[1 1 1],'FontWeight','normal','FontSize',8,'ColumnWidth',{35},'ColumnEditable',logical(ones(1,h.no_of_events)),'CellEditCallback',@Pos_control);
h.dur_table_text = uicontrol('Parent',h.panel_3,'Units','normalized','Position',[0.31 0.72 0.6 0.05],'Style','text',...
    'String','Duration of every stimulus [s]','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','center');
h.dur_table = uitable('Parent',h.panel_3,'data',h.Dur,'Units','normalized','Position',[0.31 0.6 0.6 0.12],'BackGroundColor',[1 1 1],'FontWeight','normal','FontSize',8,'ColumnWidth',{35},'ColumnEditable', logical(ones(1,h.no_of_events)),'CellEditCallback',@Dur_control);
h.ax5_text = uicontrol('Parent',h.panel_3,'Units','normalized','Position',[0.05 0.4 0.9 0.05],'Style','text',...
    'String','Stimulation vector','BackGroundColor',h.gray_1,'FontWeight','bold','FontSize',12,'HorizontalAlignment','center');
h.ax5 = axes('Units','normalized','Parent',h.panel_3,'Position',[0.05 0.1 0.9 0.28]);
h.ax5_x_label = uicontrol('Parent',h.panel_3,'Units','normalized','Position',[0.25 0 0.7 0.05],'Style','text',...
    'String','Time [s]','BackGroundColor',h.gray_1,'FontSize',10,'HorizontalAlignment','center');
subplot(h.ax5)
plot(linspace(h.M(1).dt,h.N*h.TR,h.T),h.xU)
grid on
axis([0 h.N*h.TR -0.1 1.1])

guidata(h.f,h)

function Pos_control(handle, eventdata)
% Function is used to change onset positions of events during specifying 
% stimulation vector. The unit of positions in [s]

h = guidata(handle);

time=h.N*h.TR;
pos=eventdata.Indices;
new=eventdata.NewData;
h.Pos(pos(1),pos(2))=new;
h.Pos=sort(h.Pos);

set(h.pos_table, 'Data', h.Pos)

for i=1:h.no_of_events-1
    if h.Pos(i)+h.Dur(i)>h.Pos(i+1)
        h.Dur(i)=h.Pos(i+1)-h.Pos(i);
    end
end

if h.Dur(end)>time-h.Pos(end)
     h.Dur(end)=time-h.Pos(end);
end

set(h.dur_table, 'Data', h.Dur)

h.xU=CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).dt,'secs', h.event_duration);

subplot(h.ax5)
plot(linspace(h.M(1).dt,h.N*h.TR,h.T),h.xU)
grid on
axis([0 h.N*h.TR -0.1 1.1])

guidata(h.f,h)

function equal_lengths_button(handle, eventdata)
% Function is used to activate equal lengths control. It changes duration
% of every stimulation event on one similar value (when possible)

h = guidata(handle);

pom=get(h.equal_lengths_button,'handle');
if pom=='on'
    set(h.equal_lengths_ctrl,'Enable','On')
    h=equal_lengths(handle, eventdata);
else
    set(h.equal_lengths_ctrl,'Enable','Off')
end

guidata(h.f,h)

function h=equal_lengths(handle, eventdata)
% Function changes duration of every stimululation event on one similar
% value (when possible).

h = guidata(handle);
time=h.N*h.TR;

h.equal_length=str2num(get(h.equal_lengths_ctrl,'String'));
set(h.equal_lengths_ctrl,'Value',h.equal_length)
pom=length(h.Dur);
h.Dur=h.equal_length*ones(1,pom);

for i=1:h.no_of_events-1
    if h.Pos(i)+h.Dur(i)>h.Pos(i+1)
        h.Dur(i)=h.Pos(i+1)-h.Pos(i);
    end
end

if h.Dur(end)>time-h.Pos(end)
     h.Dur(end)=time-h.Pos(end);
end

set(h.pos_table, 'Data', h.Pos)
set(h.dur_table, 'Data', h.Dur)

h.xU=CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).dt,'secs', h.event_duration);

subplot(h.ax5)
plot(linspace(h.M(1).dt,h.N*h.TR,h.T),h.xU)
grid on % udìlat ještì axis od -0.1 do 1.1
axis([0 h.N*h.TR -0.1 max(h.xU(:))+0.1])

guidata(h.f,h)

function uniform_distribution(handle, eventdata)
% Function for change of stimulation event positions. It makes uniform
% distribution of events during the time.

h = guidata(handle);

time=h.N*h.TR;
h.Pos=linspace(0,time,h.no_of_events+1);
h.Pos=h.Pos(1:end-1);

for i=1:h.no_of_events-1
    if h.Pos(i)+h.Dur(i)>h.Pos(i+1)
        h.Dur(i)=h.Pos(i+1)-h.Pos(i);
    end
end

if h.Dur(end)>time-h.Pos(end)
     h.Dur(end)=time-h.Pos(end);
end

set(h.pos_table, 'Data', h.Pos)
set(h.dur_table, 'Data', h.Dur)

h.xU=CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).dt,'secs', h.event_duration);

subplot(h.ax5)
plot(linspace(h.M(1).dt,h.N*h.TR,h.T),h.xU)
grid on
axis([0 h.N*h.TR -0.1 max(h.xU(:))+0.1])

guidata(h.f,h)


function amplitude(handle, eventdata)
% Function is used to set the amplitude of stimulation vector. Recomended
% value for this amplitude is 1

h = guidata(handle);

h.Ampl = str2double(get(handle,'string'));

if h.Ampl<=0
    h.Ampl=1;
end

set(h.amplitude_ctrl,'string',num2str(h.Ampl));

h.xU=CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).dt,'secs', h.event_duration);

subplot(h.ax5)
plot(linspace(h.M(1).dt,h.N*h.TR,h.T),h.xU)
grid on % udìlat ještì axis od -0.1 do 1.1
axis([0 h.N*h.TR -0.1 max(h.xU(:))+0.1])

guidata(h.f,h)

function save_results(handle, eventdata)
% Function saves rusults of single simulation. The name of saved file
% consists of 'Results_Connectivity_simulator_' prefix, date and time stamp

h = guidata(handle);

str=date;
clk=clock;

res.M=h.M;
res.ROI=h.ROI;
res.SNR=h.SNR;
res.TR=h.TR;
res.dt=h.dt;
res.N=h.N;
res.T=h.T;
res.inputs=h.inputs;
res.no_of_events=h.no_of_events;
res.equal_length=h.equal_length;
res.event_duration=h.event_duration;
res.Pos=h.Pos;
res.Dur=h.Dur;
res.Ampl=h.Ampl;
res.xU=h.xU;
res.xUd=h.xUd;
res.SIM=h.SIM;
res.used_method=h.used_method;
res.DCM=h.DCM;
res.GCM_ret=h.GCM_ret;
res.GC=h.GC;
res.GC2=h.GC2;
res.PR=h.PR;
res.causf=h.causf;
res.causd=h.causd;
res.model_order=h.model_order;
res.aic=h.aic;
res.bic=h.bic;

save(['Results_Connectivity_simulator_',str,'_',num2str(clk(1,4)),'_',num2str(clk(1,5))],'res')

guidata(h.f,h)

function load_results(handle, eventdata)
% Function is used for loading results of connectivity simulation. It is 
% called from enhanced graphical window of Connectivity_simulator

h = guidata(handle);

path = spm_select(1,'mat','Select Results_Connectivity_simulator *.mat file');
load(path)

h.M=res.M;
h.ROI=res.ROI;
h.SNR=res.SNR;
h.TR=res.TR;
h.dt=res.dt;
h.N=res.N;
h.T=res.T;
h.inputs=res.inputs;
h.no_of_events=res.no_of_events;
h.equal_length=res.equal_length;
h.event_duration=res.event_duration;
h.Pos=res.Pos;
h.Dur=res.Dur;
h.Ampl=res.Ampl;
h.xU=res.xU;
h.xUd=res.xUd;
h.SIM=res.SIM;
h.used_method=res.used_method;
h.DCM=res.DCM;
h.GCM_ret=res.GCM_ret;
h.GC=res.GC;
h.GC2=res.GC2;
h.PR=res.PR;
h.causf=res.causf;
h.causd=res.causd;
h.model_order=res.model_order;
h.aic=res.aic;
h.bic=res.bic;

% Basic simulation parameters setting
set(h.connectivity_matrix,'data',h.M(1).A,'ColumnFormat',repmat(cellstr('bank'),1,h.ROI),'ColumnEditable', logical(ones(1,h.ROI)));
set(h.delay_matrix,'data',h.M(1).d,'ColumnEditable', logical(ones(1,h.ROI)));

set(h.ROI_ctrl,'string',num2str(h.ROI));
set(h.SNR_ctrl,'string',num2str(h.SNR));
set(h.TR_ctrl,'string',num2str(h.TR));
set(h.time_step_ctrl,'string',num2str(h.dt));
set(h.scan_ctrl,'string',num2str(h.N));

h.xU = CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).dt,'secs', h.event_duration);
U = repmat(h.xU',h.ROI,1); % Input
h.M(1).U = U;

h.xUd = CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).Ud_dt,'secs', h.event_duration);
h.M(1).Ud = repmat(h.xUd,1,h.ROI);

subplot(h.ax1)
cca_plotcausality(h.M(1).A,[],1);

% Results setting
DCM_A=h.DCM.Ep.A;   % DCM connectivity matrix
DCM_Ap=h.DCM.Pp.A;  % DCM probability matrix
threshold=0.9;      % Threshold for DCM significan values

DCM_A2=DCM_A;       
DCM_A2(DCM_Ap<threshold)=0; % Thresholded matrix

GCM_A_all=h.GC; % GCM connectivity matrix
GCM_A=h.GC2;    % GCM thresholded connectivity matrix
GCM_Ap=h.PR;    % GCM significance matrix

subplot(h.ax2)
plot(h.SIM.y)   % Plot simulated BOLD data
grid on
for i=1:h.ROI
    leg(1,i)={['ROI ',num2str(i)]};
end
legend(leg)
axis([1 h.N min(h.SIM.y(:))-0.1 max(h.SIM.y(:))+0.1])


set(h.f,'Position',[h.fpos h.width h.height])
set(h.panel_2,'Visible','on')
set(h.panel_1,'Position',[0.01 0.01 0.3 0.97])

% ------------ DCM panel ------------
set(h.DCM_A_matrix,'Data',DCM_A)
set(h.DCM_Ap_matrix,'Data',DCM_Ap)
set(h.slider_probability_DCM,'value',0.9);
set(h.slider_edit_probability_DCM,'String',num2str(0.9));
subplot(h.ax3)
cca_plotcausality(DCM_A2,[],1);

% ------------ GCM panel ------------
set(h.GCM_A_matrix,'Data',GCM_A_all)
set(h.GCM_Ap_matrix,'Data',GCM_Ap)
set(h.slider_probability_GCM,'value',0.01);
set(h.slider_edit_probability_GCM,'String',num2str(0.01));
subplot(h.ax4)
cca_plotcausality(GCM_A,[],1);

% Comparison panel
similarity_fcn(h);

if h.used_method=='DCM';
    set(h.used_method_text,'String','Generated by DCM','Visible','on')
    set(h.method_DCM,'Value',get(h.method_DCM,'Max'))
    set(h.method_GCM,'Value',get(h.method_GCM,'Min'))
else
    set(h.used_method_text,'String','Generated by GCM','Visible','on')
    set(h.method_GCM,'Value',get(h.method_GCM,'Max'))
    set(h.method_DCM,'Value',get(h.method_DCM,'Min'))
end

set(h.dur_table, 'Data', h.Dur)
set(h.pos_table, 'Data', h.Pos)
set(h.amplitude_ctrl,'string',num2str(h.Ampl));
set(h.equal_lengths_ctrl,'Value',h.equal_length);
set(h.no_of_events_ctrl,'string',num2str(h.no_of_events));

subplot(h.ax5)
plot(linspace(h.M(1).dt,h.N*h.TR,h.T),h.xU)
grid on
axis([0 h.N*h.TR -0.1 1.1])

set(h.panel_2,'Visible','On');
set(h.panel_3,'Visible','Off');

guidata(h.f,h)

function load_results2(handle, eventdata)
% Function is used for loading results of connectivity simulation. It is 
% called from basic graphical window of Connectivity_simulator

h = guidata(handle);

path = spm_select(1,'mat','Select Results_Connectivity_simulator *.mat file');
load(path)

h.M=res.M;
h.ROI=res.ROI;
h.SNR=res.SNR;
h.TR=res.TR;
h.dt=res.dt;
h.N=res.N;
h.T=res.T;
h.inputs=res.inputs;
h.no_of_events=res.no_of_events;
h.equal_length=res.equal_length;
h.event_duration=res.event_duration;
h.Pos=res.Pos;
h.Dur=res.Dur;
h.Ampl=res.Ampl;
h.xU=res.xU;
h.xUd=res.xUd;
h.SIM=res.SIM;
h.used_method=res.used_method;
h.DCM=res.DCM;
h.GCM_ret=res.GCM_ret;
h.GC=res.GC;
h.GC2=res.GC2;
h.PR=res.PR;
h.causf=res.causf;
h.causd=res.causd;
h.model_order=res.model_order;
h.aic=res.aic;
h.bic=res.bic;

% Basic simulation parameters setting
set(h.connectivity_matrix,'data',h.M(1).A,'ColumnFormat',repmat(cellstr('bank'),1,h.ROI),'ColumnEditable', logical(ones(1,h.ROI)));
set(h.delay_matrix,'data',h.M(1).d,'ColumnEditable', logical(ones(1,h.ROI)));

set(h.ROI_ctrl,'string',num2str(h.ROI));
set(h.SNR_ctrl,'string',num2str(h.SNR));
set(h.TR_ctrl,'string',num2str(h.TR));
set(h.time_step_ctrl,'string',num2str(h.dt));
set(h.scan_ctrl,'string',num2str(h.N));

h.xU = CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).dt,'secs', h.event_duration);
U = repmat(h.xU',h.ROI,1); % Input
h.M(1).U = U;

h.xUd = CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).Ud_dt,'secs', h.event_duration);
h.M(1).Ud = repmat(h.xUd,1,h.ROI);

subplot(h.ax1)
cca_plotcausality(h.M(1).A,[],1);

% Results setting
DCM_A=h.DCM.Ep.A;   % DCM connectivity matrix
DCM_Ap=h.DCM.Pp.A;  % DCM probability matrix
threshold=0.9;      % Threshold for DCM significan values

DCM_A2=DCM_A;       
DCM_A2(DCM_Ap<threshold)=0; % Thresholded matrix

GCM_A_all=h.GC; % GCM connectivity matrix
GCM_A=h.GC2;    % GCM thresholded connectivity matrix
GCM_Ap=h.PR;    % GCM significance matrix


subplot(h.ax2)
plot(h.SIM.y)   % Plot simulated BOLD data
grid on
for i=1:h.ROI
    leg(1,i)={['ROI ',num2str(i)]};
end
legend(leg)
axis([1 h.N min(h.SIM.y(:))-0.1 max(h.SIM.y(:))+0.1])

set(h.f,'Position',[h.fpos h.width h.height])
set(h.panel_2,'Visible','on')
set(h.panel_1,'Position',[0.01 0.01 0.3 0.97])

% ------------------- DCM panel -----------------------
set(h.DCM_A_matrix,'Data',DCM_A)
set(h.DCM_Ap_matrix,'Data',DCM_Ap)

set(h.slider_probability_DCM,'value',0.9);
set(h.slider_edit_probability_DCM,'String',num2str(0.9));

subplot(h.ax3)
cca_plotcausality(DCM_A2,[],1); % Vykreslení modelu podle matice parametrù

% ------------------- GCM panel -----------------------
set(h.GCM_A_matrix,'Data',GCM_A_all)
set(h.GCM_Ap_matrix,'Data',GCM_Ap)

set(h.slider_probability_GCM,'value',0.01);
set(h.slider_edit_probability_GCM,'String',num2str(0.01));

subplot(h.ax4)
cca_plotcausality(GCM_A,[],1);

% Comparison panel
similarity_fcn(h);

if h.used_method=='DCM';
    set(h.used_method_text,'String','Generated by DCM','Visible','on')
    set(h.method_DCM,'Value',get(h.method_DCM,'Max'))
    set(h.method_GCM,'Value',get(h.method_GCM,'Min'))
else %(get(h.method_GCM,'Value') == get(h.method_GCM,'Max'))
    set(h.used_method_text,'String','Generated by GCM','Visible','on')
    set(h.method_GCM,'Value',get(h.method_GCM,'Max'))
    set(h.method_DCM,'Value',get(h.method_DCM,'Min'))
end

set(h.dur_table, 'Data', h.Dur)
set(h.pos_table, 'Data', h.Pos)
set(h.amplitude_ctrl,'string',num2str(h.Ampl));
set(h.equal_lengths_ctrl,'Value',h.equal_length);
set(h.no_of_events_ctrl,'string',num2str(h.no_of_events));

subplot(h.ax5)
plot(linspace(h.M(1).dt,h.N*h.TR,h.T),h.xU)
grid on
axis([0 h.N*h.TR -0.1 1.1])

set(h.panel_2,'Visible','On');
set(h.panel_3,'Visible','Off');

guidata(h.f,h)

function no_of_events_control(handle, eventdata)
% Function is used to control number of stimulation events in stimulation
% vector settings.

h = guidata(handle);

time=h.N*h.TR;

h.no_of_events = str2double(get(handle,'string')); % Number of stimulation events

if h.no_of_events<=0
    h.no_of_events=1;
end

set(h.no_of_events_ctrl,'string',num2str(h.no_of_events));

pom=length(h.Pos);
if pom<h.no_of_events
    h.Pos=[zeros(1,h.no_of_events-pom) h.Pos];
    h.Dur=[zeros(1,h.no_of_events-pom) h.Dur];   
end

if pom>h.no_of_events
    h.Pos=h.Pos(1:h.no_of_events);
    h.Dur=h.Dur(1:h.no_of_events);
end

h.Pos=sort(h.Pos);

set(h.pos_table, 'Data', h.Pos)

for i=1:h.no_of_events-1
    if h.Pos(i)+h.Dur(i)>h.Pos(i+1)
        h.Dur(i)=h.Pos(i+1)-h.Pos(i);
    end
end

if h.Dur(end)>time-h.Pos(end)
     h.Dur(end)=time-h.Pos(end);
end

set(h.dur_table, 'Data', h.Dur)

h.xU=CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).dt,'secs', h.event_duration);

subplot(h.ax5)
plot(linspace(h.M(1).dt,h.N*h.TR,h.T),h.xU)
grid on
axis([0 h.N*h.TR -0.1 max(h.xU(:))+0.1])

guidata(h.f,h)

function Dur_control(handle, eventdata)
% Function is used to change lengths of events during specifying 
% stimulation vector. The unit of every length is [s]

h = guidata(handle);
time=h.N*h.TR;

pos=eventdata.Indices;
new=eventdata.NewData;
h.Dur(pos(1),pos(2))=new;   % Vector of event durations

for i=1:h.no_of_events-1
    if h.Pos(i)+h.Dur(i)>h.Pos(i+1)
        h.Dur(i)=h.Pos(i+1)-h.Pos(i);
    end
end

if h.Dur(end)>time-h.Pos(end)
     h.Dur(end)=time-h.Pos(end);
end

set(h.dur_table, 'Data', h.Dur)

h.xU=CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).dt,'secs', h.event_duration);

subplot(h.ax5)
plot(linspace(h.M(1).dt,h.N*h.TR,h.T),h.xU)
grid on
axis([0 h.N*h.TR -0.1 1.1])

guidata(h.f,h)

function insert_vector(handle, eventdata)
% Fuction enhances the graphical window and enables specifying of
% stimulation vector.

h = guidata(handle);
guidata(h.f,h)

set(h.panel_2,'Visible','Off');
set(h.panel_3,'Visible','On');
set(h.f,'Position',[h.fpos h.width h.height])
set(h.panel_1,'Position',[0.01 0.01 0.3 0.97])

guidata(h.f,h)

function default_settings(handle, eventdata)
% Function sets simulation variables on initial values

h = guidata(handle);

h.M(1).A=[-0.5    0     0; ...  % Connectivity matrix  
           0.9   -0.5   0; ... 
           0.7    0.8  -0.5;];

h.M(1).d=[0   0    0;           % Delay matrix [ms]
          50  0    0;     
          50  50   0]; 

h.ROI=3;    % Number of ROIs
h.SNR=5;    % Signal to noise ratio
h.TR =2;    % Repeat time [s]
h.dt =5;    % Time step [ms]
h.N  =200;  % Number of time points (BOLD signal after downsampling with TR)
h.M(1).dt = h.dt./1000;  % Time step [s]
h.T =h.TR*h.N/h.M(1).dt; % Number of time points in fine sampling

set(h.connectivity_matrix,'data',h.M(1).A,'ColumnFormat',repmat(cellstr('bank'),1,h.ROI),'ColumnEditable', logical(ones(1,h.ROI)));
set(h.delay_matrix,'data',h.M(1).d,'ColumnEditable', logical(ones(1,h.ROI)));

set(h.ROI_ctrl,'string',num2str(h.ROI));
set(h.SNR_ctrl,'string',num2str(h.SNR));
set(h.TR_ctrl,'string',num2str(h.TR));
set(h.time_step_ctrl,'string',num2str(h.dt));
set(h.scan_ctrl,'string',num2str(h.N));
set(h.method_DCM,'Value',get(h.method_DCM,'Max'))
set(h.method_GCM,'Value',get(h.method_GCM,'Min'))

h.M(1).C = zeros(h.ROI,h.ROI); % Control matrix
h.M(1).C(1,1)=1;

% Stimulation vector parameters
event_duration = 500;    % [ms] 
h.Pos = [2 11 34 48.6 60 90 160 200 240 280 300 320 340 360 380]; % Onset positions of stimulation events [s]
h.Dur = [2 10 3 3 20 15 15 15 15 5 5 5 5 5 5 ];                   % Durations of stimulation events [s]
h.Ampl = [1];                   % Amplitude of stimulation vector
h.equal_length = 5;             % Default setting of equal lengths control
h.no_of_events = length(h.Pos); % Number of events
% Stimulation vector
h.xU = CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).dt,'secs', event_duration);
U = repmat(h.xU',h.ROI,1);      % Specifying input
h.M(1).U = U;
% Stimulation vector for DCM analysis
h.xUd = CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).Ud_dt,'secs', event_duration);
h.M(1).Ud = repmat(h.xUd,1,h.ROI);

subplot(h.ax1)
cca_plotcausality(h.M(1).A,[],1); % Model of interactions
set(h.dur_table, 'Data', h.Dur)
set(h.pos_table, 'Data', h.Pos)
set(h.amplitude_ctrl,'string',num2str(h.Ampl));
set(h.equal_lengths_ctrl,'Value',h.equal_length);
set(h.no_of_events_ctrl,'string',num2str(h.no_of_events));

subplot(h.ax5)  % Plot of stimulation vector
plot(linspace(h.M(1).dt,h.N*h.TR,h.T),h.xU)
grid on
axis([0 h.N*h.TR -0.1 1.1])

guidata(h.f,h)

function help(handle, eventdata)
% Function calls the help - informations about program

f2=figure;
ScreenSize = get(0,'ScreenSize');
set(f2,'name','Connectivity_simulator Help','Position',[100 100 ScreenSize(3)-200 ScreenSize(4)-200],'Color','White')

h.HelpText1 = uicontrol('Units','normalized','Position',[0.05 0.87 0.9 0.08],'Style','Text','String',...
                'Connectivity_simulator is tool for simulation of connectivity in fMRI data. It is used to compare two methods for estimation of connectivity – dynamic causal modelling (DCM) and Granger causal modelling (GCM). Connectivity simulator was programed in MATLAB R2009a. It is using toolbox SPM8 and GCCA. It was optimized on screen resolution 1650x1050, Windows 7 Professional, RAM 2GB and processor with clock speed 2,53 GHz.',...
                'BackGroundColor','White','FontSize',12,'HorizontalAlignment','left');
h.HelpText9 = uicontrol('Units','normalized','Position',[0.05 0.81 0.9 0.03],'Style','Text','String',...
                'To provide analysis of connectivity, setup the parameters, press „Run“ button and wait for results.',...
                'BackGroundColor','White','FontSize',12,'HorizontalAlignment','Left');
h.HelpText2 = uicontrol('Units','normalized','Position',[0.05 0.7 0.9 0.08],'Style','Text','String',...
                'The graphical interface is initialized by Connectivity_simulator_DP.m. In this window could be selected basic simulation parameters – signal to noise ratio (SNR), repetition time (TR)[s], time step for fine sampling [ms], no. of scans, no. of ROIs, method for generating data (DCG/GCM), connectivity matrix and delay matrix. „Default settings“ button is used to set all simulation variables on initial conditions. „Load results“ button is used to load saved results of completed simulation. Model of interaction is visualization of connectivity matrix.',...
                'BackGroundColor','White','FontSize',12,'HorizontalAlignment','left');
h.HelpText3 = uicontrol('Units','normalized','Position',[0.05 0.59 0.9 0.08],'Style','Text','String',...
                '„Specify stimulation“ button enhances the graphical interface and enables to specify the simulation vector. Here could the user set time position of every stimulation event (Onset) [s], duration of every stimulation event [s], no. of events and amplitude of stimulation vector. Button „Uniform distribution“ recalculates onsets of events in uniform distribution during stimulation time. Button „Equal lengths“ enables to define same durations of stimulation for every event. Stimulation vector is interactively plotted at bottom of the interface.',...
                'BackGroundColor','White','FontSize',12,'HorizontalAlignment','left');
h.HelpText4 = uicontrol('Units','normalized','Position',[0.05 0.3 0.9 0.25],'Style','Text','String',...
                '„Run” button starts the simulation. During the simulation changes the green “Ready” statement to red “Analyzing” statement. After completing the simulation the graphical interface is enhanced and results are plotted on the right. In the results panel are at the top simulated BOLD data. Middle left are results of DCM analysis: estimated connectivity matrix, probability matrix, selected threshold for significant connections estimation and estimated model of interactions. Middle down are results of GCM analysis: estimated connectivity matrix, significance matrix, selected threshold for significant connections estimation and estimated model of interactions. Middle right is comparison panel. Here are compared structural differences between simulated model and DCM estimated model and between simulated model and GCM estimated model. Also sum of square of difference between simulated and estimated model by the same method is shown here. On the right side is also “Save results” button to save results in *.mat file and “Load results” to load saved results of simulation. Saved results have prefix “Results_Connectivity_simulator_” and time and date stamp.',...
                'BackGroundColor','White','FontSize',12,'HorizontalAlignment','left');
h.HelpText5 = uicontrol('Units','normalized','Position',[0.05 0.2 0.9 0.08],'Style','Text','String',...
                '„More results” button in “DCM analysis” panel opens new figure with additional results of DCM analysis. New figure contains plot of data fitted by DCM estimation combined with stimulation vector, plot of data fitted by DCM estimation combined with simulated BOLD data, bar plot of differences between prior and posterior values of DCM estimated parameters and information about free energy (criterion) and convergence of DCM estimation.',...
                'BackGroundColor','White','FontSize',12,'HorizontalAlignment','left');
h.HelpText6 = uicontrol('Units','normalized','Position',[0.05 0.1 0.9 0.08],'Style','Text','String',...
                '„More results” button in “GCM analysis” panel opens new figure with additional results of GCM analysis. New figure contains plot of simulated BOLD data, GCM estimated model of interactions, causal flow, unit causal density, used model order for GCM analysis and Akaike and Bayesian citerion for recommending the model order.',...
                'BackGroundColor','White','FontSize',12,'HorizontalAlignment','left');
h.HelpText7 = uicontrol('Units','normalized','Position',[0.05 0.05 0.9 0.03],'Style','Text','String',...
                'For more information or bug reports contact the autor: ',...
                'BackGroundColor','White','FontSize',12,'HorizontalAlignment','Right');
h.HelpText8 = uicontrol('Units','normalized','Position',[0.05 0.02 0.9 0.03],'Style','Text','String',...
                'Martin Gajdoš; email: martin.gajdos@centrum.cz',...
                'BackGroundColor','White','FontSize',12,'HorizontalAlignment','Right');
            
function ROI_control(handle, eventdata)
% Function controls number of reions of interest (ROIs)

h = guidata(handle);
h.ROI = round(str2double(get(handle,'string'))); % Number of ROIs

if h.ROI<2
    h.ROI=2;
end

set(h.ROI_ctrl,'string',num2str(h.ROI));
temp_A=h.M(1).A;
temp_d=h.M(1).d;
temp_L=length(temp_A);
h.M(1).A = zeros(h.ROI,h.ROI);
h.M(1).d = zeros(h.ROI,h.ROI);

h.M(1).C = zeros(h.ROI,h.ROI);        % Control matrix
h.M(1).C(1,1)=1;
U = repmat(h.xU',h.ROI,1);            % Specifying input
h.M(1).U = U;
h.M(1).Ud = repmat(h.xUd,1,h.ROI);    % Input for DCM
if (get(h.method_DCM,'Value') == get(h.method_DCM,'Max'))
	h.M(1).A=h.M(1).A-0.5*eye(h.ROI); % -0.5 on diagonal of connectivity matrix
else %(get(h.method_GCM,'Value') == get(h.method_GCM,'Max'))
	h.M(1).A=h.M(1).A+0.5*eye(h.ROI);
end
% Reduction of connectivity and delay matrix (when necessary)
if temp_L>h.ROI
    h.M(1).A=temp_A(1:h.ROI,1:h.ROI);
    h.M(1).d=temp_d(1:h.ROI,1:h.ROI);
end
% Extension of connectivity and delay matrix (when necessary)
if temp_L<h.ROI
    h.M(1).A(1:temp_L,1:temp_L)=temp_A;
    h.M(1).d(1:temp_L,1:temp_L)=temp_d;   
end

set(h.connectivity_matrix,'data',h.M(1).A,'ColumnFormat',repmat(cellstr('bank'),1,h.ROI),'ColumnEditable', logical(ones(1,h.ROI)));
set(h.delay_matrix,'data',h.M(1).d,'ColumnEditable', logical(ones(1,h.ROI)));

subplot(h.ax1)  
cca_plotcausality(h.M(1).A,[],1); % Model of interactions

guidata(h.f,h)

function SNR_control(handle, eventdata)
% Function is used to control signal to noise ratio (SNR) in simulated data

h = guidata(handle);
h.SNR = str2double(get(handle,'string')); % Value of SNR

if h.SNR<0
    h.SNR=0.0001;
end

set(h.SNR_ctrl,'string',num2str(h.SNR));

guidata(h.f,h)

function TR_control(handle, eventdata)
% Function is used to control repetition time in simulated data

h = guidata(handle);
h.TR = str2double(get(handle,'string')); % Value of repetition time

if h.TR<0
    h.TR=0.0001;
end

set(h.TR_ctrl,'string',num2str(h.TR));

time = h.N*h.TR;
h.Pos(h.Pos>time)=0;
h.Pos=sort(h.Pos);
set(h.pos_table, 'Data', h.Pos)

for i=1:h.no_of_events-1
    if h.Pos(i)+h.Dur(i)>h.Pos(i+1)
        h.Dur(i)=h.Pos(i+1)-h.Pos(i);
    end
end

if h.Dur(end)>time-h.Pos(end)
     h.Dur(end)=time-h.Pos(end);
end

set(h.dur_table, 'Data', h.Dur)

set(h.scan_ctrl,'string',num2str(h.N));

h.T  =h.TR*h.N/h.M(1).dt; % Number of elements in fine sampling
h.xU=CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).dt,'secs', h.event_duration); % Stimulation vector


subplot(h.ax5) % Plot of stimulation vector
plot(linspace(h.M(1).dt,h.N*h.TR,h.T),h.xU)
grid on
axis([0 h.N*h.TR -0.1 max(h.xU(:))+0.1])

guidata(h.f,h)

function time_step_control(handle, eventdata)
% Function is used to control time step of fine sampling in simulated data

h = guidata(handle);

h.dt = str2double(get(handle,'string')); % Value of time step

if h.dt<0
    h.dt=0.0001;
end

h.M(1).dt=h.dt/1000;       % Time step [s]
h.T  = h.TR*h.N/h.M(1).dt; % Number of elements in fine sampling
h.xU=CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).dt,'secs', h.event_duration);

set(h.time_step_ctrl,'string',num2str(h.dt));

guidata(h.f,h)

function scan_control(handle, eventdata)
% Function is used to control number of scans in simulated data

h = guidata(handle);
h.N = str2double(get(handle,'string')); % Number of scans

if h.N<2
    h.N=2;
end

time = h.N*h.TR;
h.Pos(h.Pos>time)=0;
h.Pos=sort(h.Pos);
set(h.pos_table, 'Data', h.Pos)

for i=1:h.no_of_events-1
    if h.Pos(i)+h.Dur(i)>h.Pos(i+1)
        h.Dur(i)=h.Pos(i+1)-h.Pos(i);
    end
end

if h.Dur(end)>time-h.Pos(end)
     h.Dur(end)=time-h.Pos(end);
end

set(h.dur_table, 'Data', h.Dur)
set(h.scan_ctrl,'string',num2str(h.N));

h.T  =h.TR*h.N/h.M(1).dt; % Number of elements in fine sampling
h.xU=CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).dt,'secs', h.event_duration); % stimulation vector

subplot(h.ax5)
plot(linspace(h.M(1).dt,h.N*h.TR,h.T),h.xU) % Plot of stimulation vector
grid on 
axis([0 h.N*h.TR -0.1 max(h.xU(:))+0.1])

guidata(h.f,h)

function connectivity_m_fcn(handle, eventdata)
% Function controls values of connectivity matrix for data simulation
h = guidata(handle);

pos=eventdata.Indices;
new=eventdata.NewData;
h.M(1).A(pos(1),pos(2))=new; % Connectivity matrix

subplot(h.ax1)
cca_plotcausality(h.M(1).A,[],1); % Model of interactions

guidata(h.f,h)

function delay_m_fcn(handle, eventdata)
% Function controls values of delay matrix for data simulation
h = guidata(handle);

pos=eventdata.Indices;
new=eventdata.NewData;
h.M(1).d(pos(1),pos(2))=new; % Delay matrix

guidata(h.f,h)

function significance_GCM(handle, eventdata)
% Function is used for displaying significantly estimated structural 
% connections by GCM. Function uses defined p-value for finding significant
% structural connections.

h = guidata(handle);

temp = str2double(get(handle,'string'));
if ~isempty(temp)
    if temp>=1 || temp<=0
        temp=0.01;
    end
    switch handle
        case h.slider_edit_probability_GCM
            PVAL = temp;
            set(h.slider_probability_GCM,'Value',PVAL);
    end
end
PVAL = get(h.slider_probability_GCM,'Value'); % Defined p-value
set(h.slider_probability_GCM,'Value',PVAL);
set(h.slider_edit_probability_GCM,'string',num2str(PVAL));

[PR,q] = cca_findsignificance(h.GCM_ret,PVAL,1); % Finding significant connections, using function from GCCA toolbox
h.PR=PR;                % Significance matrix (GCM)
h.GC = h.GCM_ret.gc;    % Connectivity matrix (GCM)
h.GC2 = h.GC.*h.PR;     % Thresholded connectivity matrix (GCM)

set(h.GCM_Ap_matrix,'Data',h.PR)
set(h.GCM_A_matrix,'Data',h.GC)

subplot(h.ax4)
cca_plotcausality(h.GC2,[],1); % GCM estimated model of interactions

similarity_fcn(h);

guidata(h.f,h)

function threshold_DCM(handle, eventdata)
% Function is used for thresholding DCM results and for displaying them. 
% Thresholding is based on values from DCM probability matrix. These values
% are posterior probabilities of presence of estimated structural 
% connections in analyzed data.

h = guidata(handle);

temp = str2double(get(handle,'string'));
if ~isempty(temp)
    if temp>=1 || temp<=0
        temp=0.01;
    end
    switch handle
        case h.slider_edit_probability_DCM
            threshold = temp;
            set(h.slider_probability_DCM,'Value',threshold);
    end
end
threshold = get(h.slider_probability_DCM,'Value'); % DCM threshold value
set(h.slider_probability_DCM,'Value',threshold);
set(h.slider_edit_probability_DCM,'string',num2str(threshold));

temp_A=h.DCM.Ep.A;  % Connectivity matrix (DCM)
temp_Ap=h.DCM.Pp.A; % Probability matrix

temp_A(temp_Ap<threshold)=0; % Thresholding

subplot(h.ax3)
cca_plotcausality(temp_A,[],1); % Estimated model of interactions (DCM)

similarity_fcn(h)

guidata(h.f,h)

function more_results_DCM(handle, eventdata)
% Function opens new figure and plots aditional informations about DCM
% analysis. New figure contains plot of data fitted by DCM estimation 
% combined with stimulation vector, plot of data fitted by DCM estimation
% combined with simulated BOLD data, bar plot of differences between prior
% and posterior values of DCM estimated parameters and information about
% free energy (criterion) and convergence of DCM estimation.

h = guidata(handle);
figure
subplot(3,2,1:2) 
plot(h.DCM.fit-min(h.DCM.fit(:)))
hold on
plot(linspace(h.M.Ud_dt,h.N,length(h.xUd)),max(max(h.DCM.fit-min(h.DCM.fit(:))))*h.xUd,'k--')   % Stimulation vector
grid on
title('Fitted BOLD data (DCM)')
xlabel('Scans')
ylabel('BOLD')

leg_str1(1,1:1+h.ROI)={['']};
for i=1:h.ROI
    leg_str1(1,i)={['ROI ',num2str(i)]};
end
leg_str1(1,end)={['Stimulation']};
legend(leg_str1,'Location','Eastoutside')

subplot(3,2,3:4)
plot(h.DCM.fit-min(h.DCM.fit(:)))   % Plot of data fitted by DCM estimation 
hold on
% Plot of simulated BOLD data (scaled on fitted data)
plot(((h.SIM.y)./max(h.SIM.y(:))).*(max(max(h.DCM.fit-min(h.DCM.fit(:))))),':')
grid on
title('Fitted and simulated BOLD data')
xlabel('Scans')
ylabel('BOLD')
leg_str(1,1:2*h.ROI)={['']};
leg_str(1,1)={['Fitted']};
leg_str(1,h.ROI+1)={['Simulated']};
legend(leg_str,'Location','Eastoutside')

subplot(3,2,5)
% Differences between prior and posterior values of DCM estimated parameters
bar(h.DCM.parameters) 
title('Conditional minus prior expectation')
xlabel('Parameter');

subplot(3,2,6)
% Free energy
text(0.1,0.6,['Free energy F = ',num2str(h.DCM.F)],'FontSize',16)
% Convergence
if h.DCM.convergence==1
    text(0.1,0.2,['Convergence: OK'],'FontSize',16)
else
    text(0.1,0.2,['Convergence: NOT !!!'],'FontSize',16,'color','red')
end
axis off

guidata(h.f,h)

function more_results_GCM(handle, eventdata)
% Function opens new figure and plots aditional informations about GCM
% analysis. New figure contains plot of simulated BOLD data, GCM estimated
% model of interactions, causal flow, unit causal density, used model
% order for GCM analysis and Akaike and Bayesian citerion for recommending
% the model order.

h = guidata(handle);
figure

FSIZE = 8;
colormap(flipud(bone));

subplot(2,3,1:2);
set(gca,'FontSize',FSIZE);
plot(h.SIM.y); % Plot of simulated BOLD data
xlabel('Scans')
ylabel('BOLD')

grid on
for i=1:h.ROI
    leg(1,i)={['ROI ',num2str(i)]};
end
legend(leg)
axis([1 h.N min(h.SIM.y(:))-0.1 max(h.SIM.y(:))+0.1])
title('Simulated BOLD data')

subplot(233);
% plot granger causalities as a network
cca_plotcausality(h.GC2,[],1);
title(['Estimated model of interactions'])

nvar=h.ROI;

subplot(234);
% plot causal flow  (bar = unweighted, line = weighted)
set(gca,'FontSize',FSIZE);
set(gca,'Box','off');
mval1 = max(abs(h.causf.flow));
mval2 = max(abs(h.causf.wflow));
mval = max([mval1 mval2]);
bar(1:nvar,h.causf.flow,'m');
ylim([-(mval+1) mval+1]);
xlim([0.5 nvar+0.5]);
set(gca,'XTick',[1:nvar]);
set(gca,'XTickLabel',1:nvar);
title('causal flow');
ylabel('out-in');
hold on;
plot(1:nvar,h.causf.wflow);
axis('square');

subplot(235);
% plot unit causal densities  (bar = unweighted, line = weighted)
set(gca,'FontSize',FSIZE);
set(gca,'Box','off');
mval1 = max(abs(h.causd.ucd));
mval2 = max(abs(h.causd.ucdw));
mval = max([mval1 mval2]);
bar(1:nvar,h.causd.ucd,'m');
ylim([-0.25 mval+1]);
xlim([0.5 nvar+0.5]);
set(gca,'XTick',[1:nvar]);
set(gca,'XTickLabel',1:nvar);
title('unit causal density');
hold on;
plot(1:nvar,h.causd.ucdw);
axis('square');

subplot(236);
% Used model order, Akaike and Bayesian criterion 
text(0.1,0.8,['Model order: ',num2str(h.model_order)],'FontSize',16)
text(0.1,0.6,['Akaike information criterion: ',num2str(h.aic)],'FontSize',16)
text(0.1,0.4,['Bayesian information criterion: ',num2str(h.bic)],'FontSize',16)
axis off

guidata(h.f,h)

function method_change(handle, eventdata)
% Function controls change between DCM/GCM method, used for BOLD data 
% simulation.
h = guidata(handle);

if (get(h.method_DCM,'Value') == get(h.method_DCM,'Max'))
    h.M(1).A=h.M(1).A.*~eye(h.ROI)-eye(h.ROI);
else %(get(h.method_GCM,'Value') == get(h.method_GCM,'Max'))
    h.M(1).A=h.M(1).A.*~eye(h.ROI)+0.5*eye(h.ROI);    
end

set(h.connectivity_matrix,'Data',h.M(1).A);

guidata(h.f,h)

function similarity_fcn(h)
% Function is used to estimate structural difference between simulated 
% model of interactions and  model of interaction estimated by DCM and GCM

A = h.M(1).A;                       % Simulated connectivity matrix
sA=logical(A.*(~eye(length(A))));   % Diagonal values are not evaluated

DCM_A  = h.DCM.Ep.A;                         % DCM estimation of connectivity matrix
DCM_Ap = h.DCM.Pp.A;                         % DCM probability matrix
p=get(h.slider_probability_DCM,'value');     % Threshold (DCM)
DCM_A(DCM_Ap<p)=0;                           % Thresholding
DCM_sA=logical(DCM_A.*(~eye(length(DCM_A))));% Diagonal values are not evaluated

GCM_A=h.GC2;;                         % GCM estimation of connectivity matrix
GCM_sA=GCM_A;                         
GCM_sA(find(eye(length(GCM_A))))=0;   % Diagonal values are not evaluated
GCM_sA=logical(GCM_sA);

DCM_error=abs(sA-DCM_sA);
DCM_error=sum(DCM_error(:)); % Number of structural differences between DCM and model
GCM_error=sum(abs(sA-GCM_sA));
GCM_error=sum(GCM_error(:)); % Number of structural differences between GCM and model

set(h.DCM_No,'String',num2str(DCM_error));
set(h.GCM_No,'String',num2str(GCM_error));

if h.used_method == 'DCM'
    % Sum of squares between model of interactions simulated by DCM and estimated by DCM
    DCM_square_error=sum(sum(((A-DCM_A).^2)));
    set(h.sum_of_squares_text,'string','DCM model to estimation:');
    set(h.sum_of_squares_No,'string',num2str(DCM_square_error));
else
    % Sum of squares between model of interactions simulated by GCM and estimated by GCM
    ssA=A.*(~eye(length(A)));
    GCM_ssA=GCM_A;    % Diagonal values are not evaluated
    GCM_ssA(find(eye(length(GCM_A))))=0;
    GCM_square_error=sum(sum(((ssA-GCM_ssA).^2)));
    set(h.sum_of_squares_text,'string','GCM model to estimation:');
    set(h.sum_of_squares_No,'string',num2str(GCM_square_error));
end

guidata(h.f,h)


function simulation(handle, eventdata)
% Function is used for performing of simulation of connectivity. Previously
% selected parameters are used as input for simulation. Then are simulated
% data estimated, then is DCM and GCM analysis performed. At the end are 
% results plotted in graphical window.

h = guidata(handle);

set(h.state,'String','Analyzing','Backgroundcolor','red'); % Ready/Analyzing control
set(h.panel_1_3_state,'Backgroundcolor','red')
pause(0.001);

U = repmat(h.xU',h.ROI,1); % Input of stimulation vector
h.M(1).U = U;

h.xUd = CreateU(h.Pos, h.Dur, h.Ampl, h.N, h.TR, h.M(1).Ud_dt,'secs', h.event_duration); % Stimulation vector for DCM analysis
h.M(1).Ud = repmat(h.xUd,1,h.ROI); % Input of stimulation vector for DCM analysis

if (get(h.method_DCM,'Value') == get(h.method_DCM,'Max'))
    h=create_DCM(h);    % Creating simulated data by DCM
    h.used_method='DCM';
    set(h.used_method_text,'String','Generated by DCM','Visible','on')
else %(get(h.method_GCM,'Value') == get(h.method_GCM,'Max'))
    h=create_GCM(h);    % Creating simulated data by GCM
    h.used_method='GCM';
    set(h.used_method_text,'String','Generated by GCM','Visible','on')
end

h=analyze_DCM(h);   % DCM analysis

h=analyze_GCM(h);   % GCM analysis

subplot(h.ax2)      % Plot of simulated BOLD data
plot(h.SIM.y)
grid on

for i=1:h.ROI
    leg(1,i)={['ROI ',num2str(i)]};
end
legend(leg)
axis([1 h.N min(h.SIM.y(:))-0.1 max(h.SIM.y(:))+0.1])

% Results
DCM_A=h.DCM.Ep.A;  % DCM estimated connectivity matrix
DCM_Ap=h.DCM.Pp.A; % DCM probability matrix
threshold=0.9;     % DCM threshold

DCM_A2=DCM_A;      % Thresholded DCM connectivity matrix
DCM_A2(DCM_Ap<threshold)=0;

GCM_A_all=h.GC;    % GCM estimated connectivity matrix
GCM_A=h.GC2;       % thresholded GCM estimated connectivity matrix
GCM_Ap=h.PR;       % GCM significance matrix

set(h.f,'Position',[h.fpos h.width h.height]) % Graphical interface extension
set(h.panel_2,'Visible','on')                 % Panel with results initialization
set(h.panel_1,'Position',[0.01 0.01 0.3 0.97])

% ------------- Plotting DCM results ----------------------
set(h.DCM_A_matrix,'Data',DCM_A)
set(h.DCM_Ap_matrix,'Data',DCM_Ap)
set(h.slider_probability_DCM,'value',0.9);
set(h.slider_edit_probability_DCM,'String',num2str(0.9));
subplot(h.ax3)
cca_plotcausality(DCM_A2,[],1); % DCM estimated model of interactions

% -------------- Plotting GCM results ---------------------
set(h.GCM_A_matrix,'Data',GCM_A_all)
set(h.GCM_Ap_matrix,'Data',GCM_Ap)
set(h.slider_probability_GCM,'value',0.01);
set(h.slider_edit_probability_GCM,'String',num2str(0.01));
subplot(h.ax4)
cca_plotcausality(GCM_A,[],1); % GCM estimated model of interactions

similarity_fcn(h);

set(h.state,'String','Ready','Backgroundcolor','green'); % Ready/Analyzing control
set(h.panel_1_3_state,'Backgroundcolor','green');
set(h.panel_2,'Visible','On');
set(h.panel_3,'Visible','Off');

guidata(h.f,h)

function h=create_DCM(h)
% Function is used to create simulated BOLD data according DCM method

%---------------------------------------------------------------------
% Inputs 
h.M(1).n    = h.ROI;       % Number of ROIs
h.M(1).N    = h.N;         % Number of time points (BOLD signal after downsampling with TR)
h.M(1).TR   = h.TR;        % Repeat Time [s]
h.M(1).T    = (h.M(1).N*h.M(1).TR)/h.M(1).dt;   % Number of elements in fine sampling
h.M(1).SNR    = h.SNR;     % Signal to noise ratio of the simulated BOLD signal
h.M(1).G      = [ones(1,5)*0.001]';   % Varince of the additive Wiener process noise to the hemodynamic states
                           % Don't go higher than 0.02 
%---------------------------------------------------------------------
% Excitatory connections:
h.M(1).B = zeros(h.ROI,h.ROI,0); % empty so far 
%---------------------------------------------------------------------
% Nonlinear term:
h.M(1).D = zeros(h.ROI,h.ROI,0); % empty so far
%---------------------------------------------------------------------
% coding delay matrix
D = h.M(1).d/(h.M(1).dt*1000) + eye(size(h.M(1).d)); % conversion from ms to time steps; 
Db = logical(D);
Ds = [1:h.ROI];
Ds = Ds(ones(h.ROI,1),:);
Ds = Ds.*Db;
D  = D(:); D(D==0) = [];
Ds = Ds(:); Ds(Ds==0) = [];
for i = 1:length(D)
  Dc{i} = zeros(1,D(i));  % delay sequences
end
h.M(1).Db = Db;
h.M(1).Ds = Ds;
h.M(1).Dc = Dc;
%---------------------------------------------------------------------
%
% hemodynamic model parameters:
%    (a)   (b)  (c)  (d)  (e)
%---------------------------------------------------------------------
parameters = [0.64 0.32 1.00 0.32 0.32];
h.M(1).H = repmat(parameters,h.ROI,1);

% (a) Neuronal efficacy
% (b) Signal decay
% (c) Flow-dependent elimination
% (d) Hemodynamic transit time
% (e) Grubb's exponent
%--------------------------------------------------------------------------
h.SIM = DCM_generateDE(h.M); % Creating simulated data

function h=create_GCM(h)
% Function is used to create simulated BOLD data according GCM method

%---------------------------------------------------------------------
% Inputs
SNR = h.SNR;     % Signal to noise ratio
TR  = h.TR;      % Repetition time
N   = h.N;       % Number of time points (BOLD signal after downsampling with TR)
dt  = h.M(1).dt; % Time step [s]

delka_signalu=N*TR;             % Length of signal [s]
fvz=1/dt;                       % Fine sampling frequency
cas=1/fvz:1/fvz:delka_signalu;  % Time axis

stimulace = h.xU';  % Stimulation vector
A= h.M(1).A;        % Simulated connectivity matrix - diagonal elements shoulde be in range 0 - 1, ideal near 0.5
d= h.M(1).d;        % Simulated delay matrix
%--------------------------------------------------------------------------
% Creating signal on neuronal level
signal = (generate_x(stimulace,dt))';
ROI=h.ROI;
n  = ROI;
AA = repmat(A*(~eye(n)),1,max(d(:))/5); % take off-diagonal elements of A
Am = A.*(~eye(n));  
AA = diag(diag(A));   % take diagonal elements of A
for i = 1:max(d(:))/5
    Ad = Am;
    if any(any(d == i*5))
        Ad = Am.*(d == i*5);
    else
        Ad = zeros(n);
    end
    AA = [AA,Ad];
end
% i - max delay index
% main data generation:
%==========================================================================
Y     = zeros(ROI,i+1);  % Matrix of simulated neuronal signal
Y_N   = Y(:);
Y     = zeros((ROI*(i+1)),length(signal));
Y(:,1)= Y_N;
C     = eye(ROI)/5;      
U     = signal;          % Input vector to 1st node
U     = [U;zeros(ROI-1,length(U))];
Dlag  = [eye((i)*ROI),zeros((i)*ROI,ROI)];    % transformation matrix

for t = 2:length(signal) % Simulating of connectivity in data using GCM method
   Y(1:ROI,t)     = AA*Y(:,t-1) + C*U(:,t-1);
   Y(ROI+1:end,t) = Dlag*Y(:,t-1);
end
y = Y(1:ROI,:);
%==========================================================================
% Using Hemodynamic model to transform neuronal signal to BOLD signal
parameters = [0.64 0.32 1.00 0.32 0.32]; % Hemodynamic parameters
H = repmat(parameters,h.ROI,1);

[BOLD states] = generate_y(y',H,dt);% Transformation of neuronal signal to BOLD signal
y=BOLD;
%--------------------------------------------------------------------------
% Downsampling
yd=y(1:(TR/dt):end,:);
%--------------------------------------------------------------------------
% Adding noise
z=randn(size(yd));
r=diag(std(yd)/SNR);
yd=yd + z*r;

h.SIM.y=yd; % Saves downsampled simulated BOLD
h.SIM.y0=y; % Saves simulated BOLD

function h=analyze_DCM(h)
% Function is used to analyze data by DCM

DCM.a=ones(h.ROI,h.ROI);   % Possible endogenous connections
DCM.b=zeros(h.ROI,h.ROI);  % Possible bilinear modulations - zeros
DCM.c=zeros(h.ROI,h.ROI);  % Possible exogenous connections
DCM.c(1,1)=1; 
DCM.U.u=h.M(1).Ud;         % Exogenous inputs - stimulation vector
DCM.U.dt=h.M(1).Ud_dt;     % Time step of stimulation vector

DCM.Y.y=h.SIM.y;           % Downsampled BOLD with noise
DCM.Y.Q = spm_Ce(ones(1,h.ROI)*h.N);
DCM.Y.X0=[];               % Confounds - empty
DCM.Y.dt=h.TR;             % TR sampling of BOLD
DCM.n=h.ROI;               % Number of regions
DCM.v=h.N;                 % Number of scans

DCM.options.nograph=1;     % Turning off SPM interactive screen

[h.DCM] = spm_dcm_estimate_MG(DCM);

function h=analyze_GCM(h)
% Function is used to analyze data by GCM. The procedure is edit m-file
% from cca_demo.m from GCCA toolbox

% Input - simulated data
X=h.SIM.y';

N       =   length(X);  % number of observations
PVAL    =   0.01;       % probability threshold for Granger causality significance
NLAGS   =   -1;         % if -1, best model order is assessed automatically

nvar = size(X,1);

% detrend and demean data
disp('detrending and demeaning data');
X = cca_detrend(X);
X = cca_rm_temporalmean(X);

% check covariance stationarity
disp('checking for covariance stationarity ...');
uroot = cca_check_cov_stat(X,10);
inx = find(uroot);
if sum(uroot) == 0,
    disp('OK, data is covariance stationary by ADF');
else
    disp('WARNING, data is NOT covariance stationary by ADF');
    disp(['unit roots found in variables: ',num2str(inx')]);
end

% check covariance stationarity again using KPSS test
[kh,kpss] = cca_kpss(X);
inx = find(kh==0);
if isempty(inx),
    disp('OK, data is covariance stationary by KPSS');
else
    disp('WARNING, data is NOT covariance stationary by KPSS');
    disp(['unit roots found in variables: ',num2str(inx)]);
end
    
% find best model order
if NLAGS == -1,
    disp('finding best model order ...');
    [bic,aic] = cca_find_model_order(X,2,12);
    disp(['best model order by Bayesian Information Criterion = ',num2str(bic)]);
    disp(['best model order by Aikaike Information Criterion = ',num2str(aic)]);
    NLAGS = aic;
end
%--------------------------------------------------------------------------
% analyze time-domain granger
% find time-domain conditional Granger causalities [THIS IS THE KEY FUNCTION]
disp('finding conditional Granger causalities ...');
ret = cca_granger_regress(X,NLAGS,1);   % STATFLAG = 1 i.e. compute stats

% check that residuals are white
dwthresh = 0.05/nvar;    % critical threshold, Bonferroni corrected
waut = zeros(1,nvar);
for ii=1:nvar,
    if ret.waut<dwthresh,
        waut(ii)=1;
    end
end
inx = find(waut==1);
if isempty(inx),
    disp('All residuals are white by corrected Durbin-Watson test');
else
    disp(['WARNING, autocorrelated residuals in variables: ',num2str(inx)]);
end

% check model consistency, ie. proportion of correlation structure of the
% data accounted for by the MVAR model
if ret.cons>=80,
    disp(['Model consistency is OK (>80%), value=',num2str(ret.cons)]);
else
    disp(['Model consistency is <80%, value=',num2str(ret.cons)]);
end
        
% analyze adjusted r-square to check that model accounts for the data (2nd check)
rss = ret.rss_adj;
inx = find(rss<0.3);
if isempty(inx)
    disp(['Adjusted r-square is OK: >0.3 of variance is accounted for by model, val=',num2str(mean(rss))]);
else
    disp(['WARNING, low (<0.3) adjusted r-square values for variables: ',num2str(inx)]);
    disp(['corresponding values are ',num2str(rss(inx))]);
    disp('try a different model order');
end

% find significant Granger causality interactions (Bonferonni correction)
[PR,q] = cca_findsignificance(ret,PVAL,1);
disp(['testing significance at P < ',num2str(PVAL), ', corrected P-val = ',num2str(q)]);

% extract the significant causal interactions only
GC = ret.gc;
GC2 = GC.*PR;
h.GCM_ret=ret;
% calculate causal connectivity statistics
disp('calculating causal connectivity statistics');
causd = cca_causaldensity(GC,PR);
causf = cca_causalflow(GC,PR);

disp(['time-domain causal density = ',num2str(causd.cd)]);
disp(['time-domain causal density (weighted) = ',num2str(causd.cdw)]);

h.GC=GC;
h.GC2=GC2;
h.PR=PR;
h.causf=causf;
h.causd=causd;
h.model_order=NLAGS;
h.aic=aic;
h.bic=bic;


function [U] = CreateU (Pos, Dur, Ampl, nscan, TR, dt, unit_flag, event_duration);
% Create fine-sampled activation/stimulation time-series
% function [U] = CreateU (Pos, Dur, Ampl, nscan, TR, dt, unit_flag, event_duration);
%
%   Pos    =  vector of positions (begins of each trial/condition/event in scans or secs)
%   Dur    =  vector of durations  (durations of each trial in scans or secs)
%   Ampl   =  vector of amplitudes (amplitudes of each trial)
%   dt     =  fine resolution step (in seconds);
%   nscan  =  number of scans
%   TR     =  scan repetition time in seconds
%   unit_flag = definition of units for input parameters (Pos and Dur)
%                'secs' or 'seconds' for seconds
%                'scans' for scans
%   event_duration = duration for events defined with zero Dur (in miliseconds)
%                     if zero then duration = TR
%
% *************************************************************************
%  Written by Michal Mikl
%  last update 2011-03-08
%  version 1.0
% *************************************************************************

% main routine
%**************************************************************************

FSfactor = TR/dt;    % fine sampling factor
nf = FSfactor*nscan; % number of time points using fine sampling

nPos = length(Pos);
nAmpl = length(Ampl);
nDur = length(Dur);

if event_duration == 0
    event_duration = tr*1000;
end
evdur = round(event_duration/dt/1000); %event duration for zero-legnth events in fine samples

% ensure the same length of Pos, Dur and Ampl vectors
if nAmpl > nPos
    Ampl=Ampl(1:nPos);
end
if nDur > nPos
    Dur=Dur(1:nPos);
end
if nAmpl < nPos
    Ampl = Ampl*ones(1,nPos);
end;
if nDur < nPos
    Dur = Dur*ones(1,nPos);
end;

% convert seconds to scans
if (strcmp(unit_flag,'seconds')) || (strcmp(unit_flag,'secs'))
    Pos = Pos / TR;
    Dur = Dur / TR;
end


fsPos = ceil(Pos * FSfactor);        % fine sampled positions
fsDur = round(Dur * FSfactor);         % fine sampled durations
fsAmpl = Ampl;                         % fine sampled amplitudes

for i=1:nPos
    if fsDur(i)==0
        fsDur(i)=evdur;
        fsAmpl(i)=1;
    end;
end;

tmpseries=zeros(nf,1);        % raw (zero-filled) fine sampled time series

% fill in timeseries positions with amplitudes = build experimental timecourse
for i=0:nf-1
    k=find(fsPos==i);
    if (k>0)
        for l=1:length(k)
        for j=1:fsDur(k(l))
            tmpseries(i+j,1)=tmpseries(i+j,1)+fsAmpl(k(l));
        end
        end
    end
end

U = tmpseries;

function [X] = generate_x(U,dt)
% Function generates time-course of neural activation (x1 from DCM states)
% X - time-course of neural activation
% U - time-course of stimulus function (fine sampling accordint to dt, size (1,n))
% dt - step of fine sampling in seconds (e.g. 0.005 = 5 ms)
% *************************************************************************
%  Written by Michal Mikl
% *************************************************************************
           
T = size(U,2);
X        = zeros(T,1);
x = 0;

for t = 1:T,

    
        f = -x + 0.5*U(1,t);
        dx       = f * dt;
       
        x        = x + dx;
        X(t,1)   = x;            
end

function [Y X] = generate_y(X1,H,dt)
% function Y = generate_y(X,H,dt)
%
% generate bold signal from neural time-course using hemodynamic model (the same as for DCM)
% Y - BOLD signal in fine sampling (same as for neural signal)
% X - time-course of hemodynamic states
% X1 - time-course of neural signal - size T x n, 
%                         T = number of time points,
%                         n = number of regions
% H - parameters of hemodynamic model - size n x 5
% dt - time step
% *************************************************************************
%  Written by Michal Mikl
% *************************************************************************

% number of time points and regions
[T, n] = size(X1);

% initialization of hemodynamic states
x = repmat([0 0 1 1 1]',1,n);              % x(1,:) = neuronal signal
                                                % x(2,:) = vascular signal
                                                % x(3,:) = rCBF 
                                                % x(4,:) = venous volume 
                                                % x(5,:) = deoyxHb   
X        = zeros(T,n*5);
Y        = zeros(T,n);

for t = 1:T,
       
        % hemodynamic state equation
        ind = ([1:n]*5)+1-5;
        x(ind)=X1(t,:);
        f = dxdt(x,H);
        % integration by euler scheme:
        %------------------------------------------------------------------
        dx       = f * dt;
        
        x        = x + reshape(dx ,5, n);
        y        = spm_gx_fmri2(x,[],[]);    % observation equation: 
        X(t,:)   = x(:)';            
        Y(t,:)   = y(:)';
end

function y = dxdt(x,H)
% Calculate changes in hemodynamic states according to neural states and
% hemodynamic parameters
% *************************************************************************
%  Written by Michal Mikl
% *************************************************************************
x     = x';
y    = x;

% Hemodynamic motion
%==========================================================================

% hemodynamic parameters
%--------------------------------------------------------------------------
%   H(1) - signal decay                                   d(ds/dt)/ds)
%   H(2) - autoregulation                                 d(ds/dt)/df)
%   H(3) - transit time                                   (t0)
%   H(4) - exponent for Fout(v)                           (alpha)
%   H(5) - resting oxygen extraction                      (E0)
%   H(6) - ratio of intra- to extra-vascular components   (epsilon)
%          of the gradient echo signal
%--------------------------------------------------------------------------
%H        = [0.65 0.41 2.00 0.32 0.34];

% signal decay
%--------------------------------------------------------------------------
sd       = H(:,1);

% transit time
%--------------------------------------------------------------------------
tt       = H(:,3);

% Fout = f(v) - outflow
%--------------------------------------------------------------------------
fv       = x(:,4).^(1./H(:,4));
% e = f(f) - oxygen extraction
%--------------------------------------------------------------------------
ff       = (1 - (1 - H(:,5)).^(1./x(:,3)))./H(:,5);

% implement differential state equation y = dx/dt (hemodynamic)
%--------------------------------------------------------------------------
y(:,2)   = x(:,1) - sd.*x(:,2) - H(:,2).*(x(:,3) - 1);
y(:,3)   = x(:,2);
y(:,4)   = (x(:,3) - fv)./(tt);
y(:,5)   = (ff.*x(:,3) - fv.*x(:,5)./x(:,4))./(tt);

y        = y';
y        = y(:);

function [SIM] = DCM_generateDE(M)
% Generate BOLD signal according DCM model
% *************************************************************************
%  Written by Michal Mikl and Martin Havlíèek
% *************************************************************************

% initialization of hemodynamic states
x = repmat([0 0 1 1 1]',1,M(1).n);              % x(1,:) = neuronal signal
                                                % x(2,:) = vascular signal
                                                % x(3,:) = rCBF 
                                                % x(4,:) = venous volume 
                                                % x(5,:) = deoyxHb   
% Inputs
DCM(1).P.H  = M(1).H;
DCM(1).P.A  = M(1).A;
DCM(1).P.C  = M(1).C;%/sqrt(0.001);
DCM(1).P.B  = M(1).B; 
DCM(1).P.D  = M(1).D;
DCM(1).U    = M(1).U;
DCM(1).P.Ds = M(1).Ds;
DCM(1).P.Db = M(1).Db;
DCM(1).P.Dc = M(1).Dc;

EXPm   = kron(speye(M(1).n),sparse([ones(5),2*ones(5,1);zeros(1,5+1)]));
xt     = repmat([zeros(1,5) 1]',M(1).n,1);
J      = repmat([1 0 0 0 0;
                 1 1 1 0 0;
                 0 1 0 0 0;
                 0 0 1 1 0;
                 0 0 1 1 1;],1,M(1).n);
             
ww    =  zeros(5*M(1).n,1);

X        = zeros(M(1).T,M(1).n*5);
Y        = zeros(M(1).T,M(1).n);

for t = 1:M(1).T,
       if t==400
         disp('k')  
       end   
    
        [f DCM(1).P] = spm_fx_fmriD(x,DCM(1).U(:,t),DCM(1).P); % hemodynamic state equation
        % integration by euler scheme:
        %------------------------------------------------------------------
        dx       = f * M(1).dt;
        x        = x + reshape(dx + ww,5,M(1).n);
        y        = spm_gx_fmri2(x,[],DCM(1).P);    % observation equation:
        X(t,:)   = x(:)';            
        Y(t,:)   = y(:)';
end

% downsample the BOLD signal
Yd      =   Y(1:(1/(M(1).T*M(1).dt))*M(1).T*M(1).TR:end,:);

% Add noise to the BOLD signal
%--------------------------------------------------------------------------
z       = randn(M(1).N,M(1).n);
% standard deviation of additive noise, for all areas
r       = diag(std(Yd)/M(1).SNR);
Ydn     = Yd + z*r;


% save simulated data:
%-----------------------------------------------------------
SIM(1).x  = X;   % simulated hemodynamic states 
SIM(1).y  = Ydn; % similated hemodynamic response - downsampled and with additive noise 
SIM(1).y0 = Y;   % similated hemodynamic response - high resolution - no noise

function [DCM] = spm_dcm_estimate_MG(P)
% Estimates parameters of a DCM (bilinear or nonlinear) for fMRI data
% Modificated function spm_dcm_estimate.m from SPM8 toolbox
% FORMAT [DCM] = spm_dcm_estimate_MG(DCM)
%   DCM - DCM structure or its filename
%
% Expects
%--------------------------------------------------------------------------
% DCM.a                              % switch on endogenous connections
% DCM.b                              % switch on bilinear modulations
% DCM.c                              % switch on exogenous connections
% DCM.d                              % switch on nonlinear modulations
% DCM.U                              % exogenous inputs
% DCM.Y                              % responses
% DCM.Y.X0                           % confounds
% DCM.n                              % number of regions
% DCM.v                              % number of scans
%
% Options
%--------------------------------------------------------------------------
% DCM.options.two_state              % two regional populations (E and I)
% DCM.options.stochastic             % fluctuations on hidden states
% DCM.options.nonlinear              % interactions among hidden states
% DCM.options.nograph                % graphical display
% DCM.options.centre                 % mean-centre inputs
% DCM.options.P                      % Starting estimates for parameters
%
% Evaluates:
%--------------------------------------------------------------------------
% DCM.M                              % Model structure
% DCM.Ep                             % Condition means (parameter structure)
% DCM.Cp                             % Conditional covariances
% DCM.Vp                             % Conditional variances
% DCM.Pp                             % Conditional probabilities
% DCM.H1                             % 1st order hemodynamic kernels
% DCM.H2                             % 2nd order hemodynamic kernels
% DCM.K1                             % 1st order neuronal kernels
% DCM.K2                             % 2nd order neuronal kernels
% DCM.R                              % residuals
% DCM.y                              % predicted data
% DCM.T                              % Threshold for Posterior inference
% DCM.Ce                             % Error variance for each region
% DCM.F                              % Free-energy bound on log evidence
% DCM.ID                             % Data ID
% DCM.AIC                            % Akaike Information criterion
% DCM.BIC                            % Bayesian Information criterion
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
 
% Karl Friston
% $Id: spm_dcm_estimate.m 4185 2011-02-01 18:46:18Z guillaume $
 
 
%-Load DCM structure
%--------------------------------------------------------------------------
if ~nargin
 
    %-Display model details
    %----------------------------------------------------------------------
    Finter = spm_figure('GetWin','Interactive');
    set(Finter,'name','Dynamic Causal Modelling')
 
    %-Get DCM
    %----------------------------------------------------------------------
    [P, sts] = spm_select(1,'^DCM.*\.mat$','select DCM_???.mat');
    if ~sts, DCM = []; return; end
    spm('Pointer','Watch')
    spm('FigName','Estimation in progress');
 
end
if isstruct(P)
    DCM = P;
    P   = ['DCM-' date '.mat'];
else
    load(P)
end

% check options
%==========================================================================
try, DCM.options.two_state;  catch, DCM.options.two_state  = 0;     end
try, DCM.options.stochastic; catch, DCM.options.stochastic = 0;     end
try, DCM.options.nonlinear;  catch, DCM.options.nonlinear  = 0;     end
try, DCM.options.centre;     catch, DCM.options.centre     = 1;     end

try, M.nograph = DCM.options.nograph; catch, M.nograph = spm('CmdLine');end
try, M.P       = DCM.options.P ;end
 
% unpack
%--------------------------------------------------------------------------
U  = DCM.U;                             % exogenous inputs
Y  = DCM.Y;                             % responses
n  = DCM.n;                             % number of regions
v  = DCM.v;                             % number of scans
 
% detrend outputs (and inputs)  
%--------------------------------------------------------------------------
Y.y = spm_detrend(Y.y);
if DCM.options.centre
    U.u = spm_detrend(U.u);
end

% check scaling of Y (enforcing a maximum change of 4%
%--------------------------------------------------------------------------
scale   = max(max((Y.y))) - min(min((Y.y)));
scale   = 4/max(scale,4);
Y.y     = Y.y*scale;
Y.scale = scale;

% check confounds (add constant if necessary)
%--------------------------------------------------------------------------
if ~size(Y.X0,2), Y.X0 = ones(v,1); end

% fMRI slice time sampling
%--------------------------------------------------------------------------
try, M.delays = DCM.delays; end
try, M.TE     = DCM.TE;     end

% create priors
%==========================================================================

% check DCM.d (for nonlinear DCMs)
%--------------------------------------------------------------------------
try
    DCM.d;
catch
    DCM.d = zeros(n,n,0);
    DCM.options.nonlinear = 0;
end
 
% specify parameters for spm_int_D (ensuring updates every second or so)
%--------------------------------------------------------------------------
if DCM.options.nonlinear
    M.IS     = 'spm_int_D';
    M.nsteps = round(max(Y.dt,1));
    M.states = 1:n;
else
    M.IS     = 'spm_int';
end
 
% check DCM.c (for endogenous DCMs)
%--------------------------------------------------------------------------
if isempty(DCM.c)
    DCM.c  = zeros(n,1);
    DCM.b  = zeros(n,n,1);
    U.u    = zeros(v,1);
    U.name = {'null'};
    DCM.options.endogenous = 1;
end
if ~any(spm_vec(DCM.c))
    DCM.options.endogenous = 1;
else
    DCM.options.endogenous = 0;
end
if  DCM.options.endogenous
    DCM.options.two_state  = 0;
    DCM.options.stochastic = 1;
end

% priors (and initial states)
%--------------------------------------------------------------------------
[pE,pC,x] = spm_dcm_fmri_priors(DCM.a,DCM.b,DCM.c,DCM.d,DCM.options);

% complete model specification
%--------------------------------------------------------------------------
M.f   = 'spm_fx_fmri';
M.g   = 'spm_gx_fmri';
M.x   = x;
M.pE  = pE;
M.pC  = pC;
M.m   = size(U.u,2);
M.n   = size(x(:),1);
M.l   = size(x,1);
M.N   = 32;
M.dt  = 16/M.N;
M.ns  = v;
 
 
% nonlinear system identification (nlsi)
%==========================================================================
 
% nonlinear system identification (Variational EM) - deterministic DCM
%--------------------------------------------------------------------------
if ~DCM.options.endogenous
    [Ep,Cp,Eh,F,fit,parameters,convergence] = spm_nlsi_GN_MG(M,U,Y);
else
    Ep = pE;
    Eh = 3;
    F  = 0;
end

% proceed to stochastic (initialising with deterministic estimates)
%--------------------------------------------------------------------------
if DCM.options.stochastic
    
    
    % specify bilinear approximation scheme for spm_DEM_eval
    %----------------------------------------------------------------------
    if size(pE.A,3), M.E.linear = 1; end   % linear model
    if size(pE.B,3), M.E.linear = 2; end   % bilinear model
    if size(pE.D,3), M.E.linear = 3; end   % nonlinear models
 
    % Decimate U.u from micro-time
    % ---------------------------------------------------------------------
    u       = U.u;
    y       = Y.y;
    Dy      = spm_dctmtx(size(y,1),size(y,1));
    Du      = spm_dctmtx(size(u,1),size(y,1));
    Dy      = Dy*sqrt(size(y,1)/size(u,1));
    u       = Dy*(Du'*u);
 
    % DEM Structure: place model, data, input and confounds in DEM
    % ---------------------------------------------------------------------
    DEM.M   = M;
    DEM.Y   = y';
    DEM.U   = u';
    DEM.X   = Y.X0';
 
    % set inversion parameters
    % ---------------------------------------------------------------------
    DEM.M(1).E.form = 'Gaussian';         % form of random fluctuations
    DEM.M(1).E.s    = 1/2;                % smoothness of fluctuations
    DEM.M(1).E.d    = 2;                  % embedding dimension 
    DEM.M(1).E.n    = 6;                  % embedding dimension
    DEM.M(1).E.nN   = 32;                 % maximum number of iterations

 
    % adjust M.f (DEM works in time bins not seconds) and initialize M.P
    % ---------------------------------------------------------------------
    DEM.M(1).delays = M.delays/Y.dt;
    DEM.M(1).f  = inline([M.f '(x,v,P)*' num2str(Y.dt)],'x','v','P');
    DEM.M(1).P  = Ep;
    
    
    % Specify hyper-priors on (log-precision of) observation noise
    % ---------------------------------------------------------------------
    DEM.M(1).Q  = spm_Ce(ones(1,n));
    DEM.M(1).hE = Eh + 1;                 % prior expectation
    DEM.M(1).hC = 1/16;                   % prior covariance
    
    
    % Specify hyper-priors, allowing hidden causes to fluctuate
    % ---------------------------------------------------------------------
    DEM.M(1).xP = 128;                    % s.d. ~ 10%      (hidden-state)
    DEM.M(1).W  = exp(9);                 % fixed precision (hidden-motion)
    DEM.M(2).V  = exp(4);                 % fixed precision (hidden-cause)
        
    % for endgenous DCMs, allow neuronal hidden states to fluctuate
    % ---------------------------------------------------------------------
    if DCM.options.endogenous
        W           = exp(spm_vec(sparse(1:n,1,(6 - 16),n,5) + 16));
        DEM.M(1).xP = exp(6);             % fixed precision (hidden-state)
        DEM.M(1).W  = diag(W);            % fixed precision (hidden-motion)
        DEM.M(2).V  = exp(16);            % fixed precision (hidden-cause)
    end


    % Generalised filtering (under the Laplace assumption)
    % =====================================================================
    if DCM.options.stochastic == 1
            DEM  = spm_LAP(DEM);          % no mean-field assumption
    else
            DEM  = spm_DEM(DEM);          % mean-field assumption (DEM)
    end
    
    % Helpful suggestion
    %----------------------------------------------------------------------
    if length(DEM.F) < 4
        disp('This inversion may have been unstable;')
        disp('try reducing DCM.options.s to 1/4')
    end
    
    
    % Save DEM estimates
    %----------------------------------------------------------------------
    DCM.qU = DEM.qU;
    DCM.qP = DEM.qP;
    DCM.qH = DEM.qH;
    DCM.Fd = F;
    
    % unpack results
    % ---------------------------------------------------------------------
    F      = DEM.F(end);
    Ep     = DEM.qP.P{1};
    Cp     = DEM.qP.C;
 
    % predicted responses (y) and residuals (R) from DEM
    %----------------------------------------------------------------------
    y      = DEM.qU.v{1}';
    R      = DEM.qU.z{1}';
    R      = R - Y.X0*spm_inv(Y.X0'*Y.X0)*(Y.X0'*R);
    Ce     = exp(-DEM.qH.h{1});
 
else
 
    % predicted responses (y) and residuals (R) from EM
    %----------------------------------------------------------------------
    y      = feval(M.IS,Ep,M,U);
    R      = Y.y - y;
    R      = R - Y.X0*spm_inv(Y.X0'*Y.X0)*(Y.X0'*R);
    Ce     = exp(-Eh);
    
end
 
 
% Bilinear representation and first-order hemodynamic kernel
%--------------------------------------------------------------------------
[M0,M1,L1,L2] = spm_bireduce(M,Ep);
[H0,H1] = spm_kernels(M0,M1,L1,L2,M.N,M.dt);
 
% and neuronal kernels
%--------------------------------------------------------------------------
L       = sparse(1:n,[1:n] + 1,1,n,length(M0));
[K0,K1] = spm_kernels(M0,M1,L,M.N,M.dt);
 
 
% Bayesian inference and variance {threshold T = 0}
%--------------------------------------------------------------------------
T       = 0;
sw      = warning('off','SPM:negativeVariance');
Pp      = spm_unvec(1 - spm_Ncdf(T,abs(spm_vec(Ep)),diag(Cp)),Ep);
Vp      = spm_unvec(diag(Cp),Ep);
warning(sw);
 
try, M = rmfield(M,'nograph'); end

% Store parameter estimates
%--------------------------------------------------------------------------
DCM.M   = M;
DCM.Y   = Y;
DCM.U   = U;
DCM.Ce  = Ce;
DCM.Ep  = Ep;
DCM.Cp  = Cp;
DCM.Pp  = Pp;
DCM.Vp  = Vp;
DCM.H1  = H1;
DCM.K1  = K1;
DCM.R   = R;
DCM.y   = y;
DCM.T   = T;
 
 
% Data ID and log-evidence
%==========================================================================
if isfield(M,'FS')
    try
        ID = spm_data_id(feval(M.FS,Y.y,M));
    catch
        ID = spm_data_id(feval(M.FS,Y.y));
    end
else
    ID     = spm_data_id(Y.y);
end
 
% Save approximations to model evidence: negative free energy, AIC, BIC
%--------------------------------------------------------------------------
evidence   = spm_dcm_evidence(DCM);
DCM.F      = F;
DCM.ID     = ID;
DCM.AIC    = evidence.aic_overall;
DCM.BIC    = evidence.bic_overall;
DCM.fit    = fit;
DCM.parameters = parameters;
DCM.convergence = convergence;
 
%-Save DCM
%--------------------------------------------------------------------------
if ~isstruct(P)
    if spm_check_version('matlab','7') >= 0
        save(P,'-V6','DCM','F','Ep','Cp');
    else
        save(P,'DCM','F','Ep','Cp');
    end
end

if ~nargin
    spm('Pointer','Arrow');
    spm('FigName','Done');
end

function [Ep,Cp,Eh,F,fit,parameters,convergence] = spm_nlsi_GN_MG(M,U,Y)
% Bayesian inversion of a nonlinear model using a Gauss-Newton/EM algorithm
% Modificated function spm_nlsi_GN.m from SPM8 toolbox
% FORMAT [Ep,Cp,Eh,F] = spm_nlsi_GN(M,U,Y)
%
% Dynamical MIMO models
%__________________________________________________________________________
%
% M.IS - function name f(P,M,U) - generative model
%        This function specifies the nonlinear model: 
%        y = Y.y = IS(P,M,U) + X0*P0 + e
%        were e ~ N(0,C).  For dynamic systems this would be an integration
%        scheme (e.g. spm_int). spm_int expects the following:
%
%     M.f  - f(x,u,P,M)
%     M.g  - g(x,u,P,M)
%       x  - state variables
%       u  - inputs or causes
%       P  - free parameters
%       M  - fixed functional forms and parameters in M
%
% M.FS - function name f(y,M)   - feature selection
%        This [optional] function performs feature selection assuming the
%        generalized model y = FS(y,M) = FS(IS(P,M,U),M) + X0*P0 + e
%
% M.P  - starting estimates for model parameters [optional]
%
% M.pE - prior expectation  - E{P}   of model parameters
% M.pC - prior covariance   - Cov{P} of model parameters
%
% M.hE - prior expectation  - E{h}   of log-precision parameters
% M.hC - prior covariance   - Cov{h} of log-precision parameters
%
% U.u  - inputs
% U.dt - sampling interval
%
% Y.y  - outputs
% Y.dt - sampling interval for outputs
% Y.X0 - Confounds or null space      (over size(y,1) bins or all vec(y))
% Y.Q  - q error precision components (over size(y,1) bins or all vec(y))
%
%
% Parameter estimates
%--------------------------------------------------------------------------
% Ep  - (p x 1)         conditional expectation    E{P|y}
% Cp  - (p x p)         conditional covariance     Cov{P|y}
% Eh  - (q x 1)         conditional log-precisions E{h|y}
%
% log evidence
%--------------------------------------------------------------------------
% F   - [-ve] free energy F = log evidence = p(y|f,g,pE,pC) = p(y|m)
%
%__________________________________________________________________________
% Returns the moments of the posterior p.d.f. of the parameters of a
% nonlinear model specified by IS(P,M,U) under Gaussian assumptions. 
% Usually, IS is an integrator of a dynamic MIMO input-state-output model 
%
%              dx/dt = f(x,u,P)
%              y     = g(x,u,P)  + X0*P0 + e
%
% A static nonlinear observation model with fixed input or causes u
% obtains when x = []. i.e.
%
%              y     = g([],u,P) + X0*P0e + e
%
% but static nonlinear models are specified more simply using
%
%              y     = IS(P,M,U) + X0*P0 + e
%
% Priors on the free parameters P are specified in terms of expectation pE
% and covariance pC. The E-Step uses a Fisher-Scoring scheme and a Laplace
% approximation to estimate the conditional expectation and covariance of P
% If the free-energy starts to increase, a Levenberg-Marquardt scheme is
% invoked.  The M-Step estimates the precision components of e, in terms
% of [Re]ML point estimators of the log-precisions.
%
% An optional feature selection can be specified with parameters M.FS.
%
% For generic aspects of the scheme see:
% 
% Friston K, Mattout J, Trujillo-Barreto N, Ashburner J, Penny W. 
% Variational free energy and the Laplace approximation. 
% NeuroImage. 2007 Jan 1;34(1):220-34.
% 
% This scheme handels complex data along the lines originally described in:
% 
% Sehpard RJ, Lordan BP, and Grant EH. 
% Least squares analysis of complex data with applications to permittivity
% measurements.
% J. Phys. D. Appl. Phys 1970 3:1759-1764.
%
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
 
% Karl Friston
% $Id: spm_nlsi_GN.m 4509 2011-10-04 16:44:54Z guillaume $
 
% figure (unless disabled)
%--------------------------------------------------------------------------
try
    M.nograph;
catch 
    M.nograph = 0;
end
if ~M.nograph
    Fsi = spm_figure('GetWin','SI');
end

% check integrator
%--------------------------------------------------------------------------
try
    M.IS;
catch
    M.IS = 'spm_int';
end
 
% composition of feature selection and prediction (usually an integrator)
%--------------------------------------------------------------------------
try
    
    % try FS(y,M)
    %----------------------------------------------------------------------
    try
        y  = feval(M.FS,Y.y,M);
        IS = inline([M.FS '(' M.IS '(P,M,U),M)'],'P','M','U');
        
    % try FS(y)
    %----------------------------------------------------------------------
    catch
        y  = feval(M.FS,Y.y);
        IS = inline([M.FS '(' M.IS '(P,M,U))'],'P','M','U');
    end
    
catch
    
    % otherwise FS(y) = y
    %----------------------------------------------------------------------
    y   = Y.y;
    IS  = inline([M.IS '(P,M,U)'],'P','M','U');
end
 
% size of data (usually samples x channels)
%--------------------------------------------------------------------------
if iscell(y)
    ns = size(y{1},1);
else
    ns = size(y,1);
end
nr   = length(spm_vec(y))/ns;       % number of samples and responses
M.ns = ns;                          % store in M.ns for integrator
 
% initial states
%--------------------------------------------------------------------------
try
    M.x;
catch
    if ~isfield(M,'n'), M.n = 0;    end
    M.x = sparse(M.n,1);
end
 
% input
%--------------------------------------------------------------------------
try
    U;
catch
    U = [];
end
 
% initial parameters
%--------------------------------------------------------------------------
try
    spm_vec(M.P) - spm_vec(M.pE);
    fprintf('\nParameter initialisation successful\n')
catch
    M.P = M.pE;
end
 
% time-step
%--------------------------------------------------------------------------
try
    Y.dt;
catch
    Y.dt = 1;
end
 
% precision components Q
%--------------------------------------------------------------------------
try
    Q = Y.Q;
    if isnumeric(Q), Q = {Q}; end
catch
    Q = spm_Ce(ns*ones(1,nr));
end
nh    = length(Q);                  % number of precision components
nt    = length(Q{1});               % number of time bins
nq    = nr*ns/nt;                   % for compact Kronecker form of M-step
h     = sparse(nh,1);               % initialise hyperparameters
 
% prior moments
%--------------------------------------------------------------------------
pE    = M.pE;
pC    = M.pC;
 
% confounds (if specified)
%--------------------------------------------------------------------------
try
    nb   = size(Y.X0,1);            % number of bins
    nx   = nr*ns/nb;                % number of blocks
    dfdu = kron(speye(nx,nx),Y.X0);
catch
    dfdu = sparse(ns*nr,0);
end
if isempty(dfdu), dfdu = sparse(ns*nr,0); end
 
% hyperpriors - expectation
%--------------------------------------------------------------------------
try
    hE = M.hE(:);
    if length(hE) ~= nh
        hE = hE(1) + sparse(nh,1);
    end
catch
    hE = sparse(nh,1);
end
h      = hE;
 
% hyperpriors - covariance
%--------------------------------------------------------------------------
try
    ihC = inv(M.hC);
    if length(ihC) ~= nh
        ihC = ihC*speye(nh,nh);
    end
catch
    ihC = speye(nh,nh);
end
 
% unpack covariance
%--------------------------------------------------------------------------
if isstruct(pC);
    pC = spm_diag(spm_vec(pC));
end

% dimension reduction of parameter space
%--------------------------------------------------------------------------
V     = spm_svd(pC,exp(-32));
nu    = size(dfdu,2);                 % number of parameters (confounds)
np    = size(V,2);                    % number of parameters (effective)
ip    = (1:np)';
iu    = (1:nu)' + np;
 
% second-order moments (in reduced space)
%--------------------------------------------------------------------------
pC    = V'*pC*V;
uC    = speye(nu)/1e-8;
ipC   = inv(spm_cat(spm_diag({pC,uC})));
 
% initialize conditional density
%--------------------------------------------------------------------------
Eu    = spm_pinv(dfdu)*spm_vec(y);
p     = [V'*(spm_vec(M.P) - spm_vec(M.pE)); Eu];
Ep    = spm_unvec(spm_vec(pE) + V*p(ip),pE);

 
% EM
%==========================================================================
criterion = [0 0 0 0];

C.F   = -Inf;                                   % free energy
v     = -2;                                     % log ascent rate
dFdh  = zeros(nh,1);
dFdhh = zeros(nh,nh);
for k = 1:128
    
    % time
    %----------------------------------------------------------------------  
    tStart = tic;
 
    % M-Step: ReML estimator of variance components:  h = max{F(p,h)}
    %======================================================================
 
    % prediction f, and gradients; dfdp
    %----------------------------------------------------------------------
    [dfdp f] = spm_diff(IS,Ep,M,U,1,{V});
    
       
    % prediction error and full gradients
    %----------------------------------------------------------------------
    e     =  spm_vec(y) - spm_vec(f) - dfdu*p(iu);
    dfdp  =  reshape(spm_vec(dfdp),ns*nr,np);
    J     = -[dfdp dfdu];
    
 
    % M-step; Fisher scoring scheme to find h = max{F(p,h)}
    %======================================================================
    for m = 1:8
 
        % check for stability
        %------------------------------------------------------------------
        if any(isnan(J(:))) || norm(J,'inf') > exp(32), break, end
        
        % precision and conditional covariance
        %------------------------------------------------------------------
        iS    = sparse(0);
        for i = 1:nh
            iS = iS + Q{i}*(exp(-32) + exp(h(i)));
        end
        S     = spm_inv(iS);
        iS    = kron(speye(nq),iS);
        Pp    = real(J'*iS*J);
        Cp    = spm_inv(Pp + ipC);
        if any(isnan(Cp(:))) || rcond(full(Cp)) < exp(-32), break, end
 
        % precision operators for M-Step
        %------------------------------------------------------------------
        for i = 1:nh
            P{i}   = Q{i}*exp(h(i));
            PS{i}  = P{i}*S;
            P{i}   = kron(speye(nq),P{i});
            JPJ{i} = real(J'*P{i}*J);
        end
                    
        % derivatives: dLdh = dL/dh,...
        %------------------------------------------------------------------
        for i = 1:nh
            dFdh(i,1)      =   trace(PS{i})*nq/2 ...
                             - real(e'*P{i}*e)/2 ...
                             - sum(sum(Cp.*JPJ{i}))/2;
            for j = i:nh
                dFdhh(i,j) = - sum(sum(PS{i}.*PS{j}))*nq/2;
                dFdhh(j,i) =   dFdhh(i,j);
            end
        end
        
        % add hyperpriors
        %------------------------------------------------------------------
        d     = h     - hE;
        dFdh  = dFdh  - ihC*d;
        dFdhh = dFdhh - ihC;
        Ch    = spm_inv(-dFdhh); 
        
        % update ReML estimate
        %------------------------------------------------------------------
        dh    = spm_dx(dFdhh,dFdh,{4});
        dh    = min(max(dh,-1),1);
        h     = h  + dh;
 
        % convergence
        %------------------------------------------------------------------
        dF    = dFdh'*dh;
        if dF < 1e-2, break, end
 
    end

    
    % E-Step with Levenberg-Marquardt regularization
    %======================================================================
 
    % objective function: F(p) (= log evidence - divergence)
    %----------------------------------------------------------------------
    F = - real(e'*iS*e)/2 ...
        - p'*ipC*p/2 ...
        - d'*ihC*d/2 ...
        - ns*nr*log(8*atan(1))/2 ...
        - spm_logdet(S)*nq/2 ...
        + spm_logdet(ipC*Cp)/2 ...
        + spm_logdet(ihC*Ch)/2;
 
    % record increases and reference log-evidence for reporting
    %----------------------------------------------------------------------
    try
        F0; 
        fprintf(' actual: %.3e (%.2f sec)\n',full(F - C.F),toc(tStart))
    catch
        F0 = F;
    end
 
    % if F has increased, update gradients and curvatures for E-Step
    %----------------------------------------------------------------------
    if F > C.F
 
        % accept current estimates
        %------------------------------------------------------------------
        C.p   = p;
        C.h   = h;
        C.F   = F;
        C.Cp  = Cp;
        
        % E-Step: Conditional update of gradients and curvature
        %------------------------------------------------------------------
        dFdp  = -real(J'*iS*e) - ipC*p;
        dFdpp = -real(J'*iS*J) - ipC;
        
        % decrease regularization
        %------------------------------------------------------------------
        v     = min(v + 1/2,4);
        str   = 'EM:(+)';
        
    else
 
        % reset expansion point
        %------------------------------------------------------------------
        p     = C.p;
        h     = C.h;
        Cp    = C.Cp;
 
        % and increase regularization
        %------------------------------------------------------------------
        v     = min(v - 2,-4);
        str   = 'EM:(-)';
        
    end
 
    % E-Step: update
    %======================================================================
    dp    = spm_dx(dFdpp,dFdp,{v});
    p     = p + dp;
    Ep    = spm_unvec(spm_vec(pE) + V*p(ip),pE);
 
    % graphics
    %----------------------------------------------------------------------
    try
 
        % reshape prediction if necessary
        %------------------------------------------------------------------
        f  = reshape(spm_vec(f),ns,nr);
        
        % subplot prediction
        %------------------------------------------------------------------
        set(0,'CurrentFigure',Fsi)
        x    = (1:ns)*Y.dt;
        xLab = 'time (seconds)';
        try
            if length(M.Hz) == ns
                x    = Y.Hz;
                xLab = 'Frequency (Hz)';
            end
        end
        
        if isreal(e)
            
            subplot(2,1,1)
            plot(x,f), hold on
            plot(x,f + spm_unvec(e,f),':'), hold off
            xlabel(xLab)
            title(sprintf('%s: %i','prediction and response: E-Step',k))
            grid on
            
        else
            
            subplot(2,2,1)
            plot(x,real(f)), hold on
            plot(x,real(f + spm_unvec(e,f)),':'), hold off
            xlabel(xLab)
            ylabel('real')
            title(sprintf('%s: %i','prediction and response: E-Step',k))
            grid on
            
            subplot(2,2,2)
            plot(x,imag(f)), hold on
            plot(x,imag(f + spm_unvec(e,f)),':'), hold off
            xlabel(xLab)
            ylabel('imaginary')
            title(sprintf('%s: %i','prediction and response: E-Step',k))
            grid on
            
        end
        
        % subplot parameters
        %------------------------------------------------------------------
        subplot(2,1,2)
        bar(full(V*p(ip)))
        xlabel('parameter')
        title('conditional [minus prior] expectation')
        grid on
        drawnow
        
    catch
        
        if iscell(y)
            
            % subplot predictions
            %--------------------------------------------------------------
            e    = spm_unvec(e,f);
            for i = 1:length(y)
                subplot(2,length(y),i)
                plot(f{i},'ko'), hold on
                plot(f{i} + e{i},'rx'), hold off
                xlabel(sprintf('%s: %i','observation ',i))
                title(sprintf('%s: %i','prediction and response: E-Step',k))
                grid on
            end
            
            % subplot parameters
            %--------------------------------------------------------------
            subplot(2,1,2)
            bar(full(V*p(ip)))
            xlabel('parameter')
            title('conditional [minus prior] expectation')
            grid on
            drawnow
            
        end
        
    end
 
    % convergence
    %----------------------------------------------------------------------
    dF  = dFdp'*dp;
    fprintf('%-6s: %i %6s %-6.3e %6s %.3e ',str,k,'F:',full(C.F - F0),'dF predicted:',full(dF))
    
    criterion = [(dF < 1e-2) criterion(1:end - 1)];
    if all(criterion), fprintf(' convergence\n'), convergence=1; break, end
    convergence=0;
end
 
% outputs
%--------------------------------------------------------------------------
Ep     = spm_unvec(spm_vec(pE) + V*C.p(ip),pE);
Cp     = V*C.Cp(ip,ip)*V';
Eh     = C.h;
F      = C.F;
fit    = f;
parameters = full(V*p(ip));

function [y] = spm_gx_fmri2(x,u,P,M)
% Simulated BOLD response to input
% FORMAT [y] = spm_gx_fmri(x,u,P,M)
% y          - BOLD response (%)
% x          - state vector     (see spm_fx_dcm)
% P          - Parameter vector (see spm_fx_dcm)
% M          - model specification structure (see spm_nlsi)
%__________________________________________________________________________
%
% This function implements the BOLD signal model described in: 
%
% Stephan KE, Weiskopf N, Drysdale PM, Robinson PA, Friston KJ (2007)
% Comparing hemodynamic models with DCM. NeuroImage 38: 387-401.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% modified version MH 2011
% Karl Friston & Klaas Enno Stephan
% $Id: spm_gx_fmri.m 3888 2010-05-15 18:49:56Z karl $
 
 
% Biophysical constants for 1.5T
%==========================================================================
 
% time to echo (TE) (default 0.04 sec)
%--------------------------------------------------------------------------
try, TE = M.TE; catch, TE = 0.04; end
 
% resting venous volume (%)
%--------------------------------------------------------------------------
V0  = 4; % check it!!!

% estimated region-specific ratios of intra- to extra-vascular signal 
%--------------------------------------------------------------------------
ep  = 1; % ? check it!!!
 
% slope r0 of intravascular relaxation rate R_iv as a function of oxygen 
% saturation S:  R_iv = r0*[(1 - S)-(1 - S0)] (Hz)
%--------------------------------------------------------------------------
r0  = 25;
 
% frequency offset at the outer surface of magnetized vessels (Hz)
%--------------------------------------------------------------------------
nu0 = 40.3; 
 
% resting oxygen extraction fraction
%--------------------------------------------------------------------------
E0  = 0.4;
 
%-Coefficients in BOLD signal model
%==========================================================================
k1  = 4.3*nu0*E0*TE;
k2  = ep*r0*E0*TE;
k3  = 1 - ep;
 
%-Output equation of BOLD signal model
%==========================================================================
x   = x';
v   = (x(:,4));
q   = (x(:,5));
y   = V0*(k1.*(1 - q) + k2.*(1 - q./v) + k3.*(1 - v));


function [y P] = spm_fx_fmriD(x,u,P)
% state equation for a dynamic [bilinear/nonlinear/Balloon] model of fMRI
% responses
% FORMAT [y] = spm_fx_fmri(x,u,P,M)

% x      - state vector
%   x(:,1) - excitatory neuronal activity             ue
%   x(:,2) - vascular signal                          s
%   x(:,3) - rCBF                                  ln(f)
%   x(:,4) - venous volume                         ln(v)
%   x(:,5) - deoyxHb                               ln(q)
%  [x(:,6) - inhibitory neuronal activity             ui]
%
% y      - dx/dt
%
%___________________________________________________________________________
%
% References for hemodynamic & neuronal state equations:
% 1. Buxton RB, Wong EC & Frank LR. Dynamics of blood flow and oxygenation
%    changes during brain activation: The Balloon model. MRM 39:855-864,
%    1998.
% 2. Friston KJ, Mechelli A, Turner R, Price CJ. Nonlinear responses in
%    fMRI: the Balloon model, Volterra kernels, and other hemodynamics.
%    Neuroimage 12:466-477, 2000.
% 3. Stephan KE, Kasper L, Harrison LM, Daunizeau J, den Ouden HE,
%    Breakspear M, Friston KJ. Nonlinear dynamic causal models for fMRI.
%    Neuroimage 42:649-662, 2008.
% 4. Marreiros AC, Kiebel SJ, Friston KJ. Dynamic causal modelling for
%    fMRI: a two-state model.
%    Neuroimage. 2008 Jan 1;39(1):269-78.
%__________________________________________________________________________

% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Modified version MH 2011
% Karl Friston & Klaas Enno Stephan
% $Id: spm_fx_fmri.m 3888 2010-05-15 18:49:56Z karl $

xlag = zeros(size(x,2));
Dbi  = find(P.Db==1);
for i=1:length(Dbi)
    DD  =[];
    DD = [x(1,P.Ds(i)),P.Dc{i}];
    xlag(Dbi(i)) = DD(end-1);
    P.Dc{i} = DD(1:end-1);
end
% Neuronal motion
%==========================================================================
B   = full(P.B);                       % bi-linear parameters
C   = P.C/2;                           % exogenous parameters
%D   = full(P.D);                       % nonlinear parameters
x     = x';
% excitatory connections
%--------------------------------------------------------------------------
for i = 1:size(P.B,3)
    P.A = P.A + u(i)*B(:,:,i);
end

% and nonlinear (state) terms
%--------------------------------------------------------------------------
% for i = 1:size(P.D,3)
%     P.A = P.A + x(i,1)*D(:,:,i);
% end

% implement differential state equation y = dx/dt (neuronal)
%--------------------------------------------------------------------------
%try, global sigma; if isempty(sigma), sigma = 1; end; end;
y    = x;
if size(x,2) == 5
    
    % one neuronal state per region
    %----------------------------------------------------------------------
    y(:,1) = sum(P.A.*xlag,2) + C*u(:); % 

else

%     % extrinsic (two neuronal states)
%     %----------------------------------------------------------------------
%     A      = exp(P.A)/8;             % enforce positivity
%     IE     = diag(diag(A));          % inhibitory to excitatory
%     EE     = A - IE;                 % excitatory to excitatory
%     EI     = 1;                      % excitatory to inhibitory
%     SE     = 1;                      % self-inhibition (excitatory)
%     SI     = 2;                      % self-inhibition (inhibitory)
% 
%     % motion - excitatory and inhibitory: y = dx/dt
%     %----------------------------------------------------------------------
%     y(:,1) = EE*x(:,1) - SE*x(:,1) - IE*x(:,6) + P.C*u(:);
%     y(:,6) = EI*x(:,1) - SI*x(:,6);

end

% Hemodynamic motion
%==========================================================================

% hemodynamic parameters
%--------------------------------------------------------------------------
%   H(1) - signal decay                                   d(ds/dt)/ds)
%   H(2) - autoregulation                                 d(ds/dt)/df)
%   H(3) - transit time                                   (t0)
%   H(4) - exponent for Fout(v)                           (alpha)
%   H(5) - resting oxygen extraction                      (E0)
%   H(6) - ratio of intra- to extra-vascular components   (epsilon)
%          of the gradient echo signal
%--------------------------------------------------------------------------
%H        = [0.65 0.41 2.00 0.32 0.34];
try 
   H     = P.H;    
catch    
   H     = [0.64 0.32 2.00 0.32 0.32];
end

% signal decay
%--------------------------------------------------------------------------
sd       = H(:,1);

% transit time
%--------------------------------------------------------------------------
tt       = H(:,3);

% Fout = f(v) - outflow
%--------------------------------------------------------------------------
fv       = x(:,4).^(1./H(:,4));

% e = f(f) - oxygen extraction
%--------------------------------------------------------------------------
ff       = (1 - (1 - H(:,5)).^(1./x(:,3)))./H(:,5);


% implement differential state equation y = dx/dt (hemodynamic)
%--------------------------------------------------------------------------
y(:,2)   = x(:,1) - sd.*x(:,2) - H(:,2).*(x(:,3) - 1);
y(:,3)   = x(:,2);
y(:,4)   = (x(:,3) - fv)./(tt);
y(:,5)   = (ff.*x(:,3) - fv.*x(:,5)./x(:,4))./(tt);
y        = y';
y        = y(:);

return

function cca_plotcausality(M,nodenames,sfac);
%-----------------------------------------------------------------------
% FUNCTION: cca_plotcausality.m
% PURPOSE:  plot causal connectivity given a matrix of Granger causalities
% 
% INPUTS:   M:  log ratio Granger causalities, NaN for nonsig values
%           nodename: vector of node labels (optional)
%           sfac:  scaling factor for linewidth (optional)
%
% OUTPUT:   none
%          
%           Written by Anil K Seth, March 2004
%           Updated AKS August 2004
%           Updated AKS December 2005
%           Updated AKS February 2009
%           Ref: Seth, A.K. (2005) Network: Comp. Neural. Sys. 16(1):35-55
% COPYRIGHT NOTICE AT BOTTOM
%-----------------------------------------------------------------------

% This file is part of GCCAtoolbox.  
% It is Copyright (C) Anil Seth (2005-2009) 
% 
% GCCAtoolbox is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
% 
% GCCAtoolbox is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with GCCAtoolbox.  If not, see <http://www.gnu.org/licenses/>.

nvar = length(M);

M=abs(M); % MG: for DCM models - drawing attenuation connections

% check input
if nargin == 2,
    sfac = 1;
elseif nargin == 1,
    nodenames = [];
end

%   parse input
if(isempty(nodenames)),
    for i=1:nvar,
        nodenames{i} = num2str(i);
    end
end
if nargin < 3,
    sfac = 1;
end
hold off;

% replace bad values with 0
M(isnan(M)) = 0;
M(isinf(M)) = 0;
M(M<0) = 0;

% plot nodes
set(gca,'XTickLabel',[]);
angle = (2*pi)/nvar;
X = zeros(nvar,1); Y = zeros(nvar,1);
for i=1:nvar
    X(i) = cos(angle*(i-1));
    Y(i) = sin(angle*(i-1));
end
for i=1:nvar
    h=plot(X(i),Y(i),'.');
    set(h,'Color','k');
    set(h,'MarkerSize',8);
    hold on;
end
XX=X*1.1;
YY=Y*1.1;

% find reciprocal edges
temp = M.*M';
[toinx,frominx] = find(temp);
R = [toinx frominx];
for i=1:length(toinx),
    xf = X(frominx(i)); yf = Y(frominx(i));
    xt = X(toinx(i)); yt = Y(toinx(i));
    h = line([xf xt],[yf yt]);
    set(h,'LineWidth',M(toinx(i),frominx(i)).*sfac);
    set(h,'Color',[1 0 0]); % red for reciprocal edge   
end
% label nodes
for i=1:nvar
    if(isempty(nodenames))
        h=text(XX(i),YY(i),num2str(i));
    else
        h=text(XX(i),YY(i),nodenames{i});
    end
    set(h,'FontSize',14);
end

% find non-reciprocal edges
[x2,y2] = find(M);
nonR = [x2 y2];
nonR = setdiff(nonR,R,'rows');
toinx = nonR(:,1);
frominx = nonR(:,2);
for i=1:length(toinx),
    xf = X(frominx(i)); yf = Y(frominx(i));
    xt = X(toinx(i)); yt = Y(toinx(i));
    h = line([xf xt],[yf yt]);
    set(h,'LineWidth',M(toinx(i),frominx(i)).*sfac);
    set(h,'Color',[0 1 0]); % green for non-recip edges
    cca_arrowh([xf (xf+xt)/2 ],[yf (yf+yt)/2],'g',400,50);
end

set(gca,'Box','off');
axis('square');
axis off;

xlim([-1.4 1.4]);
ylim([-1.4 1.4]);

function [PR,q] = cca_findsignificance(ret,pval,CFLAG)
%-----------------------------------------------------------------------
% FUNCTION: cca_findsignificance.m
% PURPOSE:  calculate significant Granger causalities
% 
% INPUTS:   ret:    output of cca_granger_regress.m (or similar)
%           pval:   desired P-value threshold
%           CFLAG:  0 = no correction, 1 = Bonferroni, 2 = false discovery
%                   rate, 3 = approx false discovery rate (see ref below)
%
% OUTPUT:   PR:    nvar by nvar matrix of 1s (sig connections)
%           q:     threshold used (=pval for no correction) 
%
%           Written by Anil K Seth, December 2005
%           Updated May 28, 2008
%           Updated Apr 2009 AKS, to return only sig connections
%           Ref: Seth, A.K. (2005) Network: Comp. Neural. Sys. 16(1):35-55
%           Ref: Benjamini, Yoav; Hochberg, Yosef (1995):
%           "Controlling the false discovery rate: a practical and powerful approach 
%           to multiple testing". Journal of the Royal Statistical Society, 
%           Series B (Methodological) 57 (1): 289–300. 
%           FDR code adapted from: http://www.sph.umich.edu/~nichols/FDR/
% COPYRIGHT NOTICE AT BOTTOM
%-----------------------------------------------------------------------

prb = ret.prb;
nvar = length(prb);
PR = zeros(nvar);

if CFLAG == 0,   % no correction
    
    PR(prb<pval) = 1;
    q = pval;

elseif CFLAG == 1, % Bonferroni correction
    
    q = pval./(nvar*(nvar-1));
    PR(prb<q) = 1;
    
elseif CFLAG == 2, %  false discovery rate (assuming independent null hypotheses)

    temp = prb(:);
    p = sort(temp);
    q = pval;
    V = length(p);
    I = (1:V)';
    cVID = 1;
    cVN = sum(1./(1:V));
    pID = p(max(find(p<=I/V*q/cVID)));  % independent or positive dependent
    pN = p(max(find(p<=I/V*q/cVN)));    % nonparametric
    % BUG FIX Nov 02
    if isempty(pID),
        pID = -100;
    end
    % END BUG FIX Nov 02
    PR(prb<pID) = 1;    % use independent or +ve
    
    q = pID;

elseif CFLAG == 3, % 'rough' false discovery rate (RFDR), based on average alpha
    
    m = nvar*(nvar-1);
    q = ((m+1)/(2*m))*pval;
    
    PR(prb<q) = 1;
 
else
    error('CFLAG must be 1-3, in cca_findsignificance.m');
end
            
        
% This file is part of GCCAtoolbox.  
% It is Copyright (C) Anil Seth (2007-2009) and Thomas Nichols
% 
% GCCAtoolbox is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
% 
% GCCAtoolbox is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
% GNU General Public License for more details.
% 
% You should have received a copy of the GNU General Public License
% along with GCCAtoolbox.  If not, see <http://www.gnu.org/licenses/>.
     
