RsBundle  Check-in [7b213efa22]

impl<P> Evaluation<P> for SimpleEvaluation<P> {
    fn objective(&self) -> Real {
        self.objective
    }
}


/**
 * Trait for implementing a first-order problem description.
 *
 */
pub trait FirstOrderProblem<'a> {
    /// Custom error type for evaluating this oracle.
    type Error : error::Error + 'static;

     *
     * Note that `nullstep_bound` and `relprec` are usually only
     * useful if there is only a `single` subproblem.
     */
    fn evaluate(&'a mut self, i : usize, y : &DVector,
                nullstep_bound : Real,
                relprec : Real) -> Result<Self::EvalResult, Self::Error>;

    /// Aggregate primal information.
    ///
    /// This function is called from the solver when minorants are
    /// aggregated. The problem can use this information to aggregate
    /// the corresponding primal information.
    ///
    /// - `coeffs` are the coefficient for the convex combination
    /// - `primals` are the corresponding primals to be aggregated
    ///
    /// The function must return the new aggregated primal.
    ///
    /// The default implementation does nothing and simply returns the
    /// last primal. This should work if the implementing problem does
    /// not provide primal information, e.g. if `Self::Primal = ()`.
    #[allow(unused_variables)]
    fn aggregate_primals(&mut self, coeffs: &[Real], primals: Vec<Self::Primal>) -> Self::Primal {
        let mut primals = primals;
        primals.pop().unwrap()
    }
}

        debug!("  dualopt={}", self.master.dualopt());
        debug!("  etaopt={}", self.eta);
        debug!("  primoptval={}", self.primoptval);

        Ok(())
    }

    fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(Self::MinorantIndex, DVector)> {
        self.master.aggregate(fidx, mins).map_err(|err| Error::Solver(Box::new(err)))
    }

    fn get_primopt(&self) -> DVector { self.primopt.clone() }

    fn get_primoptval(&self) -> Real { self.primoptval }

                this_val = this_val.max(m.eval(y));
            }
            result += this_val;
        }
        result
    }

    fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(usize, DVector)> {
        assert!(mins.len() > 0, "No minorants specified to be aggregated");

        if mins.len() == 1 { return Ok((mins[0], dvec![1.0])) }

        // scale coefficients
        let mut sum_coeffs = 0.0;
        for &i in mins {
            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 {
            dvec![0.0; mins.len()]
        };

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

        // compute aggregated off-diagonal terms
        let mut aggr_qterm = dvec![0.0; self.qterm.len()];
        for fidx_i in 0..self.minorants.len() {
            for idx_i in 0..self.minorants[fidx_i].len() {
                let i = self.min2index[fidx_i][idx_i];
                for (idx_j, &j) in mins.iter().enumerate() {

                    aggr_qterm[i] += aggr_coeffs[idx_j] * self.qterm[i][j];
                }
            }
        }

        // aggregate the minorants
        let mut aggr = Minorant::new();
        {
            let mut aggr_mins = Vec::with_capacity(mins.len());


            for &i in mins {
                let (min_fidx, min_idx) = self.index2min[i];
                debug_assert!(min_fidx == fidx);

                let m     = self.minorants[fidx].swap_remove(min_idx);
                let idx   = self.min2index[fidx].swap_remove(min_idx);
                self.opt_mults[fidx].swap_remove(min_idx);
                self.freeinds.push(idx);
                debug_assert!(idx == i);

                aggr_mins.push(m);


                // update index2min table for moved minorant
                if min_idx < self.minorants[fidx].len() {
                    self.index2min[self.min2index[fidx][min_idx]].1 = min_idx;
                }
            }
            aggr.combine_all(&aggr_coeffs, &aggr_mins);

                let i = self.min2index[fidx_i][idx_i];
                self.qterm[i][aggr_idx] = aggr_qterm[i];
                self.qterm[aggr_idx][i] = aggr_qterm[i];
            }
        }
        self.qterm[aggr_idx][aggr_idx] = aggr_diag;

        Ok((aggr_idx, aggr_coeffs))
    }

    fn move_center(&mut self, alpha: Real, d: &DVector) {
        for mins in self.minorants.iter_mut() {
            for m in mins.iter_mut() {
                m.move_center(alpha, d);
            }

    /// subproblems) index of the minorant. This index must remain
    /// valid until the minorant is aggregated.
    fn add_minorant(&mut self, fidx: usize, minorant: Minorant) -> Result<Self::MinorantIndex>;

    /// Solve the master problem.
    fn solve(&mut self, cur_value: Real) -> Result<()>;

    /// Aggregate the given minorants according to the current
    /// solution.
    ///
    /// The (indices of the) minorants to be aggregated get invalid
    /// after this operation. The function returns the index of the
    /// aggregated minorant along with the coefficients of the convex
    /// combination. The index of the new aggregated minorant might or
    /// might not be one of indices of the original minorants.

    ///
    /// # Error The indices of the minorants `mins` must belong to
    /// subproblem `fidx`.
    fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(Self::MinorantIndex, DVector)>;

    /// Return the (primal) optimal solution $\\|d\^*\\|$.
    fn get_primopt(&self) -> DVector;

    /// Return the value of the linear model in the optimal solution.
    ///
    /// This is the term $\langle g\^*, d\^* \rangle$ where $g^\*$ is

        let mut result = NEG_INFINITY;
        for m in &self.minorants {
            result = result.max(m.eval(y));
        }
        return result;
    }

    fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(usize, DVector)> {
        assert!(fidx == 0);
        if mins.len() == 2 {
            debug!("Aggregate");
            debug!("  {} * {}", self.opt_mult[0], self.minorants[0]);
            debug!("  {} * {}", self.opt_mult[1], self.minorants[1]);
            self.minorants[0] = self.minorants[0].combine(self.opt_mult[0], self.opt_mult[1], &self.minorants[1]);
            self.minorants.truncate(1);
            self.opt_mult.truncate(1);
            let coeffs = self.opt_mult.clone();
            self.opt_mult[0] = 1.0;
            debug!("  {}", self.minorants[0]);
            Ok((0, coeffs))
        } else if mins.len() == 1 {
            Ok((mins[0], dvec![1.0]))
        } else {
            panic!("No minorants specified to be aggregated");
        }
    }

    fn move_center(&mut self, alpha: Real, d: &DVector) {
        for m in self.minorants.iter_mut() {

    /// Return the value of the current model at the given point.
    fn eval_model(&self, y: &DVector) -> Real;

    /// Aggregate the given minorants according to the current solution.
    ///
    /// The (indices of the) minorants to be aggregated get invalid
    /// after this operation. The function returns the index of the
    /// aggregated minorant along with the coefficients of the convex
    /// combination. The index of the new aggregated minorant might or
    /// might not be one of indices of the original minorants.

    ///
    /// # Error
    /// The indices of the minorants `mins` must belong to subproblem `fidx`.
    fn aggregate(&mut self, fidx: usize, mins: &[usize]) -> Result<(Self::MinorantIndex, DVector), Self::Error>;

    /// Move the center of the master problem along $\alpha \cdot d$.
    fn move_center(&mut self, alpha: Real, d: &DVector);
}

 * 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(dead_code)]

use {Real, Vector, DVector, Minorant, FirstOrderProblem, SimpleEvaluation};
use mcf;

use std::fs::File;
use std::io::{self, Read};
use std::result;

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


/// Information about a minorant.
#[derive(Debug, Clone)]
struct MinorantInfo<Pr> {
    /// The minorant's index in the master problem
    index: usize,
    /// Current multiplier.
    multiplier: usize,
    /// Primal associated with this minorant.
    primal: Option<Pr>,
}

/**
 * Implementation of a bundle method.
 */
pub struct Solver<P, Pr, E>
    where P : for <'a> FirstOrderProblem<'a,Primal=Pr,EvalResult=E>,
          E : Evaluation<Pr>,
{
    /// The first order problem description.
    problem : P,

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

     */
    start_time : Instant,

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

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


impl<P, Pr, E> Solver<P, Pr, E>
    where P : for <'a> FirstOrderProblem<'a, Primal=Pr,EvalResult=E>,
          E : Evaluation<Pr>
{
    /**
     * Create a new solver for the given problem.
     *
     * Note that the solver owns the problem, so you cannot use the
     * same problem description elsewhere as long as it is assigned to
     * the solver. However, it is possible to get a reference to the
     * internally stored problem using `Solver::problem()`.
     */
    pub fn new_params(problem : P, params : SolverParams) -> Result<Solver<P, Pr, E>> {
        Ok(Solver{
            problem: problem,
            params:  params,
            terminator: Box::new(StandardTerminator { termination_precision: 1e-3 }),
            weighter: Box::new(HKWeighter::new()),
            cur_y : dvec![],
            cur_val : 0.0,

                Err(err) => return Err(Error::Master(Box::new(err))),
            },
            minorants: vec![],
        })
    }

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

    /**
     * Set the first order problem description associated with this
     * solver.
     *

            }
        }

        try!(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 = Vec::with_capacity(m);
        for _ in 0..m { self.minorants.push(vec![]); }

        self.cur_val = 0.0;
        for i in 0..m {
            let result = match self.problem.evaluate(i, &self.cur_y, INFINITY, 0.0) {
                Ok(r) => r,
                Err(err) => return Err(Error::Eval(Box::new(err))),
            };
            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];
                self.minorants[i].push(MinorantInfo {
                    index: try!(self.master.add_minorant(i, minorant)),
                    multiplier: 0,
                    primal: Some(primal),
                });
            } else {
                return Err(Error::NoMinorant);
            }
        }

        self.cur_valid = true;

        for i in 0..self.problem.num_subproblems() {
            let n = self.master.num_minorants(i);
            if n >= self.params.max_bundle_size {
                for m in self.minorants[i].iter_mut() {
                    m.multiplier = (1e6 * self.master.multiplier(m.index)) as usize;
                }
                self.minorants[i].sort_by_key(|m| -(m.multiplier as isize));
                let aggr = self.minorants[i].split_off(self.params.max_bundle_size-2);
                let (aggr_mins, aggr_primals) : (Vec<_>, Vec<_>) = aggr.into_iter().map(|m| {
                    (m.index, m.primal.unwrap())
                }).unzip();
                let (aggr_min, aggr_coeffs) = try!(self.master.aggregate(i, &aggr_mins));
                self.minorants[i].push(MinorantInfo{
                    index: aggr_min,
                    multiplier: 0,
                    primal: Some(self.problem.aggregate_primals(&aggr_coeffs, aggr_primals)),
                });
            }
        }
        Ok(())
    }

    /// Perform a descent step.

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

    /// Perform one bundle iteration.
    pub fn step(&mut self) -> Result<Step>
    {
        if !self.cur_valid {
            // current point needs new evaluation
            try!(self.init_master());
        }

        try!(self.solve_model());
        if self.terminator.terminate(&current_state!(self, Step::Term), &self.params) {

                Ok(r) => r,
                Err(err) => return Err(Error::Eval(Box::new(err))),
            };

            let fun_ub = result.objective();

            let mut minorants = result.into_iter();
            let mut nxt_minorant;
            let nxt_primal;
            match minorants.next() {
                Some((m, p)) => {
                    nxt_minorant = m;
                    nxt_primal = p;
                },
                None => return Err(Error::NoMinorant)

            }
            let fun_lb = nxt_minorant.constant;

            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;
            self.minorants[fidx].push(MinorantInfo{
                index: try!(self.master.add_minorant(fidx, nxt_minorant)),
                multiplier: 0,
                primal: 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);
            self.cur_val = self.new_cutval;
        }