+001: # bonhomme_manresa.pyx

  __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_d, __pyx_n_s_test, __pyx_t_4) < 0) __PYX_ERR(0, 1, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

 002: # cython: profile=False

 003: # cython: language_level=3

 004: # cython: boundscheck=False, wraparound=False, nonecheck=False, cdivision=True

 005: 

 006: """

 007: Grouped-fixed-effects (GFE) estimator – Cython version.

 008: The core numerical work is still done with vectorised NumPy routines,

 009: but the tight Python loops and attribute look-ups are eliminated.

 010: """

 011: 

+012: import numpy as _np

  __pyx_t_2 = __Pyx_ImportDottedModule(__pyx_n_s_numpy, NULL); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 12, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_np, __pyx_t_2) < 0) __PYX_ERR(0, 12, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;

 013: cimport numpy as np

 014: cimport cython

 015: from libc.math cimport INFINITY

 016: 

 017: ctypedef np.float64_t DTYPE_t     # element type we will use everywhere

 018: 

 019: 

 020: # -----------------------------------------------------------------------------

 021: # Helper functions

 022: # -----------------------------------------------------------------------------

+023: cdef inline tuple _get_starting_values(

static CYTHON_INLINE PyObject *__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__get_starting_values(PyArrayObject *__pyx_v_y, PyArrayObject *__pyx_v_x, int __pyx_v_G, int __pyx_v_N, int __pyx_v_K) {
  int __pyx_v_num_start_vars;
  PyArrayObject *__pyx_v_idx_theta = 0;
  PyArrayObject *__pyx_v_X0 = 0;
  PyArrayObject *__pyx_v_Y0 = 0;
  PyArrayObject *__pyx_v_theta_init = 0;
  PyArrayObject *__pyx_v_idx_alpha = 0;
  PyArrayObject *__pyx_v_alpha_init = 0;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_X0;
  __Pyx_Buffer __pyx_pybuffer_X0;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_Y0;
  __Pyx_Buffer __pyx_pybuffer_Y0;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_alpha_init;
  __Pyx_Buffer __pyx_pybuffer_alpha_init;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_idx_alpha;
  __Pyx_Buffer __pyx_pybuffer_idx_alpha;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_idx_theta;
  __Pyx_Buffer __pyx_pybuffer_idx_theta;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_theta_init;
  __Pyx_Buffer __pyx_pybuffer_theta_init;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_x;
  __Pyx_Buffer __pyx_pybuffer_x;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_y;
  __Pyx_Buffer __pyx_pybuffer_y;
  PyObject *__pyx_r = NULL;
  __pyx_pybuffer_idx_theta.pybuffer.buf = NULL;
  __pyx_pybuffer_idx_theta.refcount = 0;
  __pyx_pybuffernd_idx_theta.data = NULL;
  __pyx_pybuffernd_idx_theta.rcbuffer = &__pyx_pybuffer_idx_theta;
  __pyx_pybuffer_X0.pybuffer.buf = NULL;
  __pyx_pybuffer_X0.refcount = 0;
  __pyx_pybuffernd_X0.data = NULL;
  __pyx_pybuffernd_X0.rcbuffer = &__pyx_pybuffer_X0;
  __pyx_pybuffer_Y0.pybuffer.buf = NULL;
  __pyx_pybuffer_Y0.refcount = 0;
  __pyx_pybuffernd_Y0.data = NULL;
  __pyx_pybuffernd_Y0.rcbuffer = &__pyx_pybuffer_Y0;
  __pyx_pybuffer_theta_init.pybuffer.buf = NULL;
  __pyx_pybuffer_theta_init.refcount = 0;
  __pyx_pybuffernd_theta_init.data = NULL;
  __pyx_pybuffernd_theta_init.rcbuffer = &__pyx_pybuffer_theta_init;
  __pyx_pybuffer_idx_alpha.pybuffer.buf = NULL;
  __pyx_pybuffer_idx_alpha.refcount = 0;
  __pyx_pybuffernd_idx_alpha.data = NULL;
  __pyx_pybuffernd_idx_alpha.rcbuffer = &__pyx_pybuffer_idx_alpha;
  __pyx_pybuffer_alpha_init.pybuffer.buf = NULL;
  __pyx_pybuffer_alpha_init.refcount = 0;
  __pyx_pybuffernd_alpha_init.data = NULL;
  __pyx_pybuffernd_alpha_init.rcbuffer = &__pyx_pybuffer_alpha_init;
  __pyx_pybuffer_y.pybuffer.buf = NULL;
  __pyx_pybuffer_y.refcount = 0;
  __pyx_pybuffernd_y.data = NULL;
  __pyx_pybuffernd_y.rcbuffer = &__pyx_pybuffer_y;
  __pyx_pybuffer_x.pybuffer.buf = NULL;
  __pyx_pybuffer_x.refcount = 0;
  __pyx_pybuffernd_x.data = NULL;
  __pyx_pybuffernd_x.rcbuffer = &__pyx_pybuffer_x;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_y.rcbuffer->pybuffer, (PyObject*)__pyx_v_y, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 3, 0, __pyx_stack) == -1)) __PYX_ERR(0, 23, __pyx_L1_error)
  }
  __pyx_pybuffernd_y.diminfo[0].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_y.diminfo[0].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_y.diminfo[1].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_y.diminfo[1].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[1]; __pyx_pybuffernd_y.diminfo[2].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[2]; __pyx_pybuffernd_y.diminfo[2].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[2];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_x.rcbuffer->pybuffer, (PyObject*)__pyx_v_x, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 3, 0, __pyx_stack) == -1)) __PYX_ERR(0, 23, __pyx_L1_error)
  }
  __pyx_pybuffernd_x.diminfo[0].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_x.diminfo[0].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_x.diminfo[1].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_x.diminfo[1].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[1]; __pyx_pybuffernd_x.diminfo[2].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[2]; __pyx_pybuffernd_x.diminfo[2].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[2];
/* … */
  /* 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_12);
  __Pyx_XDECREF(__pyx_t_13);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_X0.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_Y0.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha_init.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_idx_alpha.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_idx_theta.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_theta_init.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_y.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("groupedpaneldatamodels.utils.bonhomme_manresa._get_starting_values", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_X0.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_Y0.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha_init.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_idx_alpha.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_idx_theta.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_theta_init.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_y.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_idx_theta);
  __Pyx_XDECREF((PyObject *)__pyx_v_X0);
  __Pyx_XDECREF((PyObject *)__pyx_v_Y0);
  __Pyx_XDECREF((PyObject *)__pyx_v_theta_init);
  __Pyx_XDECREF((PyObject *)__pyx_v_idx_alpha);
  __Pyx_XDECREF((PyObject *)__pyx_v_alpha_init);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 024:         np.ndarray[DTYPE_t, ndim=3] y,

 025:         np.ndarray[DTYPE_t, ndim=3] x,

 026:         int G, int N, int K):

 027:     """

 028:     Draws random rows and solves a tiny OLS to obtain a first-step θ̂ and α̂.

 029:     """

+030:     cdef int num_start_vars = x.shape[1] + G    # same logic as original code

  __pyx_v_num_start_vars = ((__pyx_f_5numpy_7ndarray_5shape_shape(((PyArrayObject *)__pyx_v_x))[1]) + __pyx_v_G);

 031: 

 032:     # Random draws for θ

+033:     cdef np.ndarray[np.intp_t, ndim=1] idx_theta = _np.random.choice(

  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 33, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_random); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 33, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_choice); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 33, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
/* … */
  __pyx_t_5 = PyTuple_New(2); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 33, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_GIVEREF(__pyx_t_3);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_3)) __PYX_ERR(0, 33, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_4);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_5, 1, __pyx_t_4)) __PYX_ERR(0, 33, __pyx_L1_error);
  __pyx_t_3 = 0;
  __pyx_t_4 = 0;
/* … */
  __pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_5, __pyx_t_4); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 33, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

+034:         N, num_start_vars, replace=False).astype(_np.intp)

  __pyx_t_3 = __Pyx_PyInt_From_int(__pyx_v_N); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 34, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_PyInt_From_int(__pyx_v_num_start_vars); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 34, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
/* … */
  __pyx_t_4 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 34, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  if (PyDict_SetItem(__pyx_t_4, __pyx_n_s_replace, Py_False) < 0) __PYX_ERR(0, 34, __pyx_L1_error)
/* … */
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_astype); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 34, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 34, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_intp); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 34, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = NULL;
  __pyx_t_6 = 0;
  #if CYTHON_UNPACK_METHODS
  if (likely(PyMethod_Check(__pyx_t_4))) {
    __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_4);
    if (likely(__pyx_t_3)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
      __Pyx_INCREF(__pyx_t_3);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_4, function);
      __pyx_t_6 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[2] = {__pyx_t_3, __pyx_t_5};
    __pyx_t_1 = __Pyx_PyObject_FastCall(__pyx_t_4, __pyx_callargs+1-__pyx_t_6, 1+__pyx_t_6);
    __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 34, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  }
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 34, __pyx_L1_error)
  __pyx_t_7 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_idx_theta.rcbuffer->pybuffer, (PyObject*)__pyx_t_7, &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
      __pyx_v_idx_theta = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_idx_theta.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 33, __pyx_L1_error)
    } else {__pyx_pybuffernd_idx_theta.diminfo[0].strides = __pyx_pybuffernd_idx_theta.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_idx_theta.diminfo[0].shape = __pyx_pybuffernd_idx_theta.rcbuffer->pybuffer.shape[0];
    }
  }
  __pyx_t_7 = 0;
  __pyx_v_idx_theta = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

 035: 

+036:     cdef np.ndarray[DTYPE_t, ndim=2] X0 = x[idx_theta].reshape(-1, K)

  __pyx_t_4 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_x), ((PyObject *)__pyx_v_idx_theta)); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 36, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_reshape); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 36, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_PyInt_From_int(__pyx_v_K); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 36, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_3 = NULL;
  __pyx_t_6 = 0;
  #if CYTHON_UNPACK_METHODS
  if (likely(PyMethod_Check(__pyx_t_5))) {
    __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_5);
    if (likely(__pyx_t_3)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_5);
      __Pyx_INCREF(__pyx_t_3);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_5, function);
      __pyx_t_6 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[3] = {__pyx_t_3, __pyx_int_neg_1, __pyx_t_4};
    __pyx_t_1 = __Pyx_PyObject_FastCall(__pyx_t_5, __pyx_callargs+1-__pyx_t_6, 2+__pyx_t_6);
    __Pyx_XDECREF(__pyx_t_3); __pyx_t_3 = 0;
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 36, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  }
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 36, __pyx_L1_error)
  __pyx_t_8 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_X0.rcbuffer->pybuffer, (PyObject*)__pyx_t_8, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_X0 = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_X0.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 36, __pyx_L1_error)
    } else {__pyx_pybuffernd_X0.diminfo[0].strides = __pyx_pybuffernd_X0.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_X0.diminfo[0].shape = __pyx_pybuffernd_X0.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_X0.diminfo[1].strides = __pyx_pybuffernd_X0.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_X0.diminfo[1].shape = __pyx_pybuffernd_X0.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_8 = 0;
  __pyx_v_X0 = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

+037:     cdef np.ndarray[DTYPE_t, ndim=2] Y0 = y[idx_theta].reshape(-1, 1)

  __pyx_t_1 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_y), ((PyObject *)__pyx_v_idx_theta)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 37, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_reshape); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 37, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_5, __pyx_tuple__3, NULL); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 37, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 37, __pyx_L1_error)
  __pyx_t_9 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_Y0.rcbuffer->pybuffer, (PyObject*)__pyx_t_9, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_Y0 = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_Y0.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 37, __pyx_L1_error)
    } else {__pyx_pybuffernd_Y0.diminfo[0].strides = __pyx_pybuffernd_Y0.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_Y0.diminfo[0].shape = __pyx_pybuffernd_Y0.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_Y0.diminfo[1].strides = __pyx_pybuffernd_Y0.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_Y0.diminfo[1].shape = __pyx_pybuffernd_Y0.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_9 = 0;
  __pyx_v_Y0 = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;
/* … */
  __pyx_tuple__3 = PyTuple_Pack(2, __pyx_int_neg_1, __pyx_int_1); if (unlikely(!__pyx_tuple__3)) __PYX_ERR(0, 37, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple__3);
  __Pyx_GIVEREF(__pyx_tuple__3);

+038:     cdef np.ndarray[DTYPE_t, ndim=2] theta_init = _np.linalg.lstsq(X0, Y0, rcond=None)[0]

  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 38, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_linalg); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 38, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_lstsq); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 38, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = PyTuple_New(2); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 38, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_INCREF((PyObject *)__pyx_v_X0);
  __Pyx_GIVEREF((PyObject *)__pyx_v_X0);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_5, 0, ((PyObject *)__pyx_v_X0))) __PYX_ERR(0, 38, __pyx_L1_error);
  __Pyx_INCREF((PyObject *)__pyx_v_Y0);
  __Pyx_GIVEREF((PyObject *)__pyx_v_Y0);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_5, 1, ((PyObject *)__pyx_v_Y0))) __PYX_ERR(0, 38, __pyx_L1_error);
  __pyx_t_4 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 38, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  if (PyDict_SetItem(__pyx_t_4, __pyx_n_s_rcond, Py_None) < 0) __PYX_ERR(0, 38, __pyx_L1_error)
  __pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_1, __pyx_t_5, __pyx_t_4); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 38, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_GetItemInt(__pyx_t_3, 0, long, 1, __Pyx_PyInt_From_long, 0, 0, 0); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 38, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  if (!(likely(((__pyx_t_4) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_4, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 38, __pyx_L1_error)
  __pyx_t_10 = ((PyArrayObject *)__pyx_t_4);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_theta_init.rcbuffer->pybuffer, (PyObject*)__pyx_t_10, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_theta_init = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_theta_init.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 38, __pyx_L1_error)
    } else {__pyx_pybuffernd_theta_init.diminfo[0].strides = __pyx_pybuffernd_theta_init.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_theta_init.diminfo[0].shape = __pyx_pybuffernd_theta_init.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_theta_init.diminfo[1].strides = __pyx_pybuffernd_theta_init.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_theta_init.diminfo[1].shape = __pyx_pybuffernd_theta_init.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_10 = 0;
  __pyx_v_theta_init = ((PyArrayObject *)__pyx_t_4);
  __pyx_t_4 = 0;

 039: 

 040:     # Random draws for α

+041:     cdef np.ndarray[np.intp_t, ndim=1] idx_alpha = _np.random.choice(

  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 41, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_random); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 41, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_choice); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 41, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
/* … */
  __pyx_t_2 = PyTuple_New(2); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 41, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_GIVEREF(__pyx_t_5);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_t_5)) __PYX_ERR(0, 41, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_1);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_2, 1, __pyx_t_1)) __PYX_ERR(0, 41, __pyx_L1_error);
  __pyx_t_5 = 0;
  __pyx_t_1 = 0;
/* … */
  __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_2, __pyx_t_1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 41, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

+042:         N, G, replace=False).astype(_np.intp)

  __pyx_t_5 = __Pyx_PyInt_From_int(__pyx_v_N); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 42, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_G); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 42, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
/* … */
  __pyx_t_1 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 42, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_t_1, __pyx_n_s_replace, Py_False) < 0) __PYX_ERR(0, 42, __pyx_L1_error)
/* … */
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_astype); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 42, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_n_s_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 42, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_intp); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 42, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = NULL;
  __pyx_t_6 = 0;
  #if CYTHON_UNPACK_METHODS
  if (likely(PyMethod_Check(__pyx_t_1))) {
    __pyx_t_5 = PyMethod_GET_SELF(__pyx_t_1);
    if (likely(__pyx_t_5)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_1);
      __Pyx_INCREF(__pyx_t_5);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_1, function);
      __pyx_t_6 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[2] = {__pyx_t_5, __pyx_t_2};
    __pyx_t_4 = __Pyx_PyObject_FastCall(__pyx_t_1, __pyx_callargs+1-__pyx_t_6, 1+__pyx_t_6);
    __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 42, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  }
  if (!(likely(((__pyx_t_4) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_4, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 42, __pyx_L1_error)
  __pyx_t_11 = ((PyArrayObject *)__pyx_t_4);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_idx_alpha.rcbuffer->pybuffer, (PyObject*)__pyx_t_11, &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
      __pyx_v_idx_alpha = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_idx_alpha.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 41, __pyx_L1_error)
    } else {__pyx_pybuffernd_idx_alpha.diminfo[0].strides = __pyx_pybuffernd_idx_alpha.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_idx_alpha.diminfo[0].shape = __pyx_pybuffernd_idx_alpha.rcbuffer->pybuffer.shape[0];
    }
  }
  __pyx_t_11 = 0;
  __pyx_v_idx_alpha = ((PyArrayObject *)__pyx_t_4);
  __pyx_t_4 = 0;

 043: 

 044:     # Y[idx] − X[idx] @ θ

+045:     cdef np.ndarray[DTYPE_t, ndim=2] alpha_init = _np.squeeze(

  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 45, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_squeeze); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 45, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
