diff --git a/GPy/util/choleskies.py b/GPy/util/choleskies.py index 7afc3098..a12631aa 100644 --- a/GPy/util/choleskies.py +++ b/GPy/util/choleskies.py @@ -100,7 +100,7 @@ def indexes_to_fix_for_low_rank(rank, size): if config.getboolean('cython', 'working'): triang_to_flat = _triang_to_flat_cython flat_to_triang = _flat_to_triang_cython - backprop_gradient = choleskies_cython.backprop_gradient + backprop_gradient = choleskies_cython.backprop_gradient_par_c else: backprop_gradient = _backprop_gradient_pure triang_to_flat = _triang_to_flat_pure diff --git a/GPy/util/choleskies_cython.c b/GPy/util/choleskies_cython.c index c5d1867e..db7d1aed 100644 --- a/GPy/util/choleskies_cython.c +++ b/GPy/util/choleskies_cython.c @@ -1470,6 +1470,7 @@ static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_2triang_to_flat(CYTHON_ static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_4backprop_gradient(CYTHON_UNUSED PyObject *__pyx_self, PyArrayObject *__pyx_v_dL, PyArrayObject *__pyx_v_L); /* proto */ static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_6backprop_gradient_par(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_dL, __Pyx_memviewslice __pyx_v_L); /* proto */ static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_8backprop_gradient_par_c(CYTHON_UNUSED PyObject *__pyx_self, PyArrayObject *__pyx_v_dL, PyArrayObject *__pyx_v_L); /* proto */ +static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_10backprop_gradient_par2(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_dL, __Pyx_memviewslice __pyx_v_L); /* proto */ static int __pyx_pf_5numpy_7ndarray___getbuffer__(PyArrayObject *__pyx_v_self, Py_buffer *__pyx_v_info, int __pyx_v_flags); /* proto */ static void __pyx_pf_5numpy_7ndarray_2__releasebuffer__(PyArrayObject *__pyx_v_self, Py_buffer *__pyx_v_info); /* proto */ static int __pyx_array_MemoryView_5array___cinit__(struct __pyx_array_obj *__pyx_v_self, PyObject *__pyx_v_shape, Py_ssize_t __pyx_v_itemsize, PyObject *__pyx_v_format, PyObject *__pyx_v_mode, int __pyx_v_allocate_buffer); /* proto */ @@ -1533,9 +1534,12 @@ static char __pyx_k_Zd[] = "Zd"; static char __pyx_k_Zf[] = "Zf"; static char __pyx_k_Zg[] = "Zg"; static char __pyx_k_dL[] = "dL"; +static char __pyx_k_iN[] = "iN"; static char __pyx_k_id[] = "id"; +static char __pyx_k_kN[] = "kN"; static char __pyx_k_mm[] = "mm"; static char __pyx_k_np[] = "np"; +static char __pyx_k_dot[] = "dot"; static char __pyx_k_obj[] = "obj"; static char __pyx_k_ret[] = "ret"; static char __pyx_k_base[] = "base"; @@ -1568,6 +1572,7 @@ static char __pyx_k_name_2[] = "__name__"; static char __pyx_k_struct[] = "struct"; static char __pyx_k_unpack[] = "unpack"; static char __pyx_k_xrange[] = "xrange"; +static char __pyx_k_asarray[] = "asarray"; static char __pyx_k_fortran[] = "fortran"; static char __pyx_k_memview[] = "memview"; static char __pyx_k_Ellipsis[] = "Ellipsis"; @@ -1590,6 +1595,7 @@ static char __pyx_k_strided_and_indirect[] = ""; static char __pyx_k_backprop_gradient_par[] = "backprop_gradient_par"; static char __pyx_k_contiguous_and_direct[] = ""; static char __pyx_k_MemoryView_of_r_object[] = ""; +static char __pyx_k_backprop_gradient_par2[] = "backprop_gradient_par2"; static char __pyx_k_MemoryView_of_r_at_0x_x[] = ""; static char __pyx_k_backprop_gradient_par_c[] = "backprop_gradient_par_c"; static char __pyx_k_contiguous_and_indirect[] = ""; @@ -1604,7 +1610,7 @@ static char __pyx_k_itemsize_0_for_cython_array[] = "itemsize <= 0 for cython.ar static char __pyx_k_ndarray_is_not_C_contiguous[] = "ndarray is not C contiguous"; static char __pyx_k_unable_to_allocate_array_data[] = "unable to allocate array data."; static char __pyx_k_strided_and_direct_or_indirect[] = ""; -static char __pyx_k_home_james_work_GPy_GPy_util_ch[] = "/home/james/work/GPy/GPy/util/choleskies_cython.pyx"; +static char __pyx_k_Users_james_work_GPy_GPy_util_c[] = "/Users/james/work/GPy/GPy/util/choleskies_cython.pyx"; static char __pyx_k_unknown_dtype_code_in_numpy_pxd[] = "unknown dtype code in numpy.pxd (%d)"; static char __pyx_k_All_dimensions_preceding_dimensi[] = "All dimensions preceding dimension %d must be indexed and not sliced"; static char __pyx_k_Buffer_view_does_not_expose_stri[] = "Buffer view does not expose strides"; @@ -1646,10 +1652,13 @@ static PyObject *__pyx_kp_s_Out_of_bounds_on_buffer_access_a; static PyObject *__pyx_n_s_RuntimeError; static PyObject *__pyx_n_s_TypeError; static PyObject *__pyx_kp_s_Unable_to_convert_item_to_object; +static PyObject *__pyx_kp_s_Users_james_work_GPy_GPy_util_c; static PyObject *__pyx_n_s_ValueError; static PyObject *__pyx_n_s_allocate_buffer; +static PyObject *__pyx_n_s_asarray; static PyObject *__pyx_n_s_backprop_gradient; static PyObject *__pyx_n_s_backprop_gradient_par; +static PyObject *__pyx_n_s_backprop_gradient_par2; static PyObject *__pyx_n_s_backprop_gradient_par_c; static PyObject *__pyx_n_s_base; static PyObject *__pyx_n_s_c; @@ -1662,6 +1671,7 @@ static PyObject *__pyx_n_s_count; static PyObject *__pyx_n_s_d; static PyObject *__pyx_n_s_dL; static PyObject *__pyx_n_s_dL_dK; +static PyObject *__pyx_n_s_dot; static PyObject *__pyx_n_s_dtype_is_object; static PyObject *__pyx_n_s_empty; static PyObject *__pyx_n_s_enumerate; @@ -1673,14 +1683,15 @@ static PyObject *__pyx_n_s_format; static PyObject *__pyx_n_s_fortran; static PyObject *__pyx_n_u_fortran; static PyObject *__pyx_kp_s_got_differing_extents_in_dimensi; -static PyObject *__pyx_kp_s_home_james_work_GPy_GPy_util_ch; static PyObject *__pyx_n_s_i; +static PyObject *__pyx_n_s_iN; static PyObject *__pyx_n_s_id; static PyObject *__pyx_n_s_import; static PyObject *__pyx_n_s_itemsize; static PyObject *__pyx_kp_s_itemsize_0_for_cython_array; static PyObject *__pyx_n_s_j; static PyObject *__pyx_n_s_k; +static PyObject *__pyx_n_s_kN; static PyObject *__pyx_n_s_m; static PyObject *__pyx_n_s_main; static PyObject *__pyx_n_s_memview; @@ -1743,15 +1754,17 @@ static PyObject *__pyx_tuple__22; static PyObject *__pyx_tuple__24; static PyObject *__pyx_tuple__26; static PyObject *__pyx_tuple__28; -static PyObject *__pyx_tuple__29; static PyObject *__pyx_tuple__30; static PyObject *__pyx_tuple__31; static PyObject *__pyx_tuple__32; +static PyObject *__pyx_tuple__33; +static PyObject *__pyx_tuple__34; static PyObject *__pyx_codeobj__19; static PyObject *__pyx_codeobj__21; static PyObject *__pyx_codeobj__23; static PyObject *__pyx_codeobj__25; static PyObject *__pyx_codeobj__27; +static PyObject *__pyx_codeobj__29; /* "GPy/util/choleskies_cython.pyx":11 * cimport numpy as np @@ -3043,7 +3056,7 @@ static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_6backprop_gradient_par( #define unlikely(x) (x) #endif #ifdef _OPENMP - #pragma omp parallel private(__pyx_t_26, __pyx_t_14, __pyx_t_20, __pyx_t_27, __pyx_t_22, __pyx_t_10, __pyx_t_11, __pyx_t_23, __pyx_t_9, __pyx_t_25, __pyx_t_16, __pyx_t_21, __pyx_t_17, __pyx_t_19, __pyx_t_15, __pyx_t_13, __pyx_t_12, __pyx_t_18, __pyx_t_24) + #pragma omp parallel private(__pyx_t_9, __pyx_t_25, __pyx_t_20, __pyx_t_16, __pyx_t_21, __pyx_t_17, __pyx_t_14, __pyx_t_19, __pyx_t_22, __pyx_t_10, __pyx_t_11, __pyx_t_23, __pyx_t_13, __pyx_t_12, __pyx_t_18, __pyx_t_24, __pyx_t_26, __pyx_t_15, __pyx_t_27) #endif /* _OPENMP */ { @@ -3240,7 +3253,7 @@ static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_6backprop_gradient_par( * dL_dK[k, k] /= (2. * L[k, k]) * return dL_dK # <<<<<<<<<<<<<< * - * cdef extern from "cholesky_backprop.h": + * #here's a pure C version... */ __Pyx_XDECREF(__pyx_r); __pyx_t_1 = __pyx_memoryview_fromslice(__pyx_v_dL_dK, 2, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 77; __pyx_clineno = __LINE__; goto __pyx_L1_error;} @@ -3277,7 +3290,7 @@ static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_6backprop_gradient_par( return __pyx_r; } -/* "GPy/util/choleskies_cython.pyx":82 +/* "GPy/util/choleskies_cython.pyx":83 * void chol_backprop(int N, double* dL, double* L) * * def backprop_gradient_par_c(np.ndarray[double, ndim=2] dL, np.ndarray[double, ndim=2] L): # <<<<<<<<<<<<<< @@ -3317,11 +3330,11 @@ static PyObject *__pyx_pw_3GPy_4util_17choleskies_cython_9backprop_gradient_par_ case 1: if (likely((values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_L)) != 0)) kw_args--; else { - __Pyx_RaiseArgtupleInvalid("backprop_gradient_par_c", 1, 2, 2, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 82; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __Pyx_RaiseArgtupleInvalid("backprop_gradient_par_c", 1, 2, 2, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L3_error;} } } if (unlikely(kw_args > 0)) { - if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "backprop_gradient_par_c") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 82; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "backprop_gradient_par_c") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L3_error;} } } else if (PyTuple_GET_SIZE(__pyx_args) != 2) { goto __pyx_L5_argtuple_error; @@ -3334,14 +3347,14 @@ static PyObject *__pyx_pw_3GPy_4util_17choleskies_cython_9backprop_gradient_par_ } goto __pyx_L4_argument_unpacking_done; __pyx_L5_argtuple_error:; - __Pyx_RaiseArgtupleInvalid("backprop_gradient_par_c", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 82; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __Pyx_RaiseArgtupleInvalid("backprop_gradient_par_c", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L3_error;} __pyx_L3_error:; __Pyx_AddTraceback("GPy.util.choleskies_cython.backprop_gradient_par_c", __pyx_clineno, __pyx_lineno, __pyx_filename); __Pyx_RefNannyFinishContext(); return NULL; __pyx_L4_argument_unpacking_done:; - if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_dL), __pyx_ptype_5numpy_ndarray, 1, "dL", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 82; __pyx_clineno = __LINE__; goto __pyx_L1_error;} - if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_L), __pyx_ptype_5numpy_ndarray, 1, "L", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 82; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_dL), __pyx_ptype_5numpy_ndarray, 1, "dL", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_L), __pyx_ptype_5numpy_ndarray, 1, "L", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __pyx_r = __pyx_pf_3GPy_4util_17choleskies_cython_8backprop_gradient_par_c(__pyx_self, __pyx_v_dL, __pyx_v_L); /* function exit code */ @@ -3387,25 +3400,25 @@ static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_8backprop_gradient_par_ __pyx_pybuffernd_L.rcbuffer = &__pyx_pybuffer_L; { __Pyx_BufFmt_StackElem __pyx_stack[1]; - if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_dL.rcbuffer->pybuffer, (PyObject*)__pyx_v_dL, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 82; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_dL.rcbuffer->pybuffer, (PyObject*)__pyx_v_dL, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} } __pyx_pybuffernd_dL.diminfo[0].strides = __pyx_pybuffernd_dL.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_dL.diminfo[0].shape = __pyx_pybuffernd_dL.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_dL.diminfo[1].strides = __pyx_pybuffernd_dL.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_dL.diminfo[1].shape = __pyx_pybuffernd_dL.rcbuffer->pybuffer.shape[1]; { __Pyx_BufFmt_StackElem __pyx_stack[1]; - if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_L.rcbuffer->pybuffer, (PyObject*)__pyx_v_L, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 82; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_L.rcbuffer->pybuffer, (PyObject*)__pyx_v_L, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} } __pyx_pybuffernd_L.diminfo[0].strides = __pyx_pybuffernd_L.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_L.diminfo[0].shape = __pyx_pybuffernd_L.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_L.diminfo[1].strides = __pyx_pybuffernd_L.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_L.diminfo[1].shape = __pyx_pybuffernd_L.rcbuffer->pybuffer.shape[1]; - /* "GPy/util/choleskies_cython.pyx":83 + /* "GPy/util/choleskies_cython.pyx":84 * * def backprop_gradient_par_c(np.ndarray[double, ndim=2] dL, np.ndarray[double, ndim=2] L): * cdef np.ndarray[double, ndim=2] dL_dK = np.tril(dL) # makes a copy, c-contig # <<<<<<<<<<<<<< * cdef int N = L.shape[0] - * chol_backprop(N, dL_dK.data, L.data) + * with nogil: */ - __pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 84; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_t_2); - __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_tril); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_tril); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 84; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_t_3); __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0; __pyx_t_2 = NULL; @@ -3419,27 +3432,27 @@ static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_8backprop_gradient_par_ } } if (!__pyx_t_2) { - __pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_3, ((PyObject *)__pyx_v_dL)); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_3, ((PyObject *)__pyx_v_dL)); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 84; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_t_1); } else { - __pyx_t_4 = PyTuple_New(1+1); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_t_4 = PyTuple_New(1+1); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 84; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_t_4); PyTuple_SET_ITEM(__pyx_t_4, 0, __pyx_t_2); __Pyx_GIVEREF(__pyx_t_2); __pyx_t_2 = NULL; __Pyx_INCREF(((PyObject *)__pyx_v_dL)); PyTuple_SET_ITEM(__pyx_t_4, 0+1, ((PyObject *)__pyx_v_dL)); __Pyx_GIVEREF(((PyObject *)__pyx_v_dL)); - __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_4, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_4, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 84; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_t_1); __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; } __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_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 84; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __pyx_t_5 = ((PyArrayObject *)__pyx_t_1); { __Pyx_BufFmt_StackElem __pyx_stack[1]; if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_pybuffernd_dL_dK.rcbuffer->pybuffer, (PyObject*)__pyx_t_5, &__Pyx_TypeInfo_double, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) { __pyx_v_dL_dK = ((PyArrayObject *)Py_None); __Pyx_INCREF(Py_None); __pyx_pybuffernd_dL_dK.rcbuffer->pybuffer.buf = NULL; - {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + {__pyx_filename = __pyx_f[0]; __pyx_lineno = 84; __pyx_clineno = __LINE__; goto __pyx_L1_error;} } else {__pyx_pybuffernd_dL_dK.diminfo[0].strides = __pyx_pybuffernd_dL_dK.rcbuffer->pybuffer.strides[0]; __pyx_pybuffernd_dL_dK.diminfo[0].shape = __pyx_pybuffernd_dL_dK.rcbuffer->pybuffer.shape[0]; __pyx_pybuffernd_dL_dK.diminfo[1].strides = __pyx_pybuffernd_dL_dK.rcbuffer->pybuffer.strides[1]; __pyx_pybuffernd_dL_dK.diminfo[1].shape = __pyx_pybuffernd_dL_dK.rcbuffer->pybuffer.shape[1]; } } @@ -3447,36 +3460,70 @@ static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_8backprop_gradient_par_ __pyx_v_dL_dK = ((PyArrayObject *)__pyx_t_1); __pyx_t_1 = 0; - /* "GPy/util/choleskies_cython.pyx":84 + /* "GPy/util/choleskies_cython.pyx":85 * def backprop_gradient_par_c(np.ndarray[double, ndim=2] dL, np.ndarray[double, ndim=2] L): * cdef np.ndarray[double, ndim=2] dL_dK = np.tril(dL) # makes a copy, c-contig * cdef int N = L.shape[0] # <<<<<<<<<<<<<< - * chol_backprop(N, dL_dK.data, L.data) - * return dL_dK + * with nogil: + * chol_backprop(N, dL_dK.data, L.data) */ __pyx_v_N = (__pyx_v_L->dimensions[0]); - /* "GPy/util/choleskies_cython.pyx":85 + /* "GPy/util/choleskies_cython.pyx":86 * cdef np.ndarray[double, ndim=2] dL_dK = np.tril(dL) # makes a copy, c-contig * cdef int N = L.shape[0] - * chol_backprop(N, dL_dK.data, L.data) # <<<<<<<<<<<<<< + * with nogil: # <<<<<<<<<<<<<< + * chol_backprop(N, dL_dK.data, L.data) + * return dL_dK + */ + { + #ifdef WITH_THREAD + PyThreadState *_save; + Py_UNBLOCK_THREADS + #endif + /*try:*/ { + + /* "GPy/util/choleskies_cython.pyx":87 + * cdef int N = L.shape[0] + * with nogil: + * chol_backprop(N, dL_dK.data, L.data) # <<<<<<<<<<<<<< * return dL_dK * */ - chol_backprop(__pyx_v_N, ((double *)__pyx_v_dL_dK->data), ((double *)__pyx_v_L->data)); + chol_backprop(__pyx_v_N, ((double *)__pyx_v_dL_dK->data), ((double *)__pyx_v_L->data)); + } - /* "GPy/util/choleskies_cython.pyx":86 + /* "GPy/util/choleskies_cython.pyx":86 + * cdef np.ndarray[double, ndim=2] dL_dK = np.tril(dL) # makes a copy, c-contig * cdef int N = L.shape[0] - * chol_backprop(N, dL_dK.data, L.data) + * with nogil: # <<<<<<<<<<<<<< + * chol_backprop(N, dL_dK.data, L.data) + * return dL_dK + */ + /*finally:*/ { + /*normal exit:*/{ + #ifdef WITH_THREAD + Py_BLOCK_THREADS + #endif + goto __pyx_L5; + } + __pyx_L5:; + } + } + + /* "GPy/util/choleskies_cython.pyx":88 + * with nogil: + * chol_backprop(N, dL_dK.data, L.data) * return dL_dK # <<<<<<<<<<<<<< * + * */ __Pyx_XDECREF(__pyx_r); __Pyx_INCREF(((PyObject *)__pyx_v_dL_dK)); __pyx_r = ((PyObject *)__pyx_v_dL_dK); goto __pyx_L0; - /* "GPy/util/choleskies_cython.pyx":82 + /* "GPy/util/choleskies_cython.pyx":83 * void chol_backprop(int N, double* dL, double* L) * * def backprop_gradient_par_c(np.ndarray[double, ndim=2] dL, np.ndarray[double, ndim=2] L): # <<<<<<<<<<<<<< @@ -3510,6 +3557,606 @@ static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_8backprop_gradient_par_ return __pyx_r; } +/* "GPy/util/choleskies_cython.pyx":92 + * + * # TODO: with the next release of cython, cimport scipy.linalg.cython_blas as blas, then blas the hell out of this. + * def backprop_gradient_par2(double[:,:] dL, double[:,:] L): # <<<<<<<<<<<<<< + * """ + * a very slow implementation, but the clearest I hope + */ + +/* Python wrapper */ +static PyObject *__pyx_pw_3GPy_4util_17choleskies_cython_11backprop_gradient_par2(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/ +static char __pyx_doc_3GPy_4util_17choleskies_cython_10backprop_gradient_par2[] = "\n a very slow implementation, but the clearest I hope\n "; +static PyMethodDef __pyx_mdef_3GPy_4util_17choleskies_cython_11backprop_gradient_par2 = {"backprop_gradient_par2", (PyCFunction)__pyx_pw_3GPy_4util_17choleskies_cython_11backprop_gradient_par2, METH_VARARGS|METH_KEYWORDS, __pyx_doc_3GPy_4util_17choleskies_cython_10backprop_gradient_par2}; +static PyObject *__pyx_pw_3GPy_4util_17choleskies_cython_11backprop_gradient_par2(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) { + __Pyx_memviewslice __pyx_v_dL = { 0, 0, { 0 }, { 0 }, { 0 } }; + __Pyx_memviewslice __pyx_v_L = { 0, 0, { 0 }, { 0 }, { 0 } }; + int __pyx_lineno = 0; + const char *__pyx_filename = NULL; + int __pyx_clineno = 0; + PyObject *__pyx_r = 0; + __Pyx_RefNannyDeclarations + __Pyx_RefNannySetupContext("backprop_gradient_par2 (wrapper)", 0); + { + static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_dL,&__pyx_n_s_L,0}; + PyObject* values[2] = {0,0}; + if (unlikely(__pyx_kwds)) { + Py_ssize_t kw_args; + const Py_ssize_t pos_args = PyTuple_GET_SIZE(__pyx_args); + switch (pos_args) { + case 2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1); + case 1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0); + case 0: break; + default: goto __pyx_L5_argtuple_error; + } + kw_args = PyDict_Size(__pyx_kwds); + switch (pos_args) { + case 0: + if (likely((values[0] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_dL)) != 0)) kw_args--; + else goto __pyx_L5_argtuple_error; + case 1: + if (likely((values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s_L)) != 0)) kw_args--; + else { + __Pyx_RaiseArgtupleInvalid("backprop_gradient_par2", 1, 2, 2, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + } + } + if (unlikely(kw_args > 0)) { + if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "backprop_gradient_par2") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + } + } else if (PyTuple_GET_SIZE(__pyx_args) != 2) { + goto __pyx_L5_argtuple_error; + } else { + values[0] = PyTuple_GET_ITEM(__pyx_args, 0); + values[1] = PyTuple_GET_ITEM(__pyx_args, 1); + } + __pyx_v_dL = __Pyx_PyObject_to_MemoryviewSlice_dsds_double(values[0]); if (unlikely(!__pyx_v_dL.memview)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __pyx_v_L = __Pyx_PyObject_to_MemoryviewSlice_dsds_double(values[1]); if (unlikely(!__pyx_v_L.memview)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + } + goto __pyx_L4_argument_unpacking_done; + __pyx_L5_argtuple_error:; + __Pyx_RaiseArgtupleInvalid("backprop_gradient_par2", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L3_error;} + __pyx_L3_error:; + __Pyx_AddTraceback("GPy.util.choleskies_cython.backprop_gradient_par2", __pyx_clineno, __pyx_lineno, __pyx_filename); + __Pyx_RefNannyFinishContext(); + return NULL; + __pyx_L4_argument_unpacking_done:; + __pyx_r = __pyx_pf_3GPy_4util_17choleskies_cython_10backprop_gradient_par2(__pyx_self, __pyx_v_dL, __pyx_v_L); + + /* function exit code */ + __Pyx_RefNannyFinishContext(); + return __pyx_r; +} + +static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_10backprop_gradient_par2(CYTHON_UNUSED PyObject *__pyx_self, __Pyx_memviewslice __pyx_v_dL, __Pyx_memviewslice __pyx_v_L) { + __Pyx_memviewslice __pyx_v_dL_dK = { 0, 0, { 0 }, { 0 }, { 0 } }; + int __pyx_v_N; + int __pyx_v_k; + int __pyx_v_i; + PyObject *__pyx_r = NULL; + __Pyx_RefNannyDeclarations + PyObject *__pyx_t_1 = NULL; + PyObject *__pyx_t_2 = NULL; + PyObject *__pyx_t_3 = NULL; + PyObject *__pyx_t_4 = NULL; + PyObject *__pyx_t_5 = NULL; + PyObject *__pyx_t_6 = NULL; + __Pyx_memviewslice __pyx_t_7 = { 0, 0, { 0 }, { 0 }, { 0 } }; + int __pyx_t_8; + int __pyx_t_9; + int __pyx_t_10; + __Pyx_memviewslice __pyx_t_11 = { 0, 0, { 0 }, { 0 }, { 0 } }; + int __pyx_t_12; + __Pyx_memviewslice __pyx_t_13 = { 0, 0, { 0 }, { 0 }, { 0 } }; + Py_ssize_t __pyx_t_14; + double __pyx_t_15; + int __pyx_t_16; + __Pyx_memviewslice __pyx_t_17 = { 0, 0, { 0 }, { 0 }, { 0 } }; + int __pyx_t_18; + __Pyx_memviewslice __pyx_t_19 = { 0, 0, { 0 }, { 0 }, { 0 } }; + int __pyx_t_20; + int __pyx_t_21; + int __pyx_t_22; + int __pyx_t_23; + int __pyx_t_24; + int __pyx_t_25; + int __pyx_t_26; + int __pyx_t_27; + int __pyx_t_28; + int __pyx_t_29; + int __pyx_t_30; + int __pyx_t_31; + int __pyx_t_32; + int __pyx_lineno = 0; + const char *__pyx_filename = NULL; + int __pyx_clineno = 0; + __Pyx_RefNannySetupContext("backprop_gradient_par2", 0); + + /* "GPy/util/choleskies_cython.pyx":96 + * a very slow implementation, but the clearest I hope + * """ + * cdef double[:,:] dL_dK = np.tril(dL).copy() # <<<<<<<<<<<<<< + * cdef int N = L.shape[0] + * cdef int k, j, i, iN, kN + */ + __pyx_t_3 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 96; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_3); + __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_3, __pyx_n_s_tril); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 96; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_4); + __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0; + __pyx_t_3 = __pyx_memoryview_fromslice(__pyx_v_dL, 2, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 96; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_3); + __pyx_t_5 = NULL; + if (CYTHON_COMPILING_IN_CPYTHON && unlikely(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); + } + } + if (!__pyx_t_5) { + __pyx_t_2 = __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_3); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 96; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0; + __Pyx_GOTREF(__pyx_t_2); + } else { + __pyx_t_6 = PyTuple_New(1+1); if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 96; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_6); + PyTuple_SET_ITEM(__pyx_t_6, 0, __pyx_t_5); __Pyx_GIVEREF(__pyx_t_5); __pyx_t_5 = NULL; + PyTuple_SET_ITEM(__pyx_t_6, 0+1, __pyx_t_3); + __Pyx_GIVEREF(__pyx_t_3); + __pyx_t_3 = 0; + __pyx_t_2 = __Pyx_PyObject_Call(__pyx_t_4, __pyx_t_6, NULL); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 96; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_2); + __Pyx_DECREF(__pyx_t_6); __pyx_t_6 = 0; + } + __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; + __pyx_t_4 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_copy); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 96; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_4); + __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0; + __pyx_t_2 = NULL; + if (CYTHON_COMPILING_IN_CPYTHON && likely(PyMethod_Check(__pyx_t_4))) { + __pyx_t_2 = PyMethod_GET_SELF(__pyx_t_4); + if (likely(__pyx_t_2)) { + PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_4); + __Pyx_INCREF(__pyx_t_2); + __Pyx_INCREF(function); + __Pyx_DECREF_SET(__pyx_t_4, function); + } + } + if (__pyx_t_2) { + __pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_4, __pyx_t_2); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 96; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0; + } else { + __pyx_t_1 = __Pyx_PyObject_CallNoArg(__pyx_t_4); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 96; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + } + __Pyx_GOTREF(__pyx_t_1); + __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; + __pyx_t_7 = __Pyx_PyObject_to_MemoryviewSlice_dsds_double(__pyx_t_1); + if (unlikely(!__pyx_t_7.memview)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 96; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; + __pyx_v_dL_dK = __pyx_t_7; + __pyx_t_7.memview = NULL; + __pyx_t_7.data = NULL; + + /* "GPy/util/choleskies_cython.pyx":97 + * """ + * cdef double[:,:] dL_dK = np.tril(dL).copy() + * cdef int N = L.shape[0] # <<<<<<<<<<<<<< + * cdef int k, j, i, iN, kN + * for k in range(N - 1, -1, -1): + */ + __pyx_v_N = (__pyx_v_L.shape[0]); + + /* "GPy/util/choleskies_cython.pyx":99 + * cdef int N = L.shape[0] + * cdef int k, j, i, iN, kN + * for k in range(N - 1, -1, -1): # <<<<<<<<<<<<<< + * #pragma this loop: + * for i in range(k+1, N): + */ + for (__pyx_t_8 = (__pyx_v_N - 1); __pyx_t_8 > -1; __pyx_t_8-=1) { + __pyx_v_k = __pyx_t_8; + + /* "GPy/util/choleskies_cython.pyx":101 + * for k in range(N - 1, -1, -1): + * #pragma this loop: + * for i in range(k+1, N): # <<<<<<<<<<<<<< + * dL_dK[i,+k] -= np.dot(dL_dK[i,k+1:], L[k+1:,k]) + * dL_dK[i,+k] -= np.dot(dL_dK[i:,i], L[i:,k]) + */ + __pyx_t_9 = __pyx_v_N; + for (__pyx_t_10 = (__pyx_v_k + 1); __pyx_t_10 < __pyx_t_9; __pyx_t_10+=1) { + __pyx_v_i = __pyx_t_10; + + /* "GPy/util/choleskies_cython.pyx":102 + * #pragma this loop: + * for i in range(k+1, N): + * dL_dK[i,+k] -= np.dot(dL_dK[i,k+1:], L[k+1:,k]) # <<<<<<<<<<<<<< + * dL_dK[i,+k] -= np.dot(dL_dK[i:,i], L[i:,k]) + * + */ + __pyx_t_4 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_4); + __pyx_t_2 = __Pyx_PyObject_GetAttrStr(__pyx_t_4, __pyx_n_s_dot); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_2); + __Pyx_DECREF(__pyx_t_4); __pyx_t_4 = 0; + __pyx_t_12 = -1; + __pyx_t_11.data = __pyx_v_dL_dK.data; + __pyx_t_11.memview = __pyx_v_dL_dK.memview; + __PYX_INC_MEMVIEW(&__pyx_t_11, 0); + { + Py_ssize_t __pyx_tmp_idx = __pyx_v_i; + Py_ssize_t __pyx_tmp_shape = __pyx_v_dL_dK.shape[0]; + Py_ssize_t __pyx_tmp_stride = __pyx_v_dL_dK.strides[0]; + if (1 && (__pyx_tmp_idx < 0)) + __pyx_tmp_idx += __pyx_tmp_shape; + if (0 && (__pyx_tmp_idx < 0 || __pyx_tmp_idx >= __pyx_tmp_shape)) { + PyErr_SetString(PyExc_IndexError, "Index out of bounds (axis 0)"); + {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + } + __pyx_t_11.data += __pyx_tmp_idx * __pyx_tmp_stride; +} + +if (unlikely(__pyx_memoryview_slice_memviewslice( + &__pyx_t_11, + __pyx_v_dL_dK.shape[1], __pyx_v_dL_dK.strides[1], __pyx_v_dL_dK.suboffsets[1], + 1, + 0, + &__pyx_t_12, + (__pyx_v_k + 1), + 0, + 0, + 1, + 0, + 0, + 1) < 0)) +{ + {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;} +} + +__pyx_t_4 = __pyx_memoryview_fromslice(__pyx_t_11, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_4); + __PYX_XDEC_MEMVIEW(&__pyx_t_11, 1); + __pyx_t_12 = -1; + __pyx_t_13.data = __pyx_v_L.data; + __pyx_t_13.memview = __pyx_v_L.memview; + __PYX_INC_MEMVIEW(&__pyx_t_13, 0); + if (unlikely(__pyx_memoryview_slice_memviewslice( + &__pyx_t_13, + __pyx_v_L.shape[0], __pyx_v_L.strides[0], __pyx_v_L.suboffsets[0], + 0, + 0, + &__pyx_t_12, + (__pyx_v_k + 1), + 0, + 0, + 1, + 0, + 0, + 1) < 0)) +{ + {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;} +} + +{ + Py_ssize_t __pyx_tmp_idx = __pyx_v_k; + Py_ssize_t __pyx_tmp_shape = __pyx_v_L.shape[1]; + Py_ssize_t __pyx_tmp_stride = __pyx_v_L.strides[1]; + if (1 && (__pyx_tmp_idx < 0)) + __pyx_tmp_idx += __pyx_tmp_shape; + if (0 && (__pyx_tmp_idx < 0 || __pyx_tmp_idx >= __pyx_tmp_shape)) { + PyErr_SetString(PyExc_IndexError, "Index out of bounds (axis 1)"); + {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + } + __pyx_t_13.data += __pyx_tmp_idx * __pyx_tmp_stride; +} + +__pyx_t_6 = __pyx_memoryview_fromslice(__pyx_t_13, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_6); + __PYX_XDEC_MEMVIEW(&__pyx_t_13, 1); + __pyx_t_3 = NULL; + __pyx_t_14 = 0; + if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_2))) { + __pyx_t_3 = PyMethod_GET_SELF(__pyx_t_2); + if (likely(__pyx_t_3)) { + PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_2); + __Pyx_INCREF(__pyx_t_3); + __Pyx_INCREF(function); + __Pyx_DECREF_SET(__pyx_t_2, function); + __pyx_t_14 = 1; + } + } + __pyx_t_5 = PyTuple_New(2+__pyx_t_14); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_5); + if (__pyx_t_3) { + PyTuple_SET_ITEM(__pyx_t_5, 0, __pyx_t_3); __Pyx_GIVEREF(__pyx_t_3); __pyx_t_3 = NULL; + } + PyTuple_SET_ITEM(__pyx_t_5, 0+__pyx_t_14, __pyx_t_4); + __Pyx_GIVEREF(__pyx_t_4); + PyTuple_SET_ITEM(__pyx_t_5, 1+__pyx_t_14, __pyx_t_6); + __Pyx_GIVEREF(__pyx_t_6); + __pyx_t_4 = 0; + __pyx_t_6 = 0; + __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_2, __pyx_t_5, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_1); + __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0; + __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0; + __pyx_t_15 = __pyx_PyFloat_AsDouble(__pyx_t_1); if (unlikely((__pyx_t_15 == (double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 102; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; + __pyx_t_12 = __pyx_v_i; + __pyx_t_16 = __pyx_v_k; + if (__pyx_t_12 < 0) __pyx_t_12 += __pyx_v_dL_dK.shape[0]; + if (__pyx_t_16 < 0) __pyx_t_16 += __pyx_v_dL_dK.shape[1]; + *((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_dL_dK.data + __pyx_t_12 * __pyx_v_dL_dK.strides[0]) ) + __pyx_t_16 * __pyx_v_dL_dK.strides[1]) )) -= __pyx_t_15; + + /* "GPy/util/choleskies_cython.pyx":103 + * for i in range(k+1, N): + * dL_dK[i,+k] -= np.dot(dL_dK[i,k+1:], L[k+1:,k]) + * dL_dK[i,+k] -= np.dot(dL_dK[i:,i], L[i:,k]) # <<<<<<<<<<<<<< + * + * for i in range(k+1, N): + */ + __pyx_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 103; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_2); + __pyx_t_5 = __Pyx_PyObject_GetAttrStr(__pyx_t_2, __pyx_n_s_dot); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 103; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_5); + __Pyx_DECREF(__pyx_t_2); __pyx_t_2 = 0; + __pyx_t_18 = -1; + __pyx_t_17.data = __pyx_v_dL_dK.data; + __pyx_t_17.memview = __pyx_v_dL_dK.memview; + __PYX_INC_MEMVIEW(&__pyx_t_17, 0); + if (unlikely(__pyx_memoryview_slice_memviewslice( + &__pyx_t_17, + __pyx_v_dL_dK.shape[0], __pyx_v_dL_dK.strides[0], __pyx_v_dL_dK.suboffsets[0], + 0, + 0, + &__pyx_t_18, + __pyx_v_i, + 0, + 0, + 1, + 0, + 0, + 1) < 0)) +{ + {__pyx_filename = __pyx_f[0]; __pyx_lineno = 103; __pyx_clineno = __LINE__; goto __pyx_L1_error;} +} + +{ + Py_ssize_t __pyx_tmp_idx = __pyx_v_i; + Py_ssize_t __pyx_tmp_shape = __pyx_v_dL_dK.shape[1]; + Py_ssize_t __pyx_tmp_stride = __pyx_v_dL_dK.strides[1]; + if (1 && (__pyx_tmp_idx < 0)) + __pyx_tmp_idx += __pyx_tmp_shape; + if (0 && (__pyx_tmp_idx < 0 || __pyx_tmp_idx >= __pyx_tmp_shape)) { + PyErr_SetString(PyExc_IndexError, "Index out of bounds (axis 1)"); + {__pyx_filename = __pyx_f[0]; __pyx_lineno = 103; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + } + __pyx_t_17.data += __pyx_tmp_idx * __pyx_tmp_stride; +} + +__pyx_t_2 = __pyx_memoryview_fromslice(__pyx_t_17, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 103; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_2); + __PYX_XDEC_MEMVIEW(&__pyx_t_17, 1); + __pyx_t_18 = -1; + __pyx_t_19.data = __pyx_v_L.data; + __pyx_t_19.memview = __pyx_v_L.memview; + __PYX_INC_MEMVIEW(&__pyx_t_19, 0); + if (unlikely(__pyx_memoryview_slice_memviewslice( + &__pyx_t_19, + __pyx_v_L.shape[0], __pyx_v_L.strides[0], __pyx_v_L.suboffsets[0], + 0, + 0, + &__pyx_t_18, + __pyx_v_i, + 0, + 0, + 1, + 0, + 0, + 1) < 0)) +{ + {__pyx_filename = __pyx_f[0]; __pyx_lineno = 103; __pyx_clineno = __LINE__; goto __pyx_L1_error;} +} + +{ + Py_ssize_t __pyx_tmp_idx = __pyx_v_k; + Py_ssize_t __pyx_tmp_shape = __pyx_v_L.shape[1]; + Py_ssize_t __pyx_tmp_stride = __pyx_v_L.strides[1]; + if (1 && (__pyx_tmp_idx < 0)) + __pyx_tmp_idx += __pyx_tmp_shape; + if (0 && (__pyx_tmp_idx < 0 || __pyx_tmp_idx >= __pyx_tmp_shape)) { + PyErr_SetString(PyExc_IndexError, "Index out of bounds (axis 1)"); + {__pyx_filename = __pyx_f[0]; __pyx_lineno = 103; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + } + __pyx_t_19.data += __pyx_tmp_idx * __pyx_tmp_stride; +} + +__pyx_t_6 = __pyx_memoryview_fromslice(__pyx_t_19, 1, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_6)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 103; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_6); + __PYX_XDEC_MEMVIEW(&__pyx_t_19, 1); + __pyx_t_4 = NULL; + __pyx_t_14 = 0; + if (CYTHON_COMPILING_IN_CPYTHON && 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_14 = 1; + } + } + __pyx_t_3 = PyTuple_New(2+__pyx_t_14); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 103; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_3); + if (__pyx_t_4) { + PyTuple_SET_ITEM(__pyx_t_3, 0, __pyx_t_4); __Pyx_GIVEREF(__pyx_t_4); __pyx_t_4 = NULL; + } + PyTuple_SET_ITEM(__pyx_t_3, 0+__pyx_t_14, __pyx_t_2); + __Pyx_GIVEREF(__pyx_t_2); + PyTuple_SET_ITEM(__pyx_t_3, 1+__pyx_t_14, __pyx_t_6); + __Pyx_GIVEREF(__pyx_t_6); + __pyx_t_2 = 0; + __pyx_t_6 = 0; + __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_5, __pyx_t_3, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 103; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_1); + __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0; + __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0; + __pyx_t_15 = __pyx_PyFloat_AsDouble(__pyx_t_1); if (unlikely((__pyx_t_15 == (double)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 103; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; + __pyx_t_18 = __pyx_v_i; + __pyx_t_20 = __pyx_v_k; + if (__pyx_t_18 < 0) __pyx_t_18 += __pyx_v_dL_dK.shape[0]; + if (__pyx_t_20 < 0) __pyx_t_20 += __pyx_v_dL_dK.shape[1]; + *((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_dL_dK.data + __pyx_t_18 * __pyx_v_dL_dK.strides[0]) ) + __pyx_t_20 * __pyx_v_dL_dK.strides[1]) )) -= __pyx_t_15; + } + + /* "GPy/util/choleskies_cython.pyx":105 + * dL_dK[i,+k] -= np.dot(dL_dK[i:,i], L[i:,k]) + * + * for i in range(k+1, N): # <<<<<<<<<<<<<< + * dL_dK[i, k] /= L[k, k]; + * dL_dK[k, k] -= L[i, k] * dL_dK[i, k]; + */ + __pyx_t_9 = __pyx_v_N; + for (__pyx_t_10 = (__pyx_v_k + 1); __pyx_t_10 < __pyx_t_9; __pyx_t_10+=1) { + __pyx_v_i = __pyx_t_10; + + /* "GPy/util/choleskies_cython.pyx":106 + * + * for i in range(k+1, N): + * dL_dK[i, k] /= L[k, k]; # <<<<<<<<<<<<<< + * dL_dK[k, k] -= L[i, k] * dL_dK[i, k]; + * + */ + __pyx_t_21 = __pyx_v_k; + __pyx_t_22 = __pyx_v_k; + if (__pyx_t_21 < 0) __pyx_t_21 += __pyx_v_L.shape[0]; + if (__pyx_t_22 < 0) __pyx_t_22 += __pyx_v_L.shape[1]; + __pyx_t_23 = __pyx_v_i; + __pyx_t_24 = __pyx_v_k; + if (__pyx_t_23 < 0) __pyx_t_23 += __pyx_v_dL_dK.shape[0]; + if (__pyx_t_24 < 0) __pyx_t_24 += __pyx_v_dL_dK.shape[1]; + *((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_dL_dK.data + __pyx_t_23 * __pyx_v_dL_dK.strides[0]) ) + __pyx_t_24 * __pyx_v_dL_dK.strides[1]) )) /= (*((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_L.data + __pyx_t_21 * __pyx_v_L.strides[0]) ) + __pyx_t_22 * __pyx_v_L.strides[1]) ))); + + /* "GPy/util/choleskies_cython.pyx":107 + * for i in range(k+1, N): + * dL_dK[i, k] /= L[k, k]; + * dL_dK[k, k] -= L[i, k] * dL_dK[i, k]; # <<<<<<<<<<<<<< + * + * dL_dK[k, k] /= (2. * L[k, k]); + */ + __pyx_t_25 = __pyx_v_i; + __pyx_t_26 = __pyx_v_k; + if (__pyx_t_25 < 0) __pyx_t_25 += __pyx_v_L.shape[0]; + if (__pyx_t_26 < 0) __pyx_t_26 += __pyx_v_L.shape[1]; + __pyx_t_27 = __pyx_v_i; + __pyx_t_28 = __pyx_v_k; + if (__pyx_t_27 < 0) __pyx_t_27 += __pyx_v_dL_dK.shape[0]; + if (__pyx_t_28 < 0) __pyx_t_28 += __pyx_v_dL_dK.shape[1]; + __pyx_t_29 = __pyx_v_k; + __pyx_t_30 = __pyx_v_k; + if (__pyx_t_29 < 0) __pyx_t_29 += __pyx_v_dL_dK.shape[0]; + if (__pyx_t_30 < 0) __pyx_t_30 += __pyx_v_dL_dK.shape[1]; + *((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_dL_dK.data + __pyx_t_29 * __pyx_v_dL_dK.strides[0]) ) + __pyx_t_30 * __pyx_v_dL_dK.strides[1]) )) -= ((*((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_L.data + __pyx_t_25 * __pyx_v_L.strides[0]) ) + __pyx_t_26 * __pyx_v_L.strides[1]) ))) * (*((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_dL_dK.data + __pyx_t_27 * __pyx_v_dL_dK.strides[0]) ) + __pyx_t_28 * __pyx_v_dL_dK.strides[1]) )))); + } + + /* "GPy/util/choleskies_cython.pyx":109 + * dL_dK[k, k] -= L[i, k] * dL_dK[i, k]; + * + * dL_dK[k, k] /= (2. * L[k, k]); # <<<<<<<<<<<<<< + * return np.asarray(dL_dK) + * + */ + __pyx_t_9 = __pyx_v_k; + __pyx_t_10 = __pyx_v_k; + if (__pyx_t_9 < 0) __pyx_t_9 += __pyx_v_L.shape[0]; + if (__pyx_t_10 < 0) __pyx_t_10 += __pyx_v_L.shape[1]; + __pyx_t_31 = __pyx_v_k; + __pyx_t_32 = __pyx_v_k; + if (__pyx_t_31 < 0) __pyx_t_31 += __pyx_v_dL_dK.shape[0]; + if (__pyx_t_32 < 0) __pyx_t_32 += __pyx_v_dL_dK.shape[1]; + *((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_dL_dK.data + __pyx_t_31 * __pyx_v_dL_dK.strides[0]) ) + __pyx_t_32 * __pyx_v_dL_dK.strides[1]) )) /= (2. * (*((double *) ( /* dim=1 */ (( /* dim=0 */ (__pyx_v_L.data + __pyx_t_9 * __pyx_v_L.strides[0]) ) + __pyx_t_10 * __pyx_v_L.strides[1]) )))); + } + + /* "GPy/util/choleskies_cython.pyx":110 + * + * dL_dK[k, k] /= (2. * L[k, k]); + * return np.asarray(dL_dK) # <<<<<<<<<<<<<< + * + */ + __Pyx_XDECREF(__pyx_r); + __pyx_t_5 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 110; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_5); + __pyx_t_3 = __Pyx_PyObject_GetAttrStr(__pyx_t_5, __pyx_n_s_asarray); if (unlikely(!__pyx_t_3)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 110; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_3); + __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0; + __pyx_t_5 = __pyx_memoryview_fromslice(__pyx_v_dL_dK, 2, (PyObject *(*)(char *)) __pyx_memview_get_double, (int (*)(char *, PyObject *)) __pyx_memview_set_double, 0);; if (unlikely(!__pyx_t_5)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 110; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_5); + __pyx_t_6 = NULL; + if (CYTHON_COMPILING_IN_CPYTHON && unlikely(PyMethod_Check(__pyx_t_3))) { + __pyx_t_6 = PyMethod_GET_SELF(__pyx_t_3); + if (likely(__pyx_t_6)) { + PyObject* function = PyMethod_GET_FUNCTION(__pyx_t_3); + __Pyx_INCREF(__pyx_t_6); + __Pyx_INCREF(function); + __Pyx_DECREF_SET(__pyx_t_3, function); + } + } + if (!__pyx_t_6) { + __pyx_t_1 = __Pyx_PyObject_CallOneArg(__pyx_t_3, __pyx_t_5); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 110; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_DECREF(__pyx_t_5); __pyx_t_5 = 0; + __Pyx_GOTREF(__pyx_t_1); + } else { + __pyx_t_2 = PyTuple_New(1+1); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 110; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_2); + PyTuple_SET_ITEM(__pyx_t_2, 0, __pyx_t_6); __Pyx_GIVEREF(__pyx_t_6); __pyx_t_6 = NULL; + PyTuple_SET_ITEM(__pyx_t_2, 0+1, __pyx_t_5); + __Pyx_GIVEREF(__pyx_t_5); + __pyx_t_5 = 0; + __pyx_t_1 = __Pyx_PyObject_Call(__pyx_t_3, __pyx_t_2, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 110; __pyx_clineno = __LINE__; goto __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_r = __pyx_t_1; + __pyx_t_1 = 0; + goto __pyx_L0; + + /* "GPy/util/choleskies_cython.pyx":92 + * + * # TODO: with the next release of cython, cimport scipy.linalg.cython_blas as blas, then blas the hell out of this. + * def backprop_gradient_par2(double[:,:] dL, double[:,:] L): # <<<<<<<<<<<<<< + * """ + * a very slow implementation, but the clearest I hope + */ + + /* function exit code */ + __pyx_L1_error:; + __Pyx_XDECREF(__pyx_t_1); + __Pyx_XDECREF(__pyx_t_2); + __Pyx_XDECREF(__pyx_t_3); + __Pyx_XDECREF(__pyx_t_4); + __Pyx_XDECREF(__pyx_t_5); + __Pyx_XDECREF(__pyx_t_6); + __PYX_XDEC_MEMVIEW(&__pyx_t_7, 1); + __PYX_XDEC_MEMVIEW(&__pyx_t_11, 1); + __PYX_XDEC_MEMVIEW(&__pyx_t_13, 1); + __PYX_XDEC_MEMVIEW(&__pyx_t_17, 1); + __PYX_XDEC_MEMVIEW(&__pyx_t_19, 1); + __Pyx_AddTraceback("GPy.util.choleskies_cython.backprop_gradient_par2", __pyx_clineno, __pyx_lineno, __pyx_filename); + __pyx_r = NULL; + __pyx_L0:; + __PYX_XDEC_MEMVIEW(&__pyx_v_dL_dK, 1); + __PYX_XDEC_MEMVIEW(&__pyx_v_dL, 1); + __PYX_XDEC_MEMVIEW(&__pyx_v_L, 1); + __Pyx_XGIVEREF(__pyx_r); + __Pyx_RefNannyFinishContext(); + return __pyx_r; +} + /* "../../../../anaconda/lib/python2.7/site-packages/Cython/Includes/numpy/__init__.pxd":194 * # experimental exception made for __getbuffer__ and __releasebuffer__ * # -- the details of this may change. @@ -16885,10 +17532,13 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = { {&__pyx_n_s_RuntimeError, __pyx_k_RuntimeError, sizeof(__pyx_k_RuntimeError), 0, 0, 1, 1}, {&__pyx_n_s_TypeError, __pyx_k_TypeError, sizeof(__pyx_k_TypeError), 0, 0, 1, 1}, {&__pyx_kp_s_Unable_to_convert_item_to_object, __pyx_k_Unable_to_convert_item_to_object, sizeof(__pyx_k_Unable_to_convert_item_to_object), 0, 0, 1, 0}, + {&__pyx_kp_s_Users_james_work_GPy_GPy_util_c, __pyx_k_Users_james_work_GPy_GPy_util_c, sizeof(__pyx_k_Users_james_work_GPy_GPy_util_c), 0, 0, 1, 0}, {&__pyx_n_s_ValueError, __pyx_k_ValueError, sizeof(__pyx_k_ValueError), 0, 0, 1, 1}, {&__pyx_n_s_allocate_buffer, __pyx_k_allocate_buffer, sizeof(__pyx_k_allocate_buffer), 0, 0, 1, 1}, + {&__pyx_n_s_asarray, __pyx_k_asarray, sizeof(__pyx_k_asarray), 0, 0, 1, 1}, {&__pyx_n_s_backprop_gradient, __pyx_k_backprop_gradient, sizeof(__pyx_k_backprop_gradient), 0, 0, 1, 1}, {&__pyx_n_s_backprop_gradient_par, __pyx_k_backprop_gradient_par, sizeof(__pyx_k_backprop_gradient_par), 0, 0, 1, 1}, + {&__pyx_n_s_backprop_gradient_par2, __pyx_k_backprop_gradient_par2, sizeof(__pyx_k_backprop_gradient_par2), 0, 0, 1, 1}, {&__pyx_n_s_backprop_gradient_par_c, __pyx_k_backprop_gradient_par_c, sizeof(__pyx_k_backprop_gradient_par_c), 0, 0, 1, 1}, {&__pyx_n_s_base, __pyx_k_base, sizeof(__pyx_k_base), 0, 0, 1, 1}, {&__pyx_n_s_c, __pyx_k_c, sizeof(__pyx_k_c), 0, 0, 1, 1}, @@ -16901,6 +17551,7 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = { {&__pyx_n_s_d, __pyx_k_d, sizeof(__pyx_k_d), 0, 0, 1, 1}, {&__pyx_n_s_dL, __pyx_k_dL, sizeof(__pyx_k_dL), 0, 0, 1, 1}, {&__pyx_n_s_dL_dK, __pyx_k_dL_dK, sizeof(__pyx_k_dL_dK), 0, 0, 1, 1}, + {&__pyx_n_s_dot, __pyx_k_dot, sizeof(__pyx_k_dot), 0, 0, 1, 1}, {&__pyx_n_s_dtype_is_object, __pyx_k_dtype_is_object, sizeof(__pyx_k_dtype_is_object), 0, 0, 1, 1}, {&__pyx_n_s_empty, __pyx_k_empty, sizeof(__pyx_k_empty), 0, 0, 1, 1}, {&__pyx_n_s_enumerate, __pyx_k_enumerate, sizeof(__pyx_k_enumerate), 0, 0, 1, 1}, @@ -16912,14 +17563,15 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = { {&__pyx_n_s_fortran, __pyx_k_fortran, sizeof(__pyx_k_fortran), 0, 0, 1, 1}, {&__pyx_n_u_fortran, __pyx_k_fortran, sizeof(__pyx_k_fortran), 0, 1, 0, 1}, {&__pyx_kp_s_got_differing_extents_in_dimensi, __pyx_k_got_differing_extents_in_dimensi, sizeof(__pyx_k_got_differing_extents_in_dimensi), 0, 0, 1, 0}, - {&__pyx_kp_s_home_james_work_GPy_GPy_util_ch, __pyx_k_home_james_work_GPy_GPy_util_ch, sizeof(__pyx_k_home_james_work_GPy_GPy_util_ch), 0, 0, 1, 0}, {&__pyx_n_s_i, __pyx_k_i, sizeof(__pyx_k_i), 0, 0, 1, 1}, + {&__pyx_n_s_iN, __pyx_k_iN, sizeof(__pyx_k_iN), 0, 0, 1, 1}, {&__pyx_n_s_id, __pyx_k_id, sizeof(__pyx_k_id), 0, 0, 1, 1}, {&__pyx_n_s_import, __pyx_k_import, sizeof(__pyx_k_import), 0, 0, 1, 1}, {&__pyx_n_s_itemsize, __pyx_k_itemsize, sizeof(__pyx_k_itemsize), 0, 0, 1, 1}, {&__pyx_kp_s_itemsize_0_for_cython_array, __pyx_k_itemsize_0_for_cython_array, sizeof(__pyx_k_itemsize_0_for_cython_array), 0, 0, 1, 0}, {&__pyx_n_s_j, __pyx_k_j, sizeof(__pyx_k_j), 0, 0, 1, 1}, {&__pyx_n_s_k, __pyx_k_k, sizeof(__pyx_k_k), 0, 0, 1, 1}, + {&__pyx_n_s_kN, __pyx_k_kN, sizeof(__pyx_k_kN), 0, 0, 1, 1}, {&__pyx_n_s_m, __pyx_k_m, sizeof(__pyx_k_m), 0, 0, 1, 1}, {&__pyx_n_s_main, __pyx_k_main, sizeof(__pyx_k_main), 0, 0, 1, 1}, {&__pyx_n_s_memview, __pyx_k_memview, sizeof(__pyx_k_memview), 0, 0, 1, 1}, @@ -17179,7 +17831,7 @@ static int __Pyx_InitCachedConstants(void) { __pyx_tuple__18 = PyTuple_Pack(9, __pyx_n_s_flat, __pyx_n_s_M, __pyx_n_s_D, __pyx_n_s_N, __pyx_n_s_count, __pyx_n_s_ret, __pyx_n_s_d, __pyx_n_s_m, __pyx_n_s_mm); if (unlikely(!__pyx_tuple__18)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 11; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_tuple__18); __Pyx_GIVEREF(__pyx_tuple__18); - __pyx_codeobj__19 = (PyObject*)__Pyx_PyCode_New(2, 0, 9, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__18, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_home_james_work_GPy_GPy_util_ch, __pyx_n_s_flat_to_triang, 11, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__19)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 11; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_codeobj__19 = (PyObject*)__Pyx_PyCode_New(2, 0, 9, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__18, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_james_work_GPy_GPy_util_c, __pyx_n_s_flat_to_triang, 11, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__19)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 11; __pyx_clineno = __LINE__; goto __pyx_L1_error;} /* "GPy/util/choleskies_cython.pyx":31 * return ret @@ -17191,7 +17843,7 @@ static int __Pyx_InitCachedConstants(void) { __pyx_tuple__20 = PyTuple_Pack(9, __pyx_n_s_L, __pyx_n_s_D, __pyx_n_s_M, __pyx_n_s_N, __pyx_n_s_count, __pyx_n_s_flat, __pyx_n_s_d, __pyx_n_s_m, __pyx_n_s_mm); if (unlikely(!__pyx_tuple__20)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_tuple__20); __Pyx_GIVEREF(__pyx_tuple__20); - __pyx_codeobj__21 = (PyObject*)__Pyx_PyCode_New(1, 0, 9, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__20, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_home_james_work_GPy_GPy_util_ch, __pyx_n_s_triang_to_flat, 31, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__21)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_codeobj__21 = (PyObject*)__Pyx_PyCode_New(1, 0, 9, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__20, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_james_work_GPy_GPy_util_c, __pyx_n_s_triang_to_flat, 31, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__21)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 31; __pyx_clineno = __LINE__; goto __pyx_L1_error;} /* "GPy/util/choleskies_cython.pyx":47 * @@ -17203,7 +17855,7 @@ static int __Pyx_InitCachedConstants(void) { __pyx_tuple__22 = PyTuple_Pack(7, __pyx_n_s_dL, __pyx_n_s_L, __pyx_n_s_dL_dK, __pyx_n_s_N, __pyx_n_s_k, __pyx_n_s_j, __pyx_n_s_i); if (unlikely(!__pyx_tuple__22)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 47; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_tuple__22); __Pyx_GIVEREF(__pyx_tuple__22); - __pyx_codeobj__23 = (PyObject*)__Pyx_PyCode_New(2, 0, 7, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__22, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_home_james_work_GPy_GPy_util_ch, __pyx_n_s_backprop_gradient, 47, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__23)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 47; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_codeobj__23 = (PyObject*)__Pyx_PyCode_New(2, 0, 7, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__22, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_james_work_GPy_GPy_util_c, __pyx_n_s_backprop_gradient, 47, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__23)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 47; __pyx_clineno = __LINE__; goto __pyx_L1_error;} /* "GPy/util/choleskies_cython.pyx":62 * return dL_dK @@ -17215,19 +17867,31 @@ static int __Pyx_InitCachedConstants(void) { __pyx_tuple__24 = PyTuple_Pack(7, __pyx_n_s_dL, __pyx_n_s_L, __pyx_n_s_dL_dK, __pyx_n_s_N, __pyx_n_s_k, __pyx_n_s_j, __pyx_n_s_i); if (unlikely(!__pyx_tuple__24)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 62; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_tuple__24); __Pyx_GIVEREF(__pyx_tuple__24); - __pyx_codeobj__25 = (PyObject*)__Pyx_PyCode_New(2, 0, 7, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__24, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_home_james_work_GPy_GPy_util_ch, __pyx_n_s_backprop_gradient_par, 62, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__25)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 62; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_codeobj__25 = (PyObject*)__Pyx_PyCode_New(2, 0, 7, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__24, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_james_work_GPy_GPy_util_c, __pyx_n_s_backprop_gradient_par, 62, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__25)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 62; __pyx_clineno = __LINE__; goto __pyx_L1_error;} - /* "GPy/util/choleskies_cython.pyx":82 + /* "GPy/util/choleskies_cython.pyx":83 * void chol_backprop(int N, double* dL, double* L) * * def backprop_gradient_par_c(np.ndarray[double, ndim=2] dL, np.ndarray[double, ndim=2] L): # <<<<<<<<<<<<<< * cdef np.ndarray[double, ndim=2] dL_dK = np.tril(dL) # makes a copy, c-contig * cdef int N = L.shape[0] */ - __pyx_tuple__26 = PyTuple_Pack(4, __pyx_n_s_dL, __pyx_n_s_L, __pyx_n_s_dL_dK, __pyx_n_s_N); if (unlikely(!__pyx_tuple__26)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 82; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_tuple__26 = PyTuple_Pack(4, __pyx_n_s_dL, __pyx_n_s_L, __pyx_n_s_dL_dK, __pyx_n_s_N); if (unlikely(!__pyx_tuple__26)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_tuple__26); __Pyx_GIVEREF(__pyx_tuple__26); - __pyx_codeobj__27 = (PyObject*)__Pyx_PyCode_New(2, 0, 4, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__26, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_home_james_work_GPy_GPy_util_ch, __pyx_n_s_backprop_gradient_par_c, 82, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__27)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 82; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_codeobj__27 = (PyObject*)__Pyx_PyCode_New(2, 0, 4, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__26, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_james_work_GPy_GPy_util_c, __pyx_n_s_backprop_gradient_par_c, 83, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__27)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + + /* "GPy/util/choleskies_cython.pyx":92 + * + * # TODO: with the next release of cython, cimport scipy.linalg.cython_blas as blas, then blas the hell out of this. + * def backprop_gradient_par2(double[:,:] dL, double[:,:] L): # <<<<<<<<<<<<<< + * """ + * a very slow implementation, but the clearest I hope + */ + __pyx_tuple__28 = PyTuple_Pack(9, __pyx_n_s_dL, __pyx_n_s_L, __pyx_n_s_dL_dK, __pyx_n_s_N, __pyx_n_s_k, __pyx_n_s_j, __pyx_n_s_i, __pyx_n_s_iN, __pyx_n_s_kN); if (unlikely(!__pyx_tuple__28)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_tuple__28); + __Pyx_GIVEREF(__pyx_tuple__28); + __pyx_codeobj__29 = (PyObject*)__Pyx_PyCode_New(2, 0, 9, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__28, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_Users_james_work_GPy_GPy_util_c, __pyx_n_s_backprop_gradient_par2, 92, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__29)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L1_error;} /* "View.MemoryView":276 * return self.name @@ -17236,9 +17900,9 @@ static int __Pyx_InitCachedConstants(void) { * cdef strided = Enum("") # default * cdef indirect = Enum("") */ - __pyx_tuple__28 = PyTuple_Pack(1, __pyx_kp_s_strided_and_direct_or_indirect); if (unlikely(!__pyx_tuple__28)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 276; __pyx_clineno = __LINE__; goto __pyx_L1_error;} - __Pyx_GOTREF(__pyx_tuple__28); - __Pyx_GIVEREF(__pyx_tuple__28); + __pyx_tuple__30 = PyTuple_Pack(1, __pyx_kp_s_strided_and_direct_or_indirect); if (unlikely(!__pyx_tuple__30)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 276; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_tuple__30); + __Pyx_GIVEREF(__pyx_tuple__30); /* "View.MemoryView":277 * @@ -17247,9 +17911,9 @@ static int __Pyx_InitCachedConstants(void) { * cdef indirect = Enum("") * */ - __pyx_tuple__29 = PyTuple_Pack(1, __pyx_kp_s_strided_and_direct); if (unlikely(!__pyx_tuple__29)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 277; __pyx_clineno = __LINE__; goto __pyx_L1_error;} - __Pyx_GOTREF(__pyx_tuple__29); - __Pyx_GIVEREF(__pyx_tuple__29); + __pyx_tuple__31 = PyTuple_Pack(1, __pyx_kp_s_strided_and_direct); if (unlikely(!__pyx_tuple__31)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 277; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_tuple__31); + __Pyx_GIVEREF(__pyx_tuple__31); /* "View.MemoryView":278 * cdef generic = Enum("") @@ -17258,9 +17922,9 @@ static int __Pyx_InitCachedConstants(void) { * * */ - __pyx_tuple__30 = PyTuple_Pack(1, __pyx_kp_s_strided_and_indirect); if (unlikely(!__pyx_tuple__30)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 278; __pyx_clineno = __LINE__; goto __pyx_L1_error;} - __Pyx_GOTREF(__pyx_tuple__30); - __Pyx_GIVEREF(__pyx_tuple__30); + __pyx_tuple__32 = PyTuple_Pack(1, __pyx_kp_s_strided_and_indirect); if (unlikely(!__pyx_tuple__32)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 278; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_tuple__32); + __Pyx_GIVEREF(__pyx_tuple__32); /* "View.MemoryView":281 * @@ -17269,9 +17933,9 @@ static int __Pyx_InitCachedConstants(void) { * cdef indirect_contiguous = Enum("") * */ - __pyx_tuple__31 = PyTuple_Pack(1, __pyx_kp_s_contiguous_and_direct); if (unlikely(!__pyx_tuple__31)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 281; __pyx_clineno = __LINE__; goto __pyx_L1_error;} - __Pyx_GOTREF(__pyx_tuple__31); - __Pyx_GIVEREF(__pyx_tuple__31); + __pyx_tuple__33 = PyTuple_Pack(1, __pyx_kp_s_contiguous_and_direct); if (unlikely(!__pyx_tuple__33)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 281; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_tuple__33); + __Pyx_GIVEREF(__pyx_tuple__33); /* "View.MemoryView":282 * @@ -17280,9 +17944,9 @@ static int __Pyx_InitCachedConstants(void) { * * */ - __pyx_tuple__32 = PyTuple_Pack(1, __pyx_kp_s_contiguous_and_indirect); if (unlikely(!__pyx_tuple__32)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 282; __pyx_clineno = __LINE__; goto __pyx_L1_error;} - __Pyx_GOTREF(__pyx_tuple__32); - __Pyx_GIVEREF(__pyx_tuple__32); + __pyx_tuple__34 = PyTuple_Pack(1, __pyx_kp_s_contiguous_and_indirect); if (unlikely(!__pyx_tuple__34)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 282; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_tuple__34); + __Pyx_GIVEREF(__pyx_tuple__34); __Pyx_RefNannyFinishContext(); return 0; __pyx_L1_error:; @@ -17496,16 +18160,28 @@ PyMODINIT_FUNC PyInit_choleskies_cython(void) if (PyDict_SetItem(__pyx_d, __pyx_n_s_backprop_gradient_par, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 62; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; - /* "GPy/util/choleskies_cython.pyx":82 + /* "GPy/util/choleskies_cython.pyx":83 * void chol_backprop(int N, double* dL, double* L) * * def backprop_gradient_par_c(np.ndarray[double, ndim=2] dL, np.ndarray[double, ndim=2] L): # <<<<<<<<<<<<<< * cdef np.ndarray[double, ndim=2] dL_dK = np.tril(dL) # makes a copy, c-contig * cdef int N = L.shape[0] */ - __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_3GPy_4util_17choleskies_cython_9backprop_gradient_par_c, NULL, __pyx_n_s_GPy_util_choleskies_cython); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 82; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_3GPy_4util_17choleskies_cython_9backprop_gradient_par_c, NULL, __pyx_n_s_GPy_util_choleskies_cython); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_t_1); - if (PyDict_SetItem(__pyx_d, __pyx_n_s_backprop_gradient_par_c, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 82; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + if (PyDict_SetItem(__pyx_d, __pyx_n_s_backprop_gradient_par_c, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; + + /* "GPy/util/choleskies_cython.pyx":92 + * + * # TODO: with the next release of cython, cimport scipy.linalg.cython_blas as blas, then blas the hell out of this. + * def backprop_gradient_par2(double[:,:] dL, double[:,:] L): # <<<<<<<<<<<<<< + * """ + * a very slow implementation, but the clearest I hope + */ + __pyx_t_1 = PyCFunction_NewEx(&__pyx_mdef_3GPy_4util_17choleskies_cython_11backprop_gradient_par2, NULL, __pyx_n_s_GPy_util_choleskies_cython); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __Pyx_GOTREF(__pyx_t_1); + if (PyDict_SetItem(__pyx_d, __pyx_n_s_backprop_gradient_par2, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 92; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; /* "GPy/util/choleskies_cython.pyx":1 @@ -17538,7 +18214,7 @@ PyMODINIT_FUNC PyInit_choleskies_cython(void) * cdef strided = Enum("") # default * cdef indirect = Enum("") */ - __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)((PyObject *)__pyx_MemviewEnum_type)), __pyx_tuple__28, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 276; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)((PyObject *)__pyx_MemviewEnum_type)), __pyx_tuple__30, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 276; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_t_1); __Pyx_XGOTREF(generic); __Pyx_DECREF_SET(generic, __pyx_t_1); @@ -17552,7 +18228,7 @@ PyMODINIT_FUNC PyInit_choleskies_cython(void) * cdef indirect = Enum("") * */ - __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)((PyObject *)__pyx_MemviewEnum_type)), __pyx_tuple__29, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 277; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)((PyObject *)__pyx_MemviewEnum_type)), __pyx_tuple__31, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 277; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_t_1); __Pyx_XGOTREF(strided); __Pyx_DECREF_SET(strided, __pyx_t_1); @@ -17566,7 +18242,7 @@ PyMODINIT_FUNC PyInit_choleskies_cython(void) * * */ - __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)((PyObject *)__pyx_MemviewEnum_type)), __pyx_tuple__30, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 278; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)((PyObject *)__pyx_MemviewEnum_type)), __pyx_tuple__32, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 278; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_t_1); __Pyx_XGOTREF(indirect); __Pyx_DECREF_SET(indirect, __pyx_t_1); @@ -17580,7 +18256,7 @@ PyMODINIT_FUNC PyInit_choleskies_cython(void) * cdef indirect_contiguous = Enum("") * */ - __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)((PyObject *)__pyx_MemviewEnum_type)), __pyx_tuple__31, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 281; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)((PyObject *)__pyx_MemviewEnum_type)), __pyx_tuple__33, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 281; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_t_1); __Pyx_XGOTREF(contiguous); __Pyx_DECREF_SET(contiguous, __pyx_t_1); @@ -17594,7 +18270,7 @@ PyMODINIT_FUNC PyInit_choleskies_cython(void) * * */ - __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)((PyObject *)__pyx_MemviewEnum_type)), __pyx_tuple__32, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 282; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_t_1 = __Pyx_PyObject_Call(((PyObject *)((PyObject *)__pyx_MemviewEnum_type)), __pyx_tuple__34, NULL); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[2]; __pyx_lineno = 282; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_GOTREF(__pyx_t_1); __Pyx_XGOTREF(indirect_contiguous); __Pyx_DECREF_SET(indirect_contiguous, __pyx_t_1); diff --git a/GPy/util/choleskies_cython.pyx b/GPy/util/choleskies_cython.pyx index 50fa6709..bd6069b9 100644 --- a/GPy/util/choleskies_cython.pyx +++ b/GPy/util/choleskies_cython.pyx @@ -76,12 +76,36 @@ def backprop_gradient_par(double[:,:] dL, double[:,:] L): dL_dK[k, k] /= (2. * L[k, k]) return dL_dK -cdef extern from "cholesky_backprop.h": +#here's a pure C version... +cdef extern from "cholesky_backprop.h" nogil: void chol_backprop(int N, double* dL, double* L) def backprop_gradient_par_c(np.ndarray[double, ndim=2] dL, np.ndarray[double, ndim=2] L): cdef np.ndarray[double, ndim=2] dL_dK = np.tril(dL) # makes a copy, c-contig cdef int N = L.shape[0] - chol_backprop(N, dL_dK.data, L.data) + with nogil: + chol_backprop(N, dL_dK.data, L.data) return dL_dK + +# TODO: with the next release of cython, cimport scipy.linalg.cython_blas as blas, then blas the hell out of this. +def backprop_gradient_par2(double[:,:] dL, double[:,:] L): + """ + a very slow implementation, but the clearest I hope + """ + cdef double[:,:] dL_dK = np.tril(dL).copy() + cdef int N = L.shape[0] + cdef int k, j, i, iN, kN + for k in range(N - 1, -1, -1): + #pragma this loop: + for i in range(k+1, N): + dL_dK[i,+k] -= np.dot(dL_dK[i,k+1:], L[k+1:,k]) + dL_dK[i,+k] -= np.dot(dL_dK[i:,i], L[i:,k]) + + for i in range(k+1, N): + dL_dK[i, k] /= L[k, k]; + dL_dK[k, k] -= L[i, k] * dL_dK[i, k]; + + dL_dK[k, k] /= (2. * L[k, k]); + return np.asarray(dL_dK) + diff --git a/GPy/util/cholesky_backprop.c b/GPy/util/cholesky_backprop.c index e3919abb..975e63cc 100644 --- a/GPy/util/cholesky_backprop.c +++ b/GPy/util/cholesky_backprop.c @@ -1,22 +1,32 @@ +#include +double mydot(int n, double* a, int stride_a, double* b, int stride_b){ + double ret = 0; + for(int i=0; i(-1);k--){ - #pragma omp parallel for private(i,j) - for(i=k+1;i(-1);k--){ + kN = k*N; + #pragma omp parallel for private(iN) + for(int i=k+1; i