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 27 using namespace LinearProgramming; 28 using namespace std; 29 30 31 EquationSystem::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 54 EquationSystem::~EquationSystem() 55 { 56 free_matrix(fMatrix); 57 delete[] fB; 58 delete[] fRowIndices; 59 delete[] fColumnIndices; 60 } 61 62 63 void 64 EquationSystem::SetRows(int32 rows) 65 { 66 fRows = rows; 67 } 68 69 70 int32 71 EquationSystem::Rows() 72 { 73 return fRows; 74 } 75 76 77 int32 78 EquationSystem::Columns() 79 { 80 return fColumns; 81 } 82 83 84 bool 85 EquationSystem::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 95 double& 96 EquationSystem::A(int32 row, int32 column) 97 { 98 return fMatrix[fRowIndices[row]][fColumnIndices[column]]; 99 } 100 101 102 double& 103 EquationSystem::B(int32 row) 104 { 105 return fB[row]; 106 } 107 108 109 void 110 EquationSystem::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 122 void 123 EquationSystem::SwapColumn(int32 i, int32 j) 124 { 125 swap(fColumnIndices[i], fColumnIndices[j]); 126 } 127 128 129 void 130 EquationSystem::SwapRow(int32 i, int32 j) 131 { 132 swap(fRowIndices[i], fRowIndices[j]); 133 swap(fB[i], fB[j]); 134 } 135 136 137 bool 138 EquationSystem::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 177 bool 178 EquationSystem::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 190 void 191 EquationSystem::GaussJordan(int32 i) 192 { 193 _EliminateColumn(i, 0, fRows - 1); 194 } 195 196 197 void 198 EquationSystem::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 232 void 233 EquationSystem::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 254 void 255 EquationSystem::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 264 void 265 EquationSystem::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 275 void 276 EquationSystem::_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 302 namespace { 303 304 305 Constraint* 306 AddMinConstraint(LinearSpec* spec, Variable* var) 307 { 308 return spec->AddConstraint(1, var, kEQ, 0, 5, 5); 309 } 310 311 312 Constraint* 313 AddMaxConstraint(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 322 ActiveSetSolver::ActiveSetSolver(LinearSpec* linearSpec) 323 : 324 SolverInterface(linearSpec), 325 326 fVariables(linearSpec->UsedVariables()), 327 fConstraints(linearSpec->Constraints()) 328 { 329 330 } 331 332 333 ActiveSetSolver::~ActiveSetSolver() 334 { 335 336 } 337 338 339 /* Using algorithm found in: 340 Solving Inequalities and Proving Farkas's Lemma Made Easy 341 David Avis and Bohdan Kaluzny 342 The American Mathematical Monthly 343 Vol. 111, No. 2 (Feb., 2004), pp. 152-157 */ 344 static bool 345 solve(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 393 static 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 403 class SoftInequalityAdder { 404 public: 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 447 private: 448 LinearSpec* fLinearSpec; 449 ConstraintList fModifiedIneqConstraints; 450 ConstraintList fInequalitySoftConstraints; 451 }; 452 453 454 ResultType 455 ActiveSetSolver::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 528 bool 529 ActiveSetSolver::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 543 bool 544 ActiveSetSolver::VariableRemoved(Variable* variable) 545 { 546 fVariableGEConstraints.RemoveItemAt(variable->GlobalIndex()); 547 fVariableLEConstraints.RemoveItemAt(variable->GlobalIndex()); 548 return true; 549 } 550 551 552 bool 553 ActiveSetSolver::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 586 bool 587 ActiveSetSolver::ConstraintAdded(Constraint* constraint) 588 { 589 return true; 590 } 591 592 593 bool 594 ActiveSetSolver::ConstraintRemoved(Constraint* constraint) 595 { 596 return true; 597 } 598 599 600 bool 601 ActiveSetSolver::LeftSideChanged(Constraint* constraint) 602 { 603 return true; 604 } 605 606 607 bool 608 ActiveSetSolver::RightSideChanged(Constraint* constraint) 609 { 610 return true; 611 } 612 613 614 bool 615 ActiveSetSolver::OperatorChanged(Constraint* constraint) 616 { 617 return true; 618 } 619 620 621 bool 622 ActiveSetSolver::PenaltiesChanged(Constraint* constraint) 623 { 624 return true; 625 } 626 627 628 bool 629 ActiveSetSolver::SaveModel(const char* fileName) 630 { 631 return false; 632 } 633 634 635 ResultType 636 ActiveSetSolver::FindMaxs(const VariableList* variables) 637 { 638 return _FindWithConstraintsNoSoft(variables, AddMaxConstraint); 639 } 640 641 642 ResultType 643 ActiveSetSolver::FindMins(const VariableList* variables) 644 { 645 return _FindWithConstraintsNoSoft(variables, AddMinConstraint); 646 } 647 648 649 void 650 ActiveSetSolver::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 662 void 663 ActiveSetSolver::_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 678 void 679 ActiveSetSolver::_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 690 ResultType 691 ActiveSetSolver::_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