/*
 * 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 }, ())],
            })
        }
    }
}
