Added the option usecss true/false

.gitignore

@@ -0,0 +1 @@
+*.m~

README.md

@@ -1,5 +1,8 @@
 # analyzePRF
 
+___Forked by garikoitz to make edits and include it in the tool vistalab/PRFmodel___
+
+
 analyzePRF is a MATLAB toolbox for fitting population receptive field (pRF) models
 to fMRI data.  It is developed by Kendrick Kay (kendrick@post.harvard.edu).
 

analyzePRF.m

@@ -51,6 +51,9 @@
 %   <display> (optional) is 'iter' | 'final' | 'off'.  default: 'iter'.
 %   <typicalgain> (optional) is a typical value for the gain in each time-series.
 %     default: 10.
+%   <usecss> (optional) logical true/false value, if set to false: (1) changes
+%   the boundaries to [NaN 1] and (2) the seed values of the exponentials to be
+%   always 1. Default is true.
 %
 % Analyze pRF data and return the results.
 %
@@ -177,10 +180,10 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INTERNAL CONSTANTS
 
 % define
-remotedir = '/scratch/knk/input/';
-remotedir2 = '/scratch/knk/output/';
-remotelogin = 'knk@login2.chpc.wustl.edu';
-remoteuser = 'knk';
+% remotedir = '/scratch/knk/input/';
+% remotedir2 = '/scratch/knk/output/';
+% remotelogin = 'knk@login2.chpc.wustl.edu';
+% remoteuser = 'knk';
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SETUP AND PREPARATION
 
@@ -244,6 +247,10 @@
 if ~isfield(options,'typicalgain') || isempty(options.typicalgain)
   options.typicalgain = 10;
 end
+if ~isfield(options,'usecss') || isempty(options.usecss)
+  options.usecss = true;
+end
+
 
 % massage
 wantquick = isequal(options.seedmode,-2);
@@ -314,34 +321,85 @@
 [d,xx,yy] = makegaussian2d(resmx,2,2,2,2);
 
 % define the model (parameters are R C S G N)
