spm_find_pC.m 3.1 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114
  1. function [i,pC,pE,Np] = spm_find_pC(varargin)
  2. % Utility routine that finds the indices of non-zero covariance
  3. % FORMAT [i,pC,pE,Np] = spm_find_pC(pC,pE,fields)
  4. % FORMAT [i,pC,pE,Np] = spm_find_pC(DCM,fields)
  5. % FORMAT [i,pC,pE,Np] = spm_find_pC(DCM)
  6. %
  7. % pC - covariance matrix or variance stucture
  8. % pE - parameter structure
  9. % fields - desired fields of pE
  10. %
  11. % or
  12. %
  13. % DCM - DCM structure
  14. %
  15. % i - find(diag(pC) > TOL)
  16. % rC - reduced covariances
  17. % rE - reduced expectation
  18. %
  19. %__________________________________________________________________________
  20. % Copyright (C) 2015 Wellcome Trust Centre for Neuroimaging
  21. % Karl Friston
  22. % $Id: spm_find_pC.m 7271 2018-03-04 13:11:54Z karl $
  23. %-parse input arguments
  24. %--------------------------------------------------------------------------
  25. if nargin > 2
  26. pC = varargin{1};
  27. pE = varargin{2};
  28. fields = varargin{3};
  29. elseif numel(varargin) > 1
  30. DCM = varargin{1};
  31. fields = varargin{2};
  32. else
  33. DCM = varargin{1};
  34. end
  35. %-get prior density from DCM
  36. %--------------------------------------------------------------------------
  37. if nargin < 3
  38. if ischar(DCM)
  39. DCM = load(DCM,'DCM');
  40. DCM = DCM.DCM;
  41. end
  42. if any(isfield(DCM,{'options','M'}))
  43. try, [pC,pE] = spm_find_rC(DCM); end
  44. end
  45. end
  46. %-Deal with variance structures
  47. %--------------------------------------------------------------------------
  48. if isstruct(pC)
  49. q = spm_vec(pC);
  50. else
  51. q = diag(pC);
  52. end
  53. %-Get indices
  54. %--------------------------------------------------------------------------
  55. i = find(q > mean(q(q < 1024))/1024);
  56. Np = numel(q);
  57. %-subsample fields if necessary
  58. %--------------------------------------------------------------------------
  59. if nargin > 1
  60. if ischar(fields), fields = {fields}; end
  61. if isstruct(pE)
  62. j = spm_fieldindices(pE,fields{:});
  63. if isempty(j) && ~(~isempty(fields) && strcmp(fields{1},'none'))
  64. warning('%s not found. Returning all fields',...
  65. strjoin(cellstr(fields),','));
  66. else
  67. i = j(ismember(j,i));
  68. end
  69. end
  70. end
  71. return
  72. function [pC,pE] = spm_find_rC(DCM)
  73. % FORMAT [pC,pE] = spm_find_rC(DCM)
  74. % model priors
  75. %__________________________________________________________________________
  76. % Get full priors and posteriors
  77. %--------------------------------------------------------------------------
  78. try
  79. pC = DCM.M.pC;
  80. pE = DCM.M.pE;
  81. return
  82. end
  83. % get priors from model specification
  84. %--------------------------------------------------------------------------
  85. if isfield(DCM.options,'spatial')
  86. % EEG or MEG
  87. %----------------------------------------------------------------------
  88. if strcmpi(DCM.options.analysis,'IND')
  89. [pE,~,pC] = spm_ind_priors(DCM.A,DCM.B,DCM.C,DCM.Nf);
  90. else
  91. [pE, pC] = spm_dcm_neural_priors(DCM.A,DCM.B,DCM.C,DCM.options.model);
  92. end
  93. else
  94. % fMRI
  95. %----------------------------------------------------------------------
  96. if ~isfield(DCM,'d')
  97. DCM.d = zeros(size(DCM.a,1),size(DCM.a,1),0);
  98. end
  99. [pE,pC] = spm_dcm_fmri_priors(DCM.a,DCM.b,DCM.c,DCM.d,DCM.options);
  100. end