/* … */
  if (!(likely(((__pyx_t_4) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_4, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 45, __pyx_L1_error)
  __pyx_t_14 = ((PyArrayObject *)__pyx_t_4);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_alpha_init.rcbuffer->pybuffer, (PyObject*)__pyx_t_14, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_alpha_init = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_alpha_init.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 45, __pyx_L1_error)
    } else {__pyx_pybuffernd_alpha_init.diminfo[0].strides = __pyx_pybuffernd_alpha_init.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_alpha_init.diminfo[0].shape = __pyx_pybuffernd_alpha_init.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_alpha_init.diminfo[1].strides = __pyx_pybuffernd_alpha_init.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_alpha_init.diminfo[1].shape = __pyx_pybuffernd_alpha_init.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_14 = 0;
  __pyx_v_alpha_init = ((PyArrayObject *)__pyx_t_4);
  __pyx_t_4 = 0;

+046:         y[idx_alpha] - _np.matmul(x[idx_alpha], theta_init))

  __pyx_t_1 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_y), ((PyObject *)__pyx_v_idx_alpha)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_12 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_matmul); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_12);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_x), ((PyObject *)__pyx_v_idx_alpha)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_13 = NULL;
  __pyx_t_6 = 0;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_12))) {
    __pyx_t_13 = PyMethod_GET_SELF(__pyx_t_12);
    if (likely(__pyx_t_13)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_12);
      __Pyx_INCREF(__pyx_t_13);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_12, function);
      __pyx_t_6 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[3] = {__pyx_t_13, __pyx_t_3, ((PyObject *)__pyx_v_theta_init)};
    __pyx_t_5 = __Pyx_PyObject_FastCall(__pyx_t_12, __pyx_callargs+1-__pyx_t_6, 2+__pyx_t_6);
    __Pyx_XDECREF(__pyx_t_13); __pyx_t_13 = 0;
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 46, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_5);
    __Pyx_DECREF(__pyx_t_12); __pyx_t_12 = 0;
  }
  __pyx_t_12 = PyNumber_Subtract(__pyx_t_1, __pyx_t_5); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 46, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_12);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = NULL;
  __pyx_t_6 = 0;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_2))) {
    __pyx_t_5 = PyMethod_GET_SELF(__pyx_t_2);
    if (likely(__pyx_t_5)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_2);
      __Pyx_INCREF(__pyx_t_5);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_2, function);
      __pyx_t_6 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[2] = {__pyx_t_5, __pyx_t_12};
    __pyx_t_4 = __Pyx_PyObject_FastCall(__pyx_t_2, __pyx_callargs+1-__pyx_t_6, 1+__pyx_t_6);
    __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_12); __pyx_t_12 = 0;
    if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 45, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_4);
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  }

 047: 

+048:     return theta_init, alpha_init

  __Pyx_XDECREF(__pyx_r);
  __pyx_t_4 = PyTuple_New(2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 48, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_INCREF((PyObject *)__pyx_v_theta_init);
  __Pyx_GIVEREF((PyObject *)__pyx_v_theta_init);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_4, 0, ((PyObject *)__pyx_v_theta_init))) __PYX_ERR(0, 48, __pyx_L1_error);
  __Pyx_INCREF((PyObject *)__pyx_v_alpha_init);
  __Pyx_GIVEREF((PyObject *)__pyx_v_alpha_init);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_4, 1, ((PyObject *)__pyx_v_alpha_init))) __PYX_ERR(0, 48, __pyx_L1_error);
  __pyx_r = ((PyObject*)__pyx_t_4);
  __pyx_t_4 = 0;
  goto __pyx_L0;

 049: 

 050: 

+051: cdef inline np.ndarray[np.intp_t, ndim=1] _compute_groupings(

static CYTHON_INLINE PyArrayObject *__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_groupings(PyArrayObject *__pyx_v_res, PyArrayObject *__pyx_v_alpha) {
  PyArrayObject *__pyx_v_dist = 0;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_alpha;
  __Pyx_Buffer __pyx_pybuffer_alpha;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_dist;
  __Pyx_Buffer __pyx_pybuffer_dist;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_res;
  __Pyx_Buffer __pyx_pybuffer_res;
  PyArrayObject *__pyx_r = NULL;
  __pyx_pybuffer_dist.pybuffer.buf = NULL;
  __pyx_pybuffer_dist.refcount = 0;
  __pyx_pybuffernd_dist.data = NULL;
  __pyx_pybuffernd_dist.rcbuffer = &__pyx_pybuffer_dist;
  __pyx_pybuffer_res.pybuffer.buf = NULL;
  __pyx_pybuffer_res.refcount = 0;
  __pyx_pybuffernd_res.data = NULL;
  __pyx_pybuffernd_res.rcbuffer = &__pyx_pybuffer_res;
  __pyx_pybuffer_alpha.pybuffer.buf = NULL;
  __pyx_pybuffer_alpha.refcount = 0;
  __pyx_pybuffernd_alpha.data = NULL;
  __pyx_pybuffernd_alpha.rcbuffer = &__pyx_pybuffer_alpha;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_res.rcbuffer->pybuffer, (PyObject*)__pyx_v_res, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 51, __pyx_L1_error)
  }
  __pyx_pybuffernd_res.diminfo[0].strides = __pyx_pybuffernd_res.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_res.diminfo[0].shape = __pyx_pybuffernd_res.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_res.diminfo[1].strides = __pyx_pybuffernd_res.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_res.diminfo[1].shape = __pyx_pybuffernd_res.rcbuffer->pybuffer.shape[1];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer, (PyObject*)__pyx_v_alpha, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 51, __pyx_L1_error)
  }
  __pyx_pybuffernd_alpha.diminfo[0].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_alpha.diminfo[0].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_alpha.diminfo[1].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_alpha.diminfo[1].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[1];
/* … */
  /* 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_5);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_dist.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_res.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("groupedpaneldatamodels.utils.bonhomme_manresa._compute_groupings", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_dist.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_res.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_dist);
  __Pyx_XGIVEREF((PyObject *)__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 052:         np.ndarray[DTYPE_t, ndim=2] res,

 053:         np.ndarray[DTYPE_t, ndim=2] alpha):

 054:     """

 055:     For every observation choose the closest α (Euclidean distance).

 056:     """

 057:     # res shape: (N, K)  -> add singleton axis so broadcasting works

 058:     # alpha shape: (G, K)

 059:     cdef np.ndarray[DTYPE_t, ndim=2] dist = (

+060:         (res[None, :, :] - alpha[:, None, :]) ** 2).sum(axis=2)

  __pyx_slice__4 = PySlice_New(Py_None, Py_None, Py_None); if (unlikely(!__pyx_slice__4)) __PYX_ERR(0, 60, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_slice__4);
  __Pyx_GIVEREF(__pyx_slice__4);
/* … */
  __pyx_t_1 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_res), __pyx_tuple__5); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 60, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_tuple__5 = PyTuple_Pack(3, Py_None, __pyx_slice__4, __pyx_slice__4); if (unlikely(!__pyx_tuple__5)) __PYX_ERR(0, 60, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple__5);
  __Pyx_GIVEREF(__pyx_tuple__5);
  __pyx_t_2 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_alpha), __pyx_tuple__6); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 60, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = PyNumber_Subtract(__pyx_t_1, __pyx_t_2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 60, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = PyNumber_Power(__pyx_t_3, __pyx_int_2, Py_None); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 60, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_sum); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 60, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 60, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  if (PyDict_SetItem(__pyx_t_2, __pyx_n_s_axis, __pyx_int_2) < 0) __PYX_ERR(0, 60, __pyx_L1_error)
  __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_empty_tuple, __pyx_t_2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 60, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 60, __pyx_L1_error)
  __pyx_t_4 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_dist.rcbuffer->pybuffer, (PyObject*)__pyx_t_4, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_dist = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_dist.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 59, __pyx_L1_error)
    } else {__pyx_pybuffernd_dist.diminfo[0].strides = __pyx_pybuffernd_dist.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_dist.diminfo[0].shape = __pyx_pybuffernd_dist.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_dist.diminfo[1].strides = __pyx_pybuffernd_dist.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_dist.diminfo[1].shape = __pyx_pybuffernd_dist.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_4 = 0;
  __pyx_v_dist = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;
  __pyx_tuple__6 = PyTuple_Pack(3, __pyx_slice__4, Py_None, __pyx_slice__4); if (unlikely(!__pyx_tuple__6)) __PYX_ERR(0, 60, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple__6);
  __Pyx_GIVEREF(__pyx_tuple__6);

 061: 

+062:     return _np.argmin(dist, axis=0)      # (N,) int array

  __Pyx_XDECREF((PyObject *)__pyx_r);
  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 62, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_argmin); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 62, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = PyTuple_New(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 62, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_INCREF((PyObject *)__pyx_v_dist);
  __Pyx_GIVEREF((PyObject *)__pyx_v_dist);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_1, 0, ((PyObject *)__pyx_v_dist))) __PYX_ERR(0, 62, __pyx_L1_error);
  __pyx_t_3 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 62, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  if (PyDict_SetItem(__pyx_t_3, __pyx_n_s_axis, __pyx_int_0) < 0) __PYX_ERR(0, 62, __pyx_L1_error)
  __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_1, __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 62, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  if (!(likely(((__pyx_t_5) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_5, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 62, __pyx_L1_error)
  __pyx_r = ((PyArrayObject *)__pyx_t_5);
  __pyx_t_5 = 0;
  goto __pyx_L0;

 063: 

 064: 

+065: cdef inline np.ndarray[DTYPE_t, ndim=2] _compute_alpha(

static CYTHON_INLINE PyArrayObject *__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_alpha(PyArrayObject *__pyx_v_res, PyArrayObject *__pyx_v_g, int __pyx_v_G) {
  PyArrayObject *__pyx_v_counts = 0;
  PyArrayObject *__pyx_v_sums = 0;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_counts;
  __Pyx_Buffer __pyx_pybuffer_counts;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_g;
  __Pyx_Buffer __pyx_pybuffer_g;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_res;
  __Pyx_Buffer __pyx_pybuffer_res;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_sums;
  __Pyx_Buffer __pyx_pybuffer_sums;
  PyArrayObject *__pyx_r = NULL;
  __pyx_pybuffer_counts.pybuffer.buf = NULL;
  __pyx_pybuffer_counts.refcount = 0;
  __pyx_pybuffernd_counts.data = NULL;
  __pyx_pybuffernd_counts.rcbuffer = &__pyx_pybuffer_counts;
  __pyx_pybuffer_sums.pybuffer.buf = NULL;
  __pyx_pybuffer_sums.refcount = 0;
  __pyx_pybuffernd_sums.data = NULL;
  __pyx_pybuffernd_sums.rcbuffer = &__pyx_pybuffer_sums;
  __pyx_pybuffer_res.pybuffer.buf = NULL;
  __pyx_pybuffer_res.refcount = 0;
  __pyx_pybuffernd_res.data = NULL;
  __pyx_pybuffernd_res.rcbuffer = &__pyx_pybuffer_res;
  __pyx_pybuffer_g.pybuffer.buf = NULL;
  __pyx_pybuffer_g.refcount = 0;
  __pyx_pybuffernd_g.data = NULL;
  __pyx_pybuffernd_g.rcbuffer = &__pyx_pybuffer_g;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_res.rcbuffer->pybuffer, (PyObject*)__pyx_v_res, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 65, __pyx_L1_error)
  }
  __pyx_pybuffernd_res.diminfo[0].strides = __pyx_pybuffernd_res.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_res.diminfo[0].shape = __pyx_pybuffernd_res.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_res.diminfo[1].strides = __pyx_pybuffernd_res.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_res.diminfo[1].shape = __pyx_pybuffernd_res.rcbuffer->pybuffer.shape[1];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_g.rcbuffer->pybuffer, (PyObject*)__pyx_v_g, &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 65, __pyx_L1_error)
  }
  __pyx_pybuffernd_g.diminfo[0].strides = __pyx_pybuffernd_g.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_g.diminfo[0].shape = __pyx_pybuffernd_g.rcbuffer->pybuffer.shape[0];
/* … */
  /* 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);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_counts.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_res.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_sums.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("groupedpaneldatamodels.utils.bonhomme_manresa._compute_alpha", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_counts.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_res.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_sums.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_counts);
  __Pyx_XDECREF((PyObject *)__pyx_v_sums);
  __Pyx_XGIVEREF((PyObject *)__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 066:         np.ndarray[DTYPE_t, ndim=2] res,

 067:         np.ndarray[np.intp_t, ndim=1] g,

 068:         int G):

 069:     """

 070:     α̂_g = average residual in group g.

 071:     """

+072:     cdef np.ndarray[np.intp_t, ndim=1] counts = _np.bincount(

  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 72, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_bincount); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 72, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* … */
  __pyx_t_2 = PyTuple_New(1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 72, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_INCREF((PyObject *)__pyx_v_g);
  __Pyx_GIVEREF((PyObject *)__pyx_v_g);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_2, 0, ((PyObject *)__pyx_v_g))) __PYX_ERR(0, 72, __pyx_L1_error);
/* … */
  __pyx_t_5 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_2, __pyx_t_4); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 72, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

+073:         g, minlength=G).astype(_np.intp)          # (G,)

  __pyx_t_4 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 73, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_5 = __Pyx_PyInt_From_int(__pyx_v_G); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 73, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  if (PyDict_SetItem(__pyx_t_4, __pyx_n_s_minlength, __pyx_t_5) < 0) __PYX_ERR(0, 73, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
/* … */
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_astype); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 73, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_n_s_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 73, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_intp); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 73, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = NULL;
  __pyx_t_6 = 0;
  #if CYTHON_UNPACK_METHODS
  if (likely(PyMethod_Check(__pyx_t_4))) {
    __pyx_t_5 = PyMethod_GET_SELF(__pyx_t_4);
    if (likely(__pyx_t_5)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4);
      __Pyx_INCREF(__pyx_t_5);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_4, function);
      __pyx_t_6 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[2] = {__pyx_t_5, __pyx_t_2};
    __pyx_t_1 = __Pyx_PyObject_FastCall(__pyx_t_4, __pyx_callargs+1-__pyx_t_6, 1+__pyx_t_6);
    __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 73, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  }
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 73, __pyx_L1_error)
  __pyx_t_7 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_counts.rcbuffer->pybuffer, (PyObject*)__pyx_t_7, &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
      __pyx_v_counts = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_counts.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 72, __pyx_L1_error)
    } else {__pyx_pybuffernd_counts.diminfo[0].strides = __pyx_pybuffernd_counts.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_counts.diminfo[0].shape = __pyx_pybuffernd_counts.rcbuffer->pybuffer.shape[0];
    }
  }
  __pyx_t_7 = 0;
  __pyx_v_counts = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

 074: 

+075:     cdef np.ndarray[DTYPE_t, ndim=2] sums = _np.zeros(

  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 75, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 75, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
/* … */
  __pyx_t_2 = PyTuple_New(1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 75, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_GIVEREF(__pyx_t_5);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_t_5)) __PYX_ERR(0, 75, __pyx_L1_error);
  __pyx_t_5 = 0;
/* … */
  __pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_2, __pyx_t_5); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 75, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 75, __pyx_L1_error)
  __pyx_t_8 = ((PyArrayObject *)__pyx_t_3);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_sums.rcbuffer->pybuffer, (PyObject*)__pyx_t_8, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_sums = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_sums.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 75, __pyx_L1_error)
    } else {__pyx_pybuffernd_sums.diminfo[0].strides = __pyx_pybuffernd_sums.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_sums.diminfo[0].shape = __pyx_pybuffernd_sums.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_sums.diminfo[1].strides = __pyx_pybuffernd_sums.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_sums.diminfo[1].shape = __pyx_pybuffernd_sums.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_8 = 0;
  __pyx_v_sums = ((PyArrayObject *)__pyx_t_3);
  __pyx_t_3 = 0;

+076:         (G, res.shape[1]), dtype=_np.float64)              # (G, K)

  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_G); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 76, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_2 = __Pyx_PyInt_From_npy_intp((__pyx_f_5numpy_7ndarray_5shape_shape(((PyArrayObject *)__pyx_v_res))[1])); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 76, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_5 = PyTuple_New(2); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 76, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_GIVEREF(__pyx_t_1);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_1)) __PYX_ERR(0, 76, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_2);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_5, 1, __pyx_t_2)) __PYX_ERR(0, 76, __pyx_L1_error);
  __pyx_t_1 = 0;
  __pyx_t_2 = 0;
/* … */
  __pyx_t_5 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 76, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 76, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_float64); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 76, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  if (PyDict_SetItem(__pyx_t_5, __pyx_n_s_dtype, __pyx_t_3) < 0) __PYX_ERR(0, 76, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;

 077: 

+078:     _np.add.at(sums, g, res)                               # in-place accumulation

  __Pyx_GetModuleGlobalName(__pyx_t_5, __pyx_n_s_np); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 78, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_add); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 78, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_at); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 78, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = NULL;
  __pyx_t_6 = 0;
  #if CYTHON_UNPACK_METHODS
  if (likely(PyMethod_Check(__pyx_t_5))) {
    __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_5);
    if (likely(__pyx_t_2)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_5);
      __Pyx_INCREF(__pyx_t_2);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_5, function);
      __pyx_t_6 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[4] = {__pyx_t_2, ((PyObject *)__pyx_v_sums), ((PyObject *)__pyx_v_g), ((PyObject *)__pyx_v_res)};
    __pyx_t_3 = __Pyx_PyObject_FastCall(__pyx_t_5, __pyx_callargs+1-__pyx_t_6, 3+__pyx_t_6);
    __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
    if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 78, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  }
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;

