1.1 --- a/lemon/lp_base.h Fri Jun 03 12:25:23 2005 +0000
1.2 +++ b/lemon/lp_base.h Sat Jun 04 12:50:15 2005 +0000
1.3 @@ -410,6 +410,116 @@
1.4 }
1.5 };
1.6
1.7 + ///Linear expression of rows
1.8 +
1.9 + ///This data structure represents a column of the matrix,
1.10 + ///thas is it strores a linear expression of the dual variables
1.11 + ///(\ref Row "Row"s).
1.12 + ///
1.13 + ///There are several ways to access and modify the contents of this
1.14 + ///container.
1.15 + ///- Its it fully compatible with \c std::map<Row,double>, so for expamle
1.16 + ///if \c e is an DualExpr and \c v
1.17 + ///and \c w are of type \ref Row, then you can
1.18 + ///read and modify the coefficients like
1.19 + ///these.
1.20 + ///\code
1.21 + ///e[v]=5;
1.22 + ///e[v]+=12;
1.23 + ///e.erase(v);
1.24 + ///\endcode
1.25 + ///or you can also iterate through its elements.
1.26 + ///\code
1.27 + ///double s=0;
1.28 + ///for(LpSolverBase::DualExpr::iterator i=e.begin();i!=e.end();++i)
1.29 + /// s+=i->second;
1.30 + ///\endcode
1.31 + ///(This code computes the sum of all coefficients).
1.32 + ///- Numbers (<tt>double</tt>'s)
1.33 + ///and variables (\ref Row "Row"s) directly convert to an
1.34 + ///\ref DualExpr and the usual linear operations are defined so
1.35 + ///\code
1.36 + ///v+w
1.37 + ///2*v-3.12*(v-w/2)
1.38 + ///v*2.1+(3*v+(v*12+w)*3)/2
1.39 + ///\endcode
1.40 + ///are valid \ref DualExpr "DualExpr"essions.
1.41 + ///The usual assignment operations are also defined.
1.42 + ///\code
1.43 + ///e=v+w;
1.44 + ///e+=2*v-3.12*(v-w/2);
1.45 + ///e*=3.4;
1.46 + ///e/=5;
1.47 + ///\endcode
1.48 + ///
1.49 + ///\sa Expr
1.50 + ///
1.51 + class DualExpr : public std::map<Row,Value>
1.52 + {
1.53 + public:
1.54 + typedef LpSolverBase::Row Key;
1.55 + typedef LpSolverBase::Value Value;
1.56 +
1.57 + protected:
1.58 + typedef std::map<Row,Value> Base;
1.59 +
1.60 + public:
1.61 + typedef True IsLinExpression;
1.62 + ///\e
1.63 + DualExpr() : Base() { }
1.64 + ///\e
1.65 + DualExpr(const Key &v) {
1.66 + Base::insert(std::make_pair(v, 1));
1.67 + }
1.68 + ///\e
1.69 + DualExpr(const Value &v) {}
1.70 + ///\e
1.71 + void set(const Key &v,const Value &c) {
1.72 + Base::insert(std::make_pair(v, c));
1.73 + }
1.74 +
1.75 + ///Removes the components with zero coefficient.
1.76 + void simplify() {
1.77 + for (Base::iterator i=Base::begin(); i!=Base::end();) {
1.78 + Base::iterator j=i;
1.79 + ++j;
1.80 + if ((*i).second==0) Base::erase(i);
1.81 + j=i;
1.82 + }
1.83 + }
1.84 +
1.85 + ///Sets all coefficients to 0.
1.86 + void clear() {
1.87 + Base::clear();
1.88 + }
1.89 +
1.90 + ///\e
1.91 + DualExpr &operator+=(const DualExpr &e) {
1.92 + for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
1.93 + (*this)[j->first]+=j->second;
1.94 + ///\todo it might be speeded up using "hints"
1.95 + return *this;
1.96 + }
1.97 + ///\e
1.98 + DualExpr &operator-=(const DualExpr &e) {
1.99 + for (Base::const_iterator j=e.begin(); j!=e.end(); ++j)
1.100 + (*this)[j->first]-=j->second;
1.101 + return *this;
1.102 + }
1.103 + ///\e
1.104 + DualExpr &operator*=(const Value &c) {
1.105 + for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
1.106 + j->second*=c;
1.107 + return *this;
1.108 + }
1.109 + ///\e
1.110 + DualExpr &operator/=(const Value &c) {
1.111 + for (Base::iterator j=Base::begin(); j!=Base::end(); ++j)
1.112 + j->second/=c;
1.113 + return *this;
1.114 + }
1.115 + };
1.116 +
1.117
1.118 protected:
1.119 _FixId rows;
1.120 @@ -547,13 +657,113 @@
1.121 }
1.122 #endif
1.123
1.124 - ///Add a new empty row (i.e a new constaint) to the LP
1.125 + ///Set a column (i.e a dual constraint) of the LP
1.126
1.127 - ///This function adds a new empty row (i.e a new constaint) to the LP.
1.128 + ///\param c is the column to be modified
1.129 + ///\param e is a dual linear expression (see \ref DualExpr)
1.130 + ///\bug This is a temportary function. The interface will change to
1.131 + ///a better one.
1.132 + void setCol(Col c,const DualExpr &e) {
1.133 + std::vector<int> indices;
1.134 + std::vector<Value> values;
1.135 + indices.push_back(0);
1.136 + values.push_back(0);
1.137 + for(DualExpr::const_iterator i=e.begin(); i!=e.end(); ++i)
1.138 + if((*i).second!=0) { ///\bug EPSILON would be necessary here!!!
1.139 + indices.push_back(cols.floatingId((*i).first.id));
1.140 + values.push_back((*i).second);
1.141 + }
1.142 + _setColCoeffs(cols.floatingId(c.id),indices.size()-1,
1.143 + &indices[0],&values[0]);
1.144 + }
1.145 +
1.146 + ///Add a new column to the LP
1.147 +
1.148 + ///\param e is a dual linear expression (see \ref DualExpr)
1.149 + ///\param obj is the corresponding component of the objective
1.150 + ///function. It is 0 by default.
1.151 + ///\return The created column.
1.152 + ///\bug This is a temportary function. The interface will change to
1.153 + ///a better one.
1.154 + Col addCol(Value l,const DualExpr &e, Value obj=0) {
1.155 + Col c=addCol();
1.156 + setCol(c,e);
1.157 + objCoeff(c,0);
1.158 + return c;
1.159 + }
1.160 +
1.161 + ///Add a new empty row (i.e a new constraint) to the LP
1.162 +
1.163 + ///This function adds a new empty row (i.e a new constraint) to the LP.
1.164 ///\return The created row
1.165 Row addRow() { Row r; r.id=rows.insert(_addRow()); return r;}
1.166
1.167 - ///Set a row (i.e a constaint) of the LP
1.168 + ///\brief Adds several new row
1.169 + ///(i.e a variables) at once
1.170 + ///
1.171 + ///This magic function takes a container as its argument
1.172 + ///and fills its elements
1.173 + ///with new row (i.e. variables)
1.174 + ///\param t can be
1.175 + ///- a standard STL compatible iterable container with
1.176 + ///\ref Row as its \c values_type
1.177 + ///like
1.178 + ///\code
1.179 + ///std::vector<LpSolverBase::Row>
1.180 + ///std::list<LpSolverBase::Row>
1.181 + ///\endcode
1.182 + ///- a standard STL compatible iterable container with
1.183 + ///\ref Row as its \c mapped_type
1.184 + ///like
1.185 + ///\code
1.186 + ///std::map<AnyType,LpSolverBase::Row>
1.187 + ///\endcode
1.188 + ///- an iterable lemon \ref concept::WriteMap "write map" like
1.189 + ///\code
1.190 + ///ListGraph::NodeMap<LpSolverBase::Row>
1.191 + ///ListGraph::EdgeMap<LpSolverBase::Row>
1.192 + ///\endcode
1.193 + ///\return The number of rows created.
1.194 +#ifdef DOXYGEN
1.195 + template<class T>
1.196 + int addRowSet(T &t) { return 0;}
1.197 +#else
1.198 + template<class T>
1.199 + typename enable_if<typename T::value_type::LpSolverRow,int>::type
1.200 + addRowSet(T &t,dummy<0> = 0) {
1.201 + int s=0;
1.202 + for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
1.203 + return s;
1.204 + }
1.205 + template<class T>
1.206 + typename enable_if<typename T::value_type::second_type::LpSolverRow,
1.207 + int>::type
1.208 + addRowSet(T &t,dummy<1> = 1) {
1.209 + int s=0;
1.210 + for(typename T::iterator i=t.begin();i!=t.end();++i) {
1.211 + i->second=addRow();
1.212 + s++;
1.213 + }
1.214 + return s;
1.215 + }
1.216 + template<class T>
1.217 + typename enable_if<typename T::ValueSet::value_type::LpSolverRow,
1.218 + int>::type
1.219 + addRowSet(T &t,dummy<2> = 2) {
1.220 + ///\bug <tt>return addRowSet(t.valueSet());</tt> should also work.
1.221 + int s=0;
1.222 + for(typename T::ValueSet::iterator i=t.valueSet().begin();
1.223 + i!=t.valueSet().end();
1.224 + ++i)
1.225 + {
1.226 + *i=addRow();
1.227 + s++;
1.228 + }
1.229 + return s;
1.230 + }
1.231 +#endif
1.232 +
1.233 + ///Set a row (i.e a constraint) of the LP
1.234
1.235 ///\param r is the row to be modified
1.236 ///\param l is lower bound (-\ref INF means no bound)
1.237 @@ -580,7 +790,7 @@
1.238 _setRowBounds(rows.floatingId(r.id),l-e.constComp(),u-e.constComp());
1.239 }
1.240
1.241 - ///Set a row (i.e a constaint) of the LP
1.242 + ///Set a row (i.e a constraint) of the LP
1.243
1.244 ///\param r is the row to be modified
1.245 ///\param c is a linear expression (see \ref Constr)
1.246 @@ -591,7 +801,7 @@
1.247 c.upperBounded()?c.upperBound():INF);
1.248 }
1.249
1.250 - ///Add a new row (i.e a new constaint) to the LP
1.251 + ///Add a new row (i.e a new constraint) to the LP
1.252
1.253 ///\param l is the lower bound (-\ref INF means no bound)
1.254 ///\param e is a linear expression (see \ref Expr)
1.255 @@ -605,7 +815,7 @@
1.256 return r;
1.257 }
1.258
1.259 - ///Add a new row (i.e a new constaint) to the LP
1.260 + ///Add a new row (i.e a new constraint) to the LP
1.261
1.262 ///\param c is a linear expression (see \ref Constr)
1.263 ///\return The created row.
1.264 @@ -917,6 +1127,63 @@
1.265 return tmp;
1.266 }
1.267
1.268 + ///\e
1.269 +
1.270 + ///\relates LpSolverBase::DualExpr
1.271 + ///
1.272 + inline LpSolverBase::DualExpr operator+(const LpSolverBase::DualExpr &a,
1.273 + const LpSolverBase::DualExpr &b)
1.274 + {
1.275 + LpSolverBase::DualExpr tmp(a);
1.276 + tmp+=b; ///\todo Doesn't STL have some special 'merge' algorithm?
1.277 + return tmp;
1.278 + }
1.279 + ///\e
1.280 +
1.281 + ///\relates LpSolverBase::DualExpr
1.282 + ///
1.283 + inline LpSolverBase::DualExpr operator-(const LpSolverBase::DualExpr &a,
1.284 + const LpSolverBase::DualExpr &b)
1.285 + {
1.286 + LpSolverBase::DualExpr tmp(a);
1.287 + tmp-=b; ///\todo Doesn't STL have some special 'merge' algorithm?
1.288 + return tmp;
1.289 + }
1.290 + ///\e
1.291 +
1.292 + ///\relates LpSolverBase::DualExpr
1.293 + ///
1.294 + inline LpSolverBase::DualExpr operator*(const LpSolverBase::DualExpr &a,
1.295 + const LpSolverBase::Value &b)
1.296 + {
1.297 + LpSolverBase::DualExpr tmp(a);
1.298 + tmp*=b; ///\todo Doesn't STL have some special 'merge' algorithm?
1.299 + return tmp;
1.300 + }
1.301 +
1.302 + ///\e
1.303 +
1.304 + ///\relates LpSolverBase::DualExpr
1.305 + ///
1.306 + inline LpSolverBase::DualExpr operator*(const LpSolverBase::Value &a,
1.307 + const LpSolverBase::DualExpr &b)
1.308 + {
1.309 + LpSolverBase::DualExpr tmp(b);
1.310 + tmp*=a; ///\todo Doesn't STL have some special 'merge' algorithm?
1.311 + return tmp;
1.312 + }
1.313 + ///\e
1.314 +
1.315 + ///\relates LpSolverBase::DualExpr
1.316 + ///
1.317 + inline LpSolverBase::DualExpr operator/(const LpSolverBase::DualExpr &a,
1.318 + const LpSolverBase::Value &b)
1.319 + {
1.320 + LpSolverBase::DualExpr tmp(a);
1.321 + tmp/=b; ///\todo Doesn't STL have some special 'merge' algorithm?
1.322 + return tmp;
1.323 + }
1.324 +
1.325
1.326 } //namespace lemon
1.327
2.1 --- a/test/lp_test.cc Fri Jun 03 12:25:23 2005 +0000
2.2 +++ b/test/lp_test.cc Sat Jun 04 12:50:15 2005 +0000
2.3 @@ -34,108 +34,140 @@
2.4
2.5 lp.addColSet(z);
2.6
2.7 + {
2.8 + LP::Expr e,f,g;
2.9 + LP::Col p1,p2,p3,p4,p5;
2.10 + LP::Constr c;
2.11 +
2.12 + e[p1]=2;
2.13 + e.constComp()=12;
2.14 + e[p1]+=2;
2.15 + e.constComp()+=12;
2.16 + e[p1]-=2;
2.17 + e.constComp()-=12;
2.18 +
2.19 + e=2;
2.20 + e=2.2;
2.21 + e=p1;
2.22 + e=f;
2.23 +
2.24 + e+=2;
2.25 + e+=2.2;
2.26 + e+=p1;
2.27 + e+=f;
2.28 +
2.29 + e-=2;
2.30 + e-=2.2;
2.31 + e-=p1;
2.32 + e-=f;
2.33 +
2.34 + e*=2;
2.35 + e*=2.2;
2.36 + e/=2;
2.37 + e/=2.2;
2.38 +
2.39 + e=((p1+p2)+(p1-p2)+(p1+12)+(12+p1)+(p1-12)+(12-p1)+
2.40 + (f+12)+(12+f)+(p1+f)+(f+p1)+(f+g)+
2.41 + (f-12)+(12-f)+(p1-f)+(f-p1)+(f-g)+
2.42 + 2.2*f+f*2.2+f/2.2+
2.43 + 2*f+f*2+f/2+
2.44 + 2.2*p1+p1*2.2+p1/2.2+
2.45 + 2*p1+p1*2+p1/2
2.46 + );
2.47
2.48 - LP::Expr e,f,g;
2.49 - LP::Col p1,p2,p3,p4,p5;
2.50 - LP::Constr c;
2.51 +
2.52 + c = (e <= f );
2.53 + c = (e <= 2.2);
2.54 + c = (e <= 2 );
2.55 + c = (e <= p1 );
2.56 + c = (2.2<= f );
2.57 + c = (2 <= f );
2.58 + c = (p1 <= f );
2.59 + c = (p1 <= p2 );
2.60 + c = (p1 <= 2.2);
2.61 + c = (p1 <= 2 );
2.62 + c = (2.2<= p2 );
2.63 + c = (2 <= p2 );
2.64 +
2.65 + c = (e >= f );
2.66 + c = (e >= 2.2);
2.67 + c = (e >= 2 );
2.68 + c = (e >= p1 );
2.69 + c = (2.2>= f );
2.70 + c = (2 >= f );
2.71 + c = (p1 >= f );
2.72 + c = (p1 >= p2 );
2.73 + c = (p1 >= 2.2);
2.74 + c = (p1 >= 2 );
2.75 + c = (2.2>= p2 );
2.76 + c = (2 >= p2 );
2.77 +
2.78 + c = (e == f );
2.79 + c = (e == 2.2);
2.80 + c = (e == 2 );
2.81 + c = (e == p1 );
2.82 + c = (2.2== f );
2.83 + c = (2 == f );
2.84 + c = (p1 == f );
2.85 + //c = (p1 == p2 );
2.86 + c = (p1 == 2.2);
2.87 + c = (p1 == 2 );
2.88 + c = (2.2== p2 );
2.89 + c = (2 == p2 );
2.90 +
2.91 + c = (2 <= e <= 3);
2.92 + c = (2 <= p1<= 3);
2.93 +
2.94 + c = (2 >= e >= 3);
2.95 + c = (2 >= p1>= 3);
2.96 +
2.97 + e[x[3]]=2;
2.98 + e[x[3]]=4;
2.99 + e[x[3]]=1;
2.100 + e.constComp()=12;
2.101 +
2.102 + lp.addRow(LP::INF,e,23);
2.103 + lp.addRow(LP::INF,3.0*(x[1]+x[2]/2)-x[3],23);
2.104 + lp.addRow(LP::INF,3.0*(x[1]+x[2]*2-5*x[3]+12-x[4]/3)+2*x[4]-4,23);
2.105 +
2.106 + lp.addRow(x[1]+x[3]<=x[5]-3);
2.107 + lp.addRow(-7<=x[1]+x[3]-12<=3);
2.108 + lp.addRow(x[1]<=x[5]);
2.109 + }
2.110
2.111 - e[p1]=2;
2.112 - e.constComp()=12;
2.113 - e[p1]+=2;
2.114 - e.constComp()+=12;
2.115 - e[p1]-=2;
2.116 - e.constComp()-=12;
2.117 -
2.118 - e=2;
2.119 - e=2.2;
2.120 - e=p1;
2.121 - e=f;
2.122 -
2.123 - e+=2;
2.124 - e+=2.2;
2.125 - e+=p1;
2.126 - e+=f;
2.127 -
2.128 - e-=2;
2.129 - e-=2.2;
2.130 - e-=p1;
2.131 - e-=f;
2.132 -
2.133 - e*=2;
2.134 - e*=2.2;
2.135 - e/=2;
2.136 - e/=2.2;
2.137 -
2.138 - e=((p1+p2)+(p1-p2)+(p1+12)+(12+p1)+(p1-12)+(12-p1)+
2.139 - (f+12)+(12+f)+(p1+f)+(f+p1)+(f+g)+
2.140 - (f-12)+(12-f)+(p1-f)+(f-p1)+(f-g)+
2.141 - 2.2*f+f*2.2+f/2.2+
2.142 - 2*f+f*2+f/2+
2.143 - 2.2*p1+p1*2.2+p1/2.2+
2.144 - 2*p1+p1*2+p1/2
2.145 - );
2.146 + {
2.147 + LP::DualExpr e,f,g;
2.148 + LP::Row p1,p2,p3,p4,p5;
2.149 +
2.150 + e[p1]=2;
2.151 + e[p1]+=2;
2.152 + e[p1]-=2;
2.153 +
2.154 + e=p1;
2.155 + e=f;
2.156 +
2.157 + e+=p1;
2.158 + e+=f;
2.159 +
2.160 + e-=p1;
2.161 + e-=f;
2.162 +
2.163 + e*=2;
2.164 + e*=2.2;
2.165 + e/=2;
2.166 + e/=2.2;
2.167 +
2.168 + e=((p1+p2)+(p1-p2)+(p1+12)+(12+p1)+(p1-12)+(12-p1)+
2.169 + (p1+f)+(f+p1)+(f+g)+
2.170 + (p1-f)+(f-p1)+(f-g)+
2.171 + 2.2*f+f*2.2+f/2.2+
2.172 + 2*f+f*2+f/2+
2.173 + 2.2*p1+p1*2.2+p1/2.2+
2.174 + 2*p1+p1*2+p1/2
2.175 + );
2.176 + }
2.177
2.178
2.179 - c = (e <= f );
2.180 - c = (e <= 2.2);
2.181 - c = (e <= 2 );
2.182 - c = (e <= p1 );
2.183 - c = (2.2<= f );
2.184 - c = (2 <= f );
2.185 - c = (p1 <= f );
2.186 - c = (p1 <= p2 );
2.187 - c = (p1 <= 2.2);
2.188 - c = (p1 <= 2 );
2.189 - c = (2.2<= p2 );
2.190 - c = (2 <= p2 );
2.191 -
2.192 - c = (e >= f );
2.193 - c = (e >= 2.2);
2.194 - c = (e >= 2 );
2.195 - c = (e >= p1 );
2.196 - c = (2.2>= f );
2.197 - c = (2 >= f );
2.198 - c = (p1 >= f );
2.199 - c = (p1 >= p2 );
2.200 - c = (p1 >= 2.2);
2.201 - c = (p1 >= 2 );
2.202 - c = (2.2>= p2 );
2.203 - c = (2 >= p2 );
2.204 -
2.205 - c = (e == f );
2.206 - c = (e == 2.2);
2.207 - c = (e == 2 );
2.208 - c = (e == p1 );
2.209 - c = (2.2== f );
2.210 - c = (2 == f );
2.211 - c = (p1 == f );
2.212 - //c = (p1 == p2 );
2.213 - c = (p1 == 2.2);
2.214 - c = (p1 == 2 );
2.215 - c = (2.2== p2 );
2.216 - c = (2 == p2 );
2.217 -
2.218 - c = (2 <= e <= 3);
2.219 - c = (2 <= p1<= 3);
2.220 -
2.221 - c = (2 >= e >= 3);
2.222 - c = (2 >= p1>= 3);
2.223 -
2.224 - e[x[3]]=2;
2.225 - e[x[3]]=4;
2.226 - e[x[3]]=1;
2.227 - e.constComp()=12;
2.228 -
2.229 - lp.addRow(LP::INF,e,23);
2.230 - lp.addRow(LP::INF,3.0*(x[1]+x[2]/2)-x[3],23);
2.231 - lp.addRow(LP::INF,3.0*(x[1]+x[2]*2-5*x[3]+12-x[4]/3)+2*x[4]-4,23);
2.232 -
2.233 - lp.addRow(x[1]+x[3]<=x[5]-3);
2.234 - lp.addRow(-7<=x[1]+x[3]-12<=3);
2.235 - lp.addRow(x[1]<=x[5]);
2.236 -
2.237 -
2.238 -
2.239 }
2.240
2.241 int main()