gravatar
alpar (Alpar Juttner)
alpar@cs.elte.hu
Merge
0 5 0
merge default
0 files changed with 37 insertions and 14 deletions:
↑ Collapse diff ↑
Ignore white space 16777216 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_GRAPH_TO_EPS_H
20 20
#define LEMON_GRAPH_TO_EPS_H
21 21

	
22 22
#include<iostream>
23 23
#include<fstream>
24 24
#include<sstream>
25 25
#include<algorithm>
26 26
#include<vector>
27 27

	
28 28
#ifndef WIN32
29 29
#include<sys/time.h>
30 30
#include<ctime>
31 31
#else
32
#ifndef WIN32_LEAN_AND_MEAN
32 33
#define WIN32_LEAN_AND_MEAN
34
#endif
35
#ifndef NOMINMAX
33 36
#define NOMINMAX
37
#endif
34 38
#include<windows.h>
35 39
#endif
36 40

	
37 41
#include<lemon/math.h>
38 42
#include<lemon/core.h>
39 43
#include<lemon/dim2.h>
40 44
#include<lemon/maps.h>
41 45
#include<lemon/color.h>
42 46
#include<lemon/bits/bezier.h>
43 47
#include<lemon/error.h>
44 48

	
45 49

	
46 50
///\ingroup eps_io
47 51
///\file
48 52
///\brief A well configurable tool for visualizing graphs
49 53

	
50 54
namespace lemon {
51 55

	
52 56
  namespace _graph_to_eps_bits {
53 57
    template<class MT>
54 58
    class _NegY {
55 59
    public:
56 60
      typedef typename MT::Key Key;
57 61
      typedef typename MT::Value Value;
58 62
      const MT &map;
59 63
      int yscale;
60 64
      _NegY(const MT &m,bool b) : map(m), yscale(1-b*2) {}
61 65
      Value operator[](Key n) { return Value(map[n].x,map[n].y*yscale);}
62 66
    };
63 67
  }
64 68

	
65 69
///Default traits class of GraphToEps
66 70

	
67 71
///Default traits class of \ref GraphToEps.
68 72
///
69 73
///\c G is the type of the underlying graph.
70 74
template<class G>
71 75
struct DefaultGraphToEpsTraits
72 76
{
73 77
  typedef G Graph;
74 78
  typedef typename Graph::Node Node;
75 79
  typedef typename Graph::NodeIt NodeIt;
76 80
  typedef typename Graph::Arc Arc;
77 81
  typedef typename Graph::ArcIt ArcIt;
78 82
  typedef typename Graph::InArcIt InArcIt;
79 83
  typedef typename Graph::OutArcIt OutArcIt;
80 84

	
81 85

	
82 86
  const Graph &g;
83 87

	
84 88
  std::ostream& os;
85 89

	
86 90
  typedef ConstMap<typename Graph::Node,dim2::Point<double> > CoordsMapType;
87 91
  CoordsMapType _coords;
88 92
  ConstMap<typename Graph::Node,double > _nodeSizes;
89 93
  ConstMap<typename Graph::Node,int > _nodeShapes;
90 94

	
91 95
  ConstMap<typename Graph::Node,Color > _nodeColors;
92 96
  ConstMap<typename Graph::Arc,Color > _arcColors;
93 97

	
94 98
  ConstMap<typename Graph::Arc,double > _arcWidths;
95 99

	
96 100
  double _arcWidthScale;
97 101

	
98 102
  double _nodeScale;
99 103
  double _xBorder, _yBorder;
100 104
  double _scale;
101 105
  double _nodeBorderQuotient;
102 106

	
103 107
  bool _drawArrows;
104 108
  double _arrowLength, _arrowWidth;
105 109

	
106 110
  bool _showNodes, _showArcs;
107 111

	
108 112
  bool _enableParallel;
109 113
  double _parArcDist;
110 114

	
111 115
  bool _showNodeText;
112 116
  ConstMap<typename Graph::Node,bool > _nodeTexts;
113 117
  double _nodeTextSize;
114 118

	
115 119
  bool _showNodePsText;
116 120
  ConstMap<typename Graph::Node,bool > _nodePsTexts;
117 121
  char *_nodePsTextsPreamble;
118 122

	
119 123
  bool _undirected;
120 124

	
121 125
  bool _pleaseRemoveOsStream;
122 126

	
123 127
  bool _scaleToA4;
124 128

	
125 129
  std::string _title;
126 130
  std::string _copyright;
127 131

	
128 132
  enum NodeTextColorType
129 133
    { DIST_COL=0, DIST_BW=1, CUST_COL=2, SAME_COL=3 } _nodeTextColorType;
130 134
  ConstMap<typename Graph::Node,Color > _nodeTextColors;
131 135

	
132 136
  bool _autoNodeScale;
133 137
  bool _autoArcWidthScale;
134 138

	
135 139
  bool _absoluteNodeSizes;
136 140
  bool _absoluteArcWidths;
137 141

	
138 142
  bool _negY;
139 143

	
140 144
  bool _preScale;
141 145
  ///Constructor
142 146

	
143 147
  ///Constructor
144 148
  ///\param _g  Reference to the graph to be printed.
145 149
  ///\param _os Reference to the output stream.
146 150
  ///\param _os Reference to the output stream.
147 151
  ///By default it is <tt>std::cout</tt>.
148 152
  ///\param _pros If it is \c true, then the \c ostream referenced by \c _os
149 153
  ///will be explicitly deallocated by the destructor.
150 154
  DefaultGraphToEpsTraits(const G &_g,std::ostream& _os=std::cout,
151 155
                          bool _pros=false) :
152 156
    g(_g), os(_os),
153 157
    _coords(dim2::Point<double>(1,1)), _nodeSizes(1), _nodeShapes(0),
154 158
    _nodeColors(WHITE), _arcColors(BLACK),
155 159
    _arcWidths(1.0), _arcWidthScale(0.003),
156 160
    _nodeScale(.01), _xBorder(10), _yBorder(10), _scale(1.0),
157 161
    _nodeBorderQuotient(.1),
158 162
    _drawArrows(false), _arrowLength(1), _arrowWidth(0.3),
159 163
    _showNodes(true), _showArcs(true),
160 164
    _enableParallel(false), _parArcDist(1),
161 165
    _showNodeText(false), _nodeTexts(false), _nodeTextSize(1),
162 166
    _showNodePsText(false), _nodePsTexts(false), _nodePsTextsPreamble(0),
163 167
    _undirected(lemon::UndirectedTagIndicator<G>::value),
164 168
    _pleaseRemoveOsStream(_pros), _scaleToA4(false),
165 169
    _nodeTextColorType(SAME_COL), _nodeTextColors(BLACK),
166 170
    _autoNodeScale(false),
167 171
    _autoArcWidthScale(false),
168 172
    _absoluteNodeSizes(false),
169 173
    _absoluteArcWidths(false),
170 174
    _negY(false),
171 175
    _preScale(true)
172 176
  {}
173 177
};
174 178

	
175 179
///Auxiliary class to implement the named parameters of \ref graphToEps()
176 180

	
177 181
///Auxiliary class to implement the named parameters of \ref graphToEps().
178 182
///
179 183
///For detailed examples see the \ref graph_to_eps_demo.cc demo file.
180 184
template<class T> class GraphToEps : public T
181 185
{
182 186
  // Can't believe it is required by the C++ standard
183 187
  using T::g;
184 188
  using T::os;
185 189

	
186 190
  using T::_coords;
187 191
  using T::_nodeSizes;
188 192
  using T::_nodeShapes;
189 193
  using T::_nodeColors;
190 194
  using T::_arcColors;
191 195
  using T::_arcWidths;
192 196

	
193 197
  using T::_arcWidthScale;
194 198
  using T::_nodeScale;
195 199
  using T::_xBorder;
196 200
  using T::_yBorder;
197 201
  using T::_scale;
198 202
  using T::_nodeBorderQuotient;
199 203

	
200 204
  using T::_drawArrows;
201 205
  using T::_arrowLength;
202 206
  using T::_arrowWidth;
203 207

	
204 208
  using T::_showNodes;
205 209
  using T::_showArcs;
206 210

	
207 211
  using T::_enableParallel;
208 212
  using T::_parArcDist;
209 213

	
210 214
  using T::_showNodeText;
211 215
  using T::_nodeTexts;
212 216
  using T::_nodeTextSize;
213 217

	
214 218
  using T::_showNodePsText;
215 219
  using T::_nodePsTexts;
216 220
  using T::_nodePsTextsPreamble;
217 221

	
218 222
  using T::_undirected;
219 223

	
220 224
  using T::_pleaseRemoveOsStream;
221 225

	
222 226
  using T::_scaleToA4;
223 227

	
224 228
  using T::_title;
225 229
  using T::_copyright;
226 230

	
227 231
  using T::NodeTextColorType;
228 232
  using T::CUST_COL;
229 233
  using T::DIST_COL;
230 234
  using T::DIST_BW;
231 235
  using T::_nodeTextColorType;
232 236
  using T::_nodeTextColors;
233 237

	
234 238
  using T::_autoNodeScale;
235 239
  using T::_autoArcWidthScale;
236 240

	
237 241
  using T::_absoluteNodeSizes;
238 242
  using T::_absoluteArcWidths;
239 243

	
240 244

	
241 245
  using T::_negY;
242 246
  using T::_preScale;
243 247

	
244 248
  // dradnats ++C eht yb deriuqer si ti eveileb t'naC
245 249

	
246 250
  typedef typename T::Graph Graph;
247 251
  typedef typename Graph::Node Node;
248 252
  typedef typename Graph::NodeIt NodeIt;
249 253
  typedef typename Graph::Arc Arc;
250 254
  typedef typename Graph::ArcIt ArcIt;
251 255
  typedef typename Graph::InArcIt InArcIt;
252 256
  typedef typename Graph::OutArcIt OutArcIt;
253 257

	
254 258
  static const int INTERPOL_PREC;
255 259
  static const double A4HEIGHT;
256 260
  static const double A4WIDTH;
257 261
  static const double A4BORDER;
258 262

	
259 263
  bool dontPrint;
260 264

	
261 265
public:
262 266
  ///Node shapes
263 267

	
264 268
  ///Node shapes.
265 269
  ///
266 270
  enum NodeShapes {
267 271
    /// = 0
268 272
    ///\image html nodeshape_0.png
269 273
    ///\image latex nodeshape_0.eps "CIRCLE shape (0)" width=2cm
270 274
    CIRCLE=0,
271 275
    /// = 1
272 276
    ///\image html nodeshape_1.png
273 277
    ///\image latex nodeshape_1.eps "SQUARE shape (1)" width=2cm
274 278
    ///
275 279
    SQUARE=1,
276 280
    /// = 2
277 281
    ///\image html nodeshape_2.png
278 282
    ///\image latex nodeshape_2.eps "DIAMOND shape (2)" width=2cm
279 283
    ///
280 284
    DIAMOND=2,
281 285
    /// = 3
282 286
    ///\image html nodeshape_3.png
283 287
    ///\image latex nodeshape_2.eps "MALE shape (4)" width=2cm
284 288
    ///
285 289
    MALE=3,
286 290
    /// = 4
287 291
    ///\image html nodeshape_4.png
288 292
    ///\image latex nodeshape_2.eps "FEMALE shape (4)" width=2cm
289 293
    ///
290 294
    FEMALE=4
291 295
  };
292 296

	
293 297
private:
294 298
  class arcLess {
295 299
    const Graph &g;
296 300
  public:
297 301
    arcLess(const Graph &_g) : g(_g) {}
298 302
    bool operator()(Arc a,Arc b) const
299 303
    {
300 304
      Node ai=std::min(g.source(a),g.target(a));
301 305
      Node aa=std::max(g.source(a),g.target(a));
302 306
      Node bi=std::min(g.source(b),g.target(b));
303 307
      Node ba=std::max(g.source(b),g.target(b));
304 308
      return ai<bi ||
305 309
        (ai==bi && (aa < ba ||
306 310
                    (aa==ba && ai==g.source(a) && bi==g.target(b))));
307 311
    }
308 312
  };
309 313
  bool isParallel(Arc e,Arc f) const
310 314
  {
311 315
    return (g.source(e)==g.source(f)&&
312 316
            g.target(e)==g.target(f)) ||
313 317
      (g.source(e)==g.target(f)&&
314 318
       g.target(e)==g.source(f));
315 319
  }
316 320
  template<class TT>
317 321
  static std::string psOut(const dim2::Point<TT> &p)
318 322
    {
319 323
      std::ostringstream os;
320 324
      os << p.x << ' ' << p.y;
321 325
      return os.str();
322 326
    }
323 327
  static std::string psOut(const Color &c)
324 328
    {
325 329
      std::ostringstream os;
326 330
      os << c.red() << ' ' << c.green() << ' ' << c.blue();
327 331
      return os.str();
328 332
    }
329 333

	
330 334
public:
331 335
  GraphToEps(const T &t) : T(t), dontPrint(false) {};
332 336

	
333 337
  template<class X> struct CoordsTraits : public T {
334 338
  typedef X CoordsMapType;
335 339
    const X &_coords;
336 340
    CoordsTraits(const T &t,const X &x) : T(t), _coords(x) {}
337 341
  };
338 342
  ///Sets the map of the node coordinates
339 343

	
340 344
  ///Sets the map of the node coordinates.
341 345
  ///\param x must be a node map with \ref dim2::Point "dim2::Point<double>" or
342 346
  ///\ref dim2::Point "dim2::Point<int>" values.
343 347
  template<class X> GraphToEps<CoordsTraits<X> > coords(const X &x) {
344 348
    dontPrint=true;
345 349
    return GraphToEps<CoordsTraits<X> >(CoordsTraits<X>(*this,x));
346 350
  }
347 351
  template<class X> struct NodeSizesTraits : public T {
348 352
    const X &_nodeSizes;
349 353
    NodeSizesTraits(const T &t,const X &x) : T(t), _nodeSizes(x) {}
350 354
  };
351 355
  ///Sets the map of the node sizes
352 356

	
353 357
  ///Sets the map of the node sizes.
354 358
  ///\param x must be a node map with \c double (or convertible) values.
355 359
  template<class X> GraphToEps<NodeSizesTraits<X> > nodeSizes(const X &x)
356 360
  {
357 361
    dontPrint=true;
358 362
    return GraphToEps<NodeSizesTraits<X> >(NodeSizesTraits<X>(*this,x));
359 363
  }
360 364
  template<class X> struct NodeShapesTraits : public T {
361 365
    const X &_nodeShapes;
362 366
    NodeShapesTraits(const T &t,const X &x) : T(t), _nodeShapes(x) {}
363 367
  };
364 368
  ///Sets the map of the node shapes
365 369

	
366 370
  ///Sets the map of the node shapes.
367 371
  ///The available shape values
368 372
  ///can be found in \ref NodeShapes "enum NodeShapes".
369 373
  ///\param x must be a node map with \c int (or convertible) values.
370 374
  ///\sa NodeShapes
371 375
  template<class X> GraphToEps<NodeShapesTraits<X> > nodeShapes(const X &x)
372 376
  {
373 377
    dontPrint=true;
374 378
    return GraphToEps<NodeShapesTraits<X> >(NodeShapesTraits<X>(*this,x));
375 379
  }
376 380
  template<class X> struct NodeTextsTraits : public T {
377 381
    const X &_nodeTexts;
378 382
    NodeTextsTraits(const T &t,const X &x) : T(t), _nodeTexts(x) {}
379 383
  };
380 384
  ///Sets the text printed on the nodes
381 385

	
382 386
  ///Sets the text printed on the nodes.
383 387
  ///\param x must be a node map with type that can be pushed to a standard
384 388
  ///\c ostream.
385 389
  template<class X> GraphToEps<NodeTextsTraits<X> > nodeTexts(const X &x)
386 390
  {
387 391
    dontPrint=true;
388 392
    _showNodeText=true;
389 393
    return GraphToEps<NodeTextsTraits<X> >(NodeTextsTraits<X>(*this,x));
390 394
  }
391 395
  template<class X> struct NodePsTextsTraits : public T {
392 396
    const X &_nodePsTexts;
393 397
    NodePsTextsTraits(const T &t,const X &x) : T(t), _nodePsTexts(x) {}
394 398
  };
395 399
  ///Inserts a PostScript block to the nodes
396 400

	
397 401
  ///With this command it is possible to insert a verbatim PostScript
398 402
  ///block to the nodes.
399 403
  ///The PS current point will be moved to the center of the node before
400 404
  ///the PostScript block inserted.
401 405
  ///
402 406
  ///Before and after the block a newline character is inserted so you
403 407
  ///don't have to bother with the separators.
404 408
  ///
405 409
  ///\param x must be a node map with type that can be pushed to a standard
406 410
  ///\c ostream.
407 411
  ///
408 412
  ///\sa nodePsTextsPreamble()
409 413
  template<class X> GraphToEps<NodePsTextsTraits<X> > nodePsTexts(const X &x)
410 414
  {
411 415
    dontPrint=true;
412 416
    _showNodePsText=true;
413 417
    return GraphToEps<NodePsTextsTraits<X> >(NodePsTextsTraits<X>(*this,x));
414 418
  }
415 419
  template<class X> struct ArcWidthsTraits : public T {
416 420
    const X &_arcWidths;
417 421
    ArcWidthsTraits(const T &t,const X &x) : T(t), _arcWidths(x) {}
418 422
  };
419 423
  ///Sets the map of the arc widths
420 424

	
421 425
  ///Sets the map of the arc widths.
422 426
  ///\param x must be an arc map with \c double (or convertible) values.
423 427
  template<class X> GraphToEps<ArcWidthsTraits<X> > arcWidths(const X &x)
424 428
  {
425 429
    dontPrint=true;
426 430
    return GraphToEps<ArcWidthsTraits<X> >(ArcWidthsTraits<X>(*this,x));
427 431
  }
428 432

	
429 433
  template<class X> struct NodeColorsTraits : public T {
430 434
    const X &_nodeColors;
431 435
    NodeColorsTraits(const T &t,const X &x) : T(t), _nodeColors(x) {}
432 436
  };
433 437
  ///Sets the map of the node colors
434 438

	
435 439
  ///Sets the map of the node colors.
436 440
  ///\param x must be a node map with \ref Color values.
437 441
  ///
438 442
  ///\sa Palette
439 443
  template<class X> GraphToEps<NodeColorsTraits<X> >
440 444
  nodeColors(const X &x)
441 445
  {
442 446
    dontPrint=true;
443 447
    return GraphToEps<NodeColorsTraits<X> >(NodeColorsTraits<X>(*this,x));
444 448
  }
445 449
  template<class X> struct NodeTextColorsTraits : public T {
446 450
    const X &_nodeTextColors;
447 451
    NodeTextColorsTraits(const T &t,const X &x) : T(t), _nodeTextColors(x) {}
448 452
  };
449 453
  ///Sets the map of the node text colors
450 454

	
451 455
  ///Sets the map of the node text colors.
452 456
  ///\param x must be a node map with \ref Color values.
453 457
  ///
454 458
  ///\sa Palette
455 459
  template<class X> GraphToEps<NodeTextColorsTraits<X> >
456 460
  nodeTextColors(const X &x)
457 461
  {
458 462
    dontPrint=true;
459 463
    _nodeTextColorType=CUST_COL;
460 464
    return GraphToEps<NodeTextColorsTraits<X> >
461 465
      (NodeTextColorsTraits<X>(*this,x));
462 466
  }
463 467
  template<class X> struct ArcColorsTraits : public T {
464 468
    const X &_arcColors;
465 469
    ArcColorsTraits(const T &t,const X &x) : T(t), _arcColors(x) {}
466 470
  };
467 471
  ///Sets the map of the arc colors
468 472

	
469 473
  ///Sets the map of the arc colors.
470 474
  ///\param x must be an arc map with \ref Color values.
471 475
  ///
472 476
  ///\sa Palette
473 477
  template<class X> GraphToEps<ArcColorsTraits<X> >
474 478
  arcColors(const X &x)
475 479
  {
476 480
    dontPrint=true;
477 481
    return GraphToEps<ArcColorsTraits<X> >(ArcColorsTraits<X>(*this,x));
478 482
  }
479 483
  ///Sets a global scale factor for node sizes
480 484

	
481 485
  ///Sets a global scale factor for node sizes.
482 486
  ///
483 487
  /// If nodeSizes() is not given, this function simply sets the node
484 488
  /// sizes to \c d.  If nodeSizes() is given, but
485 489
  /// autoNodeScale() is not, then the node size given by
486 490
  /// nodeSizes() will be multiplied by the value \c d.
487 491
  /// If both nodeSizes() and autoNodeScale() are used, then the
488 492
  /// node sizes will be scaled in such a way that the greatest size will be
489 493
  /// equal to \c d.
490 494
  /// \sa nodeSizes()
491 495
  /// \sa autoNodeScale()
492 496
  GraphToEps<T> &nodeScale(double d=.01) {_nodeScale=d;return *this;}
493 497
  ///Turns on/off the automatic node size scaling.
494 498

	
495 499
  ///Turns on/off the automatic node size scaling.
496 500
  ///
497 501
  ///\sa nodeScale()
498 502
  ///
499 503
  GraphToEps<T> &autoNodeScale(bool b=true) {
500 504
    _autoNodeScale=b;return *this;
501 505
  }
502 506

	
503 507
  ///Turns on/off the absolutematic node size scaling.
504 508

	
505 509
  ///Turns on/off the absolutematic node size scaling.
506 510
  ///
507 511
  ///\sa nodeScale()
508 512
  ///
509 513
  GraphToEps<T> &absoluteNodeSizes(bool b=true) {
510 514
    _absoluteNodeSizes=b;return *this;
511 515
  }
512 516

	
513 517
  ///Negates the Y coordinates.
514 518
  GraphToEps<T> &negateY(bool b=true) {
515 519
    _negY=b;return *this;
516 520
  }
517 521

	
518 522
  ///Turn on/off pre-scaling
519 523

	
520 524
  ///By default graphToEps() rescales the whole image in order to avoid
521 525
  ///very big or very small bounding boxes.
522 526
  ///
523 527
  ///This (p)rescaling can be turned off with this function.
524 528
  ///
525 529
  GraphToEps<T> &preScale(bool b=true) {
526 530
    _preScale=b;return *this;
527 531
  }
528 532

	
529 533
  ///Sets a global scale factor for arc widths
530 534

	
531 535
  /// Sets a global scale factor for arc widths.
532 536
  ///
533 537
  /// If arcWidths() is not given, this function simply sets the arc
534 538
  /// widths to \c d.  If arcWidths() is given, but
535 539
  /// autoArcWidthScale() is not, then the arc withs given by
536 540
  /// arcWidths() will be multiplied by the value \c d.
537 541
  /// If both arcWidths() and autoArcWidthScale() are used, then the
538 542
  /// arc withs will be scaled in such a way that the greatest width will be
539 543
  /// equal to \c d.
540 544
  GraphToEps<T> &arcWidthScale(double d=.003) {_arcWidthScale=d;return *this;}
541 545
  ///Turns on/off the automatic arc width scaling.
542 546

	
543 547
  ///Turns on/off the automatic arc width scaling.
544 548
  ///
545 549
  ///\sa arcWidthScale()
546 550
  ///
547 551
  GraphToEps<T> &autoArcWidthScale(bool b=true) {
548 552
    _autoArcWidthScale=b;return *this;
549 553
  }
550 554
  ///Turns on/off the absolutematic arc width scaling.
551 555

	
552 556
  ///Turns on/off the absolutematic arc width scaling.
553 557
  ///
554 558
  ///\sa arcWidthScale()
555 559
  ///
556 560
  GraphToEps<T> &absoluteArcWidths(bool b=true) {
557 561
    _absoluteArcWidths=b;return *this;
558 562
  }
559 563
  ///Sets a global scale factor for the whole picture
560 564
  GraphToEps<T> &scale(double d) {_scale=d;return *this;}
561 565
  ///Sets the width of the border around the picture
562 566
  GraphToEps<T> &border(double b=10) {_xBorder=_yBorder=b;return *this;}
563 567
  ///Sets the width of the border around the picture
564 568
  GraphToEps<T> &border(double x, double y) {
565 569
    _xBorder=x;_yBorder=y;return *this;
566 570
  }
567 571
  ///Sets whether to draw arrows
568 572
  GraphToEps<T> &drawArrows(bool b=true) {_drawArrows=b;return *this;}
569 573
  ///Sets the length of the arrowheads
570 574
  GraphToEps<T> &arrowLength(double d=1.0) {_arrowLength*=d;return *this;}
571 575
  ///Sets the width of the arrowheads
572 576
  GraphToEps<T> &arrowWidth(double d=.3) {_arrowWidth*=d;return *this;}
573 577

	
574 578
  ///Scales the drawing to fit to A4 page
575 579
  GraphToEps<T> &scaleToA4() {_scaleToA4=true;return *this;}
576 580

	
577 581
  ///Enables parallel arcs
578 582
  GraphToEps<T> &enableParallel(bool b=true) {_enableParallel=b;return *this;}
579 583

	
580 584
  ///Sets the distance between parallel arcs
581 585
  GraphToEps<T> &parArcDist(double d) {_parArcDist*=d;return *this;}
582 586

	
583 587
  ///Hides the arcs
584 588
  GraphToEps<T> &hideArcs(bool b=true) {_showArcs=!b;return *this;}
585 589
  ///Hides the nodes
586 590
  GraphToEps<T> &hideNodes(bool b=true) {_showNodes=!b;return *this;}
587 591

	
588 592
  ///Sets the size of the node texts
589 593
  GraphToEps<T> &nodeTextSize(double d) {_nodeTextSize=d;return *this;}
590 594

	
591 595
  ///Sets the color of the node texts to be different from the node color
592 596

	
593 597
  ///Sets the color of the node texts to be as different from the node color
594 598
  ///as it is possible.
595 599
  GraphToEps<T> &distantColorNodeTexts()
596 600
  {_nodeTextColorType=DIST_COL;return *this;}
597 601
  ///Sets the color of the node texts to be black or white and always visible.
598 602

	
599 603
  ///Sets the color of the node texts to be black or white according to
600 604
  ///which is more different from the node color.
601 605
  GraphToEps<T> &distantBWNodeTexts()
602 606
  {_nodeTextColorType=DIST_BW;return *this;}
603 607

	
604 608
  ///Gives a preamble block for node Postscript block.
605 609

	
606 610
  ///Gives a preamble block for node Postscript block.
607 611
  ///
608 612
  ///\sa nodePsTexts()
609 613
  GraphToEps<T> & nodePsTextsPreamble(const char *str) {
610 614
    _nodePsTextsPreamble=str ;return *this;
611 615
  }
612 616
  ///Sets whether the graph is undirected
613 617

	
614 618
  ///Sets whether the graph is undirected.
615 619
  ///
616 620
  ///This setting is the default for undirected graphs.
617 621
  ///
618 622
  ///\sa directed()
619 623
   GraphToEps<T> &undirected(bool b=true) {_undirected=b;return *this;}
620 624

	
621 625
  ///Sets whether the graph is directed
622 626

	
623 627
  ///Sets whether the graph is directed.
624 628
  ///Use it to show the edges as a pair of directed ones.
625 629
  ///
626 630
  ///This setting is the default for digraphs.
627 631
  ///
628 632
  ///\sa undirected()
629 633
  GraphToEps<T> &directed(bool b=true) {_undirected=!b;return *this;}
630 634

	
631 635
  ///Sets the title.
632 636

	
633 637
  ///Sets the title of the generated image,
634 638
  ///namely it inserts a <tt>%%Title:</tt> DSC field to the header of
635 639
  ///the EPS file.
636 640
  GraphToEps<T> &title(const std::string &t) {_title=t;return *this;}
637 641
  ///Sets the copyright statement.
638 642

	
639 643
  ///Sets the copyright statement of the generated image,
640 644
  ///namely it inserts a <tt>%%Copyright:</tt> DSC field to the header of
641 645
  ///the EPS file.
642 646
  GraphToEps<T> &copyright(const std::string &t) {_copyright=t;return *this;}
643 647

	
644 648
protected:
645 649
  bool isInsideNode(dim2::Point<double> p, double r,int t)
646 650
  {
647 651
    switch(t) {
648 652
    case CIRCLE:
649 653
    case MALE:
650 654
    case FEMALE:
651 655
      return p.normSquare()<=r*r;
652 656
    case SQUARE:
653 657
      return p.x<=r&&p.x>=-r&&p.y<=r&&p.y>=-r;
654 658
    case DIAMOND:
655 659
      return p.x+p.y<=r && p.x-p.y<=r && -p.x+p.y<=r && -p.x-p.y<=r;
656 660
    }
657 661
    return false;
658 662
  }
659 663

	
660 664
public:
661 665
  ~GraphToEps() { }
662 666

	
663 667
  ///Draws the graph.
664 668

	
665 669
  ///Like other functions using
666 670
  ///\ref named-templ-func-param "named template parameters",
667 671
  ///this function calls the algorithm itself, i.e. in this case
668 672
  ///it draws the graph.
669 673
  void run() {
670 674
    const double EPSILON=1e-9;
671 675
    if(dontPrint) return;
672 676

	
673 677
    _graph_to_eps_bits::_NegY<typename T::CoordsMapType>
674 678
      mycoords(_coords,_negY);
675 679

	
676 680
    os << "%!PS-Adobe-2.0 EPSF-2.0\n";
677 681
    if(_title.size()>0) os << "%%Title: " << _title << '\n';
678 682
     if(_copyright.size()>0) os << "%%Copyright: " << _copyright << '\n';
679 683
    os << "%%Creator: LEMON, graphToEps()\n";
680 684

	
681 685
    {
682 686
#ifndef WIN32
683 687
      timeval tv;
684 688
      gettimeofday(&tv, 0);
685 689

	
686 690
      char cbuf[26];
687 691
      ctime_r(&tv.tv_sec,cbuf);
688 692
      os << "%%CreationDate: " << cbuf;
689 693
#else
690 694
      SYSTEMTIME time;
691
      char buf1[11], buf2[9], buf3[5];
692

	
693 695
      GetSystemTime(&time);
696
#if defined(_MSC_VER) && (_MSC_VER < 1500)
697
      LPWSTR buf1, buf2, buf3;
694 698
      if (GetDateFormat(LOCALE_USER_DEFAULT, 0, &time,
695
                        "ddd MMM dd", buf1, 11) &&
699
                        L"ddd MMM dd", buf1, 11) &&
696 700
          GetTimeFormat(LOCALE_USER_DEFAULT, 0, &time,
697
                        "HH':'mm':'ss", buf2, 9) &&
701
                        L"HH':'mm':'ss", buf2, 9) &&
698 702
          GetDateFormat(LOCALE_USER_DEFAULT, 0, &time,
699
                                "yyyy", buf3, 5)) {
703
                        L"yyyy", buf3, 5)) {
700 704
        os << "%%CreationDate: " << buf1 << ' '
701 705
           << buf2 << ' ' << buf3 << std::endl;
702 706
      }
707
#else
708
        char buf1[11], buf2[9], buf3[5];
709
        if (GetDateFormat(LOCALE_USER_DEFAULT, 0, &time,
710
                          "ddd MMM dd", buf1, 11) &&
711
            GetTimeFormat(LOCALE_USER_DEFAULT, 0, &time,
712
                          "HH':'mm':'ss", buf2, 9) &&
713
            GetDateFormat(LOCALE_USER_DEFAULT, 0, &time,
714
                          "yyyy", buf3, 5)) {
715
          os << "%%CreationDate: " << buf1 << ' '
716
             << buf2 << ' ' << buf3 << std::endl;
717
        }
718
#endif
703 719
#endif
704 720
    }
