Index: .fossil-settings/ignore-glob ================================================================== --- .fossil-settings/ignore-glob +++ .fossil-settings/ignore-glob @@ -1,4 +1,5 @@ target/ *.log *.cpxlog Cargo.lock +instances/ Index: Cargo.toml ================================================================== --- Cargo.toml +++ Cargo.toml @@ -1,14 +1,15 @@ [package] name = "bundle" -version = "0.5.0" +version = "0.7.0-dev" +edition = "2018" authors = ["Frank Fischer "] [dependencies] -itertools = "^0.7.2" +itertools = "^0.8" libc = "^0.2.6" log = "^0.4" c_str_macro = "^1.0" -cplex-sys = "^0.5" +cplex-sys = "^0.6" [dev-dependencies] -env_logger = "^0.5.3" +env_logger = "^0.7" Index: examples/mmcf.rs ================================================================== --- examples/mmcf.rs +++ examples/mmcf.rs @@ -13,17 +13,17 @@ * * You should have received a copy of the GNU General Public License * along with this program. If not, see */ -extern crate bundle; -extern crate env_logger; -#[macro_use] -extern crate log; +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(); @@ -41,11 +41,12 @@ max_bundle_size: 25, min_weight: 1e-3, max_weight: 100.0, ..Default::default() }, - ).unwrap(); + ) + .unwrap(); solver.terminator = Box::new(StandardTerminator { termination_precision: 1e-6, }); solver.solve().unwrap(); @@ -52,11 +53,12 @@ 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(); + }) + .sum(); info!("Primal costs: {}", costs); } else { panic!("Usage: {} FILENAME", program); } } Index: examples/quadratic.rs ================================================================== --- examples/quadratic.rs +++ examples/quadratic.rs @@ -15,15 +15,13 @@ * along with this program. If not, see */ use std::error::Error; -#[macro_use] -extern crate bundle; -extern crate env_logger; -#[macro_use] -extern crate log; +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], @@ -94,8 +92,9 @@ SolverParams { min_weight: 1.0, max_weight: 1.0, ..Default::default() }, - ).unwrap(); + ) + .unwrap(); solver.solve().unwrap(); } Index: src/firstorderproblem.rs ================================================================== --- src/firstorderproblem.rs +++ src/firstorderproblem.rs @@ -1,6 +1,6 @@ -// Copyright (c) 2016, 2017 Frank Fischer +// Copyright (c) 2016, 2017, 2019 Frank Fischer // // 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. @@ -14,12 +14,12 @@ // along with this program. If not, see // //! Problem description of a first-order convex optimization problem. -use solver::UpdateState; -use {Minorant, Real}; +use crate::solver::UpdateState; +use crate::{Minorant, Real}; use std::result::Result; use std::vec::IntoIter; /** @@ -65,12 +65,14 @@ /// 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. -#[derive(Debug, Clone, Copy)] -pub enum Update { +/// +/// - `P` is the type of primals +/// - `E` is the error type +pub enum Update { /// Add a variable with bounds. /// /// The initial value of the variable will be the feasible value /// closest to 0. AddVariable { lower: Real, upper: Real }, @@ -77,10 +79,15 @@ /// 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 Result<(), E>>, + }, } /** * Trait for implementing a first-order problem description. * @@ -175,21 +182,29 @@ } /// Return updates of the problem. /// /// The default implementation returns no updates. - fn update(&mut self, _state: &UpdateState) -> Result, Self::Err> { + fn update( + &mut self, + _state: &UpdateState, + ) -> Result>, Self::Err> { Ok(vec![]) } /// Return new components for a subgradient. /// - /// The components are typically generated by some primal - /// information. The corresponding primal is passed as a + /// 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, _primal: &Self::Primal, _vars: &[usize]) -> Result, Self::Err> { + fn extend_subgradient( + &mut self, + _i: usize, + _primal: &Self::Primal, + _vars: &[usize], + ) -> Result, Self::Err> { unimplemented!() } } Index: src/hkweighter.rs ================================================================== --- src/hkweighter.rs +++ src/hkweighter.rs @@ -20,12 +20,14 @@ //! //! > Helmberg, C. and Kiwiel, K.C. (2002): A spectral bundle method //! > with bounds, Math. Programming A 93, 173--194 //! -use Real; -use {BundleState, SolverParams, Step, Weighter}; +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; @@ -81,11 +83,12 @@ 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) + } + .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; @@ -99,11 +102,12 @@ 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(FACTOR * state.weight) .min(params.max_weight); if new_weight > state.weight { self.iter = -1 } else { self.iter = min(self.iter - 1, -1); @@ -122,11 +126,12 @@ w } else if self.iter > 3 || state.nxt_val < self.model_max { state.weight / 2.0 } else { state.weight - }.max(state.weight / FACTOR) + } + .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; Index: src/lib.rs ================================================================== --- src/lib.rs +++ src/lib.rs @@ -14,17 +14,10 @@ // along with this program. If not, see // //! Proximal bundle method implementation. -#[macro_use] -extern crate c_str_macro; -extern crate itertools; - -#[macro_use] -extern crate log; - /// Type used for real numbers throughout the library. pub type Real = f64; #[macro_export] macro_rules! dvec { @@ -31,28 +24,25 @@ ( $ elem : expr ; $ n : expr ) => { DVector(vec![$elem; $n]) }; ( $ ( $ x : expr ) , * ) => { DVector(vec![$($x),*]) }; ( $ ( $ x : expr , ) * ) => { DVector(vec![$($x,)*]) }; } -#[macro_use] -extern crate cplex_sys; - pub mod vector; -pub use vector::{DVector, Vector}; +pub use crate::vector::{DVector, Vector}; pub mod minorant; -pub use minorant::Minorant; +pub use crate::minorant::Minorant; pub mod firstorderproblem; -pub use firstorderproblem::{Evaluation, FirstOrderProblem, SimpleEvaluation, Update}; +pub use crate::firstorderproblem::{Evaluation, FirstOrderProblem, SimpleEvaluation, Update}; pub mod solver; -pub use solver::{ +pub use crate::solver::{ BundleState, IterationInfo, Solver, SolverParams, StandardTerminator, Step, Terminator, UpdateState, Weighter, }; mod hkweighter; -pub use hkweighter::HKWeighter; +pub use crate::hkweighter::HKWeighter; mod master; pub mod mcf; Index: src/master/base.rs ================================================================== --- src/master/base.rs +++ src/master/base.rs @@ -1,6 +1,6 @@ -// Copyright (c) 2016, 2017 Frank Fischer +// Copyright (c) 2016, 2017, 2019 Frank Fischer // // 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. @@ -12,11 +12,11 @@ // // You should have received a copy of the GNU General Public License // along with this program. If not, see // -use {DVector, Minorant, Real}; +use crate::{DVector, Minorant, Real}; use std::error::Error; use std::fmt; use std::result; @@ -44,11 +44,11 @@ } } } impl Error for MasterProblemError { - fn cause(&self) -> Option<&Error> { + fn cause(&self) -> Option<&dyn Error> { use self::MasterProblemError::*; match self { SubgradientExtension(err) | Solver(err) | Custom(err) => Some(err.as_ref()), NoMinorants => None, } @@ -56,10 +56,13 @@ } /// Result type of master problems. pub type Result = result::Result; +/// Callback for subgradient extensions. +pub type SubgradientExtension<'a, I> = dyn FnMut(usize, I, &[usize]) -> result::Result> + 'a; + pub trait MasterProblem { /// Unique index for a minorant. type MinorantIndex: Copy + Eq; /// Set the number of subproblems. @@ -81,18 +84,18 @@ fn set_max_updates(&mut self, max_updates: usize) -> Result<()>; /// Return the current number of inner iterations. fn cnt_updates(&self) -> usize; - /// Add or movesome variables with bounds. + /// 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, Real, Real)], - extend_subgradient: &mut FnMut(usize, Self::MinorantIndex, &[usize]) -> result::Result>, + extend_subgradient: &mut SubgradientExtension, ) -> Result<()>; /// Add a new minorant to the model. /// /// The function returns a unique (among all minorants of all Index: src/master/boxed.rs ================================================================== --- src/master/boxed.rs +++ src/master/boxed.rs @@ -1,6 +1,6 @@ -// Copyright (c) 2016, 2017, 2018 Frank Fischer +// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer // // 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. @@ -12,19 +12,19 @@ // // You should have received a copy of the GNU General Public License // along with this program. If not, see // -use master::MasterProblem; -use master::UnconstrainedMasterProblem; -use {DVector, Minorant, Real}; +use crate::master::UnconstrainedMasterProblem; +use crate::master::{MasterProblem, SubgradientExtension}; +use crate::{DVector, Minorant, Real}; use super::Result; + use itertools::multizip; -use std::error::Error; +use log::debug; use std::f64::{EPSILON, INFINITY, NEG_INFINITY}; -use std::result; /** * Turn unconstrained master problem into box-constrained one. * * This master problem adds box constraints to an unconstrainted @@ -159,22 +159,22 @@ } else if ub < INFINITY { ub * eta } else { 0.0 } - }).sum() + }) + .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 { - let dualopt = self.master.dualopt(); - dualopt.iter().zip(self.eta.iter()).map(|(x, y)| x * y).sum() + self.eta.dot(self.master.dualopt()) } } impl MasterProblem for BoxedMasterProblem { type MinorantIndex = M::MinorantIndex; @@ -208,11 +208,11 @@ } fn add_vars( &mut self, bounds: &[(Option, Real, Real)], - extend_subgradient: &mut FnMut(usize, Self::MinorantIndex, &[usize]) -> result::Result>, + extend_subgradient: &mut SubgradientExtension, ) -> 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; @@ -227,11 +227,11 @@ } else { Ok(()) } } - #[cfg_attr(feature = "cargo-clippy", allow(cyclomatic_complexity))] + #[allow(clippy::cognitive_complexity)] fn solve(&mut self, center_value: Real) -> Result<()> { debug!("Solve Master"); debug!(" lb ={}", self.lb); debug!(" ub ={}", self.ub); Index: src/master/cpx.rs ================================================================== --- src/master/cpx.rs +++ src/master/cpx.rs @@ -1,6 +1,6 @@ -// Copyright (c) 2016, 2017, 2018 Frank Fischer +// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer // // 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. @@ -16,22 +16,25 @@ //! Master problem implementation using CPLEX. #![allow(unused_unsafe)] -use master::{Error as MasterProblemError, UnconstrainedMasterProblem}; -use {DVector, Minorant, Real}; - -use cplex_sys as cpx; +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::error::Error; use std::f64::{self, NEG_INFINITY}; +use std::iter::{once, repeat}; use std::os::raw::{c_char, c_int}; use std::ptr; -use std::result; impl From for MasterProblemError { fn from(err: cpx::CplexError) -> MasterProblemError { MasterProblemError::Solver(err.into()) } @@ -155,11 +158,11 @@ fn add_vars( &mut self, nnew: usize, changed: &[usize], - extend_subgradient: &mut FnMut(usize, Self::MinorantIndex, &[usize]) -> result::Result>, + extend_subgradient: &mut SubgradientExtension, ) -> Result<()> { debug_assert!(!self.minorants[0].is_empty()); let noldvars = self.minorants[0][0].linear.len(); let nnewvars = noldvars + nnew; @@ -167,12 +170,12 @@ 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, i, &changedvars).map_err(MasterProblemError::SubgradientExtension)?; + 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()..]); } @@ -205,11 +208,11 @@ Ok(()) } fn solve(&mut self, eta: &DVector, _fbound: Real, _augbound: Real, _relprec: Real) -> Result<()> { if self.force_update || !self.updateinds.is_empty() { - try!(self.init_qp()); + self.init_qp()?; } let nvars = unsafe { cpx::getnumcols(cpx::env(), self.lp) as usize }; if nvars == 0 { return Err(MasterProblemError::NoMinorants); @@ -232,10 +235,14 @@ 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); @@ -294,11 +301,12 @@ let aggr_coeffs = if sum_coeffs != 0.0 { mins.iter() .map(|&i| self.opt_mults[fidx][self.index2min[i].1] / sum_coeffs) .collect::() } else { - dvec![0.0; mins.len()] + // 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() { Index: src/master/minimal.rs ================================================================== --- src/master/minimal.rs +++ src/master/minimal.rs @@ -1,6 +1,6 @@ -// Copyright (c) 2016, 2017 Frank Fischer +// Copyright (c) 2016, 2017, 2019 Frank Fischer // // 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. @@ -12,14 +12,17 @@ // // You should have received a copy of the GNU General Public License // along with this program. If not, see // -use master::{Error as MasterProblemError, UnconstrainedMasterProblem}; -use {DVector, Minorant, Real}; +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; @@ -123,11 +126,11 @@ fn add_vars( &mut self, nnew: usize, changed: &[usize], - extend_subgradient: &mut FnMut(usize, Self::MinorantIndex, &[usize]) -> result::Result>, + extend_subgradient: &mut SubgradientExtension, ) -> Result<()> { if !self.minorants.is_empty() { let noldvars = self.minorants[0].linear.len(); let mut changedvars = vec![]; changedvars.extend_from_slice(changed); Index: src/master/mod.rs ================================================================== --- src/master/mod.rs +++ src/master/mod.rs @@ -1,6 +1,6 @@ -// Copyright (c) 2016, 2017 Frank Fischer +// Copyright (c) 2016, 2017, 2019 Frank Fischer // // 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. @@ -12,11 +12,10 @@ // // You should have received a copy of the GNU General Public License // along with this program. If not, see // -#![cfg_attr(feature = "cargo-clippy", allow(doc_markdown))] //! Bundle master problem solver. //! //! This module contains solvers for the bundle master problem, i.e. //! for solving convex optimization problems of the form //! @@ -36,11 +35,11 @@ //! * 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}; +pub use self::base::{MasterProblem, MasterProblemError as Error, Result, SubgradientExtension}; mod boxed; pub use self::boxed::BoxedMasterProblem; mod unconstrained; @@ -51,6 +50,6 @@ // pub mod grb; // pub use master::grb::GurobiMaster; mod cpx; -pub use master::cpx::CplexMaster; +pub use crate::master::cpx::CplexMaster; Index: src/master/unconstrained.rs ================================================================== --- src/master/unconstrained.rs +++ src/master/unconstrained.rs @@ -1,6 +1,6 @@ -// Copyright (c) 2016, 2017 Frank Fischer +// Copyright (c) 2016, 2017, 2019 Frank Fischer // // 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. @@ -12,16 +12,13 @@ // // You should have received a copy of the GNU General Public License // along with this program. If not, see // -use {DVector, Minorant, Real}; +use crate::{DVector, Minorant, Real}; -use super::Result; - -use std::error::Error; -use std::result; +use super::{Result, SubgradientExtension}; /** * Trait for master problems without box constraints. * * Implementors of this trait are supposed to solve quadratic @@ -79,11 +76,11 @@ /// are added. fn add_vars( &mut self, nnew: usize, changed: &[usize], - extend_subgradient: &mut FnMut(usize, Self::MinorantIndex, &[usize]) -> result::Result>, + extend_subgradient: &mut SubgradientExtension, ) -> Result<()>; /// Solve the master problem. fn solve(&mut self, eta: &DVector, fbound: Real, augbound: Real, relprec: Real) -> Result<()>; Index: src/mcf/mod.rs ================================================================== --- src/mcf/mod.rs +++ src/mcf/mod.rs @@ -15,9 +15,9 @@ // //! Multi-commodity min-cost-flow subproblems. mod solver; -pub use mcf::solver::Solver; +pub use crate::mcf::solver::Solver; mod problem; -pub use mcf::problem::MMCFProblem; +pub use crate::mcf::problem::MMCFProblem; Index: src/mcf/problem.rs ================================================================== --- src/mcf/problem.rs +++ src/mcf/problem.rs @@ -1,6 +1,6 @@ -// Copyright (c) 2016, 2017 Frank Fischer +// Copyright (c) 2016, 2017, 2019 Frank Fischer // // 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. @@ -12,12 +12,14 @@ // // You should have received a copy of the GNU General Public License // along with this program. If not, see // -use mcf; -use {DVector, FirstOrderProblem, Minorant, Real, SimpleEvaluation}; +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; @@ -37,11 +39,11 @@ } impl Error for MMCFFormatError {} /// Result type of the MMCFProblem. -pub type Result = result::Result>; +pub type Result = result::Result>; #[derive(Clone, Copy, Debug)] struct ArcInfo { arc: usize, src: usize, @@ -67,22 +69,23 @@ impl MMCFProblem { pub fn read_mnetgen(basename: &str) -> Result { let mut buffer = String::new(); { - let mut f = try!(File::open(&format!("{}.nod", basename))); - try!(f.read_to_string(&mut buffer)); + let mut f = File::open(&format!("{}.nod", basename))?; + f.read_to_string(&mut buffer)?; } let fnod = buffer .split_whitespace() .map(|x| x.parse::().unwrap()) .collect::>(); if fnod.len() != 4 { return Err(MMCFFormatError { msg: format!("Expected 4 numbers in {}.nod, but got {}", basename, fnod.len()), - }.into()); + } + .into()); } let ncom = fnod[0]; let nnodes = fnod[1]; let narcs = fnod[2]; @@ -89,45 +92,45 @@ let ncaps = fnod[3]; // read nodes let mut nets = Vec::with_capacity(ncom); for _ in 0..ncom { - nets.push(try!(mcf::Solver::new(nnodes))) + nets.push(mcf::Solver::new(nnodes)?) } { - let mut f = try!(File::open(&format!("{}.sup", basename))); + let mut f = File::open(&format!("{}.sup", basename))?; buffer.clear(); - try!(f.read_to_string(&mut buffer)); + f.read_to_string(&mut buffer)?; } for line in buffer.lines() { let mut data = line.split_whitespace(); - let node = try!(data.next().unwrap().parse::()); - let com = try!(data.next().unwrap().parse::()); - let supply = try!(data.next().unwrap().parse::()); - try!(nets[com - 1].set_balance(node - 1, supply)); + let node = data.next().unwrap().parse::()?; + let com = data.next().unwrap().parse::()?; + let supply = data.next().unwrap().parse::()?; + 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 = try!(File::open(&format!("{}.arc", basename))); + let mut f = File::open(&format!("{}.arc", basename))?; buffer.clear(); - try!(f.read_to_string(&mut buffer)); + f.read_to_string(&mut buffer)?; } for line in buffer.lines() { let mut data = line.split_whitespace(); - let arc = try!(data.next().unwrap().parse::()) - 1; - let src = try!(data.next().unwrap().parse::()) - 1; - let snk = try!(data.next().unwrap().parse::()) - 1; - let com = try!(data.next().unwrap().parse::()) - 1; - let cost = try!(data.next().unwrap().parse::()); - let cap = try!(data.next().unwrap().parse::()); - let mt = try!(data.next().unwrap().parse::()) - 1; + let arc = data.next().unwrap().parse::()? - 1; + let src = data.next().unwrap().parse::()? - 1; + let snk = data.next().unwrap().parse::()? - 1; + let com = data.next().unwrap().parse::()? - 1; + let cost = data.next().unwrap().parse::()?; + let cap = data.next().unwrap().parse::()?; + let mt = data.next().unwrap().parse::()? - 1; assert!( arc < narcs, format!("Wrong arc number (got: {}, expected in 1..{})", arc + 1, narcs) ); // set internal coeff @@ -136,11 +139,11 @@ arc: arc + 1, src: src + 1, snk: snk + 1, }); // add arc - try!(nets[com].add_arc(src, snk, cost, if cap < 0.0 { INFINITY } else { cap })); + 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); @@ -147,19 +150,19 @@ } } // read rhs of coupling constraints { - let mut f = try!(File::open(&format!("{}.mut", basename))); + let mut f = File::open(&format!("{}.mut", basename))?; buffer.clear(); - try!(f.read_to_string(&mut buffer)); + 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 = try!(data.next().unwrap().parse::()) - 1; - let cap = try!(data.next().unwrap().parse::()); + let mt = data.next().unwrap().parse::()? - 1; + let cap = data.next().unwrap().parse::()?; rhs[mt] = cap; } // set lhs let mut lhs = vec![vec![vec![]; ncom]; ncaps]; @@ -206,11 +209,12 @@ .iter() .map(|x| { let mut r = dvec![]; r.scal(primals[0].0, x); r - }).collect::>(); + }) + .collect::>(); for &(alpha, primal) in &primals[1..] { for (j, x) in primal.iter().enumerate() { aggr[j].add_scaled(alpha, x); } @@ -266,32 +270,32 @@ } debug!("y={:?}", y); for i in 0..self.nets.len() { debug!("c[{}]={}", i, self.c[i]); - try!(self.nets[i].set_objective(&self.c[i])); + self.nets[i].set_objective(&self.c[i])?; } // solve subproblems for (i, net) in self.nets.iter_mut().enumerate() { - try!(net.solve()); - debug!("c[{}]={}", i, try!(net.objective())); + 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 - try!(self.nets[fidx].objective()); + objective = self.rhsval - self.nets[fidx].objective()?; } else { subg = dvec![0.0; self.rhs.len()]; - objective = -try!(self.nets[fidx].objective()); + objective = -self.nets[fidx].objective()?; } - let sol = try!(self.nets[fidx].get_solution()); + 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::(); } Ok(SimpleEvaluation { @@ -306,12 +310,12 @@ }) } else { let mut objective = self.rhsval; let mut sols = Vec::with_capacity(self.nets.len()); for i in 0..self.nets.len() { - objective -= try!(self.nets[i].objective()); - sols.push(try!(self.nets[i].get_solution())); + 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() { Index: src/mcf/solver.rs ================================================================== --- src/mcf/solver.rs +++ src/mcf/solver.rs @@ -1,6 +1,6 @@ -// Copyright (c) 2016, 2017, 2018 Frank Fischer +// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer // // 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. @@ -14,13 +14,15 @@ // along with this program. If not, see // #![allow(unused_unsafe)] -use {DVector, Real}; +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; @@ -46,11 +48,11 @@ pub fn new(nnodes: usize) -> Result { let mut status: c_int; let mut net = ptr::null_mut(); unsafe { - #[cfg_attr(feature = "cargo-clippy", allow(never_loop))] + #[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; } @@ -75,11 +77,12 @@ 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()); + } + .into()); } } Ok(Solver { net }) } Index: src/minorant.rs ================================================================== --- src/minorant.rs +++ src/minorant.rs @@ -14,11 +14,11 @@ // along with this program. If not, see // //! A linear minorant. -use {DVector, Real}; +use crate::{DVector, Real}; use std::fmt; /** * A linear minorant of a convex function. @@ -40,11 +40,11 @@ pub linear: DVector, } impl fmt::Display for Minorant { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - try!(write!(f, "{} + y * {}", self.constant, self.linear)); + write!(f, "{} + y * {}", self.constant, self.linear)?; Ok(()) } } impl Default for Minorant { Index: src/solver.rs ================================================================== --- src/solver.rs +++ src/solver.rs @@ -1,6 +1,6 @@ -// Copyright (c) 2016, 2017, 2018 Frank Fischer +// Copyright (c) 2016, 2017, 2018, 2019 Frank Fischer // // 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. @@ -14,15 +14,17 @@ // along with this program. If not, see // //! The main bundle method solver. -use {DVector, Real}; -use {Evaluation, FirstOrderProblem, HKWeighter, Update}; +use crate::{DVector, Real}; +use crate::{Evaluation, FirstOrderProblem, HKWeighter, Update}; -use master::{BoxedMasterProblem, Error as MasterProblemError, MasterProblem, UnconstrainedMasterProblem}; -use master::{CplexMaster, MinimalMaster}; +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; @@ -48,10 +50,12 @@ 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 fmt::Display for SolverError { @@ -71,17 +75,22 @@ 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 Error for SolverError { - fn cause(&self) -> Option<&Error> { + fn cause(&self) -> Option<&dyn Error> { match self { SolverError::Evaluation(err) => Some(err), SolverError::Update(err) => Some(err), SolverError::Master(err) => Some(err), _ => None, @@ -92,10 +101,16 @@ impl From for SolverError { fn from(err: ParameterError) -> SolverError { SolverError::Parameter(err) } } + +impl From for SolverError { + fn from(err: MasterProblemError) -> SolverError { + SolverError::Master(err) + } +} /** * The current state of the bundle method. * * Captures the current state of the bundle method during the run of @@ -321,17 +336,15 @@ Term, } /// Information about a minorant. #[derive(Debug, Clone)] -struct MinorantInfo { +struct MinorantInfo { /// The minorant's index in the master problem index: usize, /// Current multiplier. multiplier: Real, - /// Primal associated with this minorant. - primal: Option, } /// Information about the last iteration. #[derive(Debug, Clone, Copy, PartialEq)] pub enum IterationInfo { @@ -341,11 +354,13 @@ } /// State information for the update callback. pub struct UpdateState<'a, Pr: 'a> { /// Current model minorants. - minorants: &'a [Vec>], + minorants: &'a [Vec], + /// The primals. + primals: &'a Vec>, /// 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 @@ -358,19 +373,19 @@ 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, m.primal.as_ref().unwrap())) + .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| m.primal.as_ref()) + self.minorants[fidx].last().and_then(|m| self.primals[m.index].as_ref()) } } /** * Implementation of a bundle method. @@ -381,14 +396,14 @@ /// The solver parameter. pub params: SolverParams, /// Termination predicate. - pub terminator: Box, + pub terminator: Box, /// Weighter heuristic. - pub weighter: Box, + pub weighter: Box, /// Lower and upper bounds of all variables. bounds: Vec<(Real, Real)>, /// Current center of stability. @@ -453,14 +468,17 @@ * This is actually the time of the last call to `Solver::init`. */ start_time: Instant, /// The master problem. - master: Box>, + master: Box>, /// The active minorant indices for each subproblem. - minorants: Vec>>, + minorants: Vec>, + + /// The primals associated with each global minorant index. + primals: Vec>, /// Accumulated information about the last iteration. iterinfos: Vec, } @@ -501,14 +519,13 @@ sgnorm: 0.0, expected_progress: 0.0, cnt_descent: 0, cnt_null: 0, start_time: Instant::now(), - master: Box::new(BoxedMasterProblem::new( - MinimalMaster::new().map_err(SolverError::Master)?, - )), + master: Box::new(BoxedMasterProblem::new(MinimalMaster::new()?)), minorants: vec![], + primals: vec![], iterinfos: vec![], }) } /// A new solver with default parameter. @@ -574,18 +591,27 @@ self.start_time = Instant::now(); Ok(()) } - /// Solve the problem. + /// 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> { 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> { + // First initialize the internal data structures. + self.init()?; - if self.solve_iter(LIMIT)? { + if self.solve_iter(iter_limit)? { Ok(()) } else { - Err(SolverError::IterationLimit { limit: LIMIT }) + Err(SolverError::IterationLimit { limit: iter_limit }) } } /// Solve the problem but stop after `niter` iterations. /// @@ -618,10 +644,11 @@ /// separate new constraints depending on the current solution. fn update_problem(&mut self, term: Step) -> Result> { 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 @@ -677,26 +704,40 @@ 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 minorants = &self.minorants; - self.master - .add_vars( - &newvars.iter().map(|v| (v.0, v.1, v.2)).collect::>(), - &mut |fidx, minidx, vars| { - problem - .extend_subgradient(minorants[fidx][minidx].primal.as_ref().unwrap(), vars) - .map(DVector) - .map_err(|e| e.into()) - }, - ).map_err(SolverError::Master)?; + let primals = &self.primals; + self.master.add_vars( + &newvars.iter().map(|v| (v.0, v.1, v.2)).collect::>(), + &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; @@ -718,11 +759,11 @@ /// 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, m.primal.as_ref().unwrap())) + .map(|m| (m.multiplier, self.primals[m.index].as_ref().unwrap())) .collect() } fn show_info(&self, step: Step) { let time = self.start_time.elapsed(); @@ -766,18 +807,14 @@ fn init_master(&mut self) -> Result<(), SolverError> { 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().map_err(SolverError::Master)?, - )) + Box::new(BoxedMasterProblem::new(MinimalMaster::new()?)) } else { debug!("Use CPLEX master problem"); - Box::new(BoxedMasterProblem::new( - CplexMaster::new().map_err(SolverError::Master)?, - )) + Box::new(BoxedMasterProblem::new(CplexMaster::new()?)) }; let lb = self.problem.lower_bounds().map(DVector); let ub = self.problem.upper_bounds().map(DVector); @@ -794,17 +831,13 @@ .unwrap_or(false) { return Err(SolverError::Dimension); } - self.master.set_num_subproblems(m).map_err(SolverError::Master)?; - self.master - .set_vars(self.problem.num_variables(), lb, ub) - .map_err(SolverError::Master)?; - self.master - .set_max_updates(self.params.max_updates) - .map_err(SolverError::Master)?; + 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 { @@ -817,15 +850,19 @@ 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: self.master.add_minorant(i, minorant).map_err(SolverError::Master)?, + index: minidx, multiplier: 0.0, - primal: Some(primal), }); + if minidx >= self.primals.len() { + self.primals.resize_with(minidx + 1, || None); + } + self.primals[minidx] = Some(primal); } else { return Err(SolverError::NoMinorant); } } @@ -837,27 +874,27 @@ // 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).map_err(SolverError::Master)?; - self.master.solve(self.cur_val).map_err(SolverError::Master)?; + 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).map_err(SolverError::Master)?; + 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> { - self.master.solve(self.cur_val).map_err(SolverError::Master)?; + 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; @@ -883,31 +920,33 @@ 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, m.primal.unwrap())).unzip(); - let (aggr_min, aggr_coeffs) = self.master.aggregate(i, &aggr_mins).map_err(SolverError::Master)?; + 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, - primal: Some( - self.problem - .aggregate_primals(aggr_coeffs.into_iter().zip(aggr_primals.into_iter()).collect()), - ), }); + 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> { let new_weight = self.weighter.weight(¤t_state!(self, Step::Descent), &self.params); - self.master.set_weight(new_weight).map_err(SolverError::Master)?; + 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); @@ -920,18 +959,18 @@ } /// Perform a null step. fn null_step(&mut self) -> Result<(), SolverError> { let new_weight = self.weighter.weight(¤t_state!(self, Step::Null), &self.params); - self.master.set_weight(new_weight).map_err(SolverError::Master)?; + self.master.set_weight(new_weight)?; self.cnt_null += 1; debug!("Null Step"); Ok(()) } /// Perform one bundle iteration. - #[cfg_attr(feature = "cargo-clippy", allow(collapsible_if))] + #[allow(clippy::collapsible_if)] pub fn step(&mut self) -> Result> { self.iterinfos.clear(); if !self.cur_valid { // current point needs new evaluation @@ -980,18 +1019,19 @@ 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: self - .master - .add_minorant(fidx, nxt_minorant) - .map_err(SolverError::Master)?, + index: minidx, multiplier: 0.0, - primal: Some(nxt_primal), }); + 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:{}", Index: src/vector.rs ================================================================== --- src/vector.rs +++ src/vector.rs @@ -14,13 +14,13 @@ // along with this program. If not, see // //! Finite-dimensional sparse and dense vectors. +use crate::Real; use std::fmt; use std::ops::{Deref, DerefMut}; -use Real; // use std::cmp::min; use std::iter::FromIterator; use std::vec::IntoIter; /// Type of dense vectors. @@ -41,18 +41,18 @@ } } impl fmt::Display for DVector { fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { - try!(write!(f, "(")); + write!(f, "(")?; for (i, x) in self.iter().enumerate() { if i > 0 { - try!(write!(f, ", ")); + write!(f, ", ")?; } - try!(write!(f, "{}", x)) + write!(f, "{}", x)? } - try!(write!(f, ")")); + write!(f, ")")?; Ok(()) } } impl FromIterator for DVector { @@ -89,15 +89,15 @@ 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)); + write!(f, "{}:(", size)?; if let Some(&(i, x)) = it.next() { - try!(write!(f, "{}:{}", i, x)); + write!(f, "{}:{}", i, x)?; for &(i, x) in it { - try!(write!(f, ", {}:{}", i, x)); + write!(f, ", {}:{}", i, x)?; } } write!(f, ")") } } @@ -189,11 +189,11 @@ for (x, y) in self.iter_mut().zip(y.iter()) { *x += alpha * y; } let n = self.len(); if n < y.len() { - self.extend_from_slice(&y[n..]); + 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() { @@ -283,5 +283,13 @@ 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]); +}