Merge pull request #423 from MethodicalAcceleratorDesign/dev

sync master

Author: ldeniau

src/mad_tpsa.c

@@ -569,20 +569,26 @@ FUN(cycle) (const T *t, idx_t i, ssz_t n, ord_t m_[n], NUM *v_)
 {
   assert(t); DBGFUN(->); DBGTPSA(t);
   const D *d = t->d;
-  i += 1;
-  ensure(0 <= i && i <= d->nc, "index %d out of bounds", i);
+  if (++i == d->nc) { DBGFUN(<-); return -1; }
+  ensure(0 <= i && i < d->nc, "index %d out of bounds", i);
 
   const idx_t *o2i = d->ord2idx;
-  idx_t ni = o2i[MIN(t->hi,d->to)+1];
-  for (i = MAX(i,o2i[t->lo]); i < ni && t->coef[i] == 0; ++i) ;
+  ord_t lo = MAX(t->lo, d->ords[i]);
+  idx_t hi = MIN(t->hi, d->to);
+
+  for (ord_t o = lo; o <= hi ; ++o) {
+    if (!mad_bit_tst(t->nz,o)) continue;
+    for (i = MAX(i,o2i[o]); i < o2i[o+1] && !t->coef[i]; ++i) ;
+    if (i < o2i[o+1] && t->coef[i]) break;
+  }
 
-  if (i >= ni) { DBGFUN(<-); return -1; }
+  if (i >= o2i[hi+1]) { DBGFUN(<-); return -1; }
 
+  if (v_) *v_ = t->coef[i];
   if (m_) {
     ensure(0 <= n && n <= d->nn, "invalid monomial length %d", n);
     mad_mono_copy(n, d->To[i], m_);
   }
-  if (v_) *v_ = t->coef[i];
 
   DBGFUN(<-); return i;
 }
@@ -604,7 +610,7 @@ FUN(geti) (const T *t, idx_t i)
   const D *d = t->d;
   ensure(0 <= i && i < d->nc, "index %d out of bounds", i);
   ord_t o = d->ords[i];
-  NUM ret = t->lo <= o && o <= MIN(t->hi,d->to) ? t->coef[i] : 0;
+  NUM ret = t->lo <= o && o <= MIN(t->hi,d->to) && mad_bit_tst(t->nz,o) ? t->coef[i] : 0;
   DBGFUN(<-); return ret;
 }
 
@@ -617,7 +623,7 @@ FUN(gets) (const T *t, ssz_t n, str_t s)
   idx_t i = mad_desc_idxs(d,n,s);
   ensure(i >= 0, "invalid monomial");
   ord_t o = d->ords[i];
-  NUM ret = t->lo <= o && o <= MIN(t->hi,d->to) ? t->coef[i] : 0;
+  NUM ret = t->lo <= o && o <= MIN(t->hi,d->to) && mad_bit_tst(t->nz,o) ? t->coef[i] : 0;
   DBGFUN(<-); return ret;
 }
 
@@ -629,7 +635,7 @@ FUN(getm) (const T *t, ssz_t n, const ord_t m[n])
   idx_t i = mad_desc_idxm(d,n,m);
   ensure(i >= 0, "invalid monomial");
   ord_t o = d->ords[i];
-  NUM ret = t->lo <= o && o <= MIN(t->hi,d->to) ? t->coef[i] : 0;
+  NUM ret = t->lo <= o && o <= MIN(t->hi,d->to) && mad_bit_tst(t->nz,o) ? t->coef[i] : 0;
   DBGFUN(<-); return ret;
 }
 
@@ -642,7 +648,7 @@ FUN(getsm) (const T *t, ssz_t n, const idx_t m[n])
   idx_t i = mad_desc_idxsm(d,n,m);
   ensure(i >= 0, "invalid monomial");
   ord_t o = d->ords[i];
-  NUM ret = t->lo <= o && o <= MIN(t->hi,d->to) ? t->coef[i] : 0;
+  NUM ret = t->lo <= o && o <= MIN(t->hi,d->to) && mad_bit_tst(t->nz,o) ? t->coef[i] : 0;
   DBGFUN(<-); return ret;
 }
 
@@ -658,7 +664,7 @@ FUN(getv) (const T *t, idx_t i, ssz_t n, NUM v[n])
   const ord_t *ord = d->ords+i;
   const NUM  *coef = t->coef+i;
   for (idx_t j = 0; j < n; ++j)
-    v[j] = t->lo <= ord[j] && ord[j] <= hi ? coef[j] : 0;
+    v[j] = t->lo <= ord[j] && ord[j] <= hi && mad_bit_tst(t->nz,ord[j]) ? coef[j] : 0;
 
   DBGFUN(<-);
 }

src/mad_tpsa_ops.c