+079:     return sums / counts[:, None]                          # broadcast along columns

  __Pyx_XDECREF((PyObject *)__pyx_r);
  __pyx_t_3 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_counts), __pyx_tuple__7); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 79, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_5 = __Pyx_PyNumber_Divide(((PyObject *)__pyx_v_sums), __pyx_t_3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 79, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  if (!(likely(((__pyx_t_5) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_5, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 79, __pyx_L1_error)
  __pyx_r = ((PyArrayObject *)__pyx_t_5);
  __pyx_t_5 = 0;
  goto __pyx_L0;
/* … */
  __pyx_tuple__7 = PyTuple_Pack(2, __pyx_slice__4, Py_None); if (unlikely(!__pyx_tuple__7)) __PYX_ERR(0, 79, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple__7);
  __Pyx_GIVEREF(__pyx_tuple__7);

 080: 

 081: # ---------------------------------------------------------------------------

 082: # Fast α̂ update using explicit C loops (no np.add.at)

 083: # ---------------------------------------------------------------------------

 084: @cython.boundscheck(False)

 085: @cython.wraparound(False)

+086: cdef inline np.ndarray[DTYPE_t, ndim=2] _compute_alpha_fast(

static CYTHON_INLINE PyArrayObject *__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_alpha_fast(PyArrayObject *__pyx_v_res, PyArrayObject *__pyx_v_g, int __pyx_v_G) {
  int __pyx_v_N;
  int __pyx_v_L;
  PyArrayObject *__pyx_v_sums = 0;
  PyArrayObject *__pyx_v_counts = 0;
  int __pyx_v_i;
  int __pyx_v_l;
  int __pyx_v_gi;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_counts;
  __Pyx_Buffer __pyx_pybuffer_counts;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_g;
  __Pyx_Buffer __pyx_pybuffer_g;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_res;
  __Pyx_Buffer __pyx_pybuffer_res;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_sums;
  __Pyx_Buffer __pyx_pybuffer_sums;
  PyArrayObject *__pyx_r = NULL;
  __pyx_pybuffer_sums.pybuffer.buf = NULL;
  __pyx_pybuffer_sums.refcount = 0;
  __pyx_pybuffernd_sums.data = NULL;
  __pyx_pybuffernd_sums.rcbuffer = &__pyx_pybuffer_sums;
  __pyx_pybuffer_counts.pybuffer.buf = NULL;
  __pyx_pybuffer_counts.refcount = 0;
  __pyx_pybuffernd_counts.data = NULL;
  __pyx_pybuffernd_counts.rcbuffer = &__pyx_pybuffer_counts;
  __pyx_pybuffer_res.pybuffer.buf = NULL;
  __pyx_pybuffer_res.refcount = 0;
  __pyx_pybuffernd_res.data = NULL;
  __pyx_pybuffernd_res.rcbuffer = &__pyx_pybuffer_res;
  __pyx_pybuffer_g.pybuffer.buf = NULL;
  __pyx_pybuffer_g.refcount = 0;
  __pyx_pybuffernd_g.data = NULL;
  __pyx_pybuffernd_g.rcbuffer = &__pyx_pybuffer_g;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_res.rcbuffer->pybuffer, (PyObject*)__pyx_v_res, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 86, __pyx_L1_error)
  }
  __pyx_pybuffernd_res.diminfo[0].strides = __pyx_pybuffernd_res.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_res.diminfo[0].shape = __pyx_pybuffernd_res.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_res.diminfo[1].strides = __pyx_pybuffernd_res.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_res.diminfo[1].shape = __pyx_pybuffernd_res.rcbuffer->pybuffer.shape[1];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_g.rcbuffer->pybuffer, (PyObject*)__pyx_v_g, &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 86, __pyx_L1_error)
  }
  __pyx_pybuffernd_g.diminfo[0].strides = __pyx_pybuffernd_g.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_g.diminfo[0].shape = __pyx_pybuffernd_g.rcbuffer->pybuffer.shape[0];
/* … */
  /* 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_6);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_counts.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_res.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_sums.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("groupedpaneldatamodels.utils.bonhomme_manresa._compute_alpha_fast", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_counts.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_res.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_sums.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_sums);
  __Pyx_XDECREF((PyObject *)__pyx_v_counts);
  __Pyx_XGIVEREF((PyObject *)__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 087:         np.ndarray[DTYPE_t, ndim=2] res,

 088:         np.ndarray[np.intp_t, ndim=1] g,

 089:         int G):

 090:     """

 091:     Faster α̂ update using tight nested loops to accumulate group means.

 092:     """

+093:     cdef int N = res.shape[0]

  __pyx_v_N = (__pyx_f_5numpy_7ndarray_5shape_shape(((PyArrayObject *)__pyx_v_res))[0]);

+094:     cdef int L = res.shape[1]

  __pyx_v_L = (__pyx_f_5numpy_7ndarray_5shape_shape(((PyArrayObject *)__pyx_v_res))[1]);

+095:     cdef np.ndarray[DTYPE_t, ndim=2] sums = _np.zeros((G, L), dtype=res.dtype)

  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 95, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 95, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_G); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 95, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_3 = __Pyx_PyInt_From_int(__pyx_v_L); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 95, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = PyTuple_New(2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 95, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_GIVEREF(__pyx_t_1);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_1)) __PYX_ERR(0, 95, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_3);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_4, 1, __pyx_t_3)) __PYX_ERR(0, 95, __pyx_L1_error);
  __pyx_t_1 = 0;
  __pyx_t_3 = 0;
  __pyx_t_3 = PyTuple_New(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 95, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_GIVEREF(__pyx_t_4);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_4)) __PYX_ERR(0, 95, __pyx_L1_error);
  __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 95, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(((PyObject *)__pyx_v_res), __pyx_n_s_dtype); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 95, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (PyDict_SetItem(__pyx_t_4, __pyx_n_s_dtype, __pyx_t_1) < 0) __PYX_ERR(0, 95, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_3, __pyx_t_4); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 95, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 95, __pyx_L1_error)
  __pyx_t_5 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_sums.rcbuffer->pybuffer, (PyObject*)__pyx_t_5, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_sums = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_sums.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 95, __pyx_L1_error)
    } else {__pyx_pybuffernd_sums.diminfo[0].strides = __pyx_pybuffernd_sums.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_sums.diminfo[0].shape = __pyx_pybuffernd_sums.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_sums.diminfo[1].strides = __pyx_pybuffernd_sums.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_sums.diminfo[1].shape = __pyx_pybuffernd_sums.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_5 = 0;
  __pyx_v_sums = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

+096:     cdef np.ndarray[np.int64_t, ndim=1] counts = _np.zeros(G, dtype=_np.int64)

  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 96, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_zeros); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 96, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_G); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 96, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_3 = PyTuple_New(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 96, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_GIVEREF(__pyx_t_1);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_1)) __PYX_ERR(0, 96, __pyx_L1_error);
  __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 96, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 96, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_6 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_int64); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 96, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  if (PyDict_SetItem(__pyx_t_1, __pyx_n_s_dtype, __pyx_t_6) < 0) __PYX_ERR(0, 96, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
  __pyx_t_6 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_3, __pyx_t_1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 96, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_6);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  if (!(likely(((__pyx_t_6) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_6, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 96, __pyx_L1_error)
  __pyx_t_7 = ((PyArrayObject *)__pyx_t_6);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_counts.rcbuffer->pybuffer, (PyObject*)__pyx_t_7, &__Pyx_TypeInfo_nn___pyx_t_5numpy_int64_t, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 1, 0, __pyx_stack) == -1)) {
      __pyx_v_counts = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_counts.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 96, __pyx_L1_error)
    } else {__pyx_pybuffernd_counts.diminfo[0].strides = __pyx_pybuffernd_counts.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_counts.diminfo[0].shape = __pyx_pybuffernd_counts.rcbuffer->pybuffer.shape[0];
    }
  }
  __pyx_t_7 = 0;
  __pyx_v_counts = ((PyArrayObject *)__pyx_t_6);
  __pyx_t_6 = 0;

 097:     cdef int i, l, gi

+098:     for i in range(N):

  __pyx_t_8 = __pyx_v_N;
  __pyx_t_9 = __pyx_t_8;
  for (__pyx_t_10 = 0; __pyx_t_10 < __pyx_t_9; __pyx_t_10+=1) {
    __pyx_v_i = __pyx_t_10;

+099:         gi = <int>g[i]

    __pyx_t_11 = __pyx_v_i;
    __pyx_v_gi = ((int)(*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_intp_t *, __pyx_pybuffernd_g.rcbuffer->pybuffer.buf, __pyx_t_11, __pyx_pybuffernd_g.diminfo[0].strides)));

+100:         counts[gi] += 1

    __pyx_t_11 = __pyx_v_gi;
    *__Pyx_BufPtrStrided1d(__pyx_t_5numpy_int64_t *, __pyx_pybuffernd_counts.rcbuffer->pybuffer.buf, __pyx_t_11, __pyx_pybuffernd_counts.diminfo[0].strides) += 1;

+101:         for l in range(L):

    __pyx_t_12 = __pyx_v_L;
    __pyx_t_13 = __pyx_t_12;
    for (__pyx_t_14 = 0; __pyx_t_14 < __pyx_t_13; __pyx_t_14+=1) {
      __pyx_v_l = __pyx_t_14;

+102:             sums[gi, l] += res[i, l]

      __pyx_t_11 = __pyx_v_i;
      __pyx_t_15 = __pyx_v_l;
      __pyx_t_16 = __pyx_v_gi;
      __pyx_t_17 = __pyx_v_l;
      *__Pyx_BufPtrStrided2d(__pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t *, __pyx_pybuffernd_sums.rcbuffer->pybuffer.buf, __pyx_t_16, __pyx_pybuffernd_sums.diminfo[0].strides, __pyx_t_17, __pyx_pybuffernd_sums.diminfo[1].strides) += (*__Pyx_BufPtrStrided2d(__pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t *, __pyx_pybuffernd_res.rcbuffer->pybuffer.buf, __pyx_t_11, __pyx_pybuffernd_res.diminfo[0].strides, __pyx_t_15, __pyx_pybuffernd_res.diminfo[1].strides));
    }
  }

+103:     for gi in range(G):

  __pyx_t_8 = __pyx_v_G;
  __pyx_t_9 = __pyx_t_8;
  for (__pyx_t_10 = 0; __pyx_t_10 < __pyx_t_9; __pyx_t_10+=1) {
    __pyx_v_gi = __pyx_t_10;

+104:         if counts[gi]:

    __pyx_t_15 = __pyx_v_gi;
    __pyx_t_18 = ((*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_int64_t *, __pyx_pybuffernd_counts.rcbuffer->pybuffer.buf, __pyx_t_15, __pyx_pybuffernd_counts.diminfo[0].strides)) != 0);
    if (__pyx_t_18) {
/* … */
    }
  }

+105:             sums[gi, :] /= counts[gi]

      __pyx_t_6 = __Pyx_PyInt_From_int(__pyx_v_gi); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 105, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_6);
      __pyx_t_1 = PyTuple_New(2); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 105, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_1);
      __Pyx_GIVEREF(__pyx_t_6);
      if (__Pyx_PyTuple_SET_ITEM(__pyx_t_1, 0, __pyx_t_6)) __PYX_ERR(0, 105, __pyx_L1_error);
      __Pyx_INCREF(__pyx_slice__4);
      __Pyx_GIVEREF(__pyx_slice__4);
      if (__Pyx_PyTuple_SET_ITEM(__pyx_t_1, 1, __pyx_slice__4)) __PYX_ERR(0, 105, __pyx_L1_error);
      __pyx_t_6 = 0;
      __pyx_t_6 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_sums), __pyx_t_1); if (unlikely(!__pyx_t_6)) __PYX_ERR(0, 105, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_6);
      __pyx_t_15 = __pyx_v_gi;
      __pyx_t_3 = __Pyx_PyInt_From_npy_int64((*__Pyx_BufPtrStrided1d(__pyx_t_5numpy_int64_t *, __pyx_pybuffernd_counts.rcbuffer->pybuffer.buf, __pyx_t_15, __pyx_pybuffernd_counts.diminfo[0].strides))); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 105, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_3);
      __pyx_t_4 = __Pyx_PyNumber_InPlaceDivide(__pyx_t_6, __pyx_t_3); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 105, __pyx_L1_error)
      __Pyx_GOTREF(__pyx_t_4);
      __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0;
      __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
      if (unlikely((PyObject_SetItem(((PyObject *)__pyx_v_sums), __pyx_t_1, __pyx_t_4) < 0))) __PYX_ERR(0, 105, __pyx_L1_error)
      __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
      __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;

+106:     return sums

  __Pyx_XDECREF((PyObject *)__pyx_r);
  __Pyx_INCREF((PyObject *)__pyx_v_sums);
  __pyx_r = ((PyArrayObject *)__pyx_v_sums);
  goto __pyx_L0;

 107: 

+108: cdef inline np.ndarray[DTYPE_t, ndim=2] _compute_theta(

static CYTHON_INLINE PyArrayObject *__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_theta(PyArrayObject *__pyx_v_x, PyArrayObject *__pyx_v_y, PyArrayObject *__pyx_v_alpha, PyArrayObject *__pyx_v_g) {
  int __pyx_v_K;
  PyArrayObject *__pyx_v_theta = 0;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_alpha;
  __Pyx_Buffer __pyx_pybuffer_alpha;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_g;
  __Pyx_Buffer __pyx_pybuffer_g;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_theta;
  __Pyx_Buffer __pyx_pybuffer_theta;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_x;
  __Pyx_Buffer __pyx_pybuffer_x;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_y;
  __Pyx_Buffer __pyx_pybuffer_y;
  PyArrayObject *__pyx_r = NULL;
  __pyx_pybuffer_theta.pybuffer.buf = NULL;
  __pyx_pybuffer_theta.refcount = 0;
  __pyx_pybuffernd_theta.data = NULL;
  __pyx_pybuffernd_theta.rcbuffer = &__pyx_pybuffer_theta;
  __pyx_pybuffer_x.pybuffer.buf = NULL;
  __pyx_pybuffer_x.refcount = 0;
  __pyx_pybuffernd_x.data = NULL;
  __pyx_pybuffernd_x.rcbuffer = &__pyx_pybuffer_x;
  __pyx_pybuffer_y.pybuffer.buf = NULL;
  __pyx_pybuffer_y.refcount = 0;
  __pyx_pybuffernd_y.data = NULL;
  __pyx_pybuffernd_y.rcbuffer = &__pyx_pybuffer_y;
  __pyx_pybuffer_alpha.pybuffer.buf = NULL;
  __pyx_pybuffer_alpha.refcount = 0;
  __pyx_pybuffernd_alpha.data = NULL;
  __pyx_pybuffernd_alpha.rcbuffer = &__pyx_pybuffer_alpha;
  __pyx_pybuffer_g.pybuffer.buf = NULL;
  __pyx_pybuffer_g.refcount = 0;
  __pyx_pybuffernd_g.data = NULL;
  __pyx_pybuffernd_g.rcbuffer = &__pyx_pybuffer_g;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_x.rcbuffer->pybuffer, (PyObject*)__pyx_v_x, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 3, 0, __pyx_stack) == -1)) __PYX_ERR(0, 108, __pyx_L1_error)
  }
  __pyx_pybuffernd_x.diminfo[0].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_x.diminfo[0].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_x.diminfo[1].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_x.diminfo[1].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[1]; __pyx_pybuffernd_x.diminfo[2].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[2]; __pyx_pybuffernd_x.diminfo[2].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[2];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_y.rcbuffer->pybuffer, (PyObject*)__pyx_v_y, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 3, 0, __pyx_stack) == -1)) __PYX_ERR(0, 108, __pyx_L1_error)
  }
  __pyx_pybuffernd_y.diminfo[0].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_y.diminfo[0].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_y.diminfo[1].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_y.diminfo[1].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[1]; __pyx_pybuffernd_y.diminfo[2].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[2]; __pyx_pybuffernd_y.diminfo[2].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[2];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer, (PyObject*)__pyx_v_alpha, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 108, __pyx_L1_error)
  }
  __pyx_pybuffernd_alpha.diminfo[0].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_alpha.diminfo[0].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_alpha.diminfo[1].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_alpha.diminfo[1].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[1];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_g.rcbuffer->pybuffer, (PyObject*)__pyx_v_g, &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 108, __pyx_L1_error)
  }
  __pyx_pybuffernd_g.diminfo[0].strides = __pyx_pybuffernd_g.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_g.diminfo[0].shape = __pyx_pybuffernd_g.rcbuffer->pybuffer.shape[0];
