RsBundle  Changes On Branch modifyprimals

target/
*.log
*.cpxlog
Cargo.lock
instances/

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

[dependencies]
itertools = "^0.8"
libc = "^0.2.6"
log = "^0.4"
c_str_macro = "^1.0"
cplex-sys = "^0.6"

[dev-dependencies]
env_logger = "^0.7"

 * 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();
}

// Copyright (c) 2016, 2017, 2019 Frank Fischer <frank-fischer@shadow-soft.de>
//
// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// 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/>
//

//! 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.
 *

}

/// Problem update information.
///
/// The solver calls the `update` method of the problem regularly.
/// This method can modify the problem by adding (or removing)
/// variables. The possible updates are encoded in this type.

///
/// - `P` is the type of primals
/// - `E` is the error type
pub enum Update<P, E> {
    /// Add a variable with bounds.
    ///
    /// The initial value of the variable will be the feasible value
    /// closest to 0.
    AddVariable { lower: Real, upper: Real },
    /// Add a variable with bounds and initial value.
    AddVariableValue { lower: Real, upper: Real, value: Real },
    /// Change the current value of a variable. The bounds remain
    /// unchanged.
    MoveVariable { index: usize, value: Real },
    /// Update primal variables of a subproblem.
    ModifyPrimals {
        subproblem: usize,
        modify: Box<dyn Fn(&mut P) -> Result<(), E>>,
    },
}

/**
 * Trait for implementing a first-order problem description.
 *
 */
pub trait FirstOrderProblem {

        let mut primals = primals;
        primals.pop().unwrap().1
    }

    /// Return updates of the problem.
    ///
    /// The default implementation returns no updates.
    fn update(
        &mut self,
        _state: &UpdateState<Self::Primal>,
    ) -> Result<Vec<Update<Self::Primal, Self::Err>>, Self::Err> {
        Ok(vec![])
    }

    /// Return new components for a subgradient.
    ///
    /// The components are typically generated by some primal information. The
    /// corresponding primal along with its subproblem index is passed as a
    /// parameter.
    ///
    /// The default implementation fails because it should never be
    /// called.
    fn extend_subgradient(
        &mut self,
        _i: usize,
        _primal: &Self::Primal,
        _vars: &[usize],
    ) -> Result<Vec<Real>, Self::Err> {
        unimplemented!()
    }
}

//!
//! 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;

// Copyright (c) 2016, 2017, 2019 Frank Fischer <frank-fischer@shadow-soft.de>
//
// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// 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)]

            Solver(err) => write!(fmt, "Solver error: {}", err),
            Custom(err) => err.fmt(fmt),
        }
    }
}

impl Error for MasterProblemError {
    fn cause(&self) -> Option<&dyn Error> {
        use self::MasterProblemError::*;
        match self {
            SubgradientExtension(err) | Solver(err) | Custom(err) => Some(err.as_ref()),
            NoMinorants => None,
        }
    }
}

/// Result type of master problems.
pub type Result<T> = result::Result<T, MasterProblemError>;

/// Callback for subgradient extensions.
pub type SubgradientExtension<'a, I> = dyn FnMut(usize, I, &[usize]) -> result::Result<DVector, Box<dyn Error>> + 'a;

pub trait MasterProblem {
    /// Unique index for a minorant.
    type MinorantIndex: Copy + Eq;

    /// Set the number of subproblems.
    fn set_num_subproblems(&mut self, n: usize) -> Result<()>;

    /// Set the maximal number of inner iterations.
    fn set_max_updates(&mut self, max_updates: usize) -> Result<()>;

    /// Return the current number of inner iterations.
    fn cnt_updates(&self) -> usize;

    /// Add or move some variables with bounds.
    ///
    /// If an index is specified, existing variables are moved,
    /// otherwise new variables are generated.
    fn add_vars(
        &mut self,
        bounds: &[(Option<usize>, Real, Real)],
        extend_subgradient: &mut SubgradientExtension<Self::MinorantIndex>,
    ) -> Result<()>;

