123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205 |
- function spm_dcm_create(syn_model,source_model,SNR)
- % Create a DCM with simulated data (specified via GUI or an existing model)
- % FORMAT spm_dcm_create(syn_model,source_model,SNR)
- %
- % syn_model - name of the synthetic DCM to be created
- % source_model - define new model ('GUI')
- % or import existing model via file selector ('import')
- % or import existing model (directly specified by directory
- % and name)
- % [default: 'GUI']
- % SNR - signal-to-noise ratio [default: 1]
- %
- % This function allows to create DCM networks with known connectivity
- % parameters from which synthetic data are then generated by calling
- % spm_dcm_generate.
- %
- % This function is very much like spm_dcm_specify_ui but inputs etc. are
- % specified either via the user interface or from an existing model.
- % Currently, the interface provided by this function does not allow for
- % manual specification of nonlinear DCMs; however, these can be imported
- % from existing files.
- %__________________________________________________________________________
- % Copyright (C) 2002-2014 Wellcome Trust Centre for Neuroimaging
- % Will Penny & Klaas Enno Stephan & Peter Zeidman
- % $Id: spm_dcm_create.m 6037 2014-06-04 09:01:16Z peter $
- Finter = spm_figure('GetWin','Interactive');
- header = get(Finter,'Name');
- set(Finter,'Name','Dynamic Causal Modelling')
- %-Check parameters and insert default values, if necessary
- %==========================================================================
- if nargin == 0
- syn_model = spm_input('Name for target DCM_???.mat','+1','s');
- SNR = spm_input('Signal-to-noise ratio (SNR)? ','+1','r',[],1);
- source_model = 'GUI';
- else
- try, source_model; catch, source_model = 'GUI'; end
- try, SNR; catch, SNR = 1; end
- end
- %-Outputs
- %==========================================================================
- switch upper(source_model)
-
- case 'GUI'
-
- %-Define model by GUI
- %==================================================================
-
- % Load SPM for experimental timing
- %------------------------------------------------------------------
- [spm_file, sts] = spm_select(1,'^SPM\.mat$','Select SPM.mat');
- if ~sts, return; end
- try
- SPM = load(spm_file);
- SPM = SPM.SPM;
- catch
- error('Cannot read %s.',spm_file);
- end
-
- % Check for which session to create the DCM
- %------------------------------------------------------------------
- session = length(SPM.Sess);
- if session > 1
- session = spm_input('which session','!+1','n1',1,session);
- if isempty(session) || ~isnumeric(session)
- error('A session number is required');
- end
- end
-
- % Get cell array of region structures
- %------------------------------------------------------------------
- n = spm_input('Enter number of regions',1,'r',[],1);
- for i=1:n
- str = sprintf('Region %d',i);
- xY(i).name = spm_input(['Name for ',str],'+1','s');
- % Make up spurious VOI info
- % for compatibility with spm_dcm_display
- xY(i).xyz = [i i i]'*10;
- xY(i).XYZmm = [i i i]'*10;
- xY(i).s = 1;
- xY(i).spec = 1;
- % for compatibility with spm_dcm_specify
- xY(i).Sess = session;
- xY(i).u = 1;
- xY(i).X0 = [];
- end
-
- % Run through standard specification questions
- %------------------------------------------------------------------
- DCM = spm_dcm_specify_ui(SPM,xY);
-
- % Get desired number of scans
- %------------------------------------------------------------------
- DCM.v = SPM.nscan(session);
-
- % Set default connection strengths to reasonable values
- %------------------------------------------------------------------
- DEFAULT_SELF = 0;
- DEFAULT_BETWEEN = 0.3;
- DEFAULT_DRIVING = 0.8;
-
- def = struct();
- def.A = DEFAULT_BETWEEN .* DCM.a;
- def.A = def.A - diag(diag(def.A));
- def.A = def.A + diag(repmat(DEFAULT_SELF,n,1));
- def.B = DEFAULT_BETWEEN .* DCM.b;
- def.C = DEFAULT_DRIVING .* DCM.c;
-
- % Build A-matrix
- %------------------------------------------------------------------
-
- % Initialize DCM.Ep with prior values
- DCM.Ep = spm_dcm_fmri_priors(DCM.a,DCM.b,DCM.c,DCM.d);
-
- enabled = struct('A',DCM.a,'B',DCM.b,'C',DCM.c);
-
- entry_accepted = false;
- while ~entry_accepted
- % Prompt for A-matrix
- str = 'Connectivity for';
- con = spm_dcm_connectivity_ui(DCM,'A', str, def, enabled);
- DCM.Ep.A = con.A;
-
- % Check stability
- if spm_dcm_check_stability(DCM)
- entry_accepted = true;
- else
- % Query for another attempt
- qtext = ...
- 'Parameters may lead to an unstable estimate. Edit your selection?';
- options = struct('Default','Yes','Interpreter','none');
- choice = questdlg(qtext,'Stability','Yes','No',options);
- entry_accepted = ~strcmp(choice,'Yes');
-
- % Ensure any new values are re-populated
- def = struct('A',DCM.Ep.A,'B',DCM.Ep.B,'C',DCM.Ep.C);
- end
- end
-
- % Build B-matrix and C-matrix
- %------------------------------------------------------------------
- con = spm_dcm_connectivity_ui(DCM,'B', str, def, enabled);
- DCM.Ep.B = con.B;
- con = spm_dcm_connectivity_ui(DCM,'C', str, def, enabled);
- DCM.Ep.C = con.C;
-
- case 'IMPORT'
-
- % Import existing model - prompt user to choose it
- %==================================================================
- P = spm_select(1,'^DCM.*\.mat$','Select source DCM_???.mat');
- load(P)
- otherwise
-
- % Import existing model (directly specified by directory & name)
- %==================================================================
- try
- load(source_model)
- catch
- error('Cannot load source model');
- end
-
- end
- % Now set up output structure
- %--------------------------------------------------------------------------
- X0 = ones(DCM.v,1);
- switch upper(source_model)
- case 'GUI'
- Y = DCM.Y;
- DCM.Ep.decay = sparse(DCM.n,1);
- DCM.Ep.transit = sparse(DCM.n,1);
- DCM.Ep.epsilon = sparse(1,1);
- otherwise
- try
- Y.dt = DCM.Y.dt;
- catch
- Y.dt = DCM.delays(1);
- end
- end
- Y.X0 = X0;
- for i = 1:DCM.n
- Y.name{i} = DCM.xY(i).name;
- end
- Y.Q = spm_Ce(ones(1,DCM.n)*DCM.v);
- DCM.Y = Y;
- %-Save
- %--------------------------------------------------------------------------
- dcm_filename = ['DCM_' syn_model '.mat'];
- save(dcm_filename,'DCM', spm_get_defaults('mat.format'));
- % Generate synthetic output data
- %==========================================================================
- spm_dcm_generate(dcm_filename,SNR);
- spm('FigName',header);
|