705 721

	
706 722
    if (_autoArcWidthScale) {
707 723
      double max_w=0;
708 724
      for(ArcIt e(g);e!=INVALID;++e)
709 725
        max_w=std::max(double(_arcWidths[e]),max_w);
710 726
      if(max_w>EPSILON) {
711 727
        _arcWidthScale/=max_w;
712 728
      }
713 729
    }
714 730

	
715 731
    if (_autoNodeScale) {
716 732
      double max_s=0;
717 733
      for(NodeIt n(g);n!=INVALID;++n)
718 734
        max_s=std::max(double(_nodeSizes[n]),max_s);
719 735
      if(max_s>EPSILON) {
720 736
        _nodeScale/=max_s;
721 737
      }
722 738
    }
723 739

	
724 740
    double diag_len = 1;
725 741
    if(!(_absoluteNodeSizes&&_absoluteArcWidths)) {
726 742
      dim2::Box<double> bb;
727 743
      for(NodeIt n(g);n!=INVALID;++n) bb.add(mycoords[n]);
728 744
      if (bb.empty()) {
729 745
        bb = dim2::Box<double>(dim2::Point<double>(0,0));
730 746
      }
731 747
      diag_len = std::sqrt((bb.bottomLeft()-bb.topRight()).normSquare());
732 748
      if(diag_len<EPSILON) diag_len = 1;
733 749
      if(!_absoluteNodeSizes) _nodeScale*=diag_len;
734 750
      if(!_absoluteArcWidths) _arcWidthScale*=diag_len;
735 751
    }
736 752

	
737 753
    dim2::Box<double> bb;
738 754
    for(NodeIt n(g);n!=INVALID;++n) {
739 755
      double ns=_nodeSizes[n]*_nodeScale;
740 756
      dim2::Point<double> p(ns,ns);
741 757
      switch(_nodeShapes[n]) {
742 758
      case CIRCLE:
743 759
      case SQUARE:
744 760
      case DIAMOND:
745 761
        bb.add(p+mycoords[n]);
746 762
        bb.add(-p+mycoords[n]);
747 763
        break;
748 764
      case MALE:
749 765
        bb.add(-p+mycoords[n]);
750 766
        bb.add(dim2::Point<double>(1.5*ns,1.5*std::sqrt(3.0)*ns)+mycoords[n]);
751 767
        break;
752 768
      case FEMALE:
753 769
        bb.add(p+mycoords[n]);
754 770
        bb.add(dim2::Point<double>(-ns,-3.01*ns)+mycoords[n]);
755 771
        break;
756 772
      }
757 773
    }
758 774
    if (bb.empty()) {
759 775
      bb = dim2::Box<double>(dim2::Point<double>(0,0));
760 776
    }
761 777

	
762 778
    if(_scaleToA4)
763 779
      os <<"%%BoundingBox: 0 0 596 842\n%%DocumentPaperSizes: a4\n";
764 780
    else {
765 781
      if(_preScale) {
766 782
        //Rescale so that BoundingBox won't be neither to big nor too small.
767 783
        while(bb.height()*_scale>1000||bb.width()*_scale>1000) _scale/=10;
768 784
        while(bb.height()*_scale<100||bb.width()*_scale<100) _scale*=10;
769 785
      }
770 786

	
771 787
      os << "%%BoundingBox: "
772 788
         << int(floor(bb.left()   * _scale - _xBorder)) << ' '
773 789
         << int(floor(bb.bottom() * _scale - _yBorder)) << ' '
774 790
         << int(ceil(bb.right()  * _scale + _xBorder)) << ' '
775 791
         << int(ceil(bb.top()    * _scale + _yBorder)) << '\n';
776 792
    }
777 793

	
778 794
    os << "%%EndComments\n";
779 795

	
780 796
    //x1 y1 x2 y2 x3 y3 cr cg cb w
781 797
    os << "/lb { setlinewidth setrgbcolor newpath moveto\n"
782 798
       << "      4 2 roll 1 index 1 index curveto stroke } bind def\n";
783 799
    os << "/l { setlinewidth setrgbcolor newpath moveto lineto stroke }"
784 800
       << " bind def\n";
785 801
    //x y r
786 802
    os << "/c { newpath dup 3 index add 2 index moveto 0 360 arc closepath }"
787 803
       << " bind def\n";
788 804
    //x y r
789 805
    os << "/sq { newpath 2 index 1 index add 2 index 2 index add moveto\n"
790 806
       << "      2 index 1 index sub 2 index 2 index add lineto\n"
791 807
       << "      2 index 1 index sub 2 index 2 index sub lineto\n"
792 808
       << "      2 index 1 index add 2 index 2 index sub lineto\n"
793 809
       << "      closepath pop pop pop} bind def\n";
794 810
    //x y r
795 811
    os << "/di { newpath 2 index 1 index add 2 index moveto\n"
796 812
       << "      2 index             2 index 2 index add lineto\n"
797 813
       << "      2 index 1 index sub 2 index             lineto\n"
798 814
       << "      2 index             2 index 2 index sub lineto\n"
799 815
       << "      closepath pop pop pop} bind def\n";
800 816
    // x y r cr cg cb
801 817
    os << "/nc { 0 0 0 setrgbcolor 5 index 5 index 5 index c fill\n"
802 818
       << "     setrgbcolor " << 1+_nodeBorderQuotient << " div c fill\n"
803 819
       << "   } bind def\n";
804 820
    os << "/nsq { 0 0 0 setrgbcolor 5 index 5 index 5 index sq fill\n"
805 821
       << "     setrgbcolor " << 1+_nodeBorderQuotient << " div sq fill\n"
806 822
       << "   } bind def\n";
807 823
    os << "/ndi { 0 0 0 setrgbcolor 5 index 5 index 5 index di fill\n"
808 824
       << "     setrgbcolor " << 1+_nodeBorderQuotient << " div di fill\n"
809 825
       << "   } bind def\n";
810 826
    os << "/nfemale { 0 0 0 setrgbcolor 3 index "
811 827
       << _nodeBorderQuotient/(1+_nodeBorderQuotient)
812 828
       << " 1.5 mul mul setlinewidth\n"
813 829
       << "  newpath 5 index 5 index moveto "
814 830
       << "5 index 5 index 5 index 3.01 mul sub\n"
815 831
       << "  lineto 5 index 4 index .7 mul sub 5 index 5 index 2.2 mul sub"
816 832
       << " moveto\n"
817 833
       << "  5 index 4 index .7 mul add 5 index 5 index 2.2 mul sub lineto "
818 834
       << "stroke\n"
819 835
       << "  5 index 5 index 5 index c fill\n"
820 836
       << "  setrgbcolor " << 1+_nodeBorderQuotient << " div c fill\n"
821 837
       << "  } bind def\n";
822 838
    os << "/nmale {\n"
823 839
       << "  0 0 0 setrgbcolor 3 index "
824 840
       << _nodeBorderQuotient/(1+_nodeBorderQuotient)
825 841
       <<" 1.5 mul mul setlinewidth\n"
826 842
       << "  newpath 5 index 5 index moveto\n"
827 843
       << "  5 index 4 index 1 mul 1.5 mul add\n"
828 844
       << "  5 index 5 index 3 sqrt 1.5 mul mul add\n"
829 845
       << "  1 index 1 index lineto\n"
830 846
       << "  1 index 1 index 7 index sub moveto\n"
831 847
       << "  1 index 1 index lineto\n"
832 848
       << "  exch 5 index 3 sqrt .5 mul mul sub exch 5 index .5 mul sub"
833 849
       << " lineto\n"
834 850
       << "  stroke\n"
835 851
       << "  5 index 5 index 5 index c fill\n"
836 852
       << "  setrgbcolor " << 1+_nodeBorderQuotient << " div c fill\n"
837 853
       << "  } bind def\n";
838 854

	
839 855

	
840 856
    os << "/arrl " << _arrowLength << " def\n";
841 857
    os << "/arrw " << _arrowWidth << " def\n";
842 858
    // l dx_norm dy_norm
843 859
    os << "/lrl { 2 index mul exch 2 index mul exch rlineto pop} bind def\n";
844 860
    //len w dx_norm dy_norm x1 y1 cr cg cb
845 861
    os << "/arr { setrgbcolor /y1 exch def /x1 exch def /dy exch def /dx "
846 862
       << "exch def\n"
847 863
       << "       /w exch def /len exch def\n"
848 864
      //<< "0.1 setlinewidth x1 y1 moveto dx len mul dy len mul rlineto stroke"
849 865
       << "       newpath x1 dy w 2 div mul add y1 dx w 2 div mul sub moveto\n"
850 866
       << "       len w sub arrl sub dx dy lrl\n"
851 867
       << "       arrw dy dx neg lrl\n"
852 868
       << "       dx arrl w add mul dy w 2 div arrw add mul sub\n"
853 869
       << "       dy arrl w add mul dx w 2 div arrw add mul add rlineto\n"
854 870
       << "       dx arrl w add mul neg dy w 2 div arrw add mul sub\n"
855 871
       << "       dy arrl w add mul neg dx w 2 div arrw add mul add rlineto\n"
856 872
       << "       arrw dy dx neg lrl\n"
857 873
       << "       len w sub arrl sub neg dx dy lrl\n"
858 874
       << "       closepath fill } bind def\n";
859 875
    os << "/cshow { 2 index 2 index moveto dup stringwidth pop\n"
860 876
       << "         neg 2 div fosi .35 mul neg rmoveto show pop pop} def\n";
861 877

	
862 878
    os << "\ngsave\n";
863 879
    if(_scaleToA4)
864 880
      if(bb.height()>bb.width()) {
865 881
        double sc= std::min((A4HEIGHT-2*A4BORDER)/bb.height(),
866 882
                  (A4WIDTH-2*A4BORDER)/bb.width());
867 883
        os << ((A4WIDTH -2*A4BORDER)-sc*bb.width())/2 + A4BORDER << ' '
868 884
           << ((A4HEIGHT-2*A4BORDER)-sc*bb.height())/2 + A4BORDER
869 885
           << " translate\n"
870 886
           << sc << " dup scale\n"
871 887
           << -bb.left() << ' ' << -bb.bottom() << " translate\n";
872 888
      }
873 889
      else {
874 890
        double sc= std::min((A4HEIGHT-2*A4BORDER)/bb.width(),
875 891
                  (A4WIDTH-2*A4BORDER)/bb.height());
876 892
        os << ((A4WIDTH -2*A4BORDER)-sc*bb.height())/2 + A4BORDER << ' '
877 893
           << ((A4HEIGHT-2*A4BORDER)-sc*bb.width())/2 + A4BORDER
878 894
           << " translate\n"
879 895
           << sc << " dup scale\n90 rotate\n"
880 896
           << -bb.left() << ' ' << -bb.top() << " translate\n";
881 897
        }
