/*
* 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, Vector, DVector, Minorant, FirstOrderProblem, SimpleEvaluation};
use mcf;
use std::fs::File;
use std::io::{self, Read};
use std::result;
use std::num::{ParseIntError, ParseFloatError};
use std::f64::INFINITY;
quick_error! {
/// A solver error.
#[derive(Debug)]
pub enum Error {
Io(err: io::Error) {
cause(err)
description(err.description())
display("Io Error: {}", err)
from()
}
Solver(err: mcf::solver::Error) {
cause(err)
description(err.description())
display("Solver Error: {}", err)
from()
}
Format(msg: String) {
description("Format error")
display("Format error: {}", msg)
from(err: ParseIntError) -> (format!("{}", err))
from(err: ParseFloatError) -> (format!("{}", err))
}
}
}
pub type Result<T> = result::Result<T, Error>;
#[derive(Clone, Copy, Debug)]
struct ArcInfo { arc : usize, src : usize, snk : usize }
#[derive(Clone, Copy, Debug)]
struct Elem { ind : usize, val : Real }
pub struct MMCFProblem {
pub multimodel : bool,
nets : Vec<mcf::Solver>,
lhs : Vec<Vec<Vec<Elem>>>,
rhs : DVector,
rhsval : Real,
cbase : Vec<DVector>,
c : Vec<DVector>,
}
impl MMCFProblem {
pub fn read_mnetgen(basename: &str) -> Result<MMCFProblem> {
let mut buffer = String::new();
{
let mut f = try!(File::open(&format!("{}.nod", basename)));
try!(f.read_to_string(&mut buffer));
}
let fnod = buffer.split_whitespace().map(|x| x.parse::<usize>().unwrap()).collect::<Vec<_>>();
if fnod.len() != 4 {
return Err(Error::Format(format!("Expected 4 numbers in {}.nod, but got {}", basename, fnod.len())));
}
let ncom = fnod[0];
let nnodes = fnod[1];
let narcs = fnod[2];
let ncaps = fnod[3];
// read nodes
let mut nets = Vec::with_capacity(ncom);
for _ in 0..ncom { nets.push(try!(mcf::Solver::new(nnodes)))};
{
let mut f = try!(File::open(&format!("{}.sup", basename)));
buffer.clear();
try!(f.read_to_string(&mut buffer));
}
for line in buffer.lines() {
let mut data = line.split_whitespace();
let node = try!(data.next().unwrap().parse::<usize>());
let com = try!(data.next().unwrap().parse::<usize>());
let supply = try!(data.next().unwrap().parse::<Real>());
try!(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)));
buffer.clear();
try!(f.read_to_string(&mut buffer));
}
for line in buffer.lines() {
let mut data = line.split_whitespace();
let arc = try!(data.next().unwrap().parse::<usize>())-1;
let src = try!(data.next().unwrap().parse::<usize>())-1;
let snk = try!(data.next().unwrap().parse::<usize>())-1;
let com = try!(data.next().unwrap().parse::<usize>())-1;
let cost = try!(data.next().unwrap().parse::<Real>());
let cap = try!(data.next().unwrap().parse::<Real>());
let mt = try!(data.next().unwrap().parse::<isize>())-1;
// set internal coeff
let coeff = arcmap[com].len();
arcmap[com].push(ArcInfo { 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 }));
// set objective
cbase[com].push(cost); // + 1e-6 * coeff
// add to mutual capacity constraint
if mt >= 0 {
lhsidx[mt as usize][com].push(coeff);
}
}
// read rhs of coupling constraints
{
let mut f = try!(File::open(&format!("{}.mut", basename)));
buffer.clear();
try!(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::<usize>())-1;
let cap = try!(data.next().unwrap().parse::<Real>());
rhs[mt] = cap;
}
// set lhs
let mut lhs = vec![vec![vec![]; ncom]; ncaps];
for i in 0..ncaps {
for fidx in 0..ncom {
lhs[i][fidx] = lhsidx[i][fidx].iter().map(|&j| Elem{ ind: j, val: 1.0 }).collect();
}
}
Ok(MMCFProblem {
multimodel : false,
nets: nets,
lhs: lhs,
rhs: rhs,
rhsval : 0.0,
cbase: cbase,
c: vec![dvec![]; ncom],
})
}
}
impl<'a> FirstOrderProblem<'a> for MMCFProblem {
type Error = Error;
type Primal = ();
type EvalResult = SimpleEvaluation<()>;
fn num_variables(&self) -> usize { self.lhs.len() }
fn lower_bounds(&self) -> Option<Vector> {
Some(Vector::new_sparse(self.lhs.len(), &[], &[]))
}
fn upper_bounds(&self) -> Option<Vector> { None }
fn num_subproblems(&self) -> usize {
if self.multimodel {self.nets.len()} else {1}
}
#[allow(unused_variables)]
fn evaluate(&'a mut self, fidx : usize, y : &DVector,
nullstep_bound : Real,
relprec : Real) -> result::Result<Self::EvalResult, Self::Error>
{
// compute costs
self.rhsval = 0.0;
for i in 0..self.c.len() {
self.c[i].clear();
self.c[i].extend(self.cbase[i].iter());
}
for i in 0..self.lhs.len() {
if y[i] != 0.0 {
self.rhsval += self.rhs[i] * y[i];
for fidx in 0..self.nets.len() {
for elem in &self.lhs[i][fidx] {
self.c[fidx][elem.ind] += y[i] * elem.val;
}
}
}
}
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]));
}
// solve subproblems
for (i,net) in self.nets.iter_mut().enumerate() {
try!(net.solve());
debug!("c[{}]={}", i, try!(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());
} else {
subg = dvec![0.0; self.rhs.len()];
objective = -try!(self.nets[fidx].objective());
}
let sol = try!(self.nets[fidx].get_solution());
for i in 0..self.lhs.len() {
for elem in &self.lhs[i][fidx] {
subg[i] -= elem.val * sol[elem.ind];
}
}
Ok(SimpleEvaluation {
objective: objective,
minorants: vec![(Minorant { constant: objective, linear: subg }, ())],
})
} 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()));
}
let mut subg = self.rhs.clone();
for i in 0..self.lhs.len() {
for (fidx, flhs) in self.lhs[i].iter().enumerate() {
for elem in flhs {
subg[i] -= elem.val * sols[fidx][elem.ind];
}
}
}
Ok(SimpleEvaluation {
objective: objective,
minorants: vec![(Minorant { constant: objective, linear: subg }, ())],
})
}
}
}