import * as bfjs from "./bf.js";
/**
* High-precision Function Minimization using Brent's Method (similar to MATLAB's fminbnd).
*
* This algorithm finds a local minimum of a function of one variable within a fixed interval.
* It combines Golden Section Search (linear convergence) with Parabolic Interpolation
* (superlinear convergence) for reliability and speed.
*
* @param {Function} f - The objective function to minimize.
* Must accept a BigFloat argument and return a BigFloat result.
*
* @param {number|string|BigFloat} _ax - The start of the search interval.
* @param {number|string|BigFloat} _bx - The end of the search interval.
*
*
* @param {Object} [info={}] - Configuration and Status object.
* Updates in-place with statistics (iterations, execution time, error estimate).
* // --- Input Configuration Properties ---
* @param {number|string|BigFloat} [info._e=1e-30] - Absolute Error Tolerance.
*
* @param {number|string|BigFloat} [info._re=info._e] - Relative Error Tolerance.
* The convergence criteria is based on the position x, not the function value.
* tol = |x| * _re + _e
* @param {number} [info.max_step=500] - Maximum number of iterations allowed.
* If this limit is reached without convergence, the function returns null and logs a warning.
* @param {number} [info.max_time=60000] - Maximum execution time in milliseconds.
* Prevents the function from hanging in infinite loops or extremely slow computations.
* @param {Function} [info.cb] - Optional callback function.
* If defined, this function is called after every iteration. Useful for updating UI progress or logging.
* @param {boolean} [info.debug] - Optional flag to enable debug logging (implementation dependent).
*
* // --- Output Status Properties (Updated during execution) ---
* @param {BigFloat|null} info.result - The final found root.
* Returns a BigFloat if converged, or null if failed.
* @param {BigFloat} info.lastresult - The result of the last iteration (current best guess `b`).
* Even if convergence fails, this contains the value closest to the root when execution stopped.
* @param {string} info.eff_result - String representation of the result based on effective precision.
* Generated by truncating `lastresult` according to `eff_decimal_precision`.
* @param {number} info.steps - The number of iterations currently executed.
* @param {number} info.exectime - The elapsed execution time in milliseconds.
* @param {BigFloat} info.error - The estimated error bound.
* Typically represents half the width of the current search interval (`xm`).
* @param {BigFloat} info.residual - The absolute value of the function at the current best guess: `|f(x)|`.
* Ideally, this value should be close to zero.
* @param {number} info.eff_decimal_precision - Estimated number of significant decimal digits.
* Calculated as `-log10(error)`.
* @param {Function} info.toString - Helper method.
* Returns a formatted string containing steps, error, residual, and execution time.
*
* @returns {BigFloat|null}
* Returns the BigFloat x where f(x) is minimized.
* Returns `null` if max steps or time limit exceeded.
*/
export function fminbnd(f, _ax, _bx, info = {}) {
let _e = info._e ?? 1e-30;
let _re = info._re ?? _e;
// 1. Configuration & Constants
let max_step = info.max_step || 500, // Minimization often takes more steps than root finding
max_time = info.max_time || 60000;
if (typeof(_e) !== 'number' || typeof(_re) !== 'number' || typeof(info) !== "object") {
throw new Error("arguments error");
}
const start_time = new Date().getTime();
// Cache BigFloat constants
const bf_zero = bfjs.zero;
const bf_one = bfjs.one;
const bf_two = bfjs.two;
const bf_half = bfjs.half;
// Golden Section Ratio: (3 - sqrt(5)) / 2 ≈ 0.381966011...
// We calculate it at runtime to match current precision
const bf_golden = bfjs.three.sub(bfjs.bf(5).sqrt()).mul(bf_half);
// Initial Interval
let a = bfjs.bf(_ax);
let b = bfjs.bf(_bx);
// Ensure a < b
if (a.cmp(b) > 0) {
let temp = a; a = b; b = temp;
}
// Initialization of Brent's variables
// x: is the point with the very least function value found so far
// w: is the point with the second least function value
// v: is the previous value of w
// u: is the point at which the function has been evaluated most recently
// Start x, w, v at: a + golden_ratio * (b - a)
let c = b.sub(a);
let x = a.add(c.mul(bf_golden));
let w = x;
let v = x;
let fx = f(x);
let fw = fx;
let fv = fx;
let d = bf_zero; // Step taken in the iteration before last
let e = bf_zero; // Step taken in the last iteration
// Tolerances
let tol_act = bfjs.bf(_e);
let tol_rel = bfjs.bf(_re);
// Helper to sync local variables to info object
const updateInfo = (iter, errorBound, currentMinVal) => {
info.exectime = new Date().getTime() - start_time;
info.steps = iter;
info.lastresult = x; // Current best position
info.min_value = currentMinVal; // Current best function value
info.error = errorBound; // Interval width or step estimate
// Calculate effective decimal precision of the position x
if (errorBound.cmp(bf_zero) === 0) {
info.eff_decimal_precision = bfjs.decimalPrecision();
} else {
info.eff_decimal_precision = Math.floor(-errorBound.log().f64() / Math.log(10));
}
if (info.eff_decimal_precision <= 0) {
info.eff_decimal_precision = 0;
info.eff_result = '';
} else {
let limit = bfjs.decimalPrecision();
let prec = info.eff_decimal_precision > limit ? limit : info.eff_decimal_precision;
info.eff_result = x.toString(10, prec);
}
};
info.toString = function() {
return `xmin=${this.eff_result},
f(xmin)=${this.min_value ? this.min_value.toString(10, 6) : 'N/A'},
steps=${this.steps}/${max_step},
error=${this.error ? this.error.toString(10, 3) : 'N/A'},
exectime=${this.exectime}/${max_time}`;
};
// 2. Main Loop
for (let iter = 1; iter <= max_step; ++iter) {
// Midpoint of current interval
let xm = a.add(b).mul(bf_half);
// Calculate tolerance based on current x magnitude
// tol1 = |x| * tol_rel + tol_act
let tol1 = x.abs().mul(tol_rel).add(tol_act);
let tol2 = tol1.mul(bf_two);
// Check for Convergence
// Stop if the interval |b-a| is small enough (roughly < 2*tol)
// Strictly: |x - xm| <= (2*tol1 - 0.5*(b-a))
// Simplified robust check: Is max distance from x to a or b <= 2*tol1?
let dist = x.sub(xm).abs().add(b.sub(a).mul(bf_half));
if (dist.cmp(tol2) <= 0) {
updateInfo(iter, dist, fx);
info.result = x;
return x;
}
// Callback and Timeout Check
if (info.cb) {
updateInfo(iter, dist, fx);
info.cb();
}
if (new Date().getTime() - start_time > max_time) {
updateInfo(iter, dist, fx);
console.log("fminbnd: Timeout reached.");
info.result = null;
return null;
}
// --- Determine Step ---
let p = bf_zero;
let q = bf_zero;
let r = bf_zero;
let new_step = bf_zero;
// Is Parabolic Interpolation possible?
// Needs |e| > tol1 (previous step was significant)
if (e.abs().cmp(tol1) > 0) {
// Parabolic fit using x, w, v
// r = (x - w) * (fx - fv)
// q = (x - v) * (fx - fw)
// p = (x - v) * q - (x - w) * r
// q = 2.0 * (q - r)
r = x.sub(w).mul(fx.sub(fv));
q = x.sub(v).mul(fx.sub(fw));
p = x.sub(v).mul(q).sub(x.sub(w).mul(r));
q = q.sub(r).mul(bf_two);
if (q.cmp(bf_zero) > 0) {
p = p.neg();
}
q = q.abs();
let etemp = e;
e = d;
// Check if parabolic fit is acceptable
// 1. p must be within interval (a, b) specific bounds: |p| < |q * 0.5 * (a/b - x)|
// Strictly: (abs(p) < abs(0.5*q*etemp)) AND (p > q*(a-x)) AND (p < q*(b-x))
// 2. Step size must be shrinking
let is_parabolic_valid = true;
// Condition 1: Step smaller than half of previous-previous step? (Convergence check)
if (p.abs().cmp(q.mul(etemp).abs().mul(bf_half)) >= 0) {
is_parabolic_valid = false;
}
// Condition 2: Must land within (a, b)
// p/q calculates the offset from x. x + p/q must be in [a, b]
// Actually, we usually require it to be within [a+tol, b-tol]
else {
let p_bound_a = q.mul(a.sub(x));
let p_bound_b = q.mul(b.sub(x));
// Since q > 0, we can compare directly
if (p.cmp(p_bound_a) <= 0 || p.cmp(p_bound_b) >= 0) {
is_parabolic_valid = false;
}
}
if (is_parabolic_valid) {
// Accept Parabolic Step
d = p.div(q);
new_step = d;
// If step lands very close to bounds, adjust to ensure we don't violate precision
let u_tentative = x.add(d);
if (u_tentative.sub(a).cmp(tol2) < 0 || b.sub(u_tentative).cmp(tol2) < 0) {
// d = sign(xm - x) * tol1
let sign = xm.sub(x).cmp(bf_zero) >= 0 ? bf_one : bf_one.neg();
d = tol1.mul(sign);
}
} else {
// Reject Parabolic, switch to Golden Section
e = (x.cmp(xm) >= 0 ? a : b).sub(x);
d = bf_golden.mul(e);
}
} else {
// Golden Section Step (first iteration or small steps)
e = (x.cmp(xm) >= 0 ? a : b).sub(x);
d = bf_golden.mul(e);
}
// --- Apply Step ---
// Ensure we move at least tol1 to allow convergence check to pass eventually
// u = x + d (if |d| >= tol1) else x + sign(d)*tol1
if (d.abs().cmp(tol1) >= 0) {
new_step = d;
} else {
let sign = d.cmp(bf_zero) >= 0 ? bf_one : bf_one.neg();
new_step = tol1.mul(sign);
}
let u = x.add(new_step);
let fu = f(u);
// --- Update Brackets (a, b) and Points (x, w, v) ---
if (fu.cmp(fx) <= 0) {
// Found a new minimum!
if (u.cmp(x) >= 0) {
a = x;
} else {
b = x;
}
v = w; fv = fw;
w = x; fw = fx;
x = u; fx = fu;
} else {
// Not a new minimum, but use it to tighten interval
if (u.cmp(x) < 0) {
a = u;
} else {
b = u;
}
// Should we update w or v?
if (fu.cmp(fw) <= 0 || w.cmp(x) === 0) {
v = w; fv = fw;
w = u; fw = fu;
} else if (fu.cmp(fv) <= 0 || v.cmp(x) === 0 || v.cmp(w) === 0) {
v = u; fv = fu;
}
}
}
// 3. Exit (Max steps)
let final_dist = x.sub(a.add(b).mul(bf_half)).abs().add(b.sub(a).mul(bf_half));
updateInfo(max_step, final_dist, fx);
console.log(`fminbnd: Failed to converge after ${max_step} steps.`);
info.result = null;
return null;
}