[Teem-commit] SF.net SVN: teem:[7430] teem/trunk/src/ten/tendLmdemo .c

[<kin...@us...>]

Revision: 7430
          http://sourceforge.net/p/teem/code/7430
Author:   kindlmann
Date:     2025-09-14 09:28:38 +0000 (Sun, 14 Sep 2025)
Log Message:
-----------
basically done with self-contained instructive example of using the levmar functions

Modified Paths:
--------------
    teem/trunk/src/ten/tendLmdemo.c

Modified: teem/trunk/src/ten/tendLmdemo.c
===================================================================
--- teem/trunk/src/ten/tendLmdemo.c	2025-09-13 20:12:21 UTC (rev 7429)
+++ teem/trunk/src/ten/tendLmdemo.c	2025-09-14 09:28:38 UTC (rev 7430)
@@ -45,34 +45,6 @@
           "GNUmake-ing with environment variable TEEM_LEVMAR.");
 #endif
 
-/* excerpts from levmar.h:
-
--- double precision LM, with & without Jacobian --
--- unconstrained minimization --
-extern int dlevmar_der(
-      void (*func)(double *p, double *hx, int m, int n, void *adata),
-      void (*jacf)(double *p, double *j, int m, int n, void *adata),
-      double *p, double *x, int m, int n, int itmax, double *opts,
-      double *info, double *work, double *covar, void *adata);
-
-extern int dlevmar_dif(
-      void (*func)(double *p, double *hx, int m, int n, void *adata),
-      double *p, double *x, int m, int n, int itmax, double *opts,
-      double *info, double *work, double *covar, void *adata);
-
--- box-constrained minimization --
-extern int dlevmar_bc_der(
-       void (*func)(double *p, double *hx, int m, int n, void *adata),
-       void (*jacf)(double *p, double *j, int m, int n, void *adata),
-       double *p, double *x, int m, int n, double *lb, double *ub, double *dscl,
-       int itmax, double *opts, double *info, double *work, double *covar, void *adata);
-
-extern int dlevmar_bc_dif(
-       void (*func)(double *p, double *hx, int m, int n, void *adata),
-       double *p, double *x, int m, int n, double *lb, double *ub, double *dscl,
-       int itmax, double *opts, double *info, double *work, double *covar, void *adata);
-*/
-
 static double
 polyEval(const double *parmCurr, unsigned int M, double xx) {
   double ret = 0;
@@ -94,7 +66,7 @@
 
 /*
 funcPredict: the "func" callback for levmar functions
-Forward modeling or prediction callback for both dlevmar_der and dlevmar_dif, which
+Forward modeling or prediction, for all of dlevmar{_bc,}_{der,dif}, which
 the https://users.ics.forth.gr/~lourakis/levmar/ docs describe as:
    functional relation describing measurements.
    A pc[] \in R^mm yields a \hat{y} \in  R^nn
@@ -108,9 +80,9 @@
 funcPredict(double *parmCurr, double *hy, int mm, int nn, void *adata) {
   bagOstate *bag = (bagOstate *)adata;
   if (bag->verb > 2) {
-    printf("%s: called at parmCurr =", __func__);
+    printf("%s:  called at parmCurr =", __func__);
     for (int pi = 0; pi < mm; pi++) {
-      printf(" %g", parmCurr[pi]);
+      printf(" %.17g", parmCurr[pi]);
     }
     printf("\n");
   }
@@ -123,14 +95,47 @@
   return;
 }
 
+/* Analytic Jacobian J[i,j] = d(predVal[i]) / d(parm[j])
+ * e.g. polynomial p(x) = a + b*x + c*x*x => m = #parms = 3 => j=0..2
+ *  => i'th row of J = [  1    xi   xi*xi  ]
+ * levmar wants "row-major" => rows are contiguous in memory
+ * => j is faster axis, and i is slower axis
+ */
+static void
+funcJacobian(double *parmCurr, double *jac, int mm, int nn, void *adata) {
+  bagOstate *bag = (bagOstate *)adata;
+  if (bag->verb > 2) {
+    printf("%s: called at parmCurr =", __func__);
+    for (int pi = 0; pi < mm; pi++) {
+      printf(" %.17g", parmCurr[pi]);
+    }
+    printf("\n");
+  }
+  assert(AIR_UINT(mm) == bag->M);
+  assert(AIR_UINT(nn) == bag->N);
+  for (int si = 0; si < nn; si++) { // si = sample index
+    double sx = NRRD_CELL_POS(bag->xmm[0], bag->xmm[1], nn, si);
+    double pxj = 1; // pow(x,j) = x^j but computed incrementally
+    for (int jj = 0; jj < mm; jj++) {
+      jac[si * mm + jj] = pxj; // jj faster, si slower
+      /* If using this busted Jacobian: jac[jj * nn + si] = pxj; levmar still converges
+      (at least it can converge due to small Dp) but with a huge residual.  There does
+      NOT seem to be any internal levmar sanity check that the Jacobian callback is
+      giving results consistent with the data prediction callback */
+      pxj *= sx;
+    }
+  }
+  return;
+}
+
 // Synthesize data and store in ndata
 static void
-dataSynth(Nrrd *ndata, const bagOstate *bag, double rnStdv, unsigned int rnSeed,
+dataSynth(Nrrd *ndata, const bagOstate *bag, const double rndStdvSeed[2],
           airArray *mop) {
   double *data = AIR_CAST(double *, ndata->data);
   airJSFRand *jsf = NULL;
-  if (rnStdv > 0) {
-    jsf = airJSFRandNew(rnSeed);
+  if (rndStdvSeed[0] > 0) {
+    jsf = airJSFRandNew(AIR_UINT(rndStdvSeed[1]));
     airMopAdd(mop, jsf, (airMopper)airJSFRandNix, airMopAlways);
   }
   unsigned int N = AIR_UINT(ndata->axis[0].size); // N = # of data points
@@ -138,11 +143,11 @@
   for (unsigned int si = 0; si < N; si++) { // si = sample index
     double sx = NRRD_CELL_POS(bag->xmm[0], bag->xmm[1], N, si);
     double val = polyEval(bag->tp, bag->M, sx);
-    double noise = jsf ? rnStdv * airJSFRandNormal_d(jsf) : 0;
+    double noise = jsf ? rndStdvSeed[0] * airJSFRandNormal_d(jsf) : 0;
     /*
     if (bag->verb > 1) {
-      printf("%s: data[%u] = poly(%g) + noise = %g + %g = %g\n", __func__, si, sx, val,
-             noise, val + noise);
+      printf("%s: data[%u] = poly(%.17g) + noise = %.17g + %.17g = %.17g\n", __func__,
+    si, sx, val, noise, val + noise);
     }
     */
     data[si] = val + noise;
@@ -157,68 +162,162 @@
 To be instructive, this levmar wrapper function is trying to be general purpose;
 the fact that funcPredict is evaluating a polynomial and that we're here to fit a
 polynomial should not be exposed by this function.
+
+Summary of info from from https://users.ics.forth.gr/~lourakis/levmar/
+and https://users.ics.forth.gr/~lourakis/levmar/levmar.pdf
+and the levmar-2.6 source code, specifically lm_core.c and lmbc_core.c
+
+notation:
+e in R^n = epsilon vector = (given data x  -  predicted data hx or \hat{x})
+p in R^m = parameter vector
+Dp in R^m = delta p = last step taken in parameter space
+J nxm matrix = Jacobian of i'th predicted data w.r.t j'th parm; J_ij = d hx_i / d p_j
+mu in R = damping term added to diagonal of [J^T J] to normalize it
+   (mu set to tau * max[J^T J]_ii)
+lb, ub = lower (per-parameter) bounds, upper bounds
+
+double opts[5] for controlling levmar functions:
+opts[0] = (\tau) the scale factor for initial \mu  (default LM_INIT_MU)
+(GLK renamed these eps from "epsilon" to avoid confusion with epsilon vector e)
+opts[1] = (\eps1) stopping threshold for ||J^T e||_inf (default LM_STOP_THRESH)
+opts[2] = (\eps2) stopping threshold for ||Dp||_2 (default LM_STOP_THRESH)
+opts[3] = (\eps3) stopping threshold for ||e||_2 (default LM_STOP_THRESH)
+(opts[4] is only for dlevmar_dif(), and it is not a stopping threshold;
+ dlevmar_der() only reads opts[0..3])
+opts[4] = (\delta) step size used in difference approximation of Jacobian.
+If \delta<0, the Jacobian is approximated with central differences which
+are more accurate (but slower!) compared to the forward differences
+(as employed by default). (default LM_DIFF_DELTA)
+Set to opts=NULL for defaults to be used (listed above)
+
+GLK's understanding of these functions is also in the `if (bag->verb)` diagnostics
+about the info[LM_INFO_SZ] vector and interpreting the stopping info[6]
 */
+#if TEEM_LEVMAR
+// these constants come from levmar.h
+static double optsDefault[5] = {LM_INIT_MU,     // 0  (from looking at lm_core.c)
+                                LM_STOP_THRESH, // 1
+                                LM_STOP_THRESH, // 2
+                                LM_STOP_THRESH, // 3
+                                LM_DIFF_DELTA}; // 4
+#else
+static double optsDefault[5] = {0, 0, 0, 0, 0};
+#endif
+
 static void
-doLM(Nrrd *ndata, int itmax, /* const */ double opts[5], bagOstate *bag) {
-  int iters = -1;
+levmarCall(Nrrd *ndata, int ajac, int itmax,
+           /* const (levmar sadly has no const correctness) */ double *bndLo,
+           /* const */ double *bndHi,
+           /* const */ double opts[5], bagOstate *bag) {
+  int IT = -1;     // short variable name just for formatting
+  double info[10]; // == LM_INFO_SZ, but hard-coding so non-LEVMAR version compiles
+  // NaN-out the output info[] vector
+  for (unsigned int ii = 0; ii < 10; ii++) {
+    info[ii] = AIR_NAN;
+  }
+  int M = (int)(bag->M);
+  int N = (int)(bag->N);
+  assert(!!bndLo == !!bndHi);           // either have both bndLo and bndHi or neither
+  double *ym = (double *)(ndata->data); // measured y values
   AIR_UNUSED(ndata);
   AIR_UNUSED(itmax);
-  double info[LM_INFO_SZ];
-  for (unsigned int ii = 0; ii < LM_INFO_SZ; ii++) {
-    info[ii] = AIR_NAN;
-  }
+  AIR_UNUSED(M);
+  AIR_UNUSED(N);
+  AIR_UNUSED(ym);
+  char lmfname[AIR_STRLEN_SMALL + 1];
+  if (ajac) {     // analytic Jacobian
+    if (!bndLo) { // no lower bounds or upper bounds
+      strcpy(lmfname, "dlevmar_der");
+      /* extern int dlevmar_der(
+         void (*func)(double *p, double *hx, int m, int n, void *adata),
+         void (*jacf)(double *p, double *j, int m, int n, void *adata),
+         double *p, double *x, int m, int n, int itmax, double *opts,
+         double *info, double *work, double *covar, void *adata); */
 #if TEEM_LEVMAR
-  double *ym = (double *)(ndata->data); // measured y values
-  /* extern int dlevmar_dif(
+      IT = dlevmar_der(funcPredict,                          //
+                       funcJacobian,                         //
+                       bag->parmCurr, ym, M, N, itmax, opts, //
+                       info, NULL, NULL, AIR_VOIDP(bag));
+#endif
+    } else { // bndLo and bndHi form _b_ox _c_constraints -> bc
+      strcpy(lmfname, "dlevmar_bc_der");
+      /* Why do the bc functions need dscl?  GLK not sure. Quoting from lmbc_core.c:
+      "The algorithm implemented by this function employs projected gradient steps. Since
+      steepest descent is very sensitive to poor scaling, diagonal scaling has been
+      implemented through the dscl argument: Instead of minimizing f(p) for p, f(D*q) is
+      minimized for q=D^-1*p, D being a diagonal scaling matrix whose diagonal equals
+      dscl (see Nocedal-Wright p.27). dscl should contain "typical" magnitudes for the
+      parameters p. A NULL value for dscl implies no scaling. i.e. D=I. To account for
+      scaling, the code divides the starting point and box bounds pointwise by dscl.
+      Moreover, before calling func and jacf the scaling has to be undone (by
+      multiplying), as should be done with the final point. Note also that jac_q=jac_p*D,
+      where jac_q, jac_p are the jacobians w.r.t. q & p, resp."
+      (but nothing there about dscl is specific to having box constraints?) */
+      /* extern int dlevmar_bc_der(
       void (*func)(double *p, double *hx, int m, int n, void *adata),
-      double *p, double *x, int m, int n, int itmax, double *opts,
-      double *info, double *work, double *covar, void *adata); */
-  iters = dlevmar_dif(funcPredict,                                                  //
-                      bag->parmCurr, ym, (int)(bag->M), (int)(bag->N), itmax, opts, //
-                      info, NULL, NULL, AIR_VOIDP(bag));
-#else
-  printf("%s: not calling dlevmar_dif() because we don't have TEEM_LEVMAR\n", __func__);
+      void (*jacf)(double *p, double *j, int m, int n, void *adata),
+      double *p, double *x, int m, int n, double *lb, double *ub, double *dscl,
+      int itmax, double *opts, double *info, double *work, double *covar, void *adata);
+      */
+#if TEEM_LEVMAR
+      IT = dlevmar_bc_der(funcPredict,                                 //
+                          funcJacobian,                                //
+                          bag->parmCurr, ym, M, N, bndLo, bndHi, NULL, //
+                          itmax, opts, info, NULL, NULL, AIR_VOIDP(bag));
 #endif
-  printf("%s: after %d iters, ended at parmCurr =", __func__, iters);
+    }
+  } else { // !ajac: numerical derivatives of predicted data
+    if (!bndLo) {
+      strcpy(lmfname, "dlevmar_dif");
+      /* extern int dlevmar_dif(
+         void (*func)(double *p, double *hx, int m, int n, void *adata),
+         double *p, double *x, int m, int n, int itmax, double *opts,
+         double *info, double *work, double *covar, void *adata); */
+#if TEEM_LEVMAR
+      IT = dlevmar_dif(funcPredict,                          //
+                       bag->parmCurr, ym, M, N, itmax, opts, //
+                       info, NULL, NULL, AIR_VOIDP(bag));
+#endif
+    } else {
+      strcpy(lmfname, "dlevmar_bc_dif");
+      /* extern int dlevmar_bc_dif(
+      void (*func)(double *p, double *hx, int m, int n, void *adata),
+      double *p, double *x, int m, int n, double *lb, double *ub, double *dscl,
+      int itmax, double *opts, double *info, double *work, double *covar, void *adata);
+      */
+#if TEEM_LEVMAR
+      IT = dlevmar_bc_dif(funcPredict,                                 //
+                          bag->parmCurr, ym, M, N, bndLo, bndHi, NULL, //
+                          itmax, opts, info, NULL, NULL, AIR_VOIDP(bag));
+#endif
+    }
+  }
+#if !(TEEM_LEVMAR)
+  printf("%s: Did NOT call %s() because we don't have TEEM_LEVMAR\n", __func__, lmfname);
+#endif
+  printf("%s: After %d iters, %s ended at parmCurr =\n", __func__, IT, lmfname);
   for (unsigned int pi = 0; pi < bag->M; pi++) {
-    printf(" %g", bag->parmCurr[pi]);
+    printf(" %.17g", bag->parmCurr[pi]);
   }
   printf("\n");
   if (bag->verb) {
-    /*
-    Summary of info from from https://users.ics.forth.gr/~lourakis/levmar/
-    and https://users.ics.forth.gr/~lourakis/levmar/levmar.pdf
-    e in R^n = epsilon vector = given data - predicted data
-    p in R^m = parameter vector
-    dp in R^m = delta p = last step taken in parameter space
-    J nxm matrix = Jacobian; J_ij = d hx_i / d p_j
-    mu in R = damping term added to diagonal of [J^T J] to normalize it
-       (mu set to tau * max[J^T J]_ii)
-
-    double opts[5] for controlling levmar functions:
-    opts[0-4] =  [\tau, \epsilon1, \epsilon2, \epsilon3, \delta]. Respectively
-    opts[0] = (\tau) the scale factor for initial \mu,
-    opts[1] = (\epsilon1) stopping threshold for ||J^T e||_inf,
-    opts[2] = (\epsilon2) stopping threshold for ||Dp||_2
-    opts[3] = (\epsilon3) stopping threshold for ||e||_2
-    opts[4] = (\delta) stopping thresh for step used in difference approximation of
-    Jacobian. If \delta<0, the Jacobian is approximated with central differences which
-    are more accurate (but slower!) compared to the forward differences employed by
-    default. Set to opts=NULL for defaults to be used.
-    */
-    printf("info[0]: Error at initial p (||e||^2): %g\n", info[0]);
-    printf("info[1]: Error at final p (||e||^2): %g\n", info[1]);
-    printf("info[2]: Gradient norm at final p (||J^T e||_inf): %g\n", info[2]);
-    printf("info[3]: Norm of last step (||dp||_2): %g\n", info[3]);
-    printf("info[4]: Final damping scaling tau (mu/max[J^T J]_ii): %g\n", info[4]);
+    printf("info[0]: Chi^2 error at initial p (||e||^2): %.17g\n", info[0]);
+    printf("info[1]: Chi^2 error at final p (||e||^2): %.17g\n", info[1]);
+    printf("info[2]: Gradient norm at final p (||J^T e||_inf): %.17g\n", info[2]);
+    printf("info[3]: Norm of last step (||Dp||_2): %.17g\n", info[3]);
+    printf("info[4]: Final damping scaling tau (mu/max[J^T J]_ii): %.17g\n", info[4]);
     printf("info[5]: # of iterations: %.0f\n", info[5]);
     printf("info[6]: Termination reason: ");
-    switch ((int)(info[6] + 0.5)) { // round to nearest int
+    char obuff[AIR_STRLEN_SMALL + 1];
+#define OPTS(I)                                                                         \
+  (opts ? (sprintf(obuff, "%.17g", opts[I]), obuff)                                     \
+        : (sprintf(obuff, "default %.17g", optsDefault[I]), obuff))
+    switch ((int)(info[6])) {
     case 1:
-      printf("(1) stopped by small gradient J^T e\n");
+      printf("(1) stopped by small gradient ||J^T e||_inf < %s = eps1\n", OPTS(1));
       break;
     case 2:
-      printf("(2) stopped by small dp\n");
+      printf("(2) stopped by small parameter step ||Dp||_2 < %s = eps2\n", OPTS(2));
       break;
     case 3:
       printf("(3) stopped by itmax %d\n", itmax);
@@ -231,7 +330,7 @@
       printf("(5) no further error reduction is possible. Restart with increased mu\n");
       break;
     case 6:
-      printf("(6) stopped by small ||e||_2\n");
+      printf("(6) stopped by small residual ||e||_2 < %s = eps3\n", OPTS(3));
       break;
     case 7:
       printf("(7) stopped by invalid (i.e. NaN or Inf) func values; a user error\n");
@@ -255,32 +354,73 @@
   char *perr, *err;
 
   bagOstate bag[1]; // allocate one bag-o-state; can access member m with bag->m
-  hestOptAdd_Nv_Double(&hopt, "tp", "true poly coeffs", 1, -1, &(bag->tp), "1",
+  hestOptAdd_Nv_Double(&hopt, "tp", "true poly coeffs", 1, -1, &(bag->tp), NULL,
                        "coefficients of (ground-truth) polynomial to sample, to "
                        "synthsize the data to later fit. \"-tp A B C\" means "
                        "A + Bx + Cx^2.  These coefficients are the _M_ "
                        "parameters that the LM method seeks to recover",
                        &(bag->M));
+  // HEY should add hest opt for initial parameter vector
   hestOptAdd_1_UInt(&hopt, "N", "# points", &(bag->N), "42",
                     "How many times to sample polynomial to generate data _N_ data "
-                    "points.  Part of the purpose of this demo is to connect the "
+                    "points. Part of the purpose of this demo is to connect the "
                     "stupid single-letter variable names used in the levmar docs "
-                    "to this concrete example.");
+                    "to this concrete example.  NOTE that how the levmar code and docs "
+                    "use \"n\" vs \"m\" is *flipped* from many other presentations of "
+                    "Levenberg-Marquardt.");
   hestOptAdd_2_Double(&hopt, "xmm", "xmin xmax", bag->xmm, "-1 1",
                       "polynomial will be evaluated at _N_ cell-centered points along "
                       "this interval on X axis");
-  double rnStdv;
-  hestOptAdd_1_Double(&hopt, "r", "noise", &rnStdv, "0",
+  double rndStdvSeed[2];
+  hestOptAdd_2_Double(&hopt, "noise", "stdv seed", rndStdvSeed, "0 67",
                       "Gaussian noise of this stdv will be added to polynomial "
-                      "evaluations to generate data");
-  unsigned int rnSeed;
-  hestOptAdd_1_UInt(&hopt, "e", "rng sEed", &rnSeed, "67",
-                    "Seed for random number generator");
+                      "evaluations to generate data, and the seed value (cast to "
+                      "a uint) will be used to initialized the RNG");
   char *outS;
   hestOptAdd_1_String(&hopt, "o", "data out", &outS, NULL,
                       "If given a filename here, will save out synthetic data");
   int itmax;
-  hestOptAdd_1_Int(&hopt, "im", "itmax", &itmax, "100", "cap on # iterations");
+  hestOptAdd_1_Int(&hopt, "itmax", "# iters", &itmax, "100", "cap on # iterations");
+// https://gcc.gnu.org/onlinedocs/cpp/Stringizing.html#Stringizing
+#if TEEM_LEVMAR
+#  define _str(s) #s
+#  define str(s)  _str(s)
+#else
+#  define str(s) ("nan")
+#endif
+  double tau;
+  hestOptAdd_1_Double(&hopt, "tau", "tau", &tau, str(LM_INIT_MU),
+                      "Initial damping mu is found by multiplying max[J^T J]_ii by this "
+                      "number, called tau");
+  double eps1;
+  hestOptAdd_1_Double(&hopt, "eps1", "thresh", &eps1, str(LM_STOP_THRESH),
+                      "stopping thresh on (Linf of) the parm gradient [J^T (x-hx)]");
+  double eps2;
+  hestOptAdd_1_Double(&hopt, "eps2", "thresh", &eps2, str(LM_STOP_THRESH),
+                      "stopping thresh on (L2 of) the parm delta Dp");
+  double eps3;
+  hestOptAdd_1_Double(&hopt, "eps3", "thresh", &eps3, str(LM_STOP_THRESH),
+                      "stopping thresh on (L2 of) the (x-hx) residual vector");
+  int ajac;
+  hestOptAdd_Flag(&hopt, "ajac", &ajac,
+                  "use analytic Jacobian of predicted (modeled) data w.r.t parameters");
+  double bc;
+  hestOptAdd_1_Double(&hopt, "bc", "scaling", &bc, "0.0",
+                      "If > 0, then create box constraints around the parameter "
+                      "by first finding pmax the max absolute value of the "
+                      "ground truth parameters, then box is bc*[-pmax,pmax]^M. "
+                      "So to give breathing room want bc well above 1, but can "
+                      "set bc < 1 for testing purposes.");
+  double delta;
+  hestOptAdd_1_Double(&hopt, "delta", "delta", &delta, str(LM_DIFF_DELTA),
+                      "if not using analytic Jacobian, this is the "
+                      "per-parameter delta to use for numerically computing it, "
+                      "via forward (delta>0) or central (delta<0) differences");
+  int nulopt;
+  hestOptAdd_Flag(&hopt, "nulopt", &nulopt,
+                  "Instead of creating the opts[] vector from the previous "
+                  "tau,eps{1,2,3},delta options, just use NULL, which invokes "
+                  "levmar's internal defaults");
   hestOptAdd_1_Int(&hopt, "v", "verbosity", &(bag->verb), "1", "verbosity level");
 
   airArray *mop = airMopNew();
@@ -307,7 +447,7 @@
     return 1;
   }
   // synthesize the data and maybe save it out
-  dataSynth(ndata, bag, rnStdv, rnSeed, mop);
+  dataSynth(ndata, bag, rndStdvSeed, mop);
   if (airStrlen(outS)) {
     NrrdIoState *nio = nrrdIoStateNew();
     airMopAdd(mop, nio, (airMopper)nrrdIoStateNix, airMopAlways);
@@ -321,16 +461,42 @@
     }
   }
 
-  // prepare for LM call
+  // prepare for levmar call
   bag->parmCurr = AIR_MALLOC(bag->M, double);
   assert(bag->parmCurr);
   airMopAdd(mop, bag->parmCurr, airFree, airMopAlways);
+  double *bndLo = NULL, *bndHi = NULL;
+  if (bc > 0) {
+    double pmax = 0;
+    for (unsigned int pi = 0; pi < bag->M; pi++) {
+      pmax = AIR_MAX(pmax, fabs(bag->tp[pi]));
+    }
+    if (!pmax) {
+      pmax = 1;
+    }
+    bndLo = AIR_MALLOC(bag->M, double);
+    assert(bndLo);
+    airMopAdd(mop, bndLo, airFree, airMopAlways);
+    bndHi = AIR_MALLOC(bag->M, double);
+    assert(bndHi);
+    airMopAdd(mop, bndHi, airFree, airMopAlways);
+    double bnd = bc * pmax;
+    for (unsigned int pi = 0; pi < bag->M; pi++) {
+      bndLo[pi] = -bnd;
+      bndHi[pi] = bnd;
+    }
+    if (bag->verb) {
+      printf("%s: bc %g, pmax %g => box constraints [%g,%g]\n", me, bc, pmax, -bnd, bnd);
+    }
+  }
   // initialize parms that will be fitted
   for (unsigned int pi = 0; pi < bag->M; pi++) {
     bag->parmCurr[pi] = 0;
   }
-  AIR_UNUSED(funcPredict); // to quiet warnings if not levmar
-  doLM(ndata, itmax, /* opts */ NULL, bag);
+  AIR_UNUSED(funcPredict);  // to quiet warnings if no TEEM_LEVMAR
+  AIR_UNUSED(funcJacobian); // to quiet warnings if no TEEM_LEVMAR
+  double opts[5] = {tau, eps1, eps2, eps3, delta};
+  levmarCall(ndata, ajac, itmax, bndLo, bndHi, nulopt ? NULL : opts, bag);
 
   airMopOkay(mop);
   return 0;

This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site.