Merge branch 'master' into nazara-next

thirdparty/include/Newton/dgTypes.h

@@ -0,0 +1,696 @@
+/* Copyright (c) <2003-2019> <Julio Jerez, Newton Game Dynamics>
+* 
+* This software is provided 'as-is', without any express or implied
+* warranty. In no event will the authors be held liable for any damages
+* arising from the use of this software.
+* 
+* Permission is granted to anyone to use this software for any purpose,
+* including commercial applications, and to alter it and redistribute it
+* freely, subject to the following restrictions:
+* 
+* 1. The origin of this software must not be misrepresented; you must not
+* claim that you wrote the original software. If you use this software
+* in a product, an acknowledgment in the product documentation would be
+* appreciated but is not required.
+* 
+* 2. Altered source versions must be plainly marked as such, and must not be
+* misrepresented as being the original software.
+* 
+* 3. This notice may not be removed or altered from any source distribution.
+*/
+
+#ifndef __DGTYPES_H__
+#define __DGTYPES_H__
+
+#ifdef _MSC_VER 
+	#ifdef _M_ARM
+		#ifndef _ARM_VER
+			#define _ARM_VER
+		#endif
+	#elif defined (_M_X64)
+		#ifndef _WIN_64_VER
+			#define _WIN_64_VER
+		#endif
+	#elif defined (_M_IX86)
+		#ifndef _WIN_32_VER
+			#define _WIN_32_VER
+		#endif
+	#else 
+		#error target platform not defined
+	#endif
+
+	#if _MSC_VER >= 1400
+		#define HAVE_STRUCT_TIMESPEC
+	#endif
+#endif
+
+#if (defined (_WIN_32_VER) || defined (_WIN_64_VER) || (defined (_MSC_VER ) && defined (_ARM_VER)) )
+	#pragma warning (disable: 4100) //unreferenced formal parameter
+	#pragma warning (disable: 4201) //nonstandard extension used : nameless struct/union
+	#pragma warning (disable: 4324) //structure was padded due to __declspec(align())
+	#pragma warning (disable: 4514) //unreferenced inline function has been removed
+	#pragma warning (disable: 4530) //C++ exception handler used, but unwind semantics are not enabled. Specify /EHsc
+	#pragma warning (disable: 4625) //copy constructor could not be generated because a base class copy constructor is inaccessible or deleted
+	#pragma warning (disable: 4626) //assignment operator could not be generated because a base class assignment operator is inaccessible or deleted
+	#pragma warning (disable: 4640) //construction of local static object is not thread-safe
+	#pragma warning (disable: 4820) //bytes padding added after data member
+	#pragma warning (disable: 4005) //'__useHeader': macro redefinition
+	//#pragma warning (disable: 4577) // 'noexcept' used with no exception handling mode specified; termination on exception is not guaranteed. Specify /EHsc
+    
+	#include <io.h> 
+	#include <direct.h> 
+	#include <malloc.h>
+	#include <float.h>
+	#include <stdarg.h>
+	#include <process.h>
+
+	#ifdef _DEBUG
+		#pragma warning (disable: 4127)	//conditional expression is constant
+	#endif
+
+	#pragma warning (push, 3) 
+		#include <windows.h>
+		#include <crtdbg.h>
+		//#ifndef _DURANGO
+		//	#include <tlhelp32.h>
+		//#endif
+	#pragma warning (pop) 
+#endif
+
+#include <new>
+#include <stdio.h>
+#include <stdlib.h>
+#include <string.h>
+#include <stdarg.h>
+#include <time.h>
+#include <math.h>
+#include <float.h>
+#include <ctype.h>
+
+#if (defined (__MINGW32__) || defined (__MINGW64__))
+	#include <io.h> 
+	#include <direct.h> 
+	#include <malloc.h>
+	#include <float.h>
+	#include <windows.h>
+	#include <process.h>
+#endif
+
+#if (defined (_WIN_32_VER) || defined (_WIN_64_VER))
+	#include <intrin.h>
+	#include <emmintrin.h> 
+	#include <pmmintrin.h>
+#endif
+
+#if (defined (_POSIX_VER) || defined (_POSIX_VER_64) || defined (__MINGW32__) || defined (__MINGW64__))
+  /* CMake defines NDEBUG for _not_ debug builds. Therefore, set
+     Newton's _DEBUG flag only when NDEBUG is not defined.
+  */
+  #ifndef NDEBUG
+    #define _DEBUG 1
+  #endif
+
+	#include <unistd.h>
+	#include <assert.h>
+	#if (!defined(__arm__) && !defined(__aarch64__)) // it was __ARMCC_VERSION before, it should be __ARM__ or aarch64, otherwise cross compiling in gcc fails.
+		extern "C" 
+		{ 
+			// for SSE3 and up
+			#include <pmmintrin.h> 
+			#include <emmintrin.h> 
+			#include <mmintrin.h> 
+		} 
+	#endif
+#endif
+
+#ifdef _MACOSX_VER
+	#include <unistd.h>
+	#include <sys/sysctl.h>
+    #include <assert.h> 
+    #if (defined __i386__ || defined __x86_64__)
+		#include <fenv.h>
+		#include <pmmintrin.h> 
+		#include <emmintrin.h>  //sse3
+        #include <mmintrin.h> 
+    #endif
+#endif
+
+// uncomment out D_PROFILER to enable profiler frame capture profiler traces
+// alternatively the end application can use a command line option to enable this define
+//#define D_PROFILER
+
+// uncomment this for Scalar floating point 
+// alternatively the end application can use a command line option to enable this define
+//#define DG_SCALAR_VECTOR_CLASS
+
+// uncomment this for Scalar floating point 
+// alternatively the end application can use a command line option to enable this define
+//#define __ANDROID__
+
+// by default newton run on a separate thread and 
+// optionally concurrent with the calling thread,
+// it also uses a thread job pool for multi core systems.
+// define DG_USE_THREAD_EMULATION on the command line for 
+// platform that do not support hardware multi threading or 
+// if the and application want to control threading at the application level 
+//#define DG_USE_THREAD_EMULATION
+
+#if (defined (_WIN_32_VER) || defined (_WIN_64_VER))
+	#if _MSC_VER < 1700
+		#ifndef DG_USE_THREAD_EMULATION
+			#define DG_USE_THREAD_EMULATION
+		#endif
+	#endif
+#endif
+
+#ifndef DG_USE_THREAD_EMULATION
+	#include <mutex>
+	#include <thread>
+	#include <condition_variable>
+#endif
+
+
+//************************************************************
+#ifdef DG_DISABLE_ASSERT
+	#define dgAssert(x)
+#else 
+	#if defined (_WIN_32_VER) || defined (_WIN_64_VER) 
+		#define dgAssert(x) _ASSERTE(x)
+	#elif defined (_MSC_VER ) && defined (_ARM_VER) 
+		#define dgAssert(x) _ASSERTE(x)
+	#else
+		#ifdef _DEBUG
+			#define dgAssert(x) assert(x)
+		#else 
+			#define dgAssert(x)
+		#endif
+	#endif
+#endif
+
+
+#define	DG_MAX_THREADS_HIVE_COUNT		16
+
+#ifdef _DEBUG
+//#define __ENABLE_DG_CONTAINERS_SANITY_CHECK 
+#endif
+
+
+#ifdef DLL_DECLSPEC
+#undef DLL_DECLSPEC
+#endif
+
+#ifdef _DEBUG
+	#define DG_INLINE inline
+#else
+	#if defined(_MSC_VER)
+			#define DG_INLINE __forceinline 
+	#else 
+		#define DG_INLINE	inline
+		//#define DG_INLINE	 __attribute__((always_inline))
+	#endif
+#endif
+
+
+#define DG_VECTOR_SIMD_SIZE		16
+#define DG_VECTOR_AVX2_SIZE		32
+
+#if defined(_MSC_VER)
+	#define	DG_GCC_VECTOR_ALIGNMENT
+	#define	DG_MSC_VECTOR_ALIGNMENT    __declspec(align(DG_VECTOR_SIMD_SIZE))
+#else
+	#define	DG_GCC_VECTOR_ALIGNMENT    __attribute__ ((aligned (DG_VECTOR_SIMD_SIZE)))
+	#define	DG_MSC_VECTOR_ALIGNMENT
+#endif
+
+#if defined(_MSC_VER)
+	#define	DG_GCC_AVX_ALIGNMENT
+	#define	DG_MSC_AVX_ALIGNMENT       __declspec(align(DG_VECTOR_AVX2_SIZE))
+#else
+	#define	DG_GCC_AVX_ALIGNMENT       __attribute__ ((aligned (DG_VECTOR_AVX2_SIZE)))
+	#define	DG_MSC_AVX_ALIGNMENT
+#endif
+
+#if defined(_MSC_VER)
+	#define DG_LIBRARY_EXPORT __declspec(dllexport)
+	#define DG_LIBRARY_IMPORT __declspec(dllimport)
+	#define DG_LIBRARY_STATIC
+#else
+	#define DG_LIBRARY_EXPORT __attribute__((visibility("default")))
+	#define DG_LIBRARY_IMPORT __attribute__((visibility("default")))
+	#define DG_LIBRARY_STATIC
+#endif
+
+
+#if ((defined (_WIN_32_VER) || defined (_WIN_64_VER)) && (_MSC_VER  >= 1600))
+	#include <stdint.h>
+	typedef int8_t dgInt8;
+	typedef uint8_t dgUnsigned8;
+
+	typedef int16_t dgInt16;
+	typedef uint16_t dgUnsigned16;
+
+	typedef int32_t dgInt32;
+	typedef uint32_t dgUnsigned32;
+
+	typedef int64_t dgInt64;
+	typedef uint64_t dgUnsigned64;
+#else
+	typedef char dgInt8;
+	typedef unsigned char dgUnsigned8;
+
+	typedef short dgInt16;
+	typedef unsigned short dgUnsigned16;
+
+	typedef int dgInt32;
+	typedef unsigned dgUnsigned32;
+	typedef unsigned int dgUnsigned32;
+
+	typedef long long dgInt64;
+	typedef unsigned long long dgUnsigned64;
+	typedef double dgFloat64;
+#endif
+
+
+typedef double dgFloat64;
+
+#ifdef _NEWTON_USE_DOUBLE
+	typedef double dgFloat32;
+#else
+	typedef float dgFloat32;
+#endif
+
+
+class dgTriplex
+{
+	public:
+	dgFloat32 m_x;
+	dgFloat32 m_y;
+	dgFloat32 m_z;
+};
+
+#define dgPi		 	dgFloat32 (3.141592f)
+#define dgPI2		 	dgFloat32 (dgPi * 2.0f)
+#define dgEXP		 	dgFloat32 (2.71828f)
+#define dgEpsilon	  	dgFloat32 (1.0e-5f)
+#define dgDegreeToRad  	dgFloat32 (dgPi / 180.0f)
+#define dgRadToDegree  	dgFloat32 (180.0f / dgPi)
+
+class dgBigVector;
+#ifndef _NEWTON_USE_DOUBLE
+class dgVector;
+#endif 
+
+
+#if (defined (_WIN_32_VER) || defined (_WIN_64_VER))
+	#define dgApi __cdecl 	
+	#define dgStdApi __stdcall 	
+#else
+	#define dgApi 	
+	#define dgStdApi
+#endif
+
+
+
+#if (defined (_WIN_32_VER) || defined (_WIN_64_VER))
+	#define dgCheckFloat(x) (_finite(x) && !_isnan(x))
+//	#define dgCheckFloat(x) 1
+#else
+	#define dgCheckFloat(x) (isfinite(x) && !isnan(x))
+//		#define dgCheckFloat(x) 1
+#endif
+
+typedef void (dgApi *dgDeserialize) (void* const userData, void* buffer, dgInt32 size);
+typedef void (dgApi *dgSerialize) (void* const userData, const void* const buffer, dgInt32 size);
+typedef bool (dgApi *dgReportProgress) (dgFloat32 progressNormalzedPercent, void* const userData);
+
+// assume this function returns memory aligned to 16 bytes
+#define dgAlloca(type, count) (type*) alloca (sizeof (type) * (count))
+
+DG_INLINE dgInt32 dgExp2 (dgInt32 x)
+{
+	dgInt32 exp;
+	for (exp = -1; x; x >>= 1) {
+		exp ++;
+	}
+	return exp;
+}
+
+DG_INLINE dgInt32 dgBitReversal(dgInt32 v, dgInt32 base)
+{
+	dgInt32 x = 0;
+	dgInt32 power = dgExp2 (base) - 1;
+	do {
+		x += (v & 1) << power;
+		v >>= 1;
+		power--;
+	} while (v);
+	dgAssert(x < base);
+	return x;
+}
+
+
+template <class T> 
+DG_INLINE T dgMin(T A, T B)
+{
+	return (A < B) ? A : B; 
+}
+
+template <class T> 
+DG_INLINE T dgMax(T A, T B)
+{
+	return (A > B) ? A : B; 
+}
+
+template <class T>
+DG_INLINE T dgMin(T A, T B, T C)
+{
+	return dgMin(dgMin (A, B), C);
+}
+
+template <class T>
+DG_INLINE T dgMax(T A, T B, T C)
+{
+	return dgMax(dgMax (A, B), C);
+}
+
+template <class T>
+DG_INLINE T dgClamp(T val, T min, T max)
+{
+	return dgMax (min, dgMin (max, val));
+}
+
+template <class T> 
+DG_INLINE void dgSwap(T& A, T& B)
+{
+	T tmp (A);
+	A = B;
+	B = tmp;
+}	
+
+template <class T>
+DG_INLINE T dgAbs(T A)
+{
+	// according to Intel this is better because is does not read after write
+	return (A >= T(0)) ? A : -A;
+}
+
+template <class T>
+DG_INLINE T dgSign(T A)
+{
+	return (A >= T(0)) ? T(1) : T(-1);
+}
+
+template <class T> 
+DG_INLINE bool dgAreEqual(T A, T B, T tol)
+{
+	if ((dgAbs(A) < tol) && (dgAbs(B) < tol)) {
+		return true;
+	}
+/*
+	dgInt32 exp0;
+	dgFloat64 mantissa0 = frexp(dgFloat64 (A), &exp0);
+
+	dgInt32 exp1;
+	dgFloat64 mantissa1 = frexp(dgFloat64(B), &exp1);
+
+	if ((exp0 < -12) && (exp1 < -12)) {
+		return true;
+	}
+
+	if (exp0 != exp1) {
+		return false;
+	}
+	return dgAbs(mantissa0 - mantissa1) < tol;
+*/	
+	T den = dgMax(dgAbs(A), dgAbs(B)) + tol;
+	A /= den;
+	B /= den;
+	return dgAbs(A - B) < tol;
+}
+
+#ifdef _NEWTON_USE_DOUBLE
+	union dgFloatSign
+	{
+		struct {
+			dgInt32 m_dommy;
+			dgInt32 m_iVal;
+		} m_integer;
+		dgFloat64 m_fVal;
+	};
+#else
+	union dgFloatSign
+	{
+		struct {
+			dgInt32 m_iVal;
+		} m_integer;
+		dgFloat32 m_fVal;
+	};
+#endif
+
+union dgDoubleInt
+{
+	struct {
+		dgInt32 m_intL;
+		dgInt32 m_intH;
+	};
+	void* m_ptr;
+	dgInt64 m_int;
+	dgFloat64 m_float;
+};
+
+
+void dgGetMinMax (dgBigVector &Min, dgBigVector &Max, const dgFloat64* const vArray, dgInt32 vCount, dgInt32 strideInBytes);
+dgInt32 dgVertexListToIndexList (dgFloat64* const vertexList, dgInt32 strideInBytes, dgInt32 compareCount,     dgInt32 vertexCount,         dgInt32* const indexListOut, dgFloat64 tolerance = dgEpsilon);
+dgInt32 dgVertexListToIndexList (dgFloat32* const vertexList, dgInt32 strideInBytes, dgInt32 floatSizeInBytes, dgInt32 unsignedSizeInBytes, dgInt32 vertexCount, dgInt32* const indexListOut, dgFloat32 tolerance = dgEpsilon);
+
+#define PointerToInt(x) ((size_t)x)
+#define IntToPointer(x) ((void*)(size_t(x)))
+
+
+#ifndef _MSC_VER 
+	#define _stricmp(x,y) strcasecmp(x,y)
+#endif
+
+#define dgSqrt(x)			dgFloat32 (sqrt(x))	
+#define dgSin(x)			dgFloat32 (sin(x))
+#define dgCos(x)			dgFloat32 (cos(x))
+#define dgAsin(x)			dgFloat32 (asin(x))
+#define dgAcos(x)			dgFloat32 (acos(x))
+#define dgLog(x)			dgFloat32 (log(x))
+#define dgCeil(x)			dgFloat32 (ceil(x))
+#define dgFloor(x)			dgFloat32 (floor(x))	
+#define dgPow(x,y)			dgFloat32 (pow(x,y))
+#define dgFmod(x,y)			dgFloat32 (fmod(x,y))
+#define dgAtan2(x,y)		dgFloat32 (atan2(x,y))
+#define dgRsqrt(x)			(dgFloat32 (1.0f) / dgSqrt(x))
+#define dgClearFP()			_clearfp() 
+#define dgControlFP(x,y)	_controlfp(x,y)
+
+enum dgSerializeRevisionNumber
+{
+	m_firstRevision = 100,
+	// add new serialization revision number here
+	m_currentRevision 
+};
+
+dgUnsigned64 dgGetTimeInMicrosenconds();
+dgFloat64 dgRoundToFloat(dgFloat64 val);
+void dgSerializeMarker(dgSerialize serializeCallback, void* const userData);
+dgInt32 dgDeserializeMarker(dgDeserialize serializeCallback, void* const userData);
+
+class dgFloatExceptions
+{
+	public:
+	#ifdef _MSC_VER
+		#define DG_FLOAT_EXECTIONS_MASK	(EM_INVALID | EM_DENORMAL | EM_ZERODIVIDE)
+	#else 
+		#define DG_FLOAT_EXECTIONS_MASK	0
+	#endif
+
+	dgFloatExceptions(dgUnsigned32 mask = DG_FLOAT_EXECTIONS_MASK);
+	~dgFloatExceptions();
+
+	private:
+	#if (defined (_MSC_VER) && defined (_WIN_32_VER))
+	dgUnsigned32 m_mask;
+#endif
+};
+
+class dgSetPrecisionDouble 
+{
+	public:
+	dgSetPrecisionDouble();
+	~dgSetPrecisionDouble();
+	#if (defined (_MSC_VER) && defined (_WIN_32_VER))
+	dgInt32 m_mask; 
+	#endif
+};
+
+DG_INLINE dgInt32 dgAtomicExchangeAndAdd (dgInt32* const addend, dgInt32 amount)
+{
+	#if (defined (_WIN_32_VER) || defined (_WIN_64_VER))
+		return _InterlockedExchangeAdd((long*)addend, long(amount));
+	#elif defined (_MSC_VER ) && defined (_ARM_VER) 
+		return _InterlockedExchangeAdd((long*)addend, long(amount));
+	#elif (defined (__MINGW32__) || defined (__MINGW64__))
+		return InterlockedExchangeAdd((long*)addend, long(amount));
+	#elif (defined (_POSIX_VER) || defined (_POSIX_VER_64) ||defined (_MACOSX_VER)|| defined ANDROID)
+		return __sync_fetch_and_add((int32_t*)addend, amount);
+	#else
+		#error "dgAtomicExchangeAndAdd implementation required"
+	#endif
+}
+
+DG_INLINE dgInt32 dgInterlockedExchange(dgInt32* const ptr, dgInt32 value)
+{
+	#if (defined (_WIN_32_VER) || defined (_WIN_64_VER))
+		return _InterlockedExchange((long*) ptr, value);
+	#elif defined (_MSC_VER ) && defined (_ARM_VER) 
+		return _InterlockedExchange((long*)ptr, value);
+	#elif (defined (__MINGW32__) || defined (__MINGW64__))
+		return InterlockedExchange((long*) ptr, value);
+	#elif (defined (_POSIX_VER) || defined (_POSIX_VER_64) ||defined (_MACOSX_VER))
+		return __sync_lock_test_and_set((int32_t*)ptr, value);
+	#else
+		#error "dgInterlockedExchange implementation required"
+	#endif
+}
+
+/*
+DG_INLINE void* dgInterlockedExchange(void** const ptr, void* value)
+{
+	#if defined(_WIN32)
+		#ifdef _M_X64 
+			return (void*)_InterlockedExchange64((dgInt64*)ptr, dgInt64 (value));
+		#else
+			return (void*)_InterlockedExchange((long*) ptr, dgInt32(value));
+		#endif
+	#elif defined(__linux__) || defined(_MACOSX_VER)
+		#if defined(__ANDROID__)
+			return __sync_lock_test_and_set(ptr, value);
+		#else
+			#ifdef __x86_64__
+				return (void*)__sync_lock_test_and_set((dgInt64*) ptr, (dgInt64)value);
+			#else
+				return (void*)__sync_lock_test_and_set((dgInt32*) ptr, (dgInt32) value);
+			#endif
+		#endif
+	#else
+		#error "dgInterlockedExchange implementation required"
+	#endif
+}
+*/
+
+DG_INLINE dgInt32 dgInterlockedTest(dgInt32* const ptr, dgInt32 value)
+{
+	#if (defined (_WIN_32_VER) || defined (_WIN_64_VER))
+		return _InterlockedCompareExchange((long*)ptr, value, value);
+	#elif defined (_MSC_VER ) && defined (_ARM_VER) 
+		return _InterlockedCompareExchange((long*)ptr, value, value);
+	#elif (defined (__MINGW32__) || defined (__MINGW64__))
+		return InterlockedCompareExchange((long*)ptr, value, value);
+	#elif (defined (_POSIX_VER) || defined (_POSIX_VER_64) ||defined (_MACOSX_VER))
+		return __sync_lock_test_and_set((int32_t*)ptr, value);
+	#else
+		#error "dgInterlockedTest implementation required"
+	#endif
+}
+
+DG_INLINE void dgThreadYield()
+{
+#ifndef DG_USE_THREAD_EMULATION
+	std::this_thread::yield();
+#endif
+}
+
+DG_INLINE void dgThreadPause()
+{
+#ifndef DG_USE_THREAD_EMULATION
+	//#if defined (_WIN_32_VER) || defined (_WIN_64_VER) || defined (WIN32) || defined (i386_) || defined (x86_64_)
+	#if defined (_WIN_32_VER) || defined (_WIN_64_VER)
+		_mm_pause();
+	#else 
+		std::this_thread::yield();
+	#endif
+#endif
+}
+
+DG_INLINE void dgSpinLock(dgInt32* const ptr)
+{
+#ifndef DG_USE_THREAD_EMULATION 
+	while (dgInterlockedExchange(ptr, 1)) {
+		dgThreadYield();
+	}
+#endif
+}
+
+
+DG_INLINE void dgSpinUnlock (dgInt32* const ptr)
+{
+	#ifndef DG_USE_THREAD_EMULATION 
+		dgInterlockedExchange(ptr, 0);
+	#else
+		*ptr = 0;
+	#endif
+}
+
+
+class dgScopeSpinLock
+{
+	public:
+	DG_INLINE dgScopeSpinLock(dgInt32* const lock)
+		:m_atomicLock(lock)
+	{
+		dgSpinLock(m_atomicLock);
+	}
+
+	DG_INLINE ~dgScopeSpinLock()
+	{
+		dgSpinUnlock(m_atomicLock);
+	}
+
+	dgInt32* m_atomicLock;
+};
+
+class dgScopeSpinPause
+{
+	public:
+	DG_INLINE dgScopeSpinPause(dgInt32* const lock)
+		:m_atomicLock(lock)
+	{
+		while (dgInterlockedExchange(m_atomicLock, 1)) {
+			dgThreadPause();
+		}
+	}
+
+	DG_INLINE ~dgScopeSpinPause()
+	{
+		dgSpinUnlock(m_atomicLock);
+	}
+
+	dgInt32* m_atomicLock;
+};
+
+
+
+#ifdef _MACOSX_VER
+#include <sys/time.h>
+#ifndef CLOCK_REALTIME
+	#define CLOCK_REALTIME 0
+#endif
+#ifndef CLOCK_MONOTONIC
+	#define CLOCK_MONOTONIC 0
+#endif
+//clock_gettime is not implemented on OSX
+DG_INLINE int clock_gettime(int /*clk_id*/, struct timespec* t) {
+    struct timeval now;
+    int rv = gettimeofday(&now, NULL);
+    if (rv) return rv;
+    t->tv_sec  = now.tv_sec;
+    t->tv_nsec = now.tv_usec * 1000;
+    return 0;
+}
+#endif
+
+#endif
+

thirdparty/include/kiwi/AssocVector.h

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

thirdparty/include/kiwi/constraint.h

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

thirdparty/include/kiwi/debug.h

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

thirdparty/include/kiwi/errors.h

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

thirdparty/include/kiwi/expression.h

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

thirdparty/include/kiwi/kiwi.h

@@ -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"

thirdparty/include/kiwi/maptype.h

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

thirdparty/include/kiwi/row.h

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

thirdparty/include/kiwi/shareddata.h

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

thirdparty/include/kiwi/solver.h

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

thirdparty/include/kiwi/solverimpl.h

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

thirdparty/include/kiwi/strength.h

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

thirdparty/include/kiwi/symbol.h

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

thirdparty/include/kiwi/symbolics.h

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

thirdparty/include/kiwi/term.h

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

thirdparty/include/kiwi/util.h

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

thirdparty/include/kiwi/variable.h

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

thirdparty/include/kiwi/version.h

@@ -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"