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 printf("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 ActiveSetSolver::ActiveSetSolver(LinearSpec* linearSpec) 297 : 298 QPSolverInterface(linearSpec), 299 300 fVariables(linearSpec->UsedVariables()), 301 fConstraints(linearSpec->Constraints()) 302 { 303 304 } 305 306 307 ActiveSetSolver::~ActiveSetSolver() 308 { 309 310 } 311 312 313 /* Using algorithm found in: 314 Solving Inequalities and Proving Farkas's Lemma Made Easy 315 David Avis and Bohdan Kaluzny 316 The American Mathematical Monthly 317 Vol. 111, No. 2 (Feb., 2004), pp. 152-157 */ 318 static bool 319 solve(EquationSystem& system) 320 { 321 // basic solve 322 if (!system.GaussJordan()) 323 return false; 324 325 bool done = false; 326 while (!done) { 327 double smallestB = HUGE_VALF; 328 int smallestBRow = -1; 329 for (int row = 0; row < system.Rows(); row++) { 330 if (system.B(row) > 0 || fuzzy_equals(system.B(row), 0)) 331 continue; 332 333 double bValue = fabs(system.B(row)); 334 if (bValue < smallestB) { 335 smallestB = bValue; 336 smallestBRow = row; 337 } 338 } 339 if (smallestBRow == -1) { 340 done = true; 341 break; 342 } 343 344 int negValueCol = -1; 345 for (int col = system.Rows(); col < system.Columns(); col++) { 346 double value = system.A(smallestBRow, col); 347 if (value > 0 || fuzzy_equals(value, 0)) 348 continue; 349 negValueCol = col; 350 break; 351 } 352 if (negValueCol == -1) { 353 printf("can't solve\n"); 354 return false; 355 } 356 357 system.SwapColumn(smallestBRow, negValueCol); 358 359 // eliminate 360 system.GaussJordan(smallestBRow); 361 } 362 363 return true; 364 } 365 366 367 ResultType 368 ActiveSetSolver::Solve() 369 { 370 int32 nConstraints = fConstraints.CountItems(); 371 int32 nVariables = fVariables.CountItems(); 372 373 if (nVariables > nConstraints) { 374 printf("More variables then constraints! vars: %i, constraints: %i\n", 375 (int)nVariables, (int)nConstraints); 376 return kInfeasible; 377 } 378 379 /* First find an initial solution and then optimize it using the active set 380 method. */ 381 EquationSystem system(nConstraints, nVariables + nConstraints); 382 383 int32 slackIndex = nVariables; 384 // setup constraint matrix and add slack variables if necessary 385 int32 rowIndex = 0; 386 for (int32 c = 0; c < nConstraints; c++) { 387 Constraint* constraint = fConstraints.ItemAt(c); 388 if (constraint->IsSoft()) 389 continue; 390 SummandList* leftSide = constraint->LeftSide(); 391 system.B(rowIndex) = constraint->RightSide(); 392 for (int32 sIndex = 0; sIndex < leftSide->CountItems(); sIndex++ ) { 393 Summand* summand = leftSide->ItemAt(sIndex); 394 double coefficient = summand->Coeff(); 395 system.A(rowIndex, summand->VariableIndex()) = coefficient; 396 } 397 if (constraint->Op() == kLE) { 398 system.A(rowIndex, slackIndex) = 1; 399 slackIndex++; 400 } else if (constraint->Op() == kGE) { 401 system.A(rowIndex, slackIndex) = -1; 402 slackIndex++; 403 } 404 rowIndex++; 405 } 406 407 system.SetRows(rowIndex); 408 409 system.RemoveLinearlyDependentRows(); 410 system.RemoveUnusedVariables(); 411 412 if (!solve(system)) 413 return kInfeasible; 414 415 double results[nVariables + nConstraints]; 416 system.Results(results, nVariables + nConstraints); 417 TRACE("base system solved\n"); 418 419 LayoutOptimizer optimizer(fConstraints, nVariables); 420 optimizer.Solve(results); 421 422 // back to the variables 423 for (int32 i = 0; i < nVariables; i++) 424 fVariables.ItemAt(i)->SetValue(results[i]); 425 426 for (int32 i = 0; i < nVariables; i++) 427 TRACE("var %f\n", results[i]); 428 429 return kOptimal; 430 } 431 432 433 bool 434 ActiveSetSolver::VariableAdded(Variable* variable) 435 { 436 // TODO: error checks 437 fVariableGEConstraints.AddItem(NULL); 438 fVariableLEConstraints.AddItem(NULL); 439 return true; 440 } 441 442 443 bool 444 ActiveSetSolver::VariableRemoved(Variable* variable) 445 { 446 fVariableGEConstraints.RemoveItemAt(variable->Index()); 447 fVariableLEConstraints.RemoveItemAt(variable->Index()); 448 return true; 449 } 450 451 452 bool 453 ActiveSetSolver::VariableRangeChanged(Variable* variable) 454 { 455 double min = variable->Min(); 456 double max = variable->Max(); 457 int32 variableIndex = variable->Index(); 458 459 Constraint* constraintGE = fVariableGEConstraints.ItemAt(variableIndex); 460 Constraint* constraintLE = fVariableLEConstraints.ItemAt(variableIndex); 461 if (constraintGE == NULL && min > -20000) { 462 constraintGE = fLinearSpec->AddConstraint(1, variable, kGE, 0); 463 if (constraintGE == NULL) 464 return false; 465 fVariableGEConstraints.RemoveItemAt(variableIndex); 466 fVariableGEConstraints.AddItem(constraintGE, variableIndex); 467 } 468 if (constraintLE == NULL && max < 20000) { 469 constraintLE = fLinearSpec->AddConstraint(1, variable, kLE, 20000); 470 if (constraintLE == NULL) 471 return false; 472 fVariableLEConstraints.RemoveItemAt(variableIndex); 473 fVariableLEConstraints.AddItem(constraintLE, variableIndex); 474 } 475 476 if (constraintGE) 477 constraintGE->SetRightSide(min); 478 if (constraintLE) 479 constraintLE->SetRightSide(max); 480 return true; 481 } 482 483 484 bool 485 ActiveSetSolver::ConstraintAdded(Constraint* constraint) 486 { 487 return QPSolverInterface::ConstraintAdded(constraint); 488 } 489 490 491 bool 492 ActiveSetSolver::ConstraintRemoved(Constraint* constraint) 493 { 494 return QPSolverInterface::ConstraintRemoved(constraint); 495 } 496 497 498 bool 499 ActiveSetSolver::LeftSideChanged(Constraint* constraint) 500 { 501 return true; 502 } 503 504 505 bool 506 ActiveSetSolver::RightSideChanged(Constraint* constraint) 507 { 508 return true; 509 } 510 511 512 bool 513 ActiveSetSolver::OperatorChanged(Constraint* constraint) 514 { 515 return true; 516 } 517 518 519 bool 520 ActiveSetSolver::SaveModel(const char* fileName) 521 { 522 return false; 523 } 524 525 526 BSize 527 ActiveSetSolver::MinSize(Variable* width, Variable* height) 528 { 529 ConstraintList softConstraints; 530 _RemoveSoftConstraint(softConstraints); 531 532 Constraint* heightConstraint = fLinearSpec->AddConstraint(1, height, 533 kEQ, 0, 5, 5); 534 Constraint* widthConstraint = fLinearSpec->AddConstraint(1, width, 535 kEQ, 0, 5, 5); 536 ResultType result = Solve(); 537 fLinearSpec->RemoveConstraint(heightConstraint); 538 fLinearSpec->RemoveConstraint(widthConstraint); 539 540 _AddSoftConstraint(softConstraints); 541 542 if (result == kUnbounded) 543 return kMinSize; 544 if (result != kOptimal) 545 printf("Could not solve the layout specification (%d). ", result); 546 547 return BSize(width->Value(), height->Value()); 548 } 549 550 551 BSize 552 ActiveSetSolver::MaxSize(Variable* width, Variable* height) 553 { 554 ConstraintList softConstraints; 555 _RemoveSoftConstraint(softConstraints); 556 557 const double kHugeValue = 32000; 558 Constraint* heightConstraint = fLinearSpec->AddConstraint(1, height, 559 kEQ, kHugeValue, 5, 5); 560 Constraint* widthConstraint = fLinearSpec->AddConstraint(1, width, 561 kEQ, kHugeValue, 5, 5); 562 ResultType result = Solve(); 563 fLinearSpec->RemoveConstraint(heightConstraint); 564 fLinearSpec->RemoveConstraint(widthConstraint); 565 566 _AddSoftConstraint(softConstraints); 567 568 if (result == kUnbounded) 569 return kMaxSize; 570 if (result != kOptimal) 571 printf("Could not solve the layout specification (%d). ", result); 572 573 return BSize(width->Value(), height->Value()); 574 } 575 576 577 void 578 ActiveSetSolver::_RemoveSoftConstraint(ConstraintList& list) 579 { 580 ConstraintList allConstraints = fLinearSpec->Constraints(); 581 for (int i = 0; i < allConstraints.CountItems(); i++) { 582 Constraint* constraint = allConstraints.ItemAt(i); 583 if (!constraint->IsSoft()) 584 continue; 585 586 if (fLinearSpec->RemoveConstraint(constraint, false) == true) 587 list.AddItem(constraint); 588 } 589 } 590 591 592 void 593 ActiveSetSolver::_AddSoftConstraint(const ConstraintList& list) 594 { 595 for (int i = 0; i < list.CountItems(); i++) { 596 Constraint* constraint = list.ItemAt(i); 597 // at least don't leak it 598 if (fLinearSpec->AddConstraint(constraint) == false) 599 delete constraint; 600 } 601 } 602 603