add Kiwi (Cassowary implementation, as a header-only library)

This commit is contained in:
Paul Davis 2020-05-20 21:32:41 -06:00
parent 1c24e9abef
commit c8f85d6b6c
18 changed files with 3299 additions and 0 deletions

356
libs/kiwi/AssocVector.h Normal file
View file

@ -0,0 +1,356 @@
////////////////////////////////////////////////////////////////////////////////
// The Loki Library
// Copyright (c) 2001 by Andrei Alexandrescu
// This code accompanies the book:
// Alexandrescu, Andrei. "Modern C++ Design: Generic Programming and Design
// Patterns Applied". Copyright (c) 2001. Addison-Wesley.
// Permission to use, copy, modify, distribute and sell this software for any
// purpose is hereby granted without fee, provided that the above copyright
// notice appear in all copies and that both that copyright notice and this
// permission notice appear in supporting documentation.
// The author or Addison-Wesley Longman make no representations about the
// suitability of this software for any purpose. It is provided "as is"
// without express or implied warranty.
////////////////////////////////////////////////////////////////////////////////
// Updated 2019 by Matthieu Dartiailh for C++11 compliancy
////////////////////////////////////////////////////////////////////////////////
#pragma once
// $Id: AssocVector.h 765 2006-10-18 13:55:32Z syntheticpp $
#include <algorithm>
#include <functional>
#include <vector>
#include <utility>
namespace Loki
{
////////////////////////////////////////////////////////////////////////////////
// class template AssocVectorCompare
// Used by AssocVector
////////////////////////////////////////////////////////////////////////////////
namespace Private
{
template <class Value, class C>
class AssocVectorCompare : public C
{
typedef std::pair<typename C::first_argument_type, Value>
Data;
typedef typename C::first_argument_type first_argument_type;
public:
AssocVectorCompare()
{}
AssocVectorCompare(const C& src) : C(src)
{}
bool operator()(const first_argument_type& lhs,
const first_argument_type& rhs) const
{ return C::operator()(lhs, rhs); }
bool operator()(const Data& lhs, const Data& rhs) const
{ return operator()(lhs.first, rhs.first); }
bool operator()(const Data& lhs,
const first_argument_type& rhs) const
{ return operator()(lhs.first, rhs); }
bool operator()(const first_argument_type& lhs,
const Data& rhs) const
{ return operator()(lhs, rhs.first); }
};
}
////////////////////////////////////////////////////////////////////////////////
// class template AssocVector
// An associative vector built as a syntactic drop-in replacement for std::map
// BEWARE: AssocVector doesn't respect all map's guarantees, the most important
// being:
// * iterators are invalidated by insert and erase operations
// * the complexity of insert/erase is O(N) not O(log N)
// * value_type is std::pair<K, V> not std::pair<const K, V>
// * iterators are random
////////////////////////////////////////////////////////////////////////////////
template
<
class K,
class V,
class C = std::less<K>,
class A = std::allocator< std::pair<K, V> >
>
class AssocVector
: private std::vector< std::pair<K, V>, A >
, private Private::AssocVectorCompare<V, C>
{
typedef std::vector<std::pair<K, V>, A> Base;
typedef Private::AssocVectorCompare<V, C> MyCompare;
public:
typedef K key_type;
typedef V mapped_type;
typedef typename Base::value_type value_type;
typedef C key_compare;
typedef A allocator_type;
typedef typename A::reference reference;
typedef typename A::const_reference const_reference;
typedef typename Base::iterator iterator;
typedef typename Base::const_iterator const_iterator;
typedef typename Base::size_type size_type;
typedef typename Base::difference_type difference_type;
typedef typename A::pointer pointer;
typedef typename A::const_pointer const_pointer;
typedef typename Base::reverse_iterator reverse_iterator;
typedef typename Base::const_reverse_iterator const_reverse_iterator;
class value_compare
: public std::function<bool(value_type, value_type)>
, private key_compare
{
friend class AssocVector;
protected:
value_compare(key_compare pred) : key_compare(pred)
{}
public:
bool operator()(const value_type& lhs, const value_type& rhs) const
{ return key_compare::operator()(lhs.first, rhs.first); }
};
// 23.3.1.1 construct/copy/destroy
explicit AssocVector(const key_compare& comp = key_compare(),
const A& alloc = A())
: Base(alloc), MyCompare(comp)
{}
template <class InputIterator>
AssocVector(InputIterator first, InputIterator last,
const key_compare& comp = key_compare(),
const A& alloc = A())
: Base(first, last, alloc), MyCompare(comp)
{
MyCompare& me = *this;
std::sort(begin(), end(), me);
}
AssocVector& operator=(const AssocVector& rhs)
{
AssocVector(rhs).swap(*this);
return *this;
}
// iterators:
// The following are here because MWCW gets 'using' wrong
iterator begin() { return Base::begin(); }
const_iterator begin() const { return Base::begin(); }
iterator end() { return Base::end(); }
const_iterator end() const { return Base::end(); }
reverse_iterator rbegin() { return Base::rbegin(); }
const_reverse_iterator rbegin() const { return Base::rbegin(); }
reverse_iterator rend() { return Base::rend(); }
const_reverse_iterator rend() const { return Base::rend(); }
// capacity:
bool empty() const { return Base::empty(); }
size_type size() const { return Base::size(); }
size_type max_size() { return Base::max_size(); }
// 23.3.1.2 element access:
mapped_type& operator[](const key_type& key)
{ return insert(value_type(key, mapped_type())).first->second; }
// modifiers:
std::pair<iterator, bool> insert(const value_type& val)
{
bool found(true);
iterator i(lower_bound(val.first));
if (i == end() || this->operator()(val.first, i->first))
{
i = Base::insert(i, val);
found = false;
}
return std::make_pair(i, !found);
}
//Section [23.1.2], Table 69
//http://developer.apple.com/documentation/DeveloperTools/gcc-3.3/libstdc++/23_containers/howto.html#4
iterator insert(iterator pos, const value_type& val)
{
if( (pos == begin() || this->operator()(*(pos-1),val)) &&
(pos == end() || this->operator()(val, *pos)) )
{
return Base::insert(pos, val);
}
return insert(val).first;
}
template <class InputIterator>
void insert(InputIterator first, InputIterator last)
{ for (; first != last; ++first) insert(*first); }
void erase(iterator pos)
{ Base::erase(pos); }
size_type erase(const key_type& k)
{
iterator i(find(k));
if (i == end()) return 0;
erase(i);
return 1;
}
void erase(iterator first, iterator last)
{ Base::erase(first, last); }
void swap(AssocVector& other)
{
Base::swap(other);
MyCompare& me = *this;
MyCompare& rhs = other;
std::swap(me, rhs);
}
void clear()
{ Base::clear(); }
// observers:
key_compare key_comp() const
{ return *this; }
value_compare value_comp() const
{
const key_compare& comp = *this;
return value_compare(comp);
}
// 23.3.1.3 map operations:
iterator find(const key_type& k)
{
iterator i(lower_bound(k));
if (i != end() && this->operator()(k, i->first))
{
i = end();
}
return i;
}
const_iterator find(const key_type& k) const
{
const_iterator i(lower_bound(k));
if (i != end() && this->operator()(k, i->first))
{
i = end();
}
return i;
}
size_type count(const key_type& k) const
{ return find(k) != end(); }
iterator lower_bound(const key_type& k)
{
MyCompare& me = *this;
return std::lower_bound(begin(), end(), k, me);
}
const_iterator lower_bound(const key_type& k) const
{
const MyCompare& me = *this;
return std::lower_bound(begin(), end(), k, me);
}
iterator upper_bound(const key_type& k)
{
MyCompare& me = *this;
return std::upper_bound(begin(), end(), k, me);
}
const_iterator upper_bound(const key_type& k) const
{
const MyCompare& me = *this;
return std::upper_bound(begin(), end(), k, me);
}
std::pair<iterator, iterator> equal_range(const key_type& k)
{
MyCompare& me = *this;
return std::equal_range(begin(), end(), k, me);
}
std::pair<const_iterator, const_iterator> equal_range(
const key_type& k) const
{
const MyCompare& me = *this;
return std::equal_range(begin(), end(), k, me);
}
template <class K1, class V1, class C1, class A1>
friend bool operator==(const AssocVector<K1, V1, C1, A1>& lhs,
const AssocVector<K1, V1, C1, A1>& rhs);
bool operator<(const AssocVector& rhs) const
{
const Base& me = *this;
const Base& yo = rhs;
return me < yo;
}
template <class K1, class V1, class C1, class A1>
friend bool operator!=(const AssocVector<K1, V1, C1, A1>& lhs,
const AssocVector<K1, V1, C1, A1>& rhs);
template <class K1, class V1, class C1, class A1>
friend bool operator>(const AssocVector<K1, V1, C1, A1>& lhs,
const AssocVector<K1, V1, C1, A1>& rhs);
template <class K1, class V1, class C1, class A1>
friend bool operator>=(const AssocVector<K1, V1, C1, A1>& lhs,
const AssocVector<K1, V1, C1, A1>& rhs);
template <class K1, class V1, class C1, class A1>
friend bool operator<=(const AssocVector<K1, V1, C1, A1>& lhs,
const AssocVector<K1, V1, C1, A1>& rhs);
};
template <class K, class V, class C, class A>
inline bool operator==(const AssocVector<K, V, C, A>& lhs,
const AssocVector<K, V, C, A>& rhs)
{
const std::vector<std::pair<K, V>, A>& me = lhs;
return me == rhs;
}
template <class K, class V, class C, class A>
inline bool operator!=(const AssocVector<K, V, C, A>& lhs,
const AssocVector<K, V, C, A>& rhs)
{ return !(lhs == rhs); }
template <class K, class V, class C, class A>
inline bool operator>(const AssocVector<K, V, C, A>& lhs,
const AssocVector<K, V, C, A>& rhs)
{ return rhs < lhs; }
template <class K, class V, class C, class A>
inline bool operator>=(const AssocVector<K, V, C, A>& lhs,
const AssocVector<K, V, C, A>& rhs)
{ return !(lhs < rhs); }
template <class K, class V, class C, class A>
inline bool operator<=(const AssocVector<K, V, C, A>& lhs,
const AssocVector<K, V, C, A>& rhs)
{ return !(rhs < lhs); }
// specialized algorithms:
template <class K, class V, class C, class A>
void swap(AssocVector<K, V, C, A>& lhs, AssocVector<K, V, C, A>& rhs)
{ lhs.swap(rhs); }
} // namespace Loki

