1.1 --- a/lemon/Makefile.am	Wed Nov 29 17:34:29 2006 +0000
     1.2 +++ b/lemon/Makefile.am	Wed Nov 29 17:35:31 2006 +0000
     1.3 @@ -76,6 +76,7 @@
     1.4  	lemon/lp_glpk.h \
     1.5  	lemon/lp_skeleton.h \
     1.6  	lemon/lp_soplex.h \
     1.7 +	lemon/lp_utils.h \
     1.8  	lemon/map_iterator.h \
     1.9  	lemon/maps.h \
    1.10  	lemon/matrix_maps.h \

     2.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     2.2 +++ b/lemon/lp_utils.h	Wed Nov 29 17:35:31 2006 +0000
     2.3 @@ -0,0 +1,364 @@
     2.4 +/* -*- C++ -*-
     2.5 + *
     2.6 + * This file is a part of LEMON, a generic C++ optimization library
     2.7 + *
     2.8 + * Copyright (C) 2003-2006
     2.9 + * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
    2.10 + * (Egervary Research Group on Combinatorial Optimization, EGRES).
    2.11 + *
    2.12 + * Permission to use, modify and distribute this software is granted
    2.13 + * provided that this copyright notice appears in all copies. For
    2.14 + * precise terms see the accompanying LICENSE file.
    2.15 + *
    2.16 + * This software is provided "AS IS" with no warranty of any kind,
    2.17 + * express or implied, and with no claim as to its suitability for any
    2.18 + * purpose.
    2.19 + *
    2.20 + */
    2.21 +
    2.22 +#ifndef LEMON_LP_UTILS_H
    2.23 +#define LEMON_LP_UTILS_H
    2.24 +
    2.25 +#include <lemon/lp_base.h>
    2.26 +
    2.27 +#include <lemon/lemon_reader.h>
    2.28 +
    2.29 +namespace lemon {
    2.30 +
    2.31 +  class LpReader : public LemonReader::SectionReader {
    2.32 +    typedef LemonReader::SectionReader Parent;
    2.33 +  public:
    2.34 +
    2.35 +
    2.36 +    /// \brief Constructor.
    2.37 +    ///
    2.38 +    /// Constructor for LpReader. It creates the LpReader and attach
    2.39 +    /// it into the given LemonReader. The lp reader will add
    2.40 +    /// variables, constraints and objective function to the
    2.41 +    /// given lp solver.
    2.42 +    LpReader(LemonReader& _reader, LpSolverBase& _lp, 
    2.43 +             const std::string& _name = std::string())
    2.44 +      : Parent(_reader), lp(_lp), name(_name) {} 
    2.45 +
    2.46 +
    2.47 +    /// \brief Destructor.
    2.48 +    ///
    2.49 +    /// Destructor for NodeSetReader.
    2.50 +    virtual ~LpReader() {}
    2.51 +
    2.52 +  private:
    2.53 +    LpReader(const LpReader&);
    2.54 +    void operator=(const LpReader&);
    2.55 +  
    2.56 +  protected:
    2.57 +
    2.58 +    /// \brief Gives back true when the SectionReader can process 
    2.59 +    /// the section with the given header line.
    2.60 +    ///
    2.61 +    /// It gives back true when the header line starts with \c \@lp,
    2.62 +    /// and the header line's name and the nodeset's name are the same.
    2.63 +    virtual bool header(const std::string& line) {
    2.64 +      std::istringstream ls(line);
    2.65 +      std::string command;
    2.66 +      std::string id;
    2.67 +      ls >> command >> id;
    2.68 +      return command == "@lp" && name == id;
    2.69 +    }
    2.70 +
    2.71 +  private:
    2.72 +
    2.73 +    enum Relation {
    2.74 +      LE, EQ, GE
    2.75 +    };
    2.76 +
    2.77 +    std::istream& readConstraint(std::istream& is, LpSolverBase::Constr& c) {
    2.78 +      char x;
    2.79 +      LpSolverBase::Expr e1, e2;
    2.80 +      Relation op1;
    2.81 +      is >> std::ws;
    2.82 +      readExpression(is, e1);
    2.83 +      is >> std::ws;
    2.84 +      readRelation(is, op1);
    2.85 +      is >> std::ws;
    2.86 +      readExpression(is, e2);
    2.87 +      if (!is.get(x)) {
    2.88 +        if (op1 == LE) {
    2.89 +          c = e1 <= e2;
    2.90 +        } else if (op1 == GE) {
    2.91 +          c = e1 >= e2;
    2.92 +        } else {
    2.93 +          c = e1 == e2;
    2.94 +        }
    2.95 +      } else {
    2.96 +        is.putback(x);
    2.97 +        LpSolverBase::Expr e3;
    2.98 +        Relation op2;
    2.99 +        readRelation(is, op2);
   2.100 +        is >> std::ws;
   2.101 +        readExpression(is, e3);
   2.102 +        if (!e1.empty() || !e3.empty()) {
   2.103 +          throw DataFormatError("Wrong range format");
   2.104 +        }
   2.105 +        if (op2 != op1 || op1 == EQ) {
   2.106 +          throw DataFormatError("Wrong range format");
   2.107 +        }
   2.108 +        if (op1 == LE) {
   2.109 +          c = e1.constComp() <= e2 <= e3.constComp();
   2.110 +        } else {
   2.111 +          c = e1.constComp() >= e2 >= e3.constComp();
   2.112 +        }
   2.113 +      }
   2.114 +    }
   2.115 +
   2.116 +    std::istream& readExpression(std::istream& is, LpSolverBase::Expr& e) {
   2.117 +      LpSolverBase::Col c;
   2.118 +      double d;
   2.119 +      char x;
   2.120 +      readElement(is, c, d);
   2.121 +      if (c != INVALID) {
   2.122 +        e += d * c;
   2.123 +      } else {
   2.124 +        e += d;
   2.125 +      }
   2.126 +      is >> std::ws;
   2.127 +      while (is.get(x) && (x == '+' || x == '-')) {
   2.128 +        is >> std::ws;
   2.129 +        readElement(is, c, d);
   2.130 +        if (c != INVALID) {
   2.131 +          e += (x == '+' ? d : -d) * c;
   2.132 +        } else {
   2.133 +          e += (x == '+' ? d : -d);
   2.134 +        }
   2.135 +        is >> std::ws;
   2.136 +      }
   2.137 +      if (!is) {
   2.138 +        is.clear();
   2.139 +      } else {
   2.140 +        is.putback(x);
   2.141 +      }
   2.142 +      return is;
   2.143 +    }
   2.144 +
   2.145 +    std::istream& readElement(std::istream& is, 
   2.146 +                              LpSolverBase::Col& c, double& d) { 
   2.147 +      d = 1.0;
   2.148 +      c = INVALID;
   2.149 +      char x, y;
   2.150 +      if (!is.get(x)) throw DataFormatError("Cannot find lp element");
   2.151 +      if (x == '+' || x == '-') {
   2.152 +        is >> std::ws;
   2.153 +        d *= x == '-' ? -1 : 1;
   2.154 +        while (is.get(x) && (x == '+' || x == '-')) {
   2.155 +          d *= x == '-' ? -1 : 1;
   2.156 +          is >> std::ws;
   2.157 +        }
   2.158 +        if (!is) throw DataFormatError("Cannot find lp element");
   2.159 +      }
   2.160 +      if (numFirstChar(x)) {
   2.161 +        is.putback(x);
   2.162 +        double e;
   2.163 +        readNum(is, e);
   2.164 +        d *= e;
   2.165 +      } else if (varFirstChar(x)) {
   2.166 +        is.putback(x);
   2.167 +        LpSolverBase::Col f;
   2.168 +        readCol(is, f);
   2.169 +        c = f;
   2.170 +      } else {
   2.171 +        throw DataFormatError("Invalid expression format");          
   2.172 +      }
   2.173 +      is >> std::ws;
   2.174 +      while (is.get(y) && (y == '*' || y == '/')) {
   2.175 +        is >> std::ws;
   2.176 +        if (!is.get(x)) throw DataFormatError("Cannot find lp element");
   2.177 +        if (x == '+' || x == '-') {
   2.178 +          is >> std::ws;
   2.179 +          d *= x == '-' ? -1 : 1;
   2.180 +          while (is.get(x) && (x == '+' || x == '-')) {
   2.181 +            d *= x == '-' ? -1 : 1;
   2.182 +            is >> std::ws;
   2.183 +          }
   2.184 +          if (!is) throw DataFormatError("Cannot find lp element");
   2.185 +        }
   2.186 +        if (numFirstChar(x)) {
   2.187 +          is.putback(x);
   2.188 +          double e;
   2.189 +          readNum(is, e);
   2.190 +          if (y == '*') {
   2.191 +            d *= e;
   2.192 +          } else {
   2.193 +            d /= e;
   2.194 +          }
   2.195 +        } else if (varFirstChar(x)) {
   2.196 +          is.putback(x);
   2.197 +          LpSolverBase::Col f;
   2.198 +          readCol(is, f);
   2.199 +          if (y == '*') {
   2.200 +            if (c == INVALID) {
   2.201 +              c = f;
   2.202 +            } else {
   2.203 +              throw DataFormatError("Quadratic element in expression");
   2.204 +            }
   2.205 +          } else {
   2.206 +            throw DataFormatError("Division by variable");
   2.207 +          }
   2.208 +        } else {
   2.209 +          throw DataFormatError("Invalid expression format");          
   2.210 +        }
   2.211 +        is >> std::ws;
   2.212 +      }
   2.213 +      if (!is) {
   2.214 +        is.clear();
   2.215 +      } else {
   2.216 +        is.putback(y);
   2.217 +      }
   2.218 +      return is;
   2.219 +    }
   2.220 +
   2.221 +    std::istream& readCol(std::istream& is, LpSolverBase::Col& c) {
   2.222 +      char x;
   2.223 +      std::string var;
   2.224 +      while (is.get(x) && varChar(x)) {
   2.225 +        var += x;
   2.226 +      }
   2.227 +      if (!is) {
   2.228 +        is.clear();
   2.229 +      } else {
   2.230 +        is.putback(x);
   2.231 +      }
   2.232 +      ColMap::const_iterator it = cols.find(var);
   2.233 +      if (cols.find(var) != cols.end()) {
   2.234 +        c = it->second;
   2.235 +      } else {
   2.236 +        c = lp.addCol();
   2.237 +        cols.insert(std::make_pair(var, c));
   2.238 +        lp.colName(c, var);
   2.239 +      }
   2.240 +      return is;
   2.241 +    }
   2.242 +
   2.243 +    std::istream& readNum(std::istream& is, double& d) {
   2.244 +      is >> d;
   2.245 +      if (!is) throw DataFormatError("Wrong number format");
   2.246 +      return is;
   2.247 +    }
   2.248 +
   2.249 +    std::istream& readRelation(std::istream& is, Relation& op) {
   2.250 +      char x, y;
   2.251 +      if (!is.get(x) || !(x == '<' || x == '=' || x == '>')) {
   2.252 +        throw DataFormatError("Wrong relation operator");
   2.253 +      }
   2.254 +      if (!is.get(y) || y != '=') {
   2.255 +        throw DataFormatError("Wrong relation operator");
   2.256 +      }
   2.257 +      switch (x) {
   2.258 +      case '<': op = LE; 
   2.259 +        break;
   2.260 +      case '=': op = EQ; 
   2.261 +        break;
   2.262 +      case '>': op = GE; 
   2.263 +        break;
   2.264 +      }
   2.265 +      return is;
   2.266 +    }
   2.267 +
   2.268 +    static bool relationFirstChar(char c) {
   2.269 +      return c == '<' || c == '=' || c == '>';
   2.270 +    }
   2.271 +
   2.272 +    static bool varFirstChar(char c) {
   2.273 +      return (c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z');
   2.274 +    }
   2.275 +
   2.276 +    static bool numFirstChar(char c) {
   2.277 +      return (c >= '0' && c <= '9') || c == '.';
   2.278 +    }
   2.279 +
   2.280 +    static bool varChar(char c) {
   2.281 +      return (c >= '0' && c <= '9') || 
   2.282 +        (c >= 'a' && c <= 'z') || (c >= 'A' && c <= 'Z');
   2.283 +    }
   2.284 +
   2.285 +
   2.286 +    void addConstraint(const LpSolverBase::Constr& constr) {
   2.287 +      if (constr.expr().size() != 1) {
   2.288 +        lp.addRow(constr);
   2.289 +      } else {
   2.290 +        Lp::Expr e = constr.expr();
   2.291 +        LpSolverBase::Col col = e.begin()->first;
   2.292 +        double coeff = e.begin()->second;
   2.293 +        double lb = LpSolverBase::NaN;
   2.294 +        double ub = LpSolverBase::NaN;
   2.295 +        if (coeff > 0) {
   2.296 +          if (constr.upperBounded()) {
   2.297 +            lb = (constr.lowerBound() - e.constComp()) / coeff;
   2.298 +          }
   2.299 +          if (constr.lowerBounded()) {
   2.300 +            ub = (constr.upperBound() - e.constComp()) / coeff;
   2.301 +          }
   2.302 +        } else if (coeff < 0) {
   2.303 +          if (constr.upperBounded()) {
   2.304 +            lb = (constr.upperBound() - e.constComp()) / coeff;
   2.305 +          }
   2.306 +          if (constr.lowerBounded()) {
   2.307 +            ub = (constr.lowerBound() - e.constComp()) / coeff;
   2.308 +          }
   2.309 +        } else {
   2.310 +          lb = -LpSolverBase::INF;
   2.311 +          ub = LpSolverBase::INF;
   2.312 +        }
   2.313 +        lp.colLowerBound(col, lb);
   2.314 +        lp.colUpperBound(col, ub);
   2.315 +      }
   2.316 +    }
   2.317 +    
   2.318 +  protected:
   2.319 +
   2.320 +    /// \brief Reader function of the section.
   2.321 +    ///
   2.322 +    /// It reads the content of the section.
   2.323 +    virtual void read(std::istream& is) {
   2.324 +      std::string line;
   2.325 +      std::map<std::string, LpSolverBase::Col> vars;
   2.326 +      while (getline(is, line)) {
   2.327 +        std::istringstream ls(line);
   2.328 +        std::string sense;
   2.329 +        ls >> sense;
   2.330 +        if (sense == "min" || sense == "max") {
   2.331 +          LpSolverBase::Expr expr;
   2.332 +          ls >> std::ws;
   2.333 +          readExpression(ls, expr);
   2.334 +          lp.setObj(expr);
   2.335 +          if (sense == "min") {
   2.336 +            lp.min();
   2.337 +          } else {
   2.338 +            lp.max();
   2.339 +          }
   2.340 +        } else {
   2.341 +          ls.str(line);
   2.342 +          LpSolverBase::Constr constr;
   2.343 +          ls >> std::ws;
   2.344 +          readConstraint(ls, constr);
   2.345 +          addConstraint(constr);
   2.346 +        }        
   2.347 +      }
   2.348 +    }
   2.349 +      
   2.350 +    virtual void missing() {
   2.351 +      ErrorMessage msg;
   2.352 +      msg << "Lp section not found in file: @lp " << name;
   2.353 +      throw IoParameterError(msg.message());
   2.354 +    }
   2.355 +
   2.356 +  private:
   2.357 +
   2.358 +    typedef std::map<std::string, LpSolverBase::Col> ColMap;
   2.359 +      
   2.360 +    LpSolverBase& lp;
   2.361 +    std::string name;
   2.362 +    ColMap cols;
   2.363 +  };
   2.364 +
   2.365 +}
   2.366 +
   2.367 +#endif