/* … */
  /* 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);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_theta.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_y.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("groupedpaneldatamodels.utils.bonhomme_manresa._compute_theta", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_theta.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_y.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_theta);
  __Pyx_XGIVEREF((PyObject *)__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 109:         np.ndarray[DTYPE_t, ndim=3] x,

 110:         np.ndarray[DTYPE_t, ndim=3] y,

 111:         np.ndarray[DTYPE_t, ndim=2] alpha,

 112:         np.ndarray[np.intp_t, ndim=1] g):

 113:     """

 114:     θ̂ from pooled OLS of (y − α_g) on x.

 115:     """

+116:     cdef int K = x.shape[2]

  __pyx_v_K = (__pyx_f_5numpy_7ndarray_5shape_shape(((PyArrayObject *)__pyx_v_x))[2]);

 117:     # broadcast α_g so that (N,L,1) − (N,L,1) is well‑defined

+118:     cdef np.ndarray[DTYPE_t, ndim=2] theta = _np.linalg.lstsq(

  __Pyx_GetModuleGlobalName(__pyx_t_1, __pyx_n_s_np); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 118, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_1, __pyx_n_s_linalg); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 118, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_lstsq); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 118, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* … */
  __pyx_t_4 = PyTuple_New(2); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 118, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_GIVEREF(__pyx_t_2);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_2)) __PYX_ERR(0, 118, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_3);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_4, 1, __pyx_t_3)) __PYX_ERR(0, 118, __pyx_L1_error);
  __pyx_t_2 = 0;
  __pyx_t_3 = 0;
/* … */
  __pyx_t_2 = __Pyx_PyObject_Call(__pyx_t_1, __pyx_t_4, __pyx_t_3); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 118, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;

+119:         x.reshape(-1, K),

  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(((PyObject *)__pyx_v_x), __pyx_n_s_reshape); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 119, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_PyInt_From_int(__pyx_v_K); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 119, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_5 = NULL;
  __pyx_t_6 = 0;
  #if CYTHON_UNPACK_METHODS
  if (likely(PyMethod_Check(__pyx_t_3))) {
    __pyx_t_5 = PyMethod_GET_SELF(__pyx_t_3);
    if (likely(__pyx_t_5)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
      __Pyx_INCREF(__pyx_t_5);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_3, function);
      __pyx_t_6 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[3] = {__pyx_t_5, __pyx_int_neg_1, __pyx_t_4};
    __pyx_t_2 = __Pyx_PyObject_FastCall(__pyx_t_3, __pyx_callargs+1-__pyx_t_6, 2+__pyx_t_6);
    __Pyx_XDECREF(__pyx_t_5); __pyx_t_5 = 0;
    __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
    if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 119, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_2);
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  }

+120:         (y - alpha[g][:, :, None]).reshape(-1, 1),

  __pyx_t_3 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_alpha), ((PyObject *)__pyx_v_g)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 120, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = __Pyx_PyObject_GetItem(__pyx_t_3, __pyx_tuple__8); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 120, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = PyNumber_Subtract(((PyObject *)__pyx_v_y), __pyx_t_4); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 120, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_reshape); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 120, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_tuple__3, NULL); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 120, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
/* … */
  __pyx_tuple__8 = PyTuple_Pack(3, __pyx_slice__4, __pyx_slice__4, Py_None); if (unlikely(!__pyx_tuple__8)) __PYX_ERR(0, 120, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple__8);
  __Pyx_GIVEREF(__pyx_tuple__8);

+121:         rcond=None

  __pyx_t_3 = __Pyx_PyDict_NewPresized(1); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 121, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  if (PyDict_SetItem(__pyx_t_3, __pyx_n_s_rcond, Py_None) < 0) __PYX_ERR(0, 121, __pyx_L1_error)

+122:     )[0]

  __pyx_t_3 = __Pyx_GetItemInt(__pyx_t_2, 0, long, 1, __Pyx_PyInt_From_long, 0, 0, 0); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 122, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 122, __pyx_L1_error)
  __pyx_t_7 = ((PyArrayObject *)__pyx_t_3);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_theta.rcbuffer->pybuffer, (PyObject*)__pyx_t_7, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_theta = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_theta.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 118, __pyx_L1_error)
    } else {__pyx_pybuffernd_theta.diminfo[0].strides = __pyx_pybuffernd_theta.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_theta.diminfo[0].shape = __pyx_pybuffernd_theta.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_theta.diminfo[1].strides = __pyx_pybuffernd_theta.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_theta.diminfo[1].shape = __pyx_pybuffernd_theta.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_7 = 0;
  __pyx_v_theta = ((PyArrayObject *)__pyx_t_3);
  __pyx_t_3 = 0;

+123:     return theta

  __Pyx_XDECREF((PyObject *)__pyx_r);
  __Pyx_INCREF((PyObject *)__pyx_v_theta);
  __pyx_r = ((PyArrayObject *)__pyx_v_theta);
  goto __pyx_L0;

 124: 

 125: 

 126: # ---------------------------------------------------------------------------

 127: # Fast θ̂ update using precomputed X'X and X'

 128: # ---------------------------------------------------------------------------

+129: cdef inline np.ndarray[DTYPE_t, ndim=2] _compute_theta_fast(

static CYTHON_INLINE PyArrayObject *__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_theta_fast(PyArrayObject *__pyx_v_XtX, PyArrayObject *__pyx_v_Xt, PyArrayObject *__pyx_v_y, PyArrayObject *__pyx_v_alpha, PyArrayObject *__pyx_v_g) {
  PyArrayObject *__pyx_v_XtY = 0;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_Xt;
  __Pyx_Buffer __pyx_pybuffer_Xt;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_XtX;
  __Pyx_Buffer __pyx_pybuffer_XtX;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_XtY;
  __Pyx_Buffer __pyx_pybuffer_XtY;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_alpha;
  __Pyx_Buffer __pyx_pybuffer_alpha;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_g;
  __Pyx_Buffer __pyx_pybuffer_g;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_y;
  __Pyx_Buffer __pyx_pybuffer_y;
  PyArrayObject *__pyx_r = NULL;
  __pyx_pybuffer_XtY.pybuffer.buf = NULL;
  __pyx_pybuffer_XtY.refcount = 0;
  __pyx_pybuffernd_XtY.data = NULL;
  __pyx_pybuffernd_XtY.rcbuffer = &__pyx_pybuffer_XtY;
  __pyx_pybuffer_XtX.pybuffer.buf = NULL;
  __pyx_pybuffer_XtX.refcount = 0;
  __pyx_pybuffernd_XtX.data = NULL;
  __pyx_pybuffernd_XtX.rcbuffer = &__pyx_pybuffer_XtX;
  __pyx_pybuffer_Xt.pybuffer.buf = NULL;
  __pyx_pybuffer_Xt.refcount = 0;
  __pyx_pybuffernd_Xt.data = NULL;
  __pyx_pybuffernd_Xt.rcbuffer = &__pyx_pybuffer_Xt;
  __pyx_pybuffer_y.pybuffer.buf = NULL;
  __pyx_pybuffer_y.refcount = 0;
  __pyx_pybuffernd_y.data = NULL;
  __pyx_pybuffernd_y.rcbuffer = &__pyx_pybuffer_y;
  __pyx_pybuffer_alpha.pybuffer.buf = NULL;
  __pyx_pybuffer_alpha.refcount = 0;
  __pyx_pybuffernd_alpha.data = NULL;
  __pyx_pybuffernd_alpha.rcbuffer = &__pyx_pybuffer_alpha;
  __pyx_pybuffer_g.pybuffer.buf = NULL;
  __pyx_pybuffer_g.refcount = 0;
  __pyx_pybuffernd_g.data = NULL;
  __pyx_pybuffernd_g.rcbuffer = &__pyx_pybuffer_g;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_XtX.rcbuffer->pybuffer, (PyObject*)__pyx_v_XtX, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 129, __pyx_L1_error)
  }
  __pyx_pybuffernd_XtX.diminfo[0].strides = __pyx_pybuffernd_XtX.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_XtX.diminfo[0].shape = __pyx_pybuffernd_XtX.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_XtX.diminfo[1].strides = __pyx_pybuffernd_XtX.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_XtX.diminfo[1].shape = __pyx_pybuffernd_XtX.rcbuffer->pybuffer.shape[1];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_Xt.rcbuffer->pybuffer, (PyObject*)__pyx_v_Xt, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 129, __pyx_L1_error)
  }
  __pyx_pybuffernd_Xt.diminfo[0].strides = __pyx_pybuffernd_Xt.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_Xt.diminfo[0].shape = __pyx_pybuffernd_Xt.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_Xt.diminfo[1].strides = __pyx_pybuffernd_Xt.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_Xt.diminfo[1].shape = __pyx_pybuffernd_Xt.rcbuffer->pybuffer.shape[1];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_y.rcbuffer->pybuffer, (PyObject*)__pyx_v_y, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 3, 0, __pyx_stack) == -1)) __PYX_ERR(0, 129, __pyx_L1_error)
  }
  __pyx_pybuffernd_y.diminfo[0].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_y.diminfo[0].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_y.diminfo[1].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_y.diminfo[1].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[1]; __pyx_pybuffernd_y.diminfo[2].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[2]; __pyx_pybuffernd_y.diminfo[2].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[2];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer, (PyObject*)__pyx_v_alpha, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 129, __pyx_L1_error)
  }
  __pyx_pybuffernd_alpha.diminfo[0].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_alpha.diminfo[0].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_alpha.diminfo[1].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_alpha.diminfo[1].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[1];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_g.rcbuffer->pybuffer, (PyObject*)__pyx_v_g, &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 129, __pyx_L1_error)
  }
  __pyx_pybuffernd_g.diminfo[0].strides = __pyx_pybuffernd_g.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_g.diminfo[0].shape = __pyx_pybuffernd_g.rcbuffer->pybuffer.shape[0];
/* … */
  /* 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);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_Xt.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_XtX.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_XtY.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_y.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("groupedpaneldatamodels.utils.bonhomme_manresa._compute_theta_fast", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_Xt.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_XtX.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_XtY.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_y.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_XtY);
  __Pyx_XGIVEREF((PyObject *)__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 130:         np.ndarray[DTYPE_t, ndim=2] XtX,

 131:         np.ndarray[DTYPE_t, ndim=2] Xt,

 132:         np.ndarray[DTYPE_t, ndim=3] y,

 133:         np.ndarray[DTYPE_t, ndim=2] alpha,

 134:         np.ndarray[np.intp_t, ndim=1] g):

 135:     """

 136:     Faster θ̂ update using precomputed X'X and X'.

 137:     """

+138:     cdef np.ndarray[DTYPE_t, ndim=2] XtY = _np.matmul(

  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 138, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_matmul); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 138, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
/* … */
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 138, __pyx_L1_error)
  __pyx_t_6 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_XtY.rcbuffer->pybuffer, (PyObject*)__pyx_t_6, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_XtY = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_XtY.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 138, __pyx_L1_error)
    } else {__pyx_pybuffernd_XtY.diminfo[0].strides = __pyx_pybuffernd_XtY.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_XtY.diminfo[0].shape = __pyx_pybuffernd_XtY.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_XtY.diminfo[1].strides = __pyx_pybuffernd_XtY.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_XtY.diminfo[1].shape = __pyx_pybuffernd_XtY.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_6 = 0;
  __pyx_v_XtY = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

+139:         Xt, (y - alpha[g][:, :, None]).reshape(-1, 1))

  __pyx_t_2 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_alpha), ((PyObject *)__pyx_v_g)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 139, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_4 = __Pyx_PyObject_GetItem(__pyx_t_2, __pyx_tuple__8); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 139, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = PyNumber_Subtract(((PyObject *)__pyx_v_y), __pyx_t_4); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 139, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_reshape); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 139, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_tuple__3, NULL); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 139, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = NULL;
  __pyx_t_5 = 0;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_3))) {
    __pyx_t_4 = PyMethod_GET_SELF(__pyx_t_3);
    if (likely(__pyx_t_4)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
      __Pyx_INCREF(__pyx_t_4);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_3, function);
      __pyx_t_5 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[3] = {__pyx_t_4, ((PyObject *)__pyx_v_Xt), __pyx_t_2};
    __pyx_t_1 = __Pyx_PyObject_FastCall(__pyx_t_3, __pyx_callargs+1-__pyx_t_5, 2+__pyx_t_5);
    __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 138, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  }

+140:     return _np.linalg.solve(XtX, XtY)

  __Pyx_XDECREF((PyObject *)__pyx_r);
  __Pyx_GetModuleGlobalName(__pyx_t_3, __pyx_n_s_np); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 140, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_linalg); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 140, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_solve); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 140, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = NULL;
  __pyx_t_5 = 0;
  #if CYTHON_UNPACK_METHODS
  if (likely(PyMethod_Check(__pyx_t_3))) {
    __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_3);
    if (likely(__pyx_t_2)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
      __Pyx_INCREF(__pyx_t_2);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_3, function);
      __pyx_t_5 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[3] = {__pyx_t_2, ((PyObject *)__pyx_v_XtX), ((PyObject *)__pyx_v_XtY)};
    __pyx_t_1 = __Pyx_PyObject_FastCall(__pyx_t_3, __pyx_callargs+1-__pyx_t_5, 2+__pyx_t_5);
    __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 140, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  }
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 140, __pyx_L1_error)
  __pyx_r = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;
  goto __pyx_L0;

 141: 

 142: 

+143: cdef inline np.ndarray[DTYPE_t, ndim=2] _compute_residuals(

static CYTHON_INLINE PyArrayObject *__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_residuals(PyArrayObject *__pyx_v_y, PyArrayObject *__pyx_v_x, PyArrayObject *__pyx_v_theta) {
  __Pyx_LocalBuf_ND __pyx_pybuffernd_theta;
  __Pyx_Buffer __pyx_pybuffer_theta;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_x;
  __Pyx_Buffer __pyx_pybuffer_x;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_y;
  __Pyx_Buffer __pyx_pybuffer_y;
  PyArrayObject *__pyx_r = NULL;
  __pyx_pybuffer_y.pybuffer.buf = NULL;
  __pyx_pybuffer_y.refcount = 0;
  __pyx_pybuffernd_y.data = NULL;
  __pyx_pybuffernd_y.rcbuffer = &__pyx_pybuffer_y;
  __pyx_pybuffer_x.pybuffer.buf = NULL;
  __pyx_pybuffer_x.refcount = 0;
  __pyx_pybuffernd_x.data = NULL;
  __pyx_pybuffernd_x.rcbuffer = &__pyx_pybuffer_x;
  __pyx_pybuffer_theta.pybuffer.buf = NULL;
  __pyx_pybuffer_theta.refcount = 0;
  __pyx_pybuffernd_theta.data = NULL;
  __pyx_pybuffernd_theta.rcbuffer = &__pyx_pybuffer_theta;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_y.rcbuffer->pybuffer, (PyObject*)__pyx_v_y, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 3, 0, __pyx_stack) == -1)) __PYX_ERR(0, 143, __pyx_L1_error)
  }
  __pyx_pybuffernd_y.diminfo[0].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_y.diminfo[0].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_y.diminfo[1].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_y.diminfo[1].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[1]; __pyx_pybuffernd_y.diminfo[2].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[2]; __pyx_pybuffernd_y.diminfo[2].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[2];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_x.rcbuffer->pybuffer, (PyObject*)__pyx_v_x, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 3, 0, __pyx_stack) == -1)) __PYX_ERR(0, 143, __pyx_L1_error)
  }
  __pyx_pybuffernd_x.diminfo[0].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_x.diminfo[0].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_x.diminfo[1].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_x.diminfo[1].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[1]; __pyx_pybuffernd_x.diminfo[2].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[2]; __pyx_pybuffernd_x.diminfo[2].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[2];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_theta.rcbuffer->pybuffer, (PyObject*)__pyx_v_theta, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 143, __pyx_L1_error)
  }
  __pyx_pybuffernd_theta.diminfo[0].strides = __pyx_pybuffernd_theta.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_theta.diminfo[0].shape = __pyx_pybuffernd_theta.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_theta.diminfo[1].strides = __pyx_pybuffernd_theta.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_theta.diminfo[1].shape = __pyx_pybuffernd_theta.rcbuffer->pybuffer.shape[1];
