RsBundle  Check-in [51732172c2]

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"

        _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!()
    }
}

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

        }
    }
}

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

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

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

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

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

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

        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;

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

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

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

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

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

}

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

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

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

    /// 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]);
}