-modelfun = @(pp,dd) conv2run(posrect(pp(4)) * (dd*[vflatten(placematrix(zeros(res),makegaussian2d(resmx,pp(1),pp(2),abs(pp(3)),abs(pp(3)),xx,yy,0,0) / (2*pi*abs(pp(3))^2))); 0]) .^ posrect(pp(5)),options.hrf,dd(:,prod(res)+1));
-model = {{[] [1-res(1)+1 1-res(2)+1 0    0   NaN;
-              2*res(1)-1 2*res(2)-1 Inf  Inf Inf] modelfun} ...
-         {@(ss)ss [1-res(1)+1 1-res(2)+1 0    0   0;
-                   2*res(1)-1 2*res(2)-1 Inf  Inf Inf] @(ss)modelfun}};
+% Select the correct boundaries depending if we want CSS or not
+if options.usecss
+    modelfun = @(pp,dd) conv2run(posrect(pp(4)) * (dd*[vflatten(placematrix(zeros(res),makegaussian2d(resmx,pp(1),pp(2),abs(pp(3)),abs(pp(3)),xx,yy,0,0) / (2*pi*abs(pp(3))^2))); 0]) .^ posrect(pp(5)),options.hrf,dd(:,prod(res)+1));
+    model = {{[] [1-res(1)+1 1-res(2)+1 0    0   NaN;
+        2*res(1)-1 2*res(2)-1 Inf  Inf Inf] modelfun} ...
+        {@(ss)ss [1-res(1)+1 1-res(2)+1 0    0   0;
+        2*res(1)-1 2*res(2)-1 Inf  Inf Inf] @(ss)modelfun}};
+else
+    modelfun = @(pp,dd) conv2run(posrect(pp(4)) * (dd*[vflatten(placematrix(zeros(res),makegaussian2d(resmx,pp(1),pp(2),abs(pp(3)),abs(pp(3)),xx,yy,0,0))); 0]) .^ posrect(pp(5)),options.hrf,dd(:,prod(res)+1));
+    model = {{[] [1-res(1)+1 1-res(2)+1 0    0   NaN;
+        2*res(1)-1 2*res(2)-1 Inf  Inf 1] modelfun}};
+end
+
+
+%% Little effort to clarify the modelfun function. This is just a comment.
+
+% Inputs of the function
+%   pp = parameters matrix
+%   dd = data matrix
+% Components of the algorithm, from above:
+%   The model involves computing the dot-product between the stimulus and a 2D isotropic
+%   Gaussian, raising the result to an exponent, scaling the result by a gain factor,
+%   and then convolving the result with a hemodynamic response function (HRF).  Polynomial
+%   terms are included (on a run-by-run basis) to model the baseline signal level.
+
+% modelfun = @(pp,dd) ...  % Defines the inputs to the function, params (to be guessed) and data
+%            conv2run(...  % Defines the main function, a convolution. We will want to guess the params
+%               posrect(pp(4)) * ...  % FIRST part of the convolution. It has form A * B. This is A
+%                 (...                  % B starts here. Separate it as well
+%                   dd * ...               % Data matrix
+%                   [vflatten( ...         % Vector. vflatten just returns a vertical vector. Same as kk(:).
+%                         placematrix( ...    % Substitutes matrix 2 into matrix 1 depending on matrix 3
+%                              zeros(res), ...    % Matrix 1
+%                              makegaussian2d(resmx,pp(1),pp(2),abs(pp(3)),abs(pp(3)),xx,yy,0,0) / ... % Matrix 2: element1
+%                              (2*pi*abs(pp(3))^2) ...                                                 % Matrix 2: element2
+%                                     ) ...   % Close placematrix
+%                             )...          % Close vflatten 
+%                     ;0]...               % Adds a 0 to the flattened matrix (vector) at the end
+%                  ).^posrect(pp(5)),...% End of B part B of A*B part of the convolution. This is the parameter in the exponential
+%                options.hrf,...      % SECOND part of the convolution, the HRF signal
+%                dd(:,prod(res)+1) ...% THIRD part of the conv, according knk. Separates convs between runs. This third part is basically a categorization. See help conv2run
+%              );          % Close the convolution function
+% Now the model function needs to be incorporated into a model that lsqcurvefit understands
+% Create two different functions that, fitnonlinear.m will loop over the 2 functions. From above:
+%   A two-stage optimization strategy is used whereby all parameters excluding the
+%   exponent parameter are first optimized (holding the exponent parameter fixed) and 
+%   then all parameters are optimized (including the exponent parameter).  This 
+%   strategy helps avoid local minima.
+
 
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE SEEDS
 
 % init
 seeds = [];
 
+% Select the css expontial for the different seeds
+if options.usecss
+    cssexpt = 0.5;
+else
+    cssexpt = 1;
+end
+
 % generic large seed
 if ismember(0,options.seedmode)
   seeds = [seeds;
-           (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(0.5) options.typicalgain 0.5];
+           (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(cssexpt) options.typicalgain cssexpt];
 end
 
 % generic small seed
 if ismember(1,options.seedmode)
   seeds = [seeds;
-           (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(0.5)/10 options.typicalgain 0.5];
+           (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(cssexpt)/10 options.typicalgain cssexpt];
 end
 
 % super-grid seed
 if any(ismember([2 -2],options.seedmode))
   [supergridseeds,rvalues] = analyzePRFcomputesupergridseeds(res,stimulus,data,modelfun, ...
                                                    options.maxpolydeg,dimdata,dimtime, ...
-                                                   options.typicalgain,noisereg);
+                                                   options.typicalgain,noisereg,options.usecss);
 end
 
 % make a function that individualizes the seeds
@@ -499,7 +557,7 @@
     'wantremovepoly',1, ...
     'extraregressors',noiseregINPUT, ...
     'wantremoveextra',0, ...
-    'dontsave',{{'modelfit' 'opt' 'vxsfull' 'modelpred' 'testdata'}});  % 'resnorms' 'numiters' 
+    'dontsave',{{'opt' 'vxsfull'}}); % GLU: we want to see the fit and testdata that goes to lsqcurvefit
 
           %  'outputfcn',@(a,b,c,d) pause2(.1) | outputfcnsanitycheck(a,b,c,1e-6,10) | outputfcnplot(a,b,c,1,d), ...
           %'outputfcn',@(a,b,c,d) pause2(.1) | outputfcnsanitycheck(a,b,c,1e-6,10) | outputfcnplot(a,b,c,1,d));
