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

Last change on this file since 1034 was 1034, checked in by forrest, 13 years ago

moving branches/devel to trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 33.9 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4#include "CoinPragma.hpp"
5#include "CoinHelperFunctions.hpp"
6#include "CoinIndexedVector.hpp"
7#include "ClpFactorization.hpp"
8#include "ClpSimplex.hpp"
9#include "ClpQuadraticObjective.hpp"
10//#############################################################################
11// Constructors / Destructor / Assignment
12//#############################################################################
13//-------------------------------------------------------------------
14// Default Constructor
15//-------------------------------------------------------------------
16ClpQuadraticObjective::ClpQuadraticObjective () 
17: ClpObjective()
18{
19  type_=2;
20  objective_=NULL;
21  quadraticObjective_=NULL;
22  gradient_ = NULL;
23  numberColumns_=0;
24  numberExtendedColumns_=0;
25  activated_=0;
26  fullMatrix_=false;
27}
28
29//-------------------------------------------------------------------
30// Useful Constructor
31//-------------------------------------------------------------------
32ClpQuadraticObjective::ClpQuadraticObjective (const double * objective , 
33                                              int numberColumns,
34                                              const CoinBigIndex * start,
35                                              const int * column, const double * element,
36                                              int numberExtendedColumns)
37  : ClpObjective()
38{
39  type_=2;
40  numberColumns_ = numberColumns;
41  if (numberExtendedColumns>=0)
42    numberExtendedColumns_= max(numberColumns_,numberExtendedColumns);
43  else
44    numberExtendedColumns_= numberColumns_;
45  if (objective) {
46    objective_ = new double [numberExtendedColumns_];
47    memcpy(objective_,objective,numberColumns_*sizeof(double));
48    memset(objective_+numberColumns_,0,(numberExtendedColumns_-numberColumns_)*sizeof(double));
49  } else {
50    objective_ = new double [numberExtendedColumns_];
51    memset(objective_,0,numberExtendedColumns_*sizeof(double));
52  }
53  if (start) 
54    quadraticObjective_ = new CoinPackedMatrix(true,numberColumns,numberColumns,
55                                             start[numberColumns],element,column,start,NULL);
56  else
57  quadraticObjective_=NULL;
58  gradient_ = NULL;
59  activated_=1;
60  fullMatrix_=false;
61}
62
63//-------------------------------------------------------------------
64// Copy constructor
65//-------------------------------------------------------------------
66ClpQuadraticObjective::ClpQuadraticObjective (const ClpQuadraticObjective & rhs,
67                                              int type) 
68: ClpObjective(rhs)
69{ 
70  numberColumns_=rhs.numberColumns_;
71  numberExtendedColumns_=rhs.numberExtendedColumns_;
72  fullMatrix_=rhs.fullMatrix_;
73  if (rhs.objective_) {
74    objective_ = new double [numberExtendedColumns_];
75    memcpy(objective_,rhs.objective_,numberExtendedColumns_*sizeof(double));
76  } else {
77    objective_=NULL;
78  }
79  if (rhs.gradient_) {
80    gradient_ = new double [numberExtendedColumns_];
81    memcpy(gradient_,rhs.gradient_,numberExtendedColumns_*sizeof(double));
82  } else {
83    gradient_=NULL;
84  }
85  if (rhs.quadraticObjective_) {
86    // see what type of matrix wanted
87    if (type==0) {
88      // just copy
89      quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
90    } else if (type==1) {
91      // expand to full symmetric
92      fullMatrix_=true;
93      const int * columnQuadratic1 = rhs.quadraticObjective_->getIndices();
94      const CoinBigIndex * columnQuadraticStart1 = rhs.quadraticObjective_->getVectorStarts();
95      const int * columnQuadraticLength1 = rhs.quadraticObjective_->getVectorLengths();
96      const double * quadraticElement1 = rhs.quadraticObjective_->getElements();
97      CoinBigIndex * columnQuadraticStart2 = new CoinBigIndex [numberExtendedColumns_+1];
98      int * columnQuadraticLength2 = new int [numberExtendedColumns_];
99      int iColumn;
100      int numberColumns = rhs.quadraticObjective_->getNumCols();
101      int numberBelow=0;
102      int numberAbove=0;
103      int numberDiagonal=0;
104      CoinZeroN(columnQuadraticLength2,numberExtendedColumns_);
105      for (iColumn=0;iColumn<numberColumns;iColumn++) {
106        for (CoinBigIndex j=columnQuadraticStart1[iColumn];
107             j<columnQuadraticStart1[iColumn]+columnQuadraticLength1[iColumn];j++) {
108          int jColumn = columnQuadratic1[j];
109          if (jColumn>iColumn) {
110            numberBelow++;
111            columnQuadraticLength2[jColumn]++;
112            columnQuadraticLength2[iColumn]++;
113          } else if (jColumn==iColumn) {
114            numberDiagonal++;
115            columnQuadraticLength2[iColumn]++;
116          } else {
117            numberAbove++;
118          }
119        }
120      }
121      if (numberAbove>0) {
122        if (numberAbove==numberBelow) {
123          // already done
124          quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
125          delete [] columnQuadraticStart2;
126          delete [] columnQuadraticLength2;
127        } else {
128          printf("number above = %d, number below = %d, error\n",
129                 numberAbove,numberBelow);
130          abort();
131        }
132      } else {
133        int numberElements=numberDiagonal+2*numberBelow;
134        int * columnQuadratic2 = new int [numberElements];
135        double * quadraticElement2 = new double [numberElements];
136        columnQuadraticStart2[0]=0;
137        numberElements=0;
138        for (iColumn=0;iColumn<numberColumns;iColumn++) {
139          int n=columnQuadraticLength2[iColumn];
140          columnQuadraticLength2[iColumn]=0;
141          numberElements += n;
142          columnQuadraticStart2[iColumn+1]=numberElements;
143        }
144        for (iColumn=0;iColumn<numberColumns;iColumn++) {
145          for (CoinBigIndex j=columnQuadraticStart1[iColumn];
146               j<columnQuadraticStart1[iColumn]+columnQuadraticLength1[iColumn];j++) {
147            int jColumn = columnQuadratic1[j];
148            if (jColumn>iColumn) {
149              // put in two places
150              CoinBigIndex put=columnQuadraticLength2[jColumn]+columnQuadraticStart2[jColumn];
151              columnQuadraticLength2[jColumn]++;
152              quadraticElement2[put]=quadraticElement1[j];
153              columnQuadratic2[put]=iColumn;
154              put=columnQuadraticLength2[iColumn]+columnQuadraticStart2[iColumn];
155              columnQuadraticLength2[iColumn]++;
156              quadraticElement2[put]=quadraticElement1[j];
157              columnQuadratic2[put]=jColumn;
158            } else if (jColumn==iColumn) {
159              CoinBigIndex put=columnQuadraticLength2[iColumn]+columnQuadraticStart2[iColumn];
160              columnQuadraticLength2[iColumn]++;
161              quadraticElement2[put]=quadraticElement1[j];
162              columnQuadratic2[put]=iColumn;
163            } else {
164              abort();
165            }
166          }
167        }
168        // Now create
169        quadraticObjective_ = 
170          new CoinPackedMatrix (true,
171                                rhs.numberExtendedColumns_,
172                                rhs.numberExtendedColumns_,
173                                numberElements,
174                                quadraticElement2,
175                                columnQuadratic2,
176                                columnQuadraticStart2,
177                                columnQuadraticLength2,0.0,0.0);
178        delete [] columnQuadraticStart2;
179        delete [] columnQuadraticLength2;
180        delete [] columnQuadratic2;
181        delete [] quadraticElement2;
182      }
183    } else {
184      fullMatrix_=false;
185      abort(); // code when needed
186    }
187           
188  } else {
189    quadraticObjective_=NULL;
190  }
191}
192/* Subset constructor.  Duplicates are allowed
193   and order is as given.
194*/
195ClpQuadraticObjective::ClpQuadraticObjective (const ClpQuadraticObjective &rhs,
196                                        int numberColumns, 
197                                        const int * whichColumn) 
198: ClpObjective(rhs)
199{
200  fullMatrix_=rhs.fullMatrix_;
201  objective_=NULL;
202  int extra = rhs.numberExtendedColumns_-rhs.numberColumns_;
203  numberColumns_=0;
204  numberExtendedColumns_ = numberColumns+extra;
205  if (numberColumns>0) {
206    // check valid lists
207    int numberBad=0;
208    int i;
209    for (i=0;i<numberColumns;i++)
210      if (whichColumn[i]<0||whichColumn[i]>=rhs.numberColumns_)
211        numberBad++;
212    if (numberBad)
213      throw CoinError("bad column list", "subset constructor", 
214                      "ClpQuadraticObjective");
215    numberColumns_ = numberColumns;
216    objective_ = new double[numberExtendedColumns_];
217    for (i=0;i<numberColumns_;i++) 
218      objective_[i]=rhs.objective_[whichColumn[i]];
219    memcpy(objective_+numberColumns_,rhs.objective_+rhs.numberColumns_,
220           (numberExtendedColumns_-numberColumns_)*sizeof(double));
221    if (rhs.gradient_) {
222      gradient_ = new double[numberExtendedColumns_];
223      for (i=0;i<numberColumns_;i++) 
224        gradient_[i]=rhs.gradient_[whichColumn[i]];
225      memcpy(gradient_+numberColumns_,rhs.gradient_+rhs.numberColumns_,
226             (numberExtendedColumns_-numberColumns_)*sizeof(double));
227    } else {
228      gradient_=NULL;
229    }
230  } else {
231    gradient_=NULL;
232    objective_=NULL;
233  }
234  if (rhs.quadraticObjective_) {
235    quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_,
236                                               numberColumns,whichColumn,
237                                               numberColumns,whichColumn);
238  } else {
239    quadraticObjective_=NULL;
240  }
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      memcpy(objective_,rhs.objective_,numberExtendedColumns_*sizeof(double));
272    } else {
273      objective_=NULL;
274    }
275    if (rhs.gradient_) {
276      gradient_ = new double [numberExtendedColumns_];
277      memcpy(gradient_,rhs.gradient_,numberExtendedColumns_*sizeof(double));
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()||
299              model->objectiveScale()!=1.0||model->optimizationDirection()!=1.0))
300    scaling=true;
301  const double * cost = NULL;
302  if (model)
303    cost = model->costRegion();
304  if (!cost) {
305    // not in solve
306    cost = objective_;
307    scaling=false;
308  }
309  if (!scaling) {
310    if (!quadraticObjective_||!solution||!activated_) {
311      return objective_;
312    } else {
313      if (refresh||!gradient_) {
314        if (!gradient_) 
315          gradient_ = new double[numberExtendedColumns_];
316        const int * columnQuadratic = quadraticObjective_->getIndices();
317        const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
318        const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
319        const double * quadraticElement = quadraticObjective_->getElements();
320        offset=0.0;
321        // use current linear cost region
322        if (includeLinear==1)
323          memcpy(gradient_,cost,numberExtendedColumns_*sizeof(double));
324        else if (includeLinear==2)
325          memcpy(gradient_,objective_,numberExtendedColumns_*sizeof(double));
326        else
327          memset(gradient_,0,numberExtendedColumns_*sizeof(double));
328        if (activated_) {
329          if (!fullMatrix_) {
330            int iColumn;
331            for (iColumn=0;iColumn<numberColumns_;iColumn++) {
332              double valueI = solution[iColumn];
333              CoinBigIndex j;
334              for (j=columnQuadraticStart[iColumn];
335                   j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
336                int jColumn = columnQuadratic[j];
337                double valueJ = solution[jColumn];
338                double elementValue = quadraticElement[j];
339                if (iColumn!=jColumn) {
340                  offset += valueI*valueJ*elementValue;
341                  //if (fabs(valueI*valueJ*elementValue)>1.0e-12)
342                  //printf("%d %d %g %g %g -> %g\n",
343                  //       iColumn,jColumn,valueI,valueJ,elementValue,
344                  //       valueI*valueJ*elementValue);
345                  double gradientI = valueJ*elementValue;
346                  double gradientJ = valueI*elementValue;
347                  gradient_[iColumn] += gradientI;
348                  gradient_[jColumn] += gradientJ;
349                } else {
350                  offset += 0.5*valueI*valueI*elementValue;
351                  //if (fabs(valueI*valueI*elementValue)>1.0e-12)
352                  //printf("XX %d %g %g -> %g\n",
353                  //       iColumn,valueI,elementValue,
354                  //       0.5*valueI*valueI*elementValue);
355                  double gradientI = valueI*elementValue;
356                  gradient_[iColumn] += gradientI;
357                }
358              }
359            }
360          } else {
361            // full matrix
362            int iColumn;
363            offset *= 2.0;
364            for (iColumn=0;iColumn<numberColumns_;iColumn++) {
365              CoinBigIndex j;
366              double value=0.0;
367              double current = gradient_[iColumn];
368              for (j=columnQuadraticStart[iColumn];
369                   j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
370                int jColumn = columnQuadratic[j];
371                double valueJ = solution[jColumn]*quadraticElement[j];
372                value += valueJ;
373              }
374              offset += value*solution[iColumn];
375              gradient_[iColumn] = current+value;
376            }
377            offset *= 0.5;
378          }
379        }
380      }
381      if (model)
382        offset *= model->optimizationDirection()*model->objectiveScale();
383      return gradient_;
384    }
385  } else {
386    // do scaling
387    assert (solution);
388    // for now only if half
389    assert (!fullMatrix_);
390    if (refresh||!gradient_) {
391      if (!gradient_) 
392        gradient_ = new double[numberExtendedColumns_];
393      double direction = model->optimizationDirection()*model->objectiveScale();
394      // direction is actually scale out not scale in
395      //if (direction)
396      //direction = 1.0/direction;
397      const int * columnQuadratic = quadraticObjective_->getIndices();
398      const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
399      const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
400      const double * quadraticElement = quadraticObjective_->getElements();
401      int iColumn;
402      const double * columnScale = model->columnScale();
403      // use current linear cost region (already scaled)
404      if (includeLinear==1) {
405        memcpy(gradient_,model->costRegion(),numberExtendedColumns_*sizeof(double));
406      } else if (includeLinear==2) {
407        memset(gradient_+numberColumns_,0,(numberExtendedColumns_-numberColumns_)*sizeof(double));
408        if (!columnScale) {
409          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
410            gradient_[iColumn]= objective_[iColumn]*direction;
411          }
412        } else {
413          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
414            gradient_[iColumn]= objective_[iColumn]*direction*columnScale[iColumn];
415          }
416        }
417      } else {
418        memset(gradient_,0,numberExtendedColumns_*sizeof(double));
419      }
420      if (!columnScale) {
421        if (activated_) {
422          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
423            double valueI = solution[iColumn];
424            CoinBigIndex j;
425            for (j=columnQuadraticStart[iColumn];
426                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
427              int jColumn = columnQuadratic[j];
428              double valueJ = solution[jColumn];
429              double elementValue = quadraticElement[j];
430              elementValue *= direction;
431              if (iColumn!=jColumn) {
432                offset += valueI*valueJ*elementValue;
433                double gradientI = valueJ*elementValue;
434                double gradientJ = valueI*elementValue;
435                gradient_[iColumn] += gradientI;
436                gradient_[jColumn] += gradientJ;
437              } else {
438                offset += 0.5*valueI*valueI*elementValue;
439                double gradientI = valueI*elementValue;
440                gradient_[iColumn] += gradientI;
441              }
442            }
443          }
444        }
445      } else {
446        // scaling
447        if (activated_) {
448          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
449            double valueI = solution[iColumn];
450            double scaleI = columnScale[iColumn]*direction;
451            CoinBigIndex j;
452            for (j=columnQuadraticStart[iColumn];
453                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
454              int jColumn = columnQuadratic[j];
455              double valueJ = solution[jColumn];
456              double elementValue = quadraticElement[j];
457              double scaleJ = columnScale[jColumn];
458              elementValue *= scaleI*scaleJ;
459              if (iColumn!=jColumn) {
460                offset += valueI*valueJ*elementValue;
461                double gradientI = valueJ*elementValue;
462                double gradientJ = valueI*elementValue;
463                gradient_[iColumn] += gradientI;
464                gradient_[jColumn] += gradientJ;
465              } else {
466                offset += 0.5*valueI*valueI*elementValue;
467                double gradientI = valueI*elementValue;
468                gradient_[iColumn] += gradientI;
469              }
470            }
471          }
472        }
473      }
474    }
475    if (model)
476      offset *= model->optimizationDirection();
477    return gradient_;
478  }
479}
480 
481//-------------------------------------------------------------------
482// Clone
483//-------------------------------------------------------------------
484ClpObjective * ClpQuadraticObjective::clone() const
485{
486  return new ClpQuadraticObjective(*this);
487}
488/* Subset clone.  Duplicates are allowed
489   and order is as given.
490*/
491ClpObjective * 
492ClpQuadraticObjective::subsetClone (int numberColumns, 
493                           const int * whichColumns) const
494{
495  return new ClpQuadraticObjective(*this, numberColumns, whichColumns);
496}
497// Resize objective
498void 
499ClpQuadraticObjective::resize(int newNumberColumns)
500{
501  if (numberColumns_!=newNumberColumns) {
502    int newExtended = newNumberColumns + (numberExtendedColumns_-numberColumns_);
503    int i;
504    double * newArray = new double[newExtended];
505    if (objective_)
506      memcpy(newArray,objective_,
507             CoinMin(newExtended,numberExtendedColumns_)*sizeof(double));
508    delete [] objective_;
509    objective_ = newArray;
510    for (i=numberColumns_;i<newNumberColumns;i++) 
511      objective_[i]=0.0;
512    if (gradient_) {
513      newArray = new double[newExtended];
514      if (gradient_)
515        memcpy(newArray,gradient_,
516               CoinMin(newExtended,numberExtendedColumns_)*sizeof(double));
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    memcpy(objective_+newNumberColumns,objective_+numberColumns_,
571           (numberExtendedColumns_-numberColumns_)*sizeof(double));
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    memcpy(gradient_+newNumberColumns,gradient_+numberColumns_,
598           (numberExtendedColumns_-numberColumns_)*sizeof(double));
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      memcpy(newArray,objective_,numberColumns_*sizeof(double));
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      memcpy(newArray,gradient_,numberColumns_*sizeof(double));
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      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}
Note: See TracBrowser for help on using the repository browser.