882 898
    else if(_scale!=1.0) os << _scale << " dup scale\n";
883 899

	
884 900
    if(_showArcs) {
885 901
      os << "%Arcs:\ngsave\n";
886 902
      if(_enableParallel) {
887 903
        std::vector<Arc> el;
888 904
        for(ArcIt e(g);e!=INVALID;++e)
889 905
          if((!_undirected||g.source(e)<g.target(e))&&_arcWidths[e]>0
890 906
             &&g.source(e)!=g.target(e))
891 907
            el.push_back(e);
892 908
        std::sort(el.begin(),el.end(),arcLess(g));
893 909

	
894 910
        typename std::vector<Arc>::iterator j;
895 911
        for(typename std::vector<Arc>::iterator i=el.begin();i!=el.end();i=j) {
896 912
          for(j=i+1;j!=el.end()&&isParallel(*i,*j);++j) ;
897 913

	
898 914
          double sw=0;
899 915
          for(typename std::vector<Arc>::iterator e=i;e!=j;++e)
900 916
            sw+=_arcWidths[*e]*_arcWidthScale+_parArcDist;
901 917
          sw-=_parArcDist;
902 918
          sw/=-2.0;
903 919
          dim2::Point<double>
904 920
            dvec(mycoords[g.target(*i)]-mycoords[g.source(*i)]);
905 921
          double l=std::sqrt(dvec.normSquare());
906 922
          dim2::Point<double> d(dvec/std::max(l,EPSILON));
907 923
          dim2::Point<double> m;
908 924
//           m=dim2::Point<double>(mycoords[g.target(*i)]+
909 925
//                                 mycoords[g.source(*i)])/2.0;
910 926

	
911 927
//            m=dim2::Point<double>(mycoords[g.source(*i)])+
912 928
//             dvec*(double(_nodeSizes[g.source(*i)])/
913 929
//                (_nodeSizes[g.source(*i)]+_nodeSizes[g.target(*i)]));
914 930

	
915 931
          m=dim2::Point<double>(mycoords[g.source(*i)])+
916 932
            d*(l+_nodeSizes[g.source(*i)]-_nodeSizes[g.target(*i)])/2.0;
917 933

	
918 934
          for(typename std::vector<Arc>::iterator e=i;e!=j;++e) {
919 935
            sw+=_arcWidths[*e]*_arcWidthScale/2.0;
920 936
            dim2::Point<double> mm=m+rot90(d)*sw/.75;
921 937
            if(_drawArrows) {
922 938
              int node_shape;
923 939
              dim2::Point<double> s=mycoords[g.source(*e)];
924 940
              dim2::Point<double> t=mycoords[g.target(*e)];
925 941
              double rn=_nodeSizes[g.target(*e)]*_nodeScale;
926 942
              node_shape=_nodeShapes[g.target(*e)];
927 943
              dim2::Bezier3 bez(s,mm,mm,t);
928 944
              double t1=0,t2=1;
929 945
              for(int ii=0;ii<INTERPOL_PREC;++ii)
930 946
                if(isInsideNode(bez((t1+t2)/2)-t,rn,node_shape)) t2=(t1+t2)/2;
931 947
                else t1=(t1+t2)/2;
932 948
              dim2::Point<double> apoint=bez((t1+t2)/2);
933 949
              rn = _arrowLength+_arcWidths[*e]*_arcWidthScale;
934 950
              rn*=rn;
935 951
              t2=(t1+t2)/2;t1=0;
936 952
              for(int ii=0;ii<INTERPOL_PREC;++ii)
937 953
                if((bez((t1+t2)/2)-apoint).normSquare()>rn) t1=(t1+t2)/2;
938 954
                else t2=(t1+t2)/2;
939 955
              dim2::Point<double> linend=bez((t1+t2)/2);
940 956
              bez=bez.before((t1+t2)/2);
941 957
//               rn=_nodeSizes[g.source(*e)]*_nodeScale;
942 958
//               node_shape=_nodeShapes[g.source(*e)];
943 959
//               t1=0;t2=1;
944 960
//               for(int i=0;i<INTERPOL_PREC;++i)
945 961
//                 if(isInsideNode(bez((t1+t2)/2)-t,rn,node_shape))
946 962
//                   t1=(t1+t2)/2;
947 963
//                 else t2=(t1+t2)/2;
948 964
//               bez=bez.after((t1+t2)/2);
949 965
              os << _arcWidths[*e]*_arcWidthScale << " setlinewidth "
950 966
                 << _arcColors[*e].red() << ' '
951 967
                 << _arcColors[*e].green() << ' '
952 968
                 << _arcColors[*e].blue() << " setrgbcolor newpath\n"
953 969
                 << bez.p1.x << ' ' <<  bez.p1.y << " moveto\n"
954 970
                 << bez.p2.x << ' ' << bez.p2.y << ' '
955 971
                 << bez.p3.x << ' ' << bez.p3.y << ' '
956 972
                 << bez.p4.x << ' ' << bez.p4.y << " curveto stroke\n";
957 973
              dim2::Point<double> dd(rot90(linend-apoint));
958 974
              dd*=(.5*_arcWidths[*e]*_arcWidthScale+_arrowWidth)/
959 975
                std::sqrt(dd.normSquare());
960 976
              os << "newpath " << psOut(apoint) << " moveto "
961 977
                 << psOut(linend+dd) << " lineto "
962 978
                 << psOut(linend-dd) << " lineto closepath fill\n";
963 979
            }
964 980
            else {
965 981
              os << mycoords[g.source(*e)].x << ' '
966 982
                 << mycoords[g.source(*e)].y << ' '
967 983
                 << mm.x << ' ' << mm.y << ' '
968 984
                 << mycoords[g.target(*e)].x << ' '
969 985
                 << mycoords[g.target(*e)].y << ' '
970 986
                 << _arcColors[*e].red() << ' '
971 987
                 << _arcColors[*e].green() << ' '
972 988
                 << _arcColors[*e].blue() << ' '
973 989
                 << _arcWidths[*e]*_arcWidthScale << " lb\n";
974 990
            }
975 991
            sw+=_arcWidths[*e]*_arcWidthScale/2.0+_parArcDist;
976 992
          }
977 993
        }
978 994
      }
979 995
      else for(ArcIt e(g);e!=INVALID;++e)
980 996
        if((!_undirected||g.source(e)<g.target(e))&&_arcWidths[e]>0
981 997
           &&g.source(e)!=g.target(e)) {
982 998
          if(_drawArrows) {
983 999
            dim2::Point<double> d(mycoords[g.target(e)]-mycoords[g.source(e)]);
984 1000
            double rn=_nodeSizes[g.target(e)]*_nodeScale;
985 1001
            int node_shape=_nodeShapes[g.target(e)];
986 1002
            double t1=0,t2=1;
987 1003
            for(int i=0;i<INTERPOL_PREC;++i)
988 1004
              if(isInsideNode((-(t1+t2)/2)*d,rn,node_shape)) t1=(t1+t2)/2;
989 1005
              else t2=(t1+t2)/2;
990 1006
            double l=std::sqrt(d.normSquare());
991 1007
            d/=l;
992 1008

	
993 1009
            os << l*(1-(t1+t2)/2) << ' '
994 1010
               << _arcWidths[e]*_arcWidthScale << ' '
995 1011
               << d.x << ' ' << d.y << ' '
996 1012
               << mycoords[g.source(e)].x << ' '
997 1013
               << mycoords[g.source(e)].y << ' '
998 1014
               << _arcColors[e].red() << ' '
999 1015
               << _arcColors[e].green() << ' '
1000 1016
               << _arcColors[e].blue() << " arr\n";
1001 1017
          }
1002 1018
          else os << mycoords[g.source(e)].x << ' '
1003 1019
                  << mycoords[g.source(e)].y << ' '
1004 1020
                  << mycoords[g.target(e)].x << ' '
1005 1021
                  << mycoords[g.target(e)].y << ' '
1006 1022
                  << _arcColors[e].red() << ' '
1007 1023
                  << _arcColors[e].green() << ' '
1008 1024
                  << _arcColors[e].blue() << ' '
1009 1025
                  << _arcWidths[e]*_arcWidthScale << " l\n";
1010 1026
        }
1011 1027
      os << "grestore\n";
1012 1028
    }
1013 1029
    if(_showNodes) {
1014 1030
      os << "%Nodes:\ngsave\n";
1015 1031
      for(NodeIt n(g);n!=INVALID;++n) {
1016 1032
        os << mycoords[n].x << ' ' << mycoords[n].y << ' '
1017 1033
           << _nodeSizes[n]*_nodeScale << ' '
1018 1034
           << _nodeColors[n].red() << ' '
1019 1035
           << _nodeColors[n].green() << ' '
1020 1036
           << _nodeColors[n].blue() << ' ';
1021 1037
        switch(_nodeShapes[n]) {
1022 1038
        case CIRCLE:
1023 1039
          os<< "nc";break;
1024 1040
        case SQUARE:
1025 1041
          os<< "nsq";break;
1026 1042
        case DIAMOND:
1027 1043
          os<< "ndi";break;
1028 1044
        case MALE:
1029 1045
          os<< "nmale";break;
1030 1046
        case FEMALE:
1031 1047
          os<< "nfemale";break;
1032 1048
        }
1033 1049
        os<<'\n';
1034 1050
      }
1035 1051
      os << "grestore\n";
1036 1052
    }
1037 1053
    if(_showNodeText) {
1038 1054
      os << "%Node texts:\ngsave\n";
1039 1055
      os << "/fosi " << _nodeTextSize << " def\n";
1040 1056
      os << "(Helvetica) findfont fosi scalefont setfont\n";
1041 1057
      for(NodeIt n(g);n!=INVALID;++n) {
1042 1058
        switch(_nodeTextColorType) {
1043 1059
        case DIST_COL:
1044 1060
          os << psOut(distantColor(_nodeColors[n])) << " setrgbcolor\n";
1045 1061
          break;
1046 1062
        case DIST_BW:
1047 1063
          os << psOut(distantBW(_nodeColors[n])) << " setrgbcolor\n";
1048 1064
          break;
1049 1065
        case CUST_COL:
1050 1066
          os << psOut(distantColor(_nodeTextColors[n])) << " setrgbcolor\n";
1051 1067
          break;
1052 1068
        default:
1053 1069
          os << "0 0 0 setrgbcolor\n";
1054 1070
        }
1055 1071
        os << mycoords[n].x << ' ' << mycoords[n].y
1056 1072
           << " (" << _nodeTexts[n] << ") cshow\n";
1057 1073
      }
1058 1074
      os << "grestore\n";
1059 1075
    }
1060 1076
    if(_showNodePsText) {
1061 1077
      os << "%Node PS blocks:\ngsave\n";
1062 1078
      for(NodeIt n(g);n!=INVALID;++n)
1063 1079
        os << mycoords[n].x << ' ' << mycoords[n].y
1064 1080
           << " moveto\n" << _nodePsTexts[n] << "\n";
1065 1081
      os << "grestore\n";
1066 1082
    }
1067 1083

	
1068 1084
    os << "grestore\nshowpage\n";
1069 1085

	
1070 1086
    //CleanUp:
1071 1087
    if(_pleaseRemoveOsStream) {delete &os;}
1072 1088
  }
1073 1089

	
1074 1090
  ///\name Aliases
1075 1091
  ///These are just some aliases to other parameter setting functions.
1076 1092

	
1077 1093
  ///@{
1078 1094

	
1079 1095
  ///An alias for arcWidths()
1080 1096
  template<class X> GraphToEps<ArcWidthsTraits<X> > edgeWidths(const X &x)
1081 1097
  {
1082 1098
    return arcWidths(x);
1083 1099
  }
1084 1100

	
1085 1101
  ///An alias for arcColors()
1086 1102
  template<class X> GraphToEps<ArcColorsTraits<X> >
1087 1103
  edgeColors(const X &x)
1088 1104
  {
1089 1105
    return arcColors(x);
1090 1106
  }
1091 1107

	
1092 1108
  ///An alias for arcWidthScale()
1093 1109
  GraphToEps<T> &edgeWidthScale(double d) {return arcWidthScale(d);}
1094 1110

	
1095 1111
  ///An alias for autoArcWidthScale()
1096 1112
  GraphToEps<T> &autoEdgeWidthScale(bool b=true)
1097 1113
  {
1098 1114
    return autoArcWidthScale(b);
1099 1115
  }
1100 1116

	
1101 1117
  ///An alias for absoluteArcWidths()
1102 1118
  GraphToEps<T> &absoluteEdgeWidths(bool b=true)
1103 1119
  {
1104 1120
    return absoluteArcWidths(b);
1105 1121
  }
1106 1122

	
1107 1123
  ///An alias for parArcDist()
1108 1124
  GraphToEps<T> &parEdgeDist(double d) {return parArcDist(d);}
1109 1125

	
1110 1126
  ///An alias for hideArcs()
1111 1127
  GraphToEps<T> &hideEdges(bool b=true) {return hideArcs(b);}
1112 1128

	
1113 1129
  ///@}
