Add an alternative solver to lp_solve. The solver based on Ingo's active set solver but is able to handle arbitrary hard and soft constraints. The advantage to lp_solve is that the active set solver can optimize variable in respect to a quadratic objective function. This makes it possible to minimise the quadratic derivation to a desired value e.g. \Sum_i(x_i - x_{i,pref})^2 -> min.

The solver part has been refactored in this way that both solver can be used with the same layout specifications. The active set solver is default now; the performance is not as good as lp_solve, though.



git-svn-id: file:///srv/svn/repos/haiku/haiku/trunk@40285 a95241bf-73f2-0310-859d-f6bbb57e9c96

headers/libs/alm/ALMLayout.h

@@ -124,8 +124,6 @@ private:
 			/*! Add a view without initialize the Area. */
 			BLayoutItem*		_CreateLayoutItem(BView* view);
 
-			void				_SolveLayout();
-
 			void				_UpdateAreaConstraints();
 
 			BSize				_CalculateMinSize();
@@ -152,8 +150,10 @@ private:
 
 			Area*				fCurrentArea;
 
+#if USE_SCALE_VARIABLE
 			Variable*			fScaleWidth;
 			Variable*			fScaleHeight;
+#endif
 };
 
 }	// namespace BALM

headers/libs/alm/Area.h

@@ -20,6 +20,9 @@
 #include "Tab.h"
 
 
+#define USE_SCALE_VARIABLE 1
+
+
 class Constraint;
 
 
@@ -106,6 +109,7 @@ public:
 private:
 								Area(BLayoutItem* item);
 
+#if USE_SCALE_VARIABLE
 			void				_Init(LinearSpec* ls, XTab* left, YTab* top,
 									XTab* right, YTab* bottom,
 									Variable* scaleWidth,
@@ -113,6 +117,11 @@ private:
 			void				_Init(LinearSpec* ls, Row* row, Column* column,
 									Variable* scaleWidth,
 									Variable* scaleHeight);
+#else
+			void				_Init(LinearSpec* ls, XTab* left, YTab* top,
+									XTab* right, YTab* bottom);
+			void				_Init(LinearSpec* ls, Row* row, Column* column);
+#endif
 
 			void				_DoLayout();
 
@@ -152,9 +161,10 @@ private:
 			double				fContentAspectRatio;
 			Constraint*			fContentAspectRatioC;
 
+#if USE_SCALE_VARIABLE
 			Variable*			fScaleWidth;
 			Variable*			fScaleHeight;
-
+#endif
 public:
 	friend class		BALMLayout;
 

headers/libs/linprog/Constraint.h

@@ -27,7 +27,6 @@ class LinearSpec;
  * May render a specification infeasible.
  */
 class Constraint {
-
 public:
 			int32				Index() const;
 
@@ -57,24 +56,26 @@ public:
 			const char*			Label();
 			void				SetLabel(const char* label);
 
-			void				WriteXML(BFile* file);
-
 			Variable*			DNeg() const;
 			Variable*			DPos() const;
 
+			bool				IsSoft() const;
+
 			bool				IsValid();
 			void				Invalidate();
 
 								operator BString() const;
 			void				GetString(BString& string) const;
+			void				PrintToStream();
 
 								~Constraint();
 
 protected:
 								Constraint(LinearSpec* ls,
 									SummandList* summands, OperatorType op,
-									double rightSide, double penaltyNeg,
+									double rightSide,
-									double penaltyPos);
+									double penaltyNeg = -1,
+									double penaltyPos = -1);
 
 private:
 			LinearSpec*			fLS;
@@ -92,6 +93,7 @@ private:
 
 public:
 	friend class		LinearSpec;
+	friend class		LPSolveInterface;
 
 };
 

headers/libs/linprog/LinearSpec.h

@@ -11,47 +11,51 @@
 
 #include <List.h>
 #include <OS.h>
+#include <Size.h>
 #include <String.h>
 #include <SupportDefs.h>
 
 #include "Constraint.h"
 #include "LinearProgrammingTypes.h"
-#include "PenaltyFunction.h"
 #include "Summand.h"
 #include "Variable.h"
 
 
 namespace LinearProgramming {
 
+
+class LinearSpec;
+
+
+const BSize kMinSize(0, 0);
+const BSize kMaxSize(B_SIZE_UNLIMITED, B_SIZE_UNLIMITED);
+
+
 class SolverInterface {
 public:
+								SolverInterface(LinearSpec* linSpec);
+
 	virtual						~SolverInterface() {}
 
-	virtual ResultType			Solve(VariableList& variables) = 0;
+	virtual ResultType			Solve() = 0;
-	virtual double				GetObjectiveValue() = 0;
 
-	virtual	bool				AddVariable() = 0;
+	virtual	bool				VariableAdded(Variable* variable) = 0;
-	virtual	bool				RemoveVariable(int variable) = 0;
+	virtual	bool				VariableRemoved(Variable* variable) = 0;
-	virtual	bool				SetVariableRange(int variable, double min,
+	virtual	bool				VariableRangeChanged(Variable* variable) = 0;
-									double max) = 0;
 
-	virtual	bool				AddConstraint(int nElements,
+	virtual	bool				ConstraintAdded(Constraint* constraint) = 0;
-									double* coefficients, int* variableIndices,
+	virtual	bool				ConstraintRemoved(Constraint* constraint) = 0;
-									OperatorType op, double rightSide) = 0;
+	virtual	bool				LeftSideChanged(Constraint* constraint) = 0;
-	virtual	bool				RemoveConstraint(int constraint) = 0;
+	virtual	bool				RightSideChanged(Constraint* constraint) = 0;
-	virtual	bool				SetLeftSide(int constraint, int nElements,
+	virtual	bool				OperatorChanged(Constraint* constraint) = 0;
-									double* coefficients,
-									int* variableIndices) = 0;
-	virtual	bool				SetRightSide(int constraint, double value) = 0;
-	virtual	bool				SetOperator(int constraint,
-									OperatorType op) = 0;
-
-	virtual bool				SetObjectiveFunction(int nElements,
-									double* coefficients,
-									int* variableIndices) = 0;
-	virtual bool				SetOptimization(OptimizationType value) = 0;
 
 	virtual bool				SaveModel(const char* fileName) = 0;
+
+	virtual	BSize				MinSize(Variable* width, Variable* height) = 0;
+	virtual	BSize				MaxSize(Variable* width, Variable* height) = 0;
+
+protected:
+			LinearSpec*			fLinearSpec; 
 };
 
 
@@ -113,31 +117,22 @@ public:
 									OperatorType op, double rightSide,
 									double penaltyNeg, double penaltyPos);
 
-			PenaltyFunction*	AddPenaltyFunction(Variable* var, BList* xs,
+			BSize				MinSize(Variable* width, Variable* height);
-									BList* gs);
+			BSize				MaxSize(Variable* width, Variable* height);
-
-			SummandList*		ObjectiveFunction();
-			//! Caller takes ownership of the Summand's and the SummandList.
-			SummandList*		SwapObjectiveFunction(
-									SummandList* objFunction);
-			void				SetObjectiveFunction(SummandList* objFunction);
-			void				UpdateObjectiveFunction();
 
 			ResultType			Solve();
 			bool				Save(const char* fileName);
 
 			int32				CountColumns() const;
-			OptimizationType	Optimization() const;
-			void				SetOptimization(OptimizationType value);
 
 			ResultType			Result() const;
-			double				ObjectiveValue() const;
+			bigtime_t			SolvingTime() const;
-			double				SolvingTime() const;
 
 			operator BString() const;
 			void				GetString(BString& string) const;
 
 	const	ConstraintList&		Constraints() const;
+	const	VariableList&		Variables() const;
 
 protected:
 	friend class Constraint;
@@ -152,14 +147,10 @@ private:
 									OperatorType op, double rightSide,
 									double penaltyNeg, double penaltyPos);
 
-			OptimizationType	fOptimization;
-
-			SummandList*		fObjFunction;
 			VariableList		fVariables;
 			ConstraintList		fConstraints;
 			ResultType			fResult;
-			double 				fObjectiveValue;
+			bigtime_t 			fSolvingTime;
-			double 				fSolvingTime;
 
 			SolverInterface*	fSolver;
 };

headers/libs/linprog/PenaltyFunction.h

@@ -1,50 +0,0 @@
-/*
- * Copyright 2007-2008, Christof Lutteroth, [email protected]
- * Copyright 2007-2008, James Kim, [email protected]
- * Distributed under the terms of the MIT License.
- */
-
-#ifndef	PENALTY_FUNCTION_H
-#define	PENALTY_FUNCTION_H
-
-#include <List.h>
-
-
-namespace LinearProgramming {
-
-class LinearSpec;
-class Variable;
-
-/**
- * Penalty function.
- */
-class PenaltyFunction {
-
-protected:
-						PenaltyFunction(LinearSpec* ls, Variable* var, BList* xs, BList* gs);	
-
-public:
-						~PenaltyFunction();
-	const Variable*		Var() const;
-	const BList*		Xs() const;
-	const BList*		Gs() const;
-
-private:
-	LinearSpec*			fLS;
-	Variable*			fVar;
-	BList*				fXs;		// double
-	BList*				fGs;		// double
-	BList*				fConstraints;
-	BList*				fObjFunctionSummands;
-
-public:
-	friend class			LinearSpec;
-
-};
-
-}	// namespace LinearProgramming
-
-using LinearProgramming::PenaltyFunction;
-
-#endif	// PENALTY_FUNCTION_H
-

headers/libs/linprog/Summand.h

@@ -27,6 +27,7 @@ public:
 			Variable*			Var();
 			void				SetVar(Variable* var);
 
+			int32				VariableIndex();
 private:
 			double				fCoeff;
 			Variable*			fVar;

src/libs/alm/ALMLayout.cpp

@@ -19,8 +19,6 @@ using namespace LinearProgramming;
 
 
 const BSize kUnsetSize(B_SIZE_UNSET, B_SIZE_UNSET);
-const BSize kMinSize(0, 0);
-const BSize kMaxSize(B_SIZE_UNLIMITED, B_SIZE_UNLIMITED);
 
 
 /*!
@@ -42,7 +40,8 @@ BALMLayout::BALMLayout(float spacing, BALMLayout* friendLayout)
 	fTop = AddYTab();
 	fBottom = AddYTab();
 
-	// the Left tab is always at x-position 0, and the Top tab is always at y-position 0
+	// the Left tab is always at x-position 0, and the Top tab is always at
+	// y-position 0
 	fLeft->SetRange(0, 0);
 	fTop->SetRange(0, 0);
 
@@ -54,15 +53,19 @@ BALMLayout::BALMLayout(float spacing, BALMLayout* friendLayout)
 
 	fPerformancePath = NULL;
 
+#if USE_SCALE_VARIABLE
 	fScaleWidth = fSolver->AddVariable();
 	fScaleHeight = fSolver->AddVariable();
+#endif
 }
 
 
 BALMLayout::~BALMLayout()
 {
+#if USE_SCALE_VARIABLE
 	delete fScaleWidth;
 	delete fScaleHeight;
+#endif
 }
 
 
@@ -484,7 +487,11 @@ BALMLayout::AddItem(BLayoutItem* item, XTab* left, YTab* top, XTab* right,
 		return NULL;
 	fCurrentArea = area;
 
+#if USE_SCALE_VARIABLE
 	area->_Init(fSolver, left, top, right, bottom, fScaleWidth, fScaleHeight);
+#else
+	area->_Init(fSolver, left, top, right, bottom);
+#endif
 	return area;
 }
 
@@ -499,7 +506,11 @@ BALMLayout::AddItem(BLayoutItem* item, Row* row, Column* column)
 		return NULL;
 	fCurrentArea = area;
 
+#if USE_SCALE_VARIABLE
 	area->_Init(fSolver, row, column, fScaleWidth, fScaleHeight);
+#else
+	area->_Init(fSolver, row, column);
+#endif
 	return area;
 }
 
@@ -701,7 +712,7 @@ BALMLayout::DerivedLayoutItems()
 	Right()->SetRange(area.right, area.right);
 	Bottom()->SetRange(area.bottom, area.bottom);
 
-	_SolveLayout();
+	fSolver->Solve();
 
 	// if new layout is infeasible, use previous layout
 	if (fSolver->Result() == kInfeasible)
@@ -786,36 +797,6 @@ BALMLayout::_CreateLayoutItem(BView* view)
 }
 
 
-void
-BALMLayout::_SolveLayout()
-{
-	// Try to solve the layout until the result is kOptimal or kInfeasible,
-	// maximally 15 tries sometimes the solving algorithm encounters numerical
-	// problems (NUMFAILURE), and repeating the solving often helps to overcome
-	// them.
-	BFile* file = NULL;
-	if (fPerformancePath != NULL) {
-		file = new(std::nothrow) BFile(fPerformancePath,
-			B_READ_WRITE | B_CREATE_FILE | B_OPEN_AT_END);
-	}
-
-	ResultType result;
-	for (int32 tries = 0; tries < 15; tries++) {
-		result = fSolver->Solve();
-		if (fPerformancePath != NULL) {
-			/*char buffer [100];
-			file->Write(buffer, sprintf(buffer, "%d\t%fms\t#vars=%ld\t"
-				"#constraints=%ld\n", result, fSolver->SolvingTime(),
-				fSolver->Variables()->CountItems(),
-				fSolver->Constraints()->CountItems()));*/
-		}
-		if (result == kOptimal || result == kInfeasible)
-			break;
-	}
-	delete file;
-}
-
-
 /**
  * Caculates the miminum size.
  */
@@ -824,24 +805,7 @@ BALMLayout::_CalculateMinSize()
 {
 	_UpdateAreaConstraints();
 
-	SummandList* newObjFunction = new(std::nothrow) SummandList(2);
+	return fSolver->MinSize(Right(), Bottom());
-	newObjFunction->AddItem(new(std::nothrow) Summand(1.0, fRight));
-	newObjFunction->AddItem(new(std::nothrow) Summand(1.0, fBottom));
-	SummandList* oldObjFunction = fSolver->SwapObjectiveFunction(
-		newObjFunction);
-	_SolveLayout();
-	fSolver->SetObjectiveFunction(oldObjFunction);
-
-	if (fSolver->Result() == kUnbounded)
-		return kMinSize;
-	if (fSolver->Result() != kOptimal) {
-		fSolver->Save("failed-layout.txt");
-		printf("Could not solve the layout specification (%d). "
-			"Saved specification in file failed-layout.txt", fSolver->Result());
-	}
-
-	return BSize(Right()->Value() - Left()->Value(),
-		Bottom()->Value() - Top()->Value());
 }
 
 
