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::GaussJordan() 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 GaussJordan(i); 162 } 163 return true; 164 } 165 166 167 void 168 EquationSystem::GaussJordan(int32 i) 169 { 170 double value = fMatrix[fRowIndices[i]][fColumnIndices[i]]; 171 for (int j = 0; j < fColumns; j++) 172 fMatrix[fRowIndices[i]][fColumnIndices[j]] /= value; 173 fB[i] /= value; 174 175 for (int r = 0; r < fRows; r++) { 176 if (r == i) 177 continue; 178 double q = -fMatrix[fRowIndices[r]][fColumnIndices[i]]; 179 // don't need to do nothing, since matrix is typically sparse this 180 // should save some work 181 if (fuzzy_equals(q, 0)) 182 continue; 183 for (int c = 0; c < fColumns; c++) 184 fMatrix[fRowIndices[r]][fColumnIndices[c]] 185 += fMatrix[fRowIndices[i]][fColumnIndices[c]] * q; 186 187 fB[r] += fB[i] * q; 188 } 189 } 190 191 192 void 193 EquationSystem::RemoveLinearlyDependentRows() 194 { 195 double oldB[fRows]; 196 for (int r = 0; r < fRows; r++) 197 oldB[r] = fB[r]; 198 199 double** temp = allocate_matrix(fRows, fColumns); 200 bool independentRows[fRows]; 201 202 // copy to temp 203 copy_matrix(fMatrix, temp, fRows, fColumns); 204 int nIndependent = compute_dependencies(temp, fRows, fColumns, 205 independentRows); 206 if (nIndependent == fRows) 207 return; 208 209 // remove the rows 210 for (int i = 0; i < fRows; i++) { 211 if (!independentRows[i]) { 212 int lastDepRow = -1; 213 for (int d = fRows - 1; d > i; d--) { 214 if (independentRows[d]) { 215 lastDepRow = d; 216 break; 217 } 218 } 219 if (lastDepRow < 0) 220 break; 221 SwapRow(i, lastDepRow); 222 fRows--; 223 } 224 } 225 fRows = nIndependent; 226 227 free_matrix(temp); 228 } 229 230 231 void 232 EquationSystem::RemoveUnusedVariables() 233 { 234 for (int c = 0; c < fColumns; c++) { 235 bool used = false; 236 for (int r = 0; r < fRows; r++) { 237 if (!fuzzy_equals(fMatrix[r][fColumnIndices[c]], 0)) { 238 used = true; 239 break; 240 } 241 } 242 if (used) 243 continue; 244 245 //MoveColumnRight(c, fColumns - 1); 246 SwapColumn(c, fColumns - 1); 247 fColumns--; 248 c--; 249 } 250 } 251 252 253 void 254 EquationSystem::MoveColumnRight(int32 i, int32 target) 255 { 256 int32 index = fColumnIndices[i]; 257 for (int c = i; c < target; c++) 258 fColumnIndices[c] = fColumnIndices[c + 1]; 259 fColumnIndices[target] = index; 260 } 261 262 263 void 264 EquationSystem::Print() 265 { 266 for (int m = 0; m < fRows; m++) { 267 for (int n = 0; n < fColumns; n++) 268 printf("%.1f ", fMatrix[fRowIndices[m]][fColumnIndices[n]]); 269 printf("= %.1f\n", fB[m]); 270 } 271 } 272 273 274 ActiveSetSolver::ActiveSetSolver(LinearSpec* linearSpec) 275 : 276 QPSolverInterface(linearSpec), 277 278 fVariables(linearSpec->UsedVariables()), 279 fConstraints(linearSpec->Constraints()) 280 { 281 282 } 283 284 285 ActiveSetSolver::~ActiveSetSolver() 286 { 287 288 } 289 290 291 /* Using algorithm found in: 292 Solving Inequalities and Proving Farkas's Lemma Made Easy 293 David Avis and Bohdan Kaluzny 294 The American Mathematical Monthly 295 Vol. 111, No. 2 (Feb., 2004), pp. 152-157 */ 296 bool 297 solve(EquationSystem& system) 298 { 299 // basic solve 300 if (!system.GaussJordan()) 301 return false; 302 303 bool done = false; 304 while (!done) { 305 double smallestB = HUGE_VALF; 306 int smallestBRow = -1; 307 for (int row = 0; row < system.Rows(); row++) { 308 if (system.B(row) > 0 || fuzzy_equals(system.B(row), 0)) 309 continue; 310 311 double bValue = fabs(system.B(row)); 312 if (bValue < smallestB) { 313 smallestB = bValue; 314 smallestBRow = row; 315 } 316 } 317 if (smallestBRow == -1) { 318 done = true; 319 break; 320 } 321 322 int negValueCol = -1; 323 for (int col = system.Rows(); col < system.Columns(); col++) { 324 double value = system.A(smallestBRow, col); 325 if (value > 0 || fuzzy_equals(value, 0)) 326 continue; 327 negValueCol = col; 328 break; 329 } 330 if (negValueCol == -1) { 331 printf("can't solve\n"); 332 return false; 333 } 334 335 system.SwapColumn(smallestBRow, negValueCol); 336 337 // eliminate 338 system.GaussJordan(smallestBRow); 339 } 340 341 return true; 342 } 343 344 345 ResultType 346 ActiveSetSolver::Solve() 347 { 348 int32 nConstraints = fConstraints.CountItems(); 349 int32 nVariables = fVariables.CountItems(); 350 351 if (nVariables > nConstraints) { 352 printf("More variables then constraints! vars: %i, constraints: %i\n", 353 (int)nVariables, (int)nConstraints); 354 return kInfeasible; 355 } 356 357 /* First find an initial solution and then optimize it using the active set 358 method. */ 359 EquationSystem system(nConstraints, nVariables + nConstraints); 360 361 int32 slackIndex = nVariables; 362 // setup constraint matrix and add slack variables if necessary 363 int32 rowIndex = 0; 364 for (int32 c = 0; c < nConstraints; c++) { 365 Constraint* constraint = fConstraints.ItemAt(c); 366 if (constraint->IsSoft()) 367 continue; 368 SummandList* leftSide = constraint->LeftSide(); 369 system.B(rowIndex) = constraint->RightSide(); 370 for (int32 sIndex = 0; sIndex < leftSide->CountItems(); sIndex++ ) { 371 Summand* summand = leftSide->ItemAt(sIndex); 372 double coefficient = summand->Coeff(); 373 system.A(rowIndex, summand->VariableIndex()) = coefficient; 374 } 375 if (constraint->Op() == kLE) { 376 system.A(rowIndex, slackIndex) = 1; 377 slackIndex++; 378 } else if (constraint->Op() == kGE) { 379 system.A(rowIndex, slackIndex) = -1; 380 slackIndex++; 381 } 382 rowIndex++; 383 } 384 385 system.SetRows(rowIndex); 386 387 system.RemoveLinearlyDependentRows(); 388 system.RemoveUnusedVariables(); 389 390 if (!solve(system)) 391 return kInfeasible; 392 393 double results[nVariables + nConstraints]; 394 system.Results(results, nVariables + nConstraints); 395 printf("base system solved\n"); 396 397 LayoutOptimizer optimizer(fConstraints, nVariables); 398 optimizer.Solve(results); 399 400 // back to the variables 401 for (int32 i = 0; i < nVariables; i++) 402 fVariables.ItemAt(i)->SetValue(results[i]); 403 404 for (int32 i = 0; i < nVariables; i++) 405 TRACE("var %f\n", results[i]); 406 407 return kOptimal; 408 } 409 410 411 bool 412 ActiveSetSolver::VariableAdded(Variable* variable) 413 { 414 // TODO: error checks 415 fVariableGEConstraints.AddItem(NULL); 416 fVariableLEConstraints.AddItem(NULL); 417 return true; 418 } 419 420 421 bool 422 ActiveSetSolver::VariableRemoved(Variable* variable) 423 { 424 fVariableGEConstraints.RemoveItemAt(variable->Index()); 425 fVariableLEConstraints.RemoveItemAt(variable->Index()); 426 return true; 427 } 428 429 430 bool 431 ActiveSetSolver::VariableRangeChanged(Variable* variable) 432 { 433 double min = variable->Min(); 434 double max = variable->Max(); 435 int32 variableIndex = variable->Index(); 436 437 Constraint* constraintGE = fVariableGEConstraints.ItemAt(variableIndex); 438 Constraint* constraintLE = fVariableLEConstraints.ItemAt(variableIndex); 439 if (constraintGE == NULL && min > -20000) { 440 constraintGE = fLinearSpec->AddConstraint(1, variable, kGE, 0); 441 if (constraintGE == NULL) 442 return false; 443 fVariableGEConstraints.RemoveItemAt(variableIndex); 444 fVariableGEConstraints.AddItem(constraintGE, variableIndex); 445 } 446 if (constraintLE == NULL && max < 20000) { 447 constraintLE = fLinearSpec->AddConstraint(1, variable, kLE, 20000); 448 if (constraintLE == NULL) 449 return false; 450 fVariableLEConstraints.RemoveItemAt(variableIndex); 451 fVariableLEConstraints.AddItem(constraintLE, variableIndex); 452 } 453 454 if (constraintGE) 455 constraintGE->SetRightSide(min); 456 if (constraintLE) 457 constraintLE->SetRightSide(max); 458 return true; 459 } 460 461 462 bool 463 ActiveSetSolver::ConstraintAdded(Constraint* constraint) 464 { 465 return QPSolverInterface::ConstraintAdded(constraint); 466 } 467 468 469 bool 470 ActiveSetSolver::ConstraintRemoved(Constraint* constraint) 471 { 472 return QPSolverInterface::ConstraintRemoved(constraint); 473 } 474 475 476 bool 477 ActiveSetSolver::LeftSideChanged(Constraint* constraint) 478 { 479 return true; 480 } 481 482 483 bool 484 ActiveSetSolver::RightSideChanged(Constraint* constraint) 485 { 486 return true; 487 } 488 489 490 bool 491 ActiveSetSolver::OperatorChanged(Constraint* constraint) 492 { 493 return true; 494 } 495 496 497 bool 498 ActiveSetSolver::SaveModel(const char* fileName) 499 { 500 return false; 501 } 502 503 504 BSize 505 ActiveSetSolver::MinSize(Variable* width, Variable* height) 506 { 507 ConstraintList softConstraints; 508 _RemoveSoftConstraint(softConstraints); 509 510 Constraint* heightConstraint = fLinearSpec->AddConstraint(1, height, 511 kEQ, 0, 5, 5); 512 Constraint* widthConstraint = fLinearSpec->AddConstraint(1, width, 513 kEQ, 0, 5, 5); 514 ResultType result = Solve(); 515 fLinearSpec->RemoveConstraint(heightConstraint); 516 fLinearSpec->RemoveConstraint(widthConstraint); 517 518 _AddSoftConstraint(softConstraints); 519 520 if (result == kUnbounded) 521 return kMinSize; 522 if (result != kOptimal) 523 printf("Could not solve the layout specification (%d). ", result); 524 525 return BSize(width->Value(), height->Value()); 526 } 527 528 529 BSize 530 ActiveSetSolver::MaxSize(Variable* width, Variable* height) 531 { 532 ConstraintList softConstraints; 533 _RemoveSoftConstraint(softConstraints); 534 535 const double kHugeValue = 32000; 536 Constraint* heightConstraint = fLinearSpec->AddConstraint(1, height, 537 kEQ, kHugeValue, 5, 5); 538 Constraint* widthConstraint = fLinearSpec->AddConstraint(1, width, 539 kEQ, kHugeValue, 5, 5); 540 ResultType result = Solve(); 541 fLinearSpec->RemoveConstraint(heightConstraint); 542 fLinearSpec->RemoveConstraint(widthConstraint); 543 544 _AddSoftConstraint(softConstraints); 545 546 if (result == kUnbounded) 547 return kMaxSize; 548 if (result != kOptimal) 549 printf("Could not solve the layout specification (%d). ", result); 550 551 return BSize(width->Value(), height->Value()); 552 } 553 554 555 void 556 ActiveSetSolver::_RemoveSoftConstraint(ConstraintList& list) 557 { 558 ConstraintList allConstraints = fLinearSpec->Constraints(); 559 for (int i = 0; i < allConstraints.CountItems(); i++) { 560 Constraint* constraint = allConstraints.ItemAt(i); 561 if (!constraint->IsSoft()) 562 continue; 563 564 if (fLinearSpec->RemoveConstraint(constraint, false) == true) 565 list.AddItem(constraint); 566 } 567 } 568 569 570 void 571 ActiveSetSolver::_AddSoftConstraint(const ConstraintList& list) 572 { 573 for (int i = 0; i < list.CountItems(); i++) { 574 Constraint* constraint = list.ItemAt(i); 575 // at least don't leak it 576 if (fLinearSpec->AddConstraint(constraint) == false) 577 delete constraint; 578 } 579 } 580 581