119
libs/kiwi/constraint.h Normal file
View file

@ -0,0 +1,119 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#include <map>
#include <vector>
#include "expression.h"
#include "shareddata.h"
#include "strength.h"
#include "term.h"
#include "variable.h"
namespace kiwi
{
enum RelationalOperator
{
OP_LE,
OP_GE,
OP_EQ
};
class Constraint
{
public:
Constraint() : m_data(0) {}
Constraint(const Expression &expr,
RelationalOperator op,
double strength = strength::required) : m_data(new ConstraintData(expr, op, strength)) {}
Constraint(const Constraint &other, double strength) : m_data(new ConstraintData(other, strength)) {}
~Constraint() {}
const Expression &expression() const
{
return m_data->m_expression;
}
RelationalOperator op() const
{
return m_data->m_op;
}
double strength() const
{
return m_data->m_strength;
}
bool operator!() const
{
return !m_data;
}
private:
static Expression reduce(const Expression &expr)
{
std::map<Variable, double> vars;
typedef std::vector<Term>::const_iterator iter_t;
iter_t end = expr.terms().end();
for (iter_t it = expr.terms().begin(); it != end; ++it)
vars[it->variable()] += it->coefficient();
std::vector<Term> terms(vars.begin(), vars.end());
return Expression(terms, expr.constant());
}
class ConstraintData : public SharedData
{
public:
ConstraintData(const Expression &expr,
RelationalOperator op,
double strength) : SharedData(),
m_expression(reduce(expr)),
m_strength(strength::clip(strength)),
m_op(op) {}
ConstraintData(const Constraint &other, double strength) : SharedData(),
m_expression(other.expression()),
m_strength(strength::clip(strength)),
m_op(other.op()) {}
~ConstraintData() {}
Expression m_expression;
double m_strength;
RelationalOperator m_op;
private:
ConstraintData(const ConstraintData &other);
ConstraintData &operator=(const ConstraintData &other);
};
SharedDataPtr<ConstraintData> m_data;
friend bool operator<(const Constraint &lhs, const Constraint &rhs)
{
return lhs.m_data < rhs.m_data;
}
friend bool operator==(const Constraint &lhs, const Constraint &rhs)
{
return lhs.m_data == rhs.m_data;
}
friend bool operator!=(const Constraint &lhs, const Constraint &rhs)
{
return lhs.m_data != rhs.m_data;
}
};
} // namespace kiwi

200
libs/kiwi/debug.h Normal file
View file

@ -0,0 +1,200 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#include <iostream>
#include <sstream>
#include <vector>
#include "constraint.h"
#include "solverimpl.h"
#include "term.h"
namespace kiwi
{
namespace impl
{
class DebugHelper
{
public:
static void dump(const SolverImpl &solver, std::ostream &out)
{
out << "Objective" << std::endl;
out << "---------" << std::endl;
dump(*solver.m_objective, out);
out << std::endl;
out << "Tableau" << std::endl;
out << "-------" << std::endl;
dump(solver.m_rows, out);
out << std::endl;
out << "Infeasible" << std::endl;
out << "----------" << std::endl;
dump(solver.m_infeasible_rows, out);
out << std::endl;
out << "Variables" << std::endl;
out << "---------" << std::endl;
dump(solver.m_vars, out);
out << std::endl;
out << "Edit Variables" << std::endl;
out << "--------------" << std::endl;
dump(solver.m_edits, out);
out << std::endl;
out << "Constraints" << std::endl;
out << "-----------" << std::endl;
dump(solver.m_cns, out);
out << std::endl;
out << std::endl;
}
static void dump(const SolverImpl::RowMap &rows, std::ostream &out)
{
typedef SolverImpl::RowMap::const_iterator iter_t;
iter_t end = rows.end();
for (iter_t it = rows.begin(); it != end; ++it)
{
dump(it->first, out);
out << " | ";
dump(*it->second, out);
}
}
static void dump(const std::vector<Symbol> &symbols, std::ostream &out)
{
typedef std::vector<Symbol>::const_iterator iter_t;
iter_t end = symbols.end();
for (iter_t it = symbols.begin(); it != end; ++it)
{
dump(*it, out);
out << std::endl;
}
}
static void dump(const SolverImpl::VarMap &vars, std::ostream &out)
{
typedef SolverImpl::VarMap::const_iterator iter_t;
iter_t end = vars.end();
for (iter_t it = vars.begin(); it != end; ++it)
{
out << it->first.name() << " = ";
dump(it->second, out);
out << std::endl;
}
}
static void dump(const SolverImpl::CnMap &cns, std::ostream &out)
{
typedef SolverImpl::CnMap::const_iterator iter_t;
iter_t end = cns.end();
for (iter_t it = cns.begin(); it != end; ++it)
dump(it->first, out);
}
static void dump(const SolverImpl::EditMap &edits, std::ostream &out)
{
typedef SolverImpl::EditMap::const_iterator iter_t;
iter_t end = edits.end();
for (iter_t it = edits.begin(); it != end; ++it)
out << it->first.name() << std::endl;
}
static void dump(const Row &row, std::ostream &out)
{
typedef Row::CellMap::const_iterator iter_t;
out << row.constant();
iter_t end = row.cells().end();
for (iter_t it = row.cells().begin(); it != end; ++it)
{
out << " + " << it->second << " * ";
dump(it->first, out);
}
out << std::endl;
}
static void dump(const Symbol &symbol, std::ostream &out)
{
switch (symbol.type())
{
case Symbol::Invalid:
out << "i";
break;
case Symbol::External:
out << "v";
break;
case Symbol::Slack:
out << "s";
break;
case Symbol::Error:
out << "e";
break;
case Symbol::Dummy:
out << "d";
break;
default:
break;
}
out << symbol.id();
}
static void dump(const Constraint &cn, std::ostream &out)
{
typedef std::vector<Term>::const_iterator iter_t;
iter_t begin = cn.expression().terms().begin();
iter_t end = cn.expression().terms().end();
for (iter_t it = begin; it != end; ++it)
{
out << it->coefficient() << " * ";
out << it->variable().name() << " + ";
}
out << cn.expression().constant();
switch (cn.op())
{
case OP_LE:
out << " <= 0 ";
break;
case OP_GE:
out << " >= 0 ";
break;
case OP_EQ:
out << " == 0 ";
break;
default:
break;
}
out << " | strength = " << cn.strength() << std::endl;
}
};
} // namespace impl
namespace debug
{
template <typename T>
void dump(const T &value)
{
impl::DebugHelper::dump(value, std::cout);
}
template <typename T>
void dump(const T &value, std::ostream &out)
{
impl::DebugHelper::dump(value, out);
}
template <typename T>
std::string dumps(const T &value)
{
std::stringstream stream;
impl::DebugHelper::dump(value, stream);
return stream.str();
}
} // namespace debug
} // namespace kiwi

