aGrUM  0.20.3
a C++ library for (probabilistic) graphical models
LpInterface_tpl.h
Go to the documentation of this file.
1 /**
2  *
3  * Copyright (c) 2005-2021 by Pierre-Henri WUILLEMIN(@LIP6) & Christophe GONZALES(@AMU)
4  * info_at_agrum_dot_org
5  *
6  * This library is free software: you can redistribute it and/or modify
7  * it under the terms of the GNU Lesser General Public License as published by
8  * the Free Software Foundation, either version 3 of the License, or
9  * (at your option) any later version.
10  *
11  * This library is distributed in the hope that it will be useful,
12  * but WITHOUT ANY WARRANTY; without even the implied warranty of
13  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14  * GNU Lesser General Public License for more details.
15  *
16  * You should have received a copy of the GNU Lesser General Public License
17  * along with this library. If not, see <http://www.gnu.org/licenses/>.
18  *
19  */
20 
21 
22 #include "LpInterface.h"
23 
24 namespace gum {
25  namespace credal {
26  namespace lp {
27 
28  /**
29  * class LpExpr
30  */
31  template < typename SCALAR >
32  LpExpr& LpExpr::operator=(const SCALAR& rhs) {
33  clear();
34 
35  _mValue_ = rhs;
36  _imiddle_ = true;
37 
38  return *this;
39  }
40 
41  template < typename T >
42  LpExpr& LpExpr::operator+=(const T& rhs) {
43  if (_ileft_ || _iright_)
44  GUM_ERROR(OperationNotAllowed, "expr::operator+= (expr) : <= present on one side of expr")
45 
46  if (!_imiddle_) _imiddle_ = true;
47 
48  _mValue_ += rhs;
49 
50  return *this;
51  }
52 
53  template < typename T >
54  LpExpr& LpExpr::operator-=(const T& rhs) {
55  if (_ileft_ || _iright_)
56  GUM_ERROR(OperationNotAllowed, "expr::operator-= (rhs) : <= present in one of expr")
57 
58  if (!_imiddle_) _imiddle_ = true;
59 
60  _mValue_ -= rhs;
61 
62  return *this;
63  }
64 
65  template < typename SCALAR >
66  void LpExpr::_addSide_(const SCALAR& from) {
67  if (!_ileft_) {
68  _lValue_ = from;
69  _ileft_ = true;
70  } else if (!_imiddle_) {
71  _mValue_ = from;
72  _imiddle_ = true;
73  } else if (!_iright_) {
74  _rValue_ = from;
75  _iright_ = true;
76  } else
77  GUM_ERROR(OperationNotAllowed,
78  "LpExpr::setSide ( const LpCol & from "
79  ") : too many <= ; no free side");
80  }
81 
82  /**
83  * class LpInterface
84  */
85  template < typename GUM_SCALAR >
86  LpInterface< GUM_SCALAR >::LpInterface() {
87  _positivity_ = false;
88  _sumIsOne_ = false;
89  GUM_CONSTRUCTOR(LpInterface);
90  }
91 
92  template < typename GUM_SCALAR >
93  LpInterface< GUM_SCALAR >::LpInterface(const LpInterface< GUM_SCALAR >& from) :
94  _cols_(from._cols_), _positivity_(from._positivity_), _sumIsOne_(from._sumIsOne_) {
95  _rows_.resize(from._rows_.size());
96 
97  for (unsigned int i = 0, end = from._rows_.size(); i < end; i++)
98  _rows_[i] = new LpRow(*from._rows_[i]);
99 
100  GUM_CONS_CPY(LpInterface);
101  }
102 
103  template < typename GUM_SCALAR >
104  LpInterface< GUM_SCALAR >::LpInterface(LpInterface< GUM_SCALAR >&& from) :
105  _positivity_(from._positivity_), _sumIsOne_(from._sumIsOne_) {
106  _rows_.swap(from._rows_);
107  _cols_.swap(from._cols_);
108  GUM_CONS_CPY(LpInterface);
109  }
110 
111  template < typename GUM_SCALAR >
112  LpInterface< GUM_SCALAR >::~LpInterface() {
113  for (const auto row: _rows_)
114  delete row;
115 
116  GUM_DESTRUCTOR(LpInterface);
117  }
118 
119  template < typename GUM_SCALAR >
120  LpInterface< GUM_SCALAR >&
121  LpInterface< GUM_SCALAR >::operator=(const LpInterface< GUM_SCALAR >& from) {
122  /// faster than clear (), delete only rows
123  for (const auto& row: _rows_)
124  delete row;
125 
126  _rows_.clear();
127  _rows_.shrink_to_fit();
128 
129  _rows_.resize(from._rows_.size());
130 
131  for (unsigned int i = 0, end = from._rows_.size(); i < end; i++)
132  _rows_[i] = new LpRow(*from._rows_[i]);
133 
134  _cols_ = from._cols_;
135  _positivity_ = from._positivity_;
136  _sumIsOne_ = from._sumIsOne_;
137 
138  return *this;
139  }
140 
141  template < typename GUM_SCALAR >
142  LpInterface< GUM_SCALAR >&
143  LpInterface< GUM_SCALAR >::operator=(LpInterface< GUM_SCALAR >&& from) {
144  _rows_.swap(from._rows_);
145  _cols_.swap(from._cols_);
146 
147  _positivity_ = from._positivity_;
148  _sumIsOne_ = from._sumIsOne_;
149 
150  return *this;
151  }
152 
153  template < typename T >
154  std::ostream& operator<<(std::ostream& out, const LpInterface< T >& lpi) {
155  out << lpi.toString();
156  return out;
157  }
158 
159  template < typename GUM_SCALAR >
160  LpCol LpInterface< GUM_SCALAR >::addCol() {
161  LpCol col((unsigned int)_cols_.size());
162 
163  _cols_.push_back(col);
164 
165  return col;
166  }
167 
168  template < typename GUM_SCALAR >
169  std::vector< LpCol > LpInterface< GUM_SCALAR >::addCols(const unsigned int& cols) {
170  if (cols < 1)
171  GUM_ERROR(OperationNotAllowed,
172  "LpInterface::addCols ( cols ) : cols "
173  "needs must be equal or greater than 1 : "
174  << cols << " < 1");
175 
176  for (unsigned int i = 0; i < cols; i++) {
177  _cols_.push_back(LpCol((unsigned int)_cols_.size()));
178  }
179 
180  return _cols_;
181  }
182 
183  template < typename GUM_SCALAR >
184  void LpInterface< GUM_SCALAR >::addRow(const LpExpr& expr) {
185  if (!expr._ileft_ && !expr._iright_)
186  GUM_ERROR(OperationNotAllowed,
187  "addRow ( const LpExpr & expr ) : expr : " << expr.toString()
188  << "is not an inequality.");
189 
190  if ((expr._ileft_ && !expr._iright_) || (!expr._ileft_ && expr._iright_)) {
191  _rows_.push_back(new LpRow(expr, _cols_));
192  } else {
193  LpExpr lexpr(expr, true, true, false);
194  LpExpr rexpr(expr, false, true, true);
195 
196  _rows_.push_back(new LpRow(std::move(lexpr),
197  _cols_)); /// lexpr not used anymore, use move constructor
198  _rows_.push_back(new LpRow(std::move(rexpr),
199  _cols_)); /// rexpr not used anymore, use move constructor
200  }
201  }
202 
203  template < typename GUM_SCALAR >
204  void LpInterface< GUM_SCALAR >::addRow(LpExpr&& expr) {
205  if (!expr._ileft_ && !expr._iright_)
206  GUM_ERROR(OperationNotAllowed,
207  "addRow ( const LpExpr & expr ) : expr : " << expr.toString()
208  << "is not an inequality.");
209 
210  if ((expr._ileft_ && !expr._iright_) || (!expr._ileft_ && expr._iright_)) {
211  _rows_.push_back(new LpRow(std::move(expr), _cols_));
212  } else {
213  LpExpr lexpr(std::move(expr), true, true, false);
214 
215  /// expr pointers on maps now are nullptr except right side
216  LpExpr rexpr(std::move(expr), false, false, true);
217 
218  /// rexpr miss middle side, copy it from lexpr
219 
220  *rexpr._mCoeffs_ = *lexpr._mCoeffs_;
221  rexpr._mValue_ = lexpr._mValue_;
222  rexpr._imiddle_ = true;
223 
224  _rows_.push_back(new LpRow(std::move(lexpr),
225  _cols_)); /// lexpr not used anymore, use move constructor
226  _rows_.push_back(new LpRow(std::move(rexpr),
227  _cols_)); /// rexpr not used anymore, use move constructor
228  }
229  }
230 
231  template < typename GUM_SCALAR >
232  void LpInterface< GUM_SCALAR >::addPositivity() {
233  if (_positivity_) return;
234 
235  for (const auto& col: _cols_)
236  addRow(0 <= col);
237 
238  _positivity_ = true;
239  }
240 
241  template < typename GUM_SCALAR >
242  void LpInterface< GUM_SCALAR >::addSumIsOne() {
243  if (_sumIsOne_) return;
244 
245  LpExpr expr;
246 
247  for (const auto& col: _cols_)
248  expr += col;
249 
250  addRow(1 <= std::move(expr) <= 1);
251 
252  _sumIsOne_ = true;
253  }
254 
255  template < typename GUM_SCALAR >
256  void LpInterface< GUM_SCALAR >::addProba() {
257  if (_positivity_ && _sumIsOne_) {
258  return;
259  } else if (_positivity_ && !_sumIsOne_) {
260  addSumIsOne();
261  return;
262  } else if (!_positivity_ && _sumIsOne_) {
263  addPositivity();
264  return;
265  }
266 
267  // we can do both with one loop, don't call the above functions.
268  // addPositivity();
269  // addSumIsOne();
270  LpExpr expr;
271 
272  for (const auto& col: _cols_) {
273  addRow(0 <= col);
274  expr += col;
275  }
276 
277  addRow(1 <= std::move(expr) <= 1);
278 
279  _sumIsOne_ = true;
280  _positivity_ = true;
281  }
282 
283  template < typename GUM_SCALAR >
284  std::vector< std::vector< GUM_SCALAR > > LpInterface< GUM_SCALAR >::solve() {
285  LRSWrapper< GUM_SCALAR > lrs;
286 
287  lrs.setUpH((unsigned int)_cols_.size());
288 
289  std::vector< std::vector< GUM_SCALAR > > lrsMatrix;
290 
291  for (const auto& row: _rows_) {
292  std::vector< GUM_SCALAR > expandedRow(_cols_.size() + 1, 0);
293 
294  expandedRow[0] = row->_cste_;
295 
296  for (const auto& elt: *row->_coeffs_)
297  expandedRow[elt.first.id() + 1] = elt.second;
298 
299  lrsMatrix.push_back(expandedRow);
300  }
301 
302  lrs.fillMatrix(lrsMatrix);
303 
304  lrs.H2V();
305 
306  return lrs.getOutput();
307  }
308 
309  template < typename GUM_SCALAR >
310  std::vector< LpCol > LpInterface< GUM_SCALAR >::getCols() const {
311  return _cols_;
312  }
313 
314  template < typename GUM_SCALAR >
315  std::string LpInterface< GUM_SCALAR >::toString() const {
316  std::ostringstream s;
317 
318  s << std::endl << std::endl << "Variables : " << std::endl;
319 
320  for (const auto& col: _cols_)
321  s << " " << col.toString();
322 
323  s << std::endl;
324 
325  for (const auto& row: _rows_)
326  s << std::endl << row->toString();
327 
328  s << std::endl << std::endl;
329 
330  return s.str();
331  }
332 
333  template < typename GUM_SCALAR >
334  void LpInterface< GUM_SCALAR >::clear() {
335  for (const auto& row: _rows_)
336  delete row;
337 
338  _rows_.clear();
339  _rows_.shrink_to_fit(); /// to really clear content memory, otherwise
340  /// we have
341  /// to wait for (*this) destruction ???
342  /// checked with sizeof( _rows_ ) + sizeof( LpRow ) * _rows_.capacity()
343 
344  _cols_.clear();
345  _cols_.shrink_to_fit();
346 
347  _positivity_ = false;
348  _sumIsOne_ = false;
349  }
350 
351  template < typename GUM_SCALAR >
352  void LpInterface< GUM_SCALAR >::clearRows() {
353  for (const auto& row: _rows_)
354  delete row;
355 
356  _rows_.clear();
357  _rows_.shrink_to_fit();
358 
359  _positivity_ = false;
360  _sumIsOne_ = false;
361  }
362 
363  ///////////////////////////////////////////////////////
364  template < typename T2 >
365  LpExpr operator+(LpExpr&& lhs, const T2& rhs) {
366  LpExpr expr = std::move(lhs);
367  expr += rhs;
368 
369  return expr;
370  }
371 
372  template < typename T2 >
373  LpExpr operator+(const LpExpr& lhs, const T2& rhs) {
374  LpExpr expr(lhs);
375  expr += rhs;
376 
377  return expr;
378  }
379 
380  template < typename T1, forbidden_type< T1, LpExpr > >
381  LpExpr operator+(const T1& lhs, LpExpr&& rhs) {
382  LpExpr expr = std::move(rhs);
383  ;
384  expr += lhs;
385 
386  return expr;
387  }
388 
389  template < typename T1, forbidden_type< T1, LpExpr > >
390  LpExpr operator+(const T1& lhs, LpExpr& rhs) {
391  LpExpr expr(rhs);
392  expr += lhs;
393 
394  return expr;
395  }
396 
397  template < typename T2, forbidden_type< T2, LpExpr > >
398  LpExpr operator+(const LpCol& lhs, const T2& rhs) {
399  LpExpr expr;
400  expr += lhs;
401  expr += rhs;
402 
403  return expr;
404  }
405 
406  template < typename T1, forbidden_type< T1, LpExpr >, forbidden_type< T1, LpCol > >
407  LpExpr operator+(const T1& lhs, const LpCol& rhs) {
408  LpExpr expr;
409  expr += rhs;
410  expr += lhs;
411 
412  return expr;
413  }
414 
415  ///////////////////////////////////////////////////////
416  template < typename T2 >
417  LpExpr operator-(LpExpr&& lhs, const T2& rhs) {
418  LpExpr expr = std::move(lhs);
419  expr -= rhs;
420 
421  return expr;
422  }
423 
424  template < typename T2 >
425  LpExpr operator-(const LpExpr& lhs, const T2& rhs) {
426  LpExpr expr(lhs);
427  expr -= rhs;
428 
429  return expr;
430  }
431 
432  template < typename T1, forbidden_type< T1, LpExpr > >
433  LpExpr operator-(const T1& lhs, LpExpr&& rhs) {
434  LpExpr expr;
435  expr += std::move(rhs);
436  ;
437  expr -= lhs;
438 
439  return expr;
440  }
441 
442  template < typename T1, forbidden_type< T1, LpExpr > >
443  LpExpr operator-(const T1& lhs, LpExpr& rhs) {
444  LpExpr expr;
445  expr += rhs;
446  expr -= lhs;
447 
448  return expr;
449  }
450 
451  template < typename T2, forbidden_type< T2, LpExpr > >
452  LpExpr operator-(const LpCol& lhs, const T2& rhs) {
453  LpExpr expr;
454  expr += lhs;
455  expr -= rhs;
456 
457  return expr;
458  }
459 
460  template < typename T1, forbidden_type< T1, LpExpr >, forbidden_type< T1, LpCol > >
461  LpExpr operator-(const T1& lhs, const LpCol& rhs) {
462  LpExpr expr;
463  expr += rhs;
464  expr -= lhs;
465 
466  return expr;
467  }
468 
469  ///////////////////////////////////////////////////////
470  template < typename SCALAR >
471  INLINE LpExpr LpExpr::multiply(const SCALAR& lhs, const LpCol& rhs) {
472  LpExpr expr;
473  expr._mCoeffs_->insert(rhs, lhs);
474  expr._imiddle_ = true;
475  return expr;
476  }
477 
478  template < typename SCALAR >
479  LpExpr operator*(const SCALAR& lhs, const LpCol& rhs) {
480  return LpExpr::multiply(lhs, rhs);
481  }
482 
483  template < typename SCALAR >
484  LpExpr operator*(const LpCol& lhs, const SCALAR& rhs) {
485  return LpExpr::multiply(rhs, lhs);
486  }
487 
488  ///////////////////////////////////////////////////////
489  template < typename T1, typename T2 >
490  INLINE LpExpr LpExpr::lessThan(T1&& lhs, T2&& rhs) {
491  LpExpr expr;
492  expr._addSide_(std::forward< T1 >(lhs));
493  expr._addSide_(std::forward< T2 >(rhs));
494  return expr;
495  }
496 
497  // const lvalue
498  template < typename T2 >
499  LpExpr operator<=(const LpExpr& lhs, T2&& rhs) {
500  return LpExpr::lessThan(lhs, std::forward< T2 >(rhs));
501  }
502 
503  template < typename T2 >
504  LpExpr operator<=(const LpCol& lhs, T2&& rhs) {
505  return LpExpr::lessThan(lhs, std::forward< T2 >(rhs));
506  }
507 
508  template < typename T1, forbidden_type< T1, LpExpr& >, forbidden_type< T1, LpCol& > >
509  LpExpr operator<=(T1&& lhs, const LpExpr& rhs) {
510  return LpExpr::lessThan(std::forward< T1 >(lhs), rhs);
511  }
512 
513  template < typename T1, forbidden_type< T1, LpExpr& >, forbidden_type< T1, LpCol& > >
514  LpExpr operator<=(T1&& lhs, const LpCol& rhs) {
515  return LpExpr::lessThan(std::forward< T1 >(lhs), rhs);
516  }
517 
518  // rvaue
519  template < typename T2 >
520  LpExpr operator<=(LpExpr&& lhs, T2&& rhs) {
521  return LpExpr::lessThan(std::move(lhs), std::forward< T2 >(rhs));
522  }
523 
524  template < typename T2 >
525  LpExpr operator<=(LpCol&& lhs, T2&& rhs) {
526  return LpExpr::lessThan(std::move(lhs), std::forward< T2 >(rhs));
527  }
528 
529  template < typename T1, forbidden_type< T1, LpExpr >, forbidden_type< T1, LpCol > >
530  LpExpr operator<=(T1&& lhs, LpExpr&& rhs) {
531  return LpExpr::lessThan(std::forward< T1 >(lhs), std::move(rhs));
532  }
533 
534  template < typename T1, forbidden_type< T1, LpExpr >, forbidden_type< T1, LpCol > >
535  LpExpr operator<=(T1&& lhs, LpCol&& rhs) {
536  return LpExpr::lessThan(std::forward< T1 >(lhs), std::move(rhs));
537  }
538  } // namespace lp
539 
540  } // namespace credal
541 
542 } // namespace gum