diff --git a/mrLoadRet/Plugin/pRF/fminsearch/fminsearchbnd.m b/mrLoadRet/Plugin/pRF/fminsearch/fminsearchbnd.m new file mode 100755 index 000000000..0448eae0f --- /dev/null +++ b/mrLoadRet/Plugin/pRF/fminsearch/fminsearchbnd.m @@ -0,0 +1,307 @@ +function [x,fval,exitflag,output] = fminsearchbnd(fun,x0,LB,UB,options,varargin) +% FMINSEARCHBND: FMINSEARCH, but with bound constraints by transformation +% usage: x=FMINSEARCHBND(fun,x0) +% usage: x=FMINSEARCHBND(fun,x0,LB) +% usage: x=FMINSEARCHBND(fun,x0,LB,UB) +% usage: x=FMINSEARCHBND(fun,x0,LB,UB,options) +% usage: x=FMINSEARCHBND(fun,x0,LB,UB,options,p1,p2,...) +% usage: [x,fval,exitflag,output]=FMINSEARCHBND(fun,x0,...) +% +% arguments: +% fun, x0, options - see the help for FMINSEARCH +% +% LB - lower bound vector or array, must be the same size as x0 +% +% If no lower bounds exist for one of the variables, then +% supply -inf for that variable. +% +% If no lower bounds at all, then LB may be left empty. +% +% Variables may be fixed in value by setting the corresponding +% lower and upper bounds to exactly the same value. +% +% UB - upper bound vector or array, must be the same size as x0 +% +% If no upper bounds exist for one of the variables, then +% supply +inf for that variable. +% +% If no upper bounds at all, then UB may be left empty. +% +% Variables may be fixed in value by setting the corresponding +% lower and upper bounds to exactly the same value. +% +% Notes: +% +% If options is supplied, then TolX will apply to the transformed +% variables. All other FMINSEARCH parameters should be unaffected. +% +% Variables which are constrained by both a lower and an upper +% bound will use a sin transformation. Those constrained by +% only a lower or an upper bound will use a quadratic +% transformation, and unconstrained variables will be left alone. +% +% Variables may be fixed by setting their respective bounds equal. +% In this case, the problem will be reduced in size for FMINSEARCH. +% +% The bounds are inclusive inequalities, which admit the +% boundary values themselves, but will not permit ANY function +% evaluations outside the bounds. These constraints are strictly +% followed. +% +% If your problem has an EXCLUSIVE (strict) constraint which will +% not admit evaluation at the bound itself, then you must provide +% a slightly offset bound. An example of this is a function which +% contains the log of one of its parameters. If you constrain the +% variable to have a lower bound of zero, then FMINSEARCHBND may +% try to evaluate the function exactly at zero. +% +% +% Example usage: +% rosen = @(x) (1-x(1)).^2 + 105*(x(2)-x(1).^2).^2; +% +% fminsearch(rosen,[3 3]) % unconstrained +% ans = +% 1.0000 1.0000 +% +% fminsearchbnd(rosen,[3 3],[2 2],[]) % constrained +% ans = +% 2.0000 4.0000 +% +% See test_main.m for other examples of use. +% +% +% See also: fminsearch, fminspleas +% +% +% Author: John D'Errico +% E-mail: woodchips@rochester.rr.com +% Release: 4 +% Release date: 7/23/06 + +% size checks +xsize = size(x0); +x0 = x0(:); +n=length(x0); + +if (nargin<3) || isempty(LB) + LB = repmat(-inf,n,1); +else + LB = LB(:); +end +if (nargin<4) || isempty(UB) + UB = repmat(inf,n,1); +else + UB = UB(:); +end + +if (n~=length(LB)) || (n~=length(UB)) + error 'x0 is incompatible in size with either LB or UB.' +end + +% set default options if necessary +if (nargin<5) || isempty(options) + options = optimset('fminsearch'); +end + +% stuff into a struct to pass around +params.args = varargin; +params.LB = LB; +params.UB = UB; +params.fun = fun; +params.n = n; +% note that the number of parameters may actually vary if +% a user has chosen to fix one or more parameters +params.xsize = xsize; +params.OutputFcn = []; + +% 0 --> unconstrained variable +% 1 --> lower bound only +% 2 --> upper bound only +% 3 --> dual finite bounds +% 4 --> fixed variable +params.BoundClass = zeros(n,1); +for i=1:n + k = isfinite(LB(i)) + 2*isfinite(UB(i)); + params.BoundClass(i) = k; + if (k==3) && (LB(i)==UB(i)) + params.BoundClass(i) = 4; + end +end + +% transform starting values into their unconstrained +% surrogates. Check for infeasible starting guesses. +x0u = x0; +k=1; +for i = 1:n + switch params.BoundClass(i) + case 1 + % lower bound only + if x0(i)<=LB(i) + % infeasible starting value. Use bound. + x0u(k) = 0; + else + x0u(k) = sqrt(x0(i) - LB(i)); + end + + % increment k + k=k+1; + case 2 + % upper bound only + if x0(i)>=UB(i) + % infeasible starting value. use bound. + x0u(k) = 0; + else + x0u(k) = sqrt(UB(i) - x0(i)); + end + + % increment k + k=k+1; + case 3 + % lower and upper bounds + if x0(i)<=LB(i) + % infeasible starting value + x0u(k) = -pi/2; + elseif x0(i)>=UB(i) + % infeasible starting value + x0u(k) = pi/2; + else + x0u(k) = 2*(x0(i) - LB(i))/(UB(i)-LB(i)) - 1; + % shift by 2*pi to avoid problems at zero in fminsearch + % otherwise, the initial simplex is vanishingly small + x0u(k) = 2*pi+asin(max(-1,min(1,x0u(k)))); + end + + % increment k + k=k+1; + case 0 + % unconstrained variable. x0u(i) is set. + x0u(k) = x0(i); + + % increment k + k=k+1; + case 4 + % fixed variable. drop it before fminsearch sees it. + % k is not incremented for this variable. + end + +end +% if any of the unknowns were fixed, then we need to shorten +% x0u now. +if k<=n + x0u(k:n) = []; +end + +% were all the variables fixed? +if isempty(x0u) + % All variables were fixed. quit immediately, setting the + % appropriate parameters, then return. + + % undo the variable transformations into the original space + x = xtransform(x0u,params); + + % final reshape + x = reshape(x,xsize); + + % stuff fval with the final value + fval = feval(params.fun,x,params.args{:}); + + % fminsearchbnd was not called + exitflag = 0; + + output.iterations = 0; + output.funcCount = 1; + output.algorithm = 'fminsearch'; + output.message = 'All variables were held fixed by the applied bounds'; + + % return with no call at all to fminsearch + return +end + +% Check for an outputfcn. If there is any, then substitute my +% own wrapper function. +if ~isempty(options.OutputFcn) + params.OutputFcn = options.OutputFcn; + options.OutputFcn = @outfun_wrapper; +end + +% now we can call fminsearch, but with our own +% intra-objective function. +[xu,fval,exitflag,output] = fminsearch(@intrafun,x0u,options,params); + +% undo the variable transformations into the original space +x = xtransform(xu,params); + +% final reshape to make sure the result has the proper shape +x = reshape(x,xsize); + +% Use a nested function as the OutputFcn wrapper + function stop = outfun_wrapper(x,varargin); + % we need to transform x first + xtrans = xtransform(x,params); + + % then call the user supplied OutputFcn + stop = params.OutputFcn(xtrans,varargin{1:(end-1)}); + + end + +end % mainline end + +% ====================================== +% ========= begin subfunctions ========= +% ====================================== +function fval = intrafun(x,params) +% transform variables, then call original function + +% transform +xtrans = xtransform(x,params); + +% and call fun +fval = feval(params.fun,reshape(xtrans,params.xsize),params.args{:}); + +end % sub function intrafun end + +% ====================================== +function xtrans = xtransform(x,params) +% converts unconstrained variables into their original domains + +xtrans = zeros(params.xsize); +% k allows some variables to be fixed, thus dropped from the +% optimization. +k=1; +for i = 1:params.n + switch params.BoundClass(i) + case 1 + % lower bound only + xtrans(i) = params.LB(i) + x(k).^2; + + k=k+1; + case 2 + % upper bound only + xtrans(i) = params.UB(i) - x(k).^2; + + k=k+1; + case 3 + % lower and upper bounds + xtrans(i) = (sin(x(k))+1)/2; + xtrans(i) = xtrans(i)*(params.UB(i) - params.LB(i)) + params.LB(i); + % just in case of any floating point problems + xtrans(i) = max(params.LB(i),min(params.UB(i),xtrans(i))); + + k=k+1; + case 4 + % fixed variable, bounds are equal, set it at either bound + xtrans(i) = params.LB(i); + case 0 + % unconstrained variable. + xtrans(i) = x(k); + + k=k+1; + end +end + +end % sub function xtransform end + + + + + diff --git a/mrLoadRet/Plugin/pRF/models/pRFModelTemplate.m b/mrLoadRet/Plugin/pRF/models/pRFModelTemplate.m new file mode 100644 index 000000000..3e1cc9ca0 --- /dev/null +++ b/mrLoadRet/Plugin/pRF/models/pRFModelTemplate.m @@ -0,0 +1,123 @@ +% pRFModelTemplate.m +% +% $Id:$ +% usage: pRFModelTemplate(varargin) +% by: akshay jagadeesh +% date: 08/23/16 +% purpose: Template file to create new models +% +% - Allows you to create a new model type +% - Simply add the model type to pRFGUI and create a +% file like this one for your model with the appropriate +% methods filled in. +% +% This model template is set up for the standard Gaussian model. + +function output = pRFModelTemplate(varargin) + +if nargin <=2 + disp(sprintf('Not enough arguments')); + return +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%pRFModelTemplate(command, fitParams, rfModel, hrf, i)%%%%% +%%%% Called from getModelResidual %%%%% +if strcmp(varargin{1}, 'getModelResponse') + + fitParams = varargin{2}; + rfModel = varargin{3}; + hrf = varargin{4}; + i = varargin{5}; + + nFrames = fitParams.concatInfo.runTransition(i,2); + thisModelResponse = convolveModelWithStimulus(rfModel,fitParams.stim{i},nFrames); + + % and convolve in time. + thisModelResponse = convolveModelResponseWithHRF(thisModelResponse,hrf); + + % drop junk frames here + thisModelResponse = thisModelResponse(fitParams.concatInfo.totalJunkedFrames(i)+1:end); + + % return the calculated model response + output = thisModelResponse; +%end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%pRFModelTemplate(command, fitParams)%%%%% +%%%% Called from setFitParams %%%%% + +elseif strcmp(varargin{1}, 'setParams') + + fitParams = varargin{2}; + + fitParams.paramNames = {'x','y','rfWidth'}; + fitParams.paramDescriptions = {'RF x position','RF y position','RF width (std of gaussian)'}; + fitParams.paramIncDec = [1 1 1]; + fitParams.paramMin = [-inf -inf 0]; + fitParams.paramMax = [inf inf inf]; + % set min/max and init + fitParams.minParams = [fitParams.stimExtents(1) fitParams.stimExtents(2) 0]; + fitParams.maxParams = [fitParams.stimExtents(3) fitParams.stimExtents(4) inf]; + fitParams.initParams = [0 0 4]; + + % return fitParams with modified values + output = fitParams; +%end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%pRFModelTemplate(command, fitParams, params)%%%%% +%%%% Called from getFitParams %%%%% + +elseif strcmp(varargin{1}, 'getFitParams') + + fitParams = varargin{2}; + params = varargin{3}; + + p.rfType = fitParams.rfType; + + % Define your parameters here + p.x = params(1); + p.y = params(2); + p.std = params(3); + % use a fixed single gaussian + p.canonical.type = 'gamma'; + p.canonical.lengthInSeconds = 25; + p.canonical.timelag = fitParams.timelag; + p.canonical.tau = fitParams.tau; + p.canonical.exponent = fitParams.exponent; + p.canonical.offset = 0; + p.canonical.diffOfGamma = fitParams.diffOfGamma; + p.canonical.amplitudeRatio = fitParams.amplitudeRatio; + p.canonical.timelag2 = fitParams.timelag2; + p.canonical.tau2 = fitParams.tau2; + p.canonical.exponent2 = fitParams.exponent2; + p.canonical.offset2 = 0; + +end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelWithStimulus %% +function modelResponse = convolveModelWithStimulus(rfModel,stim,nFrames) + +% get number of frames +nStimFrames = size(stim.im,3); + +% preallocate memory +modelResponse = zeros(1,nFrames); + +for frameNum = 1:nStimFrames + % multipy the stimulus frame by frame with the rfModel + % and take the sum + modelResponse(frameNum) = sum(sum(rfModel.*stim.im(:,:,frameNum))); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelResponseWithHRF %% +function modelTimecourse = convolveModelResponseWithHRF(modelTimecourse,hrf) + +n = length(modelTimecourse); +modelTimecourse = conv(modelTimecourse,hrf.hrf); +modelTimecourse = modelTimecourse(1:n); + diff --git a/mrLoadRet/Plugin/pRF/models/pRF_DoG_CSS.m b/mrLoadRet/Plugin/pRF/models/pRF_DoG_CSS.m new file mode 100644 index 000000000..d0407f21b --- /dev/null +++ b/mrLoadRet/Plugin/pRF/models/pRF_DoG_CSS.m @@ -0,0 +1,135 @@ +% pRF_DoG_CSS.m +% +% $Id:$ +% usage: pRF_DoG_CSS(varargin) +% by: akshay jagadeesh +% date: 09/02/16 +% purpose: Model of Difference of Gaussians with Compressive Static Non linearity +% +% - Allows you to create a new model type +% - Simply add the model type to pRFGUI and create a +% file like this one for your model with the appropriate +% methods filled in. +% +% This model template is set up for the standard Gaussian model. + +function output = pRF_DoG_CSS(varargin) + +if nargin < 2 + disp(sprintf('Not enough arguments')); + disp(sprintf('Number of arguments: %d', nargin)); + celldisp(varargin) + return +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% pRF_gaussian('getModelResponse', fitParams, rfModel, hrf, p, i) %%%% +%%% Called from getModelResidual %%%% +if strcmp(varargin{1}, 'getModelResponse') + + fitParams = varargin{2}; + rfModel1 = varargin{3}; + hrf = varargin{4}; + p = varargin{5}; + i = varargin{6}; + + rfModel2 = exp(-(((fitParams.stimX - p.x).^2)/(2*(p.std2^2))+((fitParams.stimY-p.y).^2)/(2*(p.std2^2)))); + nFrames = fitParams.concatInfo.runTransition(i,2); + rPlus = convolveModelWithStimulus(rfModel1,fitParams.stim{i},nFrames); + rMinus = convolveModelWithStimulus(rfModel2, fitParams.stim{i}, nFrames); + + % and convolve in time. + pPlus = convolveModelResponseWithHRF(rPlus,hrf); + pMinus = convolveModelResponseWithHRF(rMinus,hrf); + + % Model response is the difference of Gaussians, weighted by Beta amplitudes + thisModelResponse = power(pPlus*p.B1 + pMinus*p.B2, p.exp); + + % drop junk frames here + thisModelResponse = thisModelResponse(fitParams.concatInfo.totalJunkedFrames(i)+1:end); + + % return the calculated model response + output = thisModelResponse; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% pRF_gaussian('setParams', fitParams) %%%% +%%% Called from setParams %%%% + +elseif strcmp(varargin{1}, 'setParams') + + fitParams = varargin{2}; + + fitParams.paramNames = {'x','y','rfWidth', 'surroundWidth', 'centerAmplitude', 'surroundAmplitude', 'exponent'}; + fitParams.paramDescriptions = {'RF x position','RF y position','RF width (std of gaussian)', 'RF width of surround pool', 'Beta weight amplitude of positive Gaussian', 'Beta weight amplitude for negative Gaussian', 'Exponent static nonlinearity'}; + fitParams.paramIncDec = [1 1 1 1 1 1 1]; + fitParams.paramMin = [-inf -inf 0 0 -inf -inf -inf]; + fitParams.paramMax = [inf inf inf inf inf inf inf]; + % set min/max and init + fitParams.minParams = [fitParams.stimExtents(1) fitParams.stimExtents(2) 0 0 0 -inf -5]; + fitParams.maxParams = [fitParams.stimExtents(3) fitParams.stimExtents(4) inf inf inf 0 5]; + fitParams.initParams = [0 0 4 4 1 0 1]; + + % return fitParams with modified values + output = fitParams; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%% pRF_gaussian(command, fitParams, params) %%%%% +%%%% Called from getFitParams %%%%% + +elseif strcmp(varargin{1}, 'getFitParams') + + fitParams = varargin{2}; + params = varargin{3}; + + p.rfType = fitParams.rfType; + + % Define your parameters here + p.x = params(1); + p.y = params(2); + p.std = params(3); + p.std2 = params(4); + p.B1 = params(5); + p.B2 = params(6); + p.exp = params(7); + % use a fixed single gaussian + p.canonical.type = 'gamma'; + p.canonical.lengthInSeconds = 25; + p.canonical.timelag = fitParams.timelag; + p.canonical.tau = fitParams.tau; + p.canonical.exponent = fitParams.exponent; + p.canonical.offset = 0; + p.canonical.diffOfGamma = fitParams.diffOfGamma; + p.canonical.amplitudeRatio = fitParams.amplitudeRatio; + p.canonical.timelag2 = fitParams.timelag2; + p.canonical.tau2 = fitParams.tau2; + p.canonical.exponent2 = fitParams.exponent2; + p.canonical.offset2 = 0; + + output = struct(p); + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelWithStimulus %% +function modelResponse = convolveModelWithStimulus(rfModel,stim,nFrames) + +% get number of frames +nStimFrames = size(stim.im,3); + +% preallocate memory +modelResponse = zeros(1,nFrames); + +for frameNum = 1:nStimFrames + % multipy the stimulus frame by frame with the rfModel + % and take the sum + modelResponse(frameNum) = sum(sum(rfModel.*stim.im(:,:,frameNum))); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelResponseWithHRF %% +function modelTimecourse = convolveModelResponseWithHRF(modelTimecourse,hrf) + +n = length(modelTimecourse); +modelTimecourse = conv(modelTimecourse,hrf.hrf); +modelTimecourse = modelTimecourse(1:n); + diff --git a/mrLoadRet/Plugin/pRF/models/pRF_diffGaussian.m b/mrLoadRet/Plugin/pRF/models/pRF_diffGaussian.m new file mode 100644 index 000000000..5436bb88c --- /dev/null +++ b/mrLoadRet/Plugin/pRF/models/pRF_diffGaussian.m @@ -0,0 +1,142 @@ +% pRF_gaussian.m +% +% $Id:$ +% usage: pRF_gaussian(varargin) +% by: akshay jagadeesh +% date: 09/02/16 +% purpose: Template file to create new models +% +% - Allows you to create a new model type +% - Simply add the model type to pRFGUI and create a +% file like this one for your model with the appropriate +% methods filled in. +% +% This model template is set up for the standard Gaussian model. + +function output = pRF_diffGaussian(varargin) + +if nargin < 2 + disp(sprintf('Not enough arguments')); + disp(sprintf('Number of arguments: %d', nargin)); + celldisp(varargin) + return +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% pRF_gaussian('getModelResponse', fitParams, rfModel, hrf, p, i) %%%% +%%% Called from getModelResidual %%%% +if strcmp(varargin{1}, 'getModelResponse') + + fitParams = varargin{2}; + rfModel1 = varargin{3}; + hrf = varargin{4}; + p = varargin{5}; + i = varargin{6}; + + rfModel2 = exp(-(((fitParams.stimX - p.x).^2)/(2*((p.std*p.stdRatio)^2))+((fitParams.stimY-p.y).^2)/(2*((p.std*p.stdRatio)^2)))); + nFrames = fitParams.concatInfo.runTransition(i,2); + rPlus = convolveModelWithStimulus(rfModel1,fitParams.stim{i},nFrames); + rMinus = convolveModelWithStimulus(rfModel2, fitParams.stim{i}, nFrames); + + thisModelResponse = rPlus*p.B1 - rMinus*p.B2; + thisModelResponse = convolveModelResponseWithHRF(thisModelResponse, hrf); + + % and convolve in time. + %pPlus = convolveModelResponseWithHRF(rPlus,hrf); + %pMinus = convolveModelResponseWithHRF(rMinus,hrf); + + % Model response is the difference of Gaussians, weighted by Beta amplitudes + %thisModelResponse = pPlus*p.B1 + pMinus*p.B2; + + % drop junk frames here + %thisModelResponse = thisModelResponse(fitParams.concatInfo.totalJunkedFrames(i)+1:end); + if fitParams.concatInfo.isConcat + thisModelResponse = thisModelResponse(fitParams.concatInfo.junkFrames(i)+1:end); + end + + + % return the calculated model response + output = thisModelResponse; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% pRF_gaussian('setParams', fitParams) %%%% +%%% Called from setParams %%%% + +elseif strcmp(varargin{1}, 'setParams') + + fitParams = varargin{2}; + + fitParams.paramNames = {'x','y','rfWidth', 'surroundRatio', 'centerGain', 'surroundGain'}; + fitParams.paramDescriptions = {'RF x position','RF y position','RF width (std of gaussian)', 'Ratio of surround width to center width', 'Center Gain', 'Surround Gain'}; + + fitParams.paramIncDec = [1 1 1 1 1 1]; + fitParams.paramMin = [-inf -inf 0 0 -inf -inf]; + fitParams.paramMax = [inf inf inf inf inf inf]; + % set min/max and init + fitParams.minParams = [fitParams.stimExtents(1) fitParams.stimExtents(2) 0 1 0 0]; + fitParams.maxParams = [fitParams.stimExtents(3) fitParams.stimExtents(4) inf inf inf inf]; + fitParams.initParams = [0 0 4 2 1 1]; + + % return fitParams with modified values + output = fitParams; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%% pRF_gaussian(command, fitParams, params) %%%%% +%%%% Called from getFitParams %%%%% + +elseif strcmp(varargin{1}, 'getFitParams') + + fitParams = varargin{2}; + params = varargin{3}; + + p.rfType = fitParams.rfType; + + % Define your parameters here + p.x = params(1); + p.y = params(2); + p.std = params(3); + p.stdRatio = params(4); + p.B1 = params(5); + p.B2 = params(6); + % use a fixed single gaussian + p.canonical.type = 'gamma'; + p.canonical.lengthInSeconds = 25; + p.canonical.timelag = fitParams.timelag; + p.canonical.tau = fitParams.tau; + p.canonical.exponent = fitParams.exponent; + p.canonical.offset = 0; + p.canonical.diffOfGamma = fitParams.diffOfGamma; + p.canonical.amplitudeRatio = fitParams.amplitudeRatio; + p.canonical.timelag2 = fitParams.timelag2; + p.canonical.tau2 = fitParams.tau2; + p.canonical.exponent2 = fitParams.exponent2; + p.canonical.offset2 = 0; + + output = struct(p); + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelWithStimulus %% +function modelResponse = convolveModelWithStimulus(rfModel,stim,nFrames) + +% get number of frames +nStimFrames = size(stim.im,3); + +% preallocate memory +modelResponse = zeros(1,nStimFrames); + +for frameNum = 1:nStimFrames + % multipy the stimulus frame by frame with the rfModel + % and take the sum + modelResponse(frameNum) = sum(sum(rfModel.*stim.im(:,:,frameNum))); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelResponseWithHRF %% +function modelTimecourse = convolveModelResponseWithHRF(modelTimecourse,hrf) + +n = length(modelTimecourse); +modelTimecourse = conv(modelTimecourse,hrf.hrf); +modelTimecourse = modelTimecourse(1:n); + diff --git a/mrLoadRet/Plugin/pRF/models/pRF_divGaussian.m b/mrLoadRet/Plugin/pRF/models/pRF_divGaussian.m new file mode 100644 index 000000000..691aab70c --- /dev/null +++ b/mrLoadRet/Plugin/pRF/models/pRF_divGaussian.m @@ -0,0 +1,137 @@ +% pRF_divGaussian.m +% +% $Id:$ +% usage: pRF_gaussian(varargin) +% by: akshay jagadeesh +% date: 09/02/16 +% purpose: Template file to create new models +% +% - Allows you to create a new model type +% - Simply add the model type to pRFGUI and create a +% file like this one for your model with the appropriate +% methods filled in. +% +% This model template is set up for the standard Gaussian model. + +function output = pRF_divGaussian(varargin) + +if nargin < 2 + disp(sprintf('Not enough arguments')); + disp(sprintf('Number of arguments: %d', nargin)); + celldisp(varargin) + return +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% pRF_gaussian('getModelResponse', fitParams, rfModel, hrf, p, i) %%%% +%%% Called from getModelResidual %%%% +if strcmp(varargin{1}, 'getModelResponse') + + fitParams = varargin{2}; + rfModel1 = varargin{3}; + hrf = varargin{4}; + p = varargin{5}; + i = varargin{6}; + + rfModel2 = exp(-(((fitParams.stimX - p.x).^2)/(2*((p.std*p.stdRatio)^2))+((fitParams.stimY-p.y).^2)/(2*((p.std*p.stdRatio)^2)))); + nFrames = fitParams.concatInfo.runTransition(i,2); + + % Convolve model with stimulus + rPlus = convolveModelWithStimulus(rfModel1,fitParams.stim{i},nFrames); + rMinus = convolveModelWithStimulus(rfModel2, fitParams.stim{i}, nFrames); + + % Apply divisive normalization + resp = rPlus ./ (1 + p.b1*rMinus); + + % Convolve response with hemodynamic response function + thisModelResponse = convolveModelResponseWithHRF(resp, hrf); + + % drop junk frames here + %thisModelResponse = thisModelResponse(fitParams.concatInfo.totalJunkedFrames(i)+1:end); + if fitParams.concatInfo.isConcat + thisModelResponse = thisModelResponse(fitParams.concatInfo.junkFrames(i)+1:end); + end + % return the calculated model response + output = thisModelResponse; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% pRF_gaussian('setParams', fitParams) %%%% +%%% Called from setParams %%%% + +elseif strcmp(varargin{1}, 'setParams') + + fitParams = varargin{2}; + + fitParams.paramNames = {'x','y','rfWidth', 'surroundRatio', 'surroundGain'}; + fitParams.paramDescriptions = {'RF x position','RF y position','RF width (std of gaussian)', 'Ratio of surround width to center width', 'Surround Gain'}; + fitParams.paramIncDec = [1 1 1 1 1]; + fitParams.paramMin = [-inf -inf 0 0 -inf]; + fitParams.paramMax = [inf inf inf inf -inf]; + % set min/max and init + fitParams.minParams = [fitParams.stimExtents(1) fitParams.stimExtents(2) 0 1 0]; + fitParams.maxParams = [fitParams.stimExtents(3) fitParams.stimExtents(4) inf inf inf]; + fitParams.initParams = [0 0 4 2 1]; + + % return fitParams with modified values + output = fitParams; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%% pRF_gaussian(command, fitParams, params) %%%%% +%%%% Called from getFitParams %%%%% + +elseif strcmp(varargin{1}, 'getFitParams') + + fitParams = varargin{2}; + params = varargin{3}; + + p.rfType = fitParams.rfType; + + % Define your parameters here + p.x = params(1); + p.y = params(2); + p.std = params(3); + p.stdRatio = params(4); + p.b1 = params(5); + p.b2 = params(6); + % use a fixed single gaussian + p.canonical.type = 'gamma'; + p.canonical.lengthInSeconds = 25; + p.canonical.timelag = fitParams.timelag; + p.canonical.tau = fitParams.tau; + p.canonical.exponent = fitParams.exponent; + p.canonical.offset = 0; + p.canonical.diffOfGamma = fitParams.diffOfGamma; + p.canonical.amplitudeRatio = fitParams.amplitudeRatio; + p.canonical.timelag2 = fitParams.timelag2; + p.canonical.tau2 = fitParams.tau2; + p.canonical.exponent2 = fitParams.exponent2; + p.canonical.offset2 = 0; + + output = struct(p); + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelWithStimulus %% +function modelResponse = convolveModelWithStimulus(rfModel,stim,nFrames) + +% get number of frames +nStimFrames = size(stim.im,3); + +% preallocate memory +modelResponse = zeros(1,nStimFrames); + +for frameNum = 1:nStimFrames + % multipy the stimulus frame by frame with the rfModel + % and take the sum + modelResponse(frameNum) = sum(sum(rfModel.*stim.im(:,:,frameNum))); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelResponseWithHRF %% +function modelTimecourse = convolveModelResponseWithHRF(modelTimecourse,hrf) + +n = length(modelTimecourse); +modelTimecourse = conv(modelTimecourse,hrf.hrf); +modelTimecourse = modelTimecourse(1:n); + diff --git a/mrLoadRet/Plugin/pRF/models/pRF_exp.m b/mrLoadRet/Plugin/pRF/models/pRF_exp.m new file mode 100644 index 000000000..a83043dca --- /dev/null +++ b/mrLoadRet/Plugin/pRF/models/pRF_exp.m @@ -0,0 +1,139 @@ +% pRF_exp.m +% +% $Id:$ +% usage: pRF_exp(varargin) +% by: akshay jagadeesh +% date: 09/01/16 +% purpose: Model file to specify a new rftype +% +% +% +% This model template is set up for the standard Gaussian model. + +function output = pRF_exp(varargin) + +if nargin < 2 + disp(sprintf('Not enough arguments')); + disp(sprintf('Number of arguments: %d', nargin)); + celldisp(varargin) + return +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%pRF_exp(command, fitParams, rfModel, hrf, i)%%%%% +%%%% Called from getModelResidual %%%%% +if strcmp(varargin{1}, 'getModelResponse') + + fitParams = varargin{2}; + rfModel = varargin{3}; + hrf = varargin{4}; + p = varargin{5}; + i = varargin{6}; + + nFrames = fitParams.concatInfo.runTransition(i,2); + thisModelResponse = convolveModelWithStimulus(rfModel,fitParams.stim{i},nFrames); + + % FOR GAUSSIAN-EXP ONLY: include exponent non-linearity + if strcmp(fitParams.rfType, 'gaussian-exp') + if ~isnan(thisModelResponse) + thisModelResponse = power(thisModelResponse, p.exp); + end + end + + % and convolve in time. + thisModelResponse = convolveModelResponseWithHRF(thisModelResponse,hrf); + + % drop junk frames here + %thisModelResponse = thisModelResponse(fitParams.concatInfo.totalJunkedFrames(i)+1:end); + if fitParams.concatInfo.isConcat + thisModelResponse = thisModelResponse(fitParams.concatInfo.junkFrames(i)+1:end); + end + % return the calculated model response + output = thisModelResponse; +%end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%pRF_exp(command, fitParams)%%%%% +%%%% Called from setFitParams %%%%% + +elseif strcmp(varargin{1}, 'setParams') + + fitParams = varargin{2}; + + fitParams.paramNames = {'x','y','rfWidth','exp'}; + fitParams.paramDescriptions = {'RF x position','RF y position','RF width (std of gaussian)', 'Spatial summation exponent'}; + fitParams.paramIncDec = [1 1 1 1]; + fitParams.paramMin = [-inf -inf 0 0]; + fitParams.paramMax = [inf inf inf inf]; + % set min/max and init + fitParams.minParams = [fitParams.stimExtents(1) fitParams.stimExtents(2) 0 0]; + fitParams.maxParams = [fitParams.stimExtents(3) fitParams.stimExtents(4) inf 5]; % Set max exponent to 5 + fitParams.initParams = [0 0 4 1]; %Initialize exponent to 1 (linear) + + % return fitParams with modified values + output = struct(fitParams); +%end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%pRFModelTemplate(command, fitParams, params)%%%%% +%%%% Called from getFitParams %%%%% + +elseif strcmp(varargin{1}, 'getFitParams') + + fitParams = varargin{2}; + params = varargin{3}; + + p.rfType = fitParams.rfType; + + % Define your parameters here + p.x = params(1); + p.y = params(2); + p.std = params(3); + p.exp = params(4); + % use a fixed single gaussian + p.canonical.type = 'gamma'; + p.canonical.lengthInSeconds = 25; + p.canonical.timelag = fitParams.timelag; + p.canonical.tau = fitParams.tau; + p.canonical.exponent = fitParams.exponent; + p.canonical.offset = 0; + p.canonical.diffOfGamma = fitParams.diffOfGamma; + p.canonical.amplitudeRatio = fitParams.amplitudeRatio; + p.canonical.timelag2 = fitParams.timelag2; + p.canonical.tau2 = fitParams.tau2; + p.canonical.exponent2 = fitParams.exponent2; + p.canonical.offset2 = 0; + + output = struct(p); + + +else + disp(sprintf('Function did not execute properly')); + return +end + + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelWithStimulus %% +function modelResponse = convolveModelWithStimulus(rfModel,stim,nFrames) + +% get number of frames +nStimFrames = size(stim.im,3); + +% preallocate memory +modelResponse = zeros(1,nStimFrames); + +for frameNum = 1:nStimFrames + % multipy the stimulus frame by frame with the rfModel + % and take the sum + modelResponse(frameNum) = sum(sum(rfModel.*stim.im(:,:,frameNum))); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelResponseWithHRF %% +function modelTimecourse = convolveModelResponseWithHRF(modelTimecourse,hrf) + +n = length(modelTimecourse); +modelTimecourse = conv(modelTimecourse,hrf.hrf); +modelTimecourse = modelTimecourse(1:n); + diff --git a/mrLoadRet/Plugin/pRF/models/pRF_gaussian.m b/mrLoadRet/Plugin/pRF/models/pRF_gaussian.m new file mode 100644 index 000000000..12e3dc7c4 --- /dev/null +++ b/mrLoadRet/Plugin/pRF/models/pRF_gaussian.m @@ -0,0 +1,129 @@ +% pRF_gaussian.m +% +% $Id:$ +% usage: pRF_gaussian(varargin) +% by: akshay jagadeesh +% date: 09/02/16 +% purpose: Template file to create new models +% +% - Allows you to create a new model type +% - Simply add the model type to pRFGUI and create a +% file like this one for your model with the appropriate +% methods filled in. +% +% This model template is set up for the standard Gaussian model. + +function output = pRF_gaussian(varargin) + +if nargin < 2 + disp(sprintf('Not enough arguments')); + disp(sprintf('Number of argumnets: %d', nargin)); + celldisp(varargin) + return +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% pRF_gaussian('getModelResponse', fitParams, rfModel, hrf, p, i) %%%% +%%% Called from getModelResidual %%%% +if strcmp(varargin{1}, 'getModelResponse') + + fitParams = varargin{2}; + rfModel = varargin{3}; + hrf = varargin{4}; + p = varargin{5}; + i = varargin{6}; + + nFrames = fitParams.concatInfo.runTransition(i,2); + thisModelResponse = convolveModelWithStimulus(rfModel,fitParams.stim{i},nFrames); + + % and convolve in time. + thisModelResponse = convolveModelResponseWithHRF(thisModelResponse,hrf); + + % drop junk frames here + if fitParams.concatInfo.isConcat + thisModelResponse = thisModelResponse(fitParams.concatInfo.junkFrames(i)+1:end); + end + %%%% Maybe uncomment this later: + %thisModelResponse = thisModelResponse(fitParams.concatInfo.totalJunkedFrames(i)+1:end); + + % return the calculated model response + output = thisModelResponse; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%% pRF_gaussian('setParams', fitParams) %%%% +%%% Called from setParams %%%% + +elseif strcmp(varargin{1}, 'setParams') + + fitParams = varargin{2}; + + fitParams.paramNames = {'x','y','rfWidth'}; + fitParams.paramDescriptions = {'RF x position','RF y position','RF width (std of gaussian)'}; + fitParams.paramIncDec = [1 1 1]; + fitParams.paramMin = [-inf -inf 0]; + fitParams.paramMax = [inf inf inf]; + % set min/max and init + fitParams.minParams = [fitParams.stimExtents(1) fitParams.stimExtents(2) 0]; + fitParams.maxParams = [fitParams.stimExtents(3) fitParams.stimExtents(4) inf]; + fitParams.initParams = [0 0 4]; + + % return fitParams with modified values + output = fitParams; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%% pRF_gaussian(command, fitParams, params) %%%%% +%%%% Called from getFitParams %%%%% + +elseif strcmp(varargin{1}, 'getFitParams') + + fitParams = varargin{2}; + params = varargin{3}; + + p.rfType = fitParams.rfType; + + % Define your parameters here + p.x = params(1); + p.y = params(2); + p.std = params(3); + % use a fixed single gaussian + p.canonical.type = 'gamma'; + p.canonical.lengthInSeconds = 25; + p.canonical.timelag = fitParams.timelag; + p.canonical.tau = fitParams.tau; + p.canonical.exponent = fitParams.exponent; + p.canonical.offset = 0; + p.canonical.diffOfGamma = fitParams.diffOfGamma; + p.canonical.amplitudeRatio = fitParams.amplitudeRatio; + p.canonical.timelag2 = fitParams.timelag2; + p.canonical.tau2 = fitParams.tau2; + p.canonical.exponent2 = fitParams.exponent2; + p.canonical.offset2 = 0; + + output = struct(p); + +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelWithStimulus %% +function modelResponse = convolveModelWithStimulus(rfModel,stim,nFrames) + +% get number of frames +nStimFrames = size(stim.im,3); + +% preallocate memory +modelResponse = zeros(1,nStimFrames); + +for frameNum = 1:nStimFrames + % multipy the stimulus frame by frame with the rfModel + % and take the sum + modelResponse(frameNum) = sum(sum(rfModel.*stim.im(:,:,frameNum))); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelResponseWithHRF %% +function modelTimecourse = convolveModelResponseWithHRF(modelTimecourse,hrf) + +n = length(modelTimecourse); +modelTimecourse = conv(modelTimecourse,hrf.hrf); +modelTimecourse = modelTimecourse(1:n); + diff --git a/mrLoadRet/Plugin/pRF/models/pRF_gaussianhdr.m b/mrLoadRet/Plugin/pRF/models/pRF_gaussianhdr.m new file mode 100644 index 000000000..3fbb58366 --- /dev/null +++ b/mrLoadRet/Plugin/pRF/models/pRF_gaussianhdr.m @@ -0,0 +1,136 @@ +% pRF_gaussian.m +% +% $Id:$ +% usage: pRF_gaussianhdr(varargin) +% by: akshay jagadeesh +% date: 09/02/16 +% purpose: Template file to create new models +% +% - Allows you to create a new model type +% - Simply add the model type to pRFGUI and create a +% file like this one for your model with the appropriate +% methods filled in. +% +% This model template is set up for the standard Gaussian model. + +function output = pRF_gaussianhdr(varargin) + +if nargin < 2 + disp(sprintf('Not enough arguments')); + return +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%pRF_gaussianhdr('getModelResponse', fitParams, rfModel, hrf, i)%%%%% +%%%% Called from getModelResidual %%%%% +if strcmp(varargin{1}, 'getModelResponse') + + fitParams = varargin{2}; + rfModel = varargin{3}; + hrf = varargin{4}; + p = varargin{5}; + i = varargin{6}; + + nFrames = fitParams.concatInfo.runTransition(i,2); + thisModelResponse = convolveModelWithStimulus(rfModel,fitParams.stim{i},nFrames); + + % and convolve in time. + thisModelResponse = convolveModelResponseWithHRF(thisModelResponse,hrf); + + % drop junk frames here + thisModelResponse = thisModelResponse(fitParams.concatInfo.totalJunkedFrames(i)+1:end); + + % return the calculated model response + output = thisModelResponse; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%%%%pRF_gaussian('setParams', fitParams)%%%%% +%%%% Called from setParams %%%%% + +elseif strcmp(varargin{1}, 'setParams') + + fitParams = varargin{2}; + + fitParams.paramNames = {'x','y','rfWidth', 'timelag', 'tau'}; + fitParams.paramDescriptions = {'RF x position','RF y position','RF width (std of gaussian)', 'Time before start of rise of hemodynamic function', 'Width of the hemodynamic function (tau parameter of gamma)'}; + fitParams.paramIncDec = [1 1 1 0.1 0.5]; + fitParams.paramMin = [-inf -inf 0 0 0]; + fitParams.paramMax = [inf inf inf inf inf]; + % set min/max and init + fitParams.minParams = [fitParams.stimExtents(1) fitParams.stimExtents(2) 0 0 0]; + fitParams.maxParams = [fitParams.stimExtents(3) fitParams.stimExtents(4) inf 3 inf]; + fitParams.initParams = [0 0 4 fitParams.timelag fitParams.tau]; + % add on parameters for difference of gamma + if fitParams.diffOfGamma + % parameter names/descriptions and other information for allowing user to set them + fitParams.paramNames = {fitParams.paramNames{:} 'amp2' 'timelag2','tau2'}; + fitParams.paramDescriptions = {fitParams.paramDescriptions{:} 'Amplitude of second gamma for HDR' 'Timelag for second gamma for HDR','tau for second gamma for HDR'}; + fitParams.paramIncDec = [fitParams.paramIncDec(:)' 0.1 0.1 0.5]; + fitParams.paramMin = [fitParams.paramMin(:)' 0 0 0]; + fitParams.paramMax = [fitParams.paramMax(:)' inf inf inf]; + % set min/max and init + fitParams.minParams = [fitParams.minParams 0 0 0]; + fitParams.maxParams = [fitParams.maxParams inf 6 inf]; + fitParams.initParams = [fitParams.initParams fitParams.amplitudeRatio fitParams.timelag2 fitParams.tau2]; + end + + % return fitParams with modified values + output = fitParams; + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%0%%%%% +%%%% pRF_gaussianhdr(command, fitParams, params)%%%%% +%%%% Called from getFitParams %%%%% + +elseif strcmp(varargin{1}, 'getFitParams') + + fitParams = varargin{2}; + params = varargin{3}; + + p.rfType = fitParams.rfType; + + % Define your parameters here + p.x = params(1); + p.y = params(2); + p.std = params(3); + % use a fixed single gaussian + p.canonical.type = 'gamma'; + p.canonical.lengthInSeconds = 25; + p.canonical.timelag = params(4); + p.canonical.tau = params(5); + p.canonical.exponent = fitParams.exponent; + p.canonical.offset = 0; + p.canonical.diffOfGamma = fitParams.diffOfGamma; + if fitParams.diffOfGamma + p.canonical.amplitudeRatio = params(6); + p.canonical.timelag2 = params(7); + p.canonical.tau2 = params(8); + p.canonical.exponent2 = fitParams.exponent2; + p.canonical.offset2 = 0; + end + output = struct(p); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelWithStimulus %% +function modelResponse = convolveModelWithStimulus(rfModel,stim,nFrames) + +% get number of frames +nStimFrames = size(stim.im,3); + +% preallocate memory +modelResponse = zeros(1,nFrames); + +for frameNum = 1:nStimFrames + % multipy the stimulus frame by frame with the rfModel + % and take the sum + modelResponse(frameNum) = sum(sum(rfModel.*stim.im(:,:,frameNum))); +end + +%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% +%% convolveModelResponseWithHRF %% +function modelTimecourse = convolveModelResponseWithHRF(modelTimecourse,hrf) + +n = length(modelTimecourse); +modelTimecourse = conv(modelTimecourse,hrf.hrf); +modelTimecourse = modelTimecourse(1:n); + diff --git a/mrLoadRet/Plugin/pRF/pRFFit.m b/mrLoadRet/Plugin/pRF/pRFFit.m index e81a57df7..4d9899f12 100644 --- a/mrLoadRet/Plugin/pRF/pRFFit.m +++ b/mrLoadRet/Plugin/pRF/pRFFit.m @@ -54,7 +54,6 @@ fit = fitParams.stim; return end - % get the tSeries if ~isempty(x) % if tSeries was not passed in then load it @@ -67,7 +66,7 @@ % but useful for raw/motionCorrected time series. Also, it is very important that % the tSeries is properly mean subtracted if ~isfield(fitParams.concatInfo,'hipassfilter') - tSeries = percentTSeries(tSeries,'detrend','Linear','spatialNormalization','Divide by mean','subtractMean', 'Yes', 'temporalNormalization', 'No'); + tSeries = percentTSeries(tSeries,'detrend','Linear', 'spatialNormalization','Divide by mean','subtractMean', 'Yes', 'temporalNormalization', 'No'); end % if there are any nans in the tSeries then don't fit @@ -80,7 +79,6 @@ else tSeries = []; end - % handle junk frames (i.e. ones that have not already been junked) if ~isempty(fitParams.junkFrames) && ~isequal(fitParams.junkFrames,0) % drop junk frames @@ -143,7 +141,9 @@ if ~isfield(fitParams.prefit,'modelResponse') % get number of workers nProcessors = mlrNumWorkers; - disppercent(-inf,sprintf('(pRFFit) Computing %i prefit model responses using %i processors',fitParams.prefit.n,nProcessors)); + if fitParams.verbose==1 + disppercent(-inf,sprintf('(pRFFit) Computing %i prefit model responses using %i processors',fitParams.prefit.n,nProcessors)); + end % first convert the x/y and width parameters into sizes % on the actual screen fitParams.prefit.x = fitParams.prefit.x*fitParams.stimWidth; @@ -184,10 +184,13 @@ if fitParams.prefitOnly % return if we are just doing a prefit fit = getFitParams(fitParams.initParams,fitParams); + fit.modelResponse = fitParams.prefit.modelResponse; + fit.tSeries = tSeries; fit.rfType = fitParams.rfType; fit.params = fitParams.initParams; fit.r2 = maxr^2; fit.r = maxr; + fit.bestFitVoxel = bestModel; [fit.polarAngle fit.eccentricity] = cart2pol(fit.x,fit.y); % display if fitParams.verbose @@ -202,6 +205,8 @@ [params resnorm residual exitflag output lambda jacobian] = lsqnonlin(@getModelResidual,fitParams.initParams,fitParams.minParams,fitParams.maxParams,fitParams.optimParams,tSeries,fitParams); elseif strcmp(lower(fitParams.algorithm),'nelder-mead') [params fval exitflag] = fminsearch(@getModelResidual,fitParams.initParams,fitParams.optimParams,(tSeries-mean(tSeries))/var(tSeries.^2),fitParams); +elseif strcmp(lower(fitParams.algorithm), 'nelder-mead-bnd') + [params fval exitflag] = fminsearchbnd(@getModelResidual, fitParams.initParams, fitParams.minParams, fitParams.maxParams, fitParams.optimParams, (tSeries-mean(tSeries))/var(tSeries.^2), fitParams); else disp(sprintf('(pRFFit) Unknown optimization algorithm: %s',fitParams.algorithm)); return @@ -218,15 +223,20 @@ fit.r2 = 1-sum((residual-mean(residual)).^2)/sum((tSeries-mean(tSeries)).^2); elseif strcmp(lower(fitParams.algorithm),'nelder-mead') fit.r2 = residual^2; +elseif strcmp(lower(fitParams.algorithm), 'nelder-mead-bnd') + fit.r2 = residual^2; end % compute polar coordinates [fit.polarAngle fit.eccentricity] = cart2pol(fit.x,fit.y); % display -if fitParams.verbose - disp(sprintf('%s[%2.f %2.f %2.f] r2=%0.2f polarAngle=%6.1f eccentricity=%6.1f rfHalfWidth=%6.1f',fitParams.dispstr,x,y,z,fit.r2,r2d(fit.polarAngle),fit.eccentricity,fit.std)); +if fitParams.verbose && strcmp(fitParams.rfType, 'gaussian-exp') + disp(sprintf('%s[%2.f %2.f %2.f] r2=%0.2f polarAngle=%6.1f eccentricity=%6.1f rfHalfWidth=%6.1f exp=%f',fitParams.dispstr,x,y,z,fit.r2,r2d(fit.polarAngle),fit.eccentricity,fit.std, fit.exp)); +elseif fitParams.verbose + disp(sprintf('%s[%2.f %2.f %2.f] r2=%0.2f polarAngle=%6.1f eccentricity=%6.1f rfHalfWidth=%6.1f', fitParams.dispstr,x,y,z,fit.r2,r2d(fit.polarAngle),fit.eccentricity,fit.std)); end +%keyboard %%%%%%%%%%%%%%%%%%%%%% % setFitParams % @@ -253,46 +263,10 @@ if ~isfield(fitParams,'initParams') % check the rfType to get the correct min/max arrays - switch (fitParams.rfType) - case 'gaussian' - % parameter names/descriptions and other information for allowing user to set them - fitParams.paramNames = {'x','y','rfWidth'}; - fitParams.paramDescriptions = {'RF x position','RF y position','RF width (std of gaussian)'}; - fitParams.paramIncDec = [1 1 1]; - fitParams.paramMin = [-inf -inf 0]; - fitParams.paramMax = [inf inf inf]; - % set min/max and init - fitParams.minParams = [fitParams.stimExtents(1) fitParams.stimExtents(2) 0]; - fitParams.maxParams = [fitParams.stimExtents(3) fitParams.stimExtents(4) inf]; - fitParams.initParams = [0 0 4]; - case 'gaussian-hdr' - % parameter names/descriptions and other information for allowing user to set them - fitParams.paramNames = {'x','y','rfWidth','timelag','tau'}; - fitParams.paramDescriptions = {'RF x position','RF y position','RF width (std of gaussian)','Time before start of rise of hemodynamic function','Width of the hemodynamic function (tau parameter of gamma)'}; - fitParams.paramIncDec = [1 1 1 0.1 0.5]; - fitParams.paramMin = [-inf -inf 0 0 0]; - fitParams.paramMax = [inf inf inf inf inf]; - % set min/max and init - fitParams.minParams = [fitParams.stimExtents(1) fitParams.stimExtents(2) 0 0 0]; - fitParams.maxParams = [fitParams.stimExtents(3) fitParams.stimExtents(4) inf 3 inf]; - fitParams.initParams = [0 0 4 fitParams.timelag fitParams.tau]; - % add on parameters for difference of gamma - if fitParams.diffOfGamma - % parameter names/descriptions and other information for allowing user to set them - fitParams.paramNames = {fitParams.paramNames{:} 'amp2' 'timelag2','tau2'}; - fitParams.paramDescriptions = {fitParams.paramDescriptions{:} 'Amplitude of second gamma for HDR' 'Timelag for second gamma for HDR','tau for second gamma for HDR'}; - fitParams.paramIncDec = [fitParams.paramIncDec(:)' 0.1 0.1 0.5]; - fitParams.paramMin = [fitParams.paramMin(:)' 0 0 0]; - fitParams.paramMax = [fitParams.paramMax(:)' inf inf inf]; - % set min/max and init - fitParams.minParams = [fitParams.minParams 0 0 0]; - fitParams.maxParams = [fitParams.maxParams inf 6 inf]; - fitParams.initParams = [fitParams.initParams fitParams.amplitudeRatio fitParams.timelag2 fitParams.tau2]; - end - otherwise - disp(sprintf('(pRFFit:setFitParams) Unknown rfType %s',rfType)); - return - end + + %%%%%%%%%%%%%%%%%%% + fitParams = prfModel('setParams', fitParams); + %%%%%%%%%%%%%%%%%%% % round constraints fitParams.minParams = round(fitParams.minParams*10)/10; @@ -300,7 +274,7 @@ % handle constraints here % Check if fit algorithm is one that allows constraints - algorithmsWithConstraints = {'levenberg-marquardt'}; + algorithmsWithConstraints = {'levenberg-marquardt', 'nelder-mead-bnd'}; if any(strcmp(fitParams.algorithm,algorithmsWithConstraints)) % if constraints allowed then allow user to adjust them here (if they set defaultConstraints) if isfield(fitParams,'defaultConstraints') && ~fitParams.defaultConstraints @@ -328,7 +302,9 @@ end else % no constraints allowed - disp(sprintf('(pRFFit) !!! Fit constraints ignored for algorithm: %s (if you want to constrain the fits, then use: %s) !!!',fitParams.algorithm,cell2mat(algorithmsWithConstraints))); + if fitParams.verbose==1 + disp(sprintf('(pRFFit) !!! Fit constraints ignored for algorithm: %s (if you want to constrain the fits, then use: %s) !!!',fitParams.algorithm,cell2mat(algorithmsWithConstraints))); + end end end @@ -336,7 +312,8 @@ % optimization parameters if ~isfield(fitParams,'algorithm') || isempty(fitParams.algorithm) - fitParams.algorithm = 'nelder-mead'; + fitParams.algorithm = 'nelder-mead-bnd'; + disp('(pRFFit) No algorithm provided. Using Default: Nelder-Mead-Bnd'); end fitParams.optimParams = optimset('MaxIter',inf,'Display',fitParams.optimDisplay); @@ -371,18 +348,18 @@ % create the model for each concat for i = 1:fitParams.concatInfo.n - % get model response - nFrames = fitParams.concatInfo.runTransition(i,2); - thisModelResponse = convolveModelWithStimulus(rfModel,fitParams.stim{i},nFrames); - % get a model hrf hrf = getCanonicalHRF(p.canonical,fitParams.framePeriod); - % and convolve in time. - thisModelResponse = convolveModelResponseWithHRF(thisModelResponse,hrf); + %%%%%%%%%%%%%%%%%%% + % Get model response, which involves convolving model with stimulus, and with HRF and dropping junk frames. + thisModelResponse = prfModel('getModelResponse', fitParams, rfModel, hrf, p, i); + %%%%%%%%%%%%%%%%%%% + + %~~~~~~~~~~ + %keyboard + %~~~~~~~~~~ - % drop junk frames here - thisModelResponse = thisModelResponse(fitParams.concatInfo.totalJunkedFrames(i)+1:end); % apply concat filtering if isfield(fitParams,'applyFiltering') && fitParams.applyFiltering @@ -435,7 +412,7 @@ end % for nelder-mead just compute correlation and return 1-4 -if strcmp(lower(fitParams.algorithm),'nelder-mead') +if strcmp(lower(fitParams.algorithm),'nelder-mead') || strcmp(lower(fitParams.algorithm), 'nelder-mead-bnd') residual = -corr(modelResponse,tSeries); % disp(sprintf('(pRFFit:getModelResidual) r: %f',residual)); end @@ -486,7 +463,7 @@ function dispModelFit(params,fitParams,modelResponse,tSeries,rfModel) designMatrix(:,2) = 1; % get beta weight for the modelResponse -if ~any(isnan(modelResponse)) +if ~any(isnan(modelResponse)) && ~any(isinf(modelResponse)) beta = pinv(designMatrix)*tSeries; beta(1) = max(beta(1),0); modelResponse = designMatrix*beta; @@ -501,47 +478,7 @@ function dispModelFit(params,fitParams,modelResponse,tSeries,rfModel) function p = getFitParams(params,fitParams) p.rfType = fitParams.rfType; - -switch (fitParams.rfType) - case 'gaussian' - p.x = params(1); - p.y = params(2); - p.std = params(3); - % use a fixed single gaussian - p.canonical.type = 'gamma'; - p.canonical.lengthInSeconds = 25; - p.canonical.timelag = fitParams.timelag; - p.canonical.tau = fitParams.tau; - p.canonical.exponent = fitParams.exponent; - p.canonical.offset = 0; - p.canonical.diffOfGamma = fitParams.diffOfGamma; - p.canonical.amplitudeRatio = fitParams.amplitudeRatio; - p.canonical.timelag2 = fitParams.timelag2; - p.canonical.tau2 = fitParams.tau2; - p.canonical.exponent2 = fitParams.exponent2; - p.canonical.offset2 = 0; - case 'gaussian-hdr' - p.x = params(1); - p.y = params(2); - p.std = params(3); - % use a fixed single gaussian - p.canonical.type = 'gamma'; - p.canonical.lengthInSeconds = 25; - p.canonical.timelag = params(4); - p.canonical.tau = params(5); - p.canonical.exponent = fitParams.exponent; - p.canonical.offset = 0; - p.canonical.diffOfGamma = fitParams.diffOfGamma; - if fitParams.diffOfGamma - p.canonical.amplitudeRatio = params(6); - p.canonical.timelag2 = params(7); - p.canonical.tau2 = params(8); - p.canonical.exponent2 = fitParams.exponent2; - p.canonical.offset2 = 0; - end -otherwise - disp(sprintf('(pRFFit) Unknown rfType %s',rfType)); -end +p = prfModel('getFitParams', fitParams, params); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% convolveModelWithStimulus %% @@ -552,7 +489,7 @@ function dispModelFit(params,fitParams,modelResponse,tSeries,rfModel) nStimFrames = size(stim.im,3); % preallocate memory -modelResponse = zeros(1,nFrames); +modelResponse = zeros(1,nStimFrames); for frameNum = 1:nStimFrames % multipy the stimulus frame by frame with the rfModel @@ -619,12 +556,34 @@ function dispModelFit(params,fitParams,modelResponse,tSeries,rfModel) rfModel = []; % now gernerate the rfModel -if any(strcmp(fitParams.rfType,{'gaussian','gaussian-hdr'})) +if any(strcmp(fitParams.rfType,{'gaussian','gaussian-hdr', 'gaussian-exp', 'gaussian-diffs', 'gaussian-divs', 'gaussian-DoG-CSS'})) rfModel = makeRFGaussian(params,fitParams); else disp(sprintf('(pRFFit:getRFModel) Unknown rfType: %s',fitParams.rfType)); end +%%%%%%%%%%%%%%%%%%% +%% prfModel %% +%%%%%%%%%%%%%%%%%%% +function output = prfModel(varargin) + +fitParams = varargin{2}; +switch (fitParams.rfType) + case 'gaussian' + output = pRF_gaussian(varargin{:}); + case 'gaussian-hdr' + output = pRF_gaussianhdr(varargin{:}); + case 'gaussian-diffs' + output = pRF_diffGaussian(varargin{:}); + case 'gaussian-divs' + output = pRF_divGaussian(varargin{:}); + case 'gaussian-exp' + output = pRF_exp(varargin{:}); + otherwise % 'gaussian-DoG-CSS' + testModel = @pRF_DoG_CSS; %%% Only need to change this line to specify a new model. + output = testModel(varargin{:}); +end + %%%%%%%%%%%%%%%%%%%%%%%% %% makeRFGaussian %% diff --git a/mrLoadRet/Plugin/pRF/pRFGUI.m b/mrLoadRet/Plugin/pRF/pRFGUI.m index 9d8d101e2..703fc162e 100644 --- a/mrLoadRet/Plugin/pRF/pRFGUI.m +++ b/mrLoadRet/Plugin/pRF/pRFGUI.m @@ -116,9 +116,9 @@ end %all of these parameters are for pRFFit -paramsInfo{end+1} = {'rfType',{'gaussian','gaussian-hdr'},'Type of pRF fit. Gaussian fits a gaussian with x,y,width as parameters to each voxel. gaussian-hdr fits also the hemodynamic response with the parameters of the hdr as below.'}; +paramsInfo{end+1} = {'rfType',{'gaussian','gaussian-hdr', 'gaussian-exp', 'gaussian-diffs', 'gaussian-divs', 'gaussian-DoG-CSS'},'Type of pRF fit. Gaussian fits a gaussian with x,y,width as parameters to each voxel. gaussian-hdr fits also the hemodynamic response with the parameters of the hdr as below. gaussian-exp incorporates an exponent non-linearity when calculating neuron response, gaussian-diffs calculates response as the difference between a center and a surround gaussians. gaussian-divs takes the quotient of center and surround gaussians.'}; paramsInfo{end+1} = {'betaEachScan',false,'type=checkbox','Compute a separate beta weight (scaling) for each scan in the concanetation. This may be useful if there is some reason to believe that different scans have different magnitude responses, this will allow the fit to scale the magnitude for each scan'}; -paramsInfo{end+1} = {'algorithm',{'nelder-mead','levenberg-marquardt'},'Which algorithm to use for optimization. Levenberg-marquardt seems to get stuck in local minimum, so the default is nelder-mead. However, levenberg-marquardt can set bounds for parameters, so may be better for when you are trying to fit the hdr along with the rf, since the hdr parameters can fly off to strange values.'}; +paramsInfo{end+1} = {'algorithm',{'nelder-mead','levenberg-marquardt', 'nelder-mead-bnd'},'Which algorithm to use for optimization. Levenberg-marquardt seems to get stuck in local minimum, so the default is nelder-mead. However, levenberg-marquardt can set bounds for parameters, so may be better for when you are trying to fit the hdr along with the rf, since the hdr parameters can fly off to strange values.'}; paramsInfo{end+1} = {'defaultConstraints',1,'type=checkbox','Sets how to constrain the search (i.e. what are the allowed range of stimulus parameters). The default is to constrain so that the x,y of the RF has to be within the stimulus extents (other parameter constrains will print to the matlab window). If you click this off a dialog box will come up after the stimulus has been calculated from the stimfiles allowing you to specify the constraints on the parameters of the model. You may want to custom constrain the parameters if you know something about the RFs you are trying to model (like how big they are) to keep the nonlinear fits from finding unlikely parameter estimates. Note that nelder-mead is an unconstrained fit so this will not do anything.'}; paramsInfo{end+1} = {'prefitOnly',false,'type=checkbox','Check this if you want to ONLY do a prefit and not optimize further. The prefit computes a preset set of model parameters (x,y,rfHalfWidth) and picks the one that produces a mdoel with the highest correlation with the time series. You may want to do this to get a quick but accurate fit so that you can draw a set of ROIs for a full analysis'}; paramsInfo{end+1} = {'quickPrefit',false,'type=checkbox','Check this if you want to do a quick prefit - this samples fewer x,y and rfWidth points. It is faster (especially if coupled with prefitOnly for a fast check), but the optimization routines may be more likely to get trapped into local minima or have to search a long time for the minimum'}; diff --git a/mrLoadRet/Plugin/pRF/pRFGetStimImageFromStimfile.m b/mrLoadRet/Plugin/pRF/pRFGetStimImageFromStimfile.m index bba7b545b..8a7b0528e 100644 --- a/mrLoadRet/Plugin/pRF/pRFGetStimImageFromStimfile.m +++ b/mrLoadRet/Plugin/pRF/pRFGetStimImageFromStimfile.m @@ -451,7 +451,7 @@ % have a task which is mglRetinotopy taskNum = []; for iTask = 1:2 - if (length(thiss.task) >= iTask) && (isequal(thiss.task{iTask}{1}.taskFilename,'mglRetinotopy.m') || isequal(thiss.task{iTask}{1}.taskFilename,'gruRetinotopy.m')) + if (length(thiss.task) >= iTask) && (isequal(thiss.task{iTask}{1}.taskFilename,'mglRetinotopy.m') || isequal(thiss.task{iTask}{1}.taskFilename,'gruRetinotopy.m') || isequal(thiss.task{iTask}{1}.taskFilename, 'mglDoubleBars.m')) taskNum = iTask; end end @@ -467,7 +467,7 @@ if ~isfield(thiss.task{taskNum}{1},'randVars') missing = 'randVars';end if ~isfield(thiss.task{taskNum}{1},'parameter') missing = 'parameter';end if ~any(strcmp('maskPhase',thiss.myscreen.traceNames)) missing = 'maskPhase';end - if ~any(strcmp('blank',thiss.myscreen.traceNames)) missing = 'blank';end + if ~any(strcmp('blank',thiss.myscreen.traceNames)) && isequal(thiss.task{iTask}{1}.taskFilename,'mglRetinotopy.m') missing = 'blank';end if ~isempty(missing) disp(sprintf('(pRFGetStimImageFromStimfile:checkStimfile) Stimfile: %s',dispstr)); disp(sprintf('(pRFGetStimImageFromStimfile:checkStimfile) The stimfile does not appear to have been created by the latest version of mglRetinotopy which contains the field %s necessary for reconstructing the stimulus. Consider running a dummy run with a newer version of mglRetinotpy with the same parameters (see mglSimulateRun to simulate backticks) and then use that stimfile instead of this one.',missing)); @@ -480,6 +480,9 @@ % now check for each variable that we need varnames = {'blank'}; + if isequal(thiss.task{taskNum}{1}.taskFilename, 'mglDoubleBars.m') % mglDoubleBars doesn't use blanks + varnames = {}; + end for i = 1:length(varnames) varval = getVarFromParameters(varnames{i},e); if isempty(varval) @@ -494,12 +497,15 @@ if ~isempty(stimulusType) && (stimulusType ~= thiss.stimulusType) disp(sprintf('(pRFGetStimImageFromStimfile:checkStimfile) !!! Stimfile %s does not match previous one !!! Have you averaged together scans with different stimulus conditions?')); end - if any(thiss.stimulus.stimulusType == [3 4]) + if ~isequal(thiss.task{taskNum}{1}.taskFilename, 'mglDoubleBars.m') && any(thiss.stimulus.stimulusType == [3 4]) varval = getVarFromParameters('barAngle',e); if ~isempty(barAngle) && ~isequal(varval,barAngle) disp(sprintf('(pRFGetStimImageFromStimfile:checkStimfile) !!! Stimfile %s does not match previous one !!! The barAngles are different! Have you averaged together scans with different stimulus conditions?')); end barAngle = varval; + elseif isequal(thiss.task{taskNum}{1}.taskFilename, 'mglDoubleBars.m') + varval = getVarFromParameters('conditionNum', e); + barAngle = thiss.stimulus.conditions(varval, 3:4); else if ~isempty(direction) && (thiss.stimulus.direction ~= direction) disp(sprintf('(pRFGetStimImageFromStimfile:checkStimfile) !!! Stimfile %s does not match previous one !!! The directions are different! Have you averaged together scans with different stimulus conditions?'));