162
libs/kiwi/errors.h Normal file
View file

@ -0,0 +1,162 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#include <exception>
#include <string>
#include "constraint.h"
#include "variable.h"
namespace kiwi
{
class UnsatisfiableConstraint : public std::exception
{
public:
UnsatisfiableConstraint(const Constraint &constraint) : m_constraint(constraint) {}
~UnsatisfiableConstraint() throw() {}
const char *what() const throw()
{
return "The constraint can not be satisfied.";
}
const Constraint &constraint() const
{
return m_constraint;
}
private:
Constraint m_constraint;
};
class UnknownConstraint : public std::exception
{
public:
UnknownConstraint(const Constraint &constraint) : m_constraint(constraint) {}
~UnknownConstraint() throw() {}
const char *what() const throw()
{
return "The constraint has not been added to the solver.";
}
const Constraint &constraint() const
{
return m_constraint;
}
private:
Constraint m_constraint;
};
class DuplicateConstraint : public std::exception
{
public:
DuplicateConstraint(const Constraint &constraint) : m_constraint(constraint) {}
~DuplicateConstraint() throw() {}
const char *what() const throw()
{
return "The constraint has already been added to the solver.";
}
const Constraint &constraint() const
{
return m_constraint;
}
private:
Constraint m_constraint;
};
class UnknownEditVariable : public std::exception
{
public:
UnknownEditVariable(const Variable &variable) : m_variable(variable) {}
~UnknownEditVariable() throw() {}
const char *what() const throw()
{
return "The edit variable has not been added to the solver.";
}
const Variable &variable() const
{
return m_variable;
}
private:
Variable m_variable;
};
class DuplicateEditVariable : public std::exception
{
public:
DuplicateEditVariable(const Variable &variable) : m_variable(variable) {}
~DuplicateEditVariable() throw() {}
const char *what() const throw()
{
return "The edit variable has already been added to the solver.";
}
const Variable &variable() const
{
return m_variable;
}
private:
Variable m_variable;
};
class BadRequiredStrength : public std::exception
{
public:
BadRequiredStrength() {}
~BadRequiredStrength() throw() {}
const char *what() const throw()
{
return "A required strength cannot be used in this context.";
}
};
class InternalSolverError : public std::exception
{
public:
InternalSolverError() : m_msg("An internal solver error ocurred.") {}
InternalSolverError(const char *msg) : m_msg(msg) {}
InternalSolverError(const std::string &msg) : m_msg(msg) {}
~InternalSolverError() throw() {}
const char *what() const throw()
{
return m_msg.c_str();
}
private:
std::string m_msg;
};
} // namespace kiwi

52
libs/kiwi/expression.h Normal file
View file

@ -0,0 +1,52 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#include <vector>
#include "term.h"
namespace kiwi
{
class Expression
{
public:
Expression(double constant = 0.0) : m_constant(constant) {}
Expression(const Term &term, double constant = 0.0) : m_terms(1, term), m_constant(constant) {}
Expression(const std::vector<Term> &terms, double constant = 0.0) : m_terms(terms), m_constant(constant) {}
~Expression() {}
const std::vector<Term> &terms() const
{
return m_terms;
}
double constant() const
{
return m_constant;
}
double value() const
{
typedef std::vector<Term>::const_iterator iter_t;
double result = m_constant;
iter_t end = m_terms.end();
for (iter_t it = m_terms.begin(); it != end; ++it)
result += it->value();
return result;
}
private:
std::vector<Term> m_terms;
double m_constant;
};
} // namespace kiwi

19
libs/kiwi/kiwi.h Normal file
View file

@ -0,0 +1,19 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#include "constraint.h"
#include "debug.h"
#include "errors.h"
#include "expression.h"
#include "shareddata.h"
#include "solver.h"
#include "strength.h"
#include "symbolics.h"
#include "term.h"
#include "variable.h"
#include "version.h"

37
libs/kiwi/maptype.h Normal file
View file

@ -0,0 +1,37 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2019, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#include <functional>
#include <map>
#include <memory>
#include <utility>
#include "AssocVector.h"
namespace kiwi
{
namespace impl
{
template <
typename K,
typename V,
typename C = std::less<K>,
typename A = std::allocator<std::pair<K, V>>>
using MapType = Loki::AssocVector<K, V, C, A>;
// template<
// typename K,
// typename V,
// typename C = std::less<K>,
// typename A = std::allocator< std::pair<const K, V> > >
// using MapType = std::map<K, V, C, A>;
} // namespace impl
} // namespace kiwi

188
libs/kiwi/row.h Normal file
View file

@ -0,0 +1,188 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#include "maptype.h"
#include "symbol.h"
#include "util.h"
namespace kiwi
{
namespace impl
{
class Row
{
public:
typedef MapType<Symbol, double> CellMap;
Row() : m_constant(0.0) {}
Row(double constant) : m_constant(constant) {}
Row(const Row &other) : m_cells(other.m_cells), m_constant(other.m_constant) {}
~Row() {}
const CellMap &cells() const
{
return m_cells;
}
double constant() const
{
return m_constant;
}
/* Add a constant value to the row constant.
The new value of the constant is returned.
*/
double add(double value)
{
return m_constant += value;
}
/* Insert a symbol into the row with a given coefficient.
If the symbol already exists in the row, the coefficient will be
added to the existing coefficient. If the resulting coefficient
is zero, the symbol will be removed from the row.
*/
void insert(const Symbol &symbol, double coefficient = 1.0)
{
if (nearZero(m_cells[symbol] += coefficient))
m_cells.erase(symbol);
}
/* Insert a row into this row with a given coefficient.
The constant and the cells of the other row will be multiplied by
the coefficient and added to this row. Any cell with a resulting
coefficient of zero will be removed from the row.
*/
void insert(const Row &other, double coefficient = 1.0)
{
typedef CellMap::const_iterator iter_t;
m_constant += other.m_constant * coefficient;
iter_t end = other.m_cells.end();
for (iter_t it = other.m_cells.begin(); it != end; ++it)
{
double coeff = it->second * coefficient;
if (nearZero(m_cells[it->first] += coeff))
m_cells.erase(it->first);
}
}
/* Remove the given symbol from the row.
*/
void remove(const Symbol &symbol)
{
CellMap::iterator it = m_cells.find(symbol);
if (it != m_cells.end())
m_cells.erase(it);
}
/* Reverse the sign of the constant and all cells in the row.
*/
void reverseSign()
{
typedef CellMap::iterator iter_t;
m_constant = -m_constant;
iter_t end = m_cells.end();
for (iter_t it = m_cells.begin(); it != end; ++it)
it->second = -it->second;
}
/* Solve the row for the given symbol.
This method assumes the row is of the form a * x + b * y + c = 0
and (assuming solve for x) will modify the row to represent the
right hand side of x = -b/a * y - c / a. The target symbol will
be removed from the row, and the constant and other cells will
be multiplied by the negative inverse of the target coefficient.
The given symbol *must* exist in the row.
*/
void solveFor(const Symbol &symbol)
{
typedef CellMap::iterator iter_t;
double coeff = -1.0 / m_cells[symbol];
m_cells.erase(symbol);
m_constant *= coeff;
iter_t end = m_cells.end();
for (iter_t it = m_cells.begin(); it != end; ++it)
it->second *= coeff;
}
/* Solve the row for the given symbols.
This method assumes the row is of the form x = b * y + c and will
solve the row such that y = x / b - c / b. The rhs symbol will be
removed from the row, the lhs added, and the result divided by the
negative inverse of the rhs coefficient.
The lhs symbol *must not* exist in the row, and the rhs symbol
*must* exist in the row.
*/
void solveFor(const Symbol &lhs, const Symbol &rhs)
{
insert(lhs, -1.0);
solveFor(rhs);
}
/* Get the coefficient for the given symbol.
If the symbol does not exist in the row, zero will be returned.
*/
double coefficientFor(const Symbol &symbol) const
{
CellMap::const_iterator it = m_cells.find(symbol);
if (it == m_cells.end())
return 0.0;
return it->second;
}
/* Substitute a symbol with the data from another row.
Given a row of the form a * x + b and a substitution of the
form x = 3 * y + c the row will be updated to reflect the
expression 3 * a * y + a * c + b.
If the symbol does not exist in the row, this is a no-op.
*/
void substitute(const Symbol &symbol, const Row &row)
{
typedef CellMap::iterator iter_t;
iter_t it = m_cells.find(symbol);
if (it != m_cells.end())
{
double coefficient = it->second;
m_cells.erase(it);
insert(row, coefficient);
}
}
private:
CellMap m_cells;
double m_constant;
};
} // namespace impl
} // namespace kiwi