1114 1130
};
1115 1131

	
1116 1132
template<class T>
1117 1133
const int GraphToEps<T>::INTERPOL_PREC = 20;
1118 1134
template<class T>
1119 1135
const double GraphToEps<T>::A4HEIGHT = 841.8897637795276;
1120 1136
template<class T>
1121 1137
const double GraphToEps<T>::A4WIDTH  = 595.275590551181;
1122 1138
template<class T>
1123 1139
const double GraphToEps<T>::A4BORDER = 15;
1124 1140

	
1125 1141

	
1126 1142
///Generates an EPS file from a graph
1127 1143

	
1128 1144
///\ingroup eps_io
1129 1145
///Generates an EPS file from a graph.
1130 1146
///\param g Reference to the graph to be printed.
1131 1147
///\param os Reference to the output stream.
1132 1148
///By default it is <tt>std::cout</tt>.
1133 1149
///
1134 1150
///This function also has a lot of
1135 1151
///\ref named-templ-func-param "named parameters",
1136 1152
///they are declared as the members of class \ref GraphToEps. The following
1137 1153
///example shows how to use these parameters.
1138 1154
///\code
1139 1155
/// graphToEps(g,os).scale(10).coords(coords)
1140 1156
///              .nodeScale(2).nodeSizes(sizes)
1141 1157
///              .arcWidthScale(.4).run();
1142 1158
///\endcode
1143 1159
///
1144 1160
///For more detailed examples see the \ref graph_to_eps_demo.cc demo file.
1145 1161
///
1146 1162
///\warning Don't forget to put the \ref GraphToEps::run() "run()"
1147 1163
///to the end of the parameter list.
1148 1164
///\sa GraphToEps
1149 1165
///\sa graphToEps(G &g, const char *file_name)
1150 1166
template<class G>
1151 1167
GraphToEps<DefaultGraphToEpsTraits<G> >
1152 1168
graphToEps(G &g, std::ostream& os=std::cout)
1153 1169
{
1154 1170
  return
1155 1171
    GraphToEps<DefaultGraphToEpsTraits<G> >(DefaultGraphToEpsTraits<G>(g,os));
1156 1172
}
1157 1173

	
1158 1174
///Generates an EPS file from a graph
1159 1175

	
1160 1176
///\ingroup eps_io
1161 1177
///This function does the same as
1162 1178
///\ref graphToEps(G &g,std::ostream& os)
1163 1179
///but it writes its output into the file \c file_name
1164 1180
///instead of a stream.
1165 1181
///\sa graphToEps(G &g, std::ostream& os)
1166 1182
template<class G>
1167 1183
GraphToEps<DefaultGraphToEpsTraits<G> >
1168 1184
graphToEps(G &g,const char *file_name)
1169 1185
{
1170 1186
  std::ostream* os = new std::ofstream(file_name);
1171 1187
  if (!(*os)) {
1172 1188
    delete os;
1173 1189
    throw IoError("Cannot write file", file_name);
1174 1190
  }
1175 1191
  return GraphToEps<DefaultGraphToEpsTraits<G> >
1176 1192
    (DefaultGraphToEpsTraits<G>(g,*os,true));
1177 1193
}
1178 1194

	
1179 1195
///Generates an EPS file from a graph
1180 1196

	
1181 1197
///\ingroup eps_io
1182 1198
///This function does the same as
1183 1199
///\ref graphToEps(G &g,std::ostream& os)
1184 1200
///but it writes its output into the file \c file_name
1185 1201
///instead of a stream.
1186 1202
///\sa graphToEps(G &g, std::ostream& os)
1187 1203
template<class G>
1188 1204
GraphToEps<DefaultGraphToEpsTraits<G> >
1189 1205
graphToEps(G &g,const std::string& file_name)
1190 1206
{
1191 1207
  std::ostream* os = new std::ofstream(file_name.c_str());
1192 1208
  if (!(*os)) {
1193 1209
    delete os;
1194 1210
    throw IoError("Cannot write file", file_name);
1195 1211
  }
1196 1212
  return GraphToEps<DefaultGraphToEpsTraits<G> >
1197 1213
    (DefaultGraphToEpsTraits<G>(g,*os,true));
1198 1214
}
1199 1215

	
1200 1216
} //END OF NAMESPACE LEMON
1201 1217

	
1202 1218
#endif // LEMON_GRAPH_TO_EPS_H
Ignore white space 16777216 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2008
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
///\file
20 20
///\brief The implementation of the LP solver interface.
21 21

	
22 22
#include <lemon/lp_base.h>
23 23
namespace lemon {
24 24

	
25
  const LpBase::Value LpBase::INF = std::numeric_limits<Value>::infinity();
26
  const LpBase::Value LpBase::NaN = std::numeric_limits<Value>::quiet_NaN();
25
  const LpBase::Value LpBase::INF =
26
    std::numeric_limits<LpBase::Value>::infinity();
27
  const LpBase::Value LpBase::NaN =
28
    std::numeric_limits<LpBase::Value>::quiet_NaN();
27 29

	
28 30
} //namespace lemon
Ignore white space 16777216 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2008
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_LP_BASE_H
20 20
#define LEMON_LP_BASE_H
21 21

	
22 22
#include<iostream>
23 23
#include<vector>
24 24
#include<map>
25 25
#include<limits>
26 26
#include<lemon/math.h>
27 27

	
28 28
#include<lemon/error.h>
29 29
#include<lemon/assert.h>
30 30

	
31 31
#include<lemon/core.h>
32 32
#include<lemon/bits/solver_bits.h>
33 33

	
34 34
///\file
35 35
///\brief The interface of the LP solver interface.
36 36
///\ingroup lp_group
37 37
namespace lemon {
38 38

	
39 39
  ///Common base class for LP and MIP solvers
40 40

	
41 41
  ///Usually this class is not used directly, please use one of the concrete
42 42
  ///implementations of the solver interface.
43 43
  ///\ingroup lp_group
44 44
  class LpBase {
45 45

	
46 46
  protected:
47 47

	
48 48
    _solver_bits::VarIndex rows;
49 49
    _solver_bits::VarIndex cols;
50 50

	
51 51
  public:
52 52

	
53 53
    ///Possible outcomes of an LP solving procedure
54 54
    enum SolveExitStatus {
55 55
      ///This means that the problem has been successfully solved: either
56 56
      ///an optimal solution has been found or infeasibility/unboundedness
57 57
      ///has been proved.
58 58
      SOLVED = 0,
59 59
      ///Any other case (including the case when some user specified
60 60
      ///limit has been exceeded)
61 61
      UNSOLVED = 1
62 62
    };
63 63

	
64 64
    ///Direction of the optimization
65 65
    enum Sense {
66 66
      /// Minimization
67 67
      MIN,
68 68
      /// Maximization
69 69
      MAX
70 70
    };
71 71

	
72 72
    ///The floating point type used by the solver
73 73
    typedef double Value;
74 74
    ///The infinity constant
75 75
    static const Value INF;
76 76
    ///The not a number constant
77 77
    static const Value NaN;
78 78

	
79 79
    friend class Col;
80 80
    friend class ColIt;
81 81
    friend class Row;
82 82
    friend class RowIt;
83 83

	
84 84
    ///Refer to a column of the LP.
85 85

	
86 86
    ///This type is used to refer to a column of the LP.
87 87
    ///
88 88
    ///Its value remains valid and correct even after the addition or erase of
89 89
    ///other columns.
90 90
    ///
91 91
    ///\note This class is similar to other Item types in LEMON, like
92 92
    ///Node and Arc types in digraph.
93 93
    class Col {
94 94
      friend class LpBase;
95 95
    protected:
96 96
      int _id;
97 97
      explicit Col(int id) : _id(id) {}
98 98
    public:
99 99
      typedef Value ExprValue;
100 100
      typedef True LpCol;
101 101
      /// Default constructor
102 102
      
103 103
      /// \warning The default constructor sets the Col to an
104 104
      /// undefined value.
105 105
      Col() {}
106 106
      /// Invalid constructor \& conversion.
107 107
      
108 108
      /// This constructor initializes the Col to be invalid.
109 109
      /// \sa Invalid for more details.      
110 110
      Col(const Invalid&) : _id(-1) {}
111 111
      /// Equality operator
112 112

	
113 113
      /// Two \ref Col "Col"s are equal if and only if they point to
114 114
      /// the same LP column or both are invalid.
115 115
      bool operator==(Col c) const  {return _id == c._id;}
116 116
      /// Inequality operator
117 117

	
118 118
      /// \sa operator==(Col c)
119 119
      ///
120 120
      bool operator!=(Col c) const  {return _id != c._id;}
121 121
      /// Artificial ordering operator.
122 122

	
123 123
      /// To allow the use of this object in std::map or similar
124 124
      /// associative container we require this.
125 125
      ///
126 126
      /// \note This operator only have to define some strict ordering of
127 127
      /// the items; this order has nothing to do with the iteration
128 128
      /// ordering of the items.
129 129
      bool operator<(Col c) const  {return _id < c._id;}
130 130
    };
131 131

	
132 132
    ///Iterator for iterate over the columns of an LP problem
133 133

	
134 134
    /// Its usage is quite simple, for example you can count the number
135 135
    /// of columns in an LP \c lp:
136 136
    ///\code
137 137
    /// int count=0;
138 138
    /// for (LpBase::ColIt c(lp); c!=INVALID; ++c) ++count;
139 139
    ///\endcode
140 140
    class ColIt : public Col {
141 141
      const LpBase *_solver;
142 142
    public:
143 143
      /// Default constructor
144 144
      
145 145
      /// \warning The default constructor sets the iterator
146 146
      /// to an undefined value.
147 147
      ColIt() {}
148 148
      /// Sets the iterator to the first Col
149 149
      
150 150
      /// Sets the iterator to the first Col.
151 151
      ///
152 152
      ColIt(const LpBase &solver) : _solver(&solver)
153 153
      {
154 154
        _solver->cols.firstItem(_id);
155 155
      }
156 156
      /// Invalid constructor \& conversion
157 157
      
158 158
      /// Initialize the iterator to be invalid.
159 159
      /// \sa Invalid for more details.
160 160
      ColIt(const Invalid&) : Col(INVALID) {}
161 161
      /// Next column
162 162
      
163 163
      /// Assign the iterator to the next column.
164 164
      ///
165 165
      ColIt &operator++()
166 166
      {
167 167
        _solver->cols.nextItem(_id);
168 168
        return *this;
169 169
      }
170 170
    };
171 171

	
172 172
    /// \brief Returns the ID of the column.
173 173
    static int id(const Col& col) { return col._id; }
174 174
    /// \brief Returns the column with the given ID.
175 175
    ///
176 176
    /// \pre The argument should be a valid column ID in the LP problem.
177 177
    static Col colFromId(int id) { return Col(id); }
178 178

	
179 179
    ///Refer to a row of the LP.
180 180

	
181 181
    ///This type is used to refer to a row of the LP.
182 182
    ///
183 183
    ///Its value remains valid and correct even after the addition or erase of
184 184
    ///other rows.
185 185
    ///
186 186
    ///\note This class is similar to other Item types in LEMON, like
187 187
    ///Node and Arc types in digraph.
188 188
    class Row {
189 189
      friend class LpBase;
190 190
    protected:
191 191
      int _id;
192 192
      explicit Row(int id) : _id(id) {}
193 193
    public:
194 194
      typedef Value ExprValue;
195 195
      typedef True LpRow;
196 196
      /// Default constructor
197 197
      
198 198
      /// \warning The default constructor sets the Row to an
199 199
      /// undefined value.
200 200
      Row() {}
201 201
      /// Invalid constructor \& conversion.
202 202
      
203 203
      /// This constructor initializes the Row to be invalid.
204 204
      /// \sa Invalid for more details.      
205 205
      Row(const Invalid&) : _id(-1) {}
206 206
      /// Equality operator
207 207

	
208 208
      /// Two \ref Row "Row"s are equal if and only if they point to
209 209
      /// the same LP row or both are invalid.
210 210
      bool operator==(Row r) const  {return _id == r._id;}
211 211
      /// Inequality operator
212 212
      
213 213
      /// \sa operator==(Row r)
214 214
      ///
215 215
      bool operator!=(Row r) const  {return _id != r._id;}
216 216
      /// Artificial ordering operator.
217 217

	
218 218
      /// To allow the use of this object in std::map or similar
219 219
      /// associative container we require this.
220 220
      ///
221 221
      /// \note This operator only have to define some strict ordering of
222 222
      /// the items; this order has nothing to do with the iteration
223 223
      /// ordering of the items.
224 224
      bool operator<(Row r) const  {return _id < r._id;}
225 225
    };
226 226

	
227 227
    ///Iterator for iterate over the rows of an LP problem
228 228

	
229 229
    /// Its usage is quite simple, for example you can count the number
230 230
    /// of rows in an LP \c lp:
231 231
    ///\code
232 232
    /// int count=0;
233 233
    /// for (LpBase::RowIt c(lp); c!=INVALID; ++c) ++count;
234 234
    ///\endcode
235 235
    class RowIt : public Row {
236 236
      const LpBase *_solver;
237 237
    public:
238 238
      /// Default constructor
239 239
      
240 240
      /// \warning The default constructor sets the iterator
241 241
      /// to an undefined value.
242 242
      RowIt() {}
243 243
      /// Sets the iterator to the first Row
244 244
      
245 245
      /// Sets the iterator to the first Row.
246 246
      ///
247 247
      RowIt(const LpBase &solver) : _solver(&solver)
248 248
      {
249 249
        _solver->rows.firstItem(_id);
250 250
      }
251 251
      /// Invalid constructor \& conversion
252 252
      
253 253
      /// Initialize the iterator to be invalid.
254 254
      /// \sa Invalid for more details.
255 255
      RowIt(const Invalid&) : Row(INVALID) {}
256 256
      /// Next row
257 257
      
258 258
      /// Assign the iterator to the next row.
259 259
      ///
260 260
      RowIt &operator++()
261 261
      {
262 262
        _solver->rows.nextItem(_id);
263 263
        return *this;
264 264
      }
265 265
    };
266 266

	
267 267
    /// \brief Returns the ID of the row.
268 268
    static int id(const Row& row) { return row._id; }
269 269
    /// \brief Returns the row with the given ID.
270 270
    ///
271 271
    /// \pre The argument should be a valid row ID in the LP problem.
272 272
    static Row rowFromId(int id) { return Row(id); }
273 273

	
274 274
  public:
275 275

	
276 276
    ///Linear expression of variables and a constant component
277 277

	
278 278
    ///This data structure stores a linear expression of the variables
279 279
    ///(\ref Col "Col"s) and also has a constant component.
280 280
    ///
281 281
    ///There are several ways to access and modify the contents of this
282 282
    ///container.
283 283
    ///\code
284 284
    ///e[v]=5;
285 285
    ///e[v]+=12;
286 286
    ///e.erase(v);
287 287
    ///\endcode
288 288
    ///or you can also iterate through its elements.
289 289
    ///\code
290 290
    ///double s=0;
291 291
    ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
292 292
    ///  s+=*i * primal(i);
293 293
    ///\endcode
294 294
    ///(This code computes the primal value of the expression).
295 295
    ///- Numbers (<tt>double</tt>'s)
296 296
    ///and variables (\ref Col "Col"s) directly convert to an
297 297
    ///\ref Expr and the usual linear operations are defined, so
298 298
    ///\code
299 299
    ///v+w
300 300
    ///2*v-3.12*(v-w/2)+2
301 301
    ///v*2.1+(3*v+(v*12+w+6)*3)/2
302 302
    ///\endcode
303 303
    ///are valid expressions.
304 304
    ///The usual assignment operations are also defined.
305 305
    ///\code
306 306
    ///e=v+w;
307 307
    ///e+=2*v-3.12*(v-w/2)+2;
308 308
    ///e*=3.4;
309 309
    ///e/=5;
310 310
    ///\endcode
311 311
    ///- The constant member can be set and read by dereference
312 312
    ///  operator (unary *)
313 313
    ///
314 314
    ///\code
315 315
    ///*e=12;
316 316
    ///double c=*e;
317 317
    ///\endcode
318 318
    ///
319 319
    ///\sa Constr
320 320
    class Expr {
321 321
      friend class LpBase;
322 322
    public:
323 323
      /// The key type of the expression
324 324
      typedef LpBase::Col Key;
325 325
      /// The value type of the expression
326 326
      typedef LpBase::Value Value;
327 327

	
328 328
    protected:
329 329
      Value const_comp;
330 330
      std::map<int, Value> comps;
331 331

	
332 332
    public:
333 333
      typedef True SolverExpr;
334 334
      /// Default constructor
335 335
      
336 336
      /// Construct an empty expression, the coefficients and
337 337
      /// the constant component are initialized to zero.
338 338
      Expr() : const_comp(0) {}
339 339
      /// Construct an expression from a column
340 340

	
341 341
      /// Construct an expression, which has a term with \c c variable
342 342
      /// and 1.0 coefficient.
343 343
      Expr(const Col &c) : const_comp(0) {
344 344
        typedef std::map<int, Value>::value_type pair_type;
345 345
        comps.insert(pair_type(id(c), 1));
346 346
      }
347 347
      /// Construct an expression from a constant
348 348

	
349 349
      /// Construct an expression, which's constant component is \c v.
350 350
      ///
351 351
      Expr(const Value &v) : const_comp(v) {}
352 352
      /// Returns the coefficient of the column
353 353
      Value operator[](const Col& c) const {
354 354
        std::map<int, Value>::const_iterator it=comps.find(id(c));
355 355
        if (it != comps.end()) {
356 356
          return it->second;
357 357
        } else {
358 358
          return 0;
359 359
        }
360 360
      }
361 361
      /// Returns the coefficient of the column
362 362
      Value& operator[](const Col& c) {
363 363
        return comps[id(c)];
364 364
      }
365 365
      /// Sets the coefficient of the column
366 366
      void set(const Col &c, const Value &v) {
367 367
        if (v != 0.0) {
368 368
          typedef std::map<int, Value>::value_type pair_type;
369 369
          comps.insert(pair_type(id(c), v));
370 370
        } else {
371 371
          comps.erase(id(c));
372 372
        }
373 373
      }
374 374
      /// Returns the constant component of the expression
375 375
      Value& operator*() { return const_comp; }
376 376
      /// Returns the constant component of the expression
377 377
      const Value& operator*() const { return const_comp; }
378 378
      /// \brief Removes the coefficients which's absolute value does
379 379
      /// not exceed \c epsilon. It also sets to zero the constant
380 380
      /// component, if it does not exceed epsilon in absolute value.
381 381
      void simplify(Value epsilon = 0.0) {
382 382
        std::map<int, Value>::iterator it=comps.begin();
383 383
        while (it != comps.end()) {
384 384
          std::map<int, Value>::iterator jt=it;
385 385
          ++jt;
386 386
          if (std::fabs((*it).second) <= epsilon) comps.erase(it);
387 387
          it=jt;
388 388
        }
389 389
        if (std::fabs(const_comp) <= epsilon) const_comp = 0;
390 390
      }
391 391

	
392 392
      void simplify(Value epsilon = 0.0) const {
393 393
        const_cast<Expr*>(this)->simplify(epsilon);
394 394
      }
395 395

	
396 396
      ///Sets all coefficients and the constant component to 0.
397 397
      void clear() {
398 398
        comps.clear();
399 399
        const_comp=0;
400 400
      }
401 401

	
402 402
      ///Compound assignment
403 403
      Expr &operator+=(const Expr &e) {
404 404
        for (std::map<int, Value>::const_iterator it=e.comps.begin();
405 405
             it!=e.comps.end(); ++it)
406 406
          comps[it->first]+=it->second;
407 407
        const_comp+=e.const_comp;
408 408
        return *this;
409 409
      }
410 410
      ///Compound assignment
411 411
      Expr &operator-=(const Expr &e) {
412 412
        for (std::map<int, Value>::const_iterator it=e.comps.begin();
413 413
             it!=e.comps.end(); ++it)
414 414
          comps[it->first]-=it->second;
415 415
        const_comp-=e.const_comp;
416 416
        return *this;
417 417
      }
418 418
      ///Multiply with a constant
419 419
      Expr &operator*=(const Value &v) {
420 420
        for (std::map<int, Value>::iterator it=comps.begin();
421 421
             it!=comps.end(); ++it)
422 422
          it->second*=v;
423 423
        const_comp*=v;
424 424
        return *this;
425 425
      }
426 426
      ///Division with a constant
427 427
      Expr &operator/=(const Value &c) {
428 428
        for (std::map<int, Value>::iterator it=comps.begin();
429 429
             it!=comps.end(); ++it)
430 430
          it->second/=c;
431 431
        const_comp/=c;
432 432
        return *this;
433 433
      }
434 434

	
435 435
      ///Iterator over the expression
436 436
      
437 437
      ///The iterator iterates over the terms of the expression. 
438 438
      /// 
439 439
      ///\code
440 440
      ///double s=0;
441 441
      ///for(LpBase::Expr::CoeffIt i(e);i!=INVALID;++i)
442 442
      ///  s+= *i * primal(i);
443 443
      ///\endcode
444 444
      class CoeffIt {
445 445
      private:
446 446

	
447 447
        std::map<int, Value>::iterator _it, _end;
448 448

	
449 449
      public:
450 450

	
451 451
        /// Sets the iterator to the first term
452 452
        
453 453
        /// Sets the iterator to the first term of the expression.
454 454
        ///
455 455
        CoeffIt(Expr& e)
456 456
          : _it(e.comps.begin()), _end(e.comps.end()){}
457 457

	
458 458
        /// Convert the iterator to the column of the term
459 459
        operator Col() const {
460 460
          return colFromId(_it->first);
461 461
        }
462 462

	
463 463
        /// Returns the coefficient of the term
464 464
        Value& operator*() { return _it->second; }
465 465

	
466 466
        /// Returns the coefficient of the term
467 467
        const Value& operator*() const { return _it->second; }
468 468
        /// Next term
469 469
        
470 470
        /// Assign the iterator to the next term.
471 471
        ///
472 472
        CoeffIt& operator++() { ++_it; return *this; }
473 473

	
474 474
        /// Equality operator
475 475
        bool operator==(Invalid) const { return _it == _end; }
476 476
        /// Inequality operator
477 477
        bool operator!=(Invalid) const { return _it != _end; }
478 478
      };
479 479

	
480 480
      /// Const iterator over the expression
481 481
      
482 482
      ///The iterator iterates over the terms of the expression. 
483 483
      /// 
484 484
      ///\code
485 485
      ///double s=0;
486 486
      ///for(LpBase::Expr::ConstCoeffIt i(e);i!=INVALID;++i)
487 487
      ///  s+=*i * primal(i);
488 488
      ///\endcode
489 489
      class ConstCoeffIt {
490 490
      private:
491 491

	
492 492
        std::map<int, Value>::const_iterator _it, _end;
493 493

	
494 494
      public:
495 495

	
496 496
        /// Sets the iterator to the first term
497 497
        
498 498
        /// Sets the iterator to the first term of the expression.
499 499
        ///
500 500
        ConstCoeffIt(const Expr& e)
501 501
          : _it(e.comps.begin()), _end(e.comps.end()){}
502 502

	
503 503
        /// Convert the iterator to the column of the term
504 504
        operator Col() const {
505 505
          return colFromId(_it->first);
506 506
        }
507 507

	
508 508
        /// Returns the coefficient of the term
509 509
        const Value& operator*() const { return _it->second; }
510 510

	
511 511
        /// Next term
512 512
        
513 513
        /// Assign the iterator to the next term.
514 514
        ///
515 515
        ConstCoeffIt& operator++() { ++_it; return *this; }
516 516

	
517 517
        /// Equality operator
518 518
        bool operator==(Invalid) const { return _it == _end; }
519 519
        /// Inequality operator
520 520
        bool operator!=(Invalid) const { return _it != _end; }
521 521
      };
522 522

	
523 523
    };
524 524

	
525 525
    ///Linear constraint
526 526

	
527 527
    ///This data stucture represents a linear constraint in the LP.
528 528
    ///Basically it is a linear expression with a lower or an upper bound
529 529
    ///(or both). These parts of the constraint can be obtained by the member
530 530
    ///functions \ref expr(), \ref lowerBound() and \ref upperBound(),
531 531
    ///respectively.
532 532
    ///There are two ways to construct a constraint.
533 533
    ///- You can set the linear expression and the bounds directly
534 534
    ///  by the functions above.
535 535
    ///- The operators <tt>\<=</tt>, <tt>==</tt> and  <tt>\>=</tt>
536 536
    ///  are defined between expressions, or even between constraints whenever
537 537
    ///  it makes sense. Therefore if \c e and \c f are linear expressions and
538 538
    ///  \c s and \c t are numbers, then the followings are valid expressions
539 539
    ///  and thus they can be used directly e.g. in \ref addRow() whenever
540 540
    ///  it makes sense.
541 541
    ///\code
542 542
    ///  e<=s
543 543
    ///  e<=f
544 544
    ///  e==f
545 545
    ///  s<=e<=t
546 546
    ///  e>=t
547 547
    ///\endcode
548 548
    ///\warning The validity of a constraint is checked only at run
549 549
    ///time, so e.g. \ref addRow(<tt>x[1]\<=x[2]<=5</tt>) will
550 550
    ///compile, but will fail an assertion.
551 551
    class Constr
552 552
    {
553 553
    public:
554 554
      typedef LpBase::Expr Expr;
555 555
      typedef Expr::Key Key;
556 556
      typedef Expr::Value Value;
557 557

	
558 558
    protected:
559 559
      Expr _expr;
560 560
      Value _lb,_ub;
561 561
    public:
562 562
      ///\e
563 563
      Constr() : _expr(), _lb(NaN), _ub(NaN) {}
564 564
      ///\e
565 565
      Constr(Value lb, const Expr &e, Value ub) :
566 566
        _expr(e), _lb(lb), _ub(ub) {}
567 567
      Constr(const Expr &e) :
568 568
        _expr(e), _lb(NaN), _ub(NaN) {}
569 569
      ///\e
570 570
      void clear()
571 571
      {
572 572
        _expr.clear();
573 573
        _lb=_ub=NaN;
574 574
      }
575 575

	
576 576
      ///Reference to the linear expression
577 577
      Expr &expr() { return _expr; }
578 578
      ///Cont reference to the linear expression
579 579
      const Expr &expr() const { return _expr; }
580 580
      ///Reference to the lower bound.
581 581

	
582 582
      ///\return
583 583
      ///- \ref INF "INF": the constraint is lower unbounded.
584 584
      ///- \ref NaN "NaN": lower bound has not been set.
585 585
      ///- finite number: the lower bound
586 586
      Value &lowerBound() { return _lb; }
587 587
      ///The const version of \ref lowerBound()
588 588
      const Value &lowerBound() const { return _lb; }
589 589
      ///Reference to the upper bound.
590 590

	
591 591
      ///\return
592 592
      ///- \ref INF "INF": the constraint is upper unbounded.
593 593
      ///- \ref NaN "NaN": upper bound has not been set.
594 594
      ///- finite number: the upper bound
595 595
      Value &upperBound() { return _ub; }
596 596
      ///The const version of \ref upperBound()
597 597
      const Value &upperBound() const { return _ub; }
598 598
      ///Is the constraint lower bounded?
599 599
      bool lowerBounded() const {
600
        return _lb != -INF && !isnan(_lb);
600
        return _lb != -INF && !isNaN(_lb);
601 601
      }
602 602
      ///Is the constraint upper bounded?
603 603
      bool upperBounded() const {
604
        return _ub != INF && !isnan(_ub);
604
        return _ub != INF && !isNaN(_ub);
605 605
      }
606 606

	
607 607
    };
608 608

	
609 609
    ///Linear expression of rows
610 610

	
611 611
    ///This data structure represents a column of the matrix,
612 612
    ///thas is it strores a linear expression of the dual variables
613 613
    ///(\ref Row "Row"s).
614 614
    ///
615 615
    ///There are several ways to access and modify the contents of this
616 616
    ///container.
617 617
    ///\code
618 618
    ///e[v]=5;
619 619
    ///e[v]+=12;
620 620
    ///e.erase(v);
621 621
    ///\endcode
622 622
    ///or you can also iterate through its elements.
623 623
    ///\code
624 624
    ///double s=0;
625 625
    ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
626 626
    ///  s+=*i;
627 627
    ///\endcode
628 628
    ///(This code computes the sum of all coefficients).
629 629
    ///- Numbers (<tt>double</tt>'s)
630 630
    ///and variables (\ref Row "Row"s) directly convert to an
631 631
    ///\ref DualExpr and the usual linear operations are defined, so
632 632
    ///\code
633 633
    ///v+w
634 634
    ///2*v-3.12*(v-w/2)
635 635
    ///v*2.1+(3*v+(v*12+w)*3)/2
636 636
    ///\endcode
637 637
    ///are valid \ref DualExpr dual expressions.
638 638
    ///The usual assignment operations are also defined.
639 639
    ///\code
640 640
    ///e=v+w;
641 641
    ///e+=2*v-3.12*(v-w/2);
642 642
    ///e*=3.4;
643 643
    ///e/=5;
644 644
    ///\endcode
645 645
    ///
646 646
    ///\sa Expr
647 647
    class DualExpr {
648 648
      friend class LpBase;
649 649
    public:
650 650
      /// The key type of the expression
651 651
      typedef LpBase::Row Key;
652 652
      /// The value type of the expression
653 653
      typedef LpBase::Value Value;
654 654

	
655 655
    protected:
656 656
      std::map<int, Value> comps;
657 657

	
658 658
    public:
659 659
      typedef True SolverExpr;
660 660
      /// Default constructor
661 661
      
662 662
      /// Construct an empty expression, the coefficients are
663 663
      /// initialized to zero.
664 664
      DualExpr() {}
665 665
      /// Construct an expression from a row
666 666

	
667 667
      /// Construct an expression, which has a term with \c r dual
668 668
      /// variable and 1.0 coefficient.
669 669
      DualExpr(const Row &r) {
670 670
        typedef std::map<int, Value>::value_type pair_type;
671 671
        comps.insert(pair_type(id(r), 1));
672 672
      }
673 673
      /// Returns the coefficient of the row
674 674
      Value operator[](const Row& r) const {
675 675
        std::map<int, Value>::const_iterator it = comps.find(id(r));
676 676
        if (it != comps.end()) {
677 677
          return it->second;
678 678
        } else {
679 679
          return 0;
680 680
        }
681 681
      }
682 682
      /// Returns the coefficient of the row
683 683
      Value& operator[](const Row& r) {
684 684
        return comps[id(r)];
685 685
      }
686 686
      /// Sets the coefficient of the row
687 687
      void set(const Row &r, const Value &v) {
688 688
        if (v != 0.0) {
689 689
          typedef std::map<int, Value>::value_type pair_type;
690 690
          comps.insert(pair_type(id(r), v));
691 691
        } else {
692 692
          comps.erase(id(r));
693 693
        }
694 694
      }
695 695
      /// \brief Removes the coefficients which's absolute value does
696 696
      /// not exceed \c epsilon. 
697 697
      void simplify(Value epsilon = 0.0) {
698 698
        std::map<int, Value>::iterator it=comps.begin();
699 699
        while (it != comps.end()) {
700 700
          std::map<int, Value>::iterator jt=it;
701 701
          ++jt;
702 702
          if (std::fabs((*it).second) <= epsilon) comps.erase(it);
703 703
          it=jt;
704 704
        }
705 705
      }
706 706

	
707 707
      void simplify(Value epsilon = 0.0) const {
708 708
        const_cast<DualExpr*>(this)->simplify(epsilon);
709 709
      }
710 710

	
711 711
      ///Sets all coefficients to 0.
712 712
      void clear() {
713 713
        comps.clear();
714 714
      }
715 715
      ///Compound assignment
716 716
      DualExpr &operator+=(const DualExpr &e) {
717 717
        for (std::map<int, Value>::const_iterator it=e.comps.begin();
718 718
             it!=e.comps.end(); ++it)
719 719
          comps[it->first]+=it->second;
720 720
        return *this;
721 721
      }
722 722
      ///Compound assignment
723 723
      DualExpr &operator-=(const DualExpr &e) {
724 724
        for (std::map<int, Value>::const_iterator it=e.comps.begin();
725 725
             it!=e.comps.end(); ++it)
726 726
          comps[it->first]-=it->second;
727 727
        return *this;
728 728
      }
729 729
      ///Multiply with a constant
730 730
      DualExpr &operator*=(const Value &v) {
731 731
        for (std::map<int, Value>::iterator it=comps.begin();
732 732
             it!=comps.end(); ++it)
733 733
          it->second*=v;
734 734
        return *this;
735 735
      }
736 736
      ///Division with a constant
737 737
      DualExpr &operator/=(const Value &v) {
738 738
        for (std::map<int, Value>::iterator it=comps.begin();
739 739
             it!=comps.end(); ++it)
740 740
          it->second/=v;
741 741
        return *this;
742 742
      }
743 743

	
744 744
      ///Iterator over the expression
745 745
      
746 746
      ///The iterator iterates over the terms of the expression. 
747 747
      /// 
748 748
      ///\code
749 749
      ///double s=0;
750 750
      ///for(LpBase::DualExpr::CoeffIt i(e);i!=INVALID;++i)
751 751
      ///  s+= *i * dual(i);
752 752
      ///\endcode
753 753
      class CoeffIt {
754 754
      private:
755 755

	
756 756
        std::map<int, Value>::iterator _it, _end;
757 757

	
758 758
      public:
759 759

	
760 760
        /// Sets the iterator to the first term
761 761
        
762 762
        /// Sets the iterator to the first term of the expression.
763 763
        ///
764 764
        CoeffIt(DualExpr& e)
765 765
          : _it(e.comps.begin()), _end(e.comps.end()){}
766 766

	
767 767
        /// Convert the iterator to the row of the term
768 768
        operator Row() const {
769 769
          return rowFromId(_it->first);
770 770
        }
771 771

	
772 772
        /// Returns the coefficient of the term
773 773
        Value& operator*() { return _it->second; }
774 774

	
775 775
        /// Returns the coefficient of the term
776 776
        const Value& operator*() const { return _it->second; }
777 777

	
778 778
        /// Next term
779 779
        
780 780
        /// Assign the iterator to the next term.
781 781
        ///
782 782
        CoeffIt& operator++() { ++_it; return *this; }
783 783

	
784 784
        /// Equality operator
785 785
        bool operator==(Invalid) const { return _it == _end; }
786 786
        /// Inequality operator
787 787
        bool operator!=(Invalid) const { return _it != _end; }
788 788
      };
789 789

	
790 790
      ///Iterator over the expression
791 791
      
792 792
      ///The iterator iterates over the terms of the expression. 
793 793
      /// 
794 794
      ///\code
795 795
      ///double s=0;
796 796
      ///for(LpBase::DualExpr::ConstCoeffIt i(e);i!=INVALID;++i)
797 797
      ///  s+= *i * dual(i);
798 798
      ///\endcode
799 799
      class ConstCoeffIt {
800 800
      private:
801 801

	
802 802
        std::map<int, Value>::const_iterator _it, _end;
803 803

	
804 804
      public:
805 805

	
806 806
        /// Sets the iterator to the first term
807 807
        
808 808
        /// Sets the iterator to the first term of the expression.
809 809
        ///
810 810
        ConstCoeffIt(const DualExpr& e)
811 811
          : _it(e.comps.begin()), _end(e.comps.end()){}
812 812

	
813 813
        /// Convert the iterator to the row of the term
814 814
        operator Row() const {
815 815
          return rowFromId(_it->first);
816 816
        }
817 817

	
818 818
        /// Returns the coefficient of the term
819 819
        const Value& operator*() const { return _it->second; }
820 820

	
821 821
        /// Next term
822 822
        
823 823
        /// Assign the iterator to the next term.
824 824
        ///
825 825
        ConstCoeffIt& operator++() { ++_it; return *this; }
826 826

	
827 827
        /// Equality operator
828 828
        bool operator==(Invalid) const { return _it == _end; }
829 829
        /// Inequality operator
830 830
        bool operator!=(Invalid) const { return _it != _end; }
831 831
      };
832 832
    };
833 833

	
834 834

	
835 835
  protected:
836 836

	
837 837
    class InsertIterator {
838 838
    private:
839 839

	
840 840
      std::map<int, Value>& _host;
841 841
      const _solver_bits::VarIndex& _index;
842 842

	
843 843
    public:
844 844

	
845 845
      typedef std::output_iterator_tag iterator_category;
846 846
      typedef void difference_type;
847 847
      typedef void value_type;
848 848
      typedef void reference;
849 849
      typedef void pointer;
850 850

	
851 851
      InsertIterator(std::map<int, Value>& host,
852 852
                   const _solver_bits::VarIndex& index)
853 853
        : _host(host), _index(index) {}
854 854

	
855 855
      InsertIterator& operator=(const std::pair<int, Value>& value) {
856 856
        typedef std::map<int, Value>::value_type pair_type;
857 857
        _host.insert(pair_type(_index[value.first], value.second));
858 858
        return *this;
859 859
      }
860 860

	
861 861
      InsertIterator& operator*() { return *this; }
862 862
      InsertIterator& operator++() { return *this; }
863 863
      InsertIterator operator++(int) { return *this; }
864 864

	
865 865
    };
866 866

	
867 867
    class ExprIterator {
868 868
    private:
869 869
      std::map<int, Value>::const_iterator _host_it;
870 870
      const _solver_bits::VarIndex& _index;
871 871
    public:
872 872

	
873 873
      typedef std::bidirectional_iterator_tag iterator_category;
874 874
      typedef std::ptrdiff_t difference_type;
875 875
      typedef const std::pair<int, Value> value_type;
876 876
      typedef value_type reference;
877 877

	
878 878
      class pointer {
879 879
      public:
880 880
        pointer(value_type& _value) : value(_value) {}
881 881
        value_type* operator->() { return &value; }
882 882
      private:
883 883
        value_type value;
884 884
      };
885 885

	
886 886
      ExprIterator(const std::map<int, Value>::const_iterator& host_it,
887 887
                   const _solver_bits::VarIndex& index)
888 888
        : _host_it(host_it), _index(index) {}
889 889

	
890 890
      reference operator*() {
891 891
        return std::make_pair(_index(_host_it->first), _host_it->second);
892 892
      }
893 893

	
894 894
      pointer operator->() {
895 895
        return pointer(operator*());
896 896
      }
897 897

	
898 898
      ExprIterator& operator++() { ++_host_it; return *this; }
899 899
      ExprIterator operator++(int) {
900 900
        ExprIterator tmp(*this); ++_host_it; return tmp;
901 901
      }
902 902

	
903 903
      ExprIterator& operator--() { --_host_it; return *this; }
904 904
      ExprIterator operator--(int) {
905 905
        ExprIterator tmp(*this); --_host_it; return tmp;
906 906
      }
907 907

	
908 908
      bool operator==(const ExprIterator& it) const {
909 909
        return _host_it == it._host_it;
910 910
      }
911 911

	
912 912
      bool operator!=(const ExprIterator& it) const {
913 913
        return _host_it != it._host_it;
914 914
      }
915 915

	
916 916
    };
917 917

	
918 918
  protected:
919 919

	
920 920
    //Abstract virtual functions
921 921
    virtual LpBase* _newSolver() const = 0;
922 922
    virtual LpBase* _cloneSolver() const = 0;
923 923

	
924 924
    virtual int _addColId(int col) { return cols.addIndex(col); }
925 925
    virtual int _addRowId(int row) { return rows.addIndex(row); }
926 926

	
927 927
    virtual void _eraseColId(int col) { cols.eraseIndex(col); }
928 928
    virtual void _eraseRowId(int row) { rows.eraseIndex(row); }
929 929

	
930 930
    virtual int _addCol() = 0;
931 931
    virtual int _addRow() = 0;
932 932

	
933 933
    virtual void _eraseCol(int col) = 0;
934 934
    virtual void _eraseRow(int row) = 0;
935 935

	
936 936
    virtual void _getColName(int col, std::string& name) const = 0;
937 937
    virtual void _setColName(int col, const std::string& name) = 0;
938 938
    virtual int _colByName(const std::string& name) const = 0;
939 939

	
940 940
    virtual void _getRowName(int row, std::string& name) const = 0;
941 941
    virtual void _setRowName(int row, const std::string& name) = 0;
942 942
    virtual int _rowByName(const std::string& name) const = 0;
943 943

	
944 944
    virtual void _setRowCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
945 945
    virtual void _getRowCoeffs(int i, InsertIterator b) const = 0;
946 946

	
947 947
    virtual void _setColCoeffs(int i, ExprIterator b, ExprIterator e) = 0;
948 948
    virtual void _getColCoeffs(int i, InsertIterator b) const = 0;
949 949

	
950 950
    virtual void _setCoeff(int row, int col, Value value) = 0;
951 951
    virtual Value _getCoeff(int row, int col) const = 0;
952 952

	
953 953
    virtual void _setColLowerBound(int i, Value value) = 0;
954 954
    virtual Value _getColLowerBound(int i) const = 0;
955 955

	
956 956
    virtual void _setColUpperBound(int i, Value value) = 0;
957 957
    virtual Value _getColUpperBound(int i) const = 0;
958 958

	
959 959
    virtual void _setRowLowerBound(int i, Value value) = 0;
960 960
    virtual Value _getRowLowerBound(int i) const = 0;
961 961

	
962 962
    virtual void _setRowUpperBound(int i, Value value) = 0;
963 963
    virtual Value _getRowUpperBound(int i) const = 0;
964 964

	
965 965
    virtual void _setObjCoeffs(ExprIterator b, ExprIterator e) = 0;
966 966
    virtual void _getObjCoeffs(InsertIterator b) const = 0;
967 967

	
968 968
    virtual void _setObjCoeff(int i, Value obj_coef) = 0;
969 969
    virtual Value _getObjCoeff(int i) const = 0;
970 970

	
971 971
    virtual void _setSense(Sense) = 0;
972 972
    virtual Sense _getSense() const = 0;
973 973

	
974 974
    virtual void _clear() = 0;
975 975

	
976 976
    virtual const char* _solverName() const = 0;
977 977

	
978 978
    //Own protected stuff
979 979

	
980 980
    //Constant component of the objective function
981 981
    Value obj_const_comp;
982 982

	
983 983
    LpBase() : rows(), cols(), obj_const_comp(0) {}
984 984

	
985 985
  public:
986 986

	
987 987
    /// Virtual destructor
988 988
    virtual ~LpBase() {}
989 989

	
990 990
    ///Creates a new LP problem
991 991
    LpBase* newSolver() {return _newSolver();}
992 992
    ///Makes a copy of the LP problem
993 993
    LpBase* cloneSolver() {return _cloneSolver();}
994 994

	
995 995
    ///Gives back the name of the solver.
996 996
    const char* solverName() const {return _solverName();}
997 997

	
998 998
    ///\name Build up and modify the LP
999 999

	
1000 1000
    ///@{
1001 1001

	
1002 1002
    ///Add a new empty column (i.e a new variable) to the LP
1003 1003
    Col addCol() { Col c; c._id = _addColId(_addCol()); return c;}
1004 1004

	
1005 1005
    ///\brief Adds several new columns (i.e variables) at once
1006 1006
    ///
1007 1007
    ///This magic function takes a container as its argument and fills
1008 1008
    ///its elements with new columns (i.e. variables)
1009 1009
    ///\param t can be
1010 1010
    ///- a standard STL compatible iterable container with
1011 1011
    ///\ref Col as its \c values_type like
1012 1012
    ///\code
1013 1013
    ///std::vector<LpBase::Col>
1014 1014
    ///std::list<LpBase::Col>
1015 1015
    ///\endcode
1016 1016
    ///- a standard STL compatible iterable container with
1017 1017
    ///\ref Col as its \c mapped_type like
1018 1018
    ///\code
1019 1019
    ///std::map<AnyType,LpBase::Col>
1020 1020
    ///\endcode
1021 1021
    ///- an iterable lemon \ref concepts::WriteMap "write map" like
1022 1022
    ///\code
1023 1023
    ///ListGraph::NodeMap<LpBase::Col>
1024 1024
    ///ListGraph::ArcMap<LpBase::Col>
1025 1025
    ///\endcode
1026 1026
    ///\return The number of the created column.
1027 1027
#ifdef DOXYGEN
1028 1028
    template<class T>
1029 1029
    int addColSet(T &t) { return 0;}
1030 1030
#else
1031 1031
    template<class T>
1032 1032
    typename enable_if<typename T::value_type::LpCol,int>::type
1033 1033
    addColSet(T &t,dummy<0> = 0) {
1034 1034
      int s=0;
1035 1035
      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addCol();s++;}
1036 1036
      return s;
1037 1037
    }
1038 1038
    template<class T>
1039 1039
    typename enable_if<typename T::value_type::second_type::LpCol,
1040 1040
                       int>::type
1041 1041
    addColSet(T &t,dummy<1> = 1) {
1042 1042
      int s=0;
1043 1043
      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1044 1044
        i->second=addCol();
1045 1045
        s++;
1046 1046
      }
1047 1047
      return s;
1048 1048
    }