/* … */
  /* 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);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_theta.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_y.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("groupedpaneldatamodels.utils.bonhomme_manresa._compute_residuals", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_theta.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_y.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XGIVEREF((PyObject *)__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 144:         np.ndarray[DTYPE_t, ndim=3] y,

 145:         np.ndarray[DTYPE_t, ndim=3] x,

 146:         np.ndarray[DTYPE_t, ndim=2] theta):

 147:     """

 148:     r_i = y_i − x_i θ̂

 149:     """

+150:     return _np.squeeze(y - _np.matmul(x, theta))

  __Pyx_XDECREF((PyObject *)__pyx_r);
  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 150, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_squeeze); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 150, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __Pyx_GetModuleGlobalName(__pyx_t_4, __pyx_n_s_np); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 150, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_matmul); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 150, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;
  __pyx_t_4 = NULL;
  __pyx_t_6 = 0;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_5))) {
    __pyx_t_4 = PyMethod_GET_SELF(__pyx_t_5);
    if (likely(__pyx_t_4)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_5);
      __Pyx_INCREF(__pyx_t_4);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_5, function);
      __pyx_t_6 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[3] = {__pyx_t_4, ((PyObject *)__pyx_v_x), ((PyObject *)__pyx_v_theta)};
    __pyx_t_2 = __Pyx_PyObject_FastCall(__pyx_t_5, __pyx_callargs+1-__pyx_t_6, 2+__pyx_t_6);
    __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
    if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 150, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_2);
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  }
  __pyx_t_5 = PyNumber_Subtract(((PyObject *)__pyx_v_y), __pyx_t_2); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 150, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = NULL;
  __pyx_t_6 = 0;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_3))) {
    __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_3);
    if (likely(__pyx_t_2)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
      __Pyx_INCREF(__pyx_t_2);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_3, function);
      __pyx_t_6 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[2] = {__pyx_t_2, __pyx_t_5};
    __pyx_t_1 = __Pyx_PyObject_FastCall(__pyx_t_3, __pyx_callargs+1-__pyx_t_6, 1+__pyx_t_6);
    __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
    __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 150, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  }
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 150, __pyx_L1_error)
  __pyx_r = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;
  goto __pyx_L0;

 151: 

 152: 

+153: cdef inline double _compute_objective_value(

static CYTHON_INLINE double __pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_objective_value(PyArrayObject *__pyx_v_res, PyArrayObject *__pyx_v_alpha, PyArrayObject *__pyx_v_g) {
  __Pyx_LocalBuf_ND __pyx_pybuffernd_alpha;
  __Pyx_Buffer __pyx_pybuffer_alpha;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_g;
  __Pyx_Buffer __pyx_pybuffer_g;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_res;
  __Pyx_Buffer __pyx_pybuffer_res;
  double __pyx_r;
  __pyx_pybuffer_res.pybuffer.buf = NULL;
  __pyx_pybuffer_res.refcount = 0;
  __pyx_pybuffernd_res.data = NULL;
  __pyx_pybuffernd_res.rcbuffer = &__pyx_pybuffer_res;
  __pyx_pybuffer_alpha.pybuffer.buf = NULL;
  __pyx_pybuffer_alpha.refcount = 0;
  __pyx_pybuffernd_alpha.data = NULL;
  __pyx_pybuffernd_alpha.rcbuffer = &__pyx_pybuffer_alpha;
  __pyx_pybuffer_g.pybuffer.buf = NULL;
  __pyx_pybuffer_g.refcount = 0;
  __pyx_pybuffernd_g.data = NULL;
  __pyx_pybuffernd_g.rcbuffer = &__pyx_pybuffer_g;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_res.rcbuffer->pybuffer, (PyObject*)__pyx_v_res, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 153, __pyx_L1_error)
  }
  __pyx_pybuffernd_res.diminfo[0].strides = __pyx_pybuffernd_res.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_res.diminfo[0].shape = __pyx_pybuffernd_res.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_res.diminfo[1].strides = __pyx_pybuffernd_res.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_res.diminfo[1].shape = __pyx_pybuffernd_res.rcbuffer->pybuffer.shape[1];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer, (PyObject*)__pyx_v_alpha, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) __PYX_ERR(0, 153, __pyx_L1_error)
  }
  __pyx_pybuffernd_alpha.diminfo[0].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_alpha.diminfo[0].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_alpha.diminfo[1].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_alpha.diminfo[1].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[1];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_g.rcbuffer->pybuffer, (PyObject*)__pyx_v_g, &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) __PYX_ERR(0, 153, __pyx_L1_error)
  }
  __pyx_pybuffernd_g.diminfo[0].strides = __pyx_pybuffernd_g.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_g.diminfo[0].shape = __pyx_pybuffernd_g.rcbuffer->pybuffer.shape[0];
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_1);
  __Pyx_XDECREF(__pyx_t_2);
  __Pyx_XDECREF(__pyx_t_3);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_res.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("groupedpaneldatamodels.utils.bonhomme_manresa._compute_objective_value", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = -1;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_res.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 154:         np.ndarray[DTYPE_t, ndim=2] res,

 155:         np.ndarray[DTYPE_t, ndim=2] alpha,

 156:         np.ndarray[np.intp_t, ndim=1] g):

 157:     """

 158:     Sum of squared residuals.

 159:     """

+160:     return float(((res - alpha[g]) ** 2).sum())

  __pyx_t_2 = __Pyx_PyObject_GetItem(((PyObject *)__pyx_v_alpha), ((PyObject *)__pyx_v_g)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 160, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = PyNumber_Subtract(((PyObject *)__pyx_v_res), __pyx_t_2); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 160, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = PyNumber_Power(__pyx_t_3, __pyx_int_2, Py_None); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 160, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_sum); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 160, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = NULL;
  __pyx_t_4 = 0;
  #if CYTHON_UNPACK_METHODS
  if (likely(PyMethod_Check(__pyx_t_3))) {
    __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_3);
    if (likely(__pyx_t_2)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
      __Pyx_INCREF(__pyx_t_2);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_3, function);
      __pyx_t_4 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[2] = {__pyx_t_2, NULL};
    __pyx_t_1 = __Pyx_PyObject_FastCall(__pyx_t_3, __pyx_callargs+1-__pyx_t_4, 0+__pyx_t_4);
    __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 160, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  }
  __pyx_t_5 = __Pyx_PyObject_AsDouble(__pyx_t_1); if (unlikely(__pyx_t_5 == ((double)((double)-1)) && PyErr_Occurred())) __PYX_ERR(0, 160, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  __pyx_r = __pyx_t_5;
  goto __pyx_L0;

 161: 

 162: 

 163: # -----------------------------------------------------------------------------

 164: # One full inner cycle of the algorithm

 165: # -----------------------------------------------------------------------------

 166: 

+167: cdef tuple _gfe_iteration(

static PyObject *__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__gfe_iteration(PyArrayObject *__pyx_v_y, PyArrayObject *__pyx_v_x, int __pyx_v_G, int __pyx_v_N, int __pyx_v_K, int __pyx_v_max_iter, double __pyx_v_tol) {
  PyArrayObject *__pyx_v_theta = 0;
  PyArrayObject *__pyx_v_alpha = 0;
  PyArrayObject *__pyx_v_g = 0;
  PyArrayObject *__pyx_v_res = 0;
  double __pyx_v_old_obj;
  double __pyx_v_new_obj;
  int __pyx_v_it_used;
  PyArrayObject *__pyx_v_X_flat = 0;
  PyArrayObject *__pyx_v_Xt = 0;
  PyArrayObject *__pyx_v_XtX = 0;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_X_flat;
  __Pyx_Buffer __pyx_pybuffer_X_flat;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_Xt;
  __Pyx_Buffer __pyx_pybuffer_Xt;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_XtX;
  __Pyx_Buffer __pyx_pybuffer_XtX;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_alpha;
  __Pyx_Buffer __pyx_pybuffer_alpha;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_g;
  __Pyx_Buffer __pyx_pybuffer_g;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_res;
  __Pyx_Buffer __pyx_pybuffer_res;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_theta;
  __Pyx_Buffer __pyx_pybuffer_theta;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_x;
  __Pyx_Buffer __pyx_pybuffer_x;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_y;
  __Pyx_Buffer __pyx_pybuffer_y;
  PyObject *__pyx_r = NULL;
  __pyx_pybuffer_theta.pybuffer.buf = NULL;
  __pyx_pybuffer_theta.refcount = 0;
  __pyx_pybuffernd_theta.data = NULL;
  __pyx_pybuffernd_theta.rcbuffer = &__pyx_pybuffer_theta;
  __pyx_pybuffer_alpha.pybuffer.buf = NULL;
  __pyx_pybuffer_alpha.refcount = 0;
  __pyx_pybuffernd_alpha.data = NULL;
  __pyx_pybuffernd_alpha.rcbuffer = &__pyx_pybuffer_alpha;
  __pyx_pybuffer_g.pybuffer.buf = NULL;
  __pyx_pybuffer_g.refcount = 0;
  __pyx_pybuffernd_g.data = NULL;
  __pyx_pybuffernd_g.rcbuffer = &__pyx_pybuffer_g;
  __pyx_pybuffer_res.pybuffer.buf = NULL;
  __pyx_pybuffer_res.refcount = 0;
  __pyx_pybuffernd_res.data = NULL;
  __pyx_pybuffernd_res.rcbuffer = &__pyx_pybuffer_res;
  __pyx_pybuffer_X_flat.pybuffer.buf = NULL;
  __pyx_pybuffer_X_flat.refcount = 0;
  __pyx_pybuffernd_X_flat.data = NULL;
  __pyx_pybuffernd_X_flat.rcbuffer = &__pyx_pybuffer_X_flat;
  __pyx_pybuffer_Xt.pybuffer.buf = NULL;
  __pyx_pybuffer_Xt.refcount = 0;
  __pyx_pybuffernd_Xt.data = NULL;
  __pyx_pybuffernd_Xt.rcbuffer = &__pyx_pybuffer_Xt;
  __pyx_pybuffer_XtX.pybuffer.buf = NULL;
  __pyx_pybuffer_XtX.refcount = 0;
  __pyx_pybuffernd_XtX.data = NULL;
  __pyx_pybuffernd_XtX.rcbuffer = &__pyx_pybuffer_XtX;
  __pyx_pybuffer_y.pybuffer.buf = NULL;
  __pyx_pybuffer_y.refcount = 0;
  __pyx_pybuffernd_y.data = NULL;
  __pyx_pybuffernd_y.rcbuffer = &__pyx_pybuffer_y;
  __pyx_pybuffer_x.pybuffer.buf = NULL;
  __pyx_pybuffer_x.refcount = 0;
  __pyx_pybuffernd_x.data = NULL;
  __pyx_pybuffernd_x.rcbuffer = &__pyx_pybuffer_x;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_y.rcbuffer->pybuffer, (PyObject*)__pyx_v_y, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 3, 0, __pyx_stack) == -1)) __PYX_ERR(0, 167, __pyx_L1_error)
  }
  __pyx_pybuffernd_y.diminfo[0].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_y.diminfo[0].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_y.diminfo[1].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_y.diminfo[1].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[1]; __pyx_pybuffernd_y.diminfo[2].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[2]; __pyx_pybuffernd_y.diminfo[2].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[2];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_x.rcbuffer->pybuffer, (PyObject*)__pyx_v_x, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 3, 0, __pyx_stack) == -1)) __PYX_ERR(0, 167, __pyx_L1_error)
  }
  __pyx_pybuffernd_x.diminfo[0].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_x.diminfo[0].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_x.diminfo[1].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_x.diminfo[1].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[1]; __pyx_pybuffernd_x.diminfo[2].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[2]; __pyx_pybuffernd_x.diminfo[2].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[2];
/* … */
  /* 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);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_X_flat.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_Xt.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_XtX.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_res.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_theta.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_y.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("groupedpaneldatamodels.utils.bonhomme_manresa._gfe_iteration", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = 0;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_X_flat.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_Xt.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_XtX.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_res.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_theta.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_y.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_theta);
  __Pyx_XDECREF((PyObject *)__pyx_v_alpha);
  __Pyx_XDECREF((PyObject *)__pyx_v_g);
  __Pyx_XDECREF((PyObject *)__pyx_v_res);
  __Pyx_XDECREF((PyObject *)__pyx_v_X_flat);
  __Pyx_XDECREF((PyObject *)__pyx_v_Xt);
  __Pyx_XDECREF((PyObject *)__pyx_v_XtX);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

 168:         np.ndarray[DTYPE_t, ndim=3] y,

 169:         np.ndarray[DTYPE_t, ndim=3] x,

 170:         int G, int N, int K,

 171:         int max_iter, double tol):

 172:     """

 173:     Core fixed-point iterations (one random initialisation).

 174:     """

 175:     cdef np.ndarray[DTYPE_t, ndim=2] theta, alpha

 176:     cdef np.ndarray[np.intp_t, ndim=1] g

 177:     cdef np.ndarray[DTYPE_t, ndim=2] res

+178:     cdef double old_obj = INFINITY, new_obj

  __pyx_v_old_obj = INFINITY;

+179:     cdef int it_used = 0

  __pyx_v_it_used = 0;

 180: 

 181:     # Pre‑compute X'X and X' once: they depend only on x

+182:     cdef np.ndarray[DTYPE_t, ndim=2] X_flat = x.reshape(-1, K)

  __pyx_t_2 = __Pyx_PyObject_GetAttrStr(((PyObject *)__pyx_v_x), __pyx_n_s_reshape); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 182, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = __Pyx_PyInt_From_int(__pyx_v_K); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 182, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __pyx_t_4 = NULL;
  __pyx_t_5 = 0;
  #if CYTHON_UNPACK_METHODS
  if (likely(PyMethod_Check(__pyx_t_2))) {
    __pyx_t_4 = PyMethod_GET_SELF(__pyx_t_2);
    if (likely(__pyx_t_4)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_2);
      __Pyx_INCREF(__pyx_t_4);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_2, function);
      __pyx_t_5 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[3] = {__pyx_t_4, __pyx_int_neg_1, __pyx_t_3};
    __pyx_t_1 = __Pyx_PyObject_FastCall(__pyx_t_2, __pyx_callargs+1-__pyx_t_5, 2+__pyx_t_5);
    __Pyx_XDECREF(__pyx_t_4); __pyx_t_4 = 0;
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 182, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  }
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 182, __pyx_L1_error)
  __pyx_t_6 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_X_flat.rcbuffer->pybuffer, (PyObject*)__pyx_t_6, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_X_flat = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_X_flat.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 182, __pyx_L1_error)
    } else {__pyx_pybuffernd_X_flat.diminfo[0].strides = __pyx_pybuffernd_X_flat.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_X_flat.diminfo[0].shape = __pyx_pybuffernd_X_flat.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_X_flat.diminfo[1].strides = __pyx_pybuffernd_X_flat.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_X_flat.diminfo[1].shape = __pyx_pybuffernd_X_flat.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_6 = 0;
  __pyx_v_X_flat = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

+183:     cdef np.ndarray[DTYPE_t, ndim=2] Xt = X_flat.T

  __pyx_t_1 = __Pyx_PyObject_GetAttrStr(((PyObject *)__pyx_v_X_flat), __pyx_n_s_T); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 183, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 183, __pyx_L1_error)
  __pyx_t_7 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_Xt.rcbuffer->pybuffer, (PyObject*)__pyx_t_7, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_Xt = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_Xt.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 183, __pyx_L1_error)
    } else {__pyx_pybuffernd_Xt.diminfo[0].strides = __pyx_pybuffernd_Xt.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_Xt.diminfo[0].shape = __pyx_pybuffernd_Xt.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_Xt.diminfo[1].strides = __pyx_pybuffernd_Xt.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_Xt.diminfo[1].shape = __pyx_pybuffernd_Xt.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_7 = 0;
  __pyx_v_Xt = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

+184:     cdef np.ndarray[DTYPE_t, ndim=2] XtX = _np.matmul(Xt, X_flat)

  __Pyx_GetModuleGlobalName(__pyx_t_2, __pyx_n_s_np); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 184, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_matmul); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 184, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0;
  __pyx_t_2 = NULL;
  __pyx_t_5 = 0;
  #if CYTHON_UNPACK_METHODS
  if (unlikely(PyMethod_Check(__pyx_t_3))) {
    __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_3);
    if (likely(__pyx_t_2)) {
      PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3);
      __Pyx_INCREF(__pyx_t_2);
      __Pyx_INCREF(function);
      __Pyx_DECREF_SET(__pyx_t_3, function);
      __pyx_t_5 = 1;
    }
  }
  #endif
  {
    PyObject *__pyx_callargs[3] = {__pyx_t_2, ((PyObject *)__pyx_v_Xt), ((PyObject *)__pyx_v_X_flat)};
    __pyx_t_1 = __Pyx_PyObject_FastCall(__pyx_t_3, __pyx_callargs+1-__pyx_t_5, 2+__pyx_t_5);
    __Pyx_XDECREF(__pyx_t_2); __pyx_t_2 = 0;
    if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 184, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0;
  }
  if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 184, __pyx_L1_error)
  __pyx_t_8 = ((PyArrayObject *)__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_XtX.rcbuffer->pybuffer, (PyObject*)__pyx_t_8, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_XtX = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_XtX.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 184, __pyx_L1_error)
    } else {__pyx_pybuffernd_XtX.diminfo[0].strides = __pyx_pybuffernd_XtX.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_XtX.diminfo[0].shape = __pyx_pybuffernd_XtX.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_XtX.diminfo[1].strides = __pyx_pybuffernd_XtX.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_XtX.diminfo[1].shape = __pyx_pybuffernd_XtX.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_8 = 0;
  __pyx_v_XtX = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

 185: 

+186:     theta, alpha = _get_starting_values(y, x, G, N, K)

  __pyx_t_1 = __pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__get_starting_values(((PyArrayObject *)__pyx_v_y), ((PyArrayObject *)__pyx_v_x), __pyx_v_G, __pyx_v_N, __pyx_v_K); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 186, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  if (likely(__pyx_t_1 != Py_None)) {
    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, 186, __pyx_L1_error)
    }
    #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS
    __pyx_t_3 = PyTuple_GET_ITEM(sequence, 0); 
    __pyx_t_2 = PyTuple_GET_ITEM(sequence, 1); 
    __Pyx_INCREF(__pyx_t_3);
    __Pyx_INCREF(__pyx_t_2);
    #else
    __pyx_t_3 = PySequence_ITEM(sequence, 0); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 186, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_3);
    __pyx_t_2 = PySequence_ITEM(sequence, 1); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 186, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_2);
    #endif
    __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0;
  } else {
    __Pyx_RaiseNoneNotIterableError(); __PYX_ERR(0, 186, __pyx_L1_error)
  }
  if (!(likely(((__pyx_t_3) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_3, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 186, __pyx_L1_error)
  if (!(likely(((__pyx_t_2) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_2, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 186, __pyx_L1_error)
  __pyx_t_9 = ((PyArrayObject *)__pyx_t_3);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_theta.rcbuffer->pybuffer);
    __pyx_t_10 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_theta.rcbuffer->pybuffer, (PyObject*)__pyx_t_9, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
    if (unlikely(__pyx_t_10 < 0)) {
      PyErr_Fetch(&__pyx_t_11, &__pyx_t_12, &__pyx_t_13);
      if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_theta.rcbuffer->pybuffer, (PyObject*)__pyx_v_theta, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
        Py_XDECREF(__pyx_t_11); Py_XDECREF(__pyx_t_12); Py_XDECREF(__pyx_t_13);
        __Pyx_RaiseBufferFallbackError();
      } else {
        PyErr_Restore(__pyx_t_11, __pyx_t_12, __pyx_t_13);
      }
      __pyx_t_11 = __pyx_t_12 = __pyx_t_13 = 0;
    }
    __pyx_pybuffernd_theta.diminfo[0].strides = __pyx_pybuffernd_theta.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_theta.diminfo[0].shape = __pyx_pybuffernd_theta.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_theta.diminfo[1].strides = __pyx_pybuffernd_theta.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_theta.diminfo[1].shape = __pyx_pybuffernd_theta.rcbuffer->pybuffer.shape[1];
    if (unlikely((__pyx_t_10 < 0))) __PYX_ERR(0, 186, __pyx_L1_error)
  }
  __pyx_t_9 = 0;
  __pyx_v_theta = ((PyArrayObject *)__pyx_t_3);
  __pyx_t_3 = 0;
  __pyx_t_9 = ((PyArrayObject *)__pyx_t_2);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer);
    __pyx_t_10 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer, (PyObject*)__pyx_t_9, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
    if (unlikely(__pyx_t_10 < 0)) {
      PyErr_Fetch(&__pyx_t_13, &__pyx_t_12, &__pyx_t_11);
      if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer, (PyObject*)__pyx_v_alpha, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
        Py_XDECREF(__pyx_t_13); Py_XDECREF(__pyx_t_12); Py_XDECREF(__pyx_t_11);
        __Pyx_RaiseBufferFallbackError();
      } else {
        PyErr_Restore(__pyx_t_13, __pyx_t_12, __pyx_t_11);
      }
      __pyx_t_13 = __pyx_t_12 = __pyx_t_11 = 0;
    }
    __pyx_pybuffernd_alpha.diminfo[0].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_alpha.diminfo[0].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_alpha.diminfo[1].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_alpha.diminfo[1].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[1];
    if (unlikely((__pyx_t_10 < 0))) __PYX_ERR(0, 186, __pyx_L1_error)
  }
  __pyx_t_9 = 0;
  __pyx_v_alpha = ((PyArrayObject *)__pyx_t_2);
  __pyx_t_2 = 0;

+187:     res = _compute_residuals(y, x, theta)

  __pyx_t_1 = ((PyObject *)__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_residuals(((PyArrayObject *)__pyx_v_y), ((PyArrayObject *)__pyx_v_x), ((PyArrayObject *)__pyx_v_theta))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 187, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_res.rcbuffer->pybuffer);
    __pyx_t_10 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_res.rcbuffer->pybuffer, (PyObject*)((PyArrayObject *)__pyx_t_1), &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
    if (unlikely(__pyx_t_10 < 0)) {
      PyErr_Fetch(&__pyx_t_11, &__pyx_t_12, &__pyx_t_13);
      if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_res.rcbuffer->pybuffer, (PyObject*)__pyx_v_res, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
        Py_XDECREF(__pyx_t_11); Py_XDECREF(__pyx_t_12); Py_XDECREF(__pyx_t_13);
        __Pyx_RaiseBufferFallbackError();
      } else {
        PyErr_Restore(__pyx_t_11, __pyx_t_12, __pyx_t_13);
      }
      __pyx_t_11 = __pyx_t_12 = __pyx_t_13 = 0;
    }
    __pyx_pybuffernd_res.diminfo[0].strides = __pyx_pybuffernd_res.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_res.diminfo[0].shape = __pyx_pybuffernd_res.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_res.diminfo[1].strides = __pyx_pybuffernd_res.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_res.diminfo[1].shape = __pyx_pybuffernd_res.rcbuffer->pybuffer.shape[1];
    if (unlikely((__pyx_t_10 < 0))) __PYX_ERR(0, 187, __pyx_L1_error)
  }
  __pyx_v_res = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

+188:     g = _compute_groupings(res, alpha)

  __pyx_t_1 = ((PyObject *)__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_groupings(((PyArrayObject *)__pyx_v_res), ((PyArrayObject *)__pyx_v_alpha))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 188, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
    __pyx_t_10 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_g.rcbuffer->pybuffer, (PyObject*)((PyArrayObject *)__pyx_t_1), &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack);
    if (unlikely(__pyx_t_10 < 0)) {
      PyErr_Fetch(&__pyx_t_13, &__pyx_t_12, &__pyx_t_11);
      if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_g.rcbuffer->pybuffer, (PyObject*)__pyx_v_g, &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
        Py_XDECREF(__pyx_t_13); Py_XDECREF(__pyx_t_12); Py_XDECREF(__pyx_t_11);
        __Pyx_RaiseBufferFallbackError();
      } else {
        PyErr_Restore(__pyx_t_13, __pyx_t_12, __pyx_t_11);
      }
      __pyx_t_13 = __pyx_t_12 = __pyx_t_11 = 0;
    }
    __pyx_pybuffernd_g.diminfo[0].strides = __pyx_pybuffernd_g.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_g.diminfo[0].shape = __pyx_pybuffernd_g.rcbuffer->pybuffer.shape[0];
    if (unlikely((__pyx_t_10 < 0))) __PYX_ERR(0, 188, __pyx_L1_error)
  }
  __pyx_v_g = ((PyArrayObject *)__pyx_t_1);
  __pyx_t_1 = 0;

 189: 

+190:     for it_used in range(max_iter):

  __pyx_t_10 = __pyx_v_max_iter;
  __pyx_t_14 = __pyx_t_10;
  for (__pyx_t_15 = 0; __pyx_t_15 < __pyx_t_14; __pyx_t_15+=1) {
    __pyx_v_it_used = __pyx_t_15;

+191:         alpha = _compute_alpha_fast(res, g, G)

    __pyx_t_1 = ((PyObject *)__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_alpha_fast(((PyArrayObject *)__pyx_v_res), ((PyArrayObject *)__pyx_v_g), __pyx_v_G)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 191, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer);
      __pyx_t_16 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer, (PyObject*)((PyArrayObject *)__pyx_t_1), &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
      if (unlikely(__pyx_t_16 < 0)) {
        PyErr_Fetch(&__pyx_t_11, &__pyx_t_12, &__pyx_t_13);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer, (PyObject*)__pyx_v_alpha, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_11); Py_XDECREF(__pyx_t_12); Py_XDECREF(__pyx_t_13);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_11, __pyx_t_12, __pyx_t_13);
        }
        __pyx_t_11 = __pyx_t_12 = __pyx_t_13 = 0;
      }
      __pyx_pybuffernd_alpha.diminfo[0].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_alpha.diminfo[0].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_alpha.diminfo[1].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_alpha.diminfo[1].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[1];
      if (unlikely((__pyx_t_16 < 0))) __PYX_ERR(0, 191, __pyx_L1_error)
    }
    __Pyx_DECREF_SET(__pyx_v_alpha, ((PyArrayObject *)__pyx_t_1));
    __pyx_t_1 = 0;

+192:         theta = _compute_theta_fast(XtX, Xt, y, alpha, g)

    __pyx_t_1 = ((PyObject *)__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_theta_fast(((PyArrayObject *)__pyx_v_XtX), ((PyArrayObject *)__pyx_v_Xt), ((PyArrayObject *)__pyx_v_y), ((PyArrayObject *)__pyx_v_alpha), ((PyArrayObject *)__pyx_v_g))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 192, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_theta.rcbuffer->pybuffer);
      __pyx_t_16 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_theta.rcbuffer->pybuffer, (PyObject*)((PyArrayObject *)__pyx_t_1), &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
      if (unlikely(__pyx_t_16 < 0)) {
        PyErr_Fetch(&__pyx_t_13, &__pyx_t_12, &__pyx_t_11);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_theta.rcbuffer->pybuffer, (PyObject*)__pyx_v_theta, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_13); Py_XDECREF(__pyx_t_12); Py_XDECREF(__pyx_t_11);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_13, __pyx_t_12, __pyx_t_11);
        }
        __pyx_t_13 = __pyx_t_12 = __pyx_t_11 = 0;
      }
      __pyx_pybuffernd_theta.diminfo[0].strides = __pyx_pybuffernd_theta.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_theta.diminfo[0].shape = __pyx_pybuffernd_theta.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_theta.diminfo[1].strides = __pyx_pybuffernd_theta.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_theta.diminfo[1].shape = __pyx_pybuffernd_theta.rcbuffer->pybuffer.shape[1];
      if (unlikely((__pyx_t_16 < 0))) __PYX_ERR(0, 192, __pyx_L1_error)
    }
    __Pyx_DECREF_SET(__pyx_v_theta, ((PyArrayObject *)__pyx_t_1));
    __pyx_t_1 = 0;

+193:         res = _compute_residuals(y, x, theta)

    __pyx_t_1 = ((PyObject *)__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_residuals(((PyArrayObject *)__pyx_v_y), ((PyArrayObject *)__pyx_v_x), ((PyArrayObject *)__pyx_v_theta))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 193, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_res.rcbuffer->pybuffer);
      __pyx_t_16 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_res.rcbuffer->pybuffer, (PyObject*)((PyArrayObject *)__pyx_t_1), &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
      if (unlikely(__pyx_t_16 < 0)) {
        PyErr_Fetch(&__pyx_t_11, &__pyx_t_12, &__pyx_t_13);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_res.rcbuffer->pybuffer, (PyObject*)__pyx_v_res, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_11); Py_XDECREF(__pyx_t_12); Py_XDECREF(__pyx_t_13);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_11, __pyx_t_12, __pyx_t_13);
        }
        __pyx_t_11 = __pyx_t_12 = __pyx_t_13 = 0;
      }
      __pyx_pybuffernd_res.diminfo[0].strides = __pyx_pybuffernd_res.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_res.diminfo[0].shape = __pyx_pybuffernd_res.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_res.diminfo[1].strides = __pyx_pybuffernd_res.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_res.diminfo[1].shape = __pyx_pybuffernd_res.rcbuffer->pybuffer.shape[1];
      if (unlikely((__pyx_t_16 < 0))) __PYX_ERR(0, 193, __pyx_L1_error)
    }
    __Pyx_DECREF_SET(__pyx_v_res, ((PyArrayObject *)__pyx_t_1));
    __pyx_t_1 = 0;

+194:         alpha = _compute_alpha_fast(res, g, G)

    __pyx_t_1 = ((PyObject *)__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_alpha_fast(((PyArrayObject *)__pyx_v_res), ((PyArrayObject *)__pyx_v_g), __pyx_v_G)); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 194, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer);
      __pyx_t_16 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer, (PyObject*)((PyArrayObject *)__pyx_t_1), &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
      if (unlikely(__pyx_t_16 < 0)) {
        PyErr_Fetch(&__pyx_t_13, &__pyx_t_12, &__pyx_t_11);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_alpha.rcbuffer->pybuffer, (PyObject*)__pyx_v_alpha, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_13); Py_XDECREF(__pyx_t_12); Py_XDECREF(__pyx_t_11);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_13, __pyx_t_12, __pyx_t_11);
        }
        __pyx_t_13 = __pyx_t_12 = __pyx_t_11 = 0;
      }
      __pyx_pybuffernd_alpha.diminfo[0].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_alpha.diminfo[0].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_alpha.diminfo[1].strides = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_alpha.diminfo[1].shape = __pyx_pybuffernd_alpha.rcbuffer->pybuffer.shape[1];
      if (unlikely((__pyx_t_16 < 0))) __PYX_ERR(0, 194, __pyx_L1_error)
    }
    __Pyx_DECREF_SET(__pyx_v_alpha, ((PyArrayObject *)__pyx_t_1));
    __pyx_t_1 = 0;

+195:         g = _compute_groupings(res, alpha)

    __pyx_t_1 = ((PyObject *)__pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_groupings(((PyArrayObject *)__pyx_v_res), ((PyArrayObject *)__pyx_v_alpha))); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 195, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_1);
    {
      __Pyx_BufFmt_StackElem __pyx_stack[1];
      __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_g.rcbuffer->pybuffer);
      __pyx_t_16 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_g.rcbuffer->pybuffer, (PyObject*)((PyArrayObject *)__pyx_t_1), &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack);
      if (unlikely(__pyx_t_16 < 0)) {
        PyErr_Fetch(&__pyx_t_11, &__pyx_t_12, &__pyx_t_13);
        if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_g.rcbuffer->pybuffer, (PyObject*)__pyx_v_g, &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
          Py_XDECREF(__pyx_t_11); Py_XDECREF(__pyx_t_12); Py_XDECREF(__pyx_t_13);
          __Pyx_RaiseBufferFallbackError();
        } else {
          PyErr_Restore(__pyx_t_11, __pyx_t_12, __pyx_t_13);
        }
        __pyx_t_11 = __pyx_t_12 = __pyx_t_13 = 0;
      }
      __pyx_pybuffernd_g.diminfo[0].strides = __pyx_pybuffernd_g.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_g.diminfo[0].shape = __pyx_pybuffernd_g.rcbuffer->pybuffer.shape[0];
      if (unlikely((__pyx_t_16 < 0))) __PYX_ERR(0, 195, __pyx_L1_error)
    }
    __Pyx_DECREF_SET(__pyx_v_g, ((PyArrayObject *)__pyx_t_1));
    __pyx_t_1 = 0;

 196: 

+197:         new_obj = _compute_objective_value(res, alpha, g)

    __pyx_t_17 = __pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__compute_objective_value(((PyArrayObject *)__pyx_v_res), ((PyArrayObject *)__pyx_v_alpha), ((PyArrayObject *)__pyx_v_g)); if (unlikely(__pyx_t_17 == ((double)-1) && PyErr_Occurred())) __PYX_ERR(0, 197, __pyx_L1_error)
    __pyx_v_new_obj = __pyx_t_17;

+198:         if abs(old_obj - new_obj) < tol:

    __pyx_t_18 = (fabs((__pyx_v_old_obj - __pyx_v_new_obj)) < __pyx_v_tol);
    if (__pyx_t_18) {
/* … */
    }