@@ -675,13 +675,16 @@ num_t
 FUN(nrm) (const T *a)
 {
   assert(a); DBGFUN(->); DBGTPSA(a);
-
-  num_t nrm = 0;
   ord_t hi  = MIN(a->hi, a->d->to);
+  num_t nrm = 0;
+
   if (mad_bit_hcut(a->nz,hi)) {
     const idx_t *o2i = a->d->ord2idx;
-    for (idx_t i = o2i[a->lo]; i < o2i[hi+1]; ++i)
-      nrm += fabs(a->coef[i]);
+    for (ord_t o = a->lo; o <= hi ; ++o) {
+      if (!mad_bit_tst(a->nz,o)) continue;
+      for (idx_t i = o2i[o]; i < o2i[o+1]; ++i)
+        nrm += fabs(a->coef[i]);
+    }
   }
 
   DBGFUN(<-); return nrm;

src/madl_lsopt.mad

@@ -265,13 +265,13 @@ local function ld_lmdif (env)
   local jstra, jiter, jtau, bisec, rcond, ncond in obj
 
   -- default setup
-  jstra, jiter, jtau, bisec, rcond = jstra or 1, max(1, jiter or 10),
-             max(eps, jtau or 1e-3), bisec or obj.exec and 0 or 3,
-                                     rcond or 1e-12
+  jstra, jiter, jtau, bisec, rcond, ncond = jstra or 1, max(1, jiter or 10),
+                       max(eps, jtau or 1e-3), bisec or obj.exec and 0 or 3,
+                                            rcond or 1e-12, ncond or 0
 
   if env.info >= 3 then
-    printf("nlopt: jstra=%d, jiter=%d, jtau=%g, bisec=%d, rcond=%g\n",
-                   jstra,    jiter,    jtau,    bisec,    rcond)
+    printf("nlopt: jstra=%d, jiter=%d, jtau=%g, bisec=%d, rcond=%g, ncond=%d\n",
+                   jstra,    jiter,    jtau,    bisec,    rcond,    ncond)
   end
 
   -- adjust variables (i.e. adjust x to fit bbox if any)
@@ -420,13 +420,13 @@ local function ld_jacobian (env)
   local jstra, jiter, bisec, rcond, ncond in obj
 
   -- default setup
-  jstra, jiter, bisec, rcond = jstra or 1, max(1, jiter or 10),
-                               bisec or obj.exec and 0 or 3,
-                               rcond or 1e-12
+  jstra, jiter, bisec, rcond, ncond = jstra or 1, max(1, jiter or 10),
+                                      bisec or obj.exec and 0 or 3,
+                                      rcond or 1e-12, ncond or 0
 
   if env.info >= 3 then
-    printf("nlopt: jstra=%d, jiter=%d, bisec=%d, rcond=%g\n",
-                   jstra,    jiter,    bisec,    rcond)
+    printf("nlopt: jstra=%d, jiter=%d, bisec=%d, rcond=%g, ncond=%d\n",
+                   jstra,    jiter,    bisec,    rcond,    ncond)
   end
 
   -- adjust variables (i.e. adjust x to fit bbox if any)

src/madl_match.mad

@@ -560,13 +560,18 @@ local function compute_fgrd (env)
   -- backup states
   backup(var, bak)
 
-  -- last step norm
+  -- estimate dh from last step or current point
   local hn = h:norm()
+  local dh = 1e-4*hn
+  if dh == 0 then dh = 1e-8*x:norm() end
+  if dh == 0 then dh = 1e-10 end
+
+  if trc then printf("nlopt: computing derivatives with dh=%-.5e\n", dh) end
 
   -- Broyden's Jacobian update
   local bro, Bc = env.objective.broyden and prv.c ~= nil and hn ~= 0
   if bro then
-    if trc then printf("nlopt: Broyden's update with |hn|=%-.8e\n", hn) end
+    if trc then printf("nlopt: Broyden's update with |hn|=%-.5e\n", hn) end
     Bc = cjac + (c - prv.c - cjac*h)*h:t()/hn^2
   end
 
@@ -575,12 +580,13 @@ local function compute_fgrd (env)
     if bro and h[iv] >= 0.8*hn then             -- Broyden's rank one update
       cjac:setcol(iv, Bc:getcol(iv))
       if trc then
-        printf("nlopt: Broyden's update for variable %d (%-.8e)\n", iv, h[iv])
+        printf("nlopt: Broyden's update for variable %d (%-.5e)\n", iv, h[iv])
       end
     else                                        -- finite difference required
-      local ih = x[iv] ~= 0 and x[iv]*xstp[iv] or 1e-12
+      local ih = x[iv]*xstp[iv]
+      if ih == 0 then ih = dh end
       x[iv] = x[iv]+ih ; var.fval = fval
-      if x[iv] < xmin[iv] or x[iv] > xmax[iv] then
+      if x[iv] < xmin[iv] or x[iv] > xmax[iv] or ih*xslp[iv] < 0 then
         x[iv], ih = x[iv]-2*ih, -ih   -- take -ih if bbox are violated
       end
 
@@ -814,7 +820,7 @@ local function exec (self)
     local v = variables[i] or {}
     if v.get then var.x0[i] = v.get(env) end
     var.x   [i] = var.x0[i]
-    var.xstp[i] = v.rstp   or variables.rstp  or  1e-8
+    var.xstp[i] = v.rstp   or variables.rstp  or  0
     var.xslp[i] = v.slope  or variables.slope or  0
     var.xtol[i] = v.tol    or variables.tol   or  0
     var.xmin[i] = v.min    or variables.min   or -inf