basin 1.10.0

Numerical optimization in pure Rust, with pluggable linear-algebra backends and WASM support.
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75

//! Linear interpolation models for COBYLA.
//!
//! At each iteration COBYLA reads the gradient of the linear model of the
//! objective, `g`, and of the constraints, `A` (n × m, one column of gradients
//! per constraint), from the simplex displacements via `simi = inv(sim[:,0..n])`.
//! These are PRIMA's `g = matprod(fval(1:n)−fval(n+1), simi)` and
//! `A = transpose((conmat(:,1:n)−conmat(:,n+1)) · simi)` (with `m_lcon = 0`,
//! since basin's COBYLA carries no separate linear-constraint block).

use crate::core::math::Scalar;

/// Reusable interpolation-model storage.
pub(crate) struct ModelWork<F> {
    pub(crate) g: Vec<F>,
    pub(crate) a: Vec<F>,
    pub(crate) b: Vec<F>,
    difference: Vec<F>,
}

impl<F: Scalar> ModelWork<F> {
    pub(crate) fn new(n: usize, m: usize) -> Self {
        Self {
            g: vec![F::zero(); n],
            a: vec![F::zero(); n * m],
            b: vec![F::zero(); m],
            difference: vec![F::zero(); n],
        }
    }

    pub(crate) fn build(&mut self, fval: &[F], conmat: &[F], simi: &[F]) {
        let n = self.g.len();
        let m = self.b.len();
        build_g_into(fval, simi, n, &mut self.g);
        build_a_into(conmat, simi, n, m, &mut self.a, &mut self.difference);
        for i in 0..m {
            self.b[i] = -conmat[i + n * m];
        }
    }
}

/// Objective-model gradient `g` (length `n`):
/// `g[l] = Σ_i (fval[i] − fval[n]) · simi[i, l]`.
fn build_g_into<F: Scalar>(fval: &[F], simi: &[F], n: usize, g: &mut [F]) {
    let fn_pole = fval[n];
    for l in 0..n {
        g[l] = (0..n).map(|i| (fval[i] - fn_pole) * simi[i + l * n]).sum();
    }
}

/// Constraint-model gradients `A` (n × m, column-major; column `i` is the
/// gradient of constraint `i`):
/// `A[l, i] = Σ_j (conmat[i, j] − conmat[i, n]) · simi[j, l]`.
fn build_a_into<F: Scalar>(
    conmat: &[F],
    simi: &[F],
    n: usize,
    m: usize,
    a: &mut [F],
    difference: &mut [F],
) {
    for i in 0..m {
        let pole = conmat[i + n * m];
        for j in 0..n {
            difference[j] = conmat[i + j * m] - pole;
        }
        for (value, column) in
            a[i * n..(i + 1) * n].iter_mut().zip(simi.chunks_exact(n))
        {
            *value = difference
                .iter()
                .zip(column)
                .fold(F::zero(), |s, (&d, &v)| s + d * v);
        }
    }
}