+199:             break

      goto __pyx_L4_break;

+200:         old_obj = new_obj

    __pyx_v_old_obj = __pyx_v_new_obj;
  }
  __pyx_L4_break:;

 201: 

+202:     return theta, alpha, g, it_used, new_obj

  __Pyx_XDECREF(__pyx_r);
  __pyx_t_1 = __Pyx_PyInt_From_int(__pyx_v_it_used); if (unlikely(!__pyx_t_1)) __PYX_ERR(0, 202, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_1);
  __pyx_t_2 = PyFloat_FromDouble(__pyx_v_new_obj); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 202, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);
  __pyx_t_3 = PyTuple_New(5); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 202, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);
  __Pyx_INCREF((PyObject *)__pyx_v_theta);
  __Pyx_GIVEREF((PyObject *)__pyx_v_theta);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_3, 0, ((PyObject *)__pyx_v_theta))) __PYX_ERR(0, 202, __pyx_L1_error);
  __Pyx_INCREF((PyObject *)__pyx_v_alpha);
  __Pyx_GIVEREF((PyObject *)__pyx_v_alpha);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_3, 1, ((PyObject *)__pyx_v_alpha))) __PYX_ERR(0, 202, __pyx_L1_error);
  __Pyx_INCREF((PyObject *)__pyx_v_g);
  __Pyx_GIVEREF((PyObject *)__pyx_v_g);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_3, 2, ((PyObject *)__pyx_v_g))) __PYX_ERR(0, 202, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_1);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_3, 3, __pyx_t_1)) __PYX_ERR(0, 202, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_2);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_3, 4, __pyx_t_2)) __PYX_ERR(0, 202, __pyx_L1_error);
  __pyx_t_1 = 0;
  __pyx_t_2 = 0;
  __pyx_r = ((PyObject*)__pyx_t_3);
  __pyx_t_3 = 0;
  goto __pyx_L0;

 203: 

 204: 

 205: # -----------------------------------------------------------------------------

 206: # Public function – exactly the same signature you had before

 207: # -----------------------------------------------------------------------------

