aGrUM  0.20.2
a C++ library for (probabilistic) graphical models
LpInterface_tpl.h
Go to the documentation of this file.
1 /**
2  *
3  * Copyright 2005-2020 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,
45  "expr::operator+= (expr) : <= present on one side of expr");
46 
47  if (!imiddle__) imiddle__ = true;
48 
49  mValue__ += rhs;
50 
51  return *this;
52  }
53 
54  template < typename T >
55  LpExpr& LpExpr::operator-=(const T& rhs) {
56  if (ileft__ || iright__)
57  GUM_ERROR(OperationNotAllowed,
58  "expr::operator-= (rhs) : <= present in one of expr");
59 
60  if (!imiddle__) imiddle__ = true;
61 
62  mValue__ -= rhs;
63 
64  return *this;
65  }
66 
67  template < typename SCALAR >
68  void LpExpr::addSide__(const SCALAR& from) {
69  if (!ileft__) {
70  lValue__ = from;
71  ileft__ = true;
72  } else if (!imiddle__) {
73  mValue__ = from;
74  imiddle__ = true;
75  } else if (!iright__) {
76  rValue__ = from;
77  iright__ = true;
78  } else
79  GUM_ERROR(OperationNotAllowed,
80  "LpExpr::setSide ( const LpCol & from "
81  ") : too many <= ; no free side");
82  }
83 
84  /**
85  * class LpInterface
86  */
87  template < typename GUM_SCALAR >
88  LpInterface< GUM_SCALAR >::LpInterface() {
89  positivity__ = false;
90  sumIsOne__ = false;
91  GUM_CONSTRUCTOR(LpInterface);
92  }
93 
94  template < typename GUM_SCALAR >
95  LpInterface< GUM_SCALAR >::LpInterface(
96  const LpInterface< GUM_SCALAR >& from) :
97  cols__(from.cols__),
98  positivity__(from.positivity__), sumIsOne__(from.sumIsOne__) {
99  rows__.resize(from.rows__.size());
100 
101  for (unsigned int i = 0, end = from.rows__.size(); i < end; i++)
102  rows__[i] = new LpRow(*from.rows__[i]);
103 
104  GUM_CONS_CPY(LpInterface);
105  }
106 
107  template < typename GUM_SCALAR >
108  LpInterface< GUM_SCALAR >::LpInterface(LpInterface< GUM_SCALAR >&& from) :
109  positivity__(from.positivity__), sumIsOne__(from.sumIsOne__) {
110  rows__.swap(from.rows__);
111  cols__.swap(from.cols__);
112  GUM_CONS_CPY(LpInterface);
113  }
114 
115  template < typename GUM_SCALAR >
116  LpInterface< GUM_SCALAR >::~LpInterface() {
117  for (const auto row: rows__)
118  delete row;
119 
120  GUM_DESTRUCTOR(LpInterface);
121  }
122 
123  template < typename GUM_SCALAR >
124  LpInterface< GUM_SCALAR >& LpInterface< GUM_SCALAR >::operator=(
125  const LpInterface< GUM_SCALAR >& from) {
126  /// faster than clear (), delete only rows
127  for (const auto& row: rows__)
128  delete row;
129 
130  rows__.clear();
131  rows__.shrink_to_fit();
132 
133  rows__.resize(from.rows__.size());
134 
135  for (unsigned int i = 0, end = from.rows__.size(); i < end; i++)
136  rows__[i] = new LpRow(*from.rows__[i]);
137 
138  cols__ = from.cols__;
139  positivity__ = from.positivity__;
140  sumIsOne__ = from.sumIsOne__;
141 
142  return *this;
143  }
144 
145  template < typename GUM_SCALAR >
146  LpInterface< GUM_SCALAR >&
147  LpInterface< GUM_SCALAR >::operator=(LpInterface< GUM_SCALAR >&& from) {
148  rows__.swap(from.rows__);
149  cols__.swap(from.cols__);
150 
151  positivity__ = from.positivity__;
152  sumIsOne__ = from.sumIsOne__;
153 
154  return *this;
155  }
156 
157  template < typename T >
158  std::ostream& operator<<(std::ostream& out, const LpInterface< T >& lpi) {
159  out << lpi.toString();
160  return out;
161  }
162 
163  template < typename GUM_SCALAR >
164  LpCol LpInterface< GUM_SCALAR >::addCol() {
165  LpCol col((unsigned int)cols__.size());
166 
167  cols__.push_back(col);
168 
169  return col;
170  }
171 
172  template < typename GUM_SCALAR >
173  std::vector< LpCol >
174  LpInterface< GUM_SCALAR >::addCols(const unsigned int& cols) {
175  if (cols < 1)
176  GUM_ERROR(OperationNotAllowed,
177  "LpInterface::addCols ( cols ) : cols "
178  "needs must be equal or greater than 1 : "
179  << cols << " < 1");
180 
181  for (unsigned int i = 0; i < cols; i++) {
182  cols__.push_back(LpCol((unsigned int)cols__.size()));
183  }
184 
185  return cols__;
186  }
187 
188  template < typename GUM_SCALAR >
189  void LpInterface< GUM_SCALAR >::addRow(const LpExpr& expr) {
190  if (!expr.ileft__ && !expr.iright__)
191  GUM_ERROR(OperationNotAllowed,
192  "addRow ( const LpExpr & expr ) : expr : "
193  << expr.toString() << "is not an inequality.");
194 
195  if ((expr.ileft__ && !expr.iright__) || (!expr.ileft__ && expr.iright__)) {
196  rows__.push_back(new LpRow(expr, cols__));
197  } else {
198  LpExpr lexpr(expr, true, true, false);
199  LpExpr rexpr(expr, false, true, true);
200 
201  rows__.push_back(
202  new LpRow(std::move(lexpr),
203  cols__)); /// lexpr not used anymore, use move constructor
204  rows__.push_back(
205  new LpRow(std::move(rexpr),
206  cols__)); /// rexpr not used anymore, use move constructor
207  }
208  }
209 
210  template < typename GUM_SCALAR >
211  void LpInterface< GUM_SCALAR >::addRow(LpExpr&& expr) {
212  if (!expr.ileft__ && !expr.iright__)
213  GUM_ERROR(OperationNotAllowed,
214  "addRow ( const LpExpr & expr ) : expr : "
215  << expr.toString() << "is not an inequality.");
216 
217  if ((expr.ileft__ && !expr.iright__) || (!expr.ileft__ && expr.iright__)) {
218  rows__.push_back(new LpRow(std::move(expr), cols__));
219  } else {
220  LpExpr lexpr(std::move(expr), true, true, false);
221 
222  /// expr pointers on maps now are nullptr except right side
223  LpExpr rexpr(std::move(expr), false, false, true);
224 
225  /// rexpr miss middle side, copy it from lexpr
226 
227  *rexpr.mCoeffs__ = *lexpr.mCoeffs__;
228  rexpr.mValue__ = lexpr.mValue__;
229  rexpr.imiddle__ = true;
230 
231  rows__.push_back(
232  new LpRow(std::move(lexpr),
233  cols__)); /// lexpr not used anymore, use move constructor
234  rows__.push_back(
235  new LpRow(std::move(rexpr),
236  cols__)); /// rexpr not used anymore, use move constructor
237  }
238  }
239 
240  template < typename GUM_SCALAR >
241  void LpInterface< GUM_SCALAR >::addPositivity() {
242  if (positivity__) return;
243 
244  for (const auto& col: cols__)
245  addRow(0 <= col);
246 
247  positivity__ = true;
248  }
249 
250  template < typename GUM_SCALAR >
251  void LpInterface< GUM_SCALAR >::addSumIsOne() {
252  if (sumIsOne__) return;
253 
254  LpExpr expr;
255 
256  for (const auto& col: cols__)
257  expr += col;
258 
259  addRow(1 <= std::move(expr) <= 1);
260 
261  sumIsOne__ = true;
262  }
263 
264  template < typename GUM_SCALAR >
265  void LpInterface< GUM_SCALAR >::addProba() {
266  if (positivity__ && sumIsOne__) {
267  return;
268  } else if (positivity__ && !sumIsOne__) {
269  addSumIsOne();
270  return;
271  } else if (!positivity__ && sumIsOne__) {
272  addPositivity();
273  return;
274  }
275 
276  // we can do both with one loop, don't call the above functions.
277  // addPositivity();
278  // addSumIsOne();
279  LpExpr expr;
280 
281  for (const auto& col: cols__) {
282  addRow(0 <= col);
283  expr += col;
284  }
285 
286  addRow(1 <= std::move(expr) <= 1);
287 
288  sumIsOne__ = true;
289  positivity__ = true;
290  }
291 
292  template < typename GUM_SCALAR >
293  std::vector< std::vector< GUM_SCALAR > > LpInterface< GUM_SCALAR >::solve() {
294  LRSWrapper< GUM_SCALAR > lrs;
295 
296  lrs.setUpH((unsigned int)cols__.size());
297 
298  std::vector< std::vector< GUM_SCALAR > > lrsMatrix;
299 
300  for (const auto& row: rows__) {
301  std::vector< GUM_SCALAR > expandedRow(cols__.size() + 1, 0);
302 
303  expandedRow[0] = row->cste__;
304 
305  for (const auto& elt: *row->coeffs__)
306  expandedRow[elt.first.id() + 1] = elt.second;
307 
308  lrsMatrix.push_back(expandedRow);
309  }
310 
311  lrs.fillMatrix(lrsMatrix);
312 
313  lrs.H2V();
314 
315  return lrs.getOutput();
316  }
317 
318  template < typename GUM_SCALAR >
319  std::vector< LpCol > LpInterface< GUM_SCALAR >::getCols() const {
320  return cols__;
321  }
322 
323  template < typename GUM_SCALAR >
324  std::string LpInterface< GUM_SCALAR >::toString() const {
325  std::ostringstream s;
326 
327  s << std::endl << std::endl << "Variables : " << std::endl;
328 
329  for (const auto& col: cols__)
330  s << " " << col.toString();
331 
332  s << std::endl;
333 
334  for (const auto& row: rows__)
335  s << std::endl << row->toString();
336 
337  s << std::endl << std::endl;
338 
339  return s.str();
340  }
341 
342  template < typename GUM_SCALAR >
343  void LpInterface< GUM_SCALAR >::clear() {
344  for (const auto& row: rows__)
345  delete row;
346 
347  rows__.clear();
348  rows__.shrink_to_fit(); /// to really clear content memory, otherwise
349  /// we have
350  /// to wait for (*this) destruction ???
351  /// checked with sizeof( rows__ ) + sizeof( LpRow ) * rows__.capacity()
352 
353  cols__.clear();
354  cols__.shrink_to_fit();
355 
356  positivity__ = false;
357  sumIsOne__ = false;
358  }
359 
360  template < typename GUM_SCALAR >
361  void LpInterface< GUM_SCALAR >::clearRows() {
362  for (const auto& row: rows__)
363  delete row;
364 
365  rows__.clear();
366  rows__.shrink_to_fit();
367 
368  positivity__ = false;
369  sumIsOne__ = false;
370  }
371 
372  ///////////////////////////////////////////////////////
373  template < typename T2 >
374  LpExpr operator+(LpExpr&& lhs, const T2& rhs) {
375  LpExpr expr = std::move(lhs);
376  expr += rhs;
377 
378  return expr;
379  }
380 
381  template < typename T2 >
382  LpExpr operator+(const LpExpr& lhs, const T2& rhs) {
383  LpExpr expr(lhs);
384  expr += rhs;
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 = std::move(rhs);
392  ;
393  expr += lhs;
394 
395  return expr;
396  }
397 
398  template < typename T1, forbidden_type< T1, LpExpr > >
399  LpExpr operator+(const T1& lhs, LpExpr& rhs) {
400  LpExpr expr(rhs);
401  expr += lhs;
402 
403  return expr;
404  }
405 
406  template < typename T2, forbidden_type< T2, LpExpr > >
407  LpExpr operator+(const LpCol& lhs, const T2& rhs) {
408  LpExpr expr;
409  expr += lhs;
410  expr += rhs;
411 
412  return expr;
413  }
414 
415  template < typename T1,
416  forbidden_type< T1, LpExpr >,
417  forbidden_type< T1, LpCol > >
418  LpExpr operator+(const T1& lhs, const LpCol& rhs) {
419  LpExpr expr;
420  expr += rhs;
421  expr += lhs;
422 
423  return expr;
424  }
425 
426  ///////////////////////////////////////////////////////
427  template < typename T2 >
428  LpExpr operator-(LpExpr&& lhs, const T2& rhs) {
429  LpExpr expr = std::move(lhs);
430  expr -= rhs;
431 
432  return expr;
433  }
434 
435  template < typename T2 >
436  LpExpr operator-(const LpExpr& lhs, const T2& rhs) {
437  LpExpr expr(lhs);
438  expr -= rhs;
439 
440  return expr;
441  }
442 
443  template < typename T1, forbidden_type< T1, LpExpr > >
444  LpExpr operator-(const T1& lhs, LpExpr&& rhs) {
445  LpExpr expr;
446  expr += std::move(rhs);
447  ;
448  expr -= lhs;
449 
450  return expr;
451  }
452 
453  template < typename T1, forbidden_type< T1, LpExpr > >
454  LpExpr operator-(const T1& lhs, LpExpr& rhs) {
455  LpExpr expr;
456  expr += rhs;
457  expr -= lhs;
458 
459  return expr;
460  }
461 
462  template < typename T2, forbidden_type< T2, LpExpr > >
463  LpExpr operator-(const LpCol& lhs, const T2& rhs) {
464  LpExpr expr;
465  expr += lhs;
466  expr -= rhs;
467 
468  return expr;
469  }
470 
471  template < typename T1,
472  forbidden_type< T1, LpExpr >,
473  forbidden_type< T1, LpCol > >
474  LpExpr operator-(const T1& lhs, const LpCol& rhs) {
475  LpExpr expr;
476  expr += rhs;
477  expr -= lhs;
478 
479  return expr;
480  }
481 
482  ///////////////////////////////////////////////////////
483  template < typename SCALAR >
484  INLINE LpExpr LpExpr::multiply(const SCALAR& lhs, const LpCol& rhs) {
485  LpExpr expr;
486  expr.mCoeffs__->insert(rhs, lhs);
487  expr.imiddle__ = true;
488  return expr;
489  }
490 
491  template < typename SCALAR >
492  LpExpr operator*(const SCALAR& lhs, const LpCol& rhs) {
493  return LpExpr::multiply(lhs, rhs);
494  }
495 
496  template < typename SCALAR >
497  LpExpr operator*(const LpCol& lhs, const SCALAR& rhs) {
498  return LpExpr::multiply(rhs, lhs);
499  }
500 
501  ///////////////////////////////////////////////////////
502  template < typename T1, typename T2 >
503  INLINE LpExpr LpExpr::lessThan(T1&& lhs, T2&& rhs) {
504  LpExpr expr;
505  expr.addSide__(std::forward< T1 >(lhs));
506  expr.addSide__(std::forward< T2 >(rhs));
507  return expr;
508  }
509 
510  // const lvalue
511  template < typename T2 >
512  LpExpr operator<=(const LpExpr& lhs, T2&& rhs) {
513  return LpExpr::lessThan(lhs, std::forward< T2 >(rhs));
514  }
515 
516  template < typename T2 >
517  LpExpr operator<=(const LpCol& lhs, T2&& rhs) {
518  return LpExpr::lessThan(lhs, std::forward< T2 >(rhs));
519  }
520 
521  template < typename T1,
522  forbidden_type< T1, LpExpr& >,
523  forbidden_type< T1, LpCol& > >
524  LpExpr operator<=(T1&& lhs, const LpExpr& rhs) {
525  return LpExpr::lessThan(std::forward< T1 >(lhs), rhs);
526  }
527 
528  template < typename T1,
529  forbidden_type< T1, LpExpr& >,
530  forbidden_type< T1, LpCol& > >
531  LpExpr operator<=(T1&& lhs, const LpCol& rhs) {
532  return LpExpr::lessThan(std::forward< T1 >(lhs), rhs);
533  }
534 
535  // rvaue
536  template < typename T2 >
537  LpExpr operator<=(LpExpr&& lhs, T2&& rhs) {
538  return LpExpr::lessThan(std::move(lhs), std::forward< T2 >(rhs));
539  }
540 
541  template < typename T2 >
542  LpExpr operator<=(LpCol&& lhs, T2&& rhs) {
543  return LpExpr::lessThan(std::move(lhs), std::forward< T2 >(rhs));
544  }
545 
546  template < typename T1,
547  forbidden_type< T1, LpExpr >,
548  forbidden_type< T1, LpCol > >
549  LpExpr operator<=(T1&& lhs, LpExpr&& rhs) {
550  return LpExpr::lessThan(std::forward< T1 >(lhs), std::move(rhs));
551  }
552 
553  template < typename T1,
554  forbidden_type< T1, LpExpr >,
555  forbidden_type< T1, LpCol > >
556  LpExpr operator<=(T1&& lhs, LpCol&& rhs) {
557  return LpExpr::lessThan(std::forward< T1 >(lhs), std::move(rhs));
558  }
559  } // namespace lp
560 
561  } // namespace credal
562 
563 } // namespace gum