From 4ab90a68252cc265ae30b9a5d11e7ffaf25c69b9 Mon Sep 17 00:00:00 2001 From: Laurent Deniau Date: Mon, 22 Jan 2024 15:26:30 +0100 Subject: [PATCH 1/5] fix in cycle to take nz into account and avoid returning garbage --- src/mad_tpsa.c | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/mad_tpsa.c b/src/mad_tpsa.c index 00c410b31..94744ac87 100644 --- a/src/mad_tpsa.c +++ b/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; } From aff813eea1e394d169b1400edd35e85c6ddb3903 Mon Sep 17 00:00:00 2001 From: Laurent Deniau Date: Mon, 22 Jan 2024 16:31:21 +0100 Subject: [PATCH 2/5] fix nrm to take into account nz --- src/mad_tpsa_ops.c | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/mad_tpsa_ops.c b/src/mad_tpsa_ops.c index 65882a423..ab33f3bac 100644 --- a/src/mad_tpsa_ops.c +++ b/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; From fe64824714f1b9aa5d06d331d04ef62d4728e4c5 Mon Sep 17 00:00:00 2001 From: Laurent Deniau Date: Tue, 23 Jan 2024 16:18:06 +0100 Subject: [PATCH 3/5] fix nz check in getters --- src/mad_tpsa.c | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/mad_tpsa.c b/src/mad_tpsa.c index 94744ac87..5e37be11f 100644 --- a/src/mad_tpsa.c +++ b/src/mad_tpsa.c @@ -610,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; } @@ -623,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; } @@ -635,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; } @@ -648,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; } @@ -664,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(<-); } From cfad68ced15375e6360836efcb82862c4dab4993 Mon Sep 17 00:00:00 2001 From: Laurent Deniau Date: Wed, 24 Jan 2024 16:07:53 +0100 Subject: [PATCH 4/5] add ncond to debug output --- src/madl_lsopt.mad | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/src/madl_lsopt.mad b/src/madl_lsopt.mad index 6a1aa80d2..070c2c1c0 100644 --- a/src/madl_lsopt.mad +++ b/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) From ead7372c2a0abab925af20485daaf51c5f58327a Mon Sep 17 00:00:00 2001 From: Laurent Deniau Date: Wed, 24 Jan 2024 16:08:31 +0100 Subject: [PATCH 5/5] restore (much) better dh estimate --- src/madl_match.mad | 18 ++++++++++++------ 1 file changed, 12 insertions(+), 6 deletions(-) diff --git a/src/madl_match.mad b/src/madl_match.mad index d30d60b7b..901ff000f 100644 --- a/src/madl_match.mad +++ b/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