1049 1049
    template<class T>
1050 1050
    typename enable_if<typename T::MapIt::Value::LpCol,
1051 1051
                       int>::type
1052 1052
    addColSet(T &t,dummy<2> = 2) {
1053 1053
      int s=0;
1054 1054
      for(typename T::MapIt i(t); i!=INVALID; ++i)
1055 1055
        {
1056 1056
          i.set(addCol());
1057 1057
          s++;
1058 1058
        }
1059 1059
      return s;
1060 1060
    }
1061 1061
#endif
1062 1062

	
1063 1063
    ///Set a column (i.e a dual constraint) of the LP
1064 1064

	
1065 1065
    ///\param c is the column to be modified
1066 1066
    ///\param e is a dual linear expression (see \ref DualExpr)
1067 1067
    ///a better one.
1068 1068
    void col(Col c, const DualExpr &e) {
1069 1069
      e.simplify();
1070 1070
      _setColCoeffs(cols(id(c)), ExprIterator(e.comps.begin(), rows),
1071 1071
                    ExprIterator(e.comps.end(), rows));
1072 1072
    }
1073 1073

	
1074 1074
    ///Get a column (i.e a dual constraint) of the LP
1075 1075

	
1076 1076
    ///\param c is the column to get
1077 1077
    ///\return the dual expression associated to the column
1078 1078
    DualExpr col(Col c) const {
1079 1079
      DualExpr e;
1080 1080
      _getColCoeffs(cols(id(c)), InsertIterator(e.comps, rows));
1081 1081
      return e;
1082 1082
    }
1083 1083

	
1084 1084
    ///Add a new column to the LP
1085 1085

	
1086 1086
    ///\param e is a dual linear expression (see \ref DualExpr)
1087 1087
    ///\param o is the corresponding component of the objective
1088 1088
    ///function. It is 0 by default.
1089 1089
    ///\return The created column.
1090 1090
    Col addCol(const DualExpr &e, Value o = 0) {
1091 1091
      Col c=addCol();
1092 1092
      col(c,e);
1093 1093
      objCoeff(c,o);
1094 1094
      return c;
1095 1095
    }
1096 1096

	
1097 1097
    ///Add a new empty row (i.e a new constraint) to the LP
1098 1098

	
1099 1099
    ///This function adds a new empty row (i.e a new constraint) to the LP.
1100 1100
    ///\return The created row
1101 1101
    Row addRow() { Row r; r._id = _addRowId(_addRow()); return r;}
1102 1102

	
1103 1103
    ///\brief Add several new rows (i.e constraints) at once
1104 1104
    ///
1105 1105
    ///This magic function takes a container as its argument and fills
1106 1106
    ///its elements with new row (i.e. variables)
1107 1107
    ///\param t can be
1108 1108
    ///- a standard STL compatible iterable container with
1109 1109
    ///\ref Row as its \c values_type like
1110 1110
    ///\code
1111 1111
    ///std::vector<LpBase::Row>
1112 1112
    ///std::list<LpBase::Row>
1113 1113
    ///\endcode
1114 1114
    ///- a standard STL compatible iterable container with
1115 1115
    ///\ref Row as its \c mapped_type like
1116 1116
    ///\code
1117 1117
    ///std::map<AnyType,LpBase::Row>
1118 1118
    ///\endcode
1119 1119
    ///- an iterable lemon \ref concepts::WriteMap "write map" like
1120 1120
    ///\code
1121 1121
    ///ListGraph::NodeMap<LpBase::Row>
1122 1122
    ///ListGraph::ArcMap<LpBase::Row>
1123 1123
    ///\endcode
1124 1124
    ///\return The number of rows created.
1125 1125
#ifdef DOXYGEN
1126 1126
    template<class T>
1127 1127
    int addRowSet(T &t) { return 0;}
1128 1128
#else
1129 1129
    template<class T>
1130 1130
    typename enable_if<typename T::value_type::LpRow,int>::type
1131 1131
    addRowSet(T &t, dummy<0> = 0) {
1132 1132
      int s=0;
1133 1133
      for(typename T::iterator i=t.begin();i!=t.end();++i) {*i=addRow();s++;}
1134 1134
      return s;
1135 1135
    }
1136 1136
    template<class T>
1137 1137
    typename enable_if<typename T::value_type::second_type::LpRow, int>::type
1138 1138
    addRowSet(T &t, dummy<1> = 1) {
1139 1139
      int s=0;
1140 1140
      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1141 1141
        i->second=addRow();
1142 1142
        s++;
1143 1143
      }
1144 1144
      return s;
1145 1145
    }
1146 1146
    template<class T>
1147 1147
    typename enable_if<typename T::MapIt::Value::LpRow, int>::type
1148 1148
    addRowSet(T &t, dummy<2> = 2) {
1149 1149
      int s=0;
1150 1150
      for(typename T::MapIt i(t); i!=INVALID; ++i)
1151 1151
        {
1152 1152
          i.set(addRow());
1153 1153
          s++;
1154 1154
        }
1155 1155
      return s;
1156 1156
    }
1157 1157
#endif
1158 1158

	
1159 1159
    ///Set a row (i.e a constraint) of the LP
1160 1160

	
1161 1161
    ///\param r is the row to be modified
1162 1162
    ///\param l is lower bound (-\ref INF means no bound)
1163 1163
    ///\param e is a linear expression (see \ref Expr)
1164 1164
    ///\param u is the upper bound (\ref INF means no bound)
