 Timestamp:
 May 3, 2007 11:45:47 AM (13 years ago)
 Location:
 branches/devel/Clp/src
 Files:

 8 edited
Legend:
 Unmodified
 Added
 Removed

branches/devel/Clp/src/ClpConstraint.hpp
r997 r998 22 22 23 23 /** Fills gradient. If Linear then solution may be NULL, 24 also returns value of function24 also returns true value of function and offset so we can use x not deltaX in constraint 25 25 If refresh is false then uses last solution 26 26 Uses model for scaling … … 31 31 double * gradient, 32 32 double & functionValue , 33 double & offset, 33 34 bool refresh=true) const =0; 34 35 /// Resize constraint … … 42 43 */ 43 44 virtual int markNonlinear(char * which) const = 0; 45 /** Given a zeroed array sets possible nonzero coefficients to 1. 46 Returns number of nonzeros 47 */ 48 virtual int markNonzero(char * which) const = 0; 44 49 //@} 45 50 
branches/devel/Clp/src/ClpConstraintLinear.cpp
r997 r998 28 28 // Useful Constructor 29 29 // 30 ClpConstraintLinear::ClpConstraintLinear (int numberCoefficents ,30 ClpConstraintLinear::ClpConstraintLinear (int row, int numberCoefficents , 31 31 int numberColumns, 32 32 const int * column, const double * coefficient) … … 34 34 { 35 35 type_=0; 36 rowNumber_=row; 36 37 numberColumns_ = numberColumns; 37 38 numberCoefficients_=numberCoefficents; … … 92 93 const double * solution, 93 94 double * gradient, 94 double & functionValue, bool refresh) const 95 double & functionValue, 96 double & offset, 97 bool refresh) const 95 98 { 96 99 if (refresh!lastGradient_) { … … 121 124 } 122 125 functionValue = functionValue_; 126 offset=0.0; 123 127 memcpy(gradient,lastGradient_,numberColumns_*sizeof(double)); 124 128 return 0; … … 178 182 return 0; 179 183 } 184 /* Given a zeroed array sets possible nonzero coefficients to 1. 185 Returns number of nonzeros 186 */ 187 int 188 ClpConstraintLinear::markNonzero(char * which) const 189 { 190 for (int i=0;i<numberCoefficients_;i++) { 191 int iColumn = column_[i]; 192 which[iColumn]=1; 193 } 194 return numberCoefficients_; 195 } 
branches/devel/Clp/src/ClpConstraintLinear.hpp
r997 r998 21 21 22 22 /** Fills gradient. If Linear then solution may be NULL, 23 also returns value of function23 also returns true value of function and offset so we can use x not deltaX in constraint 24 24 If refresh is false then uses last solution 25 25 Uses model for scaling … … 30 30 double * gradient, 31 31 double & functionValue , 32 double & offset, 32 33 bool refresh=true) const ; 33 34 /// Resize constraint … … 41 42 */ 42 43 virtual int markNonlinear(char * which) const ; 44 /** Given a zeroed array sets possible nonzero coefficients to 1. 45 Returns number of nonzeros 46 */ 47 virtual int markNonzero(char * which) const; 43 48 //@} 44 49 … … 50 55 51 56 /// Constructor from constraint 52 ClpConstraintLinear(int numberCoefficients, int numberColumns,57 ClpConstraintLinear(int row, int numberCoefficients, int numberColumns, 53 58 const int * column, const double * element); 54 59 
branches/devel/Clp/src/ClpSimplex.cpp
r990 r998 4836 4836 { 4837 4837 return ((ClpSimplexNonlinear *) this)>primalSLP(numberPasses,deltaTolerance); 4838 } 4839 /* Solves problem with nonlinear constraints using SLP  may be used as crash 4840 for other algorithms when number of iterations small. 4841 Also exits if all problematical variables are changing 4842 less than deltaTolerance 4843 */ 4844 int 4845 ClpSimplex::nonlinearSLP(int numberConstraints, ClpConstraint ** constraints, 4846 int numberPasses,double deltaTolerance) 4847 { 4848 return ((ClpSimplexNonlinear *) this)>primalSLP(numberConstraints,constraints,numberPasses,deltaTolerance); 4838 4849 } 4839 4850 // Solves nonlinear using reduced gradient 
branches/devel/Clp/src/ClpSimplex.hpp
r990 r998 26 26 class CoinWarmStartBasis; 27 27 class ClpDisasterHandler; 28 class ClpConstraint; 28 29 29 30 /** This solves LPs using the simplex method … … 249 250 */ 250 251 int nonlinearSLP(int numberPasses,double deltaTolerance); 252 /** Solves problem with nonlinear constraints using SLP  may be used as crash 253 for other algorithms when number of iterations small. 254 Also exits if all problematical variables are changing 255 less than deltaTolerance 256 */ 257 int nonlinearSLP(int numberConstraints, ClpConstraint ** constraints, 258 int numberPasses,double deltaTolerance); 251 259 /** Solves using barrier (assumes you have good cholesky factor code). 252 260 Does crossover to simplex if asked*/ 
branches/devel/Clp/src/ClpSimplexNonlinear.cpp
r913 r998 11 11 #include "ClpNonLinearCost.hpp" 12 12 #include "ClpLinearObjective.hpp" 13 #include "ClpConstraint.hpp" 13 14 #include "ClpQuadraticObjective.hpp" 14 15 #include "CoinPackedMatrix.hpp" … … 3310 3311 return 0; 3311 3312 } 3313 /* Primal algorithm for nonlinear constraints 3314 Using a semitrust region approach as for pooling problem 3315 This is in because I have it lying around 3316 3317 */ 3318 int 3319 ClpSimplexNonlinear::primalSLP(int numberConstraints, ClpConstraint ** constraints, 3320 int numberPasses, double deltaTolerance) 3321 { 3322 if (!numberConstraints) { 3323 // no nonlinear part 3324 return ClpSimplex::primal(0); 3325 } 3326 // check all matrix for odd rows is empty 3327 int iConstraint; 3328 int numberBad=0; 3329 CoinPackedMatrix * columnCopy = matrix(); 3330 // Get a row copy in standard format 3331 CoinPackedMatrix copy; 3332 copy.reverseOrderedCopyOf(*columnCopy); 3333 // get matrix data pointers 3334 //const int * column = copy.getIndices(); 3335 //const CoinBigIndex * rowStart = copy.getVectorStarts(); 3336 const int * rowLength = copy.getVectorLengths(); 3337 //const double * elementByRow = copy.getElements(); 3338 int numberArtificials=0; 3339 for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) { 3340 ClpConstraint * constraint = constraints[iConstraint]; 3341 int iRow = constraint>rowNumber(); 3342 if (iRow>=0) { 3343 if (iRow>=numberRows_) 3344 numberBad++; 3345 else if (rowLength[iRow]) 3346 numberBad++; 3347 if (rowLower_[iRow]>1.0e20) 3348 numberArtificials++; 3349 if (rowUpper_[iRow]<1.0e20) 3350 numberArtificials++; 3351 } else { 3352 // objective 3353 assert (iRow==1); 3354 ClpLinearObjective * linearObj = (dynamic_cast< ClpLinearObjective*>(objective_)); 3355 assert (linearObj); 3356 // check empty 3357 const double * obj = objective(); 3358 int iColumn; 3359 for (iColumn=0;iColumn<numberColumns_;iColumn++) { 3360 if (obj[iColumn]) 3361 numberBad++; 3362 } 3363 } 3364 } 3365 if (numberBad) 3366 return numberBad; 3367 ClpSimplex newModel(*this); 3368 int numberColumns2 = numberColumns_; 3369 double penalty=1.0e9; 3370 if (numberArtificials) { 3371 numberColumns2 += numberArtificials; 3372 int * addStarts = new int [numberArtificials+1]; 3373 int * addRow = new int[numberArtificials]; 3374 double * addElement = new double[numberArtificials]; 3375 addStarts[0]=0; 3376 double * addCost = new double [numberArtificials]; 3377 numberArtificials=0; 3378 for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) { 3379 ClpConstraint * constraint = constraints[iConstraint]; 3380 int iRow = constraint>rowNumber(); 3381 if (iRow>=0) { 3382 if (rowLower_[iRow]>1.0e20) { 3383 addRow[numberArtificials]=iRow; 3384 addElement[numberArtificials]=1.0; 3385 addCost[numberArtificials]=penalty; 3386 numberArtificials++; 3387 addStarts[numberArtificials]=numberArtificials; 3388 } 3389 if (rowUpper_[iRow]<1.0e20) { 3390 addRow[numberArtificials]=iRow; 3391 addElement[numberArtificials]=1.0; 3392 addCost[numberArtificials]=penalty; 3393 numberArtificials++; 3394 addStarts[numberArtificials]=numberArtificials; 3395 } 3396 } 3397 } 3398 newModel.addColumns(numberArtificials,NULL,NULL,addCost, 3399 addStarts,addRow,addElement); 3400 delete [] addStarts; 3401 delete [] addRow; 3402 delete [] addElement; 3403 delete [] addCost; 3404 } 3405 // find nonlinear columns 3406 int * listNonLinearColumn = new int [numberColumns_]; 3407 char * mark = new char[numberColumns_]; 3408 memset(mark,0,numberColumns_); 3409 for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) { 3410 ClpConstraint * constraint = constraints[iConstraint]; 3411 constraint>markNonlinear(mark); 3412 } 3413 int iColumn; 3414 int numberNonLinearColumns=0; 3415 for (iColumn=0;iColumn<numberColumns_;iColumn++) { 3416 if (mark[iColumn]) 3417 listNonLinearColumn[numberNonLinearColumns++]=iColumn; 3418 } 3419 // build row copy of matrix with space for nonzeros 3420 // Get column copy 3421 columnCopy = matrix(); 3422 copy.reverseOrderedCopyOf(*columnCopy); 3423 // get matrix data pointers 3424 const int * column = copy.getIndices(); 3425 const CoinBigIndex * rowStart = copy.getVectorStarts(); 3426 rowLength = copy.getVectorLengths(); 3427 const double * elementByRow = copy.getElements(); 3428 CoinBigIndex numberExtra=0; 3429 for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) { 3430 ClpConstraint * constraint = constraints[iConstraint]; 3431 numberExtra+=constraint>markNonzero(mark); 3432 } 3433 numberExtra +=copy.getNumElements(); 3434 CoinBigIndex * newStarts = new CoinBigIndex [numberRows_+1]; 3435 int * newColumn = new int[numberExtra]; 3436 double * newElement = new double[numberExtra]; 3437 newStarts[0]=0; 3438 int * backRow = new int [numberRows_]; 3439 int iRow; 3440 for (iRow=0;iRow<numberRows_;iRow++) 3441 backRow[iRow]=1; 3442 for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) { 3443 ClpConstraint * constraint = constraints[iConstraint]; 3444 iRow = constraint>rowNumber(); 3445 backRow[iRow]=iConstraint; 3446 } 3447 numberExtra=0; 3448 for (iRow=0;iRow<numberRows_;iRow++) { 3449 if (backRow[iRow]<0) { 3450 // copy normal 3451 for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow]; 3452 j++) { 3453 newColumn[numberExtra]=column[j]; 3454 newElement[numberExtra++] = elementByRow[j]; 3455 } 3456 } else { 3457 // all possible 3458 memset(mark,0,numberColumns_); 3459 ClpConstraint * constraint = constraints[backRow[iRow]]; 3460 assert(iRow == constraint>rowNumber()); 3461 constraint>markNonzero(mark); 3462 for (int k=0;k<numberColumns_;k++) { 3463 if (mark[k]) { 3464 newColumn[numberExtra]=k; 3465 newElement[numberExtra++] = 1.0; 3466 } 3467 } 3468 } 3469 newStarts[iRow+1]=numberExtra; 3470 } 3471 delete [] backRow; 3472 CoinPackedMatrix saveMatrix(false,numberColumns2,numberRows_, 3473 numberExtra,newElement,newColumn,newStarts,NULL,0.0,0.0); 3474 delete [] newStarts; 3475 delete [] newColumn; 3476 delete [] newElement; 3477 delete [] mark; 3478 int numberRows = newModel.numberRows(); 3479 double * columnLower = newModel.columnLower(); 3480 double * columnUpper = newModel.columnUpper(); 3481 double * objective = newModel.objective(); 3482 double * rowLower = newModel.rowLower(); 3483 double * rowUpper = newModel.rowUpper(); 3484 double * solution = newModel.primalColumnSolution(); 3485 int jNon; 3486 int * last[3]; 3487 3488 double * trust = new double[numberNonLinearColumns]; 3489 double * trueLower = new double[numberNonLinearColumns]; 3490 double * trueUpper = new double[numberNonLinearColumns]; 3491 double objectiveOffset; 3492 double objectiveOffset2; 3493 getDblParam(ClpObjOffset,objectiveOffset); 3494 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { 3495 iColumn=listNonLinearColumn[jNon]; 3496 double upper = columnUpper[iColumn]; 3497 double lower = columnLower[iColumn]; 3498 if (upper>1.0e10) 3499 upper = max(1000.0,10*solution[iColumn]); 3500 //columnUpper[iColumn]=upper; 3501 trust[jNon]=0.05*(1.0+upperlower); 3502 trueLower[jNon]=lower; 3503 //trueUpper[jNon]=upper; 3504 trueUpper[jNon]=columnUpper[iColumn]; 3505 } 3506 bool tryFix=false; 3507 int iPass; 3508 double lastObjective=1.0e31; 3509 double * saveSolution = new double [numberColumns2+numberRows]; 3510 char * saveStatus = new char [numberColumns2+numberRows]; 3511 double targetDrop=1.0e31; 3512 //double objectiveOffset2 = objectiveOffset; 3513 // 1 bound up, 2 up, 1 bound down, 2 down, 0 no change 3514 for (iPass=0;iPass<3;iPass++) { 3515 last[iPass]=new int[numberNonLinearColumns]; 3516 for (jNon=0;jNon<numberNonLinearColumns;jNon++) 3517 last[iPass][jNon]=0; 3518 } 3519 int numberZeroPasses=0; 3520 bool zeroTargetDrop=false; 3521 double * gradient = new double [numberColumns_]; 3522 #define SMALL_FIX 0.0 3523 for (iPass=0;iPass<numberPasses;iPass++) { 3524 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { 3525 iColumn=listNonLinearColumn[jNon]; 3526 if (solution[iColumn]<trueLower[jNon]) 3527 solution[iColumn]=trueLower[jNon]; 3528 else if (solution[iColumn]>trueUpper[jNon]) 3529 solution[iColumn]=trueUpper[jNon]; 3530 if (iPass) { 3531 columnLower[iColumn]=max(solution[iColumn] 3532 trust[jNon], 3533 trueLower[jNon]); 3534 if (!trueLower[jNon]&&columnLower[iColumn]<SMALL_FIX) 3535 columnLower[iColumn]=SMALL_FIX; 3536 columnUpper[iColumn]=min(solution[iColumn] 3537 +trust[jNon], 3538 trueUpper[jNon]); 3539 if (!trueLower[jNon]) 3540 columnUpper[iColumn] = max(columnUpper[iColumn], 3541 columnLower[iColumn]+SMALL_FIX); 3542 if (!trueLower[jNon]&&tryFix&& 3543 columnLower[iColumn]==SMALL_FIX&& 3544 columnUpper[iColumn]<3.0*SMALL_FIX) { 3545 columnLower[iColumn]=0.0; 3546 solution[iColumn]=0.0; 3547 columnUpper[iColumn]=0.0; 3548 printf("fixing %d\n",iColumn); 3549 } 3550 } else { 3551 #if 1 3552 // fix nonlinear variables 3553 columnLower[iColumn]=solution[iColumn]; 3554 columnUpper[iColumn]=solution[iColumn]; 3555 #endif 3556 } 3557 } 3558 // redo matrix 3559 double offset; 3560 CoinPackedMatrix newMatrix(saveMatrix); 3561 // get matrix data pointers 3562 column = newMatrix.getIndices(); 3563 rowStart = newMatrix.getVectorStarts(); 3564 rowLength = newMatrix.getVectorLengths(); 3565 double * changeableElement = newMatrix.getMutableElements(); 3566 for (iConstraint=0;iConstraint<numberConstraints;iConstraint++) { 3567 ClpConstraint * constraint = constraints[iConstraint]; 3568 int iRow = constraint>rowNumber(); 3569 double functionValue; 3570 int numberErrors = constraint>gradient(&newModel,solution,gradient,functionValue,offset); 3571 assert (!numberErrors); 3572 if (iRow>=0) { 3573 if (numberRows<50) 3574 printf("For row %d current value is %g (activity %g) , dual is %g\n",iRow,offset, 3575 newModel.primalRowSolution()[iRow], 3576 newModel.dualRowSolution()[iRow]); 3577 for (CoinBigIndex j=rowStart[iRow];j<rowStart[iRow]+rowLength[iRow]; 3578 j++) { 3579 int iColumn = column[j]; 3580 changeableElement[j] = gradient[iColumn]; 3581 gradient[iColumn]=0.0; 3582 } 3583 for (int k=0;k<numberColumns_;k++) 3584 assert (!gradient[k]); 3585 if (rowLower_[iRow]>1.0e20) 3586 rowLower[iRow] = rowLower_[iRow]  offset; 3587 if (rowUpper_[iRow]<1.0e20) 3588 rowUpper[iRow] = rowUpper_[iRow]  offset; 3589 } else { 3590 memcpy(objective,gradient,numberColumns_*sizeof(double)); 3591 printf("objective offset %g\n",offset); 3592 objectiveOffset2 =objectiveOffsetoffset; 3593 newModel.setObjectiveOffset(objectiveOffset2); 3594 } 3595 } 3596 // Replace matrix 3597 // Get a column copy in standard format 3598 CoinPackedMatrix * columnCopy = new CoinPackedMatrix();; 3599 columnCopy>reverseOrderedCopyOf(newMatrix); 3600 newModel.replaceMatrix(columnCopy,true); 3601 double objValue=0.0; 3602 double infValue=0.0; 3603 { 3604 // get matrix data pointers 3605 const int * row = columnCopy>getIndices(); 3606 const CoinBigIndex * columnStart = columnCopy>getVectorStarts(); 3607 const int * columnLength = columnCopy>getVectorLengths(); 3608 const double * elementByColumn = columnCopy>getElements(); 3609 int iRow,iColumn; 3610 double * accumulate = new double [numberRows]; 3611 //const double * solution = newModel.primalColumnSolution(); 3612 double * partSolution = new double [numberColumns2]; 3613 memcpy(partSolution,solution,numberColumns_*sizeof(double)); 3614 memset(partSolution+numberColumns_,0, 3615 (numberColumns2numberColumns_)*sizeof(double)); 3616 memset(accumulate,0,numberRows*sizeof(double)); 3617 //newModel.matrix()>times(partSolution,accumulate); 3618 for (iColumn=0;iColumn<numberColumns_;iColumn++) { 3619 if (solution[iColumn]) { 3620 double value=solution[iColumn]; 3621 objValue += objective[iColumn]*solution[iColumn]; 3622 int j; 3623 for (j=columnStart[iColumn]; 3624 j<columnStart[iColumn]+columnLength[iColumn];j++) { 3625 int iRow=row[j]; 3626 accumulate[iRow] += value*elementByColumn[j]; 3627 } 3628 } 3629 } 3630 for (iRow=0;iRow<numberRows;iRow++) { 3631 if (accumulate[iRow]<rowLower[iRow]1.0e5) { 3632 double under = rowLower[iRow]accumulate[iRow]; 3633 infValue += under; 3634 #if 0 3635 printf("Row %d l=%g u=%g a=%g",iRow,rowLower[iRow], 3636 rowUpper[iRow],accumulate[iRow]); 3637 printf(" ** under by %g",under); 3638 printf("\n"); 3639 #endif 3640 } else if (accumulate[iRow]>rowUpper[iRow]+1.0e5) { 3641 double over = accumulate[iRow]rowUpper[iRow]; 3642 infValue += over; 3643 #if 0 3644 printf("Row %d l=%g u=%g a=%g",iRow,rowLower[iRow], 3645 rowUpper[iRow],accumulate[iRow]); 3646 printf(" ** over by %g",over); 3647 printf("\n"); 3648 #endif 3649 } 3650 } 3651 if (infValue) 3652 printf("Sum infeasibilities %g ",infValue); 3653 if (infValue<0.1) 3654 infValue=0.0; 3655 infValue *= penalty; 3656 if (infValue) 3657 printf("Infeasible obj %g ",infValue); 3658 delete [] accumulate; 3659 delete [] partSolution; 3660 } 3661 if (objectiveOffset2) 3662 printf("offset2 %g\n",objectiveOffset2); 3663 objValue = objValue + infValue objectiveOffset2; 3664 printf("True objective %g\n",objValue); 3665 if (iPass) { 3666 double drop = lastObjectiveobjValue; 3667 std::cout<<"True drop was "<<drop<<std::endl; 3668 if (drop<0.05*fabs(objValue)1.0e4) { 3669 // pretty bad  go back and halve 3670 memcpy(solution,saveSolution,numberColumns2*sizeof(double)); 3671 memcpy(newModel.primalRowSolution(),saveSolution+numberColumns2, 3672 numberRows*sizeof(double)); 3673 memcpy(newModel.statusArray(),saveStatus, 3674 numberColumns2+numberRows); 3675 for (jNon=0;jNon<numberNonLinearColumns;jNon++) 3676 if (trust[jNon]>0.1) 3677 trust[jNon] *= 0.5; 3678 else 3679 trust[jNon] *= 0.9; 3680 3681 printf("** Halving trust\n"); 3682 objValue=lastObjective; 3683 continue; 3684 } else if ((iPass%25)==1) { 3685 for (jNon=0;jNon<numberNonLinearColumns;jNon++) 3686 trust[jNon] *= 2.0; 3687 printf("** Doubling trust\n"); 3688 } 3689 int * temp=last[2]; 3690 last[2]=last[1]; 3691 last[1]=last[0]; 3692 last[0]=temp; 3693 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { 3694 iColumn=listNonLinearColumn[jNon]; 3695 double change = solution[iColumn]saveSolution[iColumn]; 3696 if (change<1.0e5) { 3697 if (fabs(change+trust[jNon])<1.0e5) 3698 temp[jNon]=1; 3699 else 3700 temp[jNon]=2; 3701 } else if(change>1.0e5) { 3702 if (fabs(changetrust[jNon])<1.0e5) 3703 temp[jNon]=1; 3704 else 3705 temp[jNon]=2; 3706 } else { 3707 temp[jNon]=0; 3708 } 3709 } 3710 double maxDelta=0.0; 3711 double smallestTrust=1.0e31; 3712 double smallestNonLinearGap=1.0e31; 3713 double smallestGap=1.0e31; 3714 bool increasing=false; 3715 for (iColumn=0;iColumn<numberColumns_;iColumn++) { 3716 double gap = columnUpper[iColumn]columnLower[iColumn]; 3717 assert (gap>=0.0); 3718 if (gap) 3719 smallestGap = min(smallestGap,gap); 3720 } 3721 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { 3722 iColumn=listNonLinearColumn[jNon]; 3723 double gap = columnUpper[iColumn]columnLower[iColumn]; 3724 assert (gap>=0.0); 3725 if (gap) { 3726 smallestNonLinearGap = min(smallestNonLinearGap,gap); 3727 if (gap<1.0e7&&iPass==1) { 3728 printf("Small gap %d %d %g %g %g\n", 3729 jNon,iColumn,columnLower[iColumn],columnUpper[iColumn], 3730 gap); 3731 //trueUpper[jNon]=trueLower[jNon]; 3732 //columnUpper[iColumn]=columnLower[iColumn]; 3733 } 3734 } 3735 maxDelta = max(maxDelta, 3736 fabs(solution[iColumn]saveSolution[iColumn])); 3737 if (last[0][jNon]*last[1][jNon]<0) { 3738 // halve 3739 if (trust[jNon]>1.0) 3740 trust[jNon] *= 0.5; 3741 else 3742 trust[jNon] *= 0.7; 3743 } else { 3744 // ? only increase if +=1 ? 3745 if (last[0][jNon]==last[1][jNon]&& 3746 last[0][jNon]==last[2][jNon]&& 3747 last[0][jNon]) { 3748 trust[jNon] *= 1.8; 3749 increasing=true; 3750 } 3751 } 3752 smallestTrust = min(smallestTrust,trust[jNon]); 3753 } 3754 std::cout<<"largest delta is "<<maxDelta 3755 <<", smallest trust is "<<smallestTrust 3756 <<", smallest gap is "<<smallestGap 3757 <<", smallest nonlinear gap is "<<smallestNonLinearGap 3758 <<std::endl; 3759 if (maxDelta<deltaTolerance&&!increasing&&iPass>100) { 3760 numberZeroPasses++; 3761 if (numberZeroPasses==3) { 3762 if (tryFix) { 3763 std::cout<<"Exiting"<<std::endl; 3764 break; 3765 } else { 3766 tryFix=true; 3767 for (jNon=0;jNon<numberNonLinearColumns;jNon++) 3768 trust[jNon] = max(trust[jNon],1.0e1); 3769 numberZeroPasses=0; 3770 } 3771 } 3772 } else { 3773 numberZeroPasses=0; 3774 } 3775 } 3776 memcpy(saveSolution,solution,numberColumns2*sizeof(double)); 3777 memcpy(saveSolution+numberColumns2,newModel.primalRowSolution(), 3778 numberRows*sizeof(double)); 3779 memcpy(saveStatus,newModel.statusArray(), 3780 numberColumns2+numberRows); 3781 3782 targetDrop=0.0; 3783 if (iPass) { 3784 // get reduced costs 3785 const double * pi = newModel.dualRowSolution(); 3786 newModel.matrix()>transposeTimes(pi, 3787 newModel.dualColumnSolution()); 3788 double * r = newModel.dualColumnSolution(); 3789 for (iColumn=0;iColumn<numberColumns_;iColumn++) 3790 r[iColumn] = objective[iColumn]r[iColumn]; 3791 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { 3792 iColumn=listNonLinearColumn[jNon]; 3793 double dj = r[iColumn]; 3794 if (dj<1.0e6) { 3795 double drop = dj*(columnUpper[iColumn]solution[iColumn]); 3796 //double upper = min(trueUpper[jNon],solution[iColumn]+0.1); 3797 //double drop2 = dj*(uppersolution[iColumn]); 3798 #if 0 3799 if (drop>1.0e8drop2>100.0*drop(drop>1.0e2&&iPass>100)) 3800 printf("Big drop %d %g %g %g %g T %g\n", 3801 iColumn,columnLower[iColumn],solution[iColumn], 3802 columnUpper[iColumn],dj,trueUpper[jNon]); 3803 #endif 3804 targetDrop += drop; 3805 if (dj<1.0e1&&trust[jNon]<1.0e3 3806 &&trueUpper[jNon]solution[iColumn]>1.0e2) { 3807 trust[jNon] *= 1.5; 3808 //printf("Increasing trust on %d to %g\n", 3809 // iColumn,trust[jNon]); 3810 } 3811 } else if (dj>1.0e6) { 3812 double drop = dj*(columnLower[iColumn]solution[iColumn]); 3813 //double lower = max(trueLower[jNon],solution[iColumn]0.1); 3814 //double drop2 = dj*(lowersolution[iColumn]); 3815 #if 0 3816 if (drop>1.0e8drop2>100.0*drop(drop>1.0e2)) 3817 printf("Big drop %d %g %g %g %g T %g\n", 3818 iColumn,columnLower[iColumn],solution[iColumn], 3819 columnUpper[iColumn],dj,trueLower[jNon]); 3820 #endif 3821 targetDrop += drop; 3822 if (dj>1.0e1&&trust[jNon]<1.0e3 3823 &&solution[iColumn]trueLower[jNon]>1.0e2) { 3824 trust[jNon] *= 1.5; 3825 printf("Increasing trust on %d to %g\n", 3826 iColumn,trust[jNon]); 3827 } 3828 } 3829 } 3830 } 3831 std::cout<<"Pass  "<<iPass 3832 <<", target drop is "<<targetDrop 3833 <<std::endl; 3834 if (iPass>1&&targetDrop<1.0e5&&zeroTargetDrop) 3835 break; 3836 if (iPass>1&&targetDrop<1.0e5) 3837 zeroTargetDrop = true; 3838 else 3839 zeroTargetDrop = false; 3840 //if (iPass==5) 3841 //newModel.setLogLevel(63); 3842 lastObjective = objValue; 3843 // take out when ClpPackedMatrix changed 3844 //newModel.scaling(false); 3845 #if 0 3846 CoinMpsIO writer; 3847 writer.setMpsData(*newModel.matrix(), COIN_DBL_MAX, 3848 newModel.getColLower(), newModel.getColUpper(), 3849 newModel.getObjCoefficients(), 3850 (const char*) 0 /*integrality*/, 3851 newModel.getRowLower(), newModel.getRowUpper(), 3852 NULL,NULL); 3853 writer.writeMps("xx.mps"); 3854 #endif 3855 if (iPass<3) { 3856 newModel.primal(1); 3857 abort(); 3858 } else { 3859 abort(); 3860 } 3861 if (newModel.status()==1) { 3862 // Infeasible! 3863 newModel.allSlackBasis(); 3864 newModel.primal(); 3865 assert(!newModel.status()); 3866 } 3867 double sumInfeas=0.0; 3868 int numberInfeas=0; 3869 for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn++) { 3870 if (solution[iColumn]>1.0e8) { 3871 numberInfeas++; 3872 sumInfeas += solution[iColumn]; 3873 } 3874 } 3875 printf("%d infeasibilities  summing to %g\n", 3876 numberInfeas,sumInfeas); 3877 //model.writeMps("xx"); 3878 } 3879 delete [] saveSolution; 3880 delete [] saveStatus; 3881 for (iPass=0;iPass<3;iPass++) 3882 delete [] last[iPass]; 3883 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { 3884 iColumn=listNonLinearColumn[jNon]; 3885 columnLower[iColumn]=trueLower[jNon]; 3886 columnUpper[iColumn]=trueUpper[jNon]; 3887 } 3888 delete [] trust; 3889 delete [] trueUpper; 3890 delete [] trueLower; 3891 // Simplest way to get true row activity ? 3892 double * rowActivity = newModel.primalRowSolution(); 3893 for (iRow=0;iRow<numberRows;iRow++) { 3894 double difference; 3895 if (fabs(rowLower_[iRow])<fabs(rowUpper_[iRow])) 3896 difference = rowLower_[iRow]rowLower[iRow]; 3897 else 3898 difference = rowUpper_[iRow]rowUpper[iRow]; 3899 rowLower[iRow]=rowLower_[iRow]; 3900 rowUpper[iRow]=rowUpper_[iRow]; 3901 if (difference) { 3902 if (numberRows<50) 3903 printf("For row %d activity changes from %g to %g\n", 3904 iRow,rowActivity[iRow],rowActivity[iRow]+difference); 3905 rowActivity[iRow]+= difference; 3906 } 3907 } 3908 delete [] saveSolution; 3909 delete [] saveStatus; 3910 for (iPass=0;iPass<3;iPass++) 3911 delete [] last[iPass]; 3912 for (jNon=0;jNon<numberNonLinearColumns;jNon++) { 3913 iColumn=listNonLinearColumn[jNon]; 3914 columnLower[iColumn]=trueLower[jNon]; 3915 columnUpper[iColumn]=trueUpper[jNon]; 3916 } 3917 delete [] trust; 3918 delete [] trueUpper; 3919 delete [] trueLower; 3920 delete [] listNonLinearColumn; 3921 delete [] gradient; 3922 printf("solution still in newModel  do objective etc!\n"); 3923 objectiveValue_=newModel.objectiveValue(); 3924 numberIterations_ =newModel.numberIterations(); 3925 problemStatus_ =newModel.problemStatus(); 3926 secondaryStatus_ =newModel.secondaryStatus(); 3927 memcpy(columnActivity_,newModel.primalColumnSolution(),numberColumns_*sizeof(double)); 3928 // should do status region 3929 CoinZeroN(rowActivity_,numberRows_); 3930 matrix_>times(1.0,columnActivity_,rowActivity_); 3931 abort(); 3932 return 0; 3933 } 
branches/devel/Clp/src/ClpSimplexNonlinear.hpp
r754 r998 13 13 class ClpNonlinearInfo; 14 14 class ClpQuadraticObjective; 15 class ClpConstraint; 15 16 16 17 #include "ClpSimplexPrimal.hpp" … … 42 43 */ 43 44 int primalSLP(int numberPasses, double deltaTolerance); 45 /** Primal algorithm for nonlinear constraints 46 Using a semitrust region approach as for pooling problem 47 This is in because I have it lying around 48 49 */ 50 int primalSLP(int numberConstraints, ClpConstraint ** constraints, 51 int numberPasses, double deltaTolerance); 44 52 45 53 /** Creates direction vector. note longArray is long enough 
branches/devel/Clp/src/Makefile.in
r996 r998 69 69 libClp_la_LIBADD = 70 70 am_libClp_la_OBJECTS = ClpCholeskyBase.lo ClpCholeskyDense.lo \ 71 ClpCholeskyUfl.lo ClpCholeskyWssmp.lo ClpCholeskyWssmpKKT.lo \ 72 ClpConstraint.lo ClpConstraintLinear.lo Clp_C_Interface.lo \ 73 ClpDualRowDantzig.lo ClpDualRowPivot.lo ClpDualRowSteepest.lo \ 74 ClpDummyMatrix.lo ClpDynamicExampleMatrix.lo \ 75 ClpDynamicMatrix.lo ClpEventHandler.lo ClpFactorization.lo \ 76 ClpGubDynamicMatrix.lo ClpGubMatrix.lo ClpHelperFunctions.lo \ 77 ClpInterior.lo ClpLinearObjective.lo ClpMatrixBase.lo \ 78 ClpMessage.lo ClpModel.lo ClpNetworkBasis.lo \ 79 ClpNetworkMatrix.lo ClpNonLinearCost.lo ClpObjective.lo \ 80 ClpPackedMatrix.lo ClpPlusMinusOneMatrix.lo \ 81 ClpPredictorCorrector.lo ClpPresolve.lo \ 82 ClpPrimalColumnDantzig.lo ClpPrimalColumnPivot.lo \ 83 ClpPrimalColumnSteepest.lo ClpQuadraticObjective.lo \ 84 ClpSimplex.lo ClpSimplexDual.lo ClpSimplexNonlinear.lo \ 85 ClpSimplexOther.lo ClpSimplexPrimal.lo ClpSolve.lo Idiot.lo \ 86 IdiSolve.lo 71 ClpCholeskyUfl.lo ClpConstraint.lo ClpConstraintLinear.lo \ 72 Clp_C_Interface.lo ClpDualRowDantzig.lo ClpDualRowPivot.lo \ 73 ClpDualRowSteepest.lo ClpDummyMatrix.lo \ 74 ClpDynamicExampleMatrix.lo ClpDynamicMatrix.lo \ 75 ClpEventHandler.lo ClpFactorization.lo ClpGubDynamicMatrix.lo \ 76 ClpGubMatrix.lo ClpHelperFunctions.lo ClpInterior.lo \ 77 ClpLinearObjective.lo ClpMatrixBase.lo ClpMessage.lo \ 78 ClpModel.lo ClpNetworkBasis.lo ClpNetworkMatrix.lo \ 79 ClpNonLinearCost.lo ClpObjective.lo ClpPackedMatrix.lo \ 80 ClpPlusMinusOneMatrix.lo ClpPredictorCorrector.lo \ 81 ClpPresolve.lo ClpPrimalColumnDantzig.lo \ 82 ClpPrimalColumnPivot.lo ClpPrimalColumnSteepest.lo \ 83 ClpQuadraticObjective.lo ClpSimplex.lo ClpSimplexDual.lo \ 84 ClpSimplexNonlinear.lo ClpSimplexOther.lo ClpSimplexPrimal.lo \ 85 ClpSolve.lo Idiot.lo IdiSolve.lo 87 86 libClp_la_OBJECTS = $(am_libClp_la_OBJECTS) 88 87 binPROGRAMS_INSTALL = $(INSTALL_PROGRAM) … … 276 275 ClpCholeskyDense.cpp ClpCholeskyDense.hpp \ 277 276 ClpCholeskyUfl.cpp ClpCholeskyUfl.hpp \ 278 ClpCholeskyWssmp.cpp ClpCholeskyWssmp.hpp \279 ClpCholeskyWssmpKKT.cpp ClpCholeskyWssmpKKT.hpp \280 277 ClpConstraint.cpp ClpConstraint.hpp \ 281 278 ClpConstraintLinear.cpp ClpConstraintLinear.hpp \ … … 514 511 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpCholeskyDense.Plo@am__quote@ 515 512 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpCholeskyUfl.Plo@am__quote@ 516 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpCholeskyWssmp.Plo@am__quote@517 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpCholeskyWssmpKKT.Plo@am__quote@518 513 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpConstraint.Plo@am__quote@ 519 514 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/ClpConstraintLinear.Plo@am__quote@
Note: See TracChangeset
for help on using the changeset viewer.