File Coverage

/usr/local/lib/perl5/site_perl/5.26.1/x86_64-linux/auto/share/dist/Alien-Kiwisolver/include/kiwi/solverimpl.h
Criterion Covered Total %
statement 200 317 63.0
branch 192 422 45.5
condition n/a
subroutine n/a
pod n/a
total 392 739 53.0


line stmt bran cond sub pod time code
1             /*-----------------------------------------------------------------------------
2             | Copyright (c) 2013-2017, Nucleic Development Team.
3             |
4             | Distributed under the terms of the Modified BSD License.
5             |
6             | The full license is in the file LICENSE, distributed with this software.
7             |----------------------------------------------------------------------------*/
8             #pragma once
9             #include <algorithm>
10             #include <limits>
11             #include <memory>
12             #include <vector>
13             #include "constraint.h"
14             #include "errors.h"
15             #include "expression.h"
16             #include "maptype.h"
17             #include "row.h"
18             #include "symbol.h"
19             #include "term.h"
20             #include "util.h"
21             #include "variable.h"
22              
23              
24             namespace kiwi
25             {
26              
27             namespace impl
28             {
29              
30             class SolverImpl
31             {
32             friend class DebugHelper;
33              
34 17178           struct Tag
35             {
36             Symbol marker;
37             Symbol other;
38             };
39              
40 38662           struct EditInfo
41             {
42             Tag tag;
43             Constraint constraint;
44             double constant;
45             };
46              
47             using VarMap = MapType<Variable, Symbol>;
48              
49             using RowMap = MapType<Symbol, Row*>;
50              
51             using CnMap = MapType<Constraint, Tag>;
52              
53             using EditMap = MapType<Variable, EditInfo>;
54              
55             struct DualOptimizeGuard
56             {
57 619           DualOptimizeGuard( SolverImpl& impl ) : m_impl( impl ) {}
58 1238           ~DualOptimizeGuard() { m_impl.dualOptimize(); }
59             SolverImpl& m_impl;
60             };
61              
62             public:
63              
64 5 50         SolverImpl() : m_objective( new Row() ), m_id_tick( 1 ) {}
    50          
    50          
    50          
    50          
    50          
65              
66             SolverImpl( const SolverImpl& ) = delete;
67              
68             SolverImpl( SolverImpl&& ) = delete;
69              
70 10           ~SolverImpl() { clearRows(); }
71              
72             /* Add a constraint to the solver.
73              
74             Throws
75             ------
76             DuplicateConstraint
77             The given constraint has already been added to the solver.
78              
79             UnsatisfiableConstraint
80             The given constraint is required and cannot be satisfied.
81              
82             */
83 3216           void addConstraint( const Constraint& constraint )
84             {
85 3216 50         if( m_cns.find( constraint ) != m_cns.end() )
    50          
86 0 0         throw DuplicateConstraint( constraint );
87              
88             // Creating a row causes symbols to be reserved for the variables
89             // in the constraint. If this method exits with an exception,
90             // then its possible those variables will linger in the var map.
91             // Since its likely that those variables will be used in other
92             // constraints and since exceptional conditions are uncommon,
93             // i'm not too worried about aggressive cleanup of the var map.
94 3216           Tag tag;
95 6432 50         std::unique_ptr<Row> rowptr( createRow( constraint, tag ) );
96 3216 50         Symbol subject( chooseSubject( *rowptr, tag ) );
97              
98             // If chooseSubject could not find a valid entering symbol, one
99             // last option is available if the entire row is composed of
100             // dummy variables. If the constant of the row is zero, then
101             // this represents redundant constraints and the new dummy
102             // marker can enter the basis. If the constant is non-zero,
103             // then it represents an unsatisfiable constraint.
104 3216 100         if( subject.type() == Symbol::Invalid && allDummies( *rowptr ) )
    50          
    50          
    50          
105             {
106 0 0         if( !nearZero( rowptr->constant() ) )
107 0 0         throw UnsatisfiableConstraint( constraint );
108             else
109 0           subject = tag.marker;
110             }
111              
112             // If an entering symbol still isn't found, then the row must
113             // be added using an artificial variable. If that fails, then
114             // the row represents an unsatisfiable constraint.
115 3216 100         if( subject.type() == Symbol::Invalid )
116             {
117 415 50         if( !addWithArtificialVariable( *rowptr ) )
    50          
118 0 0         throw UnsatisfiableConstraint( constraint );
119             }
120             else
121             {
122 2801 50         rowptr->solveFor( subject );
123 2801 50         substitute( subject, *rowptr );
124 2801 50         m_rows[ subject ] = rowptr.release();
125             }
126              
127 3216 50         m_cns[ constraint ] = tag;
128              
129             // Optimizing after each constraint is added performs less
130             // aggregate work due to a smaller average system size. It
131             // also ensures the solver remains in a consistent state.
132 3216 50         optimize( *m_objective );
133 3216           }
134              
135             /* Remove a constraint from the solver.
136              
137             Throws
138             ------
139             UnknownConstraint
140             The given constraint has not been added to the solver.
141              
142             */
143 0           void removeConstraint( const Constraint& constraint )
144             {
145 0 0         auto cn_it = m_cns.find( constraint );
146 0 0         if( cn_it == m_cns.end() )
147 0 0         throw UnknownConstraint( constraint );
148              
149 0           Tag tag( cn_it->second );
150 0 0         m_cns.erase( cn_it );
151              
152             // Remove the error effects from the objective function
153             // *before* pivoting, or substitutions into the objective
154             // will lead to incorrect solver results.
155 0 0         removeConstraintEffects( constraint, tag );
156              
157             // If the marker is basic, simply drop the row. Otherwise,
158             // pivot the marker into the basis and then drop the row.
159 0 0         auto row_it = m_rows.find( tag.marker );
160 0 0         if( row_it != m_rows.end() )
161             {
162 0           std::unique_ptr<Row> rowptr( row_it->second );
163 0 0         m_rows.erase( row_it );
164             }
165             else
166             {
167 0 0         row_it = getMarkerLeavingRow( tag.marker );
168 0 0         if( row_it == m_rows.end() )
169 0 0         throw InternalSolverError( "failed to find leaving row" );
170 0           Symbol leaving( row_it->first );
171 0           std::unique_ptr<Row> rowptr( row_it->second );
172 0 0         m_rows.erase( row_it );
173 0 0         rowptr->solveFor( leaving, tag.marker );
174 0 0         substitute( tag.marker, *rowptr );
175             }
176              
177             // Optimizing after each constraint is removed ensures that the
178             // solver remains consistent. It makes the solver api easier to
179             // use at a small tradeoff for speed.
180 0 0         optimize( *m_objective );
181 0           }
182              
183             /* Test whether a constraint has been added to the solver.
184              
185             */
186 0           bool hasConstraint( const Constraint& constraint ) const
187             {
188 0 0         return m_cns.find( constraint ) != m_cns.end();
189             }
190              
191             /* Add an edit variable to the solver.
192              
193             This method should be called before the `suggestValue` method is
194             used to supply a suggested value for the given edit variable.
195              
196             Throws
197             ------
198             DuplicateEditVariable
199             The given edit variable has already been added to the solver.
200              
201             BadRequiredStrength
202             The given strength is >= required.
203              
204             */
205 719           void addEditVariable( const Variable& variable, double strength )
206             {
207 719 50         if( m_edits.find( variable ) != m_edits.end() )
    50          
208 0 0         throw DuplicateEditVariable( variable );
209 719           strength = strength::clip( strength );
210 719 50         if( strength == strength::required )
211 0           throw BadRequiredStrength();
212 1438 50         Constraint cn( Expression( variable ), OP_EQ, strength );
    50          
    50          
213 719 50         addConstraint( cn );
214 1438 50         EditInfo info;
215 719 50         info.tag = m_cns[ cn ];
216 719 50         info.constraint = cn;
217 719           info.constant = 0.0;
218 719 50         m_edits[ variable ] = info;
    50          
219 719           }
220              
221             /* Remove an edit variable from the solver.
222              
223             Throws
224             ------
225             UnknownEditVariable
226             The given edit variable has not been added to the solver.
227              
228             */
229 0           void removeEditVariable( const Variable& variable )
230             {
231 0 0         auto it = m_edits.find( variable );
232 0 0         if( it == m_edits.end() )
233 0 0         throw UnknownEditVariable( variable );
234 0 0         removeConstraint( it->second.constraint );
235 0 0         m_edits.erase( it );
236 0           }
237              
238             /* Test whether an edit variable has been added to the solver.
239              
240             */
241 0           bool hasEditVariable( const Variable& variable ) const
242             {
243 0 0         return m_edits.find( variable ) != m_edits.end();
244             }
245              
246             /* Suggest a value for the given edit variable.
247              
248             This method should be used after an edit variable as been added to
249             the solver in order to suggest the value for that variable.
250              
251             Throws
252             ------
253             UnknownEditVariable
254             The given edit variable has not been added to the solver.
255              
256             */
257 619           void suggestValue( const Variable& variable, double value )
258             {
259 619 50         auto it = m_edits.find( variable );
260 619 50         if( it == m_edits.end() )
261 0 0         throw UnknownEditVariable( variable );
262              
263 1238           DualOptimizeGuard guard( *this );
264 619           EditInfo& info = it->second;
265 619           double delta = value - info.constant;
266 619           info.constant = value;
267              
268             // Check first if the positive error variable is basic.
269 619 50         auto row_it = m_rows.find( info.tag.marker );
270 619 100         if( row_it != m_rows.end() )
271             {
272 3 100         if( row_it->second->add( -delta ) < 0.0 )
273 2 50         m_infeasible_rows.push_back( row_it->first );
274 3           return;
275             }
276              
277             // Check next if the negative error variable is basic.
278 616 50         row_it = m_rows.find( info.tag.other );
279 616 50         if( row_it != m_rows.end() )
280             {
281 0 0         if( row_it->second->add( delta ) < 0.0 )
282 0 0         m_infeasible_rows.push_back( row_it->first );
283 0           return;
284             }
285              
286             // Otherwise update each row where the error variables exist.
287 277841 100         for (const auto & rowPair : m_rows)
    100          
288             {
289 277222 50         double coeff = rowPair.second->coefficientFor( info.tag.marker );
290 279051 100         if( coeff != 0.0 &&
    100          
291 279051 100         rowPair.second->add( delta * coeff ) < 0.0 &&
    50          
292 3           rowPair.first.type() != Symbol::External )
293 3 50         m_infeasible_rows.push_back( rowPair.first );
294             }
295             }
296              
297             /* Update the values of the external solver variables.
298              
299             */
300 8           void updateVariables()
301             {
302 8           auto row_end = m_rows.end();
303              
304 3718 100         for (auto &varPair : m_vars)
305             {
306 3710           Variable& var = varPair.first;
307 3710 50         auto row_it = m_rows.find( varPair.second );
308 3710 50         if( row_it == row_end )
309 0 0         var.setValue( 0.0 );
310             else
311 3710 50         var.setValue( row_it->second->constant() );
312             }
313 8           }
314              
315             /* Reset the solver to the empty starting condition.
316              
317             This method resets the internal solver state to the empty starting
318             condition, as if no constraints or edit variables have been added.
319             This can be faster than deleting the solver and creating a new one
320             when the entire system must change, since it can avoid unecessary
321             heap (de)allocations.
322              
323             */
324 0           void reset()
325             {
326 0           clearRows();
327 0           m_cns.clear();
328 0           m_vars.clear();
329 0           m_edits.clear();
330 0           m_infeasible_rows.clear();
331 0 0         m_objective.reset( new Row() );
332 0           m_artificial.reset();
333 0           m_id_tick = 1;
334 0           }
335              
336             SolverImpl& operator=( const SolverImpl& ) = delete;
337              
338             SolverImpl& operator=( SolverImpl&& ) = delete;
339              
340             private:
341              
342             struct RowDeleter
343             {
344             template<typename T>
345 6432 50         void operator()( T& pair ) { delete pair.second; }
346             };
347              
348 5           void clearRows()
349             {
350 5 50         std::for_each( m_rows.begin(), m_rows.end(), RowDeleter() );
351 5           m_rows.clear();
352 5           }
353              
354             /* Get the symbol for the given variable.
355              
356             If a symbol does not exist for the variable, one will be created.
357              
358             */
359 6121           Symbol getVarSymbol( const Variable& variable )
360             {
361 6121 50         auto it = m_vars.find( variable );
362 6121 100         if( it != m_vars.end() )
363 3638           return it->second;
364 2483           Symbol symbol( Symbol::External, m_id_tick++ );
365 2483 50         m_vars[ variable ] = symbol;
366 6121           return symbol;
367             }
368              
369             /* Create a new Row object for the given constraint.
370              
371             The terms in the constraint will be converted to cells in the row.
372             Any term in the constraint with a coefficient of zero is ignored.
373             This method uses the `getVarSymbol` method to get the symbol for
374             the variables added to the row. If the symbol for a given cell
375             variable is basic, the cell variable will be substituted with the
376             basic row.
377              
378             The necessary slack and error variables will be added to the row.
379             If the constant for the row is negative, the sign for the row
380             will be inverted so the constant becomes positive.
381              
382             The tag will be updated with the marker and error symbols to use
383             for tracking the movement of the constraint in the tableau.
384              
385             */
386 3216           std::unique_ptr<Row> createRow( const Constraint& constraint, Tag& tag )
387             {
388 3216           const Expression& expr( constraint.expression() );
389 3216 50         std::unique_ptr<Row> row( new Row( expr.constant() ) );
    50          
390              
391             // Substitute the current basic variables into the row.
392 9337 100         for (const auto &term : expr.terms())
393             {
394 6121 50         if( !nearZero( term.coefficient() ) )
395             {
396 6121 50         Symbol symbol( getVarSymbol( term.variable() ) );
397 6121 50         auto row_it = m_rows.find( symbol );
398 6121 100         if( row_it != m_rows.end() )
399 2639 50         row->insert( *row_it->second, term.coefficient() );
400             else
401 6121 50         row->insert( symbol, term.coefficient() );
402             }
403             }
404              
405             // Add the necessary slack, error, and dummy variables.
406 3216 50         switch( constraint.op() )
407             {
408             case OP_LE:
409             case OP_GE:
410             {
411 759 50         double coeff = constraint.op() == OP_LE ? 1.0 : -1.0;
    100          
412 759           Symbol slack( Symbol::Slack, m_id_tick++ );
413 759           tag.marker = slack;
414 759 50         row->insert( slack, coeff );
415 759 50         if( constraint.strength() < strength::required )
    50          
416             {
417 0           Symbol error( Symbol::Error, m_id_tick++ );
418 0           tag.other = error;
419 0 0         row->insert( error, -coeff );
420 0 0         m_objective->insert( error, constraint.strength() );
    0          
421             }
422 759           break;
423             }
424             case OP_EQ:
425             {
426 2457 50         if( constraint.strength() < strength::required )
    100          
427             {
428 720           Symbol errplus( Symbol::Error, m_id_tick++ );
429 720           Symbol errminus( Symbol::Error, m_id_tick++ );
430 720           tag.marker = errplus;
431 720           tag.other = errminus;
432 720 50         row->insert( errplus, -1.0 ); // v = eplus - eminus
433 720 50         row->insert( errminus, 1.0 ); // v - eplus + eminus = 0
434 720 50         m_objective->insert( errplus, constraint.strength() );
    50          
435 720 50         m_objective->insert( errminus, constraint.strength() );
    50          
436             }
437             else
438             {
439 1737           Symbol dummy( Symbol::Dummy, m_id_tick++ );
440 1737           tag.marker = dummy;
441 1737 50         row->insert( dummy );
442             }
443 2457           break;
444             }
445             }
446              
447             // Ensure the row as a positive constant.
448 3216 100         if( row->constant() < 0.0 )
449 578 50         row->reverseSign();
450              
451 3216           return row;
452             }
453              
454             /* Choose the subject for solving for the row.
455              
456             This method will choose the best subject for using as the solve
457             target for the row. An invalid symbol will be returned if there
458             is no valid target.
459              
460             The symbols are chosen according to the following precedence:
461              
462             1) The first symbol representing an external variable.
463             2) A negative slack or error tag variable.
464              
465             If a subject cannot be found, an invalid symbol will be returned.
466              
467             */
468 3216           Symbol chooseSubject( const Row& row, const Tag& tag ) const
469             {
470 25153 100         for (const auto &cellPair : row.cells())
471             {
472 24420 100         if( cellPair.first.type() == Symbol::External )
473 2483           return cellPair.first;
474             }
475 733 100         if( tag.marker.type() == Symbol::Slack || tag.marker.type() == Symbol::Error )
    100          
    100          
476             {
477 508 100         if( row.coefficientFor( tag.marker ) < 0.0 )
478 316           return tag.marker;
479             }
480 417 50         if( tag.other.type() == Symbol::Slack || tag.other.type() == Symbol::Error )
    100          
    100          
481             {
482 2 50         if( row.coefficientFor( tag.other ) < 0.0 )
483 2           return tag.other;
484             }
485 415           return Symbol();
486             }
487              
488             /* Add the row to the tableau using an artificial variable.
489              
490             This will return false if the constraint cannot be satisfied.
491              
492             */
493 415           bool addWithArtificialVariable( const Row& row )
494             {
495             // Create and add the artificial variable to the tableau
496 415           Symbol art( Symbol::Slack, m_id_tick++ );
497 415 50         m_rows[ art ] = new Row( row );
    50          
    50          
498 415 50         m_artificial.reset( new Row( row ) );
    50          
499              
500             // Optimize the artificial objective. This is successful
501             // only if the artificial objective is optimized to zero.
502 415 50         optimize( *m_artificial );
503 415           bool success = nearZero( m_artificial->constant() );
504 415           m_artificial.reset();
505              
506             // If the artificial variable is not basic, pivot the row so that
507             // it becomes basic. If the row is constant, exit early.
508 415 50         auto it = m_rows.find( art );
509 415 50         if( it != m_rows.end() )
510             {
511 0           std::unique_ptr<Row> rowptr( it->second );
512 0 0         m_rows.erase( it );
513 0 0         if( rowptr->cells().empty() )
514 0           return success;
515 0 0         Symbol entering( anyPivotableSymbol( *rowptr ) );
516 0 0         if( entering.type() == Symbol::Invalid )
517 0           return false; // unsatisfiable (will this ever happen?)
518 0 0         rowptr->solveFor( art, entering );
519 0 0         substitute( entering, *rowptr );
520 0 0         m_rows[ entering ] = rowptr.release();
    0          
521             }
522              
523             // Remove the artificial variable from the tableau.
524 172254 100         for (auto &rowPair : m_rows)
525 171839 50         rowPair.second->remove(art);
526              
527 415 50         m_objective->remove( art );
528 415           return success;
529             }
530              
531             /* Substitute the parametric symbol with the given row.
532              
533             This method will substitute all instances of the parametric symbol
534             in the tableau and the objective function with the given row.
535              
536             */
537 10753           void substitute( const Symbol& symbol, const Row& row )
538             {
539 3075174 100         for( auto& rowPair : m_rows )
540             {
541 3064421 50         rowPair.second->substitute( symbol, row );
542 3894969 100         if( rowPair.first.type() != Symbol::External &&
543 830548           rowPair.second->constant() < 0.0 )
544 2 50         m_infeasible_rows.push_back( rowPair.first );
545             }
546 10753           m_objective->substitute( symbol, row );
547 10753 100         if( m_artificial.get() )
548 2937           m_artificial->substitute( symbol, row );
549 10753           }
550              
551             /* Optimize the system for the given objective function.
552              
553             This method performs iterations of Phase 2 of the simplex method
554             until the objective function reaches a minimum.
555              
556             Throws
557             ------
558             InternalSolverError
559             The value of the objective function is unbounded.
560              
561             */
562 3631           void optimize( const Row& objective )
563             {
564 7946           while( true )
565             {
566 11577 50         Symbol entering( getEnteringSymbol( objective ) );
567 11577 100         if( entering.type() == Symbol::Invalid )
568 3631           return;
569 7946 50         auto it = getLeavingRow( entering );
570 7946 50         if( it == m_rows.end() )
571 0 0         throw InternalSolverError( "The objective is unbounded." );
572             // pivot the entering symbol into the basis
573 7946           Symbol leaving( it->first );
574 7946           Row* row = it->second;
575 7946 50         m_rows.erase( it );
576 7946 50         row->solveFor( leaving, entering );
577 7946 50         substitute( entering, *row );
578 7946 50         m_rows[ entering ] = row;
579             }
580             }
581              
582             /* Optimize the system using the dual of the simplex method.
583              
584             The current state of the system should be such that the objective
585             function is optimal, but not feasible. This method will perform
586             an iteration of the dual simplex method to make the solution both
587             optimal and feasible.
588              
589             Throws
590             ------
591             InternalSolverError
592             The system cannot be dual optimized.
593              
594             */
595 619           void dualOptimize()
596             {
597 626 100         while( !m_infeasible_rows.empty() )
598             {
599              
600 7           Symbol leaving( m_infeasible_rows.back() );
601 7           m_infeasible_rows.pop_back();
602 7 50         auto it = m_rows.find( leaving );
603 13 50         if( it != m_rows.end() && !nearZero( it->second->constant() ) &&
    50          
    100          
604 6           it->second->constant() < 0.0 )
605             {
606 6 50         Symbol entering( getDualEnteringSymbol( *it->second ) );
607 6 50         if( entering.type() == Symbol::Invalid )
608 0 0         throw InternalSolverError( "Dual optimize failed." );
609             // pivot the entering symbol into the basis
610 6           Row* row = it->second;
611 6 50         m_rows.erase( it );
612 6 50         row->solveFor( leaving, entering );
613 6 50         substitute( entering, *row );
614 6 50         m_rows[ entering ] = row;
615             }
616             }
617 619           }
618              
619             /* Compute the entering variable for a pivot operation.
620              
621             This method will return first symbol in the objective function which
622             is non-dummy and has a coefficient less than zero. If no symbol meets
623             the criteria, it means the objective function is at a minimum, and an
624             invalid symbol is returned.
625              
626             */
627 11577           Symbol getEnteringSymbol( const Row& objective ) const
628             {
629 1606101 100         for (const auto &cellPair : objective.cells())
630             {
631 1602470 100         if( cellPair.first.type() != Symbol::Dummy && cellPair.second < 0.0 )
    100          
    100          
632 7946           return cellPair.first;
633             }
634 3631           return Symbol();
635             }
636              
637             /* Compute the entering symbol for the dual optimize operation.
638              
639             This method will return the symbol in the row which has a positive
640             coefficient and yields the minimum ratio for its respective symbol
641             in the objective function. The provided row *must* be infeasible.
642             If no symbol is found which meats the criteria, an invalid symbol
643             is returned.
644              
645             */
646 6           Symbol getDualEnteringSymbol( const Row& row ) const
647             {
648 6           Symbol entering;
649 6           double ratio = std::numeric_limits<double>::max();
650 2687 100         for (const auto &cellPair : row.cells())
651             {
652 2681 100         if( cellPair.second > 0.0 && cellPair.first.type() != Symbol::Dummy )
    100          
    100          
653             {
654 589 50         double coeff = m_objective->coefficientFor( cellPair.first );
655 589           double r = coeff / cellPair.second;
656 589 100         if( r < ratio )
657             {
658 6           ratio = r;
659 589           entering = cellPair.first;
660             }
661             }
662             }
663 6           return entering;
664             }
665              
666             /* Get the first Slack or Error symbol in the row.
667              
668             If no such symbol is present, and Invalid symbol will be returned.
669              
670             */
671 0           Symbol anyPivotableSymbol( const Row& row ) const
672             {
673 0 0         for (const auto &cellPair : row.cells())
674             {
675 0           const Symbol& sym( cellPair.first );
676 0 0         if( sym.type() == Symbol::Slack || sym.type() == Symbol::Error )
    0          
    0          
677 0           return sym;
678             }
679 0           return Symbol();
680             }
681              
682             /* Compute the row which holds the exit symbol for a pivot.
683              
684             This method will return an iterator to the row in the row map
685             which holds the exit symbol. If no appropriate exit symbol is
686             found, the end() iterator will be returned. This indicates that
687             the objective function is unbounded.
688              
689             */
690 7946           RowMap::iterator getLeavingRow( const Symbol& entering )
691             {
692 7946           double ratio = std::numeric_limits<double>::max();
693 7946           auto end = m_rows.end();
694 7946           auto found = m_rows.end();
695 1789089 100         for( auto it = m_rows.begin(); it != end; ++it )
696             {
697 1781143 100         if( it->first.type() != Symbol::External )
698             {
699 674527 50         double temp = it->second->coefficientFor( entering );
700 674527 100         if( temp < 0.0 )
701             {
702 183815           double temp_ratio = -it->second->constant() / temp;
703 183815 100         if( temp_ratio < ratio )
704             {
705 131700           ratio = temp_ratio;
706 674527           found = it;
707             }
708             }
709             }
710             }
711 7946           return found;
712             }
713              
714             /* Compute the leaving row for a marker variable.
715              
716             This method will return an iterator to the row in the row map
717             which holds the given marker variable. The row will be chosen
718             according to the following precedence:
719              
720             1) The row with a restricted basic varible and a negative coefficient
721             for the marker with the smallest ratio of -constant / coefficient.
722              
723             2) The row with a restricted basic variable and the smallest ratio
724             of constant / coefficient.
725              
726             3) The last unrestricted row which contains the marker.
727              
728             If the marker does not exist in any row, the row map end() iterator
729             will be returned. This indicates an internal solver error since
730             the marker *should* exist somewhere in the tableau.
731              
732             */
733 0           RowMap::iterator getMarkerLeavingRow( const Symbol& marker )
734             {
735 0           const double dmax = std::numeric_limits<double>::max();
736 0           double r1 = dmax;
737 0           double r2 = dmax;
738 0           auto end = m_rows.end();
739 0           auto first = end;
740 0           auto second = end;
741 0           auto third = end;
742 0 0         for( auto it = m_rows.begin(); it != end; ++it )
743             {
744 0 0         double c = it->second->coefficientFor( marker );
745 0 0         if( c == 0.0 )
746 0           continue;
747 0 0         if( it->first.type() == Symbol::External )
748             {
749 0           third = it;
750             }
751 0 0         else if( c < 0.0 )
752             {
753 0           double r = -it->second->constant() / c;
754 0 0         if( r < r1 )
755             {
756 0           r1 = r;
757 0           first = it;
758             }
759             }
760             else
761             {
762 0           double r = it->second->constant() / c;
763 0 0         if( r < r2 )
764             {
765 0           r2 = r;
766 0           second = it;
767             }
768             }
769             }
770 0 0         if( first != end )
771 0           return first;
772 0 0         if( second != end )
773 0           return second;
774 0           return third;
775             }
776              
777             /* Remove the effects of a constraint on the objective function.
778              
779             */
780 0           void removeConstraintEffects( const Constraint& cn, const Tag& tag )
781             {
782 0 0         if( tag.marker.type() == Symbol::Error )
783 0           removeMarkerEffects( tag.marker, cn.strength() );
784 0 0         if( tag.other.type() == Symbol::Error )
785 0           removeMarkerEffects( tag.other, cn.strength() );
786 0           }
787              
788             /* Remove the effects of an error marker on the objective function.
789              
790             */
791 0           void removeMarkerEffects( const Symbol& marker, double strength )
792             {
793 0 0         auto row_it = m_rows.find( marker );
794 0 0         if( row_it != m_rows.end() )
795 0 0         m_objective->insert( *row_it->second, -strength );
796             else
797 0 0         m_objective->insert( marker, -strength );
798 0           }
799              
800             /* Test whether a row is composed of all dummy variables.
801              
802             */
803 415           bool allDummies( const Row& row ) const
804             {
805 478 50         for (const auto &rowPair : row.cells())
806             {
807 478 100         if( rowPair.first.type() != Symbol::Dummy )
808 415           return false;
809             }
810 0           return true;
811             }
812              
813             CnMap m_cns;
814             RowMap m_rows;
815             VarMap m_vars;
816             EditMap m_edits;
817             std::vector<Symbol> m_infeasible_rows;
818             std::unique_ptr<Row> m_objective;
819             std::unique_ptr<Row> m_artificial;
820             Symbol::Id m_id_tick;
821             };
822              
823             } // namespace impl
824              
825             } // namespace kiwi