Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with HTTPS or Subversion.

Download ZIP
Browse files

More on LARS performance: triangular solving and cholesky deletes.

seems to be working ...
  • Loading branch information...
commit a79553e56f014819f407d71e73893fc991082aeb 1 parent c463cea
@fabianp fabianp authored
View
6 scikits/learn/glm/lars.py
@@ -221,12 +221,14 @@ def lars_path(X, y, Gram=None, max_iter=None, alpha_min=0,
beta[n_iter, active] = beta[n_iter - 1, active] + gamma_ * b
if drop:
+ arrayfuncs.cholesky_delete (L[:n_pred, :n_pred], idx)
n_pred -= 1
drop_idx = active.pop (idx)
unactive.append(drop_idx)
active_mask[drop_idx] = False
- Xa = Xt[active] # duplicate
- L[:n_pred, :n_pred] = linalg.cholesky(np.dot(Xa, Xa.T), lower=True)
+ # pdb.set_trace()
+ # Xa = Xt[active] # duplicate
+ # L[:n_pred, :n_pred] = linalg.cholesky(np.dot(Xa, Xa.T), lower=True)
sign_active = np.delete (sign_active, idx) # do an append to maintain size
sign_active = np.append (sign_active, 0.)
# should be done using cholesky deletes
View
224 scikits/learn/utils/arrayfuncs.c
@@ -1,4 +1,4 @@
-/* Generated by Cython 0.12.1 on Tue Sep 14 22:01:04 2010 */
+/* Generated by Cython 0.12.1 on Tue Sep 14 22:45:32 2010 */
#define PY_SSIZE_T_CLEAN
#include "Python.h"
@@ -158,6 +158,7 @@
#include "numpy/arrayobject.h"
#include "numpy/ufuncobject.h"
#include "cblas.h"
+#include "src/cholesky_delete.c"
#ifndef CYTHON_INLINE
#if defined(__GNUC__)
@@ -661,8 +662,8 @@ static char __pyx_k_4[] = "unknown dtype code in numpy.pxd (%d)";
static char __pyx_k_5[] = "Format string allocated too short, see comment in numpy.pxd";
static char __pyx_k_6[] = "Format string allocated too short.";
static char __pyx_k_7[] = "\nSmall collection of auxiliary functions that operate on arrays\n\n";
-static char __pyx_k_8[] = "dot_over (line 39)";
-static char __pyx_k_9[] = "solve_triangular (line 62)";
+static char __pyx_k_8[] = "dot_over (line 44)";
+static char __pyx_k_9[] = "solve_triangular (line 67)";
static char __pyx_k__B[] = "B";
static char __pyx_k__H[] = "H";
static char __pyx_k__I[] = "I";
@@ -698,6 +699,7 @@ static char __pyx_k__range[] = "range";
static char __pyx_k__shape[] = "shape";
static char __pyx_k__fields[] = "fields";
static char __pyx_k__format[] = "format";
+static char __pyx_k__go_out[] = "go_out";
static char __pyx_k__strides[] = "strides";
static char __pyx_k____main__[] = "__main__";
static char __pyx_k____test__[] = "__test__";
@@ -718,6 +720,7 @@ static PyObject *__pyx_kp_u_5;
static PyObject *__pyx_kp_u_6;
static PyObject *__pyx_kp_u_8;
static PyObject *__pyx_kp_u_9;
+static PyObject *__pyx_n_s__L;
static PyObject *__pyx_n_s__RuntimeError;
static PyObject *__pyx_n_s__ValueError;
static PyObject *__pyx_n_s__X;
@@ -731,6 +734,7 @@ static PyObject *__pyx_n_s__descr;
static PyObject *__pyx_n_s__dot_over;
static PyObject *__pyx_n_s__fields;
static PyObject *__pyx_n_s__format;
+static PyObject *__pyx_n_s__go_out;
static PyObject *__pyx_n_s__itemsize;
static PyObject *__pyx_n_s__mask;
static PyObject *__pyx_n_s__names;
@@ -750,7 +754,7 @@ static PyObject *__pyx_n_s__val;
static PyObject *__pyx_n_s__y;
static PyObject *__pyx_int_15;
-/* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":39
+/* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":44
*
*
* def dot_over ( # <<<<<<<<<<<<<<
@@ -819,34 +823,34 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over(PyObject
values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__y);
if (likely(values[1])) kw_args--;
else {
- __Pyx_RaiseArgtupleInvalid("dot_over", 1, 5, 5, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 39; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __Pyx_RaiseArgtupleInvalid("dot_over", 1, 5, 5, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 44; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
}
case 2:
values[2] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__mask);
if (likely(values[2])) kw_args--;
else {
- __Pyx_RaiseArgtupleInvalid("dot_over", 1, 5, 5, 2); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 39; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __Pyx_RaiseArgtupleInvalid("dot_over", 1, 5, 5, 2); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 44; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
}
case 3:
values[3] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__val);
if (likely(values[3])) kw_args--;
else {
- __Pyx_RaiseArgtupleInvalid("dot_over", 1, 5, 5, 3); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 39; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __Pyx_RaiseArgtupleInvalid("dot_over", 1, 5, 5, 3); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 44; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
}
case 4:
values[4] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__out);
if (likely(values[4])) kw_args--;
else {
- __Pyx_RaiseArgtupleInvalid("dot_over", 1, 5, 5, 4); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 39; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __Pyx_RaiseArgtupleInvalid("dot_over", 1, 5, 5, 4); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 44; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
}
}
if (unlikely(kw_args > 0)) {
- if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, PyTuple_GET_SIZE(__pyx_args), "dot_over") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 39; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, PyTuple_GET_SIZE(__pyx_args), "dot_over") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 44; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
}
__pyx_v_X = ((PyArrayObject *)values[0]);
__pyx_v_y = ((PyArrayObject *)values[1]);
__pyx_v_mask = ((PyArrayObject *)values[2]);
- __pyx_v_val = __Pyx_PyInt_from_py_npy_uint8(values[3]); if (unlikely((__pyx_v_val == (npy_uint8)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 43; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __pyx_v_val = __Pyx_PyInt_from_py_npy_uint8(values[3]); if (unlikely((__pyx_v_val == (npy_uint8)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 48; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
__pyx_v_out = ((PyArrayObject *)values[4]);
} else if (PyTuple_GET_SIZE(__pyx_args) != 5) {
goto __pyx_L5_argtuple_error;
@@ -854,12 +858,12 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over(PyObject
__pyx_v_X = ((PyArrayObject *)PyTuple_GET_ITEM(__pyx_args, 0));
__pyx_v_y = ((PyArrayObject *)PyTuple_GET_ITEM(__pyx_args, 1));
__pyx_v_mask = ((PyArrayObject *)PyTuple_GET_ITEM(__pyx_args, 2));
- __pyx_v_val = __Pyx_PyInt_from_py_npy_uint8(PyTuple_GET_ITEM(__pyx_args, 3)); if (unlikely((__pyx_v_val == (npy_uint8)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 43; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __pyx_v_val = __Pyx_PyInt_from_py_npy_uint8(PyTuple_GET_ITEM(__pyx_args, 3)); if (unlikely((__pyx_v_val == (npy_uint8)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 48; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
__pyx_v_out = ((PyArrayObject *)PyTuple_GET_ITEM(__pyx_args, 4));
}
goto __pyx_L4_argument_unpacking_done;
__pyx_L5_argtuple_error:;
- __Pyx_RaiseArgtupleInvalid("dot_over", 1, 5, 5, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 39; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __Pyx_RaiseArgtupleInvalid("dot_over", 1, 5, 5, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 44; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
__pyx_L3_error:;
__Pyx_AddTraceback("scikits.learn.utils.arrayfuncs.dot_over");
return NULL;
@@ -872,36 +876,36 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over(PyObject
__pyx_bstruct_y.buf = NULL;
__pyx_bstruct_mask.buf = NULL;
__pyx_bstruct_out.buf = NULL;
- if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_X), __pyx_ptype_5numpy_ndarray, 1, "X", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 40; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_y), __pyx_ptype_5numpy_ndarray, 1, "y", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 41; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_mask), __pyx_ptype_5numpy_ndarray, 1, "mask", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 42; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_out), __pyx_ptype_5numpy_ndarray, 1, "out", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 44; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_X), __pyx_ptype_5numpy_ndarray, 1, "X", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 45; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_y), __pyx_ptype_5numpy_ndarray, 1, "y", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 46; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_mask), __pyx_ptype_5numpy_ndarray, 1, "mask", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 47; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_out), __pyx_ptype_5numpy_ndarray, 1, "out", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 49; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
{
__Pyx_BufFmt_StackElem __pyx_stack[1];
- if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_X, (PyObject*)__pyx_v_X, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 39; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_X, (PyObject*)__pyx_v_X, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 44; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
__pyx_bstride_0_X = __pyx_bstruct_X.strides[0]; __pyx_bstride_1_X = __pyx_bstruct_X.strides[1];
__pyx_bshape_0_X = __pyx_bstruct_X.shape[0]; __pyx_bshape_1_X = __pyx_bstruct_X.shape[1];
{
__Pyx_BufFmt_StackElem __pyx_stack[1];
- if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_y, (PyObject*)__pyx_v_y, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 39; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_y, (PyObject*)__pyx_v_y, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 44; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
__pyx_bstride_0_y = __pyx_bstruct_y.strides[0];
__pyx_bshape_0_y = __pyx_bstruct_y.shape[0];
{
__Pyx_BufFmt_StackElem __pyx_stack[1];
- if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_mask, (PyObject*)__pyx_v_mask, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_BOOL, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 39; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_mask, (PyObject*)__pyx_v_mask, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_BOOL, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 44; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
__pyx_bstride_0_mask = __pyx_bstruct_mask.strides[0];
__pyx_bshape_0_mask = __pyx_bstruct_mask.shape[0];
{
__Pyx_BufFmt_StackElem __pyx_stack[1];
- if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_out, (PyObject*)__pyx_v_out, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 1, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 39; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_out, (PyObject*)__pyx_v_out, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES| PyBUF_WRITABLE, 1, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 44; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
__pyx_bstride_0_out = __pyx_bstruct_out.strides[0];
__pyx_bshape_0_out = __pyx_bstruct_out.shape[0];
- /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":49
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":54
* """
*
* cdef int i, j=0, N, incX, incY, incXY # <<<<<<<<<<<<<<
@@ -910,7 +914,7 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over(PyObject
*/
__pyx_v_j = 0;
- /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":50
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":55
*
* cdef int i, j=0, N, incX, incY, incXY
* cdef double *X_ptr = <double *> X.data # <<<<<<<<<<<<<<
@@ -919,7 +923,7 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over(PyObject
*/
__pyx_v_X_ptr = ((double *)__pyx_v_X->data);
- /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":51
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":56
* cdef int i, j=0, N, incX, incY, incXY
* cdef double *X_ptr = <double *> X.data
* N = X.shape[1] # <<<<<<<<<<<<<<
@@ -928,7 +932,7 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over(PyObject
*/
__pyx_v_N = (__pyx_v_X->dimensions[1]);
- /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":52
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":57
* cdef double *X_ptr = <double *> X.data
* N = X.shape[1]
* incX = X.strides[1] / sizeof (double) # <<<<<<<<<<<<<<
@@ -939,11 +943,11 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over(PyObject
__pyx_t_2 = (sizeof(double));
if (unlikely(__pyx_t_2 == 0)) {
PyErr_Format(PyExc_ZeroDivisionError, "integer division or modulo by zero");
- {__pyx_filename = __pyx_f[0]; __pyx_lineno = 52; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ {__pyx_filename = __pyx_f[0]; __pyx_lineno = 57; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
__pyx_v_incX = (__pyx_t_1 / __pyx_t_2);
- /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":53
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":58
* N = X.shape[1]
* incX = X.strides[1] / sizeof (double)
* incY = y.strides[0] / sizeof (double) # <<<<<<<<<<<<<<
@@ -954,11 +958,11 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over(PyObject
__pyx_t_2 = (sizeof(double));
if (unlikely(__pyx_t_2 == 0)) {
PyErr_Format(PyExc_ZeroDivisionError, "integer division or modulo by zero");
- {__pyx_filename = __pyx_f[0]; __pyx_lineno = 53; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ {__pyx_filename = __pyx_f[0]; __pyx_lineno = 58; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
__pyx_v_incY = (__pyx_t_1 / __pyx_t_2);
- /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":54
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":59
* incX = X.strides[1] / sizeof (double)
* incY = y.strides[0] / sizeof (double)
* incXY = X.strides[0] / sizeof (double) # <<<<<<<<<<<<<<
@@ -969,11 +973,11 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over(PyObject
__pyx_t_2 = (sizeof(double));
if (unlikely(__pyx_t_2 == 0)) {
PyErr_Format(PyExc_ZeroDivisionError, "integer division or modulo by zero");
- {__pyx_filename = __pyx_f[0]; __pyx_lineno = 54; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ {__pyx_filename = __pyx_f[0]; __pyx_lineno = 59; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
__pyx_v_incXY = (__pyx_t_1 / __pyx_t_2);
- /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":56
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":61
* incXY = X.strides[0] / sizeof (double)
*
* for i in range(X.shape[0]): # <<<<<<<<<<<<<<
@@ -984,7 +988,7 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over(PyObject
for (__pyx_t_3 = 0; __pyx_t_3 < __pyx_t_1; __pyx_t_3+=1) {
__pyx_v_i = __pyx_t_3;
- /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":57
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":62
*
* for i in range(X.shape[0]):
* if mask[i] == val: # <<<<<<<<<<<<<<
@@ -999,12 +1003,12 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over(PyObject
} else if (unlikely(__pyx_t_4 >= __pyx_bshape_0_mask)) __pyx_t_5 = 0;
if (unlikely(__pyx_t_5 != -1)) {
__Pyx_RaiseBufferIndexError(__pyx_t_5);
- {__pyx_filename = __pyx_f[0]; __pyx_lineno = 57; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ {__pyx_filename = __pyx_f[0]; __pyx_lineno = 62; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
__pyx_t_6 = ((*__Pyx_BufPtrStrided1d(__pyx_t_7scikits_5learn_5utils_10arrayfuncs_BOOL *, __pyx_bstruct_mask.buf, __pyx_t_4, __pyx_bstride_0_mask)) == __pyx_v_val);
if (__pyx_t_6) {
- /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":58
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":63
* for i in range(X.shape[0]):
* if mask[i] == val:
* out[j] = cblas_ddot (N, X_ptr + i*incXY, incX, <double *> y.data, incY) # <<<<<<<<<<<<<<
@@ -1019,11 +1023,11 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over(PyObject
} else if (unlikely(__pyx_t_5 >= __pyx_bshape_0_out)) __pyx_t_7 = 0;
if (unlikely(__pyx_t_7 != -1)) {
__Pyx_RaiseBufferIndexError(__pyx_t_7);
- {__pyx_filename = __pyx_f[0]; __pyx_lineno = 58; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ {__pyx_filename = __pyx_f[0]; __pyx_lineno = 63; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
*__Pyx_BufPtrStrided1d(__pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE *, __pyx_bstruct_out.buf, __pyx_t_5, __pyx_bstride_0_out) = cblas_ddot(__pyx_v_N, (__pyx_v_X_ptr + (__pyx_v_i * __pyx_v_incXY)), __pyx_v_incX, ((double *)__pyx_v_y->data), __pyx_v_incY);
- /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":59
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":64
* if mask[i] == val:
* out[j] = cblas_ddot (N, X_ptr + i*incXY, incX, <double *> y.data, incY)
* j += 1 # <<<<<<<<<<<<<<
@@ -1064,7 +1068,7 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over(PyObject
return __pyx_r;
}
-/* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":62
+/* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":67
*
*
* def solve_triangular ( # <<<<<<<<<<<<<<
@@ -1110,11 +1114,11 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_solve_triangular(P
values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__y);
if (likely(values[1])) kw_args--;
else {
- __Pyx_RaiseArgtupleInvalid("solve_triangular", 1, 2, 2, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 62; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __Pyx_RaiseArgtupleInvalid("solve_triangular", 1, 2, 2, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 67; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
}
}
if (unlikely(kw_args > 0)) {
- if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, PyTuple_GET_SIZE(__pyx_args), "solve_triangular") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 62; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, PyTuple_GET_SIZE(__pyx_args), "solve_triangular") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 67; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
}
__pyx_v_X = ((PyArrayObject *)values[0]);
__pyx_v_y = ((PyArrayObject *)values[1]);
@@ -1126,29 +1130,29 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_solve_triangular(P
}
goto __pyx_L4_argument_unpacking_done;
__pyx_L5_argtuple_error:;
- __Pyx_RaiseArgtupleInvalid("solve_triangular", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 62; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __Pyx_RaiseArgtupleInvalid("solve_triangular", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 67; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
__pyx_L3_error:;
__Pyx_AddTraceback("scikits.learn.utils.arrayfuncs.solve_triangular");
return NULL;
__pyx_L4_argument_unpacking_done:;
__pyx_bstruct_X.buf = NULL;
__pyx_bstruct_y.buf = NULL;
- if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_X), __pyx_ptype_5numpy_ndarray, 1, "X", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 63; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
- if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_y), __pyx_ptype_5numpy_ndarray, 1, "y", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 64; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_X), __pyx_ptype_5numpy_ndarray, 1, "X", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 68; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_y), __pyx_ptype_5numpy_ndarray, 1, "y", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 69; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
{
__Pyx_BufFmt_StackElem __pyx_stack[1];
- if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_X, (PyObject*)__pyx_v_X, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 62; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_X, (PyObject*)__pyx_v_X, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 67; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
__pyx_bstride_0_X = __pyx_bstruct_X.strides[0]; __pyx_bstride_1_X = __pyx_bstruct_X.strides[1];
__pyx_bshape_0_X = __pyx_bstruct_X.shape[0]; __pyx_bshape_1_X = __pyx_bstruct_X.shape[1];
{
__Pyx_BufFmt_StackElem __pyx_stack[1];
- if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_y, (PyObject*)__pyx_v_y, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 62; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_y, (PyObject*)__pyx_v_y, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES, 1, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 67; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
__pyx_bstride_0_y = __pyx_bstruct_y.strides[0];
__pyx_bshape_0_y = __pyx_bstruct_y.shape[0];
- /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":71
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":76
* TODO: option for upper, lower, etc.
* """
* cdef int lda = <int> X.strides[0] / sizeof(double) # <<<<<<<<<<<<<<
@@ -1159,14 +1163,16 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_solve_triangular(P
__pyx_t_2 = (sizeof(double));
if (unlikely(__pyx_t_2 == 0)) {
PyErr_Format(PyExc_ZeroDivisionError, "integer division or modulo by zero");
- {__pyx_filename = __pyx_f[0]; __pyx_lineno = 71; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ {__pyx_filename = __pyx_f[0]; __pyx_lineno = 76; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
}
__pyx_v_lda = (__pyx_t_1 / __pyx_t_2);
- /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":75
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":80
* cblas_dtrsv (CblasRowMajor, CblasLower, CblasNoTrans,
* CblasNonUnit, <int> X.shape[0], <double *> X.data,
* lda, <double *> y.data, 1); # <<<<<<<<<<<<<<
+ *
+ *
*/
cblas_dtrsv(CblasRowMajor, CblasLower, CblasNoTrans, CblasNonUnit, ((int)(__pyx_v_X->dimensions[0])), ((double *)__pyx_v_X->data), __pyx_v_lda, ((double *)__pyx_v_y->data), 1);
@@ -1190,6 +1196,127 @@ static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_solve_triangular(P
return __pyx_r;
}
+/* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":83
+ *
+ *
+ * def cholesky_delete ( # <<<<<<<<<<<<<<
+ * np.ndarray[DOUBLE, ndim=2] L,
+ * int go_out):
+ */
+
+static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_cholesky_delete(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds); /*proto*/
+static PyObject *__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_cholesky_delete(PyObject *__pyx_self, PyObject *__pyx_args, PyObject *__pyx_kwds) {
+ PyArrayObject *__pyx_v_L = 0;
+ int __pyx_v_go_out;
+ int __pyx_v_n;
+ int __pyx_v_m;
+ Py_buffer __pyx_bstruct_L;
+ Py_ssize_t __pyx_bstride_0_L = 0;
+ Py_ssize_t __pyx_bstride_1_L = 0;
+ Py_ssize_t __pyx_bshape_0_L = 0;
+ Py_ssize_t __pyx_bshape_1_L = 0;
+ PyObject *__pyx_r = NULL;
+ int __pyx_t_1;
+ size_t __pyx_t_2;
+ static PyObject **__pyx_pyargnames[] = {&__pyx_n_s__L,&__pyx_n_s__go_out,0};
+ __Pyx_RefNannySetupContext("cholesky_delete");
+ __pyx_self = __pyx_self;
+ if (unlikely(__pyx_kwds)) {
+ Py_ssize_t kw_args = PyDict_Size(__pyx_kwds);
+ PyObject* values[2] = {0,0};
+ switch (PyTuple_GET_SIZE(__pyx_args)) {
+ case 2: values[1] = PyTuple_GET_ITEM(__pyx_args, 1);
+ case 1: values[0] = PyTuple_GET_ITEM(__pyx_args, 0);
+ case 0: break;
+ default: goto __pyx_L5_argtuple_error;
+ }
+ switch (PyTuple_GET_SIZE(__pyx_args)) {
+ case 0:
+ values[0] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__L);
+ if (likely(values[0])) kw_args--;
+ else goto __pyx_L5_argtuple_error;
+ case 1:
+ values[1] = PyDict_GetItem(__pyx_kwds, __pyx_n_s__go_out);
+ if (likely(values[1])) kw_args--;
+ else {
+ __Pyx_RaiseArgtupleInvalid("cholesky_delete", 1, 2, 2, 1); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ }
+ }
+ if (unlikely(kw_args > 0)) {
+ if (unlikely(__Pyx_ParseOptionalKeywords(__pyx_kwds, __pyx_pyargnames, 0, values, PyTuple_GET_SIZE(__pyx_args), "cholesky_delete") < 0)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ }
+ __pyx_v_L = ((PyArrayObject *)values[0]);
+ __pyx_v_go_out = __Pyx_PyInt_AsInt(values[1]); if (unlikely((__pyx_v_go_out == (int)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 85; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ } else if (PyTuple_GET_SIZE(__pyx_args) != 2) {
+ goto __pyx_L5_argtuple_error;
+ } else {
+ __pyx_v_L = ((PyArrayObject *)PyTuple_GET_ITEM(__pyx_args, 0));
+ __pyx_v_go_out = __Pyx_PyInt_AsInt(PyTuple_GET_ITEM(__pyx_args, 1)); if (unlikely((__pyx_v_go_out == (int)-1) && PyErr_Occurred())) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 85; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ }
+ goto __pyx_L4_argument_unpacking_done;
+ __pyx_L5_argtuple_error:;
+ __Pyx_RaiseArgtupleInvalid("cholesky_delete", 1, 2, 2, PyTuple_GET_SIZE(__pyx_args)); {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L3_error;}
+ __pyx_L3_error:;
+ __Pyx_AddTraceback("scikits.learn.utils.arrayfuncs.cholesky_delete");
+ return NULL;
+ __pyx_L4_argument_unpacking_done:;
+ __pyx_bstruct_L.buf = NULL;
+ if (unlikely(!__Pyx_ArgTypeTest(((PyObject *)__pyx_v_L), __pyx_ptype_5numpy_ndarray, 1, "L", 0))) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 84; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ {
+ __Pyx_BufFmt_StackElem __pyx_stack[1];
+ if (unlikely(__Pyx_GetBufferAndValidate(&__pyx_bstruct_L, (PyObject*)__pyx_v_L, &__Pyx_TypeInfo_nn___pyx_t_7scikits_5learn_5utils_10arrayfuncs_DOUBLE, PyBUF_FORMAT| PyBUF_STRIDES, 2, 0, __pyx_stack) == -1)) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 83; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ }
+ __pyx_bstride_0_L = __pyx_bstruct_L.strides[0]; __pyx_bstride_1_L = __pyx_bstruct_L.strides[1];
+ __pyx_bshape_0_L = __pyx_bstruct_L.shape[0]; __pyx_bshape_1_L = __pyx_bstruct_L.shape[1];
+
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":87
+ * int go_out):
+ *
+ * cdef int n = <int> L.shape[0] # <<<<<<<<<<<<<<
+ * cdef int m = <int> L.strides[0] / sizeof (double)
+ * c_cholesky_delete (m, n, <double *> L.data, go_out)
+ */
+ __pyx_v_n = ((int)(__pyx_v_L->dimensions[0]));
+
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":88
+ *
+ * cdef int n = <int> L.shape[0]
+ * cdef int m = <int> L.strides[0] / sizeof (double) # <<<<<<<<<<<<<<
+ * c_cholesky_delete (m, n, <double *> L.data, go_out)
+ */
+ __pyx_t_1 = ((int)(__pyx_v_L->strides[0]));
+ __pyx_t_2 = (sizeof(double));
+ if (unlikely(__pyx_t_2 == 0)) {
+ PyErr_Format(PyExc_ZeroDivisionError, "integer division or modulo by zero");
+ {__pyx_filename = __pyx_f[0]; __pyx_lineno = 88; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ }
+ __pyx_v_m = (__pyx_t_1 / __pyx_t_2);
+
+ /* "/Users/fabian/dev/scikit-learn/scikits/learn/utils/arrayfuncs.pyx":89
+ * cdef int n = <int> L.shape[0]
+ * cdef int m = <int> L.strides[0] / sizeof (double)
+ * c_cholesky_delete (m, n, <double *> L.data, go_out) # <<<<<<<<<<<<<<
+ */
+ cholesky_delete(__pyx_v_m, __pyx_v_n, ((double *)__pyx_v_L->data), __pyx_v_go_out);
+
+ __pyx_r = Py_None; __Pyx_INCREF(Py_None);
+ goto __pyx_L0;
+ __pyx_L1_error:;
+ { PyObject *__pyx_type, *__pyx_value, *__pyx_tb;
+ __Pyx_ErrFetch(&__pyx_type, &__pyx_value, &__pyx_tb);
+ __Pyx_SafeReleaseBuffer(&__pyx_bstruct_L);
+ __Pyx_ErrRestore(__pyx_type, __pyx_value, __pyx_tb);}
+ __Pyx_AddTraceback("scikits.learn.utils.arrayfuncs.cholesky_delete");
+ __pyx_r = NULL;
+ goto __pyx_L2;
+ __pyx_L0:;
+ __Pyx_SafeReleaseBuffer(&__pyx_bstruct_L);
+ __pyx_L2:;
+ __Pyx_XGIVEREF(__pyx_r);
+ __Pyx_RefNannyFinishContext();
+ return __pyx_r;
+}
+
/* "/opt/local/Library/Frameworks/Python.framework/Versions/2.6/lib/python2.6/site-packages/Cython/Includes/numpy.pxd":187
* # experimental exception made for __getbuffer__ and __releasebuffer__
* # -- the details of this may change.
@@ -3125,6 +3252,7 @@ static CYTHON_INLINE PyObject *__pyx_f_5numpy_get_array_base(PyArrayObject *__py
static struct PyMethodDef __pyx_methods[] = {
{__Pyx_NAMESTR("dot_over"), (PyCFunction)__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_dot_over, METH_VARARGS|METH_KEYWORDS, __Pyx_DOCSTR(__pyx_doc_7scikits_5learn_5utils_10arrayfuncs_dot_over)},
{__Pyx_NAMESTR("solve_triangular"), (PyCFunction)__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_solve_triangular, METH_VARARGS|METH_KEYWORDS, __Pyx_DOCSTR(__pyx_doc_7scikits_5learn_5utils_10arrayfuncs_solve_triangular)},
+ {__Pyx_NAMESTR("cholesky_delete"), (PyCFunction)__pyx_pf_7scikits_5learn_5utils_10arrayfuncs_cholesky_delete, METH_VARARGS|METH_KEYWORDS, __Pyx_DOCSTR(0)},
{0, 0, 0, 0}
};
@@ -3153,6 +3281,7 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = {
{&__pyx_kp_u_6, __pyx_k_6, sizeof(__pyx_k_6), 0, 1, 0, 0},
{&__pyx_kp_u_8, __pyx_k_8, sizeof(__pyx_k_8), 0, 1, 0, 0},
{&__pyx_kp_u_9, __pyx_k_9, sizeof(__pyx_k_9), 0, 1, 0, 0},
+ {&__pyx_n_s__L, __pyx_k__L, sizeof(__pyx_k__L), 0, 0, 1, 1},
{&__pyx_n_s__RuntimeError, __pyx_k__RuntimeError, sizeof(__pyx_k__RuntimeError), 0, 0, 1, 1},
{&__pyx_n_s__ValueError, __pyx_k__ValueError, sizeof(__pyx_k__ValueError), 0, 0, 1, 1},
{&__pyx_n_s__X, __pyx_k__X, sizeof(__pyx_k__X), 0, 0, 1, 1},
@@ -3166,6 +3295,7 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = {
{&__pyx_n_s__dot_over, __pyx_k__dot_over, sizeof(__pyx_k__dot_over), 0, 0, 1, 1},
{&__pyx_n_s__fields, __pyx_k__fields, sizeof(__pyx_k__fields), 0, 0, 1, 1},
{&__pyx_n_s__format, __pyx_k__format, sizeof(__pyx_k__format), 0, 0, 1, 1},
+ {&__pyx_n_s__go_out, __pyx_k__go_out, sizeof(__pyx_k__go_out), 0, 0, 1, 1},
{&__pyx_n_s__itemsize, __pyx_k__itemsize, sizeof(__pyx_k__itemsize), 0, 0, 1, 1},
{&__pyx_n_s__mask, __pyx_k__mask, sizeof(__pyx_k__mask), 0, 0, 1, 1},
{&__pyx_n_s__names, __pyx_k__names, sizeof(__pyx_k__names), 0, 0, 1, 1},
@@ -3186,7 +3316,7 @@ static __Pyx_StringTabEntry __pyx_string_tab[] = {
{0, 0, 0, 0, 0, 0, 0}
};
static int __Pyx_InitCachedBuiltins(void) {
- __pyx_builtin_range = __Pyx_GetName(__pyx_b, __pyx_n_s__range); if (!__pyx_builtin_range) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 56; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
+ __pyx_builtin_range = __Pyx_GetName(__pyx_b, __pyx_n_s__range); if (!__pyx_builtin_range) {__pyx_filename = __pyx_f[0]; __pyx_lineno = 61; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__pyx_builtin_ValueError = __Pyx_GetName(__pyx_b, __pyx_n_s__ValueError); if (!__pyx_builtin_ValueError) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 205; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
__pyx_builtin_RuntimeError = __Pyx_GetName(__pyx_b, __pyx_n_s__RuntimeError); if (!__pyx_builtin_RuntimeError) {__pyx_filename = __pyx_f[1]; __pyx_lineno = 786; __pyx_clineno = __LINE__; goto __pyx_L1_error;}
return 0;
View
14 scikits/learn/utils/arrayfuncs.pyx
@@ -30,6 +30,11 @@ cdef extern from "cblas.h":
int incX)
+cdef extern from "src/cholesky_delete.c":
+ int c_cholesky_delete "cholesky_delete" (int m, int n, double *L, int go_out)
+
+
+
ctypedef np.float64_t DOUBLE
ctypedef np.uint8_t BOOL
# cython has no support for np.bool type. See
@@ -73,3 +78,12 @@ def solve_triangular (
cblas_dtrsv (CblasRowMajor, CblasLower, CblasNoTrans,
CblasNonUnit, <int> X.shape[0], <double *> X.data,
lda, <double *> y.data, 1);
+
+
+def cholesky_delete (
+ np.ndarray[DOUBLE, ndim=2] L,
+ int go_out):
+
+ cdef int n = <int> L.shape[0]
+ cdef int m = <int> L.strides[0] / sizeof (double)
+ c_cholesky_delete (m, n, <double *> L.data, go_out)
View
47 scikits/learn/utils/src/cholesky_delete.c
@@ -0,0 +1,47 @@
+#include <math.h>
+#include <cblas.h>
+
+#ifdef _MSC_VER
+ #define copysign _copysign
+#endif
+
+/* General Cholesky Delete.
+ * Remove an element from the cholesky factorization
+ * m = columns
+ * n = rows
+ *
+ * TODO: put transpose as an option
+ *
+ */
+int cholesky_delete (int m, int n, double *L, int go_out) {
+
+ double c, s;
+
+ /* delete row go_out */
+ double * _L = L + (go_out * m);
+ int i;
+ for (i = go_out; i < n - 1; ++i) {
+ cblas_dcopy (i + 2, _L + m , 1, _L, 1);
+ _L += m;
+ }
+
+ _L = L + (go_out * m);
+ for (i=go_out; i < n - 1; ++i) {
+
+ cblas_drotg (_L + i, _L + i + 1, &c, &s);
+ if (_L[i] < 0) {
+ /* Diagonals cannot be negative */
+ _L[i] = copysign(_L[i], 1.0);
+ c = -c;
+ s = -s;
+ }
+ _L[i+1] = 0.; /* just for cleanup */
+ _L += m;
+
+ cblas_drot (n - (i + 2), _L + i, m, _L + i + 1,
+ m, c, s);
+
+ }
+
+ return 0;
+}
Please sign in to comment.
Something went wrong with that request. Please try again.