---------------------------------------------------------------------------
LinAlgError Traceback (most recent call last)
Input In [29], in <cell line: 1>()
----> 1 sn.fit()
File /opt/conda/anaconda3/envs/pisco/lib/python3.10/site-packages/piscola/sn_class.py:281, in Supernova.fit(self, kernel1, kernel2, gp_mean)
264 def fit(self, kernel1="matern52", kernel2="squaredexp", gp_mean="mean"):
265 """Fits and corrects the multi-colour light curves.
266
267 The corrections include Milky-Way dust extinction and mangling of the SED.
(...)
279 Gaussian process mean function. Either ``mean``, ``max`` or ``min``.
280 """
--> 281 self._fit_lcs(kernel1, kernel2, gp_mean) # to get initial tmax
283 sed_lcs = self.sed.obs_lcs_fit # interpolated light curves
284 sed_times = sed_lcs.phase.values + self.init_tmax
File /opt/conda/anaconda3/envs/pisco/lib/python3.10/site-packages/piscola/sn_class.py:228, in Supernova._fit_lcs(self, kernel1, kernel2, gp_mean)
212 """Fits the multi-colour light-curve data with gaussian process.
213
214 The time of rest-frame B-band peak luminosity is estimated by finding where the derivative is equal to zero.
(...)
225 Gaussian process mean function. Either ``mean``, ``max`` or ``min``.
226 """
227 self._stack_lcs()
--> 228 timeXwave, lc_mean, lc_std, gp_pred, gp = gp_2d_fit(
229 self._stacked_time,
230 self._stacked_wave,
231 self._stacked_flux,
232 self._stacked_flux_err,
233 kernel1,
234 kernel2,
235 gp_mean,
236 )
238 self.init_fits["timeXwave"], self.init_fits["lc_mean"] = timeXwave, lc_mean
239 self.init_fits["lc_std"], self.init_fits["gp_pred"] = lc_std, gp_pred
File /opt/conda/anaconda3/envs/pisco/lib/python3.10/site-packages/piscola/gaussian_process.py:190, in gp_2d_fit(x1_data, x2_data, y_data, yerr_data, kernel1, kernel2, gp_mean)
188 # optimization routine for hyperparameters
189 p0 = gp.get_parameter_vector()
--> 190 results = scipy.optimize.minimize(
191 neg_ln_like, p0, jac=grad_neg_ln_like, method="BFGS"
192 )
193 gp.set_parameter_vector(results.x)
195 # extrapolation edges
File /opt/conda/anaconda3/envs/pisco/lib/python3.10/site-packages/scipy/optimize/_minimize.py:694, in minimize(fun, x0, args, method, jac, hess, hessp, bounds, constraints, tol, callback, options)
692 res = _minimize_cg(fun, x0, args, jac, callback, **options)
693 elif meth == 'bfgs':
--> 694 res = _minimize_bfgs(fun, x0, args, jac, callback, **options)
695 elif meth == 'newton-cg':
696 res = _minimize_newtoncg(fun, x0, args, jac, hess, hessp, callback,
697 **options)
File /opt/conda/anaconda3/envs/pisco/lib/python3.10/site-packages/scipy/optimize/_optimize.py:1309, in _minimize_bfgs(fun, x0, args, jac, callback, gtol, norm, eps, maxiter, disp, return_all, finite_diff_rel_step, **unknown_options)
1306 pk = -np.dot(Hk, gfk)
1307 try:
1308 alpha_k, fc, gc, old_fval, old_old_fval, gfkp1 = \
-> 1309 _line_search_wolfe12(f, myfprime, xk, pk, gfk,
1310 old_fval, old_old_fval, amin=1e-100, amax=1e100)
1311 except _LineSearchError:
1312 # Line search failed to find a better solution.
1313 warnflag = 2
File /opt/conda/anaconda3/envs/pisco/lib/python3.10/site-packages/scipy/optimize/_optimize.py:1087, in _line_search_wolfe12(f, fprime, xk, pk, gfk, old_fval, old_old_fval, **kwargs)
1073 """
1074 Same as line_search_wolfe1, but fall back to line_search_wolfe2 if
1075 suitable step length is not found, and raise an exception if a
(...)
1082
1083 """
1085 extra_condition = kwargs.pop('extra_condition', None)
-> 1087 ret = line_search_wolfe1(f, fprime, xk, pk, gfk,
1088 old_fval, old_old_fval,
1089 **kwargs)
1091 if ret[0] is not None and extra_condition is not None:
1092 xp1 = xk + ret[0] * pk
File /opt/conda/anaconda3/envs/pisco/lib/python3.10/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 /opt/conda/anaconda3/envs/pisco/lib/python3.10/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 /opt/conda/anaconda3/envs/pisco/lib/python3.10/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 /opt/conda/anaconda3/envs/pisco/lib/python3.10/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 /opt/conda/anaconda3/envs/pisco/lib/python3.10/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 /opt/conda/anaconda3/envs/pisco/lib/python3.10/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 /opt/conda/anaconda3/envs/pisco/lib/python3.10/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 /opt/conda/anaconda3/envs/pisco/lib/python3.10/site-packages/piscola/gaussian_process.py:142, in gp_2d_fit.<locals>.neg_ln_like(p)
140 def neg_ln_like(p):
141 gp.set_parameter_vector(p)
--> 142 return -gp.log_likelihood(y)
File /opt/conda/anaconda3/envs/pisco/lib/python3.10/site-packages/george/gp.py:385, in GP.log_likelihood(self, y, quiet)
369 def log_likelihood(self, y, quiet=False):
370 """
371 Compute the logarithm of the marginalized likelihood of a set of
372 observations under the Gaussian process model. You must call
(...)
383
384 """
--> 385 if not self.recompute(quiet=quiet):
386 return -np.inf
387 try:
File /opt/conda/anaconda3/envs/pisco/lib/python3.10/site-packages/george/gp.py:355, in GP.recompute(self, quiet, **kwargs)
351 raise RuntimeError("You need to compute the model first")
352 try:
353 # Update the model making sure that we store the original
354 # ordering of the points.
--> 355 self.compute(self._x, np.sqrt(self._yerr2), **kwargs)
356 except (ValueError, LinAlgError):
357 if quiet:
File /opt/conda/anaconda3/envs/pisco/lib/python3.10/site-packages/george/gp.py:331, in GP.compute(self, x, yerr, **kwargs)
329 # Include the white noise term.
330 yerr = np.sqrt(self._yerr2 + np.exp(self._call_white_noise(self._x)))
--> 331 self.solver.compute(self._x, yerr, **kwargs)
333 self._const = -0.5 * (
334 len(self._x) * np.log(2 * np.pi) + self.solver.log_determinant
335 )
336 self.computed = True
File /opt/conda/anaconda3/envs/pisco/lib/python3.10/site-packages/george/solvers/basic.py:68, in BasicSolver.compute(self, x, yerr)
65 K[np.diag_indices_from(K)] += yerr ** 2
67 # Factor the matrix and compute the log-determinant.
---> 68 self._factor = (cholesky(K, overwrite_a=True, lower=False), False)
69 self.log_determinant = 2 * np.sum(np.log(np.diag(self._factor[0])))
70 self.computed = True
File /opt/conda/anaconda3/envs/pisco/lib/python3.10/site-packages/scipy/linalg/_decomp_cholesky.py:88, in cholesky(a, lower, overwrite_a, check_finite)
45 def cholesky(a, lower=False, overwrite_a=False, check_finite=True):
46 """
47 Compute the Cholesky decomposition of a matrix.
48
(...)
86
87 """
---> 88 c, lower = _cholesky(a, lower=lower, overwrite_a=overwrite_a, clean=True,
89 check_finite=check_finite)
90 return c
File /opt/conda/anaconda3/envs/pisco/lib/python3.10/site-packages/scipy/linalg/_decomp_cholesky.py:37, in _cholesky(a, lower, overwrite_a, clean, check_finite)
35 c, info = potrf(a1, lower=lower, overwrite_a=overwrite_a, clean=clean)
36 if info > 0:
---> 37 raise LinAlgError("%d-th leading minor of the array is not positive "
38 "definite" % info)
39 if info < 0:
40 raise ValueError('LAPACK reported an illegal value in {}-th argument'
41 'on entry to "POTRF".'.format(-info))
LinAlgError: 104-th leading minor of the array is not positive definite