+208: def grouped_fixed_effects(np.ndarray[DTYPE_t, ndim=3] y,

/* Python wrapper */
static PyObject *__pyx_pw_22groupedpaneldatamodels_5utils_16bonhomme_manresa_1grouped_fixed_effects(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
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_22groupedpaneldatamodels_5utils_16bonhomme_manresa_grouped_fixed_effects, "\n    Run several random starts and return the best solution.\n    ");
static PyMethodDef __pyx_mdef_22groupedpaneldatamodels_5utils_16bonhomme_manresa_1grouped_fixed_effects = {"grouped_fixed_effects", (PyCFunction)(void*)(__Pyx_PyCFunction_FastCallWithKeywords)__pyx_pw_22groupedpaneldatamodels_5utils_16bonhomme_manresa_1grouped_fixed_effects, __Pyx_METH_FASTCALL|METH_KEYWORDS, __pyx_doc_22groupedpaneldatamodels_5utils_16bonhomme_manresa_grouped_fixed_effects};
static PyObject *__pyx_pw_22groupedpaneldatamodels_5utils_16bonhomme_manresa_1grouped_fixed_effects(PyObject *__pyx_self, 
#if CYTHON_METH_FASTCALL
PyObject *const *__pyx_args, Py_ssize_t __pyx_nargs, PyObject *__pyx_kwds
#else
PyObject *__pyx_args, PyObject *__pyx_kwds
#endif
) {
  PyArrayObject *__pyx_v_y = 0;
  PyArrayObject *__pyx_v_x = 0;
  int __pyx_v_G;
  int __pyx_v_max_iter;
  double __pyx_v_tol;
  int __pyx_v_gfe_iterations;
  #if !CYTHON_METH_FASTCALL
  CYTHON_UNUSED Py_ssize_t __pyx_nargs;
  #endif
  CYTHON_UNUSED PyObject *const *__pyx_kwvalues;
  PyObject *__pyx_r = 0;
  __Pyx_RefNannyDeclarations
  __Pyx_RefNannySetupContext("grouped_fixed_effects (wrapper)", 0);
  #if !CYTHON_METH_FASTCALL
  #if CYTHON_ASSUME_SAFE_MACROS
  __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 **__pyx_pyargnames[] = {&__pyx_n_s_y,&__pyx_n_s_x,&__pyx_n_s_G,&__pyx_n_s_max_iter,&__pyx_n_s_tol,&__pyx_n_s_gfe_iterations,0};
  PyObject* values[6] = {0,0,0,0,0,0};
    if (__pyx_kwds) {
      Py_ssize_t kw_args;
      switch (__pyx_nargs) {
        case  6: values[5] = __Pyx_Arg_FASTCALL(__pyx_args, 5);
        CYTHON_FALLTHROUGH;
        case  5: values[4] = __Pyx_Arg_FASTCALL(__pyx_args, 4);
        CYTHON_FALLTHROUGH;
        case  4: values[3] = __Pyx_Arg_FASTCALL(__pyx_args, 3);
        CYTHON_FALLTHROUGH;
        case  3: values[2] = __Pyx_Arg_FASTCALL(__pyx_args, 2);
        CYTHON_FALLTHROUGH;
        case  2: values[1] = __Pyx_Arg_FASTCALL(__pyx_args, 1);
        CYTHON_FALLTHROUGH;
        case  1: values[0] = __Pyx_Arg_FASTCALL(__pyx_args, 0);
        CYTHON_FALLTHROUGH;
        case  0: break;
        default: goto __pyx_L5_argtuple_error;
      }
      kw_args = __Pyx_NumKwargs_FASTCALL(__pyx_kwds);
      switch (__pyx_nargs) {
        case  0:
        if (likely((values[0] = __Pyx_GetKwValue_FASTCALL(__pyx_kwds, __pyx_kwvalues, __pyx_n_s_y)) != 0)) {
          (void)__Pyx_Arg_NewRef_FASTCALL(values[0]);
          kw_args--;
        }
        else if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 208, __pyx_L3_error)
        else goto __pyx_L5_argtuple_error;
        CYTHON_FALLTHROUGH;
        case  1:
        if (likely((values[1] = __Pyx_GetKwValue_FASTCALL(__pyx_kwds, __pyx_kwvalues, __pyx_n_s_x)) != 0)) {
          (void)__Pyx_Arg_NewRef_FASTCALL(values[1]);
          kw_args--;
        }
        else if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 208, __pyx_L3_error)
        else {
          __Pyx_RaiseArgtupleInvalid("grouped_fixed_effects", 0, 3, 6, 1); __PYX_ERR(0, 208, __pyx_L3_error)
        }
        CYTHON_FALLTHROUGH;
        case  2:
        if (likely((values[2] = __Pyx_GetKwValue_FASTCALL(__pyx_kwds, __pyx_kwvalues, __pyx_n_s_G)) != 0)) {
          (void)__Pyx_Arg_NewRef_FASTCALL(values[2]);
          kw_args--;
        }
        else if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 208, __pyx_L3_error)
        else {
          __Pyx_RaiseArgtupleInvalid("grouped_fixed_effects", 0, 3, 6, 2); __PYX_ERR(0, 208, __pyx_L3_error)
        }
        CYTHON_FALLTHROUGH;
        case  3:
        if (kw_args > 0) {
          PyObject* value = __Pyx_GetKwValue_FASTCALL(__pyx_kwds, __pyx_kwvalues, __pyx_n_s_max_iter);
          if (value) { values[3] = __Pyx_Arg_NewRef_FASTCALL(value); kw_args--; }
          else if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 208, __pyx_L3_error)
        }
        CYTHON_FALLTHROUGH;
        case  4:
        if (kw_args > 0) {
          PyObject* value = __Pyx_GetKwValue_FASTCALL(__pyx_kwds, __pyx_kwvalues, __pyx_n_s_tol);
          if (value) { values[4] = __Pyx_Arg_NewRef_FASTCALL(value); kw_args--; }
          else if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 208, __pyx_L3_error)
        }
        CYTHON_FALLTHROUGH;
        case  5:
        if (kw_args > 0) {
          PyObject* value = __Pyx_GetKwValue_FASTCALL(__pyx_kwds, __pyx_kwvalues, __pyx_n_s_gfe_iterations);
          if (value) { values[5] = __Pyx_Arg_NewRef_FASTCALL(value); kw_args--; }
          else if (unlikely(PyErr_Occurred())) __PYX_ERR(0, 208, __pyx_L3_error)
        }
      }
      if (unlikely(kw_args > 0)) {
        const Py_ssize_t kwd_pos_args = __pyx_nargs;
        if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_kwvalues, __pyx_pyargnames, 0, values + 0, kwd_pos_args, "grouped_fixed_effects") < 0)) __PYX_ERR(0, 208, __pyx_L3_error)
      }
    } else {
      switch (__pyx_nargs) {
        case  6: values[5] = __Pyx_Arg_FASTCALL(__pyx_args, 5);
        CYTHON_FALLTHROUGH;
        case  5: values[4] = __Pyx_Arg_FASTCALL(__pyx_args, 4);
        CYTHON_FALLTHROUGH;
        case  4: values[3] = __Pyx_Arg_FASTCALL(__pyx_args, 3);
        CYTHON_FALLTHROUGH;
        case  3: values[2] = __Pyx_Arg_FASTCALL(__pyx_args, 2);
        values[1] = __Pyx_Arg_FASTCALL(__pyx_args, 1);
        values[0] = __Pyx_Arg_FASTCALL(__pyx_args, 0);
        break;
        default: goto __pyx_L5_argtuple_error;
      }
    }
    __pyx_v_y = ((PyArrayObject *)values[0]);
    __pyx_v_x = ((PyArrayObject *)values[1]);
    __pyx_v_G = __Pyx_PyInt_As_int(values[2]); if (unlikely((__pyx_v_G == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 210, __pyx_L3_error)
    if (values[3]) {
      __pyx_v_max_iter = __Pyx_PyInt_As_int(values[3]); if (unlikely((__pyx_v_max_iter == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 211, __pyx_L3_error)
    } else {
      __pyx_v_max_iter = ((int)((int)0x2710));
    }
    if (values[4]) {
      __pyx_v_tol = __pyx_PyFloat_AsDouble(values[4]); if (unlikely((__pyx_v_tol == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 212, __pyx_L3_error)
    } else {
      __pyx_v_tol = ((double)((double)1e-8));
    }
    if (values[5]) {
      __pyx_v_gfe_iterations = __Pyx_PyInt_As_int(values[5]); if (unlikely((__pyx_v_gfe_iterations == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 213, __pyx_L3_error)
    } else {
      __pyx_v_gfe_iterations = ((int)((int)20));
    }
  }
  goto __pyx_L6_skip;
  __pyx_L5_argtuple_error:;
  __Pyx_RaiseArgtupleInvalid("grouped_fixed_effects", 0, 3, 6, __pyx_nargs); __PYX_ERR(0, 208, __pyx_L3_error)
  __pyx_L6_skip:;
  goto __pyx_L4_argument_unpacking_done;
  __pyx_L3_error:;
  {
    Py_ssize_t __pyx_temp;
    for (__pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
      __Pyx_Arg_XDECREF_FASTCALL(values[__pyx_temp]);
    }
  }
  __Pyx_AddTraceback("groupedpaneldatamodels.utils.bonhomme_manresa.grouped_fixed_effects", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __Pyx_RefNannyFinishContext();
  return NULL;
  __pyx_L4_argument_unpacking_done:;
  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_y), __pyx_ptype_5numpy_ndarray, 1, "y", 0))) __PYX_ERR(0, 208, __pyx_L1_error)
  if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_x), __pyx_ptype_5numpy_ndarray, 1, "x", 0))) __PYX_ERR(0, 209, __pyx_L1_error)
  __pyx_r = __pyx_pf_22groupedpaneldatamodels_5utils_16bonhomme_manresa_grouped_fixed_effects(__pyx_self, __pyx_v_y, __pyx_v_x, __pyx_v_G, __pyx_v_max_iter, __pyx_v_tol, __pyx_v_gfe_iterations);
  int __pyx_lineno = 0;
  const char *__pyx_filename = NULL;
  int __pyx_clineno = 0;

  /* function exit code */
  goto __pyx_L0;
  __pyx_L1_error:;
  __pyx_r = NULL;
  __pyx_L0:;
  {
    Py_ssize_t __pyx_temp;
    for (__pyx_temp=0; __pyx_temp < (Py_ssize_t)(sizeof(values)/sizeof(values[0])); ++__pyx_temp) {
      __Pyx_Arg_XDECREF_FASTCALL(values[__pyx_temp]);
    }
  }
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}

static PyObject *__pyx_pf_22groupedpaneldatamodels_5utils_16bonhomme_manresa_grouped_fixed_effects(CYTHON_UNUSED PyObject *__pyx_self, PyArrayObject *__pyx_v_y, PyArrayObject *__pyx_v_x, int __pyx_v_G, int __pyx_v_max_iter, double __pyx_v_tol, int __pyx_v_gfe_iterations) {
  int __pyx_v_N;
  int __pyx_v_K;
  PyArrayObject *__pyx_v_best_theta = 0;
  PyArrayObject *__pyx_v_best_alpha = 0;
  PyArrayObject *__pyx_v_best_g = 0;
  int __pyx_v_best_it;
  double __pyx_v_best_obj;
  PyObject *__pyx_v_res_tup = 0;
  CYTHON_UNUSED int __pyx_v__;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_best_alpha;
  __Pyx_Buffer __pyx_pybuffer_best_alpha;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_best_g;
  __Pyx_Buffer __pyx_pybuffer_best_g;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_best_theta;
  __Pyx_Buffer __pyx_pybuffer_best_theta;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_x;
  __Pyx_Buffer __pyx_pybuffer_x;
  __Pyx_LocalBuf_ND __pyx_pybuffernd_y;
  __Pyx_Buffer __pyx_pybuffer_y;
  PyObject *__pyx_r = NULL;
  __pyx_pybuffer_best_theta.pybuffer.buf = NULL;
  __pyx_pybuffer_best_theta.refcount = 0;
  __pyx_pybuffernd_best_theta.data = NULL;
  __pyx_pybuffernd_best_theta.rcbuffer = &__pyx_pybuffer_best_theta;
  __pyx_pybuffer_best_alpha.pybuffer.buf = NULL;
  __pyx_pybuffer_best_alpha.refcount = 0;
  __pyx_pybuffernd_best_alpha.data = NULL;
  __pyx_pybuffernd_best_alpha.rcbuffer = &__pyx_pybuffer_best_alpha;
  __pyx_pybuffer_best_g.pybuffer.buf = NULL;
  __pyx_pybuffer_best_g.refcount = 0;
  __pyx_pybuffernd_best_g.data = NULL;
  __pyx_pybuffernd_best_g.rcbuffer = &__pyx_pybuffer_best_g;
  __pyx_pybuffer_y.pybuffer.buf = NULL;
  __pyx_pybuffer_y.refcount = 0;
  __pyx_pybuffernd_y.data = NULL;
  __pyx_pybuffernd_y.rcbuffer = &__pyx_pybuffer_y;
  __pyx_pybuffer_x.pybuffer.buf = NULL;
  __pyx_pybuffer_x.refcount = 0;
  __pyx_pybuffernd_x.data = NULL;
  __pyx_pybuffernd_x.rcbuffer = &__pyx_pybuffer_x;
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_y.rcbuffer->pybuffer, (PyObject*)__pyx_v_y, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 3, 0, __pyx_stack) == -1)) __PYX_ERR(0, 208, __pyx_L1_error)
  }
  __pyx_pybuffernd_y.diminfo[0].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_y.diminfo[0].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_y.diminfo[1].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_y.diminfo[1].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[1]; __pyx_pybuffernd_y.diminfo[2].strides = __pyx_pybuffernd_y.rcbuffer->pybuffer.strides[2]; __pyx_pybuffernd_y.diminfo[2].shape = __pyx_pybuffernd_y.rcbuffer->pybuffer.shape[2];
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_x.rcbuffer->pybuffer, (PyObject*)__pyx_v_x, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 3, 0, __pyx_stack) == -1)) __PYX_ERR(0, 208, __pyx_L1_error)
  }
  __pyx_pybuffernd_x.diminfo[0].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_x.diminfo[0].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_x.diminfo[1].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_x.diminfo[1].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[1]; __pyx_pybuffernd_x.diminfo[2].strides = __pyx_pybuffernd_x.rcbuffer->pybuffer.strides[2]; __pyx_pybuffernd_x.diminfo[2].shape = __pyx_pybuffernd_x.rcbuffer->pybuffer.shape[2];
/* … */
  /* function exit code */
  __pyx_L1_error:;
  __Pyx_XDECREF(__pyx_t_7);
  __Pyx_XDECREF(__pyx_t_8);
  __Pyx_XDECREF(__pyx_t_10);
  __Pyx_XDECREF(__pyx_t_11);
  __Pyx_XDECREF(__pyx_t_12);
  { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
    __Pyx_PyThreadState_declare
    __Pyx_PyThreadState_assign
    __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_best_alpha.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_best_g.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_best_theta.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x.rcbuffer->pybuffer);
    __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_y.rcbuffer->pybuffer);
  __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
  __Pyx_AddTraceback("groupedpaneldatamodels.utils.bonhomme_manresa.grouped_fixed_effects", __pyx_clineno, __pyx_lineno, __pyx_filename);
  __pyx_r = NULL;
  goto __pyx_L2;
  __pyx_L0:;
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_best_alpha.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_best_g.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_best_theta.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_x.rcbuffer->pybuffer);
  __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_y.rcbuffer->pybuffer);
  __pyx_L2:;
  __Pyx_XDECREF((PyObject *)__pyx_v_best_theta);
  __Pyx_XDECREF((PyObject *)__pyx_v_best_alpha);
  __Pyx_XDECREF((PyObject *)__pyx_v_best_g);
  __Pyx_XDECREF(__pyx_v_res_tup);
  __Pyx_XGIVEREF(__pyx_r);
  __Pyx_RefNannyFinishContext();
  return __pyx_r;
}
/* … */
  __pyx_tuple__11 = PyTuple_Pack(15, __pyx_n_s_y, __pyx_n_s_x, __pyx_n_s_G, __pyx_n_s_max_iter, __pyx_n_s_tol, __pyx_n_s_gfe_iterations, __pyx_n_s_N, __pyx_n_s_K, __pyx_n_s_best_theta, __pyx_n_s_best_alpha, __pyx_n_s_best_g, __pyx_n_s_best_it, __pyx_n_s_best_obj, __pyx_n_s_res_tup, __pyx_n_s__10); if (unlikely(!__pyx_tuple__11)) __PYX_ERR(0, 208, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_tuple__11);
  __Pyx_GIVEREF(__pyx_tuple__11);
