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

Last change on this file since 1197 was 1197, checked in by forrest, 12 years ago

many changes to try and improve performance

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 33.8 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    CoinMemcpyN(objective,numberColumns_,objective_);
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    CoinMemcpyN(rhs.objective_,numberExtendedColumns_,objective_);
76  } else {
77    objective_=NULL;
78  }
79  if (rhs.gradient_) {
80    gradient_ = new double [numberExtendedColumns_];
81    CoinMemcpyN(rhs.gradient_,numberExtendedColumns_,gradient_);
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    CoinMemcpyN(rhs.objective_+rhs.numberColumns_,      (numberExtendedColumns_-numberColumns_),
220                 objective_+numberColumns_);
221    if (rhs.gradient_) {
222      gradient_ = new double[numberExtendedColumns_];
223      for (i=0;i<numberColumns_;i++) 
224        gradient_[i]=rhs.gradient_[whichColumn[i]];
225      CoinMemcpyN(rhs.gradient_+rhs.numberColumns_,     (numberExtendedColumns_-numberColumns_),
226                   gradient_+numberColumns_);
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      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()||
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   CoinMemcpyN(cost,numberExtendedColumns_,gradient_);
324        else if (includeLinear==2)
325   CoinMemcpyN(objective_,numberExtendedColumns_,gradient_);
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 CoinMemcpyN(model->costRegion(),numberExtendedColumns_,gradient_);
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      CoinMemcpyN(objective_,   CoinMin(newExtended,numberExtendedColumns_),newArray);
507    delete [] objective_;
508    objective_ = newArray;
509    for (i=numberColumns_;i<newNumberColumns;i++) 
510      objective_[i]=0.0;
511    if (gradient_) {
512      newArray = new double[newExtended];
513      if (gradient_)
514 CoinMemcpyN(gradient_, CoinMin(newExtended,numberExtendedColumns_),newArray);
515      delete [] gradient_;
516      gradient_ = newArray;
517      for (i=numberColumns_;i<newNumberColumns;i++) 
518        gradient_[i]=0.0;
519    }
520    if (quadraticObjective_) {
521      if (newNumberColumns<numberColumns_) {
522        int * which = new int[numberColumns_-newNumberColumns];
523        int i;
524        for (i=newNumberColumns;i<numberColumns_;i++) 
525          which[i-newNumberColumns]=i;
526        quadraticObjective_->deleteRows(numberColumns_-newNumberColumns,which);
527        quadraticObjective_->deleteCols(numberColumns_-newNumberColumns,which);
528        delete [] which;
529      } else {
530        quadraticObjective_->setDimensions(newNumberColumns,newNumberColumns);
531      }
532    }
533    numberColumns_ = newNumberColumns;
534    numberExtendedColumns_ = newExtended;
535  } 
536 
537}
538// Delete columns in  objective
539void 
540ClpQuadraticObjective::deleteSome(int numberToDelete, const int * which) 
541{
542  int newNumberColumns = numberColumns_-numberToDelete;
543  int newExtended = numberExtendedColumns_ - numberToDelete;
544  if (objective_) {
545    int i ;
546    char * deleted = new char[numberColumns_];
547    int numberDeleted=0;
548    memset(deleted,0,numberColumns_*sizeof(char));
549    for (i=0;i<numberToDelete;i++) {
550      int j = which[i];
551      if (j>=0&&j<numberColumns_&&!deleted[j]) {
552        numberDeleted++;
553        deleted[j]=1;
554      }
555    }
556    newNumberColumns = numberColumns_-numberDeleted;
557    newExtended = numberExtendedColumns_ - numberDeleted;
558    double * newArray = new double[newExtended];
559    int put=0;
560    for (i=0;i<numberColumns_;i++) {
561      if (!deleted[i]) {
562        newArray[put++]=objective_[i];
563      }
564    }
565    delete [] objective_;
566    objective_ = newArray;
567    delete [] deleted;
568    CoinMemcpyN(objective_+numberColumns_,      (numberExtendedColumns_-numberColumns_),
569                 objective_+newNumberColumns);
570  }
571  if (gradient_) {
572    int i ;
573    char * deleted = new char[numberColumns_];
574    int numberDeleted=0;
575    memset(deleted,0,numberColumns_*sizeof(char));
576    for (i=0;i<numberToDelete;i++) {
577      int j = which[i];
578      if (j>=0&&j<numberColumns_&&!deleted[j]) {
579        numberDeleted++;
580        deleted[j]=1;
581      }
582    }
583    newNumberColumns = numberColumns_-numberDeleted;
584    newExtended = numberExtendedColumns_ - numberDeleted;
585    double * newArray = new double[newExtended];
586    int put=0;
587    for (i=0;i<numberColumns_;i++) {
588      if (!deleted[i]) {
589        newArray[put++]=gradient_[i];
590      }
591    }
592    delete [] gradient_;
593    gradient_ = newArray;
594    delete [] deleted;
595    CoinMemcpyN(gradient_+numberColumns_,       (numberExtendedColumns_-numberColumns_),
596                 gradient_+newNumberColumns);
597  }
598  numberColumns_ = newNumberColumns;
599  numberExtendedColumns_ = newExtended;
600  if (quadraticObjective_) {
601    quadraticObjective_->deleteCols(numberToDelete,which);
602    quadraticObjective_->deleteRows(numberToDelete,which);
603  }
604}
605
606// Load up quadratic objective
607void 
608ClpQuadraticObjective::loadQuadraticObjective(const int numberColumns, const CoinBigIndex * start,
609                              const int * column, const double * element,int numberExtended)
610{
611  fullMatrix_=false;
612  delete quadraticObjective_;
613  quadraticObjective_ = new CoinPackedMatrix(true,numberColumns,numberColumns,
614                                             start[numberColumns],element,column,start,NULL);
615  numberColumns_=numberColumns;
616  if (numberExtended>numberExtendedColumns_) {
617    if (objective_) {
618      // make correct size
619      double * newArray = new double[numberExtended];
620      CoinMemcpyN(objective_,numberColumns_,newArray);
621      delete [] objective_;
622      objective_ = newArray;
623      memset(objective_+numberColumns_,0,(numberExtended-numberColumns_)*sizeof(double));
624    }
625    if (gradient_) {
626      // make correct size
627      double * newArray = new double[numberExtended];
628      CoinMemcpyN(gradient_,numberColumns_,newArray);
629      delete [] gradient_;
630      gradient_ = newArray;
631      memset(gradient_+numberColumns_,0,(numberExtended-numberColumns_)*sizeof(double));
632    }
633    numberExtendedColumns_ = numberExtended;
634  } else {
635    numberExtendedColumns_ = numberColumns_;
636  }
637}
638void 
639ClpQuadraticObjective::loadQuadraticObjective (  const CoinPackedMatrix& matrix)
640{
641  delete quadraticObjective_;
642  quadraticObjective_ = new CoinPackedMatrix(matrix);
643}
644// Get rid of quadratic objective
645void 
646ClpQuadraticObjective::deleteQuadraticObjective()
647{
648  delete quadraticObjective_;
649  quadraticObjective_ = NULL;
650}
651/* Returns reduced gradient.Returns an offset (to be added to current one).
652 */
653double 
654ClpQuadraticObjective::reducedGradient(ClpSimplex * model, double * region,
655                                       bool useFeasibleCosts)
656{
657  int numberRows = model->numberRows();
658  int numberColumns=model->numberColumns();
659 
660  //work space
661  CoinIndexedVector  * workSpace = model->rowArray(0);
662 
663  CoinIndexedVector arrayVector;
664  arrayVector.reserve(numberRows+1);
665 
666  int iRow;
667#ifdef CLP_DEBUG
668  workSpace->checkClear();
669#endif
670  double * array = arrayVector.denseVector();
671  int * index = arrayVector.getIndices();
672  int number=0;
673  const double * costNow = gradient(model,model->solutionRegion(),offset_,
674                                    true,useFeasibleCosts ? 2 : 1);
675  double * cost = model->costRegion();
676  const int * pivotVariable = model->pivotVariable();
677  for (iRow=0;iRow<numberRows;iRow++) {
678    int iPivot=pivotVariable[iRow];
679    double value;
680    if (iPivot<numberColumns)
681      value = costNow[iPivot];
682    else if (!useFeasibleCosts) 
683      value = cost[iPivot];
684    else 
685      value=0.0;
686    if (value) {
687      array[iRow]=value;
688      index[number++]=iRow;
689    }
690  }
691  arrayVector.setNumElements(number);
692
693  // Btran basic costs
694  model->factorization()->updateColumnTranspose(workSpace,&arrayVector);
695  double * work = workSpace->denseVector();
696  ClpFillN(work,numberRows,0.0);
697  // now look at dual solution
698  double * rowReducedCost = region+numberColumns;
699  double * dual = rowReducedCost;
700  const double * rowCost = cost+numberColumns;
701  for (iRow=0;iRow<numberRows;iRow++) {
702    dual[iRow]=array[iRow];
703  }
704  double * dj = region;
705  ClpDisjointCopyN(costNow,numberColumns,dj);
706 
707  model->transposeTimes(-1.0,dual,dj);
708  for (iRow=0;iRow<numberRows;iRow++) {
709    // slack
710    double value = dual[iRow];
711    value += rowCost[iRow];
712    rowReducedCost[iRow]=value;
713  }
714  return offset_;
715}
716/* Returns step length which gives minimum of objective for
717   solution + theta * change vector up to maximum theta.
718   
719   arrays are numberColumns+numberRows
720*/
721double 
722ClpQuadraticObjective::stepLength(ClpSimplex * model,
723                                  const double * solution,
724                                  const double * change,
725                                  double maximumTheta,
726                                  double & currentObj,
727                                  double & predictedObj,
728                                  double & thetaObj)
729{
730  const double * cost = model->costRegion();
731  bool inSolve=true;
732  if (!cost) {
733    // not in solve
734    cost = objective_;
735    inSolve=false;
736  }
737  double delta=0.0;
738  double linearCost =0.0;
739  int numberRows = model->numberRows();
740  int numberColumns = model->numberColumns();
741  int numberTotal = numberColumns;
742  if (inSolve)
743    numberTotal += numberRows;
744  currentObj=0.0;
745  thetaObj=0.0;
746  for (int iColumn=0;iColumn<numberTotal;iColumn++) {
747    delta += cost[iColumn]*change[iColumn];
748    linearCost += cost[iColumn]*solution[iColumn];
749  }
750  if (!activated_||!quadraticObjective_) {
751    currentObj=linearCost;
752    thetaObj =currentObj + delta*maximumTheta;
753    if (delta<0.0) {
754      return maximumTheta;
755    } else {
756      printf("odd linear direction %g\n",delta);
757      return 0.0;
758    }
759  }
760  assert (model);
761  bool scaling=false;
762  if ((model->rowScale()||
763       model->objectiveScale()!=1.0||model->optimizationDirection()!=1.0)&&inSolve)
764    scaling=true;
765  const int * columnQuadratic = quadraticObjective_->getIndices();
766  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
767  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
768  const double * quadraticElement = quadraticObjective_->getElements();
769  double a=0.0;
770  double b=delta;
771  double c=0.0;
772  if (!scaling) {
773    if (!fullMatrix_) {
774      int iColumn;
775      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
776        double valueI = solution[iColumn];
777        double changeI = change[iColumn];
778        CoinBigIndex j;
779        for (j=columnQuadraticStart[iColumn];
780             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
781          int jColumn = columnQuadratic[j];
782          double valueJ = solution[jColumn];
783          double changeJ = change[jColumn];
784          double elementValue = quadraticElement[j];
785          if (iColumn!=jColumn) {
786            a += changeI*changeJ*elementValue;
787            b += (changeI*valueJ+changeJ*valueI)*elementValue;
788            c += valueI*valueJ*elementValue;
789          } else {
790            a += 0.5*changeI*changeI*elementValue;
791            b += changeI*valueI*elementValue;
792            c += 0.5*valueI*valueI*elementValue;
793          }
794        }
795      }
796    } else {
797      // full matrix stored
798      int iColumn;
799      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
800        double valueI = solution[iColumn];
801        double changeI = change[iColumn];
802        CoinBigIndex j;
803        for (j=columnQuadraticStart[iColumn];
804             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
805          int jColumn = columnQuadratic[j];
806          double valueJ = solution[jColumn];
807          double changeJ = change[jColumn];
808          double elementValue = quadraticElement[j];
809          valueJ *= elementValue;
810          a += changeI*changeJ*elementValue;
811          b += changeI*valueJ;
812          c += valueI*valueJ;
813        }
814      }
815      a *= 0.5;
816      c *= 0.5;
817    }
818  } else {
819    // scaling
820    // for now only if half
821    assert (!fullMatrix_);
822    const double * columnScale = model->columnScale();
823    double direction = model->optimizationDirection()*model->objectiveScale();
824    // direction is actually scale out not scale in
825    if (direction)
826      direction = 1.0/direction;
827    if (!columnScale) {
828      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
829        double valueI = solution[iColumn];
830        double changeI = change[iColumn];
831        CoinBigIndex j;
832        for (j=columnQuadraticStart[iColumn];
833             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
834          int jColumn = columnQuadratic[j];
835          double valueJ = solution[jColumn];
836          double changeJ = change[jColumn];
837          double elementValue = quadraticElement[j];
838          elementValue *= direction;
839          if (iColumn!=jColumn) {
840            a += changeI*changeJ*elementValue;
841            b += (changeI*valueJ+changeJ*valueI)*elementValue;
842            c += valueI*valueJ*elementValue;
843          } else {
844            a += 0.5*changeI*changeI*elementValue;
845            b += changeI*valueI*elementValue;
846            c += 0.5*valueI*valueI*elementValue;
847          }
848        }
849      }
850    } else {
851      // scaling
852      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
853        double valueI = solution[iColumn];
854        double changeI = change[iColumn];
855        double scaleI = columnScale[iColumn]*direction;
856        CoinBigIndex j;
857        for (j=columnQuadraticStart[iColumn];
858             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
859          int jColumn = columnQuadratic[j];
860          double valueJ = solution[jColumn];
861          double changeJ = change[jColumn];
862          double elementValue = quadraticElement[j];
863          elementValue *= scaleI*columnScale[jColumn];
864          if (iColumn!=jColumn) {
865            a += changeI*changeJ*elementValue;
866            b += (changeI*valueJ+changeJ*valueI)*elementValue;
867            c += valueI*valueJ*elementValue;
868          } else {
869            a += 0.5*changeI*changeI*elementValue;
870            b += changeI*valueI*elementValue;
871            c += 0.5*valueI*valueI*elementValue;
872          }
873        }
874      }
875    }
876  }
877  double theta;
878  //printf("Current cost %g\n",c+linearCost);
879  currentObj = c+linearCost;
880  thetaObj = currentObj + a*maximumTheta*maximumTheta+b*maximumTheta;
881  // minimize a*x*x + b*x + c
882  if (a<=0.0) {
883    theta = maximumTheta;
884  } else {
885    theta = -0.5*b/a;
886  }
887  predictedObj = currentObj + a*theta*theta+b*theta;
888  if (b>0.0) {
889    if (model->messageHandler()->logLevel()&32)
890      printf("a %g b %g c %g => %g\n",a,b,c,theta); 
891    b=0.0;
892  }
893  return CoinMin(theta,maximumTheta);
894}
895// Return objective value (without any ClpModel offset) (model may be NULL)
896double 
897ClpQuadraticObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
898{
899  bool scaling=false;
900  if (model&&(model->rowScale()||
901              model->objectiveScale()!=1.0))
902    scaling=true;
903  const double * cost = NULL;
904  if (model)
905    cost = model->costRegion();
906  if (!cost) {
907    // not in solve
908    cost = objective_;
909    scaling=false;
910  }
911  double linearCost =0.0;
912  int numberColumns = model->numberColumns();
913  int numberTotal = numberColumns;
914  double currentObj=0.0;
915  for (int iColumn=0;iColumn<numberTotal;iColumn++) {
916    linearCost += cost[iColumn]*solution[iColumn];
917  }
918  if (!activated_||!quadraticObjective_) {
919    currentObj=linearCost;
920    return currentObj;
921  }
922  assert (model);
923  const int * columnQuadratic = quadraticObjective_->getIndices();
924  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
925  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
926  const double * quadraticElement = quadraticObjective_->getElements();
927  double c=0.0;
928  if (!scaling) {
929    if (!fullMatrix_) {
930      int iColumn;
931      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
932        double valueI = solution[iColumn];
933        CoinBigIndex j;
934        for (j=columnQuadraticStart[iColumn];
935             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
936          int jColumn = columnQuadratic[j];
937          double valueJ = solution[jColumn];
938          double elementValue = quadraticElement[j];
939          if (iColumn!=jColumn) {
940            c += valueI*valueJ*elementValue;
941          } else {
942            c += 0.5*valueI*valueI*elementValue;
943          }
944        }
945      }
946    } else {
947      // full matrix stored
948      int iColumn;
949      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
950        double valueI = solution[iColumn];
951        CoinBigIndex j;
952        for (j=columnQuadraticStart[iColumn];
953             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
954          int jColumn = columnQuadratic[j];
955          double valueJ = solution[jColumn];
956          double elementValue = quadraticElement[j];
957          valueJ *= elementValue;
958          c += valueI*valueJ;
959        }
960      }
961      c *= 0.5;
962    }
963  } else {
964    // scaling
965    // for now only if half
966    assert (!fullMatrix_);
967    const double * columnScale = model->columnScale();
968    double direction = model->objectiveScale();
969    // direction is actually scale out not scale in
970    if (direction)
971      direction = 1.0/direction;
972    if (!columnScale) {
973      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
974        double valueI = solution[iColumn];
975        CoinBigIndex j;
976        for (j=columnQuadraticStart[iColumn];
977             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
978          int jColumn = columnQuadratic[j];
979          double valueJ = solution[jColumn];
980          double elementValue = quadraticElement[j];
981          elementValue *= direction;
982          if (iColumn!=jColumn) {
983            c += valueI*valueJ*elementValue;
984          } else {
985            c += 0.5*valueI*valueI*elementValue;
986          }
987        }
988      }
989    } else {
990      // scaling
991      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
992        double valueI = solution[iColumn];
993        double scaleI = columnScale[iColumn]*direction;
994        CoinBigIndex j;
995        for (j=columnQuadraticStart[iColumn];
996             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
997          int jColumn = columnQuadratic[j];
998          double valueJ = solution[jColumn];
999          double elementValue = quadraticElement[j];
1000          elementValue *= scaleI*columnScale[jColumn];
1001          if (iColumn!=jColumn) {
1002            c += valueI*valueJ*elementValue;
1003          } else {
1004            c += 0.5*valueI*valueI*elementValue;
1005          }
1006        }
1007      }
1008    }
1009  }
1010  currentObj = c+linearCost;
1011  return currentObj;
1012}
1013// Scale objective
1014void 
1015ClpQuadraticObjective::reallyScale(const double * columnScale) 
1016{
1017  const int * columnQuadratic = quadraticObjective_->getIndices();
1018  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
1019  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
1020  double * quadraticElement = quadraticObjective_->getMutableElements();
1021  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
1022    double scaleI = columnScale[iColumn];
1023    objective_[iColumn] *= scaleI;
1024    CoinBigIndex j;
1025    for (j=columnQuadraticStart[iColumn];
1026         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1027      int jColumn = columnQuadratic[j];
1028      quadraticElement[j] *= scaleI*columnScale[jColumn];
1029    }
1030  }
1031}
1032/* Given a zeroed array sets nonlinear columns to 1.
1033   Returns number of nonlinear columns
1034*/
1035int 
1036ClpQuadraticObjective::markNonlinear(char * which)
1037{
1038  int iColumn;
1039  const int * columnQuadratic = quadraticObjective_->getIndices();
1040  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
1041  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
1042  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1043    CoinBigIndex j;
1044    for (j=columnQuadraticStart[iColumn];
1045         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1046      int jColumn = columnQuadratic[j];
1047      which[jColumn]=1;
1048      which[iColumn]=1;
1049    }
1050  }
1051  int numberNonLinearColumns = 0;
1052  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1053    if(which[iColumn])
1054      numberNonLinearColumns++;
1055  }
1056  return numberNonLinearColumns;
1057}
Note: See TracBrowser for help on using the repository browser.