1165 1165
    void row(Row r, Value l, const Expr &e, Value u) {
1166 1166
      e.simplify();
1167 1167
      _setRowCoeffs(rows(id(r)), ExprIterator(e.comps.begin(), cols),
1168 1168
                    ExprIterator(e.comps.end(), cols));
1169 1169
      _setRowLowerBound(rows(id(r)),l - *e);
1170 1170
      _setRowUpperBound(rows(id(r)),u - *e);
1171 1171
    }
1172 1172

	
1173 1173
    ///Set a row (i.e a constraint) of the LP
1174 1174

	
1175 1175
    ///\param r is the row to be modified
1176 1176
    ///\param c is a linear expression (see \ref Constr)
1177 1177
    void row(Row r, const Constr &c) {
1178 1178
      row(r, c.lowerBounded()?c.lowerBound():-INF,
1179 1179
          c.expr(), c.upperBounded()?c.upperBound():INF);
1180 1180
    }
1181 1181

	
1182 1182

	
1183 1183
    ///Get a row (i.e a constraint) of the LP
1184 1184

	
1185 1185
    ///\param r is the row to get
1186 1186
    ///\return the expression associated to the row
1187 1187
    Expr row(Row r) const {
1188 1188
      Expr e;
1189 1189
      _getRowCoeffs(rows(id(r)), InsertIterator(e.comps, cols));
1190 1190
      return e;
1191 1191
    }
1192 1192

	
1193 1193
    ///Add a new row (i.e a new constraint) to the LP
1194 1194

	
1195 1195
    ///\param l is the lower bound (-\ref INF means no bound)
1196 1196
    ///\param e is a linear expression (see \ref Expr)
1197 1197
    ///\param u is the upper bound (\ref INF means no bound)
1198 1198
    ///\return The created row.
1199 1199
    Row addRow(Value l,const Expr &e, Value u) {
1200 1200
      Row r=addRow();
1201 1201
      row(r,l,e,u);
1202 1202
      return r;
1203 1203
    }
1204 1204

	
1205 1205
    ///Add a new row (i.e a new constraint) to the LP
1206 1206

	
1207 1207
    ///\param c is a linear expression (see \ref Constr)
1208 1208
    ///\return The created row.
1209 1209
    Row addRow(const Constr &c) {
1210 1210
      Row r=addRow();
1211 1211
      row(r,c);
1212 1212
      return r;
1213 1213
    }
1214 1214
    ///Erase a column (i.e a variable) from the LP
1215 1215

	
1216 1216
    ///\param c is the column to be deleted
1217 1217
    void erase(Col c) {
1218 1218
      _eraseCol(cols(id(c)));
1219 1219
      _eraseColId(cols(id(c)));
1220 1220
    }
1221 1221
    ///Erase a row (i.e a constraint) from the LP
1222 1222

	
1223 1223
    ///\param r is the row to be deleted
1224 1224
    void erase(Row r) {
1225 1225
      _eraseRow(rows(id(r)));
1226 1226
      _eraseRowId(rows(id(r)));
1227 1227
    }
1228 1228

	
1229 1229
    /// Get the name of a column
1230 1230

	
1231 1231
    ///\param c is the coresponding column
1232 1232
    ///\return The name of the colunm
1233 1233
    std::string colName(Col c) const {
1234 1234
      std::string name;
1235 1235
      _getColName(cols(id(c)), name);
1236 1236
      return name;
1237 1237
    }
1238 1238

	
1239 1239
    /// Set the name of a column
1240 1240

	
1241 1241
    ///\param c is the coresponding column
1242 1242
    ///\param name The name to be given
1243 1243
    void colName(Col c, const std::string& name) {
1244 1244
      _setColName(cols(id(c)), name);
1245 1245
    }
1246 1246

	
1247 1247
    /// Get the column by its name
1248 1248

	
1249 1249
    ///\param name The name of the column
1250 1250
    ///\return the proper column or \c INVALID
1251 1251
    Col colByName(const std::string& name) const {
1252 1252
      int k = _colByName(name);
1253 1253
      return k != -1 ? Col(cols[k]) : Col(INVALID);
1254 1254
    }
1255 1255

	
1256 1256
    /// Get the name of a row
1257 1257

	
1258 1258
    ///\param r is the coresponding row
1259 1259
    ///\return The name of the row
1260 1260
    std::string rowName(Row r) const {
1261 1261
      std::string name;
1262 1262
      _getRowName(rows(id(r)), name);
1263 1263
      return name;
1264 1264
    }
1265 1265

	
1266 1266
    /// Set the name of a row
1267 1267

	
1268 1268
    ///\param r is the coresponding row
1269 1269
    ///\param name The name to be given
1270 1270
    void rowName(Row r, const std::string& name) {
1271 1271
      _setRowName(rows(id(r)), name);
1272 1272
    }
1273 1273

	
1274 1274
    /// Get the row by its name
1275 1275

	
1276 1276
    ///\param name The name of the row
1277 1277
    ///\return the proper row or \c INVALID
1278 1278
    Row rowByName(const std::string& name) const {
1279 1279
      int k = _rowByName(name);
1280 1280
      return k != -1 ? Row(rows[k]) : Row(INVALID);
1281 1281
    }
1282 1282

	
1283 1283
    /// Set an element of the coefficient matrix of the LP
1284 1284

	
1285 1285
    ///\param r is the row of the element to be modified
1286 1286
    ///\param c is the column of the element to be modified
1287 1287
    ///\param val is the new value of the coefficient
1288 1288
    void coeff(Row r, Col c, Value val) {
1289 1289
      _setCoeff(rows(id(r)),cols(id(c)), val);
1290 1290
    }
1291 1291

	
1292 1292
    /// Get an element of the coefficient matrix of the LP
1293 1293

	
1294 1294
    ///\param r is the row of the element
1295 1295
    ///\param c is the column of the element
1296 1296
    ///\return the corresponding coefficient
1297 1297
    Value coeff(Row r, Col c) const {
1298 1298
      return _getCoeff(rows(id(r)),cols(id(c)));
1299 1299
    }
1300 1300

	
1301 1301
    /// Set the lower bound of a column (i.e a variable)
1302 1302

	
1303 1303
    /// The lower bound of a variable (column) has to be given by an
1304 1304
    /// extended number of type Value, i.e. a finite number of type
1305 1305
    /// Value or -\ref INF.
1306 1306
    void colLowerBound(Col c, Value value) {
1307 1307
      _setColLowerBound(cols(id(c)),value);
1308 1308
    }
1309 1309

	
1310 1310
    /// Get the lower bound of a column (i.e a variable)
1311 1311

	
1312 1312
    /// This function returns the lower bound for column (variable) \c c
1313 1313
    /// (this might be -\ref INF as well).
1314 1314
    ///\return The lower bound for column \c c
1315 1315
    Value colLowerBound(Col c) const {
1316 1316
      return _getColLowerBound(cols(id(c)));
1317 1317
    }
1318 1318

	
1319 1319
    ///\brief Set the lower bound of  several columns
1320 1320
    ///(i.e variables) at once
1321 1321
    ///
1322 1322
    ///This magic function takes a container as its argument
1323 1323
    ///and applies the function on all of its elements.
1324 1324
    ///The lower bound of a variable (column) has to be given by an
1325 1325
    ///extended number of type Value, i.e. a finite number of type
1326 1326
    ///Value or -\ref INF.
1327 1327
#ifdef DOXYGEN
1328 1328
    template<class T>
1329 1329
    void colLowerBound(T &t, Value value) { return 0;}
1330 1330
#else
1331 1331
    template<class T>
1332 1332
    typename enable_if<typename T::value_type::LpCol,void>::type
1333 1333
    colLowerBound(T &t, Value value,dummy<0> = 0) {
1334 1334
      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1335 1335
        colLowerBound(*i, value);
1336 1336
      }
1337 1337
    }
1338 1338
    template<class T>
1339 1339
    typename enable_if<typename T::value_type::second_type::LpCol,
1340 1340
                       void>::type
1341 1341
    colLowerBound(T &t, Value value,dummy<1> = 1) {
1342 1342
      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1343 1343
        colLowerBound(i->second, value);
1344 1344
      }
1345 1345
    }
1346 1346
    template<class T>
1347 1347
    typename enable_if<typename T::MapIt::Value::LpCol,
1348 1348
                       void>::type
1349 1349
    colLowerBound(T &t, Value value,dummy<2> = 2) {
1350 1350
      for(typename T::MapIt i(t); i!=INVALID; ++i){
1351 1351
        colLowerBound(*i, value);
1352 1352
      }
1353 1353
    }
1354 1354
#endif
1355 1355

	
1356 1356
    /// Set the upper bound of a column (i.e a variable)
1357 1357

	
1358 1358
    /// The upper bound of a variable (column) has to be given by an
1359 1359
    /// extended number of type Value, i.e. a finite number of type
1360 1360
    /// Value or \ref INF.
1361 1361
    void colUpperBound(Col c, Value value) {
1362 1362
      _setColUpperBound(cols(id(c)),value);
1363 1363
    };
1364 1364

	
1365 1365
    /// Get the upper bound of a column (i.e a variable)
1366 1366

	
1367 1367
    /// This function returns the upper bound for column (variable) \c c
1368 1368
    /// (this might be \ref INF as well).
1369 1369
    /// \return The upper bound for column \c c
1370 1370
    Value colUpperBound(Col c) const {
1371 1371
      return _getColUpperBound(cols(id(c)));
1372 1372
    }
1373 1373

	
1374 1374
    ///\brief Set the upper bound of  several columns
1375 1375
    ///(i.e variables) at once
1376 1376
    ///
1377 1377
    ///This magic function takes a container as its argument
1378 1378
    ///and applies the function on all of its elements.
1379 1379
    ///The upper bound of a variable (column) has to be given by an
1380 1380
    ///extended number of type Value, i.e. a finite number of type
1381 1381
    ///Value or \ref INF.
1382 1382
#ifdef DOXYGEN
1383 1383
    template<class T>
1384 1384
    void colUpperBound(T &t, Value value) { return 0;}
1385 1385
#else
1386 1386
    template<class T>
1387 1387
    typename enable_if<typename T::value_type::LpCol,void>::type
1388 1388
    colUpperBound(T &t, Value value,dummy<0> = 0) {
1389 1389
      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1390 1390
        colUpperBound(*i, value);
1391 1391
      }
1392 1392
    }
1393 1393
    template<class T>
1394 1394
    typename enable_if<typename T::value_type::second_type::LpCol,
1395 1395
                       void>::type
1396 1396
    colUpperBound(T &t, Value value,dummy<1> = 1) {
1397 1397
      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1398 1398
        colUpperBound(i->second, value);
1399 1399
      }
1400 1400
    }
1401 1401
    template<class T>
1402 1402
    typename enable_if<typename T::MapIt::Value::LpCol,
1403 1403
                       void>::type
1404 1404
    colUpperBound(T &t, Value value,dummy<2> = 2) {
1405 1405
      for(typename T::MapIt i(t); i!=INVALID; ++i){
1406 1406
        colUpperBound(*i, value);
1407 1407
      }
1408 1408
    }
1409 1409
#endif
1410 1410

	
1411 1411
    /// Set the lower and the upper bounds of a column (i.e a variable)
1412 1412

	
1413 1413
    /// The lower and the upper bounds of
1414 1414
    /// a variable (column) have to be given by an
1415 1415
    /// extended number of type Value, i.e. a finite number of type
1416 1416
    /// Value, -\ref INF or \ref INF.
1417 1417
    void colBounds(Col c, Value lower, Value upper) {
1418 1418
      _setColLowerBound(cols(id(c)),lower);
1419 1419
      _setColUpperBound(cols(id(c)),upper);
1420 1420
    }
1421 1421

	
1422 1422
    ///\brief Set the lower and the upper bound of several columns
1423 1423
    ///(i.e variables) at once
1424 1424
    ///
1425 1425
    ///This magic function takes a container as its argument
1426 1426
    ///and applies the function on all of its elements.
1427 1427
    /// The lower and the upper bounds of
1428 1428
    /// a variable (column) have to be given by an
1429 1429
    /// extended number of type Value, i.e. a finite number of type
1430 1430
    /// Value, -\ref INF or \ref INF.
1431 1431
#ifdef DOXYGEN
1432 1432
    template<class T>
1433 1433
    void colBounds(T &t, Value lower, Value upper) { return 0;}
1434 1434
#else
1435 1435
    template<class T>
1436 1436
    typename enable_if<typename T::value_type::LpCol,void>::type
1437 1437
    colBounds(T &t, Value lower, Value upper,dummy<0> = 0) {
1438 1438
      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1439 1439
        colBounds(*i, lower, upper);
1440 1440
      }
1441 1441
    }
1442 1442
    template<class T>
1443 1443
    typename enable_if<typename T::value_type::second_type::LpCol, void>::type
1444 1444
    colBounds(T &t, Value lower, Value upper,dummy<1> = 1) {
1445 1445
      for(typename T::iterator i=t.begin();i!=t.end();++i) {
1446 1446
        colBounds(i->second, lower, upper);
1447 1447
      }
1448 1448
    }
1449 1449
    template<class T>
1450 1450
    typename enable_if<typename T::MapIt::Value::LpCol, void>::type
1451 1451
    colBounds(T &t, Value lower, Value upper,dummy<2> = 2) {
1452 1452
      for(typename T::MapIt i(t); i!=INVALID; ++i){
1453 1453
        colBounds(*i, lower, upper);
1454 1454
      }
1455 1455
    }
1456 1456
#endif
1457 1457

	
1458 1458
    /// Set the lower bound of a row (i.e a constraint)
1459 1459

	
1460 1460
    /// The lower bound of a constraint (row) has to be given by an
1461 1461
    /// extended number of type Value, i.e. a finite number of type
1462 1462
    /// Value or -\ref INF.
1463 1463
    void rowLowerBound(Row r, Value value) {
1464 1464
      _setRowLowerBound(rows(id(r)),value);
1465 1465
    }
1466 1466

	
1467 1467
    /// Get the lower bound of a row (i.e a constraint)
1468 1468

	
1469 1469
    /// This function returns the lower bound for row (constraint) \c c
1470 1470
    /// (this might be -\ref INF as well).
1471 1471
    ///\return The lower bound for row \c r
1472 1472
    Value rowLowerBound(Row r) const {
1473 1473
      return _getRowLowerBound(rows(id(r)));
1474 1474
    }
1475 1475

	
1476 1476
    /// Set the upper bound of a row (i.e a constraint)
1477 1477

	
1478 1478
    /// The upper bound of a constraint (row) has to be given by an
1479 1479
    /// extended number of type Value, i.e. a finite number of type
1480 1480
    /// Value or -\ref INF.
1481 1481
    void rowUpperBound(Row r, Value value) {
1482 1482
      _setRowUpperBound(rows(id(r)),value);
1483 1483
    }
1484 1484

	
1485 1485
    /// Get the upper bound of a row (i.e a constraint)
1486 1486

	
1487 1487
    /// This function returns the upper bound for row (constraint) \c c
1488 1488
    /// (this might be -\ref INF as well).
1489 1489
    ///\return The upper bound for row \c r
1490 1490
    Value rowUpperBound(Row r) const {
1491 1491
      return _getRowUpperBound(rows(id(r)));
1492 1492
    }
1493 1493

	
1494 1494
    ///Set an element of the objective function
1495 1495
    void objCoeff(Col c, Value v) {_setObjCoeff(cols(id(c)),v); };
1496 1496

	
1497 1497
    ///Get an element of the objective function
1498 1498
    Value objCoeff(Col c) const { return _getObjCoeff(cols(id(c))); };
1499 1499

	
1500 1500
    ///Set the objective function
1501 1501

	
1502 1502
    ///\param e is a linear expression of type \ref Expr.
1503 1503
    ///
1504 1504
    void obj(const Expr& e) {
1505 1505
      _setObjCoeffs(ExprIterator(e.comps.begin(), cols),
1506 1506
                    ExprIterator(e.comps.end(), cols));
1507 1507
      obj_const_comp = *e;
1508 1508
    }
1509 1509

	
1510 1510
    ///Get the objective function
1511 1511

	
1512 1512
    ///\return the objective function as a linear expression of type
1513 1513
    ///Expr.
1514 1514
    Expr obj() const {
1515 1515
      Expr e;
1516 1516
      _getObjCoeffs(InsertIterator(e.comps, cols));
1517 1517
      *e = obj_const_comp;
1518 1518
      return e;
1519 1519
    }
1520 1520

	
1521 1521

	
1522 1522
    ///Set the direction of optimization
1523 1523
    void sense(Sense sense) { _setSense(sense); }
1524 1524

	
1525 1525
    ///Query the direction of the optimization
1526 1526
    Sense sense() const {return _getSense(); }
1527 1527

	
1528 1528
    ///Set the sense to maximization
1529 1529
    void max() { _setSense(MAX); }
1530 1530

	
1531 1531
    ///Set the sense to maximization
1532 1532
    void min() { _setSense(MIN); }
1533 1533

	
1534 1534
    ///Clears the problem
1535 1535
    void clear() { _clear(); }
1536 1536

	
1537 1537
    ///@}
1538 1538

	
1539 1539
  };
1540 1540

	
1541 1541
  /// Addition
1542 1542

	
1543 1543
  ///\relates LpBase::Expr
1544 1544
  ///
1545 1545
  inline LpBase::Expr operator+(const LpBase::Expr &a, const LpBase::Expr &b) {
1546 1546
    LpBase::Expr tmp(a);
1547 1547
    tmp+=b;
1548 1548
    return tmp;
1549 1549
  }
1550 1550
  ///Substraction
1551 1551

	
1552 1552
  ///\relates LpBase::Expr
1553 1553
  ///
1554 1554
  inline LpBase::Expr operator-(const LpBase::Expr &a, const LpBase::Expr &b) {
1555 1555
    LpBase::Expr tmp(a);
1556 1556
    tmp-=b;
1557 1557
    return tmp;
1558 1558
  }
1559 1559
  ///Multiply with constant
1560 1560

	
1561 1561
  ///\relates LpBase::Expr
1562 1562
  ///
1563 1563
  inline LpBase::Expr operator*(const LpBase::Expr &a, const LpBase::Value &b) {
1564 1564
    LpBase::Expr tmp(a);
1565 1565
    tmp*=b;
1566 1566
    return tmp;
1567 1567
  }
1568 1568

	
1569 1569
  ///Multiply with constant
1570 1570

	
1571 1571
  ///\relates LpBase::Expr
1572 1572
  ///
1573 1573
  inline LpBase::Expr operator*(const LpBase::Value &a, const LpBase::Expr &b) {
1574 1574
    LpBase::Expr tmp(b);
1575 1575
    tmp*=a;
1576 1576
    return tmp;
1577 1577
  }
1578 1578
  ///Divide with constant
1579 1579

	
1580 1580
  ///\relates LpBase::Expr
1581 1581
  ///
1582 1582
  inline LpBase::Expr operator/(const LpBase::Expr &a, const LpBase::Value &b) {
1583 1583
    LpBase::Expr tmp(a);
1584 1584
    tmp/=b;
1585 1585
    return tmp;
1586 1586
  }
1587 1587

	
1588 1588
  ///Create constraint
1589 1589

	
1590 1590
  ///\relates LpBase::Constr
1591 1591
  ///
1592 1592
  inline LpBase::Constr operator<=(const LpBase::Expr &e,
1593 1593
                                   const LpBase::Expr &f) {
1594 1594
    return LpBase::Constr(0, f - e, LpBase::INF);
1595 1595
  }
1596 1596

	
1597 1597
  ///Create constraint
1598 1598

	
1599 1599
  ///\relates LpBase::Constr
1600 1600
  ///
1601 1601
  inline LpBase::Constr operator<=(const LpBase::Value &e,
1602 1602
                                   const LpBase::Expr &f) {
1603 1603
    return LpBase::Constr(e, f, LpBase::NaN);
1604 1604
  }
1605 1605

	
1606 1606
  ///Create constraint
1607 1607

	
1608 1608
  ///\relates LpBase::Constr
1609 1609
  ///
1610 1610
  inline LpBase::Constr operator<=(const LpBase::Expr &e,
1611 1611
                                   const LpBase::Value &f) {
1612 1612
    return LpBase::Constr(- LpBase::INF, e, f);
1613 1613
  }
1614 1614

	
1615 1615
  ///Create constraint
1616 1616

	
1617 1617
  ///\relates LpBase::Constr
1618 1618
  ///
1619 1619
  inline LpBase::Constr operator>=(const LpBase::Expr &e,
1620 1620
                                   const LpBase::Expr &f) {
1621 1621
    return LpBase::Constr(0, e - f, LpBase::INF);
1622 1622
  }
1623 1623

	
1624 1624

	
1625 1625
  ///Create constraint
1626 1626

	
1627 1627
  ///\relates LpBase::Constr
1628 1628
  ///
1629 1629
  inline LpBase::Constr operator>=(const LpBase::Value &e,
1630 1630
                                   const LpBase::Expr &f) {
1631 1631
    return LpBase::Constr(LpBase::NaN, f, e);
1632 1632
  }
1633 1633

	
1634 1634

	
1635 1635
  ///Create constraint
1636 1636

	
1637 1637
  ///\relates LpBase::Constr
1638 1638
  ///
