function results = analyzePRF(stimulus,data,tr,options)
<stimulus> provides the apertures as a cell vector of R x C x time.
values should be in [0,1]. the number of time points can differ across runs.
<data> provides the data as a cell vector of voxels x time. can also be
X x Y x Z x time. the number of time points should match the number of
time points in <stimulus>.
<tr> is the TR in seconds (e.g. 1.5)
<options> (optional) is a struct with the following fields:
<vxs> (optional) is a vector of voxel indices to analyze. (we automatically
sort the voxel indices and ensure uniqueness.) default is 1:V where
V is total number of voxels. to reduce computational time, you may want
to create a binary brain mask, perform find() on it, and use the result as <vxs>.
<wantglmdenoise> (optional) is whether to use GLMdenoise to determine
nuisance regressors to add into the PRF model. note that in order to use
this feature, there must be at least two runs (and conditions must repeat
across runs). we automatically determine the GLM design matrix based on
the contents of <stimulus>. special case is to pass in the noise regressors
directly (e.g. from a previous call). default: 0.
<hrf> (optional) is a column vector with the hemodynamic response function (HRF)
to use in the model. the first value of <hrf> should be coincident with the onset
of the stimulus, and the HRF should indicate the timecourse of the response to
a stimulus that lasts for one TR. default is to use a canonical HRF (calculated
using getcanonicalhrf(tr,tr)').
<maxpolydeg> (optional) is a non-negative integer indicating the maximum polynomial
degree to use for drift terms. can be a vector whose length matches the number
of runs in <data>. default is to use round(L/2) where L is the number of minutes
in the duration of a given run.
<seedmode> (optional) is a vector consisting of one or more of the
following values (we automatically sort and ensure uniqueness):
0 means use generic large PRF seed
1 means use generic small PRF seed
2 means use best seed based on super-grid
default: [0 1 2].
<xvalmode> (optional) is
0 means just fit all the data
1 means two-fold cross-validation (first half of runs; second half of runs)
2 means two-fold cross-validation (first half of each run; second half of each run)
default: 0. (note that we round when halving.)
<numperjob> (optional) is
[] means to run locally (not on the cluster)
N where N is a positive integer indicating the number of voxels to analyze in each
cluster job. this option requires a customized computational setup!
default: [].
<maxiter> (optional) is the maximum number of iterations. default: 500.
<display> (optional) is 'iter' | 'final' | 'off'. default: 'iter'.
<typicalgain> (optional) is a typical value for the gain in each time-series.
default: 10.
Analyze pRF data and return the results.
The results structure contains the following fields:
<ang> contains pRF angle estimates. Values range between 0 and 360 degrees.
0 corresponds to the right horizontal meridian, 90 corresponds to the upper vertical
meridian, and so on.
<ecc> contains pRF eccentricity estimates. Values are in pixel units with a lower
bound of 0 pixels.
<rfsize> contains pRF size estimates. pRF size is defined as sigma/sqrt(n) where
sigma is the standard of the 2D Gaussian and n is the exponent of the power-law
function. Values are in pixel units with a lower bound of 0 pixels.
<expt> contains pRF exponent estimates.
<gain> contains pRF gain estimates. Values are in the same units of the data
and are constrained to be non-negative.
<R2> contains R^2 values that indicate the goodness-of-fit of the model to the data.
Values are in percentages and generally range between 0% and 100%. The R^2 values
are computed after projecting out polynomials from both the data and the model fit.
(Because of this projection, R^2 values can sometimes drop below 0%.) Note that
if cross-validation is used (see <xvalmode>), the interpretation of <R2> changes
accordingly.
<resnorms> and <numiters> contain optimization details (residual norms and
number of iterations, respectively).
<meanvol> contains the mean volume, that is, the mean of each voxel's time-series.
<noisereg> contains a record of the noise regressors used in the model.
<params> contains a record of the raw parameter estimates that are obtained internally
in the code. These raw parameters are transformed to a more palatable format for
the user (as described above).
<options> contains a record of the options used in the call to analyzePRF.m.
Details on the pRF model:
- Before analysis, we zero out any voxel that has a non-finite value or has all zeros
in at least one of the runs. This prevents weird issues due to missing or bad data.
- The pRF model that is fit is similar to that described in Dumoulin and Wandell (2008),
except that a static power-law nonlinearity is added to the model. This new model,
called the Compressive Spatial Summation (CSS) model, is described in Kay, Winawer,
Mezer, & Wandell (2013).
- 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.
- The 2D isotropic Gaussian is scaled such that the summation of the values in the
Gaussian is equal to one. This eases the interpretation of the gain of the model.
- The exponent parameter in the model is constrained to be non-negative.
- The gain factor in the model is constrained to be non-negative; this aids the
interpretation of the model (e.g. helps avoid voxels with negative BOLD responses
to the stimuli).
- The workhorse of the analysis is fitnonlinearmodel.m, which is essentially a wrapper
around routines in the MATLAB Optimization Toolbox. We use the Levenberg-Marquardt
algorithm for optimization, minimizing squared error between the model and the data.
- 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.
Regarding GLMdenoise:
- If the <wantglmdenoise> option is specified, we derive noise regressors using
GLMdenoise prior to model fitting. This is done by creating a GLM design matrix
based on the contents of <stimulus> and then using this design matrix in conjunction
with GLMdenoise to analyze the data. The noise regressors identified by GLMdenoise
are then used in the fitting of the pRF models (the regressors enter the model
additively, just like the polynomial regressors).
Regarding seeding issues:
- To minimize the impact of local minima, the default strategy is to perform full
optimizations starting from three different initial seeds.
- The first seed is a generic large pRF that is centered with respect to the stimulus,
has a pRF size equal to 1/4th of the stimulus extent (thus, +/- 2 pRF sizes matches
the stimulus extent), and has an exponent of 0.5.
- The second seed is a generic small pRF that is just like the first seed except has
a pRF size that is 10 times smaller.
- The third seed is a "supergrid" seed that is identified by performing a quick grid
search prior to optimization (similar in spirit to methods described in Dumoulin and
Wandell, 2008). In this procedure, a list of potential seeds is constructed by
exploring a range of eccentricities, angles, and exponents. For each potential
seed, the model prediction is computed, and the seed that produces the closest
match to the data is identified. Note that the supergrid seed may be different
for different voxels.
history:
2014/06/17 - version 1.1
2014/06/15 - add inputs <hrf> and <maxpolydeg>.
2014/06/10 - version 1.0
2014/04/27 - gain seed is now set to 0; add gain to the output
2014/04/29 - use typicalgain now (default 10). allow display input.0001 function results = analyzePRF(stimulus,data,tr,options) 0002 0003 % function results = analyzePRF(stimulus,data,tr,options) 0004 % 0005 % <stimulus> provides the apertures as a cell vector of R x C x time. 0006 % values should be in [0,1]. the number of time points can differ across runs. 0007 % <data> provides the data as a cell vector of voxels x time. can also be 0008 % X x Y x Z x time. the number of time points should match the number of 0009 % time points in <stimulus>. 0010 % <tr> is the TR in seconds (e.g. 1.5) 0011 % <options> (optional) is a struct with the following fields: 0012 % <vxs> (optional) is a vector of voxel indices to analyze. (we automatically 0013 % sort the voxel indices and ensure uniqueness.) default is 1:V where 0014 % V is total number of voxels. to reduce computational time, you may want 0015 % to create a binary brain mask, perform find() on it, and use the result as <vxs>. 0016 % <wantglmdenoise> (optional) is whether to use GLMdenoise to determine 0017 % nuisance regressors to add into the PRF model. note that in order to use 0018 % this feature, there must be at least two runs (and conditions must repeat 0019 % across runs). we automatically determine the GLM design matrix based on 0020 % the contents of <stimulus>. special case is to pass in the noise regressors 0021 % directly (e.g. from a previous call). default: 0. 0022 % <hrf> (optional) is a column vector with the hemodynamic response function (HRF) 0023 % to use in the model. the first value of <hrf> should be coincident with the onset 0024 % of the stimulus, and the HRF should indicate the timecourse of the response to 0025 % a stimulus that lasts for one TR. default is to use a canonical HRF (calculated 0026 % using getcanonicalhrf(tr,tr)'). 0027 % <maxpolydeg> (optional) is a non-negative integer indicating the maximum polynomial 0028 % degree to use for drift terms. can be a vector whose length matches the number 0029 % of runs in <data>. default is to use round(L/2) where L is the number of minutes 0030 % in the duration of a given run. 0031 % <seedmode> (optional) is a vector consisting of one or more of the 0032 % following values (we automatically sort and ensure uniqueness): 0033 % 0 means use generic large PRF seed 0034 % 1 means use generic small PRF seed 0035 % 2 means use best seed based on super-grid 0036 % default: [0 1 2]. 0037 % <xvalmode> (optional) is 0038 % 0 means just fit all the data 0039 % 1 means two-fold cross-validation (first half of runs; second half of runs) 0040 % 2 means two-fold cross-validation (first half of each run; second half of each run) 0041 % default: 0. (note that we round when halving.) 0042 % <numperjob> (optional) is 0043 % [] means to run locally (not on the cluster) 0044 % N where N is a positive integer indicating the number of voxels to analyze in each 0045 % cluster job. this option requires a customized computational setup! 0046 % default: []. 0047 % <maxiter> (optional) is the maximum number of iterations. default: 500. 0048 % <display> (optional) is 'iter' | 'final' | 'off'. default: 'iter'. 0049 % <typicalgain> (optional) is a typical value for the gain in each time-series. 0050 % default: 10. 0051 % 0052 % Analyze pRF data and return the results. 0053 % 0054 % The results structure contains the following fields: 0055 % <ang> contains pRF angle estimates. Values range between 0 and 360 degrees. 0056 % 0 corresponds to the right horizontal meridian, 90 corresponds to the upper vertical 0057 % meridian, and so on. 0058 % <ecc> contains pRF eccentricity estimates. Values are in pixel units with a lower 0059 % bound of 0 pixels. 0060 % <rfsize> contains pRF size estimates. pRF size is defined as sigma/sqrt(n) where 0061 % sigma is the standard of the 2D Gaussian and n is the exponent of the power-law 0062 % function. Values are in pixel units with a lower bound of 0 pixels. 0063 % <expt> contains pRF exponent estimates. 0064 % <gain> contains pRF gain estimates. Values are in the same units of the data 0065 % and are constrained to be non-negative. 0066 % <R2> contains R^2 values that indicate the goodness-of-fit of the model to the data. 0067 % Values are in percentages and generally range between 0% and 100%. The R^2 values 0068 % are computed after projecting out polynomials from both the data and the model fit. 0069 % (Because of this projection, R^2 values can sometimes drop below 0%.) Note that 0070 % if cross-validation is used (see <xvalmode>), the interpretation of <R2> changes 0071 % accordingly. 0072 % <resnorms> and <numiters> contain optimization details (residual norms and 0073 % number of iterations, respectively). 0074 % <meanvol> contains the mean volume, that is, the mean of each voxel's time-series. 0075 % <noisereg> contains a record of the noise regressors used in the model. 0076 % <params> contains a record of the raw parameter estimates that are obtained internally 0077 % in the code. These raw parameters are transformed to a more palatable format for 0078 % the user (as described above). 0079 % <options> contains a record of the options used in the call to analyzePRF.m. 0080 % 0081 % Details on the pRF model: 0082 % - Before analysis, we zero out any voxel that has a non-finite value or has all zeros 0083 % in at least one of the runs. This prevents weird issues due to missing or bad data. 0084 % - The pRF model that is fit is similar to that described in Dumoulin and Wandell (2008), 0085 % except that a static power-law nonlinearity is added to the model. This new model, 0086 % called the Compressive Spatial Summation (CSS) model, is described in Kay, Winawer, 0087 % Mezer, & Wandell (2013). 0088 % - The model involves computing the dot-product between the stimulus and a 2D isotropic 0089 % Gaussian, raising the result to an exponent, scaling the result by a gain factor, 0090 % and then convolving the result with a hemodynamic response function (HRF). Polynomial 0091 % terms are included (on a run-by-run basis) to model the baseline signal level. 0092 % - The 2D isotropic Gaussian is scaled such that the summation of the values in the 0093 % Gaussian is equal to one. This eases the interpretation of the gain of the model. 0094 % - The exponent parameter in the model is constrained to be non-negative. 0095 % - The gain factor in the model is constrained to be non-negative; this aids the 0096 % interpretation of the model (e.g. helps avoid voxels with negative BOLD responses 0097 % to the stimuli). 0098 % - The workhorse of the analysis is fitnonlinearmodel.m, which is essentially a wrapper 0099 % around routines in the MATLAB Optimization Toolbox. We use the Levenberg-Marquardt 0100 % algorithm for optimization, minimizing squared error between the model and the data. 0101 % - A two-stage optimization strategy is used whereby all parameters excluding the 0102 % exponent parameter are first optimized (holding the exponent parameter fixed) and 0103 % then all parameters are optimized (including the exponent parameter). This 0104 % strategy helps avoid local minima. 0105 % 0106 % Regarding GLMdenoise: 0107 % - If the <wantglmdenoise> option is specified, we derive noise regressors using 0108 % GLMdenoise prior to model fitting. This is done by creating a GLM design matrix 0109 % based on the contents of <stimulus> and then using this design matrix in conjunction 0110 % with GLMdenoise to analyze the data. The noise regressors identified by GLMdenoise 0111 % are then used in the fitting of the pRF models (the regressors enter the model 0112 % additively, just like the polynomial regressors). 0113 % 0114 % Regarding seeding issues: 0115 % - To minimize the impact of local minima, the default strategy is to perform full 0116 % optimizations starting from three different initial seeds. 0117 % - The first seed is a generic large pRF that is centered with respect to the stimulus, 0118 % has a pRF size equal to 1/4th of the stimulus extent (thus, +/- 2 pRF sizes matches 0119 % the stimulus extent), and has an exponent of 0.5. 0120 % - The second seed is a generic small pRF that is just like the first seed except has 0121 % a pRF size that is 10 times smaller. 0122 % - The third seed is a "supergrid" seed that is identified by performing a quick grid 0123 % search prior to optimization (similar in spirit to methods described in Dumoulin and 0124 % Wandell, 2008). In this procedure, a list of potential seeds is constructed by 0125 % exploring a range of eccentricities, angles, and exponents. For each potential 0126 % seed, the model prediction is computed, and the seed that produces the closest 0127 % match to the data is identified. Note that the supergrid seed may be different 0128 % for different voxels. 0129 % 0130 % history: 0131 % 2014/06/17 - version 1.1 0132 % 2014/06/15 - add inputs <hrf> and <maxpolydeg>. 0133 % 2014/06/10 - version 1.0 0134 % 2014/04/27 - gain seed is now set to 0; add gain to the output 0135 % 2014/04/29 - use typicalgain now (default 10). allow display input. 0136 0137 % internal notes: 0138 % - for cluster mode, need to make sure fitnonlinearmodel is compiled (compilemcc.m) 0139 % - to check whether local minima are a problem, can look at results.resnorms 0140 0141 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% REPORT 0142 0143 fprintf('*** analyzePRF: started at %s. ***\n',datestr(now)); 0144 stime = clock; % start time 0145 0146 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% INTERNAL CONSTANTS 0147 0148 % define 0149 remotedir = '/scratch/knk/input/'; 0150 remotedir2 = '/scratch/knk/output/'; 0151 remotelogin = 'knk@login2.chpc.wustl.edu'; 0152 remoteuser = 'knk'; 0153 0154 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% SETUP AND PREPARATION 0155 0156 % massage cell inputs 0157 if ~iscell(stimulus) 0158 stimulus = {stimulus}; 0159 end 0160 if ~iscell(data) 0161 data = {data}; 0162 end 0163 0164 % calc 0165 is3d = size(data{1},4) > 1; 0166 if is3d 0167 dimdata = 3; 0168 dimtime = 4; 0169 xyzsize = sizefull(data{1},3); 0170 else 0171 dimdata = 1; 0172 dimtime = 2; 0173 xyzsize = size(data{1},1); 0174 end 0175 numvxs = prod(xyzsize); 0176 0177 % calc 0178 res = sizefull(stimulus{1},2); 0179 resmx = max(res); 0180 numruns = length(data); 0181 0182 % deal with inputs 0183 if ~exist('options','var') || isempty(options) 0184 options = struct(); 0185 end 0186 if ~isfield(options,'vxs') || isempty(options.vxs) 0187 options.vxs = 1:numvxs; 0188 end 0189 if ~isfield(options,'wantglmdenoise') || isempty(options.wantglmdenoise) 0190 options.wantglmdenoise = 0; 0191 end 0192 if ~isfield(options,'hrf') || isempty(options.hrf) 0193 options.hrf = []; 0194 end 0195 if ~isfield(options,'maxpolydeg') || isempty(options.maxpolydeg) 0196 options.maxpolydeg = []; 0197 end 0198 if ~isfield(options,'seedmode') || isempty(options.seedmode) 0199 options.seedmode = [0 1 2]; 0200 end 0201 if ~isfield(options,'xvalmode') || isempty(options.xvalmode) 0202 options.xvalmode = 0; 0203 end 0204 if ~isfield(options,'numperjob') || isempty(options.numperjob) 0205 options.numperjob = []; 0206 end 0207 if ~isfield(options,'maxiter') || isempty(options.maxiter) 0208 options.maxiter = 500; 0209 end 0210 if ~isfield(options,'display') || isempty(options.display) 0211 options.display = 'iter'; 0212 end 0213 if ~isfield(options,'typicalgain') || isempty(options.typicalgain) 0214 options.typicalgain = 10; 0215 end 0216 0217 % massage 0218 options.seedmode = union(options.seedmode(:),[]); 0219 0220 % calc 0221 usecluster = ~isempty(options.numperjob); 0222 0223 % prepare stimuli 0224 for p=1:length(stimulus) 0225 stimulus{p} = squish(stimulus{p},2)'; % frames x pixels 0226 stimulus{p} = [stimulus{p} p*ones(size(stimulus{p},1),1)]; % add a dummy column to indicate run breaks 0227 stimulus{p} = single(stimulus{p}); % make single to save memory 0228 end 0229 0230 % deal with data badness (set bad voxels to be always all 0) 0231 bad = cellfun(@(x) any(~isfinite(x),dimtime) | all(x==0,dimtime),data,'UniformOutput',0); % if non-finite or all 0 0232 bad = any(cat(dimtime,bad{:}),dimtime); % badness in ANY run 0233 for p=1:numruns 0234 data{p}(repmat(bad,[ones(1,dimdata) size(data{p},dimtime)])) = 0; 0235 end 0236 0237 % calc mean volume 0238 meanvol = mean(catcell(dimtime,data),dimtime); 0239 0240 % what HRF should we use? 0241 if isempty(options.hrf) 0242 options.hrf = getcanonicalhrf(tr,tr)'; 0243 end 0244 numinhrf = length(options.hrf); 0245 0246 % what polynomials should we use? 0247 if isempty(options.maxpolydeg) 0248 options.maxpolydeg = cellfun(@(x) round(size(x,dimtime)*tr/60/2),data); 0249 end 0250 if isscalar(options.maxpolydeg) 0251 options.maxpolydeg = repmat(options.maxpolydeg,[1 numruns]); 0252 end 0253 fprintf('using the following maximum polynomial degrees: %s\n',mat2str(options.maxpolydeg)); 0254 0255 % initialize cluster stuff 0256 if usecluster 0257 localfilestodelete = {}; 0258 remotefilestodelete = {}; 0259 end 0260 0261 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FIGURE OUT NOISE REGRESSORS 0262 0263 if isequal(options.wantglmdenoise,1) 0264 noisereg = analyzePRFcomputeGLMdenoiseregressors(stimulus,data,tr); 0265 elseif isequal(options.wantglmdenoise,0) 0266 noisereg = []; 0267 else 0268 noisereg = options.wantglmdenoise; 0269 end 0270 0271 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE MODEL 0272 0273 % pre-compute some cache 0274 [d,xx,yy] = makegaussian2d(resmx,2,2,2,2); 0275 0276 % define the model (parameters are R C S G N) 0277 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)); 0278 model = {{[] [1-res(1)+1 1-res(2)+1 0 0 NaN; 0279 2*res(1)-1 2*res(2)-1 Inf Inf Inf] modelfun} ... 0280 {@(ss)ss [1-res(1)+1 1-res(2)+1 0 0 0; 0281 2*res(1)-1 2*res(2)-1 Inf Inf Inf] @(ss)modelfun}}; 0282 0283 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE SEEDS 0284 0285 % init 0286 seeds = []; 0287 0288 % generic large seed 0289 if ismember(0,options.seedmode) 0290 seeds = [seeds; 0291 (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(0.5) options.typicalgain 0.5]; 0292 end 0293 0294 % generic small seed 0295 if ismember(1,options.seedmode) 0296 seeds = [seeds; 0297 (1+res(1))/2 (1+res(2))/2 resmx/4*sqrt(0.5)/10 options.typicalgain 0.5]; 0298 end 0299 0300 % super-grid seed 0301 if ismember(2,options.seedmode) 0302 supergridseeds = analyzePRFcomputesupergridseeds(res,stimulus,data,modelfun, ... 0303 options.maxpolydeg,dimdata,dimtime, ... 0304 options.typicalgain); 0305 end 0306 0307 % make a function that individualizes the seeds 0308 if exist('supergridseeds','var') 0309 seedfun = @(vx) [[seeds]; 0310 [subscript(squish(supergridseeds,dimdata),{vx ':'})]]; 0311 else 0312 seedfun = @(vx) [seeds]; 0313 end 0314 0315 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE RESAMPLING STUFF 0316 0317 % define wantresampleruns and resampling 0318 switch options.xvalmode 0319 case 0 0320 wantresampleruns = []; 0321 resampling = 0; 0322 case 1 0323 wantresampleruns = 1; 0324 half1 = copymatrix(zeros(1,length(data)),1:round(length(data)/2),1); 0325 half2 = ~half1; 0326 resampling = [(1)*half1 + (-1)*half2; 0327 (-1)*half1 + (1)*half2]; 0328 case 2 0329 wantresampleruns = 0; 0330 resampling = []; 0331 for p=1:length(data) 0332 half1 = copymatrix(zeros(1,size(data{p},2)),1:round(size(data{p},2)/2),1); 0333 half2 = ~half1; 0334 resampling = cat(2,resampling,[(1)*half1 + (-1)*half2; 0335 (-1)*half1 + (1)*half2]); 0336 end 0337 end 0338 0339 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE STIMULUS AND DATA 0340 0341 %%%%% CLUSTER CASE 0342 0343 if usecluster 0344 0345 % save stimulus and transport to the remote server 0346 while 1 0347 filename0 = sprintf('stim%s.mat',randomword(5)); % file name 0348 localfile0 = [tempdir '/' filename0]; % local path to file 0349 remotefile0 = [remotedir '/' filename0]; % remote path to file 0350 0351 % redo if file already exists locally or remotely 0352 if exist(localfile0) || 0==unix(sprintf('ssh %s ls %s',remotelogin,remotefile0)) 0353 continue; 0354 end 0355 0356 % save file and transport it 0357 save(localfile0,'stimulus'); 0358 assert(0==unix(sprintf('rsync -av %s %s:"%s/"',localfile0,remotelogin,remotedir))); 0359 0360 % record 0361 localfilestodelete{end+1} = localfile0; 0362 remotefilestodelete{end+1} = remotefile0; 0363 0364 % stop 0365 break; 0366 end 0367 clear stimulus; % don't let it bleed through anywhere! 0368 0369 % define stimulus 0370 stimulus = @() loadmulti(remotefile0,'stimulus'); 0371 0372 % save data and transport to the remote server 0373 while 1 0374 filename0 = sprintf('data%s',randomword(5)); % directory name that will contain 001.bin, etc. 0375 localfile0 = [tempdir '/' filename0]; % local path to dir 0376 remotefile0 = [remotedir '/' filename0]; % remote path to dir 0377 0378 % redo if dir already exists locally or remotely 0379 if exist(localfile0) || 0==unix(sprintf('ssh %s ls %s',remotelogin,remotefile0)) 0380 continue; 0381 end 0382 0383 % save files and transport them 0384 assert(mkdir(localfile0)); 0385 for p=1:numruns 0386 savebinary([localfile0 sprintf('/%03d.bin',p)],'single',squish(data{p},dimdata)'); % notice squish 0387 end 0388 assert(0==unix(sprintf('rsync -av %s %s:"%s/"',localfile0,remotelogin,remotedir))); 0389 0390 % record 0391 localfilestodelete{end+1} = localfile0; 0392 remotefilestodelete{end+1} = remotefile0; 0393 0394 % stop 0395 break; 0396 end 0397 clear data; 0398 0399 % define data 0400 binfiles = cellfun(@(x) [remotefile0 sprintf('/%03d.bin',x)],num2cell(1:numruns),'UniformOutput',0); 0401 data = @(vxs) cellfun(@(x) double(loadbinary(x,'single',[0 numvxs],-vxs)),binfiles,'UniformOutput',0); 0402 0403 % prepare the output directory name 0404 while 1 0405 filename0 = sprintf('prfresults%s',randomword(5)); 0406 localfile0 = [tempdir '/' filename0]; 0407 remotefile0 = [remotedir2 '/' filename0]; 0408 if exist(localfile0) || 0==unix(sprintf('ssh %s ls %s',remotelogin,remotefile0)) 0409 continue; 0410 end 0411 localfilestodelete{end+1} = localfile0; 0412 localfilestodelete{end+1} = [localfile0 '.mat']; % after consolidation 0413 remotefilestodelete{end+1} = remotefile0; 0414 break; 0415 end 0416 outputdirlocal = localfile0; 0417 outputdirremote = remotefile0; 0418 outputdir = outputdirremote; 0419 0420 %%%%% NON-CLUSTER CASE 0421 0422 else 0423 0424 stimulus = {stimulus}; 0425 data = @(vxs) cellfun(@(x) subscript(squish(x,dimdata),{vxs ':'})',data,'UniformOutput',0); 0426 outputdir = []; 0427 0428 end 0429 0430 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE OPTIONS 0431 0432 % last-minute prep 0433 if iscell(noisereg) 0434 noiseregINPUT = {noisereg}; 0435 else 0436 noiseregINPUT = noisereg; 0437 end 0438 0439 % construct the options struct 0440 opt = struct( ... 0441 'outputdir',outputdir, ... 0442 'stimulus',stimulus, ... 0443 'data',data, ... 0444 'vxs',options.vxs, ... 0445 'model',{model}, ... 0446 'seed',seedfun, ... 0447 'optimoptions',{{'Display' options.display 'Algorithm' 'levenberg-marquardt' 'MaxIter' options.maxiter}}, ... 0448 'wantresampleruns',wantresampleruns, ... 0449 'resampling',resampling, ... 0450 'metric',@calccod, ... 0451 'maxpolydeg',options.maxpolydeg, ... 0452 'wantremovepoly',1, ... 0453 'extraregressors',noiseregINPUT, ... 0454 'wantremoveextra',0, ... 0455 'dontsave',{{'modelfit' 'opt' 'vxsfull' 'modelpred' 'testdata'}}); % 'resnorms' 'numiters' 0456 0457 % 'outputfcn',@(a,b,c,d) pause2(.1) | outputfcnsanitycheck(a,b,c,1e-6,10) | outputfcnplot(a,b,c,1,d), ... 0458 %'outputfcn',@(a,b,c,d) pause2(.1) | outputfcnsanitycheck(a,b,c,1e-6,10) | outputfcnplot(a,b,c,1,d)); 0459 % % debugging: 0460 % chpcstimfile = '/stone/ext1/knk/HCPretinotopy/conimagesB.mat'; 0461 % chpcdatadir2 = outputdir2; % go back 0462 % opt.outputdir='~/temp1'; 0463 % profile on; 0464 % results = fitnonlinearmodel(opt,100,100); 0465 % results = fitnonlinearmodel(opt,1,715233); 0466 % profsave(profile('info'),'~/inout/profile_results'); 0467 % % modelfit = feval(modelfun,results.params,feval(stimulusINPUT)); 0468 % % thedata = feval(dataINPUT,52948); 0469 % % pmatrix = projectionmatrix(constructpolynomialmatrix(304,0:3)); 0470 % % figure; hold on; 0471 % % plot(pmatrix*thedata,'k-'); 0472 % % plot(pmatrix*modelfit,'r-'); 0473 0474 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% FIT MODEL 0475 0476 %%%%% CLUSTER CASE 0477 0478 if usecluster 0479 0480 % submit jobs 0481 jobnames = {}; 0482 jobnames = [jobnames {makedirid(opt.outputdir,1)}]; 0483 jobids = []; 0484 jobids = [jobids chpcrun(jobnames{end},'fitnonlinearmodel',options.numperjob, ... 0485 1,ceil(length(options.vxs)/options.numperjob),[], ... 0486 {'data' 'stimulus' 'bad' 'd' 'xx' 'yy' 'modelfun' 'model'})]; 0487 0488 % record additional files to delete 0489 for p=1:length(jobnames) 0490 remotefilestodelete{end+1} = sprintf('~/sgeoutput/job_%s.*',jobnames{p}); % .o and .e files 0491 remotefilestodelete{end+1} = sprintf('~/mcc/job_%s.mat',jobnames{p}); 0492 localfilestodelete{end+1} = sprintf('~/mcc/job_%s.mat',jobnames{p}); 0493 end 0494 0495 % wait for jobs to finish 0496 sgewaitjobs(jobnames,jobids,remotelogin,remoteuser); 0497 0498 % download the results 0499 assert(0==unix(sprintf('rsync -a %s:"%s" "%s/"',remotelogin,outputdirremote,tempdir))); 0500 0501 % consolidate the results 0502 fitnonlinearmodel_consolidate(outputdirlocal); 0503 0504 % load the results 0505 a1 = load([outputdirlocal '.mat']); 0506 0507 %%%%% NON-CLUSTER CASE 0508 0509 else 0510 0511 a1 = fitnonlinearmodel(opt); 0512 0513 end 0514 0515 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% PREPARE OUTPUT 0516 0517 % calc 0518 numfits = size(a1.params,1); 0519 0520 % init 0521 clear results; 0522 results.ang = NaN*zeros(numvxs,numfits); 0523 results.ecc = NaN*zeros(numvxs,numfits); 0524 results.expt = NaN*zeros(numvxs,numfits); 0525 results.rfsize = NaN*zeros(numvxs,numfits); 0526 results.R2 = NaN*zeros(numvxs,numfits); 0527 results.gain = NaN*zeros(numvxs,numfits); 0528 results.resnorms = cell(numvxs,1); 0529 results.numiters = cell(numvxs,1); 0530 0531 % massage model parameters for output and put in 'results' struct 0532 results.ang(options.vxs,:) = permute(mod(atan2((1+res(1))/2 - a1.params(:,1,:), ... 0533 a1.params(:,2,:) - (1+res(2))/2),2*pi)/pi*180,[3 1 2]); 0534 results.ecc(options.vxs,:) = permute(sqrt(((1+res(1))/2 - a1.params(:,1,:)).^2 + ... 0535 (a1.params(:,2,:) - (1+res(2))/2).^2),[3 1 2]); 0536 results.expt(options.vxs,:) = permute(posrect(a1.params(:,5,:)),[3 1 2]); 0537 results.rfsize(options.vxs,:) = permute(abs(a1.params(:,3,:)) ./ sqrt(posrect(a1.params(:,5,:))),[3 1 2]); 0538 results.R2(options.vxs,:) = permute(a1.trainperformance,[2 1]); 0539 results.gain(options.vxs,:) = permute(posrect(a1.params(:,4,:)),[3 1 2]); 0540 results.resnorms(options.vxs) = a1.resnorms; 0541 results.numiters(options.vxs) = a1.numiters; 0542 0543 % reshape 0544 results.ang = reshape(results.ang, [xyzsize numfits]); 0545 results.ecc = reshape(results.ecc, [xyzsize numfits]); 0546 results.expt = reshape(results.expt, [xyzsize numfits]); 0547 results.rfsize = reshape(results.rfsize, [xyzsize numfits]); 0548 results.R2 = reshape(results.R2, [xyzsize numfits]); 0549 results.gain = reshape(results.gain, [xyzsize numfits]); 0550 results.resnorms = reshape(results.resnorms, [xyzsize 1]); 0551 results.numiters = reshape(results.numiters, [xyzsize 1]); 0552 0553 % add some more stuff 0554 results.meanvol = meanvol; 0555 results.noisereg = noisereg; 0556 results.params = a1.params; 0557 results.options = options; 0558 0559 % save 'results' to a temporary file so we don't lose these precious results! 0560 file0 = [tempname '.mat']; 0561 fprintf('saving results to %s (just in case).\n',file0); 0562 save(file0,'results'); 0563 0564 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% CLEAN UP 0565 0566 %%%%% CLUSTER CASE 0567 0568 if usecluster 0569 0570 % delete local files and directories 0571 for p=1:length(localfilestodelete) 0572 if exist(localfilestodelete{p},'dir') % first dir, then file 0573 rmdir(localfilestodelete{p},'s'); 0574 elseif exist(localfilestodelete{p},'file') 0575 delete(localfilestodelete{p}); 0576 end 0577 end 0578 0579 % delete remote files and directories 0580 for p=1:length(remotefilestodelete) 0581 assert(0==unix(sprintf('ssh %s "rm -rf %s"',remotelogin,remotefilestodelete{p}))); 0582 end 0583 0584 %%%%% NON-CLUSTER CASE 0585 0586 else 0587 0588 end 0589 0590 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% REPORT 0591 0592 fprintf('*** analyzePRF: ended at %s (%.1f minutes). ***\n', ... 0593 datestr(now),etime(clock,stime)/60); 0594 0595 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% JUNK 0596 0597 % % define the model (parameters are R C S G N [HRF]) 0598 % 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)),pp(5+(1:numinhrf))',dd(:,prod(res)+1)); 0599 % model = {{[] [1-res(1)+1 1-res(2)+1 0 0 NaN repmat(NaN,[1 numinhrf]); 0600 % 2*res(1)-1 2*res(2)-1 Inf Inf Inf repmat(Inf,[1 numinhrf])] modelfun} ... 0601 % {@(ss)ss [1-res(1)+1 1-res(2)+1 0 0 0 repmat(NaN,[1 numinhrf]); 0602 % 2*res(1)-1 2*res(2)-1 Inf Inf Inf repmat(Inf,[1 numinhrf])] @(ss)modelfun} ... 0603 % {@(ss)ss [1-res(1)+1 1-res(2)+1 0 0 0 repmat(-Inf,[1 numinhrf]); 0604 % 2*res(1)-1 2*res(2)-1 Inf Inf Inf repmat(Inf,[1 numinhrf])] @(ss)modelfun}}; 0605 % 0606 % % if not fitting the HRF, exclude the last model step 0607 % if ~wantfithrf 0608 % model = model(1:2); 0609 % end 0610 %wantfithrf = 0; % for now, leave at 0 0611 0612 % results.hrf = NaN*zeros(numvxs,numinhrf,numfits); 0613 % results.hrf(options.vxs,:,:) = permute(a1.params(:,5+(1:numinhrf),:),[3 2 1]); 0614 % results.hrf = reshape(results.hrf, [xyzsize numinhrf numfits]);