1/* $Id: ClpQuadraticObjective.cpp 1753 2011-06-19 16:27:26Z stefan $ */
2// Copyright (C) 2003, 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#include "CoinHelperFunctions.hpp"
8#include "CoinIndexedVector.hpp"
9#include "ClpFactorization.hpp"
10#include "ClpSimplex.hpp"
11#include "ClpQuadraticObjective.hpp"
12//#############################################################################
13// Constructors / Destructor / Assignment
14//#############################################################################
15//-------------------------------------------------------------------
16// Default Constructor
17//-------------------------------------------------------------------
18ClpQuadraticObjective::ClpQuadraticObjective ()
19 : ClpObjective()
20{
21 type_ = 2;
22 objective_ = NULL;
23 quadraticObjective_ = NULL;
24 gradient_ = NULL;
25 numberColumns_ = 0;
26 numberExtendedColumns_ = 0;
27 activated_ = 0;
28 fullMatrix_ = false;
29}
30
31//-------------------------------------------------------------------
32// Useful Constructor
33//-------------------------------------------------------------------
34ClpQuadraticObjective::ClpQuadraticObjective (const double * objective ,
35 int numberColumns,
36 const CoinBigIndex * start,
37 const int * column, const double * element,
38 int numberExtendedColumns)
39 : ClpObjective()
40{
41 type_ = 2;
42 numberColumns_ = numberColumns;
43 if (numberExtendedColumns >= 0)
44 numberExtendedColumns_ = CoinMax(numberColumns_, numberExtendedColumns);
45 else
46 numberExtendedColumns_ = numberColumns_;
47 if (objective) {
48 objective_ = new double [numberExtendedColumns_];
49 CoinMemcpyN(objective, numberColumns_, objective_);
50 memset(objective_ + numberColumns_, 0, (numberExtendedColumns_ - numberColumns_)*sizeof(double));
51 } else {
52 objective_ = new double [numberExtendedColumns_];
53 memset(objective_, 0, numberExtendedColumns_ * sizeof(double));
54 }
55 if (start)
56 quadraticObjective_ = new CoinPackedMatrix(true, numberColumns, numberColumns,
57 start[numberColumns], element, column, start, NULL);
58 else
59 quadraticObjective_ = NULL;
60 gradient_ = NULL;
61 activated_ = 1;
62 fullMatrix_ = false;
63}
64
65//-------------------------------------------------------------------
66// Copy constructor
67//-------------------------------------------------------------------
68ClpQuadraticObjective::ClpQuadraticObjective (const ClpQuadraticObjective & rhs,
69 int type)
70 : ClpObjective(rhs)
71{
72 numberColumns_ = rhs.numberColumns_;
73 numberExtendedColumns_ = rhs.numberExtendedColumns_;
74 fullMatrix_ = rhs.fullMatrix_;
75 if (rhs.objective_) {
76 objective_ = new double [numberExtendedColumns_];
77 CoinMemcpyN(rhs.objective_, numberExtendedColumns_, objective_);
78 } else {
79 objective_ = NULL;
80 }
81 if (rhs.gradient_) {
82 gradient_ = new double [numberExtendedColumns_];
83 CoinMemcpyN(rhs.gradient_, numberExtendedColumns_, gradient_);
84 } else {
85 gradient_ = NULL;
86 }
87 if (rhs.quadraticObjective_) {
88 // see what type of matrix wanted
89 if (type == 0) {
90 // just copy
91 quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
92 } else if (type == 1) {
93 // expand to full symmetric
94 fullMatrix_ = true;
95 const int * columnQuadratic1 = rhs.quadraticObjective_->getIndices();
96 const CoinBigIndex * columnQuadraticStart1 = rhs.quadraticObjective_->getVectorStarts();
97 const int * columnQuadraticLength1 = rhs.quadraticObjective_->getVectorLengths();
98 const double * quadraticElement1 = rhs.quadraticObjective_->getElements();
99 CoinBigIndex * columnQuadraticStart2 = new CoinBigIndex [numberExtendedColumns_+1];
100 int * columnQuadraticLength2 = new int [numberExtendedColumns_];
101 int iColumn;
102 int numberColumns = rhs.quadraticObjective_->getNumCols();
103 int numberBelow = 0;
104 int numberAbove = 0;
105 int numberDiagonal = 0;
106 CoinZeroN(columnQuadraticLength2, numberExtendedColumns_);
107 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
108 for (CoinBigIndex j = columnQuadraticStart1[iColumn];
109 j < columnQuadraticStart1[iColumn] + columnQuadraticLength1[iColumn]; j++) {
110 int jColumn = columnQuadratic1[j];
111 if (jColumn > iColumn) {
112 numberBelow++;
113 columnQuadraticLength2[jColumn]++;
114 columnQuadraticLength2[iColumn]++;
115 } else if (jColumn == iColumn) {
116 numberDiagonal++;
117 columnQuadraticLength2[iColumn]++;
118 } else {
119 numberAbove++;
120 }
121 }
122 }
123 if (numberAbove > 0) {
124 if (numberAbove == numberBelow) {
125 // already done
126 quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
127 delete [] columnQuadraticStart2;
128 delete [] columnQuadraticLength2;
129 } else {
130 printf("number above = %d, number below = %d, error\n",
131 numberAbove, numberBelow);
132 abort();
133 }
134 } else {
135 int numberElements = numberDiagonal + 2 * numberBelow;
136 int * columnQuadratic2 = new int [numberElements];
137 double * quadraticElement2 = new double [numberElements];
138 columnQuadraticStart2[0] = 0;
139 numberElements = 0;
140 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
141 int n = columnQuadraticLength2[iColumn];
142 columnQuadraticLength2[iColumn] = 0;
143 numberElements += n;
144 columnQuadraticStart2[iColumn+1] = numberElements;
145 }
146 for (iColumn = 0; iColumn < numberColumns; iColumn++) {
147 for (CoinBigIndex j = columnQuadraticStart1[iColumn];
148 j < columnQuadraticStart1[iColumn] + columnQuadraticLength1[iColumn]; j++) {
149 int jColumn = columnQuadratic1[j];
150 if (jColumn > iColumn) {
151 // put in two places
152 CoinBigIndex put = columnQuadraticLength2[jColumn] + columnQuadraticStart2[jColumn];
153 columnQuadraticLength2[jColumn]++;
154 quadraticElement2[put] = quadraticElement1[j];
155 columnQuadratic2[put] = iColumn;
156 put = columnQuadraticLength2[iColumn] + columnQuadraticStart2[iColumn];
157 columnQuadraticLength2[iColumn]++;
158 quadraticElement2[put] = quadraticElement1[j];
159 columnQuadratic2[put] = jColumn;
160 } else if (jColumn == iColumn) {
161 CoinBigIndex put = columnQuadraticLength2[iColumn] + columnQuadraticStart2[iColumn];
162 columnQuadraticLength2[iColumn]++;
163 quadraticElement2[put] = quadraticElement1[j];
164 columnQuadratic2[put] = iColumn;
165 } else {
166 abort();
167 }
168 }
169 }
170 // Now create
171 quadraticObjective_ =
172 new CoinPackedMatrix (true,
173 rhs.numberExtendedColumns_,
174 rhs.numberExtendedColumns_,
175 numberElements,
176 quadraticElement2,
177 columnQuadratic2,
178 columnQuadraticStart2,
179 columnQuadraticLength2, 0.0, 0.0);
180 delete [] columnQuadraticStart2;
181 delete [] columnQuadraticLength2;
182 delete [] columnQuadratic2;
183 delete [] quadraticElement2;
184 }
185 } else {
186 fullMatrix_ = false;
187 abort(); // code when needed
188 }
189
190 } else {
191 quadraticObjective_ = NULL;
192 }
193}
194/* Subset constructor. Duplicates are allowed
195 and order is as given.
196*/
197ClpQuadraticObjective::ClpQuadraticObjective (const ClpQuadraticObjective &rhs,
198 int numberColumns,
199 const int * whichColumn)
200 : ClpObjective(rhs)
201{
202 fullMatrix_ = rhs.fullMatrix_;
203 objective_ = NULL;
204 int extra = rhs.numberExtendedColumns_ - rhs.numberColumns_;
205 numberColumns_ = 0;
206 numberExtendedColumns_ = numberColumns + extra;
207 if (numberColumns > 0) {
208 // check valid lists
209 int numberBad = 0;
210 int i;
211 for (i = 0; i < numberColumns; i++)
212 if (whichColumn[i] < 0 || whichColumn[i] >= rhs.numberColumns_)
213 numberBad++;
214 if (numberBad)
215 throw CoinError("bad column list", "subset constructor",
216 "ClpQuadraticObjective");
217 numberColumns_ = numberColumns;
218 objective_ = new double[numberExtendedColumns_];
219 for (i = 0; i < numberColumns_; i++)
220 objective_[i] = rhs.objective_[whichColumn[i]];
221 CoinMemcpyN(rhs.objective_ + rhs.numberColumns_, (numberExtendedColumns_ - numberColumns_),
222 objective_ + numberColumns_);
223 if (rhs.gradient_) {
224 gradient_ = new double[numberExtendedColumns_];
225 for (i = 0; i < numberColumns_; i++)
226 gradient_[i] = rhs.gradient_[whichColumn[i]];
227 CoinMemcpyN(rhs.gradient_ + rhs.numberColumns_, (numberExtendedColumns_ - numberColumns_),
228 gradient_ + numberColumns_);
229 } else {
230 gradient_ = NULL;
231 }
232 } else {
233 gradient_ = NULL;
234 objective_ = NULL;
235 }
236 if (rhs.quadraticObjective_) {
237 quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_,
238 numberColumns, whichColumn,
239 numberColumns, whichColumn);
240 } else {
241 quadraticObjective_ = NULL;
242 }
243}
244
245
246//-------------------------------------------------------------------
247// Destructor
248//-------------------------------------------------------------------
249ClpQuadraticObjective::~ClpQuadraticObjective ()
250{
251 delete [] objective_;
252 delete [] gradient_;
253 delete quadraticObjective_;
254}
255
256//----------------------------------------------------------------
257// Assignment operator
258//-------------------------------------------------------------------
259ClpQuadraticObjective &
260ClpQuadraticObjective::operator=(const ClpQuadraticObjective& rhs)
261{
262 if (this != &rhs) {
263 fullMatrix_ = rhs.fullMatrix_;
264 delete quadraticObjective_;
265 quadraticObjective_ = NULL;
266 delete [] objective_;
267 delete [] gradient_;
268 ClpObjective::operator=(rhs);
269 numberColumns_ = rhs.numberColumns_;
270 numberExtendedColumns_ = rhs.numberExtendedColumns_;
271 if (rhs.objective_) {
272 objective_ = new double [numberExtendedColumns_];
273 CoinMemcpyN(rhs.objective_, numberExtendedColumns_, objective_);
274 } else {
275 objective_ = NULL;
276 }
277 if (rhs.gradient_) {
278 gradient_ = new double [numberExtendedColumns_];
279 CoinMemcpyN(rhs.gradient_, numberExtendedColumns_, gradient_);
280 } else {
281 gradient_ = NULL;
282 }
283 if (rhs.quadraticObjective_) {
284 quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
285 } else {
286 quadraticObjective_ = NULL;
287 }
288 }
289 return *this;
290}
291
292// Returns gradient
293double *
294ClpQuadraticObjective::gradient(const ClpSimplex * model,
295 const double * solution, double & offset, bool refresh,
296 int includeLinear)
297{
298 offset = 0.0;
299 bool scaling = false;
300 if (model && (model->rowScale() ||
301 model->objectiveScale() != 1.0 || model->optimizationDirection() != 1.0))
302 scaling = true;
303 const double * cost = NULL;
304 if (model)
305 cost = model->costRegion();
306 if (!cost) {
307 // not in solve
308 cost = objective_;
309 scaling = false;
310 }
311 if (!scaling) {
312 if (!quadraticObjective_ || !solution || !activated_) {
313 return objective_;
314 } else {
315 if (refresh || !gradient_) {
316 if (!gradient_)
317 gradient_ = new double[numberExtendedColumns_];
318 const int * columnQuadratic = quadraticObjective_->getIndices();
319 const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
320 const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
321 const double * quadraticElement = quadraticObjective_->getElements();
322 offset = 0.0;
323 // use current linear cost region
324 if (includeLinear == 1)
325 CoinMemcpyN(cost, numberExtendedColumns_, gradient_);
326 else if (includeLinear == 2)
327 CoinMemcpyN(objective_, numberExtendedColumns_, gradient_);
328 else
329 memset(gradient_, 0, numberExtendedColumns_ * sizeof(double));
330 if (activated_) {
331 if (!fullMatrix_) {
332 int iColumn;
333 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
334 double valueI = solution[iColumn];
335 CoinBigIndex j;
336 for (j = columnQuadraticStart[iColumn];
337 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
338 int jColumn = columnQuadratic[j];
339 double valueJ = solution[jColumn];
340 double elementValue = quadraticElement[j];
341 if (iColumn != jColumn) {
342 offset += valueI * valueJ * elementValue;
343 //if (fabs(valueI*valueJ*elementValue)>1.0e-12)
344 //printf("%d %d %g %g %g -> %g\n",
345 // iColumn,jColumn,valueI,valueJ,elementValue,
346 // valueI*valueJ*elementValue);
347 double gradientI = valueJ * elementValue;
348 double gradientJ = valueI * elementValue;
349 gradient_[iColumn] += gradientI;
350 gradient_[jColumn] += gradientJ;
351 } else {
352 offset += 0.5 * valueI * valueI * elementValue;
353 //if (fabs(valueI*valueI*elementValue)>1.0e-12)
354 //printf("XX %d %g %g -> %g\n",
355 // iColumn,valueI,elementValue,
356 // 0.5*valueI*valueI*elementValue);
357 double gradientI = valueI * elementValue;
358 gradient_[iColumn] += gradientI;
359 }
360 }
361 }
362 } else {
363 // full matrix
364 int iColumn;
365 offset *= 2.0;
366 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
367 CoinBigIndex j;
368 double value = 0.0;
369 double current = gradient_[iColumn];
370 for (j = columnQuadraticStart[iColumn];
371 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
372 int jColumn = columnQuadratic[j];
373 double valueJ = solution[jColumn] * quadraticElement[j];
374 value += valueJ;
375 }
376 offset += value * solution[iColumn];
377 gradient_[iColumn] = current + value;
378 }
379 offset *= 0.5;
380 }
381 }
382 }
383 if (model)
384 offset *= model->optimizationDirection() * model->objectiveScale();
385 return gradient_;
386 }
387 } else {
388 // do scaling
389 assert (solution);
390 // for now only if half
391 assert (!fullMatrix_);
392 if (refresh || !gradient_) {
393 if (!gradient_)
394 gradient_ = new double[numberExtendedColumns_];
395 double direction = model->optimizationDirection() * model->objectiveScale();
396 // direction is actually scale out not scale in
397 //if (direction)
398 //direction = 1.0/direction;
399 const int * columnQuadratic = quadraticObjective_->getIndices();
400 const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
401 const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
402 const double * quadraticElement = quadraticObjective_->getElements();
403 int iColumn;
404 const double * columnScale = model->columnScale();
405 // use current linear cost region (already scaled)
406 if (includeLinear == 1) {
407 CoinMemcpyN(model->costRegion(), numberExtendedColumns_, gradient_);
408 } else if (includeLinear == 2) {
409 memset(gradient_ + numberColumns_, 0, (numberExtendedColumns_ - numberColumns_)*sizeof(double));
410 if (!columnScale) {
411 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
412 gradient_[iColumn] = objective_[iColumn] * direction;
413 }
414 } else {
415 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
416 gradient_[iColumn] = objective_[iColumn] * direction * columnScale[iColumn];
417 }
418 }
419 } else {
420 memset(gradient_, 0, numberExtendedColumns_ * sizeof(double));
421 }
422 if (!columnScale) {
423 if (activated_) {
424 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
425 double valueI = solution[iColumn];
426 CoinBigIndex j;
427 for (j = columnQuadraticStart[iColumn];
428 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
429 int jColumn = columnQuadratic[j];
430 double valueJ = solution[jColumn];
431 double elementValue = quadraticElement[j];
432 elementValue *= direction;
433 if (iColumn != jColumn) {
434 offset += valueI * valueJ * elementValue;
435 double gradientI = valueJ * elementValue;
436 double gradientJ = valueI * elementValue;
437 gradient_[iColumn] += gradientI;
438 gradient_[jColumn] += gradientJ;
439 } else {
440 offset += 0.5 * valueI * valueI * elementValue;
441 double gradientI = valueI * elementValue;
442 gradient_[iColumn] += gradientI;
443 }
444 }
445 }
446 }
447 } else {
448 // scaling
449 if (activated_) {
450 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
451 double valueI = solution[iColumn];
452 double scaleI = columnScale[iColumn] * direction;
453 CoinBigIndex j;
454 for (j = columnQuadraticStart[iColumn];
455 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
456 int jColumn = columnQuadratic[j];
457 double valueJ = solution[jColumn];
458 double elementValue = quadraticElement[j];
459 double scaleJ = columnScale[jColumn];
460 elementValue *= scaleI * scaleJ;
461 if (iColumn != jColumn) {
462 offset += valueI * valueJ * elementValue;
463 double gradientI = valueJ * elementValue;
464 double gradientJ = valueI * elementValue;
465 gradient_[iColumn] += gradientI;
466 gradient_[jColumn] += gradientJ;
467 } else {
468 offset += 0.5 * valueI * valueI * elementValue;
469 double gradientI = valueI * elementValue;
470 gradient_[iColumn] += gradientI;
471 }
472 }
473 }
474 }
475 }
476 }
477 if (model)
478 offset *= model->optimizationDirection();
479 return gradient_;
480 }
481}
482
483//-------------------------------------------------------------------
484// Clone
485//-------------------------------------------------------------------
486ClpObjective * ClpQuadraticObjective::clone() const
487{
488 return new ClpQuadraticObjective(*this);
489}
490/* Subset clone. Duplicates are allowed
491 and order is as given.
492*/
493ClpObjective *
494ClpQuadraticObjective::subsetClone (int numberColumns,
495 const int * whichColumns) const
496{
497 return new ClpQuadraticObjective(*this, numberColumns, whichColumns);
498}
499// Resize objective
500void
501ClpQuadraticObjective::resize(int newNumberColumns)
502{
503 if (numberColumns_ != newNumberColumns) {
504 int newExtended = newNumberColumns + (numberExtendedColumns_ - numberColumns_);
505 int i;
506 double * newArray = new double[newExtended];
507 if (objective_)
508 CoinMemcpyN(objective_, CoinMin(newExtended, numberExtendedColumns_), newArray);
509 delete [] objective_;
510 objective_ = newArray;
511 for (i = numberColumns_; i < newNumberColumns; i++)
512 objective_[i] = 0.0;
513 if (gradient_) {
514 newArray = new double[newExtended];
515 if (gradient_)
516 CoinMemcpyN(gradient_, CoinMin(newExtended, numberExtendedColumns_), newArray);
517 delete [] gradient_;
518 gradient_ = newArray;
519 for (i = numberColumns_; i < newNumberColumns; i++)
520 gradient_[i] = 0.0;
521 }
522 if (quadraticObjective_) {
523 if (newNumberColumns < numberColumns_) {
524 int * which = new int[numberColumns_-newNumberColumns];
525 int i;
526 for (i = newNumberColumns; i < numberColumns_; i++)
527 which[i-newNumberColumns] = i;
528 quadraticObjective_->deleteRows(numberColumns_ - newNumberColumns, which);
529 quadraticObjective_->deleteCols(numberColumns_ - newNumberColumns, which);
530 delete [] which;
531 } else {
532 quadraticObjective_->setDimensions(newNumberColumns, newNumberColumns);
533 }
534 }
535 numberColumns_ = newNumberColumns;
536 numberExtendedColumns_ = newExtended;
537 }
538
539}
540// Delete columns in objective
541void
542ClpQuadraticObjective::deleteSome(int numberToDelete, const int * which)
543{
544 int newNumberColumns = numberColumns_ - numberToDelete;
545 int newExtended = numberExtendedColumns_ - numberToDelete;
546 if (objective_) {
547 int i ;
548 char * deleted = new char[numberColumns_];
549 int numberDeleted = 0;
550 memset(deleted, 0, numberColumns_ * sizeof(char));
551 for (i = 0; i < numberToDelete; i++) {
552 int j = which[i];
553 if (j >= 0 && j < numberColumns_ && !deleted[j]) {
554 numberDeleted++;
555 deleted[j] = 1;
556 }
557 }
558 newNumberColumns = numberColumns_ - numberDeleted;
559 newExtended = numberExtendedColumns_ - numberDeleted;
560 double * newArray = new double[newExtended];
561 int put = 0;
562 for (i = 0; i < numberColumns_; i++) {
563 if (!deleted[i]) {
564 newArray[put++] = objective_[i];
565 }
566 }
567 delete [] objective_;
568 objective_ = newArray;
569 delete [] deleted;
570 CoinMemcpyN(objective_ + numberColumns_, (numberExtendedColumns_ - numberColumns_),
571 objective_ + newNumberColumns);
572 }
573 if (gradient_) {
574 int i ;
575 char * deleted = new char[numberColumns_];
576 int numberDeleted = 0;
577 memset(deleted, 0, numberColumns_ * sizeof(char));
578 for (i = 0; i < numberToDelete; i++) {
579 int j = which[i];
580 if (j >= 0 && j < numberColumns_ && !deleted[j]) {
581 numberDeleted++;
582 deleted[j] = 1;
583 }
584 }
585 newNumberColumns = numberColumns_ - numberDeleted;
586 newExtended = numberExtendedColumns_ - numberDeleted;
587 double * newArray = new double[newExtended];
588 int put = 0;
589 for (i = 0; i < numberColumns_; i++) {
590 if (!deleted[i]) {
591 newArray[put++] = gradient_[i];
592 }
593 }
594 delete [] gradient_;
595 gradient_ = newArray;
596 delete [] deleted;
597 CoinMemcpyN(gradient_ + numberColumns_, (numberExtendedColumns_ - numberColumns_),
598 gradient_ + newNumberColumns);
599 }
600 numberColumns_ = newNumberColumns;
601 numberExtendedColumns_ = newExtended;
602 if (quadraticObjective_) {
603 quadraticObjective_->deleteCols(numberToDelete, which);
604 quadraticObjective_->deleteRows(numberToDelete, which);
605 }
606}
607
608// Load up quadratic objective
609void
610ClpQuadraticObjective::loadQuadraticObjective(const int numberColumns, const CoinBigIndex * start,
611 const int * column, const double * element, int numberExtended)
612{
613 fullMatrix_ = false;
614 delete quadraticObjective_;
615 quadraticObjective_ = new CoinPackedMatrix(true, numberColumns, numberColumns,
616 start[numberColumns], element, column, start, NULL);
617 numberColumns_ = numberColumns;
618 if (numberExtended > numberExtendedColumns_) {
619 if (objective_) {
620 // make correct size
621 double * newArray = new double[numberExtended];
622 CoinMemcpyN(objective_, numberColumns_, newArray);
623 delete [] objective_;
624 objective_ = newArray;
625 memset(objective_ + numberColumns_, 0, (numberExtended - numberColumns_)*sizeof(double));
626 }
627 if (gradient_) {
628 // make correct size
629 double * newArray = new double[numberExtended];
630 CoinMemcpyN(gradient_, numberColumns_, newArray);
631 delete [] gradient_;
632 gradient_ = newArray;
633 memset(gradient_ + numberColumns_, 0, (numberExtended - numberColumns_)*sizeof(double));
634 }
635 numberExtendedColumns_ = numberExtended;
636 } else {
637 numberExtendedColumns_ = numberColumns_;
638 }
639}
640void
641ClpQuadraticObjective::loadQuadraticObjective ( const CoinPackedMatrix& matrix)
642{
643 delete quadraticObjective_;
644 quadraticObjective_ = new CoinPackedMatrix(matrix);
645}
646// Get rid of quadratic objective
647void
648ClpQuadraticObjective::deleteQuadraticObjective()
649{
650 delete quadraticObjective_;
651 quadraticObjective_ = NULL;
652}
653/* Returns reduced gradient.Returns an offset (to be added to current one).
654 */
655double
656ClpQuadraticObjective::reducedGradient(ClpSimplex * model, double * region,
657 bool useFeasibleCosts)
658{
659 int numberRows = model->numberRows();
660 int numberColumns = model->numberColumns();
661
662 //work space
663 CoinIndexedVector * workSpace = model->rowArray(0);
664
665 CoinIndexedVector arrayVector;
666 arrayVector.reserve(numberRows + 1);
667
668 int iRow;
669#ifdef CLP_DEBUG
670 workSpace->checkClear();
671#endif
672 double * array = arrayVector.denseVector();
673 int * index = arrayVector.getIndices();
674 int number = 0;
675 const double * costNow = gradient(model, model->solutionRegion(), offset_,
676 true, useFeasibleCosts ? 2 : 1);
677 double * cost = model->costRegion();
678 const int * pivotVariable = model->pivotVariable();
679 for (iRow = 0; iRow < numberRows; iRow++) {
680 int iPivot = pivotVariable[iRow];
681 double value;
682 if (iPivot < numberColumns)
683 value = costNow[iPivot];
684 else if (!useFeasibleCosts)
685 value = cost[iPivot];
686 else
687 value = 0.0;
688 if (value) {
689 array[iRow] = value;
690 index[number++] = iRow;
691 }
692 }
693 arrayVector.setNumElements(number);
694
695 // Btran basic costs
696 model->factorization()->updateColumnTranspose(workSpace, &arrayVector);
697 double * work = workSpace->denseVector();
698 ClpFillN(work, numberRows, 0.0);
699 // now look at dual solution
700 double * rowReducedCost = region + numberColumns;
701 double * dual = rowReducedCost;
702 const double * rowCost = cost + numberColumns;
703 for (iRow = 0; iRow < numberRows; iRow++) {
704 dual[iRow] = array[iRow];
705 }
706 double * dj = region;
707 ClpDisjointCopyN(costNow, numberColumns, dj);
708
709 model->transposeTimes(-1.0, dual, dj);
710 for (iRow = 0; iRow < numberRows; iRow++) {
711 // slack
712 double value = dual[iRow];
713 value += rowCost[iRow];
714 rowReducedCost[iRow] = value;
715 }
716 return offset_;
717}
718/* Returns step length which gives minimum of objective for
719 solution + theta * change vector up to maximum theta.
720
721 arrays are numberColumns+numberRows
722*/
723double
724ClpQuadraticObjective::stepLength(ClpSimplex * model,
725 const double * solution,
726 const double * change,
727 double maximumTheta,
728 double & currentObj,
729 double & predictedObj,
730 double & thetaObj)
731{
732 const double * cost = model->costRegion();
733 bool inSolve = true;
734 if (!cost) {
735 // not in solve
736 cost = objective_;
737 inSolve = false;
738 }
739 double delta = 0.0;
740 double linearCost = 0.0;
741 int numberRows = model->numberRows();
742 int numberColumns = model->numberColumns();
743 int numberTotal = numberColumns;
744 if (inSolve)
745 numberTotal += numberRows;
746 currentObj = 0.0;
747 thetaObj = 0.0;
748 for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
749 delta += cost[iColumn] * change[iColumn];
750 linearCost += cost[iColumn] * solution[iColumn];
751 }
752 if (!activated_ || !quadraticObjective_) {
753 currentObj = linearCost;
754 thetaObj = currentObj + delta * maximumTheta;
755 if (delta < 0.0) {
756 return maximumTheta;
757 } else {
758 COIN_DETAIL_PRINT(printf("odd linear direction %g\n", delta));
759 return 0.0;
760 }
761 }
762 assert (model);
763 bool scaling = false;
764 if ((model->rowScale() ||
765 model->objectiveScale() != 1.0 || model->optimizationDirection() != 1.0) && inSolve)
766 scaling = true;
767 const int * columnQuadratic = quadraticObjective_->getIndices();
768 const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
769 const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
770 const double * quadraticElement = quadraticObjective_->getElements();
771 double a = 0.0;
772 double b = delta;
773 double c = 0.0;
774 if (!scaling) {
775 if (!fullMatrix_) {
776 int iColumn;
777 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
778 double valueI = solution[iColumn];
779 double changeI = change[iColumn];
780 CoinBigIndex j;
781 for (j = columnQuadraticStart[iColumn];
782 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
783 int jColumn = columnQuadratic[j];
784 double valueJ = solution[jColumn];
785 double changeJ = change[jColumn];
786 double elementValue = quadraticElement[j];
787 if (iColumn != jColumn) {
788 a += changeI * changeJ * elementValue;
789 b += (changeI * valueJ + changeJ * valueI) * elementValue;
790 c += valueI * valueJ * elementValue;
791 } else {
792 a += 0.5 * changeI * changeI * elementValue;
793 b += changeI * valueI * elementValue;
794 c += 0.5 * valueI * valueI * elementValue;
795 }
796 }
797 }
798 } else {
799 // full matrix stored
800 int iColumn;
801 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
802 double valueI = solution[iColumn];
803 double changeI = change[iColumn];
804 CoinBigIndex j;
805 for (j = columnQuadraticStart[iColumn];
806 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
807 int jColumn = columnQuadratic[j];
808 double valueJ = solution[jColumn];
809 double changeJ = change[jColumn];
810 double elementValue = quadraticElement[j];
811 valueJ *= elementValue;
812 a += changeI * changeJ * elementValue;
813 b += changeI * valueJ;
814 c += valueI * valueJ;
815 }
816 }
817 a *= 0.5;
818 c *= 0.5;
819 }
820 } else {
821 // scaling
822 // for now only if half
823 assert (!fullMatrix_);
824 const double * columnScale = model->columnScale();
825 double direction = model->optimizationDirection() * model->objectiveScale();
826 // direction is actually scale out not scale in
827 if (direction)
828 direction = 1.0 / direction;
829 if (!columnScale) {
830 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
831 double valueI = solution[iColumn];
832 double changeI = change[iColumn];
833 CoinBigIndex j;
834 for (j = columnQuadraticStart[iColumn];
835 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
836 int jColumn = columnQuadratic[j];
837 double valueJ = solution[jColumn];
838 double changeJ = change[jColumn];
839 double elementValue = quadraticElement[j];
840 elementValue *= direction;
841 if (iColumn != jColumn) {
842 a += changeI * changeJ * elementValue;
843 b += (changeI * valueJ + changeJ * valueI) * elementValue;
844 c += valueI * valueJ * elementValue;
845 } else {
846 a += 0.5 * changeI * changeI * elementValue;
847 b += changeI * valueI * elementValue;
848 c += 0.5 * valueI * valueI * elementValue;
849 }
850 }
851 }
852 } else {
853 // scaling
854 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
855 double valueI = solution[iColumn];
856 double changeI = change[iColumn];
857 double scaleI = columnScale[iColumn] * direction;
858 CoinBigIndex j;
859 for (j = columnQuadraticStart[iColumn];
860 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
861 int jColumn = columnQuadratic[j];
862 double valueJ = solution[jColumn];
863 double changeJ = change[jColumn];
864 double elementValue = quadraticElement[j];
865 elementValue *= scaleI * columnScale[jColumn];
866 if (iColumn != jColumn) {
867 a += changeI * changeJ * elementValue;
868 b += (changeI * valueJ + changeJ * valueI) * elementValue;
869 c += valueI * valueJ * elementValue;
870 } else {
871 a += 0.5 * changeI * changeI * elementValue;
872 b += changeI * valueI * elementValue;
873 c += 0.5 * valueI * valueI * elementValue;
874 }
875 }
876 }
877 }
878 }
879 double theta;
880 //printf("Current cost %g\n",c+linearCost);
881 currentObj = c + linearCost;
882 thetaObj = currentObj + a * maximumTheta * maximumTheta + b * maximumTheta;
883 // minimize a*x*x + b*x + c
884 if (a <= 0.0) {
885 theta = maximumTheta;
886 } else {
887 theta = -0.5 * b / a;
888 }
889 predictedObj = currentObj + a * theta * theta + b * theta;
890 if (b > 0.0) {
891 if (model->messageHandler()->logLevel() & 32)
892 printf("a %g b %g c %g => %g\n", a, b, c, theta);
893 b = 0.0;
894 }
895 return CoinMin(theta, maximumTheta);
896}
897// Return objective value (without any ClpModel offset) (model may be NULL)
898double
899ClpQuadraticObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
900{
901 bool scaling = false;
902 if (model && (model->rowScale() ||
903 model->objectiveScale() != 1.0))
904 scaling = true;
905 const double * cost = NULL;
906 if (model)
907 cost = model->costRegion();
908 if (!cost) {
909 // not in solve
910 cost = objective_;
911 scaling = false;
912 }
913 double linearCost = 0.0;
914 int numberColumns = model->numberColumns();
915 int numberTotal = numberColumns;
916 double currentObj = 0.0;
917 for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
918 linearCost += cost[iColumn] * solution[iColumn];
919 }
920 if (!activated_ || !quadraticObjective_) {
921 currentObj = linearCost;
922 return currentObj;
923 }
924 assert (model);
925 const int * columnQuadratic = quadraticObjective_->getIndices();
926 const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
927 const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
928 const double * quadraticElement = quadraticObjective_->getElements();
929 double c = 0.0;
930 if (!scaling) {
931 if (!fullMatrix_) {
932 int iColumn;
933 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
934 double valueI = solution[iColumn];
935 CoinBigIndex j;
936 for (j = columnQuadraticStart[iColumn];
937 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
938 int jColumn = columnQuadratic[j];
939 double valueJ = solution[jColumn];
940 double elementValue = quadraticElement[j];
941 if (iColumn != jColumn) {
942 c += valueI * valueJ * elementValue;
943 } else {
944 c += 0.5 * valueI * valueI * elementValue;
945 }
946 }
947 }
948 } else {
949 // full matrix stored
950 int iColumn;
951 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
952 double valueI = solution[iColumn];
953 CoinBigIndex j;
954 for (j = columnQuadraticStart[iColumn];
955 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
956 int jColumn = columnQuadratic[j];
957 double valueJ = solution[jColumn];
958 double elementValue = quadraticElement[j];
959 valueJ *= elementValue;
960 c += valueI * valueJ;
961 }
962 }
963 c *= 0.5;
964 }
965 } else {
966 // scaling
967 // for now only if half
968 assert (!fullMatrix_);
969 const double * columnScale = model->columnScale();
970 double direction = model->objectiveScale();
971 // direction is actually scale out not scale in
972 if (direction)
973 direction = 1.0 / direction;
974 if (!columnScale) {
975 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
976 double valueI = solution[iColumn];
977 CoinBigIndex j;
978 for (j = columnQuadraticStart[iColumn];
979 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
980 int jColumn = columnQuadratic[j];
981 double valueJ = solution[jColumn];
982 double elementValue = quadraticElement[j];
983 elementValue *= direction;
984 if (iColumn != jColumn) {
985 c += valueI * valueJ * elementValue;
986 } else {
987 c += 0.5 * valueI * valueI * elementValue;
988 }
989 }
990 }
991 } else {
992 // scaling
993 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
994 double valueI = solution[iColumn];
995 double scaleI = columnScale[iColumn] * direction;
996 CoinBigIndex j;
997 for (j = columnQuadraticStart[iColumn];
998 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
999 int jColumn = columnQuadratic[j];
1000 double valueJ = solution[jColumn];
1001 double elementValue = quadraticElement[j];
1002 elementValue *= scaleI * columnScale[jColumn];
1003 if (iColumn != jColumn) {
1004 c += valueI * valueJ * elementValue;
1005 } else {
1006 c += 0.5 * valueI * valueI * elementValue;
1007 }
1008 }
1009 }
1010 }
1011 }
1012 currentObj = c + linearCost;
1013 return currentObj;
1014}
1015// Scale objective
1016void
1017ClpQuadraticObjective::reallyScale(const double * columnScale)
1018{
1019 const int * columnQuadratic = quadraticObjective_->getIndices();
1020 const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
1021 const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
1022 double * quadraticElement = quadraticObjective_->getMutableElements();
1023 for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1024 double scaleI = columnScale[iColumn];
1025 objective_[iColumn] *= scaleI;
1026 CoinBigIndex j;
1027 for (j = columnQuadraticStart[iColumn];
1028 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
1029 int jColumn = columnQuadratic[j];
1030 quadraticElement[j] *= scaleI * columnScale[jColumn];
1031 }
1032 }
1033}
1034/* Given a zeroed array sets nonlinear columns to 1.
1035 Returns number of nonlinear columns
1036*/
1037int
1038ClpQuadraticObjective::markNonlinear(char * which)
1039{
1040 int iColumn;
1041 const int * columnQuadratic = quadraticObjective_->getIndices();
1042 const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
1043 const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
1044 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1045 CoinBigIndex j;
1046 for (j = columnQuadraticStart[iColumn];
1047 j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
1048 int jColumn = columnQuadratic[j];
1049 which[jColumn] = 1;
1050 which[iColumn] = 1;
1051 }
1052 }
1053 int numberNonLinearColumns = 0;
1054 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1055 if(which[iColumn])
1056 numberNonLinearColumns++;
1057 }
1058 return numberNonLinearColumns;
1059}
1060