    /// Add a new minorant to the model.
    ///
    /// The function returns a unique (among all minorants of all
    /// subproblems) index of the minorant. This index must remain
    /// valid until the minorant is aggregated.

// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer <frank-fischer@shadow-soft.de>
//
// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// 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::UnconstrainedMasterProblem;
use crate::master::{MasterProblem, SubgradientExtension};
use crate::{DVector, Minorant, Real};

use super::Result;

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


/**
 * Turn unconstrained master problem into box-constrained one.
 *
 * This master problem adds box constraints to an unconstrainted
 * master problem implementation. The box constraints are enforced by
 * an additional outer optimization loop.

                        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$.
     */
    fn get_norm_subg2(&self) -> Real {
        self.eta.dot(self.master.dualopt())

    }
}

impl<M: UnconstrainedMasterProblem> MasterProblem for BoxedMasterProblem<M> {
    type MinorantIndex = M::MinorantIndex;

    fn set_num_subproblems(&mut self, n: usize) -> Result<()> {

    fn set_weight(&mut self, weight: Real) -> Result<()> {
        self.master.set_weight(weight)
    }

    fn add_vars(
        &mut self,
        bounds: &[(Option<usize>, Real, Real)],
        extend_subgradient: &mut SubgradientExtension<Self::MinorantIndex>,
    ) -> Result<()> {
        if !bounds.is_empty() {
            for (index, l, u) in bounds.iter().filter_map(|v| v.0.map(|i| (i, v.1, v.2))) {
                self.lb[index] = l;
                self.ub[index] = u;
            }
            self.lb.extend(bounds.iter().filter(|v| v.0.is_none()).map(|x| x.1));
            self.ub.extend(bounds.iter().filter(|v| v.0.is_none()).map(|x| x.2));
            self.eta.resize(self.lb.len(), 0.0);
            self.need_new_candidate = true;
            let nnew = bounds.iter().filter(|v| v.0.is_none()).count();
            let changed = bounds.iter().filter_map(|v| v.0).collect::<Vec<_>>();
            self.master.add_vars(nnew, &changed, extend_subgradient)
        } else {
            Ok(())
        }
    }

    #[allow(clippy::cognitive_complexity)]
    fn solve(&mut self, center_value: Real) -> Result<()> {
        debug!("Solve Master");
        debug!("  lb      ={}", self.lb);
        debug!("  ub      ={}", self.ub);

        if self.need_new_candidate {
            self.compute_candidate();

// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer <frank-fischer@shadow-soft.de>
//
// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// 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/>
//

//! Master problem implementation using CPLEX.

#![allow(unused_unsafe)]

use crate::master::{Error as MasterProblemError, SubgradientExtension, 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, warn};


use std;

use std::f64::{self, NEG_INFINITY};
use std::iter::{once, repeat};
use std::os::raw::{c_char, c_int};
use std::ptr;


impl From<cpx::CplexError> for MasterProblemError {
    fn from(err: cpx::CplexError) -> MasterProblemError {
        MasterProblemError::Solver(err.into())
    }
}

        }
    }

    fn add_vars(
        &mut self,
        nnew: usize,
        changed: &[usize],
        extend_subgradient: &mut SubgradientExtension<Self::MinorantIndex>,
    ) -> Result<()> {
        debug_assert!(!self.minorants[0].is_empty());
        let noldvars = self.minorants[0][0].linear.len();
        let nnewvars = noldvars + nnew;

        let mut changedvars = vec![];
        changedvars.extend_from_slice(changed);
        changedvars.extend(noldvars..nnewvars);
        for (fidx, mins) in self.minorants.iter_mut().enumerate() {
            if !mins.is_empty() {
                for (i, m) in mins.iter_mut().enumerate() {
                    let new_subg = extend_subgradient(fidx, self.min2index[fidx][i], &changedvars)
                        .map_err(MasterProblemError::SubgradientExtension)?;
                    for (&j, &g) in changed.iter().zip(new_subg.iter()) {
                        m.linear[j] = g;
                    }
                    m.linear.extend_from_slice(&new_subg[changed.len()..]);
                }
            }
        }

        // 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

                nvars as c_int,
                inds.as_ptr(),
                c.as_ptr()
            ));
        }

        trycpx!(cpx::qpopt(cpx::env(), self.lp));
        let solstat = unsafe { cpx::getstat(cpx::env(), self.lp) };
        if solstat != cpx::Stat::Optimal.to_c() && solstat != cpx::Stat::NumBest.to_c() {
            warn!("Problem solving the master QP with Cplex (status: {:?})", solstat)
        }
        let mut sol = vec![0.0; nvars];
        trycpx!(cpx::getx(cpx::env(), self.lp, sol.as_mut_ptr(), 0, nvars as c_int - 1));

        let mut idx = 0;
        let mut mults = Vec::with_capacity(nvars);
        let mut mins = Vec::with_capacity(nvars);
        for fidx in 0..self.minorants.len() {

            sum_coeffs += self.opt_mults[fidx][self.index2min[i].1];
        }
        let aggr_coeffs = if sum_coeffs != 0.0 {
            mins.iter()
                .map(|&i| self.opt_mults[fidx][self.index2min[i].1] / sum_coeffs)
                .collect::<DVector>()
        } else {
            // All coefficients are zero, we can simply remove all but one of the minorants.
            once(1.0).chain(repeat(0.0)).take(mins.len()).collect()
        };

        // compute aggregated diagonal term
        let mut aggr_diag = 0.0;
        for (idx_i, &i) in mins.iter().enumerate() {
            for (idx_j, &j) in mins.iter().enumerate() {
                aggr_diag += aggr_coeffs[idx_i] * aggr_coeffs[idx_j] * self.qterm[i][j];

// Copyright (c) 2016, 2017, 2019 Frank Fischer <frank-fischer@shadow-soft.de>
//
// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// 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, SubgradientExtension, 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)]

        Ok(self.minorants.len() - 1)
    }

    fn add_vars(
        &mut self,
        nnew: usize,
        changed: &[usize],
        extend_subgradient: &mut SubgradientExtension<Self::MinorantIndex>,
    ) -> Result<()> {
        if !self.minorants.is_empty() {
            let noldvars = self.minorants[0].linear.len();
            let mut changedvars = vec![];
            changedvars.extend_from_slice(changed);
            changedvars.extend(noldvars..noldvars + nnew);
            for (i, m) in self.minorants.iter_mut().enumerate() {

// Copyright (c) 2016, 2017, 2019 Frank Fischer <frank-fischer@shadow-soft.de>
//
// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// 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/>
//


//! Bundle master problem solver.
//!
//! This module contains solvers for the bundle master problem, i.e.
//! for solving convex optimization problems of the form
//!
//! \\[ \min \left\\{ \hat{f}(d) + \frac{w}{2} \\|d\\|\^2 \colon d \in [l,u] \right\\}, \\]
//!

//! * changing the weight parameter $w$,
//! * modifying $\hat{f}$ by adding or removing linear functions $\ell_i$,
//! * moving the center of the linear functions $\ell_i$ (and the
//!   bounds), i.e. replacing $\hat{f}$ by $d \mapsto \hat{f}(d -
//!   \hat{d})$ for some given $\hat{d} \in \mathbb{R}\^n$.

mod base;
pub use self::base::{MasterProblem, MasterProblemError as Error, Result, SubgradientExtension};

mod boxed;
pub use self::boxed::BoxedMasterProblem;

mod unconstrained;
pub use self::unconstrained::UnconstrainedMasterProblem;

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

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

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

// Copyright (c) 2016, 2017, 2019 Frank Fischer <frank-fischer@shadow-soft.de>
//
// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// 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, SubgradientExtension};




/**
 * Trait for master problems without box constraints.
 *
 * Implementors of this trait are supposed to solve quadratic
 * optimization problems of the form
 *

    /// The variables in `changed` have been changed, so the subgradient
    /// information must be updated. Furthermore, `nnew` new variables
    /// are added.
    fn add_vars(
        &mut self,
        nnew: usize,
        changed: &[usize],
        extend_subgradient: &mut SubgradientExtension<Self::MinorantIndex>,
    ) -> Result<()>;

    /// Solve the master problem.
    fn solve(&mut self, eta: &DVector, fbound: Real, augbound: Real, relprec: Real) -> Result<()>;

    /// Return the current dual optimal solution.
    fn dualopt(&self) -> &DVector;

// 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;

// Copyright (c) 2016, 2017, 2019 Frank Fischer <frank-fischer@shadow-soft.de>
//
// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// 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;

        write!(fmt, "Format error: {}", self.msg)
    }
}

impl Error for MMCFFormatError {}

/// Result type of the MMCFProblem.
pub type Result<T> = result::Result<T, Box<dyn Error>>;

#[derive(Clone, Copy, Debug)]
struct ArcInfo {
    arc: usize,
    src: usize,
    snk: usize,
}

    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>();
                }

// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer <frank-fischer@shadow-soft.de>
//
// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// 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/>
//

#![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;

impl Solver {
    pub fn new(nnodes: usize) -> Result<Solver> {
        let mut status: c_int;
        let mut net = ptr::null_mut();

        unsafe {
            #[allow(clippy::never_loop)]
            loop {
                status = cpx::setlogfilename(cpx::env(), c_str!("mcf.cpxlog").as_ptr(), c_str!("w").as_ptr());
                if status != 0 {
                    break;
                }

                net = cpx::NETcreateprob(cpx::env(), &mut status, c_str!("mcf").as_ptr());

            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 {

// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer <frank-fischer@shadow-soft.de>
//
// This program is free software: you can redistribute it and/or
// modify it under the terms of the GNU General Public License as
// published by the Free Software Foundation, either version 3 of the
// License, or (at your option) any later version.
//
// This program is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// 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/>
//

//! 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;

    Parameter(ParameterError),
    /// The lower bound of a variable is larger than the upper bound.
    InvalidBounds { lower: Real, upper: Real },
    /// The value of a variable is outside its bounds.
    ViolatedBounds { lower: Real, upper: Real, value: Real },
    /// The variable index is out of bounds.
    InvalidVariable { index: usize, nvars: usize },
    /// A subproblem index is out of bounds.
    InvalidSubproblem { subproblem: usize, nsubs: usize },
    /// Iteration limit has been reached.
    IterationLimit { limit: usize },
}

impl<E: fmt::Display> fmt::Display for SolverError<E> {
    fn fmt(&self, fmt: &mut fmt::Formatter) -> Result<(), fmt::Error> {
        use self::SolverError::*;

                fmt,
                "Violated bounds, lower:{}, upper:{}, value:{}",
                lower, upper, value
            ),
            InvalidVariable { index, nvars } => {
                write!(fmt, "Variable index out of bounds, got:{} must be < {}", index, nvars)
            }
            InvalidSubproblem { subproblem, nsubs } => write!(
                fmt,
                "Subproblem index out of bounds, got:{} must be < {}",
                subproblem, nsubs
            ),
            IterationLimit { limit } => write!(fmt, "The iteration limit of {} has been reached.", limit),
        }
    }
}

impl<E: Error> Error for SolverError<E> {
    fn cause(&self) -> Option<&dyn Error> {
        match self {
            SolverError::Evaluation(err) => Some(err),
            SolverError::Update(err) => Some(err),
            SolverError::Master(err) => Some(err),
            _ => None,
        }
    }
}

impl<E> From<ParameterError> for SolverError<E> {
    fn from(err: ParameterError) -> SolverError<E> {
        SolverError::Parameter(err)
    }
}

impl<E> From<MasterProblemError> for SolverError<E> {
    fn from(err: MasterProblemError) -> SolverError<E> {
        SolverError::Master(err)
    }
}

/**
 * The current state of the bundle method.
 *
 * Captures the current state of the bundle method during the run of
 * the algorithm. This state is passed to certain callbacks like
 * Terminator or Weighter so that they can compute their result

    Descent,
    /// No step but the algorithm has been terminated.
    Term,
}

/// Information about a minorant.
#[derive(Debug, Clone)]
struct MinorantInfo {
    /// The minorant's index in the master problem
    index: usize,
    /// Current multiplier.
    multiplier: Real,


}

/// Information about the last iteration.
#[derive(Debug, Clone, Copy, PartialEq)]
pub enum IterationInfo {
    NewMinorantTooHigh { new: Real, old: Real },
    UpperBoundNullStep,
    ShallowCut,
}

/// State information for the update callback.
pub struct UpdateState<'a, Pr: 'a> {
    /// Current model minorants.
    minorants: &'a [Vec<MinorantInfo>],
    /// The primals.
    primals: &'a Vec<Option<Pr>>,
    /// The last step type.
    pub step: Step,
    /// Iteration information.
    pub iteration_info: &'a [IterationInfo],
    /// The current candidate. If the step was a descent step, this is
    /// the new center.
    pub nxt_y: &'a DVector,
    /// The center. IF the step was a descent step, this is the old
    /// center.
    pub cur_y: &'a DVector,
}

impl<'a, Pr: 'a> UpdateState<'a, Pr> {
    pub fn aggregated_primals(&self, subproblem: usize) -> Vec<(Real, &Pr)> {
        self.minorants[subproblem]
            .iter()
            .map(|m| (m.multiplier, self.primals[m.index].as_ref().unwrap()))
            .collect()
    }

    /// Return the last primal for a given subproblem.
    ///
    /// This is the last primal generated by the oracle.
    pub fn last_primal(&self, fidx: usize) -> Option<&Pr> {
        self.minorants[fidx].last().and_then(|m| self.primals[m.index].as_ref())
    }
}

/**
 * Implementation of a bundle method.
 */
pub struct Solver<P: FirstOrderProblem> {
    /// The first order problem description.
    problem: P,

    /// The solver parameter.
    pub params: SolverParams,

    /// Termination predicate.
    pub terminator: Box<dyn Terminator>,

    /// Weighter heuristic.
    pub weighter: Box<dyn Weighter>,

    /// Lower and upper bounds of all variables.
    bounds: Vec<(Real, Real)>,

    /// Current center of stability.
    cur_y: DVector,

     * Time when the solution process started.
     *
     * This is actually the time of the last call to `Solver::init`.
     */
    start_time: Instant,

    /// The master problem.
    master: Box<dyn MasterProblem<MinorantIndex = usize>>,

    /// The active minorant indices for each subproblem.
    minorants: Vec<Vec<MinorantInfo>>,

    /// The primals associated with each global minorant index.
    primals: Vec<Option<P::Primal>>,

    /// Accumulated information about the last iteration.
    iterinfos: Vec<IterationInfo>,
}

impl<P: FirstOrderProblem> Solver<P>
where

            nxt_mods: dvec![],
            new_cutval: 0.0,
            sgnorm: 0.0,
            expected_progress: 0.0,
            cnt_descent: 0,
            cnt_null: 0,
            start_time: Instant::now(),
            master: Box::new(BoxedMasterProblem::new(MinimalMaster::new()?)),


            minorants: vec![],
            primals: vec![],
            iterinfos: vec![],
        })
    }

    /// A new solver with default parameter.
    pub fn new(problem: P) -> Result<Solver<P>, SolverError<P::Err>> {
        Solver::new_params(problem, SolverParams::default())

        self.nxt_mods.init0(m);

        self.start_time = Instant::now();

        Ok(())
    }

    /// Solve the problem with at most 10_000 iterations.
    ///
    /// Use `solve_with_limit` for an explicit iteration limit.
    pub fn solve(&mut self) -> Result<(), SolverError<P::Err>> {
        const LIMIT: usize = 10_000;
        self.solve_with_limit(LIMIT)
    }

    /// Solve the problem with explicit iteration limit.
    pub fn solve_with_limit(&mut self, iter_limit: usize) -> Result<(), SolverError<P::Err>> {
        // First initialize the internal data structures.
        self.init()?;

        if self.solve_iter(iter_limit)? {
            Ok(())
        } else {
            Err(SolverError::IterationLimit { limit: iter_limit })
        }
    }

    /// Solve the problem but stop after `niter` iterations.
    ///
    /// The function returns `Ok(true)` if the termination criterion
    /// has been satisfied. Otherwise it returns `Ok(false)` or an

    ///
    /// Calling this function typically triggers the problem to
    /// separate new constraints depending on the current solution.
    fn update_problem(&mut self, term: Step) -> Result<bool, SolverError<P::Err>> {
        let updates = {
            let state = UpdateState {
                minorants: &self.minorants,
                primals: &self.primals,
                step: term,
                iteration_info: &self.iterinfos,
                // this is a dirty trick: when updating the center, we
                // simply swapped the `cur_*` fields with the `nxt_*`
                // fields
                cur_y: if term == Step::Descent {
                    &self.nxt_y

                    }
                    let (lower, upper) = self.bounds[index];
                    if value < lower || value > upper {
                        return Err(SolverError::ViolatedBounds { lower, upper, value });
                    }
                    newvars.push((Some(index), lower - value, upper - value, value));
                }
                Update::ModifyPrimals { subproblem, modify } => {
                    if subproblem >= self.minorants.len() {
                        return Err(SolverError::InvalidSubproblem {
                            subproblem: subproblem,
                            nsubs: self.minorants.len(),
                        });
                    }
                    for m in &mut self.minorants[subproblem] {
                        if let Some(ref mut p) = self.primals[m.index] {
                            if let Err(err) = modify(p) {
                                return Err(SolverError::Update(err));
                            }
                        }
                    }
                }
            }
        }

        if !newvars.is_empty() {
            let problem = &mut self.problem;
            let primals = &self.primals;
            self.master.add_vars(

                &newvars.iter().map(|v| (v.0, v.1, v.2)).collect::<Vec<_>>(),
                &mut |fidx, minidx, vars| {
                    problem
                        .extend_subgradient(fidx, primals[minidx].as_ref().unwrap(), vars)
                        .map(DVector)
                        .map_err(|e| e.into())
                },
            )?;
            // 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

    /// This function returns all currently used minorants $x_i$ along
    /// with their coefficients $\alpha_i$. The aggregated primal can
    /// be computed by combining the minorants $\bar{x} =
    /// \sum_{i=1}\^m \alpha_i x_i$.
    pub fn aggregated_primals(&self, subproblem: usize) -> Vec<(Real, &P::Primal)> {
        self.minorants[subproblem]
            .iter()
            .map(|m| (m.multiplier, self.primals[m.index].as_ref().unwrap()))
            .collect()
    }

    fn show_info(&self, step: Step) {
        let time = self.start_time.elapsed();
        info!(
            "{} {:0>2}:{:0>2}:{:0>2}.{:0>2} {:4} {:4} {:4}{:1}  {:9.4} {:9.4} \

     * information.
     */
    fn init_master(&mut self) -> Result<(), SolverError<P::Err>> {
        let m = self.problem.num_subproblems();

        self.master = if m == 1 && self.params.max_bundle_size == 2 {
            debug!("Use minimal master problem");
            Box::new(BoxedMasterProblem::new(MinimalMaster::new()?))


        } else {
            debug!("Use CPLEX master problem");
            Box::new(BoxedMasterProblem::new(CplexMaster::new()?))


        };

        let lb = self.problem.lower_bounds().map(DVector);
        let ub = self.problem.upper_bounds().map(DVector);

        if lb
            .as_ref()
            .map(|lb| lb.len() != self.problem.num_variables())
            .unwrap_or(false)
        {
            return Err(SolverError::Dimension);
        }
        if ub
            .as_ref()
            .map(|ub| ub.len() != self.problem.num_variables())
            .unwrap_or(false)
        {
            return Err(SolverError::Dimension);
        }

        self.master.set_num_subproblems(m)?;

        self.master.set_vars(self.problem.num_variables(), lb, ub)?;


        self.master.set_max_updates(self.params.max_updates)?;


        self.minorants = (0..m).map(|_| vec![]).collect();

        self.cur_val = 0.0;
        for i in 0..m {
            let result = self
                .problem
                .evaluate(i, &self.cur_y, INFINITY, 0.0)
                .map_err(SolverError::Evaluation)?;
            self.cur_vals[i] = result.objective();
            self.cur_val += self.cur_vals[i];

            let mut minorants = result.into_iter();
            if let Some((minorant, primal)) = minorants.next() {
                self.cur_mods[i] = minorant.constant;
                self.cur_mod += self.cur_mods[i];
                let minidx = self.master.add_minorant(i, minorant)?;
                self.minorants[i].push(MinorantInfo {
                    index: minidx,
                    multiplier: 0.0,

                });
                if minidx >= self.primals.len() {
                    self.primals.resize_with(minidx + 1, || None);
                }
                self.primals[minidx] = Some(primal);
            } else {
                return Err(SolverError::NoMinorant);
            }
        }

        self.cur_valid = true;

        // Solve the master problem once to compute the initial
        // subgradient.
        //
        // We could compute that subgradient directly by
        // adding up the initial minorants, but this would not include
        // the eta terms. However, this is a heuristic anyway because
        // we assume an initial weight of 1.0, which, in general, will
        // *not* be the initial weight for the first iteration.
        self.master.set_weight(1.0)?;
        self.master.solve(self.cur_val)?;
        self.sgnorm = self.master.get_dualoptnorm2().sqrt();

        // Compute the real initial weight.
        let state = current_state!(self, Step::Term);
        let new_weight = self.weighter.weight(&state, &self.params);
        self.master.set_weight(new_weight)?;

        debug!("Init master completed");

        Ok(())
    }

    /// Solve the model (i.e. master problem) to compute the next candidate.
    fn solve_model(&mut self) -> Result<(), SolverError<P::Err>> {
        self.master.solve(self.cur_val)?;
        self.nxt_d = self.master.get_primopt();
        self.nxt_y.add(&self.cur_y, &self.nxt_d);
        self.nxt_mod = self.master.get_primoptval();
        self.sgnorm = self.master.get_dualoptnorm2().sqrt();
        self.expected_progress = self.cur_val - self.nxt_mod;

        // update multiplier from master solution

        for i in 0..self.problem.num_subproblems() {
            let n = self.master.num_minorants(i);
            if n >= self.params.max_bundle_size {
                // aggregate minorants with smallest coefficients
                self.minorants[i].sort_by_key(|m| -((1e6 * m.multiplier) as isize));
                let aggr = self.minorants[i].split_off(self.params.max_bundle_size - 2);
                let aggr_sum = aggr.iter().map(|m| m.multiplier).sum();
                let (aggr_mins, aggr_primals): (Vec<_>, Vec<_>) = aggr
                    .into_iter()
                    .map(|m| (m.index, self.primals[m.index].take().unwrap()))
                    .unzip();
                let (aggr_min, aggr_coeffs) = self.master.aggregate(i, &aggr_mins)?;
                // append aggregated minorant
                self.minorants[i].push(MinorantInfo {
                    index: aggr_min,
                    multiplier: aggr_sum,
                });
                self.primals[aggr_min] = Some(
                    self.problem
                        .aggregate_primals(aggr_coeffs.into_iter().zip(aggr_primals.into_iter()).collect()),

                );
            }
        }
        Ok(())
    }

    /// Perform a descent step.
    fn descent_step(&mut self) -> Result<(), SolverError<P::Err>> {
        let new_weight = self.weighter.weight(&current_state!(self, Step::Descent), &self.params);
        self.master.set_weight(new_weight)?;
        self.cnt_descent += 1;
        swap(&mut self.cur_y, &mut self.nxt_y);
        swap(&mut self.cur_val, &mut self.nxt_val);
        swap(&mut self.cur_mod, &mut self.nxt_mod);
        swap(&mut self.cur_vals, &mut self.nxt_vals);
        swap(&mut self.cur_mods, &mut self.nxt_mods);
        self.master.move_center(1.0, &self.nxt_d);
        debug!("Descent Step");
        debug!("  dir ={}", self.nxt_d);
        debug!("  newy={}", self.cur_y);
        Ok(())
    }

    /// Perform a null step.
    fn null_step(&mut self) -> Result<(), SolverError<P::Err>> {
        let new_weight = self.weighter.weight(&current_state!(self, Step::Null), &self.params);
        self.master.set_weight(new_weight)?;
        self.cnt_null += 1;
        debug!("Null Step");
        Ok(())
    }

    /// Perform one bundle iteration.
    #[allow(clippy::collapsible_if)]
    pub fn step(&mut self) -> Result<Step, SolverError<P::Err>> {
        self.iterinfos.clear();

        if !self.cur_valid {
            // current point needs new evaluation
            self.init_master()?;
        }

            nxt_lb += fun_lb;
            nxt_ub += fun_ub;
            self.nxt_vals[fidx] = fun_ub;

            // move center of minorant to cur_y
            nxt_minorant.move_center(-1.0, &self.nxt_d);
            self.new_cutval += nxt_minorant.constant;
            let minidx = self.master.add_minorant(fidx, nxt_minorant)?;
            self.minorants[fidx].push(MinorantInfo {
                index: minidx,



                multiplier: 0.0,

            });
            if minidx >= self.primals.len() {
                self.primals.resize_with(minidx + 1, || None);
            }
            self.primals[minidx] = Some(nxt_primal);
        }

        if self.new_cutval > self.cur_val + 1e-3 {
            warn!(
                "New minorant has higher value in center new:{} old:{}",
                self.new_cutval, self.cur_val
            );

//
// 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, ")")
            }
        }
    }
}

    /// the smaller vector are assumed to be 0.0.
    pub fn add_scaled_begin(&mut self, alpha: Real, y: &DVector) {
        for (x, y) in self.iter_mut().zip(y.iter()) {
            *x += alpha * y;
        }
        let n = self.len();
        if n < y.len() {
            self.extend(y[n..].iter().map(|y| alpha * y));
        }
        // if self.len() < y.len() {
        //     self.resize(y.len(), 0.0);
        // }
        // for i in 0..y.len() {
        //     unsafe { *self.get_unchecked_mut(i) += alpha * *y.get_unchecked(i) };
        // }

                    unsafe { *v.get_unchecked_mut(i) = x };
                }
                DVector(v)
            }
        }
    }
}

#[test]
fn test_add_scaled_begin() {
    let mut x = dvec![1.0; 5];
    let y = dvec![2.0; 7];
    x.add_scaled_begin(3.0, &y);
    assert_eq!(x, dvec![7.0, 7.0, 7.0, 7.0, 7.0, 6.0, 6.0]);
}