Cross validation error

[RobertoDF]

Hi, I run in the following error during cross validation:

for ID in ecephys_session_table.iloc[:3]["behavior_session_id"]:
    ps.process_session(ID, version=1)
    ps.plot_fit(ID)
    ps.plot_session_summary(IDS)

/alzheimer/Roberto/Allen_Institute/alexpiet/psy_fits_v1/session_fits/1044408432
Starting Fit now
Pulling Data
Loading SDK object
removing passive session stimuli
Adding stimulus annotations
Annotating lick bouts
Formating Data
Loading options for version 1
This session had 3040 licks and 127 rewards
Initial Fit
{'omissions1': 1, 'bias': 1, 'timing1D': 1, 'omissions': 1, 'task0': 1}
WARNING: sigDay has no effect, dayLength not supplied in dat

/alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/allensdk/brain_observatory/behavior/data_objects/metadata/subject_metadata/full_genotype.py:57: UserWarning: Unable to parse cre_line from full_genotype
  warnings.warn('Unable to parse cre_line from full_genotype')

Cross Validation Analysis
running fold 9

---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[10], line 2
      1 for ID in ecephys_session_table.iloc[:3]["behavior_session_id"]:
----> 2     ps.process_session(ID, version=1)
      3     ps.plot_fit(ID)
      4     ps.plot_session_summary(IDS)

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/licking_behavior_NP/psy_tools.py:90, in process_session(bsid, complete, version, format_options, refit)
     86 print("Cross Validation Analysis")
     87 cross_psydata = psy.trim(psydata, 
     88     END=int(np.floor(len(psydata['y'])/format_options['num_cv_folds'])\
     89     *format_options['num_cv_folds'])) 
---> 90 cross_results = compute_cross_validation(cross_psydata, hyp, weights,
     91     folds=format_options['num_cv_folds'])
     92 cv_pred = compute_cross_validation_ypred(cross_psydata, cross_results,ypred)
     94 if complete:

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/licking_behavior_NP/psy_tools.py:657, in compute_cross_validation(psydata, hyp, weights, folds)
    655 for k in range(folds):
    656     print("\rrunning fold " +str(k),end="") 
--> 657     _,_,wMode_K,_ = psy.hyperOpt(trainDs[k], hyp, weights, ['sigma'],hess_calc=None)
    658     logli, gw = xval_loglike(testDs[k], wMode_K, trainDs[k]['missing_trials'], 
    659         weights)
    660     res = {'logli' : np.sum(logli), 'gw' : gw, 'test_inds' : testDs[k]['test_inds']}

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/psytrack/hyperOpt.py:177, in hyperOpt(dat, hyper, weights, optList, method, showOpt, jump, hess_calc)
    174     else:
    175         optVals += np.log2(current_hyper[val]).tolist()
--> 177 result = minimize(
    178     hyperOpt_lossfun,
    179     optVals,
    180     args=opt_keywords,
    181     method='BFGS',
    182     options=opts,
    183     callback=callback,
    184 )
    186 diff = np.linalg.norm((optVals - result.x) / optVals)
    187 if showOpt:

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/scipy/optimize/_minimize.py:691, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
    689     res = _minimize_cg(fun, x0, args, jac, callback, **options)
    690 elif meth == 'bfgs':
--> 691     res = _minimize_bfgs(fun, x0, args, jac, callback, **options)
    692 elif meth == 'newton-cg':
    693     res = _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,
    694                              **options)

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/scipy/optimize/_optimize.py:1388, in _minimize_bfgs(fun, x0, args, jac, callback, gtol, norm, eps, maxiter, disp, return_all, finite_diff_rel_step, xrtol, **unknown_options)
   1385 pk = -np.dot(Hk, gfk)
   1386 try:
   1387     alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \
-> 1388              _line_search_wolfe12(f, myfprime, xk, pk, gfk,
   1389                                   old_fval, old_old_fval, amin=1e-100, amax=1e100)
   1390 except _LineSearchError:
   1391     # Line search failed to find a better solution.
   1392     warnflag = 2

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/scipy/optimize/_optimize.py:1160, in _line_search_wolfe12(f, fprime, xk, pk, gfk, old_fval, old_old_fval, **kwargs)
   1146 """
   1147 Same as line_search_wolfe1, but fall back to line_search_wolfe2 if
   1148 suitable step length is not found, and raise an exception if a
   (...)
   1155 
   1156 """
   1158 extra_condition = kwargs.pop('extra_condition', None)
