lots of work on cython choleskies

This commit is contained in:
James Hensman 2015-06-11 18:00:47 +01:00
parent 93076df259
commit 60fabe41e2
5 changed files with 243 additions and 570 deletions

View file

@ -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[] = "<strided and indirect>";
static char __pyx_k_backprop_gradient_par[] = "backprop_gradient_par";
static char __pyx_k_contiguous_and_direct[] = "<contiguous and direct>";
static char __pyx_k_MemoryView_of_r_object[] = "<MemoryView of %r object>";
static char __pyx_k_backprop_gradient_par2[] = "backprop_gradient_par2";
static char __pyx_k_MemoryView_of_r_at_0x_x[] = "<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[] = "<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[] = "<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, <double*> dL_dK.data, <double*> 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, <double*> dL_dK.data, <double*> 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, <double*> dL_dK.data, <double*> 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, <double*> dL_dK.data, <double*> 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, <double*> dL_dK.data, <double*> 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, <double*> dL_dK.data, <double*> 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

View file

@ -87,25 +87,21 @@ def backprop_gradient_par_c(np.ndarray[double, ndim=2] dL, np.ndarray[double, nd
chol_backprop(N, <double*> dL_dK.data, <double*> 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, <double*> dL_dK.data, <double*> 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)

View file

@ -1,4 +1,25 @@
#include <omp.h>
#include <cblas.h>
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<n; i++){
@ -6,26 +27,24 @@ double mydot(int n, double* a, int stride_a, double* b, int stride_b){
}
return ret;
}
void chol_backprop(int N, double* dL, double* U){
void old_chol_backprop(int N, double* dL, double* U){
//at the input to this fn, dL is df_dL. after this fn is complet, dL is df_dK
int iN, kN;
for(int k=N-1;k>(-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<N; i++){
#pragma omp parallel for private(i,iN)
for(i=k+1; i<N; i++){
iN = i*N;
dL[iN+k] -= mydot(i-k, &dL[iN+k+1], 1, &U[kN+k+1], 1);
dL[iN+k] -= mydot(N-i, &dL[iN+i], N, &U[kN+i], 1);
}
for(int i=(k + 1); i<N; i++){
for(i=(k + 1); i<N; i++){
iN = i*N;
dL[iN + k] /= U[kN + k];
dL[kN + k] -= U[kN + i] * dL[iN + k];
}
dL[kN + k] /= (2. * U[kN + k]);
}
}

View file

@ -1,2 +1,5 @@
#include <cblas.h>
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);

View file

@ -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'],