Optimization of time constants when solving the AR2 Model

[ehennestad]

I am trying to solve the 2nd order AR model and I want to optimize the time constants. The only method that seem to provide this option is the thresholded. However, when setting 'optimize_pars' and 'optimize_smin' to true (to get into the right case in deconvolveCa) the code exits with an error.

The first error happens on line 136 in thresholded_oasisAR2.m. The function 'update_smin' takes 8 inputs, but only 7 is given, (s_max is missing).

I changed line 136-137 to the following: [smin, solution, spks, active_set] = update_smin(y-b, g, smin,... solution, spks, active_set, sqrt(thresh), max(spks)); And tried again.

This time, I get an error on line 143 in the same function. The variable 'lam' is passed to function update_g, but lam is not previously defined. There is a comment on line 108 to update lam, but there is no reference to lam there.

[epnev]

@ehennestad I'll look more into it but I believe that thresholded_OASIS is basically run with lam=0. Instead of the regularization you put a hard threshold on the minimum spike amplitude s_min. Can you replace lam with 0 on line 143. The AR1 version of thresholded_OASIS is better documented.

[ehennestad]

Yes, now I get the following:

Not enough input arguments.

Error in thresholded_oasisAR2>update_g (line 237)
[c,s,active_set] = oasisAR2(y, g, 0, smin);

Error in thresholded_oasisAR2 (line 144)
                [solution, active_set, g, spks] = update_g(y-b, active_set, lam);

Error in deconvolveCa (line 145)
                [c, s, options.b, options.pars, options.smin] = thresholded_oasisAR2(y,...

The problem is that 'smin' which is passed to oasisAR2 is not present, because it in turn is not passed to update_g in line 144 of thresholded_oasisAR2.

First I set smin to zero before line 237, which is the default value given to it within oasisAR2, and then I get a final solution, but it consist of complex numbers.

Second, I change line 144 to this: [solution, active_set, g, spks] = update_g(y-b, active_set, lam, smin);

The initial value of smin is 0. This was a necessary condition to be able to call thresholded_oasisAR2 from deconvolveCa in the first place. Now, on line 137-138 in thresholded_oasisAR2 (right before the call to update_g), 'update_smin' is called. However, update_smin returns smin=0, so the solution of the deconvolution is still complex.

Next I bypass the condition of deconvolveCa, and pass a non-zero value of smin to thresholded_oasisAR2, this time 1.5*getSn(y) (which makes smin=0.056). The value of s_min which is returned from update_smin is still identical to what was passed to update_smin.

Looking at the function update_smin found here. It seems like smin is only updated when the error is below a threshold value. So even if the error gets smaller, if it doesn't fall below this threshold, smin is not updated, even if it continues to be zero. Is this intentional? I checked out update_smin from thresholded_oasisAR1, and this function is identical.

The second thing to note is that the time constants do change from the initial estimate, but they become much lower. tau_rise equals 1ms and tau_decay equals 30ms. This happens regardless of whether I start with smin = 0 or smin = 1.5*getSn(y).

[epnev]

@ehennestad The bug should be fixed with f2eb9cb6960d4f963451873923c76f62448a17e8 For some reason the update_g function is called in a different way depending on whether you are in AR1 or AR2. I haven't written this function so not sure why this happening. Perhaps @zhoupc might be able to help if needed.

[ehennestad]

Thank you for fixing this! I tried to run it, but it ended up taking a very long time. Running the thresholded method (AR2) without updating time constants usually take less than a second for one roi (~50000 samples). Setting optimize_pars to true, I was waiting around 15 minutes. When I paused the script it was on the 2nd out of the 10th iteration. Is this expected?

[epnev]

@ehennestad Unfortunately that might be the case for this particular combination (thresholded, AR2, 50k samples, optimize_pars=true, matlab). The python implementation relies on some cython code and is considerably faster.

One thing you can do, is run the deconvolution on a shorter (contiguous) segment (e.g., 5k samples) in this fully general setup and once it converges, grab all the parameters and run it on the full trace with these parameters and set optimize_pars=false. This should give you some decent speed-up.

[ehennestad]

Thanks, thats a good idea, I will definitely try that. Do you know how long it would take in python? I could also use python for the deconvolution.