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

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version
  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 30.3 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                  //if (fabs(valueI*valueJ*elementValue)>1.0e-12)
340                  //printf("%d %d %g %g %g -> %g\n",
341                  //       iColumn,jColumn,valueI,valueJ,elementValue,
342                  //       valueI*valueJ*elementValue);
343                  double gradientI = valueJ*elementValue;
344                  double gradientJ = valueI*elementValue;
345                  gradient_[iColumn] += gradientI;
346                  gradient_[jColumn] += gradientJ;
347                } else {
348                  offset += 0.5*valueI*valueI*elementValue;
349                  //if (fabs(valueI*valueI*elementValue)>1.0e-12)
350                  //printf("XX %d %g %g -> %g\n",
351                  //       iColumn,valueI,elementValue,
352                  //       0.5*valueI*valueI*elementValue);
353                  double gradientI = valueI*elementValue;
354                  gradient_[iColumn] += gradientI;
355                }
356              }
357            }
358          } else {
359            // full matrix
360            int iColumn;
361            offset *= 2.0;
362            for (iColumn=0;iColumn<numberColumns_;iColumn++) {
363              CoinBigIndex j;
364              double value=0.0;
365              double current = gradient_[iColumn];
366              for (j=columnQuadraticStart[iColumn];
367                   j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
368                int jColumn = columnQuadratic[j];
369                double valueJ = solution[jColumn]*quadraticElement[j];
370                value += valueJ;
371              }
372              offset += value*solution[iColumn];
373              gradient_[iColumn] = current+value;
374            }
375            offset *= 0.5;
376          }
377        }
378      }
379      if (model)
380        offset *= model->optimizationDirection()*model->objectiveScale();
381      return gradient_;
382    }
383  } else {
384    // do scaling
385    assert (solution);
386    // for now only if half
387    assert (!fullMatrix_);
388    if (refresh||!gradient_) {
389      if (!gradient_) 
390        gradient_ = new double[numberExtendedColumns_];
391      double direction = model->optimizationDirection()*model->objectiveScale();
392      // direction is actually scale out not scale in
393      //if (direction)
394      //direction = 1.0/direction;
395      const int * columnQuadratic = quadraticObjective_->getIndices();
396      const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
397      const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
398      const double * quadraticElement = quadraticObjective_->getElements();
399      int iColumn;
400      const double * columnScale = model->columnScale();
401      // use current linear cost region (already scaled)
402      if (includeLinear==1) {
403        memcpy(gradient_,model->costRegion(),numberExtendedColumns_*sizeof(double));
404      } else if (includeLinear==2) {
405        memset(gradient_+numberColumns_,0,(numberExtendedColumns_-numberColumns_)*sizeof(double));
406        if (!columnScale) {
407          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
408            gradient_[iColumn]= objective_[iColumn]*direction;
409          }
410        } else {
411          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
412            gradient_[iColumn]= objective_[iColumn]*direction*columnScale[iColumn];
413          }
414        }
415      } else {
416        memset(gradient_,0,numberExtendedColumns_*sizeof(double));
417      }
418      if (!columnScale) {
419        if (activated_) {
420          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
421            double valueI = solution[iColumn];
422            CoinBigIndex j;
423            for (j=columnQuadraticStart[iColumn];
424                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
425              int jColumn = columnQuadratic[j];
426              double valueJ = solution[jColumn];
427              double elementValue = quadraticElement[j];
428              elementValue *= direction;
429              if (iColumn!=jColumn) {
430                offset += valueI*valueJ*elementValue;
431                double gradientI = valueJ*elementValue;
432                double gradientJ = valueI*elementValue;
433                gradient_[iColumn] += gradientI;
434                gradient_[jColumn] += gradientJ;
435              } else {
436                offset += 0.5*valueI*valueI*elementValue;
437                double gradientI = valueI*elementValue;
438                gradient_[iColumn] += gradientI;
439              }
440            }
441          }
442        }
443      } else {
444        // scaling
445        if (activated_) {
446          for (iColumn=0;iColumn<numberColumns_;iColumn++) {
447            double valueI = solution[iColumn];
448            double scaleI = columnScale[iColumn]*direction;
449            CoinBigIndex j;
450            for (j=columnQuadraticStart[iColumn];
451                 j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
452              int jColumn = columnQuadratic[j];
453              double valueJ = solution[jColumn];
454              double elementValue = quadraticElement[j];
455              double scaleJ = columnScale[jColumn];
456              elementValue *= scaleI*scaleJ;
457              if (iColumn!=jColumn) {
458                offset += valueI*valueJ*elementValue;
459                double gradientI = valueJ*elementValue;
460                double gradientJ = valueI*elementValue;
461                gradient_[iColumn] += gradientI;
462                gradient_[jColumn] += gradientJ;
463              } else {
464                offset += 0.5*valueI*valueI*elementValue;
465                double gradientI = valueI*elementValue;
466                gradient_[iColumn] += gradientI;
467              }
468            }
469          }
470        }
471      }
472    }
473    if (model)
474      offset *= model->optimizationDirection();
475    return gradient_;
476  }
477}
478 
479//-------------------------------------------------------------------
480// Clone
481//-------------------------------------------------------------------
482ClpObjective * ClpQuadraticObjective::clone() const
483{
484  return new ClpQuadraticObjective(*this);
485}
486/* Subset clone.  Duplicates are allowed
487   and order is as given.
488*/
489ClpObjective * 
490ClpQuadraticObjective::subsetClone (int numberColumns, 
491                           const int * whichColumns) const
492{
493  return new ClpQuadraticObjective(*this, numberColumns, whichColumns);
494}
495// Resize objective
496void 
497ClpQuadraticObjective::resize(int newNumberColumns)
498{
499  if (numberColumns_!=newNumberColumns) {
500    int newExtended = newNumberColumns + (numberExtendedColumns_-numberColumns_);
501    int i;
502    double * newArray = new double[newExtended];
503    if (objective_)
504      memcpy(newArray,objective_,
505             CoinMin(newExtended,numberExtendedColumns_)*sizeof(double));
506    delete [] objective_;
507    objective_ = newArray;
508    for (i=numberColumns_;i<newNumberColumns;i++) 
509      objective_[i]=0.0;
510    if (gradient_) {
511      newArray = new double[newExtended];
512      if (gradient_)
513        memcpy(newArray,gradient_,
514               CoinMin(newExtended,numberExtendedColumns_)*sizeof(double));
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    memcpy(objective_+newNumberColumns,objective_+numberColumns_,
569           (numberExtendedColumns_-numberColumns_)*sizeof(double));
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    memcpy(gradient_+newNumberColumns,gradient_+numberColumns_,
596           (numberExtendedColumns_-numberColumns_)*sizeof(double));
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      memcpy(newArray,objective_,numberColumns_*sizeof(double));
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      memcpy(newArray,gradient_,numberColumns_*sizeof(double));
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  double * work = workSpace->denseVector();
695  model->factorization()->updateColumnTranspose(workSpace,&arrayVector);
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// Scale objective
896void 
897ClpQuadraticObjective::reallyScale(const double * columnScale) 
898{
899  const int * columnQuadratic = quadraticObjective_->getIndices();
900  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
901  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
902  double * quadraticElement = quadraticObjective_->getMutableElements();
903  for (int iColumn=0;iColumn<numberColumns_;iColumn++) {
904    double scaleI = columnScale[iColumn];
905    objective_[iColumn] *= scaleI;
906    CoinBigIndex j;
907    for (j=columnQuadraticStart[iColumn];
908         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
909      int jColumn = columnQuadratic[j];
910      quadraticElement[j] *= scaleI*columnScale[jColumn];
911    }
912  }
913}
914/* Given a zeroed array sets nonlinear columns to 1.
915   Returns number of nonlinear columns
916*/
917int 
918ClpQuadraticObjective::markNonlinear(char * which)
919{
920  int iColumn;
921  const int * columnQuadratic = quadraticObjective_->getIndices();
922  const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
923  const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
924  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
925    CoinBigIndex j;
926    for (j=columnQuadraticStart[iColumn];
927         j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
928      int jColumn = columnQuadratic[j];
929      which[jColumn]=1;
930      which[iColumn]=1;
931    }
932  }
933  int numberNonLinearColumns = 0;
934  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
935    if(which[iColumn])
936      numberNonLinearColumns++;
937  }
938  return numberNonLinearColumns;
939}
Note: See TracBrowser for help on using the repository browser.