/* … */
  __pyx_t_5 = PyTuple_New(3); if (unlikely(!__pyx_t_5)) __PYX_ERR(0, 208, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_5);
  __Pyx_GIVEREF(__pyx_t_2);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_2)) __PYX_ERR(0, 208, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_3);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_5, 1, __pyx_t_3)) __PYX_ERR(0, 208, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_4);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_5, 2, __pyx_t_4)) __PYX_ERR(0, 208, __pyx_L1_error);
  __pyx_t_2 = 0;
  __pyx_t_3 = 0;
  __pyx_t_4 = 0;
  __pyx_t_4 = __Pyx_CyFunction_New(&__pyx_mdef_22groupedpaneldatamodels_5utils_16bonhomme_manresa_1grouped_fixed_effects, 0, __pyx_n_s_grouped_fixed_effects, NULL, __pyx_n_s_groupedpaneldatamodels_utils_bon, __pyx_d, ((PyObject *)__pyx_codeobj__12)); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 208, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);
  __Pyx_CyFunction_SetDefaultsTuple(__pyx_t_4, __pyx_t_5);
  __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0;
  if (PyDict_SetItem(__pyx_d, __pyx_n_s_grouped_fixed_effects, __pyx_t_4) < 0) __PYX_ERR(0, 208, __pyx_L1_error)
  __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0;

 209:                           np.ndarray[DTYPE_t, ndim=3] x,

 210:                           int G,

+211:                           int max_iter=10000,

  __pyx_t_2 = __Pyx_PyInt_From_int(((int)0x2710)); if (unlikely(!__pyx_t_2)) __PYX_ERR(0, 211, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_2);

+212:                           double tol=1e-8,

  __pyx_t_3 = PyFloat_FromDouble(((double)1e-8)); if (unlikely(!__pyx_t_3)) __PYX_ERR(0, 212, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_3);

+213:                           int gfe_iterations=20):

  __pyx_t_4 = __Pyx_PyInt_From_int(((int)20)); if (unlikely(!__pyx_t_4)) __PYX_ERR(0, 213, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_4);

 214:     """

 215:     Run several random starts and return the best solution.

 216:     """

+217:     cdef int N = x.shape[0]

  __pyx_v_N = (__pyx_f_5numpy_7ndarray_5shape_shape(((PyArrayObject *)__pyx_v_x))[0]);

+218:     cdef int K = x.shape[2]

  __pyx_v_K = (__pyx_f_5numpy_7ndarray_5shape_shape(((PyArrayObject *)__pyx_v_x))[2]);

 219: 

+220:     cdef np.ndarray[DTYPE_t, ndim=2] best_theta = None

  __pyx_t_1 = ((PyArrayObject *)Py_None);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_best_theta.rcbuffer->pybuffer, (PyObject*)__pyx_t_1, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_best_theta = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_best_theta.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 220, __pyx_L1_error)
    } else {__pyx_pybuffernd_best_theta.diminfo[0].strides = __pyx_pybuffernd_best_theta.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_best_theta.diminfo[0].shape = __pyx_pybuffernd_best_theta.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_best_theta.diminfo[1].strides = __pyx_pybuffernd_best_theta.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_best_theta.diminfo[1].shape = __pyx_pybuffernd_best_theta.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_1 = 0;
  __Pyx_INCREF(Py_None);
  __pyx_v_best_theta = ((PyArrayObject *)Py_None);

+221:     cdef np.ndarray[DTYPE_t, ndim=2] best_alpha = None

  __pyx_t_2 = ((PyArrayObject *)Py_None);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_best_alpha.rcbuffer->pybuffer, (PyObject*)__pyx_t_2, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
      __pyx_v_best_alpha = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_best_alpha.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 221, __pyx_L1_error)
    } else {__pyx_pybuffernd_best_alpha.diminfo[0].strides = __pyx_pybuffernd_best_alpha.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_best_alpha.diminfo[0].shape = __pyx_pybuffernd_best_alpha.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_best_alpha.diminfo[1].strides = __pyx_pybuffernd_best_alpha.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_best_alpha.diminfo[1].shape = __pyx_pybuffernd_best_alpha.rcbuffer->pybuffer.shape[1];
    }
  }
  __pyx_t_2 = 0;
  __Pyx_INCREF(Py_None);
  __pyx_v_best_alpha = ((PyArrayObject *)Py_None);

+222:     cdef np.ndarray[np.intp_t, ndim=1] best_g = None

  __pyx_t_3 = ((PyArrayObject *)Py_None);
  {
    __Pyx_BufFmt_StackElem __pyx_stack[1];
    if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_best_g.rcbuffer->pybuffer, (PyObject*)__pyx_t_3, &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
      __pyx_v_best_g = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_best_g.rcbuffer->pybuffer.buf = NULL;
      __PYX_ERR(0, 222, __pyx_L1_error)
    } else {__pyx_pybuffernd_best_g.diminfo[0].strides = __pyx_pybuffernd_best_g.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_best_g.diminfo[0].shape = __pyx_pybuffernd_best_g.rcbuffer->pybuffer.shape[0];
    }
  }
  __pyx_t_3 = 0;
  __Pyx_INCREF(Py_None);
  __pyx_v_best_g = ((PyArrayObject *)Py_None);

+223:     cdef int best_it = -1

  __pyx_v_best_it = -1;

+224:     cdef double best_obj = INFINITY

  __pyx_v_best_obj = INFINITY;

 225: 

 226:     cdef tuple res_tup

+227:     for _ in range(gfe_iterations):

  __pyx_t_4 = __pyx_v_gfe_iterations;
  __pyx_t_5 = __pyx_t_4;
  for (__pyx_t_6 = 0; __pyx_t_6 < __pyx_t_5; __pyx_t_6+=1) {
    __pyx_v__ = __pyx_t_6;

+228:         res_tup = _gfe_iteration(y, x, G, N, K, max_iter, tol)

    __pyx_t_7 = __pyx_f_22groupedpaneldatamodels_5utils_16bonhomme_manresa__gfe_iteration(((PyArrayObject *)__pyx_v_y), ((PyArrayObject *)__pyx_v_x), __pyx_v_G, __pyx_v_N, __pyx_v_K, __pyx_v_max_iter, __pyx_v_tol); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 228, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_7);
    __Pyx_XDECREF_SET(__pyx_v_res_tup, ((PyObject*)__pyx_t_7));
    __pyx_t_7 = 0;

+229:         if res_tup[4] < best_obj:

    if (unlikely(__pyx_v_res_tup == Py_None)) {
      PyErr_SetString(PyExc_TypeError, "'NoneType' object is not subscriptable");
      __PYX_ERR(0, 229, __pyx_L1_error)
    }
    __pyx_t_7 = PyFloat_FromDouble(__pyx_v_best_obj); if (unlikely(!__pyx_t_7)) __PYX_ERR(0, 229, __pyx_L1_error)
    __Pyx_GOTREF(__pyx_t_7);
    __pyx_t_8 = PyObject_RichCompare(PyTuple_GET_ITEM(__pyx_v_res_tup, 4), __pyx_t_7, Py_LT); __Pyx_XGOTREF(__pyx_t_8); if (unlikely(!__pyx_t_8)) __PYX_ERR(0, 229, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_7); __pyx_t_7 = 0;
    __pyx_t_9 = __Pyx_PyObject_IsTrue(__pyx_t_8); if (unlikely((__pyx_t_9 < 0))) __PYX_ERR(0, 229, __pyx_L1_error)
    __Pyx_DECREF(__pyx_t_8); __pyx_t_8 = 0;
    if (__pyx_t_9) {
/* … */
    }
  }

+230:             best_theta, best_alpha, best_g, best_it, best_obj = res_tup

      if (likely(__pyx_v_res_tup != Py_None)) {
        PyObject* sequence = __pyx_v_res_tup;
        Py_ssize_t size = __Pyx_PySequence_SIZE(sequence);
        if (unlikely(size != 5)) {
          if (size > 5) __Pyx_RaiseTooManyValuesError(5);
          else if (size >= 0) __Pyx_RaiseNeedMoreValuesError(size);
          __PYX_ERR(0, 230, __pyx_L1_error)
        }
        #if CYTHON_ASSUME_SAFE_MACROS && !CYTHON_AVOID_BORROWED_REFS
        __pyx_t_8 = PyTuple_GET_ITEM(sequence, 0); 
        __pyx_t_7 = PyTuple_GET_ITEM(sequence, 1); 
        __pyx_t_10 = PyTuple_GET_ITEM(sequence, 2); 
        __pyx_t_11 = PyTuple_GET_ITEM(sequence, 3); 
        __pyx_t_12 = PyTuple_GET_ITEM(sequence, 4); 
        __Pyx_INCREF(__pyx_t_8);
        __Pyx_INCREF(__pyx_t_7);
        __Pyx_INCREF(__pyx_t_10);
        __Pyx_INCREF(__pyx_t_11);
        __Pyx_INCREF(__pyx_t_12);
        #else
        {
          Py_ssize_t i;
          PyObject** temps[5] = {&__pyx_t_8,&__pyx_t_7,&__pyx_t_10,&__pyx_t_11,&__pyx_t_12};
          for (i=0; i < 5; i++) {
            PyObject* item = PySequence_ITEM(sequence, i); if (unlikely(!item)) __PYX_ERR(0, 230, __pyx_L1_error)
            __Pyx_GOTREF(item);
            *(temps[i]) = item;
          }
        }
        #endif
      } else {
        __Pyx_RaiseNoneNotIterableError(); __PYX_ERR(0, 230, __pyx_L1_error)
      }
      if (!(likely(((__pyx_t_8) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_8, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 230, __pyx_L1_error)
      if (!(likely(((__pyx_t_7) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_7, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 230, __pyx_L1_error)
      if (!(likely(((__pyx_t_10) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_10, __pyx_ptype_5numpy_ndarray))))) __PYX_ERR(0, 230, __pyx_L1_error)
      __pyx_t_13 = __Pyx_PyInt_As_int(__pyx_t_11); if (unlikely((__pyx_t_13 == (int)-1) && PyErr_Occurred())) __PYX_ERR(0, 230, __pyx_L1_error)
      __Pyx_DECREF(__pyx_t_11); __pyx_t_11 = 0;
      __pyx_t_14 = __pyx_PyFloat_AsDouble(__pyx_t_12); if (unlikely((__pyx_t_14 == (double)-1) && PyErr_Occurred())) __PYX_ERR(0, 230, __pyx_L1_error)
      __Pyx_DECREF(__pyx_t_12); __pyx_t_12 = 0;
      __pyx_t_1 = ((PyArrayObject *)__pyx_t_8);
      {
        __Pyx_BufFmt_StackElem __pyx_stack[1];
        __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_best_theta.rcbuffer->pybuffer);
        __pyx_t_15 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_best_theta.rcbuffer->pybuffer, (PyObject*)__pyx_t_1, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
        if (unlikely(__pyx_t_15 < 0)) {
          PyErr_Fetch(&__pyx_t_16, &__pyx_t_17, &__pyx_t_18);
          if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_best_theta.rcbuffer->pybuffer, (PyObject*)__pyx_v_best_theta, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
            Py_XDECREF(__pyx_t_16); Py_XDECREF(__pyx_t_17); Py_XDECREF(__pyx_t_18);
            __Pyx_RaiseBufferFallbackError();
          } else {
            PyErr_Restore(__pyx_t_16, __pyx_t_17, __pyx_t_18);
          }
          __pyx_t_16 = __pyx_t_17 = __pyx_t_18 = 0;
        }
        __pyx_pybuffernd_best_theta.diminfo[0].strides = __pyx_pybuffernd_best_theta.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_best_theta.diminfo[0].shape = __pyx_pybuffernd_best_theta.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_best_theta.diminfo[1].strides = __pyx_pybuffernd_best_theta.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_best_theta.diminfo[1].shape = __pyx_pybuffernd_best_theta.rcbuffer->pybuffer.shape[1];
        if (unlikely((__pyx_t_15 < 0))) __PYX_ERR(0, 230, __pyx_L1_error)
      }
      __pyx_t_1 = 0;
      __Pyx_DECREF_SET(__pyx_v_best_theta, ((PyArrayObject *)__pyx_t_8));
      __pyx_t_8 = 0;
      __pyx_t_2 = ((PyArrayObject *)__pyx_t_7);
      {
        __Pyx_BufFmt_StackElem __pyx_stack[1];
        __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_best_alpha.rcbuffer->pybuffer);
        __pyx_t_15 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_best_alpha.rcbuffer->pybuffer, (PyObject*)__pyx_t_2, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack);
        if (unlikely(__pyx_t_15 < 0)) {
          PyErr_Fetch(&__pyx_t_18, &__pyx_t_17, &__pyx_t_16);
          if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_best_alpha.rcbuffer->pybuffer, (PyObject*)__pyx_v_best_alpha, &__Pyx_TypeInfo_nn___pyx_t_22groupedpaneldatamodels_5utils_16bonhomme_manresa_DTYPE_t, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {
            Py_XDECREF(__pyx_t_18); Py_XDECREF(__pyx_t_17); Py_XDECREF(__pyx_t_16);
            __Pyx_RaiseBufferFallbackError();
          } else {
            PyErr_Restore(__pyx_t_18, __pyx_t_17, __pyx_t_16);
          }
          __pyx_t_18 = __pyx_t_17 = __pyx_t_16 = 0;
        }
        __pyx_pybuffernd_best_alpha.diminfo[0].strides = __pyx_pybuffernd_best_alpha.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_best_alpha.diminfo[0].shape = __pyx_pybuffernd_best_alpha.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_best_alpha.diminfo[1].strides = __pyx_pybuffernd_best_alpha.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_best_alpha.diminfo[1].shape = __pyx_pybuffernd_best_alpha.rcbuffer->pybuffer.shape[1];
        if (unlikely((__pyx_t_15 < 0))) __PYX_ERR(0, 230, __pyx_L1_error)
      }
      __pyx_t_2 = 0;
      __Pyx_DECREF_SET(__pyx_v_best_alpha, ((PyArrayObject *)__pyx_t_7));
      __pyx_t_7 = 0;
      __pyx_t_3 = ((PyArrayObject *)__pyx_t_10);
      {
        __Pyx_BufFmt_StackElem __pyx_stack[1];
        __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_best_g.rcbuffer->pybuffer);
        __pyx_t_15 = __Pyx_GetBufferAndValidate(&__pyx_pybuffernd_best_g.rcbuffer->pybuffer, (PyObject*)__pyx_t_3, &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack);
        if (unlikely(__pyx_t_15 < 0)) {
          PyErr_Fetch(&__pyx_t_16, &__pyx_t_17, &__pyx_t_18);
          if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_best_g.rcbuffer->pybuffer, (PyObject*)__pyx_v_best_g, &__Pyx_TypeInfo_nn___pyx_t_5numpy_intp_t, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {
            Py_XDECREF(__pyx_t_16); Py_XDECREF(__pyx_t_17); Py_XDECREF(__pyx_t_18);
            __Pyx_RaiseBufferFallbackError();
          } else {
            PyErr_Restore(__pyx_t_16, __pyx_t_17, __pyx_t_18);
          }
          __pyx_t_16 = __pyx_t_17 = __pyx_t_18 = 0;
        }
        __pyx_pybuffernd_best_g.diminfo[0].strides = __pyx_pybuffernd_best_g.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_best_g.diminfo[0].shape = __pyx_pybuffernd_best_g.rcbuffer->pybuffer.shape[0];
        if (unlikely((__pyx_t_15 < 0))) __PYX_ERR(0, 230, __pyx_L1_error)
      }
      __pyx_t_3 = 0;
      __Pyx_DECREF_SET(__pyx_v_best_g, ((PyArrayObject *)__pyx_t_10));
      __pyx_t_10 = 0;
      __pyx_v_best_it = __pyx_t_13;
      __pyx_v_best_obj = __pyx_t_14;

 231: 

+232:     return best_theta, best_alpha, best_g, best_it, best_obj

  __Pyx_XDECREF(__pyx_r);
  __pyx_t_12 = __Pyx_PyInt_From_int(__pyx_v_best_it); if (unlikely(!__pyx_t_12)) __PYX_ERR(0, 232, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_12);
  __pyx_t_11 = PyFloat_FromDouble(__pyx_v_best_obj); if (unlikely(!__pyx_t_11)) __PYX_ERR(0, 232, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_11);
  __pyx_t_10 = PyTuple_New(5); if (unlikely(!__pyx_t_10)) __PYX_ERR(0, 232, __pyx_L1_error)
  __Pyx_GOTREF(__pyx_t_10);
  __Pyx_INCREF((PyObject *)__pyx_v_best_theta);
  __Pyx_GIVEREF((PyObject *)__pyx_v_best_theta);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_10, 0, ((PyObject *)__pyx_v_best_theta))) __PYX_ERR(0, 232, __pyx_L1_error);
  __Pyx_INCREF((PyObject *)__pyx_v_best_alpha);
  __Pyx_GIVEREF((PyObject *)__pyx_v_best_alpha);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_10, 1, ((PyObject *)__pyx_v_best_alpha))) __PYX_ERR(0, 232, __pyx_L1_error);
  __Pyx_INCREF((PyObject *)__pyx_v_best_g);
  __Pyx_GIVEREF((PyObject *)__pyx_v_best_g);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_10, 2, ((PyObject *)__pyx_v_best_g))) __PYX_ERR(0, 232, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_12);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_10, 3, __pyx_t_12)) __PYX_ERR(0, 232, __pyx_L1_error);
  __Pyx_GIVEREF(__pyx_t_11);
  if (__Pyx_PyTuple_SET_ITEM(__pyx_t_10, 4, __pyx_t_11)) __PYX_ERR(0, 232, __pyx_L1_error);
  __pyx_t_12 = 0;
  __pyx_t_11 = 0;
  __pyx_r = __pyx_t_10;
  __pyx_t_10 = 0;
  goto __pyx_L0;