From 3aea3b9d91fe178fa22a8de99ca47935cb7f7dbb Mon Sep 17 00:00:00 2001 From: Garikoitz Lerma-Usabiaga Date: Tue, 9 Jul 2019 18:04:51 -0400 Subject: [PATCH 1/8] Update README.md --- README.md | 3 +++ 1 file changed, 3 insertions(+) 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). From e4aa73dea8b109ed0e6c9df079fb4cdcf6a71575 Mon Sep 17 00:00:00 2001 From: Garikoitz Date: Wed, 10 Jul 2019 01:34:15 -0400 Subject: [PATCH 2/8] pushing changes to test with Jon --- .gitignore | 1 + analyzePRF.m | 54 +++++++++++++++++++++++++++++++++++++++++++++------- 2 files changed, 48 insertions(+), 7 deletions(-) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..6d725d9 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +*.m~ diff --git a/analyzePRF.m b/analyzePRF.m index 3d8666d..46224bc 100644 --- a/analyzePRF.m +++ b/analyzePRF.m @@ -171,10 +171,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 @@ -257,7 +257,7 @@ for p=1:length(stimulus) stimulus{p} = squish(stimulus{p},2)'; % frames x pixels stimulus{p} = [stimulus{p} p*ones(size(stimulus{p},1),1)]; % add a dummy column to indicate run breaks - stimulus{p} = single(stimulus{p}); % make single to save memory + % stimulus{p} = single(stimulus{p}); % make single to save memory end % deal with data badness (set bad voxels to be always all 0) @@ -313,6 +313,42 @@ {@(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}}; +%% 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 @@ -492,7 +528,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)); @@ -610,6 +646,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); @@ -670,4 +710,4 @@ % results.hrf = NaN*zeros(numvxs,numinhrf,numfits); % results.hrf(options.vxs,:,:) = permute(a1.params(:,5+(1:numinhrf),:),[3 2 1]); -% results.hrf = reshape(results.hrf, [xyzsize numinhrf numfits]); +% results.hrf = reshape(results.hrf, [xyzsize numinhrf numfits]); \ No newline at end of file From d201f5783a85e54cd5261f1e084628e8d025e940 Mon Sep 17 00:00:00 2001 From: Garikoitz Date: Wed, 10 Jul 2019 16:26:25 -0400 Subject: [PATCH 3/8] latest fitnonlinear from kendrick/analyzePRF and the fix to have the testdata --- utilities/fitnonlinearmodel.m | 57 +++++++++++++++++++++++++---------- 1 file changed, 41 insertions(+), 16 deletions(-) diff --git a/utilities/fitnonlinearmodel.m b/utilities/fitnonlinearmodel.m index 743b256..bc1e05d 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) @@ -621,6 +642,10 @@ % loop over voxels clear results0; +% GLU: he was using parfor here, but for testing with few voxels is anoying. +warning('We removed parfor for testing. To put it back go to function fitnonlinearmodel.m and put it back there') + +% parfor p=1:vnum parfor p=1:vnum % report From 74590bf3f93c51f2593a11a7e21be45262a00275 Mon Sep 17 00:00:00 2001 From: Garikoitz Date: Thu, 11 Jul 2019 09:28:27 -0400 Subject: [PATCH 4/8] default is now without parfor --- utilities/fitnonlinearmodel.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/utilities/fitnonlinearmodel.m b/utilities/fitnonlinearmodel.m index bc1e05d..9af68f9 100644 --- a/utilities/fitnonlinearmodel.m +++ b/utilities/fitnonlinearmodel.m @@ -646,7 +646,7 @@ warning('We removed parfor for testing. To put it back go to function fitnonlinearmodel.m and put it back there') % parfor p=1:vnum -parfor p=1:vnum +for p=1:vnum % report fprintf('*** fitnonlinearmodel: processing voxel %d (%d of %d). ***\n',vxs(p),p,vnum); From b4cf76d1b94a9a5db05240e74568e29e2d4ed9cb Mon Sep 17 00:00:00 2001 From: Garikoitz Date: Mon, 9 Sep 2019 06:28:43 -0700 Subject: [PATCH 5/8] changed param(5) boundaries to [1 1] --- analyzePRF.m | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/analyzePRF.m b/analyzePRF.m index a7386aa..11aa7fd 100644 --- a/analyzePRF.m +++ b/analyzePRF.m @@ -313,13 +313,33 @@ % pre-compute some cache [d,xx,yy] = makegaussian2d(resmx,2,2,2,2); + + +% This is the old version +%{ % 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}}; +%} + + +% This is the new version +% { +% 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 1; + 2*res(1)-1 2*res(2)-1 Inf Inf 1] modelfun} ... + {@(ss)ss [1-res(1)+1 1-res(2)+1 0 0 1; + 2*res(1)-1 2*res(2)-1 Inf Inf 1] @(ss)modelfun}}; +%} + + + + %% Little effort to clarify the modelfun function. This is just a comment. % Inputs of the function From 56b89c4bd377ff2b5034fc5fa95127fe21fdd78e Mon Sep 17 00:00:00 2001 From: Garikoitz Date: Wed, 11 Sep 2019 01:32:21 -0700 Subject: [PATCH 6/8] removed CSS to expnt=1 --- analyzePRF.m | 28 ++++++++++++++++++++++------ analyzePRFcomputesupergridseeds.m | 6 ++++-- 2 files changed, 26 insertions(+), 8 deletions(-) diff --git a/analyzePRF.m b/analyzePRF.m index 11aa7fd..179eab0 100644 --- a/analyzePRF.m +++ b/analyzePRF.m @@ -330,15 +330,25 @@ % { % 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 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} ... - {@(ss)ss [1-res(1)+1 1-res(2)+1 0 0 1; + {@(ss)ss [1-res(1)+1 1-res(2)+1 0 0 NaN; 2*res(1)-1 2*res(2)-1 Inf Inf 1] @(ss)modelfun}}; %} +% After Noah's suggestion +% The NaN in the top row actually indicates that the parameter should not be fit. In the original version there are two rows in the model cell-array because the first row corresponds to the "fit only R, C, S, and G" step of the fits and the second row corresponds to the "start at the previously found point and fit all parameters including N" step of the fitting (this is how analyzePRF works I believe). So you should be able to tell it not to fit the n parameter by setting its lower-bound value to NaN. +%{ +% 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 1] ... + modelfun}}; +%} - - +% To fix it to 1, I think you would need to modify the analyzePRF_computesupergridseeds +% function--the row `expts = [0.5 0.25 0.125];` should be just `expts = [1];` %% Little effort to clarify the modelfun function. This is just a comment. @@ -383,14 +393,20 @@ % generic large seed if ismember(0,options.seedmode) + % edit GLU, we want CSS=1 + % seeds = [seeds; + % (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(0.5) options.typicalgain 0.5]; 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(0.5) options.typicalgain 1]; end % generic small seed if ismember(1,options.seedmode) + % edit GLU, we want CSS=1 + % seeds = [seeds; + % (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(0.5)/10 options.typicalgain 0.5]; 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(0.5)/10 options.typicalgain 1]; end % super-grid seed diff --git a/analyzePRFcomputesupergridseeds.m b/analyzePRFcomputesupergridseeds.m index 6ff0b7e..daab440 100644 --- a/analyzePRFcomputesupergridseeds.m +++ b/analyzePRFcomputesupergridseeds.m @@ -1,4 +1,4 @@ -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) % function [seeds,rvalues] = analyzePRF_computesupergridseeds(res,stimulus,data,modelfun,maxpolydeg,dimdata,dimtime,typicalgain,noisereg) % @@ -27,7 +27,9 @@ % 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]; +% GLU EDIT +% expts = [0.5 0.25 0.125]; +expts = [1]; % calc numvxs = prod(sizefull(data{1},dimdata)); % total number of voxels From c146442bdf48ac559da3125dcb62efae40109fe5 Mon Sep 17 00:00:00 2001 From: Garikoitz Date: Thu, 19 Sep 2019 18:22:41 -0700 Subject: [PATCH 7/8] added option.usecss, now we can choose if CSS implemented or not --- analyzePRF.m | 87 ++++++++++++++++--------------- analyzePRFcomputesupergridseeds.m | 16 ++++-- 2 files changed, 55 insertions(+), 48 deletions(-) diff --git a/analyzePRF.m b/analyzePRF.m index 179eab0..ca5009c 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. % @@ -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); @@ -313,42 +320,22 @@ % pre-compute some cache [d,xx,yy] = makegaussian2d(resmx,2,2,2,2); - - -% This is the old version -%{ % 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) / (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 1] modelfun} ... + {@(ss)ss [1-res(1)+1 1-res(2)+1 0 0 NaN; + 2*res(1)-1 2*res(2)-1 Inf Inf 1] @(ss)modelfun}}; +end -% This is the new version -% { -% 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 1] modelfun} ... - {@(ss)ss [1-res(1)+1 1-res(2)+1 0 0 NaN; - 2*res(1)-1 2*res(2)-1 Inf Inf 1] @(ss)modelfun}}; -%} - -% After Noah's suggestion -% The NaN in the top row actually indicates that the parameter should not be fit. In the original version there are two rows in the model cell-array because the first row corresponds to the "fit only R, C, S, and G" step of the fits and the second row corresponds to the "start at the previously found point and fit all parameters including N" step of the fitting (this is how analyzePRF works I believe). So you should be able to tell it not to fit the n parameter by setting its lower-bound value to NaN. -%{ -% 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 1] ... - modelfun}}; -%} - -% To fix it to 1, I think you would need to modify the analyzePRF_computesupergridseeds -% function--the row `expts = [0.5 0.25 0.125];` should be just `expts = [1];` %% Little effort to clarify the modelfun function. This is just a comment. @@ -391,29 +378,30 @@ % 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) - % edit GLU, we want CSS=1 - % seeds = [seeds; - % (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(0.5) options.typicalgain 0.5]; seeds = [seeds; - (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(0.5) options.typicalgain 1]; + (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(cssexpt) options.typicalgain cssexpt]; end % generic small seed if ismember(1,options.seedmode) - % edit GLU, we want CSS=1 - % seeds = [seeds; - % (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(0.5)/10 options.typicalgain 0.5]; seeds = [seeds; - (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(0.5)/10 options.typicalgain 1]; + (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 @@ -673,6 +661,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]); diff --git a/analyzePRFcomputesupergridseeds.m b/analyzePRFcomputesupergridseeds.m index daab440..6c5fa55 100644 --- a/analyzePRFcomputesupergridseeds.m +++ b/analyzePRFcomputesupergridseeds.m @@ -1,6 +1,6 @@ -function [seeds,rvalues] = analyzePRFcomputesupergridseeds(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,9 +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); -% GLU EDIT -% expts = [0.5 0.25 0.125]; -expts = [1]; + + +% 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 From 871c533c9e6925cbef923a0937afddb96dd3d4a4 Mon Sep 17 00:00:00 2001 From: Garikoitz Date: Sun, 29 Sep 2019 00:40:24 -0700 Subject: [PATCH 8/8] removed second part of model when usecss=false --- analyzePRF.m | 6 ++---- utilities/fitnonlinearmodel.m | 7 ++----- 2 files changed, 4 insertions(+), 9 deletions(-) diff --git a/analyzePRF.m b/analyzePRF.m index ca5009c..db5ab35 100644 --- a/analyzePRF.m +++ b/analyzePRF.m @@ -329,11 +329,9 @@ {@(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) / (2*pi*abs(pp(3))^2))); 0]) .^ posrect(pp(5)),options.hrf,dd(:,prod(res)+1)); + 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} ... - {@(ss)ss [1-res(1)+1 1-res(2)+1 0 0 NaN; - 2*res(1)-1 2*res(2)-1 Inf Inf 1] @(ss)modelfun}}; + 2*res(1)-1 2*res(2)-1 Inf Inf 1] modelfun}}; end diff --git a/utilities/fitnonlinearmodel.m b/utilities/fitnonlinearmodel.m index deadb9f..768eaad 100644 --- a/utilities/fitnonlinearmodel.m +++ b/utilities/fitnonlinearmodel.m @@ -642,11 +642,6 @@ % loop over voxels clear results0; -% GLU: he was using parfor here, but for testing with few voxels is anoying. -% warning('We removed parfor for testing. To put it back go to function fitnonlinearmodel.m and put it back there') -% for p=1:vnum - -% for or parfor here parfor p=1:vnum % report fprintf('*** fitnonlinearmodel: processing voxel %d (%d of %d). ***\n',vxs(p),p,vnum); @@ -682,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);