diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..6d725d9 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +*.m~ diff --git a/README.md b/README.md index 9fd18b9..acf00b4 100644 --- a/README.md +++ b/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). diff --git a/analyzePRF.m b/analyzePRF.m index 3c6f159..db5ab35 100644 --- a/analyzePRF.m +++ b/analyzePRF.m @@ -51,6 +51,9 @@ % (optional) is 'iter' | 'final' | 'off'. default: 'iter'. % (optional) is a typical value for the gain in each time-series. % default: 10. +% (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); diff --git a/analyzePRFcomputesupergridseeds.m b/analyzePRFcomputesupergridseeds.m index 6ff0b7e..6c5fa55 100644 --- a/analyzePRFcomputesupergridseeds.m +++ b/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) % % is [R C] with the resolution of the stimuli % is a cell vector of time x (pixels+1) @@ -11,6 +11,7 @@ % is the dimension that is the time dimension % is a typical value for the gain in each time-series % is [] or a set of noise regressors (cell vector of matrices) +% logical, true use CSS, false make the exponential 1 % % this is an internal function called by analyzePRF.m. this function returns % 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 diff --git a/utilities/fitnonlinearmodel.m b/utilities/fitnonlinearmodel.m index 743b256..768eaad 100644 --- a/utilities/fitnonlinearmodel.m +++ b/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);