RsBundle  Check-in [f0ef5ba392]

//!   trycpx!(CPXNETchgobj(cplex::env(), net, CPX_MAX))
//! }
//! ```

#![allow(dead_code)]

use libc::{c_int, c_double, c_char};
use std::ptr;
use std::fmt;
use std::error;

pub use std::result::Result as cplex_std_Result;
pub use std::convert::From as cplex_std_From;

pub const CPXMSGBUFSIZE: usize = 1024;

pub const CPX_ALG_NET: CPXINT = 3;
pub const CPX_ALG_BARRIER: CPXINT = 4;

/// Globally unique environment.
static mut ENV: *mut CPXenv = 0 as *mut CPXenv;

pub unsafe fn env() -> *mut CPXenv {
    if ENV == ptr::null_mut() {
        let mut status: c_int = 0;
        ENV = CPXopenCPLEX(&mut status);
        if status != 0 {
            panic!("Can't open CPLEX environment");
        }
    }
    ENV

            eps_weight: 1e30,
            m_r: m_r,
            iter: 0,
            model_max: NEG_INFINITY,
        }
    }
}







impl Weighter for HKWeighter {
    fn weight(&mut self, state: &BundleState, params: &SolverParams) -> Real {
        assert!(params.min_weight > 0.0);
        assert!(params.max_weight >= params.min_weight);

        debug!("HKWeighter {:?} iter:{}", state.step, self.iter);

                   sgnorm,
                   state.cur_val,
                   state.new_cutval,
                   lin_err,
                   self.eps_weight);
            debug!("  new_weight={}", new_weight);

            return new_weight;
        } else {
            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.weight / 2.0
                } else if 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);
            }

            debug!("  model_max={}", self.model_max);
            debug!("  new_weight={}", new_weight);

            return new_weight;
        }
    }
}

// Copyright (c) 2016 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 {Real, DVector, Minorant};
use master::master::{MasterProblem, Error, Result};
use master::UnconstrainedMasterProblem;

