1 | /* $Id: ClpSimplexNonlinear.cpp 2385 2019-01-06 19:43:06Z unxusr $ */ |
---|
2 | // Copyright (C) 2004, International Business Machines |
---|
3 | // Corporation and others. All Rights Reserved. |
---|
4 | // This code is licensed under the terms of the Eclipse Public License (EPL). |
---|
5 | |
---|
6 | #include "CoinPragma.hpp" |
---|
7 | |
---|
8 | #include <math.h> |
---|
9 | #include "CoinHelperFunctions.hpp" |
---|
10 | #include "ClpHelperFunctions.hpp" |
---|
11 | #include "ClpSimplexNonlinear.hpp" |
---|
12 | #include "ClpSimplexOther.hpp" |
---|
13 | #include "ClpSimplexDual.hpp" |
---|
14 | #include "ClpFactorization.hpp" |
---|
15 | #include "ClpNonLinearCost.hpp" |
---|
16 | #include "ClpLinearObjective.hpp" |
---|
17 | #include "ClpConstraint.hpp" |
---|
18 | #include "ClpPresolve.hpp" |
---|
19 | #include "ClpQuadraticObjective.hpp" |
---|
20 | #include "CoinPackedMatrix.hpp" |
---|
21 | #include "CoinIndexedVector.hpp" |
---|
22 | #include "ClpPrimalColumnPivot.hpp" |
---|
23 | #include "ClpMessage.hpp" |
---|
24 | #include "ClpEventHandler.hpp" |
---|
25 | #include <cfloat> |
---|
26 | #include <cassert> |
---|
27 | #include <string> |
---|
28 | #include <stdio.h> |
---|
29 | #include <iostream> |
---|
30 | #ifndef NDEBUG |
---|
31 | #define CLP_DEBUG 1 |
---|
32 | #endif |
---|
33 | // primal |
---|
34 | int ClpSimplexNonlinear::primal() |
---|
35 | { |
---|
36 | |
---|
37 | int ifValuesPass = 1; |
---|
38 | algorithm_ = +3; |
---|
39 | |
---|
40 | // save data |
---|
41 | ClpDataSave data = saveData(); |
---|
42 | matrix_->refresh(this); // make sure matrix okay |
---|
43 | |
---|
44 | // Save objective |
---|
45 | ClpObjective *saveObjective = NULL; |
---|
46 | if (objective_->type() > 1) { |
---|
47 | // expand to full if quadratic |
---|
48 | #ifndef NO_RTTI |
---|
49 | ClpQuadraticObjective *quadraticObj = (dynamic_cast< ClpQuadraticObjective * >(objective_)); |
---|
50 | #else |
---|
51 | ClpQuadraticObjective *quadraticObj = NULL; |
---|
52 | if (objective_->type() == 2) |
---|
53 | quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_)); |
---|
54 | #endif |
---|
55 | // for moment only if no scaling |
---|
56 | // May be faster if switched off - but can't see why |
---|
57 | if (!quadraticObj->fullMatrix() && (!rowScale_ && !scalingFlag_) && objectiveScale_ == 1.0) { |
---|
58 | saveObjective = objective_; |
---|
59 | objective_ = new ClpQuadraticObjective(*quadraticObj, 1); |
---|
60 | } |
---|
61 | } |
---|
62 | double bestObjectiveWhenFlagged = COIN_DBL_MAX; |
---|
63 | int pivotMode = 15; |
---|
64 | //pivotMode=20; |
---|
65 | |
---|
66 | // initialize - maybe values pass and algorithm_ is +1 |
---|
67 | // true does something (? perturbs) |
---|
68 | if (!startup(true)) { |
---|
69 | |
---|
70 | // Set average theta |
---|
71 | nonLinearCost_->setAverageTheta(1.0e3); |
---|
72 | int lastCleaned = 0; // last time objective or bounds cleaned up |
---|
73 | |
---|
74 | // Say no pivot has occurred (for steepest edge and updates) |
---|
75 | pivotRow_ = -2; |
---|
76 | |
---|
77 | // This says whether to restore things etc |
---|
78 | int factorType = 0; |
---|
79 | // Start check for cycles |
---|
80 | progress_.startCheck(); |
---|
81 | /* |
---|
82 | Status of problem: |
---|
83 | 0 - optimal |
---|
84 | 1 - infeasible |
---|
85 | 2 - unbounded |
---|
86 | -1 - iterating |
---|
87 | -2 - factorization wanted |
---|
88 | -3 - redo checking without factorization |
---|
89 | -4 - looks infeasible |
---|
90 | -5 - looks unbounded |
---|
91 | */ |
---|
92 | while (problemStatus_ < 0) { |
---|
93 | int iRow, iColumn; |
---|
94 | // clear |
---|
95 | for (iRow = 0; iRow < 4; iRow++) { |
---|
96 | rowArray_[iRow]->clear(); |
---|
97 | } |
---|
98 | |
---|
99 | for (iColumn = 0; iColumn < SHORT_REGION; iColumn++) { |
---|
100 | columnArray_[iColumn]->clear(); |
---|
101 | } |
---|
102 | |
---|
103 | // give matrix (and model costs and bounds a chance to be |
---|
104 | // refreshed (normally null) |
---|
105 | matrix_->refresh(this); |
---|
106 | // If getting nowhere - why not give it a kick |
---|
107 | // If we have done no iterations - special |
---|
108 | if (lastGoodIteration_ == numberIterations_ && factorType) |
---|
109 | factorType = 3; |
---|
110 | |
---|
111 | // may factorize, checks if problem finished |
---|
112 | if (objective_->type() > 1 && lastFlaggedIteration_ >= 0 && numberIterations_ > lastFlaggedIteration_ + 507) { |
---|
113 | unflag(); |
---|
114 | lastFlaggedIteration_ = numberIterations_; |
---|
115 | if (pivotMode >= 10) { |
---|
116 | pivotMode--; |
---|
117 | #ifdef CLP_DEBUG |
---|
118 | if (handler_->logLevel() & 32) |
---|
119 | printf("pivot mode now %d\n", pivotMode); |
---|
120 | #endif |
---|
121 | if (pivotMode == 9) |
---|
122 | pivotMode = 0; // switch off fast attempt |
---|
123 | } |
---|
124 | } |
---|
125 | statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, |
---|
126 | bestObjectiveWhenFlagged); |
---|
127 | |
---|
128 | // Say good factorization |
---|
129 | factorType = 1; |
---|
130 | |
---|
131 | // Say no pivot has occurred (for steepest edge and updates) |
---|
132 | pivotRow_ = -2; |
---|
133 | |
---|
134 | // exit if victory declared |
---|
135 | if (problemStatus_ >= 0) |
---|
136 | break; |
---|
137 | |
---|
138 | // test for maximum iterations |
---|
139 | if (hitMaximumIterations() || (ifValuesPass == 2 && firstFree_ < 0)) { |
---|
140 | problemStatus_ = 3; |
---|
141 | break; |
---|
142 | } |
---|
143 | |
---|
144 | if (firstFree_ < 0) { |
---|
145 | if (ifValuesPass) { |
---|
146 | // end of values pass |
---|
147 | ifValuesPass = 0; |
---|
148 | int status = eventHandler_->event(ClpEventHandler::endOfValuesPass); |
---|
149 | if (status >= 0) { |
---|
150 | problemStatus_ = 5; |
---|
151 | secondaryStatus_ = ClpEventHandler::endOfValuesPass; |
---|
152 | break; |
---|
153 | } |
---|
154 | } |
---|
155 | } |
---|
156 | // Check event |
---|
157 | { |
---|
158 | int status = eventHandler_->event(ClpEventHandler::endOfFactorization); |
---|
159 | if (status >= 0) { |
---|
160 | problemStatus_ = 5; |
---|
161 | secondaryStatus_ = ClpEventHandler::endOfFactorization; |
---|
162 | break; |
---|
163 | } |
---|
164 | } |
---|
165 | // Iterate |
---|
166 | whileIterating(pivotMode); |
---|
167 | } |
---|
168 | } |
---|
169 | // if infeasible get real values |
---|
170 | if (problemStatus_ == 1) { |
---|
171 | infeasibilityCost_ = 0.0; |
---|
172 | createRim(1 + 4); |
---|
173 | delete nonLinearCost_; |
---|
174 | nonLinearCost_ = new ClpNonLinearCost(this); |
---|
175 | nonLinearCost_->checkInfeasibilities(0.0); |
---|
176 | sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities(); |
---|
177 | numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities(); |
---|
178 | // and get good feasible duals |
---|
179 | computeDuals(NULL); |
---|
180 | } |
---|
181 | // correct objective value |
---|
182 | if (numberColumns_) |
---|
183 | objectiveValue_ = nonLinearCost_->feasibleCost() + objective_->nonlinearOffset(); |
---|
184 | objectiveValue_ /= (objectiveScale_ * rhsScale_); |
---|
185 | // clean up |
---|
186 | unflag(); |
---|
187 | finish(); |
---|
188 | restoreData(data); |
---|
189 | // restore objective if full |
---|
190 | if (saveObjective) { |
---|
191 | delete objective_; |
---|
192 | objective_ = saveObjective; |
---|
193 | } |
---|
194 | return problemStatus_; |
---|
195 | } |
---|
196 | /* Refactorizes if necessary |
---|
197 | Checks if finished. Updates status. |
---|
198 | lastCleaned refers to iteration at which some objective/feasibility |
---|
199 | cleaning too place. |
---|
200 | |
---|
201 | type - 0 initial so set up save arrays etc |
---|
202 | - 1 normal -if good update save |
---|
203 | - 2 restoring from saved |
---|
204 | */ |
---|
205 | void ClpSimplexNonlinear::statusOfProblemInPrimal(int &lastCleaned, int type, |
---|
206 | ClpSimplexProgress *progress, |
---|
207 | bool doFactorization, |
---|
208 | double &bestObjectiveWhenFlagged) |
---|
209 | { |
---|
210 | int dummy; // for use in generalExpanded |
---|
211 | if (type == 2) { |
---|
212 | // trouble - restore solution |
---|
213 | CoinMemcpyN(saveStatus_, (numberColumns_ + numberRows_), status_); |
---|
214 | CoinMemcpyN(savedSolution_ + numberColumns_, numberRows_, rowActivityWork_); |
---|
215 | CoinMemcpyN(savedSolution_, numberColumns_, columnActivityWork_); |
---|
216 | // restore extra stuff |
---|
217 | matrix_->generalExpanded(this, 6, dummy); |
---|
218 | forceFactorization_ = 1; // a bit drastic but .. |
---|
219 | pivotRow_ = -1; // say no weights update |
---|
220 | changeMade_++; // say change made |
---|
221 | } |
---|
222 | int saveThreshold = factorization_->sparseThreshold(); |
---|
223 | int tentativeStatus = problemStatus_; |
---|
224 | int numberThrownOut = 1; // to loop round on bad factorization in values pass |
---|
225 | while (numberThrownOut) { |
---|
226 | if (problemStatus_ > -3 || problemStatus_ == -4) { |
---|
227 | // factorize |
---|
228 | // later on we will need to recover from singularities |
---|
229 | // also we could skip if first time |
---|
230 | // do weights |
---|
231 | // This may save pivotRow_ for use |
---|
232 | if (doFactorization) |
---|
233 | primalColumnPivot_->saveWeights(this, 1); |
---|
234 | |
---|
235 | if (type && doFactorization) { |
---|
236 | // is factorization okay? |
---|
237 | int factorStatus = internalFactorize(1); |
---|
238 | if (factorStatus) { |
---|
239 | if (type != 1 || largestPrimalError_ > 1.0e3 |
---|
240 | || largestDualError_ > 1.0e3) { |
---|
241 | // was ||largestDualError_>1.0e3||objective_->type()>1) { |
---|
242 | // switch off dense |
---|
243 | int saveDense = factorization_->denseThreshold(); |
---|
244 | factorization_->setDenseThreshold(0); |
---|
245 | // make sure will do safe factorization |
---|
246 | pivotVariable_[0] = -1; |
---|
247 | internalFactorize(2); |
---|
248 | factorization_->setDenseThreshold(saveDense); |
---|
249 | // Go to safe |
---|
250 | factorization_->pivotTolerance(0.99); |
---|
251 | // restore extra stuff |
---|
252 | matrix_->generalExpanded(this, 6, dummy); |
---|
253 | } else { |
---|
254 | // no - restore previous basis |
---|
255 | CoinMemcpyN(saveStatus_, (numberColumns_ + numberRows_), status_); |
---|
256 | CoinMemcpyN(savedSolution_ + numberColumns_, numberRows_, rowActivityWork_); |
---|
257 | CoinMemcpyN(savedSolution_, numberColumns_, columnActivityWork_); |
---|
258 | // restore extra stuff |
---|
259 | matrix_->generalExpanded(this, 6, dummy); |
---|
260 | matrix_->generalExpanded(this, 5, dummy); |
---|
261 | forceFactorization_ = 1; // a bit drastic but .. |
---|
262 | type = 2; |
---|
263 | // Go to safe |
---|
264 | factorization_->pivotTolerance(0.99); |
---|
265 | if (internalFactorize(1) != 0) |
---|
266 | largestPrimalError_ = 1.0e4; // force other type |
---|
267 | } |
---|
268 | // flag incoming |
---|
269 | if (sequenceIn_ >= 0 && getStatus(sequenceIn_) != basic) { |
---|
270 | setFlagged(sequenceIn_); |
---|
271 | saveStatus_[sequenceIn_] = status_[sequenceIn_]; |
---|
272 | } |
---|
273 | changeMade_++; // say change made |
---|
274 | } |
---|
275 | } |
---|
276 | if (problemStatus_ != -4) |
---|
277 | problemStatus_ = -3; |
---|
278 | } |
---|
279 | // at this stage status is -3 or -5 if looks unbounded |
---|
280 | // get primal and dual solutions |
---|
281 | // put back original costs and then check |
---|
282 | createRim(4); |
---|
283 | // May need to do more if column generation |
---|
284 | dummy = 4; |
---|
285 | matrix_->generalExpanded(this, 9, dummy); |
---|
286 | numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= 0)); |
---|
287 | if (numberThrownOut) { |
---|
288 | problemStatus_ = tentativeStatus; |
---|
289 | doFactorization = true; |
---|
290 | } |
---|
291 | } |
---|
292 | // Double check reduced costs if no action |
---|
293 | if (progress->lastIterationNumber(0) == numberIterations_) { |
---|
294 | if (primalColumnPivot_->looksOptimal()) { |
---|
295 | numberDualInfeasibilities_ = 0; |
---|
296 | sumDualInfeasibilities_ = 0.0; |
---|
297 | } |
---|
298 | } |
---|
299 | // Check if looping |
---|
300 | int loop; |
---|
301 | if (type != 2) |
---|
302 | loop = progress->looping(); |
---|
303 | else |
---|
304 | loop = -1; |
---|
305 | if (loop >= 0) { |
---|
306 | if (!problemStatus_) { |
---|
307 | // declaring victory |
---|
308 | numberPrimalInfeasibilities_ = 0; |
---|
309 | sumPrimalInfeasibilities_ = 0.0; |
---|
310 | } else { |
---|
311 | problemStatus_ = loop; //exit if in loop |
---|
312 | problemStatus_ = 10; // instead - try other algorithm |
---|
313 | } |
---|
314 | problemStatus_ = 10; // instead - try other algorithm |
---|
315 | return; |
---|
316 | } else if (loop < -1) { |
---|
317 | // Is it time for drastic measures |
---|
318 | if (nonLinearCost_->numberInfeasibilities() && progress->badTimes() > 5 && progress->oddState() < 10 && progress->oddState() >= 0) { |
---|
319 | progress->newOddState(); |
---|
320 | nonLinearCost_->zapCosts(); |
---|
321 | } |
---|
322 | // something may have changed |
---|
323 | gutsOfSolution(NULL, NULL, true); |
---|
324 | } |
---|
325 | // If progress then reset costs |
---|
326 | if (loop == -1 && !nonLinearCost_->numberInfeasibilities() && progress->oddState() < 0) { |
---|
327 | createRim(4, false); // costs back |
---|
328 | delete nonLinearCost_; |
---|
329 | nonLinearCost_ = new ClpNonLinearCost(this); |
---|
330 | progress->endOddState(); |
---|
331 | gutsOfSolution(NULL, NULL, true); |
---|
332 | } |
---|
333 | // Flag to say whether to go to dual to clean up |
---|
334 | //bool goToDual = false; |
---|
335 | // really for free variables in |
---|
336 | //if((progressFlag_&2)!=0) |
---|
337 | //problemStatus_=-1;; |
---|
338 | progressFlag_ = 0; //reset progress flag |
---|
339 | |
---|
340 | handler_->message(CLP_SIMPLEX_STATUS, messages_) |
---|
341 | << numberIterations_ << nonLinearCost_->feasibleReportCost(); |
---|
342 | handler_->printing(nonLinearCost_->numberInfeasibilities() > 0) |
---|
343 | << nonLinearCost_->sumInfeasibilities() << nonLinearCost_->numberInfeasibilities(); |
---|
344 | handler_->printing(sumDualInfeasibilities_ > 0.0) |
---|
345 | << sumDualInfeasibilities_ << numberDualInfeasibilities_; |
---|
346 | handler_->printing(numberDualInfeasibilitiesWithoutFree_ |
---|
347 | < numberDualInfeasibilities_) |
---|
348 | << numberDualInfeasibilitiesWithoutFree_; |
---|
349 | handler_->message() << CoinMessageEol; |
---|
350 | if (!primalFeasible()) { |
---|
351 | nonLinearCost_->checkInfeasibilities(primalTolerance_); |
---|
352 | gutsOfSolution(NULL, NULL, true); |
---|
353 | nonLinearCost_->checkInfeasibilities(primalTolerance_); |
---|
354 | } |
---|
355 | double trueInfeasibility = nonLinearCost_->sumInfeasibilities(); |
---|
356 | if (trueInfeasibility > 1.0) { |
---|
357 | // If infeasibility going up may change weights |
---|
358 | double testValue = trueInfeasibility - 1.0e-4 * (10.0 + trueInfeasibility); |
---|
359 | if (progress->lastInfeasibility() < testValue) { |
---|
360 | if (infeasibilityCost_ < 1.0e14) { |
---|
361 | infeasibilityCost_ *= 1.5; |
---|
362 | if (handler_->logLevel() == 63) |
---|
363 | printf("increasing weight to %g\n", infeasibilityCost_); |
---|
364 | gutsOfSolution(NULL, NULL, true); |
---|
365 | } |
---|
366 | } |
---|
367 | } |
---|
368 | // we may wish to say it is optimal even if infeasible |
---|
369 | bool alwaysOptimal = (specialOptions_ & 1) != 0; |
---|
370 | // give code benefit of doubt |
---|
371 | if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) { |
---|
372 | // say optimal (with these bounds etc) |
---|
373 | numberDualInfeasibilities_ = 0; |
---|
374 | sumDualInfeasibilities_ = 0.0; |
---|
375 | numberPrimalInfeasibilities_ = 0; |
---|
376 | sumPrimalInfeasibilities_ = 0.0; |
---|
377 | } |
---|
378 | // had ||(type==3&&problemStatus_!=-5) -- ??? why ???? |
---|
379 | if (dualFeasible() || problemStatus_ == -4) { |
---|
380 | // see if extra helps |
---|
381 | if (nonLinearCost_->numberInfeasibilities() && (nonLinearCost_->sumInfeasibilities() > 1.0e-3 || sumOfRelaxedPrimalInfeasibilities_) |
---|
382 | && !alwaysOptimal) { |
---|
383 | //may need infeasiblity cost changed |
---|
384 | // we can see if we can construct a ray |
---|
385 | // make up a new objective |
---|
386 | double saveWeight = infeasibilityCost_; |
---|
387 | // save nonlinear cost as we are going to switch off costs |
---|
388 | ClpNonLinearCost *nonLinear = nonLinearCost_; |
---|
389 | // do twice to make sure Primal solution has settled |
---|
390 | // put non-basics to bounds in case tolerance moved |
---|
391 | // put back original costs |
---|
392 | createRim(4); |
---|
393 | nonLinearCost_->checkInfeasibilities(0.0); |
---|
394 | gutsOfSolution(NULL, NULL, true); |
---|
395 | |
---|
396 | infeasibilityCost_ = 1.0e100; |
---|
397 | // put back original costs |
---|
398 | createRim(4); |
---|
399 | nonLinearCost_->checkInfeasibilities(primalTolerance_); |
---|
400 | // may have fixed infeasibilities - double check |
---|
401 | if (nonLinearCost_->numberInfeasibilities() == 0) { |
---|
402 | // carry on |
---|
403 | problemStatus_ = -1; |
---|
404 | infeasibilityCost_ = saveWeight; |
---|
405 | nonLinearCost_->checkInfeasibilities(primalTolerance_); |
---|
406 | } else { |
---|
407 | nonLinearCost_ = NULL; |
---|
408 | // scale |
---|
409 | int i; |
---|
410 | for (i = 0; i < numberRows_ + numberColumns_; i++) |
---|
411 | cost_[i] *= 1.0e-95; |
---|
412 | gutsOfSolution(NULL, NULL, false); |
---|
413 | nonLinearCost_ = nonLinear; |
---|
414 | infeasibilityCost_ = saveWeight; |
---|
415 | if ((infeasibilityCost_ >= 1.0e18 || numberDualInfeasibilities_ == 0) && perturbation_ == 101) { |
---|
416 | /* goToDual = unPerturb(); // stop any further perturbation |
---|
417 | if (nonLinearCost_->sumInfeasibilities() > 1.0e-1) |
---|
418 | goToDual = false; |
---|
419 | */ |
---|
420 | unPerturb(); |
---|
421 | nonLinearCost_->checkInfeasibilities(primalTolerance_); |
---|
422 | numberDualInfeasibilities_ = 1; // carry on |
---|
423 | problemStatus_ = -1; |
---|
424 | } |
---|
425 | if (infeasibilityCost_ >= 1.0e20 || numberDualInfeasibilities_ == 0) { |
---|
426 | // we are infeasible - use as ray |
---|
427 | delete[] ray_; |
---|
428 | ray_ = new double[numberRows_]; |
---|
429 | CoinMemcpyN(dual_, numberRows_, ray_); |
---|
430 | // and get feasible duals |
---|
431 | infeasibilityCost_ = 0.0; |
---|
432 | createRim(4); |
---|
433 | nonLinearCost_->checkInfeasibilities(primalTolerance_); |
---|
434 | gutsOfSolution(NULL, NULL, true); |
---|
435 | // so will exit |
---|
436 | infeasibilityCost_ = 1.0e30; |
---|
437 | // reset infeasibilities |
---|
438 | sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities(); |
---|
439 | ; |
---|
440 | numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities(); |
---|
441 | } |
---|
442 | if (infeasibilityCost_ < 1.0e20) { |
---|
443 | infeasibilityCost_ *= 5.0; |
---|
444 | changeMade_++; // say change made |
---|
445 | unflag(); |
---|
446 | handler_->message(CLP_PRIMAL_WEIGHT, messages_) |
---|
447 | << infeasibilityCost_ |
---|
448 | << CoinMessageEol; |
---|
449 | // put back original costs and then check |
---|
450 | createRim(4); |
---|
451 | nonLinearCost_->checkInfeasibilities(0.0); |
---|
452 | gutsOfSolution(NULL, NULL, true); |
---|
453 | problemStatus_ = -1; //continue |
---|
454 | //goToDual = false; |
---|
455 | } else { |
---|
456 | // say infeasible |
---|
457 | problemStatus_ = 1; |
---|
458 | } |
---|
459 | } |
---|
460 | } else { |
---|
461 | // may be optimal |
---|
462 | if (perturbation_ == 101) { |
---|
463 | /* goToDual =*/unPerturb(); // stop any further perturbation |
---|
464 | lastCleaned = -1; // carry on |
---|
465 | } |
---|
466 | bool unflagged = unflag() != 0; |
---|
467 | if (lastCleaned != numberIterations_ || unflagged) { |
---|
468 | handler_->message(CLP_PRIMAL_OPTIMAL, messages_) |
---|
469 | << primalTolerance_ |
---|
470 | << CoinMessageEol; |
---|
471 | double current = nonLinearCost_->feasibleReportCost(); |
---|
472 | if (numberTimesOptimal_ < 4) { |
---|
473 | if (bestObjectiveWhenFlagged <= current) { |
---|
474 | numberTimesOptimal_++; |
---|
475 | #ifdef CLP_DEBUG |
---|
476 | if (handler_->logLevel() & 32) |
---|
477 | printf("%d times optimal, current %g best %g\n", numberTimesOptimal_, |
---|
478 | current, bestObjectiveWhenFlagged); |
---|
479 | #endif |
---|
480 | } else { |
---|
481 | bestObjectiveWhenFlagged = current; |
---|
482 | } |
---|
483 | changeMade_++; // say change made |
---|
484 | if (numberTimesOptimal_ == 1) { |
---|
485 | // better to have small tolerance even if slower |
---|
486 | factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15)); |
---|
487 | } |
---|
488 | lastCleaned = numberIterations_; |
---|
489 | if (primalTolerance_ != dblParam_[ClpPrimalTolerance]) |
---|
490 | handler_->message(CLP_PRIMAL_ORIGINAL, messages_) |
---|
491 | << CoinMessageEol; |
---|
492 | double oldTolerance = primalTolerance_; |
---|
493 | primalTolerance_ = dblParam_[ClpPrimalTolerance]; |
---|
494 | // put back original costs and then check |
---|
495 | createRim(4); |
---|
496 | nonLinearCost_->checkInfeasibilities(oldTolerance); |
---|
497 | gutsOfSolution(NULL, NULL, true); |
---|
498 | if (sumOfRelaxedDualInfeasibilities_ == 0.0 && sumOfRelaxedPrimalInfeasibilities_ == 0.0) { |
---|
499 | // say optimal (with these bounds etc) |
---|
500 | numberDualInfeasibilities_ = 0; |
---|
501 | sumDualInfeasibilities_ = 0.0; |
---|
502 | numberPrimalInfeasibilities_ = 0; |
---|
503 | sumPrimalInfeasibilities_ = 0.0; |
---|
504 | } |
---|
505 | if (dualFeasible() && !nonLinearCost_->numberInfeasibilities() && lastCleaned >= 0) |
---|
506 | problemStatus_ = 0; |
---|
507 | else |
---|
508 | problemStatus_ = -1; |
---|
509 | } else { |
---|
510 | problemStatus_ = 0; // optimal |
---|
511 | if (lastCleaned < numberIterations_) { |
---|
512 | handler_->message(CLP_SIMPLEX_GIVINGUP, messages_) |
---|
513 | << CoinMessageEol; |
---|
514 | } |
---|
515 | } |
---|
516 | } else { |
---|
517 | problemStatus_ = 0; // optimal |
---|
518 | } |
---|
519 | } |
---|
520 | } else { |
---|
521 | // see if looks unbounded |
---|
522 | if (problemStatus_ == -5) { |
---|
523 | if (nonLinearCost_->numberInfeasibilities()) { |
---|
524 | if (infeasibilityCost_ > 1.0e18 && perturbation_ == 101) { |
---|
525 | // back off weight |
---|
526 | infeasibilityCost_ = 1.0e13; |
---|
527 | unPerturb(); // stop any further perturbation |
---|
528 | } |
---|
529 | //we need infeasiblity cost changed |
---|
530 | if (infeasibilityCost_ < 1.0e20) { |
---|
531 | infeasibilityCost_ *= 5.0; |
---|
532 | changeMade_++; // say change made |
---|
533 | unflag(); |
---|
534 | handler_->message(CLP_PRIMAL_WEIGHT, messages_) |
---|
535 | << infeasibilityCost_ |
---|
536 | << CoinMessageEol; |
---|
537 | // put back original costs and then check |
---|
538 | createRim(4); |
---|
539 | gutsOfSolution(NULL, NULL, true); |
---|
540 | problemStatus_ = -1; //continue |
---|
541 | } else { |
---|
542 | // say unbounded |
---|
543 | problemStatus_ = 2; |
---|
544 | } |
---|
545 | } else { |
---|
546 | // say unbounded |
---|
547 | problemStatus_ = 2; |
---|
548 | } |
---|
549 | } else { |
---|
550 | if (type == 3 && problemStatus_ != -5) |
---|
551 | unflag(); // odd |
---|
552 | // carry on |
---|
553 | problemStatus_ = -1; |
---|
554 | } |
---|
555 | } |
---|
556 | // save extra stuff |
---|
557 | matrix_->generalExpanded(this, 5, dummy); |
---|
558 | if (type == 0 || type == 1) { |
---|
559 | if (type != 1 || !saveStatus_) { |
---|
560 | // create save arrays |
---|
561 | delete[] saveStatus_; |
---|
562 | delete[] savedSolution_; |
---|
563 | saveStatus_ = new unsigned char[numberRows_ + numberColumns_]; |
---|
564 | savedSolution_ = new double[numberRows_ + numberColumns_]; |
---|
565 | } |
---|
566 | // save arrays |
---|
567 | CoinMemcpyN(status_, (numberColumns_ + numberRows_), saveStatus_); |
---|
568 | CoinMemcpyN(rowActivityWork_, numberRows_, savedSolution_ + numberColumns_); |
---|
569 | CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_); |
---|
570 | } |
---|
571 | if (doFactorization) { |
---|
572 | // restore weights (if saved) - also recompute infeasibility list |
---|
573 | if (tentativeStatus > -3) |
---|
574 | primalColumnPivot_->saveWeights(this, (type < 2) ? 2 : 4); |
---|
575 | else |
---|
576 | primalColumnPivot_->saveWeights(this, 3); |
---|
577 | if (saveThreshold) { |
---|
578 | // use default at present |
---|
579 | factorization_->sparseThreshold(0); |
---|
580 | factorization_->goSparse(); |
---|
581 | } |
---|
582 | } |
---|
583 | if (problemStatus_ < 0 && !changeMade_) { |
---|
584 | problemStatus_ = 4; // unknown |
---|
585 | } |
---|
586 | lastGoodIteration_ = numberIterations_; |
---|
587 | //if (goToDual) |
---|
588 | //problemStatus_=10; // try dual |
---|
589 | // Allow matrices to be sorted etc |
---|
590 | int fake = -999; // signal sort |
---|
591 | matrix_->correctSequence(this, fake, fake); |
---|
592 | } |
---|
593 | /* |
---|
594 | Reasons to come out: |
---|
595 | -1 iterations etc |
---|
596 | -2 inaccuracy |
---|
597 | -3 slight inaccuracy (and done iterations) |
---|
598 | -4 end of values pass and done iterations |
---|
599 | +0 looks optimal (might be infeasible - but we will investigate) |
---|
600 | +2 looks unbounded |
---|
601 | +3 max iterations |
---|
602 | */ |
---|
603 | int ClpSimplexNonlinear::whileIterating(int &pivotMode) |
---|
604 | { |
---|
605 | // Say if values pass |
---|
606 | //int ifValuesPass=(firstFree_>=0) ? 1 : 0; |
---|
607 | int ifValuesPass = 1; |
---|
608 | int returnCode = -1; |
---|
609 | // status stays at -1 while iterating, >=0 finished, -2 to invert |
---|
610 | // status -3 to go to top without an invert |
---|
611 | int numberInterior = 0; |
---|
612 | int nextUnflag = 10; |
---|
613 | int nextUnflagIteration = numberIterations_ + 10; |
---|
614 | // get two arrays |
---|
615 | double *array1 = new double[2 * (numberRows_ + numberColumns_)]; |
---|
616 | double solutionError = -1.0; |
---|
617 | while (problemStatus_ == -1) { |
---|
618 | int result; |
---|
619 | rowArray_[1]->clear(); |
---|
620 | //#define CLP_DEBUG |
---|
621 | #if CLP_DEBUG > 1 |
---|
622 | rowArray_[0]->checkClear(); |
---|
623 | rowArray_[1]->checkClear(); |
---|
624 | rowArray_[2]->checkClear(); |
---|
625 | rowArray_[3]->checkClear(); |
---|
626 | columnArray_[0]->checkClear(); |
---|
627 | #endif |
---|
628 | if (numberInterior >= 5) { |
---|
629 | // this can go when we have better minimization |
---|
630 | if (pivotMode < 10) |
---|
631 | pivotMode = 1; |
---|
632 | unflag(); |
---|
633 | #ifdef CLP_DEBUG |
---|
634 | if (handler_->logLevel() & 32) |
---|
635 | printf("interior unflag\n"); |
---|
636 | #endif |
---|
637 | numberInterior = 0; |
---|
638 | nextUnflag = 10; |
---|
639 | nextUnflagIteration = numberIterations_ + 10; |
---|
640 | } else { |
---|
641 | if (numberInterior > nextUnflag && numberIterations_ > nextUnflagIteration) { |
---|
642 | nextUnflagIteration = numberIterations_ + 10; |
---|
643 | nextUnflag += 10; |
---|
644 | unflag(); |
---|
645 | #ifdef CLP_DEBUG |
---|
646 | if (handler_->logLevel() & 32) |
---|
647 | printf("unflagging as interior\n"); |
---|
648 | #endif |
---|
649 | } |
---|
650 | } |
---|
651 | pivotRow_ = -1; |
---|
652 | result = pivotColumn(rowArray_[3], rowArray_[0], |
---|
653 | columnArray_[0], rowArray_[1], pivotMode, solutionError, |
---|
654 | array1); |
---|
655 | if (result) { |
---|
656 | if (result == 2 && sequenceIn_ < 0) { |
---|
657 | // does not look good |
---|
658 | double currentObj; |
---|
659 | double thetaObj; |
---|
660 | double predictedObj; |
---|
661 | objective_->stepLength(this, solution_, solution_, 0.0, |
---|
662 | currentObj, thetaObj, predictedObj); |
---|
663 | if (currentObj == predictedObj) { |
---|
664 | #ifdef CLP_INVESTIGATE |
---|
665 | printf("looks bad - no change in obj %g\n", currentObj); |
---|
666 | #endif |
---|
667 | if (factorization_->pivots()) |
---|
668 | result = 3; |
---|
669 | else |
---|
670 | problemStatus_ = 0; |
---|
671 | } |
---|
672 | } |
---|
673 | if (result == 3) |
---|
674 | break; // null vector not accurate |
---|
675 | #ifdef CLP_DEBUG |
---|
676 | if (handler_->logLevel() & 32) { |
---|
677 | double currentObj; |
---|
678 | double thetaObj; |
---|
679 | double predictedObj; |
---|
680 | objective_->stepLength(this, solution_, solution_, 0.0, |
---|
681 | currentObj, thetaObj, predictedObj); |
---|
682 | printf("obj %g after interior move\n", currentObj); |
---|
683 | } |
---|
684 | #endif |
---|
685 | // just move and try again |
---|
686 | if (pivotMode < 10) { |
---|
687 | pivotMode = result - 1; |
---|
688 | numberInterior++; |
---|
689 | } |
---|
690 | continue; |
---|
691 | } else { |
---|
692 | if (pivotMode < 10) { |
---|
693 | if (theta_ > 0.001) |
---|
694 | pivotMode = 0; |
---|
695 | else if (pivotMode == 2) |
---|
696 | pivotMode = 1; |
---|
697 | } |
---|
698 | numberInterior = 0; |
---|
699 | nextUnflag = 10; |
---|
700 | nextUnflagIteration = numberIterations_ + 10; |
---|
701 | } |
---|
702 | sequenceOut_ = -1; |
---|
703 | rowArray_[1]->clear(); |
---|
704 | if (sequenceIn_ >= 0) { |
---|
705 | // we found a pivot column |
---|
706 | assert(!flagged(sequenceIn_)); |
---|
707 | #ifdef CLP_DEBUG |
---|
708 | if ((handler_->logLevel() & 32)) { |
---|
709 | char x = isColumn(sequenceIn_) ? 'C' : 'R'; |
---|
710 | std::cout << "pivot column " << x << sequenceWithin(sequenceIn_) << std::endl; |
---|
711 | } |
---|
712 | #endif |
---|
713 | // do second half of iteration |
---|
714 | if (pivotRow_ < 0 && theta_ < 1.0e-8) { |
---|
715 | assert(sequenceIn_ >= 0); |
---|
716 | returnCode = pivotResult(ifValuesPass); |
---|
717 | } else { |
---|
718 | // specialized code |
---|
719 | returnCode = pivotNonlinearResult(); |
---|
720 | //printf("odd pivrow %d\n",sequenceOut_); |
---|
721 | if (sequenceOut_ >= 0 && theta_ < 1.0e-5) { |
---|
722 | if (getStatus(sequenceOut_) != isFixed) { |
---|
723 | if (getStatus(sequenceOut_) == atUpperBound) |
---|
724 | solution_[sequenceOut_] = upper_[sequenceOut_]; |
---|
725 | else if (getStatus(sequenceOut_) == atLowerBound) |
---|
726 | solution_[sequenceOut_] = lower_[sequenceOut_]; |
---|
727 | setFlagged(sequenceOut_); |
---|
728 | } |
---|
729 | } |
---|
730 | } |
---|
731 | if (returnCode < -1 && returnCode > -5) { |
---|
732 | problemStatus_ = -2; // |
---|
733 | } else if (returnCode == -5) { |
---|
734 | // something flagged - continue; |
---|
735 | } else if (returnCode == 2) { |
---|
736 | problemStatus_ = -5; // looks unbounded |
---|
737 | } else if (returnCode == 4) { |
---|
738 | problemStatus_ = -2; // looks unbounded but has iterated |
---|
739 | } else if (returnCode != -1) { |
---|
740 | assert(returnCode == 3); |
---|
741 | problemStatus_ = 3; |
---|
742 | } |
---|
743 | } else { |
---|
744 | // no pivot column |
---|
745 | #ifdef CLP_DEBUG |
---|
746 | if (handler_->logLevel() & 32) |
---|
747 | printf("** no column pivot\n"); |
---|
748 | #endif |
---|
749 | if (pivotMode < 10) { |
---|
750 | // looks optimal |
---|
751 | primalColumnPivot_->setLooksOptimal(true); |
---|
752 | } else { |
---|
753 | pivotMode--; |
---|
754 | #ifdef CLP_DEBUG |
---|
755 | if (handler_->logLevel() & 32) |
---|
756 | printf("pivot mode now %d\n", pivotMode); |
---|
757 | #endif |
---|
758 | if (pivotMode == 9) |
---|
759 | pivotMode = 0; // switch off fast attempt |
---|
760 | unflag(); |
---|
761 | } |
---|
762 | if (nonLinearCost_->numberInfeasibilities()) |
---|
763 | problemStatus_ = -4; // might be infeasible |
---|
764 | returnCode = 0; |
---|
765 | break; |
---|
766 | } |
---|
767 | } |
---|
768 | delete[] array1; |
---|
769 | return returnCode; |
---|
770 | } |
---|
771 | // Creates direction vector |
---|
772 | void ClpSimplexNonlinear::directionVector(CoinIndexedVector *vectorArray, |
---|
773 | CoinIndexedVector *spare1, CoinIndexedVector *spare2, |
---|
774 | int pivotMode2, |
---|
775 | double &normFlagged, double &normUnflagged, |
---|
776 | int &numberNonBasic) |
---|
777 | { |
---|
778 | #if CLP_DEBUG > 1 |
---|
779 | vectorArray->checkClear(); |
---|
780 | spare1->checkClear(); |
---|
781 | spare2->checkClear(); |
---|
782 | #endif |
---|
783 | double *array = vectorArray->denseVector(); |
---|
784 | int *index = vectorArray->getIndices(); |
---|
785 | int number = 0; |
---|
786 | sequenceIn_ = -1; |
---|
787 | normFlagged = 0.0; |
---|
788 | normUnflagged = 1.0; |
---|
789 | double dualTolerance2 = CoinMin(1.0e-8, 1.0e-2 * dualTolerance_); |
---|
790 | double dualTolerance3 = CoinMin(1.0e-2, 1.0e3 * dualTolerance_); |
---|
791 | if (!numberNonBasic) { |
---|
792 | //if (nonLinearCost_->sumInfeasibilities()>1.0e-4) |
---|
793 | //printf("infeasible\n"); |
---|
794 | if (!pivotMode2 || pivotMode2 >= 10) { |
---|
795 | normUnflagged = 0.0; |
---|
796 | double bestDj = 0.0; |
---|
797 | double bestSuper = 0.0; |
---|
798 | double sumSuper = 0.0; |
---|
799 | sequenceIn_ = -1; |
---|
800 | int nSuper = 0; |
---|
801 | for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) { |
---|
802 | array[iSequence] = 0.0; |
---|
803 | if (flagged(iSequence)) { |
---|
804 | // accumulate norm |
---|
805 | switch (getStatus(iSequence)) { |
---|
806 | |
---|
807 | case basic: |
---|
808 | case ClpSimplex::isFixed: |
---|
809 | break; |
---|
810 | case atUpperBound: |
---|
811 | if (dj_[iSequence] > dualTolerance3) |
---|
812 | normFlagged += dj_[iSequence] * dj_[iSequence]; |
---|
813 | break; |
---|
814 | case atLowerBound: |
---|
815 | if (dj_[iSequence] < -dualTolerance3) |
---|
816 | normFlagged += dj_[iSequence] * dj_[iSequence]; |
---|
817 | break; |
---|
818 | case isFree: |
---|
819 | case superBasic: |
---|
820 | if (fabs(dj_[iSequence]) > dualTolerance3) |
---|
821 | normFlagged += dj_[iSequence] * dj_[iSequence]; |
---|
822 | break; |
---|
823 | } |
---|
824 | continue; |
---|
825 | } |
---|
826 | switch (getStatus(iSequence)) { |
---|
827 | |
---|
828 | case basic: |
---|
829 | case ClpSimplex::isFixed: |
---|
830 | break; |
---|
831 | case atUpperBound: |
---|
832 | if (dj_[iSequence] > dualTolerance_) { |
---|
833 | if (dj_[iSequence] > dualTolerance3) |
---|
834 | normUnflagged += dj_[iSequence] * dj_[iSequence]; |
---|
835 | if (pivotMode2 < 10) { |
---|
836 | array[iSequence] = -dj_[iSequence]; |
---|
837 | index[number++] = iSequence; |
---|
838 | } else { |
---|
839 | if (dj_[iSequence] > bestDj) { |
---|
840 | bestDj = dj_[iSequence]; |
---|
841 | sequenceIn_ = iSequence; |
---|
842 | } |
---|
843 | } |
---|
844 | } |
---|
845 | break; |
---|
846 | case atLowerBound: |
---|
847 | if (dj_[iSequence] < -dualTolerance_) { |
---|
848 | if (dj_[iSequence] < -dualTolerance3) |
---|
849 | normUnflagged += dj_[iSequence] * dj_[iSequence]; |
---|
850 | if (pivotMode2 < 10) { |
---|
851 | array[iSequence] = -dj_[iSequence]; |
---|
852 | index[number++] = iSequence; |
---|
853 | } else { |
---|
854 | if (-dj_[iSequence] > bestDj) { |
---|
855 | bestDj = -dj_[iSequence]; |
---|
856 | sequenceIn_ = iSequence; |
---|
857 | } |
---|
858 | } |
---|
859 | } |
---|
860 | break; |
---|
861 | case isFree: |
---|
862 | case superBasic: |
---|
863 | //#define ALLSUPER |
---|
864 | #define NOSUPER |
---|
865 | #ifndef ALLSUPER |
---|
866 | if (fabs(dj_[iSequence]) > dualTolerance_) { |
---|
867 | if (fabs(dj_[iSequence]) > dualTolerance3) |
---|
868 | normUnflagged += dj_[iSequence] * dj_[iSequence]; |
---|
869 | nSuper++; |
---|
870 | bestSuper = CoinMax(fabs(dj_[iSequence]), bestSuper); |
---|
871 | sumSuper += fabs(dj_[iSequence]); |
---|
872 | } |
---|
873 | if (fabs(dj_[iSequence]) > dualTolerance2) { |
---|
874 | array[iSequence] = -dj_[iSequence]; |
---|
875 | index[number++] = iSequence; |
---|
876 | } |
---|
877 | #else |
---|
878 | array[iSequence] = -dj_[iSequence]; |
---|
879 | index[number++] = iSequence; |
---|
880 | if (pivotMode2 >= 10) |
---|
881 | bestSuper = CoinMax(fabs(dj_[iSequence]), bestSuper); |
---|
882 | #endif |
---|
883 | break; |
---|
884 | } |
---|
885 | } |
---|
886 | #ifdef NOSUPER |
---|
887 | // redo |
---|
888 | bestSuper = sumSuper; |
---|
889 | if (sequenceIn_ >= 0 && bestDj > bestSuper) { |
---|
890 | int j; |
---|
891 | // get rid of superbasics |
---|
892 | for (j = 0; j < number; j++) { |
---|
893 | int iSequence = index[j]; |
---|
894 | array[iSequence] = 0.0; |
---|
895 | } |
---|
896 | number = 0; |
---|
897 | array[sequenceIn_] = -dj_[sequenceIn_]; |
---|
898 | index[number++] = sequenceIn_; |
---|
899 | } else { |
---|
900 | sequenceIn_ = -1; |
---|
901 | } |
---|
902 | #else |
---|
903 | if (pivotMode2 >= 10 || !nSuper) { |
---|
904 | bool takeBest = true; |
---|
905 | if (pivotMode2 == 100 && nSuper > 1) |
---|
906 | takeBest = false; |
---|
907 | if (sequenceIn_ >= 0 && takeBest) { |
---|
908 | if (fabs(dj_[sequenceIn_]) > bestSuper) { |
---|
909 | array[sequenceIn_] = -dj_[sequenceIn_]; |
---|
910 | index[number++] = sequenceIn_; |
---|
911 | } else { |
---|
912 | sequenceIn_ = -1; |
---|
913 | } |
---|
914 | } else { |
---|
915 | sequenceIn_ = -1; |
---|
916 | } |
---|
917 | } |
---|
918 | #endif |
---|
919 | #ifdef CLP_DEBUG |
---|
920 | if (handler_->logLevel() & 32) { |
---|
921 | if (sequenceIn_ >= 0) |
---|
922 | printf("%d superBasic, chosen %d - dj %g\n", nSuper, sequenceIn_, |
---|
923 | dj_[sequenceIn_]); |
---|
924 | else |
---|
925 | printf("%d superBasic - none chosen\n", nSuper); |
---|
926 | } |
---|
927 | #endif |
---|
928 | } else { |
---|
929 | double bestDj = 0.0; |
---|
930 | double saveDj = 0.0; |
---|
931 | if (sequenceOut_ >= 0) { |
---|
932 | saveDj = dj_[sequenceOut_]; |
---|
933 | dj_[sequenceOut_] = 0.0; |
---|
934 | switch (getStatus(sequenceOut_)) { |
---|
935 | |
---|
936 | case basic: |
---|
937 | sequenceOut_ = -1; |
---|
938 | case ClpSimplex::isFixed: |
---|
939 | break; |
---|
940 | case atUpperBound: |
---|
941 | if (dj_[sequenceOut_] > dualTolerance_) { |
---|
942 | #ifdef CLP_DEBUG |
---|
943 | if (handler_->logLevel() & 32) |
---|
944 | printf("after pivot out %d values %g %g %g, dj %g\n", |
---|
945 | sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_], |
---|
946 | upper_[sequenceOut_], dj_[sequenceOut_]); |
---|
947 | #endif |
---|
948 | } |
---|
949 | break; |
---|
950 | case atLowerBound: |
---|
951 | if (dj_[sequenceOut_] < -dualTolerance_) { |
---|
952 | #ifdef CLP_DEBUG |
---|
953 | if (handler_->logLevel() & 32) |
---|
954 | printf("after pivot out %d values %g %g %g, dj %g\n", |
---|
955 | sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_], |
---|
956 | upper_[sequenceOut_], dj_[sequenceOut_]); |
---|
957 | #endif |
---|
958 | } |
---|
959 | break; |
---|
960 | case isFree: |
---|
961 | case superBasic: |
---|
962 | if (dj_[sequenceOut_] > dualTolerance_) { |
---|
963 | #ifdef CLP_DEBUG |
---|
964 | if (handler_->logLevel() & 32) |
---|
965 | printf("after pivot out %d values %g %g %g, dj %g\n", |
---|
966 | sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_], |
---|
967 | upper_[sequenceOut_], dj_[sequenceOut_]); |
---|
968 | #endif |
---|
969 | } else if (dj_[sequenceOut_] < -dualTolerance_) { |
---|
970 | #ifdef CLP_DEBUG |
---|
971 | if (handler_->logLevel() & 32) |
---|
972 | printf("after pivot out %d values %g %g %g, dj %g\n", |
---|
973 | sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_], |
---|
974 | upper_[sequenceOut_], dj_[sequenceOut_]); |
---|
975 | #endif |
---|
976 | } |
---|
977 | break; |
---|
978 | } |
---|
979 | } |
---|
980 | // Go for dj |
---|
981 | pivotMode2 = 3; |
---|
982 | for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) { |
---|
983 | array[iSequence] = 0.0; |
---|
984 | if (flagged(iSequence)) |
---|
985 | continue; |
---|
986 | switch (getStatus(iSequence)) { |
---|
987 | |
---|
988 | case basic: |
---|
989 | case ClpSimplex::isFixed: |
---|
990 | break; |
---|
991 | case atUpperBound: |
---|
992 | if (dj_[iSequence] > dualTolerance_) { |
---|
993 | double distance = CoinMin(1.0e-2, solution_[iSequence] - lower_[iSequence]); |
---|
994 | double merit = distance * dj_[iSequence]; |
---|
995 | if (pivotMode2 == 1) |
---|
996 | merit *= 1.0e-20; // discourage |
---|
997 | if (pivotMode2 == 3) |
---|
998 | merit = fabs(dj_[iSequence]); |
---|
999 | if (merit > bestDj) { |
---|
1000 | sequenceIn_ = iSequence; |
---|
1001 | bestDj = merit; |
---|
1002 | } |
---|
1003 | } |
---|
1004 | break; |
---|
1005 | case atLowerBound: |
---|
1006 | if (dj_[iSequence] < -dualTolerance_) { |
---|
1007 | double distance = CoinMin(1.0e-2, upper_[iSequence] - solution_[iSequence]); |
---|
1008 | double merit = -distance * dj_[iSequence]; |
---|
1009 | if (pivotMode2 == 1) |
---|
1010 | merit *= 1.0e-20; // discourage |
---|
1011 | if (pivotMode2 == 3) |
---|
1012 | merit = fabs(dj_[iSequence]); |
---|
1013 | if (merit > bestDj) { |
---|
1014 | sequenceIn_ = iSequence; |
---|
1015 | bestDj = merit; |
---|
1016 | } |
---|
1017 | } |
---|
1018 | break; |
---|
1019 | case isFree: |
---|
1020 | case superBasic: |
---|
1021 | if (dj_[iSequence] > dualTolerance_) { |
---|
1022 | double distance = CoinMin(1.0e-2, solution_[iSequence] - lower_[iSequence]); |
---|
1023 | distance = CoinMin(solution_[iSequence] - lower_[iSequence], |
---|
1024 | upper_[iSequence] - solution_[iSequence]); |
---|
1025 | double merit = distance * dj_[iSequence]; |
---|
1026 | if (pivotMode2 == 1) |
---|
1027 | merit = distance; |
---|
1028 | if (pivotMode2 == 3) |
---|
1029 | merit = fabs(dj_[iSequence]); |
---|
1030 | if (merit > bestDj) { |
---|
1031 | sequenceIn_ = iSequence; |
---|
1032 | bestDj = merit; |
---|
1033 | } |
---|
1034 | } else if (dj_[iSequence] < -dualTolerance_) { |
---|
1035 | double distance = CoinMin(1.0e-2, upper_[iSequence] - solution_[iSequence]); |
---|
1036 | distance = CoinMin(solution_[iSequence] - lower_[iSequence], |
---|
1037 | upper_[iSequence] - solution_[iSequence]); |
---|
1038 | double merit = -distance * dj_[iSequence]; |
---|
1039 | if (pivotMode2 == 1) |
---|
1040 | merit = distance; |
---|
1041 | if (pivotMode2 == 3) |
---|
1042 | merit = fabs(dj_[iSequence]); |
---|
1043 | if (merit > bestDj) { |
---|
1044 | sequenceIn_ = iSequence; |
---|
1045 | bestDj = merit; |
---|
1046 | } |
---|
1047 | } |
---|
1048 | break; |
---|
1049 | } |
---|
1050 | } |
---|
1051 | if (sequenceOut_ >= 0) { |
---|
1052 | dj_[sequenceOut_] = saveDj; |
---|
1053 | sequenceOut_ = -1; |
---|
1054 | } |
---|
1055 | if (sequenceIn_ >= 0) { |
---|
1056 | array[sequenceIn_] = -dj_[sequenceIn_]; |
---|
1057 | index[number++] = sequenceIn_; |
---|
1058 | } |
---|
1059 | } |
---|
1060 | numberNonBasic = number; |
---|
1061 | } else { |
---|
1062 | // compute norms |
---|
1063 | normUnflagged = 0.0; |
---|
1064 | for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) { |
---|
1065 | if (flagged(iSequence)) { |
---|
1066 | // accumulate norm |
---|
1067 | switch (getStatus(iSequence)) { |
---|
1068 | |
---|
1069 | case basic: |
---|
1070 | case ClpSimplex::isFixed: |
---|
1071 | break; |
---|
1072 | case atUpperBound: |
---|
1073 | if (dj_[iSequence] > dualTolerance_) |
---|
1074 | normFlagged += dj_[iSequence] * dj_[iSequence]; |
---|
1075 | break; |
---|
1076 | case atLowerBound: |
---|
1077 | if (dj_[iSequence] < -dualTolerance_) |
---|
1078 | normFlagged += dj_[iSequence] * dj_[iSequence]; |
---|
1079 | break; |
---|
1080 | case isFree: |
---|
1081 | case superBasic: |
---|
1082 | if (fabs(dj_[iSequence]) > dualTolerance_) |
---|
1083 | normFlagged += dj_[iSequence] * dj_[iSequence]; |
---|
1084 | break; |
---|
1085 | } |
---|
1086 | } |
---|
1087 | } |
---|
1088 | // re-use list |
---|
1089 | number = 0; |
---|
1090 | int j; |
---|
1091 | for (j = 0; j < numberNonBasic; j++) { |
---|
1092 | int iSequence = index[j]; |
---|
1093 | if (flagged(iSequence)) |
---|
1094 | continue; |
---|
1095 | switch (getStatus(iSequence)) { |
---|
1096 | |
---|
1097 | case basic: |
---|
1098 | case ClpSimplex::isFixed: |
---|
1099 | continue; //abort(); |
---|
1100 | break; |
---|
1101 | case atUpperBound: |
---|
1102 | if (dj_[iSequence] > dualTolerance_) { |
---|
1103 | number++; |
---|
1104 | normUnflagged += dj_[iSequence] * dj_[iSequence]; |
---|
1105 | } |
---|
1106 | break; |
---|
1107 | case atLowerBound: |
---|
1108 | if (dj_[iSequence] < -dualTolerance_) { |
---|
1109 | number++; |
---|
1110 | normUnflagged += dj_[iSequence] * dj_[iSequence]; |
---|
1111 | } |
---|
1112 | break; |
---|
1113 | case isFree: |
---|
1114 | case superBasic: |
---|
1115 | if (fabs(dj_[iSequence]) > dualTolerance_) { |
---|
1116 | number++; |
---|
1117 | normUnflagged += dj_[iSequence] * dj_[iSequence]; |
---|
1118 | } |
---|
1119 | break; |
---|
1120 | } |
---|
1121 | array[iSequence] = -dj_[iSequence]; |
---|
1122 | } |
---|
1123 | // switch to large |
---|
1124 | normUnflagged = 1.0; |
---|
1125 | if (!number) { |
---|
1126 | for (j = 0; j < numberNonBasic; j++) { |
---|
1127 | int iSequence = index[j]; |
---|
1128 | array[iSequence] = 0.0; |
---|
1129 | } |
---|
1130 | numberNonBasic = 0; |
---|
1131 | } |
---|
1132 | number = numberNonBasic; |
---|
1133 | } |
---|
1134 | if (number) { |
---|
1135 | int j; |
---|
1136 | // Now do basic ones - how do I compensate for small basic infeasibilities? |
---|
1137 | int iRow; |
---|
1138 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1139 | int iPivot = pivotVariable_[iRow]; |
---|
1140 | double value = 0.0; |
---|
1141 | if (solution_[iPivot] > upper_[iPivot]) { |
---|
1142 | value = upper_[iPivot] - solution_[iPivot]; |
---|
1143 | } else if (solution_[iPivot] < lower_[iPivot]) { |
---|
1144 | value = lower_[iPivot] - solution_[iPivot]; |
---|
1145 | } |
---|
1146 | //if (value) |
---|
1147 | //printf("inf %d %g %g %g\n",iPivot,lower_[iPivot],solution_[iPivot], |
---|
1148 | // upper_[iPivot]); |
---|
1149 | //value=0.0; |
---|
1150 | value *= -1.0e0; |
---|
1151 | if (value) { |
---|
1152 | array[iPivot] = value; |
---|
1153 | index[number++] = iPivot; |
---|
1154 | } |
---|
1155 | } |
---|
1156 | double *array2 = spare1->denseVector(); |
---|
1157 | int *index2 = spare1->getIndices(); |
---|
1158 | int number2 = 0; |
---|
1159 | times(-1.0, array, array2); |
---|
1160 | array = array + numberColumns_; |
---|
1161 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1162 | double value = array2[iRow] + array[iRow]; |
---|
1163 | if (value) { |
---|
1164 | array2[iRow] = value; |
---|
1165 | index2[number2++] = iRow; |
---|
1166 | } else { |
---|
1167 | array2[iRow] = 0.0; |
---|
1168 | } |
---|
1169 | } |
---|
1170 | array -= numberColumns_; |
---|
1171 | spare1->setNumElements(number2); |
---|
1172 | // Ftran |
---|
1173 | factorization_->updateColumn(spare2, spare1); |
---|
1174 | number2 = spare1->getNumElements(); |
---|
1175 | for (j = 0; j < number2; j++) { |
---|
1176 | int iSequence = index2[j]; |
---|
1177 | double value = array2[iSequence]; |
---|
1178 | array2[iSequence] = 0.0; |
---|
1179 | if (value) { |
---|
1180 | int iPivot = pivotVariable_[iSequence]; |
---|
1181 | double oldValue = array[iPivot]; |
---|
1182 | if (!oldValue) { |
---|
1183 | array[iPivot] = value; |
---|
1184 | index[number++] = iPivot; |
---|
1185 | } else { |
---|
1186 | // something already there |
---|
1187 | array[iPivot] = value + oldValue; |
---|
1188 | } |
---|
1189 | } |
---|
1190 | } |
---|
1191 | spare1->setNumElements(0); |
---|
1192 | } |
---|
1193 | vectorArray->setNumElements(number); |
---|
1194 | } |
---|
1195 | #define MINTYPE 1 |
---|
1196 | #if MINTYPE == 2 |
---|
1197 | static double |
---|
1198 | innerProductIndexed(const double *region1, int size, const double *region2, const int *which) |
---|
1199 | { |
---|
1200 | int i; |
---|
1201 | double value = 0.0; |
---|
1202 | for (i = 0; i < size; i++) { |
---|
1203 | int j = which[i]; |
---|
1204 | value += region1[j] * region2[j]; |
---|
1205 | } |
---|
1206 | return value; |
---|
1207 | } |
---|
1208 | #endif |
---|
1209 | /* |
---|
1210 | Row array and column array have direction |
---|
1211 | Returns 0 - can do normal iteration (basis change) |
---|
1212 | 1 - no basis change |
---|
1213 | */ |
---|
1214 | int ClpSimplexNonlinear::pivotColumn(CoinIndexedVector *longArray, |
---|
1215 | CoinIndexedVector *rowArray, |
---|
1216 | CoinIndexedVector *columnArray, |
---|
1217 | CoinIndexedVector *spare, |
---|
1218 | int &pivotMode, |
---|
1219 | double &solutionError, |
---|
1220 | double *dArray) |
---|
1221 | { |
---|
1222 | // say not optimal |
---|
1223 | primalColumnPivot_->setLooksOptimal(false); |
---|
1224 | double acceptablePivot = 1.0e-10; |
---|
1225 | int lastSequenceIn = -1; |
---|
1226 | if (pivotMode && pivotMode < 10) { |
---|
1227 | acceptablePivot = 1.0e-6; |
---|
1228 | if (factorization_->pivots()) |
---|
1229 | acceptablePivot = 1.0e-5; // if we have iterated be more strict |
---|
1230 | } |
---|
1231 | double acceptableBasic = 1.0e-7; |
---|
1232 | |
---|
1233 | int number = longArray->getNumElements(); |
---|
1234 | int numberTotal = numberRows_ + numberColumns_; |
---|
1235 | int bestSequence = -1; |
---|
1236 | int bestBasicSequence = -1; |
---|
1237 | double eps = 1.0e-1; |
---|
1238 | eps = 1.0e-6; |
---|
1239 | double basicTheta = 1.0e30; |
---|
1240 | double objTheta = 0.0; |
---|
1241 | bool finished = false; |
---|
1242 | sequenceIn_ = -1; |
---|
1243 | int nPasses = 0; |
---|
1244 | int nTotalPasses = 0; |
---|
1245 | int nBigPasses = 0; |
---|
1246 | double djNorm0 = 0.0; |
---|
1247 | double djNorm = 0.0; |
---|
1248 | double normFlagged = 0.0; |
---|
1249 | double normUnflagged = 0.0; |
---|
1250 | int localPivotMode = pivotMode; |
---|
1251 | bool allFinished = false; |
---|
1252 | bool justOne = false; |
---|
1253 | int returnCode = 1; |
---|
1254 | double currentObj; |
---|
1255 | double predictedObj; |
---|
1256 | double thetaObj; |
---|
1257 | objective_->stepLength(this, solution_, solution_, 0.0, |
---|
1258 | currentObj, predictedObj, thetaObj); |
---|
1259 | double saveObj = currentObj; |
---|
1260 | #if MINTYPE == 2 |
---|
1261 | // try Shanno's method |
---|
1262 | //would be memory leak |
---|
1263 | //double * saveY=new double[numberTotal]; |
---|
1264 | //double * saveS=new double[numberTotal]; |
---|
1265 | //double * saveY2=new double[numberTotal]; |
---|
1266 | //double * saveS2=new double[numberTotal]; |
---|
1267 | double saveY[100]; |
---|
1268 | double saveS[100]; |
---|
1269 | double saveY2[100]; |
---|
1270 | double saveS2[100]; |
---|
1271 | double zz[10000]; |
---|
1272 | #endif |
---|
1273 | double *dArray2 = dArray + numberTotal; |
---|
1274 | // big big loop |
---|
1275 | while (!allFinished) { |
---|
1276 | double *work = longArray->denseVector(); |
---|
1277 | int *which = longArray->getIndices(); |
---|
1278 | allFinished = true; |
---|
1279 | // CONJUGATE 0 - never, 1 as pivotMode, 2 as localPivotMode, 3 always |
---|
1280 | //#define SMALLTHETA1 1.0e-25 |
---|
1281 | //#define SMALLTHETA2 1.0e-25 |
---|
1282 | #define SMALLTHETA1 1.0e-10 |
---|
1283 | #define SMALLTHETA2 1.0e-10 |
---|
1284 | #define CONJUGATE 2 |
---|
1285 | #if CONJUGATE == 0 |
---|
1286 | int conjugate = 0; |
---|
1287 | #elif CONJUGATE == 1 |
---|
1288 | int conjugate = (pivotMode < 10) ? MINTYPE : 0; |
---|
1289 | #elif CONJUGATE == 2 |
---|
1290 | int conjugate = MINTYPE; |
---|
1291 | #else |
---|
1292 | int conjugate = MINTYPE; |
---|
1293 | #endif |
---|
1294 | |
---|
1295 | //if (pivotMode==1) |
---|
1296 | //localPivotMode=11;; |
---|
1297 | #if CLP_DEBUG > 1 |
---|
1298 | longArray->checkClear(); |
---|
1299 | #endif |
---|
1300 | int numberNonBasic = 0; |
---|
1301 | int startLocalMode = -1; |
---|
1302 | while (!finished) { |
---|
1303 | double simpleObjective = COIN_DBL_MAX; |
---|
1304 | returnCode = 1; |
---|
1305 | int iSequence; |
---|
1306 | objective_->reducedGradient(this, dj_, false); |
---|
1307 | // get direction vector in longArray |
---|
1308 | longArray->clear(); |
---|
1309 | // take out comment to try slightly different idea |
---|
1310 | if (nPasses + nTotalPasses > 3000 || nBigPasses > 100) { |
---|
1311 | if (factorization_->pivots()) |
---|
1312 | returnCode = 3; |
---|
1313 | break; |
---|
1314 | } |
---|
1315 | if (!nPasses) { |
---|
1316 | numberNonBasic = 0; |
---|
1317 | nBigPasses++; |
---|
1318 | } |
---|
1319 | // just do superbasic if in cleanup mode |
---|
1320 | int local = localPivotMode; |
---|
1321 | if (!local && pivotMode >= 10 && nBigPasses < 10) { |
---|
1322 | local = 100; |
---|
1323 | } else if (local == -1 || nBigPasses >= 10) { |
---|
1324 | local = 0; |
---|
1325 | localPivotMode = 0; |
---|
1326 | } |
---|
1327 | if (justOne) { |
---|
1328 | local = 2; |
---|
1329 | //local=100; |
---|
1330 | justOne = false; |
---|
1331 | } |
---|
1332 | if (!nPasses) |
---|
1333 | startLocalMode = local; |
---|
1334 | directionVector(longArray, spare, rowArray, local, |
---|
1335 | normFlagged, normUnflagged, numberNonBasic); |
---|
1336 | { |
---|
1337 | // check null vector |
---|
1338 | double *rhs = spare->denseVector(); |
---|
1339 | #if CLP_DEBUG > 1 |
---|
1340 | spare->checkClear(); |
---|
1341 | #endif |
---|
1342 | int iRow; |
---|
1343 | multiplyAdd(solution_ + numberColumns_, numberRows_, -1.0, rhs, 0.0); |
---|
1344 | matrix_->times(1.0, solution_, rhs, rowScale_, columnScale_); |
---|
1345 | double largest = 0.0; |
---|
1346 | #if CLP_DEBUG > 0 |
---|
1347 | int iLargest = -1; |
---|
1348 | #endif |
---|
1349 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1350 | double value = fabs(rhs[iRow]); |
---|
1351 | rhs[iRow] = 0.0; |
---|
1352 | if (value > largest) { |
---|
1353 | largest = value; |
---|
1354 | #if CLP_DEBUG > 0 |
---|
1355 | iLargest = iRow; |
---|
1356 | #endif |
---|
1357 | } |
---|
1358 | } |
---|
1359 | #if CLP_DEBUG > 0 |
---|
1360 | if ((handler_->logLevel() & 32) && largest > 1.0e-8) |
---|
1361 | printf("largest rhs error %g on row %d\n", largest, iLargest); |
---|
1362 | #endif |
---|
1363 | if (solutionError < 0.0) { |
---|
1364 | solutionError = largest; |
---|
1365 | } else if (largest > CoinMax(1.0e-8, 1.0e2 * solutionError) && factorization_->pivots()) { |
---|
1366 | longArray->clear(); |
---|
1367 | pivotRow_ = -1; |
---|
1368 | theta_ = 0.0; |
---|
1369 | return 3; |
---|
1370 | } |
---|
1371 | } |
---|
1372 | if (sequenceIn_ >= 0) |
---|
1373 | lastSequenceIn = sequenceIn_; |
---|
1374 | #if MINTYPE != 2 |
---|
1375 | double djNormSave = djNorm; |
---|
1376 | #endif |
---|
1377 | djNorm = 0.0; |
---|
1378 | int iIndex; |
---|
1379 | for (iIndex = 0; iIndex < numberNonBasic; iIndex++) { |
---|
1380 | int iSequence = which[iIndex]; |
---|
1381 | double alpha = work[iSequence]; |
---|
1382 | djNorm += alpha * alpha; |
---|
1383 | } |
---|
1384 | // go to conjugate gradient if necessary |
---|
1385 | if (numberNonBasic && localPivotMode >= 10 && (nPasses > 4 || sequenceIn_ < 0)) { |
---|
1386 | localPivotMode = 0; |
---|
1387 | nTotalPasses += nPasses; |
---|
1388 | nPasses = 0; |
---|
1389 | } |
---|
1390 | #if CONJUGATE == 2 |
---|
1391 | conjugate = (localPivotMode < 10) ? MINTYPE : 0; |
---|
1392 | #endif |
---|
1393 | //printf("bigP %d pass %d nBasic %d norm %g normI %g normF %g\n", |
---|
1394 | // nBigPasses,nPasses,numberNonBasic,normUnflagged,normFlagged); |
---|
1395 | if (!nPasses) { |
---|
1396 | //printf("numberNon %d\n",numberNonBasic); |
---|
1397 | #if MINTYPE == 2 |
---|
1398 | assert(numberNonBasic < 100); |
---|
1399 | memset(zz, 0, numberNonBasic * numberNonBasic * sizeof(double)); |
---|
1400 | int put = 0; |
---|
1401 | for (int iVariable = 0; iVariable < numberNonBasic; iVariable++) { |
---|
1402 | zz[put] = 1.0; |
---|
1403 | put += numberNonBasic + 1; |
---|
1404 | } |
---|
1405 | #endif |
---|
1406 | djNorm0 = CoinMax(djNorm, 1.0e-20); |
---|
1407 | CoinMemcpyN(work, numberTotal, dArray); |
---|
1408 | CoinMemcpyN(work, numberTotal, dArray2); |
---|
1409 | if (sequenceIn_ >= 0 && numberNonBasic == 1) { |
---|
1410 | // see if simple move |
---|
1411 | double objTheta2 = objective_->stepLength(this, solution_, work, 1.0e30, |
---|
1412 | currentObj, predictedObj, thetaObj); |
---|
1413 | rowArray->clear(); |
---|
1414 | |
---|
1415 | // update the incoming column |
---|
1416 | unpackPacked(rowArray); |
---|
1417 | factorization_->updateColumnFT(spare, rowArray); |
---|
1418 | theta_ = 0.0; |
---|
1419 | double *work2 = rowArray->denseVector(); |
---|
1420 | int number = rowArray->getNumElements(); |
---|
1421 | int *which2 = rowArray->getIndices(); |
---|
1422 | int iIndex; |
---|
1423 | bool easyMove = false; |
---|
1424 | double way; |
---|
1425 | if (dj_[sequenceIn_] > 0.0) |
---|
1426 | way = -1.0; |
---|
1427 | else |
---|
1428 | way = 1.0; |
---|
1429 | double largest = COIN_DBL_MAX; |
---|
1430 | #ifdef CLP_DEBUG |
---|
1431 | int kPivot = -1; |
---|
1432 | #endif |
---|
1433 | for (iIndex = 0; iIndex < number; iIndex++) { |
---|
1434 | int iRow = which2[iIndex]; |
---|
1435 | double alpha = way * work2[iIndex]; |
---|
1436 | int iPivot = pivotVariable_[iRow]; |
---|
1437 | if (alpha < -1.0e-5) { |
---|
1438 | double distance = upper_[iPivot] - solution_[iPivot]; |
---|
1439 | if (distance < -largest * alpha) { |
---|
1440 | #ifdef CLP_DEBUG |
---|
1441 | kPivot = iPivot; |
---|
1442 | #endif |
---|
1443 | largest = CoinMax(0.0, -distance / alpha); |
---|
1444 | } |
---|
1445 | if (distance < -1.0e-12 * alpha) { |
---|
1446 | easyMove = true; |
---|
1447 | break; |
---|
1448 | } |
---|
1449 | } else if (alpha > 1.0e-5) { |
---|
1450 | double distance = solution_[iPivot] - lower_[iPivot]; |
---|
1451 | if (distance < largest * alpha) { |
---|
1452 | #ifdef CLP_DEBUG |
---|
1453 | kPivot = iPivot; |
---|
1454 | #endif |
---|
1455 | largest = CoinMax(0.0, distance / alpha); |
---|
1456 | } |
---|
1457 | if (distance < 1.0e-12 * alpha) { |
---|
1458 | easyMove = true; |
---|
1459 | break; |
---|
1460 | } |
---|
1461 | } |
---|
1462 | } |
---|
1463 | // But largest has to match up with change |
---|
1464 | assert(work[sequenceIn_]); |
---|
1465 | largest /= fabs(work[sequenceIn_]); |
---|
1466 | if (largest < objTheta2) { |
---|
1467 | easyMove = true; |
---|
1468 | } else if (!easyMove) { |
---|
1469 | double objDrop = currentObj - predictedObj; |
---|
1470 | double th = objective_->stepLength(this, solution_, work, largest, |
---|
1471 | currentObj, predictedObj, simpleObjective); |
---|
1472 | simpleObjective = CoinMax(simpleObjective, predictedObj); |
---|
1473 | double easyDrop = currentObj - simpleObjective; |
---|
1474 | if (easyDrop > 1.0e-8 && easyDrop > 0.5 * objDrop) { |
---|
1475 | easyMove = true; |
---|
1476 | #ifdef CLP_DEBUG |
---|
1477 | if (handler_->logLevel() & 32) |
---|
1478 | printf("easy - obj drop %g, easy drop %g\n", objDrop, easyDrop); |
---|
1479 | #endif |
---|
1480 | if (easyDrop > objDrop) { |
---|
1481 | // debug |
---|
1482 | printf("****** th %g simple %g\n", th, simpleObjective); |
---|
1483 | objective_->stepLength(this, solution_, work, 1.0e30, |
---|
1484 | currentObj, predictedObj, simpleObjective); |
---|
1485 | objective_->stepLength(this, solution_, work, largest, |
---|
1486 | currentObj, predictedObj, simpleObjective); |
---|
1487 | } |
---|
1488 | } |
---|
1489 | } |
---|
1490 | rowArray->clear(); |
---|
1491 | #ifdef CLP_DEBUG |
---|
1492 | if (handler_->logLevel() & 32) |
---|
1493 | printf("largest %g piv %d\n", largest, kPivot); |
---|
1494 | #endif |
---|
1495 | if (easyMove) { |
---|
1496 | valueIn_ = solution_[sequenceIn_]; |
---|
1497 | dualIn_ = dj_[sequenceIn_]; |
---|
1498 | lowerIn_ = lower_[sequenceIn_]; |
---|
1499 | upperIn_ = upper_[sequenceIn_]; |
---|
1500 | if (dualIn_ > 0.0) |
---|
1501 | directionIn_ = -1; |
---|
1502 | else |
---|
1503 | directionIn_ = 1; |
---|
1504 | longArray->clear(); |
---|
1505 | pivotRow_ = -1; |
---|
1506 | theta_ = 0.0; |
---|
1507 | return 0; |
---|
1508 | } |
---|
1509 | } |
---|
1510 | } else { |
---|
1511 | #if MINTYPE == 1 |
---|
1512 | if (conjugate) { |
---|
1513 | double djNorm2 = djNorm; |
---|
1514 | if (numberNonBasic && false) { |
---|
1515 | int iIndex; |
---|
1516 | djNorm2 = 0.0; |
---|
1517 | for (iIndex = 0; iIndex < numberNonBasic; iIndex++) { |
---|
1518 | int iSequence = which[iIndex]; |
---|
1519 | double alpha = work[iSequence]; |
---|
1520 | //djNorm2 += alpha*alpha; |
---|
1521 | double alpha2 = work[iSequence] - dArray2[iSequence]; |
---|
1522 | djNorm2 += alpha * alpha2; |
---|
1523 | } |
---|
1524 | //printf("a %.18g b %.18g\n",djNorm,djNorm2); |
---|
1525 | } |
---|
1526 | djNorm = djNorm2; |
---|
1527 | double beta = djNorm2 / djNormSave; |
---|
1528 | // reset beta every so often |
---|
1529 | //if (numberNonBasic&&nPasses>numberNonBasic&&(nPasses%(3*numberNonBasic))==1) |
---|
1530 | //beta=0.0; |
---|
1531 | //printf("current norm %g, old %g - beta %g\n", |
---|
1532 | // djNorm,djNormSave,beta); |
---|
1533 | for (iSequence = 0; iSequence < numberTotal; iSequence++) { |
---|
1534 | dArray[iSequence] = work[iSequence] + beta * dArray[iSequence]; |
---|
1535 | dArray2[iSequence] = work[iSequence]; |
---|
1536 | } |
---|
1537 | } else { |
---|
1538 | for (iSequence = 0; iSequence < numberTotal; iSequence++) |
---|
1539 | dArray[iSequence] = work[iSequence]; |
---|
1540 | } |
---|
1541 | #else |
---|
1542 | int number2 = numberNonBasic; |
---|
1543 | if (number2) { |
---|
1544 | // pack down into dArray |
---|
1545 | int jLast = -1; |
---|
1546 | for (iSequence = 0; iSequence < numberNonBasic; iSequence++) { |
---|
1547 | int j = which[iSequence]; |
---|
1548 | assert(j > jLast); |
---|
1549 | jLast = j; |
---|
1550 | double value = work[j]; |
---|
1551 | dArray[iSequence] = -value; |
---|
1552 | } |
---|
1553 | // see whether to restart |
---|
1554 | // check signs - gradient |
---|
1555 | double g1 = innerProduct(dArray, number2, dArray); |
---|
1556 | double g2 = innerProduct(dArray, number2, saveY2); |
---|
1557 | // Get differences |
---|
1558 | for (iSequence = 0; iSequence < numberNonBasic; iSequence++) { |
---|
1559 | saveY2[iSequence] = dArray[iSequence] - saveY2[iSequence]; |
---|
1560 | saveS2[iSequence] = solution_[iSequence] - saveS2[iSequence]; |
---|
1561 | } |
---|
1562 | double g3 = innerProduct(saveS2, number2, saveY2); |
---|
1563 | printf("inner %g\n", g3); |
---|
1564 | //assert(g3>0); |
---|
1565 | double zzz[10000]; |
---|
1566 | int iVariable; |
---|
1567 | g2 = 1.0e50; // temp |
---|
1568 | if (fabs(g2) >= 0.2 * fabs(g1)) { |
---|
1569 | // restart |
---|
1570 | double delta = innerProduct(saveY2, number2, saveS2) / innerProduct(saveY2, number2, saveY2); |
---|
1571 | delta = 1.0; //temp |
---|
1572 | memset(zz, 0, number2 * sizeof(double)); |
---|
1573 | int put = 0; |
---|
1574 | for (iVariable = 0; iVariable < number2; iVariable++) { |
---|
1575 | zz[put] = delta; |
---|
1576 | put += number2 + 1; |
---|
1577 | } |
---|
1578 | } else { |
---|
1579 | } |
---|
1580 | CoinMemcpyN(zz, number2 * number2, zzz); |
---|
1581 | double ww[100]; |
---|
1582 | // get sk -Hkyk |
---|
1583 | for (iVariable = 0; iVariable < number2; iVariable++) { |
---|
1584 | double value = 0.0; |
---|
1585 | for (int jVariable = 0; jVariable < number2; jVariable++) { |
---|
1586 | value += saveY2[jVariable] * zzz[iVariable + number2 * jVariable]; |
---|
1587 | } |
---|
1588 | ww[iVariable] = saveS2[iVariable] - value; |
---|
1589 | } |
---|
1590 | double ys = innerProduct(saveY2, number2, saveS2); |
---|
1591 | double multiplier1 = 1.0 / ys; |
---|
1592 | double multiplier2 = innerProduct(saveY2, number2, ww) / (ys * ys); |
---|
1593 | #if 1 |
---|
1594 | // and second way |
---|
1595 | // Hy |
---|
1596 | double h[100]; |
---|
1597 | for (iVariable = 0; iVariable < number2; iVariable++) { |
---|
1598 | double value = 0.0; |
---|
1599 | for (int jVariable = 0; jVariable < number2; jVariable++) { |
---|
1600 | value += saveY2[jVariable] * zzz[iVariable + number2 * jVariable]; |
---|
1601 | } |
---|
1602 | h[iVariable] = value; |
---|
1603 | } |
---|
1604 | double hh[10000]; |
---|
1605 | double yhy1 = innerProduct(h, number2, saveY2) * multiplier1 + 1.0; |
---|
1606 | yhy1 *= multiplier1; |
---|
1607 | for (iVariable = 0; iVariable < number2; iVariable++) { |
---|
1608 | for (int jVariable = 0; jVariable < number2; jVariable++) { |
---|
1609 | int put = iVariable + number2 * jVariable; |
---|
1610 | double value = zzz[put]; |
---|
1611 | value += yhy1 * saveS2[iVariable] * saveS2[jVariable]; |
---|
1612 | hh[put] = value; |
---|
1613 | } |
---|
1614 | } |
---|
1615 | for (iVariable = 0; iVariable < number2; iVariable++) { |
---|
1616 | for (int jVariable = 0; jVariable < number2; jVariable++) { |
---|
1617 | int put = iVariable + number2 * jVariable; |
---|
1618 | double value = hh[put]; |
---|
1619 | value -= multiplier1 * (saveS2[iVariable] * h[jVariable]); |
---|
1620 | value -= multiplier1 * (saveS2[jVariable] * h[iVariable]); |
---|
1621 | hh[put] = value; |
---|
1622 | } |
---|
1623 | } |
---|
1624 | #endif |
---|
1625 | // now update H |
---|
1626 | for (iVariable = 0; iVariable < number2; iVariable++) { |
---|
1627 | for (int jVariable = 0; jVariable < number2; jVariable++) { |
---|
1628 | int put = iVariable + number2 * jVariable; |
---|
1629 | double value = zzz[put]; |
---|
1630 | value += multiplier1 * (ww[iVariable] * saveS2[jVariable] + ww[jVariable] * saveS2[iVariable]); |
---|
1631 | value -= multiplier2 * saveS2[iVariable] * saveS2[jVariable]; |
---|
1632 | zzz[put] = value; |
---|
1633 | } |
---|
1634 | } |
---|
1635 | //memcpy(zzz,hh,size*sizeof(double)); |
---|
1636 | // do search direction |
---|
1637 | memset(dArray, 0, numberTotal * sizeof(double)); |
---|
1638 | for (iVariable = 0; iVariable < numberNonBasic; iVariable++) { |
---|
1639 | double value = 0.0; |
---|
1640 | for (int jVariable = 0; jVariable < number2; jVariable++) { |
---|
1641 | int k = which[jVariable]; |
---|
1642 | value += work[k] * zzz[iVariable + number2 * jVariable]; |
---|
1643 | } |
---|
1644 | int i = which[iVariable]; |
---|
1645 | dArray[i] = value; |
---|
1646 | } |
---|
1647 | // Now fill out dArray |
---|
1648 | { |
---|
1649 | int j; |
---|
1650 | // Now do basic ones |
---|
1651 | int iRow; |
---|
1652 | CoinIndexedVector *spare1 = spare; |
---|
1653 | CoinIndexedVector *spare2 = rowArray; |
---|
1654 | #if CLP_DEBUG > 1 |
---|
1655 | spare1->checkClear(); |
---|
1656 | spare2->checkClear(); |
---|
1657 | #endif |
---|
1658 | double *array2 = spare1->denseVector(); |
---|
1659 | int *index2 = spare1->getIndices(); |
---|
1660 | int number2 = 0; |
---|
1661 | times(-1.0, dArray, array2); |
---|
1662 | dArray = dArray + numberColumns_; |
---|
1663 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1664 | double value = array2[iRow] + dArray[iRow]; |
---|
1665 | if (value) { |
---|
1666 | array2[iRow] = value; |
---|
1667 | index2[number2++] = iRow; |
---|
1668 | } else { |
---|
1669 | array2[iRow] = 0.0; |
---|
1670 | } |
---|
1671 | } |
---|
1672 | dArray -= numberColumns_; |
---|
1673 | spare1->setNumElements(number2); |
---|
1674 | // Ftran |
---|
1675 | factorization_->updateColumn(spare2, spare1); |
---|
1676 | number2 = spare1->getNumElements(); |
---|
1677 | for (j = 0; j < number2; j++) { |
---|
1678 | int iSequence = index2[j]; |
---|
1679 | double value = array2[iSequence]; |
---|
1680 | array2[iSequence] = 0.0; |
---|
1681 | if (value) { |
---|
1682 | int iPivot = pivotVariable_[iSequence]; |
---|
1683 | double oldValue = dArray[iPivot]; |
---|
1684 | dArray[iPivot] = value + oldValue; |
---|
1685 | } |
---|
1686 | } |
---|
1687 | spare1->setNumElements(0); |
---|
1688 | } |
---|
1689 | //assert (innerProductIndexed(dArray,number2,work,which)>0); |
---|
1690 | //printf ("innerP1 %g\n",innerProduct(dArray,numberTotal,work)); |
---|
1691 | printf("innerP1 %g innerP2 %g\n", innerProduct(dArray, numberTotal, work), |
---|
1692 | innerProductIndexed(dArray, numberNonBasic, work, which)); |
---|
1693 | assert(innerProduct(dArray, numberTotal, work) > 0); |
---|
1694 | #if 1 |
---|
1695 | { |
---|
1696 | // check null vector |
---|
1697 | int iRow; |
---|
1698 | double qq[10]; |
---|
1699 | memset(qq, 0, numberRows_ * sizeof(double)); |
---|
1700 | multiplyAdd(dArray + numberColumns_, numberRows_, -1.0, qq, 0.0); |
---|
1701 | matrix_->times(1.0, dArray, qq, rowScale_, columnScale_); |
---|
1702 | double largest = 0.0; |
---|
1703 | int iLargest = -1; |
---|
1704 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
1705 | double value = fabs(qq[iRow]); |
---|
1706 | if (value > largest) { |
---|
1707 | largest = value; |
---|
1708 | iLargest = iRow; |
---|
1709 | } |
---|
1710 | } |
---|
1711 | printf("largest null error %g on row %d\n", largest, iLargest); |
---|
1712 | for (iSequence = 0; iSequence < numberTotal; iSequence++) { |
---|
1713 | if (getStatus(iSequence) == basic) |
---|
1714 | assert(fabs(dj_[iSequence]) < 1.0e-3); |
---|
1715 | } |
---|
1716 | } |
---|
1717 | #endif |
---|
1718 | CoinMemcpyN(saveY2, numberNonBasic, saveY); |
---|
1719 | CoinMemcpyN(saveS2, numberNonBasic, saveS); |
---|
1720 | } |
---|
1721 | #endif |
---|
1722 | } |
---|
1723 | #if MINTYPE == 2 |
---|
1724 | for (iSequence = 0; iSequence < numberNonBasic; iSequence++) { |
---|
1725 | int j = which[iSequence]; |
---|
1726 | saveY2[iSequence] = -work[j]; |
---|
1727 | saveS2[iSequence] = solution_[j]; |
---|
1728 | } |
---|
1729 | #endif |
---|
1730 | if (djNorm < eps * djNorm0 || (nPasses > 100 && djNorm < CoinMin(1.0e-1 * djNorm0, 1.0e-12))) { |
---|
1731 | #ifdef CLP_DEBUG |
---|
1732 | if (handler_->logLevel() & 32) |
---|
1733 | printf("dj norm reduced from %g to %g\n", djNorm0, djNorm); |
---|
1734 | #endif |
---|
1735 | if (pivotMode < 10 || !numberNonBasic) { |
---|
1736 | finished = true; |
---|
1737 | } else { |
---|
1738 | localPivotMode = pivotMode; |
---|
1739 | nTotalPasses += nPasses; |
---|
1740 | nPasses = 0; |
---|
1741 | continue; |
---|
1742 | } |
---|
1743 | } |
---|
1744 | //if (nPasses==100||nPasses==50) |
---|
1745 | //printf("pass %d dj norm reduced from %g to %g - flagged norm %g\n",nPasses,djNorm0,djNorm, |
---|
1746 | // normFlagged); |
---|
1747 | if (nPasses > 100 && djNorm < 1.0e-2 * normFlagged && !startLocalMode) { |
---|
1748 | #ifdef CLP_DEBUG |
---|
1749 | if (handler_->logLevel() & 32) |
---|
1750 | printf("dj norm reduced from %g to %g - flagged norm %g - unflagging\n", djNorm0, djNorm, |
---|
1751 | normFlagged); |
---|
1752 | #endif |
---|
1753 | unflag(); |
---|
1754 | localPivotMode = 0; |
---|
1755 | nTotalPasses += nPasses; |
---|
1756 | nPasses = 0; |
---|
1757 | continue; |
---|
1758 | } |
---|
1759 | if (djNorm > 0.99 * djNorm0 && nPasses > 1500) { |
---|
1760 | finished = true; |
---|
1761 | #ifdef CLP_DEBUG |
---|
1762 | if (handler_->logLevel() & 32) |
---|
1763 | printf("dj norm NOT reduced from %g to %g\n", djNorm0, djNorm); |
---|
1764 | #endif |
---|
1765 | djNorm = 1.2345e-20; |
---|
1766 | } |
---|
1767 | number = longArray->getNumElements(); |
---|
1768 | if (!numberNonBasic) { |
---|
1769 | // looks optimal |
---|
1770 | // check if any dj look plausible |
---|
1771 | int nSuper = 0; |
---|
1772 | int nFlagged = 0; |
---|
1773 | int nNormal = 0; |
---|
1774 | for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) { |
---|
1775 | if (flagged(iSequence)) { |
---|
1776 | switch (getStatus(iSequence)) { |
---|
1777 | |
---|
1778 | case basic: |
---|
1779 | case ClpSimplex::isFixed: |
---|
1780 | break; |
---|
1781 | case atUpperBound: |
---|
1782 | if (dj_[iSequence] > dualTolerance_) |
---|
1783 | nFlagged++; |
---|
1784 | break; |
---|
1785 | case atLowerBound: |
---|
1786 | if (dj_[iSequence] < -dualTolerance_) |
---|
1787 | nFlagged++; |
---|
1788 | break; |
---|
1789 | case isFree: |
---|
1790 | case superBasic: |
---|
1791 | if (fabs(dj_[iSequence]) > dualTolerance_) |
---|
1792 | nFlagged++; |
---|
1793 | break; |
---|
1794 | } |
---|
1795 | continue; |
---|
1796 | } |
---|
1797 | switch (getStatus(iSequence)) { |
---|
1798 | |
---|
1799 | case basic: |
---|
1800 | case ClpSimplex::isFixed: |
---|
1801 | break; |
---|
1802 | case atUpperBound: |
---|
1803 | if (dj_[iSequence] > dualTolerance_) |
---|
1804 | nNormal++; |
---|
1805 | break; |
---|
1806 | case atLowerBound: |
---|
1807 | if (dj_[iSequence] < -dualTolerance_) |
---|
1808 | nNormal++; |
---|
1809 | break; |
---|
1810 | case isFree: |
---|
1811 | case superBasic: |
---|
1812 | if (fabs(dj_[iSequence]) > dualTolerance_) |
---|
1813 | nSuper++; |
---|
1814 | break; |
---|
1815 | } |
---|
1816 | } |
---|
1817 | #ifdef CLP_DEBUG |
---|
1818 | if (handler_->logLevel() & 32) |
---|
1819 | printf("%d super, %d normal, %d flagged\n", |
---|
1820 | nSuper, nNormal, nFlagged); |
---|
1821 | #endif |
---|
1822 | |
---|
1823 | int nFlagged2 = 1; |
---|
1824 | if (lastSequenceIn < 0 && !nNormal && !nSuper) { |
---|
1825 | nFlagged2 = unflag(); |
---|
1826 | if (pivotMode >= 10) { |
---|
1827 | pivotMode--; |
---|
1828 | #ifdef CLP_DEBUG |
---|
1829 | if (handler_->logLevel() & 32) |
---|
1830 | printf("pivot mode now %d\n", pivotMode); |
---|
1831 | #endif |
---|
1832 | if (pivotMode == 9) |
---|
1833 | pivotMode = 0; // switch off fast attempt |
---|
1834 | } |
---|
1835 | } else { |
---|
1836 | lastSequenceIn = -1; |
---|
1837 | } |
---|
1838 | if (!nFlagged2 || (normFlagged + normUnflagged < 1.0e-8)) { |
---|
1839 | primalColumnPivot_->setLooksOptimal(true); |
---|
1840 | return 0; |
---|
1841 | } else { |
---|
1842 | localPivotMode = -1; |
---|
1843 | nTotalPasses += nPasses; |
---|
1844 | nPasses = 0; |
---|
1845 | finished = false; |
---|
1846 | continue; |
---|
1847 | } |
---|
1848 | } |
---|
1849 | bestSequence = -1; |
---|
1850 | bestBasicSequence = -1; |
---|
1851 | |
---|
1852 | // temp |
---|
1853 | nPasses++; |
---|
1854 | if (nPasses > 2000) |
---|
1855 | finished = true; |
---|
1856 | double theta = 1.0e30; |
---|
1857 | basicTheta = 1.0e30; |
---|
1858 | theta_ = 1.0e30; |
---|
1859 | double basicTolerance = 1.0e-4 * primalTolerance_; |
---|
1860 | for (iSequence = 0; iSequence < numberTotal; iSequence++) { |
---|
1861 | //if (flagged(iSequence) |
---|
1862 | // continue; |
---|
1863 | double alpha = dArray[iSequence]; |
---|
1864 | Status thisStatus = getStatus(iSequence); |
---|
1865 | double oldValue = solution_[iSequence]; |
---|
1866 | if (thisStatus != basic) { |
---|
1867 | if (fabs(alpha) >= acceptablePivot) { |
---|
1868 | if (alpha < 0.0) { |
---|
1869 | // variable going towards lower bound |
---|
1870 | double bound = lower_[iSequence]; |
---|
1871 | oldValue -= bound; |
---|
1872 | if (oldValue + theta * alpha < 0.0) { |
---|
1873 | bestSequence = iSequence; |
---|
1874 | theta = CoinMax(0.0, oldValue / (-alpha)); |
---|
1875 | } |
---|
1876 | } else { |
---|
1877 | // variable going towards upper bound |
---|
1878 | double bound = upper_[iSequence]; |
---|
1879 | oldValue = bound - oldValue; |
---|
1880 | if (oldValue - theta * alpha < 0.0) { |
---|
1881 | bestSequence = iSequence; |
---|
1882 | theta = CoinMax(0.0, oldValue / alpha); |
---|
1883 | } |
---|
1884 | } |
---|
1885 | } |
---|
1886 | } else { |
---|
1887 | if (fabs(alpha) >= acceptableBasic) { |
---|
1888 | if (alpha < 0.0) { |
---|
1889 | // variable going towards lower bound |
---|
1890 | double bound = lower_[iSequence]; |
---|
1891 | oldValue -= bound; |
---|
1892 | if (oldValue + basicTheta * alpha < -basicTolerance) { |
---|
1893 | bestBasicSequence = iSequence; |
---|
1894 | basicTheta = CoinMax(0.0, (oldValue + basicTolerance) / (-alpha)); |
---|
1895 | } |
---|
1896 | } else { |
---|
1897 | // variable going towards upper bound |
---|
1898 | double bound = upper_[iSequence]; |
---|
1899 | oldValue = bound - oldValue; |
---|
1900 | if (oldValue - basicTheta * alpha < -basicTolerance) { |
---|
1901 | bestBasicSequence = iSequence; |
---|
1902 | basicTheta = CoinMax(0.0, (oldValue + basicTolerance) / alpha); |
---|
1903 | } |
---|
1904 | } |
---|
1905 | } |
---|
1906 | } |
---|
1907 | } |
---|
1908 | theta_ = CoinMin(theta, basicTheta); |
---|
1909 | // Now find minimum of function |
---|
1910 | double objTheta2 = objective_->stepLength(this, solution_, dArray, CoinMin(theta, basicTheta), |
---|
1911 | currentObj, predictedObj, thetaObj); |
---|
1912 | #ifdef CLP_DEBUG |
---|
1913 | if (handler_->logLevel() & 32) |
---|
1914 | printf("current obj %g thetaObj %g, predictedObj %g\n", currentObj, thetaObj, predictedObj); |
---|
1915 | #endif |
---|
1916 | objTheta2 = CoinMin(objTheta2, 1.0e29); |
---|
1917 | #if MINTYPE == 1 |
---|
1918 | if (conjugate) { |
---|
1919 | double offset; |
---|
1920 | const double *gradient = objective_->gradient(this, |
---|
1921 | dArray, offset, |
---|
1922 | true, 0); |
---|
1923 | double product = 0.0; |
---|
1924 | for (iSequence = 0; iSequence < numberColumns_; iSequence++) { |
---|
1925 | double alpha = dArray[iSequence]; |
---|
1926 | double value = alpha * gradient[iSequence]; |
---|
1927 | product += value; |
---|
1928 | } |
---|
1929 | //#define INCLUDESLACK |
---|
1930 | #ifdef INCLUDESLACK |
---|
1931 | for (; iSequence < numberColumns_ + numberRows_; iSequence++) { |
---|
1932 | double alpha = dArray[iSequence]; |
---|
1933 | double value = alpha * cost_[iSequence]; |
---|
1934 | product += value; |
---|
1935 | } |
---|
1936 | #endif |
---|
1937 | if (product > 0.0) |
---|
1938 | objTheta = djNorm / product; |
---|
1939 | else |
---|
1940 | objTheta = COIN_DBL_MAX; |
---|
1941 | #ifndef NDEBUG |
---|
1942 | if (product < -1.0e-8 && handler_->logLevel() > 1) |
---|
1943 | printf("bad product %g\n", product); |
---|
1944 | #endif |
---|
1945 | product = CoinMax(product, 0.0); |
---|
1946 | } else { |
---|
1947 | objTheta = objTheta2; |
---|
1948 | } |
---|
1949 | #else |
---|
1950 | objTheta = objTheta2; |
---|
1951 | #endif |
---|
1952 | // if very small difference then take pivot (depends on djNorm?) |
---|
1953 | // distinguish between basic and non-basic |
---|
1954 | bool chooseObjTheta = objTheta < theta_; |
---|
1955 | if (chooseObjTheta) { |
---|
1956 | if (thetaObj < currentObj - 1.0e-12 && objTheta + 1.0e-10 > theta_) |
---|
1957 | chooseObjTheta = false; |
---|
1958 | //if (thetaObj<currentObj+1.0e-12&&objTheta+1.0e-5>theta_) |
---|
1959 | //chooseObjTheta=false; |
---|
1960 | } |
---|
1961 | //if (objTheta+SMALLTHETA1<theta_||(thetaObj>currentObj+difference&&objTheta<theta_)) { |
---|
1962 | if (chooseObjTheta) { |
---|
1963 | theta_ = objTheta; |
---|
1964 | } else { |
---|
1965 | objTheta = CoinMax(objTheta, 1.00000001 * theta_ + 1.0e-12); |
---|
1966 | //if (theta+1.0e-13>basicTheta) { |
---|
1967 | //theta = CoinMax(theta,1.00000001*basicTheta); |
---|
1968 | //theta_ = basicTheta; |
---|
1969 | //} |
---|
1970 | } |
---|
1971 | // always out if one variable in and zero move |
---|
1972 | if (theta_ == basicTheta || (sequenceIn_ >= 0 && theta_ < 1.0e-10)) |
---|
1973 | finished = true; // come out |
---|
1974 | #ifdef CLP_DEBUG |
---|
1975 | if (handler_->logLevel() & 32) { |
---|
1976 | printf("%d pass,", nPasses); |
---|
1977 | if (sequenceIn_ >= 0) |
---|
1978 | printf(" Sin %d,", sequenceIn_); |
---|
1979 | if (basicTheta == theta_) |
---|
1980 | printf(" X(%d) basicTheta %g", bestBasicSequence, basicTheta); |
---|
1981 | else |
---|
1982 | printf(" basicTheta %g", basicTheta); |
---|
1983 | if (theta == theta_) |
---|
1984 | printf(" X(%d) non-basicTheta %g", bestSequence, theta); |
---|
1985 | else |
---|
1986 | printf(" non-basicTheta %g", theta); |
---|
1987 | printf(" %s objTheta %g", objTheta == theta_ ? "X" : "", objTheta); |
---|
1988 | printf(" djNorm %g\n", djNorm); |
---|
1989 | } |
---|
1990 | #endif |
---|
1991 | if (handler_->logLevel() > 3 && objTheta != theta_) { |
---|
1992 | printf("%d pass obj %g,", nPasses, currentObj); |
---|
1993 | if (sequenceIn_ >= 0) |
---|
1994 | printf(" Sin %d,", sequenceIn_); |
---|
1995 | if (basicTheta == theta_) |
---|
1996 | printf(" X(%d) basicTheta %g", bestBasicSequence, basicTheta); |
---|
1997 | else |
---|
1998 | printf(" basicTheta %g", basicTheta); |
---|
1999 | if (theta == theta_) |
---|
2000 | printf(" X(%d) non-basicTheta %g", bestSequence, theta); |
---|
2001 | else |
---|
2002 | printf(" non-basicTheta %g", theta); |
---|
2003 | printf(" %s objTheta %g", objTheta == theta_ ? "X" : "", objTheta); |
---|
2004 | printf(" djNorm %g\n", djNorm); |
---|
2005 | } |
---|
2006 | if (objTheta != theta_) { |
---|
2007 | //printf("hit boundary after %d passes\n",nPasses); |
---|
2008 | nTotalPasses += nPasses; |
---|
2009 | nPasses = 0; // start again |
---|
2010 | } |
---|
2011 | if (localPivotMode < 10 || lastSequenceIn == bestSequence) { |
---|
2012 | if (theta_ == theta && theta_ < basicTheta && theta_ < 1.0e-5) |
---|
2013 | setFlagged(bestSequence); // out of active set |
---|
2014 | } |
---|
2015 | // Update solution |
---|
2016 | for (iSequence = 0; iSequence < numberTotal; iSequence++) { |
---|
2017 | //for (iIndex=0;iIndex<number;iIndex++) { |
---|
2018 | |
---|
2019 | //int iSequence = which[iIndex]; |
---|
2020 | double alpha = dArray[iSequence]; |
---|
2021 | if (alpha) { |
---|
2022 | double value = solution_[iSequence] + theta_ * alpha; |
---|
2023 | solution_[iSequence] = value; |
---|
2024 | switch (getStatus(iSequence)) { |
---|
2025 | |
---|
2026 | case basic: |
---|
2027 | case isFixed: |
---|
2028 | case isFree: |
---|
2029 | case atUpperBound: |
---|
2030 | case atLowerBound: |
---|
2031 | nonLinearCost_->setOne(iSequence, value); |
---|
2032 | break; |
---|
2033 | case superBasic: |
---|
2034 | // To get correct action |
---|
2035 | setStatus(iSequence, isFixed); |
---|
2036 | nonLinearCost_->setOne(iSequence, value); |
---|
2037 | //assert (getStatus(iSequence)!=isFixed); |
---|
2038 | break; |
---|
2039 | } |
---|
2040 | } |
---|
2041 | } |
---|
2042 | if (objTheta < SMALLTHETA2 && objTheta == theta_) { |
---|
2043 | if (sequenceIn_ >= 0 && basicTheta < 1.0e-9) { |
---|
2044 | // try just one |
---|
2045 | localPivotMode = 0; |
---|
2046 | sequenceIn_ = -1; |
---|
2047 | nTotalPasses += nPasses; |
---|
2048 | nPasses = 0; |
---|
2049 | //finished=true; |
---|
2050 | //objTheta=0.0; |
---|
2051 | //theta_=0.0; |
---|
2052 | } else if (sequenceIn_ < 0 && nTotalPasses > 10) { |
---|
2053 | if (objTheta < 1.0e-10) { |
---|
2054 | finished = true; |
---|
2055 | //printf("zero move\n"); |
---|
2056 | break; |
---|
2057 | } |
---|
2058 | } |
---|
2059 | } |
---|
2060 | #ifdef CLP_DEBUG |
---|
2061 | if (handler_->logLevel() & 32) { |
---|
2062 | objective_->stepLength(this, solution_, work, 0.0, |
---|
2063 | currentObj, predictedObj, thetaObj); |
---|
2064 | printf("current obj %g after update - simple was %g\n", currentObj, simpleObjective); |
---|
2065 | } |
---|
2066 | #endif |
---|
2067 | if (sequenceIn_ >= 0 && !finished && objTheta > 1.0e-4) { |
---|
2068 | // we made some progress - back to normal |
---|
2069 | if (localPivotMode < 10) { |
---|
2070 | localPivotMode = 0; |
---|
2071 | sequenceIn_ = -1; |
---|
2072 | nTotalPasses += nPasses; |
---|
2073 | nPasses = 0; |
---|
2074 | } |
---|
2075 | #ifdef CLP_DEBUG |
---|
2076 | if (handler_->logLevel() & 32) |
---|
2077 | printf("some progress?\n"); |
---|
2078 | #endif |
---|
2079 | } |
---|
2080 | #if CLP_DEBUG > 1 |
---|
2081 | longArray->checkClean(); |
---|
2082 | #endif |
---|
2083 | } |
---|
2084 | #ifdef CLP_DEBUG |
---|
2085 | if (handler_->logLevel() & 32) |
---|
2086 | printf("out of loop after %d (%d) passes\n", nPasses, nTotalPasses); |
---|
2087 | #endif |
---|
2088 | if (nTotalPasses >= 1000 || (nTotalPasses > 10 && sequenceIn_ < 0 && theta_ < 1.0e-10)) |
---|
2089 | returnCode = 2; |
---|
2090 | bool ordinaryDj = false; |
---|
2091 | //if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-7&&theta_==basicTheta) |
---|
2092 | //printf("could try ordinary iteration %g\n",theta_); |
---|
2093 | if (sequenceIn_ >= 0 && numberNonBasic == 1 && theta_ < 1.0e-15) { |
---|
2094 | //printf("trying ordinary iteration\n"); |
---|
2095 | ordinaryDj = true; |
---|
2096 | } |
---|
2097 | if (!basicTheta && !ordinaryDj) { |
---|
2098 | //returnCode=2; |
---|
2099 | //objTheta=-1.0; // so we fall through |
---|
2100 | } |
---|
2101 | assert(theta_ < 1.0e30); // for now |
---|
2102 | // See if we need to pivot |
---|
2103 | if (theta_ == basicTheta || ordinaryDj) { |
---|
2104 | if (!ordinaryDj) { |
---|
2105 | // find an incoming column which will force pivot |
---|
2106 | int iRow; |
---|
2107 | pivotRow_ = -1; |
---|
2108 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
2109 | if (pivotVariable_[iRow] == bestBasicSequence) { |
---|
2110 | pivotRow_ = iRow; |
---|
2111 | break; |
---|
2112 | } |
---|
2113 | } |
---|
2114 | assert(pivotRow_ >= 0); |
---|
2115 | // Get good size for pivot |
---|
2116 | double acceptablePivot = 1.0e-7; |
---|
2117 | if (factorization_->pivots() > 10) |
---|
2118 | acceptablePivot = 1.0e-5; // if we have iterated be more strict |
---|
2119 | else if (factorization_->pivots() > 5) |
---|
2120 | acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict |
---|
2121 | // should be dArray but seems better this way! |
---|
2122 | double direction = work[bestBasicSequence] > 0.0 ? -1.0 : 1.0; |
---|
2123 | // create as packed |
---|
2124 | rowArray->createPacked(1, &pivotRow_, &direction); |
---|
2125 | factorization_->updateColumnTranspose(spare, rowArray); |
---|
2126 | // put row of tableau in rowArray and columnArray |
---|
2127 | matrix_->transposeTimes(this, -1.0, |
---|
2128 | rowArray, spare, columnArray); |
---|
2129 | // choose one futhest away from bound which has reasonable pivot |
---|
2130 | // If increasing we want negative alpha |
---|
2131 | |
---|
2132 | double *work2; |
---|
2133 | int iSection; |
---|
2134 | |
---|
2135 | sequenceIn_ = -1; |
---|
2136 | double bestValue = -1.0; |
---|
2137 | double bestDirection = 0.0; |
---|
2138 | // First pass we take correct direction and large pivots |
---|
2139 | // then correct direction |
---|
2140 | // then any |
---|
2141 | double check[] = { 1.0e-8, -1.0e-12, -1.0e30 }; |
---|
2142 | double mult[] = { 100.0, 1.0, 1.0 }; |
---|
2143 | for (int iPass = 0; iPass < 3; iPass++) { |
---|
2144 | //if (!bestValue&&iPass==2) |
---|
2145 | //bestValue=-1.0; |
---|
2146 | double acceptable = acceptablePivot * mult[iPass]; |
---|
2147 | double checkValue = check[iPass]; |
---|
2148 | for (iSection = 0; iSection < 2; iSection++) { |
---|
2149 | |
---|
2150 | int addSequence; |
---|
2151 | |
---|
2152 | if (!iSection) { |
---|
2153 | work2 = rowArray->denseVector(); |
---|
2154 | number = rowArray->getNumElements(); |
---|
2155 | which = rowArray->getIndices(); |
---|
2156 | addSequence = numberColumns_; |
---|
2157 | } else { |
---|
2158 | work2 = columnArray->denseVector(); |
---|
2159 | number = columnArray->getNumElements(); |
---|
2160 | which = columnArray->getIndices(); |
---|
2161 | addSequence = 0; |
---|
2162 | } |
---|
2163 | int i; |
---|
2164 | |
---|
2165 | for (i = 0; i < number; i++) { |
---|
2166 | int iSequence = which[i] + addSequence; |
---|
2167 | if (flagged(iSequence)) |
---|
2168 | continue; |
---|
2169 | //double distance = CoinMin(solution_[iSequence]-lower_[iSequence], |
---|
2170 | // upper_[iSequence]-solution_[iSequence]); |
---|
2171 | double alpha = work2[i]; |
---|
2172 | // should be dArray but seems better this way! |
---|
2173 | double change = work[iSequence]; |
---|
2174 | Status thisStatus = getStatus(iSequence); |
---|
2175 | double direction = 0; |
---|
2176 | ; |
---|
2177 | switch (thisStatus) { |
---|
2178 | |
---|
2179 | case basic: |
---|
2180 | case ClpSimplex::isFixed: |
---|
2181 | break; |
---|
2182 | case isFree: |
---|
2183 | case superBasic: |
---|
2184 | if (alpha < -acceptable && change > checkValue) |
---|
2185 | direction = 1.0; |
---|
2186 | else if (alpha > acceptable && change < -checkValue) |
---|
2187 | direction = -1.0; |
---|
2188 | break; |
---|
2189 | case atUpperBound: |
---|
2190 | if (alpha > acceptable && change < -checkValue) |
---|
2191 | direction = -1.0; |
---|
2192 | else if (iPass == 2 && alpha < -acceptable && change < -checkValue) |
---|
2193 | direction = 1.0; |
---|
2194 | break; |
---|
2195 | case atLowerBound: |
---|
2196 | if (alpha < -acceptable && change > checkValue) |
---|
2197 | direction = 1.0; |
---|
2198 | else if (iPass == 2 && alpha > acceptable && change > checkValue) |
---|
2199 | direction = -1.0; |
---|
2200 | break; |
---|
2201 | } |
---|
2202 | if (direction) { |
---|
2203 | if (sequenceIn_ != lastSequenceIn || localPivotMode < 10) { |
---|
2204 | if (CoinMin(solution_[iSequence] - lower_[iSequence], |
---|
2205 | upper_[iSequence] - solution_[iSequence]) |
---|
2206 | > bestValue) { |
---|
2207 | bestValue = CoinMin(solution_[iSequence] - lower_[iSequence], |
---|
2208 | upper_[iSequence] - solution_[iSequence]); |
---|
2209 | sequenceIn_ = iSequence; |
---|
2210 | bestDirection = direction; |
---|
2211 | } |
---|
2212 | } else { |
---|
2213 | // choose |
---|
2214 | bestValue = COIN_DBL_MAX; |
---|
2215 | sequenceIn_ = iSequence; |
---|
2216 | bestDirection = direction; |
---|
2217 | } |
---|
2218 | } |
---|
2219 | } |
---|
2220 | } |
---|
2221 | if (sequenceIn_ >= 0 && bestValue > 0.0) |
---|
2222 | break; |
---|
2223 | } |
---|
2224 | if (sequenceIn_ >= 0) { |
---|
2225 | valueIn_ = solution_[sequenceIn_]; |
---|
2226 | dualIn_ = dj_[sequenceIn_]; |
---|
2227 | if (bestDirection < 0.0) { |
---|
2228 | // we want positive dj |
---|
2229 | if (dualIn_ <= 0.0) { |
---|
2230 | //printf("bad dj - xx %g\n",dualIn_); |
---|
2231 | dualIn_ = 1.0; |
---|
2232 | } |
---|
2233 | } else { |
---|
2234 | // we want negative dj |
---|
2235 | if (dualIn_ >= 0.0) { |
---|
2236 | //printf("bad dj - xx %g\n",dualIn_); |
---|
2237 | dualIn_ = -1.0; |
---|
2238 | } |
---|
2239 | } |
---|
2240 | lowerIn_ = lower_[sequenceIn_]; |
---|
2241 | upperIn_ = upper_[sequenceIn_]; |
---|
2242 | if (dualIn_ > 0.0) |
---|
2243 | directionIn_ = -1; |
---|
2244 | else |
---|
2245 | directionIn_ = 1; |
---|
2246 | } else { |
---|
2247 | //ordinaryDj=true; |
---|
2248 | #ifdef CLP_DEBUG |
---|
2249 | if (handler_->logLevel() & 32) { |
---|
2250 | printf("no easy pivot - norm %g mode %d\n", djNorm, localPivotMode); |
---|
2251 | if (rowArray->getNumElements() + columnArray->getNumElements() < 12) { |
---|
2252 | for (iSection = 0; iSection < 2; iSection++) { |
---|
2253 | |
---|
2254 | int addSequence; |
---|
2255 | |
---|
2256 | if (!iSection) { |
---|
2257 | work2 = rowArray->denseVector(); |
---|
2258 | number = rowArray->getNumElements(); |
---|
2259 | which = rowArray->getIndices(); |
---|
2260 | addSequence = numberColumns_; |
---|
2261 | } else { |
---|
2262 | work2 = columnArray->denseVector(); |
---|
2263 | number = columnArray->getNumElements(); |
---|
2264 | which = columnArray->getIndices(); |
---|
2265 | addSequence = 0; |
---|
2266 | } |
---|
2267 | int i; |
---|
2268 | char section[] = { 'R', 'C' }; |
---|
2269 | for (i = 0; i < number; i++) { |
---|
2270 | int iSequence = which[i] + addSequence; |
---|
2271 | if (flagged(iSequence)) { |
---|
2272 | printf("%c%d flagged\n", section[iSection], which[i]); |
---|
2273 | continue; |
---|
2274 | } else { |
---|
2275 | printf("%c%d status %d sol %g %g %g alpha %g change %g\n", |
---|
2276 | section[iSection], which[i], status_[iSequence], |
---|
2277 | lower_[iSequence], solution_[iSequence], upper_[iSequence], |
---|
2278 | work2[i], work[iSequence]); |
---|
2279 | } |
---|
2280 | } |
---|
2281 | } |
---|
2282 | } |
---|
2283 | } |
---|
2284 | #endif |
---|
2285 | assert(sequenceIn_ < 0); |
---|
2286 | justOne = true; |
---|
2287 | allFinished = false; // Round again |
---|
2288 | finished = false; |
---|
2289 | nTotalPasses += nPasses; |
---|
2290 | nPasses = 0; |
---|
2291 | if (djNorm < 0.9 * djNorm0 && djNorm < 1.0e-3 && !localPivotMode) { |
---|
2292 | #ifdef CLP_DEBUG |
---|
2293 | if (handler_->logLevel() & 32) |
---|
2294 | printf("no pivot - mode %d norms %g %g - unflagging\n", |
---|
2295 | localPivotMode, djNorm0, djNorm); |
---|
2296 | #endif |
---|
2297 | unflag(); //unflagging |
---|
2298 | returnCode = 1; |
---|
2299 | } else { |
---|
2300 | returnCode = 2; // do single incoming |
---|
2301 | returnCode = 1; |
---|
2302 | } |
---|
2303 | } |
---|
2304 | } |
---|
2305 | rowArray->clear(); |
---|
2306 | columnArray->clear(); |
---|
2307 | longArray->clear(); |
---|
2308 | if (ordinaryDj) { |
---|
2309 | valueIn_ = solution_[sequenceIn_]; |
---|
2310 | dualIn_ = dj_[sequenceIn_]; |
---|
2311 | lowerIn_ = lower_[sequenceIn_]; |
---|
2312 | upperIn_ = upper_[sequenceIn_]; |
---|
2313 | if (dualIn_ > 0.0) |
---|
2314 | directionIn_ = -1; |
---|
2315 | else |
---|
2316 | directionIn_ = 1; |
---|
2317 | } |
---|
2318 | if (returnCode == 1) |
---|
2319 | returnCode = 0; |
---|
2320 | } else { |
---|
2321 | // round again |
---|
2322 | longArray->clear(); |
---|
2323 | if (djNorm < 1.0e-3 && !localPivotMode) { |
---|
2324 | if (djNorm == 1.2345e-20 && djNorm0 > 1.0e-4) { |
---|
2325 | #ifdef CLP_DEBUG |
---|
2326 | if (handler_->logLevel() & 32) |
---|
2327 | printf("slow convergence djNorm0 %g, %d passes, mode %d, result %d\n", djNorm0, nPasses, |
---|
2328 | localPivotMode, returnCode); |
---|
2329 | #endif |
---|
2330 | //if (!localPivotMode) |
---|
2331 | //returnCode=2; // force singleton |
---|
2332 | } else { |
---|
2333 | #ifdef CLP_DEBUG |
---|
2334 | if (handler_->logLevel() & 32) |
---|
2335 | printf("unflagging as djNorm %g %g, %d passes\n", djNorm, djNorm0, nPasses); |
---|
2336 | #endif |
---|
2337 | if (pivotMode >= 10) { |
---|
2338 | pivotMode--; |
---|
2339 | #ifdef CLP_DEBUG |
---|
2340 | if (handler_->logLevel() & 32) |
---|
2341 | printf("pivot mode now %d\n", pivotMode); |
---|
2342 | #endif |
---|
2343 | if (pivotMode == 9) |
---|
2344 | pivotMode = 0; // switch off fast attempt |
---|
2345 | } |
---|
2346 | bool unflagged = unflag() != 0; |
---|
2347 | if (!unflagged && djNorm < 1.0e-10) { |
---|
2348 | // ? declare victory |
---|
2349 | sequenceIn_ = -1; |
---|
2350 | returnCode = 0; |
---|
2351 | } |
---|
2352 | } |
---|
2353 | } |
---|
2354 | } |
---|
2355 | } |
---|
2356 | if (djNorm0 < 1.0e-12 * normFlagged) { |
---|
2357 | #ifdef CLP_DEBUG |
---|
2358 | if (handler_->logLevel() & 32) |
---|
2359 | printf("unflagging as djNorm %g %g and flagged norm %g\n", djNorm, djNorm0, normFlagged); |
---|
2360 | #endif |
---|
2361 | unflag(); |
---|
2362 | } |
---|
2363 | if (saveObj - currentObj < 1.0e-5 && nTotalPasses > 2000) { |
---|
2364 | normUnflagged = 0.0; |
---|
2365 | double dualTolerance3 = CoinMin(1.0e-2, 1.0e3 * dualTolerance_); |
---|
2366 | for (int iSequence = 0; iSequence < numberColumns_ + numberRows_; iSequence++) { |
---|
2367 | switch (getStatus(iSequence)) { |
---|
2368 | |
---|
2369 | case basic: |
---|
2370 | case ClpSimplex::isFixed: |
---|
2371 | break; |
---|
2372 | case atUpperBound: |
---|
2373 | if (dj_[iSequence] > dualTolerance3) |
---|
2374 | normFlagged += dj_[iSequence] * dj_[iSequence]; |
---|
2375 | break; |
---|
2376 | case atLowerBound: |
---|
2377 | if (dj_[iSequence] < -dualTolerance3) |
---|
2378 | normFlagged += dj_[iSequence] * dj_[iSequence]; |
---|
2379 | break; |
---|
2380 | case isFree: |
---|
2381 | case superBasic: |
---|
2382 | if (fabs(dj_[iSequence]) > dualTolerance3) |
---|
2383 | normFlagged += dj_[iSequence] * dj_[iSequence]; |
---|
2384 | break; |
---|
2385 | } |
---|
2386 | } |
---|
2387 | if (handler_->logLevel() > 2) |
---|
2388 | printf("possible optimal %d %d %g %g\n", |
---|
2389 | nBigPasses, nTotalPasses, saveObj - currentObj, normFlagged); |
---|
2390 | if (normFlagged < 1.0e-5) { |
---|
2391 | unflag(); |
---|
2392 | primalColumnPivot_->setLooksOptimal(true); |
---|
2393 | returnCode = 0; |
---|
2394 | } |
---|
2395 | } |
---|
2396 | return returnCode; |
---|
2397 | } |
---|
2398 | /* Do last half of an iteration. |
---|
2399 | Return codes |
---|
2400 | Reasons to come out normal mode |
---|
2401 | -1 normal |
---|
2402 | -2 factorize now - good iteration |
---|
2403 | -3 slight inaccuracy - refactorize - iteration done |
---|
2404 | -4 inaccuracy - refactorize - no iteration |
---|
2405 | -5 something flagged - go round again |
---|
2406 | +2 looks unbounded |
---|
2407 | +3 max iterations (iteration done) |
---|
2408 | |
---|
2409 | */ |
---|
2410 | int ClpSimplexNonlinear::pivotNonlinearResult() |
---|
2411 | { |
---|
2412 | |
---|
2413 | int returnCode = -1; |
---|
2414 | |
---|
2415 | rowArray_[1]->clear(); |
---|
2416 | |
---|
2417 | // we found a pivot column |
---|
2418 | // update the incoming column |
---|
2419 | unpackPacked(rowArray_[1]); |
---|
2420 | factorization_->updateColumnFT(rowArray_[2], rowArray_[1]); |
---|
2421 | theta_ = 0.0; |
---|
2422 | double *work = rowArray_[1]->denseVector(); |
---|
2423 | int number = rowArray_[1]->getNumElements(); |
---|
2424 | int *which = rowArray_[1]->getIndices(); |
---|
2425 | bool keepValue = false; |
---|
2426 | double saveValue = 0.0; |
---|
2427 | if (pivotRow_ >= 0) { |
---|
2428 | sequenceOut_ = pivotVariable_[pivotRow_]; |
---|
2429 | ; |
---|
2430 | valueOut_ = solution(sequenceOut_); |
---|
2431 | keepValue = true; |
---|
2432 | saveValue = valueOut_; |
---|
2433 | lowerOut_ = lower_[sequenceOut_]; |
---|
2434 | upperOut_ = upper_[sequenceOut_]; |
---|
2435 | int iIndex; |
---|
2436 | for (iIndex = 0; iIndex < number; iIndex++) { |
---|
2437 | |
---|
2438 | int iRow = which[iIndex]; |
---|
2439 | if (iRow == pivotRow_) { |
---|
2440 | alpha_ = work[iIndex]; |
---|
2441 | break; |
---|
2442 | } |
---|
2443 | } |
---|
2444 | } else { |
---|
2445 | int iIndex; |
---|
2446 | double smallest = COIN_DBL_MAX; |
---|
2447 | for (iIndex = 0; iIndex < number; iIndex++) { |
---|
2448 | int iRow = which[iIndex]; |
---|
2449 | double alpha = work[iIndex]; |
---|
2450 | if (fabs(alpha) > 1.0e-6) { |
---|
2451 | int iPivot = pivotVariable_[iRow]; |
---|
2452 | double distance = CoinMin(upper_[iPivot] - solution_[iPivot], |
---|
2453 | solution_[iPivot] - lower_[iPivot]); |
---|
2454 | if (distance < smallest) { |
---|
2455 | pivotRow_ = iRow; |
---|
2456 | alpha_ = alpha; |
---|
2457 | smallest = distance; |
---|
2458 | } |
---|
2459 | } |
---|
2460 | } |
---|
2461 | if (smallest > primalTolerance_) { |
---|
2462 | smallest = COIN_DBL_MAX; |
---|
2463 | for (iIndex = 0; iIndex < number; iIndex++) { |
---|
2464 | int iRow = which[iIndex]; |
---|
2465 | double alpha = work[iIndex]; |
---|
2466 | if (fabs(alpha) > 1.0e-6) { |
---|
2467 | double distance = randomNumberGenerator_.randomDouble(); |
---|
2468 | if (distance < smallest) { |
---|
2469 | pivotRow_ = iRow; |
---|
2470 | alpha_ = alpha; |
---|
2471 | smallest = distance; |
---|
2472 | } |
---|
2473 | } |
---|
2474 | } |
---|
2475 | } |
---|
2476 | assert(pivotRow_ >= 0); |
---|
2477 | sequenceOut_ = pivotVariable_[pivotRow_]; |
---|
2478 | ; |
---|
2479 | valueOut_ = solution(sequenceOut_); |
---|
2480 | lowerOut_ = lower_[sequenceOut_]; |
---|
2481 | upperOut_ = upper_[sequenceOut_]; |
---|
2482 | } |
---|
2483 | double newValue = valueOut_ - theta_ * alpha_; |
---|
2484 | bool isSuperBasic = false; |
---|
2485 | if (valueOut_ >= upperOut_ - primalTolerance_) { |
---|
2486 | directionOut_ = -1; // to upper bound |
---|
2487 | upperOut_ = nonLinearCost_->nearest(sequenceOut_, newValue); |
---|
2488 | upperOut_ = newValue; |
---|
2489 | } else if (valueOut_ <= lowerOut_ + primalTolerance_) { |
---|
2490 | directionOut_ = 1; // to lower bound |
---|
2491 | lowerOut_ = nonLinearCost_->nearest(sequenceOut_, newValue); |
---|
2492 | } else { |
---|
2493 | lowerOut_ = valueOut_; |
---|
2494 | upperOut_ = valueOut_; |
---|
2495 | isSuperBasic = true; |
---|
2496 | //printf("XX superbasic out\n"); |
---|
2497 | } |
---|
2498 | dualOut_ = dj_[sequenceOut_]; |
---|
2499 | // if stable replace in basis |
---|
2500 | |
---|
2501 | int updateStatus = factorization_->replaceColumn(this, |
---|
2502 | rowArray_[2], |
---|
2503 | rowArray_[1], |
---|
2504 | pivotRow_, |
---|
2505 | alpha_); |
---|
2506 | |
---|
2507 | // if no pivots, bad update but reasonable alpha - take and invert |
---|
2508 | if (updateStatus == 2 && lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e-5) |
---|
2509 | updateStatus = 4; |
---|
2510 | if (updateStatus == 1 || updateStatus == 4) { |
---|
2511 | // slight error |
---|
2512 | if (factorization_->pivots() > 5 || updateStatus == 4) { |
---|
2513 | returnCode = -3; |
---|
2514 | } |
---|
2515 | } else if (updateStatus == 2) { |
---|
2516 | // major error |
---|
2517 | // better to have small tolerance even if slower |
---|
2518 | factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15)); |
---|
2519 | int maxFactor = factorization_->maximumPivots(); |
---|
2520 | if (maxFactor > 10) { |
---|
2521 | if (forceFactorization_ < 0) |
---|
2522 | forceFactorization_ = maxFactor; |
---|
2523 | forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1)); |
---|
2524 | } |
---|
2525 | // later we may need to unwind more e.g. fake bounds |
---|
2526 | if (lastGoodIteration_ != numberIterations_) { |
---|
2527 | clearAll(); |
---|
2528 | pivotRow_ = -1; |
---|
2529 | returnCode = -4; |
---|
2530 | } else { |
---|
2531 | // need to reject something |
---|
2532 | char x = isColumn(sequenceIn_) ? 'C' : 'R'; |
---|
2533 | handler_->message(CLP_SIMPLEX_FLAG, messages_) |
---|
2534 | << x << sequenceWithin(sequenceIn_) |
---|
2535 | << CoinMessageEol; |
---|
2536 | setFlagged(sequenceIn_); |
---|
2537 | progress_.clearBadTimes(); |
---|
2538 | lastBadIteration_ = numberIterations_; // say be more cautious |
---|
2539 | clearAll(); |
---|
2540 | pivotRow_ = -1; |
---|
2541 | sequenceOut_ = -1; |
---|
2542 | returnCode = -5; |
---|
2543 | } |
---|
2544 | return returnCode; |
---|
2545 | } else if (updateStatus == 3) { |
---|
2546 | // out of memory |
---|
2547 | // increase space if not many iterations |
---|
2548 | if (factorization_->pivots() < 0.5 * factorization_->maximumPivots() && factorization_->pivots() < 200) |
---|
2549 | factorization_->areaFactor( |
---|
2550 | factorization_->areaFactor() * 1.1); |
---|
2551 | returnCode = -2; // factorize now |
---|
2552 | } else if (updateStatus == 5) { |
---|
2553 | problemStatus_ = -2; // factorize now |
---|
2554 | } |
---|
2555 | |
---|
2556 | // update primal solution |
---|
2557 | |
---|
2558 | double objectiveChange = 0.0; |
---|
2559 | // after this rowArray_[1] is not empty - used to update djs |
---|
2560 | // If pivot row >= numberRows then may be gub |
---|
2561 | updatePrimalsInPrimal(rowArray_[1], theta_, objectiveChange, 1); |
---|
2562 | |
---|
2563 | double oldValue = valueIn_; |
---|
2564 | if (directionIn_ == -1) { |
---|
2565 | // as if from upper bound |
---|
2566 | if (sequenceIn_ != sequenceOut_) { |
---|
2567 | // variable becoming basic |
---|
2568 | valueIn_ -= fabs(theta_); |
---|
2569 | } else { |
---|
2570 | valueIn_ = lowerIn_; |
---|
2571 | } |
---|
2572 | } else { |
---|
2573 | // as if from lower bound |
---|
2574 | if (sequenceIn_ != sequenceOut_) { |
---|
2575 | // variable becoming basic |
---|
2576 | valueIn_ += fabs(theta_); |
---|
2577 | } else { |
---|
2578 | valueIn_ = upperIn_; |
---|
2579 | } |
---|
2580 | } |
---|
2581 | objectiveChange += dualIn_ * (valueIn_ - oldValue); |
---|
2582 | // outgoing |
---|
2583 | if (sequenceIn_ != sequenceOut_) { |
---|
2584 | if (directionOut_ > 0) { |
---|
2585 | valueOut_ = lowerOut_; |
---|
2586 | } else { |
---|
2587 | valueOut_ = upperOut_; |
---|
2588 | } |
---|
2589 | if (valueOut_ < lower_[sequenceOut_] - primalTolerance_) |
---|
2590 | valueOut_ = lower_[sequenceOut_] - 0.9 * primalTolerance_; |
---|
2591 | else if (valueOut_ > upper_[sequenceOut_] + primalTolerance_) |
---|
2592 | valueOut_ = upper_[sequenceOut_] + 0.9 * primalTolerance_; |
---|
2593 | // may not be exactly at bound and bounds may have changed |
---|
2594 | // Make sure outgoing looks feasible |
---|
2595 | if (!isSuperBasic) |
---|
2596 | directionOut_ = nonLinearCost_->setOneOutgoing(sequenceOut_, valueOut_); |
---|
2597 | solution_[sequenceOut_] = valueOut_; |
---|
2598 | } |
---|
2599 | // change cost and bounds on incoming if primal |
---|
2600 | nonLinearCost_->setOne(sequenceIn_, valueIn_); |
---|
2601 | int whatNext = housekeeping(objectiveChange); |
---|
2602 | if (keepValue) |
---|
2603 | solution_[sequenceOut_] = saveValue; |
---|
2604 | if (isSuperBasic) |
---|
2605 | setStatus(sequenceOut_, superBasic); |
---|
2606 | //#define CLP_DEBUG |
---|
2607 | #if CLP_DEBUG > 1 |
---|
2608 | { |
---|
2609 | int ninf = matrix_->checkFeasible(this); |
---|
2610 | if (ninf) |
---|
2611 | printf("infeas %d\n", ninf); |
---|
2612 | } |
---|
2613 | #endif |
---|
2614 | if (whatNext == 1) { |
---|
2615 | returnCode = -2; // refactorize |
---|
2616 | } else if (whatNext == 2) { |
---|
2617 | // maximum iterations or equivalent |
---|
2618 | returnCode = 3; |
---|
2619 | } else if (numberIterations_ == lastGoodIteration_ + 2 * factorization_->maximumPivots()) { |
---|
2620 | // done a lot of flips - be safe |
---|
2621 | returnCode = -2; // refactorize |
---|
2622 | } |
---|
2623 | // Check event |
---|
2624 | { |
---|
2625 | int status = eventHandler_->event(ClpEventHandler::endOfIteration); |
---|
2626 | if (status >= 0) { |
---|
2627 | problemStatus_ = 5; |
---|
2628 | secondaryStatus_ = ClpEventHandler::endOfIteration; |
---|
2629 | returnCode = 4; |
---|
2630 | } |
---|
2631 | } |
---|
2632 | return returnCode; |
---|
2633 | } |
---|
2634 | // May use a cut approach for solving any LP |
---|
2635 | int ClpSimplexNonlinear::primalDualCuts(char *rowsIn, int startUp, int algorithm) |
---|
2636 | { |
---|
2637 | if (!rowsIn) { |
---|
2638 | if (algorithm > 0) |
---|
2639 | return ClpSimplex::primal(startUp); |
---|
2640 | else |
---|
2641 | //return static_cast<ClpSimplexDual *>(static_cast<ClpSimplex *>(this))->dual(startUp); |
---|
2642 | return ClpSimplex::dual(startUp); |
---|
2643 | } else { |
---|
2644 | int numberUsed = 0; |
---|
2645 | int rowsThreshold = CoinMax(100, numberRows_ / 2); |
---|
2646 | //int rowsTry=CoinMax(50,numberRows_/4); |
---|
2647 | // Just add this number of rows each time in small problem |
---|
2648 | int smallNumberRows = 2 * numberColumns_; |
---|
2649 | smallNumberRows = CoinMin(smallNumberRows, numberRows_ / 20); |
---|
2650 | // We will need arrays to choose rows to add |
---|
2651 | double *weight = new double[numberRows_]; |
---|
2652 | int *sort = new int[numberRows_ + numberColumns_]; |
---|
2653 | int *whichColumns = sort + numberRows_; |
---|
2654 | int numberSort = 0; |
---|
2655 | for (int i = 0; i < numberRows_; i++) { |
---|
2656 | if (rowsIn[i]) |
---|
2657 | numberUsed++; |
---|
2658 | } |
---|
2659 | if (numberUsed) { |
---|
2660 | // normal |
---|
2661 | } else { |
---|
2662 | // If non slack use that information |
---|
2663 | int numberBinding = 0; |
---|
2664 | numberPrimalInfeasibilities_ = 0; |
---|
2665 | sumPrimalInfeasibilities_ = 0.0; |
---|
2666 | memset(rowActivity_, 0, numberRows_ * sizeof(double)); |
---|
2667 | times(1.0, columnActivity_, rowActivity_); |
---|
2668 | for (int i = 0; i < numberRows_; i++) { |
---|
2669 | double lowerDifference = rowActivity_[i] - rowLower_[i]; |
---|
2670 | double upperDifference = rowActivity_[i] - rowUpper_[i]; |
---|
2671 | if (lowerDifference < -10 * primalTolerance_ || upperDifference > 10.0 * primalTolerance_) { |
---|
2672 | numberPrimalInfeasibilities_++; |
---|
2673 | if (lowerDifference < 0.0) |
---|
2674 | sumPrimalInfeasibilities_ -= lowerDifference; |
---|
2675 | else |
---|
2676 | sumPrimalInfeasibilities_ += upperDifference; |
---|
2677 | rowsIn[i] = 1; |
---|
2678 | } else if (getRowStatus(i) != basic) { |
---|
2679 | numberBinding++; |
---|
2680 | rowsIn[i] = 1; |
---|
2681 | } |
---|
2682 | } |
---|
2683 | if (numberBinding < rowsThreshold) { |
---|
2684 | // try |
---|
2685 | } else { |
---|
2686 | // random?? - or give up |
---|
2687 | // Set up initial list |
---|
2688 | numberSort = 0; |
---|
2689 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
2690 | weight[iRow] = 1.123e50; |
---|
2691 | if (rowLower_[iRow] == rowUpper_[iRow]) { |
---|
2692 | sort[numberSort++] = iRow; |
---|
2693 | weight[iRow] = 0.0; |
---|
2694 | } |
---|
2695 | } |
---|
2696 | numberSort /= 2; |
---|
2697 | // and pad out with random rows |
---|
2698 | double ratio = ((double)(smallNumberRows - numberSort)) / ((double)numberRows_); |
---|
2699 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
2700 | if (weight[iRow] == 1.123e50 && CoinDrand48() < ratio) |
---|
2701 | sort[numberSort++] = iRow; |
---|
2702 | } |
---|
2703 | // sort |
---|
2704 | CoinSort_2(weight, weight + numberRows_, sort); |
---|
2705 | numberSort = CoinMin(numberRows_, smallNumberRows); |
---|
2706 | memset(rowsIn, 0, numberRows_); |
---|
2707 | for (int iRow = 0; iRow < numberSort; iRow++) |
---|
2708 | rowsIn[sort[iRow]] = 1; |
---|
2709 | } |
---|
2710 | } |
---|
2711 | // see if feasible |
---|
2712 | numberPrimalInfeasibilities_ = 0; |
---|
2713 | memset(rowActivity_, 0, numberRows_ * sizeof(double)); |
---|
2714 | times(1.0, columnActivity_, rowActivity_); |
---|
2715 | for (int i = 0; i < numberRows_; i++) { |
---|
2716 | double lowerDifference = rowActivity_[i] - rowLower_[i]; |
---|
2717 | double upperDifference = rowActivity_[i] - rowUpper_[i]; |
---|
2718 | if (lowerDifference < -10 * primalTolerance_ || upperDifference > 10.0 * primalTolerance_) { |
---|
2719 | if (lowerDifference < 0.0) |
---|
2720 | sumPrimalInfeasibilities_ -= lowerDifference; |
---|
2721 | else |
---|
2722 | sumPrimalInfeasibilities_ += upperDifference; |
---|
2723 | numberPrimalInfeasibilities_++; |
---|
2724 | } |
---|
2725 | } |
---|
2726 | printf("Initial infeasibilities - %g (%d)\n", |
---|
2727 | sumPrimalInfeasibilities_, numberPrimalInfeasibilities_); |
---|
2728 | // Just do this number of passes |
---|
2729 | int maxPass = 50; |
---|
2730 | // And take out slack rows until this pass |
---|
2731 | int takeOutPass = 30; |
---|
2732 | int iPass; |
---|
2733 | |
---|
2734 | const CoinBigIndex *start = this->clpMatrix()->getVectorStarts(); |
---|
2735 | const int *length = this->clpMatrix()->getVectorLengths(); |
---|
2736 | const int *row = this->clpMatrix()->getIndices(); |
---|
2737 | problemStatus_ = 1; |
---|
2738 | for (iPass = 0; iPass < maxPass; iPass++) { |
---|
2739 | printf("Start of pass %d\n", iPass); |
---|
2740 | int numberSmallColumns = 0; |
---|
2741 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
2742 | int n = 0; |
---|
2743 | for (CoinBigIndex j = start[iColumn]; j < start[iColumn] + length[iColumn]; j++) { |
---|
2744 | int iRow = row[j]; |
---|
2745 | if (rowsIn[iRow]) |
---|
2746 | n++; |
---|
2747 | } |
---|
2748 | if (n) |
---|
2749 | whichColumns[numberSmallColumns++] = iColumn; |
---|
2750 | } |
---|
2751 | numberSort = 0; |
---|
2752 | for (int i = 0; i < numberRows_; i++) { |
---|
2753 | if (rowsIn[i]) |
---|
2754 | sort[numberSort++] = i; |
---|
2755 | } |
---|
2756 | // Create small problem |
---|
2757 | ClpSimplex small(this, numberSort, sort, numberSmallColumns, whichColumns); |
---|
2758 | printf("Small model has %d rows, %d columns and %d elements\n", |
---|
2759 | small.numberRows(), small.numberColumns(), small.getNumElements()); |
---|
2760 | small.setFactorizationFrequency(100 + numberSort / 200); |
---|
2761 | // Solve |
---|
2762 | small.setLogLevel(CoinMax(0, logLevel() - 1)); |
---|
2763 | if (iPass > 20) { |
---|
2764 | if (sumPrimalInfeasibilities_ > 1.0e-1) { |
---|
2765 | small.dual(); |
---|
2766 | } else { |
---|
2767 | small.primal(1); |
---|
2768 | } |
---|
2769 | } else { |
---|
2770 | // presolve |
---|
2771 | #if 0 |
---|
2772 | ClpSolve::SolveType method; |
---|
2773 | ClpSolve::PresolveType presolveType = ClpSolve::presolveOn; |
---|
2774 | ClpSolve solveOptions; |
---|
2775 | solveOptions.setPresolveType(presolveType, 5); |
---|
2776 | if (sumPrimalInfeasibilities_>1.0e-1) |
---|
2777 | method = ClpSolve::useDual; |
---|
2778 | else |
---|
2779 | method = ClpSolve::usePrimalorSprint; |
---|
2780 | solveOptions.setSolveType(method); |
---|
2781 | small.initialSolve(solveOptions); |
---|
2782 | #else |
---|
2783 | #if 1 |
---|
2784 | ClpPresolve *pinfo = new ClpPresolve(); |
---|
2785 | ClpSimplex *small2 = pinfo->presolvedModel(small, 1.0e-5); |
---|
2786 | assert(small2); |
---|
2787 | if (sumPrimalInfeasibilities_ > 1.0e-1) { |
---|
2788 | small2->dual(); |
---|
2789 | } else { |
---|
2790 | small2->primal(1); |
---|
2791 | } |
---|
2792 | pinfo->postsolve(true); |
---|
2793 | delete pinfo; |
---|
2794 | #else |
---|
2795 | char *types = new char[small.numberRows() + small.numberColumns()]; |
---|
2796 | memset(types, 0, small.numberRows() + small.numberColumns()); |
---|
2797 | if (sumPrimalInfeasibilities_ > 1.0e-1) { |
---|
2798 | small.miniSolve(types, types + small.numberRows(), |
---|
2799 | -1, 0); |
---|
2800 | } else { |
---|
2801 | small.miniSolve(types, types + small.numberRows(), |
---|
2802 | 1, 1); |
---|
2803 | } |
---|
2804 | delete[] types; |
---|
2805 | #endif |
---|
2806 | if (small.sumPrimalInfeasibilities() > 1.0) |
---|
2807 | small.primal(1); |
---|
2808 | #endif |
---|
2809 | } |
---|
2810 | bool dualInfeasible = (small.status() == 2); |
---|
2811 | // move solution back |
---|
2812 | const double *smallSolution = small.primalColumnSolution(); |
---|
2813 | for (int j = 0; j < numberSmallColumns; j++) { |
---|
2814 | int iColumn = whichColumns[j]; |
---|
2815 | columnActivity_[iColumn] = smallSolution[j]; |
---|
2816 | this->setColumnStatus(iColumn, small.getColumnStatus(j)); |
---|
2817 | } |
---|
2818 | for (int iRow = 0; iRow < numberSort; iRow++) { |
---|
2819 | int kRow = sort[iRow]; |
---|
2820 | this->setRowStatus(kRow, small.getRowStatus(iRow)); |
---|
2821 | } |
---|
2822 | // compute full solution |
---|
2823 | memset(rowActivity_, 0, numberRows_ * sizeof(double)); |
---|
2824 | times(1.0, columnActivity_, rowActivity_); |
---|
2825 | if (iPass != maxPass - 1) { |
---|
2826 | // Mark row as not looked at |
---|
2827 | for (int iRow = 0; iRow < numberRows_; iRow++) |
---|
2828 | weight[iRow] = 1.123e50; |
---|
2829 | // Look at rows already in small problem |
---|
2830 | int iSort; |
---|
2831 | int numberDropped = 0; |
---|
2832 | int numberKept = 0; |
---|
2833 | int numberBinding = 0; |
---|
2834 | numberPrimalInfeasibilities_ = 0; |
---|
2835 | sumPrimalInfeasibilities_ = 0.0; |
---|
2836 | bool allFeasible = small.numberIterations() == 0; |
---|
2837 | for (iSort = 0; iSort < numberSort; iSort++) { |
---|
2838 | int iRow = sort[iSort]; |
---|
2839 | //printf("%d %g %g\n",iRow,rowActivity_[iRow],small.primalRowSolution()[iSort]); |
---|
2840 | if (getRowStatus(iRow) == ClpSimplex::basic) { |
---|
2841 | // Basic - we can get rid of if early on |
---|
2842 | if (iPass < takeOutPass && !dualInfeasible) { |
---|
2843 | double infeasibility = CoinMax(rowActivity_[iRow] - rowUpper_[iRow], |
---|
2844 | rowLower_[iRow] - rowActivity_[iRow]); |
---|
2845 | weight[iRow] = -infeasibility; |
---|
2846 | if (infeasibility > primalTolerance_ && !allFeasible) { |
---|
2847 | numberPrimalInfeasibilities_++; |
---|
2848 | sumPrimalInfeasibilities_ += infeasibility; |
---|
2849 | } else { |
---|
2850 | weight[iRow] = 1.0; |
---|
2851 | numberDropped++; |
---|
2852 | } |
---|
2853 | } else { |
---|
2854 | // keep |
---|
2855 | weight[iRow] = -1.0e40; |
---|
2856 | numberKept++; |
---|
2857 | } |
---|
2858 | } else { |
---|
2859 | // keep |
---|
2860 | weight[iRow] = -1.0e50; |
---|
2861 | numberKept++; |
---|
2862 | numberBinding++; |
---|
2863 | } |
---|
2864 | } |
---|
2865 | // Now rest |
---|
2866 | for (int iRow = 0; iRow < numberRows_; iRow++) { |
---|
2867 | sort[iRow] = iRow; |
---|
2868 | if (weight[iRow] == 1.123e50) { |
---|
2869 | // not looked at yet |
---|
2870 | double infeasibility = CoinMax(rowActivity_[iRow] - rowUpper_[iRow], |
---|
2871 | rowLower_[iRow] - rowActivity_[iRow]); |
---|
2872 | weight[iRow] = -infeasibility; |
---|
2873 | if (infeasibility > primalTolerance_) { |
---|
2874 | numberPrimalInfeasibilities_++; |
---|
2875 | sumPrimalInfeasibilities_ += infeasibility; |
---|
2876 | } |
---|
2877 | } |
---|
2878 | } |
---|
2879 | // sort |
---|
2880 | CoinSort_2(weight, weight + numberRows_, sort); |
---|
2881 | numberSort = CoinMin(numberRows_, smallNumberRows + numberKept); |
---|
2882 | memset(rowsIn, 0, numberRows_); |
---|
2883 | for (int iRow = 0; iRow < numberSort; iRow++) |
---|
2884 | rowsIn[sort[iRow]] = 1; |
---|
2885 | printf("%d rows binding, %d rows kept, %d rows dropped - new size %d rows, %d columns\n", |
---|
2886 | numberBinding, numberKept, numberDropped, numberSort, numberSmallColumns); |
---|
2887 | printf("%d rows are infeasible - sum is %g\n", numberPrimalInfeasibilities_, |
---|
2888 | sumPrimalInfeasibilities_); |
---|
2889 | if (!numberPrimalInfeasibilities_) { |
---|
2890 | problemStatus_ = 0; |
---|
2891 | printf("Exiting as looks optimal\n"); |
---|
2892 | break; |
---|
2893 | } |
---|
2894 | numberPrimalInfeasibilities_ = 0; |
---|
2895 | sumPrimalInfeasibilities_ = 0.0; |
---|
2896 | for (iSort = 0; iSort < numberSort; iSort++) { |
---|
2897 | if (weight[iSort] > -1.0e30 && weight[iSort] < -1.0e-8) { |
---|
2898 | numberPrimalInfeasibilities_++; |
---|
2899 | sumPrimalInfeasibilities_ += -weight[iSort]; |
---|
2900 | } |
---|
2901 | } |
---|
2902 | printf("in small model %d rows are infeasible - sum is %g\n", numberPrimalInfeasibilities_, |
---|
2903 | sumPrimalInfeasibilities_); |
---|
2904 | } else { |
---|
2905 | // out of passes |
---|
2906 | problemStatus_ = -1; |
---|
2907 | } |
---|
2908 | } |
---|
2909 | delete[] weight; |
---|
2910 | delete[] sort; |
---|
2911 | return 0; |
---|
2912 | } |
---|
2913 | } |
---|
2914 | // A sequential LP method |
---|
2915 | int ClpSimplexNonlinear::primalSLP(int numberPasses, double deltaTolerance, |
---|
2916 | int otherOptions) |
---|
2917 | { |
---|
2918 | // Are we minimizing or maximizing |
---|
2919 | double whichWay = optimizationDirection(); |
---|
2920 | if (whichWay < 0.0) |
---|
2921 | whichWay = -1.0; |
---|
2922 | else if (whichWay > 0.0) |
---|
2923 | whichWay = 1.0; |
---|
2924 | |
---|
2925 | int numberColumns = this->numberColumns(); |
---|
2926 | int numberRows = this->numberRows(); |
---|
2927 | double *columnLower = this->columnLower(); |
---|
2928 | double *columnUpper = this->columnUpper(); |
---|
2929 | double *solution = this->primalColumnSolution(); |
---|
2930 | |
---|
2931 | if (objective_->type() < 2) { |
---|
2932 | // no nonlinear part |
---|
2933 | return ClpSimplex::primal(0); |
---|
2934 | } |
---|
2935 | // Get list of non linear columns |
---|
2936 | char *markNonlinear = new char[numberColumns]; |
---|
2937 | memset(markNonlinear, 0, numberColumns); |
---|
2938 | int numberNonLinearColumns = objective_->markNonlinear(markNonlinear); |
---|
2939 | |
---|
2940 | if (!numberNonLinearColumns) { |
---|
2941 | delete[] markNonlinear; |
---|
2942 | // no nonlinear part |
---|
2943 | return ClpSimplex::primal(0); |
---|
2944 | } |
---|
2945 | int iColumn; |
---|
2946 | int *listNonLinearColumn = new int[numberNonLinearColumns]; |
---|
2947 | numberNonLinearColumns = 0; |
---|
2948 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
2949 | if (markNonlinear[iColumn]) |
---|
2950 | listNonLinearColumn[numberNonLinearColumns++] = iColumn; |
---|
2951 | } |
---|
2952 | // Replace objective |
---|
2953 | ClpObjective *trueObjective = objective_; |
---|
2954 | objective_ = new ClpLinearObjective(NULL, numberColumns); |
---|
2955 | double *objective = this->objective(); |
---|
2956 | // See if user wants to use cuts |
---|
2957 | char *rowsIn = NULL; |
---|
2958 | if ((otherOptions & 1) != 0 || numberPasses < 0) { |
---|
2959 | numberPasses = abs(numberPasses); |
---|
2960 | rowsIn = new char[numberRows_]; |
---|
2961 | memset(rowsIn, 0, numberRows_); |
---|
2962 | } |
---|
2963 | // get feasible |
---|
2964 | if (this->status() < 0 || numberPrimalInfeasibilities()) |
---|
2965 | primalDualCuts(rowsIn, 1, 1); |
---|
2966 | // still infeasible |
---|
2967 | if (problemStatus_) { |
---|
2968 | delete[] listNonLinearColumn; |
---|
2969 | return 0; |
---|
2970 | } |
---|
2971 | algorithm_ = 1; |
---|
2972 | int jNon; |
---|
2973 | int *last[3]; |
---|
2974 | |
---|
2975 | double *trust = new double[numberNonLinearColumns]; |
---|
2976 | double *trueLower = new double[numberNonLinearColumns]; |
---|
2977 | double *trueUpper = new double[numberNonLinearColumns]; |
---|
2978 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
2979 | iColumn = listNonLinearColumn[jNon]; |
---|
2980 | trust[jNon] = 0.5; |
---|
2981 | trueLower[jNon] = columnLower[iColumn]; |
---|
2982 | trueUpper[jNon] = columnUpper[iColumn]; |
---|
2983 | if (solution[iColumn] < trueLower[jNon]) |
---|
2984 | solution[iColumn] = trueLower[jNon]; |
---|
2985 | else if (solution[iColumn] > trueUpper[jNon]) |
---|
2986 | solution[iColumn] = trueUpper[jNon]; |
---|
2987 | } |
---|
2988 | int saveLogLevel = logLevel(); |
---|
2989 | int iPass; |
---|
2990 | double lastObjective = 1.0e31; |
---|
2991 | double *saveSolution = new double[numberColumns]; |
---|
2992 | double *saveRowSolution = new double[numberRows]; |
---|
2993 | memset(saveRowSolution, 0, numberRows * sizeof(double)); |
---|
2994 | double *savePi = new double[numberRows]; |
---|
2995 | double *safeSolution = new double[numberColumns]; |
---|
2996 | unsigned char *saveStatus = new unsigned char[numberRows + numberColumns]; |
---|
2997 | #define MULTIPLE 0 |
---|
2998 | #if MULTIPLE > 2 |
---|
2999 | // Duplication but doesn't really matter |
---|
3000 | double * saveSolutionM[MULTIPLE |
---|
3001 | }; |
---|
3002 | for (jNon = 0; jNon < MULTIPLE; jNon++) { |
---|
3003 | saveSolutionM[jNon] = new double[numberColumns]; |
---|
3004 | CoinMemcpyN(solution, numberColumns, saveSolutionM); |
---|
3005 | } |
---|
3006 | #endif |
---|
3007 | double targetDrop = 1.0e31; |
---|
3008 | double objectiveOffset; |
---|
3009 | getDblParam(ClpObjOffset, objectiveOffset); |
---|
3010 | // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change |
---|
3011 | for (iPass = 0; iPass < 3; iPass++) { |
---|
3012 | last[iPass] = new int[numberNonLinearColumns]; |
---|
3013 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) |
---|
3014 | last[iPass][jNon] = 0; |
---|
3015 | } |
---|
3016 | // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing |
---|
3017 | int goodMove = -2; |
---|
3018 | char *statusCheck = new char[numberColumns]; |
---|
3019 | double *changeRegion = new double[numberColumns]; |
---|
3020 | double offset = 0.0; |
---|
3021 | double objValue = 0.0; |
---|
3022 | int exitPass = 2 * numberPasses + 10; |
---|
3023 | for (iPass = 0; iPass < numberPasses; iPass++) { |
---|
3024 | exitPass--; |
---|
3025 | // redo objective |
---|
3026 | offset = 0.0; |
---|
3027 | objValue = -objectiveOffset; |
---|
3028 | // make sure x updated |
---|
3029 | trueObjective->newXValues(); |
---|
3030 | double theta = -1.0; |
---|
3031 | double maxTheta = COIN_DBL_MAX; |
---|
3032 | //maxTheta=1.0; |
---|
3033 | if (iPass) { |
---|
3034 | int jNon = 0; |
---|
3035 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
3036 | changeRegion[iColumn] = solution[iColumn] - saveSolution[iColumn]; |
---|
3037 | double alpha = changeRegion[iColumn]; |
---|
3038 | double oldValue = saveSolution[iColumn]; |
---|
3039 | if (markNonlinear[iColumn] == 0) { |
---|
3040 | // linear |
---|
3041 | if (alpha < -1.0e-15) { |
---|
3042 | // variable going towards lower bound |
---|
3043 | double bound = columnLower[iColumn]; |
---|
3044 | oldValue -= bound; |
---|
3045 | if (oldValue + maxTheta * alpha < 0.0) { |
---|
3046 | maxTheta = CoinMax(0.0, oldValue / (-alpha)); |
---|
3047 | } |
---|
3048 | } else if (alpha > 1.0e-15) { |
---|
3049 | // variable going towards upper bound |
---|
3050 | double bound = columnUpper[iColumn]; |
---|
3051 | oldValue = bound - oldValue; |
---|
3052 | if (oldValue - maxTheta * alpha < 0.0) { |
---|
3053 | maxTheta = CoinMax(0.0, oldValue / alpha); |
---|
3054 | } |
---|
3055 | } |
---|
3056 | } else { |
---|
3057 | // nonlinear |
---|
3058 | if (alpha < -1.0e-15) { |
---|
3059 | // variable going towards lower bound |
---|
3060 | double bound = trueLower[jNon]; |
---|
3061 | oldValue -= bound; |
---|
3062 | if (oldValue + maxTheta * alpha < 0.0) { |
---|
3063 | maxTheta = CoinMax(0.0, oldValue / (-alpha)); |
---|
3064 | } |
---|
3065 | } else if (alpha > 1.0e-15) { |
---|
3066 | // variable going towards upper bound |
---|
3067 | double bound = trueUpper[jNon]; |
---|
3068 | oldValue = bound - oldValue; |
---|
3069 | if (oldValue - maxTheta * alpha < 0.0) { |
---|
3070 | maxTheta = CoinMax(0.0, oldValue / alpha); |
---|
3071 | } |
---|
3072 | } |
---|
3073 | jNon++; |
---|
3074 | } |
---|
3075 | } |
---|
3076 | // make sure both accurate |
---|
3077 | memset(rowActivity_, 0, numberRows_ * sizeof(double)); |
---|
3078 | times(1.0, solution, rowActivity_); |
---|
3079 | memset(saveRowSolution, 0, numberRows_ * sizeof(double)); |
---|
3080 | times(1.0, saveSolution, saveRowSolution); |
---|
3081 | for (int iRow = 0; iRow < numberRows; iRow++) { |
---|
3082 | double alpha = rowActivity_[iRow] - saveRowSolution[iRow]; |
---|
3083 | double oldValue = saveRowSolution[iRow]; |
---|
3084 | if (alpha < -1.0e-15) { |
---|
3085 | // variable going towards lower bound |
---|
3086 | double bound = rowLower_[iRow]; |
---|
3087 | oldValue -= bound; |
---|
3088 | if (oldValue + maxTheta * alpha < 0.0) { |
---|
3089 | maxTheta = CoinMax(0.0, oldValue / (-alpha)); |
---|
3090 | } |
---|
3091 | } else if (alpha > 1.0e-15) { |
---|
3092 | // variable going towards upper bound |
---|
3093 | double bound = rowUpper_[iRow]; |
---|
3094 | oldValue = bound - oldValue; |
---|
3095 | if (oldValue - maxTheta * alpha < 0.0) { |
---|
3096 | maxTheta = CoinMax(0.0, oldValue / alpha); |
---|
3097 | } |
---|
3098 | } |
---|
3099 | } |
---|
3100 | } else { |
---|
3101 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
---|
3102 | changeRegion[iColumn] = 0.0; |
---|
3103 | saveSolution[iColumn] = solution[iColumn]; |
---|
3104 | } |
---|
3105 | CoinMemcpyN(rowActivity_, numberRows, saveRowSolution); |
---|
3106 | } |
---|
3107 | // get current value anyway |
---|
3108 | double predictedObj, thetaObj; |
---|
3109 | double maxTheta2 = 2.0; // to work out a b c |
---|
3110 | double theta2 = trueObjective->stepLength(this, saveSolution, changeRegion, maxTheta2, |
---|
3111 | objValue, predictedObj, thetaObj); |
---|
3112 | int lastMoveStatus = goodMove; |
---|
3113 | if (goodMove >= 0) { |
---|
3114 | theta = CoinMin(theta2, maxTheta); |
---|
3115 | #ifdef CLP_DEBUG |
---|
3116 | if (handler_->logLevel() & 32) |
---|
3117 | printf("theta %g, current %g, at maxtheta %g, predicted %g\n", |
---|
3118 | theta, objValue, thetaObj, predictedObj); |
---|
3119 | #endif |
---|
3120 | if (theta > 0.0 && theta <= 1.0) { |
---|
3121 | // update solution |
---|
3122 | double lambda = 1.0 - theta; |
---|
3123 | for (iColumn = 0; iColumn < numberColumns; iColumn++) |
---|
3124 | solution[iColumn] = lambda * saveSolution[iColumn] |
---|
3125 | + theta * solution[iColumn]; |
---|
3126 | memset(rowActivity_, 0, numberRows_ * sizeof(double)); |
---|
3127 | times(1.0, solution, rowActivity_); |
---|
3128 | if (lambda > 0.999) { |
---|
3129 | CoinMemcpyN(savePi, numberRows, this->dualRowSolution()); |
---|
3130 | CoinMemcpyN(saveStatus, numberRows + numberColumns, status_); |
---|
3131 | } |
---|
3132 | // Do local minimization |
---|
3133 | #define LOCAL |
---|
3134 | #ifdef LOCAL |
---|
3135 | bool absolutelyOptimal = false; |
---|
3136 | int saveScaling = scalingFlag(); |
---|
3137 | scaling(0); |
---|
3138 | int savePerturbation = perturbation_; |
---|
3139 | perturbation_ = 100; |
---|
3140 | if (saveLogLevel == 1) |
---|
3141 | setLogLevel(0); |
---|
3142 | int status = startup(1); |
---|
3143 | if (!status) { |
---|
3144 | int numberTotal = numberRows_ + numberColumns_; |
---|
3145 | // resize arrays |
---|
3146 | for (int i = 0; i < 4; i++) { |
---|
3147 | rowArray_[i]->reserve(CoinMax(numberRows_ + numberColumns_, rowArray_[i]->capacity())); |
---|
3148 | } |
---|
3149 | CoinIndexedVector *longArray = rowArray_[3]; |
---|
3150 | CoinIndexedVector *rowArray = rowArray_[0]; |
---|
3151 | //CoinIndexedVector * columnArray = columnArray_[0]; |
---|
3152 | CoinIndexedVector *spare = rowArray_[1]; |
---|
3153 | double *work = longArray->denseVector(); |
---|
3154 | int *which = longArray->getIndices(); |
---|
3155 | int nPass = 100; |
---|
3156 | //bool conjugate=false; |
---|
3157 | // Put back true bounds |
---|
3158 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
3159 | int iColumn = listNonLinearColumn[jNon]; |
---|
3160 | double value; |
---|
3161 | value = trueLower[jNon]; |
---|
3162 | trueLower[jNon] = lower_[iColumn]; |
---|
3163 | lower_[iColumn] = value; |
---|
3164 | value = trueUpper[jNon]; |
---|
3165 | trueUpper[jNon] = upper_[iColumn]; |
---|
3166 | upper_[iColumn] = value; |
---|
3167 | switch (getStatus(iColumn)) { |
---|
3168 | |
---|
3169 | case basic: |
---|
3170 | case isFree: |
---|
3171 | case superBasic: |
---|
3172 | break; |
---|
3173 | case isFixed: |
---|
3174 | case atUpperBound: |
---|
3175 | case atLowerBound: |
---|
3176 | if (solution_[iColumn] > lower_[iColumn] + primalTolerance_) { |
---|
3177 | if (solution_[iColumn] < upper_[iColumn] - primalTolerance_) { |
---|
3178 | setStatus(iColumn, superBasic); |
---|
3179 | } else { |
---|
3180 | setStatus(iColumn, atUpperBound); |
---|
3181 | } |
---|
3182 | } else { |
---|
3183 | if (solution_[iColumn] < upper_[iColumn] - primalTolerance_) { |
---|
3184 | setStatus(iColumn, atLowerBound); |
---|
3185 | } else { |
---|
3186 | setStatus(iColumn, isFixed); |
---|
3187 | } |
---|
3188 | } |
---|
3189 | break; |
---|
3190 | } |
---|
3191 | } |
---|
3192 | for (int jPass = 0; jPass < nPass; jPass++) { |
---|
3193 | trueObjective->reducedGradient(this, dj_, true); |
---|
3194 | // get direction vector in longArray |
---|
3195 | longArray->clear(); |
---|
3196 | double normFlagged = 0.0; |
---|
3197 | double normUnflagged = 0.0; |
---|
3198 | int numberNonBasic = 0; |
---|
3199 | directionVector(longArray, spare, rowArray, 0, |
---|
3200 | normFlagged, normUnflagged, numberNonBasic); |
---|
3201 | if (normFlagged + normUnflagged < 1.0e-8) { |
---|
3202 | absolutelyOptimal = true; |
---|
3203 | break; //looks optimal |
---|
3204 | } |
---|
3205 | double djNorm = 0.0; |
---|
3206 | int iIndex; |
---|
3207 | for (iIndex = 0; iIndex < numberNonBasic; iIndex++) { |
---|
3208 | int iSequence = which[iIndex]; |
---|
3209 | double alpha = work[iSequence]; |
---|
3210 | djNorm += alpha * alpha; |
---|
3211 | } |
---|
3212 | //if (!jPass) |
---|
3213 | //printf("dj norm %g - %d \n",djNorm,numberNonBasic); |
---|
3214 | //int number=longArray->getNumElements(); |
---|
3215 | if (!numberNonBasic) { |
---|
3216 | // looks optimal |
---|
3217 | absolutelyOptimal = true; |
---|
3218 | break; |
---|
3219 | } |
---|
3220 | double theta = 1.0e30; |
---|
3221 | int iSequence; |
---|
3222 | for (iSequence = 0; iSequence < numberTotal; iSequence++) { |
---|
3223 | double alpha = work[iSequence]; |
---|
3224 | double oldValue = solution_[iSequence]; |
---|
3225 | if (alpha < -1.0e-15) { |
---|
3226 | // variable going towards lower bound |
---|
3227 | double bound = lower_[iSequence]; |
---|
3228 | oldValue -= bound; |
---|
3229 | if (oldValue + theta * alpha < 0.0) { |
---|
3230 | theta = CoinMax(0.0, oldValue / (-alpha)); |
---|
3231 | } |
---|
3232 | } else if (alpha > 1.0e-15) { |
---|
3233 | // variable going towards upper bound |
---|
3234 | double bound = upper_[iSequence]; |
---|
3235 | oldValue = bound - oldValue; |
---|
3236 | if (oldValue - theta * alpha < 0.0) { |
---|
3237 | theta = CoinMax(0.0, oldValue / alpha); |
---|
3238 | } |
---|
3239 | } |
---|
3240 | } |
---|
3241 | // Now find minimum of function |
---|
3242 | double currentObj; |
---|
3243 | double predictedObj; |
---|
3244 | double thetaObj; |
---|
3245 | // need to use true objective |
---|
3246 | double *saveCost = cost_; |
---|
3247 | cost_ = NULL; |
---|
3248 | double objTheta = trueObjective->stepLength(this, solution_, work, theta, |
---|
3249 | currentObj, predictedObj, thetaObj); |
---|
3250 | cost_ = saveCost; |
---|
3251 | #ifdef CLP_DEBUG |
---|
3252 | if (handler_->logLevel() & 32) |
---|
3253 | printf("current obj %g thetaObj %g, predictedObj %g\n", currentObj, thetaObj, predictedObj); |
---|
3254 | #endif |
---|
3255 | //printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj); |
---|
3256 | //printf("objTheta %g theta %g\n",objTheta,theta); |
---|
3257 | if (theta > objTheta) { |
---|
3258 | theta = objTheta; |
---|
3259 | thetaObj = predictedObj; |
---|
3260 | } |
---|
3261 | // update one used outside |
---|
3262 | objValue = currentObj; |
---|
3263 | if (theta > 1.0e-9 && (currentObj - thetaObj < -CoinMax(1.0e-8, 1.0e-15 * fabs(currentObj)) || jPass < 5)) { |
---|
3264 | // Update solution |
---|
3265 | for (iSequence = 0; iSequence < numberTotal; iSequence++) { |
---|
3266 | double alpha = work[iSequence]; |
---|
3267 | if (alpha) { |
---|
3268 | double value = solution_[iSequence] + theta * alpha; |
---|
3269 | solution_[iSequence] = value; |
---|
3270 | switch (getStatus(iSequence)) { |
---|
3271 | |
---|
3272 | case basic: |
---|
3273 | case isFixed: |
---|
3274 | case isFree: |
---|
3275 | break; |
---|
3276 | case atUpperBound: |
---|
3277 | case atLowerBound: |
---|
3278 | case superBasic: |
---|
3279 | if (value > lower_[iSequence] + primalTolerance_) { |
---|
3280 | if (value < upper_[iSequence] - primalTolerance_) { |
---|
3281 | setStatus(iSequence, superBasic); |
---|
3282 | } else { |
---|
3283 | setStatus(iSequence, atUpperBound); |
---|
3284 | } |
---|
3285 | } else { |
---|
3286 | if (value < upper_[iSequence] - primalTolerance_) { |
---|
3287 | setStatus(iSequence, atLowerBound); |
---|
3288 | } else { |
---|
3289 | setStatus(iSequence, isFixed); |
---|
3290 | } |
---|
3291 | } |
---|
3292 | break; |
---|
3293 | } |
---|
3294 | } |
---|
3295 | } |
---|
3296 | } else { |
---|
3297 | break; |
---|
3298 | } |
---|
3299 | } |
---|
3300 | // Put back fake bounds |
---|
3301 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
3302 | int iColumn = listNonLinearColumn[jNon]; |
---|
3303 | double value; |
---|
3304 | value = trueLower[jNon]; |
---|
3305 | trueLower[jNon] = lower_[iColumn]; |
---|
3306 | lower_[iColumn] = value; |
---|
3307 | value = trueUpper[jNon]; |
---|
3308 | trueUpper[jNon] = upper_[iColumn]; |
---|
3309 | upper_[iColumn] = value; |
---|
3310 | } |
---|
3311 | } |
---|
3312 | problemStatus_ = 0; |
---|
3313 | finish(); |
---|
3314 | scaling(saveScaling); |
---|
3315 | perturbation_ = savePerturbation; |
---|
3316 | setLogLevel(saveLogLevel); |
---|
3317 | #endif |
---|
3318 | // redo rowActivity |
---|
3319 | memset(rowActivity_, 0, numberRows_ * sizeof(double)); |
---|
3320 | times(1.0, solution, rowActivity_); |
---|
3321 | if (theta > 0.99999 && theta2 < 1.9 && !absolutelyOptimal) { |
---|
3322 | // If big changes then tighten |
---|
3323 | /* thetaObj is objvalue + a*2*2 +b*2 |
---|
3324 | predictedObj is objvalue + a*theta2*theta2 +b*theta2 |
---|
3325 | */ |
---|
3326 | double rhs1 = thetaObj - objValue; |
---|
3327 | double rhs2 = predictedObj - objValue; |
---|
3328 | double subtractB = theta2 * 0.5; |
---|
3329 | double a = (rhs2 - subtractB * rhs1) / (theta2 * theta2 - 4.0 * subtractB); |
---|
3330 | double b = 0.5 * (rhs1 - 4.0 * a); |
---|
3331 | if (fabs(a + b) > 1.0e-2) { |
---|
3332 | // tighten all |
---|
3333 | goodMove = -1; |
---|
3334 | } |
---|
3335 | } |
---|
3336 | } |
---|
3337 | } |
---|
3338 | CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, 2), numberColumns, |
---|
3339 | objective); |
---|
3340 | //printf("offset comp %g orig %g\n",offset,objectiveOffset); |
---|
3341 | // could do faster |
---|
3342 | trueObjective->stepLength(this, solution, changeRegion, 0.0, |
---|
3343 | objValue, predictedObj, thetaObj); |
---|
3344 | #ifdef CLP_INVESTIGATE |
---|
3345 | printf("offset comp %g orig %g - obj from stepLength %g\n", offset, objectiveOffset, objValue); |
---|
3346 | #endif |
---|
3347 | setDblParam(ClpObjOffset, objectiveOffset + offset); |
---|
3348 | int *temp = last[2]; |
---|
3349 | last[2] = last[1]; |
---|
3350 | last[1] = last[0]; |
---|
3351 | last[0] = temp; |
---|
3352 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
3353 | iColumn = listNonLinearColumn[jNon]; |
---|
3354 | double change = solution[iColumn] - saveSolution[iColumn]; |
---|
3355 | if (change < -1.0e-5) { |
---|
3356 | if (fabs(change + trust[jNon]) < 1.0e-5) |
---|
3357 | temp[jNon] = -1; |
---|
3358 | else |
---|
3359 | temp[jNon] = -2; |
---|
3360 | } else if (change > 1.0e-5) { |
---|
3361 | if (fabs(change - trust[jNon]) < 1.0e-5) |
---|
3362 | temp[jNon] = 1; |
---|
3363 | else |
---|
3364 | temp[jNon] = 2; |
---|
3365 | } else { |
---|
3366 | temp[jNon] = 0; |
---|
3367 | } |
---|
3368 | } |
---|
3369 | // goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing |
---|
3370 | double maxDelta = 0.0; |
---|
3371 | if (goodMove >= 0) { |
---|
3372 | if (objValue - lastObjective <= 1.0e-15 * fabs(lastObjective)) |
---|
3373 | goodMove = 1; |
---|
3374 | else |
---|
3375 | goodMove = 0; |
---|
3376 | } else { |
---|
3377 | maxDelta = 1.0e10; |
---|
3378 | } |
---|
3379 | double maxGap = 0.0; |
---|
3380 | int numberSmaller = 0; |
---|
3381 | int numberSmaller2 = 0; |
---|
3382 | int numberLarger = 0; |
---|
3383 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
3384 | iColumn = listNonLinearColumn[jNon]; |
---|
3385 | maxDelta = CoinMax(maxDelta, |
---|
3386 | fabs(solution[iColumn] - saveSolution[iColumn])); |
---|
3387 | if (goodMove > 0) { |
---|
3388 | if (last[0][jNon] * last[1][jNon] < 0) { |
---|
3389 | // halve |
---|
3390 | trust[jNon] *= 0.5; |
---|
3391 | numberSmaller2++; |
---|
3392 | } else { |
---|
3393 | if (last[0][jNon] == last[1][jNon] && last[0][jNon] == last[2][jNon]) |
---|
3394 | trust[jNon] = CoinMin(1.5 * trust[jNon], 1.0e6); |
---|
3395 | numberLarger++; |
---|
3396 | } |
---|
3397 | } else if (goodMove != -2 && trust[jNon] > 10.0 * deltaTolerance) { |
---|
3398 | trust[jNon] *= 0.2; |
---|
3399 | numberSmaller++; |
---|
3400 | } |
---|
3401 | maxGap = CoinMax(maxGap, trust[jNon]); |
---|
3402 | } |
---|
3403 | #ifdef CLP_DEBUG |
---|
3404 | if (handler_->logLevel() & 32) |
---|
3405 | std::cout << "largest gap is " << maxGap << " " |
---|
3406 | << numberSmaller + numberSmaller2 << " reduced (" |
---|
3407 | << numberSmaller << " badMove ), " |
---|
3408 | << numberLarger << " increased" << std::endl; |
---|
3409 | #endif |
---|
3410 | if (iPass > 10000) { |
---|
3411 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) |
---|
3412 | trust[jNon] *= 0.0001; |
---|
3413 | } |
---|
3414 | if (lastMoveStatus == -1 && goodMove == -1) |
---|
3415 | goodMove = 1; // to force solve |
---|
3416 | if (goodMove > 0) { |
---|
3417 | double drop = lastObjective - objValue; |
---|
3418 | handler_->message(CLP_SLP_ITER, messages_) |
---|
3419 | << iPass << objValue - objectiveOffset |
---|
3420 | << drop << maxDelta |
---|
3421 | << CoinMessageEol; |
---|
3422 | if (iPass > 20 && drop < 1.0e-12 * fabs(objValue)) |
---|
3423 | drop = 0.999e-4; // so will exit |
---|
3424 | if (maxDelta < deltaTolerance && drop < 1.0e-4 && goodMove && theta < 0.99999) { |
---|
3425 | if (handler_->logLevel() > 1) |
---|
3426 | std::cout << "Exiting as maxDelta < tolerance and small drop" << std::endl; |
---|
3427 | break; |
---|
3428 | } |
---|
3429 | } else if (!numberSmaller && iPass > 1) { |
---|
3430 | if (handler_->logLevel() > 1) |
---|
3431 | std::cout << "Exiting as all gaps small" << std::endl; |
---|
3432 | break; |
---|
3433 | } |
---|
3434 | if (!iPass) |
---|
3435 | goodMove = 1; |
---|
3436 | targetDrop = 0.0; |
---|
3437 | double *r = this->dualColumnSolution(); |
---|
3438 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
3439 | iColumn = listNonLinearColumn[jNon]; |
---|
3440 | columnLower[iColumn] = CoinMax(solution[iColumn] |
---|
3441 | - trust[jNon], |
---|
3442 | trueLower[jNon]); |
---|
3443 | columnUpper[iColumn] = CoinMin(solution[iColumn] |
---|
3444 | + trust[jNon], |
---|
3445 | trueUpper[jNon]); |
---|
3446 | } |
---|
3447 | if (iPass) { |
---|
3448 | // get reduced costs |
---|
3449 | this->matrix()->transposeTimes(savePi, |
---|
3450 | this->dualColumnSolution()); |
---|
3451 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
3452 | iColumn = listNonLinearColumn[jNon]; |
---|
3453 | double dj = objective[iColumn] - r[iColumn]; |
---|
3454 | r[iColumn] = dj; |
---|
3455 | if (dj < -dualTolerance_) |
---|
3456 | targetDrop -= dj * (columnUpper[iColumn] - solution[iColumn]); |
---|
3457 | else if (dj > dualTolerance_) |
---|
3458 | targetDrop -= dj * (columnLower[iColumn] - solution[iColumn]); |
---|
3459 | } |
---|
3460 | } else { |
---|
3461 | memset(r, 0, numberColumns * sizeof(double)); |
---|
3462 | } |
---|
3463 | #if 0 |
---|
3464 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
3465 | iColumn = listNonLinearColumn[jNon]; |
---|
3466 | if (statusCheck[iColumn] == 'L' && r[iColumn] < -1.0e-4) { |
---|
3467 | columnLower[iColumn] = CoinMax(solution[iColumn], |
---|
3468 | trueLower[jNon]); |
---|
3469 | columnUpper[iColumn] = CoinMin(solution[iColumn] |
---|
3470 | + trust[jNon], |
---|
3471 | trueUpper[jNon]); |
---|
3472 | } else if (statusCheck[iColumn] == 'U' && r[iColumn] > 1.0e-4) { |
---|
3473 | columnLower[iColumn] = CoinMax(solution[iColumn] |
---|
3474 | - trust[jNon], |
---|
3475 | trueLower[jNon]); |
---|
3476 | columnUpper[iColumn] = CoinMin(solution[iColumn], |
---|
3477 | trueUpper[jNon]); |
---|
3478 | } else { |
---|
3479 | columnLower[iColumn] = CoinMax(solution[iColumn] |
---|
3480 | - trust[jNon], |
---|
3481 | trueLower[jNon]); |
---|
3482 | columnUpper[iColumn] = CoinMin(solution[iColumn] |
---|
3483 | + trust[jNon], |
---|
3484 | trueUpper[jNon]); |
---|
3485 | } |
---|
3486 | } |
---|
3487 | #endif |
---|
3488 | if (goodMove > 0) { |
---|
3489 | CoinMemcpyN(solution, numberColumns, saveSolution); |
---|
3490 | CoinMemcpyN(rowActivity_, numberRows, saveRowSolution); |
---|
3491 | CoinMemcpyN(this->dualRowSolution(), numberRows, savePi); |
---|
3492 | CoinMemcpyN(status_, numberRows + numberColumns, saveStatus); |
---|
3493 | #if MULTIPLE > 2 |
---|
3494 | double *tempSol = saveSolutionM[0]; |
---|
3495 | for (jNon = 0; jNon < MULTIPLE - 1; jNon++) { |
---|
3496 | saveSolutionM[jNon] = saveSolutionM[jNon + 1]; |
---|
3497 | } |
---|
3498 | saveSolutionM[MULTIPLE - 1] = tempSol; |
---|
3499 | CoinMemcpyN(solution, numberColumns, tempSol); |
---|
3500 | |
---|
3501 | #endif |
---|
3502 | |
---|
3503 | #ifdef CLP_DEBUG |
---|
3504 | if (handler_->logLevel() & 32) |
---|
3505 | std::cout << "Pass - " << iPass |
---|
3506 | << ", target drop is " << targetDrop |
---|
3507 | << std::endl; |
---|
3508 | #endif |
---|
3509 | lastObjective = objValue; |
---|
3510 | if (targetDrop < CoinMax(1.0e-8, CoinMin(1.0e-6, 1.0e-6 * fabs(objValue))) && goodMove && iPass > 3) { |
---|
3511 | if (handler_->logLevel() > 1) |
---|
3512 | printf("Exiting on target drop %g\n", targetDrop); |
---|
3513 | break; |
---|
3514 | } |
---|
3515 | #ifdef CLP_DEBUG |
---|
3516 | { |
---|
3517 | double *r = this->dualColumnSolution(); |
---|
3518 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
3519 | iColumn = listNonLinearColumn[jNon]; |
---|
3520 | if (handler_->logLevel() & 32) |
---|
3521 | printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n", |
---|
3522 | jNon, trust[jNon], iColumn, solution[iColumn], objective[iColumn], |
---|
3523 | r[iColumn], statusCheck[iColumn], columnLower[iColumn], |
---|
3524 | columnUpper[iColumn]); |
---|
3525 | } |
---|
3526 | } |
---|
3527 | #endif |
---|
3528 | //setLogLevel(63); |
---|
3529 | this->scaling(false); |
---|
3530 | if (saveLogLevel == 1) |
---|
3531 | setLogLevel(0); |
---|
3532 | primalDualCuts(rowsIn, 1, 1); |
---|
3533 | algorithm_ = 1; |
---|
3534 | setLogLevel(saveLogLevel); |
---|
3535 | #ifdef CLP_DEBUG |
---|
3536 | if (this->status()) { |
---|
3537 | writeMps("xx.mps"); |
---|
3538 | } |
---|
3539 | #endif |
---|
3540 | if (this->status() == 1) { |
---|
3541 | // not feasible ! - backtrack and exit |
---|
3542 | // use safe solution |
---|
3543 | CoinMemcpyN(safeSolution, numberColumns, solution); |
---|
3544 | CoinMemcpyN(solution, numberColumns, saveSolution); |
---|
3545 | memset(rowActivity_, 0, numberRows_ * sizeof(double)); |
---|
3546 | times(1.0, solution, rowActivity_); |
---|
3547 | CoinMemcpyN(rowActivity_, numberRows, saveRowSolution); |
---|
3548 | CoinMemcpyN(savePi, numberRows, this->dualRowSolution()); |
---|
3549 | CoinMemcpyN(saveStatus, numberRows + numberColumns, status_); |
---|
3550 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
3551 | iColumn = listNonLinearColumn[jNon]; |
---|
3552 | columnLower[iColumn] = CoinMax(solution[iColumn] |
---|
3553 | - trust[jNon], |
---|
3554 | trueLower[jNon]); |
---|
3555 | columnUpper[iColumn] = CoinMin(solution[iColumn] |
---|
3556 | + trust[jNon], |
---|
3557 | trueUpper[jNon]); |
---|
3558 | } |
---|
3559 | break; |
---|
3560 | } else { |
---|
3561 | // save in case problems |
---|
3562 | CoinMemcpyN(solution, numberColumns, safeSolution); |
---|
3563 | } |
---|
3564 | if (problemStatus_ == 3) |
---|
3565 | break; // probably user interrupt |
---|
3566 | goodMove = 1; |
---|
3567 | } else { |
---|
3568 | // bad pass - restore solution |
---|
3569 | #ifdef CLP_DEBUG |
---|
3570 | if (handler_->logLevel() & 32) |
---|
3571 | printf("Backtracking\n"); |
---|
3572 | #endif |
---|
3573 | CoinMemcpyN(saveSolution, numberColumns, solution); |
---|
3574 | CoinMemcpyN(saveRowSolution, numberRows, rowActivity_); |
---|
3575 | CoinMemcpyN(savePi, numberRows, this->dualRowSolution()); |
---|
3576 | CoinMemcpyN(saveStatus, numberRows + numberColumns, status_); |
---|
3577 | iPass--; |
---|
3578 | assert(exitPass > 0); |
---|
3579 | goodMove = -1; |
---|
3580 | } |
---|
3581 | } |
---|
3582 | #if MULTIPLE > 2 |
---|
3583 | for (jNon = 0; jNon < MULTIPLE; jNon++) { |
---|
3584 | delete[] saveSolutionM[jNon]; |
---|
3585 | } |
---|
3586 | #endif |
---|
3587 | // restore solution |
---|
3588 | CoinMemcpyN(saveSolution, numberColumns, solution); |
---|
3589 | CoinMemcpyN(saveRowSolution, numberRows, rowActivity_); |
---|
3590 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
3591 | iColumn = listNonLinearColumn[jNon]; |
---|
3592 | columnLower[iColumn] = CoinMax(solution[iColumn], |
---|
3593 | trueLower[jNon]); |
---|
3594 | columnUpper[iColumn] = CoinMin(solution[iColumn], |
---|
3595 | trueUpper[jNon]); |
---|
3596 | } |
---|
3597 | delete[] markNonlinear; |
---|
3598 | delete[] statusCheck; |
---|
3599 | delete[] savePi; |
---|
3600 | delete[] saveStatus; |
---|
3601 | // redo objective |
---|
3602 | CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, 2), numberColumns, |
---|
3603 | objective); |
---|
3604 | primalDualCuts(rowsIn, 1, 1); |
---|
3605 | delete objective_; |
---|
3606 | delete[] rowsIn; |
---|
3607 | objective_ = trueObjective; |
---|
3608 | // redo values |
---|
3609 | setDblParam(ClpObjOffset, objectiveOffset); |
---|
3610 | objectiveValue_ += whichWay * offset; |
---|
3611 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
3612 | iColumn = listNonLinearColumn[jNon]; |
---|
3613 | columnLower[iColumn] = trueLower[jNon]; |
---|
3614 | columnUpper[iColumn] = trueUpper[jNon]; |
---|
3615 | } |
---|
3616 | delete[] saveSolution; |
---|
3617 | delete[] safeSolution; |
---|
3618 | delete[] saveRowSolution; |
---|
3619 | for (iPass = 0; iPass < 3; iPass++) |
---|
3620 | delete[] last[iPass]; |
---|
3621 | delete[] trust; |
---|
3622 | delete[] trueUpper; |
---|
3623 | delete[] trueLower; |
---|
3624 | delete[] listNonLinearColumn; |
---|
3625 | delete[] changeRegion; |
---|
3626 | // temp |
---|
3627 | //setLogLevel(63); |
---|
3628 | return 0; |
---|
3629 | } |
---|
3630 | /* Primal algorithm for nonlinear constraints |
---|
3631 | Using a semi-trust region approach as for pooling problem |
---|
3632 | This is in because I have it lying around |
---|
3633 | |
---|
3634 | */ |
---|
3635 | int ClpSimplexNonlinear::primalSLP(int numberConstraints, ClpConstraint **constraints, |
---|
3636 | int numberPasses, double deltaTolerance) |
---|
3637 | { |
---|
3638 | if (!numberConstraints) { |
---|
3639 | // no nonlinear constraints - may be nonlinear objective |
---|
3640 | return primalSLP(numberPasses, deltaTolerance); |
---|
3641 | } |
---|
3642 | // Are we minimizing or maximizing |
---|
3643 | double whichWay = optimizationDirection(); |
---|
3644 | if (whichWay < 0.0) |
---|
3645 | whichWay = -1.0; |
---|
3646 | else if (whichWay > 0.0) |
---|
3647 | whichWay = 1.0; |
---|
3648 | // check all matrix for odd rows is empty |
---|
3649 | int iConstraint; |
---|
3650 | int numberBad = 0; |
---|
3651 | CoinPackedMatrix *columnCopy = matrix(); |
---|
3652 | // Get a row copy in standard format |
---|
3653 | CoinPackedMatrix copy; |
---|
3654 | copy.reverseOrderedCopyOf(*columnCopy); |
---|
3655 | // get matrix data pointers |
---|
3656 | //const int * column = copy.getIndices(); |
---|
3657 | //const CoinBigIndex * rowStart = copy.getVectorStarts(); |
---|
3658 | const int *rowLength = copy.getVectorLengths(); |
---|
3659 | //const double * elementByRow = copy.getElements(); |
---|
3660 | int numberArtificials = 0; |
---|
3661 | // We could use nonlinearcost to do segments - maybe later |
---|
3662 | #define SEGMENTS 3 |
---|
3663 | // Penalties may be adjusted by duals |
---|
3664 | // Both these should be modified depending on problem |
---|
3665 | // Possibly start with big bounds |
---|
3666 | //double penalties[]={1.0e-3,1.0e7,1.0e9}; |
---|
3667 | double penalties[] = { 1.0e7, 1.0e8, 1.0e9 }; |
---|
3668 | double bounds[] = { 1.0e-2, 1.0e2, COIN_DBL_MAX }; |
---|
3669 | // see how many extra we need |
---|
3670 | CoinBigIndex numberExtra = 0; |
---|
3671 | for (iConstraint = 0; iConstraint < numberConstraints; iConstraint++) { |
---|
3672 | ClpConstraint *constraint = constraints[iConstraint]; |
---|
3673 | int iRow = constraint->rowNumber(); |
---|
3674 | assert(iRow >= 0); |
---|
3675 | int n = constraint->numberCoefficients() - rowLength[iRow]; |
---|
3676 | numberExtra += n; |
---|
3677 | if (iRow >= numberRows_) |
---|
3678 | numberBad++; |
---|
3679 | else if (rowLength[iRow] && n) |
---|
3680 | numberBad++; |
---|
3681 | if (rowLower_[iRow] > -1.0e20) |
---|
3682 | numberArtificials++; |
---|
3683 | if (rowUpper_[iRow] < 1.0e20) |
---|
3684 | numberArtificials++; |
---|
3685 | } |
---|
3686 | if (numberBad) |
---|
3687 | return numberBad; |
---|
3688 | ClpObjective *trueObjective = NULL; |
---|
3689 | if (objective_->type() >= 2) { |
---|
3690 | // Replace objective |
---|
3691 | trueObjective = objective_; |
---|
3692 | objective_ = new ClpLinearObjective(NULL, numberColumns_); |
---|
3693 | } |
---|
3694 | ClpSimplex newModel(*this); |
---|
3695 | // we can put back true objective |
---|
3696 | if (trueObjective) { |
---|
3697 | // Replace objective |
---|
3698 | delete objective_; |
---|
3699 | objective_ = trueObjective; |
---|
3700 | } |
---|
3701 | int numberColumns2 = numberColumns_; |
---|
3702 | int numberSmallGap = numberArtificials; |
---|
3703 | if (numberArtificials) { |
---|
3704 | numberArtificials *= SEGMENTS; |
---|
3705 | numberColumns2 += numberArtificials; |
---|
3706 | CoinBigIndex *addStarts = new CoinBigIndex[numberArtificials + 1]; |
---|
3707 | int *addRow = new int[numberArtificials]; |
---|
3708 | double *addElement = new double[numberArtificials]; |
---|
3709 | double *addUpper = new double[numberArtificials]; |
---|
3710 | addStarts[0] = 0; |
---|
3711 | double *addCost = new double[numberArtificials]; |
---|
3712 | numberArtificials = 0; |
---|
3713 | for (iConstraint = 0; iConstraint < numberConstraints; iConstraint++) { |
---|
3714 | ClpConstraint *constraint = constraints[iConstraint]; |
---|
3715 | int iRow = constraint->rowNumber(); |
---|
3716 | if (rowLower_[iRow] > -1.0e20) { |
---|
3717 | for (int k = 0; k < SEGMENTS; k++) { |
---|
3718 | addRow[numberArtificials] = iRow; |
---|
3719 | addElement[numberArtificials] = 1.0; |
---|
3720 | addCost[numberArtificials] = penalties[k]; |
---|
3721 | addUpper[numberArtificials] = bounds[k]; |
---|
3722 | numberArtificials++; |
---|
3723 | addStarts[numberArtificials] = numberArtificials; |
---|
3724 | } |
---|
3725 | } |
---|
3726 | if (rowUpper_[iRow] < 1.0e20) { |
---|
3727 | for (int k = 0; k < SEGMENTS; k++) { |
---|
3728 | addRow[numberArtificials] = iRow; |
---|
3729 | addElement[numberArtificials] = -1.0; |
---|
3730 | addCost[numberArtificials] = penalties[k]; |
---|
3731 | addUpper[numberArtificials] = bounds[k]; |
---|
3732 | numberArtificials++; |
---|
3733 | addStarts[numberArtificials] = numberArtificials; |
---|
3734 | } |
---|
3735 | } |
---|
3736 | } |
---|
3737 | newModel.addColumns(numberArtificials, NULL, addUpper, addCost, |
---|
3738 | addStarts, addRow, addElement); |
---|
3739 | delete[] addStarts; |
---|
3740 | delete[] addRow; |
---|
3741 | delete[] addElement; |
---|
3742 | delete[] addUpper; |
---|
3743 | delete[] addCost; |
---|
3744 | // newModel.primal(1); |
---|
3745 | } |
---|
3746 | // find nonlinear columns |
---|
3747 | int *listNonLinearColumn = new int[numberColumns_ + numberSmallGap]; |
---|
3748 | char *mark = new char[numberColumns_]; |
---|
3749 | memset(mark, 0, numberColumns_); |
---|
3750 | for (iConstraint = 0; iConstraint < numberConstraints; iConstraint++) { |
---|
3751 | ClpConstraint *constraint = constraints[iConstraint]; |
---|
3752 | constraint->markNonlinear(mark); |
---|
3753 | } |
---|
3754 | if (trueObjective) |
---|
3755 | trueObjective->markNonlinear(mark); |
---|
3756 | int iColumn; |
---|
3757 | int numberNonLinearColumns = 0; |
---|
3758 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
3759 | if (mark[iColumn]) |
---|
3760 | listNonLinearColumn[numberNonLinearColumns++] = iColumn; |
---|
3761 | } |
---|
3762 | // and small gap artificials |
---|
3763 | for (iColumn = numberColumns_; iColumn < numberColumns2; iColumn += SEGMENTS) { |
---|
3764 | listNonLinearColumn[numberNonLinearColumns++] = iColumn; |
---|
3765 | } |
---|
3766 | // build row copy of matrix with space for nonzeros |
---|
3767 | // Get column copy |
---|
3768 | columnCopy = newModel.matrix(); |
---|
3769 | copy.reverseOrderedCopyOf(*columnCopy); |
---|
3770 | // get matrix data pointers |
---|
3771 | const int *column = copy.getIndices(); |
---|
3772 | const CoinBigIndex *rowStart = copy.getVectorStarts(); |
---|
3773 | rowLength = copy.getVectorLengths(); |
---|
3774 | const double *elementByRow = copy.getElements(); |
---|
3775 | numberExtra += copy.getNumElements(); |
---|
3776 | CoinBigIndex *newStarts = new CoinBigIndex[numberRows_ + 1]; |
---|
3777 | int *newColumn = new int[numberExtra]; |
---|
3778 | double *newElement = new double[numberExtra]; |
---|
3779 | newStarts[0] = 0; |
---|
3780 | int *backRow = new int[numberRows_]; |
---|
3781 | int iRow; |
---|
3782 | for (iRow = 0; iRow < numberRows_; iRow++) |
---|
3783 | backRow[iRow] = -1; |
---|
3784 | for (iConstraint = 0; iConstraint < numberConstraints; iConstraint++) { |
---|
3785 | ClpConstraint *constraint = constraints[iConstraint]; |
---|
3786 | iRow = constraint->rowNumber(); |
---|
3787 | backRow[iRow] = iConstraint; |
---|
3788 | } |
---|
3789 | numberExtra = 0; |
---|
3790 | for (iRow = 0; iRow < numberRows_; iRow++) { |
---|
3791 | if (backRow[iRow] < 0) { |
---|
3792 | // copy normal |
---|
3793 | for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; |
---|
3794 | j++) { |
---|
3795 | newColumn[numberExtra] = column[j]; |
---|
3796 | newElement[numberExtra++] = elementByRow[j]; |
---|
3797 | } |
---|
3798 | } else { |
---|
3799 | ClpConstraint *constraint = constraints[backRow[iRow]]; |
---|
3800 | assert(iRow == constraint->rowNumber()); |
---|
3801 | int numberArtificials = 0; |
---|
3802 | if (rowLower_[iRow] > -1.0e20) |
---|
3803 | numberArtificials += SEGMENTS; |
---|
3804 | if (rowUpper_[iRow] < 1.0e20) |
---|
3805 | numberArtificials += SEGMENTS; |
---|
3806 | if (numberArtificials == rowLength[iRow]) { |
---|
3807 | // all possible |
---|
3808 | memset(mark, 0, numberColumns_); |
---|
3809 | constraint->markNonzero(mark); |
---|
3810 | for (int k = 0; k < numberColumns_; k++) { |
---|
3811 | if (mark[k]) { |
---|
3812 | newColumn[numberExtra] = k; |
---|
3813 | newElement[numberExtra++] = 1.0; |
---|
3814 | } |
---|
3815 | } |
---|
3816 | // copy artificials |
---|
3817 | for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; |
---|
3818 | j++) { |
---|
3819 | newColumn[numberExtra] = column[j]; |
---|
3820 | newElement[numberExtra++] = elementByRow[j]; |
---|
3821 | } |
---|
3822 | } else { |
---|
3823 | // there already |
---|
3824 | // copy |
---|
3825 | for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; |
---|
3826 | j++) { |
---|
3827 | newColumn[numberExtra] = column[j]; |
---|
3828 | assert(elementByRow[j]); |
---|
3829 | newElement[numberExtra++] = elementByRow[j]; |
---|
3830 | } |
---|
3831 | } |
---|
3832 | } |
---|
3833 | newStarts[iRow + 1] = numberExtra; |
---|
3834 | } |
---|
3835 | delete[] backRow; |
---|
3836 | CoinPackedMatrix saveMatrix(false, numberColumns2, numberRows_, |
---|
3837 | numberExtra, newElement, newColumn, newStarts, NULL, 0.0, 0.0); |
---|
3838 | delete[] newStarts; |
---|
3839 | delete[] newColumn; |
---|
3840 | delete[] newElement; |
---|
3841 | delete[] mark; |
---|
3842 | // get feasible |
---|
3843 | if (whichWay < 0.0) { |
---|
3844 | newModel.setOptimizationDirection(1.0); |
---|
3845 | double *objective = newModel.objective(); |
---|
3846 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) |
---|
3847 | objective[iColumn] = -objective[iColumn]; |
---|
3848 | } |
---|
3849 | newModel.primal(1); |
---|
3850 | // still infeasible |
---|
3851 | if (newModel.problemStatus() == 1) { |
---|
3852 | delete[] listNonLinearColumn; |
---|
3853 | return 0; |
---|
3854 | } else if (newModel.problemStatus() == 2) { |
---|
3855 | // unbounded - add bounds |
---|
3856 | double *columnLower = newModel.columnLower(); |
---|
3857 | double *columnUpper = newModel.columnUpper(); |
---|
3858 | for (int i = 0; i < numberColumns_; i++) { |
---|
3859 | columnLower[i] = CoinMax(-1.0e8, columnLower[i]); |
---|
3860 | columnUpper[i] = CoinMin(1.0e8, columnUpper[i]); |
---|
3861 | } |
---|
3862 | newModel.primal(1); |
---|
3863 | } |
---|
3864 | int numberRows = newModel.numberRows(); |
---|
3865 | double *columnLower = newModel.columnLower(); |
---|
3866 | double *columnUpper = newModel.columnUpper(); |
---|
3867 | double *objective = newModel.objective(); |
---|
3868 | double *rowLower = newModel.rowLower(); |
---|
3869 | double *rowUpper = newModel.rowUpper(); |
---|
3870 | double *solution = newModel.primalColumnSolution(); |
---|
3871 | int jNon; |
---|
3872 | int *last[3]; |
---|
3873 | |
---|
3874 | double *trust = new double[numberNonLinearColumns]; |
---|
3875 | double *trueLower = new double[numberNonLinearColumns]; |
---|
3876 | double *trueUpper = new double[numberNonLinearColumns]; |
---|
3877 | double objectiveOffset; |
---|
3878 | double objectiveOffset2; |
---|
3879 | getDblParam(ClpObjOffset, objectiveOffset); |
---|
3880 | objectiveOffset *= whichWay; |
---|
3881 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
3882 | iColumn = listNonLinearColumn[jNon]; |
---|
3883 | double upper = columnUpper[iColumn]; |
---|
3884 | double lower = columnLower[iColumn]; |
---|
3885 | if (solution[iColumn] < lower) |
---|
3886 | solution[iColumn] = lower; |
---|
3887 | else if (solution[iColumn] > upper) |
---|
3888 | solution[iColumn] = upper; |
---|
3889 | #if 0 |
---|
3890 | double large = CoinMax(1000.0, 10.0 * fabs(solution[iColumn])); |
---|
3891 | if (upper > 1.0e10) |
---|
3892 | upper = solution[iColumn] + large; |
---|
3893 | if (lower < -1.0e10) |
---|
3894 | lower = solution[iColumn] - large; |
---|
3895 | #else |
---|
3896 | upper = solution[iColumn] + 0.5; |
---|
3897 | lower = solution[iColumn] - 0.5; |
---|
3898 | #endif |
---|
3899 | //columnUpper[iColumn]=upper; |
---|
3900 | trust[jNon] = 0.05 * (1.0 + upper - lower); |
---|
3901 | trueLower[jNon] = columnLower[iColumn]; |
---|
3902 | //trueUpper[jNon]=upper; |
---|
3903 | trueUpper[jNon] = columnUpper[iColumn]; |
---|
3904 | } |
---|
3905 | bool tryFix = false; |
---|
3906 | int iPass; |
---|
3907 | double lastObjective = 1.0e31; |
---|
3908 | double lastGoodObjective = 1.0e31; |
---|
3909 | double *bestSolution = NULL; |
---|
3910 | double *saveSolution = new double[numberColumns2 + numberRows]; |
---|
3911 | char *saveStatus = new char[numberColumns2 + numberRows]; |
---|
3912 | double targetDrop = 1.0e31; |
---|
3913 | // 1 bound up, 2 up, -1 bound down, -2 down, 0 no change |
---|
3914 | for (iPass = 0; iPass < 3; iPass++) { |
---|
3915 | last[iPass] = new int[numberNonLinearColumns]; |
---|
3916 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) |
---|
3917 | last[iPass][jNon] = 0; |
---|
3918 | } |
---|
3919 | int numberZeroPasses = 0; |
---|
3920 | bool zeroTargetDrop = false; |
---|
3921 | double *gradient = new double[numberColumns_]; |
---|
3922 | bool goneFeasible = false; |
---|
3923 | // keep sum of artificials |
---|
3924 | #define KEEP_SUM 5 |
---|
3925 | double sumArt[KEEP_SUM]; |
---|
3926 | for (jNon = 0; jNon < KEEP_SUM; jNon++) |
---|
3927 | sumArt[jNon] = COIN_DBL_MAX; |
---|
3928 | #define SMALL_FIX 0.0 |
---|
3929 | for (iPass = 0; iPass < numberPasses; iPass++) { |
---|
3930 | objectiveOffset2 = objectiveOffset; |
---|
3931 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
3932 | iColumn = listNonLinearColumn[jNon]; |
---|
3933 | if (solution[iColumn] < trueLower[jNon]) |
---|
3934 | solution[iColumn] = trueLower[jNon]; |
---|
3935 | else if (solution[iColumn] > trueUpper[jNon]) |
---|
3936 | solution[iColumn] = trueUpper[jNon]; |
---|
3937 | columnLower[iColumn] = CoinMax(solution[iColumn] |
---|
3938 | - trust[jNon], |
---|
3939 | trueLower[jNon]); |
---|
3940 | if (!trueLower[jNon] && columnLower[iColumn] < SMALL_FIX) |
---|
3941 | columnLower[iColumn] = SMALL_FIX; |
---|
3942 | columnUpper[iColumn] = CoinMin(solution[iColumn] |
---|
3943 | + trust[jNon], |
---|
3944 | trueUpper[jNon]); |
---|
3945 | if (!trueLower[jNon]) |
---|
3946 | columnUpper[iColumn] = CoinMax(columnUpper[iColumn], |
---|
3947 | columnLower[iColumn] + SMALL_FIX); |
---|
3948 | if (!trueLower[jNon] && tryFix && columnLower[iColumn] == SMALL_FIX && columnUpper[iColumn] < 3.0 * SMALL_FIX) { |
---|
3949 | columnLower[iColumn] = 0.0; |
---|
3950 | solution[iColumn] = 0.0; |
---|
3951 | columnUpper[iColumn] = 0.0; |
---|
3952 | printf("fixing %d\n", iColumn); |
---|
3953 | } |
---|
3954 | } |
---|
3955 | // redo matrix |
---|
3956 | double offset; |
---|
3957 | CoinPackedMatrix newMatrix(saveMatrix); |
---|
3958 | // get matrix data pointers |
---|
3959 | column = newMatrix.getIndices(); |
---|
3960 | rowStart = newMatrix.getVectorStarts(); |
---|
3961 | rowLength = newMatrix.getVectorLengths(); |
---|
3962 | // make sure x updated |
---|
3963 | if (numberConstraints) |
---|
3964 | constraints[0]->newXValues(); |
---|
3965 | else |
---|
3966 | trueObjective->newXValues(); |
---|
3967 | double *changeableElement = newMatrix.getMutableElements(); |
---|
3968 | if (trueObjective) { |
---|
3969 | CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, 2), numberColumns_, |
---|
3970 | objective); |
---|
3971 | } else { |
---|
3972 | CoinMemcpyN(objective_->gradient(this, solution, offset, true, 2), numberColumns_, |
---|
3973 | objective); |
---|
3974 | } |
---|
3975 | if (whichWay < 0.0) { |
---|
3976 | for (int iColumn = 0; iColumn < numberColumns_; iColumn++) |
---|
3977 | objective[iColumn] = -objective[iColumn]; |
---|
3978 | } |
---|
3979 | for (iConstraint = 0; iConstraint < numberConstraints; iConstraint++) { |
---|
3980 | ClpConstraint *constraint = constraints[iConstraint]; |
---|
3981 | int iRow = constraint->rowNumber(); |
---|
3982 | double functionValue; |
---|
3983 | #ifndef NDEBUG |
---|
3984 | int numberErrors = |
---|
3985 | #endif |
---|
3986 | constraint->gradient(&newModel, solution, gradient, functionValue, offset); |
---|
3987 | assert(!numberErrors); |
---|
3988 | // double dualValue = newModel.dualRowSolution()[iRow]; |
---|
3989 | int numberCoefficients = constraint->numberCoefficients(); |
---|
3990 | for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + numberCoefficients; j++) { |
---|
3991 | int iColumn = column[j]; |
---|
3992 | changeableElement[j] = gradient[iColumn]; |
---|
3993 | //objective[iColumn] -= dualValue*gradient[iColumn]; |
---|
3994 | gradient[iColumn] = 0.0; |
---|
3995 | } |
---|
3996 | for (int k = 0; k < numberColumns_; k++) |
---|
3997 | assert(!gradient[k]); |
---|
3998 | if (rowLower_[iRow] > -1.0e20) |
---|
3999 | rowLower[iRow] = rowLower_[iRow] - offset; |
---|
4000 | if (rowUpper_[iRow] < 1.0e20) |
---|
4001 | rowUpper[iRow] = rowUpper_[iRow] - offset; |
---|
4002 | } |
---|
4003 | // Replace matrix |
---|
4004 | // Get a column copy in standard format |
---|
4005 | CoinPackedMatrix *columnCopy = new CoinPackedMatrix(); |
---|
4006 | ; |
---|
4007 | columnCopy->reverseOrderedCopyOf(newMatrix); |
---|
4008 | newModel.replaceMatrix(columnCopy, true); |
---|
4009 | // solve |
---|
4010 | newModel.primal(1); |
---|
4011 | if (newModel.status() == 1) { |
---|
4012 | // Infeasible! |
---|
4013 | newModel.allSlackBasis(); |
---|
4014 | newModel.primal(); |
---|
4015 | newModel.writeMps("infeas.mps"); |
---|
4016 | assert(!newModel.status()); |
---|
4017 | } |
---|
4018 | double sumInfeas = 0.0; |
---|
4019 | int numberInfeas = 0; |
---|
4020 | for (iColumn = numberColumns_; iColumn < numberColumns2; iColumn++) { |
---|
4021 | if (solution[iColumn] > 1.0e-8) { |
---|
4022 | numberInfeas++; |
---|
4023 | sumInfeas += solution[iColumn]; |
---|
4024 | } |
---|
4025 | } |
---|
4026 | printf("%d artificial infeasibilities - summing to %g\n", |
---|
4027 | numberInfeas, sumInfeas); |
---|
4028 | for (jNon = 0; jNon < KEEP_SUM - 1; jNon++) |
---|
4029 | sumArt[jNon] = sumArt[jNon + 1]; |
---|
4030 | sumArt[KEEP_SUM - 1] = sumInfeas; |
---|
4031 | if (sumInfeas > 0.01 && sumInfeas * 1.1 > sumArt[0] && penalties[1] < 1.0e7) { |
---|
4032 | // not doing very well - increase - be more sophisticated later |
---|
4033 | lastObjective = COIN_DBL_MAX; |
---|
4034 | for (jNon = 0; jNon < KEEP_SUM; jNon++) |
---|
4035 | sumArt[jNon] = COIN_DBL_MAX; |
---|
4036 | //for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn+=SEGMENTS) { |
---|
4037 | //objective[iColumn+1] *= 1.5; |
---|
4038 | //} |
---|
4039 | penalties[1] *= 1.5; |
---|
4040 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) |
---|
4041 | if (trust[jNon] > 0.1) |
---|
4042 | trust[jNon] *= 2.0; |
---|
4043 | else |
---|
4044 | trust[jNon] = 0.1; |
---|
4045 | } |
---|
4046 | if (sumInfeas < 0.001 && !goneFeasible) { |
---|
4047 | goneFeasible = true; |
---|
4048 | penalties[0] = 1.0e-3; |
---|
4049 | penalties[1] = 1.0e6; |
---|
4050 | printf("Got feasible\n"); |
---|
4051 | } |
---|
4052 | double infValue = 0.0; |
---|
4053 | double objValue = 0.0; |
---|
4054 | // make sure x updated |
---|
4055 | if (numberConstraints) |
---|
4056 | constraints[0]->newXValues(); |
---|
4057 | else |
---|
4058 | trueObjective->newXValues(); |
---|
4059 | if (trueObjective) { |
---|
4060 | objValue = trueObjective->objectiveValue(this, solution); |
---|
4061 | printf("objective offset %g\n", offset); |
---|
4062 | objectiveOffset2 = objectiveOffset + offset; // ? sign |
---|
4063 | newModel.setObjectiveOffset(objectiveOffset2); |
---|
4064 | } else { |
---|
4065 | objValue = objective_->objectiveValue(this, solution); |
---|
4066 | } |
---|
4067 | objValue *= whichWay; |
---|
4068 | double infPenalty = 0.0; |
---|
4069 | // This penalty is for target drop |
---|
4070 | double infPenalty2 = 0.0; |
---|
4071 | //const int * row = columnCopy->getIndices(); |
---|
4072 | //const CoinBigIndex * columnStart = columnCopy->getVectorStarts(); |
---|
4073 | //const int * columnLength = columnCopy->getVectorLengths(); |
---|
4074 | //const double * element = columnCopy->getElements(); |
---|
4075 | double *cost = newModel.objective(); |
---|
4076 | column = newMatrix.getIndices(); |
---|
4077 | rowStart = newMatrix.getVectorStarts(); |
---|
4078 | rowLength = newMatrix.getVectorLengths(); |
---|
4079 | elementByRow = newMatrix.getElements(); |
---|
4080 | int jColumn = numberColumns_; |
---|
4081 | double objectiveAdjustment = 0.0; |
---|
4082 | for (iConstraint = 0; iConstraint < numberConstraints; iConstraint++) { |
---|
4083 | ClpConstraint *constraint = constraints[iConstraint]; |
---|
4084 | int iRow = constraint->rowNumber(); |
---|
4085 | double functionValue = constraint->functionValue(this, solution); |
---|
4086 | double dualValue = newModel.dualRowSolution()[iRow]; |
---|
4087 | if (numberConstraints < -50) |
---|
4088 | printf("For row %d current value is %g (row activity %g) , dual is %g\n", iRow, functionValue, |
---|
4089 | newModel.primalRowSolution()[iRow], |
---|
4090 | dualValue); |
---|
4091 | double movement = newModel.primalRowSolution()[iRow] + constraint->offset(); |
---|
4092 | movement = fabs((movement - functionValue) * dualValue); |
---|
4093 | infPenalty2 += movement; |
---|
4094 | double sumOfActivities = 0.0; |
---|
4095 | for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) { |
---|
4096 | int iColumn = column[j]; |
---|
4097 | sumOfActivities += fabs(solution[iColumn] * elementByRow[j]); |
---|
4098 | } |
---|
4099 | if (rowLower_[iRow] > -1.0e20) { |
---|
4100 | if (functionValue < rowLower_[iRow] - 1.0e-5) { |
---|
4101 | double infeasibility = rowLower_[iRow] - functionValue; |
---|
4102 | double thisPenalty = 0.0; |
---|
4103 | infValue += infeasibility; |
---|
4104 | int k; |
---|
4105 | assert(dualValue >= -1.0e-5); |
---|
4106 | dualValue = CoinMax(dualValue, 0.0); |
---|
4107 | for (k = 0; k < SEGMENTS; k++) { |
---|
4108 | if (infeasibility <= 0) |
---|
4109 | break; |
---|
4110 | double thisPart = CoinMin(infeasibility, bounds[k]); |
---|
4111 | thisPenalty += thisPart * cost[jColumn + k]; |
---|
4112 | infeasibility -= thisPart; |
---|
4113 | } |
---|
4114 | infeasibility = functionValue - rowUpper_[iRow]; |
---|
4115 | double newPenalty = 0.0; |
---|
4116 | for (k = 0; k < SEGMENTS; k++) { |
---|
4117 | double thisPart = CoinMin(infeasibility, bounds[k]); |
---|
4118 | cost[jColumn + k] = CoinMax(penalties[k], dualValue + 1.0e-3); |
---|
4119 | newPenalty += thisPart * cost[jColumn + k]; |
---|
4120 | infeasibility -= thisPart; |
---|
4121 | } |
---|
4122 | infPenalty += thisPenalty; |
---|
4123 | objectiveAdjustment += CoinMax(0.0, newPenalty - thisPenalty); |
---|
4124 | } |
---|
4125 | jColumn += SEGMENTS; |
---|
4126 | } |
---|
4127 | if (rowUpper_[iRow] < 1.0e20) { |
---|
4128 | if (functionValue > rowUpper_[iRow] + 1.0e-5) { |
---|
4129 | double infeasibility = functionValue - rowUpper_[iRow]; |
---|
4130 | double thisPenalty = 0.0; |
---|
4131 | infValue += infeasibility; |
---|
4132 | int k; |
---|
4133 | dualValue = -dualValue; |
---|
4134 | assert(dualValue >= -1.0e-5); |
---|
4135 | dualValue = CoinMax(dualValue, 0.0); |
---|
4136 | for (k = 0; k < SEGMENTS; k++) { |
---|
4137 | if (infeasibility <= 0) |
---|
4138 | break; |
---|
4139 | double thisPart = CoinMin(infeasibility, bounds[k]); |
---|
4140 | thisPenalty += thisPart * cost[jColumn + k]; |
---|
4141 | infeasibility -= thisPart; |
---|
4142 | } |
---|
4143 | infeasibility = functionValue - rowUpper_[iRow]; |
---|
4144 | double newPenalty = 0.0; |
---|
4145 | for (k = 0; k < SEGMENTS; k++) { |
---|
4146 | double thisPart = CoinMin(infeasibility, bounds[k]); |
---|
4147 | cost[jColumn + k] = CoinMax(penalties[k], dualValue + 1.0e-3); |
---|
4148 | newPenalty += thisPart * cost[jColumn + k]; |
---|
4149 | infeasibility -= thisPart; |
---|
4150 | } |
---|
4151 | infPenalty += thisPenalty; |
---|
4152 | objectiveAdjustment += CoinMax(0.0, newPenalty - thisPenalty); |
---|
4153 | } |
---|
4154 | jColumn += SEGMENTS; |
---|
4155 | } |
---|
4156 | } |
---|
4157 | // adjust last objective value |
---|
4158 | lastObjective += objectiveAdjustment; |
---|
4159 | if (infValue) |
---|
4160 | printf("Sum infeasibilities %g - penalty %g ", infValue, infPenalty); |
---|
4161 | if (objectiveOffset2) |
---|
4162 | printf("offset2 %g ", objectiveOffset2); |
---|
4163 | objValue -= objectiveOffset2; |
---|
4164 | printf("True objective %g or maybe %g (with penalty %g) -pen2 %g %g\n", objValue, |
---|
4165 | objValue + objectiveOffset2, objValue + objectiveOffset2 + infPenalty, infPenalty2, penalties[1]); |
---|
4166 | double useObjValue = objValue + objectiveOffset2 + infPenalty; |
---|
4167 | objValue += infPenalty + infPenalty2; |
---|
4168 | objValue = useObjValue; |
---|
4169 | if (iPass) { |
---|
4170 | double drop = lastObjective - objValue; |
---|
4171 | std::cout << "True drop was " << drop << std::endl; |
---|
4172 | if (drop < -0.05 * fabs(objValue) - 1.0e-4) { |
---|
4173 | // pretty bad - go back and halve |
---|
4174 | CoinMemcpyN(saveSolution, numberColumns2, solution); |
---|
4175 | CoinMemcpyN(saveSolution + numberColumns2, |
---|
4176 | numberRows, newModel.primalRowSolution()); |
---|
4177 | CoinMemcpyN(reinterpret_cast< unsigned char * >(saveStatus), |
---|
4178 | numberColumns2 + numberRows, newModel.statusArray()); |
---|
4179 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) |
---|
4180 | if (trust[jNon] > 0.1) |
---|
4181 | trust[jNon] *= 0.5; |
---|
4182 | else |
---|
4183 | trust[jNon] *= 0.9; |
---|
4184 | |
---|
4185 | printf("** Halving trust\n"); |
---|
4186 | objValue = lastObjective; |
---|
4187 | continue; |
---|
4188 | } else if ((iPass % 25) == -1) { |
---|
4189 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) |
---|
4190 | trust[jNon] *= 2.0; |
---|
4191 | printf("** Doubling trust\n"); |
---|
4192 | } |
---|
4193 | int *temp = last[2]; |
---|
4194 | last[2] = last[1]; |
---|
4195 | last[1] = last[0]; |
---|
4196 | last[0] = temp; |
---|
4197 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
4198 | iColumn = listNonLinearColumn[jNon]; |
---|
4199 | double change = solution[iColumn] - saveSolution[iColumn]; |
---|
4200 | if (change < -1.0e-5) { |
---|
4201 | if (fabs(change + trust[jNon]) < 1.0e-5) |
---|
4202 | temp[jNon] = -1; |
---|
4203 | else |
---|
4204 | temp[jNon] = -2; |
---|
4205 | } else if (change > 1.0e-5) { |
---|
4206 | if (fabs(change - trust[jNon]) < 1.0e-5) |
---|
4207 | temp[jNon] = 1; |
---|
4208 | else |
---|
4209 | temp[jNon] = 2; |
---|
4210 | } else { |
---|
4211 | temp[jNon] = 0; |
---|
4212 | } |
---|
4213 | } |
---|
4214 | double maxDelta = 0.0; |
---|
4215 | double smallestTrust = 1.0e31; |
---|
4216 | double smallestNonLinearGap = 1.0e31; |
---|
4217 | double smallestGap = 1.0e31; |
---|
4218 | bool increasing = false; |
---|
4219 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
---|
4220 | double gap = columnUpper[iColumn] - columnLower[iColumn]; |
---|
4221 | assert(gap >= 0.0); |
---|
4222 | if (gap) |
---|
4223 | smallestGap = CoinMin(smallestGap, gap); |
---|
4224 | } |
---|
4225 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
4226 | iColumn = listNonLinearColumn[jNon]; |
---|
4227 | double gap = columnUpper[iColumn] - columnLower[iColumn]; |
---|
4228 | assert(gap >= 0.0); |
---|
4229 | if (gap) { |
---|
4230 | smallestNonLinearGap = CoinMin(smallestNonLinearGap, gap); |
---|
4231 | if (gap < 1.0e-7 && iPass == 1) { |
---|
4232 | printf("Small gap %d %d %g %g %g\n", |
---|
4233 | jNon, iColumn, columnLower[iColumn], columnUpper[iColumn], |
---|
4234 | gap); |
---|
4235 | //trueUpper[jNon]=trueLower[jNon]; |
---|
4236 | //columnUpper[iColumn]=columnLower[iColumn]; |
---|
4237 | } |
---|
4238 | } |
---|
4239 | maxDelta = CoinMax(maxDelta, |
---|
4240 | fabs(solution[iColumn] - saveSolution[iColumn])); |
---|
4241 | if (last[0][jNon] * last[1][jNon] < 0) { |
---|
4242 | // halve |
---|
4243 | if (trust[jNon] > 1.0) |
---|
4244 | trust[jNon] *= 0.5; |
---|
4245 | else |
---|
4246 | trust[jNon] *= 0.7; |
---|
4247 | } else { |
---|
4248 | // ? only increase if +=1 ? |
---|
4249 | if (last[0][jNon] == last[1][jNon] && last[0][jNon] == last[2][jNon] && last[0][jNon]) { |
---|
4250 | trust[jNon] *= 1.8; |
---|
4251 | increasing = true; |
---|
4252 | } |
---|
4253 | } |
---|
4254 | smallestTrust = CoinMin(smallestTrust, trust[jNon]); |
---|
4255 | } |
---|
4256 | std::cout << "largest delta is " << maxDelta |
---|
4257 | << ", smallest trust is " << smallestTrust |
---|
4258 | << ", smallest gap is " << smallestGap |
---|
4259 | << ", smallest nonlinear gap is " << smallestNonLinearGap |
---|
4260 | << std::endl; |
---|
4261 | if (iPass > 200) { |
---|
4262 | //double useObjValue = objValue+objectiveOffset2+infPenalty; |
---|
4263 | if (useObjValue + 1.0e-4 > lastGoodObjective && iPass > 250) { |
---|
4264 | std::cout << "Exiting as objective not changing much" << std::endl; |
---|
4265 | break; |
---|
4266 | } else if (useObjValue < lastGoodObjective) { |
---|
4267 | lastGoodObjective = useObjValue; |
---|
4268 | if (!bestSolution) |
---|
4269 | bestSolution = new double[numberColumns2]; |
---|
4270 | CoinMemcpyN(solution, numberColumns2, bestSolution); |
---|
4271 | } |
---|
4272 | } |
---|
4273 | if (maxDelta < deltaTolerance && !increasing && iPass > 100) { |
---|
4274 | numberZeroPasses++; |
---|
4275 | if (numberZeroPasses == 3) { |
---|
4276 | if (tryFix) { |
---|
4277 | std::cout << "Exiting" << std::endl; |
---|
4278 | break; |
---|
4279 | } else { |
---|
4280 | tryFix = true; |
---|
4281 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) |
---|
4282 | trust[jNon] = CoinMax(trust[jNon], 1.0e-1); |
---|
4283 | numberZeroPasses = 0; |
---|
4284 | } |
---|
4285 | } |
---|
4286 | } else { |
---|
4287 | numberZeroPasses = 0; |
---|
4288 | } |
---|
4289 | } |
---|
4290 | CoinMemcpyN(solution, numberColumns2, saveSolution); |
---|
4291 | CoinMemcpyN(newModel.primalRowSolution(), |
---|
4292 | numberRows, saveSolution + numberColumns2); |
---|
4293 | CoinMemcpyN(newModel.statusArray(), |
---|
4294 | numberColumns2 + numberRows, |
---|
4295 | reinterpret_cast< unsigned char * >(saveStatus)); |
---|
4296 | |
---|
4297 | targetDrop = infPenalty + infPenalty2; |
---|
4298 | if (iPass) { |
---|
4299 | // get reduced costs |
---|
4300 | const double *pi = newModel.dualRowSolution(); |
---|
4301 | newModel.matrix()->transposeTimes(pi, |
---|
4302 | newModel.dualColumnSolution()); |
---|
4303 | double *r = newModel.dualColumnSolution(); |
---|
4304 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) |
---|
4305 | r[iColumn] = objective[iColumn] - r[iColumn]; |
---|
4306 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
4307 | iColumn = listNonLinearColumn[jNon]; |
---|
4308 | double dj = r[iColumn]; |
---|
4309 | if (dj < -1.0e-6) { |
---|
4310 | double drop = -dj * (columnUpper[iColumn] - solution[iColumn]); |
---|
4311 | //double upper = CoinMin(trueUpper[jNon],solution[iColumn]+0.1); |
---|
4312 | //double drop2 = -dj*(upper-solution[iColumn]); |
---|
4313 | #if 0 |
---|
4314 | if (drop > 1.0e8 || drop2 > 100.0 * drop || (drop > 1.0e-2 && iPass > 100)) |
---|
4315 | printf("Big drop %d %g %g %g %g T %g\n", |
---|
4316 | iColumn, columnLower[iColumn], solution[iColumn], |
---|
4317 | columnUpper[iColumn], dj, trueUpper[jNon]); |
---|
4318 | #endif |
---|
4319 | targetDrop += drop; |
---|
4320 | if (dj < -1.0e-1 && trust[jNon] < 1.0e-3 |
---|
4321 | && trueUpper[jNon] - solution[iColumn] > 1.0e-2) { |
---|
4322 | trust[jNon] *= 1.5; |
---|
4323 | //printf("Increasing trust on %d to %g\n", |
---|
4324 | // iColumn,trust[jNon]); |
---|
4325 | } |
---|
4326 | } else if (dj > 1.0e-6) { |
---|
4327 | double drop = -dj * (columnLower[iColumn] - solution[iColumn]); |
---|
4328 | //double lower = CoinMax(trueLower[jNon],solution[iColumn]-0.1); |
---|
4329 | //double drop2 = -dj*(lower-solution[iColumn]); |
---|
4330 | #if 0 |
---|
4331 | if (drop > 1.0e8 || drop2 > 100.0 * drop || (drop > 1.0e-2)) |
---|
4332 | printf("Big drop %d %g %g %g %g T %g\n", |
---|
4333 | iColumn, columnLower[iColumn], solution[iColumn], |
---|
4334 | columnUpper[iColumn], dj, trueLower[jNon]); |
---|
4335 | #endif |
---|
4336 | targetDrop += drop; |
---|
4337 | if (dj > 1.0e-1 && trust[jNon] < 1.0e-3 |
---|
4338 | && solution[iColumn] - trueLower[jNon] > 1.0e-2) { |
---|
4339 | trust[jNon] *= 1.5; |
---|
4340 | printf("Increasing trust on %d to %g\n", |
---|
4341 | iColumn, trust[jNon]); |
---|
4342 | } |
---|
4343 | } |
---|
4344 | } |
---|
4345 | } |
---|
4346 | std::cout << "Pass - " << iPass |
---|
4347 | << ", target drop is " << targetDrop |
---|
4348 | << std::endl; |
---|
4349 | if (iPass > 1 && targetDrop < 1.0e-5 && zeroTargetDrop) |
---|
4350 | break; |
---|
4351 | if (iPass > 1 && targetDrop < 1.0e-5) |
---|
4352 | zeroTargetDrop = true; |
---|
4353 | else |
---|
4354 | zeroTargetDrop = false; |
---|
4355 | //if (iPass==5) |
---|
4356 | //newModel.setLogLevel(63); |
---|
4357 | lastObjective = objValue; |
---|
4358 | // take out when ClpPackedMatrix changed |
---|
4359 | //newModel.scaling(false); |
---|
4360 | #if 0 |
---|
4361 | CoinMpsIO writer; |
---|
4362 | writer.setMpsData(*newModel.matrix(), COIN_DBL_MAX, |
---|
4363 | newModel.getColLower(), newModel.getColUpper(), |
---|
4364 | newModel.getObjCoefficients(), |
---|
4365 | (const char*) 0 /*integrality*/, |
---|
4366 | newModel.getRowLower(), newModel.getRowUpper(), |
---|
4367 | NULL, NULL); |
---|
4368 | writer.writeMps("xx.mps"); |
---|
4369 | #endif |
---|
4370 | } |
---|
4371 | delete[] saveSolution; |
---|
4372 | delete[] saveStatus; |
---|
4373 | for (iPass = 0; iPass < 3; iPass++) |
---|
4374 | delete[] last[iPass]; |
---|
4375 | for (jNon = 0; jNon < numberNonLinearColumns; jNon++) { |
---|
4376 | iColumn = listNonLinearColumn[jNon]; |
---|
4377 | columnLower[iColumn] = trueLower[jNon]; |
---|
4378 | columnUpper[iColumn] = trueUpper[jNon]; |
---|
4379 | } |
---|
4380 | delete[] trust; |
---|
4381 | delete[] trueUpper; |
---|
4382 | delete[] trueLower; |
---|
4383 | objectiveValue_ = newModel.objectiveValue(); |
---|
4384 | if (bestSolution) { |
---|
4385 | CoinMemcpyN(bestSolution, numberColumns2, solution); |
---|
4386 | delete[] bestSolution; |
---|
4387 | printf("restoring objective of %g\n", lastGoodObjective); |
---|
4388 | objectiveValue_ = lastGoodObjective; |
---|
4389 | } |
---|
4390 | // Simplest way to get true row activity ? |
---|
4391 | double *rowActivity = newModel.primalRowSolution(); |
---|
4392 | for (iRow = 0; iRow < numberRows; iRow++) { |
---|
4393 | double difference; |
---|
4394 | if (fabs(rowLower_[iRow]) < fabs(rowUpper_[iRow])) |
---|
4395 | difference = rowLower_[iRow] - rowLower[iRow]; |
---|
4396 | else |
---|
4397 | difference = rowUpper_[iRow] - rowUpper[iRow]; |
---|
4398 | rowLower[iRow] = rowLower_[iRow]; |
---|
4399 | rowUpper[iRow] = rowUpper_[iRow]; |
---|
4400 | if (difference) { |
---|
4401 | if (numberRows < 50) |
---|
4402 | printf("For row %d activity changes from %g to %g\n", |
---|
4403 | iRow, rowActivity[iRow], rowActivity[iRow] + difference); |
---|
4404 | rowActivity[iRow] += difference; |
---|
4405 | } |
---|
4406 | } |
---|
4407 | delete[] listNonLinearColumn; |
---|
4408 | delete[] gradient; |
---|
4409 | printf("solution still in newModel - do objective etc!\n"); |
---|
4410 | numberIterations_ = newModel.numberIterations(); |
---|
4411 | problemStatus_ = newModel.problemStatus(); |
---|
4412 | secondaryStatus_ = newModel.secondaryStatus(); |
---|
4413 | CoinMemcpyN(newModel.primalColumnSolution(), numberColumns_, columnActivity_); |
---|
4414 | // should do status region |
---|
4415 | CoinZeroN(rowActivity_, numberRows_); |
---|
4416 | matrix_->times(1.0, columnActivity_, rowActivity_); |
---|
4417 | return 0; |
---|
4418 | } |
---|
4419 | |
---|
4420 | /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 |
---|
4421 | */ |
---|