151
libs/kiwi/shareddata.h Normal file
View file

@ -0,0 +1,151 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
namespace kiwi
{
class SharedData
{
public:
SharedData() : m_refcount(0) {}
SharedData(const SharedData &other) : m_refcount(0) {}
int m_refcount;
private:
SharedData &operator=(const SharedData &other);
};
template <typename T>
class SharedDataPtr
{
public:
typedef T Type;
SharedDataPtr() : m_data(0) {}
explicit SharedDataPtr(T *data) : m_data(data)
{
incref(m_data);
}
~SharedDataPtr()
{
decref(m_data);
}
T *data()
{
return m_data;
}
const T *data() const
{
return m_data;
}
operator T *()
{
return m_data;
}
operator const T *() const
{
return m_data;
}
T *operator->()
{
return m_data;
}
const T *operator->() const
{
return m_data;
}
T &operator*()
{
return *m_data;
}
const T &operator*() const
{
return *m_data;
}
bool operator!() const
{
return !m_data;
}
bool operator<(const SharedDataPtr<T> &other) const
{
return m_data < other.m_data;
}
bool operator==(const SharedDataPtr<T> &other) const
{
return m_data == other.m_data;
}
bool operator!=(const SharedDataPtr<T> &other) const
{
return m_data != other.m_data;
}
SharedDataPtr(const SharedDataPtr<T> &other) : m_data(other.m_data)
{
incref(m_data);
}
SharedDataPtr<T> &operator=(const SharedDataPtr<T> &other)
{
if (m_data != other.m_data)
{
T *temp = m_data;
m_data = other.m_data;
incref(m_data);
decref(temp);
}
return *this;
}
SharedDataPtr<T> &operator=(T *other)
{
if (m_data != other)
{
T *temp = m_data;
m_data = other;
incref(m_data);
decref(temp);
}
return *this;
}
private:
static void incref(T *data)
{
if (data)
++data->m_refcount;
}
static void decref(T *data)
{
if (data && --data->m_refcount == 0)
delete data;
}
T *m_data;
};
} // namespace kiwi

178
libs/kiwi/solver.h Normal file
View file

@ -0,0 +1,178 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#include "constraint.h"
#include "debug.h"
#include "solverimpl.h"
#include "strength.h"
#include "variable.h"
namespace kiwi
{
class Solver
{
public:
Solver() {}
~Solver() {}
/* Add a constraint to the solver.
Throws
------
DuplicateConstraint
The given constraint has already been added to the solver.
UnsatisfiableConstraint
The given constraint is required and cannot be satisfied.
*/
void addConstraint( const Constraint& constraint )
{
m_impl.addConstraint( constraint );
}
/* Remove a constraint from the solver.
Throws
------
UnknownConstraint
The given constraint has not been added to the solver.
*/
void removeConstraint( const Constraint& constraint )
{
m_impl.removeConstraint( constraint );
}
/* Test whether a constraint has been added to the solver.
*/
bool hasConstraint( const Constraint& constraint ) const
{
return m_impl.hasConstraint( constraint );
}
/* Add an edit variable to the solver.
This method should be called before the `suggestValue` method is
used to supply a suggested value for the given edit variable.
Throws
------
DuplicateEditVariable
The given edit variable has already been added to the solver.
BadRequiredStrength
The given strength is >= required.
*/
void addEditVariable( const Variable& variable, double strength )
{
m_impl.addEditVariable( variable, strength );
}
/* Remove an edit variable from the solver.
Throws
------
UnknownEditVariable
The given edit variable has not been added to the solver.
*/
void removeEditVariable( const Variable& variable )
{
m_impl.removeEditVariable( variable );
}
/* Test whether an edit variable has been added to the solver.
*/
bool hasEditVariable( const Variable& variable ) const
{
return m_impl.hasEditVariable( variable );
}
/* Suggest a value for the given edit variable.
This method should be used after an edit variable as been added to
the solver in order to suggest the value for that variable. After
all suggestions have been made, the `solve` method can be used to
update the values of all variables.
Throws
------
UnknownEditVariable
The given edit variable has not been added to the solver.
*/
void suggestValue( const Variable& variable, double value )
{
m_impl.suggestValue( variable, value );
}
/* Update the values of the external solver variables.
*/
void updateVariables()
{
m_impl.updateVariables();
}
/* Reset the solver to the empty starting condition.
This method resets the internal solver state to the empty starting
condition, as if no constraints or edit variables have been added.
This can be faster than deleting the solver and creating a new one
when the entire system must change, since it can avoid unecessary
heap (de)allocations.
*/
void reset()
{
m_impl.reset();
}
/* Dump a representation of the solver internals to stdout.
*/
void dump()
{
debug::dump( m_impl );
}
/* Dump a representation of the solver internals to a stream.
*/
void dump( std::ostream& out )
{
debug::dump( m_impl, out );
}
/* Dump a representation of the solver internals to a string.
*/
std::string dumps()
{
return debug::dumps( m_impl );
}
private:
Solver( const Solver& );
Solver& operator=( const Solver& );
impl::SolverImpl m_impl;
};
} // namespace kiwi

840
libs/kiwi/solverimpl.h Normal file
View file