use std::f64::{INFINITY, NEG_INFINITY};

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

            self.eta[i] = neweta;
        }

        debug!("Eta update");
        debug!("  primopt={}", self.primopt);
        debug!("  eta    ={}", self.eta);

        return updated_eta;
    }

    // Compute the new candidate point.
    //
    // This consists of two steps:
    //
    // 1. the new point is computed as $-\tfrac{1}{u}\bar{g}$, where $\bar{g}$

    fn eta_cutval(&self) -> Real {
        let mut val = 0.0;
        for i in 0..self.lb.len() {
            if self.eta[i] >= 0.0 {
                if self.lb[i] > NEG_INFINITY {
                    val += self.lb[i] * self.eta[i];
                }
            } else {
                if self.ub[i] < INFINITY {
                    val += self.ub[i] * self.eta[i];
                }
            }
        }
        return val;
    }


    /**
     * 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 {
        let dualopt = self.master.dualopt();
        let mut norm2 = 0.0;
        for i in 0..self.eta.len() {
            let x = dualopt[i] - self.eta[i];
            norm2 += x * x;
        }
        return norm2;
    }
}


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

    fn multiplier(&self, min: Self::MinorantIndex) -> Real {
        self.master.multiplier(min)
    }

    fn move_center(&mut self, alpha: Real, d: &DVector) {
        self.need_new_candidate = true;
        self.master.move_center(alpha, d);
        for i in 0..self.primopt.len() {
            self.lb[i] -= alpha * d[i];
            self.ub[i] -= alpha * d[i];
        }
    }

    fn set_max_updates(&mut self, max_updates: usize) {
        BoxedMasterProblem::set_max_updates(self, max_updates);
    }

    fn cnt_updates(&self) -> usize {
        self.cnt_updates
    }
}

// Copyright (c) 2016 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

        if nvars == 0 {
            return Err(Error::NoMinorants);
        }
        // update linear costs
        {
            let mut c = Vec::with_capacity(nvars);
            let mut inds = Vec::with_capacity(nvars);
            for mins in self.minorants.iter() {
                for m in mins {
                    inds.push(c.len() as c_int);
                    c.push(-m.constant * self.weight - m.linear.dot(eta));
                }
            }
            trycpx!(CPXchgobj(env(), self.lp, nvars as c_int, inds.as_ptr(), c.as_ptr()));
        }

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


impl CplexMaster {
    fn init_qp(&mut self) -> Result<()> {
        if self.lp != ptr::null_mut() {
            trycpx!(CPXfreeprob(env(), &mut self.lp));
        }
        trycpx!({
            let mut status = 0;
            self.lp = CPXcreateprob(env(), &mut status, CString::new("mcf").unwrap().as_ptr());
            status
        });

                        self.qterm[idx_i][idx_j] = x;
                        self.qterm[idx_j][idx_i] = x;
                    }
                }
            }

            if cfg!(debug_assertions) {
                for &idx_i in activeinds.iter() {
                    for &idx_j in activeinds.iter() {
                        let (fidx_i, i) = self.index2min[idx_i];
                        let (fidx_j, j) = self.index2min[idx_j];
                        let m_i = &self.minorants[fidx_i][i];
                        let m_j = &self.minorants[fidx_j][j];
                        debug_assert!((self.qterm[idx_i][idx_j] - m_i.linear.dot(&m_j.linear)).abs() < 1e-6);
                    }
                }

    }

    fn eval_model(&self, y: &DVector) -> Real {
        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]);

            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() {
            m.move_center(alpha, d);
        }
    }
}

// Copyright (c) 2016 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\\}, \\]
//!

//! * solving the problem,
//! * 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 master;
pub use self::master::{MasterProblem, Error, Result};

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

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

// Copyright (c) 2016 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(dead_code)]

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

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

            rhsval: 0.0,
            cbase: cbase,
            c: vec![dvec![]; ncom],
        })
    }

    /// Compute costs for a primal solution.
    pub fn get_primal_costs(&self, fidx: usize, primals: &Vec<DVector>) -> Real {
        if self.multimodel {
            primals[0].iter().enumerate().map(|(i, x)| x * self.cbase[fidx][i]).sum()
        } else {
            let mut sum = 0.0;
            for (fidx, p) in primals.iter().enumerate() {
                for (i, x) in p.iter().enumerate() {
                    sum += x * self.cbase[fidx][i];

        let mut net = ptr::null_mut();
        let logfile;

        unsafe {
            loop {
                logfile = CPXfopen(CString::new("mcf.cpxlog").unwrap().as_ptr(),
                                   CString::new("w").unwrap().as_ptr());
                if logfile == ptr::null_mut() {
                    return Err(Error::Cplex(CplexError {
                        code: 0,
                        msg: "Can't open log-file".to_string(),
                    }));
                }
                status = CPXsetlogfile(env(), logfile);
                if status != 0 {

                return Err(Error::Cplex(CplexError {
                    code: status,
                    msg: CString::from_raw(msg).to_string_lossy().into_owned(),
                }));
            }
        }

        return Ok(Solver {
            net: net,
            logfile: logfile,
        });
    }

    pub fn num_nodes(&self) -> usize {
        unsafe { CPXNETgetnumnodes(env(), self.net) as usize }
    }

    pub fn num_arcs(&self) -> usize {

                               self.net,
                               &mut stat as *mut c_int,
                               &mut objval as *mut c_double,
                               sol.as_mut_ptr(),
                               ptr::null_mut(),
                               ptr::null_mut(),
                               ptr::null_mut()));
        return Ok(sol);
    }

    pub fn writelp(&self, filename: &str) -> Result<()> {
        Ok(trycpx!(CPXNETwriteprob(env(),
                                   self.net,
                                   CString::new(filename).unwrap().as_ptr(),
                                   ptr::null_mut())))
    }
}

// Copyright (c) 2016 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

        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() {
            for m in self.minorants[i].iter_mut() {
                m.multiplier = self.master.multiplier(m.index);
            }
        }

        debug!("Model result");
        debug!("  cur_val ={}", self.cur_val);
        debug!("  nxt_mod ={}", self.nxt_mod);

            if nxt_ub > descent_bnd {
                // upper bound will produce null-step
                if self.cur_val - nxt_lb > (self.cur_val - self.nxt_mod) * self.params.nullstep_factor.max(0.5) {
                    warn!("Relative precision of returned objective interval enforces null-step.");
                    self.iterinfos.push(IterationInfo::UpperBoundNullStep);
                }
            }
        } else {

            // lower bound gives already a null step
            if self.cur_val - nxt_lb > 0.8 * (self.cur_val - self.nxt_mod) {
                // subgradient won't yield much improvement
                warn!("Shallow cut (subgradient won't yield much improvement)");
                self.iterinfos.push(IterationInfo::ShallowCut);
            }
        }

        debug!("Step");
        debug!("  cur_val    ={}", self.cur_val);
        debug!("  nxt_mod    ={}", self.nxt_mod);
        debug!("  nxt_ub     ={}", self.nxt_val);
        debug!("  descent_bnd={}", descent_bnd);

        // do a descent step or null step
        if nxt_ub <= descent_bnd {
            self.descent_step();
            return Ok(Step::Descent);
        } else {
            self.null_step();
            return Ok(Step::Null);
        }
    }

    /**
     * Return the bound on the function value that enforces a
     * nullstep.
     *

// Copyright (c) 2016 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

        elems: Vec<(usize, Real)>,
    },
}


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();
                try!(write!(f, "{}:(", size));
                if let Some(&(i, x)) = it.next() {
                    try!(write!(f, "{}:{}", i, x));
                    for &(i, x) in it {
                        try!(write!(f, ", {}:{}", i, x));
                    }

    /**
     * Convert vector to a dense vector.
     *
     * This function always returns a copy of the vector.
     */
    pub fn to_dense(&self) -> DVector {
        match self {
            &Vector::Dense(ref x) => x.clone(),
            &Vector::Sparse { size: n, elems: ref xs } => {
                let mut v = vec![0.0; n];
                for &(i, x) in xs {
                    unsafe { *v.get_unchecked_mut(i) = x };
                }
                DVector(v)
            }
        }
    }
}