Skip to content

Commit

Permalink
refactor!: always use dynamic dimension
Browse files Browse the repository at this point in the history
  • Loading branch information
pnevyk committed Nov 12, 2023
1 parent 606da4f commit 4f8f2ce
Show file tree
Hide file tree
Showing 15 changed files with 368 additions and 539 deletions.
19 changes: 9 additions & 10 deletions examples/rosenbrock.rs
Original file line number Diff line number Diff line change
@@ -1,7 +1,7 @@
use gomez::nalgebra as na;
use gomez::prelude::*;
use gomez::solver::TrustRegion;
use na::{DimName, IsContiguous};
use na::{Dynamic, IsContiguous};

// https://en.wikipedia.org/wiki/Rosenbrock_function
struct Rosenbrock {
Expand All @@ -11,22 +11,21 @@ struct Rosenbrock {

impl Problem for Rosenbrock {
type Scalar = f64;
type Dim = na::U2;

fn dim(&self) -> Self::Dim {
na::U2::name()
fn domain(&self) -> Domain<Self::Scalar> {
Domain::unconstrained(2)
}
}

impl System for Rosenbrock {
fn eval<Sx, Sfx>(
&self,
x: &na::Vector<Self::Scalar, Self::Dim, Sx>,
fx: &mut na::Vector<Self::Scalar, Self::Dim, Sfx>,
x: &na::Vector<Self::Scalar, Dynamic, Sx>,
fx: &mut na::Vector<Self::Scalar, Dynamic, Sfx>,
) -> Result<(), ProblemError>
where
Sx: na::storage::Storage<Self::Scalar, Self::Dim> + IsContiguous,
Sfx: na::storage::StorageMut<Self::Scalar, Self::Dim>,
Sx: na::storage::Storage<Self::Scalar, Dynamic> + IsContiguous,
Sfx: na::storage::StorageMut<Self::Scalar, Dynamic>,
{
fx[0] = (self.a - x[0]).powi(2);
fx[1] = self.b * (x[1] - x[0].powi(2)).powi(2);
Expand All @@ -41,9 +40,9 @@ fn main() -> Result<(), String> {
let mut solver = TrustRegion::new(&f, &dom);

// Initial guess.
let mut x = na::vector![-10.0, -5.0];
let mut x = na::dvector![-10.0, -5.0];

let mut fx = na::vector![0.0, 0.0];
let mut fx = na::dvector![0.0, 0.0];

for i in 1..=100 {
solver
Expand Down
50 changes: 16 additions & 34 deletions gomez-bench/benches/main.rs
Original file line number Diff line number Diff line change
Expand Up @@ -4,7 +4,7 @@ fn main() {

use gomez::{
nalgebra as na,
nalgebra::{allocator::Allocator, DefaultAllocator},
nalgebra::Dynamic,
prelude::*,
solver::{NelderMead, TrustRegion},
testing::*,
Expand Down Expand Up @@ -168,14 +168,9 @@ mod bullard_biegler {
fn bench_solve<F, S, GF, GS>(bencher: divan::Bencher, with_system: GF, with_solver: GS)
where
GF: Fn() -> (F, usize),
GS: Fn(&F, &Domain<F::Scalar>, &na::OVector<F::Scalar, F::Dim>) -> S,
GS: Fn(&F, &Domain<F::Scalar>, &na::OVector<F::Scalar, Dynamic>) -> S,
F: TestSystem,
S: Solver<F>,
DefaultAllocator: Allocator<F::Scalar, F::Dim>,
DefaultAllocator: Allocator<F::Scalar, F::Dim, F::Dim>,
F::Dim: na::DimMin<F::Dim, Output = F::Dim>,
DefaultAllocator: Allocator<F::Scalar, <F::Dim as na::DimMin<F::Dim>>::Output>,
DefaultAllocator: na::allocator::Reallocator<F::Scalar, F::Dim, F::Dim, F::Dim, F::Dim>,
{
bencher
.with_inputs(move || {
Expand Down Expand Up @@ -209,36 +204,32 @@ where
fn with_trust_region<F>(
f: &F,
dom: &Domain<F::Scalar>,
_: &na::OVector<F::Scalar, F::Dim>,
_: &na::OVector<F::Scalar, Dynamic>,
) -> TrustRegion<F>
where
F: Problem,
DefaultAllocator: Allocator<F::Scalar, F::Dim>,
DefaultAllocator: Allocator<F::Scalar, F::Dim, F::Dim>,
{
TrustRegion::new(f, dom)
}

fn with_nelder_mead<F>(
f: &F,
dom: &Domain<F::Scalar>,
_: &na::OVector<F::Scalar, F::Dim>,
_: &na::OVector<F::Scalar, Dynamic>,
) -> NelderMead<F>
where
F: Problem,
DefaultAllocator: Allocator<F::Scalar, F::Dim>,
{
NelderMead::new(f, dom)
}

fn with_gsl_hybrids<F>(
f: &F,
_: &Domain<F::Scalar>,
x: &na::OVector<F::Scalar, F::Dim>,
x: &na::OVector<F::Scalar, Dynamic>,
) -> GslSolverWrapper<GslFunctionWrapper<F>>
where
F: TestSystem<Scalar = f64> + Clone,
DefaultAllocator: Allocator<F::Scalar, F::Dim>,
{
GslSolverWrapper::new(GslFunctionWrapper::new(
f.clone(),
Expand All @@ -257,21 +248,19 @@ impl<F> GslFunctionWrapper<F> {
}
}

impl<F: TestSystem<Scalar = f64>> GslFunction for GslFunctionWrapper<F>
where
DefaultAllocator: Allocator<F::Scalar, F::Dim>,
{
impl<F: TestSystem<Scalar = f64>> GslFunction for GslFunctionWrapper<F> {
fn eval(&self, x: &GslVec, f: &mut GslVec) -> GslStatus {
use na::DimName;
let dim = Dynamic::new(x.len());

let x = na::MatrixSlice::<f64, F::Dim, na::U1>::from_slice_generic(
let x = na::MatrixSlice::<f64, Dynamic, na::U1>::from_slice_generic(
x.as_slice(),
self.f.dim(),
dim,
na::U1::name(),
);
let mut fx = na::MatrixSliceMut::<f64, F::Dim, na::U1>::from_slice_generic(
let mut fx = na::MatrixSliceMut::<f64, Dynamic, na::U1>::from_slice_generic(
f.as_mut_slice(),
self.f.dim(),
dim,
na::U1::name(),
);

Expand All @@ -298,14 +287,7 @@ impl<F: GslFunction> GslSolverWrapper<F> {
}
}

impl<F: TestSystem<Scalar = f64>> Solver<F> for GslSolverWrapper<GslFunctionWrapper<F>>
where
DefaultAllocator: Allocator<F::Scalar, F::Dim>,
DefaultAllocator: Allocator<F::Scalar, F::Dim, F::Dim>,
F::Dim: na::DimMin<F::Dim, Output = F::Dim>,
DefaultAllocator: Allocator<F::Scalar, <F::Dim as na::DimMin<F::Dim>>::Output>,
DefaultAllocator: na::allocator::Reallocator<F::Scalar, F::Dim, F::Dim, F::Dim, F::Dim>,
{
impl<F: TestSystem<Scalar = f64>> Solver<F> for GslSolverWrapper<GslFunctionWrapper<F>> {
const NAME: &'static str = "GSL hybrids";

type Error = String;
Expand All @@ -314,12 +296,12 @@ where
&mut self,
_f: &F,
_dom: &Domain<F::Scalar>,
x: &mut na::Vector<F::Scalar, F::Dim, Sx>,
fx: &mut na::Vector<F::Scalar, F::Dim, Sfx>,
x: &mut na::Vector<F::Scalar, Dynamic, Sx>,
fx: &mut na::Vector<F::Scalar, Dynamic, Sfx>,
) -> Result<(), Self::Error>
where
Sx: na::storage::StorageMut<F::Scalar, F::Dim> + IsContiguous,
Sfx: na::storage::StorageMut<F::Scalar, F::Dim>,
Sx: na::storage::StorageMut<F::Scalar, Dynamic> + IsContiguous,
Sfx: na::storage::StorageMut<F::Scalar, Dynamic>,
{
let result = self.solver.step().to_result();
x.copy_from_slice(self.solver.root());
Expand Down
24 changes: 9 additions & 15 deletions src/analysis/initial.rs
Original file line number Diff line number Diff line change
Expand Up @@ -10,8 +10,7 @@
use std::marker::PhantomData;

use nalgebra::{
allocator::Allocator, convert, storage::StorageMut, ComplexField, DefaultAllocator, DimMin,
DimName, IsContiguous, OVector, Vector, U1,
convert, storage::StorageMut, ComplexField, DimName, Dynamic, IsContiguous, OVector, Vector, U1,
};
use thiserror::Error;

Expand All @@ -37,28 +36,23 @@ pub struct InitialGuessAnalysis<F: Problem> {
ty: PhantomData<F>,
}

impl<F: System> InitialGuessAnalysis<F>
where
DefaultAllocator: Allocator<F::Scalar, F::Dim>,
DefaultAllocator: Allocator<F::Scalar, F::Dim, F::Dim>,
F::Dim: DimMin<F::Dim, Output = F::Dim>,
DefaultAllocator: Allocator<F::Scalar, <F::Dim as DimMin<F::Dim>>::Output>,
{
impl<F: System> InitialGuessAnalysis<F> {
/// Analyze the system in given point.
pub fn analyze<Sx, Sfx>(
f: &F,
dom: &Domain<F::Scalar>,
x: &mut Vector<F::Scalar, F::Dim, Sx>,
fx: &mut Vector<F::Scalar, F::Dim, Sfx>,
x: &mut Vector<F::Scalar, Dynamic, Sx>,
fx: &mut Vector<F::Scalar, Dynamic, Sfx>,
) -> Result<Self, InitialGuessAnalysisError>
where
Sx: StorageMut<F::Scalar, F::Dim> + IsContiguous,
Sfx: StorageMut<F::Scalar, F::Dim>,
Sx: StorageMut<F::Scalar, Dynamic> + IsContiguous,
Sfx: StorageMut<F::Scalar, Dynamic>,
{
let dim = Dynamic::new(dom.dim());
let scale = dom
.scale()
.map(|scale| OVector::from_iterator_generic(f.dim(), U1::name(), scale.iter().copied()))
.unwrap_or_else(|| OVector::from_element_generic(f.dim(), U1::name(), convert(1.0)));
.map(|scale| OVector::from_iterator_generic(dim, U1::name(), scale.iter().copied()))
.unwrap_or_else(|| OVector::from_element_generic(dim, U1::name(), convert(1.0)));

// Compute F'(x) in the initial point.
f.eval(x, fx)?;
Expand Down
19 changes: 4 additions & 15 deletions src/core/base.rs
Original file line number Diff line number Diff line change
@@ -1,4 +1,4 @@
use nalgebra::{Dim, RealField};
use nalgebra::RealField;
use thiserror::Error;

use super::domain::Domain;
Expand All @@ -9,27 +9,16 @@ pub trait Problem {
/// Type of the scalar, usually f32 or f64.
type Scalar: RealField + Copy;

/// Dimension of the system. Can be fixed
/// ([`Const`](nalgebra::base::dimension::Const)) or dynamic
/// ([`Dynamic`](nalgebra::base::dimension::Dynamic)).
type Dim: Dim;

/// Return the actual dimension of the system. This is needed for dynamic
/// systems.
fn dim(&self) -> Self::Dim;

/// Get the domain (bound constraints) of the system. If not overridden, the
/// system is unconstrained.
fn domain(&self) -> Domain<Self::Scalar> {
Domain::unconstrained(self.dim().value())
}
fn domain(&self) -> Domain<Self::Scalar>;
}

/// Error encountered while applying variables to the problem.
#[derive(Debug, Error)]
pub enum ProblemError {
/// The number of variables does not match the dimensionality
/// ([`Problem::dim`]) of the problem.
/// The number of variables does not match the dimensionality of the problem
/// domain.
#[error("invalid dimensionality")]
InvalidDimensionality,
/// An invalid value (NaN, positive or negative infinity) of a residual or
Expand Down
51 changes: 23 additions & 28 deletions src/core/function.rs
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
use nalgebra::{
allocator::Allocator, storage::Storage, storage::StorageMut, DefaultAllocator, IsContiguous,
Vector,
allocator::Allocator, storage::Storage, storage::StorageMut, DefaultAllocator, Dynamic,
IsContiguous, Vector,
};
use num_traits::Zero;

Expand All @@ -14,14 +14,14 @@ use super::{
/// ## Defining a function
///
/// A function is any type that implements [`Function`] and [`Problem`] traits.
/// There are two required associated types (scalar type and dimension type) and
/// two required methods: [`apply`](Function::apply) and [`dim`](Problem::dim).
/// There is one required associated type (scalar) and one required method
/// ([`apply`](Function::apply)).
///
/// ```rust
/// use gomez::core::Function;
/// use gomez::nalgebra as na;
/// use gomez::prelude::*;
/// use na::{Dim, DimName, IsContiguous};
/// use na::{Dynamic, IsContiguous};
///
/// // A problem is represented by a type.
/// struct Rosenbrock {
Expand All @@ -32,23 +32,21 @@ use super::{
/// impl Problem for Rosenbrock {
/// // The numeric type. Usually f64 or f32.
/// type Scalar = f64;
/// // The dimension of the problem. Can be either statically known or dynamic.
/// type Dim = na::U2;
///
/// // Return the actual dimension of the system.
/// fn dim(&self) -> Self::Dim {
/// na::U2::name()
/// // The domain of the problem (could be bound-constrained).
/// fn domain(&self) -> Domain<Self::Scalar> {
/// Domain::unconstrained(2)
/// }
/// }
///
/// impl Function for Rosenbrock {
/// // Apply trial values of variables to the function.
/// fn apply<Sx>(
/// &self,
/// x: &na::Vector<Self::Scalar, Self::Dim, Sx>,
/// x: &na::Vector<Self::Scalar, Dynamic, Sx>,
/// ) -> Result<Self::Scalar, ProblemError>
/// where
/// Sx: na::storage::Storage<Self::Scalar, Self::Dim> + IsContiguous,
/// Sx: na::storage::Storage<Self::Scalar, Dynamic> + IsContiguous,
/// {
/// // Compute the function value.
/// Ok((self.a - x[0]).powi(2) + self.b * (x[1] - x[0].powi(2)).powi(2))
Expand All @@ -59,10 +57,10 @@ pub trait Function: Problem {
/// Calculate the function value given values of the variables.
fn apply<Sx>(
&self,
x: &Vector<Self::Scalar, Self::Dim, Sx>,
x: &Vector<Self::Scalar, Dynamic, Sx>,
) -> Result<Self::Scalar, ProblemError>
where
Sx: Storage<Self::Scalar, Self::Dim> + IsContiguous;
Sx: Storage<Self::Scalar, Dynamic> + IsContiguous;

/// Calculate the norm of residuals of the system given values of the
/// variable for cases when the function is actually a system of equations.
Expand All @@ -72,12 +70,12 @@ pub trait Function: Problem {
/// do not make an unnecessary allocation for it.
fn apply_eval<Sx, Sfx>(
&self,
x: &Vector<Self::Scalar, Self::Dim, Sx>,
fx: &mut Vector<Self::Scalar, Self::Dim, Sfx>,
x: &Vector<Self::Scalar, Dynamic, Sx>,
fx: &mut Vector<Self::Scalar, Dynamic, Sfx>,
) -> Result<Self::Scalar, ProblemError>
where
Sx: Storage<Self::Scalar, Self::Dim> + IsContiguous,
Sfx: StorageMut<Self::Scalar, Self::Dim>,
Sx: Storage<Self::Scalar, Dynamic> + IsContiguous,
Sfx: StorageMut<Self::Scalar, Dynamic>,
{
let norm = self.apply(x)?;
fx.fill(Self::Scalar::zero());
Expand All @@ -88,27 +86,24 @@ pub trait Function: Problem {

impl<F: System> Function for F
where
DefaultAllocator: Allocator<F::Scalar, F::Dim>,
DefaultAllocator: Allocator<F::Scalar, Dynamic>,
{
fn apply<Sx>(
&self,
x: &Vector<Self::Scalar, Self::Dim, Sx>,
) -> Result<Self::Scalar, ProblemError>
fn apply<Sx>(&self, x: &Vector<Self::Scalar, Dynamic, Sx>) -> Result<Self::Scalar, ProblemError>
where
Sx: Storage<Self::Scalar, Self::Dim> + IsContiguous,
Sx: Storage<Self::Scalar, Dynamic> + IsContiguous,
{
let mut fx = x.clone_owned();
self.apply_eval(x, &mut fx)
}

fn apply_eval<Sx, Sfx>(
&self,
x: &Vector<Self::Scalar, Self::Dim, Sx>,
fx: &mut Vector<Self::Scalar, Self::Dim, Sfx>,
x: &Vector<Self::Scalar, Dynamic, Sx>,
fx: &mut Vector<Self::Scalar, Dynamic, Sfx>,
) -> Result<Self::Scalar, ProblemError>
where
Sx: Storage<Self::Scalar, Self::Dim> + IsContiguous,
Sfx: StorageMut<Self::Scalar, Self::Dim>,
Sx: Storage<Self::Scalar, Dynamic> + IsContiguous,
Sfx: StorageMut<Self::Scalar, Dynamic>,
{
self.eval(x, fx)?;
let norm = fx.norm();
Expand Down
Loading

0 comments on commit 4f8f2ce

Please sign in to comment.