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