Merge pull request #35 from jdtuck/rbgs_c++

A c++ version of RLBFGS

.github/workflows/python-package.yml

@@ -21,6 +21,7 @@ jobs:
         python-version: ${{ matrix.python-version }}
     - name: Install dependencies
       run: |
+        sudo apt-get install -y libopenblas-dev
         python -m pip install --upgrade pip
         python -m pip install flake8 pytest pytest-runner coverage
         if [ -f requirements.txt ]; then pip install -r requirements.txt; fi

.gitignore

@@ -25,3 +25,6 @@ src/imagecpp.cpp
 fdasrsf_python.code-workspace
 *.so
 .coverage
+src/*.dSYM/
+src/rbfgs
+src/crbfgs.cpp

CHANGES.txt

@@ -1,4 +1,4 @@
-v2.5.X,  2023-12-05 -- add projection and variance explained to fpca
+v2.5.X,  2023-12-05 -- add projection and variance explained to fpca, add c version of rbfgs
 v2.5.2,  2023-11-18 -- bugfixes
 v2.5.0,  2023-11-17 -- bugfixes, add elastic changepoint
 v2.4.3,  2023-08-21 -- bugfixes, expose lam in curve, add example

fdasrsf/utility_functions.py

@@ -21,6 +21,7 @@
 import optimum_reparamN2 as orN2
 import optimum_reparam_N as orN
 import cbayesian as bay
+import crbfgs as cr
 import fdasrsf.geometry as geo
 from fdasrsf.rbfgs import rlbfgs
 import sys
@@ -43,8 +44,8 @@ def smooth_data(f, sparam=1):
     fo = f.copy()
     for k in range(0, sparam):
         for r in range(0, N):
-            fo[1 : (M - 2), r] = (
-                fo[0 : (M - 3), r] + 2 * fo[1 : (M - 2), r] + fo[2 : (M - 1), r]
+            fo[1: (M - 2), r] = (
+                fo[0: (M - 3), r] + 2 * fo[1: (M - 2), r] + fo[2: (M - 1), r]
             ) / 4
     return fo
 
@@ -135,11 +136,13 @@ def optimum_reparam(
     :param q1: vector of size N or array of NxM samples of first SRSF
     :param time: vector of size N describing the sample points
     :param q2: vector of size N or array of NxM samples samples of second SRSF
-    :param method: method to apply optimization (default="DP2") options are "DP","DP2","RBFGS"
+    :param method: method to apply optimization (default="DP2") options are 
+                   "DP","DP2","RBFGS","cRBFGS"
     :param lam: controls the amount of elasticity (default = 0.0)
-    :param penalty: penalty type (default="roughness") options are "roughness", "l2gam",
-                    "l2psi", "geodesic". Only roughness implemented in all methods. To use
-                    others method needs to be "RBFGS"
+    :param penalty: penalty type (default="roughness") options are "roughness", 
+                    "l2gam", "l2psi", "geodesic". Only roughness implemented 
+                    in all methods. To use others method needs to be "RBFGS"
+                    or "cRBFGS"
     :param grid_dim: size of the grid, for the DP2 method only (default = 7)
 
     :rtype: vector
@@ -172,17 +175,20 @@ def optimum_reparam(
     elif method == "DP2":
         if q1.ndim == 1 and q2.ndim == 1:
             gam = orN2.coptimum_reparam(
-                ascontiguousarray(q1), time, ascontiguousarray(q2), lam, grid_dim
+                ascontiguousarray(q1), time, ascontiguousarray(q2), lam,
+                grid_dim
             )
 
         if q1.ndim == 1 and q2.ndim == 2:
             gam = orN2.coptimum_reparamN(
-                ascontiguousarray(q1), time, ascontiguousarray(q2), lam, grid_dim
+                ascontiguousarray(q1), time, ascontiguousarray(q2), lam,
+                grid_dim
             )
 
         if q1.ndim == 2 and q2.ndim == 2:
             gam = orN2.coptimum_reparamN2(
-                ascontiguousarray(q1), time, ascontiguousarray(q2), lam, grid_dim
+                ascontiguousarray(q1), time, ascontiguousarray(q2), lam,
+                grid_dim
             )
     elif method == "RBFGS":
         if q1.ndim == 1 and q2.ndim == 1:
@@ -206,6 +212,38 @@ def optimum_reparam(
                 obj = rlbfgs(q1[:, i], q2[:, i], time)
                 obj.solve(lam=lam, penalty=penalty)
                 gam[:, i] = obj.gammaOpt
+    elif method == "cRBFGS":
+        if penalty == "roughness":
+            pen = 0
+        elif penalty == "l2gam":
+            pen = 1
+        elif penalty == "l2psi":
+            pen = 2
+        elif penalty == "geodesic":
+            pen = 3
+        else:
+            raise Exception("penalty not implemented")
+
+        if q1.ndim == 1 and q2.ndim == 1:
+            time = linspace(0, 1, q1.shape[0])
+            gam = cr.rlbfgs(ascontiguousarray(q1), ascontiguousarray(q2),
+                            ascontiguousarray(time), 30, lam, pen)
+
+        if q1.ndim == 1 and q2.ndim == 2:
+            gam = zeros(q2.shape)
+            time = linspace(0, 1, q1.shape[0])
+            for i in range(0, q2.shape[1]):
+                gam[:, i] = cr.rlbfgs(ascontiguousarray(q1),
+                                      ascontiguousarray(q2[:, i]),
+                                      ascontiguousarray(time), 30, lam, pen)
+
+        if q1.ndim == 2 and q2.ndim == 2:
+            gam = zeros(q2.shape)
+            time = linspace(0, 1, q1.shape[0])
+            for i in range(0, q2.shape[1]):
+                gam[:, i] = cr.rlbfgs(ascontiguousarray(q1[:, i]),
+                                      ascontiguousarray(q2[:, i]),
+                                      ascontiguousarray(time), 30, lam, pen)
 
     else:
         raise Exception("Invalid Optimization Method")
@@ -265,7 +303,8 @@ def elastic_depth(f, time, method="DP2", lam=0.0, parallel=True):
 
     :param f: matrix of size MxN (M time points for N functions)
     :param time: vector of size M describing the sample points
-    :param method: method to apply optimization (default="DP2") options are "DP","DP2","RBFGS"
+    :param method: method to apply optimization (default="DP2") 
+                   options are "DP","DP2","RBFGS","cRBFGS"
     :param lam: controls the elasticity (default = 0.0)
 
     :rtype: scalar
@@ -288,7 +327,8 @@ def elastic_depth(f, time, method="DP2", lam=0.0, parallel=True):
             phs_dist[i, :] = out[i][1]
     else:
         for i in range(0, fns):
-            amp_dist[i, :], phs_dist[i, :] = distmat(f, f[:, i], time, i, method)
+            amp_dist[i, :], phs_dist[i, :] = distmat(f, f[:, i], time, i,
+                                                     method)
 
     amp_dist = amp_dist + amp_dist.T
     phs_dist = phs_dist + phs_dist.T
@@ -309,7 +349,8 @@ def elastic_distance(f1, f2, time, method="DP2", lam=0.0, alpha=None):
     :param f1: vector of size N
     :param f2: vector of size N
     :param time: vector of size N describing the sample points
-    :param method: method to apply optimization (default="DP2") options are "DP","DP2","RBFGS"
+    :param method: method to apply optimization (default="DP2") 
+                   options are "DP","DP2","RBFGS","cRBFGS"
     :param lam: controls the elasticity (default = 0.0)
     :param alpha: makes alpha * dx + (1-alpha) * dy
 
@@ -340,7 +381,7 @@ def elastic_distance(f1, f2, time, method="DP2", lam=0.0, alpha=None):
     Dx = real(arccos(q1dotq2))
 
     if alpha is not None:
-        Dt = alpha * Dx + (1-alpha) * Dy
+        Dt = alpha * Dx + (1 - alpha) * Dy
         return Dy, Dx, Dt
     else:
         return Dy, Dx
@@ -423,7 +464,8 @@ def SqrtMeanInverse(gam):
 
 def SqrtMean(gam, parallel=False, cores=-1):
     """
-    calculates the srsf of warping functions with corresponding shooting vectors
+    calculates the srsf of warping functions with corresponding shooting 
+    vectors
 
     :param gam: numpy ndarray of shape (M,N) of M warping functions
                 with N samples
@@ -432,8 +474,10 @@ def SqrtMean(gam, parallel=False, cores=-1):
 
     :rtype: 2 numpy ndarray and vector
     :return mu: Karcher mean psi function
-    :return gam_mu: vector of dim N which is the Karcher mean warping function
-    :return psi: numpy ndarray of shape (M,N) of M SRSF of the warping functions
+    :return gam_mu: vector of dim N which is the Karcher mean warping 
+                    function
+    :return psi: numpy ndarray of shape (M,N) of M SRSF of the warping 
+                 functions
     :return vec: numpy ndarray of shape (M,N) of M shooting vectors
 
     """
@@ -462,7 +506,6 @@ def SqrtMean(gam, parallel=False, cores=-1):
     min_ind = dqq.argmin()
     mu = psi[:, min_ind]
     maxiter = 501
-    tt = 1
     lvm = zeros(maxiter)
     vec = zeros((T, n))
     stp = 0.3
@@ -512,15 +555,18 @@ def inv_exp_map_sub(mu, psi):
 
 def SqrtMedian(gam):
     """
-    calculates the median srsf of warping functions with corresponding shooting vectors
+    calculates the median srsf of warping functions with corresponding 
+    shooting vectors
 
     :param gam: numpy ndarray of shape (M,N) of M warping functions
                 with N samples
 
     :rtype: 2 numpy ndarray and vector
     :return gam_median: Karcher median warping function
-    :return psi_meidan: vector of dim N which is the Karcher median srsf function
-    :return psi: numpy ndarray of shape (M,N) of M SRSF of the warping functions
+    :return psi_meidan: vector of dim N which is the Karcher median srsf 
+                        function
+    :return psi: numpy ndarray of shape (M,N) of M SRSF of the warping 
+                 functions
     :return vec: numpy ndarray of shape (M,N) of M shooting vectors
 
     """
@@ -594,7 +640,7 @@ def cumtrapzmid(x, y, c, mid):
     fa[0:mid] = tmp[::-1]
 
     # case >= mid
-    fa[mid:a] = c + cumtrapz(y[mid - 1 : a - 1], x[mid - 1 : a - 1], initial=0)
+    fa[mid:a] = c + cumtrapz(y[mid - 1: a - 1], x[mid - 1: a - 1], initial=0)
 
     return fa
 
@@ -782,7 +828,6 @@ def geigen(Amat, Bmat, Cmat):
     p = Bmat.shape[0]
     q = Cmat.shape[0]
 
-    s = min(p, q)
     tmp = fabs(Bmat - Bmat.transpose())
     tmp1 = fabs(Bmat)
     if tmp.max() / tmp1.max() > 1e-10:
@@ -866,7 +911,6 @@ def warp_f_gamma(time, f, gam):
     :return f_temp: warped srsf
 
     """
-    M = gam.size
     f_temp = interp((time[-1] - time[0]) * gam + time[0], time, f)
 
     return f_temp