@@ -601,6 +659,19 @@
   results.numiters(options.vxs) = a1.numiters;
 end
 
+% Small sanity check. If usecss=false, results.expt need to be always one. 
+% Otherwise, throw an error:
+if ~(options.usecss)
+    if length(unique(results.expt))==1
+        if unique(results.expt) ~= 1
+            error('usecss is set to false, but the value of the exponential is not 1')
+        end
+    else
+        error('usecss is set to false, and there are more than one solutions for the exponential')
+    end
+end
+
+
 % reshape
 results.ang =      reshape(results.ang,      [xyzsize numfits]);
 results.ecc =      reshape(results.ecc,      [xyzsize numfits]);
@@ -617,6 +688,10 @@
 results.params =   paramsA;
 results.options = options;
 
+% added some more, to see model prediction and data that goes into lsqcurvefit
+results.modelpred = a1.modelpred';
+results.testdata  = a1.testdata';
+
 % save 'results' to a temporary file so we don't lose these precious results!
 file0 = [tempname '.mat'];
 fprintf('saving results to %s (just in case).\n',file0);

analyzePRFcomputesupergridseeds.m

@@ -1,6 +1,6 @@
-function [seeds,rvalues] = analyzePRF_computesupergridseeds(res,stimulus,data,modelfun,maxpolydeg,dimdata,dimtime,typicalgain,noisereg)
+function [seeds,rvalues] = analyzePRFcomputesupergridseeds(res,stimulus,data,modelfun,maxpolydeg,dimdata,dimtime,typicalgain,noisereg,usecss)
 
-% function [seeds,rvalues] = analyzePRF_computesupergridseeds(res,stimulus,data,modelfun,maxpolydeg,dimdata,dimtime,typicalgain,noisereg)
+% function [seeds,rvalues] = analyzePRF_computesupergridseeds(res,stimulus,data,modelfun,maxpolydeg,dimdata,dimtime,typicalgain,noisereg,usecss)
 %
 % <res> is [R C] with the resolution of the stimuli
 % <stimulus> is a cell vector of time x (pixels+1)
@@ -11,6 +11,7 @@
 % <dimtime> is the dimension that is the time dimension
 % <typicalgain> is a typical value for the gain in each time-series
 % <noisereg> is [] or a set of noise regressors (cell vector of matrices)
+% <usecss> logical, true use CSS, false make the exponential 1
 %
 % this is an internal function called by analyzePRF.m.  this function returns <seeds>
 % as a matrix of dimensions X x Y x Z x parameters (or XYZ x parameters)
@@ -27,7 +28,14 @@
 % internal constants
 eccs = [0 0.00551 0.014 0.0269 0.0459 0.0731 0.112 0.166 0.242 0.348 0.498 0.707 1];
 angs = linspacecircular(0,2*pi,16);
-expts = [0.5 0.25 0.125];
+
+
+% Select the css expontials for the different seeds
+if usecss
+    expts = [0.5 0.25 0.125];
+else
+    expts = [1];
+end
 
 % calc
 numvxs = prod(sizefull(data{1},dimdata));  % total number of voxels

utilities/fitnonlinearmodel.m

