source: trunk/ClpQuadraticObjective.cpp @ 468

Last change on this file since 468 was 468, checked in by forrest, 16 years ago

to make sbb faster

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 29.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//-------------------------------------------------------------------
15// Default Constructor
16//-------------------------------------------------------------------
17ClpQuadraticObjective::ClpQuadraticObjective () 
18: ClpObjective()
19{
20  type_=2;
21  objective_=NULL;
22  quadraticObjective_=NULL;
23  gradient_ = NULL;
24  numberColumns_=0;
25  numberExtendedColumns_=0;
26  activated_=0;
27  fullMatrix_=false;
28}
29
30//-------------------------------------------------------------------
31// Useful Constructor
32//-------------------------------------------------------------------
33ClpQuadraticObjective::ClpQuadraticObjective (const double * objective , 
34                                              int numberColumns,
35                                              const CoinBigIndex * start,
36                                              const int * column, const double * element,
37                                              int numberExtendedColumns)
38  : ClpObjective()
39{
40  type_=2;
41  numberColumns_ = numberColumns;
42  if (numberExtendedColumns>=0)
43    numberExtendedColumns_= max(numberColumns_,numberExtendedColumns);
44  else
45    numberExtendedColumns_= numberColumns_;
46  if (objective) {
47    objective_ = new double [numberExtendedColumns_];
48    memcpy(objective_,objective,numberColumns_*sizeof(double));
49    memset(objective_+numberColumns_,0,(numberExtendedColumns_-numberColumns_)*sizeof(double));
50  } else {
51    objective_ = new double [numberExtendedColumns_];
52    memset(objective_,0,numberExtendedColumns_*sizeof(double));
53  }
54  if (start) 
55    quadraticObjective_ = new CoinPackedMatrix(true,numberColumns,numberColumns,
56                                             start[numberColumns],element,column,start,NULL);
57  else
58  quadraticObjective_=NULL;
59  gradient_ = NULL;
60  activated_=1;
61  fullMatrix_=false;
62}
63
64//-------------------------------------------------------------------
65// Copy constructor
66//-------------------------------------------------------------------
67ClpQuadraticObjective::ClpQuadraticObjective (const ClpQuadraticObjective & rhs,
68                                              int type) 
69: ClpObjective(rhs)
70{ 
71  numberColumns_=rhs.numberColumns_;
72  numberExtendedColumns_=rhs.numberExtendedColumns_;
73  fullMatrix_=rhs.fullMatrix_;
74  if (rhs.objective_) {
75    objective_ = new double [numberExtendedColumns_];
76    memcpy(objective_,rhs.objective_,numberExtendedColumns_*sizeof(double));
77  } else {
78    objective_=NULL;
79  }
80  if (rhs.gradient_) {
81    gradient_ = new double [numberExtendedColumns_];
82    memcpy(gradient_,rhs.gradient_,numberExtendedColumns_*sizeof(double));
83  } else {
84    gradient_=NULL;
85  }
86  if (rhs.quadraticObjective_) {
87    // see what type of matrix wanted
88    if (type==0) {
89      // just copy
90      quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
91    } else if (type==1) {
92      // expand to full symmetric
93      fullMatrix_=true;
94      const int * columnQuadratic1 = rhs.quadraticObjective_->getIndices();
95      const CoinBigIndex * columnQuadraticStart1 = rhs.quadraticObjective_->getVectorStarts();
96      const int * columnQuadraticLength1 = rhs.quadraticObjective_->getVectorLengths();
97      const double * quadraticElement1 = rhs.quadraticObjective_->getElements();
98      CoinBigIndex * columnQuadraticStart2 = new CoinBigIndex [numberExtendedColumns_+1];
99      int * columnQuadraticLength2 = new int [numberExtendedColumns_];
100      int iColumn;
101      int numberColumns = rhs.quadraticObjective_->getNumCols();
102      int numberBelow=0;
103      int numberAbove=0;
104      int numberDiagonal=0;
105      CoinZeroN(columnQuadraticLength2,numberExtendedColumns_);
106      for (iColumn=0;iColumn<numberColumns;iColumn++) {
107        for (CoinBigIndex j=columnQuadraticStart1[iColumn];
108             j<columnQuadraticStart1[iColumn]+columnQuadraticLength1[iColumn];j++) {
109          int jColumn = columnQuadratic1[j];
110          if (jColumn>iColumn) {
111            numberBelow++;
112            columnQuadraticLength2[jColumn]++;
113            columnQuadraticLength2[iColumn]++;
114          } else if (jColumn==iColumn) {
115            numberDiagonal++;
116            columnQuadraticLength2[iColumn]++;
117          } else {
118            numberAbove++;
119          }
120        }
121      }
122      if (numberAbove>0) {
123        if (numberAbove==numberBelow) {
124          // already done
125          quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
126          delete [] columnQuadraticStart2;
127          delete [] columnQuadraticLength2;
128        } else {
129          printf("number above = %d, number below = %d, error\n",
130                 numberAbove,numberBelow);
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    ClpObjective::operator=(rhs);
265    numberColumns_=rhs.numberColumns_;
266    numberExtendedColumns_=rhs.numberExtendedColumns_;
267    if (rhs.objective_) {
268      objective_ = new double [numberExtendedColumns_];
269      memcpy(objective_,rhs.objective_,numberExtendedColumns_*sizeof(double));
270    } else {
271      objective_=NULL;
272    }
273    if (rhs.gradient_) {
274      gradient_ = new double [numberExtendedColumns_];
275      memcpy(gradient_,rhs.gradient_,numberExtendedColumns_*sizeof(double));
276    } else {
277      gradient_=NULL;
278    }
279    if (rhs.quadraticObjective_) {
280      quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
281    } else {
282      quadraticObjective_=NULL;
283    }
284  }
285  return *this;
286}
287
288// Returns gradient
289double * 
290ClpQuadraticObjective::gradient(const ClpSimplex * model,
291                                const double * solution, double & offset,bool refresh,
292                                int includeLinear)
293{
294  offset=0.0;
295  bool scaling=false;
296  if (model&&(model->rowScale()||
297              model->objectiveScale()!=1.0||model->optimizationDirection()!=1.0))
298    scaling=true;
299  const double * cost = NULL;
300  if (model)
301    cost = model->costRegion();
302  if (!cost) {
303    // not in solve
304    cost = objective_;
305    scaling=false;
306  }
307  if (!scaling) {
308    if (!quadraticObjective_||!solution||!activated_) {
309      return objective_;
310    } else {
311      if (refresh||!gradient_) {
312        if (!gradient_) 
313          gradient_ = new double[numberExtendedColumns_];
314        const int * columnQuadratic = quadraticObjective_->getIndices();
315        const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
316        const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
317        const double * quadraticElement = quadraticObjective_->getElements();
318        offset=0.0;
319        // use current linear cost region
320        if (includeLinear==1)
321          memcpy(gradient_,cost,numberExtendedColumns_*sizeof(double));
322        else if (includeLinear==2)
323          memcpy(gradient_,objective_,numberExtendedColumns_*sizeof(double));
324        else
325          memset(gradient_,0,numberExtendedColumns_*sizeof(double));
326        if (activated_) {
327          if (!fullMatrix_) {
328            int iColumn;
329            for (iColumn=0;iColumn<numberColumns_;iColumn++) {
330              double valueI = solution[iColumn];
331              CoinBigIndex j;
332              for (j=columnQuadraticStart[iColumn];
333                   j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
334                int jColumn = columnQuadratic[j];
335                double valueJ = solution[jColumn];
336                double elementValue = quadraticElement[j];
337                if (iColumn!=jColumn) {
338                  offset += valueI*valueJ*elementValue;
339                  double gradientI = valueJ*elementValue;
340                  double gradientJ = valueI*elementValue;
341                  gradient_[iColumn] += gradientI;
342                  gradient_[jColumn] += gradientJ;
343                } else {
344                  offset += 0.5*valueI*valueI*elementValue;
345                  double gradientI = valueI*elementValue;
346                  gradient_[iColumn] += gradientI;
347                }
348              }
349            }
350          } else {
351            // full matrix
352            int iColumn;
353            offset *= 2.0;
354            for (iColumn=0;iColumn<numberColumns_;iColumn++) {
355              CoinBigIndex j;
356              double value=0.0;
357              double current = gradient_[iColumn];
358              for (j=columnQuadraticStart[iColumn];
359                   j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
360                int jColumn = columnQuadratic[j];
361                double valueJ = solution[jColumn]*quadraticElement[j];
362                value += valueJ;
363              }
364              offset += value*solution[iColumn];
365              gradient_[iColumn] = current+value;
366            }
367            offset *= 0.5;
368          }
369        }
370      }
371      offset *= model->optimizationDirection()*model->objectiveScale();
372      return gradient_;
373    }
374  } else {
375    // do scaling
376    assert (solution);
377    // for now only if half
378    assert (!fullMatrix_);
379    if (refresh||!gradient_) {
380      if (!gradient_) 
381        gradient_ = new double[numberExtendedColumns_];
382      double direction = model->optimizationDirection()*model->objectiveScale();
383      // direction is actually scale out not scale in
384      //if (direction)
385      //direction = 1.0/direction;
386      const int * columnQuadratic = quadraticObjective_->getIndices();
387      const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
388      const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
389      const double * quadraticElement = quadraticObjective_->getElements();
390      int iColumn;
391      const double * columnScale = model->columnScale();
392      // use current linear cost region (already scaled)
393      if (includeLinear==1) {
394        memcpy(gradient_,model->costRegion(),numberExtendedColumns_*sizeof(double));
395      } else if (includeLinear==2) {
396        memset(gradient_+numberColumns_,0,(numberExtendedColumns_-numberColumns_)*sizeof(double));
397        if (!columnScale) {
398          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
399            gradient_[iColumn]= objective_[iColumn]*direction;
400          }
401        } else {
402          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
403            gradient_[iColumn]= objective_[iColumn]*direction*columnScale[iColumn];
404          }
405        }
406      } else {
407        memset(gradient_,0,numberExtendedColumns_*sizeof(double));
408      }
409      if (!columnScale) {
410        if (activated_) {
411          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
412            double valueI = solution[iColumn];
413            CoinBigIndex j;
414            for (j=columnQuadraticStart[iColumn];
415                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
416              int jColumn = columnQuadratic[j];
417              double valueJ = solution[jColumn];
418              double elementValue = quadraticElement[j];
419              elementValue *= direction;
420              if (iColumn!=jColumn) {
421                offset += valueI*valueJ*elementValue;
422                double gradientI = valueJ*elementValue;
423                double gradientJ = valueI*elementValue;
424                gradient_[iColumn] += gradientI;
425                gradient_[jColumn] += gradientJ;
426              } else {
427                offset += 0.5*valueI*valueI*elementValue;
428                double gradientI = valueI*elementValue;
429                gradient_[iColumn] += gradientI;
430              }
431            }
432          }
433        }
434      } else {
435        // scaling
436        if (activated_) {
437          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
438            double valueI = solution[iColumn];
439            double scaleI = columnScale[iColumn]*direction;
440            CoinBigIndex j;
441            for (j=columnQuadraticStart[iColumn];
442                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
443              int jColumn = columnQuadratic[j];
444              double valueJ = solution[jColumn];
445              double elementValue = quadraticElement[j];
446              double scaleJ = columnScale[jColumn];
447              elementValue *= scaleI*scaleJ;
448              if (iColumn!=jColumn) {
449                offset += valueI*valueJ*elementValue;
450                double gradientI = valueJ*elementValue;
451                double gradientJ = valueI*elementValue;
452                gradient_[iColumn] += gradientI;
453                gradient_[jColumn] += gradientJ;
454              } else {
455                offset += 0.5*valueI*valueI*elementValue;
456                double gradientI = valueI*elementValue;
457                gradient_[iColumn] += gradientI;
458              }
459            }
460          }
461        }
462      }
463    }
464    offset *= model->optimizationDirection();
465    return gradient_;
466  }
467}
468 
469//-------------------------------------------------------------------
470// Clone
471//-------------------------------------------------------------------
472ClpObjective * ClpQuadraticObjective::clone() const
473{
474  return new ClpQuadraticObjective(*this);
475}
476/* Subset clone.  Duplicates are allowed
477   and order is as given.
478*/
479ClpObjective * 
480ClpQuadraticObjective::subsetClone (int numberColumns, 
481                           const int * whichColumns) const
482{
483  return new ClpQuadraticObjective(*this, numberColumns, whichColumns);
484}
485// Resize objective
486void 
487ClpQuadraticObjective::resize(int newNumberColumns)
488{
489  if (numberColumns_!=newNumberColumns) {
490    int newExtended = newNumberColumns + (numberExtendedColumns_-numberColumns_);
491    int i;
492    double * newArray = new double[newExtended];
493    if (objective_)
494      memcpy(newArray,objective_,
495             CoinMin(newExtended,numberExtendedColumns_)*sizeof(double));
496    delete [] objective_;
497    objective_ = newArray;
498    for (i=numberColumns_;i<newNumberColumns;i++) 
499      objective_[i]=0.0;
500    if (gradient_) {
501      newArray = new double[newExtended];
502      if (gradient_)
503        memcpy(newArray,gradient_,
504               CoinMin(newExtended,numberExtendedColumns_)*sizeof(double));
505      delete [] gradient_;
506      gradient_ = newArray;
507      for (i=numberColumns_;i<newNumberColumns;i++) 
508        gradient_[i]=0.0;
509    }
510    if (quadraticObjective_) {
511      if (newNumberColumns<numberColumns_) {
512        int * which = new int[numberColumns_-newNumberColumns];
513        int i;
514        for (i=newNumberColumns;i<numberColumns_;i++) 
515          which[i-newNumberColumns]=i;
516        quadraticObjective_->deleteRows(numberColumns_-newNumberColumns,which);
517        quadraticObjective_->deleteCols(numberColumns_-newNumberColumns,which);
518        delete [] which;
519      } else {
520        quadraticObjective_->setDimensions(newNumberColumns,newNumberColumns);
521      }
522    }
523    numberColumns_ = newNumberColumns;
524    numberExtendedColumns_ = newExtended;
525  } 
526 
527}
528// Delete columns in  objective
529void 
530ClpQuadraticObjective::deleteSome(int numberToDelete, const int * which) 
531{
532  int newNumberColumns = numberColumns_-numberToDelete;
533  int newExtended = numberExtendedColumns_ - numberToDelete;
534  if (objective_) {
535    int i ;
536    char * deleted = new char[numberColumns_];
537    int numberDeleted=0;
538    memset(deleted,0,numberColumns_*sizeof(char));
539    for (i=0;i<numberToDelete;i++) {
540      int j = which[i];
541      if (j>=0&&j<numberColumns_&&!deleted[j]) {
542        numberDeleted++;
543        deleted[j]=1;
544      }
545    }
546    newNumberColumns = numberColumns_-numberDeleted;
547    newExtended = numberExtendedColumns_ - numberDeleted;
548    double * newArray = new double[newExtended];
549    int put=0;
550    for (i=0;i<numberColumns_;i++) {
551      if (!deleted[i]) {
552        newArray[put++]=objective_[i];
553      }
554    }
555    delete [] objective_;
556    objective_ = newArray;
557    delete [] deleted;
558    memcpy(objective_+newNumberColumns,objective_+numberColumns_,
559           (numberExtendedColumns_-numberColumns_)*sizeof(double));
560  }
561  if (gradient_) {
562    int i ;
563    char * deleted = new char[numberColumns_];
564    int numberDeleted=0;
565    memset(deleted,0,numberColumns_*sizeof(char));
566    for (i=0;i<numberToDelete;i++) {
567      int j = which[i];
568      if (j>=0&&j<numberColumns_&&!deleted[j]) {
569        numberDeleted++;
570        deleted[j]=1;
571      }
572    }
573    newNumberColumns = numberColumns_-numberDeleted;
574    newExtended = numberExtendedColumns_ - numberDeleted;
575    double * newArray = new double[newExtended];
576    int put=0;
577    for (i=0;i<numberColumns_;i++) {
578      if (!deleted[i]) {
579        newArray[put++]=gradient_[i];
580      }
581    }
582    delete [] gradient_;
583    gradient_ = newArray;
584    delete [] deleted;
585    memcpy(gradient_+newNumberColumns,gradient_+numberColumns_,
586           (numberExtendedColumns_-numberColumns_)*sizeof(double));
587  }
588  numberColumns_ = newNumberColumns;
589  numberExtendedColumns_ = newExtended;
590  if (quadraticObjective_) {
591    quadraticObjective_->deleteCols(numberToDelete,which);
592    quadraticObjective_->deleteRows(numberToDelete,which);
593  }
594}
595
596// Load up quadratic objective
597void 
598ClpQuadraticObjective::loadQuadraticObjective(const int numberColumns, const CoinBigIndex * start,
599                              const int * column, const double * element,int numberExtended)
600{
601  fullMatrix_=false;
602  delete quadraticObjective_;
603  quadraticObjective_ = new CoinPackedMatrix(true,numberColumns,numberColumns,
604                                             start[numberColumns],element,column,start,NULL);
605  numberColumns_=numberColumns;
606  if (numberExtended>numberExtendedColumns_) {
607    if (objective_) {
608      // make correct size
609      double * newArray = new double[numberExtended];
610      memcpy(newArray,objective_,numberColumns_*sizeof(double));
611      delete [] objective_;
612      objective_ = newArray;
613      memset(objective_+numberColumns_,0,(numberExtended-numberColumns_)*sizeof(double));
614    }
615    if (gradient_) {
616      // make correct size
617      double * newArray = new double[numberExtended];
618      memcpy(newArray,gradient_,numberColumns_*sizeof(double));
619      delete [] gradient_;
620      gradient_ = newArray;
621      memset(gradient_+numberColumns_,0,(numberExtended-numberColumns_)*sizeof(double));
622    }
623    numberExtendedColumns_ = numberExtended;
624  } else {
625    numberExtendedColumns_ = numberColumns_;
626  }
627}
628void 
629ClpQuadraticObjective::loadQuadraticObjective (  const CoinPackedMatrix& matrix)
630{
631  delete quadraticObjective_;
632  quadraticObjective_ = new CoinPackedMatrix(matrix);
633}
634// Get rid of quadratic objective
635void 
636ClpQuadraticObjective::deleteQuadraticObjective()
637{
638  delete quadraticObjective_;
639  quadraticObjective_ = NULL;
640}
641/* Returns reduced gradient.Returns an offset (to be added to current one).
642 */
643double 
644ClpQuadraticObjective::reducedGradient(ClpSimplex * model, double * region,
645                                       bool useFeasibleCosts)
646{
647  int numberRows = model->numberRows();
648  int numberColumns=model->numberColumns();
649 
650  //work space
651  CoinIndexedVector  * workSpace = model->rowArray(0);
652 
653  CoinIndexedVector arrayVector;
654  arrayVector.reserve(numberRows+1);
655 
656  int iRow;
657#ifdef CLP_DEBUG
658  workSpace->checkClear();
659#endif
660  double * array = arrayVector.denseVector();
661  int * index = arrayVector.getIndices();
662  int number=0;
663  const double * costNow = gradient(model,model->solutionRegion(),offset_,
664                                    true,useFeasibleCosts ? 2 : 1);
665  double * cost = model->costRegion();
666  const int * pivotVariable = model->pivotVariable();
667  for (iRow=0;iRow<numberRows;iRow++) {
668    int iPivot=pivotVariable[iRow];
669    double value;
670    if (iPivot<numberColumns)
671      value = costNow[iPivot];
672    else if (!useFeasibleCosts) 
673      value = cost[iPivot];
674    else 
675      value=0.0;
676    if (value) {
677      array[iRow]=value;
678      index[number++]=iRow;
679    }
680  }
681  arrayVector.setNumElements(number);
682
683  // Btran basic costs
684  double * work = workSpace->denseVector();
685  model->factorization()->updateColumnTranspose(workSpace,&arrayVector);
686  ClpFillN(work,numberRows,0.0);
687  // now look at dual solution
688  double * rowReducedCost = region+numberColumns;
689  double * dual = rowReducedCost;
690  const double * rowCost = cost+numberColumns;
691  for (iRow=0;iRow<numberRows;iRow++) {
692    dual[iRow]=array[iRow];
693  }
694  double * dj = region;
695  ClpDisjointCopyN(costNow,numberColumns,dj);
696 
697  model->transposeTimes(-1.0,dual,dj);
698  for (iRow=0;iRow<numberRows;iRow++) {
699    // slack
700    double value = dual[iRow];
701    value += rowCost[iRow];
702    rowReducedCost[iRow]=value;
703  }
704  return offset_;
705}
706/* Returns step length which gives minimum of objective for
707   solution + theta * change vector up to maximum theta.
708   
709   arrays are numberColumns+numberRows
710*/
711double 
712ClpQuadraticObjective::stepLength(ClpSimplex * model,
713                                  const double * solution,
714                                  const double * change,
715                                  double maximumTheta,
716                                  double & currentObj,
717                                  double & predictedObj,
718                                  double & thetaObj)
719{
720  const double * cost = model->costRegion();
721  bool inSolve=true;
722  if (!cost) {
723    // not in solve
724    cost = objective_;
725    inSolve=false;
726  }
727  double delta=0.0;
728  double linearCost =0.0;
729  int numberRows = model->numberRows();
730  int numberColumns = model->numberColumns();
731  int numberTotal = numberColumns;
732  if (inSolve)
733    numberTotal += numberRows;
734  currentObj=0.0;
735  thetaObj=0.0;
736  for (int iColumn=0;iColumn<numberTotal;iColumn++) {
737    delta += cost[iColumn]*change[iColumn];
738    linearCost += cost[iColumn]*solution[iColumn];
739  }
740  if (!activated_||!quadraticObjective_) {
741    currentObj=linearCost;
742    thetaObj =currentObj + delta*maximumTheta;
743    if (delta<0.0) {
744      return maximumTheta;
745    } else {
746      printf("odd linear direction %g\n",delta);
747      return 0.0;
748    }
749  }
750  assert (model);
751  bool scaling=false;
752  if ((model->rowScale()||
753       model->objectiveScale()!=1.0||model->optimizationDirection()!=1.0)&&inSolve)
754    scaling=true;
755  const int * columnQuadratic = quadraticObjective_->getIndices();
756  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
757  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
758  const double * quadraticElement = quadraticObjective_->getElements();
759  double a=0.0;
760  double b=delta;
761  double c=0.0;
762  if (!scaling) {
763    if (!fullMatrix_) {
764      int iColumn;
765      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
766        double valueI = solution[iColumn];
767        double changeI = change[iColumn];
768        CoinBigIndex j;
769        for (j=columnQuadraticStart[iColumn];
770             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
771          int jColumn = columnQuadratic[j];
772          double valueJ = solution[jColumn];
773          double changeJ = change[jColumn];
774          double elementValue = quadraticElement[j];
775          if (iColumn!=jColumn) {
776            a += changeI*changeJ*elementValue;
777            b += (changeI*valueJ+changeJ*valueI)*elementValue;
778            c += valueI*valueJ*elementValue;
779          } else {
780            a += 0.5*changeI*changeI*elementValue;
781            b += changeI*valueI*elementValue;
782            c += 0.5*valueI*valueI*elementValue;
783          }
784        }
785      }
786    } else {
787      // full matrix stored
788      int iColumn;
789      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
790        double valueI = solution[iColumn];
791        double changeI = change[iColumn];
792        CoinBigIndex j;
793        for (j=columnQuadraticStart[iColumn];
794             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
795          int jColumn = columnQuadratic[j];
796          double valueJ = solution[jColumn];
797          double changeJ = change[jColumn];
798          double elementValue = quadraticElement[j];
799          valueJ *= elementValue;
800          a += changeI*changeJ*elementValue;
801          b += changeI*valueJ;
802          c += valueI*valueJ;
803        }
804      }
805      a *= 0.5;
806      c *= 0.5;
807    }
808  } else {
809    // scaling
810    // for now only if half
811    assert (!fullMatrix_);
812    const double * columnScale = model->columnScale();
813    double direction = model->optimizationDirection()*model->objectiveScale();
814    // direction is actually scale out not scale in
815    if (direction)
816      direction = 1.0/direction;
817    if (!columnScale) {
818      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
819        double valueI = solution[iColumn];
820        double changeI = change[iColumn];
821        CoinBigIndex j;
822        for (j=columnQuadraticStart[iColumn];
823             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
824          int jColumn = columnQuadratic[j];
825          double valueJ = solution[jColumn];
826          double changeJ = change[jColumn];
827          double elementValue = quadraticElement[j];
828          elementValue *= direction;
829          if (iColumn!=jColumn) {
830            a += changeI*changeJ*elementValue;
831            b += (changeI*valueJ+changeJ*valueI)*elementValue;
832            c += valueI*valueJ*elementValue;
833          } else {
834            a += 0.5*changeI*changeI*elementValue;
835            b += changeI*valueI*elementValue;
836            c += 0.5*valueI*valueI*elementValue;
837          }
838        }
839      }
840    } else {
841      // scaling
842      for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
843        double valueI = solution[iColumn];
844        double changeI = change[iColumn];
845        double scaleI = columnScale[iColumn]*direction;
846        CoinBigIndex j;
847        for (j=columnQuadraticStart[iColumn];
848             j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
849          int jColumn = columnQuadratic[j];
850          double valueJ = solution[jColumn];
851          double changeJ = change[jColumn];
852          double elementValue = quadraticElement[j];
853          elementValue *= scaleI*columnScale[jColumn];
854          if (iColumn!=jColumn) {
855            a += changeI*changeJ*elementValue;
856            b += (changeI*valueJ+changeJ*valueI)*elementValue;
857            c += valueI*valueJ*elementValue;
858          } else {
859            a += 0.5*changeI*changeI*elementValue;
860            b += changeI*valueI*elementValue;
861            c += 0.5*valueI*valueI*elementValue;
862          }
863        }
864      }
865    }
866  }
867  double theta;
868  //printf("Current cost %g\n",c+linearCost);
869  currentObj = c+linearCost;
870  thetaObj = currentObj + a*maximumTheta*maximumTheta+b*maximumTheta;
871  // minimize a*x*x + b*x + c
872  if (a<=0.0) {
873    theta = maximumTheta;
874  } else {
875    theta = -0.5*b/a;
876  }
877  predictedObj = currentObj + a*theta*theta+b*theta;
878  if (b>0.0) {
879    if (model->messageHandler()->logLevel()&32)
880      printf("a %g b %g c %g => %g\n",a,b,c,theta); 
881    b=0.0;
882  }
883  return CoinMin(theta,maximumTheta);
884}
885// Scale objective
886void 
887ClpQuadraticObjective::reallyScale(const double * columnScale) 
888{
889  const int * columnQuadratic = quadraticObjective_->getIndices();
890  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
891  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
892  double * quadraticElement = quadraticObjective_->getMutableElements();
893  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
894    double scaleI = columnScale[iColumn];
895    objective_[iColumn] *= scaleI;
896    CoinBigIndex j;
897    for (j=columnQuadraticStart[iColumn];
898         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
899      int jColumn = columnQuadratic[j];
900      quadraticElement[j] *= scaleI*columnScale[jColumn];
901    }
902  }
903}
904/* Given a zeroed array sets nonlinear columns to 1.
905   Returns number of nonlinear columns
906*/
907int 
908ClpQuadraticObjective::markNonlinear(char * which)
909{
910  int iColumn;
911  const int * columnQuadratic = quadraticObjective_->getIndices();
912  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
913  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
914  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
915    CoinBigIndex j;
916    for (j=columnQuadraticStart[iColumn];
917         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
918      int jColumn = columnQuadratic[j];
919      which[jColumn]=1;
920      which[iColumn]=1;
921    }
922  }
923  int numberNonLinearColumns = 0;
924  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
925    if(which[iColumn])
926      numberNonLinearColumns++;
927  }
928  return numberNonLinearColumns;
929}
Note: See TracBrowser for help on using the repository browser.