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
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
use crate::params::CmaesParams;
use anyhow::Result;
use nalgebra::{DMatrix, DVector, Dyn, SymmetricEigen};
use statrs::distribution::MultivariateNormal;
/// Structure to hold state for CMA-ES.
#[derive(Debug, Clone)]
pub struct CmaesState {
pub normal_distr: MultivariateNormal<Dyn>, // Normal distributnio function
pub z: DMatrix<f32>, // Matrix of standard normal random variables.
pub y: DMatrix<f32>, // Matrix of candidate solutions.
pub best_y: DVector<f32>, // Best candidate.
pub best_y_fit: DVector<f32>, // Fitness value of the best candidate.
pub best_y_hist: Vec<f32>, // Historical fitness values of the best candidate.
pub cov: DMatrix<f32>, // Covariance matrix of the population.
pub eig_vecs: DMatrix<f32>, // Eigenvectors of the covariance matrix.
pub eig_vals: DVector<f32>, // Eigenvalues of the covariance matrix.
pub inv_sqrt: DMatrix<f32>, // Matrix for the inverse square root of the covariance matrix.
pub mean: DVector<f32>, // Mean of the population.
pub sigma: f32, // Step-size (standard deviation).
pub g: i32, // Curren generation.
pub evals_count: i32, // Number of evaluations performed.
pub ps: DVector<f32>, // Evolution path for step-size adaptation.
pub pc: DVector<f32>, // Evolution path for covariance matrix adaptation.
////////////////
// TODO
// Allow flag for verbose state, maybe with tracing
////////////////
}
/// Trait for CMA-ES State
pub trait CmaesStateLogic {
type NewState;
fn init_state(params: &CmaesParams) -> Result<Self::NewState>;
fn prepare_ask(&mut self, params: &CmaesParams) -> Result<()>;
fn eigen_decomposition(&mut self, params: &CmaesParams) -> Result<()>;
fn get_best(&self) -> Result<(DVector<f32>, DVector<f32>)>;
}
/// Implementing Trait for CMA-ES State
impl CmaesStateLogic for CmaesState {
type NewState = CmaesState;
/// Initiates a new CMA-ES state
///
/// Example
///
/// ```rust
/// use haru_cmaes::params::{CmaesParams, CmaesParamsValidator};
/// use haru_cmaes::strategy::{CmaesAlgo, CmaesAlgoOptimizer};
/// use haru_cmaes::state::{CmaesState, CmaesStateLogic};
///
/// let params = CmaesParams::new().unwrap();
/// let cmaes = CmaesAlgo::new(params).unwrap();
/// let state = CmaesState::init_state(&cmaes.params);
/// assert!(state.is_ok());
/// ```
fn init_state(params: &CmaesParams) -> Result<Self::NewState> {
// Create initial values for the state
let vec = vec![0.0; params.xstart.len()];
let normal_mean = DVector::from_vec(vec);
let normal_cov = DMatrix::identity(params.xstart.len(), params.xstart.len());
let normal_distr = MultivariateNormal::new_from_nalgebra(normal_mean, normal_cov).unwrap();
let z: DMatrix<f32> = DMatrix::zeros(params.popsize as usize, params.xstart.len());
let y: DMatrix<f32> = DMatrix::zeros(params.popsize as usize, params.xstart.len());
let best_y: DVector<f32> = DVector::zeros(params.xstart.len());
let best_y_fit: DVector<f32> = DVector::from_element(1, f32::MAX);
let best_y_hist: Vec<f32> = Vec::with_capacity(50);
let cov: DMatrix<f32> = DMatrix::identity(params.xstart.len(), params.xstart.len());
let inv_sqrt: DMatrix<f32> = DMatrix::identity(params.xstart.len(), params.xstart.len());
let eig_vecs: DMatrix<f32> = DMatrix::identity(params.xstart.len(), params.xstart.len());
let eig_vals: DVector<f32> = DVector::identity(params.xstart.len());
let mean: DVector<f32> = DVector::from_vec(params.xstart.clone());
let sigma: f32 = params.sigma;
let g: i32 = 0;
let evals_count = 0;
let ps: DVector<f32> = DVector::zeros(params.xstart.len());
let pc: DVector<f32> = DVector::zeros(params.xstart.len());
Ok(CmaesState {
normal_distr,
z,
y,
best_y,
best_y_fit,
best_y_hist,
cov,
eig_vecs,
eig_vals,
inv_sqrt,
mean,
sigma,
g,
evals_count,
ps,
pc,
})
}
/// Prepares covariance, eignevalues and eigenvectors.
fn prepare_ask(&mut self, params: &CmaesParams) -> Result<()> {
let _ = self.eigen_decomposition(params);
Ok(())
}
/// Performs eigen decomposition on the covariance matrix.
fn eigen_decomposition(&mut self, params: &CmaesParams) -> Result<()> {
// Ensure symmetric covariance
self.cov = (&self.cov + &self.cov.transpose()) / 2.0;
// For matrix eigen computation efficiency
// non-diag -> enforce sparsity
// diag -> ensure positive
// current fast approach, i.e. for_each (.map is slow)
for i in 0..self.cov.nrows() {
for j in 0..self.cov.nrows() {
let x = self.cov.index_mut((i, j));
if i == j {
// Diagonal entry
if *x < 0.0 {
*x = 0.1 // + (fastrand::f32() / 100.0);
}
} else if params.only_diag {
// Non-diagonal entry
*x = 0.0;
}
}
}
// println!("cov {:?}", &self.cov);
////////////////
// dbg!(&self);
// dbg!(&self);
// println!("");
////////////////
// Perform eigen decomposition: C = B * Λ * B^T
#[cfg(any(
feature = "openblas",
feature = "netlib",
feature = "accelerate",
feature = "intel-mkl"
))]
let mut eigen = nalgebra_lapack::SymmetricEigen::try_new(self.cov.clone())?;
#[cfg(not(any(
feature = "openblas",
feature = "netlib",
feature = "accelerate",
feature = "intel-mkl"
)))]
let eigen = SymmetricEigen::try_new(self.cov.clone(), 1e-20, 0).unwrap();
let mut eig_vals: DVector<f32> = eigen.eigenvalues;
let eig_vecs: DMatrix<f32> = eigen.eigenvectors;
// Ensure positive eigenvalues
eig_vals.iter_mut().for_each(|val| {
if *val < 0.0 {
*val = 0.1; // Adjust negative valenvalues
} else if *val > 10.0 {
*val = 10.0; // Clamp to a maximum value
}
});
// // Clamp eigenvector elements to a range to avoid extreme values
// eig_vecs.iter_mut().for_each(|vec| {
// if *vec < -10.0 {
// *vec = -10.0; // Minimum limit
// } else if *vec > 10.0 {
// *vec = 10.0; // Maximum limit
// }
// });
// Calculate the inverse square root of eigenvalues
let inv_sqrt_diag = DMatrix::from_diagonal(&eig_vals.map(|eig| eig.powf(-0.5)));
self.inv_sqrt = &eig_vecs * &inv_sqrt_diag * eig_vecs.transpose();
// print!("{:?} ", &self.cov.mean());
// print!("{:?} ", &eig_vals.mean());
// print!("{:?} ", &eig_vecs.mean());
// println!("{:?}", &inv_sqrt_diag.mean());
// println!("{}", &eig_vecs);
// Store
self.eig_vecs = eig_vecs;
self.eig_vals = eig_vals;
Ok(())
}
fn get_best(&self) -> Result<(DVector<f32>, DVector<f32>)> {
Ok((self.best_y.clone(), self.best_y_fit.clone()))
}
}