ceres_solver/nlls_problem.rs
1//! Non-Linear Least Squares problem builder and solver.
2//!
3//! The diagram shows the lifecycle of a [NllsProblem]:
4//! ```text
5//! x
6//! │ NllsProblem::new()
7//! │
8//! ┌──▼────────┐ .solve(self, options) ┌───────────────────┐
9//! ┌──►│NllsProblem├──────────────────────►│NllsProblemSolution│
10//! │ └──┬────────┘ └───────────────────┘
11//! │ │ .residual_block_builder(self)
12//! │ ┌──▼─────────────────┐
13//! │ │ResidualBlockBuilder│
14//! │ └──▲─┬───────────────┘
15//! │ │ │.set_cost(self, func, num_residuals)
16//! │ └─┤
17//! │ ▲ │
18//! │ └─┤.set_loss(self, loss)
19//! │ │
20//! │ ▲ │
21//! │ └─┤.set_parameters(self,
22//! │ │
23//! └────────┘.build_into_problem(self)
24//! ```
25//! <!-- https://asciiflow.com/#/share/eJytU1tqg0AU3cpwvyKIpPlq%2FcwCSml%2FB0TNDUivM2EerSFkF8WF5LN0NV1JR502NRpbSIajnEHuOXfOHXcg0hIhFpYoBEq3qCCGHYeKQ3wzny9CDltHF7d3jhmsjNtwYN2qOBeefr59sHsi%2FaBkRljGscDXWdD77jeOedTvRz4Ai7SkF5xppHXI5MYUUuiATVT8B66HE%2F9fTV%2BofUb1SZJtmv9x72UwCTa%2BrpLBcWwsUqiLlU0pyUjmz0lmC1qhaqMPRnquLzV270dvuWwcl53heEL14TrHdE%2Bk0SS51MbfqrUVeciELZPvBHRwUjaiVR%2FYUL5Da0BSa2%2FQ0J4iG5T%2BpbZJlftDDSqveUZtAlE7z6QQRvqRwh72XxPrJEg%3D) -->
26//!
27//! We start with [NllsProblem] with no residual blocks and cannot be solved. Next we should add a
28//! residual block by calling [NllsProblem::residual_block_builder] which is a destructive method
29//! which consumes the problem and returns a [ResidualBlockBuilder] which can be used to build a new
30//! residual block. Here we add mandatory cost function [crate::cost::CostFunctionType] and
31//! parameter blocks [crate::parameter_block::ParameterBlock]. We can also set optional loss
32//! function [crate::loss::LossFunction]. Once we are done, we call
33//! [ResidualBlockBuilder::build_into_problem] which returns previously consumed [NllsProblem].
34//! Now we can optionally add more residual blocks repeating the process: call
35//! [NllsProblem::residual_block_builder] consuming [NllsProblem], add what we need and rebuild the
36//! problem. The only difference that now we can re-use parameter blocks used in the previous
37//! residual blocks, adding them by their indexes. Once we are done, we can call
38//! [NllsProblem::solve] which consumes the problem, solves it and returns [NllsProblemSolution]
39//! which contains the solution and summary of the solver run. It returns an error if the problem
40//! has no residual blocks.
41//!
42//! # Examples
43//!
44//! ## Multiple residual blocks with shared parameters
45//!
46//! Let's solve a problem of fitting a family of functions `y_ij = a + b_i * exp(c_i * x_ij)`:
47//! all of them have the same offset `a`, but different scale parameters `b_i` and `c_i`,
48//! `i in 0..=k-1` for `k` (`N_CURVES` bellow) different sets of data.
49//!
50//! ```rust
51//! use ceres_solver::parameter_block::ParameterBlockOrIndex;
52//! use ceres_solver::{CostFunctionType, NllsProblem, SolverOptions};
53//!
54//! // Get parameters, x, y and return tuple of function value and its derivatives
55//! fn target_function(parameters: &[f64; 3], x: f64) -> (f64, [f64; 3]) {
56//! let [a, b, c] = parameters;
57//! let y = a + b * f64::exp(c * x);
58//! let dy_da = 1.0;
59//! let dy_db = f64::exp(c * x);
60//! let dy_dc = b * x * f64::exp(c * x);
61//! (y, [dy_da, dy_db, dy_dc])
62//! }
63//!
64//! const N_OBS_PER_CURVE: usize = 100;
65//! const N_CURVES: usize = 3;
66//!
67//! // True parameters
68//! let a_true = -2.0;
69//! let b_true: [_; N_CURVES] = [2.0, 2.0, -1.0];
70//! let c_true: [_; N_CURVES] = [3.0, -1.0, 3.0];
71//!
72//! // Initial parameter guesses
73//! let a_init = 0.0;
74//! let b_init = 1.0;
75//! let c_init = 1.0;
76//!
77//! // Generate data
78//! let x = vec![
79//! (0..N_OBS_PER_CURVE)
80//! .map(|i| (i as f64) / (N_OBS_PER_CURVE as f64))
81//! .collect::<Vec<_>>();
82//! 3
83//! ];
84//! let y: Vec<Vec<_>> = x
85//! .iter()
86//! .zip(b_true.iter().zip(c_true.iter()))
87//! .map(|(x, (&b, &c))| {
88//! x.iter()
89//! .map(|&x| {
90//! let (y, _) = target_function(&[a_true, b, c], x);
91//! // True value + "noise"
92//! y + 0.001 + f64::sin(1e6 * x)
93//! })
94//! .collect()
95//! })
96//! .collect();
97//!
98//! // Build the problem
99//! let mut problem = NllsProblem::new();
100//! for (i, (x, y)) in x.into_iter().zip(y.into_iter()).enumerate() {
101//! let cost: CostFunctionType = Box::new(
102//! move |parameters: &[&[f64]],
103//! residuals: &mut [f64],
104//! mut jacobians: Option<&mut [Option<&mut [&mut [f64]]>]>| {
105//! assert_eq!(parameters.len(), 3);
106//! let a = parameters[0][0];
107//! let b = parameters[1][0];
108//! let c = parameters[2][0];
109//! // Number of residuls equal to the number of observations
110//! assert_eq!(residuals.len(), N_OBS_PER_CURVE);
111//! for (j, (&x, &y)) in x.iter().zip(y.iter()).enumerate() {
112//! let (y_model, derivatives) = target_function(&[a, b, c], x);
113//! residuals[j] = y - y_model;
114//! // jacobians can be None, then you don't need to provide them
115//! if let Some(jacobians) = jacobians.as_mut() {
116//! // The size of the jacobians array is equal to the number of parameters,
117//! // each element is Option<&mut [&mut [f64]]>
118//! for (mut jacobian, &derivative) in jacobians.iter_mut().zip(&derivatives) {
119//! if let Some(jacobian) = &mut jacobian {
120//! // Each element in the jacobians array is slice of slices:
121//! // the first index is for different residuals components,
122//! // the second index is for different components of the parameter vector
123//! jacobian[j][0] = -derivative;
124//! }
125//! }
126//! }
127//! }
128//! true
129//! },
130//! );
131//! let a_parameter: ParameterBlockOrIndex = if i == 0 {
132//! vec![c_init].into()
133//! } else {
134//! 0.into()
135//! };
136//! problem = problem
137//! .residual_block_builder()
138//! .set_cost(cost, N_OBS_PER_CURVE)
139//! .add_parameter(a_parameter)
140//! .add_parameter(vec![b_init])
141//! .add_parameter(vec![c_init])
142//! .build_into_problem()
143//! .unwrap()
144//! .0;
145//! }
146//!
147//! // Solve the problem
148//! let solution = problem.solve(&SolverOptions::default()).unwrap();
149//! println!("Brief summary: {:?}", solution.summary);
150//! // Getting parameter values
151//! let a = solution.parameters[0][0];
152//! assert!((a - a_true).abs() < 0.03);
153//! let (b, c): (Vec<_>, Vec<_>) = solution.parameters[1..]
154//! .chunks(2)
155//! .map(|sl| (sl[0][0], sl[1][0]))
156//! .unzip();
157//! for (b, &b_true) in b.iter().zip(b_true.iter()) {
158//! assert!((b - b_true).abs() < 0.03);
159//! }
160//! for (c, &c_true) in c.iter().zip(c_true.iter()) {
161//! assert!((c - c_true).abs() < 0.03);
162//! }
163//! ```
164//!
165//! ## Parameter constraints
166//!
167//! Let's find a minimum of the Himmelblau's function:
168//! `f(x, y) = (x^2 + y - 11)^2 + (x + y^2 - 7)^2` with boundaries `x ∈ [0; 3.5], y ∈ [-1.8; 3.5]`
169//! and initial guess `x = 3.45, y = -1.8`. This function have four global minima, all having the
170//! the same value `f(x, y) = 0`, one of them is within the boundaries and another one is just
171//! outside of them, near the initial guess. The solver converges to the corner of the boundary.
172//!
173//! ```rust
174//! use ceres_solver::{CostFunctionType, NllsProblem, ParameterBlock, SolverOptions};
175//!
176//! const LOWER_X: f64 = 0.0;
177//! const UPPER_X: f64 = 3.5;
178//! const LOWER_Y: f64 = -1.8;
179//! const UPPER_Y: f64 = 3.5;
180//!
181//! fn solve_himmelblau(initial_x: f64, initial_y: f64) -> (f64, f64) {
182//! let x_block = {
183//! let mut block = ParameterBlock::new(vec![initial_x]);
184//! block.set_all_lower_bounds(vec![LOWER_X]);
185//! block.set_all_upper_bounds(vec![UPPER_X]);
186//! block
187//! };
188//! let y_block = {
189//! let mut block = ParameterBlock::new(vec![initial_y]);
190//! block.set_all_lower_bounds(vec![LOWER_Y]);
191//! block.set_all_upper_bounds(vec![UPPER_Y]);
192//! block
193//! };
194//!
195//! // You can skip type annotations in the closure definition, we use them for verbosity only.
196//! let cost: CostFunctionType = Box::new(
197//! move |parameters: &[&[f64]],
198//! residuals: &mut [f64],
199//! mut jacobians: Option<&mut [Option<&mut [&mut [f64]]>]>| {
200//! let x = parameters[0][0];
201//! let y = parameters[1][0];
202//! // residuals have the size of your data set, in our case it is two
203//! residuals[0] = x.powi(2) + y - 11.0;
204//! residuals[1] = x + y.powi(2) - 7.0;
205//! // jacobians can be None, then you don't need to provide them
206//! if let Some(jacobians) = jacobians {
207//! // The size of the jacobians array is equal to the number of parameters,
208//! // each element is Option<&mut [&mut [f64]]>
209//! if let Some(d_dx) = &mut jacobians[0] {
210//! // Each element in the jacobians array is slice of slices:
211//! // the first index is for different residuals components,
212//! // the second index is for different components of the parameter vector
213//! d_dx[0][0] = 2.0 * x;
214//! d_dx[1][0] = 1.0;
215//! }
216//! if let Some(d_dy) = &mut jacobians[1] {
217//! d_dy[0][0] = 1.0;
218//! d_dy[1][0] = 2.0 * y;
219//! }
220//! }
221//! true
222//! },
223//! );
224//!
225//! let solution = NllsProblem::new()
226//! .residual_block_builder() // create a builder for residual block
227//! .set_cost(cost, 2) // 2 is the number of residuals
228//! .set_parameters([x_block, y_block])
229//! .build_into_problem()
230//! .unwrap()
231//! .0 // build_into_problem returns a tuple (NllsProblem, ResidualBlockId)
232//! .solve(&SolverOptions::default()) // SolverOptions can be customized
233//! .unwrap(); // Err should happen only if we added no residual blocks
234//!
235//! // Print the full solver report
236//! println!("{}", solution.summary.full_report());
237//!
238//! (solution.parameters[0][0], solution.parameters[1][0])
239//! }
240//!
241//! // The solver converges to the corner of the boundary rectangle.
242//! let (x, y) = solve_himmelblau(3.4, -1.0);
243//! assert_eq!(UPPER_X, x);
244//! assert_eq!(LOWER_Y, y);
245//!
246//! // The solver converges to the global minimum inside the boundaries.
247//! let (x, y) = solve_himmelblau(1.0, 1.0);
248//! assert!((3.0 - x).abs() < 1e-8);
249//! assert!((2.0 - y).abs() < 1e-8);
250//! ```
251
252use crate::cost::CostFunction;
253use crate::cost::CostFunctionType;
254use crate::error::{NllsProblemError, ParameterBlockStorageError, ResidualBlockBuildingError};
255use crate::loss::LossFunction;
256use crate::parameter_block::{ParameterBlockOrIndex, ParameterBlockStorage};
257use crate::residual_block::{ResidualBlock, ResidualBlockId};
258use crate::solver::{SolverOptions, SolverSummary};
259
260use ceres_solver_sys::cxx::UniquePtr;
261use ceres_solver_sys::ffi;
262use std::pin::Pin;
263
264/// Non-Linear Least Squares problem.
265///
266/// See [module-level documentation](crate::nlls_problem) building the instance of this type.
267pub struct NllsProblem<'cost> {
268 inner: UniquePtr<ffi::Problem<'cost>>,
269 parameter_storage: ParameterBlockStorage,
270 residual_blocks: Vec<ResidualBlock>,
271}
272
273impl<'cost> NllsProblem<'cost> {
274 /// Crate a new non-linear least squares problem with no residual blocks.
275 pub fn new() -> Self {
276 Self {
277 inner: ffi::new_problem(),
278 parameter_storage: ParameterBlockStorage::new(),
279 residual_blocks: Vec::new(),
280 }
281 }
282
283 /// Capture this problem into a builder for a new residual block.
284 pub fn residual_block_builder(self) -> ResidualBlockBuilder<'cost> {
285 ResidualBlockBuilder {
286 problem: self,
287 cost: None,
288 loss: None,
289 parameters: Vec::new(),
290 }
291 }
292
293 #[inline]
294 fn inner(&self) -> &ffi::Problem<'cost> {
295 self.inner
296 .as_ref()
297 .expect("Underlying C++ unique_ptr<Problem> must hold non-null pointer")
298 }
299
300 #[inline]
301 fn inner_mut(&mut self) -> Pin<&mut ffi::Problem<'cost>> {
302 self.inner
303 .as_mut()
304 .expect("Underlying C++ unique_ptr<Problem> must hold non-null pointer")
305 }
306
307 /// Set parameter block to be constant during the optimization. Parameter block must be already
308 /// added to the problem, otherwise [ParameterBlockStorageError] returned.
309 pub fn set_parameter_block_constant(
310 &mut self,
311 block_index: usize,
312 ) -> Result<(), ParameterBlockStorageError> {
313 let block_pointer = self.parameter_storage.get_block(block_index)?.pointer_mut();
314 unsafe {
315 self.inner_mut().SetParameterBlockConstant(block_pointer);
316 }
317 Ok(())
318 }
319
320 /// Set parameter block to be variable during the optimization. Parameter block must be already
321 /// added to the problem, otherwise [ParameterBlockStorageError] returned.
322 pub fn set_parameter_block_variable(
323 &mut self,
324 block_index: usize,
325 ) -> Result<(), ParameterBlockStorageError> {
326 let block_pointer = self.parameter_storage.get_block(block_index)?.pointer_mut();
327 unsafe {
328 self.inner_mut().SetParameterBlockVariable(block_pointer);
329 }
330 Ok(())
331 }
332
333 /// Check if parameter block is constant. Parameter block must be already added to the problem,
334 /// otherwise [ParameterBlockStorageError] returned.
335 pub fn is_parameter_block_constant(
336 &self,
337 block_index: usize,
338 ) -> Result<bool, ParameterBlockStorageError> {
339 let block_pointer = self.parameter_storage.get_block(block_index)?.pointer_mut();
340 unsafe { Ok(self.inner().IsParameterBlockConstant(block_pointer)) }
341 }
342
343 /// Solve the problem.
344 pub fn solve(
345 mut self,
346 options: &SolverOptions,
347 ) -> Result<NllsProblemSolution, NllsProblemError> {
348 if self.residual_blocks.is_empty() {
349 return Err(NllsProblemError::NoResidualBlocks);
350 }
351 let mut summary = SolverSummary::new();
352 ffi::solve(
353 options
354 .0
355 .as_ref()
356 .expect("Underlying C++ SolverOptions must hold non-null pointer"),
357 self.inner_mut(),
358 summary
359 .0
360 .as_mut()
361 .expect("Underlying C++ unique_ptr<SolverSummary> must hold non-null pointer"),
362 );
363 Ok(NllsProblemSolution {
364 parameters: self.parameter_storage.to_values(),
365 summary,
366 })
367 }
368}
369
370impl Default for NllsProblem<'_> {
371 fn default() -> Self {
372 Self::new()
373 }
374}
375
376/// Solution of a non-linear least squares problem [NllsProblem].
377pub struct NllsProblemSolution {
378 /// Values of the parameters, in the same order as they were added to the problem.
379 pub parameters: Vec<Vec<f64>>,
380 /// Summary of the solver run.
381 pub summary: SolverSummary,
382}
383
384/// Builder for a new residual block. It captures [NllsProblem] and returns it back with
385/// [ResidualBlockBuilder::build_into_problem] call.
386pub struct ResidualBlockBuilder<'cost> {
387 problem: NllsProblem<'cost>,
388 cost: Option<(CostFunctionType<'cost>, usize)>,
389 loss: Option<LossFunction>,
390 parameters: Vec<ParameterBlockOrIndex>,
391}
392
393impl<'cost> ResidualBlockBuilder<'cost> {
394 /// Set cost function for the residual block.
395 ///
396 /// Arguments:
397 /// * `func` - cost function, see [CostFunction] for details on how to implement it,
398 /// * `num_residuals` - number of residuals, typically the same as the number of experiments.
399 pub fn set_cost(
400 mut self,
401 func: impl Into<CostFunctionType<'cost>>,
402 num_residuals: usize,
403 ) -> Self {
404 self.cost = Some((func.into(), num_residuals));
405 self
406 }
407
408 /// Set loss function for the residual block.
409 pub fn set_loss(mut self, loss: LossFunction) -> Self {
410 self.loss = Some(loss);
411 self
412 }
413
414 /// Set parameters for the residual block.
415 ///
416 /// The argument is an iterator over [ParameterBlockOrIndex] which can be either a new parameter
417 /// block or an index of an existing parameter block.
418 pub fn set_parameters<P>(mut self, parameters: impl IntoIterator<Item = P>) -> Self
419 where
420 P: Into<ParameterBlockOrIndex>,
421 {
422 self.parameters = parameters.into_iter().map(|p| p.into()).collect();
423 self
424 }
425
426 /// Add a new parameter block to the residual block.
427 ///
428 /// The argument is either a new parameter block or an index of an existing parameter block.
429 pub fn add_parameter<P>(mut self, parameter_block: P) -> Self
430 where
431 P: Into<ParameterBlockOrIndex>,
432 {
433 self.parameters.push(parameter_block.into());
434 self
435 }
436
437 /// Build the residual block, add to the problem and return the problem back.
438 ///
439 /// Returns [ResidualBlockBuildingError] if:
440 /// * cost function is not set,
441 /// * no parameters are set,
442 /// * any of the parameters is not a new parameter block or an index of an existing parameter.
443 ///
444 /// Otherwise returns the problem and the residual block id.
445 pub fn build_into_problem(
446 self,
447 ) -> Result<(NllsProblem<'cost>, ResidualBlockId), ResidualBlockBuildingError> {
448 let Self {
449 mut problem,
450 cost,
451 loss,
452 parameters,
453 } = self;
454 if parameters.is_empty() {
455 return Err(ResidualBlockBuildingError::MissingParameters);
456 }
457 let parameter_indices = problem.parameter_storage.extend(parameters)?;
458 let parameter_sizes: Vec<_> = parameter_indices
459 .iter()
460 // At this point we know that all parameter indices are valid.
461 .map(|&index| problem.parameter_storage.blocks()[index].len())
462 .collect();
463 let parameter_pointers: Pin<Vec<_>> = Pin::new(
464 parameter_indices
465 .iter()
466 // At this point we know that all parameter indices are valid.
467 .map(|&index| problem.parameter_storage.blocks()[index].pointer_mut())
468 .collect(),
469 );
470
471 // Create cost function
472 let cost = if let Some((func, num_redisuals)) = cost {
473 CostFunction::new(func, parameter_sizes, num_redisuals)
474 } else {
475 return Err(ResidualBlockBuildingError::MissingCost);
476 };
477
478 // Set residual block
479 let residual_block_id = unsafe {
480 ffi::add_residual_block(
481 problem
482 .inner
483 .as_mut()
484 .expect("Underlying C++ unique_ptr<Problem> must hold non-null pointer"),
485 cost.into_inner(),
486 loss.map(|loss| loss.into_inner())
487 .unwrap_or_else(UniquePtr::null),
488 parameter_pointers.as_ptr(),
489 parameter_indices.len() as i32,
490 )
491 };
492 problem.residual_blocks.push(ResidualBlock {
493 id: residual_block_id.clone(),
494 parameter_pointers,
495 });
496
497 // Set parameter bounds
498 for &index in parameter_indices.iter() {
499 let block = &problem.parameter_storage.blocks()[index];
500 if let Some(lower_bound) = block.lower_bounds() {
501 for (i, lower_bound) in lower_bound.iter().enumerate() {
502 if let Some(lower_bound) = lower_bound {
503 unsafe {
504 problem
505 .inner
506 .as_mut()
507 .expect(
508 "Underlying C++ unique_ptr<Problem> must hold non-null pointer",
509 )
510 .SetParameterLowerBound(block.pointer_mut(), i as i32, *lower_bound)
511 }
512 }
513 }
514 }
515 }
516 for &index in parameter_indices.iter() {
517 let block = &problem.parameter_storage.blocks()[index];
518 if let Some(upper_bound) = block.upper_bounds() {
519 for (i, upper_bound) in upper_bound.iter().enumerate() {
520 if let Some(upper_bound) = upper_bound {
521 unsafe {
522 problem
523 .inner
524 .as_mut()
525 .expect(
526 "Underlying C++ unique_ptr<Problem> must hold non-null pointer",
527 )
528 .SetParameterUpperBound(block.pointer_mut(), i as i32, *upper_bound)
529 }
530 }
531 }
532 }
533 }
534
535 Ok((problem, residual_block_id))
536 }
537}
538
539#[cfg(test)]
540mod tests {
541 use super::*;
542
543 use crate::cost::CostFunctionType;
544 use crate::loss::{LossFunction, LossFunctionType};
545
546 use approx::assert_abs_diff_eq;
547
548 /// Adopted from c_api_tests.cc, ceres-solver version 2.1.0
549 fn simple_end_to_end_test_with_loss(loss: LossFunction) {
550 const NUM_OBSERVATIONS: usize = 67;
551 const NDIM: usize = 2;
552 let data: [[f64; NDIM]; NUM_OBSERVATIONS] = [
553 0.000000e+00,
554 1.133898e+00,
555 7.500000e-02,
556 1.334902e+00,
557 1.500000e-01,
558 1.213546e+00,
559 2.250000e-01,
560 1.252016e+00,
561 3.000000e-01,
562 1.392265e+00,
563 3.750000e-01,
564 1.314458e+00,
565 4.500000e-01,
566 1.472541e+00,
567 5.250000e-01,
568 1.536218e+00,
569 6.000000e-01,
570 1.355679e+00,
571 6.750000e-01,
572 1.463566e+00,
573 7.500000e-01,
574 1.490201e+00,
575 8.250000e-01,
576 1.658699e+00,
577 9.000000e-01,
578 1.067574e+00,
579 9.750000e-01,
580 1.464629e+00,
581 1.050000e+00,
582 1.402653e+00,
583 1.125000e+00,
584 1.713141e+00,
585 1.200000e+00,
586 1.527021e+00,
587 1.275000e+00,
588 1.702632e+00,
589 1.350000e+00,
590 1.423899e+00,
591 1.425000e+00,
592 1.543078e+00,
593 1.500000e+00,
594 1.664015e+00,
595 1.575000e+00,
596 1.732484e+00,
597 1.650000e+00,
598 1.543296e+00,
599 1.725000e+00,
600 1.959523e+00,
601 1.800000e+00,
602 1.685132e+00,
603 1.875000e+00,
604 1.951791e+00,
605 1.950000e+00,
606 2.095346e+00,
607 2.025000e+00,
608 2.361460e+00,
609 2.100000e+00,
610 2.169119e+00,
611 2.175000e+00,
612 2.061745e+00,
613 2.250000e+00,
614 2.178641e+00,
615 2.325000e+00,
616 2.104346e+00,
617 2.400000e+00,
618 2.584470e+00,
619 2.475000e+00,
620 1.914158e+00,
621 2.550000e+00,
622 2.368375e+00,
623 2.625000e+00,
624 2.686125e+00,
625 2.700000e+00,
626 2.712395e+00,
627 2.775000e+00,
628 2.499511e+00,
629 2.850000e+00,
630 2.558897e+00,
631 2.925000e+00,
632 2.309154e+00,
633 3.000000e+00,
634 2.869503e+00,
635 3.075000e+00,
636 3.116645e+00,
637 3.150000e+00,
638 3.094907e+00,
639 3.225000e+00,
640 2.471759e+00,
641 3.300000e+00,
642 3.017131e+00,
643 3.375000e+00,
644 3.232381e+00,
645 3.450000e+00,
646 2.944596e+00,
647 3.525000e+00,
648 3.385343e+00,
649 3.600000e+00,
650 3.199826e+00,
651 3.675000e+00,
652 3.423039e+00,
653 3.750000e+00,
654 3.621552e+00,
655 3.825000e+00,
656 3.559255e+00,
657 3.900000e+00,
658 3.530713e+00,
659 3.975000e+00,
660 3.561766e+00,
661 4.050000e+00,
662 3.544574e+00,
663 4.125000e+00,
664 3.867945e+00,
665 4.200000e+00,
666 4.049776e+00,
667 4.275000e+00,
668 3.885601e+00,
669 4.350000e+00,
670 4.110505e+00,
671 4.425000e+00,
672 4.345320e+00,
673 4.500000e+00,
674 4.161241e+00,
675 4.575000e+00,
676 4.363407e+00,
677 4.650000e+00,
678 4.161576e+00,
679 4.725000e+00,
680 4.619728e+00,
681 4.800000e+00,
682 4.737410e+00,
683 4.875000e+00,
684 4.727863e+00,
685 4.950000e+00,
686 4.669206e+00,
687 ]
688 .chunks_exact(NDIM)
689 .map(|chunk| chunk.try_into().unwrap())
690 .collect::<Vec<_>>()
691 .try_into()
692 .unwrap();
693
694 let cost: CostFunctionType = Box::new(move |parameters, residuals, mut jacobians| {
695 let m = parameters[0][0];
696 let c = parameters[1][0];
697 for ((i, row), residual) in data.into_iter().enumerate().zip(residuals.iter_mut()) {
698 let x = row[0];
699 let y = row[1];
700 *residual = y - f64::exp(m * x + c);
701 if let Some(jacobians) = jacobians.as_mut() {
702 if let Some(d_dm) = jacobians[0].as_mut() {
703 d_dm[i][0] = -x * f64::exp(m * x + c);
704 }
705 if let Some(d_dc) = jacobians[1].as_mut() {
706 d_dc[i][0] = -f64::exp(m * x + c);
707 }
708 }
709 }
710 true
711 });
712
713 let initial_guess = vec![vec![0.0], vec![0.0]];
714
715 let NllsProblemSolution {
716 parameters: solution,
717 summary,
718 } = NllsProblem::new()
719 .residual_block_builder()
720 .set_cost(cost, NUM_OBSERVATIONS)
721 .set_parameters(initial_guess)
722 .set_loss(loss)
723 .build_into_problem()
724 .unwrap()
725 .0
726 .solve(&SolverOptions::default())
727 .unwrap();
728
729 assert!(summary.is_solution_usable());
730 println!("{}", summary.full_report());
731
732 let m = solution[0][0];
733 let c = solution[1][0];
734
735 assert_abs_diff_eq!(0.3, m, epsilon = 0.02);
736 assert_abs_diff_eq!(0.1, c, epsilon = 0.04);
737 }
738
739 #[test]
740 fn simple_end_to_end_test_trivial_custom_loss() {
741 let loss: LossFunctionType = Box::new(|squared_norm: f64, out: &mut [f64; 3]| {
742 out[0] = squared_norm;
743 out[1] = 1.0;
744 out[2] = 0.0;
745 });
746 simple_end_to_end_test_with_loss(LossFunction::custom(loss));
747 }
748
749 #[test]
750 fn simple_end_to_end_test_arctan_stock_loss() {
751 simple_end_to_end_test_with_loss(LossFunction::arctan(1.0));
752 }
753}