Update proxTV to easily call it from C

.gitignore

@@ -14,3 +14,11 @@ prox_tv.egg-info
 .coverage
 .idea
 wheelhouse
+*.so
+*.dSYM/
+main
+main.cpp   
+*.o            
+*.DS_Store
+*.dylib
+

src/LPopt.cpp

@@ -12,6 +12,10 @@
 #include <limits.h>
 #include "LPopt.h"
 
+#ifdef __cplusplus
+extern "C" {
+#endif
+
 /*** Internal functions headers ***/
 double PN_LPpGap(double *x, double *y, double *diff, int n, double q, double lambda, double norm);
 
@@ -208,8 +212,10 @@ int PN_LPinf(double *y,double lambda,double *x,double *info,int n,Workspace *ws)
         - ws: workspace of allocated memory to use. If NULL, any needed memory is locally managed.
         - positive: 1 if all inputs y >= 0, 0 else.
         - objGap: desired quality of the solution in terms of duality gap.
+        - ctx_ptr: pointer to context data to be passed to the callback function.
+        - callback: callback function to be called at each iteration.
 */
-int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspace *ws,int positive,double objGap){
+int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspace *ws,int positive,double objGap, void* ctx_ptr, int (*callback)(const double* s_ptr, size_t s_length, double delta_k, void* ctx_ptr)){
     double *g=NULL,*d=NULL,*xnorm=NULL,*auxv=NULL;
     double stop,stop2,q,nx,f,fupdate,aux,c,den,xp1vGrad,gRd,delta,prevDelta,improve,rhs,grad0,gap,epsilon;
     int *inactive=NULL,*signs=NULL;
@@ -243,6 +249,10 @@ int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspa
     /* Compute dual norm q */
     q = 1/(1-1/p);
 
+    /* initialize stopping condition */
+    bool stopping_condition = false;
+
+
     /* Special case where the solution is the trivial x = 0 */
     /* This is bound to happen if ||y||_q <= lambda */
     if ( LPnorm(y, n, q) <= lambda ) {
@@ -285,38 +295,50 @@ int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspa
     if(!g || !d || !xnorm || !auxv || !inactive)
         {CANCEL("out of memory",info)}
 
-    /* As initial points, use an approximation for the closed-form solution */
-
-    /* Find close-form solution for TV-L2 */
-    PN_LP2(y , lambda , x, NULL ,n);
-    /* Measure gap of this solution in terms of TV-Lp norm */
+    /* ADDED NOV 15 BY Nathan Allaire - we need x=0 as a first point, not an approximation from another norm. */
+    memset(x, 0, sizeof(double)*n); // Fill x with zeros
     for ( i = 0 ; i < n ; i++ )
         auxv[i] = x[i] - y[i];
     nx = LPnorm( x, n, p );
-    stop = PN_LPpGap(x, y, auxv, n, q, lambda, nx);
-
-    /* Find close-form solution for TV-L1 or TV-Linf, whichever is closer in norm value */
-    if ( p < 2 )
-        PN_LP1(y , lambda , xnorm, NULL ,n);
-    else
-        PN_LPinf(y , lambda , xnorm, NULL ,n, ws);
-    /* Measure gap of this solution in terms of TV-Lp norm */
-    for ( i = 0 ; i < n ; i++ )
-        auxv[i] = xnorm[i] - y[i];
-    nx = LPnorm( xnorm, n, p );
-    stop2 = PN_LPpGap(xnorm, y, auxv, n, q, lambda, nx);
-
-    #ifdef DEBUG
-        fprintf(DEBUG_FILE,"Starting candidates: p2 = %lf, p%s = %lf\n", stop, p<2 ? "1" : "inf", stop2);
-    #endif
-
-    /* Use as starting point the closed-form solution with smaller dual gap */
-    if ( stop2 < stop ) {
-        /* Copy starting point */
-        memcpy(x, xnorm, sizeof(double)*n);
-        stop = stop2;
-    }
+    stop = PN_LPpGap(x, y, auxv, n, q, lambda, nx); // Compute gap of this solution in terms of TV-Lp norm
+
+    // ADDED NOV 15 - BEGINNING OF COMMENTING LINES
+
+    // /* As initial points, use an approximation for the closed-form solution */
+
+    // /* Find close-form solution for TV-L2 */
+    // PN_LP2(y , lambda , x, NULL ,n);
+    // /* Measure gap of this solution in terms of TV-Lp norm */
+    // for ( i = 0 ; i < n ; i++ )
+    //     auxv[i] = x[i] - y[i];
+    // nx = LPnorm( x, n, p );
+    // stop = PN_LPpGap(x, y, auxv, n, q, lambda, nx);
+
+    // /* Find close-form solution for TV-L1 or TV-Linf, whichever is closer in norm value */
+    // if ( p < 2 )
+    //     PN_LP1(y , lambda , xnorm, NULL ,n);
+    // else
+    //     PN_LPinf(y , lambda , xnorm, NULL ,n, ws);
+    // /* Measure gap of this solution in terms of TV-Lp norm */
+    // for ( i = 0 ; i < n ; i++ )
+    //     auxv[i] = xnorm[i] - y[i];
+    // nx = LPnorm( xnorm, n, p );
+    // stop2 = PN_LPpGap(xnorm, y, auxv, n, q, lambda, nx);
+
+    // #ifdef DEBUG
+    //     fprintf(DEBUG_FILE,"Starting candidates: p2 = %lf, p%s = %lf\n", stop, p<2 ? "1" : "inf", stop2);
+    // #endif
+
+    // /* Use as starting point the closed-form solution with smaller dual gap */
+    // if ( stop2 < stop ) {
+    //     /* Copy starting point */
+    //     memcpy(x, xnorm, sizeof(double)*n);
+    //     stop = stop2;
+    // }
+
+    // ADDED NOV 15 - END OF COMMENTING LINES
     /* If dual gap is good enough, we are finished */
+
     if ( stop < objGap ) {
         FREE
         if (info) {
@@ -355,7 +377,6 @@ int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspa
     for ( i = 0 ; i < n ; i++ )
         if ( x[i] < epsilon )
             x[i] = epsilon;
-
     /* Initial value of the point norm */
     nx = LPnorm(x,n,p);
     /* Compute differences x - y */
@@ -364,7 +385,7 @@ int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspa
 
     /* Projected Newton loop */
     stop = gap = DBL_MAX; iters = 0;
-    for(iters=0 ; stop > STOP_PNLP && iters < MAX_ITERS_PNLP && gap > objGap ; iters++){
+    for(iters=0 ;iters < MAX_ITERS_PNLP && !stopping_condition; iters++){ //  stop > STOP_PNLP && 
         #ifdef DEBUG
             fprintf(DEBUG_FILE,"Iter %d, x=[ ",iters);
             for(i=0;i<n && i<DEBUG_N;i++) fprintf(DEBUG_FILE,"%g ",x[i]);
@@ -448,15 +469,24 @@ int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspa
                 break;
 
             /* Use minimum norm subgradient updating direction */
-            case PNLP_MNSG:
+            case PNLP_MNSG: {
                 /* Compute minimum norm subgradient at x=0 */
                 /* Since the gradient at x=0 is not to be trusted, the inactive constraint detection is ignored */
                 nI = n;
+                double MAX_D = 1e9;
                 for(i=0;i<n;i++){
                     d[i] = - pow(y[i] / lambda, 1 / (p-1));
+
+                    if (isinf(d[i])) {
+                        d[i] = (d[i] > 0) ? MAX_D : -MAX_D;
+                    } else if (fabs(d[i]) > MAX_D) {
+                        d[i] = (d[i] > 0) ? MAX_D : -MAX_D;
+                    }
                     inactive[i] = i;
                 }
                 break;
+            }
+
 
             /* Hessian (Newton) updating direction */
             case PNLP_HESSIAN:
@@ -482,6 +512,7 @@ int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspa
                     xp1vGrad += auxv[j] * g[j];
                 }
 
+
                 #ifdef DEBUG
                     fprintf(DEBUG_FILE,"Iter %d, dDiag=[ ",iters);
                     for(i=0;i<n && i<DEBUG_N;i++) fprintf(DEBUG_FILE,"%lg ",d[i]);
@@ -659,6 +690,35 @@ int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspa
 
         /* Compute dual gap */
         gap = PN_LPpGap(x, y, auxv, n, q, lambda, nx);
+
+        /* ADDED NOV 26 NATHAN ALLAIRE */
+        /* If callback and ctx_ptr are not defined, switch to regular update on dual gap. Else, call the callback function from Julia */
+        if (callback && ctx_ptr) {
+            double *x_signed = (double*)malloc(sizeof(double) * n);
+
+            /* Copy x into x_signed and adjust almost zero entries */
+            for (i = 0; i < n; i++) {
+                x_signed[i] = x[i];  // Copy the value from x
+                if (fabs(x_signed[i]) <= epsilon) {
+                    x_signed[i] = 0;  // Set near-zero values to exact zero
+                }
+            }
+
+            /* Apply signs if input was not all positive */
+            if (!positive) {
+                for (i = 0; i < n; i++) {
+                    x_signed[i] = (signs[i] == -1) ? -x_signed[i] : x_signed[i];
+                }
+            }
+
+            stopping_condition = callback(x_signed, n, gap, ctx_ptr);
+            free(x_signed);
+        } else {
+            stopping_condition = (gap < objGap);
+        } /* END OF NOV 26 ADDITION */
+
+
+
         #ifdef DEBUG
             fprintf(DEBUG_FILE,"Iter %d, stop=%lg, gap=%lg\n",iters,stop,gap);
         #endif
@@ -729,8 +789,8 @@ int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspa
 
     The proximity problem is solved to a default level of accuracy, as given by STOP_GAP_PNLP.
 */
-int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspace *ws,int positive) {
-    return PN_LPp(y, lambda, x, info, n, p, ws, positive, STOP_GAP_PNLP);
+int PN_LPp_v2(double *y,double lambda,double *x,double *info,int n,double p,Workspace *ws,int positive) {
+    return PN_LPp(y, lambda, x, info, n, p, ws, positive, STOP_GAP_PNLP, NULL, NULL);
 }
 
 /** PN_LPpGap
@@ -948,7 +1008,7 @@ int LPp_project(double *y,double lambda,double *x,double *info,int n,double p,Wo
         #endif
 
     /* Invoke Lp prox solver on dual norm */
-    if(!PN_LPp(y,lambda,x,info,n,q,ws,1))
+    if(!PN_LPp(y,lambda,x,info,n,q,ws,1, STOP_GAP_PNLP))
         {CANCEL("error in internal Lp prox solver",info)}
 
     /* Apply Moreau's decomposition to recover primal problem solution */
@@ -1057,3 +1117,7 @@ void solveLinearLP(double *z, int n, double p, double lambda, double *s) {
     for ( i = 0 ; i < n ; i++ )
         s[i] = s[i] / sNorm * lambda;
 }
+
+#ifdef __cplusplus
+}
+#endif

src/LPopt.h

@@ -38,21 +38,27 @@
 #define MIN_STEP_PNLP 1e-10
 
 /*** Module functions ***/
-
+#ifdef __cplusplus
+extern "C" {
+#endif
 /* Lp norm value */
 double LPnorm(double *x, int n, double p);
 
 /* Lp proximity solver */
 int PN_LP1(double *y,double lambda,double *x,double *info,int n);
 int PN_LP2(double *y,double lambda,double *x,double *info,int n);
 int PN_LPinf(double *y,double lambda,double *x,double *info,int n,Workspace *ws);
-int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspace *ws,int positive,double objGap);
-int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspace *ws,int positive);
+int PN_LPp(double *y,double lambda,double *x,double *info,int n,double p,Workspace *ws,int positive,double objGap, void* ctx_ptr = NULL,int (*callback)(const double* s_ptr, size_t s_length, double delta_k, void* ctx_ptr)=NULL);
+int PN_LPp_v2(double *y,double lambda,double *x,double *info,int n,double p,Workspace *ws,int positive);
 /* Lp-ball projections */
 int LP1_project(double *y,double lambda,double *x,int n,Workspace *ws);
 int LPp_project(double *y,double lambda,double *x,double *info,int n,double p,Workspace *ws);
 
 /* Linear LP constrained solver */
 void solveLinearLP(double *z, int n, double p, double lambda, double *s);
 
+#ifdef __cplusplus
+}
+#endif
+
 #endif

src/TV2DWopt.cpp

@@ -20,9 +20,13 @@
 #define LAPACK_ILP64
 
 /* Internal functions */
-void DR_proxDiff(size_t n, double* input, double* output, double* W, Workspace *ws);
-void DR_columnsPass(size_t M, size_t N, double* input, double* output, double* W, Workspace **ws);
-void DR_rowsPass(size_t M, size_t N, double* input, double* output, double* ref, double* W, Workspace **ws);
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+void DR_proxDiff_2DW(size_t n, double* input, double* output, double* W, Workspace *ws);
+void DR_columnsPass_2DW(size_t M, size_t N, double* input, double* output, double* W, Workspace **ws);
+void DR_rowsPass_2DW(size_t M, size_t N, double* input, double* output, double* ref, double* W, Workspace **ws);
 
 /* DR2L1W_TV
  *
@@ -102,7 +106,7 @@ int DR2L1W_TV(size_t M, size_t N, double*unary, double*W1, double*W2, double*s,
         fprintf(DEBUG_FILE,"Dual projection along columns\n"); fflush(DEBUG_FILE);
     #endif
     // Projection (prox) step
-    DR_columnsPass(M, N, t, s, W1, ws);
+    DR_columnsPass_2DW(M, N, t, s, W1, ws);
     // Reflection
     for (i=0; i < M*N; i++) s[i] = 2*s[i] - t[i];
 
@@ -112,7 +116,7 @@ int DR2L1W_TV(size_t M, size_t N, double*unary, double*W1, double*W2, double*s,
         fprintf(DEBUG_FILE,"Dual projection along rows\n"); fflush(DEBUG_FILE);
     #endif
     // Projection (prox) step, taking into account displacemente from reference unary signal
-    DR_rowsPass(M, N, s, tb, unary, W2, ws);
+    DR_rowsPass_2DW(M, N, s, tb, unary, W2, ws);
     // Reflection
     for (i=0; i < M*N; i++) tb[i] = -2*tb[i] - s[i];
 
@@ -121,8 +125,8 @@ int DR2L1W_TV(size_t M, size_t N, double*unary, double*W1, double*W2, double*s,
   }
 
   // DR is divergent, but with an additional projection we can recover valid solutions
-  DR_columnsPass(M, N, t, s, W1, ws);
-  DR_rowsPass(M, N, s, tb, unary, W2, ws);
+  DR_columnsPass_2DW(M, N, t, s, W1, ws);
+  DR_rowsPass_2DW(M, N, s, tb, unary, W2, ws);
   for (i = 0; i < M*N; i++) s[i] = - s[i] - tb[i];
 
     /* Gather output information */
@@ -149,7 +153,7 @@ int DR2L1W_TV(size_t M, size_t N, double*unary, double*W1, double*W2, double*s,
     @param W regularization weight along cols, as a matrix of size (M-1)xN
     @param ws array of Workspaces to use for the computation
 */
-void DR_columnsPass(size_t M, size_t N, double* input, double* output, double* W, Workspace **ws) {
+void DR_columnsPass_2DW(size_t M, size_t N, double* input, double* output, double* W, Workspace **ws) {
     #pragma omp parallel shared(M,N,input,output,W,ws) default(none)
     {
         int i,j;
@@ -170,7 +174,7 @@ void DR_columnsPass(size_t M, size_t N, double* input, double* output, double* W
             // Prepare inputs
             memcpy(wsi->in, input+(M*j), sizeof(double)*M);
             // Compute prox difference for this column
-            DR_proxDiff(M, wsi->in, wsi->out, wline, wsi);
+            DR_proxDiff_2DW(M, wsi->in, wsi->out, wline, wsi);
             // Save output
             memcpy(output+(M*j), wsi->out, sizeof(double)*M);
         }
@@ -188,7 +192,7 @@ void DR_columnsPass(size_t M, size_t N, double* input, double* output, double* W
     @param W regularization weight along rows, as a matrix of size Mx(N-1)
     @param ws array of Workspaces to use for the computation
 */
-void DR_rowsPass(size_t M, size_t N, double* input, double* output, double* ref, double* W, Workspace **ws) {
+void DR_rowsPass_2DW(size_t M, size_t N, double* input, double* output, double* ref, double* W, Workspace **ws) {
     #pragma omp parallel shared(M,N,input,ref,output,W,ws) default(none)
     {
         int i,j;
@@ -212,7 +216,7 @@ void DR_rowsPass(size_t M, size_t N, double* input, double* output, double* ref,
             for ( idx = j, i = 0 ; i < N ; i++, idx+=M )
                 wsi->in[i] = ref[idx]-input[idx];
             // Compute prox difference for this row
-            DR_proxDiff(N, wsi->in, wsi->out, wline, wsi);
+            DR_proxDiff_2DW(N, wsi->in, wsi->out, wline, wsi);
             // Save output, recovering displacement from reference signal
             for ( idx = j, i = 0 ; i < N ; i++, idx+=M )
                 output[idx] = wsi->out[i] - ref[idx];
@@ -229,7 +233,7 @@ void DR_rowsPass(size_t M, size_t N, double* input, double* output, double* ref,
  @param W weights of the TV regularization
  @param ws Workspace to use for the computation
  */
-void DR_proxDiff(size_t n, double* input, double* output, double* W, Workspace *ws) {
+void DR_proxDiff_2DW(size_t n, double* input, double* output, double* W, Workspace *ws) {
     int i;
 
     // Compute proximity
@@ -239,3 +243,6 @@ void DR_proxDiff(size_t n, double* input, double* output, double* W, Workspace *
       output[i] = input[i] - output[i];
 }
 
+#ifdef __cplusplus
+}
+#endif