lib/minpack/rosenbrock_figures.viv

# The two figures the Read view of lmder's run on Rosenbrock shows beside
# its core path (/explore/minpack-rosenbrock):
#
#   rf_path_figure  the iterates on Rosenbrock's level curves, with the
#                   trust region each accepted step was taken in and the
#                   rejected first trial;
#   rf_phi_figure   the secular function phi(lambda) that lmpar solves on the
#                   first pass that needs lambda > 0, with its bracket and
#                   the values of lambda it tries.
#
# Everything drawn is computed here, from lmder's hist and from Rosenbrock's
# residuals and Jacobian (mgh_fcn); nothing is typed in. The tests at the end
# check the numbers against the gfortran oracle's record of the same run
# (lm_rosenbrock_oracle in lmder.viv).
#
# Used with `import minpack.rosenbrock_figures`.

from minpack.lmder import (
    H_ACCEPT, H_DELTA_IN, H_FNORM1, H_PAR_LM, H_PNORM, H_X, lm_rosenbrock_oracle, lmder, lmder1,
    LMDER1_FACTOR, lmder1_maxfev,
)
from minpack.lmpar import LMPAR_MAXITER, LMPAR_TOL
from minpack.enorm import enorm, EPSMCH
from minpack.mgh import mgh_data, mgh_fcn, mgh_start

# Rosenbrock's function in MINPACK's test set: problem 4, 2 residuals in 2
# unknowns
RF_NPROB = 4.0
RF_M = 2.0
RF_N = 2.0
# Rosenbrock's first residual is RF_WEIGHT (x2 - x1^2)
RF_WEIGHT = 10.0
# a full turn, in radians
RF_TURN = 2.0 * 3.141592653589793
# points on each closed curve drawn (a level curve or a trust region)
RF_CURVE_PTS = 181
# the sums of squares the level curves are drawn at; the start's is 24.2
RF_LEVELS = [0.25, 1.0, 4.0, 24.2, 100.0]
# the window of the path figure
RF_X1_LO = -2.2
RF_X1_HI = 2.2
RF_X2_LO = -4.2
RF_X2_HI = 2.2
# the values of lambda the secular function is drawn at: log-spaced from
# RF_LAM_LO to the upper bound of lmpar's bracket
RF_PHI_PTS = 161
RF_LAM_LO = 1.0e-4
# the phi axis runs up to this many trust radii, so the band where lmpar
# stops (a tenth of a radius either side of zero) can be seen; phi at
# lambda = 0 is far above it
RF_PHI_TOP = 2.0
# the label of the upper bound sits this fraction of it to the left, inside
# the window (the axis is logarithmic)
RF_LABEL_IN = 0.3

# Rosenbrock's residuals at x, into f, and its Jacobian, into J.
def rf_eval(x: Arr, f: Arr, J: Arr) {
    data = mgh_data(RF_NPROB, RF_M, RF_N)
    mgh_fcn(data, x, f, J, 1.0)
    mgh_fcn(data, x, f, J, 2.0)