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