From 3ce67164c0210db482b9d97dd2a39727b7f97763 Mon Sep 17 00:00:00 2001 From: Jonathan Pillow Date: Fri, 22 Jan 2016 01:57:59 -0500 Subject: [PATCH] #DONE fixing GLMbi code (bilinear filter); demo 2 works!! --- demos/demo2_GLM_spatialStim.m | 54 +-- glmtools_fitting/Loss_GLM_logli.m | 10 +- glmtools_fitting/Loss_GLMbi_logli.m | 350 +++++++++---------- glmtools_fitting/MLfit_GLM.m | 2 +- glmtools_fitting/MLfit_GLMbi.m | 33 +- glmtools_fitting/initfit_stimDesignMat.m | 8 +- glmtools_fitting/initfit_stimDesignMat_bi.m | 38 ++ glmtools_fitting/initfit_stimMatrix_GLMbi.m | 36 -- glmtools_fitting/makeFittingStruct_GLM.m | 6 +- glmtools_fitting/makeFittingStruct_GLMbi.m | 52 +++ glmtools_fitting/reinsertFitPrs_GLMbi.m | 14 +- glmtools_fitting/setupfitting_GLM.m | 2 +- glmtools_misc/blkdiagrep.m | 7 - glmtools_misc/interp_spikeFilt.m | 27 -- glmtools_misc/makeFittingStruct_GLM_core.m | 58 --- glmtools_misc/makeFittingStruct_GLMbi.m | 42 --- glmtools_misc/makeInterpMatrix2.m | 24 -- unittests/unit_condIntensityConsistency_bi.m | 38 ++ 18 files changed, 350 insertions(+), 451 deletions(-) create mode 100644 glmtools_fitting/initfit_stimDesignMat_bi.m delete mode 100644 glmtools_fitting/initfit_stimMatrix_GLMbi.m create mode 100644 glmtools_fitting/makeFittingStruct_GLMbi.m delete mode 100644 glmtools_misc/blkdiagrep.m delete mode 100644 glmtools_misc/interp_spikeFilt.m delete mode 100644 glmtools_misc/makeFittingStruct_GLM_core.m delete mode 100644 glmtools_misc/makeFittingStruct_GLMbi.m delete mode 100644 glmtools_misc/makeInterpMatrix2.m create mode 100644 unittests/unit_condIntensityConsistency_bi.m diff --git a/demos/demo2_GLM_spatialStim.m b/demos/demo2_GLM_spatialStim.m index 7827a69..f76ecbc 100644 --- a/demos/demo2_GLM_spatialStim.m +++ b/demos/demo2_GLM_spatialStim.m @@ -23,11 +23,12 @@ ttk = dtStim*(-nkt+1:0)'; % time bins for filter kx = 1./sqrt(2*pi*4).*exp(-(xxk-nkx/2).^2/5); k = kt*kx'; % Make space-time separable filter -k = k./norm(k(:))*2; +k = k./norm(k(:))/1.5; % Insert into glm structure (created with default history filter) ggsim = makeSimStruct_GLM(nkt,dtStim,dtSp); % Create GLM structure with default params -ggsim.k = k./norm(k(:))*3; % Insert into simulation struct +ggsim.k = k; % Insert into simulation struct +ggsim.dc = 2; % === Make Fig: model params ======================= clf; subplot(3,3,[1 4]); % ------------------------------------------ @@ -46,7 +47,7 @@ %% 2. Generate some training data ======================================== -slen = 5000; % Stimulus length (frames); More samples gives better fit +slen = 10000; % Stimulus length (frames); More samples gives better fit Stim = round(rand(slen,nkx))*2-1; % Run model on long, binary stimulus [tsp,sps,Itot,Isp] = simGLM(ggsim,Stim); % run model nsp = length(tsp); @@ -86,7 +87,7 @@ nkbasis = 8; % number of basis vectors for representing k nhbasis = 8; % number of basis vectors for representing h hpeakFinal = .1; % time of peak of last basis vector for h -gg0 = makeFittingStruct_GLM(dtStim,dtSp,nkt,nkbasis,nhbasis,hpeakFinal,sta); +gg0 = makeFittingStruct_GLM(dtStim,dtSp,nkt,nkbasis,sta,nhbasis,hpeakFinal); gg0.sps = sps; % Insert binned spike train into fitting struct gg0.mask = exptmask; % insert mask (optional) gg0.ihw = randn(size(gg0.ihw))*1; % initialize spike-history weights randomly @@ -103,16 +104,16 @@ %% 5. Fit GLM ("bilinear stim filter version") via max likelihood % Initialize params for fitting -------------- -Filter_rank = 1; % Number of column/row vector pairs to use -gg0 = makeFittingStruct_GLMbi(sta,dtSp,Filter_rank); -gg0.tsp = tsp; -gg0.mask = exptmask; -[logli1,rr1,tt] = neglogli_GLM(gg0,Stim); % Compute logli of initial params -fprintf('Initial value of negative log-li (GLMbi): %.3f\n', logli1); +k_rank = 1; % Number of column/row vector pairs to use +gg0b = makeFittingStruct_GLMbi(k_rank,dtStim,dtSp,nkt,nkbasis,sta,nhbasis,hpeakFinal); +gg0b.sps = sps; +gg0b.mask = exptmask; +logli0b = neglogli_GLM(gg0b,Stim); % Compute logli of initial params +fprintf('Initial value of negative log-li (GLMbi): %.3f\n', logli0b); % Do ML estimation of model params -opts = {'display', 'iter', 'maxiter', 500}; -[gg2, negloglival2] = MLfit_GLMbi(gg0,Stim,opts); % do ML (requires optimization toolbox) +opts = {'display', 'iter'}; +[gg2, negloglival2] = MLfit_GLMbi(gg0b,Stim,opts); % do ML (requires optimization toolbox) %% 6. Plot results ==================== @@ -121,35 +122,34 @@ subplot(231); % True filter % --------------- imagesc(ggsim.k); colormap gray; title('True Filter');ylabel('time'); - subplot(232); % sta % ------------------------ imagesc(sta); -title('raw STA'); -ylabel('time'); +title('raw STA'); ylabel('time'); subplot(233); % sta-projection % --------------- -imagesc(gg0.k) -title('low-rank STA'); +imagesc(gg0.k); title('low-rank STA'); subplot(234); % estimated filter % --------------- -imagesc(gg1.k) -title('ML estimate: full filter'); xlabel('space'); ylabel('time'); +imagesc(gg1.k); title('ML estimate: full filter'); xlabel('space'); ylabel('time'); subplot(235); % estimated filter % --------------- -imagesc(gg2.k) -title('ML estimate: bilinear filter'); xlabel('space'); +imagesc(gg2.k); title('ML estimate: bilinear filter'); xlabel('space'); subplot(236); % ---------------------------------- plot(ggsim.iht,exp(ggsim.ih),'k', gg1.iht,exp(gg1.ihbas*gg1.ihw),'b',... gg2.iht, exp(gg2.ihbas*gg2.ihw), 'r'); -title('post-spike kernel'); -axis tight; +title('post-spike kernel'); axis tight; +xlabel('time after spike (s)'); % Errors in STA and ML estimate ktmu = normalizecols([mean(ggsim.k,2),mean(gg1.k,2),mean(gg2.k,2)]); kxmu = normalizecols([mean(ggsim.k)',mean(gg1.k)',mean(gg2.k)']); -Errs_T = [subspace(ktmu(:,1),ktmu(:,2)), subspace(ktmu(:,1),ktmu(:,3))] -Errs_X = [subspace(kxmu(:,1),kxmu(:,2)), subspace(kxmu(:,1),kxmu(:,3))] -errfun = @(x,y)(sum((x(:)-y(:)).^2)); -Errs_Total = [errfun(ggsim.k,gg1.k), errfun(ggsim.k, gg2.k)] +msefun = @(k)(sum((k(:)-ggsim.k(:)).^2)); % error function +fprintf(['K-filter errors (GLM vs. GLMbilinear):\n', ... + 'Temporal error: %.3f %.3f\n', ... + ' Spatial error: %.3f %.3f\n', ... + ' Total error: %.3f %.3f\n'], ... + subspace(ktmu(:,1),ktmu(:,2)), subspace(ktmu(:,1),ktmu(:,3)), ... + subspace(kxmu(:,1),kxmu(:,2)), subspace(kxmu(:,1),kxmu(:,3)), ... + msefun(gg1.k), msefun(gg2.k)); diff --git a/glmtools_fitting/Loss_GLM_logli.m b/glmtools_fitting/Loss_GLM_logli.m index 5341468..2029cc2 100644 --- a/glmtools_fitting/Loss_GLM_logli.m +++ b/glmtools_fitting/Loss_GLM_logli.m @@ -9,17 +9,16 @@ % Inputs: % prs = [kprs - weights for stimulus kernel % dc - dc current injection -% ihprs - weights on post-spike current]; +% ihprs - weights on post-spike current] +% % Outputs: % logli = negative log likelihood of spike train % dL = gradient with respect to prs % H = hessian -% Extract some vals from Xstruct (Opt Prs); -nktot = Xstruct.nkx*Xstruct.nkt; % total # params for k -dt = Xstruct.dtSp; % absolute bin size for spike train (in sec) % Unpack GLM prs; +nktot = Xstruct.nkx*Xstruct.nkt; % total # params for k kprs = prs(1:nktot); dc = prs(nktot+1); ihprs = prs(nktot+2:end); @@ -32,6 +31,7 @@ ihflag = Xstruct.ihflag; % flag rlen = Xstruct.rlen; % number of bins in spike train vector nsp = (sum(bsps)); % number of spikes +dt = Xstruct.dtSp; % absolute bin size for spike train (in sec) % -------- Compute sum of filter reponses ----------------------- if Xstruct.ihflag @@ -65,7 +65,7 @@ if ihflag, dLdh1 = (frac1'*XXsp(bsps,:))'; end - % Combine terms + % Combine Term 1 and Term 2 dLdk = dLdk0*dt - dLdk1; dLdb = dLdb0*dt - dLdb1; if ihflag, dLdh = dLdh0*dt - dLdh1; diff --git a/glmtools_fitting/Loss_GLMbi_logli.m b/glmtools_fitting/Loss_GLMbi_logli.m index 5d4d13d..3cfec2c 100644 --- a/glmtools_fitting/Loss_GLMbi_logli.m +++ b/glmtools_fitting/Loss_GLMbi_logli.m @@ -1,7 +1,7 @@ -function [logli, dL, H] = Loss_GLMbi_logli(prs) -% [neglogli, dL, H] = Loss_GLMbi_logli(prs) +function [logli, dL, H] = Loss_GLMbi_logli(prs,Xstruct) +% [logli, dL, H] = Loss_GLMbi_logli(prs,Xstruct) % -% Compute negative log-likelihood of data undr the GLM model, with bilinear +% Compute negative log-likelXXspood of data undr the GLM model, with bilinear % parametrization of the input current % % Uses arbitrary nonlinearity 'nlfun' instead of exponential @@ -10,216 +10,180 @@ % prs = [ktprs - weights for stimulus time kernel % kxprs - weights for stim space kernel % dc - dc current injection -% ihprs - weights on post-spike current]; +% XXspprs - weights on post-spike current]; +% % Outputs: -% logli = negative log likelihood of spike train +% logli = negative log likelXXspood of spike train % dL = gradient with respect to prs -% H = hessian - +% H = hessian -etol = 1e-100; -% global vars for optimization -global SPNDS SPNDS2 OPTprs RefreshRate -%global OPTprs.MSTM MMntrp SPNDS SPNDS2 OPTprs RefreshRate - -% Extract some vals from OPTprs (Opt Prs); -nprs = length(prs); % # of params in param vector -nkt = OPTprs.nkt; % # of basis vectors for spatial k basis -nkx = OPTprs.nkx; % # of basis vectors for spatial k basis -krank = OPTprs.krank; % rank of stim filter +% Extract some size information from Xstruct +nkt = Xstruct.nkt; % # of basis vectors for spatial k basis +nkx = Xstruct.nkx; % # of basis vectors for spatial k basis +krank = Xstruct.krank; % rank of stim filter nktprs = nkt*krank; % total # params for temporal filters nkxprs = nkx*krank; % total # params for spatial filters -nh = OPTprs.nh; % # basis vectors for history current -nh2 = OPTprs.nh2; % # basis vectors for coupling currents -nCoupled = length(SPNDS2); % # of neurons coupled to this one -nhtot = nh+nh2*nCoupled; % % total # of params for history/coupling kernels -dt = OPTprs.dt; % bin size (in frames) -ndt = round(1/dt); % # of time bins per frame -slen = size(OPTprs.MSTM,1); % length of stimulus (low resolution) -% rlen = slen*ndt; % length of stimulus (high resolution) % Unpack GLM prs; ktprs = prs(1:nktprs); kxprs = prs(nktprs+1:nktprs+nkxprs); dc = prs(nktprs+nkxprs+1); -ihprs = prs(nktprs+nkxprs+2:end); -if length(ihprs)~=nhtot - error('Problem with number of params passed in (Loss_GLMbi_logli.m)'); -end +XXspprs = prs(nktprs+nkxprs+2:end); + +% Extract some other stuff we'll use a lot +XXstm = Xstruct.Xstim; % stimulus design matrix +XXsp = Xstruct.Xsp; % spike history design matrix +bsps = Xstruct.bsps; % binary spike vector +M = Xstruct.Minterp; % matrix for interpolating from stimulus bins to spike train bins +ihflag = Xstruct.ihflag; % flag +rlen = Xstruct.rlen; % number of bins in spike train vector +nsp = (sum(bsps)); % number of spikes +dt = Xstruct.dtSp; % absolute bin size for spike train (in sec) % -- Make block-diag matrix containing nkt params ----- -Mkt = blkdiagrep(reshape(ktprs,nkt,krank),nkx); +Mkt = kron(speye(nkx),reshape(ktprs,nkt,krank)); % block-diagonal matrix for kt params inds = reshape((1:nkx*krank),krank,nkx)'; inds = inds(:); -Mkt = Mkt(:,inds); - -% Initialize likelihood, gradient and Hessian ----------- -logli = 0; -dL = zeros(nprs,1); -H = zeros(nprs,nprs); - -for jch = 1:OPTprs.nchunks - %------- Compute indices relevant to this chunk --------------- - iwin = OPTprs.ichunk(jch,:); % abs fine bin win - iwin1 = [iwin(1)+1,iwin(2)]; - jwin = OPTprs.jchunk(jch,:); % abs window of stimulus, coarse - jj = [1,diff(jwin)]; - ilen = diff(iwin); - ii = [(iwin(1)-jwin(1)*ndt+1),(iwin(1)-jwin(1)*ndt+ilen)]; - % if (iwin(1)/ndt-1) == jwin(1) %--- evenly-dividing bins ---- - % ii = [ndt+1, ndt+ilen]; +Mkt = Mkt(:,inds); % turns out we need it reshaped this way + +% -------- Compute bilinear stim filter reponse ----------------------- +dSSdx = XXstm*Mkt; % stim filtered with temporal filter +ystm = dSSdx*kxprs; + +% -------- Compute sum of filter reponses ----------------------- +if ihflag + Itot = M*ystm + XXsp*XXspprs + dc; % stim-dependent + spikehist-dependent inputs +else + Itot = M*ystm + dc; % stim-dependent input only +end + +% --------- Compute output of nonlinearity ------------------------ +[rr,drr,ddrr] = Xstruct.nlfun(Itot); + +% Check for (approximately) zero-values of rr +etol = 1e-100; % cutoff for small values of conditional intensity +iiz = find(rr <= etol); +rr(iiz) = etol; % Set value to small +drr(iiz) = 0; ddrr(iiz) = 0; % Set derivs here to 0 + +% --------- Compute log-likelXXspood --------------------------------- +Trm1 = sum(rr)*dt; % non-spike term +Trm2 = -sum(log(rr(bsps))); % spike term +logli = Trm1 + Trm2; + +% ============================================= +% --------- Compute Gradient ----------------- +% ============================================= +if (nargout > 1) + + % Compute dSSdt - stim filtered with spatial kernel + Mkx = kron(speye(nkt), reshape(kxprs,nkx,krank)); + inds = reshape((1:nkt*krank),krank,nkt)'; inds = inds(:); + Mkx = Mkx(:,inds); + inds = reshape((1:nkt*nkx),nkt,nkx)'; inds = inds(:); + dSSdt = XXstm(:,inds)*Mkx; - % -------- Compute stim filter reponse ----------------------- - SS = OPTprs.MSTM(jwin(1)+1:jwin(2),:); % Relevant stim - dSSdx = SS*Mkt; % stim filtered with temporal filter - MM = OPTprs.MMintrp(ii(1):ii(2),jj(1):jj(2)); - ystm = dSSdx*kxprs; - ystmhi = MM*ystm + dc; - - % -------------- Compute net h current ----------------------- - if OPTprs.ihflag - Ih = zeros(ilen,nh+nh2*nCoupled); - if nh>0 - Ih(:,1:nh) = spikefilt_mex(SPNDS, OPTprs.ihbas, iwin1); - end - for jcpl = 1:nCoupled - Ih(:,nh+nh2*(jcpl-1)+1:nh+nh2*jcpl) = ... - spikefilt_mex(SPNDS2{jcpl}, OPTprs.ihbas2, iwin1); - end - yh = Ih*ihprs; - Iinj = ystmhi+yh; - else - Iinj = ystmhi; + % Non-spiking terms (Term 1) + dqq = drr'*M; + dLdkx0 = (dqq*dSSdx)'; + dLdkt0 = (dqq*dSSdt)'; + dLdb0 = sum(drr); + if ihflag, dLdh0 = (drr'*XXsp)'; end - % --------- Compute likelihood itself ------------------------ - i_sp = in(SPNDS, iwin+.5)-iwin(1); - [rr,drr,ddrr] = OPTprs.nlfun(Iinj); - spflag = ~isempty(i_sp); - - % Check for zero-values of rr - iiz = find(rr <= etol); - rr(iiz) = etol; % Set value to small - drr(iiz) = 0; ddrr(iiz) = 0; % Set derivs here to 0 - - Trm1 = sum(rr)*dt/RefreshRate; % non-spike term - if spflag, - Trm2 = -sum(log(rr(i_sp))); % spike term - else - Trm2 = 0; - end - logli = logli + (Trm1 + Trm2); - - % --------- Compute Gradient ----------------- - if (nargout > 1) - % Compute dSSdt - stim filtered with spatial kernel - Mkx = blkdiagrep(reshape(kxprs,nkx,krank), nkt); - inds = reshape((1:nkt*krank),krank,nkt)'; inds = inds(:); - Mkx = Mkx(:,inds); - inds = reshape((1:nkt*nkx),nkt,nkx)'; inds = inds(:); - dSSdt = SS(:,inds)*Mkx; - - % Non-spiking terms - dqq = drr'*MM; - dLdkx0 = (dqq*dSSdx)'; - dLdkt0 = (dqq*dSSdt)'; - dLdb0 = sum(drr); - if OPTprs.ihflag - dLdh0 = (drr'*Ih)'; - end - if spflag - MMsp = MM(i_sp,:); - %isprel = i_sp+ii(1)-1; - frac1 = drr(i_sp)./rr(i_sp); - fracsp = frac1'*MMsp; - dLdkx1 = (fracsp*dSSdx)'; - dLdkt1 = (fracsp*dSSdt)'; - dLdb1 = sum(frac1); - if OPTprs.ihflag - dLdh1 = (frac1'*Ih(i_sp,:))'; - end - else - dLdkx1=0; dLdkt1=0; dLdb1=0; dLdh1=0; - end - dLdkx = dLdkx0*dt/RefreshRate - dLdkx1; - dLdkt = dLdkt0*dt/RefreshRate - dLdkt1; - dLdk = [dLdkt; dLdkx]; - dLdb = dLdb0*dt/RefreshRate - dLdb1; - if OPTprs.ihflag - dLdh = dLdh0*dt/RefreshRate - dLdh1; - else - dLdh = []; - end - dL = dL + [dLdk; dLdb; dLdh]; - + % Non-spiking terms (Term 2) + Msp = M(bsps,:); % interpolation matrix just for spike bins + frac1 = drr(bsps)./rr(bsps); + fracsp = frac1'*Msp; + dLdkx1 = (fracsp*dSSdx)'; + dLdkt1 = (fracsp*dSSdt)'; + dLdb1 = sum(frac1); + if ihflag, dLdh1 = (frac1'*XXsp(bsps,:))'; end + + % Combine Term 1 and Term 2 + dLdkx = dLdkx0*dt - dLdkx1; + dLdkt = dLdkt0*dt - dLdkt1; + dLdk = [dLdkt; dLdkx]; + dLdb = dLdb0*dt - dLdb1; + if ihflag, dLdh = dLdh0*dt - dLdh1; + else dLdh = []; + end + dL = [dLdk; dLdb; dLdh]; + +end + +% ============================================= +% --------- Compute Hessian ----------------- +% ============================================= +if nargout > 2 + + % --- Non-spiking terms ----- + % multiply each row of M with drr - % --------- Compute Hessian ----------------- - if nargout > 2 - ddrrdiag = spdiags(ddrr,0,ilen,ilen); - ddqqIntrp = ddrrdiag*MM; - ddqq = MM'*ddqqIntrp; - %Hxx and Htt - Hkt = (dSSdt'*ddqq*dSSdt)*dt/RefreshRate; - Hkx = (dSSdx'*ddqq*dSSdx)*dt/RefreshRate; - %Hxt - Hktx0a= reshape((drr'*MM)*SS,nkt,nkx); - Hktx0a = blkdiagrep(Hktx0a,krank); - Hktx = (dSSdx'*ddqq*dSSdt + Hktx0a')*dt/RefreshRate; + ddrrdiag = spdiags(ddrr,0,rlen,rlen); + ddqqIntrp = ddrrdiag*M; + ddqq = M'*ddqqIntrp; % this is MUCH faster than using bsxfun, due to sparsity! + + % Hx and Ht + Hkt = (dSSdt'*ddqq*dSSdt)*dt; + Hkx = (dSSdx'*ddqq*dSSdx)*dt; + % Hxt + Hktx0a= reshape((drr'*M)*XXstm,nkt,nkx); + Hktx0a = kron(speye(krank),Hktx0a); + Hktx = (dSSdx'*ddqq*dSSdt + Hktx0a')*dt; + % Hb + Hb = sum(ddrr)*dt; + % Hkb + sumqq = sum(ddqqIntrp,1); + Hxb = (sumqq*dSSdx)'; + Htb = (sumqq*dSSdt)'; + Hkb = [Htb;Hxb]*dt; + if ihflag + % Hh + Hh = (XXsp'*ddrrdiag*XXsp)*dt; + % Hhk + Hth = (XXsp'*ddqqIntrp)*dSSdt; + Hxh = (XXsp'*ddqqIntrp)*dSSdx; + Hkh = [Hth'; Hxh']*dt; % Hb - Hb = sum(ddrr)*dt/RefreshRate; - % Hkb - sumqq = sum(ddqqIntrp,1); - Hxb = (sumqq*dSSdx)'; - Htb = (sumqq*dSSdt)'; - Hkb = [Htb;Hxb]*dt/RefreshRate; - if OPTprs.ihflag - % Hh - Hh = (Ih'*ddrrdiag*Ih)*dt/RefreshRate; - % Hhk - Hth = (Ih'*ddqqIntrp)*dSSdt; - Hxh = (Ih'*ddqqIntrp)*dSSdx; - Hkh = [Hth'; Hxh']*dt/RefreshRate; - % Hb - Hhb = (ddrr'*Ih)'*dt/RefreshRate; - else - Hkh = []; Hhb=[]; Hh=[]; - end - if spflag % ------- If spikes in window ------ - frac2 = (rr(i_sp).*ddrr(i_sp) - drr(i_sp).^2)./rr(i_sp).^2; - fr2diag = spdiags(frac2,0,length(i_sp),length(i_sp)); - fr2Interp = fr2diag*MMsp; - fr2ddqq = MMsp'*fr2Interp; - % Spiking terms, Hxx and Htt - Hkt= Hkt - dSSdt'*fr2ddqq*dSSdt; - Hkx= Hkx - dSSdx'*fr2ddqq*dSSdx; - % Spiking terms, Hxt - Hktx1a= reshape(frac1'*MMsp*SS,nkt,nkx)'; - Hktx1a= blkdiagrep(Hktx1a,krank); - Hktx1b= dSSdx'*fr2ddqq*dSSdt; - Hktx1= Hktx1a+Hktx1b; - Hktx= Hktx - Hktx1; - % Const term - Hb = Hb-sum(frac2); - % Const by k - sumfr2 = sum(fr2Interp,1); - Hxb1 = sumfr2*dSSdx; - Htb1 = sumfr2*dSSdt; - Hkb = Hkb - [Htb1'; Hxb1']; - if OPTprs.ihflag - Ihrrsp = fr2diag*Ih(i_sp,:); - Hh= Hh - Ih(i_sp,:)'*Ihrrsp; - % Const by h - Hhb1 = sum(Ihrrsp,1)'; - Hhb = Hhb - Hhb1; - % k by h term - Hth0 = Ih(i_sp,:)'*(fr2Interp*dSSdt); - Hxh0 = Ih(i_sp,:)'*(fr2Interp*dSSdx); - Hkh = Hkh-[Hth0';Hxh0']; - end - end - Hk = [[Hkt; Hktx] [Hktx'; Hkx]]; - H = H + [[Hk Hkb Hkh]; [Hkb' Hb Hhb']; [Hkh' Hhb Hh]]; + Hhb = (ddrr'*XXsp)'*dt; + else + Hkh = []; Hhb=[]; Hh=[]; end + % --- Add in spiking terms ---- + frac2 = (rr(bsps).*ddrr(bsps) - drr(bsps).^2)./rr(bsps).^2; + fr2diag = spdiags(frac2,0,nsp,nsp); + fr2Interp = fr2diag*Msp; + fr2ddqq = Msp'*fr2Interp; + % Spiking terms, Hxx and Htt + Hkt= Hkt - dSSdt'*fr2ddqq*dSSdt; + Hkx= Hkx - dSSdx'*fr2ddqq*dSSdx; + % Spiking terms, Hxt + Hktx1a= reshape(frac1'*Msp*XXstm,nkt,nkx)'; + Hktx1a= kron(speye(krank),Hktx1a); + Hktx1b= dSSdx'*fr2ddqq*dSSdt; + Hktx1= Hktx1a+Hktx1b; + Hktx= Hktx - Hktx1; + % Const term + Hb = Hb-sum(frac2); + % Const by k + sumfr2 = sum(fr2Interp,1); + Hxb1 = sumfr2*dSSdx; + Htb1 = sumfr2*dSSdt; + Hkb = Hkb - [Htb1'; Hxb1']; + if ihflag + XXsprrsp = fr2diag*XXsp(bsps,:); + Hh= Hh - XXsp(bsps,:)'*XXsprrsp; + % Const by h + Hhb1 = sum(XXsprrsp,1)'; + Hhb = Hhb - Hhb1; + % k by h term + Hth0 = XXsp(bsps,:)'*(fr2Interp*dSSdt); + Hxh0 = XXsp(bsps,:)'*(fr2Interp*dSSdx); + Hkh = Hkh-[Hth0';Hxh0']; + end + Hk = [[Hkt; Hktx] [Hktx'; Hkx]]; + H = [[Hk Hkb Hkh]; [Hkb' Hb Hhb']; [Hkh' Hhb Hh]]; end diff --git a/glmtools_fitting/MLfit_GLM.m b/glmtools_fitting/MLfit_GLM.m index a677a99..e720a8a 100644 --- a/glmtools_fitting/MLfit_GLM.m +++ b/glmtools_fitting/MLfit_GLM.m @@ -49,5 +49,5 @@ % % ------ Check analytic gradients and Hessians ------- % HessCheck(lfunc,prs0,opts); % HessCheck_Elts(@Loss_GLM_logli, [1 12],prs0,opts); -% tic; [lival,J,H]=Loss_GLM_logli(prs0); toc; +% tic; [lival,J,H]=lfunc(prs0); toc; diff --git a/glmtools_fitting/MLfit_GLMbi.m b/glmtools_fitting/MLfit_GLMbi.m index c4a6e1b..b60935e 100644 --- a/glmtools_fitting/MLfit_GLMbi.m +++ b/glmtools_fitting/MLfit_GLMbi.m @@ -1,4 +1,4 @@ -function [gg, fval,H] = MLfit_GLMbi(gg,Stim,optimArgs); +function [gg,fval,H,Xstruct] = MLfit_GLMbi(gg,Stim,optimArgs) % [gg,fval,H] = MLfit_GLMbi(gg,Stim,optimArgs); % % Computes the ML estimate for GLM params, using grad and hessians. @@ -12,10 +12,8 @@ % Outputs: % ggnew = new param struct (with ML params); % fval = negative log-likelihood at ML estimate -% -% Created: April 08. (from maxli_GLM) - -MAXSIZE = 1e7; % Maximum amount to be held in memory at once; +% H = Hessian of negative log-likelihood at ML estimate +% Xstruct = structure with design matrices for spike-hist and stim terms % Set optimization parameters if nargin > 2 @@ -24,21 +22,28 @@ opts = optimset('Gradobj','on','Hessian','on','display','iter'); end -% Set initial params ---------------------------------- -prs0 = setupfitting_GLM(gg,Stim,MAXSIZE); +% --- Create design matrix extract initial params from gg ---------------- +[prs0,Xstruct] = setupfitting_GLM(gg,Stim); % minimize negative log likelihood -[prs,fval] = fminunc(@Loss_GLMbi_logli,prs0,opts); -if nargout > 2 % Compute Hessian at maximum +lfunc = @(prs)Loss_GLMbi_logli(prs,Xstruct); % loss function for exponential nonlinearity +[prs,fval] = fminunc(lfunc,prs0,opts); % optimize negative log-likelihood for prs + +% Compute Hessian at maximum, if requested +if nargout > 2 [fval,~,H] = Loss_GLMbi_logli(prs); end % Put returned vals back into param structure ------ -gg = reinsertFitPrs_GLMbi(gg,prs); +gg = reinsertFitPrs_GLMbi(gg,prs,Xstruct); -%---------------------------------------------------- + +% %---------------------------------------------------- +% Optional debugging code +% %---------------------------------------------------- +% % % ------ Check analytic gradients, Hessians ------- -% HessCheck(@Loss_GLMbi_logli,prs0,opts); -% HessCheck_Elts(@Loss_GLMbi_logli, [1 12],prs0,opts); -% tic; [lival,J,H]=Loss_GLMbi_logli(prs0); toc; +% DerivCheck(lfunc,prs0,opts); +% HessCheck_Elts(lfunc, [1 12],prs0,opts); +% tic; [lival,J,H]=lfunc(prs0,Xstruct); toc; diff --git a/glmtools_fitting/initfit_stimDesignMat.m b/glmtools_fitting/initfit_stimDesignMat.m index 38e5b5d..0fc62fb 100644 --- a/glmtools_fitting/initfit_stimDesignMat.m +++ b/glmtools_fitting/initfit_stimDesignMat.m @@ -8,14 +8,12 @@ nkt = size(gg.ktbas,2); % # time params per stim pixel (# t params) ncols = nkx*nkt; % total number of columns in stim design matrix -[slen,swid] = size(Stim); +[slen,swid] = size(Stim); % size of stimulus upsampfactor = (gg.dtStim/gg.dtSp); % number of times by which spikes more finely sampled than stim -rlen = slen*upsampfactor; +rlen = slen*upsampfactor; % length of spike train vector % ---- Check size of filter and width of stimulus ---------- -if (nkx ~= swid) - error('Mismatch between stim width and kernel width'); -end +assert(nkx == swid,'Mismatch between stim width and kernel width'); % ---- Convolve stimulus with spatial and temporal bases ----- Xstruct.Xstim = zeros(slen,ncols); diff --git a/glmtools_fitting/initfit_stimDesignMat_bi.m b/glmtools_fitting/initfit_stimDesignMat_bi.m new file mode 100644 index 0000000..37cd110 --- /dev/null +++ b/glmtools_fitting/initfit_stimDesignMat_bi.m @@ -0,0 +1,38 @@ +function Xstruct = initfit_stimDesignMat_bi(gg,Stim) +% Xstruct = initfit_stimDesignMat_bi(gg,Stim) +% +% Initialize parameters relating to stimulus design matrix for bilinearly +% parametrized GLM + + +% ---- Set up filter and stim processing params ------------------- +nkx = size(gg.kxbas,2); % # params per spatial vector +nkt = size(gg.ktbas,2); % # time params per stim pixel (# t params) +ncols = nkx*nkt; % number of columns in design matrix + +[slen,swid] = size(Stim); % size of stimulus +upsampfactor = (gg.dtStim/gg.dtSp); % number of times by which spikes more finely sampled than stim +rlen = slen*upsampfactor; % length of spike train vector + +% ---- Check size of filter and width of stimulus ---------- +assert(nkx == swid,'Mismatch between stim width and kernel width'); + +% ---- Convolve stimulus with spatial and temporal bases ----- +xfltStim = Stim*gg.kxbas; % stimulus filtered with spatial basis +Xstruct.Xstim = zeros(slen,ncols); % +for i = 1:nkx + for j = 1:nkt + Xstruct.Xstim(:,(i-1)*nkt+j) = sameconv(xfltStim(:,i),gg.ktbas(:,j)); + end +end + +% ---- Set up stim processing params ---------- +Xstruct.nkx = nkx; % # stimulus spatial stimulus pixels +Xstruct.nkt = nkt; % # time bins in stim filter +Xstruct.krank = gg.krank; % stim filter rank +Xstruct.slen = slen; % Total stimulus length (coarse bins) +Xstruct.rlen = rlen; % Total spike-train bins (fine bins) +Xstruct.upsampfactor = upsampfactor; % rlen / slen +Xstruct.Minterp = kron(speye(slen),ones(upsampfactor,1)); % rlen x slen matrix for upsampling and downsampling +Xstruct.dtStim = gg.dtStim; % time bin size for stimulus +Xstruct.dtSp = gg.dtSp; % time bins size for spike train diff --git a/glmtools_fitting/initfit_stimMatrix_GLMbi.m b/glmtools_fitting/initfit_stimMatrix_GLMbi.m deleted file mode 100644 index 7cc5da8..0000000 --- a/glmtools_fitting/initfit_stimMatrix_GLMbi.m +++ /dev/null @@ -1,36 +0,0 @@ -function initfit_stimMatrix_GLMbi(gg, Stim) -% initfit_stimMatrix_GLMbi(gg, Stim); -% -% Initialize parameters relating to stimulus design matrix - -global OPTprs - -% ---- Set up filter and stim processing params ------------------- -nkx = size(gg.kxbas,2); % # params per spatial vector -nkt = size(gg.ktbas,2); % # time params per stim pixel (# t params) -krank = gg.krank; % rank of filter -ncols = nkx*nkt; - -[slen,swid] = size(Stim); -rlen = round(slen/gg.dt); % total length (fine bins) - -% ---- Filter stimulus with spatial and temporal bases ----- -xfltStim = Stim*gg.kxbas; -OPTprs.MSTM = zeros(slen,ncols); -for i = 1:nkx - for j = 1:nkt - OPTprs.MSTM(:,(i-1)*nkt+j) = sameconv(xfltStim(:,i),gg.ktbas(:,j)); - end -end - -% ---- Set up stim processing params ---------- -OPTprs.nkx = nkx; -OPTprs.nkt = nkt; -OPTprs.krank = krank; -OPTprs.ktbas = gg.ktbas; -OPTprs.kxbas = gg.kxbas; -OPTprs.slen = slen; % Total stimulus length (course bins) -OPTprs.rlen = rlen; -OPTprs.dt = gg.dt; - - diff --git a/glmtools_fitting/makeFittingStruct_GLM.m b/glmtools_fitting/makeFittingStruct_GLM.m index 86a9f54..d3ad884 100644 --- a/glmtools_fitting/makeFittingStruct_GLM.m +++ b/glmtools_fitting/makeFittingStruct_GLM.m @@ -1,8 +1,8 @@ function gg = makeFittingStruct_GLM(dtStim,dtSp,klength,nkbasis,k0,nhbasis,lasthpeak) -% gg = makeFittingStruct_GLM(dtStim,dtSp,klength,nkbasis,nhbasis,lasthpeak,k0) +% gg = makeFittingStruct_GLM(dtStim,dtSp,klength,nkbasis,k0,nhbasis,lasthpeak) % -% Initialize parameter structure for fitting of GLM model, -% normal parametrization of stim kernel +% Initialize parameter structure for fitting GLM, +% with normal parametrization of stim kernel % % Inputs: % dtStim = bin size of stimulus (in s) diff --git a/glmtools_fitting/makeFittingStruct_GLMbi.m b/glmtools_fitting/makeFittingStruct_GLMbi.m new file mode 100644 index 0000000..6bcb5dd --- /dev/null +++ b/glmtools_fitting/makeFittingStruct_GLMbi.m @@ -0,0 +1,52 @@ +function gg = makeFittingStruct_GLMbi(krank,varargin) +% gg = makeFittingStruct_GLM(krank,dtStim,dtSp,klength,nkbasis,k0,nhbasis,lasthpeak) +% +% Initialize parameter struct for fitting generalized linear model (GLM), +% with bilinear (i.e., low-rank) parametrization stimulus filter +% +% Inputs: +% krank = rank of stim filter +% dtStim = bin size of stimulus (in s) +% dtSp = bin size for spike train (in s) +% klength = temporal length of stimulus filter, in # of bins (optional) +% nkbasis = # temporal basis vectors for stim filter (optional) +% nhbasis = # temporal basis vectors for spike hist filter h (optional) +% lasthpeak = time of peak of last h basis vector, in s (optional) +% k0 = initial estimate of stimulus filter (optional) +% +% Outputs: +% gg = GLM-fitting param structure + + +% ==================================================================== +% Set up structure +gg = makeFittingStruct_GLM(varargin{:}); + +% Set additional fields needed by bilinear GLM +gg.kx = []; +gg.kxbas=[]; +gg.kbasprs = []; +gg.krank = krank; +gg.ktype = 'bilinear'; + +% if initial filter passed in, use svd to set up initial K params (bilinearly parametrized) +if (length(varargin)>4) && (~isempty(varargin{3})) + k0 = varargin{5}; + [u,s,v] = svd(k0); + mux = mean(k0)'; + nkt = size(gg.ktbas,2); + nkx = size(k0,2); + kt = zeros(nkt,krank); + kx = zeros(nkx,krank); + for j = 1:krank; + if v(:,j)'*mux < 0 % Flip sign if shape looks opposite to k0 + u(:,j) = -u(:,j); v(:,j) = -v(:,j); + end + kt(:,j) = (gg.ktbas'*gg.ktbas)\(gg.ktbas'*(u(:,j)*s(j,j))); + kx(:,j) = v(:,j); + end + gg.kxbas = speye(nkx); + gg.kt = kt; + gg.kx = kx; + gg.k = (gg.ktbas*gg.kt)*(gg.kxbas*gg.kx)'; +end diff --git a/glmtools_fitting/reinsertFitPrs_GLMbi.m b/glmtools_fitting/reinsertFitPrs_GLMbi.m index 0a1ecd7..bd7a774 100644 --- a/glmtools_fitting/reinsertFitPrs_GLMbi.m +++ b/glmtools_fitting/reinsertFitPrs_GLMbi.m @@ -1,17 +1,15 @@ -function gg = reinsertFitPrs_GLMbi(gg,prs); +function gg = reinsertFitPrs_GLMbi(gg,prs,Xstruct) % gg = reinsertFitPrs_GLMbi(gg,prs); % % After fitting, reinsert params into param structure -global OPTprs - % Put returned vals back into param structure ------ -krank = OPTprs.krank; -nktprs = OPTprs.nkt*krank; -nkxprs = OPTprs.nkx*krank; +krank = Xstruct.krank; +nktprs = Xstruct.nkt*krank; +nkxprs = Xstruct.nkx*krank; nktot = nktprs+nkxprs; -nh = OPTprs.nh; -nh2 = OPTprs.nh2; +nh = Xstruct.nh; +nh2 = Xstruct.nh2; % Insert params into struct gg.kt = reshape(prs(1:nktprs),[],krank); diff --git a/glmtools_fitting/setupfitting_GLM.m b/glmtools_fitting/setupfitting_GLM.m index 077323b..e4a1b85 100644 --- a/glmtools_fitting/setupfitting_GLM.m +++ b/glmtools_fitting/setupfitting_GLM.m @@ -23,7 +23,7 @@ prs0 = [gg.kt(:); gg.dc; gg.ihw(:); gg.ihw2(:)]; elseif strcmp(gg.ktype, 'bilinear') % bilinearly-parametrized stim filter - Xstruct = initfit_stimMatrix_GLMbi(gg,Stim); % create design matrix structure + Xstruct = initfit_stimDesignMat_bi(gg,Stim); % create design matrix structure % extract parameter vector prs0 = [gg.kt(:); gg.kx(:); gg.dc; gg.ihw(:); gg.ihw2(:)]; diff --git a/glmtools_misc/blkdiagrep.m b/glmtools_misc/blkdiagrep.m deleted file mode 100644 index 5a643f2..0000000 --- a/glmtools_misc/blkdiagrep.m +++ /dev/null @@ -1,7 +0,0 @@ -function B = blkdiagrep(A, mtimes) -% B = blkdiagrep(A, mtimes) -% -% Forms a block-diagonal matrix with m blocks formed from copies of A -% (uses kron) - -B = kron(speye(mtimes),A); \ No newline at end of file diff --git a/glmtools_misc/interp_spikeFilt.m b/glmtools_misc/interp_spikeFilt.m deleted file mode 100644 index 6ef41ea..0000000 --- a/glmtools_misc/interp_spikeFilt.m +++ /dev/null @@ -1,27 +0,0 @@ -function ih1 = interp_spikeFilt(ih0,iht,dt) -% ih1 = interp_spikeFilt(ih0,iht,dt) -% -% If necessary, interpolate post-spike filter to correct time binning and -% cut off any bins <= 0 (i.e., the spike time itself). -% -% Inputs: ih0 = post-spike waveform (or matrix of post-spike waveforms) -% iht = time binning of ih0 -% dt = desired time binning -% -% Outputs ih1 = interpolated, possibly-truncated version of ih0 - -if ~isempty(ih0) - htest = round(10*diff(iht)./dt); % Check if time bins for h same as dt - if all(htest == 10) - ih1 = ih0(iht>0,:); % cut off times <= 0 - else % Interpololate for binning at dt - iht1 = [dt:dt:max(iht)]'; % time points for sampling - ih1 = interp1(iht, ih0, repmat(iht1,1,size(ih0,2)), 'linear', 0); - end - -else % ih0 is empty! - ih1 = zeros(1,size(ih0,2)); - if ~isempty(iht) - fprintf('Warning: no post-spike kernel supplied but iht nonempty!!\n'); - end -end diff --git a/glmtools_misc/makeFittingStruct_GLM_core.m b/glmtools_misc/makeFittingStruct_GLM_core.m deleted file mode 100644 index ca18e23..0000000 --- a/glmtools_misc/makeFittingStruct_GLM_core.m +++ /dev/null @@ -1,58 +0,0 @@ -function gg = makeFittingStruct_GLM_core(sta,dtStim,dtSp) -% gg = makeFittingStruct_GLM_core(sta,dtStim,dtSp) -% -% Core common to both makeFittingStruct_GLM and makeFittingStruct_GLMbi -% -% Inputs: sta = initial guess at kernel (or use all zeros if unknown) -% DTsim = bin size for simulations / likelihood calculation -% glmstruct (optional) = glm sim structure to extract from - - -% Set up structure ===================================== -gg.k = []; % Actual stim filter k -gg.ih = []; % Actual post-spike filter ih -gg.dc = 0; -gg.ihw = []; % ih weights -gg.iht = []; % ih time points -gg.ihbas = []; -gg.ihbas2 = []; -gg.ihbasprs = []; -gg.kt = []; -gg.ktbas = []; -gg.ktbasprs = []; -gg.nlfun = @expfun; % default nonlinearity: exponential -gg.tsp = []; -gg.mask = []; -gg.dtStim = dtStim; -gg.dtSp = dtSp; -gg.ihw2 = []; -gg.ihbas2 = []; -gg.ihbasprs2 = []; -gg.tsp2 = []; -gg.couplednums = []; - -nkt = size(sta,1); - -% % ----- Set up temporal basis for stimulus kernel ----------- -ktbasprs.neye = min(5,nkt); % Number of "identity" basis vectors near time of spike; -ktbasprs.ncos = min(5,nkt); % Number of raised-cosine vectors to use -ktbasprs.kpeaks = [0 ((nkt-ktbasprs.neye)/2)]; % Position of 1st and last bump -ktbasprs.b = 1; % Offset for nonlinear scaling (larger -> more linear) -ktbas = makeBasis_StimKernel(ktbasprs,nkt); -gg.ktbas = ktbas; -gg.ktbasprs = ktbasprs; - -% ====================================================================== -% Make default basis for post-spike filter - -ihbasprs.ncols = 5; % Number of basis vectors for post-spike kernel -ihbasprs.hpeaks = [0 max(dtSp*10,4)]; % Peak location for first and last vectors -ihbasprs.b = 1; % How nonlinear to make spacings -ihbasprs.absref = []; % absolute refractory period (optional) -[iht,ihbas] = makeBasis_PostSpike(ihbasprs,dtSp); -gg.iht = iht; -gg.ihbas = ihbas; -gg.ihbasprs = ihbasprs; -gg.ihw = zeros(size(ihbas,2),1); -gg.ih = gg.ihbas*gg.ihw; % Set ih to be the current value of ihbas*ihw - diff --git a/glmtools_misc/makeFittingStruct_GLMbi.m b/glmtools_misc/makeFittingStruct_GLMbi.m deleted file mode 100644 index b423f45..0000000 --- a/glmtools_misc/makeFittingStruct_GLMbi.m +++ /dev/null @@ -1,42 +0,0 @@ -function gg = makeFittingStruct_GLMbi(sta,DTsim,krank,varargin) -% gg = makeFittingStruct_GLMbi(sta,DTsim,krank,glmstruct); -% -% Initialize parameter structure for fitting of GLM model -% with bilinear parametrization of stim kernel -% -% Inputs: sta = initial guess at kernel (or use all zeros if unknown) -% DTsim = bin size for simulations / likelihood calculation -% krank = rank of stim kernel "k" -% glmstruct (optional) = glm sim structure to extract from - - -% ==================================================================== -% Set up structure -gg = makeFittingStruct_GLM(sta,DTsim,varargin{:}); -gg.kx = []; -gg.kxbas=[]; -gg.kbasprs = []; -gg.krank = krank; - -% % ================================================================== -% use svd to set up initial K params (bilinearly parametrized) -[u,s,v] = svd(sta); -mux = mean(sta)'; -nkt = size(gg.ktbas,2); -nkx = size(sta,2); - -kt = zeros(nkt,krank); -kx = zeros(nkx,krank); -for j = 1:krank; - if v(:,j)'*mux < 0 % Flip sign if match is better - u(:,j) = -u(:,j); v(:,j) = -v(:,j); - end - kt(:,j) = (gg.ktbas'*gg.ktbas)\(gg.ktbas'*(u(:,j)*s(j,j))); - kx(:,j) = v(:,j); -end -gg.kxbas = speye(nkx); -gg.kt = kt; -gg.kx = kx; -gg.k = (gg.ktbas*gg.kt)*(gg.kxbas*gg.kx)'; - -gg.ktype = 'bilinear'; \ No newline at end of file diff --git a/glmtools_misc/makeInterpMatrix2.m b/glmtools_misc/makeInterpMatrix2.m deleted file mode 100644 index dfb7ec2..0000000 --- a/glmtools_misc/makeInterpMatrix2.m +++ /dev/null @@ -1,24 +0,0 @@ -function M = makeInterpMatrix2(slen, nbn) -% M = makeInterpMatrix2(slen, dt); -% -% Make (sparse matrix) for interpolation. -% -% (second arg can be # bins or bin size). -% -% Matrix performs the convolution equivalent to "fastinterp2.c"; - -if nbn<1 - nbn = round(1./nbn); -end - -c1 = [1/nbn:1/nbn:1]'; -cc = [c1; flipud(c1)-1/nbn]; -M = spalloc(nbn*slen,slen, 2*nbn*slen); -phi = floor(nbn/2); -M(1:nbn+phi)= cc(phi+1:end); -for j = 2:slen-1 - M((j-1)*nbn-phi+1:(j+1)*nbn-phi,j) = cc; -end -M(nbn*slen-nbn-phi+1:nbn*slen,slen) = cc(1:nbn+phi); - - diff --git a/unittests/unit_condIntensityConsistency_bi.m b/unittests/unit_condIntensityConsistency_bi.m new file mode 100644 index 0000000..f631a2f --- /dev/null +++ b/unittests/unit_condIntensityConsistency_bi.m @@ -0,0 +1,38 @@ + +% --------------- +% Make param object with "true" params (if desired). +% --------------- +Filter_rank = 1; +ggTrue = makeFittingStruct_GLMbi(ggsim.k,dtSp,Filter_rank,ggsim); +% Set true kernel in this struct (i.e. not represented by default basis). +[u,s,v] = svd(ggsim.k); +ggTrue.k = ggsim.k; +ggTrue.ktbas = eye(nkt); +ggTrue.kt = u(:,1); +ggTrue.kx = v(:,1)*s(1,1); +% Insert spike times +ggTrue.tsp = tsp; + +% --------------- +% Check that conditional intensity calc is correct +% (if desired, compare to vmem returned by simGLM above) +[logliTrue, rrTrue,tt] = neglogli_GLM(ggTrue,Stim); +subplot(211); plot(tt,vmem,tt,log(rrTrue)); +title('total filter output (computed 2 ways)'); +subplot(212); plot(tt,log(rrTrue)-vmem); +title('difference'); + +% --------------- +[logliTrue, rrT,tt,Iinj,Istm,Ih] = neglogli_GLM(ggTrue,Stim); + +%% 5. Unit tests + +% On total log-conditional intensity +assert(max(abs(vmem-Iinj))<1e-8,'Unit failed: log-conditional intensity consistency'); + +% On just the spike-history component of the conditional intensity +assert(max(abs(ispk-Ih))<1e-8,'Unit failed: spike-history filter output consistency'); + +fprintf('Unit passed: condIntensityConsistency\n'); + +