1639 1639
  inline LpBase::Constr operator>=(const LpBase::Expr &e,
1640 1640
                                   const LpBase::Value &f) {
1641 1641
    return LpBase::Constr(f, e, LpBase::INF);
1642 1642
  }
1643 1643

	
1644 1644
  ///Create constraint
1645 1645

	
1646 1646
  ///\relates LpBase::Constr
1647 1647
  ///
1648 1648
  inline LpBase::Constr operator==(const LpBase::Expr &e,
1649 1649
                                   const LpBase::Value &f) {
1650 1650
    return LpBase::Constr(f, e, f);
1651 1651
  }
1652 1652

	
1653 1653
  ///Create constraint
1654 1654

	
1655 1655
  ///\relates LpBase::Constr
1656 1656
  ///
1657 1657
  inline LpBase::Constr operator==(const LpBase::Expr &e,
1658 1658
                                   const LpBase::Expr &f) {
1659 1659
    return LpBase::Constr(0, f - e, 0);
1660 1660
  }
1661 1661

	
1662 1662
  ///Create constraint
1663 1663

	
1664 1664
  ///\relates LpBase::Constr
1665 1665
  ///
1666 1666
  inline LpBase::Constr operator<=(const LpBase::Value &n,
1667 1667
                                   const LpBase::Constr &c) {
1668 1668
    LpBase::Constr tmp(c);
1669
    LEMON_ASSERT(isnan(tmp.lowerBound()), "Wrong LP constraint");
1669
    LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
1670 1670
    tmp.lowerBound()=n;
1671 1671
    return tmp;
1672 1672
  }
1673 1673
  ///Create constraint
1674 1674

	
1675 1675
  ///\relates LpBase::Constr
1676 1676
  ///
1677 1677
  inline LpBase::Constr operator<=(const LpBase::Constr &c,
1678 1678
                                   const LpBase::Value &n)
1679 1679
  {
1680 1680
    LpBase::Constr tmp(c);
1681
    LEMON_ASSERT(isnan(tmp.upperBound()), "Wrong LP constraint");
1681
    LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
1682 1682
    tmp.upperBound()=n;
1683 1683
    return tmp;
1684 1684
  }
1685 1685

	
1686 1686
  ///Create constraint
1687 1687

	
1688 1688
  ///\relates LpBase::Constr
1689 1689
  ///
1690 1690
  inline LpBase::Constr operator>=(const LpBase::Value &n,
1691 1691
                                   const LpBase::Constr &c) {
1692 1692
    LpBase::Constr tmp(c);
1693
    LEMON_ASSERT(isnan(tmp.upperBound()), "Wrong LP constraint");
1693
    LEMON_ASSERT(isNaN(tmp.upperBound()), "Wrong LP constraint");
1694 1694
    tmp.upperBound()=n;
1695 1695
    return tmp;
1696 1696
  }
1697 1697
  ///Create constraint
1698 1698

	
1699 1699
  ///\relates LpBase::Constr
1700 1700
  ///
1701 1701
  inline LpBase::Constr operator>=(const LpBase::Constr &c,
1702 1702
                                   const LpBase::Value &n)
1703 1703
  {
1704 1704
    LpBase::Constr tmp(c);
1705
    LEMON_ASSERT(isnan(tmp.lowerBound()), "Wrong LP constraint");
1705
    LEMON_ASSERT(isNaN(tmp.lowerBound()), "Wrong LP constraint");
1706 1706
    tmp.lowerBound()=n;
1707 1707
    return tmp;
1708 1708
  }
1709 1709

	
1710 1710
  ///Addition
1711 1711

	
1712 1712
  ///\relates LpBase::DualExpr
1713 1713
  ///
1714 1714
  inline LpBase::DualExpr operator+(const LpBase::DualExpr &a,
1715 1715
                                    const LpBase::DualExpr &b) {
1716 1716
    LpBase::DualExpr tmp(a);
1717 1717
    tmp+=b;
1718 1718
    return tmp;
1719 1719
  }
1720 1720
  ///Substraction
1721 1721

	
1722 1722
  ///\relates LpBase::DualExpr
1723 1723
  ///
1724 1724
  inline LpBase::DualExpr operator-(const LpBase::DualExpr &a,
1725 1725
                                    const LpBase::DualExpr &b) {
1726 1726
    LpBase::DualExpr tmp(a);
1727 1727
    tmp-=b;
1728 1728
    return tmp;
1729 1729
  }
1730 1730
  ///Multiply with constant
1731 1731

	
1732 1732
  ///\relates LpBase::DualExpr
1733 1733
  ///
1734 1734
  inline LpBase::DualExpr operator*(const LpBase::DualExpr &a,
1735 1735
                                    const LpBase::Value &b) {
1736 1736
    LpBase::DualExpr tmp(a);
1737 1737
    tmp*=b;
1738 1738
    return tmp;
1739 1739
  }
1740 1740

	
1741 1741
  ///Multiply with constant
1742 1742

	
1743 1743
  ///\relates LpBase::DualExpr
1744 1744
  ///
1745 1745
  inline LpBase::DualExpr operator*(const LpBase::Value &a,
1746 1746
                                    const LpBase::DualExpr &b) {
1747 1747
    LpBase::DualExpr tmp(b);
1748 1748
    tmp*=a;
1749 1749
    return tmp;
1750 1750
  }
1751 1751
  ///Divide with constant
1752 1752

	
1753 1753
  ///\relates LpBase::DualExpr
1754 1754
  ///
1755 1755
  inline LpBase::DualExpr operator/(const LpBase::DualExpr &a,
1756 1756
                                    const LpBase::Value &b) {
1757 1757
    LpBase::DualExpr tmp(a);
1758 1758
    tmp/=b;
1759 1759
    return tmp;
1760 1760
  }
1761 1761

	
1762 1762
  /// \ingroup lp_group
1763 1763
  ///
1764 1764
  /// \brief Common base class for LP solvers
1765 1765
  ///
1766 1766
  /// This class is an abstract base class for LP solvers. This class
1767 1767
  /// provides a full interface for set and modify an LP problem,
1768 1768
  /// solve it and retrieve the solution. You can use one of the
1769 1769
  /// descendants as a concrete implementation, or the \c Lp
1770 1770
  /// default LP solver. However, if you would like to handle LP
1771 1771
  /// solvers as reference or pointer in a generic way, you can use
1772 1772
  /// this class directly.
1773 1773
  class LpSolver : virtual public LpBase {
1774 1774
  public:
1775 1775

	
1776 1776
    /// The problem types for primal and dual problems
1777 1777
    enum ProblemType {
1778 1778
      ///Feasible solution hasn't been found (but may exist).
1779 1779
      UNDEFINED = 0,
1780 1780
      ///The problem has no feasible solution
1781 1781
      INFEASIBLE = 1,
1782 1782
      ///Feasible solution found
1783 1783
      FEASIBLE = 2,
1784 1784
      ///Optimal solution exists and found
1785 1785
      OPTIMAL = 3,
1786 1786
      ///The cost function is unbounded
1787 1787
      UNBOUNDED = 4
1788 1788
    };
1789 1789

	
1790 1790
    ///The basis status of variables
1791 1791
    enum VarStatus {
1792 1792
      /// The variable is in the basis
1793 1793
      BASIC, 
1794 1794
      /// The variable is free, but not basic
1795 1795
      FREE,
1796 1796
      /// The variable has active lower bound 
1797 1797
      LOWER,
1798 1798
      /// The variable has active upper bound
1799 1799
      UPPER,
1800 1800
      /// The variable is non-basic and fixed
1801 1801
      FIXED
1802 1802
    };
1803 1803

	
1804 1804
  protected:
1805 1805

	
1806 1806
    virtual SolveExitStatus _solve() = 0;
1807 1807

	
1808 1808
    virtual Value _getPrimal(int i) const = 0;
1809 1809
    virtual Value _getDual(int i) const = 0;
1810 1810

	
1811 1811
    virtual Value _getPrimalRay(int i) const = 0;
1812 1812
    virtual Value _getDualRay(int i) const = 0;
1813 1813

	
1814 1814
    virtual Value _getPrimalValue() const = 0;
1815 1815

	
1816 1816
    virtual VarStatus _getColStatus(int i) const = 0;
1817 1817
    virtual VarStatus _getRowStatus(int i) const = 0;
1818 1818

	
1819 1819
    virtual ProblemType _getPrimalType() const = 0;
1820 1820
    virtual ProblemType _getDualType() const = 0;
1821 1821

	
1822 1822
  public:
1823 1823

	
1824 1824
    ///\name Solve the LP
1825 1825

	
1826 1826
    ///@{
1827 1827

	
1828 1828
    ///\e Solve the LP problem at hand
1829 1829
    ///
1830 1830
    ///\return The result of the optimization procedure. Possible
1831 1831
    ///values and their meanings can be found in the documentation of
1832 1832
    ///\ref SolveExitStatus.
1833 1833
    SolveExitStatus solve() { return _solve(); }
1834 1834

	
1835 1835
    ///@}
1836 1836

	
1837 1837
    ///\name Obtain the solution
1838 1838

	
1839 1839
    ///@{
1840 1840

	
1841 1841
    /// The type of the primal problem
1842 1842
    ProblemType primalType() const {
1843 1843
      return _getPrimalType();
1844 1844
    }
1845 1845

	
1846 1846
    /// The type of the dual problem
1847 1847
    ProblemType dualType() const {
1848 1848
      return _getDualType();
1849 1849
    }
1850 1850

	
1851 1851
    /// Return the primal value of the column
1852 1852

	
1853 1853
    /// Return the primal value of the column.
1854 1854
    /// \pre The problem is solved.
1855 1855
    Value primal(Col c) const { return _getPrimal(cols(id(c))); }
1856 1856

	
1857 1857
    /// Return the primal value of the expression
1858 1858

	
1859 1859
    /// Return the primal value of the expression, i.e. the dot
1860 1860
    /// product of the primal solution and the expression.
1861 1861
    /// \pre The problem is solved.
1862 1862
    Value primal(const Expr& e) const {
1863 1863
      double res = *e;
1864 1864
      for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
1865 1865
        res += *c * primal(c);
1866 1866
      }
1867 1867
      return res;
1868 1868
    }
1869 1869
    /// Returns a component of the primal ray
1870 1870
    
1871 1871
    /// The primal ray is solution of the modified primal problem,
1872 1872
    /// where we change each finite bound to 0, and we looking for a
1873 1873
    /// negative objective value in case of minimization, and positive
1874 1874
    /// objective value for maximization. If there is such solution,
1875 1875
    /// that proofs the unsolvability of the dual problem, and if a
1876 1876
    /// feasible primal solution exists, then the unboundness of
1877 1877
    /// primal problem.
1878 1878
    ///
1879 1879
    /// \pre The problem is solved and the dual problem is infeasible.
1880 1880
    /// \note Some solvers does not provide primal ray calculation
1881 1881
    /// functions.
1882 1882
    Value primalRay(Col c) const { return _getPrimalRay(cols(id(c))); }
1883 1883

	
1884 1884
    /// Return the dual value of the row
1885 1885

	
1886 1886
    /// Return the dual value of the row.
1887 1887
    /// \pre The problem is solved.
1888 1888
    Value dual(Row r) const { return _getDual(rows(id(r))); }
1889 1889

	
1890 1890
    /// Return the dual value of the dual expression
1891 1891

	
1892 1892
    /// Return the dual value of the dual expression, i.e. the dot
1893 1893
    /// product of the dual solution and the dual expression.
1894 1894
    /// \pre The problem is solved.
1895 1895
    Value dual(const DualExpr& e) const {
1896 1896
      double res = 0.0;
1897 1897
      for (DualExpr::ConstCoeffIt r(e); r != INVALID; ++r) {
1898 1898
        res += *r * dual(r);
1899 1899
      }
1900 1900
      return res;
1901 1901
    }
1902 1902

	
1903 1903
    /// Returns a component of the dual ray
1904 1904
    
1905 1905
    /// The dual ray is solution of the modified primal problem, where
1906 1906
    /// we change each finite bound to 0 (i.e. the objective function
1907 1907
    /// coefficients in the primal problem), and we looking for a
1908 1908
    /// ositive objective value. If there is such solution, that
1909 1909
    /// proofs the unsolvability of the primal problem, and if a
1910 1910
    /// feasible dual solution exists, then the unboundness of
1911 1911
    /// dual problem.
1912 1912
    ///
1913 1913
    /// \pre The problem is solved and the primal problem is infeasible.
1914 1914
    /// \note Some solvers does not provide dual ray calculation
1915 1915
    /// functions.
1916 1916
    Value dualRay(Row r) const { return _getDualRay(rows(id(r))); }
1917 1917

	
1918 1918
    /// Return the basis status of the column
1919 1919

	
1920 1920
    /// \see VarStatus
1921 1921
    VarStatus colStatus(Col c) const { return _getColStatus(cols(id(c))); }
1922 1922

	
1923 1923
    /// Return the basis status of the row
1924 1924

	
1925 1925
    /// \see VarStatus
1926 1926
    VarStatus rowStatus(Row r) const { return _getRowStatus(rows(id(r))); }
1927 1927

	
1928 1928
    ///The value of the objective function
1929 1929

	
1930 1930
    ///\return
1931 1931
    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
1932 1932
    /// of the primal problem, depending on whether we minimize or maximize.
1933 1933
    ///- \ref NaN if no primal solution is found.
1934 1934
    ///- The (finite) objective value if an optimal solution is found.
1935 1935
    Value primal() const { return _getPrimalValue()+obj_const_comp;}
1936 1936
    ///@}
1937 1937

	
1938 1938
    LpSolver* newSolver() {return _newSolver();}
1939 1939
    LpSolver* cloneSolver() {return _cloneSolver();}
1940 1940

	
1941 1941
  protected:
1942 1942

	
1943 1943
    virtual LpSolver* _newSolver() const = 0;
1944 1944
    virtual LpSolver* _cloneSolver() const = 0;
1945 1945
  };
1946 1946

	
1947 1947

	
1948 1948
  /// \ingroup lp_group
1949 1949
  ///
1950 1950
  /// \brief Common base class for MIP solvers
1951 1951
  ///
1952 1952
  /// This class is an abstract base class for MIP solvers. This class
1953 1953
  /// provides a full interface for set and modify an MIP problem,
1954 1954
  /// solve it and retrieve the solution. You can use one of the
1955 1955
  /// descendants as a concrete implementation, or the \c Lp
1956 1956
  /// default MIP solver. However, if you would like to handle MIP
1957 1957
  /// solvers as reference or pointer in a generic way, you can use
1958 1958
  /// this class directly.
1959 1959
  class MipSolver : virtual public LpBase {
1960 1960
  public:
1961 1961

	
1962 1962
    /// The problem types for MIP problems
1963 1963
    enum ProblemType {
1964 1964
      ///Feasible solution hasn't been found (but may exist).
1965 1965
      UNDEFINED = 0,
1966 1966
      ///The problem has no feasible solution
1967 1967
      INFEASIBLE = 1,
1968 1968
      ///Feasible solution found
1969 1969
      FEASIBLE = 2,
1970 1970
      ///Optimal solution exists and found
1971 1971
      OPTIMAL = 3,
1972 1972
      ///The cost function is unbounded
1973 1973
      ///
1974 1974
      ///The Mip or at least the relaxed problem is unbounded
1975 1975
      UNBOUNDED = 4
1976 1976
    };
1977 1977

	
1978 1978
    ///\name Solve the MIP
1979 1979

	
1980 1980
    ///@{
1981 1981

	
1982 1982
    /// Solve the MIP problem at hand
1983 1983
    ///
1984 1984
    ///\return The result of the optimization procedure. Possible
1985 1985
    ///values and their meanings can be found in the documentation of
1986 1986
    ///\ref SolveExitStatus.
1987 1987
    SolveExitStatus solve() { return _solve(); }
1988 1988

	
1989 1989
    ///@}
1990 1990

	
1991 1991
    ///\name Setting column type
1992 1992
    ///@{
1993 1993

	
1994 1994
    ///Possible variable (column) types (e.g. real, integer, binary etc.)
1995 1995
    enum ColTypes {
1996 1996
      ///Continuous variable (default)
1997 1997
      REAL = 0,
1998 1998
      ///Integer variable
1999 1999
      INTEGER = 1
2000 2000
    };
2001 2001

	
2002 2002
    ///Sets the type of the given column to the given type
2003 2003

	
2004 2004
    ///Sets the type of the given column to the given type.
2005 2005
    ///
2006 2006
    void colType(Col c, ColTypes col_type) {
2007 2007
      _setColType(cols(id(c)),col_type);
2008 2008
    }
2009 2009

	
2010 2010
    ///Gives back the type of the column.
2011 2011

	
2012 2012
    ///Gives back the type of the column.
2013 2013
    ///
2014 2014
    ColTypes colType(Col c) const {
2015 2015
      return _getColType(cols(id(c)));
2016 2016
    }
2017 2017
    ///@}
2018 2018

	
2019 2019
    ///\name Obtain the solution
2020 2020

	
2021 2021
    ///@{
2022 2022

	
2023 2023
    /// The type of the MIP problem
2024 2024
    ProblemType type() const {
2025 2025
      return _getType();
2026 2026
    }
2027 2027

	
2028 2028
    /// Return the value of the row in the solution
2029 2029

	
2030 2030
    ///  Return the value of the row in the solution.
2031 2031
    /// \pre The problem is solved.
2032 2032
    Value sol(Col c) const { return _getSol(cols(id(c))); }
2033 2033

	
2034 2034
    /// Return the value of the expression in the solution
2035 2035

	
2036 2036
    /// Return the value of the expression in the solution, i.e. the
2037 2037
    /// dot product of the solution and the expression.
2038 2038
    /// \pre The problem is solved.
2039 2039
    Value sol(const Expr& e) const {
2040 2040
      double res = *e;
2041 2041
      for (Expr::ConstCoeffIt c(e); c != INVALID; ++c) {
2042 2042
        res += *c * sol(c);
2043 2043
      }
2044 2044
      return res;
2045 2045
    }
2046 2046
    ///The value of the objective function
2047 2047
    
2048 2048
    ///\return
2049 2049
    ///- \ref INF or -\ref INF means either infeasibility or unboundedness
2050 2050
    /// of the problem, depending on whether we minimize or maximize.
2051 2051
    ///- \ref NaN if no primal solution is found.
2052 2052
    ///- The (finite) objective value if an optimal solution is found.
2053 2053
    Value solValue() const { return _getSolValue()+obj_const_comp;}
2054 2054
    ///@}
2055 2055

	
2056 2056
  protected:
2057 2057

	
2058 2058
    virtual SolveExitStatus _solve() = 0;
2059 2059
    virtual ColTypes _getColType(int col) const = 0;
2060 2060
    virtual void _setColType(int col, ColTypes col_type) = 0;
2061 2061
    virtual ProblemType _getType() const = 0;
2062 2062
    virtual Value _getSol(int i) const = 0;
2063 2063
    virtual Value _getSolValue() const = 0;
2064 2064

	
2065 2065
  public:
2066 2066

	
2067 2067
    MipSolver* newSolver() {return _newSolver();}
2068 2068
    MipSolver* cloneSolver() {return _cloneSolver();}
2069 2069

	
2070 2070
  protected:
2071 2071

	
2072 2072
    virtual MipSolver* _newSolver() const = 0;
2073 2073
    virtual MipSolver* _cloneSolver() const = 0;
2074 2074
  };
2075 2075

	
2076 2076

	
2077 2077

	
2078 2078
} //namespace lemon
2079 2079

	
2080 2080
#endif //LEMON_LP_BASE_H
Ignore white space 16777216 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_MATH_H
20 20
#define LEMON_MATH_H
21 21

	
22 22
///\ingroup misc
23 23
///\file
24 24
///\brief Some extensions to the standard \c cmath library.
25 25
///
26 26
///Some extensions to the standard \c cmath library.
27 27
///
28 28
///This file includes the standard math library (cmath).
29 29

	
30 30
#include<cmath>
31 31

	
32 32
namespace lemon {
33 33

	
34 34
  /// \addtogroup misc
35 35
  /// @{
36 36

	
37 37
  /// The Euler constant
38 38
  const long double E       = 2.7182818284590452353602874713526625L;
39 39
  /// log_2(e)
40 40
  const long double LOG2E   = 1.4426950408889634073599246810018921L;
41 41
  /// log_10(e)
42 42
  const long double LOG10E  = 0.4342944819032518276511289189166051L;
43 43
  /// ln(2)
44 44
  const long double LN2     = 0.6931471805599453094172321214581766L;
45 45
  /// ln(10)
46 46
  const long double LN10    = 2.3025850929940456840179914546843642L;
47 47
  /// pi
48 48
  const long double PI      = 3.1415926535897932384626433832795029L;
49 49
  /// pi/2
50 50
  const long double PI_2    = 1.5707963267948966192313216916397514L;
51 51
  /// pi/4
52 52
  const long double PI_4    = 0.7853981633974483096156608458198757L;
53 53
  /// sqrt(2)
54 54
  const long double SQRT2   = 1.4142135623730950488016887242096981L;
55 55
  /// 1/sqrt(2)
56 56
  const long double SQRT1_2 = 0.7071067811865475244008443621048490L;
57 57

	
58 58
  ///Check whether the parameter is NaN or not
59 59
  
60 60
  ///This function checks whether the parameter is NaN or not.
61 61
  ///Is should be equivalent with std::isnan(), but it is not
62 62
  ///provided by all compilers.
63
  inline bool isnan(double v)
63
  inline bool isNaN(double v)
64 64
    {
65 65
      return v!=v;
66 66
    }
67 67

	
68 68
  /// @}
69 69

	
70 70
} //namespace lemon
71 71

	
72 72
#endif //LEMON_TOLERANCE_H
Ignore white space 16777216 line context
1 1
/* -*- mode: C++; indent-tabs-mode: nil; -*-
2 2
 *
3 3
 * This file is a part of LEMON, a generic C++ optimization library.
4 4
 *
5 5
 * Copyright (C) 2003-2009
6 6
 * Egervary Jeno Kombinatorikus Optimalizalasi Kutatocsoport
7 7
 * (Egervary Research Group on Combinatorial Optimization, EGRES).
8 8
 *
9 9
 * Permission to use, modify and distribute this software is granted
10 10
 * provided that this copyright notice appears in all copies. For
11 11
 * precise terms see the accompanying LICENSE file.
12 12
 *
13 13
 * This software is provided "AS IS" with no warranty of any kind,
14 14
 * express or implied, and with no claim as to its suitability for any
15 15
 * purpose.
16 16
 *
17 17
 */