@ -0,0 +1,840 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#include <algorithm>
#include <limits>
#include <memory>
#include <vector>
#include "constraint.h"
#include "errors.h"
#include "expression.h"
#include "maptype.h"
#include "row.h"
#include "symbol.h"
#include "term.h"
#include "util.h"
#include "variable.h"
namespace kiwi
{
namespace impl
{
class SolverImpl
{
friend class DebugHelper;
struct Tag
{
Symbol marker;
Symbol other;
};
struct EditInfo
{
Tag tag;
Constraint constraint;
double constant;
};
typedef MapType<Variable, Symbol> VarMap;
typedef MapType<Symbol, Row*> RowMap;
typedef MapType<Constraint, Tag> CnMap;
typedef MapType<Variable, EditInfo> EditMap;
struct DualOptimizeGuard
{
DualOptimizeGuard( SolverImpl& impl ) : m_impl( impl ) {}
~DualOptimizeGuard() { m_impl.dualOptimize(); }
SolverImpl& m_impl;
};
public:
SolverImpl() : m_objective( new Row() ), m_id_tick( 1 ) {}
~SolverImpl() { clearRows(); }
/* Add a constraint to the solver.
Throws
------
DuplicateConstraint
The given constraint has already been added to the solver.
UnsatisfiableConstraint
The given constraint is required and cannot be satisfied.
*/
void addConstraint( const Constraint& constraint )
{
if( m_cns.find( constraint ) != m_cns.end() )
throw DuplicateConstraint( constraint );
// Creating a row causes symbols to be reserved for the variables
// in the constraint. If this method exits with an exception,
// then its possible those variables will linger in the var map.
// Since its likely that those variables will be used in other
// constraints and since exceptional conditions are uncommon,
// i'm not too worried about aggressive cleanup of the var map.
Tag tag;
std::unique_ptr<Row> rowptr( createRow( constraint, tag ) );
Symbol subject( chooseSubject( *rowptr, tag ) );
// If chooseSubject could not find a valid entering symbol, one
// last option is available if the entire row is composed of
// dummy variables. If the constant of the row is zero, then
// this represents redundant constraints and the new dummy
// marker can enter the basis. If the constant is non-zero,
// then it represents an unsatisfiable constraint.
if( subject.type() == Symbol::Invalid && allDummies( *rowptr ) )
{
if( !nearZero( rowptr->constant() ) )
throw UnsatisfiableConstraint( constraint );
else
subject = tag.marker;
}
// If an entering symbol still isn't found, then the row must
// be added using an artificial variable. If that fails, then
// the row represents an unsatisfiable constraint.
if( subject.type() == Symbol::Invalid )
{
if( !addWithArtificialVariable( *rowptr ) )
throw UnsatisfiableConstraint( constraint );
}
else
{
rowptr->solveFor( subject );
substitute( subject, *rowptr );
m_rows[ subject ] = rowptr.release();
}
m_cns[ constraint ] = tag;
// Optimizing after each constraint is added performs less
// aggregate work due to a smaller average system size. It
// also ensures the solver remains in a consistent state.
optimize( *m_objective );
}
/* Remove a constraint from the solver.
Throws
------
UnknownConstraint
The given constraint has not been added to the solver.
*/
void removeConstraint( const Constraint& constraint )
{
CnMap::iterator cn_it = m_cns.find( constraint );
if( cn_it == m_cns.end() )
throw UnknownConstraint( constraint );
Tag tag( cn_it->second );
m_cns.erase( cn_it );
// Remove the error effects from the objective function
// *before* pivoting, or substitutions into the objective
// will lead to incorrect solver results.
removeConstraintEffects( constraint, tag );
// If the marker is basic, simply drop the row. Otherwise,
// pivot the marker into the basis and then drop the row.
RowMap::iterator row_it = m_rows.find( tag.marker );
if( row_it != m_rows.end() )
{
std::unique_ptr<Row> rowptr( row_it->second );
m_rows.erase( row_it );
}
else
{
row_it = getMarkerLeavingRow( tag.marker );
if( row_it == m_rows.end() )
throw InternalSolverError( "failed to find leaving row" );
Symbol leaving( row_it->first );
std::unique_ptr<Row> rowptr( row_it->second );
m_rows.erase( row_it );
rowptr->solveFor( leaving, tag.marker );
substitute( tag.marker, *rowptr );
}
// Optimizing after each constraint is removed ensures that the
// solver remains consistent. It makes the solver api easier to
// use at a small tradeoff for speed.
optimize( *m_objective );
}
/* Test whether a constraint has been added to the solver.
*/
bool hasConstraint( const Constraint& constraint ) const
{
return m_cns.find( constraint ) != m_cns.end();
}
/* Add an edit variable to the solver.
This method should be called before the `suggestValue` method is
used to supply a suggested value for the given edit variable.
Throws
------
DuplicateEditVariable
The given edit variable has already been added to the solver.
BadRequiredStrength
The given strength is >= required.
*/
void addEditVariable( const Variable& variable, double strength )
{
if( m_edits.find( variable ) != m_edits.end() )
throw DuplicateEditVariable( variable );
strength = strength::clip( strength );
if( strength == strength::required )
throw BadRequiredStrength();
Constraint cn( Expression( variable ), OP_EQ, strength );
addConstraint( cn );
EditInfo info;
info.tag = m_cns[ cn ];
info.constraint = cn;
info.constant = 0.0;
m_edits[ variable ] = info;
}
/* Remove an edit variable from the solver.
Throws
------
UnknownEditVariable
The given edit variable has not been added to the solver.
*/
void removeEditVariable( const Variable& variable )
{
EditMap::iterator it = m_edits.find( variable );
if( it == m_edits.end() )
throw UnknownEditVariable( variable );
removeConstraint( it->second.constraint );
m_edits.erase( it );
}
/* Test whether an edit variable has been added to the solver.
*/
bool hasEditVariable( const Variable& variable ) const
{
return m_edits.find( variable ) != m_edits.end();
}
/* Suggest a value for the given edit variable.
This method should be used after an edit variable as been added to
the solver in order to suggest the value for that variable.
Throws
------
UnknownEditVariable
The given edit variable has not been added to the solver.
*/
void suggestValue( const Variable& variable, double value )
{
EditMap::iterator it = m_edits.find( variable );
if( it == m_edits.end() )
throw UnknownEditVariable( variable );
DualOptimizeGuard guard( *this );
EditInfo& info = it->second;
double delta = value - info.constant;
info.constant = value;
// Check first if the positive error variable is basic.
RowMap::iterator row_it = m_rows.find( info.tag.marker );
if( row_it != m_rows.end() )
{
if( row_it->second->add( -delta ) < 0.0 )
m_infeasible_rows.push_back( row_it->first );
return;
}
// Check next if the negative error variable is basic.
row_it = m_rows.find( info.tag.other );
if( row_it != m_rows.end() )
{
if( row_it->second->add( delta ) < 0.0 )
m_infeasible_rows.push_back( row_it->first );
return;
}
// Otherwise update each row where the error variables exist.
RowMap::iterator end = m_rows.end();
for( row_it = m_rows.begin(); row_it != end; ++row_it )
{
double coeff = row_it->second->coefficientFor( info.tag.marker );
if( coeff != 0.0 &&
row_it->second->add( delta * coeff ) < 0.0 &&
row_it->first.type() != Symbol::External )
m_infeasible_rows.push_back( row_it->first );
}
}
/* Update the values of the external solver variables.
*/
void updateVariables()
{
typedef RowMap::iterator row_iter_t;
typedef VarMap::iterator var_iter_t;
row_iter_t row_end = m_rows.end();
var_iter_t var_end = m_vars.end();
for( var_iter_t var_it = m_vars.begin(); var_it != var_end; ++var_it )
{
Variable& var( const_cast<Variable&>( var_it->first ) );
row_iter_t row_it = m_rows.find( var_it->second );
if( row_it == row_end )
var.setValue( 0.0 );
else
var.setValue( row_it->second->constant() );
}
}
/* Reset the solver to the empty starting condition.
This method resets the internal solver state to the empty starting
condition, as if no constraints or edit variables have been added.
This can be faster than deleting the solver and creating a new one
when the entire system must change, since it can avoid unecessary
heap (de)allocations.
*/
void reset()
{
clearRows();
m_cns.clear();
m_vars.clear();
m_edits.clear();
m_infeasible_rows.clear();
m_objective.reset( new Row() );
m_artificial.reset();
m_id_tick = 1;
}
private:
SolverImpl( const SolverImpl& );
SolverImpl& operator=( const SolverImpl& );
struct RowDeleter
{
template<typename T>
void operator()( T& pair ) { delete pair.second; }
};
void clearRows()
{
std::for_each( m_rows.begin(), m_rows.end(), RowDeleter() );
m_rows.clear();
}
/* Get the symbol for the given variable.
If a symbol does not exist for the variable, one will be created.
*/
Symbol getVarSymbol( const Variable& variable )
{
VarMap::iterator it = m_vars.find( variable );
if( it != m_vars.end() )
return it->second;
Symbol symbol( Symbol::External, m_id_tick++ );
m_vars[ variable ] = symbol;
return symbol;
}
/* Create a new Row object for the given constraint.
The terms in the constraint will be converted to cells in the row.
Any term in the constraint with a coefficient of zero is ignored.
This method uses the `getVarSymbol` method to get the symbol for
the variables added to the row. If the symbol for a given cell
variable is basic, the cell variable will be substituted with the
basic row.
The necessary slack and error variables will be added to the row.
If the constant for the row is negative, the sign for the row
will be inverted so the constant becomes positive.
The tag will be updated with the marker and error symbols to use
for tracking the movement of the constraint in the tableau.
*/
Row* createRow( const Constraint& constraint, Tag& tag )
{
typedef std::vector<Term>::const_iterator iter_t;
const Expression& expr( constraint.expression() );
Row* row = new Row( expr.constant() );
// Substitute the current basic variables into the row.
iter_t end = expr.terms().end();
for( iter_t it = expr.terms().begin(); it != end; ++it )
{
if( !nearZero( it->coefficient() ) )
{
Symbol symbol( getVarSymbol( it->variable() ) );
RowMap::const_iterator row_it = m_rows.find( symbol );
if( row_it != m_rows.end() )
row->insert( *row_it->second, it->coefficient() );
else
row->insert( symbol, it->coefficient() );
}
}
// Add the necessary slack, error, and dummy variables.
switch( constraint.op() )
{
case OP_LE:
case OP_GE:
{
double coeff = constraint.op() == OP_LE ? 1.0 : -1.0;
Symbol slack( Symbol::Slack, m_id_tick++ );
tag.marker = slack;
row->insert( slack, coeff );
if( constraint.strength() < strength::required )
{
Symbol error( Symbol::Error, m_id_tick++ );
tag.other = error;
row->insert( error, -coeff );
m_objective->insert( error, constraint.strength() );
}
break;
}
case OP_EQ:
{
if( constraint.strength() < strength::required )
{
Symbol errplus( Symbol::Error, m_id_tick++ );
Symbol errminus( Symbol::Error, m_id_tick++ );
tag.marker = errplus;
tag.other = errminus;
row->insert( errplus, -1.0 ); // v = eplus - eminus
row->insert( errminus, 1.0 ); // v - eplus + eminus = 0
m_objective->insert( errplus, constraint.strength() );
m_objective->insert( errminus, constraint.strength() );
}
else
{
Symbol dummy( Symbol::Dummy, m_id_tick++ );
tag.marker = dummy;
row->insert( dummy );
}
break;
}
}
// Ensure the row as a positive constant.
if( row->constant() < 0.0 )
row->reverseSign();
return row;
}
/* Choose the subject for solving for the row.
This method will choose the best subject for using as the solve
target for the row. An invalid symbol will be returned if there
is no valid target.
The symbols are chosen according to the following precedence:
1) The first symbol representing an external variable.
2) A negative slack or error tag variable.
If a subject cannot be found, an invalid symbol will be returned.
*/
Symbol chooseSubject( const Row& row, const Tag& tag )
{
typedef Row::CellMap::const_iterator iter_t;
iter_t end = row.cells().end();
for( iter_t it = row.cells().begin(); it != end; ++it )
{
if( it->first.type() == Symbol::External )
return it->first;
}
if( tag.marker.type() == Symbol::Slack || tag.marker.type() == Symbol::Error )
{
if( row.coefficientFor( tag.marker ) < 0.0 )
return tag.marker;
}
if( tag.other.type() == Symbol::Slack || tag.other.type() == Symbol::Error )
{
if( row.coefficientFor( tag.other ) < 0.0 )
return tag.other;
}
return Symbol();
}
/* Add the row to the tableau using an artificial variable.
This will return false if the constraint cannot be satisfied.
*/
bool addWithArtificialVariable( const Row& row )
{
// Create and add the artificial variable to the tableau
Symbol art( Symbol::Slack, m_id_tick++ );
m_rows[ art ] = new Row( row );
m_artificial.reset( new Row( row ) );
// Optimize the artificial objective. This is successful
// only if the artificial objective is optimized to zero.
optimize( *m_artificial );
bool success = nearZero( m_artificial->constant() );
m_artificial.reset();
// If the artificial variable is not basic, pivot the row so that
// it becomes basic. If the row is constant, exit early.
RowMap::iterator it = m_rows.find( art );
if( it != m_rows.end() )
{
std::unique_ptr<Row> rowptr( it->second );
m_rows.erase( it );
if( rowptr->cells().empty() )
return success;
Symbol entering( anyPivotableSymbol( *rowptr ) );
if( entering.type() == Symbol::Invalid )
return false; // unsatisfiable (will this ever happen?)
rowptr->solveFor( art, entering );
substitute( entering, *rowptr );
m_rows[ entering ] = rowptr.release();
}
// Remove the artificial variable from the tableau.
RowMap::iterator end = m_rows.end();
for( it = m_rows.begin(); it != end; ++it )
it->second->remove( art );
m_objective->remove( art );
return success;
}
/* Substitute the parametric symbol with the given row.
This method will substitute all instances of the parametric symbol
in the tableau and the objective function with the given row.
*/
void substitute( const Symbol& symbol, const Row& row )
{
typedef RowMap::iterator iter_t;
iter_t end = m_rows.end();
for( iter_t it = m_rows.begin(); it != end; ++it )
{
it->second->substitute( symbol, row );
if( it->first.type() != Symbol::External &&
it->second->constant() < 0.0 )
m_infeasible_rows.push_back( it->first );
}
m_objective->substitute( symbol, row );
if( m_artificial.get() )
m_artificial->substitute( symbol, row );
}
/* Optimize the system for the given objective function.
This method performs iterations of Phase 2 of the simplex method
until the objective function reaches a minimum.
Throws
------
InternalSolverError
The value of the objective function is unbounded.
*/
void optimize( const Row& objective )
{
while( true )
{
Symbol entering( getEnteringSymbol( objective ) );
if( entering.type() == Symbol::Invalid )
return;
RowMap::iterator it = getLeavingRow( entering );
if( it == m_rows.end() )
throw InternalSolverError( "The objective is unbounded." );
// pivot the entering symbol into the basis
Symbol leaving( it->first );
Row* row = it->second;
m_rows.erase( it );
row->solveFor( leaving, entering );
substitute( entering, *row );
m_rows[ entering ] = row;
}
}
/* Optimize the system using the dual of the simplex method.
The current state of the system should be such that the objective
function is optimal, but not feasible. This method will perform
an iteration of the dual simplex method to make the solution both
optimal and feasible.
Throws
------
InternalSolverError
The system cannot be dual optimized.
*/
void dualOptimize()
{
while( !m_infeasible_rows.empty() )
{
Symbol leaving( m_infeasible_rows.back() );
m_infeasible_rows.pop_back();
RowMap::iterator it = m_rows.find( leaving );
if( it != m_rows.end() && !nearZero( it->second->constant() ) &&
it->second->constant() < 0.0 )
{
Symbol entering( getDualEnteringSymbol( *it->second ) );
if( entering.type() == Symbol::Invalid )
throw InternalSolverError( "Dual optimize failed." );
// pivot the entering symbol into the basis
Row* row = it->second;
m_rows.erase( it );
row->solveFor( leaving, entering );
substitute( entering, *row );
m_rows[ entering ] = row;
}
}
}
/* Compute the entering variable for a pivot operation.
This method will return first symbol in the objective function which
is non-dummy and has a coefficient less than zero. If no symbol meets
the criteria, it means the objective function is at a minimum, and an
invalid symbol is returned.
*/
Symbol getEnteringSymbol( const Row& objective )
{
typedef Row::CellMap::const_iterator iter_t;
iter_t end = objective.cells().end();
for( iter_t it = objective.cells().begin(); it != end; ++it )
{
if( it->first.type() != Symbol::Dummy && it->second < 0.0 )
return it->first;
}
return Symbol();
}
/* Compute the entering symbol for the dual optimize operation.
This method will return the symbol in the row which has a positive
coefficient and yields the minimum ratio for its respective symbol
in the objective function. The provided row *must* be infeasible.
If no symbol is found which meats the criteria, an invalid symbol
is returned.
*/
Symbol getDualEnteringSymbol( const Row& row )
{
typedef Row::CellMap::const_iterator iter_t;
Symbol entering;
double ratio = std::numeric_limits<double>::max();
iter_t end = row.cells().end();
for( iter_t it = row.cells().begin(); it != end; ++it )
{
if( it->second > 0.0 && it->first.type() != Symbol::Dummy )
{
double coeff = m_objective->coefficientFor( it->first );
double r = coeff / it->second;
if( r < ratio )
{
ratio = r;
entering = it->first;
}
}
}
return entering;
}
/* Get the first Slack or Error symbol in the row.
If no such symbol is present, and Invalid symbol will be returned.
*/
Symbol anyPivotableSymbol( const Row& row )
{
typedef Row::CellMap::const_iterator iter_t;
iter_t end = row.cells().end();
for( iter_t it = row.cells().begin(); it != end; ++it )
{
const Symbol& sym( it->first );
if( sym.type() == Symbol::Slack || sym.type() == Symbol::Error )
return sym;
}
return Symbol();
}
/* Compute the row which holds the exit symbol for a pivot.
This method will return an iterator to the row in the row map
which holds the exit symbol. If no appropriate exit symbol is
found, the end() iterator will be returned. This indicates that
the objective function is unbounded.
*/
RowMap::iterator getLeavingRow( const Symbol& entering )
{
typedef RowMap::iterator iter_t;
double ratio = std::numeric_limits<double>::max();
iter_t end = m_rows.end();
iter_t found = m_rows.end();
for( iter_t it = m_rows.begin(); it != end; ++it )
{
if( it->first.type() != Symbol::External )
{
double temp = it->second->coefficientFor( entering );
if( temp < 0.0 )
{
double temp_ratio = -it->second->constant() / temp;
if( temp_ratio < ratio )
{
ratio = temp_ratio;
found = it;
}
}
}
}
return found;
}
/* Compute the leaving row for a marker variable.
This method will return an iterator to the row in the row map
which holds the given marker variable. The row will be chosen
according to the following precedence:
1) The row with a restricted basic varible and a negative coefficient
for the marker with the smallest ratio of -constant / coefficient.
2) The row with a restricted basic variable and the smallest ratio
of constant / coefficient.
3) The last unrestricted row which contains the marker.
If the marker does not exist in any row, the row map end() iterator
will be returned. This indicates an internal solver error since
the marker *should* exist somewhere in the tableau.
*/
RowMap::iterator getMarkerLeavingRow( const Symbol& marker )
{
const double dmax = std::numeric_limits<double>::max();
typedef RowMap::iterator iter_t;
double r1 = dmax;
double r2 = dmax;
iter_t end = m_rows.end();
iter_t first = end;
iter_t second = end;
iter_t third = end;
for( iter_t it = m_rows.begin(); it != end; ++it )
{
double c = it->second->coefficientFor( marker );
if( c == 0.0 )
continue;
if( it->first.type() == Symbol::External )
{
third = it;
}
else if( c < 0.0 )
{
double r = -it->second->constant() / c;
if( r < r1 )
{
r1 = r;
first = it;
}
}
else
{
double r = it->second->constant() / c;
if( r < r2 )
{
r2 = r;
second = it;
}
}
}
if( first != end )
return first;
if( second != end )
return second;
return third;
}
/* Remove the effects of a constraint on the objective function.
*/
void removeConstraintEffects( const Constraint& cn, const Tag& tag )
{
if( tag.marker.type() == Symbol::Error )
removeMarkerEffects( tag.marker, cn.strength() );
if( tag.other.type() == Symbol::Error )
removeMarkerEffects( tag.other, cn.strength() );
}
/* Remove the effects of an error marker on the objective function.
*/
void removeMarkerEffects( const Symbol& marker, double strength )
{
RowMap::iterator row_it = m_rows.find( marker );
if( row_it != m_rows.end() )
m_objective->insert( *row_it->second, -strength );
else
m_objective->insert( marker, -strength );
}
/* Test whether a row is composed of all dummy variables.
*/
bool allDummies( const Row& row )
{
typedef Row::CellMap::const_iterator iter_t;
iter_t end = row.cells().end();
for( iter_t it = row.cells().begin(); it != end; ++it )
{
if( it->first.type() != Symbol::Dummy )
return false;
}
return true;
}
CnMap m_cns;
RowMap m_rows;
VarMap m_vars;
EditMap m_edits;
std::vector<Symbol> m_infeasible_rows;
std::unique_ptr<Row> m_objective;
std::unique_ptr<Row> m_artificial;
Symbol::Id m_id_tick;
};
} // namespace impl
} // namespace kiwi

