From 7583db5a1ec81636829331e5fc7cae81ca7a7ba7 Mon Sep 17 00:00:00 2001 From: Clemens Zeidler Date: Tue, 25 Jan 2011 04:59:40 +0000 Subject: [PATCH] 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 | 4 +- headers/libs/alm/Area.h | 12 +- headers/libs/linprog/Constraint.h | 12 +- headers/libs/linprog/LinearSpec.h | 71 +- headers/libs/linprog/PenaltyFunction.h | 50 -- headers/libs/linprog/Summand.h | 1 + src/libs/alm/ALMLayout.cpp | 89 +-- src/libs/alm/Area.cpp | 42 +- src/libs/linprog/ActiveSetSolver.cpp | 542 ++++++++++++++ src/libs/linprog/ActiveSetSolver.h | 82 +++ src/libs/linprog/Constraint.cpp | 106 +-- src/libs/linprog/Jamfile | 4 +- src/libs/linprog/LPSolveInterface.cpp | 357 ++++++++-- src/libs/linprog/LPSolveInterface.h | 60 +- src/libs/linprog/LayoutOptimizer.cpp | 940 +++++++++++++++++++++++++ src/libs/linprog/LayoutOptimizer.h | 85 +++ src/libs/linprog/LinearSpec.cpp | 313 ++------ src/libs/linprog/PenaltyFunction.cpp | 70 -- src/libs/linprog/Summand.cpp | 7 + 19 files changed, 2207 insertions(+), 640 deletions(-) delete mode 100644 headers/libs/linprog/PenaltyFunction.h create mode 100644 src/libs/linprog/ActiveSetSolver.cpp create mode 100644 src/libs/linprog/ActiveSetSolver.h create mode 100644 src/libs/linprog/LayoutOptimizer.cpp create mode 100644 src/libs/linprog/LayoutOptimizer.h delete mode 100644 src/libs/linprog/PenaltyFunction.cpp diff --git a/headers/libs/alm/ALMLayout.h b/headers/libs/alm/ALMLayout.h index c477400e61..67371494ae 100644 --- a/headers/libs/alm/ALMLayout.h +++ b/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 diff --git a/headers/libs/alm/Area.h b/headers/libs/alm/Area.h index 7f2cebb572..226269fab9 100644 --- a/headers/libs/alm/Area.h +++ b/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; diff --git a/headers/libs/linprog/Constraint.h b/headers/libs/linprog/Constraint.h index d69c8a8bea..0ab9be4289 100644 --- a/headers/libs/linprog/Constraint.h +++ b/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 penaltyPos); + double rightSide, + double penaltyNeg = -1, + double penaltyPos = -1); private: LinearSpec* fLS; @@ -92,6 +93,7 @@ private: public: friend class LinearSpec; + friend class LPSolveInterface; }; diff --git a/headers/libs/linprog/LinearSpec.h b/headers/libs/linprog/LinearSpec.h index 7970f6c543..7cc2db1d94 100644 --- a/headers/libs/linprog/LinearSpec.h +++ b/headers/libs/linprog/LinearSpec.h @@ -11,47 +11,51 @@ #include #include +#include #include #include #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 double GetObjectiveValue() = 0; + virtual ResultType Solve() = 0; - virtual bool AddVariable() = 0; - virtual bool RemoveVariable(int variable) = 0; - virtual bool SetVariableRange(int variable, double min, - double max) = 0; + virtual bool VariableAdded(Variable* variable) = 0; + virtual bool VariableRemoved(Variable* variable) = 0; + virtual bool VariableRangeChanged(Variable* variable) = 0; - virtual bool AddConstraint(int nElements, - double* coefficients, int* variableIndices, - OperatorType op, double rightSide) = 0; - virtual bool RemoveConstraint(int constraint) = 0; - virtual bool SetLeftSide(int constraint, int nElements, - 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 ConstraintAdded(Constraint* constraint) = 0; + virtual bool ConstraintRemoved(Constraint* constraint) = 0; + virtual bool LeftSideChanged(Constraint* constraint) = 0; + virtual bool RightSideChanged(Constraint* constraint) = 0; + virtual bool OperatorChanged(Constraint* constraint) = 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, - BList* gs); - - SummandList* ObjectiveFunction(); - //! Caller takes ownership of the Summand's and the SummandList. - SummandList* SwapObjectiveFunction( - SummandList* objFunction); - void SetObjectiveFunction(SummandList* objFunction); - void UpdateObjectiveFunction(); + BSize MinSize(Variable* width, Variable* height); + BSize MaxSize(Variable* width, Variable* height); ResultType Solve(); bool Save(const char* fileName); int32 CountColumns() const; - OptimizationType Optimization() const; - void SetOptimization(OptimizationType value); ResultType Result() const; - double ObjectiveValue() const; - double SolvingTime() const; + bigtime_t 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; - double fSolvingTime; + bigtime_t fSolvingTime; SolverInterface* fSolver; }; diff --git a/headers/libs/linprog/PenaltyFunction.h b/headers/libs/linprog/PenaltyFunction.h deleted file mode 100644 index c18d6efa90..0000000000 --- a/headers/libs/linprog/PenaltyFunction.h +++ /dev/null @@ -1,50 +0,0 @@ -/* - * Copyright 2007-2008, Christof Lutteroth, lutteroth@cs.auckland.ac.nz - * Copyright 2007-2008, James Kim, jkim202@ec.auckland.ac.nz - * Distributed under the terms of the MIT License. - */ - -#ifndef PENALTY_FUNCTION_H -#define PENALTY_FUNCTION_H - -#include - - -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 - diff --git a/headers/libs/linprog/Summand.h b/headers/libs/linprog/Summand.h index 6370455bf0..92efe1cb06 100644 --- a/headers/libs/linprog/Summand.h +++ b/headers/libs/linprog/Summand.h @@ -27,6 +27,7 @@ public: Variable* Var(); void SetVar(Variable* var); + int32 VariableIndex(); private: double fCoeff; Variable* fVar; diff --git a/src/libs/alm/ALMLayout.cpp b/src/libs/alm/ALMLayout.cpp index 97e4711400..58183090cf 100644 --- a/src/libs/alm/ALMLayout.cpp +++ b/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); - 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()); + return fSolver->MinSize(Right(), Bottom()); } @@ -853,24 +817,7 @@ BALMLayout::_CalculateMaxSize() { _UpdateAreaConstraints(); - SummandList* newObjFunction = new(std::nothrow) SummandList(2); - 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()); + return fSolver->MaxSize(Right(), Bottom()); } @@ -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). " diff --git a/src/libs/alm/Area.cpp b/src/libs/alm/Area.cpp index 051302f3d8..4d3c2a6ede 100644 --- a/src/libs/alm/Area.cpp +++ b/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 } diff --git a/src/libs/linprog/ActiveSetSolver.cpp b/src/libs/linprog/ActiveSetSolver.cpp new file mode 100644 index 0000000000..b75256e0ce --- /dev/null +++ b/src/libs/linprog/ActiveSetSolver.cpp @@ -0,0 +1,542 @@ +#include "ActiveSetSolver.h" + +#include + +#include "LayoutOptimizer.h" + + +//#define DEBUG_ACTIVE_SOLVER + +#ifdef DEBUG_ACTIVE_SOLVER +#include +#define TRACE(x...) printf(x) +#else +#define TRACE(x...) /* nothing */ +#endif + + +using namespace LinearProgramming; +using namespace BPrivate::Layout; + + +template +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()); +} diff --git a/src/libs/linprog/ActiveSetSolver.h b/src/libs/linprog/ActiveSetSolver.h new file mode 100644 index 0000000000..aee37b0319 --- /dev/null +++ b/src/libs/linprog/ActiveSetSolver.h @@ -0,0 +1,82 @@ +/* + * Copyright 2010, Clemens Zeidler + * 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 diff --git a/src/libs/linprog/Constraint.cpp b/src/libs/linprog/Constraint.cpp index a87ed524b4..982b5e364b 100644 --- a/src/libs/linprog/Constraint.cpp +++ b/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 + 1; + + return i; } @@ -222,22 +221,7 @@ Constraint::SetPenaltyNeg(double value) { fPenaltyNeg = value; - if (!fIsValid) - 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(); + fLS->UpdateLeftSide(this); } @@ -263,22 +247,7 @@ Constraint::SetPenaltyPos(double value) { fPenaltyPos = value; - if (!fIsValid) - 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(); + fLS->UpdateLeftSide(this); } @@ -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\n")); - file->Write(buffer, sprintf(buffer, "\t\t\n")); - - Summand* summand; - for (int32 i = 0; i < fLeftSide->CountItems(); i++) { - summand = (Summand*)fLeftSide->ItemAt(i); - file->Write(buffer, sprintf(buffer, "\t\t\t\n")); - file->Write(buffer, sprintf(buffer, "\t\t\t\t%f\n", - summand->Coeff())); - BString varStr = *(summand->Var()); - file->Write(buffer, sprintf(buffer, "\t\t\t\t%s\n", - varStr.String())); - file->Write(buffer, sprintf(buffer, "\t\t\t\n")); - } - - file->Write(buffer, sprintf(buffer, "\t\t\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%s\n", op)); - file->Write(buffer, sprintf(buffer, "\t\t%f\n", fRightSide)); - //~ file->Write(buffer, sprintf(buffer, "\t\t%s\n", PenaltyNeg())); - //~ file->Write(buffer, sprintf(buffer, "\t\t%s\n", PenaltyPos())); - file->Write(buffer, sprintf(buffer, "\t\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(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. */ diff --git a/src/libs/linprog/Jamfile b/src/libs/linprog/Jamfile index 8ac0101090..b11568af18 100644 --- a/src/libs/linprog/Jamfile +++ b/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 ; diff --git a/src/libs/linprog/LPSolveInterface.cpp b/src/libs/linprog/LPSolveInterface.cpp index e22a715246..e9f46d14ac 100644 --- a/src/libs/linprog/LPSolveInterface.cpp +++ b/src/libs/linprog/LPSolveInterface.cpp @@ -8,14 +8,20 @@ #include "LPSolveInterface.h" +#include + 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) { - int32 size = variables.CountItems(); - double x[size]; - if (!get_variables(fLP, &x[0])) - printf("Error in get_variables.\n"); + if (result == OPTIMAL) { + int32 size = variables.CountItems(); + double x[size]; + if (!get_variables(fLP, &x[0])) + printf("Error in get_variables.\n"); - for (int32 i = 0; i < size; i++) - variables.ItemAt(i)->SetValue(x[i]); + for (int32 i = 0; i < size; 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, - int* variableIndices, OperatorType op, double rightSide) +LPSolveInterface::ConstraintAdded(Constraint* constraint) { - 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, - double* coefficients, int* variableIndices) +LPSolveInterface::OperatorChanged(Constraint* constraint) { - if (!set_rowex(fLP, constraint, nElements, coefficients, variableIndices)) + int32 index = constraint->Index() + 1; + if (index <= 0) return false; - _RemovePresolved(); - return true; -} - -bool -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) { + OperatorType op = constraint->Op(); + if (!set_constr_type(fLP, index, 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)); } } diff --git a/src/libs/linprog/LPSolveInterface.h b/src/libs/linprog/LPSolveInterface.h index 0e3a85e13e..874f205590 100644 --- a/src/libs/linprog/LPSolveInterface.h +++ b/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); - double GetObjectiveValue(); + ResultType Solve(); - bool AddVariable(); - bool RemoveVariable(int variable); - bool SetVariableRange(int variable, double min, - double max); + bool VariableAdded(Variable* variable); + bool VariableRemoved(Variable* variable); + bool VariableRangeChanged(Variable* variable); - bool AddConstraint(int nElements, - double* coefficients, int* variableIndices, - OperatorType op, double rightSide); - bool RemoveConstraint(int constraint); - bool SetLeftSide(int constraint, int nElements, - 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 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); + + 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 diff --git a/src/libs/linprog/LayoutOptimizer.cpp b/src/libs/linprog/LayoutOptimizer.cpp new file mode 100644 index 0000000000..1f36cad608 --- /dev/null +++ b/src/libs/linprog/LayoutOptimizer.cpp @@ -0,0 +1,940 @@ +/* + * Copyright 2007, Ingo Weinhold . + * Copyright 2010, Clemens Zeidler + * Distributed under the terms of the MIT License. + */ + + +#include "LayoutOptimizer.h" + +#include +#include +#include + +#include + + +//#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 +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 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]; +} diff --git a/src/libs/linprog/LayoutOptimizer.h b/src/libs/linprog/LayoutOptimizer.h new file mode 100644 index 0000000000..1dc1aff494 --- /dev/null +++ b/src/libs/linprog/LayoutOptimizer.h @@ -0,0 +1,85 @@ +/* + * Copyright 2007, Ingo Weinhold . + * All rights reserved. Distributed under the terms of the MIT License. + */ +#ifndef LAYOUT_OPTIMIZER_H +#define LAYOUT_OPTIMIZER_H + +#include +#include + +#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 diff --git a/src/libs/linprog/LinearSpec.cpp b/src/libs/linprog/LinearSpec.cpp index fb80b50440..74a819eed6 100644 --- a/src/libs/linprog/LinearSpec.cpp +++ b/src/libs/linprog/LinearSpec.cpp @@ -12,6 +12,28 @@ #include #include "LPSolveInterface.h" +#include "ActiveSetSolver.h" + + +using namespace LinearProgramming; + + +#define DEBUG_LINEAR_SPECIFICATIONS + +#ifdef DEBUG_LINEAR_SPECIFICATIONS +#include +#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(NAN) + fSolvingTime(0) { - 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); - if (index < 0) + // must be called first otherwise the index is invalid + if (!fSolver->VariableRemoved(variable)) return false; - - if (!fSolver->RemoveVariable(index)) - return false; - fVariables.RemoveItemAt(index - 1); + fVariables.RemoveItem(variable); 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); - if (i == -1) { - printf("Variable 0x%p not part of fLS->Variables().\n", variable); - return -1; - } - return i + 1; + return fVariables.IndexOf(variable); } bool LinearSpec::UpdateRange(Variable* variable) { - if (!fSolver->SetVariableRange(IndexOf(variable), variable->Min(), - variable->Max())) + if (!fSolver->VariableRangeChanged(variable)) return false; return true; } @@ -157,49 +165,13 @@ LinearSpec::UpdateRange(Variable* variable) bool LinearSpec::AddConstraint(Constraint* constraint) { - SummandList* summands = constraint->LeftSide(); - OperatorType op = constraint->Op(); - double rightSide = constraint->RightSide(); - double penaltyNeg = constraint->PenaltyNeg(); - double penaltyPos = constraint->PenaltyPos(); + if (!fConstraints.AddItem(constraint)) + return false; - double coeffs[summands->CountItems() + 2]; - int varIndexes[summands->CountItems() + 2]; - 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)) { + if (!fSolver->ConstraintAdded(constraint)) { + fConstraints.RemoveItem(constraint); 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) { - 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()); + fSolver->ConstraintRemoved(constraint); 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, - INFINITY); + return AddConstraint(coeff1, var1, coeff2, var2, op, rightSide, -1, + -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, } -/** - * Adds a new penalty function to the specification. - * - * @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) +BSize +LinearSpec::MinSize(Variable* width, Variable* height) { - return new(std::nothrow) PenaltyFunction(this, var, xs, gs); + return fSolver->MinSize(width, height); } -/** - * Gets the objective function. - * - * @return SummandList containing the objective function's summands - */ -SummandList* -LinearSpec::ObjectiveFunction() +BSize +LinearSpec::MaxSize(Variable* width, Variable* height) { - return fObjFunction; -} - - -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"); + return fSolver->MaxSize(width, height); } @@ -638,23 +505,27 @@ LinearSpec::_AddConstraint(SummandList* leftSide, OperatorType op, } -/** - * Tries to solve the linear programming problem. - * If a cached simplified version of the problem exists, it is used instead. - * - * @return the result of the solving attempt - */ +#ifdef DEBUG_LINEAR_SPECIFICATIONS +static bigtime_t sAverageSolvingTime = 0; +static int32 sSolvedCount = 0; +#endif + + ResultType LinearSpec::Solve() { - bigtime_t start, end; - start = system_time(); + bigtime_t startTime = system_time(); - fResult = fSolver->Solve(fVariables); - fObjectiveValue = fSolver->GetObjectiveValue(); + fResult = fSolver->Solve(); - end = system_time(); - fSolvingTime = (end - start) / 1000.0; + fSolvingTime = system_time() - startTime; + +#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"; } diff --git a/src/libs/linprog/PenaltyFunction.cpp b/src/libs/linprog/PenaltyFunction.cpp deleted file mode 100644 index ac6071f2ed..0000000000 --- a/src/libs/linprog/PenaltyFunction.cpp +++ /dev/null @@ -1,70 +0,0 @@ -/* - * Copyright 2007-2008, Christof Lutteroth, lutteroth@cs.auckland.ac.nz - * Copyright 2007-2008, James Kim, jkim202@ec.auckland.ac.nz - * Distributed under the terms of the MIT License. - */ - - -#include "PenaltyFunction.h" - -#include - -#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); -} - diff --git a/src/libs/linprog/Summand.cpp b/src/libs/linprog/Summand.cpp index 749bc72d9b..1e98e10219 100644 --- a/src/libs/linprog/Summand.cpp +++ b/src/libs/linprog/Summand.cpp @@ -57,6 +57,13 @@ Summand::SetVar(Variable* var) } +int32 +Summand::VariableIndex() +{ + return fVar->Index(); +} + + /** * Destructor. */