Main function for optimization
(&mut self)
| 240 | |
| 241 | /// Main function for optimization |
| 242 | pub fn optimize(&mut self) -> Vec<f64> { |
| 243 | // Receive initial data |
| 244 | let (x_vec, y_vec) = (self.domain.clone(), self.observed.clone()); |
| 245 | let (p_init, max_iter) = (self.param.clone(), self.max_iter); |
| 246 | let safe_f = |p: &Vec<AD>| (self.func)(&x_vec, p.clone()).unwrap(); |
| 247 | let unsafe_f = |p: Vec<AD>| (self.func)(&x_vec, p); |
| 248 | |
| 249 | // Take various form of initial data |
| 250 | let p_init_vec = p_init.to_f64_vec(); |
| 251 | let y = y_vec.to_col(); |
| 252 | |
| 253 | // Declare mutable values |
| 254 | let mut p: Matrix = p_init_vec.clone().into(); |
| 255 | let mut j = jacobian(safe_f, &p_init_vec); |
| 256 | let mut y_hat: Matrix = safe_f(&p_init).to_f64_vec().into(); |
| 257 | let mut jtj = &j.t() * &j; |
| 258 | let mut valid_p = p.clone(); |
| 259 | let mut err_stack = 0usize; |
| 260 | |
| 261 | match self.method { |
| 262 | GradientDescent => { |
| 263 | let alpha = *self.hyperparams.get("lr").unwrap_or(&1e-3); |
| 264 | for i in 0..max_iter { |
| 265 | let h = alpha * j.t() * (&y - &y_hat); |
| 266 | let p_cand = &p + &h; |
| 267 | match unsafe_f(p_cand.data.to_ad_vec()) { |
| 268 | Some(value) => { |
| 269 | p = p_cand; |
| 270 | valid_p = p.clone(); |
| 271 | err_stack = 0; |
| 272 | j = jacobian(safe_f, &p.data); |
| 273 | y_hat = value.to_f64_vec().into(); |
| 274 | } |
| 275 | None => { |
| 276 | if i < max_iter - 1 && err_stack < 3 { |
| 277 | p = p_cand; |
| 278 | err_stack += 1; |
| 279 | } else { |
| 280 | p = valid_p; |
| 281 | break; |
| 282 | } |
| 283 | } |
| 284 | } |
| 285 | } |
| 286 | } |
| 287 | |
| 288 | GaussNewton => unimplemented!(), |
| 289 | |
| 290 | LevenbergMarquardt => { |
| 291 | let mut chi2 = ((&y - &y_hat).t() * (&y - &y_hat))[(0, 0)]; |
| 292 | let mut nu = 2f64; |
| 293 | let lambda_0 = *self.hyperparams.get("lambda_init").unwrap_or(&1e-3); |
| 294 | let lambda_max = *self |
| 295 | .hyperparams |
| 296 | .get("lambda_max") |
| 297 | .unwrap_or(&f64::MAX.sqrt()); |
| 298 | |
| 299 | let mut lambda = lambda_0 * max(jtj.diag()); |