diff --git a/GPy/util/choleskies_cython.c b/GPy/util/choleskies_cython.c index db7d1aed..21492a88 100644 --- a/GPy/util/choleskies_cython.c +++ b/GPy/util/choleskies_cython.c @@ -1470,7 +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 PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_10backprop_gradient_par_c_old(CYTHON_UNUSED PyObject *__pyx_self, PyArrayObject *__pyx_v_dL, PyArrayObject *__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 */ @@ -1534,12 +1534,9 @@ 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"; @@ -1572,7 +1569,6 @@ 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"; @@ -1595,7 +1591,6 @@ 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[] = ""; @@ -1606,11 +1601,12 @@ static char __pyx_k_Invalid_shape_in_axis_d_d[] = "Invalid shape in axis %d: %d. static char __pyx_k_GPy_util_choleskies_cython[] = "GPy.util.choleskies_cython"; static char __pyx_k_Index_out_of_bounds_axis_d[] = "Index out of bounds (axis %d)"; static char __pyx_k_Step_may_not_be_zero_axis_d[] = "Step may not be zero (axis %d)"; +static char __pyx_k_backprop_gradient_par_c_old[] = "backprop_gradient_par_c_old"; static char __pyx_k_itemsize_0_for_cython_array[] = "itemsize <= 0 for cython.array"; 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_Users_james_work_GPy_GPy_util_c[] = "/Users/james/work/GPy/GPy/util/choleskies_cython.pyx"; +static char __pyx_k_home_james_work_GPy_GPy_util_ch[] = "/home/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"; @@ -1652,14 +1648,12 @@ 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_backprop_gradient_par_c_old; static PyObject *__pyx_n_s_base; static PyObject *__pyx_n_s_c; static PyObject *__pyx_n_u_c; @@ -1671,7 +1665,6 @@ 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; @@ -1683,15 +1676,14 @@ 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; @@ -3056,7 +3048,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_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) + #pragma omp parallel private(__pyx_t_17, __pyx_t_14, __pyx_t_20, __pyx_t_26, __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_19, __pyx_t_27, __pyx_t_15, __pyx_t_13, __pyx_t_12, __pyx_t_18, __pyx_t_24) #endif /* _OPENMP */ { @@ -3516,7 +3508,7 @@ static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_8backprop_gradient_par_ * chol_backprop(N, dL_dK.data, L.data) * return dL_dK # <<<<<<<<<<<<<< * - * + * cdef extern from "cholesky_backprop.h" nogil: */ __Pyx_XDECREF(__pyx_r); __Pyx_INCREF(((PyObject *)__pyx_v_dL_dK)); @@ -3557,27 +3549,26 @@ static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_8backprop_gradient_par_ return __pyx_r; } -/* "GPy/util/choleskies_cython.pyx":92 +/* "GPy/util/choleskies_cython.pyx":93 + * void old_chol_backprop(int N, double* dL, double* L) * - * # 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 + * def backprop_gradient_par_c_old(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] */ /* 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 } }; +static PyObject *__pyx_pw_3GPy_4util_17choleskies_cython_11backprop_gradient_par_c_old(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/ +static PyMethodDef __pyx_mdef_3GPy_4util_17choleskies_cython_11backprop_gradient_par_c_old = {"backprop_gradient_par_c_old", (PyCFunction)__pyx_pw_3GPy_4util_17choleskies_cython_11backprop_gradient_par_c_old, METH_VARARGS|METH_KEYWORDS, 0}; +static PyObject *__pyx_pw_3GPy_4util_17choleskies_cython_11backprop_gradient_par_c_old(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) { + PyArrayObject *__pyx_v_dL = 0; + PyArrayObject *__pyx_v_L = 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); + __Pyx_RefNannySetupContext("backprop_gradient_par_c_old (wrapper)", 0); { static PyObject **__pyx_pyargnames[] = {&__pyx_n_s_dL,&__pyx_n_s_L,0}; PyObject* values[2] = {0,0}; @@ -3598,11 +3589,11 @@ static PyObject *__pyx_pw_3GPy_4util_17choleskies_cython_11backprop_gradient_par 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;} + __Pyx_RaiseArgtupleInvalid("backprop_gradient_par_c_old", 1, 2, 2, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __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;} + if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, pos_args, "backprop_gradient_par_c_old") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __pyx_clineno = __LINE__; goto __pyx_L3_error;} } } else if (PyTuple_GET_SIZE(__pyx_args) != 2) { goto __pyx_L5_argtuple_error; @@ -3610,527 +3601,193 @@ static PyObject *__pyx_pw_3GPy_4util_17choleskies_cython_11backprop_gradient_par 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;} + __pyx_v_dL = ((PyArrayObject *)values[0]); + __pyx_v_L = ((PyArrayObject *)values[1]); } 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_RaiseArgtupleInvalid("backprop_gradient_par_c_old", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __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_AddTraceback("GPy.util.choleskies_cython.backprop_gradient_par_c_old", __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); + if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_dL), __pyx_ptype_5numpy_ndarray, 1, "dL", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __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 = 93; __pyx_clineno = __LINE__; goto __pyx_L1_error;} + __pyx_r = __pyx_pf_3GPy_4util_17choleskies_cython_10backprop_gradient_par_c_old(__pyx_self, __pyx_v_dL, __pyx_v_L); /* function exit code */ + goto __pyx_L0; + __pyx_L1_error:; + __pyx_r = NULL; + __pyx_L0:; __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 } }; +static PyObject *__pyx_pf_3GPy_4util_17choleskies_cython_10backprop_gradient_par_c_old(CYTHON_UNUSED PyObject *__pyx_self, PyArrayObject *__pyx_v_dL, PyArrayObject *__pyx_v_L) { + PyArrayObject *__pyx_v_dL_dK = 0; int __pyx_v_N; - int __pyx_v_k; - int __pyx_v_i; + __Pyx_LocalBuf_ND __pyx_pybuffernd_L; + __Pyx_Buffer __pyx_pybuffer_L; + __Pyx_LocalBuf_ND __pyx_pybuffernd_dL; + __Pyx_Buffer __pyx_pybuffer_dL; + __Pyx_LocalBuf_ND __pyx_pybuffernd_dL_dK; + __Pyx_Buffer __pyx_pybuffer_dL_dK; 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; + PyArrayObject *__pyx_t_5 = NULL; int __pyx_lineno = 0; const char *__pyx_filename = NULL; int __pyx_clineno = 0; - __Pyx_RefNannySetupContext("backprop_gradient_par2", 0); + __Pyx_RefNannySetupContext("backprop_gradient_par_c_old", 0); + __pyx_pybuffer_dL_dK.pybuffer.buf = NULL; + __pyx_pybuffer_dL_dK.refcount = 0; + __pyx_pybuffernd_dL_dK.data = NULL; + __pyx_pybuffernd_dL_dK.rcbuffer = &__pyx_pybuffer_dL_dK; + __pyx_pybuffer_dL.pybuffer.buf = NULL; + __pyx_pybuffer_dL.refcount = 0; + __pyx_pybuffernd_dL.data = NULL; + __pyx_pybuffernd_dL.rcbuffer = &__pyx_pybuffer_dL; + __pyx_pybuffer_L.pybuffer.buf = NULL; + __pyx_pybuffer_L.refcount = 0; + __pyx_pybuffernd_L.data = NULL; + __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 = 93; __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 = 93; __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":96 - * a very slow implementation, but the clearest I hope - * """ - * cdef double[:,:] dL_dK = np.tril(dL).copy() # <<<<<<<<<<<<<< + /* "GPy/util/choleskies_cython.pyx":94 + * + * def backprop_gradient_par_c_old(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] - * cdef int k, j, i, iN, kN + * with nogil: */ - __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_t_2 = __Pyx_GetModuleGlobalName(__pyx_n_s_np); if (unlikely(!__pyx_t_2)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __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 = 94; __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)) { + __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_6); + __Pyx_INCREF(__pyx_t_2); __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; + 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 = 94; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __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_t_4 = PyTuple_New(1+1); if (unlikely(!__pyx_t_4)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __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 = 94; __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_4); __pyx_t_4 = 0; } __Pyx_DECREF(__pyx_t_3); __pyx_t_3 = 0; - __pyx_r = __pyx_t_1; + if (!(likely(((__pyx_t_1) == Py_None) || likely(__Pyx_TypeTest(__pyx_t_1, __pyx_ptype_5numpy_ndarray))))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 94; __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 = 94; __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]; + } + } + __pyx_t_5 = 0; + __pyx_v_dL_dK = ((PyArrayObject *)__pyx_t_1); __pyx_t_1 = 0; + + /* "GPy/util/choleskies_cython.pyx":95 + * def backprop_gradient_par_c_old(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] # <<<<<<<<<<<<<< + * with nogil: + * old_chol_backprop(N, dL_dK.data, L.data) + */ + __pyx_v_N = (__pyx_v_L->dimensions[0]); + + /* "GPy/util/choleskies_cython.pyx":96 + * cdef np.ndarray[double, ndim=2] dL_dK = np.tril(dL) # makes a copy, c-contig + * cdef int N = L.shape[0] + * with nogil: # <<<<<<<<<<<<<< + * old_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":97 + * cdef int N = L.shape[0] + * with nogil: + * old_chol_backprop(N, dL_dK.data, L.data) # <<<<<<<<<<<<<< + * return dL_dK + * + */ + old_chol_backprop(__pyx_v_N, ((double *)__pyx_v_dL_dK->data), ((double *)__pyx_v_L->data)); + } + + /* "GPy/util/choleskies_cython.pyx":96 + * cdef np.ndarray[double, ndim=2] dL_dK = np.tril(dL) # makes a copy, c-contig + * cdef int N = L.shape[0] + * with nogil: # <<<<<<<<<<<<<< + * old_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":98 + * with nogil: + * old_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":92 + /* "GPy/util/choleskies_cython.pyx":93 + * void old_chol_backprop(int N, double* dL, double* L) * - * # 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 + * def backprop_gradient_par_c_old(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] */ /* function exit code */ @@ -4139,19 +3796,21 @@ __pyx_t_6 = __pyx_memoryview_fromslice(__pyx_t_19, 1, (PyObject *(*)(char *)) __ __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); + { PyObject *__pyx_type, *__pyx_value, *__pyx_tb; + __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb); + __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_L.rcbuffer->pybuffer); + __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_dL.rcbuffer->pybuffer); + __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_dL_dK.rcbuffer->pybuffer); + __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);} + __Pyx_AddTraceback("GPy.util.choleskies_cython.backprop_gradient_par_c_old", __pyx_clineno, __pyx_lineno, __pyx_filename); __pyx_r = NULL; + goto __pyx_L2; __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_SafeReleaseBuffer(&__pyx_pybuffernd_L.rcbuffer->pybuffer); + __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_dL.rcbuffer->pybuffer); + __Pyx_SafeReleaseBuffer(&__pyx_pybuffernd_dL_dK.rcbuffer->pybuffer); + __pyx_L2:; + __Pyx_XDECREF((PyObject *)__pyx_v_dL_dK); __Pyx_XGIVEREF(__pyx_r); __Pyx_RefNannyFinishContext(); return __pyx_r; @@ -17532,14 +17191,12 @@ 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_backprop_gradient_par_c_old, __pyx_k_backprop_gradient_par_c_old, sizeof(__pyx_k_backprop_gradient_par_c_old), 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}, {&__pyx_n_u_c, __pyx_k_c, sizeof(__pyx_k_c), 0, 1, 0, 1}, @@ -17551,7 +17208,6 @@ 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}, @@ -17563,15 +17219,14 @@ 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}, @@ -17831,7 +17486,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_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;} + __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;} /* "GPy/util/choleskies_cython.pyx":31 * return ret @@ -17843,7 +17498,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_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;} + __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;} /* "GPy/util/choleskies_cython.pyx":47 * @@ -17855,7 +17510,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_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;} + __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;} /* "GPy/util/choleskies_cython.pyx":62 * return dL_dK @@ -17867,7 +17522,7 @@ 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_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;} + __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;} /* "GPy/util/choleskies_cython.pyx":83 * void chol_backprop(int N, double* dL, double* L) @@ -17879,19 +17534,19 @@ static int __Pyx_InitCachedConstants(void) { __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_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;} + __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, 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 + /* "GPy/util/choleskies_cython.pyx":93 + * void old_chol_backprop(int N, double* dL, double* L) * - * # 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 + * def backprop_gradient_par_c_old(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__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_tuple__28 = PyTuple_Pack(4, __pyx_n_s_dL, __pyx_n_s_L, __pyx_n_s_dL_dK, __pyx_n_s_N); if (unlikely(!__pyx_tuple__28)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __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;} + __pyx_codeobj__29 = (PyObject*)__Pyx_PyCode_New(2, 0, 4, 0, 0, __pyx_empty_bytes, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_tuple__28, __pyx_empty_tuple, __pyx_empty_tuple, __pyx_kp_s_home_james_work_GPy_GPy_util_ch, __pyx_n_s_backprop_gradient_par_c_old, 93, __pyx_empty_bytes); if (unlikely(!__pyx_codeobj__29)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __pyx_clineno = __LINE__; goto __pyx_L1_error;} /* "View.MemoryView":276 * return self.name @@ -18172,16 +17827,16 @@ PyMODINIT_FUNC PyInit_choleskies_cython(void) 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 + /* "GPy/util/choleskies_cython.pyx":93 + * void old_chol_backprop(int N, double* dL, double* L) * - * # 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 + * def backprop_gradient_par_c_old(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_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_t_1 = PyCFunction_NewEx(&__pyx_mdef_3GPy_4util_17choleskies_cython_11backprop_gradient_par_c_old, NULL, __pyx_n_s_GPy_util_choleskies_cython); if (unlikely(!__pyx_t_1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __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;} + if (PyDict_SetItem(__pyx_d, __pyx_n_s_backprop_gradient_par_c_old, __pyx_t_1) < 0) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 93; __pyx_clineno = __LINE__; goto __pyx_L1_error;} __Pyx_DECREF(__pyx_t_1); __pyx_t_1 = 0; /* "GPy/util/choleskies_cython.pyx":1 diff --git a/GPy/util/choleskies_cython.pyx b/GPy/util/choleskies_cython.pyx index bd6069b9..d85b19cf 100644 --- a/GPy/util/choleskies_cython.pyx +++ b/GPy/util/choleskies_cython.pyx @@ -87,25 +87,21 @@ def backprop_gradient_par_c(np.ndarray[double, ndim=2] dL, np.ndarray[double, nd chol_backprop(N, dL_dK.data, L.data) return dL_dK +cdef extern from "cholesky_backprop.h" nogil: + void old_chol_backprop(int N, double* dL, double* L) -# 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() +def backprop_gradient_par_c_old(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] - 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]) + with nogil: + old_chol_backprop(N, dL_dK.data, L.data) + return dL_dK + + + + + + - 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 975e63cc..0a4a2f03 100644 --- a/GPy/util/cholesky_backprop.c +++ b/GPy/util/cholesky_backprop.c @@ -1,4 +1,25 @@ -#include +#include +void chol_backprop(int N, double* dL, double* L){ + //at the input to this fn, dL is df_dL. after this fn is complet, dL is df_dK + int i,k; + + dL[N*N - 1] /= (2. * L[N*N - 1]); + for(k=N-2;k>(-1);k--){ + cblas_dsymv(CblasRowMajor, CblasLower, + N-k-1, -1, + &dL[(N*(k+1) + k+1)],N, + &L[k*N+k+1],1, + 1, &dL[N*(k+1)+k], N); + for(i=0;i<(N-k-1); i++){ + dL[N*(k+1+i)+k] -= dL[N*(k+1)+k+i*(N+1)+1] * L[k*N+k+1+i]; + } + + cblas_dscal(N-k-1, 1.0/L[k*N+k], &dL[(k+1)*N+k], N); + dL[k*N + k] -= cblas_ddot(N-k-1, &dL[(k+1)*N+k], N, &L[k*N+k+1], 1); + dL[k*N + k] /= (2.0 * L[k*N + k]); + } +} + double mydot(int n, double* a, int stride_a, double* b, int stride_b){ double ret = 0; for(int i=0; i(-1);k--){ + int iN, kN,i,j,k; + dL[N*N-1] /= (2. * U[N*N-1]); + for(k=N-2;k>(-1);k--){ kN = k*N; - #pragma omp parallel for private(iN) - for(int i=k+1; i + +void dsymv(int N, double*A, double*b, double*y); double mydot(int n, double* a, int stride_a, double* b, int stride_b); void chol_backprop(int N, double* dL, double* L); diff --git a/setup.py b/setup.py index e12cbaa4..8025cd5d 100644 --- a/setup.py +++ b/setup.py @@ -22,7 +22,7 @@ ext_mods = [Extension(name='GPy.kern._src.stationary_cython', Extension(name='GPy.util.choleskies_cython', sources=['GPy/util/choleskies_cython.c', 'GPy/util/cholesky_backprop.c'], include_dirs=[np.get_include()], - extra_link_args = ['-lgomp'], + extra_link_args = ['-lgomp', '-lblas'], extra_compile_args=compile_flags+['-std=c99']), Extension(name='GPy.util.linalg_cython', sources=['GPy/util/linalg_cython.c'],