-> 1160 ret = line_search_wolfe1(f, fprime, xk, pk, gfk,
   1161                          old_fval, old_old_fval,
   1162                          **kwargs)
   1164 if ret[0] is not None and extra_condition is not None:
   1165     xp1 = xk + ret[0] * pk

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/scipy/optimize/_linesearch.py:84, in line_search_wolfe1(f, fprime, xk, pk, gfk, old_fval, old_old_fval, args, c1, c2, amax, amin, xtol)
     80     return np.dot(gval[0], pk)
     82 derphi0 = np.dot(gfk, pk)
---> 84 stp, fval, old_fval = scalar_search_wolfe1(
     85         phi, derphi, old_fval, old_old_fval, derphi0,
     86         c1=c1, c2=c2, amax=amax, amin=amin, xtol=xtol)
     88 return stp, fc[0], gc[0], fval, old_fval, gval[0]

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/scipy/optimize/_linesearch.py:160, in scalar_search_wolfe1(phi, derphi, phi0, old_phi0, derphi0, c1, c2, amax, amin, xtol)
    158 if task[:2] == b'FG':
    159     alpha1 = stp
--> 160     phi1 = phi(stp)
    161     derphi1 = derphi(stp)
    162 else:

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/scipy/optimize/_linesearch.py:75, in line_search_wolfe1.<locals>.phi(s)
     73 def phi(s):
     74     fc[0] += 1
---> 75     return f(xk + s*pk, *args)

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/scipy/optimize/_differentiable_functions.py:267, in ScalarFunction.fun(self, x)
    265 if not np.array_equal(x, self.x):
    266     self._update_x_impl(x)
--> 267 self._update_fun()
    268 return self.f

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/scipy/optimize/_differentiable_functions.py:251, in ScalarFunction._update_fun(self)
    249 def _update_fun(self):
    250     if not self.f_updated:
--> 251         self._update_fun_impl()
    252         self.f_updated = True

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/scipy/optimize/_differentiable_functions.py:155, in ScalarFunction.__init__.<locals>.update_fun()
    154 def update_fun():
--> 155     self.f = fun_wrapped(self.x)

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/scipy/optimize/_differentiable_functions.py:137, in ScalarFunction.__init__.<locals>.fun_wrapped(x)
    133 self.nfev += 1
    134 # Send a copy because the user may overwrite it.
    135 # Overwriting results in undefined behaviour because
    136 # fun(self.x) will change self.x, with the two no longer linked.
--> 137 fx = fun(np.copy(x), *args)
    138 # Make sure the function returns a true scalar
    139 if not np.isscalar(fx):

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/psytrack/hyperOpt.py:291, in hyperOpt_lossfun(optVals, keywords)
    289 # Calculate posterior term, then approximate evidence for new sigma
    290 center = DL_2 + lT['ddlogli']['H']
--> 291 logterm_post = (1 / 2) * sparse_logdet(center)
    293 evd = pT['logprior'] + lT['logli'] - logterm_post
    295 return -evd

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/psytrack/helper/helperFunctions.py:49, in sparse_logdet(A)
     45 if not isspmatrix_csc(A):  # needed for splu() decomposition to work
     46     raise Exception(
     47         'sparse_logdet: matrix passed is not in sparse csc form')
---> 49 aux = splu(A)
     50 return np.sum(
     51     np.log(np.abs(aux.L.diagonal())) + np.log(np.abs(aux.U.diagonal())))

File /alzheimer/Roberto/Software/mambaforge/envs/icking_behavior/lib/python3.9/site-packages/scipy/sparse/linalg/_dsolve/linsolve.py:413, in splu(A, permc_spec, diag_pivot_thresh, relax, panel_size, options)
    410 if (_options["ColPerm"] == "NATURAL"):
    411     _options["SymmetricMode"] = True
--> 413 return _superlu.gstrf(N, A.nnz, A.data, A.indices, A.indptr,
    414                       csc_construct_func=csc_construct_func,
    415                       ilu=False, options=_options)

RuntimeError: Factor is exactly singular

[alexpiet]

Yes, this happens occasionally. The reason it happens is if during a random split of the data one of the cross validation chunks happens to not have any timepoints that contain a specific model regressor (like no omissions, or no changes, or no licks), then the model is singular and we can't fit the model. In my hands, this happens pretty infrequently, so I simply ignored the sessions where it happened. It also tends to happen the most in sessions where the mice lick very little, so I would probably exclude those sessions anyways. If you set the cross validation to be larger (like 10x validation), it should happen even less frequently.