RsBundle  Check-in [e1f5bd7920]

[package]
name = "bundle"
version = "0.6.0-dev"
edition = "2018"
authors = ["Frank Fischer <frank-fischer@shadow-soft.de>"]

[dependencies]
itertools = "^0.7.2"
libc = "^0.2.6"
log = "^0.4"
c_str_macro = "^1.0"

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see  <http://www.gnu.org/licenses/>
 */

use bundle;
use env_logger;

use log::info;

use bundle::mcf;
use bundle::{FirstOrderProblem, Solver, SolverParams, StandardTerminator};

use std::env;

fn main() {
    env_logger::init();

    let mut args = env::args();
    let program = args.next().unwrap();

            mmcf,
            SolverParams {
                max_bundle_size: 25,
                min_weight: 1e-3,
                max_weight: 100.0,
                ..Default::default()
            },
        )
        .unwrap();
        solver.terminator = Box::new(StandardTerminator {
            termination_precision: 1e-6,
        });
        solver.solve().unwrap();

        let costs: f64 = (0..solver.problem().num_subproblems())
            .map(|i| {
                let primals = solver.aggregated_primals(i);
                let aggr_primals = solver.problem().aggregate_primals_ref(&primals);
                solver.problem().get_primal_costs(i, &aggr_primals)
            })
            .sum();
        info!("Primal costs: {}", costs);
    } else {
        panic!("Usage: {} FILENAME", program);
    }
}

 *
 * You should have received a copy of the GNU General Public License
 * along with this program.  If not, see  <http://www.gnu.org/licenses/>
 */

use std::error::Error;


use bundle::{self, dvec};
use env_logger;

use log::debug;

use bundle::{DVector, FirstOrderProblem, Minorant, Real, SimpleEvaluation, Solver, SolverParams};

struct QuadraticProblem {
    a: [[Real; 2]; 2],
    b: [Real; 2],
    c: Real,

    let mut solver = Solver::new_params(
        f,
        SolverParams {
            min_weight: 1.0,
            max_weight: 1.0,
            ..Default::default()
        },
    )
    .unwrap();
    solver.solve().unwrap();
}

//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

//! Problem description of a first-order convex optimization problem.

use crate::solver::UpdateState;
use crate::{Minorant, Real};

use std::result::Result;
use std::vec::IntoIter;