@@ -215,22 +215,39 @@
 % Example 1:
 % 
 % % first, a simple example
-% x = randn(100,1);
-% y = 2*x + 3 + randn(100,1);
-% opt = struct( ...
-%   'stimulus',[x ones(100,1)], ...
-%   'data',y, ...
-%   'model',{{[1 1] [-Inf -Inf; Inf Inf] @(pp,dd) dd*pp'}});
-% results = fitnonlinearmodel(opt);
+%{
+        x = randn(100,1);
+        y = 2*x + 3 + randn(100,1);
+        % Create stimulus, data and model
+        stimulus = [x ones(100,1)];
+        data     = y;
+        % model has the data to feed the parameters in lsqcurvefit
+        % model = {{[1 1]    [-Inf -Inf; Inf Inf]     @(pp,dd) dd*pp'}};
+        % model = {{ x0    [lowerBound; upperBound]       fun}};
+        % x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub)
+        model    = {{[1 1] [-Inf -Inf; Inf Inf] @(pp,dd) dd*pp'}};
+        opt = struct( ...
+                     'stimulus',stimulus, ...
+                     'data',data, ...
+                     'model',model);
+        results = fitnonlinearmodel(opt);
+        % Check the prediction
+        isequal(results.modelpred, (results.params*stimulus')')
+
+%} 
+
 % 
-% % now, try 100 bootstraps
-% opt.resampling = 100;
-% opt.optimoptions = {'Display' 'off'};  % turn off reporting
-% results = fitnonlinearmodel(opt);
+%{
+        opt.resampling = 100;
+        opt.optimoptions = {'Display' 'off'};  % turn off reporting
+        results = fitnonlinearmodel(opt);
+%}
 % 
 % % now, try leave-one-out cross-validation
-% opt.resampling = -(2*(eye(100) - 0.5));
-% results = fitnonlinearmodel(opt);
+%{
+opt.resampling = -(2*(eye(100) - 0.5));
+results = fitnonlinearmodel(opt);
+%} 
 % 
 % Example 2:
 % 
@@ -468,7 +485,8 @@
   opt.wantremoveextra = 1;
 end
 if ~isfield(opt,'dontsave') || (isempty(opt.dontsave) && ~iscell(opt.dontsave))
-  opt.dontsave = {'modelfit' 'numiters' 'resnorms'};
+  % opt.dontsave = {'modelfit' 'numiters' 'resnorms'};
+  opt.dontsave = {'numiters' 'resnorms'};
 end
 if ~iscell(opt.dontsave)
   opt.dontsave = {opt.dontsave};
@@ -582,8 +600,11 @@
 switch resamplingmode
 case 'full'
   trainfun = {@(x) catcell(1,x)};
-  testfun =  {@(x) []};
-case 'xval'
+  % GLU: If testfun is empty, then it does not give the modelpred values back.
+  % Even though it is not correct, let's make testfun the same as trainfun to
+  % have some values. 
+  % testfun =  {@(x) []};
+  testfun = {@(x) catcell(1,x)};case 'xval'
   trainfun = {};
   testfun =  {};
   for p=1:size(opt.resampling,1)
@@ -622,7 +643,6 @@
 % loop over voxels
 clear results0;
 parfor p=1:vnum
-
   % report
   fprintf('*** fitnonlinearmodel: processing voxel %d (%d of %d). ***\n',vxs(p),p,vnum);
   vtime = clock;  % start time for current voxel
@@ -657,6 +677,8 @@
 results.testdata = cat(2,results0.testdata);
 results.modelpred = cat(2,results0.modelpred);
 results.modelfit = cat(3,results0.modelfit);
+results.modelfit = results.modelpred;
+
 results.trainperformance = cat(1,results0.trainperformance).';
 results.testperformance  = cat(1,results0.testperformance).';
 results.aggregatedtestperformance = cat(2,results0.aggregatedtestperformance);