jkl1337/ljkiwi
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Squashed 'kiwi/' content from commit 268028e

Browse Source 
git-subtree-dir: kiwi
git-subtree-split: 268028ee4a694dcd89e4b1e683bf2f9ac48c08d9
This commit is contained in:
John K. Luebs
2024-02-11 15:32:50 -06:00
commit 81396a5322
76 changed files with 13184 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

352
kiwi/AssocVector.h Normal file
View File

@@ -0,0 +1,352 @@
////////////////////////////////////////////////////////////////////////////////
// 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 K>
class AssocVectorCompare : public C
{
typedef std::pair<K, Value>
Data;
typedef K 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, K>
{
typedef std::vector<std::pair<K, V>, A> Base;
typedef Private::AssocVectorCompare<V, C, K> 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 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 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

140
kiwi/constraint.h Normal file
View File

@@ -0,0 +1,140 @@
/*-----------------------------------------------------------------------------
| 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 <cstdlib>
#include <map>
#include <vector>
#include "expression.h"
#include "shareddata.h"
#include "strength.h"
#include "term.h"
#include "variable.h"
#include "util.h"
namespace kiwi
{
enum RelationalOperator
{
OP_LE,
OP_GE,
OP_EQ
};
class Constraint
{
public:
Constraint() = default;
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 Constraint &) = default;
Constraint(Constraint &&) noexcept = default;
~Constraint() = default;
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 violated() const
{
switch (m_data->m_op)
{
case OP_EQ: return !impl::nearZero(m_data->m_expression.value());
case OP_GE: return m_data->m_expression.value() < 0.0;
case OP_LE: return m_data->m_expression.value() > 0.0;
}
std::abort();
}
bool operator!() const
{
return !m_data;
}
Constraint& operator=(const Constraint &) = default;
Constraint& operator=(Constraint &&) noexcept = default;
private:
static Expression reduce(const Expression &expr)
{
std::map<Variable, double> vars;
for (const auto & term : expr.terms())
vars[term.variable()] += term.coefficient();
std::vector<Term> terms(vars.begin(), vars.end());
return Expression(std::move(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() = default;
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

184
kiwi/debug.h Normal file
View File

@@ -0,0 +1,184 @@
/*-----------------------------------------------------------------------------
| 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)
{
for (const auto &rowPair : rows)
{
dump(rowPair.first, out);
out << " | ";
dump(*rowPair.second, out);
}
}
static void dump(const std::vector<Symbol> &symbols, std::ostream &out)
{
for (const auto &symbol : symbols)
{
dump(symbol, out);
out << std::endl;
}
}
static void dump(const SolverImpl::VarMap &vars, std::ostream &out)
{
for (const auto &varPair : vars)
{
out << varPair.first.name() << " = ";
dump(varPair.second, out);
out << std::endl;
}
}
static void dump(const SolverImpl::CnMap &cns, std::ostream &out)
{
for (const auto &cnPair : cns)
dump(cnPair.first, out);
}
static void dump(const SolverImpl::EditMap &edits, std::ostream &out)
{
for (const auto &editPair : edits)
out << editPair.first.name() << std::endl;
}
static void dump(const Row &row, std::ostream &out)
{
for (const auto &rowPair : row.cells())
{
out << " + " << rowPair.second << " * ";
dump(rowPair.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)
{
for (const auto &term : cn.expression().terms())
{
out << term.coefficient() << " * ";
out << term.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
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(Constraint constraint) : m_constraint(std::move(constraint)) {}
~UnsatisfiableConstraint() noexcept {}
const char *what() const noexcept
{
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(Constraint constraint) : m_constraint(std::move(constraint)) {}
~UnknownConstraint() noexcept {}
const char *what() const noexcept
{
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(Constraint constraint) : m_constraint(std::move(constraint)) {}
~DuplicateConstraint() noexcept {}
const char *what() const noexcept
{
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(Variable variable) : m_variable(std::move(variable)) {}
~UnknownEditVariable() noexcept {}
const char *what() const noexcept
{
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(Variable variable) : m_variable(std::move(variable)) {}
~DuplicateEditVariable() noexcept {}
const char *what() const noexcept
{
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() noexcept {}
const char *what() const noexcept
{
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(std::string msg) : m_msg(std::move(msg)) {}
~InternalSolverError() noexcept {}
const char *what() const noexcept
{
return m_msg.c_str();
}
private:
std::string m_msg;
};
} // namespace kiwi

62
kiwi/expression.h Normal file
View File

@@ -0,0 +1,62 @@
/*-----------------------------------------------------------------------------
| 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(std::vector<Term> terms, double constant = 0.0) : m_terms(std::move(terms)), m_constant(constant) {}
Expression(const Expression&) = default;
// Could be marked noexcept but for a bug in the GCC of the manylinux1 image
Expression(Expression&&) = default;
~Expression() = default;
const std::vector<Term> &terms() const
{
return m_terms;
}
double constant() const
{
return m_constant;
}
double value() const
{
double result = m_constant;
for (const Term &term : m_terms)
result += term.value();
return result;
}
Expression& operator=(const Expression&) = default;
// Could be marked noexcept but for a bug in the GCC of the manylinux1 image
Expression& operator=(Expression&&) = default;
private:
std::vector<Term> m_terms;
double m_constant;
};
} // namespace kiwi

19
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
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

182
kiwi/row.h Normal file
View File

@@ -0,0 +1,182 @@
/*-----------------------------------------------------------------------------
| 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:
using CellMap = MapType<Symbol, double>;
Row() : Row(0.0) {}
Row(double constant) : m_constant(constant) {}
Row(const Row &other) = default;
~Row() = default;
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)
{
m_constant += other.m_constant * coefficient;
for (const auto & cellPair : other.m_cells)
{
double coeff = cellPair.second * coefficient;
if (nearZero(m_cells[cellPair.first] += coeff))
m_cells.erase(cellPair.first);
}
}
/* Remove the given symbol from the row.
*/
void remove(const Symbol &symbol)
{
auto 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()
{
m_constant = -m_constant;
for (auto &cellPair : m_cells)
cellPair.second = -cellPair.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)
{
double coeff = -1.0 / m_cells[symbol];
m_cells.erase(symbol);
m_constant *= coeff;
for (auto &cellPair : m_cells)
cellPair.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)
{
auto 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

181
kiwi/shareddata.h Normal file
View File

@@ -0,0 +1,181 @@
/*-----------------------------------------------------------------------------
| 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
/*
Implementation note
===================
SharedDataPtr/SharedData offer the same basic functionality as std::shared_ptr,
but do not use atomic counters under the hood.
Since kiwi operates within a single thread context, atomic counters are not necessary,
especially given the extra CPU cost.
Therefore the use of SharedDataPtr/SharedData is preferred over std::shared_ptr.
*/
namespace kiwi
{
class SharedData
{
public:
SharedData() : m_refcount(0) {}
SharedData(const SharedData &other) = delete;
SharedData(SharedData&& other) = delete;
int m_refcount;
SharedData &operator=(const SharedData &other) = delete;
SharedData &operator=(SharedData&& other) = delete;
};
template <typename T>
class SharedDataPtr
{
public:
using Type = T;
SharedDataPtr() : m_data(nullptr) {}
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(SharedDataPtr&& other) noexcept : m_data(other.m_data)
{
other.m_data = nullptr;
}
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=(SharedDataPtr<T>&& other) noexcept
{
if (m_data != other.m_data)
{
T *temp = m_data;
m_data = other.m_data;
other.m_data = nullptr;
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
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() = default;
~Solver() = default;
/* 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

825
kiwi/solverimpl.h Normal file
View File

@@ -0,0 +1,825 @@
/*-----------------------------------------------------------------------------
| 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;
};
using VarMap = MapType<Variable, Symbol>;
using RowMap = MapType<Symbol, Row*>;
using CnMap = MapType<Constraint, Tag>;
using EditMap = MapType<Variable, EditInfo>;
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( const SolverImpl& ) = delete;
SolverImpl( SolverImpl&& ) = delete;
~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 )
{
auto 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.
auto 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 )
{
auto 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 )
{
auto 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.
auto 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.
for (const auto & rowPair : m_rows)
{
double coeff = rowPair.second->coefficientFor( info.tag.marker );
if( coeff != 0.0 &&
rowPair.second->add( delta * coeff ) < 0.0 &&
rowPair.first.type() != Symbol::External )
m_infeasible_rows.push_back( rowPair.first );
}
}
/* Update the values of the external solver variables.
*/
void updateVariables()
{
auto row_end = m_rows.end();
for (auto &varPair : m_vars)
{
Variable& var = varPair.first;
auto row_it = m_rows.find( varPair.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;
}
SolverImpl& operator=( const SolverImpl& ) = delete;
SolverImpl& operator=( SolverImpl&& ) = delete;
private:
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 )
{
auto 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.
*/
std::unique_ptr<Row> createRow( const Constraint& constraint, Tag& tag )
{
const Expression& expr( constraint.expression() );
std::unique_ptr<Row> row( new Row( expr.constant() ) );
// Substitute the current basic variables into the row.
for (const auto &term : expr.terms())
{
if( !nearZero( term.coefficient() ) )
{
Symbol symbol( getVarSymbol( term.variable() ) );
auto row_it = m_rows.find( symbol );
if( row_it != m_rows.end() )
row->insert( *row_it->second, term.coefficient() );
else
row->insert( symbol, term.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 ) const
{
for (const auto &cellPair : row.cells())
{
if( cellPair.first.type() == Symbol::External )
return cellPair.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.
auto 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.
for (auto &rowPair : m_rows)
rowPair.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 )
{
for( auto& rowPair : m_rows )
{
rowPair.second->substitute( symbol, row );
if( rowPair.first.type() != Symbol::External &&
rowPair.second->constant() < 0.0 )
m_infeasible_rows.push_back( rowPair.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;
auto 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();
auto 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 ) const
{
for (const auto &cellPair : objective.cells())
{
if( cellPair.first.type() != Symbol::Dummy && cellPair.second < 0.0 )
return cellPair.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 ) const
{
Symbol entering;
double ratio = std::numeric_limits<double>::max();
for (const auto &cellPair : row.cells())
{
if( cellPair.second > 0.0 && cellPair.first.type() != Symbol::Dummy )
{
double coeff = m_objective->coefficientFor( cellPair.first );
double r = coeff / cellPair.second;
if( r < ratio )
{
ratio = r;
entering = cellPair.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 ) const
{
for (const auto &cellPair : row.cells())
{
const Symbol& sym( cellPair.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 )
{
double ratio = std::numeric_limits<double>::max();
auto end = m_rows.end();
auto found = m_rows.end();
for( auto 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();
double r1 = dmax;
double r2 = dmax;
auto end = m_rows.end();
auto first = end;
auto second = end;
auto third = end;
for( auto 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 )
{
auto 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 ) const
{
for (const auto &rowPair : row.cells())
{
if( rowPair.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
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
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:
using Id = unsigned long long;
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() = default;
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

680
kiwi/symbolics.h Normal file
View File

@@ -0,0 +1,680 @@
/*-----------------------------------------------------------------------------
| 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() );
for (const Term &term : expression.terms())
terms.push_back(term * coefficient);
return Expression( std::move(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( std::move(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( std::move(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 )
{
return Expression( { first, second } );
}
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

59
kiwi/term.h Normal file
View File

@@ -0,0 +1,59 @@
/*-----------------------------------------------------------------------------
| 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( Variable variable, double coefficient = 1.0 ) :
m_variable( std::move(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 Term&) = default;
Term(Term&&) noexcept = default;
~Term() = default;
const Variable& variable() const
{
return m_variable;
}
double coefficient() const
{
return m_coefficient;
}
double value() const
{
return m_coefficient * m_variable.value();
}
Term& operator=(const Term&) = default;
Term& operator=(Term&&) noexcept = default;
private:
Variable m_variable;
double m_coefficient;
};
} // namespace kiwi

24
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

119
kiwi/variable.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 <memory>
#include <string>
#include "shareddata.h"
namespace kiwi
{
class Variable
{
public:
class Context
{
public:
Context() = default;
virtual ~Context() {} // LCOV_EXCL_LINE
};
Variable(Context *context = 0) : m_data(new VariableData("", context)) {}
Variable(std::string name, Context *context = 0) : m_data(new VariableData(std::move(name), context)) {}
Variable(const char *name, Context *context = 0) : m_data(new VariableData(name, context)) {}
Variable(const Variable&) = default;
Variable(Variable&&) noexcept = default;
~Variable() = default;
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;
}
Variable& operator=(const Variable&) = default;
Variable& operator=(Variable&&) noexcept = default;
private:
class VariableData : public SharedData
{
public:
VariableData(std::string name, Context *context) : SharedData(),
m_name(std::move(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() = default;
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
kiwi/version.h Normal file
View File

@@ -0,0 +1,14 @@
/*-----------------------------------------------------------------------------
| Copyright (c) 2013-2022, 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 4
#define KIWI_MICRO_VERSION 2
#define KIWI_VERSION_HEX 0x010402
#define KIWI_VERSION "1.4.2"