18 18

	
19 19
#ifndef LEMON_TIME_MEASURE_H
20 20
#define LEMON_TIME_MEASURE_H
21 21

	
22 22
///\ingroup timecount
23 23
///\file
24 24
///\brief Tools for measuring cpu usage
25 25

	
26 26
#ifdef WIN32
27
#ifndef WIN32_LEAN_AND_MEAN
27 28
#define WIN32_LEAN_AND_MEAN
29
#endif
30
#ifndef NOMINMAX
28 31
#define NOMINMAX
32
#endif
29 33
#include <windows.h>
30 34
#include <cmath>
31 35
#else
36
#include <unistd.h>
32 37
#include <sys/times.h>
33 38
#include <sys/time.h>
34 39
#endif
35 40

	
36 41
#include <string>
37 42
#include <fstream>
38 43
#include <iostream>
39 44

	
40 45
namespace lemon {
41 46

	
42 47
  /// \addtogroup timecount
43 48
  /// @{
44 49

	
45 50
  /// A class to store (cpu)time instances.
46 51

	
47 52
  /// This class stores five time values.
48 53
  /// - a real time
49 54
  /// - a user cpu time
50 55
  /// - a system cpu time
51 56
  /// - a user cpu time of children
52 57
  /// - a system cpu time of children
53 58
  ///
54 59
  /// TimeStamp's can be added to or substracted from each other and
55 60
  /// they can be pushed to a stream.
56 61
  ///
57 62
  /// In most cases, perhaps the \ref Timer or the \ref TimeReport
58 63
  /// class is what you want to use instead.
59 64

	
60 65
  class TimeStamp
61 66
  {
62 67
    double utime;
63 68
    double stime;
64 69
    double cutime;
65 70
    double cstime;
66 71
    double rtime;
67 72

	
68 73
    void _reset() {
69 74
      utime = stime = cutime = cstime = rtime = 0;
70 75
    }
71 76

	
72 77
  public:
73 78

	
74 79
    ///Read the current time values of the process
75 80
    void stamp()
76 81
    {
77 82
#ifndef WIN32
78 83
      timeval tv;
79 84
      gettimeofday(&tv, 0);
80 85
      rtime=tv.tv_sec+double(tv.tv_usec)/1e6;
81 86

	
82 87
      tms ts;
83 88
      double tck=sysconf(_SC_CLK_TCK);
84 89
      times(&ts);
85 90
      utime=ts.tms_utime/tck;
86 91
      stime=ts.tms_stime/tck;
87 92
      cutime=ts.tms_cutime/tck;
88 93
      cstime=ts.tms_cstime/tck;
89 94
#else
90 95
      static const double ch = 4294967296.0e-7;
91 96
      static const double cl = 1.0e-7;
92 97

	
93 98
      FILETIME system;
94 99
      GetSystemTimeAsFileTime(&system);
95 100
      rtime = ch * system.dwHighDateTime + cl * system.dwLowDateTime;
96 101

	
97 102
      FILETIME create, exit, kernel, user;
98 103
      if (GetProcessTimes(GetCurrentProcess(),&create, &exit, &kernel, &user)) {
99 104
        utime = ch * user.dwHighDateTime + cl * user.dwLowDateTime;
100 105
        stime = ch * kernel.dwHighDateTime + cl * kernel.dwLowDateTime;
101 106
        cutime = 0;
102 107
        cstime = 0;
103 108
      } else {
104 109
        rtime = 0;
105 110
        utime = 0;
106 111
        stime = 0;
107 112
        cutime = 0;
108 113
        cstime = 0;
109 114
      }
110 115
#endif
111 116
    }
112 117

	
113 118
    /// Constructor initializing with zero
114 119
    TimeStamp()
115 120
    { _reset(); }
116 121
    ///Constructor initializing with the current time values of the process
117 122
    TimeStamp(void *) { stamp();}
118 123

	
119 124
    ///Set every time value to zero
120 125
    TimeStamp &reset() {_reset();return *this;}
121 126

	
122 127
    ///\e
123 128
    TimeStamp &operator+=(const TimeStamp &b)
124 129
    {
125 130
      utime+=b.utime;
126 131
      stime+=b.stime;
127 132
      cutime+=b.cutime;
128 133
      cstime+=b.cstime;
129 134
      rtime+=b.rtime;
130 135
      return *this;
131 136
    }
132 137
    ///\e
133 138
    TimeStamp operator+(const TimeStamp &b) const
134 139
    {
135 140
      TimeStamp t(*this);
136 141
      return t+=b;
137 142
    }
138 143
    ///\e
139 144
    TimeStamp &operator-=(const TimeStamp &b)
140 145
    {
141 146
      utime-=b.utime;
142 147
      stime-=b.stime;
143 148
      cutime-=b.cutime;
144 149
      cstime-=b.cstime;
145 150
      rtime-=b.rtime;
146 151
      return *this;
147 152
    }
148 153
    ///\e
149 154
    TimeStamp operator-(const TimeStamp &b) const
150 155
    {
151 156
      TimeStamp t(*this);
152 157
      return t-=b;
153 158
    }
154 159
    ///\e
155 160
    TimeStamp &operator*=(double b)
156 161
    {
157 162
      utime*=b;
158 163
      stime*=b;
159 164
      cutime*=b;
160 165
      cstime*=b;
161 166
      rtime*=b;
162 167
      return *this;
163 168
    }
164 169
    ///\e
165 170
    TimeStamp operator*(double b) const
166 171
    {
167 172
      TimeStamp t(*this);
168 173
      return t*=b;
169 174
    }
170 175
    friend TimeStamp operator*(double b,const TimeStamp &t);
171 176
    ///\e
172 177
    TimeStamp &operator/=(double b)
173 178
    {
174 179
      utime/=b;
175 180
      stime/=b;
176 181
      cutime/=b;
177 182
      cstime/=b;
178 183
      rtime/=b;
179 184
      return *this;
180 185
    }
181 186
    ///\e
182 187
    TimeStamp operator/(double b) const
183 188
    {
184 189
      TimeStamp t(*this);
185 190
      return t/=b;
186 191
    }
187 192
    ///The time ellapsed since the last call of stamp()
188 193
    TimeStamp ellapsed() const
189 194
    {
190 195
      TimeStamp t(NULL);
191 196
      return t-*this;
192 197
    }
193 198

	
194 199
    friend std::ostream& operator<<(std::ostream& os,const TimeStamp &t);
195 200

	
196 201
    ///Gives back the user time of the process
197 202
    double userTime() const
198 203
    {
199 204
      return utime;
200 205
    }
201 206
    ///Gives back the system time of the process
202 207
    double systemTime() const
203 208
    {
204 209
      return stime;
205 210
    }
206 211
    ///Gives back the user time of the process' children
207 212

	
208 213
    ///\note On <tt>WIN32</tt> platform this value is not calculated.
209 214
    ///
210 215
    double cUserTime() const
211 216
    {
212 217
      return cutime;
213 218
    }
214 219
    ///Gives back the user time of the process' children
215 220

	
216 221
    ///\note On <tt>WIN32</tt> platform this value is not calculated.
217 222
    ///
218 223
    double cSystemTime() const
219 224
    {
220 225
      return cstime;
221 226
    }
222 227
    ///Gives back the real time
223 228
    double realTime() const {return rtime;}
224 229
  };
225 230

	
226 231
  TimeStamp operator*(double b,const TimeStamp &t)
227 232
  {
228 233
    return t*b;
229 234
  }
230 235

	
231 236
  ///Prints the time counters
232 237

	
233 238
  ///Prints the time counters in the following form:
234 239
  ///
235 240
  /// <tt>u: XX.XXs s: XX.XXs cu: XX.XXs cs: XX.XXs real: XX.XXs</tt>
236 241
  ///
237 242
  /// where the values are the
238 243
  /// \li \c u: user cpu time,
239 244
  /// \li \c s: system cpu time,
240 245
  /// \li \c cu: user cpu time of children,
241 246
  /// \li \c cs: system cpu time of children,
242 247
  /// \li \c real: real time.
243 248
  /// \relates TimeStamp
244 249
  /// \note On <tt>WIN32</tt> platform the cummulative values are not
245 250
  /// calculated.
246 251
  inline std::ostream& operator<<(std::ostream& os,const TimeStamp &t)
247 252
  {
248 253
    os << "u: " << t.userTime() <<
249 254
      "s, s: " << t.systemTime() <<
250 255
      "s, cu: " << t.cUserTime() <<
251 256
      "s, cs: " << t.cSystemTime() <<
252 257
      "s, real: " << t.realTime() << "s";
253 258
    return os;
254 259
  }
255 260

	
256 261
  ///Class for measuring the cpu time and real time usage of the process
257 262

	
258 263
  ///Class for measuring the cpu time and real time usage of the process.
259 264
  ///It is quite easy-to-use, here is a short example.
260 265
  ///\code
261 266
  /// #include<lemon/time_measure.h>
262 267
  /// #include<iostream>
263 268
  ///
264 269
  /// int main()
265 270
  /// {
266 271
  ///
267 272
  ///   ...
268 273
  ///
269 274
  ///   Timer t;
270 275
  ///   doSomething();
271 276
  ///   std::cout << t << '\n';
272 277
  ///   t.restart();
273 278
  ///   doSomethingElse();
274 279
  ///   std::cout << t << '\n';
275 280
  ///
276 281
  ///   ...
277 282
  ///
278 283
  /// }
279 284
  ///\endcode
280 285
  ///
281 286
  ///The \ref Timer can also be \ref stop() "stopped" and
282 287
  ///\ref start() "started" again, so it is possible to compute collected
283 288
  ///running times.
284 289
  ///
285 290
  ///\warning Depending on the operation system and its actual configuration
286 291
  ///the time counters have a certain (10ms on a typical Linux system)
287 292
  ///granularity.
288 293
  ///Therefore this tool is not appropriate to measure very short times.
289 294
  ///Also, if you start and stop the timer very frequently, it could lead to
290 295
  ///distorted results.
291 296
  ///
292 297
  ///\note If you want to measure the running time of the execution of a certain
293 298
  ///function, consider the usage of \ref TimeReport instead.
294 299
  ///
295 300
  ///\sa TimeReport
296 301
  class Timer
297 302
  {
298 303
    int _running; //Timer is running iff _running>0; (_running>=0 always holds)
299 304
    TimeStamp start_time; //This is the relativ start-time if the timer
300 305
                          //is _running, the collected _running time otherwise.
301 306

	
302 307
    void _reset() {if(_running) start_time.stamp(); else start_time.reset();}
303 308

	
304 309
  public:
305 310
    ///Constructor.
306 311

	
307 312
    ///\param run indicates whether or not the timer starts immediately.
308 313
    ///
309 314
    Timer(bool run=true) :_running(run) {_reset();}
310 315

	
311 316
    ///\name Control the state of the timer
312 317
    ///Basically a Timer can be either running or stopped,
313 318
    ///but it provides a bit finer control on the execution.
314 319
    ///The \ref lemon::Timer "Timer" also counts the number of
315 320
    ///\ref lemon::Timer::start() "start()" executions, and it stops
316 321
    ///only after the same amount (or more) \ref lemon::Timer::stop()
317 322
    ///"stop()"s. This can be useful e.g. to compute the running time
318 323
    ///of recursive functions.
319 324

	
320 325
    ///@{
321 326

	
322 327
    ///Reset and stop the time counters
323 328

	
324 329
    ///This function resets and stops the time counters
325 330
    ///\sa restart()
326 331
    void reset()
327 332
    {
328 333
      _running=0;
329 334
      _reset();
330 335
    }
331 336

	
332 337
    ///Start the time counters
333 338

	
334 339
    ///This function starts the time counters.
335 340
    ///
336 341
    ///If the timer is started more than ones, it will remain running
337 342
    ///until the same amount of \ref stop() is called.
338 343
    ///\sa stop()
339 344
    void start()
340 345
    {
341 346
      if(_running) _running++;
342 347
      else {
343 348
        _running=1;
344 349
        TimeStamp t;
345 350
        t.stamp();
346 351
        start_time=t-start_time;
347 352
      }
348 353
    }
349 354

	
350 355

	
351 356
    ///Stop the time counters
352 357

	
353 358
    ///This function stops the time counters. If start() was executed more than
354 359
    ///once, then the same number of stop() execution is necessary the really
355 360
    ///stop the timer.
356 361
    ///
357 362
    ///\sa halt()
358 363
    ///\sa start()
359 364
    ///\sa restart()
360 365
    ///\sa reset()
361 366

	
362 367
    void stop()
363 368
    {
364 369
      if(_running && !--_running) {
365 370
        TimeStamp t;
366 371
        t.stamp();
367 372
        start_time=t-start_time;
368 373
      }
369 374
    }
370 375

	
371 376
    ///Halt (i.e stop immediately) the time counters
372 377

	
373 378
    ///This function stops immediately the time counters, i.e. <tt>t.halt()</tt>
374 379
    ///is a faster
375 380
    ///equivalent of the following.
376 381
    ///\code
377 382
    ///  while(t.running()) t.stop()
378 383
    ///\endcode
379 384
    ///
380 385
    ///
381 386
    ///\sa stop()
382 387
    ///\sa restart()
383 388
    ///\sa reset()
384 389

	
385 390
    void halt()
386 391
    {
387 392
      if(_running) {
388 393
        _running=0;
389 394
        TimeStamp t;
390 395
        t.stamp();
391 396
        start_time=t-start_time;
392 397
      }
393 398
    }
394 399

	
395 400
    ///Returns the running state of the timer
396 401

	
397 402
    ///This function returns the number of stop() exections that is
398 403
    ///necessary to really stop the timer.
399 404
    ///For example the timer
400 405
    ///is running if and only if the return value is \c true
401 406
    ///(i.e. greater than
402 407
    ///zero).
403 408
    int running()  { return _running; }
404 409

	
405 410

	
406 411
    ///Restart the time counters
407 412

	
408 413
    ///This function is a shorthand for
409 414
    ///a reset() and a start() calls.
410 415
    ///
411 416
    void restart()
412 417
    {
413 418
      reset();
414 419
      start();
415 420
    }
416 421

	
417 422
    ///@}
418 423

	
419 424
    ///\name Query Functions for the ellapsed time
420 425

	
421 426
    ///@{
422 427

	
423 428
    ///Gives back the ellapsed user time of the process
424 429
    double userTime() const
425 430
    {
426 431
      return operator TimeStamp().userTime();
427 432
    }
428 433
    ///Gives back the ellapsed system time of the process
429 434
    double systemTime() const
430 435
    {
431 436
      return operator TimeStamp().systemTime();
432 437
    }
433 438
    ///Gives back the ellapsed user time of the process' children
434 439

	
435 440
    ///\note On <tt>WIN32</tt> platform this value is not calculated.
436 441
    ///
437 442
    double cUserTime() const
438 443
    {
439 444
      return operator TimeStamp().cUserTime();
440 445
    }
441 446
    ///Gives back the ellapsed user time of the process' children
442 447

	
443 448
    ///\note On <tt>WIN32</tt> platform this value is not calculated.
444 449
    ///
445 450
    double cSystemTime() const
446 451
    {
447 452
      return operator TimeStamp().cSystemTime();
448 453
    }
449 454
    ///Gives back the ellapsed real time
450 455
    double realTime() const
451 456
    {
452 457
      return operator TimeStamp().realTime();
453 458
    }
454 459
    ///Computes the ellapsed time
455 460

	
456 461
    ///This conversion computes the ellapsed time, therefore you can print
457 462
    ///the ellapsed time like this.
458 463
    ///\code
459 464
    ///  Timer t;
460 465
    ///  doSomething();
461 466
    ///  std::cout << t << '\n';
462 467
    ///\endcode
463 468
    operator TimeStamp () const
464 469
    {
465 470
      TimeStamp t;
466 471
      t.stamp();
467 472
      return _running?t-start_time:start_time;
468 473
    }
469 474

	
470 475

	
471 476
    ///@}
472 477
  };
473 478

	
474 479
  ///Same as Timer but prints a report on destruction.
475 480

	
476 481
  ///Same as \ref Timer but prints a report on destruction.
477 482
  ///This example shows its usage.
478 483
  ///\code
479 484
  ///  void myAlg(ListGraph &g,int n)
480 485
  ///  {
481 486
  ///    TimeReport tr("Running time of myAlg: ");
482 487
  ///    ... //Here comes the algorithm
483 488
  ///  }
484 489
  ///\endcode
485 490
  ///
486 491
  ///\sa Timer
487 492
  ///\sa NoTimeReport
488 493
  class TimeReport : public Timer
489 494
  {
490 495
    std::string _title;
491 496
    std::ostream &_os;
492 497
  public:
493 498
    ///Constructor
494 499

	
495 500
    ///Constructor.
496 501
    ///\param title This text will be printed before the ellapsed time.
497 502
    ///\param os The stream to print the report to.
498 503
    ///\param run Sets whether the timer should start immediately.
499 504
    TimeReport(std::string title,std::ostream &os=std::cerr,bool run=true)
500 505
      : Timer(run), _title(title), _os(os){}
501 506
    ///Destructor that prints the ellapsed time
502 507
    ~TimeReport()
503 508
    {
504 509
      _os << _title << *this << std::endl;
505 510
    }
506 511
  };
507 512

	
508 513
  ///'Do nothing' version of TimeReport
509 514

	
510 515
  ///\sa TimeReport
511 516
  ///
512 517
  class NoTimeReport
513 518
  {
514 519
  public:
515 520
    ///\e
516 521
    NoTimeReport(std::string,std::ostream &,bool) {}
517 522
    ///\e
518 523
    NoTimeReport(std::string,std::ostream &) {}
519 524
    ///\e
520 525
    NoTimeReport(std::string) {}
521 526
    ///\e Do nothing.
522 527
    ~NoTimeReport() {}
523 528

	
524 529
    operator TimeStamp () const { return TimeStamp(); }
525 530
    void reset() {}
526 531
    void start() {}
527 532
    void stop() {}
528 533
    void halt() {}
529 534
    int running() { return 0; }
530 535
    void restart() {}
531 536
    double userTime() const { return 0; }
532 537
    double systemTime() const { return 0; }
533 538
    double cUserTime() const { return 0; }
534 539
    double cSystemTime() const { return 0; }
535 540
    double realTime() const { return 0; }
536 541
  };
537 542

	
538 543
  ///Tool to measure the running time more exactly.
539 544

	
540 545
  ///This function calls \c f several times and returns the average
541 546
  ///running time. The number of the executions will be choosen in such a way
542 547
  ///that the full real running time will be roughly between \c min_time
543 548
  ///and <tt>2*min_time</tt>.
544 549
  ///\param f the function object to be measured.
545 550
  ///\param min_time the minimum total running time.
546 551
  ///\retval num if it is not \c NULL, then the actual
547 552
  ///        number of execution of \c f will be written into <tt>*num</tt>.
548 553
  ///\retval full_time if it is not \c NULL, then the actual
549 554
  ///        total running time will be written into <tt>*full_time</tt>.
550 555
  ///\return The average running time of \c f.
551 556

	
552 557
  template<class F>
553 558
  TimeStamp runningTimeTest(F f,double min_time=10,unsigned int *num = NULL,
554 559
                            TimeStamp *full_time=NULL)
555 560
  {
556 561
    TimeStamp full;
557 562
    unsigned int total=0;
558 563
    Timer t;
559 564
    for(unsigned int tn=1;tn <= 1U<<31 && full.realTime()<=min_time; tn*=2) {
560 565
      for(;total<tn;total++) f();
561 566
      full=t;
562 567
    }
563 568
    if(num) *num=total;
564 569
    if(full_time) *full_time=full;
565 570
    return full/total;
566 571
  }
567 572

	
568 573
  /// @}
569 574

	
570 575

	
571 576
} //namespace lemon
572 577

	
573 578
#endif //LEMON_TIME_MEASURE_H
0 comments (0 inline)