diff --git a/+basisFactory/deltaStim.m b/+basisFactory/deltaStim.m index f8669e5..8c25e98 100644 --- a/+basisFactory/deltaStim.m +++ b/+basisFactory/deltaStim.m @@ -1,6 +1,7 @@ function stim = deltaStim(bt, nT, v) % Returns a sparse vector with events at binned timings +bt(bt>nT) = []; o = ones(numel(bt), 1); if nargin < 3 diff --git a/+basisFactory/makeNonlinearRaisedCos.m b/+basisFactory/makeNonlinearRaisedCos.m index 4e76019..16982f0 100644 --- a/+basisFactory/makeNonlinearRaisedCos.m +++ b/+basisFactory/makeNonlinearRaisedCos.m @@ -1,4 +1,4 @@ -function bases = makeNonlinearRaisedCos(nBases, binSize, endPoints, nlOffset) +function bases = makeNonlinearRaisedCos(nBases, binSize, endPoints, nlOffset, varargin) % Make nonlinearly stretched basis consisting of raised cosines. % Nonlinear stretching allows faster changes near the event. % @@ -8,6 +8,10 @@ % (i.e. center) of the last raised cosine basis vectors % nlOffset: [1] offset for nonlinear stretching of x axis: y = log(t+nlOffset) % (larger nlOffset -> more nearly linear stretching) + +% Optional Arguments (entered as argument pairs) +% 'Normalize' (default = false) normalizes basis vectors to sum to 1 +% 'Orthogonalize' (default = false) orthogonalizes the basis so B'*B=1 % % Outputs: iht = time lattice on which basis is defined % ihbasis = basis itself @@ -16,6 +20,11 @@ % Example call % bases = basisFactory.makeNonlinearRaisedCos(10, 1, [0 500], 2); +p = inputParser(); +p.addOptional('Normalize', false); +p.addOptional('Orthogonalize', false); +p.parse(varargin{:}); + % nonlinearity for stretching x axis (and its inverse) nlin = @(x)(log(x + 1e-20)); invnl = @(x)(exp(x) - 1e-20); @@ -33,6 +42,14 @@ ihbasis = ff(repmat(nlin(iht + nlOffset), 1, nBases), repmat(ctrs, numel(iht), 1), db); ihctrs = invnl(ctrs); +if p.Results.Normalize + ihbasis = bsxfun(@rdivide, ihbasis, sum(ihbasis)); +end + +if p.Results.Orthogonalize + ihbasis = orth(ihbasis); +end + bases.type = mfilename; bases.param.nBases = nBases; bases.param.binSize = binSize; @@ -40,5 +57,5 @@ bases.param.nlOffset = nlOffset; bases.B = ihbasis; bases.edim = size(bases.B, 2); -bases.tr = iht; +bases.tr = iht/binSize; bases.centers = ihctrs; \ No newline at end of file diff --git a/+basisFactory/makeSmoothTemporalBasis.m b/+basisFactory/makeSmoothTemporalBasis.m index 25cf952..aad22d4 100644 --- a/+basisFactory/makeSmoothTemporalBasis.m +++ b/+basisFactory/makeSmoothTemporalBasis.m @@ -1,47 +1,64 @@ -function bases = makeSmoothTemporalBasis(shape, duration, nBases, binfun) +function bases = makeSmoothTemporalBasis(shape, duration, nBases, binfun, varargin) % % Input % shape: 'raised cosine' or 'boxcar' % duration: the time that needs to be covered % nBases: number of basis vectors to fill the duration -% binfun: +% binfun: +% +% Optional Arguments (entered as argument pairs) +% 'Normalize' (default = false) normalizes basis vectors to sum to 1 +% 'Orthogonalize' (default = false) orthogonalizes the basis so B'*B=1 % % Output % BBstm: basis vectors +p = inputParser(); +p.addOptional('Normalize', false); +p.addOptional('Orthogonalize', false); +p.parse(varargin{:}); + nkbins = binfun(duration); % number of bins for the basis functions ttb = repmat((1:nkbins)', 1, nBases); % time indices for basis -if strcmpi(shape, 'raised cosine') - % ^ - % / \ - % / \______ - % ^ - % / \ - % ___/ \___ - % ^ - % / \ - % ______/ \ - % For raised cosine, the spacing between the centers must be 1/4 of the - % width of the cosine - dbcenter = nkbins / (3 + nBases); % spacing between bumps - width = 4 * dbcenter; % width of each bump - bcenters = 2 * dbcenter + dbcenter*(0:nBases-1); - % location of each bump centers - % bcenters = round(bcenters); % round to closest ms <-- BAD!!! NEVER DO THIS - bfun = @(x,period)((abs(x/period)<0.5).*(cos(x*2*pi/period)*.5+.5)); - BBstm = bfun(ttb-repmat(bcenters,nkbins,1), width); -elseif strcmpi(shape, 'boxcar') - width = nkbins / nBases; - BBstm = zeros(size(ttb)); - bcenters = width * (1:nBases) - width/2; - for k = 1:nBases - idx = ttb(:, k) > ceil(width * (k-1)) & ttb(:, k) <= ceil(width * k); - BBstm(idx, k) = 1 / sum(idx); - end -else - error('Unknown basis shape'); +switch shape + case 'raised cosine' + % ^ + % / \ + % / \______ + % ^ + % / \ + % ___/ \___ + % ^ + % / \ + % ______/ \ + % For raised cosine, the spacing between the centers must be 1/4 of the + % width of the cosine + dbcenter = nkbins / (3 + nBases); % spacing between bumps + width = 4 * dbcenter; % width of each bump + bcenters = 2 * dbcenter + dbcenter*(0:nBases-1); + % location of each bump centers + bfun = @(x,period)((abs(x/period)<0.5).*(cos(x*2*pi/period)*.5+.5)); + BBstm = bfun(ttb-repmat(bcenters,nkbins,1), width); + case 'boxcar' + width = nkbins / nBases; + BBstm = zeros(size(ttb)); + bcenters = width * (1:nBases) - width/2; + for k = 1:nBases + idx = ttb(:, k) > ceil(width * (k-1)) & ttb(:, k) <= ceil(width * k); + BBstm(idx, k) = 1 / sum(idx); + end + otherwise + error('Unknown basis shape'); +end + +if p.Results.Normalize + BBstm = bsxfun(@rdivide, BBstm, sum(BBstm)); +end + +if p.Results.Orthogonalize + BBstm = orth(BBstm); end bases.type = [shape '@' mfilename]; diff --git a/+buildGLM/addBiasColumn.m b/+buildGLM/addBiasColumn.m index d391761..7ec9192 100644 --- a/+buildGLM/addBiasColumn.m +++ b/+buildGLM/addBiasColumn.m @@ -1,10 +1,14 @@ -function dm = addBiasColumn(dm) -% Add a column of ones as the first column to estimate the bias (DC term) +function dm = addBiasColumn(dm, flag) +% Add a column of ones as the first (or last) column to estimate the bias +% (DC term) % dm = addBiasColumn(dm); % % Some regression packages do not allow separate bias estimation, and % a constant column of ones must be added to the design matrix. The weight % that corresponds to this column will be the bias later. +% +% To put the bias column on the right of the design matrix, call: +% dm = addBiasColumn(dm, 'right') if nargout ~= 1 error('Must assign output back to a design matrix!'); @@ -15,5 +19,17 @@ return end -dm.X = [ones(size(dm.X, 1), 1), dm.X]; -dm.biasCol = 1; % indicating that the first column is the bias column +if nargin < 2 + flag = 'left'; % bias column defaults to the left of the design matrix +end + +switch flag + case 'left' + dm.X = [ones(size(dm.X, 1), 1), dm.X]; + dm.biasCol = 1; % indicating that the first column is the bias column + case 'right' + n = dm.dspec.edim+1; + dm.X = [dm.X, ones(size(dm.X, 1), 1)]; + dm.biasCol = n; % indicating that the last column is the bias column +end + diff --git a/+buildGLM/addCovariate.m b/+buildGLM/addCovariate.m index 4d8140e..aa33511 100644 --- a/+buildGLM/addCovariate.m +++ b/+buildGLM/addCovariate.m @@ -64,7 +64,7 @@ end if nargin >= 6 - dspec.covar(newIdx).offset = offset; + dspec.covar(newIdx).offset = dspec.expt.binfun(offset); else dspec.covar(newIdx).offset = 0; end diff --git a/+buildGLM/addCovariateSpiketrain.m b/+buildGLM/addCovariateSpiketrain.m index 7d66cd6..85dee95 100644 --- a/+buildGLM/addCovariateSpiketrain.m +++ b/+buildGLM/addCovariateSpiketrain.m @@ -33,8 +33,7 @@ assert(ischar(desc), 'Description must be a string'); offset = 1; % Make sure to be causal. No instantaneous interaction allowed. - binfun = dspec.expt.binfun; -stimHandle = @(trial, expt) basisFactory.deltaStim(binfun(trial.(stimLabel)), binfun(trial.duration)); +stimHandle = @(trial, expt) basisFactory.deltaStim(binfun(trial.(stimLabel)+expt.binSize), binfun(trial.duration)); dspec = buildGLM.addCovariate(dspec, covLabel, desc, stimHandle, basisStruct, offset, varargin{:}); \ No newline at end of file diff --git a/+buildGLM/combineWeights.m b/+buildGLM/combineWeights.m index ef63c75..4f817be 100644 --- a/+buildGLM/combineWeights.m +++ b/+buildGLM/combineWeights.m @@ -11,12 +11,20 @@ dspec = dm.dspec; binSize = dspec.expt.binSize; +wout = struct(); if isfield(dm, 'biasCol') % undo z-score operation - wout.bias = w(dm.biasCol); + if isfield(dm, 'zscore') && numel(dm.zscore.mu) == numel(w) % remove bias from zscore + zmu = dm.zscore.sigma(dm.biasCol); + zsig = dm.zscore.mu(dm.biasCol); + dm.zscore.sigma(dm.biasCol) = []; + dm.zscore.mu(dm.biasCol) = []; + else + zmu = 0; + zsig = 1; + end + wout.bias = w(dm.biasCol)*zsig + zmu; % un-z-transform the bias w(dm.biasCol) = []; - dm.zscore.sigma(dm.biasCol) = []; - dm.zscore.mu(dm.biasCol) = []; end if isfield(dm, 'zscore') % undo z-score operation @@ -35,7 +43,6 @@ end startIdx = [1 (cumsum([dspec.covar(:).edim]) + 1)]; -wout = struct(); for kCov = 1:numel(dspec.covar) covar = dspec.covar(kCov); diff --git a/+buildGLM/compileDenseDesignMatrix.m b/+buildGLM/compileDenseDesignMatrix.m index e7a06ed..59fb701 100644 --- a/+buildGLM/compileDenseDesignMatrix.m +++ b/+buildGLM/compileDenseDesignMatrix.m @@ -8,10 +8,9 @@ totalT = sum(trialT); X = zeros(totalT, dspec.edim); -trialIndices = trialIndices(:)'; - -for kTrial = trialIndices - ndx = sum(trialT(1:kTrial))-(trialT(kTrial)-1):sum(trialT(1:kTrial)); +for k = 1:numel(trialIndices) + kTrial = trialIndices(k); + ndx = sum(trialT(1:k))-(trialT(k)-1):sum(trialT(1:k)); for kCov = 1:numel(dspec.covar) % for each covariate covar = dspec.covar(kCov); diff --git a/+buildGLM/getBinnedSpikeTrain.m b/+buildGLM/getBinnedSpikeTrain.m index 178941f..c7c85f1 100644 --- a/+buildGLM/getBinnedSpikeTrain.m +++ b/+buildGLM/getBinnedSpikeTrain.m @@ -6,11 +6,12 @@ endTrialIndices = [0 cumsum(binfun([expt.trial(trialIdx).duration])) + 1]; nT = endTrialIndices(end) - 1; % how many bins total? -for kTrial = trialIdx(:)' - bst = endTrialIndices(kTrial) + binfun(expt.trial(kTrial).(spLabel)); - sts{kTrial} = bst(:); +for k = 1:numel(trialIdx) + kTrial = trialIdx(k); + bst = endTrialIndices(k) + binfun(expt.trial(kTrial).(spLabel)); + sts{k} = bst(:); end sts = cell2mat(sts); - +sts(sts>nT) = []; y = sparse(sts, 1, 1, nT, 1);