@@ -853,24 +817,7 @@ BALMLayout::_CalculateMaxSize()
 {
 	_UpdateAreaConstraints();
 
-	SummandList* newObjFunction = new(std::nothrow) SummandList(2);
+	return fSolver->MaxSize(Right(), Bottom());
-	newObjFunction->AddItem(new(std::nothrow) Summand(-1.0, fRight));
-	newObjFunction->AddItem(new(std::nothrow) Summand(-1.0, fBottom));
-	SummandList* oldObjFunction = fSolver->SwapObjectiveFunction(
-		newObjFunction);
-	_SolveLayout();
-	fSolver->SetObjectiveFunction(oldObjFunction);
-
-	if (fSolver->Result() == kUnbounded)
-		return kMaxSize;
-	if (fSolver->Result() != kOptimal) {
-		fSolver->Save("failed-layout.txt");
-		printf("Could not solve the layout specification (%d). "
-			"Saved specification in file failed-layout.txt", fSolver->Result());
-	}
-
-	return BSize(Right()->Value() - Left()->Value(),
-		Bottom()->Value() - Top()->Value());
 }
 
 
@@ -882,7 +829,7 @@ BALMLayout::_CalculatePreferredSize()
 {
 	_UpdateAreaConstraints();
 
-	_SolveLayout();
+	fSolver->Solve();
 	if (fSolver->Result() != kOptimal) {
 		fSolver->Save("failed-layout.txt");
 		printf("Could not solve the layout specification (%d). "

src/libs/alm/Area.cpp

@@ -599,19 +599,24 @@ Area::Area(BLayoutItem* item)
 /**
  * Initialize variables.
  */
+#if USE_SCALE_VARIABLE
 void
 Area::_Init(LinearSpec* ls, XTab* left, YTab* top, XTab* right, YTab* bottom,
 	Variable* scaleWidth, Variable* scaleHeight)
 {
+	fScaleWidth = scaleWidth;
+	fScaleHeight = scaleHeight;
+#else
+void
+Area::_Init(LinearSpec* ls, XTab* left, YTab* top, XTab* right, YTab* bottom)
+{
+#endif
 	fLS = ls;
 	fLeft = left;
 	fRight = right;
 	fTop = top;
 	fBottom = bottom;
 
-	fScaleWidth = scaleWidth;
-	fScaleHeight = scaleHeight;
-
 	// adds the two essential constraints of the area that make sure that the
 	// left x-tab is really to the left of the right x-tab, and the top y-tab
 	// really above the bottom y-tab
@@ -621,6 +626,7 @@ Area::_Init(LinearSpec* ls, XTab* left, YTab* top, XTab* right, YTab* bottom,
 	fConstraints.AddItem(fMinContentWidth);
 	fConstraints.AddItem(fMinContentHeight);
 
+#if USE_SCALE_VARIABLE
 	fPreferredContentWidth = fLS->AddConstraint(-1.0, fLeft, 1.0, fRight, -1.0,
 		fScaleWidth, kEQ, 0, fShrinkPenalties.Width(),
 		fGrowPenalties.Width());
@@ -628,18 +634,34 @@ Area::_Init(LinearSpec* ls, XTab* left, YTab* top, XTab* right, YTab* bottom,
 	fPreferredContentHeight = fLS->AddConstraint(-1.0, fTop, 1.0, fBottom, -1.0,
 		fScaleHeight, kEQ, 0, fShrinkPenalties.Height(),
 		fGrowPenalties.Height());
+#else
+	BSize preferredSize = fLayoutItem->PreferredSize();
+	fPreferredContentWidth = fLS->AddConstraint(-1.0, fLeft, 1.0, fRight, kEQ,
+		0, fShrinkPenalties.Width(), fGrowPenalties.Width());
+	_UpdatePreferredWidthConstraint(preferredSize);
+	fPreferredContentHeight = fLS->AddConstraint(-1.0, fTop, 1.0, fBottom, kEQ,
+		0, fShrinkPenalties.Height(), fGrowPenalties.Height());
+	_UpdatePreferredHeightConstraint(preferredSize);
+#endif
 
 	fConstraints.AddItem(fPreferredContentWidth);
 	fConstraints.AddItem(fPreferredContentHeight);
 }
 
 
+#if USE_SCALE_VARIABLE
 void
 Area::_Init(LinearSpec* ls, Row* row, Column* column, Variable* scaleWidth,
 	Variable* scaleHeight)
 {
 	_Init(ls, column->Left(), row->Top(), column->Right(), row->Bottom(),
 		scaleWidth, scaleHeight);
+#else
+void
+Area::_Init(LinearSpec* ls, Row* row, Column* column)
+{
+	_Init(ls, column->Left(), row->Top(), column->Right(), row->Bottom());
+#endif
 	fRow = row;
 	fColumn = column;
 }
@@ -723,22 +745,28 @@ Area::_UpdateMaxSizeConstraint(BSize max)
 void
 Area::_UpdatePreferredWidthConstraint(BSize& preferred)
 {
-	float width = 32000;
+	float width = 0;
 	if (preferred.width > 0)
 		width = preferred.Width() + LeftInset() + RightInset();
-	
+#if USE_SCALE_VARIABLE
 	fPreferredContentWidth->SetLeftSide(-1.0, fLeft, 1.0, fRight, -width,
 		fScaleWidth);
+#else
+	fPreferredContentWidth->SetRightSide(width);
+#endif
 }
 
 
 void
 Area::_UpdatePreferredHeightConstraint(BSize& preferred)
 {
-	float height = 32000;
+	float height = 0;
 	if (preferred.height > 0)
 		height = preferred.Height() + TopInset() + BottomInset();
-
+#if USE_SCALE_VARIABLE
 	fPreferredContentHeight->SetLeftSide(-1.0, fTop, 1.0, fBottom, -height,
 		fScaleHeight);
+#else
+	fPreferredContentHeight->SetRightSide(height);
+#endif
 }

src/libs/linprog/ActiveSetSolver.cpp

@@ -0,0 +1,542 @@
+#include "ActiveSetSolver.h"
+
+#include <stdio.h>
+
+#include "LayoutOptimizer.h"
+
+
+//#define DEBUG_ACTIVE_SOLVER
+
+#ifdef DEBUG_ACTIVE_SOLVER
+#include <stdio.h>
+#define TRACE(x...) printf(x)
+#else
+#define TRACE(x...) /* nothing */
+#endif
+
+
+using namespace LinearProgramming;
+using namespace BPrivate::Layout;
+
+
+template<typename Type>
+static inline void
+swap(Type& a, Type& b)
+{
+	Type c = a;
+	a = b;
+	b = c;
+}
+
+
+EquationSystem::EquationSystem(int32 rows, int32 columns)
+	:
+	fRows(rows),
+	fColumns(columns)
+{
+	fMatrix = allocate_matrix(fRows, fColumns);
+	fB = new double[fColumns];
+	// better init all values to prevent side cases where not all variables
+	// needed to solve the problem, coping theses values to the results could
+	// cause problems
+	for (int i = 0; i < fColumns; i++)
+		fB[i] = 0;
+	zero_matrix(fMatrix, fRows, fColumns);
+
+	fRowIndices = new int32[fRows];
+	fColumnIndices = new int32[fColumns];
+	for (int i = 0; i < fRows; i++)
+		fRowIndices[i] = i;
+	for (int i = 0; i < fColumns; i++)
+		fColumnIndices[i] = i;
+}
+
+
+EquationSystem::~EquationSystem()
+{
+	free_matrix(fMatrix);
+	delete[] fB;
+	delete[] fRowIndices;
+	delete[] fColumnIndices;
+}
+
+
+void
+EquationSystem::SetRows(int32 rows)
+{
+	fRows = rows;
+}
+
+
+int32
+EquationSystem::Rows()
+{
+	return fRows;
+}
+
+
+int32
+EquationSystem::Columns()
+{
+	return fColumns;
+}
+
+
+double&
+EquationSystem::A(int32 row, int32 column)
+{
+	return fMatrix[fRowIndices[row]][fColumnIndices[column]];
+}
+
+
+double&
+EquationSystem::B(int32 row)
+{
+	return fB[row];
+}
+
+
+void
+EquationSystem::Results(double* results, int32 size)
+{
+	for (int i = 0; i < size; i++)
+		results[i] = 0;
+	for (int i = 0; i < fColumns; i++) {
+		int32 index = fColumnIndices[i];
+		if (index < fRows)
+			results[index] = fB[i];
+	}
+}
+
+
+void
+EquationSystem::SwapColumn(int32 i, int32 j)
+{
+	swap(fColumnIndices[i], fColumnIndices[j]);
+}
+
+
+void
+EquationSystem::SwapRow(int32 i, int32 j)
+{
+	swap(fRowIndices[i], fRowIndices[j]);
+	swap(fB[i], fB[j]);
+}
+
+
+bool
+EquationSystem::GaussJordan()
+{
+	// basic solve
+	for (int i = 0; i < fRows; i++) {
+		// find none zero element
+		int swapRow = -1;
+		for (int r = i; r < fRows; r++) {
+			double& value = fMatrix[fRowIndices[r]][fColumnIndices[i]];
+			if (fuzzy_equals(value, 0))
+				continue;
+			swapRow = r;
+			break;
+		}
+		if (swapRow == -1) {
+			int swapColumn = -1;
+			for (int c = i + 1; c < fColumns; c++) {
+				double& value = fMatrix[fRowIndices[i]][fColumnIndices[c]];
+				if (fuzzy_equals(value, 0))
+					continue;
+				swapRow = i;
+				swapColumn = c;
+				break;
+			}
+			if (swapColumn == -1) {
+				printf("can't solve column %i\n", i);
+				return false;
+			}
+			SwapColumn(i, swapColumn);
+		}
+		if (i != swapRow)
+			SwapRow(i, swapRow);
+
+		// normalize
+		GaussJordan(i);
+	}
+	return true;
+}
+
+
+void
+EquationSystem::GaussJordan(int32 i)
+{
+	double value = fMatrix[fRowIndices[i]][fColumnIndices[i]];
+	for (int j = 0; j < fColumns; j++)
+		fMatrix[fRowIndices[i]][fColumnIndices[j]] /= value;
+	fB[i] /= value;
+
+	for (int r = 0; r < fRows; r++) {
+		if (r == i)
+			continue;
+		double q = -fMatrix[fRowIndices[r]][fColumnIndices[i]];
+		// don't need to do nothing, since matrix is typically sparse this
+		// should save some work
+		if (fuzzy_equals(q, 0))
+			continue;
+		for (int c = 0; c < fColumns; c++)
+			fMatrix[fRowIndices[r]][fColumnIndices[c]]
+				+= fMatrix[fRowIndices[i]][fColumnIndices[c]] * q;
+
+		fB[r] += fB[i] * q;
+	}
+}
+
+
+void
+EquationSystem::RemoveLinearlyDependentRows()
+{
+	double oldB[fRows];
+	for (int r = 0; r < fRows; r++)
+		oldB[r] = fB[r];
+
+	double** temp = allocate_matrix(fRows, fColumns);
+	bool independentRows[fRows];
+
+	// copy to temp
+	copy_matrix(fMatrix, temp, fRows, fColumns);
+	int nIndependent = compute_dependencies(temp, fRows, fColumns,
+		independentRows);
+	if (nIndependent == fRows)
+		return;
+
+	// remove the rows
+	for (int i = 0; i < fRows; i++) {
+		if (!independentRows[i]) {
+			int lastDepRow = -1;
+			for (int d = fRows - 1; d > i; d--) {
+				if (independentRows[d]) {
+					lastDepRow = d;
+					break;
+				}
+			}
+			if (lastDepRow < 0)
+				break;
+			SwapRow(i, lastDepRow);
+			fRows--;
+		}
+	}
+	fRows = nIndependent;
+
+	free_matrix(temp);
+}
+
+
+void
+EquationSystem::RemoveUnusedVariables()
+{
+	for (int c = 0; c < fColumns; c++) {
+		bool used = false;
+		for (int r = 0; r < fRows; r++) {
+			if (!fuzzy_equals(fMatrix[r][fColumnIndices[c]], 0)) {
+				used = true;
+				break;
+			}
+		}
+		if (used)
+			continue;
+
+		//MoveColumnRight(c, fColumns - 1);
+		SwapColumn(c, fColumns - 1);
+		fColumns--;
+		c--;
+	}
+}
+
+
+void
+EquationSystem::MoveColumnRight(int32 i, int32 target)
+{
+	int32 index = fColumnIndices[i];
+	for (int c = i; c < target; c++)
+		fColumnIndices[c] = fColumnIndices[c + 1];
+	fColumnIndices[target] = index;
+}
+
+
+void
+EquationSystem::Print()
+{
+	for (int m = 0; m < fRows; m++) {
+		for (int n = 0; n < fColumns; n++)
+			printf("%.1f ", fMatrix[fRowIndices[m]][fColumnIndices[n]]);
+		printf("= %.1f\n", fB[m]);
+	}
+}
+
+
+ActiveSetSolver::ActiveSetSolver(LinearSpec* linearSpec)
+	:
+	SolverInterface(linearSpec),
+
+	fVariables(linearSpec->Variables()),
+	fConstraints(linearSpec->Constraints())
+{
+
+}
+
+
+ActiveSetSolver::~ActiveSetSolver()
+{
+
+}
+
+
+/* Using algorithm found in:
+Solving Inequalities and Proving Farkas's Lemma Made Easy
+David Avis and Bohdan Kaluzny
+The American Mathematical Monthly
+Vol. 111, No. 2 (Feb., 2004), pp. 152-157 */
+bool
+solve(EquationSystem& system)
+{
+	// basic solve
+	if (!system.GaussJordan())
+		return false;
+
+	bool done = false;
+	while (!done) {
+		double smallestB = HUGE_VALF;
+		int smallestBRow = -1;
+		for (int row = 0; row < system.Rows(); row++) {
+			if (system.B(row) > 0 || fuzzy_equals(system.B(row), 0))
+				continue;
+
+			double bValue = fabs(system.B(row));
+			if (bValue < smallestB) {
+				smallestB = bValue;
+				smallestBRow = row;
+			}
+		}
+		if (smallestBRow == -1) {
+			done = true;
+			break;
+		}
+
+		int negValueCol = -1;
+		for (int col = system.Rows(); col < system.Columns(); col++) {
+			double value = system.A(smallestBRow, col);
+			if (value > 0 || fuzzy_equals(value, 0))
+				continue;
+			negValueCol = col;
+			break;
+		}
+		if (negValueCol == -1) {
+			printf("can't solve\n");
+			return false;
+		}
+
+		system.SwapColumn(smallestBRow, negValueCol);
+
+		// eliminate
+		system.GaussJordan(smallestBRow);
+	}
+
+	return true;
+}
+
+
+ResultType
+ActiveSetSolver::Solve()
+{
+	int32 nConstraints = fConstraints.CountItems();
+	int32 nVariables = fVariables.CountItems();
+
+	if (nVariables > nConstraints) {
+		printf("More variables then constraints! vars: %i, constraints: %i\n",
+			(int)nVariables, (int)nConstraints);
+		return kInfeasible;
+	}
+
+	/* First find an initial solution and then optimize it using the active set
+	method. */
+	EquationSystem system(nConstraints, nVariables + nConstraints);
+
+	int32 slackIndex = nVariables;
+	// setup constraint matrix and add slack variables if necessary
+	int32 rowIndex = 0;
+	for (int32 c = 0; c < nConstraints; c++) {
+		Constraint* constraint = fConstraints.ItemAt(c);
+		if (constraint->IsSoft())
+			continue;
+		SummandList* leftSide = constraint->LeftSide();
+		system.B(rowIndex) = constraint->RightSide();
+		for (int32 sIndex = 0; sIndex < leftSide->CountItems(); sIndex++ ) {
+			Summand* summand = leftSide->ItemAt(sIndex);
+			int32 coefficient = summand->Coeff();
+			system.A(rowIndex, summand->VariableIndex()) = coefficient;
+		}
+		if (constraint->Op() == kLE) {
+			system.A(rowIndex, slackIndex) = 1;
+			slackIndex++;
+		} else if (constraint->Op() == kGE) {
+			system.A(rowIndex, slackIndex) = -1;
+			slackIndex++;
+		}
+		rowIndex++;
+	}
+
+	system.SetRows(rowIndex);
+
+	system.RemoveLinearlyDependentRows();
+	system.RemoveUnusedVariables();
+
+	if (!solve(system))
+		return kInfeasible;
+
+	double results[nVariables + nConstraints];
+	system.Results(results, nVariables + nConstraints);
+	printf("base system solved\n");
+
+	LayoutOptimizer optimizer(fConstraints, nVariables);
+	optimizer.Solve(results);
+
+	// back to the variables
+	for (int32 i = 0; i < nVariables; i++)
+		fVariables.ItemAt(i)->SetValue(results[i]);
+
+	for (int32 i = 0; i < nVariables; i++)
+		TRACE("var %f\n", results[i]);
+
+	return kOptimal;
+}
+
+
+bool
+ActiveSetSolver::VariableAdded(Variable* variable)
+{
+	// TODO: error checks
+	fVariableGEConstraints.AddItem(NULL);
+	fVariableLEConstraints.AddItem(NULL);
+	return true;
+}
+
+
+bool
+ActiveSetSolver::VariableRemoved(Variable* variable)
+{
+	fVariableGEConstraints.RemoveItemAt(variable->Index());
+	fVariableLEConstraints.RemoveItemAt(variable->Index());
+	return true;
+}
+
+
+bool
+ActiveSetSolver::VariableRangeChanged(Variable* variable)
+{
+	double min = variable->Min();
+	double max = variable->Max();
+	int32 variableIndex = variable->Index();
+
+	Constraint* constraintGE = fVariableGEConstraints.ItemAt(variableIndex);
+	Constraint* constraintLE = fVariableLEConstraints.ItemAt(variableIndex);
+	if (constraintGE == NULL && min > -20000) {
+		constraintGE = fLinearSpec->AddConstraint(1, variable, kGE, 0);
+		if (constraintGE == NULL)
+			return false;
+		fVariableGEConstraints.RemoveItemAt(variableIndex);
+		fVariableGEConstraints.AddItem(constraintGE, variableIndex);
+	}
+	if (constraintLE == NULL && max < 20000) {
+		constraintLE = fLinearSpec->AddConstraint(1, variable, kLE, 20000);
+		if (constraintLE == NULL)
+			return false;
+		fVariableLEConstraints.RemoveItemAt(variableIndex);
+		fVariableLEConstraints.AddItem(constraintLE, variableIndex);
+	}
+
+	if (constraintGE)
+		constraintGE->SetRightSide(min);
+	if (constraintLE)
+		constraintLE->SetRightSide(max);
+	return true;
+}
+
+
+bool
+ActiveSetSolver::ConstraintAdded(Constraint* constraint)
+{
+	return true;
+}
+
+
+bool
+ActiveSetSolver::ConstraintRemoved(Constraint* constraint)
+{
+	return true;
+}
+
+
+bool
+ActiveSetSolver::LeftSideChanged(Constraint* constraint)
+{
+	return true;
+}
+
+
+bool
+ActiveSetSolver::RightSideChanged(Constraint* constraint)
+{
+	return true;
+}
+
+
+bool
+ActiveSetSolver::OperatorChanged(Constraint* constraint)
+{
+	return true;
+}
+
+
+bool
+ActiveSetSolver::SaveModel(const char* fileName)
+{
+	return false;
+}
+
+
+BSize
+ActiveSetSolver::MinSize(Variable* width, Variable* height)
+{
+	Constraint* heightConstraint = fLinearSpec->AddConstraint(1, height,
+		kEQ, 0, 5, 5);
+	Constraint* widthConstraint = fLinearSpec->AddConstraint(1, width,
+		kEQ, 0, 5, 5);
+	ResultType result = Solve();
+	fLinearSpec->RemoveConstraint(heightConstraint);
+	fLinearSpec->RemoveConstraint(widthConstraint);
+
+	if (result == kUnbounded)
+		return kMinSize;
+	if (result != kOptimal)
+		printf("Could not solve the layout specification (%d). ", result);
+
+	return BSize(width->Value(), height->Value());
+}
+
+
+BSize
+ActiveSetSolver::MaxSize(Variable* width, Variable* height)
+{
+	const double kHugeValue = 32000;
+	Constraint* heightConstraint = fLinearSpec->AddConstraint(1, height,
+		kEQ, kHugeValue, 5, 5);
+	Constraint* widthConstraint = fLinearSpec->AddConstraint(1, width,
+		kEQ, kHugeValue, 5, 5);
+	ResultType result = Solve();
+	fLinearSpec->RemoveConstraint(heightConstraint);
+	fLinearSpec->RemoveConstraint(widthConstraint);
+
+	if (result == kUnbounded)
+		return kMinSize;
+	if (result != kOptimal)
+		printf("Could not solve the layout specification (%d). ", result);
+
+	return BSize(width->Value(), height->Value());
+}

src/libs/linprog/ActiveSetSolver.h

@@ -0,0 +1,82 @@
+/*
+ * Copyright 2010, Clemens Zeidler <[email protected]>
+ * Distributed under the terms of the MIT License.
+ */
+#ifndef ACTICE_SET_SOLVER_H
+#define ACTICE_SET_SOLVER_H
+
+
+#include "LinearSpec.h"
+
+
+class EquationSystem {
+public:
+								EquationSystem(int32 rows, int32 columns);
+								~EquationSystem();
+
+			void				SetRows(int32 rows);
+			int32				Rows();
+			int32				Columns();
+
+	inline	double&				A(int32 row, int32 column);
+	inline	double&				B(int32 row);
+			/*! Copy the solved variables into results, keeping the original
+			variable order. */
+	inline	void				Results(double* results, int32 size);
+
+	inline	void				SwapColumn(int32 i, int32 j);
+	inline	void				SwapRow(int32 i, int32 j);
+
+			bool				GaussJordan();
+			/*! Gauss Jordan elimination just for one column, the diagonal
+			element must be none zero. */
+			void				GaussJordan(int32 column);
+
+			void				RemoveLinearlyDependentRows();
+			void				RemoveUnusedVariables();
+
+			void				MoveColumnRight(int32 i, int32 target);
+
+			void				Print();
+private:
+			int32*				fRowIndices;
+			int32*				fColumnIndices;
+			double**			fMatrix;
+			double*				fB;
+			int32				fRows;
+			int32				fColumns;
+};
+
+
+class ActiveSetSolver : public LinearProgramming::SolverInterface {
+public:
+								ActiveSetSolver(LinearSpec* linearSpec);
+								~ActiveSetSolver();
+
+			ResultType			Solve();
+
+			bool				VariableAdded(Variable* variable);
+			bool				VariableRemoved(Variable* variable);
+			bool				VariableRangeChanged(Variable* variable);
+
+			bool				ConstraintAdded(Constraint* constraint);
+			bool				ConstraintRemoved(Constraint* constraint);
+			bool				LeftSideChanged(Constraint* constraint);
+			bool				RightSideChanged(Constraint* constraint);
+			bool				OperatorChanged(Constraint* constraint);
+
+			bool				SaveModel(const char* fileName);
+
+			BSize				MinSize(Variable* width, Variable* height);
+			BSize				MaxSize(Variable* width, Variable* height);
+
+public:
+	const	VariableList&		fVariables;
+	const	ConstraintList&		fConstraints;
+
+			ConstraintList		fVariableGEConstraints;
+			ConstraintList		fVariableLEConstraints;
+};
+
+
+#endif // ACTICE_SET_SOLVER_H

src/libs/linprog/Constraint.cpp

@@ -34,11 +34,10 @@ int32
 Constraint::Index() const
 {
 	int32 i = fLS->Constraints().IndexOf(this);
-	if (i == -1) {
+	if (i == -1)
 		STRACE(("Constraint not part of fLS->Constraints()."));
-		return -1;
+
-	}
+	return i;
-	return i + 1;
 }
 
 
@@ -222,22 +221,7 @@ Constraint::SetPenaltyNeg(double value)
 {
 	fPenaltyNeg = value;
 
-	if (!fIsValid)
+	fLS->UpdateLeftSide(this);
-		return;
-
-	if (fDNegObjSummand == NULL) {
-		fDNegObjSummand = new(std::nothrow) Summand(value, fLS->AddVariable());
-		fLS->ObjectiveFunction()->AddItem(fDNegObjSummand);
-		fLS->UpdateLeftSide(this);
-		fLS->UpdateObjectiveFunction();
-		return;
-	}
-
-	if (value == fDNegObjSummand->Coeff())
-		return;
-
-	fDNegObjSummand->SetCoeff(value);
-	fLS->UpdateObjectiveFunction();
 }
 
 
@@ -263,22 +247,7 @@ Constraint::SetPenaltyPos(double value)
 {
 	fPenaltyPos = value;
 
-	if (!fIsValid)
+	fLS->UpdateLeftSide(this);
-		return;
-
-	if (fDPosObjSummand == NULL) {
-		fDPosObjSummand = new(std::nothrow) Summand(value, fLS->AddVariable());
-		fLS->ObjectiveFunction()->AddItem(fDPosObjSummand);
-		fLS->UpdateLeftSide(this);
-		fLS->UpdateObjectiveFunction();
-		return;
-	}
-
-	if (value == fDPosObjSummand->Coeff())
-		return;
-
-	fDPosObjSummand->SetCoeff(value);
-	fLS->UpdateObjectiveFunction();
 }
 
 
@@ -296,47 +265,6 @@ Constraint::SetLabel(const char* label)
 }
 
 
-void
-Constraint::WriteXML(BFile* file)
-{
-	if (!file->IsWritable())
-		return;
-
-	char buffer[200];
-
-	file->Write(buffer, sprintf(buffer, "\t<constraint>\n"));
-	file->Write(buffer, sprintf(buffer, "\t\t<leftside>\n"));
-
-	Summand* summand;
-	for (int32 i = 0; i < fLeftSide->CountItems(); i++) {
-		summand = (Summand*)fLeftSide->ItemAt(i);
-		file->Write(buffer, sprintf(buffer, "\t\t\t<summand>\n"));
-		file->Write(buffer, sprintf(buffer, "\t\t\t\t<coeff>%f</coeff>\n",
-			summand->Coeff()));
-		BString varStr = *(summand->Var());
-		file->Write(buffer, sprintf(buffer, "\t\t\t\t<var>%s</var>\n",
-			varStr.String()));
-		file->Write(buffer, sprintf(buffer, "\t\t\t</summand>\n"));
-	}
-
-	file->Write(buffer, sprintf(buffer, "\t\t</leftside>\n"));
-
-	const char* op = "??";
-	if (fOp == kEQ)
-		op = "EQ";
-	else if (fOp == kLE)
-		op = "LE";
-	else if (fOp == kGE)
-		op = "GE";
-
-	file->Write(buffer, sprintf(buffer, "\t\t<op>%s</op>\n", op));
-	file->Write(buffer, sprintf(buffer, "\t\t<rightside>%f</rightside>\n", fRightSide));
-	//~ file->Write(buffer, sprintf(buffer, "\t\t<penaltyneg>%s</penaltyneg>\n", PenaltyNeg()));
-	//~ file->Write(buffer, sprintf(buffer, "\t\t<penaltypos>%s</penaltypos>\n", PenaltyPos()));
-	file->Write(buffer, sprintf(buffer, "\t</constraint>\n"));
-}
-
-
 /**
  * Gets the slack variable for the negative variations.
  *
@@ -365,6 +293,18 @@ Constraint::DPos() const
 }
 
 
+bool
+Constraint::IsSoft() const
+{
+	if (fPenaltyNeg > 0. && fOp != kLE)
+		return true;
+		
+	if (fPenaltyPos > 0. && fOp != kGE)
+		return true;
+	return false;
+}
+
+
 bool
 Constraint::IsValid()
 {
@@ -404,7 +344,8 @@ Constraint::GetString(BString& string) const
 		for (int i = 0; i < fLeftSide->CountItems(); i++) {
 			Summand* s = static_cast<Summand*>(fLeftSide->ItemAt(i));
 			string << (float)s->Coeff() << "*";
-			s->Var()->GetString(string);
+			string << "x";
+			string << s->Var()->Index() - 1;
 			string << " ";
 		}
 		string << ((fOp == kEQ) ? "== "
@@ -419,6 +360,15 @@ Constraint::GetString(BString& string) const
 }
 
 
+void
+Constraint::PrintToStream()
+{
+	BString string;
+	GetString(string);
+	printf("%s\n", string.String());
+}
+
+
 /**
  * Constructor.
  */

src/libs/linprog/Jamfile

@@ -5,12 +5,14 @@ SetSubDirSupportedPlatformsBeOSCompatible ;
 UseLibraryHeaders lp_solve linprog ;
 UsePrivateHeaders shared ;
 
+
 StaticLibrary liblinprog.a : 
 	Constraint.cpp
 	LinearSpec.cpp
 	LPSolveInterface.cpp
 	Summand.cpp
-	PenaltyFunction.cpp
 	Variable.cpp
+	LayoutOptimizer.cpp
+	ActiveSetSolver.cpp
 ;
 

src/libs/linprog/LPSolveInterface.cpp

@@ -8,14 +8,20 @@
 
 #include "LPSolveInterface.h"
 
+#include <new>
+
 
 using namespace LinearProgramming;
 
 
-LPSolveInterface::LPSolveInterface()
+LPSolveInterface::LPSolveInterface(LinearSpec* linearSpec)
 	:
+	SolverInterface(linearSpec),
+
 	fLpPresolved(NULL),
-	fLP(NULL)
+	fLP(NULL),
+	fOptimization(kMinimize),
+	fObjFunction(new(std::nothrow) SummandList())
 {
 	fLP = make_lp(0, 0);
 	if (fLP == NULL)
@@ -33,41 +39,47 @@ LPSolveInterface::~LPSolveInterface()
 {
 	_RemovePresolved();
 	delete_lp(fLP);
+
+	for (int32 i = 0; i < fObjFunction->CountItems(); i++)
+		delete (Summand*)fObjFunction->ItemAt(i);
+	delete fObjFunction;
 }
 
 
 ResultType
-LPSolveInterface::Solve(VariableList& variables)
+LPSolveInterface::Solve()
 {
+	const VariableList& variables = fLinearSpec->Variables();
 	if (fLpPresolved != NULL)
 		return _Presolve(variables);
 
-	ResultType result = (ResultType)solve(fLP);
+	// Try to solve the layout until the result is kOptimal or kInfeasible,
+	// maximally 15 tries sometimes the solving algorithm encounters numerical
+	// problems (NUMFAILURE), and repeating the solving often helps to overcome
+	// them.
+	ResultType result = kInfeasible;
+	for (int32 tries = 0; tries < 15; tries++) {
+		result = (ResultType)solve(fLP);
 
-	if (result == OPTIMAL) {
+		if (result == OPTIMAL) {
-		int32 size = variables.CountItems();
+			int32 size = variables.CountItems();
-		double x[size];
+			double x[size];
-		if (!get_variables(fLP, &x[0]))
+			if (!get_variables(fLP, &x[0]))
-			printf("Error in get_variables.\n");
+				printf("Error in get_variables.\n");
 
-		for (int32 i = 0; i < size; i++)
+			for (int32 i = 0; i < size; i++)
-			variables.ItemAt(i)->SetValue(x[i]);
+				variables.ItemAt(i)->SetValue(x[i]);
+			break;
+		} else if (result == kInfeasible)
+			break;
 	}
+
 	return result;
 }
 
 
-double
-LPSolveInterface::GetObjectiveValue()
-{
-	if (fLpPresolved)
-		return get_objective(fLpPresolved);
-	return get_objective(fLP);
-}
-
-
 bool
-LPSolveInterface::AddVariable()
+LPSolveInterface::VariableAdded(Variable* variable)
 {
 	double d = 0;
 	int i = 0;
@@ -79,9 +91,9 @@ LPSolveInterface::AddVariable()
 
 
 bool
-LPSolveInterface::RemoveVariable(int variable)
+LPSolveInterface::VariableRemoved(Variable* variable)
 {
-	if (!del_column(fLP, variable))
+	if (!del_column(fLP, variable->Index() + 1))
 		return false;
 	_RemovePresolved();
 	return true;
@@ -89,9 +101,11 @@ LPSolveInterface::RemoveVariable(int variable)
 
 
 bool
-LPSolveInterface::SetVariableRange(int variable, double min, double max)
+LPSolveInterface::VariableRangeChanged(Variable* variable)
 {
-	if (!set_bounds(fLP, variable, min, max))
+	double min = variable->Min();
+	double max = variable->Max();
+	if (!set_bounds(fLP, variable->Index() + 1, min, max))
 		return false;
 	_RemovePresolved();
 	return true;
@@ -99,22 +113,144 @@ LPSolveInterface::SetVariableRange(int variable, double min, double max)
 
 
 bool
-LPSolveInterface::AddConstraint(int nElements, double* coefficients,
+LPSolveInterface::ConstraintAdded(Constraint* constraint)
-	int* variableIndices, OperatorType op, double rightSide)
 {
-	if (!add_constraintex(fLP, nElements, coefficients, variableIndices,
+	OperatorType op = constraint->Op();
+	SummandList* summands = constraint->LeftSide();
+
+	double coeffs[summands->CountItems() + 2];
+	int variableIndices[summands->CountItems() + 2];
+	int32 nCoefficient = 0;
+	for (; nCoefficient < summands->CountItems(); nCoefficient++) {
+		Summand* s = summands->ItemAt(nCoefficient);
+		coeffs[nCoefficient] = s->Coeff();
+		variableIndices[nCoefficient] = s->Var()->Index() + 1;
+	}
+
+	double penaltyNeg = constraint->PenaltyNeg();
+	if (penaltyNeg > 0. && op != kLE) {
+		constraint->fDNegObjSummand = new(std::nothrow) Summand(
+			constraint->PenaltyNeg(), fLinearSpec->AddVariable());
+		fObjFunction->AddItem(constraint->fDNegObjSummand);
+		variableIndices[nCoefficient]
+			= constraint->fDNegObjSummand->Var()->Index() + 1;
+		coeffs[nCoefficient] = 1.0;
+		nCoefficient++;
+	}
+
+	double penaltyPos = constraint->PenaltyPos();
+	if (penaltyPos > 0. && op != kGE) {
+		constraint->fDPosObjSummand = new(std::nothrow) Summand(
+			constraint->PenaltyPos(), fLinearSpec->AddVariable());
+		fObjFunction->AddItem(constraint->fDPosObjSummand);
+		variableIndices[nCoefficient]
+			= constraint->fDPosObjSummand->Var()->Index() + 1;
+		coeffs[nCoefficient] = -1.0;
+		nCoefficient++;
+	}
+
+	double rightSide = constraint->RightSide();
+	if (!add_constraintex(fLP, nCoefficient, coeffs, variableIndices,
 		(op == kEQ ? EQ : (op == kGE) ? GE : LE), rightSide)) {
 		return false;
 	}
+
+	_UpdateObjectiveFunction();
 	_RemovePresolved();
 	return true;
 }
 
 
 bool
-LPSolveInterface::RemoveConstraint(int constraint)
+LPSolveInterface::ConstraintRemoved(Constraint* constraint)
 {
-	if (!del_constraint(fLP, constraint))
+	if (constraint->fDNegObjSummand) {
+		fObjFunction->RemoveItem(constraint->fDNegObjSummand);
+		delete constraint->fDNegObjSummand->Var();
+		delete constraint->fDNegObjSummand;
+		constraint->fDNegObjSummand = NULL;
+	}
+	if (constraint->fDPosObjSummand) {
+		fObjFunction->RemoveItem(constraint->fDPosObjSummand);
+		delete constraint->fDPosObjSummand->Var();
+		delete constraint->fDPosObjSummand;
+		constraint->fDPosObjSummand = NULL;
+	}
+
+	if (!del_constraint(fLP, constraint->Index() + 1))
+		return false;
+	_UpdateObjectiveFunction();
+	_RemovePresolved();
+	return true;
+}
+
+
+bool
+LPSolveInterface::LeftSideChanged(Constraint* constraint)
+{
+	if (!constraint->IsValid())
+		return false;
+
+	int32 index = constraint->Index() + 1;
+	if (index <= 0)
+		return false;
+
+	SummandList* leftSide = constraint->LeftSide();
+	OperatorType op = constraint->Op();
+
+	double coeffs[leftSide->CountItems() + 2];
+	int variableIndices[leftSide->CountItems() + 2];
+	int32 i;
+	for (i = 0; i < leftSide->CountItems(); i++) {
+		Summand* s = leftSide->ItemAt(i);
+		coeffs[i] = s->Coeff();
+		variableIndices[i] = s->Var()->Index() + 1;
+	}
+
+	double penaltyNeg = constraint->PenaltyNeg();
+	if (penaltyNeg > 0. && op != kLE) {
+		if (!constraint->fDNegObjSummand) {
+			constraint->fDNegObjSummand = new(std::nothrow) Summand(
+				constraint->PenaltyNeg(), fLinearSpec->AddVariable());
+			fObjFunction->AddItem(constraint->fDNegObjSummand);
+		}
+		variableIndices[i] = constraint->fDNegObjSummand->Var()->Index() + 1;
+		coeffs[i] = 1.0;
+		i++;
+	} else {
+		fObjFunction->RemoveItem(constraint->fDNegObjSummand);
+		delete constraint->fDNegObjSummand;
+		constraint->fDNegObjSummand = NULL;
+	}
+
+	double penaltyPos = constraint->PenaltyPos();
+	if (penaltyPos > 0. && op != kGE) {
+		if (constraint->fDPosObjSummand == NULL) {
+			constraint->fDPosObjSummand = new(std::nothrow) Summand(penaltyPos,
+				fLinearSpec->AddVariable());
+			fObjFunction->AddItem(constraint->fDPosObjSummand);
+		}
+		variableIndices[i] = constraint->fDPosObjSummand->Var()->Index() + 1;
+		coeffs[i] = -1.0;
+		i++;
+	} else {
+		fObjFunction->RemoveItem(constraint->fDPosObjSummand);
+		delete constraint->fDPosObjSummand;
+		constraint->fDPosObjSummand = NULL;
+	}
+
+	if (!set_rowex(fLP, index, i, coeffs, variableIndices))
+		return false;
+	_UpdateObjectiveFunction();
+	_RemovePresolved();
+	return true;
+}
+
+
+bool
+LPSolveInterface::RightSideChanged(Constraint* constraint)
+{
+	if (!set_rh(fLP, constraint->Index() + 1, constraint->RightSide()))
 		return false;
 	_RemovePresolved();
 	return true;
@@ -122,31 +258,14 @@ LPSolveInterface::RemoveConstraint(int constraint)
 
 
 bool
-LPSolveInterface::SetLeftSide(int constraint, int nElements,
+LPSolveInterface::OperatorChanged(Constraint* constraint)
-	double* coefficients, int* variableIndices)
 {
-	if (!set_rowex(fLP, constraint, nElements, coefficients, variableIndices))
+	int32 index = constraint->Index() + 1;
+	if (index <= 0)
 		return false;
-	_RemovePresolved();
-	return true;
-}
 
-
+	OperatorType op = constraint->Op();
-bool
+	if (!set_constr_type(fLP, index, op == kEQ) ? EQ : (op == kGE) ? GE : LE) {
-LPSolveInterface::SetRightSide(int constraint, double value)
-{
-	if (!set_rh(fLP, constraint, value))
-		return false;
-	_RemovePresolved();
-	return true;
-}
-
-
-bool
-LPSolveInterface::SetOperator(int constraint, OperatorType op)
-{
-	if (!set_constr_type(fLP, constraint, op == kEQ) ? EQ : (op == kGE) ? GE
-		: LE) {
 		return false;
 	}
 	_RemovePresolved();
@@ -168,7 +287,8 @@ LPSolveInterface::SetObjectiveFunction(int nElements, double* coefficients,
 bool
 LPSolveInterface::SetOptimization(OptimizationType value)
 {
-	if (value == kMinimize)
+	fOptimization = value;
+	if (fOptimization == kMinimize)
 		set_minim(fLP);
 	else
 		set_maxim(fLP);
@@ -176,6 +296,19 @@ LPSolveInterface::SetOptimization(OptimizationType value)
 }
 
 
+/**
+ * Gets the current optimization.
+ * The default is minimization.
+ *
+ * @return the current optimization
+ */
+OptimizationType
+LPSolveInterface::Optimization() const
+{
+	return fOptimization;
+}
+
+
 bool
 LPSolveInterface::SaveModel(const char* fileName)
 {
@@ -186,6 +319,120 @@ LPSolveInterface::SaveModel(const char* fileName)
 }
 
 
+BSize
+LPSolveInterface::MinSize(Variable* width, Variable* height)
+{
+	SummandList* newObjFunction = new(std::nothrow) SummandList(2);
+	newObjFunction->AddItem(new(std::nothrow) Summand(1.0, width));
+	newObjFunction->AddItem(new(std::nothrow) Summand(1.0, height));
+	SummandList* oldObjFunction = SwapObjectiveFunction(newObjFunction);
+
+	ResultType result = Solve();
+
+	SetObjectiveFunction(oldObjFunction);
+
+	if (result == kUnbounded)
+		return kMinSize;
+	if (result != kOptimal)
+		printf("Could not solve the layout specification (%d). ", result);
+
+	return BSize(width->Value(), height->Value());
+}
+
+
+BSize
+LPSolveInterface::MaxSize(Variable* width, Variable* height)
+{
+	SummandList* newObjFunction = new(std::nothrow) SummandList(2);
+	newObjFunction->AddItem(new(std::nothrow) Summand(-1.0, width));
+	newObjFunction->AddItem(new(std::nothrow) Summand(-1.0, height));
+	SummandList* oldObjFunction = SwapObjectiveFunction(
+		newObjFunction);
+
+	ResultType result = Solve();
+
+	SetObjectiveFunction(oldObjFunction);
+
+	if (result == kUnbounded)
+		return kMinSize;
+	if (result != kOptimal)
+		printf("Could not solve the layout specification (%d). ", result);
+
+	return BSize(width->Value(), height->Value());
+}
+
+
+SummandList*
+LPSolveInterface::SwapObjectiveFunction(SummandList* objFunction)
+{
+	SummandList* list = fObjFunction;
+	fObjFunction = objFunction;
+	_UpdateObjectiveFunction();
+	return list;
+}
+
+
+/**
+ * Sets a new objective function.
+ *
+ * @param summands	SummandList containing the objective function's summands
+ */
+void
+LPSolveInterface::SetObjectiveFunction(SummandList* objFunction)
+{
+	for (int32 i = 0; i < fObjFunction->CountItems(); i++)
+		delete (Summand*)fObjFunction->ItemAt(i);
+	delete fObjFunction;
+
+	fObjFunction = objFunction;
+	_UpdateObjectiveFunction();
+}
+
+
+/**
+ * Gets the objective function.
+ *
+ * @return SummandList containing the objective function's summands
+ */
+SummandList*
+LPSolveInterface::ObjectiveFunction()
+{
+	return fObjFunction;
+}
+
+
+double
+LPSolveInterface::GetObjectiveValue()
+{
+	if (fLpPresolved)
+		return get_objective(fLpPresolved);
+	return get_objective(fLP);
+}
+
+
+
+/**
+ * Updates the internal representation of the objective function.
+ * Must be called whenever the summands of the objective function are changed.
+ */
+void
+LPSolveInterface::_UpdateObjectiveFunction()
+{
+	int32 size = fObjFunction->CountItems();
+	double coeffs[size];
+	int varIndexes[size];
+	Summand* current;
+	for (int32 i = 0; i < size; i++) {
+		current = (Summand*)fObjFunction->ItemAt(i);
+		coeffs[i] = current->Coeff();
+		varIndexes[i] = current->Var()->Index() + 1;
+	}
+
+	if (!SetObjectiveFunction(size, &coeffs[0], &varIndexes[0]))
+		printf("Error in set_obj_fnex.\n");
+}
+
+
 /**
  * Remove a cached presolved model, if existent.
  * This is automatically done each time after the model has been changed,
@@ -210,7 +457,7 @@ LPSolveInterface::_RemovePresolved()
  * @return the result of the solving attempt
  */
 ResultType
-LPSolveInterface::_Presolve(VariableList& variables)
+LPSolveInterface::_Presolve(const VariableList& variables)
 {
 	if (fLpPresolved == NULL) {
 		fLpPresolved = copy_lp(fLP);
@@ -225,7 +472,7 @@ LPSolveInterface::_Presolve(VariableList& variables)
 		for (int32 i = 0; i < size; i++) {
 			Variable* current = variables.ItemAt(i);
 			current->SetValue(get_var_primalresult(fLpPresolved,
-				get_Norig_rows(fLpPresolved) + current->Index()));
+				get_Norig_rows(fLpPresolved) + current->Index() + 1));
 		}
 	}
 

src/libs/linprog/LPSolveInterface.h

@@ -11,43 +11,57 @@
 #include "lp_lib.h"
 
 
-class LPSolveInterface : public LinearProgramming::SolverInterface {
+namespace LinearProgramming {
+
+
+class LPSolveInterface : public SolverInterface {
 public:
-								LPSolveInterface();
+								LPSolveInterface(LinearSpec* linearSpec);
 								~LPSolveInterface();
 
-			ResultType			Solve(VariableList& variables);
+			ResultType			Solve();
-			double				GetObjectiveValue();
 
-			bool				AddVariable();
+			bool				VariableAdded(Variable* variable);
-			bool				RemoveVariable(int variable);
+			bool				VariableRemoved(Variable* variable);
-			bool				SetVariableRange(int variable, double min,
+			bool				VariableRangeChanged(Variable* variable);
-									double max);
 
-			bool				AddConstraint(int nElements,
+			bool				ConstraintAdded(Constraint* constraint);
-									double* coefficients, int* variableIndices,
+			bool				ConstraintRemoved(Constraint* constraint);
-									OperatorType op, double rightSide);
+			bool				LeftSideChanged(Constraint* constraint);
-			bool				RemoveConstraint(int constraint);
+			bool				RightSideChanged(Constraint* constraint);
-			bool				SetLeftSide(int constraint, int nElements,
+			bool				OperatorChanged(Constraint* constraint);
-									double* coefficients, int* variableIndices);
-			bool				SetRightSide(int constraint, double value);
-			bool				SetOperator(int constraint,
-									OperatorType op);
-
-			bool				SetObjectiveFunction(int nElements,
-									double* coefficients,
-									int* variableIndices);
-			bool				SetOptimization(OptimizationType value);
 
 			bool				SaveModel(const char* fileName);
 
+			BSize				MinSize(Variable* width, Variable* height);
+			BSize				MaxSize(Variable* width, Variable* height);
+
+			bool				SetOptimization(OptimizationType value);
+			OptimizationType	Optimization() const;
+			bool				SetObjectiveFunction(int nElements,
+									double* coefficients,
+									int* variableIndices);
+			SummandList*		ObjectiveFunction();
+			double				GetObjectiveValue();
+			//! Caller takes ownership of the Summand's and the SummandList.
+			SummandList*		SwapObjectiveFunction(
+									SummandList* objFunction);
+			void				SetObjectiveFunction(SummandList* objFunction);
 private:
-			ResultType			_Presolve(VariableList& variables);
+			void				_UpdateObjectiveFunction();
+
+			ResultType			_Presolve(const VariableList& variables);
 			void				_RemovePresolved();
 
 			lprec*				fLpPresolved;
 			lprec*				fLP;
+
+			OptimizationType	fOptimization;
+			SummandList*		fObjFunction;
 };
 
 
+}	// namespace LinearProgramming
+
+
 #endif // LP_SOLVE_INTERFACE_H

src/libs/linprog/LayoutOptimizer.cpp

@@ -0,0 +1,940 @@
+/*
+ * Copyright 2007, Ingo Weinhold <[email protected]>.
+ * Copyright 2010, Clemens Zeidler <[email protected]>
+ * Distributed under the terms of the MIT License.
+ */
+
+
+#include "LayoutOptimizer.h"
+
+#include <new>
+#include <stdio.h>
+#include <string.h>
+
+#include <AutoDeleter.h>
+
+
+//#define TRACE_LAYOUT_OPTIMIZER	1
+#if TRACE_LAYOUT_OPTIMIZER
+#	define TRACE(format...)	printf(format)
+#	define TRACE_ONLY(x)	x
+#else
+#	define TRACE(format...)
+#	define TRACE_ONLY(x)
+#endif
+#define TRACE_ERROR(format...)	fprintf(stderr, format)
+
+using std::nothrow;
+
+
+/*!	\class BPrivate::Layout::LayoutOptimizer
+
+	Given a set of layout constraints, a feasible solution, and a desired
+	(non-)solution this class finds an optimal solution. The optimization
+	criterion is to minimize the norm of the difference to the desired
+	(non-)solution.
+
+	It does so by implementing an active set method algorithm. The basic idea
+	is to start with the subset of the constraints that are barely satisfied by
+	the feasible solution, i.e. including all equality constraints and those
+	inequality constraints that are still satisfied, if restricted to equality
+	constraints. This set is called active set, the contained constraints active
+	constraints.
+
+	Considering all of the active constraints equality constraints a new
+	solution is computed, which still satisfies all those equality constraints
+	and is optimal with respect to the optimization criterion.
+
+	If the new solution equals the previous one, we find the inequality
+	constraint that, by keeping it in the active set, prevents us most from
+	further optimizing the solution. If none really does, we're done, having
+	found the globally optimal solution. Otherwise we remove the found
+	constraint from the active set and try again.
+
+	If the new solution does not equal the previous one, it might violate one
+	or more of the inactive constraints. If that is the case, we add the
+	most-violated constraint to the active set and adjust the new solution such
+	that barely satisfies that constraint. Otherwise, we don't adjust the
+	computed solution. With the adjusted respectively unadjusted solution
+	we enter the next iteration, i.e. by computing a new optimal solution with
+	respect to the active set.
+*/
+
+
+// #pragma mark - vector and matrix operations
+
+
+// is_zero
+static inline bool
+is_zero(double* x, int n)
+{
+	for (int i = 0; i < n; i++) {
+		if (!fuzzy_equals(x[i], 0))
+			return false;
+	}
+
+	return true;
+}
+
+
+// add_vectors
+static inline void
+add_vectors(double* x, const double* y, int n)
+{
+	for (int i = 0; i < n; i++)
+		x[i] += y[i];
+}
+
+
+// add_vectors_scaled
+static inline void
+add_vectors_scaled(double* x, const double* y, double scalar, int n)
+{
+	for (int i = 0; i < n; i++)
+		x[i] += y[i] * scalar;
+}
+
+
+// negate_vector
+static inline void
+negate_vector(double* x, int n)
+{
+	for (int i = 0; i < n; i++)
+		x[i] = -x[i];
+}
+
+
+// allocate_matrix
+double**
+BPrivate::Layout::allocate_matrix(int m, int n)
+{
+	double** matrix = new(nothrow) double*[m];
+	if (!matrix)
+		return NULL;
+
+	double* values = new(nothrow) double[m * n];
+	if (!values) {
+		delete[] matrix;
+		return NULL;
+	}
+
+	double* row = values;
+	for (int i = 0; i < m; i++, row += n)
+		matrix[i] = row;
+
+	return matrix;
+}
+
+
+// free_matrix
+void
+BPrivate::Layout::free_matrix(double** matrix)
+{
+	if (matrix) {
+		delete[] *matrix;
+		delete[] matrix;
+	}
+}
+
+
+// multiply_matrix_vector
+/*!	y = Ax
+	A: m x n matrix
+*/
+static inline void
+multiply_matrix_vector(const double* const* A, const double* x, int m, int n,
+	double* y)
+{
+	for (int i = 0; i < m; i++) {
+		double sum = 0;
+		for (int k = 0; k < n; k++)
+			sum += A[i][k] * x[k];
+		y[i] = sum;
+	}
+}
+
+
+// multiply_matrices
+/*!	c = a*b
+*/
+static void
+multiply_matrices(const double* const* a, const double* const* b, double** c,
+	int m, int n, int l)
+{
+	for (int i = 0; i < m; i++) {
+		for (int j = 0; j < l; j++) {
+			double sum = 0;
+			for (int k = 0; k < n; k++)
+				sum += a[i][k] * b[k][j];
+			c[i][j] = sum;
+		}	
+	}
+}
+
+
+// transpose_matrix
+static inline void
+transpose_matrix(const double* const* A, double** Atrans, int m, int n)
+{
+	for (int i = 0; i < m; i++) {
+		for (int k = 0; k < n; k++)
+			Atrans[k][i] = A[i][k];
+	}
+}
+
+
+// zero_matrix
+void
+BPrivate::Layout::zero_matrix(double** A, int m, int n)
+{
+	for (int i = 0; i < m; i++) {
+		for (int k = 0; k < n; k++)
+			A[i][k] = 0;
+	}
+}
+
+
+// copy_matrix
+void
+BPrivate::Layout::copy_matrix(const double* const* A, double** B, int m, int n)
+{
+	for (int i = 0; i < m; i++) {
+		for (int k = 0; k < n; k++)
+			B[i][k] = A[i][k];
+	}
+}
+
+
+static inline void
+multiply_optimization_matrix_vector(const double* x, int n, double* y)
+{
+	// The matrix has the form:
+	//  2 -1  0     ...   0  0
+	// -1  2 -1  0  ...   .  .
+	//  0 -1  2           .  .
+	//  .  0     .        .  .
+	//  .           .     0  0
+	//  .              . -1  0
+	//  0    ...    0 -1  2 -1
+	//  0    ...         -1  1
+	if (n == 1) {
+		y[0] = x[0];
+		return;
+	}
+
+	y[0] = 2 * x[0] - x[1];
+	for (int i = 1; i < n - 1; i++)
+		y[i] = 2 * x[i] - x[i - 1] - x[i + 1];
+	y[n - 1] = x[n - 1] - x[n - 2];
+}
+
+
+static inline void
+multiply_optimization_matrix_matrix(const double* const* A, int m, int n,
+	double** B)
+{
+	if (m == 1) {
+		memcpy(B[0], A[0], n * sizeof(double));
+		return;
+	}
+
+	for (int k = 0; k < n; k++) {
+		B[0][k] = 2 * A[0][k] - A[1][k];
+		for (int i = 1; i < m - 1; i++)
+			B[i][k] = 2 * A[i][k] - A[i - 1][k] - A[i + 1][k];
+		B[m - 1][k] = A[m - 1][k] - A[m - 2][k];
+	}
+}
+
+
+template<typename Type>
+static inline void
+swap(Type& a, Type& b)
+{
+	Type c = a;
+	a = b;
+	b = c;
+}
+
+
+// #pragma mark - algorithms
+
+
+bool
+BPrivate::Layout::solve(double** a, int n, double* b)
+{
+	// index array for row permutation
+	// Note: We could eliminate it, if we would permutate the row pointers of a.
+	int indices[n];
+	for (int i = 0; i < n; i++)
+		indices[i] = i;
+
+	// forward elimination
+	for (int i = 0; i < n - 1; i++) {
+		// find pivot
+		int pivot = i;
+		double pivotValue = fabs(a[indices[i]][i]);
+		for (int j = i + 1; j < n; j++) {
+			int index = indices[j];
+			double value = fabs(a[index][i]);
+			if (value > pivotValue) {
+				pivot = j;
+				pivotValue = value;
+			}
+		}
+
+		if (fuzzy_equals(pivotValue, 0)) {
+			TRACE_ERROR("solve(): matrix is not regular\n");
+			return false;
+		}
+
+		if (pivot != i) {
+			swap(indices[i], indices[pivot]);
+			swap(b[i], b[pivot]);
+		}
+		pivot = indices[i];
+
+		// eliminate
+		for (int j = i + 1; j < n; j++) {
+			int index = indices[j];
+			double q = -a[index][i] / a[pivot][i];
+			a[index][i] = 0;
+			for (int k = i + 1; k < n; k++)
+				a[index][k] += a[pivot][k] * q;
+			b[j] += b[i] * q;
+		}
+	}
+
+	// backwards substitution
+	for (int i = n - 1; i >= 0; i--) {
+		int index = indices[i];
+		double sum = b[i];
+		for (int j = i + 1; j < n; j++)
+			sum -= a[index][j] * b[j];
+
+		b[i] = sum / a[index][i];
+	}
+
+	return true;
+}
+
+
+int
+BPrivate::Layout::compute_dependencies(double** a, int m, int n,
+	bool* independent)
+{
+	// index array for row permutation
+	// Note: We could eliminate it, if we would permutate the row pointers of a.
+	int indices[m];
+	for (int i = 0; i < m; i++)
+		indices[i] = i;
+
+	// forward elimination
+	int iterations = (m > n ? n : m);
+	int i = 0;
+	int column = 0;
+	for (; i < iterations && column < n; i++) {
+		// find next pivot
+		int pivot = i;
+		do {
+			double pivotValue = fabs(a[indices[i]][column]);
+			for (int j = i + 1; j < m; j++) {
+				int index = indices[j];
+				double value = fabs(a[index][column]);
+				if (value > pivotValue) {
+					pivot = j;
+					pivotValue = value;
+				}
+			}
+
+			if (!fuzzy_equals(pivotValue, 0))
+				break;
+
+			column++;
+		} while (column < n);
+
+		if (column == n)
+			break;
+
+		if (pivot != i)
+			swap(indices[i], indices[pivot]);
+		pivot = indices[i];
+
+		independent[pivot] = true;
+
+		// eliminate
+		for (int j = i + 1; j < m; j++) {
+			int index = indices[j];
+			double q = -a[index][column] / a[pivot][column];
+			a[index][column] = 0;
+			for (int k = column + 1; k < n; k++)
+				a[index][k] += a[pivot][k] * q;
+		}
+
+		column++;
+	}
+
+	for (int j = i; j < m; j++)
+		independent[indices[j]] = false;
+
+	return i;
+}
+
+
+// remove_linearly_dependent_rows
+int
+BPrivate::Layout::remove_linearly_dependent_rows(double** A, double** temp,
+	bool* independentRows, int m, int n)
+{
+	// copy to temp
+	copy_matrix(A, temp, m, n);
+
+	int count = compute_dependencies(temp, m, n, independentRows);
+	if (count == m)
+		return count;
+
+	// remove the rows
+	int index = 0;
+	for (int i = 0; i < m; i++) {
+		if (independentRows[i]) {
+			if (index < i) {
+				for (int k = 0; k < n; k++)
+					A[index][k] = A[i][k];
+			}
+			index++;
+		}
+	}
+
+	return count;
+}
+
+
+/*!	QR decomposition using Householder transformations.
+*/
+bool
+qr_decomposition(double** a, int m, int n, double* d, double** q)
+{
+	if (m < n)
+		return false;
+
+	for (int j = 0; j < n; j++) {
+		// inner product of the first vector x of the (j,j) minor
+		double innerProductU = 0;
+		for (int i = j + 1; i < m; i++)
+			innerProductU = innerProductU + a[i][j] * a[i][j];
+		double innerProduct = innerProductU + a[j][j] * a[j][j];
+		if (fuzzy_equals(innerProduct, 0)) {
+			TRACE_ERROR("qr_decomposition(): 0 column %d\n", j);
+			return false;
+		}
+
+		// alpha (norm of x with opposite signedness of x_1) and thus r_{j,j}
+		double alpha = (a[j][j] < 0 ? sqrt(innerProduct) : -sqrt(innerProduct));
+		d[j] = alpha;
+
+		double beta = 1 / (alpha * a[j][j] - innerProduct);
+
+		// u = x - alpha * e_1
+		// (u is a[j..n][j])
+		a[j][j] -= alpha;
+
+		// left-multiply A_k with Q_k, thus obtaining a row of R and the A_{k+1}
+		// for the next iteration
+		for (int k = j + 1; k < n; k++) {
+			double sum = 0;
+			for (int i = j; i < m; i++)
+				sum += a[i][j] * a[i][k];
+			sum *= beta;
+
+			for (int i = j; i < m; i++)
+				a[i][k] += a[i][j] * sum;
+		}
+
+		// v = u/|u|
+		innerProductU += a[j][j] * a[j][j];
+		double beta2 = -2 / innerProductU;
+
+		// right-multiply Q with Q_k
+		// Q_k = I - 2vv^T
+		// Q * Q_k = Q - 2 * Q * vv^T
+		if (j == 0) {
+			for (int k = 0; k < m; k++) {
+				for (int i = 0; i < m; i++)
+					q[k][i] = beta2 * a[k][0] * a[i][0];
+
+				q[k][k] += 1;
+			}
+		} else {
+			for (int k = 0; k < m; k++) {
+				double sum = 0;
+				for (int i = j; i < m; i++)
+					sum += q[k][i] * a[i][j];
+				sum *= beta2;
+	
+				for (int i = j; i < m; i++)
+					q[k][i] += sum * a[i][j];
+			}
+		}
+	}
+
+	return true;
+}
+
+
+// MatrixDeleter
+struct MatrixDelete {
+	inline void operator()(double** matrix)
+	{
+		BPrivate::Layout::free_matrix(matrix);
+	}
+};
+typedef BPrivate::AutoDeleter<double*, MatrixDelete> MatrixDeleter;
+
+
+// #pragma mark - LayoutOptimizer
+
+
+// constructor
+LayoutOptimizer::LayoutOptimizer(const ConstraintList& list,
+	int32 variableCount)
+	:
+	fTemp1(NULL),
+	fTemp2(NULL),
+	fZtrans(NULL),
+	fQ(NULL),
+	fSoftConstraints(NULL),
+	fG(NULL),
+	fDesired(NULL)
+{
+	SetConstraints(list, variableCount);
+}
+
+
+// destructor
+LayoutOptimizer::~LayoutOptimizer()
+{
+	_MakeEmpty();
+}
+
+
+bool
+LayoutOptimizer::SetConstraints(const ConstraintList& list, int32 variableCount)
+{
+	fConstraints = list;
+	int32 constraintCount = fConstraints.CountItems();
+
+	if (fVariableCount != variableCount) {
+		_MakeEmpty();
+		_Init(variableCount, constraintCount);
+	}
+
+	zero_matrix(fSoftConstraints, constraintCount, fVariableCount);
+	double rightSide[constraintCount];
+	// set up soft constraint matrix
+	for (int32 c = 0; c < fConstraints.CountItems(); c++) {
+		Constraint* constraint = fConstraints.ItemAt(c);
+		if (!constraint->IsSoft()) {
+			rightSide[c] = 0;
+			continue;
+		}
+		rightSide[c] = _RightSide(constraint);
+		SummandList* summands = constraint->LeftSide();
+		for (int32 s = 0; s < summands->CountItems(); s++) {
+			Summand* summand = summands->ItemAt(s);
+			int32 variable = summand->Var()->Index();
+			if (constraint->Op() == LinearProgramming::kLE)
+				fSoftConstraints[c][variable] = -summand->Coeff();
+			else
+				fSoftConstraints[c][variable] = summand->Coeff();
+		}
+	}
+
+	// create G
+	transpose_matrix(fSoftConstraints, fTemp1, constraintCount, fVariableCount);
+	multiply_matrices(fTemp1, fSoftConstraints, fG, fVariableCount,
+		constraintCount, constraintCount);
+
+	// create d
+	multiply_matrix_vector(fTemp1, rightSide, fVariableCount, constraintCount,
+		fDesired);
+	negate_vector(fDesired, fVariableCount);
+
+	return true;
+}
+
+
+// InitCheck
+status_t
+LayoutOptimizer::InitCheck() const
+{
+	if (!fTemp1 || !fTemp2 || !fZtrans || !fQ || !fSoftConstraints || !fG
+		|| !fDesired)
+		return B_NO_MEMORY;
+	return B_OK;
+}
+
+
+double
+LayoutOptimizer::_ActualValue(Constraint* constraint, double* values) const
+{
+	SummandList* summands = constraint->LeftSide();
+	double value = 0;
+	for (int32 s = 0; s < summands->CountItems(); s++) {
+		Summand* summand = summands->ItemAt(s);
+		int32 variable = summand->Var()->Index();
+		value += values[variable] * summand->Coeff();
+	}
+	if (constraint->Op() == LinearProgramming::kLE)
+		return -value;
+	return value;
+}
+
+
+double
+LayoutOptimizer::_RightSide(Constraint* constraint)
+{
+	if (constraint->Op() == LinearProgramming::kLE)
+		return -constraint->RightSide();
+	return constraint->RightSide();
+}
+
+
+void
+LayoutOptimizer::_MakeEmpty()
+{
+	free_matrix(fTemp1);
+	free_matrix(fTemp2);
+	free_matrix(fZtrans);
+	free_matrix(fSoftConstraints);
+	free_matrix(fQ);
+	free_matrix(fG);
+
+	delete[] fDesired;
+}
+
+
+void
+LayoutOptimizer::_Init(int32 variableCount, int32 nConstraints)
+{
+	fVariableCount = variableCount;
+
+	fTemp1 = allocate_matrix(nConstraints, nConstraints);
+	fTemp2 = allocate_matrix(nConstraints, nConstraints);
+	fZtrans = allocate_matrix(nConstraints, fVariableCount);
+	fSoftConstraints = allocate_matrix(nConstraints, fVariableCount);
+	fQ = allocate_matrix(nConstraints, fVariableCount);
+	fG = allocate_matrix(nConstraints, nConstraints);
+
+	fDesired = new(std::nothrow) double[fVariableCount];
+}
+
+
+// Solve
+/*!	Solves the quadratic program (QP) given by the constraints added via
+	AddConstraint(), the additional constraint \sum_{i=0}^{n-1} x_i = size,
+	and the optimization criterion to minimize
+	\sum_{i=0}^{n-1} (x_i - desired[i])^2.
+	The \a values array must contain a feasible solution when called and will
+	be overwritten with the optimial solution the method computes.
+*/
+bool
+LayoutOptimizer::Solve(double* values)
+{
+	if (values == NULL)
+		return false;
+
+	int32 constraintCount = fConstraints.CountItems();
+
+	// allocate the active constraint matrix and its transposed matrix
+	fActiveMatrix = allocate_matrix(constraintCount, fVariableCount);
+	fActiveMatrixTemp = allocate_matrix(constraintCount, fVariableCount);
+	MatrixDeleter _(fActiveMatrix);
+	MatrixDeleter _2(fActiveMatrixTemp);
+	if (!fActiveMatrix || !fActiveMatrixTemp)
+		return false;
+
+	bool success = _Solve(values);
+	return success;
+}
+
+
+// _Solve
+bool
+LayoutOptimizer::_Solve(double* values)
+{
+	int32 constraintCount = fConstraints.CountItems();
+
+TRACE_ONLY(
+	TRACE("constraints:\n");
+	for (int32 i = 0; i < constraintCount; i++) {
+		TRACE(" %-2ld:  ", i);
+		fConstraints.ItemAt(i)->PrintToStream();
+	}
+)
+
+	// our QP is supposed to be in this form:
+	//   min_x 1/2x^TGx + x^Td
+	//   s.t. a_i^Tx = b_i,  i \in E
+	//        a_i^Tx >= b_i, i \in I
+
+	// init our initial x
+	double x[fVariableCount];
+	for (int i = 0; i < fVariableCount; i++)
+		x[i] = values[i];
+
+	// init d
+	// Note that the values of d and of G result from rewriting the
+	// ||x - desired|| we actually want to minimize.
+	double d[fVariableCount];
+	for (int i = 0; i < fVariableCount; i++)
+		d[i] = fDesired[i];
+
+	// init active set
+	ConstraintList activeConstraints(constraintCount);
+
+	for (int32 i = 0; i < constraintCount; i++) {
+		Constraint* constraint = (Constraint*)fConstraints.ItemAt(i);
+		if (constraint->IsSoft())
+			continue;
+		double actualValue = _ActualValue(constraint, x);
+		TRACE("constraint %ld: actual: %f  constraint: %f\n", i, actualValue,
+			_RightSide(constraint));
+		if (fuzzy_equals(actualValue, _RightSide(constraint)))
+			activeConstraints.AddItem(constraint);
+	}
+
+	// The main loop: Each iteration we try to get closer to the optimum
+	// solution. We compute a vector p that brings our x closer to the optimum.
+	// We do that by computing the QP resulting from our active constraint set,
+	// W^k. Afterward each iteration we adjust the active set.
+TRACE_ONLY(int iteration = 0;)
+	while (true) {
+TRACE_ONLY(
+		TRACE("\n[iteration %d]\n", iteration++);
+		TRACE("active set:\n");
+		for (int32 i = 0; i < activeConstraints.CountItems(); i++) {
+			TRACE("  ");
+			activeConstraints.ItemAt(i)->PrintToStream();
+		}
+)
+
+		// solve the QP:
+		//   min_p 1/2p^TGp + g_k^Tp
+		//   s.t. a_i^Tp = 0
+		//   with a_i \in activeConstraints
+		//        g_k = Gx_k + d
+		//        p = x - x_k
+
+		int32 activeCount = activeConstraints.CountItems();
+		if (activeCount == 0) {
+			TRACE_ERROR("Solve(): Error: No more active constraints!\n");
+			return false;
+		}
+
+		// construct a matrix from the active constraints
+		int am = activeCount;
+		const int an = fVariableCount;
+		bool independentRows[activeCount];
+		zero_matrix(fActiveMatrix, am, an);
+
+		for (int32 i = 0; i < activeCount; i++) {
+			Constraint* constraint = activeConstraints.ItemAt(i);
+			SummandList* summands = constraint->LeftSide();
+			for (int32 s = 0; s < summands->CountItems(); s++) {
+				Summand* summand = summands->ItemAt(s);
+				int32 variable = summand->Var()->Index();
+				if (constraint->Op() == LinearProgramming::kLE)
+					fActiveMatrix[i][variable] = -summand->Coeff();
+				else
+					fActiveMatrix[i][variable] = summand->Coeff();
+			}
+		}
+
+// TODO: The fActiveMatrix is sparse (max 2 entries per row). There should be
+// some room for optimization.
+		am = remove_linearly_dependent_rows(fActiveMatrix, fActiveMatrixTemp,
+			independentRows, am, an);
+
+		// gxd = G * x + d
+		double gxd[fVariableCount];
+		multiply_matrix_vector(fG, x, fVariableCount, fVariableCount, gxd);
+		add_vectors(gxd, d, fVariableCount);
+
+		double p[fVariableCount];
+		if (!_SolveSubProblem(gxd, am, p))
+			return false;
+
+		if (is_zero(p, fVariableCount)) {
+			// compute Lagrange multipliers lambda_i
+			// if lambda_i >= 0 for all i \in W^k \union inequality constraints,
+			// then we're done.
+			// Otherwise remove the constraint with the smallest lambda_i
+			// from the active set.
+			// The derivation of the Lagrangian yields:
+			//   \sum_{i \in W^k}(lambda_ia_i) = Gx_k + d
+			// Which is an system we can solve:
+			//   A^Tlambda = Gx_k + d
+
+			// A^T is over-determined, hence we need to reduce the number of
+			// rows before we can solve it.
+			bool independentColumns[an];
+			double** aa = fTemp1;
+			transpose_matrix(fActiveMatrix, aa, am, an);
+			const int aam = remove_linearly_dependent_rows(aa, fTemp2,
+				independentColumns, an, am);
+			const int aan = am;
+			if (aam != aan) {
+				// This should not happen, since A has full row rank.
+				TRACE_ERROR("Solve(): Transposed A has less linear independent "
+					"rows than it has columns!\n");
+				return false;
+			}
+
+			// also reduce the number of rows on the right hand side
+			double lambda[aam];
+			int index = 0;
+			for (int i = 0; i < an; i++) {
+				if (independentColumns[i])
+					lambda[index++] = gxd[i];
+			}
+
+			bool success = solve(aa, aam, lambda);
+			if (!success) {
+				// Impossible, since we've removed all linearly dependent rows.
+				TRACE_ERROR("Solve(): Failed to compute lambda!\n");
+				return false;
+			}
+
+			// find min lambda_i (only, if it's < 0, though)
+			double minLambda = 0;
+			int minIndex = -1;
+			index = 0;
+			for (int i = 0; i < activeCount; i++) {
+				if (independentRows[i]) {
+					Constraint* constraint
+						= (Constraint*)activeConstraints.ItemAt(i);
+					if (constraint->Op() != LinearProgramming::kEQ) {
+						if (lambda[index] < minLambda) {
+							minLambda = lambda[index];
+							minIndex = i;
+						}
+					}
+
+					index++;
+				}
+			}
+
+			// if the min lambda is >= 0, we're done
+			if (minIndex < 0 || fuzzy_equals(minLambda, 0)) {
+				_SetResult(x, values);
+				return true;
+			}
+
+			// remove i from the active set
+			activeConstraints.RemoveItemAt(minIndex);
+		} else {
+			// compute alpha_k
+			double alpha = 1;
+			int barrier = -1;
+			// if alpha_k < 1, add a barrier constraint to W^k
+			for (int32 i = 0; i < constraintCount; i++) {
+				Constraint* constraint = (Constraint*)fConstraints.ItemAt(i);
+				if (activeConstraints.HasItem(constraint))
+					continue;
+
+				double divider = _ActualValue(constraint, p);
+				if (divider > 0 || fuzzy_equals(divider, 0))
+					continue;
+
+				// (b_i - a_i^Tx_k) / a_i^Tp_k
+				double alphaI = _RightSide(constraint)
+					- _ActualValue(constraint, x);
+				alphaI /= divider;
+				if (alphaI < alpha) {
+					alpha = alphaI;
+					barrier = i;
+				}
+			}
+			TRACE("alpha: %f, barrier: %d\n", alpha, barrier);
+
+			if (alpha < 1)
+				activeConstraints.AddItem(fConstraints.ItemAt(barrier));
+
+			// x += p * alpha;
+			add_vectors_scaled(x, p, alpha, fVariableCount);
+		}
+	}
+}
+
+
+bool
+LayoutOptimizer::_SolveSubProblem(const double* d, int am, double* p)
+{
+	// We have to solve the QP subproblem:
+	//   min_p 1/2p^TGp + d^Tp
+	//   s.t. a_i^Tp = 0
+	//   with a_i \in activeConstraints
+	//
+	// We use the null space method, i.e. we find matrices Y and Z, such that
+	// AZ = 0 and [Y Z] is regular. Then with
+	//   p = Yp_Y + Zp_z
+	// we get
+	//   p_Y = 0
+	// and
+	//  (Z^TGZ)p_Z = -(Z^TYp_Y + Z^Tg) = -Z^Td
+	// which is a linear equation system, which we can solve.
+
+	const int an = fVariableCount;
+
+	// we get Y and Z by QR decomposition of A^T
+	double tempD[am];
+	double** const Q = fQ;
+	transpose_matrix(fActiveMatrix, fTemp1, am, an);
+	bool success = qr_decomposition(fTemp1, an, am, tempD, Q);
+	if (!success) {
+		TRACE_ERROR("Solve(): QR decomposition failed!\n");
+		return false;
+	}
+
+	// Z is the (1, m + 1) minor of Q
+	const int zm = an;
+	const int zn = an - am;
+
+	double* Z[zm];
+	for (int i = 0; i < zm; i++)
+		Z[i] = Q[i] + am;
+
+	// solve (Z^TGZ)p_Z = -Z^Td
+
+	// Z^T
+	transpose_matrix(Z, fZtrans, zm, zn);
+	// rhs: -Z^T * d;
+	double pz[zm];
+	multiply_matrix_vector(fZtrans, d, zn, zm, pz);
+	negate_vector(pz, zn);
+
+	// fTemp2 = Ztrans * G * Z
+	//multiply_optimization_matrix_matrix(Z, an, zn, fTemp1);
+	multiply_matrices(fG, Z, fTemp1, zm, fVariableCount, zn);
+	multiply_matrices(fZtrans, fTemp1, fTemp2, zn, zm, zn);
+
+	success = solve(fTemp2, zn, pz);
+	if (!success) {
+		TRACE_ERROR("Solve(): Failed to solve() system for p_Z\n");
+		return false;
+	}
+
+	// p = Z * pz;
+	multiply_matrix_vector(Z, pz, zm, zn, p);
+
+	return true;
+}
+
+
+// _SetResult
+void
+LayoutOptimizer::_SetResult(const double* x, double* values)
+{
+	for (int i = 1; i < fVariableCount; i++)
+		values[i] = x[i];
+}

src/libs/linprog/LayoutOptimizer.h

@@ -0,0 +1,85 @@
+/*
+ * Copyright 2007, Ingo Weinhold <[email protected]>.
+ * All rights reserved. Distributed under the terms of the MIT License.
+ */
+#ifndef LAYOUT_OPTIMIZER_H
+#define LAYOUT_OPTIMIZER_H
+
+#include <List.h>
+#include <math.h>
+
+#include "LinearSpec.h"
+
+
+static const double kEqualsEpsilon = 0.000001;
+
+
+namespace BPrivate {
+namespace Layout {
+
+
+double** allocate_matrix(int m, int n);
+void free_matrix(double** matrix);
+void copy_matrix(const double* const* A, double** B, int m, int n);
+void zero_matrix(double** A, int m, int n);
+int compute_dependencies(double** a, int m, int n, bool* independent);
+int remove_linearly_dependent_rows(double** A, double** temp,
+	bool* independentRows, int m, int n);
+bool solve(double** a, int n, double* b);
+
+
+class LayoutOptimizer {
+public:
+								LayoutOptimizer(const ConstraintList& list,
+									int32 variableCount);
+								~LayoutOptimizer();
+
+			bool				SetConstraints(const ConstraintList& list,
+									int32 variableCount);
+
+			status_t			InitCheck() const;
+
+			bool				Solve(double* initialSolution);
+
+private:
+			double				_ActualValue(Constraint* constraint,
+									double* values) const;
+			double				_RightSide(Constraint* constraint);
+
+			void				_MakeEmpty();
+			void				_Init(int32 variableCount, int32 nConstraints);
+
+			bool				_Solve(double* values);
+			bool				_SolveSubProblem(const double* d, int am,
+									double* p);
+			void				_SetResult(const double* x, double* values);
+
+
+			int32				fVariableCount;
+			ConstraintList		fConstraints;
+
+			double**			fTemp1;
+			double**			fTemp2;
+			double**			fZtrans;
+			double**			fQ;
+			double**			fActiveMatrix;
+			double**			fActiveMatrixTemp;
+			double**			fSoftConstraints;
+			double**			fG;
+			double*				fDesired;
+};
+
+}	// namespace Layout
+}	// namespace BPrivate
+
+using BPrivate::Layout::LayoutOptimizer;
+
+
+inline bool
+fuzzy_equals(double a, double b)
+{
+	return fabs(a - b) < kEqualsEpsilon;
+}
+
+
+#endif	// LAYOUT_OPTIMIZER_H

src/libs/linprog/LinearSpec.cpp

@@ -12,6 +12,28 @@
 #include <stdio.h>
 
 #include "LPSolveInterface.h"
+#include "ActiveSetSolver.h"
+
+
+using namespace LinearProgramming;
+
+
+#define DEBUG_LINEAR_SPECIFICATIONS
+
+#ifdef DEBUG_LINEAR_SPECIFICATIONS
+#include <stdio.h>
+#define TRACE(x...) printf(x)
+#else
+#define TRACE(x...) /* nothing */
+#endif
+
+
+SolverInterface::SolverInterface(LinearSpec* linSpec)
+	:
+	fLinearSpec(linSpec)
+{
+
+}
 
 
 /**
@@ -20,13 +42,11 @@
  */
 LinearSpec::LinearSpec()
 	:
-	fOptimization(kMinimize),
-	fObjFunction(new(std::nothrow) SummandList()),
 	fResult(kError),
-	fObjectiveValue(NAN),
+	fSolvingTime(0)
-	fSolvingTime(NAN)
 {
-	fSolver = new LPSolveInterface;
+	//fSolver = new LPSolveInterface(this);
+	fSolver = new ActiveSetSolver(this);
 }
 
 
@@ -39,12 +59,9 @@ LinearSpec::~LinearSpec()
 {
 	for (int32 i = 0; i < fConstraints.CountItems(); i++)
 		delete (Constraint*)fConstraints.ItemAt(i);
-	for (int32 i = 0; i < fObjFunction->CountItems(); i++)
-		delete (Summand*)fObjFunction->ItemAt(i);
 	while (fVariables.CountItems() > 0)
 		RemoveVariable(fVariables.ItemAt(0));
 
-	delete fObjFunction;
 	delete fSolver;
 }
 
@@ -77,7 +94,7 @@ LinearSpec::AddVariable(Variable* variable)
 
 	if (!fVariables.AddItem(variable))
 		return false;
-	if (!fSolver->AddVariable()) {
+	if (!fSolver->VariableAdded(variable)) {
 		fVariables.RemoveItem(variable);
 		return false;
 	}
@@ -95,13 +112,10 @@ LinearSpec::AddVariable(Variable* variable)
 bool
 LinearSpec::RemoveVariable(Variable* variable, bool deleteVariable)
 {
-	int32 index = IndexOf(variable);
+	// must be called first otherwise the index is invalid
-	if (index < 0)
+	if (!fSolver->VariableRemoved(variable))
 		return false;
-
+	fVariables.RemoveItem(variable);
-	if (!fSolver->RemoveVariable(index))
-		return false;
-	fVariables.RemoveItemAt(index - 1);
 	variable->fIsValid = false;
 
 	// invalidate all constraints that use this variable
@@ -135,20 +149,14 @@ LinearSpec::RemoveVariable(Variable* variable, bool deleteVariable)
 int32
 LinearSpec::IndexOf(const Variable* variable) const
 {
-	int32 i = fVariables.IndexOf(variable);
+	return fVariables.IndexOf(variable);
-	if (i == -1) {
-		printf("Variable 0x%p not part of fLS->Variables().\n", variable);
-		return -1;
-	}
-	return i + 1;
 }
 
 
 bool
 LinearSpec::UpdateRange(Variable* variable)
 {
-	if (!fSolver->SetVariableRange(IndexOf(variable), variable->Min(),
+	if (!fSolver->VariableRangeChanged(variable))
-		variable->Max()))
 		return false;
 	return true;
 }
@@ -157,49 +165,13 @@ LinearSpec::UpdateRange(Variable* variable)
 bool
 LinearSpec::AddConstraint(Constraint* constraint)
 {
-	SummandList* summands = constraint->LeftSide();
+	if (!fConstraints.AddItem(constraint))
-	OperatorType op = constraint->Op();
+		return false;
-	double rightSide = constraint->RightSide();
-	double penaltyNeg = constraint->PenaltyNeg();
-	double penaltyPos = constraint->PenaltyPos();
 
-	double coeffs[summands->CountItems() + 2];
+	if (!fSolver->ConstraintAdded(constraint)) {
-	int varIndexes[summands->CountItems() + 2];
+		fConstraints.RemoveItem(constraint);
-	int32 nCoefficient = 0;
-	for (; nCoefficient < summands->CountItems(); nCoefficient++) {
-		Summand* s = summands->ItemAt(nCoefficient);
-		coeffs[nCoefficient] = s->Coeff();
-		varIndexes[nCoefficient] = s->Var()->Index();
-	}
-
-	if (penaltyNeg != INFINITY && penaltyNeg != 0. && op != LE) {
-		constraint->fDNegObjSummand
-			= new(std::nothrow) Summand(constraint->PenaltyNeg(),
-				AddVariable());
-		fObjFunction->AddItem(constraint->fDNegObjSummand);
-		varIndexes[nCoefficient] = constraint->fDNegObjSummand->Var()->Index();
-		coeffs[nCoefficient] = 1.0;
-		nCoefficient++;
-	}
-
-	if (penaltyPos != INFINITY && penaltyPos != 0. && op != GE) {
-		constraint->fDPosObjSummand
-			= new(std::nothrow) Summand(constraint->PenaltyPos(),
-				AddVariable());
-		fObjFunction->AddItem(constraint->fDPosObjSummand);
-		varIndexes[nCoefficient] = constraint->fDPosObjSummand->Var()->Index();
-		coeffs[nCoefficient] = -1.0;
-		nCoefficient++;
-	}
-
-	if (!fSolver->AddConstraint(nCoefficient, &coeffs[0], &varIndexes[0], op,
-		rightSide)) {
 		return false;
 	}
-
-	UpdateObjectiveFunction();
-	fConstraints.AddItem(constraint);
-
 	return true;
 }
 
@@ -207,20 +179,7 @@ LinearSpec::AddConstraint(Constraint* constraint)
 bool
 LinearSpec::RemoveConstraint(Constraint* constraint, bool deleteConstraint)
 {
-	if (constraint->fDNegObjSummand) {
+	fSolver->ConstraintRemoved(constraint);
-		fObjFunction->RemoveItem(constraint->fDNegObjSummand);
-		delete constraint->fDNegObjSummand->Var();
-		delete constraint->fDNegObjSummand;
-		constraint->fDNegObjSummand = NULL;
-	}
-	if (constraint->fDPosObjSummand) {
-		fObjFunction->RemoveItem(constraint->fDPosObjSummand);
-		delete constraint->fDPosObjSummand->Var();
-		delete constraint->fDPosObjSummand;
-		constraint->fDPosObjSummand = NULL;
-	}
-
-	fSolver->RemoveConstraint(constraint->Index());
 	fConstraints.RemoveItem(constraint);
 	constraint->fIsValid = false;
 
@@ -233,38 +192,8 @@ LinearSpec::RemoveConstraint(Constraint* constraint, bool deleteConstraint)
 bool
 LinearSpec::UpdateLeftSide(Constraint* constraint)
 {
-	if (!constraint->IsValid())
+	if (!fSolver->LeftSideChanged(constraint))
 		return false;
-
-	SummandList* leftSide = constraint->LeftSide();
-	OperatorType op = constraint->Op();
-
-	double coeffs[leftSide->CountItems() + 2];
-	int varIndexes[leftSide->CountItems() + 2];
-	int32 i;
-	for (i = 0; i < leftSide->CountItems(); i++) {
-		Summand* s = leftSide->ItemAt(i);
-		coeffs[i] = s->Coeff();
-		varIndexes[i] = s->Var()->Index();
-	}
-
-	if (constraint->fDNegObjSummand != NULL && op != OperatorType(LE)) {
-		varIndexes[i] = constraint->fDNegObjSummand->Var()->Index();
-		coeffs[i] = 1.0;
-		i++;
-	}
-
-	if (constraint->fDPosObjSummand != NULL && op != OperatorType(GE)) {
-		varIndexes[i] = constraint->fDPosObjSummand->Var()->Index();
-		coeffs[i] = -1.0;
-		i++;
-	}
-
-	if (!fSolver->SetLeftSide(constraint->Index(), i, &coeffs[0],
-		&varIndexes[0]))
-		return false;
-
-	UpdateObjectiveFunction();
 	return true;
 }
 
@@ -272,7 +201,7 @@ LinearSpec::UpdateLeftSide(Constraint* constraint)
 bool
 LinearSpec::UpdateRightSide(Constraint* constraint)
 {
-	if (!fSolver->SetRightSide(constraint->Index(), constraint->RightSide()))
+	if (!fSolver->RightSideChanged(constraint))
 		return false;
 	return true;
 }
@@ -281,7 +210,7 @@ LinearSpec::UpdateRightSide(Constraint* constraint)
 bool
 LinearSpec::UpdateOperator(Constraint* constraint)
 {
-	if (!fSolver->SetOperator(constraint->Index(), constraint->Op()))
+	if (!fSolver->OperatorChanged(constraint))
 		return false;
 	return true;
 }
@@ -300,7 +229,7 @@ Constraint*
 LinearSpec::AddConstraint(SummandList* summands, OperatorType op,
 	double rightSide)
 {
-	return AddConstraint(summands, op, rightSide, INFINITY, INFINITY);
+	return AddConstraint(summands, op, rightSide, -1, -1);
 }
 
 
@@ -317,7 +246,7 @@ Constraint*
 LinearSpec::AddConstraint(double coeff1, Variable* var1,
 	OperatorType op, double rightSide)
 {
-	return AddConstraint(coeff1, var1, op, rightSide, INFINITY, INFINITY);
+	return AddConstraint(coeff1, var1, op, rightSide, -1, -1);
 }
 
 
@@ -336,8 +265,8 @@ Constraint*
 LinearSpec::AddConstraint(double coeff1, Variable* var1,
 	double coeff2, Variable* var2, OperatorType op, double rightSide)
 {
-	return AddConstraint(coeff1, var1, coeff2, var2, op, rightSide, INFINITY,
+	return AddConstraint(coeff1, var1, coeff2, var2, op, rightSide, -1,
-		INFINITY);
+		-1);
 }
 
 
@@ -360,7 +289,7 @@ LinearSpec::AddConstraint(double coeff1, Variable* var1,
 		OperatorType op, double rightSide)
 {
 	return AddConstraint(coeff1, var1, coeff2, var2, coeff3, var3, op,
-		rightSide, INFINITY, INFINITY);
+		rightSide, -1, -1);
 }
 
 
@@ -385,7 +314,7 @@ LinearSpec::AddConstraint(double coeff1, Variable* var1,
 		double coeff4, Variable* var4, OperatorType op, double rightSide)
 {
 	return AddConstraint(coeff1, var1, coeff2, var2, coeff3, var3, coeff4, var4,
-		op, rightSide, INFINITY, INFINITY);
+		op, rightSide, -1, -1);
 }
 
 
@@ -526,79 +455,17 @@ LinearSpec::AddConstraint(double coeff1, Variable* var1,
 }
 
 
-/**
+BSize
- * Adds a new penalty function to the specification.
+LinearSpec::MinSize(Variable* width, Variable* height)
- *
- * @param var	the penalty function's variable
- * @param xs	the penalty function's sampling points
- * @param gs	the penalty function's gradients
- * @return the new penalty function
- */
-PenaltyFunction*
-LinearSpec::AddPenaltyFunction(Variable* var, BList* xs, BList* gs)
 {
-	return new(std::nothrow) PenaltyFunction(this, var, xs, gs);
+	return fSolver->MinSize(width, height);
 }
 
 
-/**
+BSize
- * Gets the objective function.
+LinearSpec::MaxSize(Variable* width, Variable* height)
- *
- * @return SummandList containing the objective function's summands
- */
-SummandList*
-LinearSpec::ObjectiveFunction()
 {
-	return fObjFunction;
+	return fSolver->MaxSize(width, height);
-}
-
-
-SummandList*
-LinearSpec::SwapObjectiveFunction(SummandList* objFunction)
-{
-	SummandList* list = fObjFunction;
-	fObjFunction = objFunction;
-	UpdateObjectiveFunction();
-	return list;
-}
-
-
-/**
- * Sets a new objective function.
- *
- * @param summands	SummandList containing the objective function's summands
- */
-void
-LinearSpec::SetObjectiveFunction(SummandList* objFunction)
-{
-	for (int32 i = 0; i < fObjFunction->CountItems(); i++)
-		delete (Summand*)fObjFunction->ItemAt(i);
-	delete fObjFunction;
-
-	fObjFunction = objFunction;
-	UpdateObjectiveFunction();
-}
-
-
-/**
- * Updates the internal representation of the objective function.
- * Must be called whenever the summands of the objective function are changed.
- */
-void
-LinearSpec::UpdateObjectiveFunction()
-{
-	int32 size = fObjFunction->CountItems();
-	double coeffs[size];
-	int varIndexes[size];
-	Summand* current;
-	for (int32 i = 0; i < size; i++) {
-		current = (Summand*)fObjFunction->ItemAt(i);
-		coeffs[i] = current->Coeff();
-		varIndexes[i] = current->Var()->Index();
-	}
-
-	if (!fSolver->SetObjectiveFunction(size, &coeffs[0], &varIndexes[0]))
-		printf("Error in set_obj_fnex.\n");
 }
 
 
@@ -638,23 +505,27 @@ LinearSpec::_AddConstraint(SummandList* leftSide, OperatorType op,
 }
 
 
-/**
+#ifdef DEBUG_LINEAR_SPECIFICATIONS
- * Tries to solve the linear programming problem.
+static bigtime_t sAverageSolvingTime = 0;
- * If a cached simplified version of the problem exists, it is used instead.
+static int32 sSolvedCount = 0;
- *
+#endif
- * @return the result of the solving attempt
+
- */
+
 ResultType
 LinearSpec::Solve()
 {
-	bigtime_t start, end;
+	bigtime_t startTime = system_time();
-	start = system_time();
 
-	fResult = fSolver->Solve(fVariables);
+	fResult = fSolver->Solve();
-	fObjectiveValue = fSolver->GetObjectiveValue();
 
-	end = system_time();
+	fSolvingTime = system_time() - startTime;
-	fSolvingTime = (end - start) / 1000.0;
+
+#ifdef DEBUG_LINEAR_SPECIFICATIONS
+	sAverageSolvingTime += fSolvingTime;
+	sSolvedCount++;
+	TRACE("Solving time %i average %i [micro s]\n", (int)fSolvingTime,
+		int(sAverageSolvingTime / sSolvedCount));
+#endif
 
 	return fResult;
 }
@@ -673,33 +544,6 @@ LinearSpec::Save(const char* fileName)
 }
 
 
-/**
- * Gets the current optimization.
- * The default is minimization.
- *
- * @return the current optimization
- */
-OptimizationType
-LinearSpec::Optimization() const
-{
-	return fOptimization;
-}
-
-
-/**
- * Sets whether the solver should minimize or maximize the objective function.
- * The default is minimization.
- *
- * @param optimization	the optimization type
- */
-void
-LinearSpec::SetOptimization(OptimizationType value)
-{
-	fOptimization = value;
-	fSolver->SetOptimization(value);
-}
-
-
 /**
  * Gets the constraints.
  *
@@ -712,6 +556,13 @@ LinearSpec::Constraints() const
 }
 
 
+const VariableList&
+LinearSpec::Variables() const
+{
+	return fVariables;
+}
+
+
 /**
  * Gets the result type.
  *
@@ -724,24 +575,12 @@ LinearSpec::Result() const
 }
 
 
-/**
- * Gets the objective value.
- *
- * @return the objective value
- */
-double
-LinearSpec::ObjectiveValue() const
-{
-	return fObjectiveValue;
-}
-
-
 /**
  * Gets the solving time.
  *
  * @return the solving time
  */
-double
+bigtime_t
 LinearSpec::SolvingTime() const
 {
 	return fSolvingTime;
@@ -789,6 +628,6 @@ LinearSpec::GetString(BString& string) const
 		string << "kNumFailure";
 	else
 		string << fResult;
-	string << " SolvingTime=" << (float)fSolvingTime << "ms";
+	string << " SolvingTime=" << fSolvingTime << "micro s";
 }
 

src/libs/linprog/PenaltyFunction.cpp

@@ -1,70 +0,0 @@
-/*
- * Copyright 2007-2008, Christof Lutteroth, [email protected]
- * Copyright 2007-2008, James Kim, [email protected]
- * Distributed under the terms of the MIT License.
- */
-
-
-#include "PenaltyFunction.h"
-
-#include <stdio.h>
-
-#include "Constraint.h"
-#include "LinearSpec.h"
-#include "Summand.h"
-#include "Variable.h"
-
-
-/**
- * Constructor.
- */
-PenaltyFunction::PenaltyFunction(LinearSpec* ls, Variable* var, BList* xs, BList* gs)
-{
-	int32 sizeXs = xs->CountItems();
-	int32 sizeGs = gs->CountItems();
-	if (var->LS() != ls)
-		printf("The variable must belong to the same linear specification as the penalty function.");
-	if (sizeXs + 1 != sizeGs)
-		printf("The number of sampling points must be exactly one less than the number of gradients.");
-	for (int32 i = 1; i < sizeGs; i++) {
-		if (*(double*)(gs->ItemAt(i - 1)) > *(double*)(gs->ItemAt(i)))
-			printf("Penalty function must be concave.");
-	}
-	
-	fVar = var;
-	fXs = xs;
-	fGs = gs;
-	fConstraints = new BList(sizeGs + 1);
-	fObjFunctionSummands = new BList(sizeGs);
-		
-	fConstraints->AddItem(ls->AddConstraint(1.0, var, kEQ, 
-		*(double*)(xs->ItemAt(0)), -*(double*)(gs->ItemAt(0)), 
-		*(double*)(gs->ItemAt(1))));
-	
-	for (int32 i = 1; i < sizeGs; i++) {
-		Variable* dPos = ls->AddVariable();
-			
-		fConstraints->AddItem(ls->AddConstraint(1.0, var, -1.0, dPos, kLE, 
-			*(double*)(xs->ItemAt(i))));
-		
-		Summand* objSummand =  new Summand(*(double*)(gs->ItemAt(i + 1)) - *(double*)(gs->ItemAt(i)), dPos);
-		ls->ObjectiveFunction()->AddItem(objSummand);
-		fObjFunctionSummands->AddItem(objSummand);
-	}
-	ls->UpdateObjectiveFunction();
-}
-
-
-/**
- * Destructor.
- * Removes all constraints and summands from the penalty function.
- */
-PenaltyFunction::~PenaltyFunction()
-{
-	for (int32 i = 0; i < fConstraints->CountItems(); i++)
-		delete (Constraint*)fConstraints->ItemAt(i);
-	
-	for (int32 i = 0; i < fObjFunctionSummands->CountItems(); i++)
-		delete (Summand*)fObjFunctionSummands->ItemAt(i);
-}
-

src/libs/linprog/Summand.cpp

@@ -57,6 +57,13 @@ Summand::SetVar(Variable* var)
 }
 
 
+int32
+Summand::VariableIndex()
+{
+	return fVar->Index();
+}
+
+
 /**
  * Destructor.
  */