44
libs/kiwi/strength.h Normal file
View file

@ -0,0 +1,44 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#include <algorithm>
namespace kiwi
{
namespace strength
{
inline double create( double a, double b, double c, double w = 1.0 )
{
double result = 0.0;
result += std::max( 0.0, std::min( 1000.0, a * w ) ) * 1000000.0;
result += std::max( 0.0, std::min( 1000.0, b * w ) ) * 1000.0;
result += std::max( 0.0, std::min( 1000.0, c * w ) );
return result;
}
const double required = create( 1000.0, 1000.0, 1000.0 );
const double strong = create( 1.0, 0.0, 0.0 );
const double medium = create( 0.0, 1.0, 0.0 );
const double weak = create( 0.0, 0.0, 1.0 );
inline double clip( double value )
{
return std::max( 0.0, std::min( required, value ) );
}
} // namespace strength
} // namespace kiwi

68
libs/kiwi/symbol.h Normal file
View file

@ -0,0 +1,68 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
namespace kiwi
{
namespace impl
{
class Symbol
{
public:
typedef unsigned long long Id;
enum Type
{
Invalid,
External,
Slack,
Error,
Dummy
};
Symbol() : m_id( 0 ), m_type( Invalid ) {}
Symbol( Type type, Id id ) : m_id( id ), m_type( type ) {}
~Symbol() {}
Id id() const
{
return m_id;
}
Type type() const
{
return m_type;
}
private:
Id m_id;
Type m_type;
friend bool operator<( const Symbol& lhs, const Symbol& rhs )
{
return lhs.m_id < rhs.m_id;
}
friend bool operator==( const Symbol& lhs, const Symbol& rhs )
{
return lhs.m_id == rhs.m_id;
}
};
} // namespace impl
} // namespace kiwi