/**
 * Trait for results of an evaluation.
 *

//!
//! The procedure is described in
//!
//! > Helmberg, C. and Kiwiel, K.C. (2002): A spectral bundle method
//! > with bounds, Math. Programming A 93, 173--194
//!

use crate::Real;
use crate::{BundleState, SolverParams, Step, Weighter};

use log::debug;

use std::cmp::{max, min};
use std::f64::NEG_INFINITY;

const FACTOR: Real = 2.0;

/**

        if state.step == Step::Term {
            self.eps_weight = 1e30;
            self.iter = 0;
            return if state.cur_y.len() == 0 || state.sgnorm < state.cur_y.len() as Real * 1e-10 {
                1.0
            } else {
                state.sgnorm.max(1e-4)
            }
            .max(params.min_weight)
            .min(params.max_weight);
        }

        let cur_nxt = state.cur_val - state.nxt_val;
        let cur_mod = state.cur_val - state.nxt_mod;
        let w = 2.0 * state.weight * (1.0 - cur_nxt / cur_mod);

        debug!("  cur_nxt={} cur_mod={} w={}", cur_nxt, cur_mod, w);

        if state.step == Step::Null {
            let sgnorm = state.sgnorm;
            let lin_err = state.cur_val - state.new_cutval;
            self.eps_weight = self.eps_weight.min(sgnorm + cur_mod - sgnorm * sgnorm / state.weight);
            let new_weight = if self.iter < -3 && lin_err > self.eps_weight.max(FACTOR * cur_mod) {
                w
            } else {
                state.weight
            }
            .min(FACTOR * state.weight)
            .min(params.max_weight);
            if new_weight > state.weight {
                self.iter = -1
            } else {
                self.iter = min(self.iter - 1, -1);
            }

            self.model_max = self.model_max.max(state.nxt_mod);
            let new_weight = if self.iter > 0 && cur_nxt > self.m_r * cur_mod {
                w
            } else if self.iter > 3 || state.nxt_val < self.model_max {
                state.weight / 2.0
            } else {
                state.weight
            }
            .max(state.weight / FACTOR)
            .max(params.min_weight);
            self.eps_weight = self.eps_weight.max(2.0 * cur_mod);
            if new_weight < state.weight {
                self.iter = 1;
                self.model_max = NEG_INFINITY;
            } else {
                self.iter = max(self.iter + 1, 1);

//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

//! Proximal bundle method implementation.








/// Type used for real numbers throughout the library.
pub type Real = f64;

#[macro_export]
macro_rules! dvec {
    ( $ elem : expr ; $ n : expr ) => { DVector(vec![$elem; $n]) };
    ( $ ( $ x : expr ) , * ) => { DVector(vec![$($x),*]) };
    ( $ ( $ x : expr , ) * ) => { DVector(vec![$($x,)*]) };
}




pub mod vector;
pub use crate::vector::{DVector, Vector};

pub mod minorant;
pub use crate::minorant::Minorant;

pub mod firstorderproblem;
pub use crate::firstorderproblem::{Evaluation, FirstOrderProblem, SimpleEvaluation, Update};

pub mod solver;
pub use crate::solver::{
    BundleState, IterationInfo, Solver, SolverParams, StandardTerminator, Step, Terminator, UpdateState, Weighter,
};

mod hkweighter;
pub use crate::hkweighter::HKWeighter;

mod master;

pub mod mcf;

// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

use crate::{DVector, Minorant, Real};

use std::error::Error;
use std::fmt;
use std::result;

/// Error type for master problems.
#[derive(Debug)]

// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

use crate::master::MasterProblem;
use crate::master::UnconstrainedMasterProblem;
use crate::{DVector, Minorant, Real};

use super::Result;

use itertools::multizip;
use log::debug;
use std::error::Error;
use std::f64::{EPSILON, INFINITY, NEG_INFINITY};
use std::result;

/**
 * Turn unconstrained master problem into box-constrained one.
 *

                        0.0
                    }
                } else if ub < INFINITY {
                    ub * eta
                } else {
                    0.0
                }
            })
            .sum()
    }

    /**
     * Return $\\|G \alpha - \eta\\|_2\^2$.
     *
     * This is the norm-square of the dual optimal solution including
     * the current box-multipliers $\eta$.

// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

//! Master problem implementation using CPLEX.

#![allow(unused_unsafe)]

use crate::master::{Error as MasterProblemError, UnconstrainedMasterProblem};
use crate::{DVector, Minorant, Real};

use super::Result;

use c_str_macro::c_str;
use cplex_sys as cpx;
use cplex_sys::trycpx;
use log::debug;


use std;
use std::error::Error;
use std::f64::{self, NEG_INFINITY};
use std::os::raw::{c_char, c_int};
use std::ptr;
use std::result;

        // self.force_update = true;

        Ok(())
    }

    fn solve(&mut self, eta: &DVector, _fbound: Real, _augbound: Real, _relprec: Real) -> Result<()> {
        if self.force_update || !self.updateinds.is_empty() {
            self.init_qp()?;
        }

        let nvars = unsafe { cpx::getnumcols(cpx::env(), self.lp) as usize };
        if nvars == 0 {
            return Err(MasterProblemError::NoMinorants);
        }
        // update linear costs

// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

use crate::master::{Error as MasterProblemError, UnconstrainedMasterProblem};
use crate::{DVector, Minorant, Real};

use super::Result;

use log::debug;

use std::error::Error;
use std::f64::NEG_INFINITY;
use std::fmt;
use std::result;

/// Minimal master problem error.
#[derive(Debug)]

pub mod minimal;
pub use self::minimal::MinimalMaster;

// pub mod grb;
// pub use master::grb::GurobiMaster;

mod cpx;
pub use crate::master::cpx::CplexMaster;

// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

use crate::{DVector, Minorant, Real};

use super::Result;

use std::error::Error;
use std::result;

/**

// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

//! Multi-commodity min-cost-flow subproblems.

mod solver;
pub use crate::mcf::solver::Solver;

mod problem;
pub use crate::mcf::problem::MMCFProblem;

// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// General Public License for more details.
//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

use crate::mcf;
use crate::{DVector, FirstOrderProblem, Minorant, Real, SimpleEvaluation};

use log::debug;

use std::error::Error;
use std::f64::INFINITY;
use std::fmt;
use std::fs::File;
use std::io::Read;
use std::result;

    c: Vec<DVector>,
}

impl MMCFProblem {
    pub fn read_mnetgen(basename: &str) -> Result<MMCFProblem> {
        let mut buffer = String::new();
        {
            let mut f = File::open(&format!("{}.nod", basename))?;
            f.read_to_string(&mut buffer)?;
        }
        let fnod = buffer
            .split_whitespace()
            .map(|x| x.parse::<usize>().unwrap())
            .collect::<Vec<_>>();

        if fnod.len() != 4 {
            return Err(MMCFFormatError {
                msg: format!("Expected 4 numbers in {}.nod, but got {}", basename, fnod.len()),
            }
            .into());
        }

        let ncom = fnod[0];
        let nnodes = fnod[1];
        let narcs = fnod[2];
        let ncaps = fnod[3];

        // read nodes
        let mut nets = Vec::with_capacity(ncom);
        for _ in 0..ncom {
            nets.push(mcf::Solver::new(nnodes)?)
        }
        {
            let mut f = File::open(&format!("{}.sup", basename))?;
            buffer.clear();
            f.read_to_string(&mut buffer)?;
        }
        for line in buffer.lines() {
            let mut data = line.split_whitespace();
            let node = data.next().unwrap().parse::<usize>()?;
            let com = data.next().unwrap().parse::<usize>()?;
            let supply = data.next().unwrap().parse::<Real>()?;
            nets[com - 1].set_balance(node - 1, supply)?;
        }

        // read arcs
        let mut arcmap = vec![vec![]; ncom];
        let mut cbase = vec![dvec![]; ncom];

        // lhs nonzeros
        let mut lhsidx = vec![vec![vec![]; ncom]; ncaps];
        {
            let mut f = File::open(&format!("{}.arc", basename))?;
            buffer.clear();
            f.read_to_string(&mut buffer)?;
        }
        for line in buffer.lines() {
            let mut data = line.split_whitespace();
            let arc = data.next().unwrap().parse::<usize>()? - 1;
            let src = data.next().unwrap().parse::<usize>()? - 1;
            let snk = data.next().unwrap().parse::<usize>()? - 1;
            let com = data.next().unwrap().parse::<usize>()? - 1;
            let cost = data.next().unwrap().parse::<Real>()?;
            let cap = data.next().unwrap().parse::<Real>()?;
            let mt = data.next().unwrap().parse::<isize>()? - 1;
            assert!(
                arc < narcs,
                format!("Wrong arc number (got: {}, expected in 1..{})", arc + 1, narcs)
            );
            // set internal coeff
            let coeff = arcmap[com].len();
            arcmap[com].push(ArcInfo {
                arc: arc + 1,
                src: src + 1,
                snk: snk + 1,
            });
            // add arc
            nets[com].add_arc(src, snk, cost, if cap < 0.0 { INFINITY } else { cap })?;
            // set objective
            cbase[com].push(cost); // + 1e-6 * coeff
                                   // add to mutual capacity constraint
            if mt >= 0 {
                lhsidx[mt as usize][com].push(coeff);
            }
        }

        // read rhs of coupling constraints
        {
            let mut f = File::open(&format!("{}.mut", basename))?;
            buffer.clear();
            f.read_to_string(&mut buffer)?;
        }
        let mut rhs = dvec![0.0; ncaps];
        for line in buffer.lines() {
            let mut data = line.split_whitespace();
            let mt = data.next().unwrap().parse::<usize>()? - 1;
            let cap = data.next().unwrap().parse::<Real>()?;
            rhs[mt] = cap;
        }

        // set lhs
        let mut lhs = vec![vec![vec![]; ncom]; ncaps];
        for i in 0..ncaps {
            for fidx in 0..ncom {

        let mut aggr = primals[0]
            .1
            .iter()
            .map(|x| {
                let mut r = dvec![];
                r.scal(primals[0].0, x);
                r
            })
            .collect::<Vec<_>>();

        for &(alpha, primal) in &primals[1..] {
            for (j, x) in primal.iter().enumerate() {
                aggr[j].add_scaled(alpha, x);
            }
        }

                }
            }
        }

        debug!("y={:?}", y);
        for i in 0..self.nets.len() {
            debug!("c[{}]={}", i, self.c[i]);
            self.nets[i].set_objective(&self.c[i])?;
        }

        // solve subproblems
        for (i, net) in self.nets.iter_mut().enumerate() {
            net.solve()?;
            debug!("c[{}]={}", i, net.objective()?);
        }

        // compute minorant
        if self.multimodel {
            let objective;
            let mut subg;
            if fidx == 0 {
                subg = self.rhs.clone();
                objective = self.rhsval - self.nets[fidx].objective()?;
            } else {
                subg = dvec![0.0; self.rhs.len()];
                objective = -self.nets[fidx].objective()?;
            }

            let sol = self.nets[fidx].get_solution()?;
            for (i, lhs) in self.lhs.iter().enumerate() {
                subg[i] -= lhs[fidx].iter().map(|elem| elem.val * sol[elem.ind]).sum::<Real>();
            }

            Ok(SimpleEvaluation {
                objective,
                minorants: vec![(
                    Minorant {
                        constant: objective,
                        linear: subg,
                    },
                    vec![sol],
                )],
            })
        } else {
            let mut objective = self.rhsval;
            let mut sols = Vec::with_capacity(self.nets.len());
            for i in 0..self.nets.len() {
                objective -= self.nets[i].objective()?;
                sols.push(self.nets[i].get_solution()?);
            }

            let mut subg = self.rhs.clone();
            for (i, lhs) in self.lhs.iter().enumerate() {
                for (fidx, flhs) in lhs.iter().enumerate() {
                    subg[i] -= flhs.iter().map(|elem| elem.val * sols[fidx][elem.ind]).sum::<Real>();
                }

//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

#![allow(unused_unsafe)]

use crate::{DVector, Real};

use c_str_macro::c_str;
use cplex_sys as cpx;
use cplex_sys::trycpx;

use std;
use std::error::Error;
use std::ffi::CString;
use std::ptr;
use std::result;

            if status != 0 {
                let msg = CString::new(vec![0; cpx::MESSAGE_BUF_SIZE]).unwrap().into_raw();
                cpx::geterrorstring(cpx::env(), status, msg);
                cpx::NETfreeprob(cpx::env(), &mut net);
                return Err(cpx::CplexError {
                    code: status,
                    msg: CString::from_raw(msg).to_string_lossy().into_owned(),
                }
                .into());
            }
        }

        Ok(Solver { net })
    }

    pub fn num_nodes(&self) -> usize {

//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

//! A linear minorant.

use crate::{DVector, Real};

use std::fmt;

/**
 * A linear minorant of a convex function.
 *
 * A linear minorant of a convex function $f \colon \mathbb{R}\^n \to

    /// The linear term.
    pub linear: DVector,
}

impl fmt::Display for Minorant {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        write!(f, "{} + y * {}", self.constant, self.linear)?;
        Ok(())
    }
}

impl Default for Minorant {
    fn default() -> Minorant {
        Minorant {

//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

//! The main bundle method solver.

use crate::{DVector, Real};
use crate::{Evaluation, FirstOrderProblem, HKWeighter, Update};

use crate::master::{BoxedMasterProblem, Error as MasterProblemError, MasterProblem, UnconstrainedMasterProblem};
use crate::master::{CplexMaster, MinimalMaster};

use log::{debug, info, warn};

use std::error::Error;
use std::f64::{INFINITY, NEG_INFINITY};
use std::fmt;
use std::mem::swap;
use std::result::Result;
use std::time::Instant;

                    &newvars.iter().map(|v| (v.0, v.1, v.2)).collect::<Vec<_>>(),
                    &mut |fidx, minidx, vars| {
                        problem
                            .extend_subgradient(minorants[fidx][minidx].primal.as_ref().unwrap(), vars)
                            .map(DVector)
                            .map_err(|e| e.into())
                    },
                )
                .map_err(SolverError::Master)?;
            // modify moved variables
            for (index, val) in newvars.iter().filter_map(|v| v.0.map(|i| (i, v.3))) {
                self.cur_y[index] = val;
                self.nxt_y[index] = val;
                self.nxt_d[index] = 0.0;
            }
            // add new variables

//
// You should have received a copy of the GNU General Public License
// along with this program.  If not, see  <http://www.gnu.org/licenses/>
//

//! Finite-dimensional sparse and dense vectors.

use crate::Real;
use std::fmt;
use std::ops::{Deref, DerefMut};

// use std::cmp::min;
use std::iter::FromIterator;
use std::vec::IntoIter;

/// Type of dense vectors.
#[derive(Debug, Clone, PartialEq, Default)]
pub struct DVector(pub Vec<Real>);

    fn deref_mut(&mut self) -> &mut Vec<Real> {
        &mut self.0
    }
}

impl fmt::Display for DVector {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        write!(f, "(")?;
        for (i, x) in self.iter().enumerate() {
            if i > 0 {
                write!(f, ", ")?;
            }
            write!(f, "{}", x)?
        }
        write!(f, ")")?;
        Ok(())
    }
}

impl FromIterator<Real> for DVector {
    fn from_iter<I: IntoIterator<Item = Real>>(iter: I) -> Self {
        DVector(Vec::from_iter(iter))

impl fmt::Display for Vector {
    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
        match *self {
            Vector::Dense(ref v) => write!(f, "{}", v),
            Vector::Sparse { size, ref elems } => {
                let mut it = elems.iter();
                write!(f, "{}:(", size)?;
                if let Some(&(i, x)) = it.next() {
                    write!(f, "{}:{}", i, x)?;
                    for &(i, x) in it {
                        write!(f, ", {}:{}", i, x)?;
                    }
                }
                write!(f, ")")
            }
        }
    }
}