/// @module std::math::optimize
/// Numerical Optimization — Nelder-Mead (downhill simplex) method
///
/// Minimizes a scalar objective function over n-dimensional space
/// with optional box constraints.
/// Options for the optimizer.
pub type OptimizeOptions {
tol: number,
max_iter: int,
bounds: Array<Array<number>>?
}
/// Result of an optimization run.
pub type OptimizeResult {
x: Array<number>,
fun: number,
converged: bool,
iterations: int
}
/// Compute the centroid of all simplex vertices except one.
fn simplex_centroid(vertices: Array<Array<number>>, exclude_idx: int) {
let n: int = vertices[0].length
let count: number = (vertices.length - 1) as number
let mut sums: Array<number> = []
let mut j = 0
while j < n {
sums = sums.push(0.0)
j = j + 1
}
let mut i = 0
while i < vertices.length {
if i != exclude_idx {
let mut j = 0
while j < n {
sums = vec_set_n(sums, j, sums[j] + vertices[i][j])
j = j + 1
}
}
i = i + 1
}
sums.map(|s| s / count)
}
/// Set element at index in an Array<number> (returns new array).
fn vec_set_n(arr: Array<number>, idx: int, val: number) -> Array<number> {
let mut result: Array<number> = [val]
let mut i = 0
while i < arr.length {
let entry: number = if i == idx { val } else { arr[i] }
if i == 0 { result = [entry] } else { result = result.push(entry) }
i = i + 1
}
result
}
/// Set element at index in an Array<Array<number>> (returns new array).
fn vec_set_an(arr: Array<Array<number>>, idx: int, val: Array<number>) -> Array<Array<number>> {
let mut result = [val]
let mut i = 0
while i < arr.length {
let entry: Array<number> = if i == idx { val } else { arr[i] }
if i == 0 { result = [entry] } else { result = result.push(entry) }
i = i + 1
}
result
}
/// Reflect the worst point through the centroid.
fn simplex_reflect(centroid: Array<number>, worst: Array<number>, alpha: number) {
let n = centroid.length
let mut result: Array<number> = []
let mut i = 0
while i < n {
result = result.push(centroid[i] + alpha * (centroid[i] - worst[i]))
i = i + 1
}
result
}
/// Expand beyond the reflected point.
fn simplex_expand(centroid: Array<number>, reflected: Array<number>, gamma: number) {
let n = centroid.length
let mut result: Array<number> = []
let mut i = 0
while i < n {
result = result.push(centroid[i] + gamma * (reflected[i] - centroid[i]))
i = i + 1
}
result
}
/// Contract toward the centroid from the worst point.
fn simplex_contract(centroid: Array<number>, worst: Array<number>, rho: number) {
let n = centroid.length
let mut result: Array<number> = []
let mut i = 0
while i < n {
result = result.push(centroid[i] + rho * (worst[i] - centroid[i]))
i = i + 1
}
result
}
/// Shrink all vertices toward the best vertex.
fn simplex_shrink(vertices: Array<Array<number>>, best_idx: int, sigma: number) -> Array<Array<number>> {
let best: Array<number> = vertices[best_idx]
let n: int = best.length
// Seed with the best vertex; then overwrite element 0 on first iter
let mut result: Array<Array<number>> = [best]
let mut i: int = 0
while i < vertices.length {
let entry: Array<number> = if i == best_idx {
best
} else {
let mut v: Array<number> = []
let mut j: int = 0
while j < n {
v = v.push(best[j] + sigma * (vertices[i][j] - best[j]))
j = j + 1
}
v
}
if i == 0 {
result = [entry]
} else {
result = result.push(entry)
}
i = i + 1
}
result
}
/// Apply box constraints to a point.
fn apply_bounds(x: Array<number>, bounds) {
if bounds == None {
x
} else {
let mut result: Array<number> = []
let mut i = 0
while i < x.length {
if i < bounds.length {
let lo = bounds[i][0]
let hi = bounds[i][1]
let v = if x[i] < lo { lo } else if x[i] > hi { hi } else { x[i] }
result = result.push(v)
} else {
result = result.push(x[i])
}
i = i + 1
}
result
}
}
/// Compute the standard deviation of an array of values.
fn vec_std(values: Array<number>) -> number {
let n: int = values.length
if n == 0 { 0.0 }
else {
let nf: number = n as number
let mut sum: number = 0.0
let mut i = 0
while i < n {
sum = sum + values[i]
i = i + 1
}
let mean: number = sum / nf
let mut sq_sum: number = 0.0
let mut i = 0
while i < n {
let d: number = values[i] - mean
sq_sum = sq_sum + d * d
i = i + 1
}
sqrt(sq_sum / nf)
}
}
/// Sort simplex vertices and their function values by ascending value.
/// Returns [sorted_vertices, sorted_values] as a 2-element array.
fn sort_simplex(vertices: Array<Array<number>>, values: Array<number>) {
// Simple insertion sort — simplex is small (n+1 elements)
let n = vertices.length
let mut verts = vertices
let mut vals: Array<number> = values
let mut i = 1
while i < n {
let mut j = i
while j > 0 {
if vals[j] < vals[j - 1] {
// Swap values
let tmp_val = vals[j]
vals = vec_set_n(vals, j, vals[j - 1])
vals = vec_set_n(vals, j - 1, tmp_val)
// Swap vertices
let tmp_vert: Array<number> = verts[j]
verts = vec_set_an(verts, j, verts[j - 1])
verts = vec_set_an(verts, j - 1, tmp_vert)
j = j - 1
} else {
j = 0 // break
}
}
i = i + 1
}
[verts, vals]
}
pub fn minimize(f: (Array<number>) => number, x0: Array<number>, options: OptimizeOptions) {
let tol: number = options.tol
let max_iter: int = options.max_iter
let bounds: Array<Array<number>>? = options.bounds
let n: int = x0.length
let alpha = 1.0
let gamma = 2.0
let rho = 0.5
let sigma = 0.5
let v0 = apply_bounds(x0, bounds)
let f0: number = f(v0)
let mut vertices = [v0]
let mut values = [f0]
// Perturb each dimension to create n more vertices
let mut i: int = 0
while i < n {
let mut v: Array<number> = []
let mut j: int = 0
while j < n {
if j == i {
let delta = if abs(x0[j]) > 0.0001 { x0[j] * 0.05 } else { 0.00025 }
v = v.push(x0[j] + delta)
} else {
v = v.push(x0[j])
}
j = j + 1
}
let vb: Array<number> = apply_bounds(v, bounds)
vertices = vertices.push(vb)
values = values.push(f(vb))
i = i + 1
}
let mut sorted = sort_simplex(vertices, values)
vertices = sorted[0]
values = sorted[1]
let mut iter: int = 0
let mut converged = false
while iter < max_iter {
if vec_std(values) < tol {
converged = true
iter = max_iter
} else {
let worst_idx: int = n
let worst: Array<number> = vertices[worst_idx]
let f_worst: number = values[worst_idx]
let f_best: number = values[0]
let f_second_worst: number = values[n - 1]
let centroid: Array<number> = simplex_centroid(vertices, worst_idx)
let x_r: Array<number> = apply_bounds(simplex_reflect(centroid, worst, alpha), bounds)
let f_r: number = f(x_r)
if f_r < f_second_worst {
if f_r < f_best {
let x_e: Array<number> = apply_bounds(simplex_expand(centroid, x_r, gamma), bounds)
let f_e: number = f(x_e)
if f_e < f_r {
vertices = vec_set_an(vertices, worst_idx, x_e)
values = vec_set_n(values, worst_idx, f_e)
} else {
vertices = vec_set_an(vertices, worst_idx, x_r)
values = vec_set_n(values, worst_idx, f_r)
}
} else {
vertices = vec_set_an(vertices, worst_idx, x_r)
values = vec_set_n(values, worst_idx, f_r)
}
} else {
let x_c: Array<number> = apply_bounds(simplex_contract(centroid, worst, rho), bounds)
let f_c: number = f(x_c)
if f_c < f_worst {
vertices = vec_set_an(vertices, worst_idx, x_c)
values = vec_set_n(values, worst_idx, f_c)
} else {
vertices = simplex_shrink(vertices, 0, sigma)
let mut new_values: Array<number> = [values[0]]
let mut k: int = 1
while k < vertices.length {
let vb: Array<number> = apply_bounds(vertices[k], bounds)
vertices = vec_set_an(vertices, k, vb)
new_values = new_values.push(f(vb))
k = k + 1
}
values = new_values
}
}
sorted = sort_simplex(vertices, values)
vertices = sorted[0]
values = sorted[1]
iter = iter + 1
}
}
OptimizeResult {
x: vertices[0],
fun: values[0],
converged: converged,
iterations: iter
}
}