685
libs/kiwi/symbolics.h Normal file
View file

@ -0,0 +1,685 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#include <vector>
#include "constraint.h"
#include "expression.h"
#include "term.h"
#include "variable.h"
namespace kiwi
{
// Variable multiply, divide, and unary invert
inline
Term operator*( const Variable& variable, double coefficient )
{
return Term( variable, coefficient );
}
inline
Term operator/( const Variable& variable, double denominator )
{
return variable * ( 1.0 / denominator );
}
inline
Term operator-( const Variable& variable )
{
return variable * -1.0;
}
// Term multiply, divide, and unary invert
inline
Term operator*( const Term& term, double coefficient )
{
return Term( term.variable(), term.coefficient() * coefficient );
}
inline
Term operator/( const Term& term, double denominator )
{
return term * ( 1.0 / denominator );
}
inline
Term operator-( const Term& term )
{
return term * -1.0;
}
// Expression multiply, divide, and unary invert
inline
Expression operator*( const Expression& expression, double coefficient )
{
std::vector<Term> terms;
terms.reserve( expression.terms().size() );
typedef std::vector<Term>::const_iterator iter_t;
iter_t begin = expression.terms().begin();
iter_t end = expression.terms().end();
for( iter_t it = begin; it != end; ++it )
terms.push_back( ( *it ) * coefficient );
return Expression( terms, expression.constant() * coefficient );
}
inline
Expression operator/( const Expression& expression, double denominator )
{
return expression * ( 1.0 / denominator );
}
inline
Expression operator-( const Expression& expression )
{
return expression * -1.0;
}
// Double multiply
inline
Expression operator*( double coefficient, const Expression& expression )
{
return expression * coefficient;
}
inline
Term operator*( double coefficient, const Term& term )
{
return term * coefficient;
}
inline
Term operator*( double coefficient, const Variable& variable )
{
return variable * coefficient;
}
// Expression add and subtract
inline
Expression operator+( const Expression& first, const Expression& second )
{
std::vector<Term> terms;
terms.reserve( first.terms().size() + second.terms().size() );
terms.insert( terms.begin(), first.terms().begin(), first.terms().end() );
terms.insert( terms.end(), second.terms().begin(), second.terms().end() );
return Expression( terms, first.constant() + second.constant() );
}
inline
Expression operator+( const Expression& first, const Term& second )
{
std::vector<Term> terms;
terms.reserve( first.terms().size() + 1 );
terms.insert( terms.begin(), first.terms().begin(), first.terms().end() );
terms.push_back( second );
return Expression( terms, first.constant() );
}
inline
Expression operator+( const Expression& expression, const Variable& variable )
{
return expression + Term( variable );
}
inline
Expression operator+( const Expression& expression, double constant )
{
return Expression( expression.terms(), expression.constant() + constant );
}
inline
Expression operator-( const Expression& first, const Expression& second )
{
return first + -second;
}
inline
Expression operator-( const Expression& expression, const Term& term )
{
return expression + -term;
}
inline
Expression operator-( const Expression& expression, const Variable& variable )
{
return expression + -variable;
}
inline
Expression operator-( const Expression& expression, double constant )
{
return expression + -constant;
}
// Term add and subtract
inline
Expression operator+( const Term& term, const Expression& expression )
{
return expression + term;
}
inline
Expression operator+( const Term& first, const Term& second )
{
std::vector<Term> terms;
terms.reserve( 2 );
terms.push_back( first );
terms.push_back( second );
return Expression( terms );
}
inline
Expression operator+( const Term& term, const Variable& variable )
{
return term + Term( variable );
}
inline
Expression operator+( const Term& term, double constant )
{
return Expression( term, constant );
}
inline
Expression operator-( const Term& term, const Expression& expression )
{
return -expression + term;
}
inline
Expression operator-( const Term& first, const Term& second )
{
return first + -second;
}
inline
Expression operator-( const Term& term, const Variable& variable )
{
return term + -variable;
}
inline
Expression operator-( const Term& term, double constant )
{
return term + -constant;
}
// Variable add and subtract
inline
Expression operator+( const Variable& variable, const Expression& expression )
{
return expression + variable;
}
inline
Expression operator+( const Variable& variable, const Term& term )
{
return term + variable;
}
inline
Expression operator+( const Variable& first, const Variable& second )
{
return Term( first ) + second;
}
inline
Expression operator+( const Variable& variable, double constant )
{
return Term( variable ) + constant;
}
inline
Expression operator-( const Variable& variable, const Expression& expression )
{
return variable + -expression;
}
inline
Expression operator-( const Variable& variable, const Term& term )
{
return variable + -term;
}
inline
Expression operator-( const Variable& first, const Variable& second )
{
return first + -second;
}
inline
Expression operator-( const Variable& variable, double constant )
{
return variable + -constant;
}
// Double add and subtract
inline
Expression operator+( double constant, const Expression& expression )
{
return expression + constant;
}
inline
Expression operator+( double constant, const Term& term )
{
return term + constant;
}
inline
Expression operator+( double constant, const Variable& variable )
{
return variable + constant;
}
inline
Expression operator-( double constant, const Expression& expression )
{
return -expression + constant;
}
inline
Expression operator-( double constant, const Term& term )
{
return -term + constant;
}
inline
Expression operator-( double constant, const Variable& variable )
{
return -variable + constant;
}
// Expression relations
inline
Constraint operator==( const Expression& first, const Expression& second )
{
return Constraint( first - second, OP_EQ );
}
inline
Constraint operator==( const Expression& expression, const Term& term )
{
return expression == Expression( term );
}
inline
Constraint operator==( const Expression& expression, const Variable& variable )
{
return expression == Term( variable );
}
inline
Constraint operator==( const Expression& expression, double constant )
{
return expression == Expression( constant );
}
inline
Constraint operator<=( const Expression& first, const Expression& second )
{
return Constraint( first - second, OP_LE );
}
inline
Constraint operator<=( const Expression& expression, const Term& term )
{
return expression <= Expression( term );
}
inline
Constraint operator<=( const Expression& expression, const Variable& variable )
{
return expression <= Term( variable );
}
inline
Constraint operator<=( const Expression& expression, double constant )
{
return expression <= Expression( constant );
}
inline
Constraint operator>=( const Expression& first, const Expression& second )
{
return Constraint( first - second, OP_GE );
}
inline
Constraint operator>=( const Expression& expression, const Term& term )
{
return expression >= Expression( term );
}
inline
Constraint operator>=( const Expression& expression, const Variable& variable )
{
return expression >= Term( variable );
}
inline
Constraint operator>=( const Expression& expression, double constant )
{
return expression >= Expression( constant );
}
// Term relations
inline
Constraint operator==( const Term& term, const Expression& expression )
{
return expression == term;
}
inline
Constraint operator==( const Term& first, const Term& second )
{
return Expression( first ) == second;
}
inline
Constraint operator==( const Term& term, const Variable& variable )
{
return Expression( term ) == variable;
}
inline
Constraint operator==( const Term& term, double constant )
{
return Expression( term ) == constant;
}
inline
Constraint operator<=( const Term& term, const Expression& expression )
{
return expression >= term;
}
inline
Constraint operator<=( const Term& first, const Term& second )
{
return Expression( first ) <= second;
}
inline
Constraint operator<=( const Term& term, const Variable& variable )
{
return Expression( term ) <= variable;
}
inline
Constraint operator<=( const Term& term, double constant )
{
return Expression( term ) <= constant;
}
inline
Constraint operator>=( const Term& term, const Expression& expression )
{
return expression <= term;
}
inline
Constraint operator>=( const Term& first, const Term& second )
{
return Expression( first ) >= second;
}
inline
Constraint operator>=( const Term& term, const Variable& variable )
{
return Expression( term ) >= variable;
}
inline
Constraint operator>=( const Term& term, double constant )
{
return Expression( term ) >= constant;
}
// Variable relations
inline
Constraint operator==( const Variable& variable, const Expression& expression )
{
return expression == variable;
}
inline
Constraint operator==( const Variable& variable, const Term& term )
{
return term == variable;
}
inline
Constraint operator==( const Variable& first, const Variable& second )
{
return Term( first ) == second;
}
inline
Constraint operator==( const Variable& variable, double constant )
{
return Term( variable ) == constant;
}
inline
Constraint operator<=( const Variable& variable, const Expression& expression )
{
return expression >= variable;
}
inline
Constraint operator<=( const Variable& variable, const Term& term )
{
return term >= variable;
}
inline
Constraint operator<=( const Variable& first, const Variable& second )
{
return Term( first ) <= second;
}
inline
Constraint operator<=( const Variable& variable, double constant )
{
return Term( variable ) <= constant;
}
inline
Constraint operator>=( const Variable& variable, const Expression& expression )
{
return expression <= variable;
}
inline
Constraint operator>=( const Variable& variable, const Term& term )
{
return term <= variable;
}
inline
Constraint operator>=( const Variable& first, const Variable& second )
{
return Term( first ) >= second;
}
inline
Constraint operator>=( const Variable& variable, double constant )
{
return Term( variable ) >= constant;
}
// Double relations
inline
Constraint operator==( double constant, const Expression& expression )
{
return expression == constant;
}
inline
Constraint operator==( double constant, const Term& term )
{
return term == constant;
}
inline
Constraint operator==( double constant, const Variable& variable )
{
return variable == constant;
}
inline
Constraint operator<=( double constant, const Expression& expression )
{
return expression >= constant;
}
inline
Constraint operator<=( double constant, const Term& term )
{
return term >= constant;
}
inline
Constraint operator<=( double constant, const Variable& variable )
{
return variable >= constant;
}
inline
Constraint operator>=( double constant, const Expression& expression )
{
return expression <= constant;
}
inline
Constraint operator>=( double constant, const Term& term )
{
return term <= constant;
}
inline
Constraint operator>=( double constant, const Variable& variable )
{
return variable <= constant;
}
// Constraint strength modifier
inline
Constraint operator|( const Constraint& constraint, double strength )
{
return Constraint( constraint, strength );
}
inline
Constraint operator|( double strength, const Constraint& constraint )
{
return constraint | strength;
}
} // namespace kiwi

