1/*
2 * Copyright 2010 - 2011, Haiku.
3 * Distributed under the terms of the MIT License.
4 *
5 * Authors:
6 *		Clemens Zeidler <haiku@clemens-zeidler.de>
7 */
8
9
10#include "ActiveSetSolver.h"
11
12#include <algorithm>
13#include <stdio.h>
14
15#include "LayoutOptimizer.h"
16
17
18// #define DEBUG_ACTIVE_SOLVER
19
20#ifdef DEBUG_ACTIVE_SOLVER
21#define TRACE(x...) printf(x)
22#else
23#define TRACE(x...) /* nothing */
24#endif
25
26
27using namespace LinearProgramming;
28using namespace std;
29
30
31EquationSystem::EquationSystem(int32 rows, int32 columns)
32	:
33	fRows(rows),
34	fColumns(columns)
35{
36	fMatrix = allocate_matrix(fRows, fColumns);
37	fB = new double[fColumns];
38	// better init all values to prevent side cases where not all variables
39	// needed to solve the problem, coping theses values to the results could
40	// cause problems
41	for (int i = 0; i < fColumns; i++)
42		fB[i] = 0;
43	zero_matrix(fMatrix, fRows, fColumns);
44
45	fRowIndices = new int32[fRows];
46	fColumnIndices = new int32[fColumns];
47	for (int i = 0; i < fRows; i++)
48		fRowIndices[i] = i;
49	for (int i = 0; i < fColumns; i++)
50		fColumnIndices[i] = i;
51}
52
53
54EquationSystem::~EquationSystem()
55{
56	free_matrix(fMatrix);
57	delete[] fB;
58	delete[] fRowIndices;
59	delete[] fColumnIndices;
60}
61
62
63void
64EquationSystem::SetRows(int32 rows)
65{
66	fRows = rows;
67}
68
69
70int32
71EquationSystem::Rows()
72{
73	return fRows;
74}
75
76
77int32
78EquationSystem::Columns()
79{
80	return fColumns;
81}
82
83
84bool
85EquationSystem::InRange(int32 row, int32 column)
86{
87	if (row < 0 || row >= fRows)
88		return false;
89	if (column < 0 || column >= fColumns)
90		return false;
91	return true;
92}
93
94
95double&
96EquationSystem::A(int32 row, int32 column)
97{
98	return fMatrix[fRowIndices[row]][fColumnIndices[column]];
99}
100
101
102double&
103EquationSystem::B(int32 row)
104{
105	return fB[row];
106}
107
108
109void
110EquationSystem::Results(double* results, int32 size)
111{
112	for (int i = 0; i < size; i++)
113		results[i] = 0;
114	for (int i = 0; i < fColumns; i++) {
115		int32 index = fColumnIndices[i];
116		if (index < fRows)
117			results[index] = fB[i];
118	}
119}
120
121
122void
123EquationSystem::SwapColumn(int32 i, int32 j)
124{
125	swap(fColumnIndices[i], fColumnIndices[j]);
126}
127
128
129void
130EquationSystem::SwapRow(int32 i, int32 j)
131{
132	swap(fRowIndices[i], fRowIndices[j]);
133	swap(fB[i], fB[j]);
134}
135
136
137bool
138EquationSystem::GaussianElimination()
139{
140	// basic solve
141	for (int i = 0; i < fRows; i++) {
142		// find none zero element
143		int swapRow = -1;
144		for (int r = i; r < fRows; r++) {
145			double& value = fMatrix[fRowIndices[r]][fColumnIndices[i]];
146			if (fuzzy_equals(value, 0))
147				continue;
148			swapRow = r;
149			break;
150		}
151		if (swapRow == -1) {
152			int swapColumn = -1;
153			for (int c = i + 1; c < fColumns; c++) {
154				double& value = fMatrix[fRowIndices[i]][fColumnIndices[c]];
155				if (fuzzy_equals(value, 0))
156					continue;
157				swapRow = i;
158				swapColumn = c;
159				break;
160			}
161			if (swapColumn == -1) {
162				TRACE("can't solve column %i\n", i);
163				return false;
164			}
165			SwapColumn(i, swapColumn);
166		}
167		if (i != swapRow)
168			SwapRow(i, swapRow);
169
170		// normalize
171		_EliminateColumn(i, i + 1, fRows - 1);
172	}
173	return true;
174}
175
176
177bool
178EquationSystem::GaussJordan()
179{
180	if (!GaussianElimination())
181		return false;
182
183	for (int32 i = fRows - 1; i >= 0; i--)
184		_EliminateColumn(i, 0, i - 1);
185
186	return true;
187}
188
189
190void
191EquationSystem::GaussJordan(int32 i)
192{
193	_EliminateColumn(i, 0, fRows - 1);
194}
195
196
197void
198EquationSystem::RemoveLinearlyDependentRows()
199{
200	double** temp = allocate_matrix(fRows, fColumns);
201	bool independentRows[fRows];
202
203	// copy to temp
204	copy_matrix(fMatrix, temp, fRows, fColumns);
205	int nIndependent = compute_dependencies(temp, fRows, fColumns,
206		independentRows);
207	if (nIndependent == fRows)
208		return;
209
210	// remove the rows
211	for (int i = 0; i < fRows; i++) {
212		if (!independentRows[i]) {
213			int lastDepRow = -1;
214			for (int d = fRows - 1; d > i; d--) {
215				if (independentRows[d]) {
216					lastDepRow = d;
217					break;
218				}
219			}
220			if (lastDepRow < 0)
221				break;
222			SwapRow(i, lastDepRow);
223			fRows--;
224		}
225	}
226	fRows = nIndependent;
227
228	free_matrix(temp);
229}
230
231
232void
233EquationSystem::RemoveUnusedVariables()
234{
235	for (int c = 0; c < fColumns; c++) {
236		bool used = false;
237		for (int r = 0; r < fRows; r++) {
238			if (!fuzzy_equals(fMatrix[r][fColumnIndices[c]], 0)) {
239				used = true;
240				break;
241			}
242		}
243		if (used)
244			continue;
245
246		//MoveColumnRight(c, fColumns - 1);
247		SwapColumn(c, fColumns - 1);
248		fColumns--;
249		c--;
250	}
251}
252
253
254void
255EquationSystem::MoveColumnRight(int32 i, int32 target)
256{
257	int32 index = fColumnIndices[i];
258	for (int c = i; c < target; c++)
259		fColumnIndices[c] = fColumnIndices[c + 1];
260	fColumnIndices[target] = index;
261}
262
263
264void
265EquationSystem::Print()
266{
267	for (int m = 0; m < fRows; m++) {
268		for (int n = 0; n < fColumns; n++)
269			printf("%.1f ", fMatrix[fRowIndices[m]][fColumnIndices[n]]);
270		printf("= %.1f\n", fB[m]);
271	}
272}
273
274
275void
276EquationSystem::_EliminateColumn(int32 column, int32 startRow, int32 endRow)
277{
278	double value = fMatrix[fRowIndices[column]][fColumnIndices[column]];
279	if (value != 1.) {
280		for (int j = column; j < fColumns; j++)
281			fMatrix[fRowIndices[column]][fColumnIndices[j]] /= value;
282		fB[column] /= value;
283	}
284
285	for (int r = startRow; r < endRow + 1; r++) {
286		if (r == column)
287			continue;
288		double q = -fMatrix[fRowIndices[r]][fColumnIndices[column]];
289		// don't need to do anything, since matrix is typically sparse this
290		// should save some work
291		if (fuzzy_equals(q, 0))
292			continue;
293		for (int c = column; c < fColumns; c++)
294			fMatrix[fRowIndices[r]][fColumnIndices[c]]
295				+= fMatrix[fRowIndices[column]][fColumnIndices[c]] * q;
296
297		fB[r] += fB[column] * q;
298	}
299}
300
301
302namespace {
303
304
305Constraint*
306AddMinConstraint(LinearSpec* spec, Variable* var)
307{
308	return spec->AddConstraint(1, var, kEQ, 0, 5, 5);
309}
310
311
312Constraint*
313AddMaxConstraint(LinearSpec* spec, Variable* var)
314{
315	static const double kHugeValue = 32000;
316	return spec->AddConstraint(1, var, kEQ, kHugeValue, 5, 5);
317}
318
319};
320
321
322ActiveSetSolver::ActiveSetSolver(LinearSpec* linearSpec)
323	:
324	SolverInterface(linearSpec),
325
326	fVariables(linearSpec->UsedVariables()),
327	fConstraints(linearSpec->Constraints())
328{
329
330}
331
332
333ActiveSetSolver::~ActiveSetSolver()
334{
335
336}
337
338
339/* Using algorithm found in:
340Solving Inequalities and Proving Farkas's Lemma Made Easy
341David Avis and Bohdan Kaluzny
342The American Mathematical Monthly
343Vol. 111, No. 2 (Feb., 2004), pp. 152-157 */
344static bool
345solve(EquationSystem& system)
346{
347	// basic solve
348	if (!system.GaussJordan())
349		return false;
350
351	bool done = false;
352	while (!done) {
353		double smallestB = HUGE_VALF;
354		int smallestBRow = -1;
355		for (int row = 0; row < system.Rows(); row++) {
356			if (system.B(row) > 0 || fuzzy_equals(system.B(row), 0))
357				continue;
358
359			double bValue = fabs(system.B(row));
360			if (bValue < smallestB) {
361				smallestB = bValue;
362				smallestBRow = row;
363			}
364		}
365		if (smallestBRow == -1) {
366			done = true;
367			break;
368		}
369
370		int negValueCol = -1;
371		for (int col = system.Rows(); col < system.Columns(); col++) {
372			double value = system.A(smallestBRow, col);
373			if (value > 0 || fuzzy_equals(value, 0))
374				continue;
375			negValueCol = col;
376			break;
377		}
378		if (negValueCol == -1) {
379			TRACE("can't solve\n");
380			return false;
381		}
382
383		system.SwapColumn(smallestBRow, negValueCol);
384
385		// eliminate
386		system.GaussJordan(smallestBRow);
387	}
388
389	return true;
390}
391
392
393static bool is_soft_inequality(Constraint* constraint)
394{
395	if (constraint->PenaltyNeg() <= 0. && constraint->PenaltyPos() <= 0.)
396		return false;
397	if (constraint->Op() != kEQ)
398		return true;
399	return false;
400}
401
402
403class SoftInequalityAdder {
404public:
405	SoftInequalityAdder(LinearSpec* linSpec, ConstraintList& allConstraints)
406		:
407		fLinearSpec(linSpec)
408	{
409		for (int32 c = 0; c < allConstraints.CountItems(); c++) {
410			Constraint* constraint = allConstraints.ItemAt(c);
411			if (!is_soft_inequality(constraint))
412				continue;
413
414			Variable* variable = fLinearSpec->AddVariable();
415			variable->SetRange(0, 20000);
416
417			Constraint* helperConst = fLinearSpec->AddConstraint(1, variable,
418				kEQ, 0, constraint->PenaltyNeg(), constraint->PenaltyPos());
419			fInequalitySoftConstraints.AddItem(helperConst);
420
421			double coeff = -1;
422			if (constraint->Op() == kGE)
423				coeff = 1;
424
425			Constraint* modifiedConstraint = new Constraint(constraint);
426			allConstraints.AddItem(modifiedConstraint, c);
427			allConstraints.RemoveItemAt(c + 1);
428			fModifiedIneqConstraints.AddItem(modifiedConstraint);
429			modifiedConstraint->LeftSide()->AddItem(
430				new Summand(coeff, variable));
431		}
432		for (int32 c = 0; c < fInequalitySoftConstraints.CountItems(); c++)
433			allConstraints.AddItem(fInequalitySoftConstraints.ItemAt(c));
434	}
435
436	~SoftInequalityAdder()
437	{
438		for (int32 c = 0; c < fModifiedIneqConstraints.CountItems(); c++)
439			delete fModifiedIneqConstraints.ItemAt(c);
440		for (int32 c = 0; c < fInequalitySoftConstraints.CountItems(); c++) {
441			Constraint* con = fInequalitySoftConstraints.ItemAt(c);
442			fLinearSpec->RemoveVariable(con->LeftSide()->ItemAt(0)->Var());
443				// this also deletes the constraint
444		}
445	}
446
447private:
448	LinearSpec*		fLinearSpec;
449	ConstraintList	fModifiedIneqConstraints;
450	ConstraintList	fInequalitySoftConstraints;
451};
452
453
454ResultType
455ActiveSetSolver::Solve()
456{
457
458	// make a copy of the original constraints and create soft inequality
459	// constraints
460	ConstraintList allConstraints(fConstraints);
461	SoftInequalityAdder adder(fLinearSpec, allConstraints);
462
463	int32 nConstraints = allConstraints.CountItems();
464	int32 nVariables = fVariables.CountItems();
465
466	if (nVariables > nConstraints) {
467		TRACE("More variables then constraints! vars: %i, constraints: %i\n",
468			(int)nVariables, (int)nConstraints);
469		return kInfeasible;
470	}
471
472	/* First find an initial solution and then optimize it using the active set
473	method. */
474	EquationSystem system(nConstraints, nVariables + nConstraints);
475
476	int32 slackIndex = nVariables;
477	// setup constraint matrix and add slack variables if necessary
478	int32 rowIndex = 0;
479	for (int32 c = 0; c < nConstraints; c++) {
480		Constraint* constraint = allConstraints.ItemAt(c);
481		if (is_soft(constraint))
482			continue;
483		SummandList* leftSide = constraint->LeftSide();
484		system.B(rowIndex) = constraint->RightSide();
485		for (int32 sIndex = 0; sIndex < leftSide->CountItems(); sIndex++ ) {
486			Summand* summand = leftSide->ItemAt(sIndex);
487			double coefficient = summand->Coeff();
488			int32 columnIndex = summand->VariableIndex();
489			system.A(rowIndex, columnIndex) = coefficient;
490		}
491		if (constraint->Op() == kLE) {
492			system.A(rowIndex, slackIndex) = 1;
493			slackIndex++;
494		} else if (constraint->Op() == kGE) {
495			system.A(rowIndex, slackIndex) = -1;
496			slackIndex++;
497		}
498		rowIndex++;
499	}
500
501	system.SetRows(rowIndex);
502
503	system.RemoveLinearlyDependentRows();
504	system.RemoveUnusedVariables();
505
506	if (!solve(system))
507		return kInfeasible;
508
509	double results[nVariables + nConstraints];
510	system.Results(results, nVariables + nConstraints);
511	TRACE("base system solved\n");
512
513	LayoutOptimizer optimizer(allConstraints, nVariables);
514	if (!optimizer.Solve(results))
515		return kInfeasible;
516
517	// back to the variables
518	for (int32 i = 0; i < nVariables; i++)
519		fVariables.ItemAt(i)->SetValue(results[i]);
520
521	for (int32 i = 0; i < nVariables; i++)
522		TRACE("var %f\n", results[i]);
523
524	return kOptimal;
525}
526
527
528bool
529ActiveSetSolver::VariableAdded(Variable* variable)
530{
531	if (fVariableGEConstraints.AddItem(NULL) == false)
532		return false;
533	if (fVariableLEConstraints.AddItem(NULL) == false) {
534		// clean up
535		int32 count = fVariableGEConstraints.CountItems();
536		fVariableGEConstraints.RemoveItemAt(count - 1);
537		return false;
538	}
539	return true;
540}
541
542
543bool
544ActiveSetSolver::VariableRemoved(Variable* variable)
545{
546	fVariableGEConstraints.RemoveItemAt(variable->GlobalIndex());
547	fVariableLEConstraints.RemoveItemAt(variable->GlobalIndex());
548	return true;
549}
550
551
552bool
553ActiveSetSolver::VariableRangeChanged(Variable* variable)
554{
555	double min = variable->Min();
556	double max = variable->Max();
557	int32 variableIndex = variable->GlobalIndex();
558
559	Constraint* constraintGE = fVariableGEConstraints.ItemAt(variableIndex);
560	Constraint* constraintLE = fVariableLEConstraints.ItemAt(variableIndex);
561	if (constraintGE == NULL && min > -20000) {
562		constraintGE = fLinearSpec->AddConstraint(1, variable, kGE, 0);
563		if (constraintGE == NULL)
564			return false;
565		constraintGE->SetLabel("Var Min");
566		fVariableGEConstraints.RemoveItemAt(variableIndex);
567		fVariableGEConstraints.AddItem(constraintGE, variableIndex);
568	}
569	if (constraintLE == NULL && max < 20000) {
570		constraintLE = fLinearSpec->AddConstraint(1, variable, kLE, 20000);
571		if (constraintLE == NULL)
572			return false;
573		constraintLE->SetLabel("Var Max");
574		fVariableLEConstraints.RemoveItemAt(variableIndex);
575		fVariableLEConstraints.AddItem(constraintLE, variableIndex);
576	}
577
578	if (constraintGE)
579		constraintGE->SetRightSide(min);
580	if (constraintLE)
581		constraintLE->SetRightSide(max);
582	return true;
583}
584
585
586bool
587ActiveSetSolver::ConstraintAdded(Constraint* constraint)
588{
589	return true;
590}
591
592
593bool
594ActiveSetSolver::ConstraintRemoved(Constraint* constraint)
595{
596	return true;
597}
598
599
600bool
601ActiveSetSolver::LeftSideChanged(Constraint* constraint)
602{
603	return true;
604}
605
606
607bool
608ActiveSetSolver::RightSideChanged(Constraint* constraint)
609{
610	return true;
611}
612
613
614bool
615ActiveSetSolver::OperatorChanged(Constraint* constraint)
616{
617	return true;
618}
619
620
621bool
622ActiveSetSolver::PenaltiesChanged(Constraint* constraint)
623{
624	return true;
625}
626
627
628bool
629ActiveSetSolver::SaveModel(const char* fileName)
630{
631	return false;
632}
633
634
635ResultType
636ActiveSetSolver::FindMaxs(const VariableList* variables)
637{
638	return _FindWithConstraintsNoSoft(variables, AddMaxConstraint);
639}
640
641
642ResultType
643ActiveSetSolver::FindMins(const VariableList* variables)
644{
645	return _FindWithConstraintsNoSoft(variables, AddMinConstraint);
646}
647
648
649void
650ActiveSetSolver::GetRangeConstraints(const Variable* var,
651	const Constraint** _min, const Constraint** _max) const
652{
653	int32 variableIndex = var->GlobalIndex();
654
655	if (_min)
656		*_min = fVariableGEConstraints.ItemAt(variableIndex);
657	if (_max)
658		*_max = fVariableLEConstraints.ItemAt(variableIndex);
659}
660
661
662void
663ActiveSetSolver::_RemoveSoftConstraint(ConstraintList& list)
664{
665	ConstraintList allConstraints = fLinearSpec->Constraints();
666	for (int i = 0; i < allConstraints.CountItems(); i++) {
667		Constraint* constraint = allConstraints.ItemAt(i);
668		// soft eq an ineq constraint?
669		if (constraint->PenaltyNeg() <= 0. && constraint->PenaltyPos() <= 0.)
670			continue;
671
672		if (RemoveConstraint(constraint, false, false) == true)
673			list.AddItem(constraint);
674	}
675}
676
677
678void
679ActiveSetSolver::_AddSoftConstraint(const ConstraintList& list)
680{
681	for (int i = 0; i < list.CountItems(); i++) {
682		Constraint* constraint = list.ItemAt(i);
683		// at least don't leak it
684		if (AddConstraint(constraint, false) == false)
685			delete constraint;
686	}
687}
688
689
690ResultType
691ActiveSetSolver::_FindWithConstraintsNoSoft(const VariableList* variables,
692	AddConstraintFunc constraintFunc)
693{
694	ConstraintList softConstraints;
695	_RemoveSoftConstraint(softConstraints);
696
697	ConstraintList constraints;
698	for (int32 i = 0; i < variables->CountItems(); i++)
699		constraints.AddItem(constraintFunc(fLinearSpec, variables->ItemAt(i)));
700
701	ResultType result = Solve();
702	for (int32 i = 0; i < constraints.CountItems(); i++)
703		fLinearSpec->RemoveConstraint(constraints.ItemAt(i));
704
705	_AddSoftConstraint(softConstraints);
706
707	if (result != kOptimal)
708		TRACE("Could not solve the layout specification (%d). ", result);
709
710	return result;
711}
712