From 34ec885a4993327b4c5bc84de4b96cd200a276eb Mon Sep 17 00:00:00 2001 From: Jacob Yates Date: Tue, 3 Mar 2015 00:08:01 -0600 Subject: [PATCH 01/11] add catch for rounding errors that push events outside of the trial --- +basisFactory/deltaStim.m | 1 + 1 file changed, 1 insertion(+) 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 From 1c6b667e31dab9ac8996a0ab769d4ce3ddaa146e Mon Sep 17 00:00:00 2001 From: Jacob Yates Date: Tue, 3 Mar 2015 00:12:54 -0600 Subject: [PATCH 02/11] need to offset the time of the event before binning... adding an offset after abolishes reractory periods for cells. We should follow up on this. --- +buildGLM/addCovariateSpiketrain.m | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) 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 From 4c8be40aafa9fead1a83885d226220809d5e29f1 Mon Sep 17 00:00:00 2001 From: Jacob Yates Date: Tue, 3 Mar 2015 00:14:32 -0600 Subject: [PATCH 03/11] correct the bias parameter for zscoring --- +buildGLM/combineWeights.m | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/+buildGLM/combineWeights.m b/+buildGLM/combineWeights.m index ef63c75..b6b78d7 100644 --- a/+buildGLM/combineWeights.m +++ b/+buildGLM/combineWeights.m @@ -13,10 +13,17 @@ binSize = dspec.expt.binSize; if isfield(dm, 'biasCol') % undo z-score operation - wout.bias = w(dm.biasCol); + if isfield(dm, 'zscore') % 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 = 1; + 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 From b1c0830a10f135910c577787ae0c7ccd78dfc632 Mon Sep 17 00:00:00 2001 From: Jacob Yates Date: Wed, 4 Mar 2015 09:29:24 -0600 Subject: [PATCH 04/11] basis factory functions take optional arguments to normalize or orthogonalize the basis --- +basisFactory/makeNonlinearRaisedCos.m | 19 +++++- +basisFactory/makeSmoothTemporalBasis.m | 79 +++++++++++++++---------- 2 files changed, 66 insertions(+), 32 deletions(-) diff --git a/+basisFactory/makeNonlinearRaisedCos.m b/+basisFactory/makeNonlinearRaisedCos.m index 4e76019..008a501 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; 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]; From 0cbf5e92f871b9c4de9f37a199be839d455be4be Mon Sep 17 00:00:00 2001 From: Jacob Yates Date: Wed, 4 Mar 2015 09:31:27 -0600 Subject: [PATCH 05/11] fixed bug that TrialIndices must start from 1 and be sequential --- +buildGLM/compileDenseDesignMatrix.m | 2 +- +buildGLM/getBinnedSpikeTrain.m | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/+buildGLM/compileDenseDesignMatrix.m b/+buildGLM/compileDenseDesignMatrix.m index e7a06ed..06d1e82 100644 --- a/+buildGLM/compileDenseDesignMatrix.m +++ b/+buildGLM/compileDenseDesignMatrix.m @@ -4,7 +4,7 @@ expt = dspec.expt; subIdxs = buildGLM.getGroupIndicesFromDesignSpec(dspec); -trialT = ceil([expt.trial(trialIndices).duration]/expt.binSize); +trialT = ceil([expt.trial(:).duration]/expt.binSize); totalT = sum(trialT); X = zeros(totalT, dspec.edim); diff --git a/+buildGLM/getBinnedSpikeTrain.m b/+buildGLM/getBinnedSpikeTrain.m index 178941f..0b027c2 100644 --- a/+buildGLM/getBinnedSpikeTrain.m +++ b/+buildGLM/getBinnedSpikeTrain.m @@ -3,7 +3,7 @@ sts = cell(numel(trialIdx), 1); binfun = expt.binfun; -endTrialIndices = [0 cumsum(binfun([expt.trial(trialIdx).duration])) + 1]; +endTrialIndices = [0 cumsum(binfun([expt.trial(:).duration])) + 1]; nT = endTrialIndices(end) - 1; % how many bins total? for kTrial = trialIdx(:)' From bf7ef9cbea0dc81647b54e2958a9ffb3a9cc9c0b Mon Sep 17 00:00:00 2001 From: Jacob Yates Date: Mon, 9 Mar 2015 17:56:27 -0500 Subject: [PATCH 06/11] fixed bug with zscored weights --- +buildGLM/combineWeights.m | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/+buildGLM/combineWeights.m b/+buildGLM/combineWeights.m index b6b78d7..4f817be 100644 --- a/+buildGLM/combineWeights.m +++ b/+buildGLM/combineWeights.m @@ -11,15 +11,16 @@ dspec = dm.dspec; binSize = dspec.expt.binSize; +wout = struct(); if isfield(dm, 'biasCol') % undo z-score operation - if isfield(dm, 'zscore') % remove bias from zscore + 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 = 1; + zmu = 0; zsig = 1; end wout.bias = w(dm.biasCol)*zsig + zmu; % un-z-transform the bias @@ -42,7 +43,6 @@ end startIdx = [1 (cumsum([dspec.covar(:).edim]) + 1)]; -wout = struct(); for kCov = 1:numel(dspec.covar) covar = dspec.covar(kCov); From cb65c34c78d0b56f30a5c0b56023327bc75326e1 Mon Sep 17 00:00:00 2001 From: Jacob Yates Date: Tue, 10 Mar 2015 00:22:06 -0500 Subject: [PATCH 07/11] reduce memory load by only using bins from trialIndices --- +buildGLM/compileDenseDesignMatrix.m | 9 ++++----- +buildGLM/getBinnedSpikeTrain.m | 9 +++++---- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/+buildGLM/compileDenseDesignMatrix.m b/+buildGLM/compileDenseDesignMatrix.m index 06d1e82..59fb701 100644 --- a/+buildGLM/compileDenseDesignMatrix.m +++ b/+buildGLM/compileDenseDesignMatrix.m @@ -4,14 +4,13 @@ expt = dspec.expt; subIdxs = buildGLM.getGroupIndicesFromDesignSpec(dspec); -trialT = ceil([expt.trial(:).duration]/expt.binSize); +trialT = ceil([expt.trial(trialIndices).duration]/expt.binSize); 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 0b027c2..49ada50 100644 --- a/+buildGLM/getBinnedSpikeTrain.m +++ b/+buildGLM/getBinnedSpikeTrain.m @@ -3,12 +3,13 @@ sts = cell(numel(trialIdx), 1); binfun = expt.binfun; -endTrialIndices = [0 cumsum(binfun([expt.trial(:).duration])) + 1]; +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); From 62e25f34ba2d2ced0d668ebeeb0961282e8f1ca4 Mon Sep 17 00:00:00 2001 From: Jacob Yates Date: Thu, 7 May 2015 11:15:52 -0500 Subject: [PATCH 08/11] bs.tr should be in bins, not time... divide by binsize --- +basisFactory/makeNonlinearRaisedCos.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/+basisFactory/makeNonlinearRaisedCos.m b/+basisFactory/makeNonlinearRaisedCos.m index 008a501..16982f0 100644 --- a/+basisFactory/makeNonlinearRaisedCos.m +++ b/+basisFactory/makeNonlinearRaisedCos.m @@ -57,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 From 9b6b51de280c1c4c5264a836b37b2df7da5482fb Mon Sep 17 00:00:00 2001 From: Jacob Yates Date: Thu, 7 May 2015 11:17:59 -0500 Subject: [PATCH 09/11] offset should be in bins --- +buildGLM/addCovariate.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From 9233f0fdd8876e5fb2cf528508c23551f0d77a94 Mon Sep 17 00:00:00 2001 From: Jacob Yates Date: Thu, 28 May 2015 20:33:55 -0500 Subject: [PATCH 10/11] can't have spikes that happen after the trial ends --- +buildGLM/getBinnedSpikeTrain.m | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/+buildGLM/getBinnedSpikeTrain.m b/+buildGLM/getBinnedSpikeTrain.m index 49ada50..c7c85f1 100644 --- a/+buildGLM/getBinnedSpikeTrain.m +++ b/+buildGLM/getBinnedSpikeTrain.m @@ -13,5 +13,5 @@ end sts = cell2mat(sts); - +sts(sts>nT) = []; y = sparse(sts, 1, 1, nT, 1); From 38dc94b295aebfb2d2bf8333b4b3439ea29cd4d5 Mon Sep 17 00:00:00 2001 From: Jacob Yates Date: Thu, 28 May 2015 20:44:36 -0500 Subject: [PATCH 11/11] add option for bias column on the right --- +buildGLM/addBiasColumn.m | 24 ++++++++++++++++++++---- 1 file changed, 20 insertions(+), 4 deletions(-) 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 +