51
libs/kiwi/term.h Normal file
View file

@ -0,0 +1,51 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#include <utility>
#include "variable.h"
namespace kiwi
{
class Term
{
public:
Term( const Variable& variable, double coefficient = 1.0 ) :
m_variable( variable ), m_coefficient( coefficient ) {}
// to facilitate efficient map -> vector copies
Term( const std::pair<const Variable, double>& pair ) :
m_variable( pair.first ), m_coefficient( pair.second ) {}
~Term() {}
const Variable& variable() const
{
return m_variable;
}
double coefficient() const
{
return m_coefficient;
}
double value() const
{
return m_coefficient * m_variable.value();
}
private:
Variable m_variable;
double m_coefficient;
};
} // namespace kiwi

24
libs/kiwi/util.h Normal file
View file

@ -0,0 +1,24 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
namespace kiwi
{
namespace impl
{
inline bool nearZero(double value)
{
const double eps = 1.0e-8;
return value < 0.0 ? -value < eps : value < eps;
}
} // namespace impl
} // namespace kiwi

111
libs/kiwi/variable.h Normal file
View file

@ -0,0 +1,111 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2017, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#include <memory>
#include <string>
#include "shareddata.h"
namespace kiwi
{
class Variable
{
public:
class Context
{
public:
Context() {}
virtual ~Context() {} // LCOV_EXCL_LINE
};
Variable(Context *context = 0) : m_data(new VariableData("", context)) {}
Variable(const std::string &name, Context *context = 0) : m_data(new VariableData(name, context)) {}
Variable(const char *name, Context *context = 0) : m_data(new VariableData(name, context)) {}
~Variable() {}
const std::string &name() const
{
return m_data->m_name;
}
void setName(const char *name)
{
m_data->m_name = name;
}
void setName(const std::string &name)
{
m_data->m_name = name;
}
Context *context() const
{
return m_data->m_context.get();
}
void setContext(Context *context)
{
m_data->m_context.reset(context);
}
double value() const
{
return m_data->m_value;
}
void setValue(double value)
{
m_data->m_value = value;
}
// operator== is used for symbolics
bool equals(const Variable &other)
{
return m_data == other.m_data;
}
private:
class VariableData : public SharedData
{
public:
VariableData(const std::string &name, Context *context) : SharedData(),
m_name(name),
m_context(context),
m_value(0.0) {}
VariableData(const char *name, Context *context) : SharedData(),
m_name(name),
m_context(context),
m_value(0.0) {}
~VariableData() {}
std::string m_name;
std::unique_ptr<Context> m_context;
double m_value;
private:
VariableData(const VariableData &other);
VariableData &operator=(const VariableData &other);
};
SharedDataPtr<VariableData> m_data;
friend bool operator<(const Variable &lhs, const Variable &rhs)
{
return lhs.m_data < rhs.m_data;
}
};
} // namespace kiwi

14
libs/kiwi/version.h Normal file
View file

@ -0,0 +1,14 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2020, Nucleic Development Team.
|
| Distributed under the terms of the Modified BSD License.
|
| The full license is in the file LICENSE, distributed with this software.
|----------------------------------------------------------------------------*/
#pragma once
#define KIWI_MAJOR_VERSION 1
#define KIWI_MINOR_VERSION 2
#define KIWI_MICRO_VERSION 0
#define KIWI_VERSION_HEX 0x010200
#define KIWI_VERSION "1.2.0"