ENH: refactor: replace naked mallocs with cython arrays

While slightly slower, allocation only takes place once, outside
of critical loops. The advantage of not having to worry about
leaks (since cython arrays are refcounted) is worth it.

Use cvarray instead of np.empty

Author: vene

polylearn/cd_direct_fast.cpp

@@ -7,8 +7,8 @@
 # Author: Vlad Niculae
 # License: BSD
 
-from libc.stdlib cimport malloc, free
 from libc.math cimport fabs
+from cython.view cimport array
 
 from lightning.impl.dataset_fast cimport ColumnDataset
 
@@ -19,7 +19,7 @@ from .cd_linear_fast cimport _cd_linear_epoch
 cdef void _precompute(ColumnDataset X,
                       double[:, :, ::1] P,
                       Py_ssize_t order,
-                      double* out,
+                      double[:, ::1] out,
                       Py_ssize_t s,
                       unsigned int degree):
 
@@ -32,15 +32,17 @@ cdef void _precompute(ColumnDataset X,
     cdef int n_nz
 
     cdef Py_ssize_t i, j, ii
+    cdef unsigned int d
+    cdef double tmp
 
     for i in range(n_samples):
-        out[i] = 0
-    
+        out[degree - 1, i] = 0
+
     for j in range(n_features):
         X.get_column_ptr(j, &indices, &data, &n_nz)
         for ii in range(n_nz):
             i = indices[ii]
-            out[i] += (data[ii] * P[order, s, j]) ** degree
+            out[degree - 1, i] += (data[ii] * P[order, s, j]) ** degree
 
 
 cdef inline double _update(int* indices,
@@ -50,12 +52,11 @@ cdef inline double _update(int* indices,
                            double[:] y,
                            double[:] y_pred,
                            LossFunction loss,
-                           double* d1,
-                           double* d2,
                            unsigned int degree,
                            double lam,
                            double beta,
-                           double* cache_kp):
+                           double[:, ::1] D,
+                           double[:] cache_kp):
 
     cdef double l1_reg = 2 * beta * fabs(lam)
 
@@ -70,10 +71,10 @@ cdef inline double _update(int* indices,
         i = indices[ii]
 
         if degree == 2:
-            kp = d1[i] - p_js * data[ii]
-        elif degree == 3:
-            kp = 0.5 * (d1[i] ** 2 - d2[i])
-            kp -= p_js * data[ii] * d1[i]
+            kp = D[0, i] - p_js * data[ii]
+        else:  #  degree == 3:
+            kp = 0.5 * (D[0, i] ** 2 - D[1, i])
+            kp -= p_js * data[ii] * D[0, i]
             kp += p_js ** 2 * data[ii] ** 2
 
         kp *= lam * data[ii]
@@ -97,12 +98,11 @@ cdef inline double _cd_direct_epoch(double[:, :, ::1] P,
                                     double[:] y,
                                     double[:] y_pred,
                                     double[:] lams,
-                                    double* d1,
-                                    double* d2,
                                     unsigned int degree,
                                     double beta,
                                     LossFunction loss,
-                                    double* cache_kp):
+                                    double[:, ::1] D,
+                                    double[:] cache_kp):
 
     cdef Py_ssize_t s, j
     cdef double p_old, update, offset
@@ -118,9 +118,9 @@ cdef inline double _cd_direct_epoch(double[:, :, ::1] P,
     for s in range(n_components):
 
         # initialize the cached ds for this s
-        _precompute(X, P, order, d1, s, 1)
+        _precompute(X, P, order, D, s, 1)
         if degree == 3:
-            _precompute(X, P, order, d2, s, 2)
+            _precompute(X, P, order, D, s, 2)
 
         for j in range(n_features):
 
@@ -129,7 +129,7 @@ cdef inline double _cd_direct_epoch(double[:, :, ::1] P,
             # compute coordinate update
             p_old = P[order, s, j]
             update = _update(indices, data, n_nz, p_old, y, y_pred,
-                             loss, d1, d2, degree, lams[s], beta, cache_kp)
+                             loss, degree, lams[s], beta, D, cache_kp)
             P[order, s, j] -= update
             sum_viol += fabs(update)
 
@@ -138,9 +138,10 @@ cdef inline double _cd_direct_epoch(double[:, :, ::1] P,
                 i = indices[ii]
 
                 if degree == 3:
-                    d2[i] -= (p_old ** 2 - P[order, s, j] ** 2) * data[ii] ** 2
+                    D[1, i] -= ((p_old ** 2 - P[order, s, j] ** 2) *
+                                data[ii] ** 2)
 
-                d1[i] -= update * data[ii]
+                D[0, i] -= update * data[ii]
                 y_pred[i] -= update * cache_kp[ii]
     return sum_viol
 
@@ -169,11 +170,8 @@ def _cd_direct_ho(double[:, :, ::1] P not None,
     cdef bint converged = False
 
     # precomputed values
-    cdef double *d1 = <double *> malloc(n_samples * sizeof(double))
-    cdef double *d2
-    if degree == 3:
-        d2 = <double *> malloc(n_samples * sizeof(double))
-    cdef double *cache_kp = <double *> malloc(n_samples * sizeof(double))
+    cdef double[:, ::1] D = array((degree - 1, n_samples), sizeof(double), 'd')
+    cdef double[:] cache_kp = array((n_samples,), sizeof(double), 'd')
 
     for it in range(max_iter):
         viol = 0
@@ -182,11 +180,11 @@ def _cd_direct_ho(double[:, :, ::1] P not None,
             viol += _cd_linear_epoch(w, X, y, y_pred, col_norm_sq, alpha, loss)
 
         if fit_lower and degree == 3:  # fit degree 2. Will be looped later.
-            viol += _cd_direct_epoch(P, 1, X, y, y_pred, lams, d1, d2,
-                                     2, beta, loss, cache_kp)
+            viol += _cd_direct_epoch(P, 1, X, y, y_pred, lams, 2, beta, loss,
+                                     D, cache_kp)
 
-        viol += _cd_direct_epoch(P, 0, X, y, y_pred, lams, d1, d2,
-                                 degree, beta, loss, cache_kp)
+        viol += _cd_direct_epoch(P, 0, X, y, y_pred, lams, degree, beta, loss,
+                                 D, cache_kp)
 
         if verbose:
             print("Iteration", it + 1, "violation sum", viol)
@@ -197,10 +195,4 @@ def _cd_direct_ho(double[:, :, ::1] P not None,
             converged = True
             break
 
-    # Free up cache
-    free(d1)
-    free(cache_kp)
-    if degree == 3:
-        free(d2)
-
     return converged