source: trunk/Clp/src/ClpQuadraticObjective.cpp @ 2385

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago

formatting

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