#[macro_export]
macro_rules! __eq_helper {
($f:path $(,)?) => {
|genes: &ndarray::Array2<f64>| -> ndarray::Array1<f64> {
$f(genes).mapv(|v| v.abs() - 1e-6)
}
};
}
#[macro_export]
macro_rules! __constraints_helper {
($($c:expr),+ $(,)?) => {
|genes: &ndarray::Array2<f64>| -> ndarray::Array2<f64> {
let cols = vec![ $( ($c)(genes) ),+ ];
let views: Vec<_> = cols.iter()
.map(|v| v.view().insert_axis(ndarray::Axis(1)))
.collect();
ndarray::concatenate(ndarray::Axis(1), &views)
.expect("Failed to concatenate constraints along axis 1")
}
};
}
#[macro_export]
macro_rules! impl_constraints_fn {
(
$name:ident
$(, ineq = [ $($ineq:path),* $(,)? ] )?
$(, eq = [ $($eq:path),* $(,)? ] )?
$(, lower_bound = $lb:expr )?
$(, upper_bound = $ub:expr )?
$(,)?
) => {
#[derive(Debug, Clone, Copy)]
pub struct $name;
impl $crate::ConstraintsFn for $name {
type Dim = ndarray::Ix2;
fn call(&self, genes: &ndarray::Array2<f64>) -> ndarray::Array2<f64> {
use ndarray::{concatenate, Axis};
let mut mats: Vec<ndarray::Array2<f64>> = Vec::new();
$( mats.push($crate::__constraints_helper!($($ineq),*)(genes)); )?
$( mats.push($crate::__constraints_helper!( $($crate::__eq_helper!($eq)),* )(genes)); )?
$( mats.push({ let lb_mat = genes.mapv(|_| $lb); lb_mat - genes }); )?
$( mats.push({ let ub_mat = genes.mapv(|_| $ub); genes - ub_mat }); )?
if mats.is_empty() {
ndarray::Array2::zeros((genes.nrows(), 0))
} else {
let views: Vec<_> = mats.iter().map(|m| m.view()).collect();
concatenate(Axis(1), &views)
.expect("Failed to concatenate constraints along axis 1")
}
}
$( fn lower_bound(&self) -> Option<f64> { Some($lb) } )?
$( fn upper_bound(&self) -> Option<f64> { Some($ub) } )?
}
};
}
#[cfg(test)]
mod tests {
use ndarray::{Array1, Array2, Axis, array};
use crate::ConstraintsFn;
fn g1(genes: &Array2<f64>) -> Array1<f64> {
genes.map_axis(Axis(1), |row| row.sum() - 1.0)
}
fn g2(genes: &Array2<f64>) -> Array1<f64> {
genes.map_axis(Axis(1), |row| row.dot(&row) - 1.0)
}
fn g3(genes: &Array2<f64>) -> Array1<f64> {
genes.map_axis(Axis(1), |row| row[0] - row[1])
}
#[test]
fn inequalities_only() {
let genes = array![[0.0, 0.0], [1.0, 1.0]]; impl_constraints_fn!(MyConstr, ineq = [g1, g2]); let res = MyConstr.call(&genes);
let expect = array![[-1.0, -1.0], [1.0, 1.0]];
assert_eq!(res.shape(), &[2, 2]);
assert_eq!(res, expect);
}
#[test]
fn mixed_ineq_and_eq() {
const EPS: f64 = 1e-6;
let genes = array![
[0.5, 0.5], [2.0, 1.0], ];
impl_constraints_fn!(MyConstr, ineq = [g1, g2], eq = [g3]);
let res = MyConstr.call(&genes);
let mut exp = Array2::<f64>::zeros((2, 3));
exp[[0, 0]] = 0.0; exp[[0, 1]] = -0.5; exp[[0, 2]] = -EPS;
exp[[1, 0]] = 2.0; exp[[1, 1]] = 4.0; exp[[1, 2]] = 1.0 - EPS;
assert_eq!(res, exp);
}
#[test]
fn equalities_only() {
let genes = array![[2.0, 2.0], [0.5, 1.5]]; impl_constraints_fn!(EqOnly, eq = [g3]);
let res = EqOnly.call(&genes);
const EPS: f64 = 1e-6;
let mut exp = Array2::<f64>::zeros((2, 1));
exp[[0, 0]] = -EPS; exp[[1, 0]] = 1.0 - EPS;
assert_eq!(res, exp);
}
#[test]
fn lower_bound_only() {
let genes = array![[1.0, 3.0], [0.0, 2.0]];
impl_constraints_fn!(LowOnly, lower_bound = 2.0);
let res = LowOnly.call(&genes);
let mut exp = Array2::<f64>::zeros((2, 2));
exp[[0, 0]] = 2.0 - 1.0;
exp[[0, 1]] = 2.0 - 3.0;
exp[[1, 0]] = 2.0 - 0.0;
exp[[1, 1]] = 2.0 - 2.0;
assert_eq!(res, exp);
}
#[test]
fn upper_bound_only() {
let genes = array![[1.0, 3.0], [0.0, 2.0]];
impl_constraints_fn!(UpOnly, upper_bound = 3.0);
let res = UpOnly.call(&genes);
let mut exp = Array2::<f64>::zeros((2, 2));
exp[[0, 0]] = 1.0 - 3.0;
exp[[0, 1]] = 3.0 - 3.0;
exp[[1, 0]] = 0.0 - 3.0;
exp[[1, 1]] = 2.0 - 3.0;
assert_eq!(res, exp);
}
#[test]
fn all_constraints_combined() {
const EPS: f64 = 1e-6;
let genes = array![[0.0, 1.0], [2.0, 1.0]];
impl_constraints_fn!(
AllC,
ineq = [g1],
eq = [g3],
lower_bound = 1.0,
upper_bound = 2.0
);
let res = AllC.call(&genes); let mut exp = Array2::<f64>::zeros((2, 6));
exp[[0, 0]] = 0.0; exp[[0, 1]] = 1.0 - EPS; exp[[0, 2]] = 1.0 - 0.0; exp[[0, 3]] = 1.0 - 1.0; exp[[0, 4]] = 0.0 - 2.0; exp[[0, 5]] = 1.0 - 2.0;
exp[[1, 0]] = 2.0; exp[[1, 1]] = 1.0 - EPS; exp[[1, 2]] = 1.0 - 2.0; exp[[1, 3]] = 1.0 - 1.0; exp[[1, 4]] = 2.0 - 2.0; exp[[1, 5]] = 1.0 - 2.0; assert_eq!(res, exp);
assert_eq!(res, exp);
}
}