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

// Copyright 2018 Paul Scott
//
// Licensed under the Apache License, Version 2.0 <LICENSE-APACHE or
// https://www.apache.org/licenses/LICENSE-2.0> or the MIT license
// <LICENSE-MIT or https://opensource.org/licenses/MIT>, at your
// option. This file may not be copied, modified, or distributed
// except according to those terms.

//! Model interface for solvers.

use crate::expr::{Expression, Par, Retrieve, Store, Var, ID};

//pub enum VarType {
//    Continuous,
//    Integer,
//    Binary,
//}

/// Constraint identifier.
#[derive(Debug, Clone, Copy)]
pub struct Con(pub ID);

/// Interface for a mathematical program with continuous variables.
///
/// # Panics
///
/// Expect a panic, either our of bounds or otherwise, if a variable or
/// parameter is used in this model that was not directly returned by the
/// `add_var` or `add_par` methods. Check individual implementations of trait
/// for details of when these panics could occur.
pub trait Model {
    /// Add variable to model with lower / upper bounds and initial value.
    fn add_var(&mut self, lb: f64, ub: f64, init: f64) -> Var;
    /// Add parameter to model with starting value.
    fn add_par(&mut self, val: f64) -> Par;
    /// Add a constraint to the model with lower and upper bounds.
    ///
    /// To have no lower / upper bounds set them to `std::f64::NEG_INFINITY` /
    /// `std::f64::INFINITY` respectively.
    fn add_con<E: Into<Expression>>(&mut self, expr: E, lb: f64, ub: f64) -> Con;
    /// Set objective of model.
    fn set_obj<E: Into<Expression>>(&mut self, expr: E);
    /// Change a parameter's value.
    fn set_par(&mut self, par: Par, val: f64);
    /// Change the initial value of a variable.
    fn set_init(&mut self, var: Var, init: f64);
    /// Solve the model.
    fn solve(&mut self) -> (SolutionStatus, Option<Solution>);
    /// Solve the model using a previous solution as a warm start.
    fn warm_solve(&mut self, sol: Solution) -> (SolutionStatus, Option<Solution>);
}

// Not used yet
//pub trait MIModel {
//    fn add_ivar(&mut self, lb: f64, ub: f64) -> Var;
//    fn add_bvar(&mut self, lb: f64, ub: f64) -> Var;
//}

/// Status of the solution.
#[derive(PartialEq, Debug)]
pub enum SolutionStatus {
    /// Problem successfully solved.
    Solved,
    /// Problem appears to be infeasible.
    Infeasible,
    /// An error occurred during the solve process.
    Error,
    /// Other currrently unhandled solver status returned.
    Other,
}

/// Data for a valid solution.
#[derive(Default)]
pub struct Solution {
    /// Objective value.
    pub obj_val: f64,
    /// Store of variable and parameter values.
    pub store: Store,
    /// Constraint Lagrange / KKT multipliers.
    pub con_mult: Vec<f64>,
    /// Variable lower bound KKT multipliers.
    pub var_lb_mult: Vec<f64>,
    /// Variable upper bound KKT multipliers.
    pub var_ub_mult: Vec<f64>,
}

impl Solution {
    pub fn new() -> Self {
        Self::default()
    }

    /// Calculate the value of an expression using the solution.
    pub fn value(&self, expr: &Expression) -> f64 {
        match expr {
            Expression::ExprFix(e) => (e.f)(&self.store.vars, &self.store.pars),
            Expression::ExprFixSum(es) => {
                let mut val = 0.0;
                for e in es {
                    val += (e.f)(&self.store.vars, &self.store.pars);
                }
                val
            }
            Expression::ExprDyn(e) => {
                let mut ns = Vec::new();
                e.expr.eval(&self.store, &mut ns)
            }
            Expression::ExprDynSum(es) => {
                let mut ns = Vec::new();
                let mut val = 0.0;
                for e in es {
                    val += e.expr.eval(&self.store, &mut ns);
                }
                val
            }
        }
    }

    /// Get the value of variable for solution.
    pub fn var(&self, v: Var) -> f64 {
        self.store.var(v)
    }

    /// Get the constraint KKT / Lagrange multiplier.
    pub fn con_mult(&self, Con(cid): Con) -> f64 {
        self.con_mult[cid]
    }

    // Could write versions that take Expr, and try and match ops[0] to Var
    /// Get the variable lower bound constraint KKT / Lagrange multiplier.
    pub fn var_lb_mult(&self, Var(vid): Var) -> f64 {
        self.var_lb_mult[vid]
    }

    /// Get the variable upper bound constraint KKT / Lagrange multiplier.
    pub fn var_ub_mult(&self, Var(vid): Var) -> f64 {
        self.var_ub_mult[vid]
    }
}