From cccfcd134f9ddfc1deda56655fa0d11ddb44dc88 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Tue, 30 Sep 2025 09:33:09 -0500 Subject: [PATCH 01/20] Wrap --- py/dfovec.py | 104 +-------------------------------------------------- 1 file changed, 2 insertions(+), 102 deletions(-) diff --git a/py/dfovec.py b/py/dfovec.py index ac0f537..18a9ff3 100644 --- a/py/dfovec.py +++ b/py/dfovec.py @@ -10,108 +10,8 @@ def set_constants(): y1 = [1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34e0, 2.1e0, 4.39e0] y2 = [1.957e-1, 1.947e-1, 1.735e-1, 1.6e-1, 8.44e-2, 6.27e-2, 4.56e-2, 3.42e-2, 3.23e-2, 2.35e-2, 2.46e-2] y3 = [3.478e4, 2.861e4, 2.365e4, 1.963e4, 1.637e4, 1.372e4, 1.154e4, 9.744e3, 8.261e3, 7.03e3, 6.005e3, 5.147e3, 4.427e3, 3.82e3, 3.307e3, 2.872e3] - y4 = [ - 8.44e-1, - 9.08e-1, - 9.32e-1, - 9.36e-1, - 9.25e-1, - 9.08e-1, - 8.81e-1, - 8.5e-1, - 8.18e-1, - 7.84e-1, - 7.51e-1, - 7.18e-1, - 6.85e-1, - 6.58e-1, - 6.28e-1, - 6.03e-1, - 5.8e-1, - 5.58e-1, - 5.38e-1, - 5.22e-1, - 5.06e-1, - 4.9e-1, - 4.78e-1, - 4.67e-1, - 4.57e-1, - 4.48e-1, - 4.38e-1, - 4.31e-1, - 4.24e-1, - 4.2e-1, - 4.14e-1, - 4.11e-1, - 4.06e-1, - ] - y5 = [ - 1.366e0, - 1.191e0, - 1.112e0, - 1.013e0, - 9.91e-1, - 8.85e-1, - 8.31e-1, - 8.47e-1, - 7.86e-1, - 7.25e-1, - 7.46e-1, - 6.79e-1, - 6.08e-1, - 6.55e-1, - 6.16e-1, - 6.06e-1, - 6.02e-1, - 6.26e-1, - 6.51e-1, - 7.24e-1, - 6.49e-1, - 6.49e-1, - 6.94e-1, - 6.44e-1, - 6.24e-1, - 6.61e-1, - 6.12e-1, - 5.58e-1, - 5.33e-1, - 4.95e-1, - 5.0e-1, - 4.23e-1, - 3.95e-1, - 3.75e-1, - 3.72e-1, - 3.91e-1, - 3.96e-1, - 4.05e-1, - 4.28e-1, - 4.29e-1, - 5.23e-1, - 5.62e-1, - 6.07e-1, - 6.53e-1, - 6.72e-1, - 7.08e-1, - 6.33e-1, - 6.68e-1, - 6.45e-1, - 6.32e-1, - 5.91e-1, - 5.59e-1, - 5.97e-1, - 6.25e-1, - 7.39e-1, - 7.1e-1, - 7.29e-1, - 7.2e-1, - 6.36e-1, - 5.81e-1, - 4.28e-1, - 2.92e-1, - 1.62e-1, - 9.8e-2, - 5.4e-2, - ] + y4 = [ 8.44e-1, 9.08e-1, 9.32e-1, 9.36e-1, 9.25e-1, 9.08e-1, 8.81e-1, 8.5e-1, 8.18e-1, 7.84e-1, 7.51e-1, 7.18e-1, 6.85e-1, 6.58e-1, 6.28e-1, 6.03e-1, 5.8e-1, 5.58e-1, 5.38e-1, 5.22e-1, 5.06e-1, 4.9e-1, 4.78e-1, 4.67e-1, 4.57e-1, 4.48e-1, 4.38e-1, 4.31e-1, 4.24e-1, 4.2e-1, 4.14e-1, 4.11e-1, 4.06e-1] + y5 = [ 1.366e0, 1.191e0, 1.112e0, 1.013e0, 9.91e-1, 8.85e-1, 8.31e-1, 8.47e-1, 7.86e-1, 7.25e-1, 7.46e-1, 6.79e-1, 6.08e-1, 6.55e-1, 6.16e-1, 6.06e-1, 6.02e-1, 6.26e-1, 6.51e-1, 7.24e-1, 6.49e-1, 6.49e-1, 6.94e-1, 6.44e-1, 6.24e-1, 6.61e-1, 6.12e-1, 5.58e-1, 5.33e-1, 4.95e-1, 5.0e-1, 4.23e-1, 3.95e-1, 3.75e-1, 3.72e-1, 3.91e-1, 3.96e-1, 4.05e-1, 4.28e-1, 4.29e-1, 5.23e-1, 5.62e-1, 6.07e-1, 6.53e-1, 6.72e-1, 7.08e-1, 6.33e-1, 6.68e-1, 6.45e-1, 6.32e-1, 5.91e-1, 5.59e-1, 5.97e-1, 6.25e-1, 7.39e-1, 7.1e-1, 7.29e-1, 7.2e-1, 6.36e-1, 5.81e-1, 4.28e-1, 2.92e-1, 1.62e-1, 9.8e-2, 5.4e-2] return c13, c14, c29, c45, v, y1, y2, y3, y4, y5 From cedabb31af47928be965c638adf64903318dcbe1 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Wed, 1 Oct 2025 09:22:39 -0500 Subject: [PATCH 02/20] A first pass, but with errors --- py/dfovec.py | 104 +++++++++++++++- py/dfovec_jax.py | 226 ++++++++++++++++++++++++++++++++++ py/test_jacobians_with_jax.py | 51 ++++++++ 3 files changed, 379 insertions(+), 2 deletions(-) create mode 100644 py/dfovec_jax.py create mode 100644 py/test_jacobians_with_jax.py diff --git a/py/dfovec.py b/py/dfovec.py index 18a9ff3..ac0f537 100644 --- a/py/dfovec.py +++ b/py/dfovec.py @@ -10,8 +10,108 @@ def set_constants(): y1 = [1.4e-1, 1.8e-1, 2.2e-1, 2.5e-1, 2.9e-1, 3.2e-1, 3.5e-1, 3.9e-1, 3.7e-1, 5.8e-1, 7.3e-1, 9.6e-1, 1.34e0, 2.1e0, 4.39e0] y2 = [1.957e-1, 1.947e-1, 1.735e-1, 1.6e-1, 8.44e-2, 6.27e-2, 4.56e-2, 3.42e-2, 3.23e-2, 2.35e-2, 2.46e-2] y3 = [3.478e4, 2.861e4, 2.365e4, 1.963e4, 1.637e4, 1.372e4, 1.154e4, 9.744e3, 8.261e3, 7.03e3, 6.005e3, 5.147e3, 4.427e3, 3.82e3, 3.307e3, 2.872e3] - y4 = [ 8.44e-1, 9.08e-1, 9.32e-1, 9.36e-1, 9.25e-1, 9.08e-1, 8.81e-1, 8.5e-1, 8.18e-1, 7.84e-1, 7.51e-1, 7.18e-1, 6.85e-1, 6.58e-1, 6.28e-1, 6.03e-1, 5.8e-1, 5.58e-1, 5.38e-1, 5.22e-1, 5.06e-1, 4.9e-1, 4.78e-1, 4.67e-1, 4.57e-1, 4.48e-1, 4.38e-1, 4.31e-1, 4.24e-1, 4.2e-1, 4.14e-1, 4.11e-1, 4.06e-1] - y5 = [ 1.366e0, 1.191e0, 1.112e0, 1.013e0, 9.91e-1, 8.85e-1, 8.31e-1, 8.47e-1, 7.86e-1, 7.25e-1, 7.46e-1, 6.79e-1, 6.08e-1, 6.55e-1, 6.16e-1, 6.06e-1, 6.02e-1, 6.26e-1, 6.51e-1, 7.24e-1, 6.49e-1, 6.49e-1, 6.94e-1, 6.44e-1, 6.24e-1, 6.61e-1, 6.12e-1, 5.58e-1, 5.33e-1, 4.95e-1, 5.0e-1, 4.23e-1, 3.95e-1, 3.75e-1, 3.72e-1, 3.91e-1, 3.96e-1, 4.05e-1, 4.28e-1, 4.29e-1, 5.23e-1, 5.62e-1, 6.07e-1, 6.53e-1, 6.72e-1, 7.08e-1, 6.33e-1, 6.68e-1, 6.45e-1, 6.32e-1, 5.91e-1, 5.59e-1, 5.97e-1, 6.25e-1, 7.39e-1, 7.1e-1, 7.29e-1, 7.2e-1, 6.36e-1, 5.81e-1, 4.28e-1, 2.92e-1, 1.62e-1, 9.8e-2, 5.4e-2] + y4 = [ + 8.44e-1, + 9.08e-1, + 9.32e-1, + 9.36e-1, + 9.25e-1, + 9.08e-1, + 8.81e-1, + 8.5e-1, + 8.18e-1, + 7.84e-1, + 7.51e-1, + 7.18e-1, + 6.85e-1, + 6.58e-1, + 6.28e-1, + 6.03e-1, + 5.8e-1, + 5.58e-1, + 5.38e-1, + 5.22e-1, + 5.06e-1, + 4.9e-1, + 4.78e-1, + 4.67e-1, + 4.57e-1, + 4.48e-1, + 4.38e-1, + 4.31e-1, + 4.24e-1, + 4.2e-1, + 4.14e-1, + 4.11e-1, + 4.06e-1, + ] + y5 = [ + 1.366e0, + 1.191e0, + 1.112e0, + 1.013e0, + 9.91e-1, + 8.85e-1, + 8.31e-1, + 8.47e-1, + 7.86e-1, + 7.25e-1, + 7.46e-1, + 6.79e-1, + 6.08e-1, + 6.55e-1, + 6.16e-1, + 6.06e-1, + 6.02e-1, + 6.26e-1, + 6.51e-1, + 7.24e-1, + 6.49e-1, + 6.49e-1, + 6.94e-1, + 6.44e-1, + 6.24e-1, + 6.61e-1, + 6.12e-1, + 5.58e-1, + 5.33e-1, + 4.95e-1, + 5.0e-1, + 4.23e-1, + 3.95e-1, + 3.75e-1, + 3.72e-1, + 3.91e-1, + 3.96e-1, + 4.05e-1, + 4.28e-1, + 4.29e-1, + 5.23e-1, + 5.62e-1, + 6.07e-1, + 6.53e-1, + 6.72e-1, + 7.08e-1, + 6.33e-1, + 6.68e-1, + 6.45e-1, + 6.32e-1, + 5.91e-1, + 5.59e-1, + 5.97e-1, + 6.25e-1, + 7.39e-1, + 7.1e-1, + 7.29e-1, + 7.2e-1, + 6.36e-1, + 5.81e-1, + 4.28e-1, + 2.92e-1, + 1.62e-1, + 9.8e-2, + 5.4e-2, + ] return c13, c14, c29, c45, v, y1, y2, y3, y4, y5 diff --git a/py/dfovec_jax.py b/py/dfovec_jax.py new file mode 100644 index 0000000..6f62c1d --- /dev/null +++ b/py/dfovec_jax.py @@ -0,0 +1,226 @@ +import jax.numpy as jnp + + +def set_constants(): + c13 = 1.3e1 + c14 = 1.4e1 + c29 = 2.9e1 + c45 = 4.5e1 + v = jnp.array([4.0, 2.0, 1.0, 0.5, 0.25, 0.167, 0.125, 0.1, 0.0833, 0.0714, 0.0625]) + y1 = jnp.array([0.14, 0.18, 0.22, 0.25, 0.29, 0.32, 0.35, 0.39, 0.37, 0.58, 0.73, 0.96, 1.34, 2.1, 4.39]) + y2 = jnp.array([0.1957, 0.1947, 0.1735, 0.16, 0.0844, 0.0627, 0.0456, 0.0342, 0.0323, 0.0235, 0.0246]) + y3 = jnp.array([34780., 28610., 23650., 19630., 16370., 13720., 11540., 9744., 8261., 7030., + 6005., 5147., 4427., 3820., 3307., 2872.]) + y4 = jnp.array([0.844, 0.908, 0.932, 0.936, 0.925, 0.908, 0.881, 0.85, 0.818, 0.784, 0.751, + 0.718, 0.685, 0.658, 0.628, 0.603, 0.58, 0.558, 0.538, 0.522, 0.506, 0.49, + 0.478, 0.467, 0.457, 0.448, 0.438, 0.431, 0.424, 0.42, 0.414, 0.411, 0.406]) + y5 = jnp.array([ + 1.366, 1.191, 1.112, 1.013, 0.991, 0.885, 0.831, 0.847, 0.786, 0.725, + 0.746, 0.679, 0.608, 0.655, 0.616, 0.606, 0.602, 0.626, 0.651, 0.724, + 0.649, 0.649, 0.694, 0.644, 0.624, 0.661, 0.612, 0.558, 0.533, 0.495, + 0.5, 0.423, 0.395, 0.375, 0.372, 0.391, 0.396, 0.405, 0.428, 0.429, + 0.523, 0.562, 0.607, 0.653, 0.672, 0.708, 0.633, 0.668, 0.645, 0.632, + 0.591, 0.559, 0.597, 0.625, 0.739, 0.71, 0.729, 0.72, 0.636, 0.581, + 0.428, 0.292, 0.162, 0.098, 0.054 + ]) + return c13, c14, c29, c45, v, y1, y2, y3, y4, y5 + + +def dfovec_jax(m, n, x, nprob): + """ + This is a Python translation of the Matlab version of the subroutine dfovec.f + This subroutine specifies the nonlinear benchmark problems in + + Benchmarking Derivative-Free Optimization Algorithms + Jorge J. More' and Stefan M. Wild + SIAM J. Optimization, Vol. 20 (1), pp.172-191, 2009. + + The latest version of this subroutine is always available at + https://github.com/POptUS/BenDFO/ + The authors would appreciate feedback and experiences from numerical + studies conducted using this subroutine. + + The data file dfo.dat defines suitable values of m and n + for each problem number nprob. + + This subroutine defines the functions of 22 nonlinear + least squares problems. The allowable values of (m,n) for + functions 1,2 and 3 are variable but with m .ge. n. + For functions 4,5,6,7,8,9 and 10 the values of (m,n) are + (2,2),(3,3),(4,4),(2,2),(15,3),(11,4) and (16,3), respectively. + Function 11 (Watson) has m = 31 with n usually 6 or 9. + However, any n, n = 2,...,31, is permitted. + Functions 12,13 and 14 have n = 3,2 and 4, respectively, but + allow any m .ge. n, with the usual choices being 10,10 and 20. + Function 15 (Chebyquad) allows m and n variable with m .ge. n. + Function 16 (Brown) allows n variable with m = n. + For functions 17 and 18, the values of (m,n) are + (33,5) and (65,11), respectively. + + fvec = dfovec(m, n, x, nprob) + fvec is an output array of length m which contains the nprob + function evaluated at x. + m and n are positive integer input variables. n must not + exceed m. + x is an input array of length n. + nprob is a positive integer input variable which defines the + number of the problem. nprob must not exceed 22. + + Argonne National Laboratory + Jorge More' and Stefan Wild. January 2008. + """ + + c13, c14, c29, c45, v, y1, y2, y3, y4, y5 = set_constants() + + # Initialize things + fvec = jnp.zeros(m) + total = 0 + + if nprob == 1: # Linear function - full rank. + total = jnp.sum(x[:n]) + temp = 2 * total / m + 1 + fvec = -temp * jnp.ones(m) + fvec = fvec.at[:n].add(x[:n]) + + elif nprob == 2: + weights = jnp.arange(1, n + 1) + total = jnp.sum(weights * x) + fvec = (jnp.arange(1, m + 1) * total) - 1 + + elif nprob == 3: + weights = jnp.arange(2, n) + total = jnp.sum(weights * x[1:-1]) + fvec = jnp.zeros(m).at[:m - 1].set(jnp.arange(m - 1) * total - 1).at[m - 1].set(-1.0) + + elif nprob == 4: + fvec = jnp.array([10 * (x[1] - x[0] ** 2), 1 - x[0]]) + + elif nprob == 5: + def theta(x0, x1): + angle = jnp.arctan2(x1, x0) / (2 * jnp.pi) + return jnp.where(x0 > 0, angle, jnp.where(x0 < 0, angle + 0.5, jnp.where((x0 == 0) & (x1 == 0), 0.0, 0.25))) + th = theta(x[0], x[1]) + r = jnp.sqrt(x[0]**2 + x[1]**2) + fvec = jnp.array([10 * (x[2] - 10 * th), 10 * (r - 1), x[2]]) + + elif nprob == 6: + fvec = jnp.array([ + x[0] + 10 * x[1], + jnp.sqrt(5.0) * (x[2] - x[3]), + (x[1] - 2 * x[2]) ** 2, + jnp.sqrt(10.0) * (x[0] - x[3]) ** 2 + ]) + + elif nprob == 7: + fvec = jnp.array([ + -c13 + x[0] + ((5 - x[1]) * x[1] - 2) * x[1], + -c29 + x[0] + ((1 + x[1]) * x[1] - c14) * x[1] + ]) + + elif nprob == 8: + i = jnp.arange(15) + tmp1 = i + 1 + tmp2 = 15 - i + tmp3 = jnp.where(i > 7, tmp2, tmp1) + denom = x[1] * tmp2 + x[2] * tmp3 + fvec = y1 - (x[0] + tmp1 / denom) + + elif nprob == 9: + tmp1 = v * (v + x[1]) + tmp2 = v * (v + x[2]) + x[3] + fvec = y2 - x[0] * tmp1 / tmp2 + + elif nprob == 10: + i = jnp.arange(16) + temp = 5 * (i + 1) + c45 + x[2] + tmp1 = x[1] / temp + tmp2 = jnp.exp(tmp1) + fvec = x[0] * tmp2 - y3 + + elif nprob == 11: + def watson_term(i): + div = (i + 1) / c29 + dx = div ** jnp.arange(n) + dx1 = jnp.arange(1, n) * div ** jnp.arange(1, n) + s1 = jnp.sum(dx1 * x[1:]) + s2 = jnp.sum(dx * x) + return s1 - s2**2 - 1 + fvec = jnp.array([watson_term(i) for i in range(29)] + [x[0], x[1] - x[0] ** 2 - 1]) + + elif nprob == 12: + i = jnp.arange(1, m + 1) + tmp1 = i / 10.0 + term = jnp.exp(-tmp1[:, None] * x[:2]) + const = jnp.exp(-i) - jnp.exp(-tmp1) + fvec = term[:, 0] - term[:, 1] + const * x[2] + + elif nprob == 13: + i = jnp.arange(1, m + 1) + fvec = 2 + 2 * i - jnp.exp(i * x[0]) - jnp.exp(i * x[1]) + + elif nprob == 14: + i = jnp.arange(1, m + 1) + temp = i / 5.0 + tmp1 = x[0] + temp * x[1] - jnp.exp(temp) + tmp2 = x[2] + jnp.sin(temp) * x[3] - jnp.cos(temp) + fvec = tmp1 ** 2 + tmp2 ** 2 + + elif nprob == 15: + t = 2 * x - 1 + T = jnp.polynomial.chebyshev.chebvander(t, m - 1).T + coeffs = jnp.mean(T, axis=1) + correction = jnp.array([0. if i % 2 == 0 else 1 / ((i + 1) ** 2 - 1) for i in range(m)]) + fvec = coeffs + correction + + elif nprob == 16: + total = jnp.sum(x) - (n + 1) + prod = jnp.prod(x) + fvec = jnp.concatenate([x[:-1] + total, jnp.array([prod - 1])]) + + elif nprob == 17: + i = jnp.arange(33) + tmp1 = jnp.exp(-x[3] * 10 * i) + tmp2 = jnp.exp(-x[4] * 10 * i) + fvec = y4 - (x[0] + x[1] * tmp1 + x[2] * tmp2) + + elif nprob == 18: + i = jnp.arange(65) + t = i / 10.0 + tmp1 = jnp.exp(-x[4] * t) + tmp2 = jnp.exp(-x[5] * (t - x[8]) ** 2) + tmp3 = jnp.exp(-x[6] * (t - x[9]) ** 2) + tmp4 = jnp.exp(-x[7] * (t - x[10]) ** 2) + fvec = y5 - (x[0] * tmp1 + x[1] * tmp2 + x[2] * tmp3 + x[3] * tmp4) + + elif nprob == 19: + f1 = -4 * x[:n - 4] + 3 + f2 = sum((i + 1) * x[i + j] ** 2 for j, i in enumerate(range(n - 4))) + fvec = jnp.concatenate([f1, f2.reshape(-1)]) + + elif nprob == 20: + fvec = jnp.zeros(n) + fvec = fvec.at[0].set(x[0] - 1.0) + fvec = fvec.at[1:].set(10 * (x[1:] - x[:-1] ** 3)) + + elif nprob == 21: + def mancino_term(i): + j = jnp.arange(n) + v2 = jnp.sqrt(x[i] ** 2 + (i + 1) / (j + 1)) + return 1400 * x[i] + (i - 49) ** 3 + jnp.sum(v2 * (jnp.sin(jnp.log(v2)) ** 5 + jnp.cos(jnp.log(v2)) ** 5)) + fvec = jnp.array([mancino_term(i) for i in range(n)]) + + elif nprob == 22: + fvec = jnp.zeros(8) + fvec = fvec.at[0].set(x[0] + x[1] + 0.69) + fvec = fvec.at[1].set(x[2] + x[3] + 0.044) + fvec = fvec.at[2].set(x[4]*x[0] + x[5]*x[1] - x[6]*x[2] - x[7]*x[3] + 1.57) + fvec = fvec.at[3].set(x[6]*x[0] + x[7]*x[1] + x[4]*x[2] + x[5]*x[3] + 1.31) + fvec = fvec.at[4].set(x[0]*(x[4]**2 - x[6]**2) - 2*x[2]*x[4]*x[6] + x[1]*(x[5]**2 - x[7]**2) - 2*x[3]*x[5]*x[7] + 2.65) + fvec = fvec.at[5].set(x[2]*(x[4]**2 - x[6]**2) + 2*x[0]*x[4]*x[6] + x[3]*(x[5]**2 - x[7]**2) + 2*x[1]*x[5]*x[7] - 2.0) + fvec = fvec.at[6].set(x[0]*x[4]*(x[4]**2 - 3*x[6]**2) + x[2]*x[6]*(x[6]**2 - 3*x[4]**2) + x[1]*x[5]*(x[5]**2 - 3*x[7]**2) + x[3]*x[7]*(x[7]**2 - 3*x[5]**2) + 12.6) + fvec = fvec.at[7].set(x[2]*x[4]*(x[4]**2 - 3*x[6]**2) - x[0]*x[6]*(x[6]**2 - 3*x[4]**2) + x[3]*x[5]*(x[5]**2 - 3*x[7]**2) - x[1]*x[7]*(x[7]**2 - 3*x[5]**2) - 9.48) + + else: + raise NotImplementedError(f"nprob={nprob} not implemented") + + return fvec diff --git a/py/test_jacobians_with_jax.py b/py/test_jacobians_with_jax.py new file mode 100644 index 0000000..c9dd269 --- /dev/null +++ b/py/test_jacobians_with_jax.py @@ -0,0 +1,51 @@ +import numpy as np +import jax +import jax.numpy as jnp + +from dfovec import dfovec +from dfovec_jax import dfovec_jax +from jacobian import jacobian +from dfoxs import dfoxs # for generating X0s + +# Load benchmark problem definitions +dfo_table = np.loadtxt("../data/dfo.dat") + +def compare_jacobians(nprob, m, n, x0): + x0 = np.array(x0) + fvec_np = dfovec(m, n, x0, nprob) + + # Get FD-based Jacobian + J_fd, _ = jacobian(m, n, x0, nprob) # shape (m, n) + + # Get JAX-based Jacobian + def f(x): + return dfovec_jax(m, n, x, nprob) + + J_jax = jax.jacfwd(f)(jnp.array(x0)) # shape (m, n) + + # Compare + diff = np.linalg.norm(J_fd - np.array(J_jax)) + rel_diff = diff / (np.linalg.norm(J_fd) + 1e-12) + + print(f"nprob={nprob:2d}, m={m:3d}, n={n:2d} | abs diff: {diff:.2e}, rel diff: {rel_diff:.2e}") + +# Loop over all benchmark problems and 3 starting points +for (nprob, n, m, factor_power) in dfo_table: + n = int(n) + m = int(m) + nprob = int(nprob) + scale = int(10 ** factor_power) + + for pt in range(3): + if pt == 0: + x0 = dfoxs(n, nprob, scale) + elif pt == 1: + x0 = 0.1 * np.ones(n) + elif pt == 2: + x0 = 0.1 * np.arange(1, n + 1) + + try: + compare_jacobians(nprob, m, n, x0) + except Exception as e: + print(f"nprob={nprob:2d} failed at pt={pt}: {e}") + From 564564a568a5c46fb281d9578ac6891bc4044e2e Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Wed, 1 Oct 2025 09:27:07 -0500 Subject: [PATCH 03/20] black --- py/dfovec_jax.py | 165 ++++++++++++++++++++++++++-------- py/test_jacobians_with_jax.py | 7 +- 2 files changed, 134 insertions(+), 38 deletions(-) diff --git a/py/dfovec_jax.py b/py/dfovec_jax.py index 6f62c1d..488176e 100644 --- a/py/dfovec_jax.py +++ b/py/dfovec_jax.py @@ -9,20 +9,113 @@ def set_constants(): v = jnp.array([4.0, 2.0, 1.0, 0.5, 0.25, 0.167, 0.125, 0.1, 0.0833, 0.0714, 0.0625]) y1 = jnp.array([0.14, 0.18, 0.22, 0.25, 0.29, 0.32, 0.35, 0.39, 0.37, 0.58, 0.73, 0.96, 1.34, 2.1, 4.39]) y2 = jnp.array([0.1957, 0.1947, 0.1735, 0.16, 0.0844, 0.0627, 0.0456, 0.0342, 0.0323, 0.0235, 0.0246]) - y3 = jnp.array([34780., 28610., 23650., 19630., 16370., 13720., 11540., 9744., 8261., 7030., - 6005., 5147., 4427., 3820., 3307., 2872.]) - y4 = jnp.array([0.844, 0.908, 0.932, 0.936, 0.925, 0.908, 0.881, 0.85, 0.818, 0.784, 0.751, - 0.718, 0.685, 0.658, 0.628, 0.603, 0.58, 0.558, 0.538, 0.522, 0.506, 0.49, - 0.478, 0.467, 0.457, 0.448, 0.438, 0.431, 0.424, 0.42, 0.414, 0.411, 0.406]) - y5 = jnp.array([ - 1.366, 1.191, 1.112, 1.013, 0.991, 0.885, 0.831, 0.847, 0.786, 0.725, - 0.746, 0.679, 0.608, 0.655, 0.616, 0.606, 0.602, 0.626, 0.651, 0.724, - 0.649, 0.649, 0.694, 0.644, 0.624, 0.661, 0.612, 0.558, 0.533, 0.495, - 0.5, 0.423, 0.395, 0.375, 0.372, 0.391, 0.396, 0.405, 0.428, 0.429, - 0.523, 0.562, 0.607, 0.653, 0.672, 0.708, 0.633, 0.668, 0.645, 0.632, - 0.591, 0.559, 0.597, 0.625, 0.739, 0.71, 0.729, 0.72, 0.636, 0.581, - 0.428, 0.292, 0.162, 0.098, 0.054 - ]) + y3 = jnp.array([34780.0, 28610.0, 23650.0, 19630.0, 16370.0, 13720.0, 11540.0, 9744.0, 8261.0, 7030.0, 6005.0, 5147.0, 4427.0, 3820.0, 3307.0, 2872.0]) + y4 = jnp.array( + [ + 0.844, + 0.908, + 0.932, + 0.936, + 0.925, + 0.908, + 0.881, + 0.85, + 0.818, + 0.784, + 0.751, + 0.718, + 0.685, + 0.658, + 0.628, + 0.603, + 0.58, + 0.558, + 0.538, + 0.522, + 0.506, + 0.49, + 0.478, + 0.467, + 0.457, + 0.448, + 0.438, + 0.431, + 0.424, + 0.42, + 0.414, + 0.411, + 0.406, + ] + ) + y5 = jnp.array( + [ + 1.366, + 1.191, + 1.112, + 1.013, + 0.991, + 0.885, + 0.831, + 0.847, + 0.786, + 0.725, + 0.746, + 0.679, + 0.608, + 0.655, + 0.616, + 0.606, + 0.602, + 0.626, + 0.651, + 0.724, + 0.649, + 0.649, + 0.694, + 0.644, + 0.624, + 0.661, + 0.612, + 0.558, + 0.533, + 0.495, + 0.5, + 0.423, + 0.395, + 0.375, + 0.372, + 0.391, + 0.396, + 0.405, + 0.428, + 0.429, + 0.523, + 0.562, + 0.607, + 0.653, + 0.672, + 0.708, + 0.633, + 0.668, + 0.645, + 0.632, + 0.591, + 0.559, + 0.597, + 0.625, + 0.739, + 0.71, + 0.729, + 0.72, + 0.636, + 0.581, + 0.428, + 0.292, + 0.162, + 0.098, + 0.054, + ] + ) return c13, c14, c29, c45, v, y1, y2, y3, y4, y5 @@ -90,32 +183,26 @@ def dfovec_jax(m, n, x, nprob): elif nprob == 3: weights = jnp.arange(2, n) total = jnp.sum(weights * x[1:-1]) - fvec = jnp.zeros(m).at[:m - 1].set(jnp.arange(m - 1) * total - 1).at[m - 1].set(-1.0) + fvec = jnp.zeros(m).at[: m - 1].set(jnp.arange(m - 1) * total - 1).at[m - 1].set(-1.0) elif nprob == 4: fvec = jnp.array([10 * (x[1] - x[0] ** 2), 1 - x[0]]) elif nprob == 5: + def theta(x0, x1): angle = jnp.arctan2(x1, x0) / (2 * jnp.pi) return jnp.where(x0 > 0, angle, jnp.where(x0 < 0, angle + 0.5, jnp.where((x0 == 0) & (x1 == 0), 0.0, 0.25))) + th = theta(x[0], x[1]) - r = jnp.sqrt(x[0]**2 + x[1]**2) + r = jnp.sqrt(x[0] ** 2 + x[1] ** 2) fvec = jnp.array([10 * (x[2] - 10 * th), 10 * (r - 1), x[2]]) elif nprob == 6: - fvec = jnp.array([ - x[0] + 10 * x[1], - jnp.sqrt(5.0) * (x[2] - x[3]), - (x[1] - 2 * x[2]) ** 2, - jnp.sqrt(10.0) * (x[0] - x[3]) ** 2 - ]) + fvec = jnp.array([x[0] + 10 * x[1], jnp.sqrt(5.0) * (x[2] - x[3]), (x[1] - 2 * x[2]) ** 2, jnp.sqrt(10.0) * (x[0] - x[3]) ** 2]) elif nprob == 7: - fvec = jnp.array([ - -c13 + x[0] + ((5 - x[1]) * x[1] - 2) * x[1], - -c29 + x[0] + ((1 + x[1]) * x[1] - c14) * x[1] - ]) + fvec = jnp.array([-c13 + x[0] + ((5 - x[1]) * x[1] - 2) * x[1], -c29 + x[0] + ((1 + x[1]) * x[1] - c14) * x[1]]) elif nprob == 8: i = jnp.arange(15) @@ -138,6 +225,7 @@ def theta(x0, x1): fvec = x[0] * tmp2 - y3 elif nprob == 11: + def watson_term(i): div = (i + 1) / c29 dx = div ** jnp.arange(n) @@ -145,6 +233,7 @@ def watson_term(i): s1 = jnp.sum(dx1 * x[1:]) s2 = jnp.sum(dx * x) return s1 - s2**2 - 1 + fvec = jnp.array([watson_term(i) for i in range(29)] + [x[0], x[1] - x[0] ** 2 - 1]) elif nprob == 12: @@ -163,13 +252,13 @@ def watson_term(i): temp = i / 5.0 tmp1 = x[0] + temp * x[1] - jnp.exp(temp) tmp2 = x[2] + jnp.sin(temp) * x[3] - jnp.cos(temp) - fvec = tmp1 ** 2 + tmp2 ** 2 + fvec = tmp1**2 + tmp2**2 elif nprob == 15: t = 2 * x - 1 T = jnp.polynomial.chebyshev.chebvander(t, m - 1).T coeffs = jnp.mean(T, axis=1) - correction = jnp.array([0. if i % 2 == 0 else 1 / ((i + 1) ** 2 - 1) for i in range(m)]) + correction = jnp.array([0.0 if i % 2 == 0 else 1 / ((i + 1) ** 2 - 1) for i in range(m)]) fvec = coeffs + correction elif nprob == 16: @@ -193,7 +282,7 @@ def watson_term(i): fvec = y5 - (x[0] * tmp1 + x[1] * tmp2 + x[2] * tmp3 + x[3] * tmp4) elif nprob == 19: - f1 = -4 * x[:n - 4] + 3 + f1 = -4 * x[: n - 4] + 3 f2 = sum((i + 1) * x[i + j] ** 2 for j, i in enumerate(range(n - 4))) fvec = jnp.concatenate([f1, f2.reshape(-1)]) @@ -203,22 +292,28 @@ def watson_term(i): fvec = fvec.at[1:].set(10 * (x[1:] - x[:-1] ** 3)) elif nprob == 21: + def mancino_term(i): j = jnp.arange(n) v2 = jnp.sqrt(x[i] ** 2 + (i + 1) / (j + 1)) return 1400 * x[i] + (i - 49) ** 3 + jnp.sum(v2 * (jnp.sin(jnp.log(v2)) ** 5 + jnp.cos(jnp.log(v2)) ** 5)) + fvec = jnp.array([mancino_term(i) for i in range(n)]) elif nprob == 22: fvec = jnp.zeros(8) fvec = fvec.at[0].set(x[0] + x[1] + 0.69) fvec = fvec.at[1].set(x[2] + x[3] + 0.044) - fvec = fvec.at[2].set(x[4]*x[0] + x[5]*x[1] - x[6]*x[2] - x[7]*x[3] + 1.57) - fvec = fvec.at[3].set(x[6]*x[0] + x[7]*x[1] + x[4]*x[2] + x[5]*x[3] + 1.31) - fvec = fvec.at[4].set(x[0]*(x[4]**2 - x[6]**2) - 2*x[2]*x[4]*x[6] + x[1]*(x[5]**2 - x[7]**2) - 2*x[3]*x[5]*x[7] + 2.65) - fvec = fvec.at[5].set(x[2]*(x[4]**2 - x[6]**2) + 2*x[0]*x[4]*x[6] + x[3]*(x[5]**2 - x[7]**2) + 2*x[1]*x[5]*x[7] - 2.0) - fvec = fvec.at[6].set(x[0]*x[4]*(x[4]**2 - 3*x[6]**2) + x[2]*x[6]*(x[6]**2 - 3*x[4]**2) + x[1]*x[5]*(x[5]**2 - 3*x[7]**2) + x[3]*x[7]*(x[7]**2 - 3*x[5]**2) + 12.6) - fvec = fvec.at[7].set(x[2]*x[4]*(x[4]**2 - 3*x[6]**2) - x[0]*x[6]*(x[6]**2 - 3*x[4]**2) + x[3]*x[5]*(x[5]**2 - 3*x[7]**2) - x[1]*x[7]*(x[7]**2 - 3*x[5]**2) - 9.48) + fvec = fvec.at[2].set(x[4] * x[0] + x[5] * x[1] - x[6] * x[2] - x[7] * x[3] + 1.57) + fvec = fvec.at[3].set(x[6] * x[0] + x[7] * x[1] + x[4] * x[2] + x[5] * x[3] + 1.31) + fvec = fvec.at[4].set(x[0] * (x[4] ** 2 - x[6] ** 2) - 2 * x[2] * x[4] * x[6] + x[1] * (x[5] ** 2 - x[7] ** 2) - 2 * x[3] * x[5] * x[7] + 2.65) + fvec = fvec.at[5].set(x[2] * (x[4] ** 2 - x[6] ** 2) + 2 * x[0] * x[4] * x[6] + x[3] * (x[5] ** 2 - x[7] ** 2) + 2 * x[1] * x[5] * x[7] - 2.0) + fvec = fvec.at[6].set( + x[0] * x[4] * (x[4] ** 2 - 3 * x[6] ** 2) + x[2] * x[6] * (x[6] ** 2 - 3 * x[4] ** 2) + x[1] * x[5] * (x[5] ** 2 - 3 * x[7] ** 2) + x[3] * x[7] * (x[7] ** 2 - 3 * x[5] ** 2) + 12.6 + ) + fvec = fvec.at[7].set( + x[2] * x[4] * (x[4] ** 2 - 3 * x[6] ** 2) - x[0] * x[6] * (x[6] ** 2 - 3 * x[4] ** 2) + x[3] * x[5] * (x[5] ** 2 - 3 * x[7] ** 2) - x[1] * x[7] * (x[7] ** 2 - 3 * x[5] ** 2) - 9.48 + ) else: raise NotImplementedError(f"nprob={nprob} not implemented") diff --git a/py/test_jacobians_with_jax.py b/py/test_jacobians_with_jax.py index c9dd269..f9e9d9b 100644 --- a/py/test_jacobians_with_jax.py +++ b/py/test_jacobians_with_jax.py @@ -10,6 +10,7 @@ # Load benchmark problem definitions dfo_table = np.loadtxt("../data/dfo.dat") + def compare_jacobians(nprob, m, n, x0): x0 = np.array(x0) fvec_np = dfovec(m, n, x0, nprob) @@ -29,12 +30,13 @@ def f(x): print(f"nprob={nprob:2d}, m={m:3d}, n={n:2d} | abs diff: {diff:.2e}, rel diff: {rel_diff:.2e}") + # Loop over all benchmark problems and 3 starting points -for (nprob, n, m, factor_power) in dfo_table: +for nprob, n, m, factor_power in dfo_table: n = int(n) m = int(m) nprob = int(nprob) - scale = int(10 ** factor_power) + scale = int(10**factor_power) for pt in range(3): if pt == 0: @@ -48,4 +50,3 @@ def f(x): compare_jacobians(nprob, m, n, x0) except Exception as e: print(f"nprob={nprob:2d} failed at pt={pt}: {e}") - From 797fc07d5301ec13c38bccf30769d2cc85d4d6ed Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Wed, 1 Oct 2025 10:39:04 -0500 Subject: [PATCH 04/20] Adding comments --- py/dfovec_jax.py | 42 +++++++++++++++++------------------ py/test_jacobians_with_jax.py | 2 ++ 2 files changed, 23 insertions(+), 21 deletions(-) diff --git a/py/dfovec_jax.py b/py/dfovec_jax.py index 488176e..f5aca0d 100644 --- a/py/dfovec_jax.py +++ b/py/dfovec_jax.py @@ -175,20 +175,20 @@ def dfovec_jax(m, n, x, nprob): fvec = -temp * jnp.ones(m) fvec = fvec.at[:n].add(x[:n]) - elif nprob == 2: + elif nprob == 2: # Linear function - rank 1. weights = jnp.arange(1, n + 1) total = jnp.sum(weights * x) fvec = (jnp.arange(1, m + 1) * total) - 1 - elif nprob == 3: + elif nprob == 3: # Linear function - rank 1 with zero columns and rows. weights = jnp.arange(2, n) total = jnp.sum(weights * x[1:-1]) fvec = jnp.zeros(m).at[: m - 1].set(jnp.arange(m - 1) * total - 1).at[m - 1].set(-1.0) - elif nprob == 4: + elif nprob == 4: # Rosenbrock function. fvec = jnp.array([10 * (x[1] - x[0] ** 2), 1 - x[0]]) - elif nprob == 5: + elif nprob == 5: # Helical valley function. def theta(x0, x1): angle = jnp.arctan2(x1, x0) / (2 * jnp.pi) @@ -198,13 +198,13 @@ def theta(x0, x1): r = jnp.sqrt(x[0] ** 2 + x[1] ** 2) fvec = jnp.array([10 * (x[2] - 10 * th), 10 * (r - 1), x[2]]) - elif nprob == 6: + elif nprob == 6: # Powell singular function. fvec = jnp.array([x[0] + 10 * x[1], jnp.sqrt(5.0) * (x[2] - x[3]), (x[1] - 2 * x[2]) ** 2, jnp.sqrt(10.0) * (x[0] - x[3]) ** 2]) - elif nprob == 7: + elif nprob == 7: # Freudenstein and Roth function. fvec = jnp.array([-c13 + x[0] + ((5 - x[1]) * x[1] - 2) * x[1], -c29 + x[0] + ((1 + x[1]) * x[1] - c14) * x[1]]) - elif nprob == 8: + elif nprob == 8: # Bard function. i = jnp.arange(15) tmp1 = i + 1 tmp2 = 15 - i @@ -212,19 +212,19 @@ def theta(x0, x1): denom = x[1] * tmp2 + x[2] * tmp3 fvec = y1 - (x[0] + tmp1 / denom) - elif nprob == 9: + elif nprob == 9: # Kowalik and Osborne function. tmp1 = v * (v + x[1]) tmp2 = v * (v + x[2]) + x[3] fvec = y2 - x[0] * tmp1 / tmp2 - elif nprob == 10: + elif nprob == 10: # Meyer function. i = jnp.arange(16) temp = 5 * (i + 1) + c45 + x[2] tmp1 = x[1] / temp tmp2 = jnp.exp(tmp1) fvec = x[0] * tmp2 - y3 - elif nprob == 11: + elif nprob == 11: # Watson function. def watson_term(i): div = (i + 1) / c29 @@ -236,43 +236,43 @@ def watson_term(i): fvec = jnp.array([watson_term(i) for i in range(29)] + [x[0], x[1] - x[0] ** 2 - 1]) - elif nprob == 12: + elif nprob == 12: # Box 3-dimensional function. i = jnp.arange(1, m + 1) tmp1 = i / 10.0 term = jnp.exp(-tmp1[:, None] * x[:2]) const = jnp.exp(-i) - jnp.exp(-tmp1) fvec = term[:, 0] - term[:, 1] + const * x[2] - elif nprob == 13: + elif nprob == 13: # Jennrich and Sampson function. i = jnp.arange(1, m + 1) fvec = 2 + 2 * i - jnp.exp(i * x[0]) - jnp.exp(i * x[1]) - elif nprob == 14: + elif nprob == 14: # Brown and Dennis function. i = jnp.arange(1, m + 1) temp = i / 5.0 tmp1 = x[0] + temp * x[1] - jnp.exp(temp) tmp2 = x[2] + jnp.sin(temp) * x[3] - jnp.cos(temp) fvec = tmp1**2 + tmp2**2 - elif nprob == 15: + elif nprob == 15: # Chebyquad function. t = 2 * x - 1 T = jnp.polynomial.chebyshev.chebvander(t, m - 1).T coeffs = jnp.mean(T, axis=1) correction = jnp.array([0.0 if i % 2 == 0 else 1 / ((i + 1) ** 2 - 1) for i in range(m)]) fvec = coeffs + correction - elif nprob == 16: + elif nprob == 16: # Brown almost-linear function. total = jnp.sum(x) - (n + 1) prod = jnp.prod(x) fvec = jnp.concatenate([x[:-1] + total, jnp.array([prod - 1])]) - elif nprob == 17: + elif nprob == 17: # Osborne 1 function. i = jnp.arange(33) tmp1 = jnp.exp(-x[3] * 10 * i) tmp2 = jnp.exp(-x[4] * 10 * i) fvec = y4 - (x[0] + x[1] * tmp1 + x[2] * tmp2) - elif nprob == 18: + elif nprob == 18: # Osborne 2 function. i = jnp.arange(65) t = i / 10.0 tmp1 = jnp.exp(-x[4] * t) @@ -281,17 +281,17 @@ def watson_term(i): tmp4 = jnp.exp(-x[7] * (t - x[10]) ** 2) fvec = y5 - (x[0] * tmp1 + x[1] * tmp2 + x[2] * tmp3 + x[3] * tmp4) - elif nprob == 19: + elif nprob == 19: # Bdqrtic f1 = -4 * x[: n - 4] + 3 f2 = sum((i + 1) * x[i + j] ** 2 for j, i in enumerate(range(n - 4))) fvec = jnp.concatenate([f1, f2.reshape(-1)]) - elif nprob == 20: + elif nprob == 20: # Cube fvec = jnp.zeros(n) fvec = fvec.at[0].set(x[0] - 1.0) fvec = fvec.at[1:].set(10 * (x[1:] - x[:-1] ** 3)) - elif nprob == 21: + elif nprob == 21: # Mancino def mancino_term(i): j = jnp.arange(n) @@ -300,7 +300,7 @@ def mancino_term(i): fvec = jnp.array([mancino_term(i) for i in range(n)]) - elif nprob == 22: + elif nprob == 22: # Heart8ls fvec = jnp.zeros(8) fvec = fvec.at[0].set(x[0] + x[1] + 0.69) fvec = fvec.at[1].set(x[2] + x[3] + 0.044) diff --git a/py/test_jacobians_with_jax.py b/py/test_jacobians_with_jax.py index f9e9d9b..b5b2c48 100644 --- a/py/test_jacobians_with_jax.py +++ b/py/test_jacobians_with_jax.py @@ -29,6 +29,8 @@ def f(x): rel_diff = diff / (np.linalg.norm(J_fd) + 1e-12) print(f"nprob={nprob:2d}, m={m:3d}, n={n:2d} | abs diff: {diff:.2e}, rel diff: {rel_diff:.2e}") + if diff > 1e-7 and rel_diff > 1e-7: + raise Exception("Difference is too large") # Loop over all benchmark problems and 3 starting points From f5caa2e8d3114ead7696a8a76370852e41e5897b Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Wed, 1 Oct 2025 10:41:00 -0500 Subject: [PATCH 05/20] Improving readability --- py/dfovec.py | 23 +++++++++++++++++++++++ 1 file changed, 23 insertions(+) diff --git a/py/dfovec.py b/py/dfovec.py index ac0f537..bf96585 100644 --- a/py/dfovec.py +++ b/py/dfovec.py @@ -174,20 +174,24 @@ def dfovec(m, n, x, nprob): fvec[i] = -temp if i < n: fvec[i] = fvec[i] + x[i] + elif nprob == 2: # Linear function - rank 1. for j in range(n): total = total + (j + 1) * x[j] for i in range(m): fvec[i] = (i + 1) * total - 1 + elif nprob == 3: # Linear function - rank 1 with zero columns and rows. for j in range(1, n - 1): total = total + (j + 1) * x[j] for i in range(m - 1): fvec[i] = i * total - 1 fvec[m - 1] = -1 + elif nprob == 4: # Rosenbrock function. fvec[0] = 10 * (x[1] - x[0] * x[0]) fvec[1] = 1 - x[0] + elif nprob == 5: # Helical valley function. if x[0] > 0: th = np.arctan(x[1] / x[0]) / (2 * np.pi) @@ -201,14 +205,17 @@ def dfovec(m, n, x, nprob): fvec[0] = 10 * (x[2] - 10 * th) fvec[1] = 10 * (r - 1) fvec[2] = x[2] + elif nprob == 6: # Powell singular function. fvec[0] = x[0] + 10 * x[1] fvec[1] = np.sqrt(5) * (x[2] - x[3]) fvec[2] = (x[1] - 2 * x[2]) ** 2 fvec[3] = np.sqrt(10) * (x[0] - x[3]) ** 2 + elif nprob == 7: # Freudenstein and Roth function. fvec[0] = -c13 + x[0] + ((5 - x[1]) * x[1] - 2) * x[1] fvec[1] = -c29 + x[0] + ((1 + x[1]) * x[1] - c14) * x[1] + elif nprob == 8: # Bard function. for i in range(15): tmp1 = i + 1 @@ -217,17 +224,20 @@ def dfovec(m, n, x, nprob): if i > 7: tmp3 = tmp2 fvec[i] = y1[i] - (x[0] + tmp1 / (x[1] * tmp2 + x[2] * tmp3)) + elif nprob == 9: # Kowalik and Osborne function. for i in range(11): tmp1 = v[i] * (v[i] + x[1]) tmp2 = v[i] * (v[i] + x[2]) + x[3] fvec[i] = y2[i] - x[0] * tmp1 / tmp2 + elif nprob == 10: # Meyer function. for i in range(16): temp = 5 * (i + 1) + c45 + x[2] tmp1 = x[1] / temp tmp2 = np.exp(tmp1) fvec[i] = x[0] * tmp2 - y3[i] + elif nprob == 11: # Watson function. for i in range(29): div = (i + 1) / c29 @@ -244,21 +254,25 @@ def dfovec(m, n, x, nprob): fvec[i] = s1 - s2 * s2 - 1 fvec[29] = x[0] fvec[30] = x[1] - x[0] * x[0] - 1 + elif nprob == 12: # Box 3-dimensional function. for i in range(m): temp = i + 1 tmp1 = temp / 10 fvec[i] = np.exp(-tmp1 * x[0]) - np.exp(-tmp1 * x[1]) + (np.exp(-temp) - np.exp(-tmp1)) * x[2] + elif nprob == 13: # Jennrich and Sampson function. for i in range(m): temp = i + 1 fvec[i] = 2 + 2 * temp - np.exp(temp * x[0]) - np.exp(temp * x[1]) + elif nprob == 14: # Brown and Dennis function. for i in range(m): temp = (i + 1) / 5 tmp1 = x[0] + temp * x[1] - np.exp(temp) tmp2 = x[2] + np.sin(temp) * x[3] - np.cos(temp) fvec[i] = tmp1 * tmp1 + tmp2 * tmp2 + elif nprob == 15: # Chebyquad function. for j in range(n): t1 = 1 @@ -275,6 +289,7 @@ def dfovec(m, n, x, nprob): if iev > 0: fvec[i] = fvec[i] + 1 / ((i + 1) ** 2 - 1) iev = -iev + elif nprob == 16: # Brown almost-linear function. total1 = -(n + 1) prod1 = 1 @@ -284,12 +299,14 @@ def dfovec(m, n, x, nprob): for i in range(n - 1): fvec[i] = x[i] + total1 fvec[n - 1] = prod1 - 1 + elif nprob == 17: # Osborne 1 function. for i in range(33): temp = 10 * i tmp1 = np.exp(-x[3] * temp) tmp2 = np.exp(-x[4] * temp) fvec[i] = y4[i] - (x[0] + x[1] * tmp1 + x[2] * tmp2) + elif nprob == 18: # Osborne 2 function. for i in range(65): temp = i / 10 @@ -298,16 +315,19 @@ def dfovec(m, n, x, nprob): tmp3 = np.exp(-x[6] * (temp - x[9]) ** 2) tmp4 = np.exp(-x[7] * (temp - x[10]) ** 2) fvec[i] = y5[i] - (x[0] * tmp1 + x[1] * tmp2 + x[2] * tmp3 + x[3] * tmp4) # noqa + elif nprob == 19: # Bdqrtic # n >= 5, m = (n-4)*2 for i in range(n - 4): fvec[i] = -4 * x[i] + 3.0 fvec[n - 4 + i] = x[i] ** 2 + 2 * x[i + 1] ** 2 + 3 * x[i + 2] ** 2 + 4 * x[i + 3] ** 2 + 5 * x[n - 1] ** 2 + elif nprob == 20: # Cube # n = 2, m = n fvec[0] = x[0] - 1.0 for i in range(1, n): fvec[i] = 10 * (x[i] - x[i - 1] ** 3) + elif nprob == 21: # Mancino # n = 2, m = n for i in range(n): @@ -316,6 +336,7 @@ def dfovec(m, n, x, nprob): v2 = np.sqrt(x[i] ** 2 + (i + 1) / (j + 1)) ss = ss + v2 * ((np.sin(np.log(v2))) ** 5 + (np.cos(np.log(v2))) ** 5) # noqa fvec[i] = 1400 * x[i] + (i - 49) ** 3 + ss + elif nprob == 22: # Heart8ls # m = n = 8 fvec[0] = x[0] + x[1] + 0.69 @@ -330,7 +351,9 @@ def dfovec(m, n, x, nprob): fvec[7] = ( x[2] * x[4] * (x[4] ** 2 - 3.0 * x[6] ** 2) - x[0] * x[6] * (x[6] ** 2 - 3.0 * x[4] ** 2) + x[3] * x[5] * (x[5] ** 2 - 3.0 * x[7] ** 2) - x[1] * x[7] * (x[7] ** 2 - 3.0 * x[5] ** 2) - 9.48 ) + else: print(f"unrecognized function number {nprob}") return None + return fvec From 4375aa41c1fa210ed499d09f6d246e48752f9bf8 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Wed, 1 Oct 2025 10:50:54 -0500 Subject: [PATCH 06/20] jnp to np for easier diffs --- py/dfovec_jax.py | 122 +++++++++++++++++++++++------------------------ 1 file changed, 61 insertions(+), 61 deletions(-) diff --git a/py/dfovec_jax.py b/py/dfovec_jax.py index f5aca0d..11e1272 100644 --- a/py/dfovec_jax.py +++ b/py/dfovec_jax.py @@ -1,4 +1,4 @@ -import jax.numpy as jnp +import jax.numpy as np def set_constants(): @@ -6,11 +6,11 @@ def set_constants(): c14 = 1.4e1 c29 = 2.9e1 c45 = 4.5e1 - v = jnp.array([4.0, 2.0, 1.0, 0.5, 0.25, 0.167, 0.125, 0.1, 0.0833, 0.0714, 0.0625]) - y1 = jnp.array([0.14, 0.18, 0.22, 0.25, 0.29, 0.32, 0.35, 0.39, 0.37, 0.58, 0.73, 0.96, 1.34, 2.1, 4.39]) - y2 = jnp.array([0.1957, 0.1947, 0.1735, 0.16, 0.0844, 0.0627, 0.0456, 0.0342, 0.0323, 0.0235, 0.0246]) - y3 = jnp.array([34780.0, 28610.0, 23650.0, 19630.0, 16370.0, 13720.0, 11540.0, 9744.0, 8261.0, 7030.0, 6005.0, 5147.0, 4427.0, 3820.0, 3307.0, 2872.0]) - y4 = jnp.array( + v = np.array([4.0, 2.0, 1.0, 0.5, 0.25, 0.167, 0.125, 0.1, 0.0833, 0.0714, 0.0625]) + y1 = np.array([0.14, 0.18, 0.22, 0.25, 0.29, 0.32, 0.35, 0.39, 0.37, 0.58, 0.73, 0.96, 1.34, 2.1, 4.39]) + y2 = np.array([0.1957, 0.1947, 0.1735, 0.16, 0.0844, 0.0627, 0.0456, 0.0342, 0.0323, 0.0235, 0.0246]) + y3 = np.array([34780.0, 28610.0, 23650.0, 19630.0, 16370.0, 13720.0, 11540.0, 9744.0, 8261.0, 7030.0, 6005.0, 5147.0, 4427.0, 3820.0, 3307.0, 2872.0]) + y4 = np.array( [ 0.844, 0.908, @@ -47,7 +47,7 @@ def set_constants(): 0.406, ] ) - y5 = jnp.array( + y5 = np.array( [ 1.366, 1.191, @@ -166,49 +166,49 @@ def dfovec_jax(m, n, x, nprob): c13, c14, c29, c45, v, y1, y2, y3, y4, y5 = set_constants() # Initialize things - fvec = jnp.zeros(m) + fvec = np.zeros(m) total = 0 if nprob == 1: # Linear function - full rank. - total = jnp.sum(x[:n]) + total = np.sum(x[:n]) temp = 2 * total / m + 1 - fvec = -temp * jnp.ones(m) + fvec = -temp * np.ones(m) fvec = fvec.at[:n].add(x[:n]) elif nprob == 2: # Linear function - rank 1. - weights = jnp.arange(1, n + 1) - total = jnp.sum(weights * x) - fvec = (jnp.arange(1, m + 1) * total) - 1 + weights = np.arange(1, n + 1) + total = np.sum(weights * x) + fvec = (np.arange(1, m + 1) * total) - 1 elif nprob == 3: # Linear function - rank 1 with zero columns and rows. - weights = jnp.arange(2, n) - total = jnp.sum(weights * x[1:-1]) - fvec = jnp.zeros(m).at[: m - 1].set(jnp.arange(m - 1) * total - 1).at[m - 1].set(-1.0) + weights = np.arange(2, n) + total = np.sum(weights * x[1:-1]) + fvec = np.zeros(m).at[: m - 1].set(np.arange(m - 1) * total - 1).at[m - 1].set(-1.0) elif nprob == 4: # Rosenbrock function. - fvec = jnp.array([10 * (x[1] - x[0] ** 2), 1 - x[0]]) + fvec = np.array([10 * (x[1] - x[0] ** 2), 1 - x[0]]) elif nprob == 5: # Helical valley function. def theta(x0, x1): - angle = jnp.arctan2(x1, x0) / (2 * jnp.pi) - return jnp.where(x0 > 0, angle, jnp.where(x0 < 0, angle + 0.5, jnp.where((x0 == 0) & (x1 == 0), 0.0, 0.25))) + angle = np.arctan2(x1, x0) / (2 * np.pi) + return np.where(x0 > 0, angle, np.where(x0 < 0, angle + 0.5, np.where((x0 == 0) & (x1 == 0), 0.0, 0.25))) th = theta(x[0], x[1]) - r = jnp.sqrt(x[0] ** 2 + x[1] ** 2) - fvec = jnp.array([10 * (x[2] - 10 * th), 10 * (r - 1), x[2]]) + r = np.sqrt(x[0] ** 2 + x[1] ** 2) + fvec = np.array([10 * (x[2] - 10 * th), 10 * (r - 1), x[2]]) elif nprob == 6: # Powell singular function. - fvec = jnp.array([x[0] + 10 * x[1], jnp.sqrt(5.0) * (x[2] - x[3]), (x[1] - 2 * x[2]) ** 2, jnp.sqrt(10.0) * (x[0] - x[3]) ** 2]) + fvec = np.array([x[0] + 10 * x[1], np.sqrt(5.0) * (x[2] - x[3]), (x[1] - 2 * x[2]) ** 2, np.sqrt(10.0) * (x[0] - x[3]) ** 2]) elif nprob == 7: # Freudenstein and Roth function. - fvec = jnp.array([-c13 + x[0] + ((5 - x[1]) * x[1] - 2) * x[1], -c29 + x[0] + ((1 + x[1]) * x[1] - c14) * x[1]]) + fvec = np.array([-c13 + x[0] + ((5 - x[1]) * x[1] - 2) * x[1], -c29 + x[0] + ((1 + x[1]) * x[1] - c14) * x[1]]) elif nprob == 8: # Bard function. - i = jnp.arange(15) + i = np.arange(15) tmp1 = i + 1 tmp2 = 15 - i - tmp3 = jnp.where(i > 7, tmp2, tmp1) + tmp3 = np.where(i > 7, tmp2, tmp1) denom = x[1] * tmp2 + x[2] * tmp3 fvec = y1 - (x[0] + tmp1 / denom) @@ -218,90 +218,90 @@ def theta(x0, x1): fvec = y2 - x[0] * tmp1 / tmp2 elif nprob == 10: # Meyer function. - i = jnp.arange(16) + i = np.arange(16) temp = 5 * (i + 1) + c45 + x[2] tmp1 = x[1] / temp - tmp2 = jnp.exp(tmp1) + tmp2 = np.exp(tmp1) fvec = x[0] * tmp2 - y3 elif nprob == 11: # Watson function. def watson_term(i): div = (i + 1) / c29 - dx = div ** jnp.arange(n) - dx1 = jnp.arange(1, n) * div ** jnp.arange(1, n) - s1 = jnp.sum(dx1 * x[1:]) - s2 = jnp.sum(dx * x) + dx = div ** np.arange(n) + dx1 = np.arange(1, n) * div ** np.arange(1, n) + s1 = np.sum(dx1 * x[1:]) + s2 = np.sum(dx * x) return s1 - s2**2 - 1 - fvec = jnp.array([watson_term(i) for i in range(29)] + [x[0], x[1] - x[0] ** 2 - 1]) + fvec = np.array([watson_term(i) for i in range(29)] + [x[0], x[1] - x[0] ** 2 - 1]) elif nprob == 12: # Box 3-dimensional function. - i = jnp.arange(1, m + 1) + i = np.arange(1, m + 1) tmp1 = i / 10.0 - term = jnp.exp(-tmp1[:, None] * x[:2]) - const = jnp.exp(-i) - jnp.exp(-tmp1) + term = np.exp(-tmp1[:, None] * x[:2]) + const = np.exp(-i) - np.exp(-tmp1) fvec = term[:, 0] - term[:, 1] + const * x[2] elif nprob == 13: # Jennrich and Sampson function. - i = jnp.arange(1, m + 1) - fvec = 2 + 2 * i - jnp.exp(i * x[0]) - jnp.exp(i * x[1]) + i = np.arange(1, m + 1) + fvec = 2 + 2 * i - np.exp(i * x[0]) - np.exp(i * x[1]) elif nprob == 14: # Brown and Dennis function. - i = jnp.arange(1, m + 1) + i = np.arange(1, m + 1) temp = i / 5.0 - tmp1 = x[0] + temp * x[1] - jnp.exp(temp) - tmp2 = x[2] + jnp.sin(temp) * x[3] - jnp.cos(temp) + tmp1 = x[0] + temp * x[1] - np.exp(temp) + tmp2 = x[2] + np.sin(temp) * x[3] - np.cos(temp) fvec = tmp1**2 + tmp2**2 elif nprob == 15: # Chebyquad function. t = 2 * x - 1 - T = jnp.polynomial.chebyshev.chebvander(t, m - 1).T - coeffs = jnp.mean(T, axis=1) - correction = jnp.array([0.0 if i % 2 == 0 else 1 / ((i + 1) ** 2 - 1) for i in range(m)]) + T = np.polynomial.chebyshev.chebvander(t, m - 1).T + coeffs = np.mean(T, axis=1) + correction = np.array([0.0 if i % 2 == 0 else 1 / ((i + 1) ** 2 - 1) for i in range(m)]) fvec = coeffs + correction elif nprob == 16: # Brown almost-linear function. - total = jnp.sum(x) - (n + 1) - prod = jnp.prod(x) - fvec = jnp.concatenate([x[:-1] + total, jnp.array([prod - 1])]) + total = np.sum(x) - (n + 1) + prod = np.prod(x) + fvec = np.concatenate([x[:-1] + total, np.array([prod - 1])]) elif nprob == 17: # Osborne 1 function. - i = jnp.arange(33) - tmp1 = jnp.exp(-x[3] * 10 * i) - tmp2 = jnp.exp(-x[4] * 10 * i) + i = np.arange(33) + tmp1 = np.exp(-x[3] * 10 * i) + tmp2 = np.exp(-x[4] * 10 * i) fvec = y4 - (x[0] + x[1] * tmp1 + x[2] * tmp2) elif nprob == 18: # Osborne 2 function. - i = jnp.arange(65) + i = np.arange(65) t = i / 10.0 - tmp1 = jnp.exp(-x[4] * t) - tmp2 = jnp.exp(-x[5] * (t - x[8]) ** 2) - tmp3 = jnp.exp(-x[6] * (t - x[9]) ** 2) - tmp4 = jnp.exp(-x[7] * (t - x[10]) ** 2) + tmp1 = np.exp(-x[4] * t) + tmp2 = np.exp(-x[5] * (t - x[8]) ** 2) + tmp3 = np.exp(-x[6] * (t - x[9]) ** 2) + tmp4 = np.exp(-x[7] * (t - x[10]) ** 2) fvec = y5 - (x[0] * tmp1 + x[1] * tmp2 + x[2] * tmp3 + x[3] * tmp4) elif nprob == 19: # Bdqrtic f1 = -4 * x[: n - 4] + 3 f2 = sum((i + 1) * x[i + j] ** 2 for j, i in enumerate(range(n - 4))) - fvec = jnp.concatenate([f1, f2.reshape(-1)]) + fvec = np.concatenate([f1, f2.reshape(-1)]) elif nprob == 20: # Cube - fvec = jnp.zeros(n) + fvec = np.zeros(n) fvec = fvec.at[0].set(x[0] - 1.0) fvec = fvec.at[1:].set(10 * (x[1:] - x[:-1] ** 3)) elif nprob == 21: # Mancino def mancino_term(i): - j = jnp.arange(n) - v2 = jnp.sqrt(x[i] ** 2 + (i + 1) / (j + 1)) - return 1400 * x[i] + (i - 49) ** 3 + jnp.sum(v2 * (jnp.sin(jnp.log(v2)) ** 5 + jnp.cos(jnp.log(v2)) ** 5)) + j = np.arange(n) + v2 = np.sqrt(x[i] ** 2 + (i + 1) / (j + 1)) + return 1400 * x[i] + (i - 49) ** 3 + np.sum(v2 * (np.sin(np.log(v2)) ** 5 + np.cos(np.log(v2)) ** 5)) - fvec = jnp.array([mancino_term(i) for i in range(n)]) + fvec = np.array([mancino_term(i) for i in range(n)]) elif nprob == 22: # Heart8ls - fvec = jnp.zeros(8) + fvec = np.zeros(8) fvec = fvec.at[0].set(x[0] + x[1] + 0.69) fvec = fvec.at[1].set(x[2] + x[3] + 0.044) fvec = fvec.at[2].set(x[4] * x[0] + x[5] * x[1] - x[6] * x[2] - x[7] * x[3] + 1.57) From ba6788f0613d03ca08f401835f3c3758f79ee5e4 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Wed, 1 Oct 2025 11:22:43 -0500 Subject: [PATCH 07/20] Fixing nprob=11 case --- py/dfovec_jax.py | 26 +++++++++++++++++--------- 1 file changed, 17 insertions(+), 9 deletions(-) diff --git a/py/dfovec_jax.py b/py/dfovec_jax.py index 11e1272..81007d0 100644 --- a/py/dfovec_jax.py +++ b/py/dfovec_jax.py @@ -1,3 +1,4 @@ +import jax import jax.numpy as np @@ -225,16 +226,23 @@ def theta(x0, x1): fvec = x[0] * tmp2 - y3 elif nprob == 11: # Watson function. + fvec_body = [] + for i in range(29): + div = (i + 1.0) / c29 + s1 = 0.0 + dx = 1.0 + for j in range(1, n): + s1 += j * dx * x[j] + dx *= div + s2 = 0.0 + dx = 1.0 + for j in range(n): + s2 += dx * x[j] + dx *= div + fvec_body.append(s1 - s2 * s2 - 1.0) + fvec_tail = [x[0], x[1] - x[0]**2 - 1.0] + fvec = np.array(fvec_body + fvec_tail) - def watson_term(i): - div = (i + 1) / c29 - dx = div ** np.arange(n) - dx1 = np.arange(1, n) * div ** np.arange(1, n) - s1 = np.sum(dx1 * x[1:]) - s2 = np.sum(dx * x) - return s1 - s2**2 - 1 - - fvec = np.array([watson_term(i) for i in range(29)] + [x[0], x[1] - x[0] ** 2 - 1]) elif nprob == 12: # Box 3-dimensional function. i = np.arange(1, m + 1) From 89408cc6c2a35e73c622f67d9f886526179bb7f3 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Wed, 1 Oct 2025 11:40:49 -0500 Subject: [PATCH 08/20] Fixing two more cases --- py/dfovec_jax.py | 29 +++++++++++++++++++++-------- py/test_jacobians_with_jax.py | 2 +- 2 files changed, 22 insertions(+), 9 deletions(-) diff --git a/py/dfovec_jax.py b/py/dfovec_jax.py index 81007d0..e43afb4 100644 --- a/py/dfovec_jax.py +++ b/py/dfovec_jax.py @@ -263,11 +263,23 @@ def theta(x0, x1): fvec = tmp1**2 + tmp2**2 elif nprob == 15: # Chebyquad function. - t = 2 * x - 1 - T = np.polynomial.chebyshev.chebvander(t, m - 1).T - coeffs = np.mean(T, axis=1) - correction = np.array([0.0 if i % 2 == 0 else 1 / ((i + 1) ** 2 - 1) for i in range(m)]) - fvec = coeffs + correction + for j in range(n): + t1 = 1.0 + t2 = 2.0 * x[j] - 1.0 + t = 2.0 * t2 + for i in range(m): + fvec = fvec.at[i].add(t2) + th = t * t2 - t1 + t1 = t2 + t2 = th + + iev = -1 + for i in range(m): + val = fvec[i] / n + if iev > 0: + val += 1.0 / ((i + 1)**2 - 1.0) + fvec = fvec.at[i].set(val) + iev = -iev elif nprob == 16: # Brown almost-linear function. total = np.sum(x) - (n + 1) @@ -290,9 +302,10 @@ def theta(x0, x1): fvec = y5 - (x[0] * tmp1 + x[1] * tmp2 + x[2] * tmp3 + x[3] * tmp4) elif nprob == 19: # Bdqrtic - f1 = -4 * x[: n - 4] + 3 - f2 = sum((i + 1) * x[i + j] ** 2 for j, i in enumerate(range(n - 4))) - fvec = np.concatenate([f1, f2.reshape(-1)]) + for i in range(n - 4): + fvec = fvec.at[i].set(-4.0 * x[i] + 3.0) + quad = x[i] ** 2 + 2.0 * x[i + 1] ** 2 + 3.0 * x[i + 2] ** 2 + 4.0 * x[i + 3] ** 2 + 5.0 * x[n - 1] ** 2 + fvec = fvec.at[n - 4 + i].set(quad) elif nprob == 20: # Cube fvec = np.zeros(n) diff --git a/py/test_jacobians_with_jax.py b/py/test_jacobians_with_jax.py index b5b2c48..4f4a900 100644 --- a/py/test_jacobians_with_jax.py +++ b/py/test_jacobians_with_jax.py @@ -29,7 +29,7 @@ def f(x): rel_diff = diff / (np.linalg.norm(J_fd) + 1e-12) print(f"nprob={nprob:2d}, m={m:3d}, n={n:2d} | abs diff: {diff:.2e}, rel diff: {rel_diff:.2e}") - if diff > 1e-7 and rel_diff > 1e-7: + if diff > 2e-6 and rel_diff > 2e-6: raise Exception("Difference is too large") From 6ad88c3c279d03aea0d241e936bbd4fd468164e4 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Wed, 1 Oct 2025 11:41:03 -0500 Subject: [PATCH 09/20] black --- py/dfovec_jax.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/py/dfovec_jax.py b/py/dfovec_jax.py index e43afb4..5b4e3c5 100644 --- a/py/dfovec_jax.py +++ b/py/dfovec_jax.py @@ -240,10 +240,9 @@ def theta(x0, x1): s2 += dx * x[j] dx *= div fvec_body.append(s1 - s2 * s2 - 1.0) - fvec_tail = [x[0], x[1] - x[0]**2 - 1.0] + fvec_tail = [x[0], x[1] - x[0] ** 2 - 1.0] fvec = np.array(fvec_body + fvec_tail) - elif nprob == 12: # Box 3-dimensional function. i = np.arange(1, m + 1) tmp1 = i / 10.0 @@ -277,7 +276,7 @@ def theta(x0, x1): for i in range(m): val = fvec[i] / n if iev > 0: - val += 1.0 / ((i + 1)**2 - 1.0) + val += 1.0 / ((i + 1) ** 2 - 1.0) fvec = fvec.at[i].set(val) iev = -iev From 1b2206607e884a68aa2f343c466a7a109ec46070 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Wed, 1 Oct 2025 14:12:53 -0500 Subject: [PATCH 10/20] Testing all outputs with/without jax --- py/test_all_outputs_with_jax.py | 78 +++++++++++++++++++++++++++++++++ py/test_jacobians_with_jax.py | 2 +- 2 files changed, 79 insertions(+), 1 deletion(-) create mode 100644 py/test_all_outputs_with_jax.py diff --git a/py/test_all_outputs_with_jax.py b/py/test_all_outputs_with_jax.py new file mode 100644 index 0000000..13b0a71 --- /dev/null +++ b/py/test_all_outputs_with_jax.py @@ -0,0 +1,78 @@ +import numpy as np +import jax +import jax.numpy as jnp + +from calfun import calfun +from dfovec_jax import dfovec_jax +from dfoxs import dfoxs +from pathlib import Path + +# Load problem definitions: nprob, n, m, scale_power +dfo_table = np.loadtxt("../data/dfo.dat") # Adjust path as needed + + +def compare_outputs(nprob, m, n, x0): + x0 = np.array(x0) + + # --- Reference values from calfun --- + fval_ref, fvec_ref, grad_ref, J_ref = calfun(x0, m, nprob, probtype="smooth", num_outs=4) + J_ref = J_ref.T # Transpose back: calfun gives (n x m), we want (m x n) + + # --- JAX values --- + x_jax = jnp.array(x0) + + def fvec_fun(x): + return dfovec_jax(m, n, x, nprob) + + def scalar_fun(x): + return jnp.sum(fvec_fun(x) ** 2) + + fvec_jax = fvec_fun(x_jax) + fval_jax = scalar_fun(x_jax) + grad_jax = jax.grad(scalar_fun)(x_jax) + J_jax = jax.jacfwd(fvec_fun)(x_jax) + + # --- Differences --- + fval_diff = np.abs(fval_ref - float(fval_jax)) + fvec_diff = np.linalg.norm(fvec_ref - np.array(fvec_jax)) + grad_diff = np.linalg.norm(grad_ref - np.array(grad_jax)) + jac_diff = np.linalg.norm(J_ref - np.array(J_jax)) + + fval_rel = fval_diff / (np.abs(fval_ref) + 1e-14) + fvec_rel = fvec_diff / (np.linalg.norm(fvec_ref) + 1e-14) + grad_rel = grad_diff / (np.linalg.norm(grad_ref) + 1e-14) + jac_rel = jac_diff / (np.linalg.norm(J_ref) + 1e-14) + + # --- Print Summary --- + print( + f"nprob={nprob:2d}, m={m:3d}, n={n:2d} | " + f"fval Δ={fval_diff:.2e} ({fval_rel:.2e}), " + f"fvec Δ={fvec_diff:.2e} ({fvec_rel:.2e}), " + f"grad Δ={grad_diff:.2e} ({grad_rel:.2e}), " + f"J Δ={jac_diff:.2e} ({jac_rel:.2e})" + ) + + # --- Fail if any relative diff is large --- + if any(val > 6e-7 for val in [fval_rel, fvec_rel, grad_rel, jac_rel]): + raise Exception("Large differences detected") + + +# Loop over all benchmark problems and 3 starting points +for nprob, n, m, factor_power in dfo_table: + n = int(n) + m = int(m) + nprob = int(nprob) + scale = int(10**factor_power) + + for pt in range(3): + if pt == 0: + x0 = dfoxs(n, nprob, scale) + elif pt == 1: + x0 = 0.1 * np.ones(n) + elif pt == 2: + x0 = 0.1 * np.arange(1, n + 1) + + try: + compare_outputs(nprob, m, n, x0) + except Exception as e: + print(f"nprob={nprob:2d}, pt={pt} failed: {e}") diff --git a/py/test_jacobians_with_jax.py b/py/test_jacobians_with_jax.py index 4f4a900..c88d603 100644 --- a/py/test_jacobians_with_jax.py +++ b/py/test_jacobians_with_jax.py @@ -29,7 +29,7 @@ def f(x): rel_diff = diff / (np.linalg.norm(J_fd) + 1e-12) print(f"nprob={nprob:2d}, m={m:3d}, n={n:2d} | abs diff: {diff:.2e}, rel diff: {rel_diff:.2e}") - if diff > 2e-6 and rel_diff > 2e-6: + if diff > 6e-7 and rel_diff > 6e-7: raise Exception("Difference is too large") From 523646c33542181a41b81cae16d0e498aa7ba7e9 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Thu, 2 Oct 2025 14:48:38 -0500 Subject: [PATCH 11/20] Updating nprob=5 case --- py/dfovec_jax.py | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/py/dfovec_jax.py b/py/dfovec_jax.py index 5b4e3c5..5c07798 100644 --- a/py/dfovec_jax.py +++ b/py/dfovec_jax.py @@ -191,13 +191,12 @@ def dfovec_jax(m, n, x, nprob): elif nprob == 5: # Helical valley function. - def theta(x0, x1): - angle = np.arctan2(x1, x0) / (2 * np.pi) - return np.where(x0 > 0, angle, np.where(x0 < 0, angle + 0.5, np.where((x0 == 0) & (x1 == 0), 0.0, 0.25))) - - th = theta(x[0], x[1]) + th = np.arctan2(x[1], x[0]) / (2.0 * np.pi) r = np.sqrt(x[0] ** 2 + x[1] ** 2) - fvec = np.array([10 * (x[2] - 10 * th), 10 * (r - 1), x[2]]) + + fvec = fvec.at[0].set(10.0 * (x[2] - 10.0 * th)) + fvec = fvec.at[1].set(10.0 * (r - 1.0)) + fvec = fvec.at[2].set(x[2]) elif nprob == 6: # Powell singular function. fvec = np.array([x[0] + 10 * x[1], np.sqrt(5.0) * (x[2] - x[3]), (x[1] - 2 * x[2]) ** 2, np.sqrt(10.0) * (x[0] - x[3]) ** 2]) From 1ebfb314e26cfd87d499b39e725be77c6882d74e Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Tue, 7 Oct 2025 12:23:01 -0500 Subject: [PATCH 12/20] Trying to debug black --- .github/workflows/black.yml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/.github/workflows/black.yml b/.github/workflows/black.yml index 92d338f..11af541 100644 --- a/.github/workflows/black.yml +++ b/.github/workflows/black.yml @@ -13,7 +13,7 @@ jobs: uses: rickstaa/action-black@v1 id: action_black with: - black_args: "--check --config=.black ." + black_args: "--check --diff --verbose --color --config=.black ." - name: Annotate diff changes using reviewdog if: steps.action_black.outputs.is_formatted == 'true' uses: reviewdog/action-suggester@v1 From 14b9e67e8da021286eee9b170e21f4497a35aa32 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Tue, 7 Oct 2025 12:24:33 -0500 Subject: [PATCH 13/20] Whitespace --- py/dfovec_jax.py | 1 - 1 file changed, 1 deletion(-) diff --git a/py/dfovec_jax.py b/py/dfovec_jax.py index 5c07798..65ab9c3 100644 --- a/py/dfovec_jax.py +++ b/py/dfovec_jax.py @@ -190,7 +190,6 @@ def dfovec_jax(m, n, x, nprob): fvec = np.array([10 * (x[1] - x[0] ** 2), 1 - x[0]]) elif nprob == 5: # Helical valley function. - th = np.arctan2(x[1], x[0]) / (2.0 * np.pi) r = np.sqrt(x[0] ** 2 + x[1] ** 2) From 7a8c377851fa1634b2a2aba96ffdbf7fa9f16264 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Tue, 7 Oct 2025 12:29:18 -0500 Subject: [PATCH 14/20] Moving test to testing --- py/test_jacobians_with_jax.py | 54 ------------------- .../compare_numpy_and_jax.py | 5 +- 2 files changed, 4 insertions(+), 55 deletions(-) delete mode 100644 py/test_jacobians_with_jax.py rename py/test_all_outputs_with_jax.py => testing/compare_numpy_and_jax.py (98%) diff --git a/py/test_jacobians_with_jax.py b/py/test_jacobians_with_jax.py deleted file mode 100644 index c88d603..0000000 --- a/py/test_jacobians_with_jax.py +++ /dev/null @@ -1,54 +0,0 @@ -import numpy as np -import jax -import jax.numpy as jnp - -from dfovec import dfovec -from dfovec_jax import dfovec_jax -from jacobian import jacobian -from dfoxs import dfoxs # for generating X0s - -# Load benchmark problem definitions -dfo_table = np.loadtxt("../data/dfo.dat") - - -def compare_jacobians(nprob, m, n, x0): - x0 = np.array(x0) - fvec_np = dfovec(m, n, x0, nprob) - - # Get FD-based Jacobian - J_fd, _ = jacobian(m, n, x0, nprob) # shape (m, n) - - # Get JAX-based Jacobian - def f(x): - return dfovec_jax(m, n, x, nprob) - - J_jax = jax.jacfwd(f)(jnp.array(x0)) # shape (m, n) - - # Compare - diff = np.linalg.norm(J_fd - np.array(J_jax)) - rel_diff = diff / (np.linalg.norm(J_fd) + 1e-12) - - print(f"nprob={nprob:2d}, m={m:3d}, n={n:2d} | abs diff: {diff:.2e}, rel diff: {rel_diff:.2e}") - if diff > 6e-7 and rel_diff > 6e-7: - raise Exception("Difference is too large") - - -# Loop over all benchmark problems and 3 starting points -for nprob, n, m, factor_power in dfo_table: - n = int(n) - m = int(m) - nprob = int(nprob) - scale = int(10**factor_power) - - for pt in range(3): - if pt == 0: - x0 = dfoxs(n, nprob, scale) - elif pt == 1: - x0 = 0.1 * np.ones(n) - elif pt == 2: - x0 = 0.1 * np.arange(1, n + 1) - - try: - compare_jacobians(nprob, m, n, x0) - except Exception as e: - print(f"nprob={nprob:2d} failed at pt={pt}: {e}") diff --git a/py/test_all_outputs_with_jax.py b/testing/compare_numpy_and_jax.py similarity index 98% rename from py/test_all_outputs_with_jax.py rename to testing/compare_numpy_and_jax.py index 13b0a71..5bcee17 100644 --- a/py/test_all_outputs_with_jax.py +++ b/testing/compare_numpy_and_jax.py @@ -1,11 +1,14 @@ import numpy as np +import sys import jax import jax.numpy as jnp +from pathlib import Path + +sys.path.append('../py') from calfun import calfun from dfovec_jax import dfovec_jax from dfoxs import dfoxs -from pathlib import Path # Load problem definitions: nprob, n, m, scale_power dfo_table = np.loadtxt("../data/dfo.dat") # Adjust path as needed From 62c19e46f99007fc53487a3fa4a1dec12b2900f4 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Tue, 7 Oct 2025 12:29:33 -0500 Subject: [PATCH 15/20] isort --- testing/compare_numpy_and_jax.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/testing/compare_numpy_and_jax.py b/testing/compare_numpy_and_jax.py index 5bcee17..db52ad3 100644 --- a/testing/compare_numpy_and_jax.py +++ b/testing/compare_numpy_and_jax.py @@ -1,9 +1,9 @@ -import numpy as np import sys +from pathlib import Path + import jax import jax.numpy as jnp - -from pathlib import Path +import numpy as np sys.path.append('../py') from calfun import calfun From 54f31948327f7eb1b6e557cb02c35e9da0eedd22 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Tue, 7 Oct 2025 12:29:54 -0500 Subject: [PATCH 16/20] isort --- py/calfun.py | 1 + py/evaluate_all_calfuns_points.py | 1 + py/jacobian.py | 1 + 3 files changed, 3 insertions(+) diff --git a/py/calfun.py b/py/calfun.py index c5f1afe..542b723 100644 --- a/py/calfun.py +++ b/py/calfun.py @@ -1,4 +1,5 @@ import numpy as np + from dfovec import dfovec from jacobian import jacobian diff --git a/py/evaluate_all_calfuns_points.py b/py/evaluate_all_calfuns_points.py index 39117a4..df55aed 100644 --- a/py/evaluate_all_calfuns_points.py +++ b/py/evaluate_all_calfuns_points.py @@ -1,5 +1,6 @@ import numpy as np import scipy as sp + from calfun import calfun from dfoxs import dfoxs diff --git a/py/jacobian.py b/py/jacobian.py index e3727df..f3c5236 100644 --- a/py/jacobian.py +++ b/py/jacobian.py @@ -1,4 +1,5 @@ import numpy as np + from g_dfovec_1d import g_dfovec_1d From ede54dc496bd2a585ebbe17d939fe465d02d6979 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Tue, 7 Oct 2025 13:17:21 -0500 Subject: [PATCH 17/20] black --- testing/compare_numpy_and_jax.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testing/compare_numpy_and_jax.py b/testing/compare_numpy_and_jax.py index db52ad3..e01abe4 100644 --- a/testing/compare_numpy_and_jax.py +++ b/testing/compare_numpy_and_jax.py @@ -5,7 +5,7 @@ import jax.numpy as jnp import numpy as np -sys.path.append('../py') +sys.path.append("../py") from calfun import calfun from dfovec_jax import dfovec_jax from dfoxs import dfoxs From 2eb7d4183abe07f295ae6fba9489c9177f0ffb51 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Wed, 29 Oct 2025 13:17:02 -0500 Subject: [PATCH 18/20] Jax is single precision by default --- testing/compare_numpy_and_jax.py | 1 + 1 file changed, 1 insertion(+) diff --git a/testing/compare_numpy_and_jax.py b/testing/compare_numpy_and_jax.py index e01abe4..9c945bb 100644 --- a/testing/compare_numpy_and_jax.py +++ b/testing/compare_numpy_and_jax.py @@ -4,6 +4,7 @@ import jax import jax.numpy as jnp import numpy as np +jax.config.update("jax_enable_x64", True) sys.path.append("../py") from calfun import calfun From 06f4d413366677550b9fdfca24b778e7b9e46724 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Wed, 29 Oct 2025 13:21:00 -0500 Subject: [PATCH 19/20] black --- testing/compare_numpy_and_jax.py | 1 + 1 file changed, 1 insertion(+) diff --git a/testing/compare_numpy_and_jax.py b/testing/compare_numpy_and_jax.py index 9c945bb..ac3c233 100644 --- a/testing/compare_numpy_and_jax.py +++ b/testing/compare_numpy_and_jax.py @@ -4,6 +4,7 @@ import jax import jax.numpy as jnp import numpy as np + jax.config.update("jax_enable_x64", True) sys.path.append("../py") From 152a7d62e62006ac331f58868b5b1b42ea3206e4 Mon Sep 17 00:00:00 2001 From: Jeffrey Larson Date: Mon, 3 Nov 2025 12:43:34 -0600 Subject: [PATCH 20/20] Noting that the errors are much much smaller now --- testing/compare_numpy_and_jax.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/testing/compare_numpy_and_jax.py b/testing/compare_numpy_and_jax.py index ac3c233..0145d6a 100644 --- a/testing/compare_numpy_and_jax.py +++ b/testing/compare_numpy_and_jax.py @@ -58,7 +58,7 @@ def scalar_fun(x): ) # --- Fail if any relative diff is large --- - if any(val > 6e-7 for val in [fval_rel, fvec_rel, grad_rel, jac_rel]): + if any(val > 5e-15 for val in [fval_rel, fvec_rel, grad_rel, jac_rel]): raise Exception("Large differences detected")