source: branches/devel/Clp/src/ClpQuadraticObjective.cpp @ 920

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

minor stuff

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