Generated by Cython 3.3.0
Yellow lines hint at Python interaction.
Click on a line that starts with a "+" to see the C code that Cython generated for it.
Raw output: _core.c
+001: # cython: boundscheck=False, wraparound=False, cdivision=True, language_level=3
__pyx_t_4 = __Pyx_PyDict_NewPresized(0); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 1, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_test, __pyx_t_4) < (0)) __PYX_ERR(0, 1, __pyx_L1_error) __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
002: # cython: nonecheck=False, initializedcheck=False
003: """
004: Cython extensions for mfe.multivariate recursion hot paths.
005:
006: Functions
007: ---------
008: _dcc_q_recursion DCC Q-process: (T, K, K) array of Q_t matrices
009: _bekk_scalar_recursion Scalar BEKK: (T, K, K) conditional covariances + log-lik
010: _bekk_diagonal_recursion Diagonal BEKK: same signature
011: _mvnorm_loglik Multivariate normal log-likelihood (vectorized, no alloc)
012: """
013:
+014: import numpy as np
__pyx_t_1 = __Pyx_Import(__pyx_mstate_global->__pyx_n_u_numpy, 0, 0, NULL, 0); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 14, __pyx_L1_error) __pyx_t_4 = __pyx_t_1; __Pyx_GOTREF(__pyx_t_4); if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_np, __pyx_t_4) < (0)) __PYX_ERR(0, 14, __pyx_L1_error) __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
015: cimport numpy as np
016: from numpy cimport ndarray, float64_t
017: from libc.math cimport log, sqrt, fabs
018:
019:
020: # ────────────────────────────────────────────────────────────────────────────
021: # Helpers
022: # ────────────────────────────────────────────────────────────────────────────
023:
+024: cdef double _logdet_2x2(double a, double b, double c, double d) nogil:
static double __pyx_f_3mfe_12multivariate_5_core__logdet_2x2(double __pyx_v_a, double __pyx_v_b, double __pyx_v_c, double __pyx_v_d) {
double __pyx_v_det;
double __pyx_r;
/* … */
/* function exit code */
__pyx_L0:;
return __pyx_r;
}
025: """log|[[a,b],[c,d]]| = log(ad - bc). Returns -inf on singularity."""
+026: cdef double det = a * d - b * c
__pyx_v_det = ((__pyx_v_a * __pyx_v_d) - (__pyx_v_b * __pyx_v_c));
+027: if det <= 0.0:
__pyx_t_1 = (__pyx_v_det <= 0.0);
if (__pyx_t_1) {
/* … */
}
+028: return -1e300
{
__pyx_r = -1e300;
}
goto __pyx_L0;
+029: return log(det)
{
__pyx_r = log(__pyx_v_det);
}
goto __pyx_L0;
030:
031:
+032: cdef double _quad_2x2(
static double __pyx_f_3mfe_12multivariate_5_core__quad_2x2(double __pyx_v_a, double __pyx_v_b, CYTHON_UNUSED double __pyx_v_c, double __pyx_v_d, double __pyx_v_x0, double __pyx_v_x1) {
double __pyx_v_det;
double __pyx_v_inv_det;
double __pyx_r;
/* … */
/* function exit code */
__pyx_L0:;
return __pyx_r;
}
033: double a, double b, double c, double d, # matrix [[a,b],[c,d]]
034: double x0, double x1, # vector
035: ) nogil:
036: """x' * inv([[a,b],[c,d]]) * x for 2x2 symmetric matrix (b == c)."""
+037: cdef double det = a * d - b * b
__pyx_v_det = ((__pyx_v_a * __pyx_v_d) - (__pyx_v_b * __pyx_v_b));
+038: if fabs(det) < 1e-300:
__pyx_t_1 = (fabs(__pyx_v_det) < 1e-300);
if (__pyx_t_1) {
/* … */
}
+039: return 1e300
{
__pyx_r = 1e300;
}
goto __pyx_L0;
+040: cdef double inv_det = 1.0 / det
__pyx_v_inv_det = (1.0 / __pyx_v_det);
041: # inv = [[d, -b],[-b, a]] / det
+042: return inv_det * (d * x0 * x0 - 2.0 * b * x0 * x1 + a * x1 * x1)
{
__pyx_r = (__pyx_v_inv_det * ((((__pyx_v_d * __pyx_v_x0) * __pyx_v_x0) - (((2.0 * __pyx_v_b) * __pyx_v_x0) * __pyx_v_x1)) + ((__pyx_v_a * __pyx_v_x1) * __pyx_v_x1)));
}
goto __pyx_L0;
043:
044:
045: # ────────────────────────────────────────────────────────────────────────────
046: # 1. DCC Q-recursion
047: # ────────────────────────────────────────────────────────────────────────────
048:
+049: def _dcc_q_recursion(
/* Python wrapper */ static PyObject *__pyx_pw_3mfe_12multivariate_5_core_1_dcc_q_recursion(PyObject *__pyx_self, #if CYTHON_VECTORCALL PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds #else PyObject *__pyx_args, PyObject *__pyx_kwds #endif ); /*proto*/ PyDoc_STRVAR(__pyx_doc_3mfe_12multivariate_5_core__dcc_q_recursion, "\n Q_t = (1 - a - b) * Q_bar + a * z_{t-1} z_{t-1}\047 + b * Q_{t-1}\n\n Returns (T, K, K) array of Q_t matrices.\n ~15x faster than the pure-Python loop for K=5, T=2000.\n "); static PyMethodDef __pyx_mdef_3mfe_12multivariate_5_core_1_dcc_q_recursion = {"_dcc_q_recursion", (PyCFunction)(void(*)(void))(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_3mfe_12multivariate_5_core_1_dcc_q_recursion, __Pyx_METH_FASTCALL|METH_KEYWORDS, __pyx_doc_3mfe_12multivariate_5_core__dcc_q_recursion}; static PyObject *__pyx_pw_3mfe_12multivariate_5_core_1_dcc_q_recursion(PyObject *__pyx_self, #if CYTHON_VECTORCALL PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds #else PyObject *__pyx_args, PyObject *__pyx_kwds #endif ) { __Pyx_memviewslice __pyx_v_z = { 0, 0, { 0 }, { 0 }, { 0 } }; __Pyx_memviewslice __pyx_v_q_bar = { 0, 0, { 0 }, { 0 }, { 0 } }; double __pyx_v_a; double __pyx_v_b; #if !CYTHON_VECTORCALL CYTHON_UNUSED Py_ssize_t __pyx_nargs; #endif CYTHON_UNUSED PyObject *const *__pyx_kwvalues; PyObject *__pyx_r = 0; __Pyx_RefNannyDeclarations __Pyx_RefNannySetupContext("_dcc_q_recursion (wrapper)", 0); #if !CYTHON_VECTORCALL #if CYTHON_ASSUME_SAFE_SIZE __pyx_nargs = PyTuple_GET_SIZE(__pyx_args); #else __pyx_nargs = PyTuple_Size(__pyx_args); if (unlikely(__pyx_nargs < 0)) return NULL; #endif #endif __pyx_kwvalues = __Pyx_KwValues_FASTCALL(__pyx_args, __pyx_nargs); { PyObject ** const __pyx_pyargnames[] = {&__pyx_mstate_global->__pyx_n_u_z,&__pyx_mstate_global->__pyx_n_u_q_bar,&__pyx_mstate_global->__pyx_n_u_a,&__pyx_mstate_global->__pyx_n_u_b,0}; PyObject* values[4] = {0,0,0,0}; const Py_ssize_t __pyx_kwds_len = (__pyx_kwds) ? __Pyx_NumKwargs_FASTCALL(__pyx_kwds) : 0; if (unlikely(__pyx_kwds_len < 0)) __PYX_ERR(0, 49, __pyx_L3_error) if (__pyx_kwds_len > 0) { switch (__pyx_nargs) { case 4: values[3] = __Pyx_ArgRef_FASTCALL(__pyx_args, 3); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[3])) __PYX_ERR(0, 49, __pyx_L3_error) CYTHON_FALLTHROUGH; case 3: values[2] = __Pyx_ArgRef_FASTCALL(__pyx_args, 2); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[2])) __PYX_ERR(0, 49, __pyx_L3_error) CYTHON_FALLTHROUGH; case 2: values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 49, __pyx_L3_error) CYTHON_FALLTHROUGH; case 1: values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 49, __pyx_L3_error) CYTHON_FALLTHROUGH; case 0: break; default: goto __pyx_L5_argtuple_error; } const Py_ssize_t kwd_pos_args = __pyx_nargs; if (__Pyx_ParseKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values, kwd_pos_args, __pyx_kwds_len, "_dcc_q_recursion", 0) < (0)) __PYX_ERR(0, 49, __pyx_L3_error) for (Py_ssize_t i = __pyx_nargs; i < 4; i++) { if (unlikely(!values[i])) { __Pyx_RaiseArgtupleInvalid("_dcc_q_recursion", 1, 4, 4, i); __PYX_ERR(0, 49, __pyx_L3_error) } } } else if (unlikely(__pyx_nargs != 4)) { goto __pyx_L5_argtuple_error; } else { values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 49, __pyx_L3_error) values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 49, __pyx_L3_error) values[2] = __Pyx_ArgRef_FASTCALL(__pyx_args, 2); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[2])) __PYX_ERR(0, 49, __pyx_L3_error) values[3] = __Pyx_ArgRef_FASTCALL(__pyx_args, 3); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[3])) __PYX_ERR(0, 49, __pyx_L3_error) } __pyx_v_z = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_z.memview)) __PYX_ERR(0, 50, __pyx_L3_error) __pyx_v_q_bar = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[1], PyBUF_WRITABLE); if (unlikely(!__pyx_v_q_bar.memview)) __PYX_ERR(0, 51, __pyx_L3_error) __pyx_v_a = __Pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_a == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 52, __pyx_L3_error) __pyx_v_b = __Pyx_PyFloat_AsDouble(values[3]); if (unlikely((__pyx_v_b == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 53, __pyx_L3_error) } goto __pyx_L6_skip; __pyx_L5_argtuple_error:; __Pyx_RaiseArgtupleInvalid("_dcc_q_recursion", 1, 4, 4, __pyx_nargs); __PYX_ERR(0, 49, __pyx_L3_error) __pyx_L6_skip:; goto __pyx_L4_argument_unpacking_done; __pyx_L3_error:; for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) { Py_XDECREF(values[__pyx_temp]); } __PYX_XCLEAR_MEMVIEW(&__pyx_v_z, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_q_bar, 1); __Pyx_AddTraceback("mfe.multivariate._core._dcc_q_recursion", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); return NULL; __pyx_L4_argument_unpacking_done:; __pyx_r = __pyx_pf_3mfe_12multivariate_5_core__dcc_q_recursion(__pyx_self, __pyx_v_z, __pyx_v_q_bar, __pyx_v_a, __pyx_v_b); int __pyx_lineno = 0; const char *__pyx_filename = NULL; int __pyx_clineno = 0; /* function exit code */ for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) { Py_XDECREF(values[__pyx_temp]); } __PYX_XCLEAR_MEMVIEW(&__pyx_v_z, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_q_bar, 1); __Pyx_RefNannyFinishContext(); return __pyx_r; } static PyObject *__pyx_pf_3mfe_12multivariate_5_core__dcc_q_recursion(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_z, __Pyx_memviewslice __pyx_v_q_bar, double __pyx_v_a, double __pyx_v_b) { int __pyx_v_T; int __pyx_v_K; int __pyx_v_t; int __pyx_v_i; int __pyx_v_j; double __pyx_v_c; PyObject *__pyx_v_Q = NULL; __Pyx_memviewslice __pyx_v_Qv = { 0, 0, { 0 }, { 0 }, { 0 } }; PyObject *__pyx_r = NULL; /* … */ /* function exit code */ __pyx_L1_error:; __Pyx_XDECREF(__pyx_t_1); __Pyx_XDECREF(__pyx_t_2); __Pyx_XDECREF(__pyx_t_3); __Pyx_XDECREF(__pyx_t_4); __Pyx_XDECREF(__pyx_t_5); __Pyx_XDECREF(__pyx_t_6); __Pyx_XDECREF(__pyx_t_7); __PYX_XCLEAR_MEMVIEW(&__pyx_t_9, 1); __Pyx_AddTraceback("mfe.multivariate._core._dcc_q_recursion", __pyx_clineno, __pyx_lineno, __pyx_filename); __pyx_r = NULL; __pyx_L0:; __Pyx_XDECREF(__pyx_v_Q); __PYX_XCLEAR_MEMVIEW(&__pyx_v_Qv, 1); __Pyx_XGIVEREF(__pyx_r); __Pyx_RefNannyFinishContext(); return __pyx_r; } /* … */ __pyx_t_4 = __Pyx_CyFunction_New(&__pyx_mdef_3mfe_12multivariate_5_core_1_dcc_q_recursion, 0, __pyx_mstate_global->__pyx_n_u_dcc_q_recursion, NULL, __pyx_mstate_global->__pyx_n_u_mfe_multivariate__core, __pyx_mstate_global->__pyx_d, ((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[0])); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 49, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); #if CYTHON_COMPILING_IN_CPYTHON && PY_VERSION_HEX >= 0x030E0000 PyUnstable_Object_EnableDeferredRefcount(__pyx_t_4); #endif if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_dcc_q_recursion, __pyx_t_4) < (0)) __PYX_ERR(0, 49, __pyx_L1_error) __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
050: double[:, ::1] z, # (T, K) standardized residuals
051: double[:, ::1] q_bar, # (K, K) unconditional covariance
052: double a,
053: double b,
054: ):
055: """
056: Q_t = (1 - a - b) * Q_bar + a * z_{t-1} z_{t-1}' + b * Q_{t-1}
057:
058: Returns (T, K, K) array of Q_t matrices.
059: ~15x faster than the pure-Python loop for K=5, T=2000.
060: """
+061: cdef int T = z.shape[0]
__pyx_v_T = (__pyx_v_z.shape[0]);
+062: cdef int K = z.shape[1]
__pyx_v_K = (__pyx_v_z.shape[1]);
063: cdef int t, i, j
+064: cdef double c = 1.0 - a - b
__pyx_v_c = ((1.0 - __pyx_v_a) - __pyx_v_b);
065:
+066: Q = np.empty((T, K, K), dtype=np.float64)
__pyx_t_2 = NULL; __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 66, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_3); __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_empty); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 66, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0; __pyx_t_3 = __Pyx_PyLong_From_int(__pyx_v_T); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 66, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_3); __pyx_t_5 = __Pyx_PyLong_From_int(__pyx_v_K); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 66, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_5); __pyx_t_6 = __Pyx_PyLong_From_int(__pyx_v_K); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 66, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); __pyx_t_7 = PyTuple_New(3); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 66, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_7); __Pyx_GIVEREF(__pyx_t_3); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_7, 0, __pyx_t_3) != (0)) __PYX_ERR(0, 66, __pyx_L1_error); __Pyx_GIVEREF(__pyx_t_5); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_7, 1, __pyx_t_5) != (0)) __PYX_ERR(0, 66, __pyx_L1_error); __Pyx_GIVEREF(__pyx_t_6); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_7, 2, __pyx_t_6) != (0)) __PYX_ERR(0, 66, __pyx_L1_error); __pyx_t_3 = 0; __pyx_t_5 = 0; __pyx_t_6 = 0; __Pyx_GetModuleGlobalName(__pyx_t_6, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 66, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_mstate_global->__pyx_n_u_float64); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 66, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_5); __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; __pyx_t_8 = 1; #if CYTHON_UNPACK_METHODS if (unlikely(PyMethod_Check(__pyx_t_4))) { __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_4); assert(__pyx_t_2); PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_4); __Pyx_INCREF(__pyx_t_2); __Pyx_INCREF(__pyx__function); __Pyx_DECREF_SET(__pyx_t_4, __pyx__function); __pyx_t_8 = 0; } #endif { PyObject *__pyx_callargs[3] = {__pyx_t_2, __pyx_t_7, __pyx_t_5}; #if CYTHON_VECTORCALL __pyx_t_6 = __pyx_mstate_global->__pyx_tuple[2]; if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 66, __pyx_L1_error) __Pyx_INCREF(__pyx_t_6); #else { PyObject *__pyx_temp[1] = {__pyx_mstate_global->__pyx_n_u_dtype}; __pyx_t_6 = __Pyx_MakeKwargDict(__pyx_temp, __pyx_callargs+2, 1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 66, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); } #endif __pyx_t_1 = __Pyx_Object_VectorcallKwds((PyObject*)__pyx_t_4, __pyx_callargs+__pyx_t_8, (2-__pyx_t_8) | (__pyx_t_8*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET), __pyx_t_6); __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0; __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0; __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0; __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 66, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); } __pyx_v_Q = __pyx_t_1; __pyx_t_1 = 0; /* … */ { PyObject* __pyx_temp[1] = {__pyx_mstate_global->__pyx_n_u_dtype}; __pyx_mstate_global->__pyx_tuple[2] = __Pyx_PyTuple_FromArray(__pyx_temp, 1); if (unlikely(!__pyx_mstate_global->__pyx_tuple[2])) __PYX_ERR(0, 66, __pyx_L1_error) __Pyx_GOTREF(__pyx_mstate_global->__pyx_tuple[2]); } __Pyx_GIVEREF(__pyx_mstate_global->__pyx_tuple[2]);
+067: cdef double[:, :, ::1] Qv = Q
__pyx_t_9 = __Pyx_PyObject_to_MemoryviewSlice_d_d_dc_double(__pyx_v_Q, PyBUF_WRITABLE); if (unlikely(!__pyx_t_9.memview)) __PYX_ERR(0, 67, __pyx_L1_error) __pyx_v_Qv = __pyx_t_9; __pyx_t_9.memview = NULL; __pyx_t_9.data = NULL;
068:
069: # Q[0] = q_bar
+070: for i in range(K):
__pyx_t_10 = __pyx_v_K;
__pyx_t_11 = __pyx_t_10;
for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
__pyx_v_i = __pyx_t_12;
+071: for j in range(K):
__pyx_t_13 = __pyx_v_K;
__pyx_t_14 = __pyx_t_13;
for (__pyx_t_15 = 0; __pyx_t_15 < __pyx_t_14; __pyx_t_15+=1) {
__pyx_v_j = __pyx_t_15;
+072: Qv[0, i, j] = q_bar[i, j]
__pyx_t_16 = __pyx_v_i;
__pyx_t_17 = __pyx_v_j;
__pyx_t_18 = 0;
__pyx_t_19 = __pyx_v_i;
__pyx_t_20 = __pyx_v_j;
*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Qv.data + __pyx_t_18 * __pyx_v_Qv.strides[0]) ) + __pyx_t_19 * __pyx_v_Qv.strides[1]) )) + __pyx_t_20)) )) = (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_q_bar.data + __pyx_t_16 * __pyx_v_q_bar.strides[0]) )) + __pyx_t_17)) )));
}
}
073:
+074: for t in range(1, T):
__pyx_t_10 = __pyx_v_T;
__pyx_t_11 = __pyx_t_10;
for (__pyx_t_12 = 1; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
__pyx_v_t = __pyx_t_12;
+075: for i in range(K):
__pyx_t_13 = __pyx_v_K;
__pyx_t_14 = __pyx_t_13;
for (__pyx_t_15 = 0; __pyx_t_15 < __pyx_t_14; __pyx_t_15+=1) {
__pyx_v_i = __pyx_t_15;
+076: for j in range(K):
__pyx_t_21 = __pyx_v_K;
__pyx_t_22 = __pyx_t_21;
for (__pyx_t_23 = 0; __pyx_t_23 < __pyx_t_22; __pyx_t_23+=1) {
__pyx_v_j = __pyx_t_23;
+077: Qv[t, i, j] = (
__pyx_t_28 = __pyx_v_t;
__pyx_t_29 = __pyx_v_i;
__pyx_t_30 = __pyx_v_j;
*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Qv.data + __pyx_t_28 * __pyx_v_Qv.strides[0]) ) + __pyx_t_29 * __pyx_v_Qv.strides[1]) )) + __pyx_t_30)) )) = (((__pyx_v_c * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_q_bar.data + __pyx_t_17 * __pyx_v_q_bar.strides[0]) )) + __pyx_t_16)) )))) + ((__pyx_v_a * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_z.data + __pyx_t_20 * __pyx_v_z.strides[0]) )) + __pyx_t_19)) )))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_z.data + __pyx_t_18 * __pyx_v_z.strides[0]) )) + __pyx_t_24)) ))))) + (__pyx_v_b * (*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Qv.data + __pyx_t_25 * __pyx_v_Qv.strides[0]) ) + __pyx_t_26 * __pyx_v_Qv.strides[1]) )) + __pyx_t_27)) )))));
}
}
}
+078: c * q_bar[i, j]
__pyx_t_17 = __pyx_v_i;
__pyx_t_16 = __pyx_v_j;
+079: + a * z[t - 1, i] * z[t - 1, j]
__pyx_t_20 = (__pyx_v_t - 1);
__pyx_t_19 = __pyx_v_i;
__pyx_t_18 = (__pyx_v_t - 1);
__pyx_t_24 = __pyx_v_j;
+080: + b * Qv[t - 1, i, j]
__pyx_t_25 = (__pyx_v_t - 1);
__pyx_t_26 = __pyx_v_i;
__pyx_t_27 = __pyx_v_j;
081: )
082:
+083: return Q
{
PyObject *__pyx_temp;
{
__pyx_temp = __pyx_r;
__Pyx_INCREF(__pyx_v_Q);
__pyx_r = __pyx_v_Q;
}
__Pyx_XDECREF(__pyx_temp);
}
goto __pyx_L0;
084:
085:
086: # ────────────────────────────────────────────────────────────────────────────
087: # 2. Scalar BEKK recursion + log-likelihood
088: # ────────────────────────────────────────────────────────────────────────────
089:
+090: def _bekk_scalar_recursion(
/* Python wrapper */ static PyObject *__pyx_pw_3mfe_12multivariate_5_core_3_bekk_scalar_recursion(PyObject *__pyx_self, #if CYTHON_VECTORCALL PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds #else PyObject *__pyx_args, PyObject *__pyx_kwds #endif ); /*proto*/ PyDoc_STRVAR(__pyx_doc_3mfe_12multivariate_5_core_2_bekk_scalar_recursion, "\n H_t = CC + a2 * eps_{t-1} eps_{t-1}\047 + b2 * H_{t-1}\n\n Returns (H_series, log_likelihood) where H_series is (T, K, K).\n Log-likelihood is Gaussian QML: sum_t [-0.5*(K*log2pi + log|H_t| + eps_t\047H_t^{-1}eps_t)]\n\n For K=2 uses an analytic 2x2 inverse/det (no LAPACK call).\n For K>2 falls back to numpy for the per-step inverse.\n "); static PyMethodDef __pyx_mdef_3mfe_12multivariate_5_core_3_bekk_scalar_recursion = {"_bekk_scalar_recursion", (PyCFunction)(void(*)(void))(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_3mfe_12multivariate_5_core_3_bekk_scalar_recursion, __Pyx_METH_FASTCALL|METH_KEYWORDS, __pyx_doc_3mfe_12multivariate_5_core_2_bekk_scalar_recursion}; static PyObject *__pyx_pw_3mfe_12multivariate_5_core_3_bekk_scalar_recursion(PyObject *__pyx_self, #if CYTHON_VECTORCALL PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds #else PyObject *__pyx_args, PyObject *__pyx_kwds #endif ) { __Pyx_memviewslice __pyx_v_eps = { 0, 0, { 0 }, { 0 }, { 0 } }; __Pyx_memviewslice __pyx_v_CC = { 0, 0, { 0 }, { 0 }, { 0 } }; double __pyx_v_a2; double __pyx_v_b2; __Pyx_memviewslice __pyx_v_H0 = { 0, 0, { 0 }, { 0 }, { 0 } }; #if !CYTHON_VECTORCALL CYTHON_UNUSED Py_ssize_t __pyx_nargs; #endif CYTHON_UNUSED PyObject *const *__pyx_kwvalues; PyObject *__pyx_r = 0; __Pyx_RefNannyDeclarations __Pyx_RefNannySetupContext("_bekk_scalar_recursion (wrapper)", 0); #if !CYTHON_VECTORCALL #if CYTHON_ASSUME_SAFE_SIZE __pyx_nargs = PyTuple_GET_SIZE(__pyx_args); #else __pyx_nargs = PyTuple_Size(__pyx_args); if (unlikely(__pyx_nargs < 0)) return NULL; #endif #endif __pyx_kwvalues = __Pyx_KwValues_FASTCALL(__pyx_args, __pyx_nargs); { PyObject ** const __pyx_pyargnames[] = {&__pyx_mstate_global->__pyx_n_u_eps,&__pyx_mstate_global->__pyx_n_u_CC,&__pyx_mstate_global->__pyx_n_u_a2,&__pyx_mstate_global->__pyx_n_u_b2,&__pyx_mstate_global->__pyx_n_u_H0,0}; PyObject* values[5] = {0,0,0,0,0}; const Py_ssize_t __pyx_kwds_len = (__pyx_kwds) ? __Pyx_NumKwargs_FASTCALL(__pyx_kwds) : 0; if (unlikely(__pyx_kwds_len < 0)) __PYX_ERR(0, 90, __pyx_L3_error) if (__pyx_kwds_len > 0) { switch (__pyx_nargs) { case 5: values[4] = __Pyx_ArgRef_FASTCALL(__pyx_args, 4); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[4])) __PYX_ERR(0, 90, __pyx_L3_error) CYTHON_FALLTHROUGH; case 4: values[3] = __Pyx_ArgRef_FASTCALL(__pyx_args, 3); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[3])) __PYX_ERR(0, 90, __pyx_L3_error) CYTHON_FALLTHROUGH; case 3: values[2] = __Pyx_ArgRef_FASTCALL(__pyx_args, 2); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[2])) __PYX_ERR(0, 90, __pyx_L3_error) CYTHON_FALLTHROUGH; case 2: values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 90, __pyx_L3_error) CYTHON_FALLTHROUGH; case 1: values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 90, __pyx_L3_error) CYTHON_FALLTHROUGH; case 0: break; default: goto __pyx_L5_argtuple_error; } const Py_ssize_t kwd_pos_args = __pyx_nargs; if (__Pyx_ParseKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values, kwd_pos_args, __pyx_kwds_len, "_bekk_scalar_recursion", 0) < (0)) __PYX_ERR(0, 90, __pyx_L3_error) for (Py_ssize_t i = __pyx_nargs; i < 5; i++) { if (unlikely(!values[i])) { __Pyx_RaiseArgtupleInvalid("_bekk_scalar_recursion", 1, 5, 5, i); __PYX_ERR(0, 90, __pyx_L3_error) } } } else if (unlikely(__pyx_nargs != 5)) { goto __pyx_L5_argtuple_error; } else { values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 90, __pyx_L3_error) values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 90, __pyx_L3_error) values[2] = __Pyx_ArgRef_FASTCALL(__pyx_args, 2); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[2])) __PYX_ERR(0, 90, __pyx_L3_error) values[3] = __Pyx_ArgRef_FASTCALL(__pyx_args, 3); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[3])) __PYX_ERR(0, 90, __pyx_L3_error) values[4] = __Pyx_ArgRef_FASTCALL(__pyx_args, 4); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[4])) __PYX_ERR(0, 90, __pyx_L3_error) } __pyx_v_eps = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_eps.memview)) __PYX_ERR(0, 91, __pyx_L3_error) __pyx_v_CC = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[1], PyBUF_WRITABLE); if (unlikely(!__pyx_v_CC.memview)) __PYX_ERR(0, 92, __pyx_L3_error) __pyx_v_a2 = __Pyx_PyFloat_AsDouble(values[2]); if (unlikely((__pyx_v_a2 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 93, __pyx_L3_error) __pyx_v_b2 = __Pyx_PyFloat_AsDouble(values[3]); if (unlikely((__pyx_v_b2 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 94, __pyx_L3_error) __pyx_v_H0 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[4], PyBUF_WRITABLE); if (unlikely(!__pyx_v_H0.memview)) __PYX_ERR(0, 95, __pyx_L3_error) } goto __pyx_L6_skip; __pyx_L5_argtuple_error:; __Pyx_RaiseArgtupleInvalid("_bekk_scalar_recursion", 1, 5, 5, __pyx_nargs); __PYX_ERR(0, 90, __pyx_L3_error) __pyx_L6_skip:; goto __pyx_L4_argument_unpacking_done; __pyx_L3_error:; for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) { Py_XDECREF(values[__pyx_temp]); } __PYX_XCLEAR_MEMVIEW(&__pyx_v_eps, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_CC, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_H0, 1); __Pyx_AddTraceback("mfe.multivariate._core._bekk_scalar_recursion", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); return NULL; __pyx_L4_argument_unpacking_done:; __pyx_r = __pyx_pf_3mfe_12multivariate_5_core_2_bekk_scalar_recursion(__pyx_self, __pyx_v_eps, __pyx_v_CC, __pyx_v_a2, __pyx_v_b2, __pyx_v_H0); int __pyx_lineno = 0; const char *__pyx_filename = NULL; int __pyx_clineno = 0; /* function exit code */ for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) { Py_XDECREF(values[__pyx_temp]); } __PYX_XCLEAR_MEMVIEW(&__pyx_v_eps, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_CC, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_H0, 1); __Pyx_RefNannyFinishContext(); return __pyx_r; } static PyObject *__pyx_pf_3mfe_12multivariate_5_core_2_bekk_scalar_recursion(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_eps, __Pyx_memviewslice __pyx_v_CC, double __pyx_v_a2, double __pyx_v_b2, __Pyx_memviewslice __pyx_v_H0) { int __pyx_v_T; int __pyx_v_K; int __pyx_v_t; int __pyx_v_i; int __pyx_v_j; PyObject *__pyx_v_H = NULL; __Pyx_memviewslice __pyx_v_Hv = { 0, 0, { 0 }, { 0 }, { 0 } }; double __pyx_v_ll; double __pyx_v_LOG2PI; double __pyx_v_h00; double __pyx_v_h01; double __pyx_v_h11; double __pyx_v_det; double __pyx_v_e0; double __pyx_v_e1; double __pyx_v_logdet; double __pyx_v_quad; PyObject *__pyx_v_H_np = NULL; PyObject *__pyx_v_eps_np = NULL; PyObject *__pyx_v_Ht = NULL; PyObject *__pyx_v_sign = NULL; PyObject *__pyx_v_ldet = NULL; PyObject *__pyx_v_Hinv = NULL; PyObject *__pyx_v_et = NULL; PyObject *__pyx_r = NULL; /* … */ /* function exit code */ __pyx_L1_error:; __Pyx_XDECREF(__pyx_t_1); __Pyx_XDECREF(__pyx_t_2); __Pyx_XDECREF(__pyx_t_3); __Pyx_XDECREF(__pyx_t_4); __Pyx_XDECREF(__pyx_t_5); __Pyx_XDECREF(__pyx_t_6); __Pyx_XDECREF(__pyx_t_7); __PYX_XCLEAR_MEMVIEW(&__pyx_t_9, 1); __Pyx_AddTraceback("mfe.multivariate._core._bekk_scalar_recursion", __pyx_clineno, __pyx_lineno, __pyx_filename); __pyx_r = NULL; __pyx_L0:; __Pyx_XDECREF(__pyx_v_H); __PYX_XCLEAR_MEMVIEW(&__pyx_v_Hv, 1); __Pyx_XDECREF(__pyx_v_H_np); __Pyx_XDECREF(__pyx_v_eps_np); __Pyx_XDECREF(__pyx_v_Ht); __Pyx_XDECREF(__pyx_v_sign); __Pyx_XDECREF(__pyx_v_ldet); __Pyx_XDECREF(__pyx_v_Hinv); __Pyx_XDECREF(__pyx_v_et); __Pyx_XGIVEREF(__pyx_r); __Pyx_RefNannyFinishContext(); return __pyx_r; } /* … */ __pyx_t_4 = __Pyx_CyFunction_New(&__pyx_mdef_3mfe_12multivariate_5_core_3_bekk_scalar_recursion, 0, __pyx_mstate_global->__pyx_n_u_bekk_scalar_recursion, NULL, __pyx_mstate_global->__pyx_n_u_mfe_multivariate__core, __pyx_mstate_global->__pyx_d, ((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[1])); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 90, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); #if CYTHON_COMPILING_IN_CPYTHON && PY_VERSION_HEX >= 0x030E0000 PyUnstable_Object_EnableDeferredRefcount(__pyx_t_4); #endif if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_bekk_scalar_recursion, __pyx_t_4) < (0)) __PYX_ERR(0, 90, __pyx_L1_error) __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
091: double[:, ::1] eps, # (T, K) residuals
092: double[:, ::1] CC, # (K, K) precomputed C'C
093: double a2, # alpha^2
094: double b2, # beta^2
095: double[:, ::1] H0, # (K, K) initial covariance
096: ):
097: """
098: H_t = CC + a2 * eps_{t-1} eps_{t-1}' + b2 * H_{t-1}
099:
100: Returns (H_series, log_likelihood) where H_series is (T, K, K).
101: Log-likelihood is Gaussian QML: sum_t [-0.5*(K*log2pi + log|H_t| + eps_t'H_t^{-1}eps_t)]
102:
103: For K=2 uses an analytic 2x2 inverse/det (no LAPACK call).
104: For K>2 falls back to numpy for the per-step inverse.
105: """
+106: cdef int T = eps.shape[0]
__pyx_v_T = (__pyx_v_eps.shape[0]);
+107: cdef int K = eps.shape[1]
__pyx_v_K = (__pyx_v_eps.shape[1]);
108: cdef int t, i, j
109:
+110: H = np.empty((T, K, K), dtype=np.float64)
__pyx_t_2 = NULL; __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 110, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_3); __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_empty); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 110, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0; __pyx_t_3 = __Pyx_PyLong_From_int(__pyx_v_T); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 110, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_3); __pyx_t_5 = __Pyx_PyLong_From_int(__pyx_v_K); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 110, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_5); __pyx_t_6 = __Pyx_PyLong_From_int(__pyx_v_K); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 110, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); __pyx_t_7 = PyTuple_New(3); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 110, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_7); __Pyx_GIVEREF(__pyx_t_3); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_7, 0, __pyx_t_3) != (0)) __PYX_ERR(0, 110, __pyx_L1_error); __Pyx_GIVEREF(__pyx_t_5); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_7, 1, __pyx_t_5) != (0)) __PYX_ERR(0, 110, __pyx_L1_error); __Pyx_GIVEREF(__pyx_t_6); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_7, 2, __pyx_t_6) != (0)) __PYX_ERR(0, 110, __pyx_L1_error); __pyx_t_3 = 0; __pyx_t_5 = 0; __pyx_t_6 = 0; __Pyx_GetModuleGlobalName(__pyx_t_6, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 110, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_mstate_global->__pyx_n_u_float64); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 110, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_5); __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; __pyx_t_8 = 1; #if CYTHON_UNPACK_METHODS if (unlikely(PyMethod_Check(__pyx_t_4))) { __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_4); assert(__pyx_t_2); PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_4); __Pyx_INCREF(__pyx_t_2); __Pyx_INCREF(__pyx__function); __Pyx_DECREF_SET(__pyx_t_4, __pyx__function); __pyx_t_8 = 0; } #endif { PyObject *__pyx_callargs[3] = {__pyx_t_2, __pyx_t_7, __pyx_t_5}; #if CYTHON_VECTORCALL __pyx_t_6 = __pyx_mstate_global->__pyx_tuple[2]; if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 110, __pyx_L1_error) __Pyx_INCREF(__pyx_t_6); #else { PyObject *__pyx_temp[1] = {__pyx_mstate_global->__pyx_n_u_dtype}; __pyx_t_6 = __Pyx_MakeKwargDict(__pyx_temp, __pyx_callargs+2, 1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 110, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); } #endif __pyx_t_1 = __Pyx_Object_VectorcallKwds((PyObject*)__pyx_t_4, __pyx_callargs+__pyx_t_8, (2-__pyx_t_8) | (__pyx_t_8*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET), __pyx_t_6); __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0; __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0; __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0; __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 110, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); } __pyx_v_H = __pyx_t_1; __pyx_t_1 = 0;
+111: cdef double[:, :, ::1] Hv = H
__pyx_t_9 = __Pyx_PyObject_to_MemoryviewSlice_d_d_dc_double(__pyx_v_H, PyBUF_WRITABLE); if (unlikely(!__pyx_t_9.memview)) __PYX_ERR(0, 111, __pyx_L1_error) __pyx_v_Hv = __pyx_t_9; __pyx_t_9.memview = NULL; __pyx_t_9.data = NULL;
112:
113: # H[0] = H0
+114: for i in range(K):
__pyx_t_10 = __pyx_v_K;
__pyx_t_11 = __pyx_t_10;
for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
__pyx_v_i = __pyx_t_12;
+115: for j in range(K):
__pyx_t_13 = __pyx_v_K;
__pyx_t_14 = __pyx_t_13;
for (__pyx_t_15 = 0; __pyx_t_15 < __pyx_t_14; __pyx_t_15+=1) {
__pyx_v_j = __pyx_t_15;
+116: Hv[0, i, j] = H0[i, j]
__pyx_t_16 = __pyx_v_i;
__pyx_t_17 = __pyx_v_j;
__pyx_t_18 = 0;
__pyx_t_19 = __pyx_v_i;
__pyx_t_20 = __pyx_v_j;
*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Hv.data + __pyx_t_18 * __pyx_v_Hv.strides[0]) ) + __pyx_t_19 * __pyx_v_Hv.strides[1]) )) + __pyx_t_20)) )) = (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_H0.data + __pyx_t_16 * __pyx_v_H0.strides[0]) )) + __pyx_t_17)) )));
}
}
117:
+118: cdef double ll = 0.0
__pyx_v_ll = 0.0;
+119: cdef double LOG2PI = 1.8378770664093453 # log(2*pi)
__pyx_v_LOG2PI = 1.8378770664093453;
120:
121: # Specialised 2x2 path avoids per-step numpy calls
122: cdef double h00, h01, h11, det, inv_det, q00, q01, q10, q11
123: cdef double e0, e1, logdet, quad
124:
+125: for t in range(1, T):
__pyx_t_10 = __pyx_v_T;
__pyx_t_11 = __pyx_t_10;
for (__pyx_t_12 = 1; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
__pyx_v_t = __pyx_t_12;
126: # H_t = CC + a2 * eps[t-1] eps[t-1]' + b2 * H[t-1]
+127: for i in range(K):
__pyx_t_13 = __pyx_v_K;
__pyx_t_14 = __pyx_t_13;
for (__pyx_t_15 = 0; __pyx_t_15 < __pyx_t_14; __pyx_t_15+=1) {
__pyx_v_i = __pyx_t_15;
+128: for j in range(K):
__pyx_t_21 = __pyx_v_K;
__pyx_t_22 = __pyx_t_21;
for (__pyx_t_23 = 0; __pyx_t_23 < __pyx_t_22; __pyx_t_23+=1) {
__pyx_v_j = __pyx_t_23;
+129: Hv[t, i, j] = (
__pyx_t_28 = __pyx_v_t;
__pyx_t_29 = __pyx_v_i;
__pyx_t_30 = __pyx_v_j;
*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Hv.data + __pyx_t_28 * __pyx_v_Hv.strides[0]) ) + __pyx_t_29 * __pyx_v_Hv.strides[1]) )) + __pyx_t_30)) )) = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_CC.data + __pyx_t_17 * __pyx_v_CC.strides[0]) )) + __pyx_t_16)) ))) + ((__pyx_v_a2 * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_eps.data + __pyx_t_20 * __pyx_v_eps.strides[0]) )) + __pyx_t_19)) )))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_eps.data + __pyx_t_18 * __pyx_v_eps.strides[0]) )) + __pyx_t_24)) ))))) + (__pyx_v_b2 * (*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Hv.data + __pyx_t_25 * __pyx_v_Hv.strides[0]) ) + __pyx_t_26 * __pyx_v_Hv.strides[1]) )) + __pyx_t_27)) )))));
}
}
}
+130: CC[i, j]
__pyx_t_17 = __pyx_v_i;
__pyx_t_16 = __pyx_v_j;
+131: + a2 * eps[t - 1, i] * eps[t - 1, j]
__pyx_t_20 = (__pyx_v_t - 1);
__pyx_t_19 = __pyx_v_i;
__pyx_t_18 = (__pyx_v_t - 1);
__pyx_t_24 = __pyx_v_j;
+132: + b2 * Hv[t - 1, i, j]
__pyx_t_25 = (__pyx_v_t - 1);
__pyx_t_26 = __pyx_v_i;
__pyx_t_27 = __pyx_v_j;
133: )
134:
135: # Log-likelihood loop — K==2 fast path
+136: if K == 2:
__pyx_t_31 = (__pyx_v_K == 2);
if (__pyx_t_31) {
/* … */
goto __pyx_L13;
}
+137: for t in range(T):
__pyx_t_10 = __pyx_v_T;
__pyx_t_11 = __pyx_t_10;
for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
__pyx_v_t = __pyx_t_12;
+138: h00 = Hv[t, 0, 0]
__pyx_t_27 = __pyx_v_t;
__pyx_t_26 = 0;
__pyx_t_25 = 0;
__pyx_v_h00 = (*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Hv.data + __pyx_t_27 * __pyx_v_Hv.strides[0]) ) + __pyx_t_26 * __pyx_v_Hv.strides[1]) )) + __pyx_t_25)) )));
+139: h01 = Hv[t, 0, 1]
__pyx_t_25 = __pyx_v_t;
__pyx_t_26 = 0;
__pyx_t_27 = 1;
__pyx_v_h01 = (*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Hv.data + __pyx_t_25 * __pyx_v_Hv.strides[0]) ) + __pyx_t_26 * __pyx_v_Hv.strides[1]) )) + __pyx_t_27)) )));
+140: h11 = Hv[t, 1, 1]
__pyx_t_27 = __pyx_v_t;
__pyx_t_26 = 1;
__pyx_t_25 = 1;
__pyx_v_h11 = (*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Hv.data + __pyx_t_27 * __pyx_v_Hv.strides[0]) ) + __pyx_t_26 * __pyx_v_Hv.strides[1]) )) + __pyx_t_25)) )));
+141: det = h00 * h11 - h01 * h01
__pyx_v_det = ((__pyx_v_h00 * __pyx_v_h11) - (__pyx_v_h01 * __pyx_v_h01));
+142: if det <= 0.0:
__pyx_t_31 = (__pyx_v_det <= 0.0);
if (__pyx_t_31) {
/* … */
}
+143: return H, 1e10
__pyx_t_1 = PyTuple_New(2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 143, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __Pyx_INCREF(__pyx_v_H); __Pyx_GIVEREF(__pyx_v_H); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_v_H) != (0)) __PYX_ERR(0, 143, __pyx_L1_error); __Pyx_INCREF(__pyx_mstate_global->__pyx_float_1e10); __Pyx_GIVEREF(__pyx_mstate_global->__pyx_float_1e10); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_1, 1, __pyx_mstate_global->__pyx_float_1e10) != (0)) __PYX_ERR(0, 143, __pyx_L1_error); { PyObject *__pyx_temp; { __pyx_temp = __pyx_r; __pyx_r = __pyx_t_1; } __Pyx_XDECREF(__pyx_temp); } __pyx_t_1 = 0; goto __pyx_L0;
+144: e0 = eps[t, 0]
__pyx_t_25 = __pyx_v_t;
__pyx_t_26 = 0;
__pyx_v_e0 = (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_eps.data + __pyx_t_25 * __pyx_v_eps.strides[0]) )) + __pyx_t_26)) )));
+145: e1 = eps[t, 1]
__pyx_t_26 = __pyx_v_t;
__pyx_t_25 = 1;
__pyx_v_e1 = (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_eps.data + __pyx_t_26 * __pyx_v_eps.strides[0]) )) + __pyx_t_25)) )));
+146: logdet = log(det)
__pyx_v_logdet = log(__pyx_v_det);
147: # quad = e' H^{-1} e = (h11*e0^2 - 2*h01*e0*e1 + h00*e1^2) / det
+148: quad = (h11 * e0 * e0 - 2.0 * h01 * e0 * e1 + h00 * e1 * e1) / det
__pyx_v_quad = (((((__pyx_v_h11 * __pyx_v_e0) * __pyx_v_e0) - (((2.0 * __pyx_v_h01) * __pyx_v_e0) * __pyx_v_e1)) + ((__pyx_v_h00 * __pyx_v_e1) * __pyx_v_e1)) / __pyx_v_det);
+149: ll += logdet + quad
__pyx_v_ll = (__pyx_v_ll + (__pyx_v_logdet + __pyx_v_quad));
}
+150: ll = 0.5 * (T * K * LOG2PI + ll)
__pyx_v_ll = (0.5 * (((__pyx_v_T * __pyx_v_K) * __pyx_v_LOG2PI) + __pyx_v_ll));
151: else:
152: # General path: use numpy for inversion
+153: H_np = np.asarray(H)
/*else*/ {
__pyx_t_4 = NULL;
__Pyx_GetModuleGlobalName(__pyx_t_6, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 153, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_mstate_global->__pyx_n_u_asarray); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 153, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__pyx_t_8 = 1;
#if CYTHON_UNPACK_METHODS
if (unlikely(PyMethod_Check(__pyx_t_5))) {
__pyx_t_4 = PyMethod_GET_SELF(__pyx_t_5);
assert(__pyx_t_4);
PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_5);
__Pyx_INCREF(__pyx_t_4);
__Pyx_INCREF(__pyx__function);
__Pyx_DECREF_SET(__pyx_t_5, __pyx__function);
__pyx_t_8 = 0;
}
#endif
{
PyObject *__pyx_callargs[2] = {__pyx_t_4, __pyx_v_H};
__pyx_t_1 = __Pyx_PyObject_FastCall((PyObject*)__pyx_t_5, __pyx_callargs+__pyx_t_8, (2-__pyx_t_8) | (__pyx_t_8*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
__Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 153, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
}
__pyx_v_H_np = __pyx_t_1;
__pyx_t_1 = 0;
+154: eps_np = np.asarray(eps)
__pyx_t_5 = NULL;
__Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 154, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_mstate_global->__pyx_n_u_asarray); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 154, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_t_4 = __pyx_memoryview_fromslice(__pyx_v_eps, 2, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 154, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__pyx_t_8 = 1;
#if CYTHON_UNPACK_METHODS
if (unlikely(PyMethod_Check(__pyx_t_6))) {
__pyx_t_5 = PyMethod_GET_SELF(__pyx_t_6);
assert(__pyx_t_5);
PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_6);
__Pyx_INCREF(__pyx_t_5);
__Pyx_INCREF(__pyx__function);
__Pyx_DECREF_SET(__pyx_t_6, __pyx__function);
__pyx_t_8 = 0;
}
#endif
{
PyObject *__pyx_callargs[2] = {__pyx_t_5, __pyx_t_4};
__pyx_t_1 = __Pyx_PyObject_FastCall((PyObject*)__pyx_t_6, __pyx_callargs+__pyx_t_8, (2-__pyx_t_8) | (__pyx_t_8*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
__Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 154, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
}
__pyx_v_eps_np = __pyx_t_1;
__pyx_t_1 = 0;
+155: for t in range(T):
__pyx_t_10 = __pyx_v_T;
__pyx_t_11 = __pyx_t_10;
for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
__pyx_v_t = __pyx_t_12;
+156: Ht = H_np[t]
__pyx_t_1 = __Pyx_GetItemInt(__pyx_v_H_np, __pyx_v_t, int, 1, __Pyx_PyLong_From_int, 0, 0, 1, __Pyx_ReferenceSharing_OwnStrongReference); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 156, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __Pyx_XDECREF_SET(__pyx_v_Ht, __pyx_t_1); __pyx_t_1 = 0;
+157: sign, ldet = np.linalg.slogdet(Ht)
__Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 157, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_mstate_global->__pyx_n_u_linalg); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 157, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_5); __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; __pyx_t_6 = __pyx_t_5; __Pyx_INCREF(__pyx_t_6); __pyx_t_8 = 0; { PyObject *__pyx_callargs[2] = {__pyx_t_6, __pyx_v_Ht}; __pyx_t_1 = __Pyx_PyObject_FastCallMethod((PyObject*)__pyx_mstate_global->__pyx_n_u_slogdet, __pyx_callargs+__pyx_t_8, (2-__pyx_t_8) | (1*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET)); __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0; __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 157, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); } if ((likely(PyTuple_CheckExact(__pyx_t_1))) || (PyList_CheckExact(__pyx_t_1))) { PyObject* sequence = __pyx_t_1; Py_ssize_t size = __Pyx_PySequence_SIZE(sequence); if (unlikely(size != 2)) { if (size > 2) __Pyx_RaiseTooManyValuesError(2); else if (size >= 0) __Pyx_RaiseNeedMoreValuesError(size); __PYX_ERR(0, 157, __pyx_L1_error) } #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS if (likely(PyTuple_CheckExact(sequence))) { __pyx_t_5 = PyTuple_GET_ITEM(sequence, 0); __Pyx_INCREF(__pyx_t_5); __pyx_t_6 = PyTuple_GET_ITEM(sequence, 1); __Pyx_INCREF(__pyx_t_6); } else { __pyx_t_5 = __Pyx_PyList_GET_ITEM_REF(sequence, 0, __Pyx_ReferenceSharing_SharedReference); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 157, __pyx_L1_error) __Pyx_XGOTREF(__pyx_t_5); __pyx_t_6 = __Pyx_PyList_GET_ITEM_REF(sequence, 1, __Pyx_ReferenceSharing_SharedReference); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 157, __pyx_L1_error) __Pyx_XGOTREF(__pyx_t_6); } #else __pyx_t_5 = __Pyx_PySequence_ITEM(sequence, 0); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 157, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_5); __pyx_t_6 = __Pyx_PySequence_ITEM(sequence, 1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 157, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); #endif __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; } else { Py_ssize_t index = -1; __pyx_t_4 = PyObject_GetIter(__pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 157, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_t_32 = (CYTHON_COMPILING_IN_LIMITED_API) ? PyIter_Next : __Pyx_PyObject_GetIterNextFunc(__pyx_t_4); index = 0; __pyx_t_5 = __pyx_t_32(__pyx_t_4); if (unlikely(!__pyx_t_5)) goto __pyx_L19_unpacking_failed; __Pyx_GOTREF(__pyx_t_5); index = 1; __pyx_t_6 = __pyx_t_32(__pyx_t_4); if (unlikely(!__pyx_t_6)) goto __pyx_L19_unpacking_failed; __Pyx_GOTREF(__pyx_t_6); if (__Pyx_IternextUnpackEndCheck(__pyx_t_32(__pyx_t_4), 2) < (0)) __PYX_ERR(0, 157, __pyx_L1_error) __pyx_t_32 = NULL; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; goto __pyx_L20_unpacking_done; __pyx_L19_unpacking_failed:; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; __pyx_t_32 = NULL; if (__Pyx_IterFinish() == 0) __Pyx_RaiseNeedMoreValuesError(index); __PYX_ERR(0, 157, __pyx_L1_error) __pyx_L20_unpacking_done:; } __Pyx_XDECREF_SET(__pyx_v_sign, __pyx_t_5); __pyx_t_5 = 0; __Pyx_XDECREF_SET(__pyx_v_ldet, __pyx_t_6); __pyx_t_6 = 0;
+158: if sign <= 0:
__pyx_t_31 = __Pyx_PyObject_CompareBoolLe_object_int(__pyx_v_sign, __pyx_mstate_global->__pyx_int_0, Py_LE); if (unlikely((__pyx_t_31 < 0))) __PYX_ERR(0, 158, __pyx_L1_error) if (__pyx_t_31) { /* … */ }
+159: return H, 1e10
__pyx_t_1 = PyTuple_New(2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 159, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __Pyx_INCREF(__pyx_v_H); __Pyx_GIVEREF(__pyx_v_H); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_v_H) != (0)) __PYX_ERR(0, 159, __pyx_L1_error); __Pyx_INCREF(__pyx_mstate_global->__pyx_float_1e10); __Pyx_GIVEREF(__pyx_mstate_global->__pyx_float_1e10); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_1, 1, __pyx_mstate_global->__pyx_float_1e10) != (0)) __PYX_ERR(0, 159, __pyx_L1_error); { PyObject *__pyx_temp; { __pyx_temp = __pyx_r; __pyx_r = __pyx_t_1; } __Pyx_XDECREF(__pyx_temp); } __pyx_t_1 = 0; goto __pyx_L0;
+160: try:
{
/*try:*/ {
/* … */
}
__Pyx_XDECREF(__pyx_t_33); __pyx_t_33 = 0;
__Pyx_XDECREF(__pyx_t_34); __pyx_t_34 = 0;
__Pyx_XDECREF(__pyx_t_35); __pyx_t_35 = 0;
goto __pyx_L29_try_end;
__pyx_L22_error:;
__Pyx_XDECREF(__pyx_t_1); __pyx_t_1 = 0;
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
__Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0;
__Pyx_XDECREF(__pyx_t_7); __pyx_t_7 = 0;
__PYX_XCLEAR_MEMVIEW(&__pyx_t_9, 1);; __pyx_t_9.memview = NULL; __pyx_t_9.data = NULL;
/* … */
__pyx_L24_except_error:;
__Pyx_XGIVEREF(__pyx_t_33);
__Pyx_XGIVEREF(__pyx_t_34);
__Pyx_XGIVEREF(__pyx_t_35);
__Pyx_ExceptionReset(__pyx_t_33, __pyx_t_34, __pyx_t_35);
goto __pyx_L1_error;
__pyx_L25_except_return:;
__Pyx_XGIVEREF(__pyx_t_33);
__Pyx_XGIVEREF(__pyx_t_34);
__Pyx_XGIVEREF(__pyx_t_35);
__Pyx_ExceptionReset(__pyx_t_33, __pyx_t_34, __pyx_t_35);
goto __pyx_L0;
__pyx_L29_try_end:;
}
+161: Hinv = np.linalg.inv(Ht)
__Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 161, __pyx_L22_error) __Pyx_GOTREF(__pyx_t_5); __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_linalg); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 161, __pyx_L22_error) __Pyx_GOTREF(__pyx_t_4); __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0; __pyx_t_6 = __pyx_t_4; __Pyx_INCREF(__pyx_t_6); __pyx_t_8 = 0; { PyObject *__pyx_callargs[2] = {__pyx_t_6, __pyx_v_Ht}; __pyx_t_1 = __Pyx_PyObject_FastCallMethod((PyObject*)__pyx_mstate_global->__pyx_n_u_inv, __pyx_callargs+__pyx_t_8, (2-__pyx_t_8) | (1*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET)); __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 161, __pyx_L22_error) __Pyx_GOTREF(__pyx_t_1); } __Pyx_XDECREF_SET(__pyx_v_Hinv, __pyx_t_1); __pyx_t_1 = 0;
+162: except Exception:
__pyx_t_13 = __Pyx_PyErr_ExceptionMatches(((PyObject *)(((PyTypeObject*)PyExc_Exception)))); if (__pyx_t_13) { __Pyx_AddTraceback("mfe.multivariate._core._bekk_scalar_recursion", __pyx_clineno, __pyx_lineno, __pyx_filename); if (__Pyx_GetException(&__pyx_t_1, &__pyx_t_4, &__pyx_t_6) < 0) __PYX_ERR(0, 162, __pyx_L24_except_error) __Pyx_XGOTREF(__pyx_t_1); __Pyx_XGOTREF(__pyx_t_4); __Pyx_XGOTREF(__pyx_t_6);
+163: return H, 1e10
__pyx_t_5 = PyTuple_New(2); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 163, __pyx_L24_except_error) __Pyx_GOTREF(__pyx_t_5); __Pyx_INCREF(__pyx_v_H); __Pyx_GIVEREF(__pyx_v_H); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_v_H) != (0)) __PYX_ERR(0, 163, __pyx_L24_except_error); __Pyx_INCREF(__pyx_mstate_global->__pyx_float_1e10); __Pyx_GIVEREF(__pyx_mstate_global->__pyx_float_1e10); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_5, 1, __pyx_mstate_global->__pyx_float_1e10) != (0)) __PYX_ERR(0, 163, __pyx_L24_except_error); { PyObject *__pyx_temp; { __pyx_temp = __pyx_r; __pyx_r = __pyx_t_5; } __Pyx_XDECREF(__pyx_temp); } __pyx_t_5 = 0; __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; goto __pyx_L25_except_return; } goto __pyx_L24_except_error;
+164: et = eps_np[t]
__pyx_t_6 = __Pyx_GetItemInt(__pyx_v_eps_np, __pyx_v_t, int, 1, __Pyx_PyLong_From_int, 0, 0, 1, __Pyx_ReferenceSharing_OwnStrongReference); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 164, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); __Pyx_XDECREF_SET(__pyx_v_et, __pyx_t_6); __pyx_t_6 = 0;
+165: ll += float(ldet) + float(et @ Hinv @ et)
__pyx_t_36 = __Pyx_PyObject_AsDouble(__pyx_v_ldet); if (unlikely(__PYX_CHECK_FLOAT_EXCEPTION(__pyx_t_36, ((double)((double)-1))) && PyErr_Occurred())) __PYX_ERR(0, 165, __pyx_L1_error) __pyx_t_6 = __Pyx_PyNumber_MatrixMultiply(__pyx_v_et, __pyx_v_Hinv); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 165, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); __pyx_t_4 = __Pyx_PyNumber_MatrixMultiply(__pyx_t_6, __pyx_v_et); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 165, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; __pyx_t_37 = __Pyx_PyObject_AsDouble(__pyx_t_4); if (unlikely(__PYX_CHECK_FLOAT_EXCEPTION(__pyx_t_37, ((double)((double)-1))) && PyErr_Occurred())) __PYX_ERR(0, 165, __pyx_L1_error) __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; __pyx_v_ll = (__pyx_v_ll + (__pyx_t_36 + __pyx_t_37)); }
+166: ll = 0.5 * (T * K * LOG2PI + ll)
__pyx_v_ll = (0.5 * (((__pyx_v_T * __pyx_v_K) * __pyx_v_LOG2PI) + __pyx_v_ll)); } __pyx_L13:;
167:
+168: return H, ll
__pyx_t_4 = PyFloat_FromDouble(__pyx_v_ll); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 168, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); __pyx_t_6 = PyTuple_New(2); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 168, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); __Pyx_INCREF(__pyx_v_H); __Pyx_GIVEREF(__pyx_v_H); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_6, 0, __pyx_v_H) != (0)) __PYX_ERR(0, 168, __pyx_L1_error); __Pyx_GIVEREF(__pyx_t_4); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_6, 1, __pyx_t_4) != (0)) __PYX_ERR(0, 168, __pyx_L1_error); __pyx_t_4 = 0; { PyObject *__pyx_temp; { __pyx_temp = __pyx_r; __pyx_r = __pyx_t_6; } __Pyx_XDECREF(__pyx_temp); } __pyx_t_6 = 0; goto __pyx_L0;
169:
170:
171: # ────────────────────────────────────────────────────────────────────────────
172: # 3. Diagonal BEKK recursion + log-likelihood
173: # ────────────────────────────────────────────────────────────────────────────
174:
+175: def _bekk_diagonal_recursion(
/* Python wrapper */ static PyObject *__pyx_pw_3mfe_12multivariate_5_core_5_bekk_diagonal_recursion(PyObject *__pyx_self, #if CYTHON_VECTORCALL PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds #else PyObject *__pyx_args, PyObject *__pyx_kwds #endif ); /*proto*/ PyDoc_STRVAR(__pyx_doc_3mfe_12multivariate_5_core_4_bekk_diagonal_recursion, "\n H_t = CC + AtA * (eps_{t-1} eps_{t-1}\047) + BtB * H_{t-1}\n\n (element-wise products \342\200\224 A, B are diagonal so AtA[i,j] = A[i]*A[j])\n "); static PyMethodDef __pyx_mdef_3mfe_12multivariate_5_core_5_bekk_diagonal_recursion = {"_bekk_diagonal_recursion", (PyCFunction)(void(*)(void))(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_3mfe_12multivariate_5_core_5_bekk_diagonal_recursion, __Pyx_METH_FASTCALL|METH_KEYWORDS, __pyx_doc_3mfe_12multivariate_5_core_4_bekk_diagonal_recursion}; static PyObject *__pyx_pw_3mfe_12multivariate_5_core_5_bekk_diagonal_recursion(PyObject *__pyx_self, #if CYTHON_VECTORCALL PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds #else PyObject *__pyx_args, PyObject *__pyx_kwds #endif ) { __Pyx_memviewslice __pyx_v_eps = { 0, 0, { 0 }, { 0 }, { 0 } }; __Pyx_memviewslice __pyx_v_CC = { 0, 0, { 0 }, { 0 }, { 0 } }; __Pyx_memviewslice __pyx_v_AtA = { 0, 0, { 0 }, { 0 }, { 0 } }; __Pyx_memviewslice __pyx_v_BtB = { 0, 0, { 0 }, { 0 }, { 0 } }; __Pyx_memviewslice __pyx_v_H0 = { 0, 0, { 0 }, { 0 }, { 0 } }; #if !CYTHON_VECTORCALL CYTHON_UNUSED Py_ssize_t __pyx_nargs; #endif CYTHON_UNUSED PyObject *const *__pyx_kwvalues; PyObject *__pyx_r = 0; __Pyx_RefNannyDeclarations __Pyx_RefNannySetupContext("_bekk_diagonal_recursion (wrapper)", 0); #if !CYTHON_VECTORCALL #if CYTHON_ASSUME_SAFE_SIZE __pyx_nargs = PyTuple_GET_SIZE(__pyx_args); #else __pyx_nargs = PyTuple_Size(__pyx_args); if (unlikely(__pyx_nargs < 0)) return NULL; #endif #endif __pyx_kwvalues = __Pyx_KwValues_FASTCALL(__pyx_args, __pyx_nargs); { PyObject ** const __pyx_pyargnames[] = {&__pyx_mstate_global->__pyx_n_u_eps,&__pyx_mstate_global->__pyx_n_u_CC,&__pyx_mstate_global->__pyx_n_u_AtA,&__pyx_mstate_global->__pyx_n_u_BtB,&__pyx_mstate_global->__pyx_n_u_H0,0}; PyObject* values[5] = {0,0,0,0,0}; const Py_ssize_t __pyx_kwds_len = (__pyx_kwds) ? __Pyx_NumKwargs_FASTCALL(__pyx_kwds) : 0; if (unlikely(__pyx_kwds_len < 0)) __PYX_ERR(0, 175, __pyx_L3_error) if (__pyx_kwds_len > 0) { switch (__pyx_nargs) { case 5: values[4] = __Pyx_ArgRef_FASTCALL(__pyx_args, 4); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[4])) __PYX_ERR(0, 175, __pyx_L3_error) CYTHON_FALLTHROUGH; case 4: values[3] = __Pyx_ArgRef_FASTCALL(__pyx_args, 3); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[3])) __PYX_ERR(0, 175, __pyx_L3_error) CYTHON_FALLTHROUGH; case 3: values[2] = __Pyx_ArgRef_FASTCALL(__pyx_args, 2); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[2])) __PYX_ERR(0, 175, __pyx_L3_error) CYTHON_FALLTHROUGH; case 2: values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 175, __pyx_L3_error) CYTHON_FALLTHROUGH; case 1: values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 175, __pyx_L3_error) CYTHON_FALLTHROUGH; case 0: break; default: goto __pyx_L5_argtuple_error; } const Py_ssize_t kwd_pos_args = __pyx_nargs; if (__Pyx_ParseKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values, kwd_pos_args, __pyx_kwds_len, "_bekk_diagonal_recursion", 0) < (0)) __PYX_ERR(0, 175, __pyx_L3_error) for (Py_ssize_t i = __pyx_nargs; i < 5; i++) { if (unlikely(!values[i])) { __Pyx_RaiseArgtupleInvalid("_bekk_diagonal_recursion", 1, 5, 5, i); __PYX_ERR(0, 175, __pyx_L3_error) } } } else if (unlikely(__pyx_nargs != 5)) { goto __pyx_L5_argtuple_error; } else { values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 175, __pyx_L3_error) values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 175, __pyx_L3_error) values[2] = __Pyx_ArgRef_FASTCALL(__pyx_args, 2); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[2])) __PYX_ERR(0, 175, __pyx_L3_error) values[3] = __Pyx_ArgRef_FASTCALL(__pyx_args, 3); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[3])) __PYX_ERR(0, 175, __pyx_L3_error) values[4] = __Pyx_ArgRef_FASTCALL(__pyx_args, 4); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[4])) __PYX_ERR(0, 175, __pyx_L3_error) } __pyx_v_eps = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_eps.memview)) __PYX_ERR(0, 176, __pyx_L3_error) __pyx_v_CC = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[1], PyBUF_WRITABLE); if (unlikely(!__pyx_v_CC.memview)) __PYX_ERR(0, 177, __pyx_L3_error) __pyx_v_AtA = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[2], PyBUF_WRITABLE); if (unlikely(!__pyx_v_AtA.memview)) __PYX_ERR(0, 178, __pyx_L3_error) __pyx_v_BtB = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[3], PyBUF_WRITABLE); if (unlikely(!__pyx_v_BtB.memview)) __PYX_ERR(0, 179, __pyx_L3_error) __pyx_v_H0 = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[4], PyBUF_WRITABLE); if (unlikely(!__pyx_v_H0.memview)) __PYX_ERR(0, 180, __pyx_L3_error) } goto __pyx_L6_skip; __pyx_L5_argtuple_error:; __Pyx_RaiseArgtupleInvalid("_bekk_diagonal_recursion", 1, 5, 5, __pyx_nargs); __PYX_ERR(0, 175, __pyx_L3_error) __pyx_L6_skip:; goto __pyx_L4_argument_unpacking_done; __pyx_L3_error:; for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) { Py_XDECREF(values[__pyx_temp]); } __PYX_XCLEAR_MEMVIEW(&__pyx_v_eps, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_CC, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_AtA, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_BtB, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_H0, 1); __Pyx_AddTraceback("mfe.multivariate._core._bekk_diagonal_recursion", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); return NULL; __pyx_L4_argument_unpacking_done:; __pyx_r = __pyx_pf_3mfe_12multivariate_5_core_4_bekk_diagonal_recursion(__pyx_self, __pyx_v_eps, __pyx_v_CC, __pyx_v_AtA, __pyx_v_BtB, __pyx_v_H0); int __pyx_lineno = 0; const char *__pyx_filename = NULL; int __pyx_clineno = 0; /* function exit code */ for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) { Py_XDECREF(values[__pyx_temp]); } __PYX_XCLEAR_MEMVIEW(&__pyx_v_eps, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_CC, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_AtA, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_BtB, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_H0, 1); __Pyx_RefNannyFinishContext(); return __pyx_r; } static PyObject *__pyx_pf_3mfe_12multivariate_5_core_4_bekk_diagonal_recursion(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_eps, __Pyx_memviewslice __pyx_v_CC, __Pyx_memviewslice __pyx_v_AtA, __Pyx_memviewslice __pyx_v_BtB, __Pyx_memviewslice __pyx_v_H0) { int __pyx_v_T; int __pyx_v_K; int __pyx_v_t; int __pyx_v_i; int __pyx_v_j; PyObject *__pyx_v_H = NULL; __Pyx_memviewslice __pyx_v_Hv = { 0, 0, { 0 }, { 0 }, { 0 } }; double __pyx_v_ll; double __pyx_v_LOG2PI; double __pyx_v_h00; double __pyx_v_h01; double __pyx_v_h11; double __pyx_v_det; double __pyx_v_e0; double __pyx_v_e1; double __pyx_v_logdet; double __pyx_v_quad; PyObject *__pyx_v_H_np = NULL; PyObject *__pyx_v_eps_np = NULL; PyObject *__pyx_v_Ht = NULL; PyObject *__pyx_v_sign = NULL; PyObject *__pyx_v_ldet = NULL; PyObject *__pyx_v_Hinv = NULL; PyObject *__pyx_v_et = NULL; PyObject *__pyx_r = NULL; /* … */ /* function exit code */ __pyx_L1_error:; __Pyx_XDECREF(__pyx_t_1); __Pyx_XDECREF(__pyx_t_2); __Pyx_XDECREF(__pyx_t_3); __Pyx_XDECREF(__pyx_t_4); __Pyx_XDECREF(__pyx_t_5); __Pyx_XDECREF(__pyx_t_6); __Pyx_XDECREF(__pyx_t_7); __PYX_XCLEAR_MEMVIEW(&__pyx_t_9, 1); __Pyx_AddTraceback("mfe.multivariate._core._bekk_diagonal_recursion", __pyx_clineno, __pyx_lineno, __pyx_filename); __pyx_r = NULL; __pyx_L0:; __Pyx_XDECREF(__pyx_v_H); __PYX_XCLEAR_MEMVIEW(&__pyx_v_Hv, 1); __Pyx_XDECREF(__pyx_v_H_np); __Pyx_XDECREF(__pyx_v_eps_np); __Pyx_XDECREF(__pyx_v_Ht); __Pyx_XDECREF(__pyx_v_sign); __Pyx_XDECREF(__pyx_v_ldet); __Pyx_XDECREF(__pyx_v_Hinv); __Pyx_XDECREF(__pyx_v_et); __Pyx_XGIVEREF(__pyx_r); __Pyx_RefNannyFinishContext(); return __pyx_r; } /* … */ __pyx_t_4 = __Pyx_CyFunction_New(&__pyx_mdef_3mfe_12multivariate_5_core_5_bekk_diagonal_recursion, 0, __pyx_mstate_global->__pyx_n_u_bekk_diagonal_recursion, NULL, __pyx_mstate_global->__pyx_n_u_mfe_multivariate__core, __pyx_mstate_global->__pyx_d, ((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[2])); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 175, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); #if CYTHON_COMPILING_IN_CPYTHON && PY_VERSION_HEX >= 0x030E0000 PyUnstable_Object_EnableDeferredRefcount(__pyx_t_4); #endif if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_bekk_diagonal_recursion, __pyx_t_4) < (0)) __PYX_ERR(0, 175, __pyx_L1_error) __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
176: double[:, ::1] eps, # (T, K)
177: double[:, ::1] CC, # (K, K) precomputed C'C
178: double[:, ::1] AtA, # (K, K) precomputed A_diag^2 as outer product (diag * diag)
179: double[:, ::1] BtB, # (K, K) precomputed B_diag^2 as outer product
180: double[:, ::1] H0, # (K, K)
181: ):
182: """
183: H_t = CC + AtA * (eps_{t-1} eps_{t-1}') + BtB * H_{t-1}
184:
185: (element-wise products — A, B are diagonal so AtA[i,j] = A[i]*A[j])
186: """
+187: cdef int T = eps.shape[0]
__pyx_v_T = (__pyx_v_eps.shape[0]);
+188: cdef int K = eps.shape[1]
__pyx_v_K = (__pyx_v_eps.shape[1]);
189: cdef int t, i, j
190:
+191: H = np.empty((T, K, K), dtype=np.float64)
__pyx_t_2 = NULL; __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 191, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_3); __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_mstate_global->__pyx_n_u_empty); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 191, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0; __pyx_t_3 = __Pyx_PyLong_From_int(__pyx_v_T); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 191, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_3); __pyx_t_5 = __Pyx_PyLong_From_int(__pyx_v_K); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 191, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_5); __pyx_t_6 = __Pyx_PyLong_From_int(__pyx_v_K); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 191, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); __pyx_t_7 = PyTuple_New(3); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 191, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_7); __Pyx_GIVEREF(__pyx_t_3); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_7, 0, __pyx_t_3) != (0)) __PYX_ERR(0, 191, __pyx_L1_error); __Pyx_GIVEREF(__pyx_t_5); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_7, 1, __pyx_t_5) != (0)) __PYX_ERR(0, 191, __pyx_L1_error); __Pyx_GIVEREF(__pyx_t_6); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_7, 2, __pyx_t_6) != (0)) __PYX_ERR(0, 191, __pyx_L1_error); __pyx_t_3 = 0; __pyx_t_5 = 0; __pyx_t_6 = 0; __Pyx_GetModuleGlobalName(__pyx_t_6, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 191, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_mstate_global->__pyx_n_u_float64); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 191, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_5); __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; __pyx_t_8 = 1; #if CYTHON_UNPACK_METHODS if (unlikely(PyMethod_Check(__pyx_t_4))) { __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_4); assert(__pyx_t_2); PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_4); __Pyx_INCREF(__pyx_t_2); __Pyx_INCREF(__pyx__function); __Pyx_DECREF_SET(__pyx_t_4, __pyx__function); __pyx_t_8 = 0; } #endif { PyObject *__pyx_callargs[3] = {__pyx_t_2, __pyx_t_7, __pyx_t_5}; #if CYTHON_VECTORCALL __pyx_t_6 = __pyx_mstate_global->__pyx_tuple[2]; if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 191, __pyx_L1_error) __Pyx_INCREF(__pyx_t_6); #else { PyObject *__pyx_temp[1] = {__pyx_mstate_global->__pyx_n_u_dtype}; __pyx_t_6 = __Pyx_MakeKwargDict(__pyx_temp, __pyx_callargs+2, 1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 191, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); } #endif __pyx_t_1 = __Pyx_Object_VectorcallKwds((PyObject*)__pyx_t_4, __pyx_callargs+__pyx_t_8, (2-__pyx_t_8) | (__pyx_t_8*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET), __pyx_t_6); __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0; __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0; __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0; __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 191, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); } __pyx_v_H = __pyx_t_1; __pyx_t_1 = 0;
+192: cdef double[:, :, ::1] Hv = H
__pyx_t_9 = __Pyx_PyObject_to_MemoryviewSlice_d_d_dc_double(__pyx_v_H, PyBUF_WRITABLE); if (unlikely(!__pyx_t_9.memview)) __PYX_ERR(0, 192, __pyx_L1_error) __pyx_v_Hv = __pyx_t_9; __pyx_t_9.memview = NULL; __pyx_t_9.data = NULL;
193:
+194: for i in range(K):
__pyx_t_10 = __pyx_v_K;
__pyx_t_11 = __pyx_t_10;
for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
__pyx_v_i = __pyx_t_12;
+195: for j in range(K):
__pyx_t_13 = __pyx_v_K;
__pyx_t_14 = __pyx_t_13;
for (__pyx_t_15 = 0; __pyx_t_15 < __pyx_t_14; __pyx_t_15+=1) {
__pyx_v_j = __pyx_t_15;
+196: Hv[0, i, j] = H0[i, j]
__pyx_t_16 = __pyx_v_i;
__pyx_t_17 = __pyx_v_j;
__pyx_t_18 = 0;
__pyx_t_19 = __pyx_v_i;
__pyx_t_20 = __pyx_v_j;
*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Hv.data + __pyx_t_18 * __pyx_v_Hv.strides[0]) ) + __pyx_t_19 * __pyx_v_Hv.strides[1]) )) + __pyx_t_20)) )) = (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_H0.data + __pyx_t_16 * __pyx_v_H0.strides[0]) )) + __pyx_t_17)) )));
}
}
197:
+198: cdef double ll = 0.0
__pyx_v_ll = 0.0;
+199: cdef double LOG2PI = 1.8378770664093453
__pyx_v_LOG2PI = 1.8378770664093453;
200:
+201: for t in range(1, T):
__pyx_t_10 = __pyx_v_T;
__pyx_t_11 = __pyx_t_10;
for (__pyx_t_12 = 1; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
__pyx_v_t = __pyx_t_12;
+202: for i in range(K):
__pyx_t_13 = __pyx_v_K;
__pyx_t_14 = __pyx_t_13;
for (__pyx_t_15 = 0; __pyx_t_15 < __pyx_t_14; __pyx_t_15+=1) {
__pyx_v_i = __pyx_t_15;
+203: for j in range(K):
__pyx_t_21 = __pyx_v_K;
__pyx_t_22 = __pyx_t_21;
for (__pyx_t_23 = 0; __pyx_t_23 < __pyx_t_22; __pyx_t_23+=1) {
__pyx_v_j = __pyx_t_23;
+204: Hv[t, i, j] = (
__pyx_t_32 = __pyx_v_t;
__pyx_t_33 = __pyx_v_i;
__pyx_t_34 = __pyx_v_j;
*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Hv.data + __pyx_t_32 * __pyx_v_Hv.strides[0]) ) + __pyx_t_33 * __pyx_v_Hv.strides[1]) )) + __pyx_t_34)) )) = (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_CC.data + __pyx_t_17 * __pyx_v_CC.strides[0]) )) + __pyx_t_16)) ))) + (((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_AtA.data + __pyx_t_20 * __pyx_v_AtA.strides[0]) )) + __pyx_t_19)) ))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_eps.data + __pyx_t_18 * __pyx_v_eps.strides[0]) )) + __pyx_t_24)) )))) * (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_eps.data + __pyx_t_25 * __pyx_v_eps.strides[0]) )) + __pyx_t_26)) ))))) + ((*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_BtB.data + __pyx_t_27 * __pyx_v_BtB.strides[0]) )) + __pyx_t_28)) ))) * (*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Hv.data + __pyx_t_29 * __pyx_v_Hv.strides[0]) ) + __pyx_t_30 * __pyx_v_Hv.strides[1]) )) + __pyx_t_31)) )))));
}
}
}
+205: CC[i, j]
__pyx_t_17 = __pyx_v_i;
__pyx_t_16 = __pyx_v_j;
+206: + AtA[i, j] * eps[t - 1, i] * eps[t - 1, j]
__pyx_t_20 = __pyx_v_i;
__pyx_t_19 = __pyx_v_j;
__pyx_t_18 = (__pyx_v_t - 1);
__pyx_t_24 = __pyx_v_i;
__pyx_t_25 = (__pyx_v_t - 1);
__pyx_t_26 = __pyx_v_j;
+207: + BtB[i, j] * Hv[t - 1, i, j]
__pyx_t_27 = __pyx_v_i;
__pyx_t_28 = __pyx_v_j;
__pyx_t_29 = (__pyx_v_t - 1);
__pyx_t_30 = __pyx_v_i;
__pyx_t_31 = __pyx_v_j;
208: )
209:
210: # Log-likelihood (2x2 fast path + general)
211: cdef double h00, h01, h11, det, e0, e1, logdet, quad
+212: if K == 2:
__pyx_t_35 = (__pyx_v_K == 2);
if (__pyx_t_35) {
/* … */
goto __pyx_L13;
}
+213: for t in range(T):
__pyx_t_10 = __pyx_v_T;
__pyx_t_11 = __pyx_t_10;
for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
__pyx_v_t = __pyx_t_12;
+214: h00 = Hv[t, 0, 0]; h01 = Hv[t, 0, 1]; h11 = Hv[t, 1, 1]
__pyx_t_31 = __pyx_v_t;
__pyx_t_30 = 0;
__pyx_t_29 = 0;
__pyx_v_h00 = (*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Hv.data + __pyx_t_31 * __pyx_v_Hv.strides[0]) ) + __pyx_t_30 * __pyx_v_Hv.strides[1]) )) + __pyx_t_29)) )));
__pyx_t_29 = __pyx_v_t;
__pyx_t_30 = 0;
__pyx_t_31 = 1;
__pyx_v_h01 = (*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Hv.data + __pyx_t_29 * __pyx_v_Hv.strides[0]) ) + __pyx_t_30 * __pyx_v_Hv.strides[1]) )) + __pyx_t_31)) )));
__pyx_t_31 = __pyx_v_t;
__pyx_t_30 = 1;
__pyx_t_29 = 1;
__pyx_v_h11 = (*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Hv.data + __pyx_t_31 * __pyx_v_Hv.strides[0]) ) + __pyx_t_30 * __pyx_v_Hv.strides[1]) )) + __pyx_t_29)) )));
+215: det = h00 * h11 - h01 * h01
__pyx_v_det = ((__pyx_v_h00 * __pyx_v_h11) - (__pyx_v_h01 * __pyx_v_h01));
+216: if det <= 0.0:
__pyx_t_35 = (__pyx_v_det <= 0.0);
if (__pyx_t_35) {
/* … */
}
+217: return H, 1e10
__pyx_t_1 = PyTuple_New(2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 217, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __Pyx_INCREF(__pyx_v_H); __Pyx_GIVEREF(__pyx_v_H); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_v_H) != (0)) __PYX_ERR(0, 217, __pyx_L1_error); __Pyx_INCREF(__pyx_mstate_global->__pyx_float_1e10); __Pyx_GIVEREF(__pyx_mstate_global->__pyx_float_1e10); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_1, 1, __pyx_mstate_global->__pyx_float_1e10) != (0)) __PYX_ERR(0, 217, __pyx_L1_error); { PyObject *__pyx_temp; { __pyx_temp = __pyx_r; __pyx_r = __pyx_t_1; } __Pyx_XDECREF(__pyx_temp); } __pyx_t_1 = 0; goto __pyx_L0;
+218: e0 = eps[t, 0]; e1 = eps[t, 1]
__pyx_t_29 = __pyx_v_t;
__pyx_t_30 = 0;
__pyx_v_e0 = (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_eps.data + __pyx_t_29 * __pyx_v_eps.strides[0]) )) + __pyx_t_30)) )));
__pyx_t_30 = __pyx_v_t;
__pyx_t_29 = 1;
__pyx_v_e1 = (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_eps.data + __pyx_t_30 * __pyx_v_eps.strides[0]) )) + __pyx_t_29)) )));
+219: logdet = log(det)
__pyx_v_logdet = log(__pyx_v_det);
+220: quad = (h11 * e0 * e0 - 2.0 * h01 * e0 * e1 + h00 * e1 * e1) / det
__pyx_v_quad = (((((__pyx_v_h11 * __pyx_v_e0) * __pyx_v_e0) - (((2.0 * __pyx_v_h01) * __pyx_v_e0) * __pyx_v_e1)) + ((__pyx_v_h00 * __pyx_v_e1) * __pyx_v_e1)) / __pyx_v_det);
+221: ll += logdet + quad
__pyx_v_ll = (__pyx_v_ll + (__pyx_v_logdet + __pyx_v_quad));
}
+222: ll = 0.5 * (T * K * LOG2PI + ll)
__pyx_v_ll = (0.5 * (((__pyx_v_T * __pyx_v_K) * __pyx_v_LOG2PI) + __pyx_v_ll));
223: else:
+224: H_np = np.asarray(H)
/*else*/ {
__pyx_t_4 = NULL;
__Pyx_GetModuleGlobalName(__pyx_t_6, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 224, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_6, __pyx_mstate_global->__pyx_n_u_asarray); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 224, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_5);
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
__pyx_t_8 = 1;
#if CYTHON_UNPACK_METHODS
if (unlikely(PyMethod_Check(__pyx_t_5))) {
__pyx_t_4 = PyMethod_GET_SELF(__pyx_t_5);
assert(__pyx_t_4);
PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_5);
__Pyx_INCREF(__pyx_t_4);
__Pyx_INCREF(__pyx__function);
__Pyx_DECREF_SET(__pyx_t_5, __pyx__function);
__pyx_t_8 = 0;
}
#endif
{
PyObject *__pyx_callargs[2] = {__pyx_t_4, __pyx_v_H};
__pyx_t_1 = __Pyx_PyObject_FastCall((PyObject*)__pyx_t_5, __pyx_callargs+__pyx_t_8, (2-__pyx_t_8) | (__pyx_t_8*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
__Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 224, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
}
__pyx_v_H_np = __pyx_t_1;
__pyx_t_1 = 0;
+225: eps_np = np.asarray(eps)
__pyx_t_5 = NULL;
__Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 225, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_mstate_global->__pyx_n_u_asarray); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 225, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_6);
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__pyx_t_4 = __pyx_memoryview_fromslice(__pyx_v_eps, 2, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 225, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_4);
__pyx_t_8 = 1;
#if CYTHON_UNPACK_METHODS
if (unlikely(PyMethod_Check(__pyx_t_6))) {
__pyx_t_5 = PyMethod_GET_SELF(__pyx_t_6);
assert(__pyx_t_5);
PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_6);
__Pyx_INCREF(__pyx_t_5);
__Pyx_INCREF(__pyx__function);
__Pyx_DECREF_SET(__pyx_t_6, __pyx__function);
__pyx_t_8 = 0;
}
#endif
{
PyObject *__pyx_callargs[2] = {__pyx_t_5, __pyx_t_4};
__pyx_t_1 = __Pyx_PyObject_FastCall((PyObject*)__pyx_t_6, __pyx_callargs+__pyx_t_8, (2-__pyx_t_8) | (__pyx_t_8*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
__Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 225, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_1);
}
__pyx_v_eps_np = __pyx_t_1;
__pyx_t_1 = 0;
+226: for t in range(T):
__pyx_t_10 = __pyx_v_T;
__pyx_t_11 = __pyx_t_10;
for (__pyx_t_12 = 0; __pyx_t_12 < __pyx_t_11; __pyx_t_12+=1) {
__pyx_v_t = __pyx_t_12;
+227: Ht = H_np[t]
__pyx_t_1 = __Pyx_GetItemInt(__pyx_v_H_np, __pyx_v_t, int, 1, __Pyx_PyLong_From_int, 0, 0, 1, __Pyx_ReferenceSharing_OwnStrongReference); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 227, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __Pyx_XDECREF_SET(__pyx_v_Ht, __pyx_t_1); __pyx_t_1 = 0;
+228: sign, ldet = np.linalg.slogdet(Ht)
__Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 228, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_mstate_global->__pyx_n_u_linalg); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 228, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_5); __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; __pyx_t_6 = __pyx_t_5; __Pyx_INCREF(__pyx_t_6); __pyx_t_8 = 0; { PyObject *__pyx_callargs[2] = {__pyx_t_6, __pyx_v_Ht}; __pyx_t_1 = __Pyx_PyObject_FastCallMethod((PyObject*)__pyx_mstate_global->__pyx_n_u_slogdet, __pyx_callargs+__pyx_t_8, (2-__pyx_t_8) | (1*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET)); __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0; __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 228, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); } if ((likely(PyTuple_CheckExact(__pyx_t_1))) || (PyList_CheckExact(__pyx_t_1))) { PyObject* sequence = __pyx_t_1; Py_ssize_t size = __Pyx_PySequence_SIZE(sequence); if (unlikely(size != 2)) { if (size > 2) __Pyx_RaiseTooManyValuesError(2); else if (size >= 0) __Pyx_RaiseNeedMoreValuesError(size); __PYX_ERR(0, 228, __pyx_L1_error) } #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS if (likely(PyTuple_CheckExact(sequence))) { __pyx_t_5 = PyTuple_GET_ITEM(sequence, 0); __Pyx_INCREF(__pyx_t_5); __pyx_t_6 = PyTuple_GET_ITEM(sequence, 1); __Pyx_INCREF(__pyx_t_6); } else { __pyx_t_5 = __Pyx_PyList_GET_ITEM_REF(sequence, 0, __Pyx_ReferenceSharing_SharedReference); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 228, __pyx_L1_error) __Pyx_XGOTREF(__pyx_t_5); __pyx_t_6 = __Pyx_PyList_GET_ITEM_REF(sequence, 1, __Pyx_ReferenceSharing_SharedReference); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 228, __pyx_L1_error) __Pyx_XGOTREF(__pyx_t_6); } #else __pyx_t_5 = __Pyx_PySequence_ITEM(sequence, 0); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 228, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_5); __pyx_t_6 = __Pyx_PySequence_ITEM(sequence, 1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 228, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); #endif __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; } else { Py_ssize_t index = -1; __pyx_t_4 = PyObject_GetIter(__pyx_t_1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 228, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __pyx_t_36 = (CYTHON_COMPILING_IN_LIMITED_API) ? PyIter_Next : __Pyx_PyObject_GetIterNextFunc(__pyx_t_4); index = 0; __pyx_t_5 = __pyx_t_36(__pyx_t_4); if (unlikely(!__pyx_t_5)) goto __pyx_L19_unpacking_failed; __Pyx_GOTREF(__pyx_t_5); index = 1; __pyx_t_6 = __pyx_t_36(__pyx_t_4); if (unlikely(!__pyx_t_6)) goto __pyx_L19_unpacking_failed; __Pyx_GOTREF(__pyx_t_6); if (__Pyx_IternextUnpackEndCheck(__pyx_t_36(__pyx_t_4), 2) < (0)) __PYX_ERR(0, 228, __pyx_L1_error) __pyx_t_36 = NULL; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; goto __pyx_L20_unpacking_done; __pyx_L19_unpacking_failed:; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; __pyx_t_36 = NULL; if (__Pyx_IterFinish() == 0) __Pyx_RaiseNeedMoreValuesError(index); __PYX_ERR(0, 228, __pyx_L1_error) __pyx_L20_unpacking_done:; } __Pyx_XDECREF_SET(__pyx_v_sign, __pyx_t_5); __pyx_t_5 = 0; __Pyx_XDECREF_SET(__pyx_v_ldet, __pyx_t_6); __pyx_t_6 = 0;
+229: if sign <= 0:
__pyx_t_35 = __Pyx_PyObject_CompareBoolLe_object_int(__pyx_v_sign, __pyx_mstate_global->__pyx_int_0, Py_LE); if (unlikely((__pyx_t_35 < 0))) __PYX_ERR(0, 229, __pyx_L1_error) if (__pyx_t_35) { /* … */ }
+230: return H, 1e10
__pyx_t_1 = PyTuple_New(2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 230, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_1); __Pyx_INCREF(__pyx_v_H); __Pyx_GIVEREF(__pyx_v_H); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_v_H) != (0)) __PYX_ERR(0, 230, __pyx_L1_error); __Pyx_INCREF(__pyx_mstate_global->__pyx_float_1e10); __Pyx_GIVEREF(__pyx_mstate_global->__pyx_float_1e10); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_1, 1, __pyx_mstate_global->__pyx_float_1e10) != (0)) __PYX_ERR(0, 230, __pyx_L1_error); { PyObject *__pyx_temp; { __pyx_temp = __pyx_r; __pyx_r = __pyx_t_1; } __Pyx_XDECREF(__pyx_temp); } __pyx_t_1 = 0; goto __pyx_L0;
+231: try:
{
/*try:*/ {
/* … */
}
__Pyx_XDECREF(__pyx_t_37); __pyx_t_37 = 0;
__Pyx_XDECREF(__pyx_t_38); __pyx_t_38 = 0;
__Pyx_XDECREF(__pyx_t_39); __pyx_t_39 = 0;
goto __pyx_L29_try_end;
__pyx_L22_error:;
__Pyx_XDECREF(__pyx_t_1); __pyx_t_1 = 0;
__Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
__Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
__Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
__Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
__Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0;
__Pyx_XDECREF(__pyx_t_7); __pyx_t_7 = 0;
__PYX_XCLEAR_MEMVIEW(&__pyx_t_9, 1);; __pyx_t_9.memview = NULL; __pyx_t_9.data = NULL;
/* … */
__pyx_L24_except_error:;
__Pyx_XGIVEREF(__pyx_t_37);
__Pyx_XGIVEREF(__pyx_t_38);
__Pyx_XGIVEREF(__pyx_t_39);
__Pyx_ExceptionReset(__pyx_t_37, __pyx_t_38, __pyx_t_39);
goto __pyx_L1_error;
__pyx_L25_except_return:;
__Pyx_XGIVEREF(__pyx_t_37);
__Pyx_XGIVEREF(__pyx_t_38);
__Pyx_XGIVEREF(__pyx_t_39);
__Pyx_ExceptionReset(__pyx_t_37, __pyx_t_38, __pyx_t_39);
goto __pyx_L0;
__pyx_L29_try_end:;
}
+232: Hinv = np.linalg.inv(Ht)
__Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 232, __pyx_L22_error) __Pyx_GOTREF(__pyx_t_5); __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_mstate_global->__pyx_n_u_linalg); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 232, __pyx_L22_error) __Pyx_GOTREF(__pyx_t_4); __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0; __pyx_t_6 = __pyx_t_4; __Pyx_INCREF(__pyx_t_6); __pyx_t_8 = 0; { PyObject *__pyx_callargs[2] = {__pyx_t_6, __pyx_v_Ht}; __pyx_t_1 = __Pyx_PyObject_FastCallMethod((PyObject*)__pyx_mstate_global->__pyx_n_u_inv, __pyx_callargs+__pyx_t_8, (2-__pyx_t_8) | (1*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET)); __Pyx_XDECREF(__pyx_t_6); __pyx_t_6 = 0; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 232, __pyx_L22_error) __Pyx_GOTREF(__pyx_t_1); } __Pyx_XDECREF_SET(__pyx_v_Hinv, __pyx_t_1); __pyx_t_1 = 0;
+233: except Exception:
__pyx_t_13 = __Pyx_PyErr_ExceptionMatches(((PyObject *)(((PyTypeObject*)PyExc_Exception)))); if (__pyx_t_13) { __Pyx_AddTraceback("mfe.multivariate._core._bekk_diagonal_recursion", __pyx_clineno, __pyx_lineno, __pyx_filename); if (__Pyx_GetException(&__pyx_t_1, &__pyx_t_4, &__pyx_t_6) < 0) __PYX_ERR(0, 233, __pyx_L24_except_error) __Pyx_XGOTREF(__pyx_t_1); __Pyx_XGOTREF(__pyx_t_4); __Pyx_XGOTREF(__pyx_t_6);
+234: return H, 1e10
__pyx_t_5 = PyTuple_New(2); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 234, __pyx_L24_except_error) __Pyx_GOTREF(__pyx_t_5); __Pyx_INCREF(__pyx_v_H); __Pyx_GIVEREF(__pyx_v_H); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_v_H) != (0)) __PYX_ERR(0, 234, __pyx_L24_except_error); __Pyx_INCREF(__pyx_mstate_global->__pyx_float_1e10); __Pyx_GIVEREF(__pyx_mstate_global->__pyx_float_1e10); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_5, 1, __pyx_mstate_global->__pyx_float_1e10) != (0)) __PYX_ERR(0, 234, __pyx_L24_except_error); { PyObject *__pyx_temp; { __pyx_temp = __pyx_r; __pyx_r = __pyx_t_5; } __Pyx_XDECREF(__pyx_temp); } __pyx_t_5 = 0; __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; goto __pyx_L25_except_return; } goto __pyx_L24_except_error;
+235: et = eps_np[t]
__pyx_t_6 = __Pyx_GetItemInt(__pyx_v_eps_np, __pyx_v_t, int, 1, __Pyx_PyLong_From_int, 0, 0, 1, __Pyx_ReferenceSharing_OwnStrongReference); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 235, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); __Pyx_XDECREF_SET(__pyx_v_et, __pyx_t_6); __pyx_t_6 = 0;
+236: ll += float(ldet) + float(et @ Hinv @ et)
__pyx_t_40 = __Pyx_PyObject_AsDouble(__pyx_v_ldet); if (unlikely(__PYX_CHECK_FLOAT_EXCEPTION(__pyx_t_40, ((double)((double)-1))) && PyErr_Occurred())) __PYX_ERR(0, 236, __pyx_L1_error) __pyx_t_6 = __Pyx_PyNumber_MatrixMultiply(__pyx_v_et, __pyx_v_Hinv); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 236, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); __pyx_t_4 = __Pyx_PyNumber_MatrixMultiply(__pyx_t_6, __pyx_v_et); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 236, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; __pyx_t_41 = __Pyx_PyObject_AsDouble(__pyx_t_4); if (unlikely(__PYX_CHECK_FLOAT_EXCEPTION(__pyx_t_41, ((double)((double)-1))) && PyErr_Occurred())) __PYX_ERR(0, 236, __pyx_L1_error) __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; __pyx_v_ll = (__pyx_v_ll + (__pyx_t_40 + __pyx_t_41)); }
+237: ll = 0.5 * (T * K * LOG2PI + ll)
__pyx_v_ll = (0.5 * (((__pyx_v_T * __pyx_v_K) * __pyx_v_LOG2PI) + __pyx_v_ll)); } __pyx_L13:;
238:
+239: return H, ll
__pyx_t_4 = PyFloat_FromDouble(__pyx_v_ll); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 239, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); __pyx_t_6 = PyTuple_New(2); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 239, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_6); __Pyx_INCREF(__pyx_v_H); __Pyx_GIVEREF(__pyx_v_H); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_6, 0, __pyx_v_H) != (0)) __PYX_ERR(0, 239, __pyx_L1_error); __Pyx_GIVEREF(__pyx_t_4); if (__Pyx_PyTuple_SET_ITEM(__pyx_t_6, 1, __pyx_t_4) != (0)) __PYX_ERR(0, 239, __pyx_L1_error); __pyx_t_4 = 0; { PyObject *__pyx_temp; { __pyx_temp = __pyx_r; __pyx_r = __pyx_t_6; } __Pyx_XDECREF(__pyx_temp); } __pyx_t_6 = 0; goto __pyx_L0;
240:
241:
242: # ────────────────────────────────────────────────────────────────────────────
243: # 4. DCC correlation log-likelihood (step 2 objective)
244: # ────────────────────────────────────────────────────────────────────────────
245:
+246: def _dcc_corr_loglik(
/* Python wrapper */ static PyObject *__pyx_pw_3mfe_12multivariate_5_core_7_dcc_corr_loglik(PyObject *__pyx_self, #if CYTHON_VECTORCALL PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds #else PyObject *__pyx_args, PyObject *__pyx_kwds #endif ); /*proto*/ PyDoc_STRVAR(__pyx_doc_3mfe_12multivariate_5_core_6_dcc_corr_loglik, "\n DCC correlation log-likelihood (step 2, Engle 2002 eq. 12):\n\n L2 = 0.5 * sum_t [log|R_t| + z_t\047 R_t^{-1} z_t - z_t\047 z_t]\n\n where R_t = diag(Q_t)^{-1/2} Q_t diag(Q_t)^{-1/2}\n\n 2x2 fast path avoids per-step numpy inversion.\n "); static PyMethodDef __pyx_mdef_3mfe_12multivariate_5_core_7_dcc_corr_loglik = {"_dcc_corr_loglik", (PyCFunction)(void(*)(void))(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_3mfe_12multivariate_5_core_7_dcc_corr_loglik, __Pyx_METH_FASTCALL|METH_KEYWORDS, __pyx_doc_3mfe_12multivariate_5_core_6_dcc_corr_loglik}; static PyObject *__pyx_pw_3mfe_12multivariate_5_core_7_dcc_corr_loglik(PyObject *__pyx_self, #if CYTHON_VECTORCALL PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds #else PyObject *__pyx_args, PyObject *__pyx_kwds #endif ) { __Pyx_memviewslice __pyx_v_z = { 0, 0, { 0 }, { 0 }, { 0 } }; __Pyx_memviewslice __pyx_v_Q = { 0, 0, { 0 }, { 0 }, { 0 } }; #if !CYTHON_VECTORCALL CYTHON_UNUSED Py_ssize_t __pyx_nargs; #endif CYTHON_UNUSED PyObject *const *__pyx_kwvalues; PyObject *__pyx_r = 0; __Pyx_RefNannyDeclarations __Pyx_RefNannySetupContext("_dcc_corr_loglik (wrapper)", 0); #if !CYTHON_VECTORCALL #if CYTHON_ASSUME_SAFE_SIZE __pyx_nargs = PyTuple_GET_SIZE(__pyx_args); #else __pyx_nargs = PyTuple_Size(__pyx_args); if (unlikely(__pyx_nargs < 0)) return NULL; #endif #endif __pyx_kwvalues = __Pyx_KwValues_FASTCALL(__pyx_args, __pyx_nargs); { PyObject ** const __pyx_pyargnames[] = {&__pyx_mstate_global->__pyx_n_u_z,&__pyx_mstate_global->__pyx_n_u_Q,0}; PyObject* values[2] = {0,0}; const Py_ssize_t __pyx_kwds_len = (__pyx_kwds) ? __Pyx_NumKwargs_FASTCALL(__pyx_kwds) : 0; if (unlikely(__pyx_kwds_len < 0)) __PYX_ERR(0, 246, __pyx_L3_error) if (__pyx_kwds_len > 0) { switch (__pyx_nargs) { case 2: values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 246, __pyx_L3_error) CYTHON_FALLTHROUGH; case 1: values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 246, __pyx_L3_error) CYTHON_FALLTHROUGH; case 0: break; default: goto __pyx_L5_argtuple_error; } const Py_ssize_t kwd_pos_args = __pyx_nargs; if (__Pyx_ParseKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values, kwd_pos_args, __pyx_kwds_len, "_dcc_corr_loglik", 0) < (0)) __PYX_ERR(0, 246, __pyx_L3_error) for (Py_ssize_t i = __pyx_nargs; i < 2; i++) { if (unlikely(!values[i])) { __Pyx_RaiseArgtupleInvalid("_dcc_corr_loglik", 1, 2, 2, i); __PYX_ERR(0, 246, __pyx_L3_error) } } } else if (unlikely(__pyx_nargs != 2)) { goto __pyx_L5_argtuple_error; } else { values[0] = __Pyx_ArgRef_FASTCALL(__pyx_args, 0); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[0])) __PYX_ERR(0, 246, __pyx_L3_error) values[1] = __Pyx_ArgRef_FASTCALL(__pyx_args, 1); if (!CYTHON_ASSUME_SAFE_MACROS && unlikely(!values[1])) __PYX_ERR(0, 246, __pyx_L3_error) } __pyx_v_z = __Pyx_PyObject_to_MemoryviewSlice_d_dc_double(values[0], PyBUF_WRITABLE); if (unlikely(!__pyx_v_z.memview)) __PYX_ERR(0, 247, __pyx_L3_error) __pyx_v_Q = __Pyx_PyObject_to_MemoryviewSlice_d_d_dc_double(values[1], PyBUF_WRITABLE); if (unlikely(!__pyx_v_Q.memview)) __PYX_ERR(0, 248, __pyx_L3_error) } goto __pyx_L6_skip; __pyx_L5_argtuple_error:; __Pyx_RaiseArgtupleInvalid("_dcc_corr_loglik", 1, 2, 2, __pyx_nargs); __PYX_ERR(0, 246, __pyx_L3_error) __pyx_L6_skip:; goto __pyx_L4_argument_unpacking_done; __pyx_L3_error:; for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) { Py_XDECREF(values[__pyx_temp]); } __PYX_XCLEAR_MEMVIEW(&__pyx_v_z, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_Q, 1); __Pyx_AddTraceback("mfe.multivariate._core._dcc_corr_loglik", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); return NULL; __pyx_L4_argument_unpacking_done:; __pyx_r = __pyx_pf_3mfe_12multivariate_5_core_6_dcc_corr_loglik(__pyx_self, __pyx_v_z, __pyx_v_Q); int __pyx_lineno = 0; const char *__pyx_filename = NULL; int __pyx_clineno = 0; /* function exit code */ for (Py_ssize_t __pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) { Py_XDECREF(values[__pyx_temp]); } __PYX_XCLEAR_MEMVIEW(&__pyx_v_z, 1); __PYX_XCLEAR_MEMVIEW(&__pyx_v_Q, 1); __Pyx_RefNannyFinishContext(); return __pyx_r; } static PyObject *__pyx_pf_3mfe_12multivariate_5_core_6_dcc_corr_loglik(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_z, __Pyx_memviewslice __pyx_v_Q) { int __pyx_v_T; int __pyx_v_K; int __pyx_v_t; double __pyx_v_ll; double __pyx_v_q00; double __pyx_v_q01; double __pyx_v_q11; double __pyx_v_r01; double __pyx_v_det_r; double __pyx_v_z0; double __pyx_v_z1; double __pyx_v_quad_r; double __pyx_v_quad_z; PyObject *__pyx_v_Q_np = NULL; PyObject *__pyx_v_z_np = NULL; double __pyx_v_ll_py; PyObject *__pyx_v_Qt = NULL; PyObject *__pyx_v_d_inv = NULL; PyObject *__pyx_v_Rt = NULL; PyObject *__pyx_v_sign = NULL; PyObject *__pyx_v_ldet = NULL; PyObject *__pyx_v_zt = NULL; PyObject *__pyx_v_Rinv = NULL; PyObject *__pyx_r = NULL; /* … */ __pyx_t_4 = __Pyx_CyFunction_New(&__pyx_mdef_3mfe_12multivariate_5_core_7_dcc_corr_loglik, 0, __pyx_mstate_global->__pyx_n_u_dcc_corr_loglik, NULL, __pyx_mstate_global->__pyx_n_u_mfe_multivariate__core, __pyx_mstate_global->__pyx_d, ((PyObject *)__pyx_mstate_global->__pyx_codeobj_tab[3])); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 246, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_4); #if CYTHON_COMPILING_IN_CPYTHON && PY_VERSION_HEX >= 0x030E0000 PyUnstable_Object_EnableDeferredRefcount(__pyx_t_4); #endif if (PyDict_SetItem(__pyx_mstate_global->__pyx_d, __pyx_mstate_global->__pyx_n_u_dcc_corr_loglik, __pyx_t_4) < (0)) __PYX_ERR(0, 246, __pyx_L1_error) __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
247: double[:, ::1] z, # (T, K) standardized residuals
248: double[:, :, ::1] Q, # (T, K, K) Q matrices from _dcc_q_recursion
249: ):
250: """
251: DCC correlation log-likelihood (step 2, Engle 2002 eq. 12):
252:
253: L2 = 0.5 * sum_t [log|R_t| + z_t' R_t^{-1} z_t - z_t' z_t]
254:
255: where R_t = diag(Q_t)^{-1/2} Q_t diag(Q_t)^{-1/2}
256:
257: 2x2 fast path avoids per-step numpy inversion.
258: """
+259: cdef int T = z.shape[0]
__pyx_v_T = (__pyx_v_z.shape[0]);
+260: cdef int K = z.shape[1]
__pyx_v_K = (__pyx_v_z.shape[1]);
261: cdef int t, i
262:
+263: cdef double ll = 0.0
__pyx_v_ll = 0.0;
264: cdef double q00, q01, q11, r00, r01, r11, det_r
265: cdef double z0, z1, quad_r, quad_z
266:
+267: if K == 2:
__pyx_t_1 = (__pyx_v_K == 2);
if (__pyx_t_1) {
/* … */
}
+268: for t in range(T):
__pyx_t_2 = __pyx_v_T;
__pyx_t_3 = __pyx_t_2;
for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) {
__pyx_v_t = __pyx_t_4;
+269: q00 = Q[t, 0, 0]; q01 = Q[t, 0, 1]; q11 = Q[t, 1, 1]
__pyx_t_5 = __pyx_v_t;
__pyx_t_6 = 0;
__pyx_t_7 = 0;
__pyx_v_q00 = (*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Q.data + __pyx_t_5 * __pyx_v_Q.strides[0]) ) + __pyx_t_6 * __pyx_v_Q.strides[1]) )) + __pyx_t_7)) )));
__pyx_t_7 = __pyx_v_t;
__pyx_t_6 = 0;
__pyx_t_5 = 1;
__pyx_v_q01 = (*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Q.data + __pyx_t_7 * __pyx_v_Q.strides[0]) ) + __pyx_t_6 * __pyx_v_Q.strides[1]) )) + __pyx_t_5)) )));
__pyx_t_5 = __pyx_v_t;
__pyx_t_6 = 1;
__pyx_t_7 = 1;
__pyx_v_q11 = (*((double *) ( /* dim=2 */ ((char *) (((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_Q.data + __pyx_t_5 * __pyx_v_Q.strides[0]) ) + __pyx_t_6 * __pyx_v_Q.strides[1]) )) + __pyx_t_7)) )));
270: # R_t = diag(Q)^{-1/2} Q diag(Q)^{-1/2}
271: # r[i,j] = q[i,j] / sqrt(q[i,i] * q[j,j])
+272: if q00 <= 0.0 or q11 <= 0.0:
__pyx_t_8 = (__pyx_v_q00 <= 0.0);
if (!__pyx_t_8) {
} else {
__pyx_t_1 = __pyx_t_8;
goto __pyx_L7_bool_binop_done;
}
__pyx_t_8 = (__pyx_v_q11 <= 0.0);
__pyx_t_1 = __pyx_t_8;
__pyx_L7_bool_binop_done:;
if (__pyx_t_1) {
/* … */
}
+273: return 1e10
{
PyObject *__pyx_temp;
{
__pyx_temp = __pyx_r;
__Pyx_INCREF(__pyx_mstate_global->__pyx_float_1e10);
__pyx_r = __pyx_mstate_global->__pyx_float_1e10;
}
__Pyx_XDECREF(__pyx_temp);
}
goto __pyx_L0;
+274: r01 = q01 / sqrt(q00 * q11)
__pyx_v_r01 = (__pyx_v_q01 / sqrt((__pyx_v_q00 * __pyx_v_q11)));
275: # R is [[1, r01],[r01, 1]], det = 1 - r01^2
+276: det_r = 1.0 - r01 * r01
__pyx_v_det_r = (1.0 - (__pyx_v_r01 * __pyx_v_r01));
+277: if det_r <= 1e-14:
__pyx_t_1 = (__pyx_v_det_r <= 1e-14);
if (__pyx_t_1) {
/* … */
}
+278: return 1e10
{
PyObject *__pyx_temp;
{
__pyx_temp = __pyx_r;
__Pyx_INCREF(__pyx_mstate_global->__pyx_float_1e10);
__pyx_r = __pyx_mstate_global->__pyx_float_1e10;
}
__Pyx_XDECREF(__pyx_temp);
}
goto __pyx_L0;
+279: z0 = z[t, 0]; z1 = z[t, 1]
__pyx_t_7 = __pyx_v_t;
__pyx_t_6 = 0;
__pyx_v_z0 = (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_z.data + __pyx_t_7 * __pyx_v_z.strides[0]) )) + __pyx_t_6)) )));
__pyx_t_6 = __pyx_v_t;
__pyx_t_7 = 1;
__pyx_v_z1 = (*((double *) ( /* dim=1 */ ((char *) (((double *) ( /* dim=0 */ (__pyx_v_z.data + __pyx_t_6 * __pyx_v_z.strides[0]) )) + __pyx_t_7)) )));
280: # z' R^{-1} z = (z0^2 - 2*r01*z0*z1 + z1^2) / det_r
+281: quad_r = (z0 * z0 - 2.0 * r01 * z0 * z1 + z1 * z1) / det_r
__pyx_v_quad_r = ((((__pyx_v_z0 * __pyx_v_z0) - (((2.0 * __pyx_v_r01) * __pyx_v_z0) * __pyx_v_z1)) + (__pyx_v_z1 * __pyx_v_z1)) / __pyx_v_det_r);
+282: quad_z = z0 * z0 + z1 * z1
__pyx_v_quad_z = ((__pyx_v_z0 * __pyx_v_z0) + (__pyx_v_z1 * __pyx_v_z1));
+283: ll += log(det_r) + quad_r - quad_z
__pyx_v_ll = (__pyx_v_ll + ((log(__pyx_v_det_r) + __pyx_v_quad_r) - __pyx_v_quad_z));
}
+284: return 0.5 * ll
__pyx_t_9 = PyFloat_FromDouble((0.5 * __pyx_v_ll)); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 284, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_9); { PyObject *__pyx_temp; { __pyx_temp = __pyx_r; __pyx_r = __pyx_t_9; } __Pyx_XDECREF(__pyx_temp); } __pyx_t_9 = 0; goto __pyx_L0;
285: else:
286: # General K: numpy fallback
+287: Q_np = np.asarray(Q)
/*else*/ {
__pyx_t_10 = NULL;
__Pyx_GetModuleGlobalName(__pyx_t_11, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 287, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_11);
__pyx_t_12 = __Pyx_PyObject_GetAttrStr(__pyx_t_11, __pyx_mstate_global->__pyx_n_u_asarray); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 287, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_12);
__Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
__pyx_t_11 = __pyx_memoryview_fromslice(__pyx_v_Q, 3, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 287, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_11);
__pyx_t_13 = 1;
#if CYTHON_UNPACK_METHODS
if (unlikely(PyMethod_Check(__pyx_t_12))) {
__pyx_t_10 = PyMethod_GET_SELF(__pyx_t_12);
assert(__pyx_t_10);
PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_12);
__Pyx_INCREF(__pyx_t_10);
__Pyx_INCREF(__pyx__function);
__Pyx_DECREF_SET(__pyx_t_12, __pyx__function);
__pyx_t_13 = 0;
}
#endif
{
PyObject *__pyx_callargs[2] = {__pyx_t_10, __pyx_t_11};
__pyx_t_9 = __Pyx_PyObject_FastCall((PyObject*)__pyx_t_12, __pyx_callargs+__pyx_t_13, (2-__pyx_t_13) | (__pyx_t_13*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
__Pyx_XDECREF(__pyx_t_10); __pyx_t_10 = 0;
__Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
__Pyx_DECREF(__pyx_t_12); __pyx_t_12 = 0;
if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 287, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_9);
}
__pyx_v_Q_np = __pyx_t_9;
__pyx_t_9 = 0;
+288: z_np = np.asarray(z)
__pyx_t_12 = NULL;
__Pyx_GetModuleGlobalName(__pyx_t_11, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 288, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_11);
__pyx_t_10 = __Pyx_PyObject_GetAttrStr(__pyx_t_11, __pyx_mstate_global->__pyx_n_u_asarray); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 288, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_10);
__Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
__pyx_t_11 = __pyx_memoryview_fromslice(__pyx_v_z, 2, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 288, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_11);
__pyx_t_13 = 1;
#if CYTHON_UNPACK_METHODS
if (unlikely(PyMethod_Check(__pyx_t_10))) {
__pyx_t_12 = PyMethod_GET_SELF(__pyx_t_10);
assert(__pyx_t_12);
PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_10);
__Pyx_INCREF(__pyx_t_12);
__Pyx_INCREF(__pyx__function);
__Pyx_DECREF_SET(__pyx_t_10, __pyx__function);
__pyx_t_13 = 0;
}
#endif
{
PyObject *__pyx_callargs[2] = {__pyx_t_12, __pyx_t_11};
__pyx_t_9 = __Pyx_PyObject_FastCall((PyObject*)__pyx_t_10, __pyx_callargs+__pyx_t_13, (2-__pyx_t_13) | (__pyx_t_13*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
__Pyx_XDECREF(__pyx_t_12); __pyx_t_12 = 0;
__Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
__Pyx_DECREF(__pyx_t_10); __pyx_t_10 = 0;
if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 288, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_9);
}
__pyx_v_z_np = __pyx_t_9;
__pyx_t_9 = 0;
+289: ll_py = 0.0
__pyx_v_ll_py = 0.0;
+290: for t in range(T):
__pyx_t_2 = __pyx_v_T;
__pyx_t_3 = __pyx_t_2;
for (__pyx_t_4 = 0; __pyx_t_4 < __pyx_t_3; __pyx_t_4+=1) {
__pyx_v_t = __pyx_t_4;
+291: Qt = Q_np[t]
__pyx_t_9 = __Pyx_GetItemInt(__pyx_v_Q_np, __pyx_v_t, int, 1, __Pyx_PyLong_From_int, 0, 0, 1, __Pyx_ReferenceSharing_OwnStrongReference); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 291, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_9); __Pyx_XDECREF_SET(__pyx_v_Qt, __pyx_t_9); __pyx_t_9 = 0;
+292: d_inv = 1.0 / np.sqrt(np.diag(Qt))
__pyx_t_10 = NULL;
__Pyx_GetModuleGlobalName(__pyx_t_11, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 292, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_11);
__pyx_t_12 = __Pyx_PyObject_GetAttrStr(__pyx_t_11, __pyx_mstate_global->__pyx_n_u_sqrt); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 292, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_12);
__Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
__pyx_t_14 = NULL;
__Pyx_GetModuleGlobalName(__pyx_t_15, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_15)) __PYX_ERR(0, 292, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_15);
__pyx_t_16 = __Pyx_PyObject_GetAttrStr(__pyx_t_15, __pyx_mstate_global->__pyx_n_u_diag); if (unlikely(!__pyx_t_16)) __PYX_ERR(0, 292, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_16);
__Pyx_DECREF(__pyx_t_15); __pyx_t_15 = 0;
__pyx_t_13 = 1;
#if CYTHON_UNPACK_METHODS
if (unlikely(PyMethod_Check(__pyx_t_16))) {
__pyx_t_14 = PyMethod_GET_SELF(__pyx_t_16);
assert(__pyx_t_14);
PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_16);
__Pyx_INCREF(__pyx_t_14);
__Pyx_INCREF(__pyx__function);
__Pyx_DECREF_SET(__pyx_t_16, __pyx__function);
__pyx_t_13 = 0;
}
#endif
{
PyObject *__pyx_callargs[2] = {__pyx_t_14, __pyx_v_Qt};
__pyx_t_11 = __Pyx_PyObject_FastCall((PyObject*)__pyx_t_16, __pyx_callargs+__pyx_t_13, (2-__pyx_t_13) | (__pyx_t_13*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
__Pyx_XDECREF(__pyx_t_14); __pyx_t_14 = 0;
__Pyx_DECREF(__pyx_t_16); __pyx_t_16 = 0;
if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 292, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_11);
}
__pyx_t_13 = 1;
#if CYTHON_UNPACK_METHODS
if (unlikely(PyMethod_Check(__pyx_t_12))) {
__pyx_t_10 = PyMethod_GET_SELF(__pyx_t_12);
assert(__pyx_t_10);
PyObject* __pyx__function = PyMethod_GET_FUNCTION(__pyx_t_12);
__Pyx_INCREF(__pyx_t_10);
__Pyx_INCREF(__pyx__function);
__Pyx_DECREF_SET(__pyx_t_12, __pyx__function);
__pyx_t_13 = 0;
}
#endif
{
PyObject *__pyx_callargs[2] = {__pyx_t_10, __pyx_t_11};
__pyx_t_9 = __Pyx_PyObject_FastCall((PyObject*)__pyx_t_12, __pyx_callargs+__pyx_t_13, (2-__pyx_t_13) | (__pyx_t_13*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET));
__Pyx_XDECREF(__pyx_t_10); __pyx_t_10 = 0;
__Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
__Pyx_DECREF(__pyx_t_12); __pyx_t_12 = 0;
if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 292, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_9);
}
__pyx_t_12 = __Pyx_PyFloat_TrueDivideCObj(__pyx_mstate_global->__pyx_float_1_0, __pyx_t_9, 1.0, 0, 1); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 292, __pyx_L1_error)
__Pyx_GOTREF(__pyx_t_12);
__Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0;
__Pyx_XDECREF_SET(__pyx_v_d_inv, __pyx_t_12);
__pyx_t_12 = 0;
+293: Rt = d_inv[:, None] * Qt * d_inv[None, :]
__pyx_t_12 = __Pyx_PyObject_GetItem(__pyx_v_d_inv, __pyx_mstate_global->__pyx_tuple[3]); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 293, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_12); __pyx_t_9 = __Pyx_PyNumber_Multiply_object_object(__pyx_t_12, __pyx_v_Qt); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 293, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_9); __Pyx_DECREF(__pyx_t_12); __pyx_t_12 = 0; /* … */ { PyObject* __pyx_temp[2] = {__pyx_mstate_global->__pyx_slice[0], Py_None}; __pyx_mstate_global->__pyx_tuple[3] = __Pyx_PyTuple_FromArray(__pyx_temp, 2); if (unlikely(!__pyx_mstate_global->__pyx_tuple[3])) __PYX_ERR(0, 293, __pyx_L1_error) __Pyx_GOTREF(__pyx_mstate_global->__pyx_tuple[3]); } __Pyx_GIVEREF(__pyx_mstate_global->__pyx_tuple[3]); __pyx_t_12 = __Pyx_PyObject_GetItem(__pyx_v_d_inv, __pyx_mstate_global->__pyx_tuple[4]); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 293, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_12); __pyx_t_11 = __Pyx_PyNumber_Multiply_object_object(__pyx_t_9, __pyx_t_12); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 293, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_11); __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0; __Pyx_DECREF(__pyx_t_12); __pyx_t_12 = 0; __Pyx_XDECREF_SET(__pyx_v_Rt, __pyx_t_11); __pyx_t_11 = 0;
+294: sign, ldet = np.linalg.slogdet(Rt)
__Pyx_GetModuleGlobalName(__pyx_t_9, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 294, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_9); __pyx_t_10 = __Pyx_PyObject_GetAttrStr(__pyx_t_9, __pyx_mstate_global->__pyx_n_u_linalg); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 294, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0; __pyx_t_12 = __pyx_t_10; __Pyx_INCREF(__pyx_t_12); __pyx_t_13 = 0; { PyObject *__pyx_callargs[2] = {__pyx_t_12, __pyx_v_Rt}; __pyx_t_11 = __Pyx_PyObject_FastCallMethod((PyObject*)__pyx_mstate_global->__pyx_n_u_slogdet, __pyx_callargs+__pyx_t_13, (2-__pyx_t_13) | (1*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET)); __Pyx_XDECREF(__pyx_t_12); __pyx_t_12 = 0; __Pyx_DECREF(__pyx_t_10); __pyx_t_10 = 0; if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 294, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_11); } if ((likely(PyTuple_CheckExact(__pyx_t_11))) || (PyList_CheckExact(__pyx_t_11))) { PyObject* sequence = __pyx_t_11; Py_ssize_t size = __Pyx_PySequence_SIZE(sequence); if (unlikely(size != 2)) { if (size > 2) __Pyx_RaiseTooManyValuesError(2); else if (size >= 0) __Pyx_RaiseNeedMoreValuesError(size); __PYX_ERR(0, 294, __pyx_L1_error) } #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS if (likely(PyTuple_CheckExact(sequence))) { __pyx_t_10 = PyTuple_GET_ITEM(sequence, 0); __Pyx_INCREF(__pyx_t_10); __pyx_t_12 = PyTuple_GET_ITEM(sequence, 1); __Pyx_INCREF(__pyx_t_12); } else { __pyx_t_10 = __Pyx_PyList_GET_ITEM_REF(sequence, 0, __Pyx_ReferenceSharing_SharedReference); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 294, __pyx_L1_error) __Pyx_XGOTREF(__pyx_t_10); __pyx_t_12 = __Pyx_PyList_GET_ITEM_REF(sequence, 1, __Pyx_ReferenceSharing_SharedReference); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 294, __pyx_L1_error) __Pyx_XGOTREF(__pyx_t_12); } #else __pyx_t_10 = __Pyx_PySequence_ITEM(sequence, 0); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 294, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_10); __pyx_t_12 = __Pyx_PySequence_ITEM(sequence, 1); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 294, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_12); #endif __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0; } else { Py_ssize_t index = -1; __pyx_t_9 = PyObject_GetIter(__pyx_t_11); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 294, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_9); __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0; __pyx_t_17 = (CYTHON_COMPILING_IN_LIMITED_API) ? PyIter_Next : __Pyx_PyObject_GetIterNextFunc(__pyx_t_9); index = 0; __pyx_t_10 = __pyx_t_17(__pyx_t_9); if (unlikely(!__pyx_t_10)) goto __pyx_L12_unpacking_failed; __Pyx_GOTREF(__pyx_t_10); index = 1; __pyx_t_12 = __pyx_t_17(__pyx_t_9); if (unlikely(!__pyx_t_12)) goto __pyx_L12_unpacking_failed; __Pyx_GOTREF(__pyx_t_12); if (__Pyx_IternextUnpackEndCheck(__pyx_t_17(__pyx_t_9), 2) < (0)) __PYX_ERR(0, 294, __pyx_L1_error) __pyx_t_17 = NULL; __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0; goto __pyx_L13_unpacking_done; __pyx_L12_unpacking_failed:; __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0; __pyx_t_17 = NULL; if (__Pyx_IterFinish() == 0) __Pyx_RaiseNeedMoreValuesError(index); __PYX_ERR(0, 294, __pyx_L1_error) __pyx_L13_unpacking_done:; } __Pyx_XDECREF_SET(__pyx_v_sign, __pyx_t_10); __pyx_t_10 = 0; __Pyx_XDECREF_SET(__pyx_v_ldet, __pyx_t_12); __pyx_t_12 = 0;
+295: if sign <= 0:
__pyx_t_1 = __Pyx_PyObject_CompareBoolLe_object_int(__pyx_v_sign, __pyx_mstate_global->__pyx_int_0, Py_LE); if (unlikely((__pyx_t_1 < 0))) __PYX_ERR(0, 295, __pyx_L1_error) if (__pyx_t_1) { /* … */ }
+296: return 1e10
{
PyObject *__pyx_temp;
{
__pyx_temp = __pyx_r;
__Pyx_INCREF(__pyx_mstate_global->__pyx_float_1e10);
__pyx_r = __pyx_mstate_global->__pyx_float_1e10;
}
__Pyx_XDECREF(__pyx_temp);
}
goto __pyx_L0;
+297: zt = z_np[t]
__pyx_t_11 = __Pyx_GetItemInt(__pyx_v_z_np, __pyx_v_t, int, 1, __Pyx_PyLong_From_int, 0, 0, 1, __Pyx_ReferenceSharing_OwnStrongReference); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 297, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_11); __Pyx_XDECREF_SET(__pyx_v_zt, __pyx_t_11); __pyx_t_11 = 0;
+298: try:
{
/*try:*/ {
/* … */
}
__Pyx_XDECREF(__pyx_t_18); __pyx_t_18 = 0;
__Pyx_XDECREF(__pyx_t_19); __pyx_t_19 = 0;
__Pyx_XDECREF(__pyx_t_20); __pyx_t_20 = 0;
goto __pyx_L22_try_end;
__pyx_L15_error:;
__Pyx_XDECREF(__pyx_t_10); __pyx_t_10 = 0;
__Pyx_XDECREF(__pyx_t_11); __pyx_t_11 = 0;
__Pyx_XDECREF(__pyx_t_12); __pyx_t_12 = 0;
__Pyx_XDECREF(__pyx_t_14); __pyx_t_14 = 0;
__Pyx_XDECREF(__pyx_t_15); __pyx_t_15 = 0;
__Pyx_XDECREF(__pyx_t_16); __pyx_t_16 = 0;
__Pyx_XDECREF(__pyx_t_9); __pyx_t_9 = 0;
/* … */
__pyx_L17_except_error:;
__Pyx_XGIVEREF(__pyx_t_18);
__Pyx_XGIVEREF(__pyx_t_19);
__Pyx_XGIVEREF(__pyx_t_20);
__Pyx_ExceptionReset(__pyx_t_18, __pyx_t_19, __pyx_t_20);
goto __pyx_L1_error;
__pyx_L18_except_return:;
__Pyx_XGIVEREF(__pyx_t_18);
__Pyx_XGIVEREF(__pyx_t_19);
__Pyx_XGIVEREF(__pyx_t_20);
__Pyx_ExceptionReset(__pyx_t_18, __pyx_t_19, __pyx_t_20);
goto __pyx_L0;
__pyx_L22_try_end:;
}
+299: Rinv = np.linalg.inv(Rt)
__Pyx_GetModuleGlobalName(__pyx_t_10, __pyx_mstate_global->__pyx_n_u_np); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 299, __pyx_L15_error) __Pyx_GOTREF(__pyx_t_10); __pyx_t_9 = __Pyx_PyObject_GetAttrStr(__pyx_t_10, __pyx_mstate_global->__pyx_n_u_linalg); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 299, __pyx_L15_error) __Pyx_GOTREF(__pyx_t_9); __Pyx_DECREF(__pyx_t_10); __pyx_t_10 = 0; __pyx_t_12 = __pyx_t_9; __Pyx_INCREF(__pyx_t_12); __pyx_t_13 = 0; { PyObject *__pyx_callargs[2] = {__pyx_t_12, __pyx_v_Rt}; __pyx_t_11 = __Pyx_PyObject_FastCallMethod((PyObject*)__pyx_mstate_global->__pyx_n_u_inv, __pyx_callargs+__pyx_t_13, (2-__pyx_t_13) | (1*__Pyx_PY_VECTORCALL_ARGUMENTS_OFFSET)); __Pyx_XDECREF(__pyx_t_12); __pyx_t_12 = 0; __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0; if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 299, __pyx_L15_error) __Pyx_GOTREF(__pyx_t_11); } __Pyx_XDECREF_SET(__pyx_v_Rinv, __pyx_t_11); __pyx_t_11 = 0;
+300: except Exception:
__pyx_t_21 = __Pyx_PyErr_ExceptionMatches(((PyObject *)(((PyTypeObject*)PyExc_Exception)))); if (__pyx_t_21) { __Pyx_ErrRestore(0,0,0);
+301: return 1e10
{
PyObject *__pyx_temp;
{
__pyx_temp = __pyx_r;
__Pyx_INCREF(__pyx_mstate_global->__pyx_float_1e10);
__pyx_r = __pyx_mstate_global->__pyx_float_1e10;
}
__Pyx_XDECREF(__pyx_temp);
}
goto __pyx_L18_except_return;
}
goto __pyx_L17_except_error;
+302: ll_py += float(ldet) + float(zt @ Rinv @ zt) - float(zt @ zt)
__pyx_t_22 = __Pyx_PyObject_AsDouble(__pyx_v_ldet); if (unlikely(__PYX_CHECK_FLOAT_EXCEPTION(__pyx_t_22, ((double)((double)-1))) && PyErr_Occurred())) __PYX_ERR(0, 302, __pyx_L1_error) __pyx_t_11 = __Pyx_PyNumber_MatrixMultiply(__pyx_v_zt, __pyx_v_Rinv); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 302, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_11); __pyx_t_9 = __Pyx_PyNumber_MatrixMultiply(__pyx_t_11, __pyx_v_zt); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 302, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_9); __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0; __pyx_t_23 = __Pyx_PyObject_AsDouble(__pyx_t_9); if (unlikely(__PYX_CHECK_FLOAT_EXCEPTION(__pyx_t_23, ((double)((double)-1))) && PyErr_Occurred())) __PYX_ERR(0, 302, __pyx_L1_error) __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0; __pyx_t_9 = __Pyx_PyNumber_MatrixMultiply(__pyx_v_zt, __pyx_v_zt); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 302, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_9); __pyx_t_24 = __Pyx_PyObject_AsDouble(__pyx_t_9); if (unlikely(__PYX_CHECK_FLOAT_EXCEPTION(__pyx_t_24, ((double)((double)-1))) && PyErr_Occurred())) __PYX_ERR(0, 302, __pyx_L1_error) __Pyx_DECREF(__pyx_t_9); __pyx_t_9 = 0; __pyx_v_ll_py = (__pyx_v_ll_py + ((__pyx_t_22 + __pyx_t_23) - __pyx_t_24)); }
+303: return 0.5 * ll_py
__pyx_t_9 = PyFloat_FromDouble((0.5 * __pyx_v_ll_py)); if (unlikely(!__pyx_t_9)) __PYX_ERR(0, 303, __pyx_L1_error) __Pyx_GOTREF(__pyx_t_9); { PyObject *__pyx_temp; { __pyx_temp = __pyx_r; __pyx_r = __pyx_t_9; } __Pyx_XDECREF(__pyx_temp); } __pyx_t_9 = 0; goto __pyx_L0; }