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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 46.6 KB
Line 
1/* $Id: ClpQuadraticObjective.cpp 1525 2010-02-26 17:27:59Z mjs $ */
2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5#include "CoinPragma.hpp"
6#include "CoinHelperFunctions.hpp"
7#include "CoinIndexedVector.hpp"
8#include "ClpFactorization.hpp"
9#include "ClpSimplex.hpp"
10#include "ClpQuadraticObjective.hpp"
11//#############################################################################
12// Constructors / Destructor / Assignment
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_ = CoinMax(numberColumns_, numberExtendedColumns);
44     else
45          numberExtendedColumns_ = numberColumns_;
46     if (objective) {
47          objective_ = new double [numberExtendedColumns_];
48          CoinMemcpyN(objective, numberColumns_, objective_);
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          CoinMemcpyN(rhs.objective_, numberExtendedColumns_, objective_);
77     } else {
78          objective_ = NULL;
79     }
80     if (rhs.gradient_) {
81          gradient_ = new double [numberExtendedColumns_];
82          CoinMemcpyN(rhs.gradient_, numberExtendedColumns_, gradient_);
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          CoinMemcpyN(rhs.objective_ + rhs.numberColumns_,      (numberExtendedColumns_ - numberColumns_),
221                      objective_ + numberColumns_);
222          if (rhs.gradient_) {
223               gradient_ = new double[numberExtendedColumns_];
224               for (i = 0; i < numberColumns_; i++)
225                    gradient_[i] = rhs.gradient_[whichColumn[i]];
226               CoinMemcpyN(rhs.gradient_ + rhs.numberColumns_,  (numberExtendedColumns_ - numberColumns_),
227                           gradient_ + numberColumns_);
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          delete [] objective_;
266          delete [] gradient_;
267          ClpObjective::operator=(rhs);
268          numberColumns_ = rhs.numberColumns_;
269          numberExtendedColumns_ = rhs.numberExtendedColumns_;
270          if (rhs.objective_) {
271               objective_ = new double [numberExtendedColumns_];
272               CoinMemcpyN(rhs.objective_, numberExtendedColumns_, objective_);
273          } else {
274               objective_ = NULL;
275          }
276          if (rhs.gradient_) {
277               gradient_ = new double [numberExtendedColumns_];
278               CoinMemcpyN(rhs.gradient_, numberExtendedColumns_, gradient_);
279          } else {
280               gradient_ = NULL;
281          }
282          if (rhs.quadraticObjective_) {
283               quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
284          } else {
285               quadraticObjective_ = NULL;
286          }
287     }
288     return *this;
289}
290
291// Returns gradient
292double *
293ClpQuadraticObjective::gradient(const ClpSimplex * model,
294                                const double * solution, double & offset, bool refresh,
295                                int includeLinear)
296{
297     offset = 0.0;
298     bool scaling = false;
299     if (model && (model->rowScale() ||
300                   model->objectiveScale() != 1.0 || model->optimizationDirection() != 1.0))
301          scaling = true;
302     const double * cost = NULL;
303     if (model)
304          cost = model->costRegion();
305     if (!cost) {
306          // not in solve
307          cost = objective_;
308          scaling = false;
309     }
310     if (!scaling) {
311          if (!quadraticObjective_ || !solution || !activated_) {
312               return objective_;
313          } else {
314               if (refresh || !gradient_) {
315                    if (!gradient_)
316                         gradient_ = new double[numberExtendedColumns_];
317                    const int * columnQuadratic = quadraticObjective_->getIndices();
318                    const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
319                    const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
320                    const double * quadraticElement = quadraticObjective_->getElements();
321                    offset = 0.0;
322                    // use current linear cost region
323                    if (includeLinear == 1)
324                         CoinMemcpyN(cost, numberExtendedColumns_, gradient_);
325                    else if (includeLinear == 2)
326                         CoinMemcpyN(objective_, numberExtendedColumns_, gradient_);
327                    else
328                         memset(gradient_, 0, numberExtendedColumns_ * sizeof(double));
329                    if (activated_) {
330                         if (!fullMatrix_) {
331                              int iColumn;
332                              for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
333                                   double valueI = solution[iColumn];
334                                   CoinBigIndex j;
335                                   for (j = columnQuadraticStart[iColumn];
336                                             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
337                                        int jColumn = columnQuadratic[j];
338                                        double valueJ = solution[jColumn];
339                                        double elementValue = quadraticElement[j];
340                                        if (iColumn != jColumn) {
341                                             offset += valueI * valueJ * elementValue;
342                                             //if (fabs(valueI*valueJ*elementValue)>1.0e-12)
343                                             //printf("%d %d %g %g %g -> %g\n",
344                                             //       iColumn,jColumn,valueI,valueJ,elementValue,
345                                             //       valueI*valueJ*elementValue);
346                                             double gradientI = valueJ * elementValue;
347                                             double gradientJ = valueI * elementValue;
348                                             gradient_[iColumn] += gradientI;
349                                             gradient_[jColumn] += gradientJ;
350                                        } else {
351                                             offset += 0.5 * valueI * valueI * elementValue;
352                                             //if (fabs(valueI*valueI*elementValue)>1.0e-12)
353                                             //printf("XX %d %g %g -> %g\n",
354                                             //       iColumn,valueI,elementValue,
355                                             //       0.5*valueI*valueI*elementValue);
356                                             double gradientI = valueI * elementValue;
357                                             gradient_[iColumn] += gradientI;
358                                        }
359                                   }
360                              }
361                         } else {
362                              // full matrix
363                              int iColumn;
364                              offset *= 2.0;
365                              for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
366                                   CoinBigIndex j;
367                                   double value = 0.0;
368                                   double current = gradient_[iColumn];
369                                   for (j = columnQuadraticStart[iColumn];
370                                             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
371                                        int jColumn = columnQuadratic[j];
372                                        double valueJ = solution[jColumn] * quadraticElement[j];
373                                        value += valueJ;
374                                   }
375                                   offset += value * solution[iColumn];
376                                   gradient_[iColumn] = current + value;
377                              }
378                              offset *= 0.5;
379                         }
380                    }
381               }
382               if (model)
383                    offset *= model->optimizationDirection() * model->objectiveScale();
384               return gradient_;
385          }
386     } else {
387          // do scaling
388          assert (solution);
389          // for now only if half
390          assert (!fullMatrix_);
391          if (refresh || !gradient_) {
392               if (!gradient_)
393                    gradient_ = new double[numberExtendedColumns_];
394               double direction = model->optimizationDirection() * model->objectiveScale();
395               // direction is actually scale out not scale in
396               //if (direction)
397               //direction = 1.0/direction;
398               const int * columnQuadratic = quadraticObjective_->getIndices();
399               const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
400               const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
401               const double * quadraticElement = quadraticObjective_->getElements();
402               int iColumn;
403               const double * columnScale = model->columnScale();
404               // use current linear cost region (already scaled)
405               if (includeLinear == 1) {
406                    CoinMemcpyN(model->costRegion(), numberExtendedColumns_, gradient_);
407               }        else if (includeLinear == 2) {
408                    memset(gradient_ + numberColumns_, 0, (numberExtendedColumns_ - numberColumns_)*sizeof(double));
409                    if (!columnScale) {
410                         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
411                              gradient_[iColumn] = objective_[iColumn] * direction;
412                         }
413                    } else {
414                         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
415                              gradient_[iColumn] = objective_[iColumn] * direction * columnScale[iColumn];
416                         }
417                    }
418               } else {
419                    memset(gradient_, 0, numberExtendedColumns_ * sizeof(double));
420               }
421               if (!columnScale) {
422                    if (activated_) {
423                         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
424                              double valueI = solution[iColumn];
425                              CoinBigIndex j;
426                              for (j = columnQuadraticStart[iColumn];
427                                        j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
428                                   int jColumn = columnQuadratic[j];
429                                   double valueJ = solution[jColumn];
430                                   double elementValue = quadraticElement[j];
431                                   elementValue *= direction;
432                                   if (iColumn != jColumn) {
433                                        offset += valueI * valueJ * elementValue;
434                                        double gradientI = valueJ * elementValue;
435                                        double gradientJ = valueI * elementValue;
436                                        gradient_[iColumn] += gradientI;
437                                        gradient_[jColumn] += gradientJ;
438                                   } else {
439                                        offset += 0.5 * valueI * valueI * elementValue;
440                                        double gradientI = valueI * elementValue;
441                                        gradient_[iColumn] += gradientI;
442                                   }
443                              }
444                         }
445                    }
446               } else {
447                    // scaling
448                    if (activated_) {
449                         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
450                              double valueI = solution[iColumn];
451                              double scaleI = columnScale[iColumn] * direction;
452                              CoinBigIndex j;
453                              for (j = columnQuadraticStart[iColumn];
454                                        j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
455                                   int jColumn = columnQuadratic[j];
456                                   double valueJ = solution[jColumn];
457                                   double elementValue = quadraticElement[j];
458                                   double scaleJ = columnScale[jColumn];
459                                   elementValue *= scaleI * scaleJ;
460                                   if (iColumn != jColumn) {
461                                        offset += valueI * valueJ * elementValue;
462                                        double gradientI = valueJ * elementValue;
463                                        double gradientJ = valueI * elementValue;
464                                        gradient_[iColumn] += gradientI;
465                                        gradient_[jColumn] += gradientJ;
466                                   } else {
467                                        offset += 0.5 * valueI * valueI * elementValue;
468                                        double gradientI = valueI * elementValue;
469                                        gradient_[iColumn] += gradientI;
470                                   }
471                              }
472                         }
473                    }
474               }
475          }
476          if (model)
477               offset *= model->optimizationDirection();
478          return gradient_;
479     }
480}
481
482//-------------------------------------------------------------------
483// Clone
484//-------------------------------------------------------------------
485ClpObjective * ClpQuadraticObjective::clone() const
486{
487     return new ClpQuadraticObjective(*this);
488}
489/* Subset clone.  Duplicates are allowed
490   and order is as given.
491*/
492ClpObjective *
493ClpQuadraticObjective::subsetClone (int numberColumns,
494                                    const int * whichColumns) const
495{
496     return new ClpQuadraticObjective(*this, numberColumns, whichColumns);
497}
498// Resize objective
499void
500ClpQuadraticObjective::resize(int newNumberColumns)
501{
502     if (numberColumns_ != newNumberColumns) {
503          int newExtended = newNumberColumns + (numberExtendedColumns_ - numberColumns_);
504          int i;
505          double * newArray = new double[newExtended];
506          if (objective_)
507               CoinMemcpyN(objective_,  CoinMin(newExtended, numberExtendedColumns_), newArray);
508          delete [] objective_;
509          objective_ = newArray;
510          for (i = numberColumns_; i < newNumberColumns; i++)
511               objective_[i] = 0.0;
512          if (gradient_) {
513               newArray = new double[newExtended];
514               if (gradient_)
515                    CoinMemcpyN(gradient_,      CoinMin(newExtended, numberExtendedColumns_), newArray);
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          CoinMemcpyN(objective_ + numberColumns_,      (numberExtendedColumns_ - numberColumns_),
570                      objective_ + newNumberColumns);
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          CoinMemcpyN(gradient_ + numberColumns_,       (numberExtendedColumns_ - numberColumns_),
597                      gradient_ + newNumberColumns);
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               CoinMemcpyN(objective_, numberColumns_, newArray);
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               CoinMemcpyN(gradient_, numberColumns_, newArray);
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// Return objective value (without any ClpModel offset) (model may be NULL)
897double
898ClpQuadraticObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
899{
900     bool scaling = false;
901     if (model && (model->rowScale() ||
902                   model->objectiveScale() != 1.0))
903          scaling = true;
904     const double * cost = NULL;
905     if (model)
906          cost = model->costRegion();
907     if (!cost) {
908          // not in solve
909          cost = objective_;
910          scaling = false;
911     }
912     double linearCost = 0.0;
913     int numberColumns = model->numberColumns();
914     int numberTotal = numberColumns;
915     double currentObj = 0.0;
916     for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
917          linearCost += cost[iColumn] * solution[iColumn];
918     }
919     if (!activated_ || !quadraticObjective_) {
920          currentObj = linearCost;
921          return currentObj;
922     }
923     assert (model);
924     const int * columnQuadratic = quadraticObjective_->getIndices();
925     const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
926     const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
927     const double * quadraticElement = quadraticObjective_->getElements();
928     double c = 0.0;
929     if (!scaling) {
930          if (!fullMatrix_) {
931               int iColumn;
932               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
933                    double valueI = solution[iColumn];
934                    CoinBigIndex j;
935                    for (j = columnQuadraticStart[iColumn];
936                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
937                         int jColumn = columnQuadratic[j];
938                         double valueJ = solution[jColumn];
939                         double elementValue = quadraticElement[j];
940                         if (iColumn != jColumn) {
941                              c += valueI * valueJ * elementValue;
942                         } else {
943                              c += 0.5 * valueI * valueI * elementValue;
944                         }
945                    }
946               }
947          } else {
948               // full matrix stored
949               int iColumn;
950               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
951                    double valueI = solution[iColumn];
952                    CoinBigIndex j;
953                    for (j = columnQuadraticStart[iColumn];
954                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
955                         int jColumn = columnQuadratic[j];
956                         double valueJ = solution[jColumn];
957                         double elementValue = quadraticElement[j];
958                         valueJ *= elementValue;
959                         c += valueI * valueJ;
960                    }
961               }
962               c *= 0.5;
963          }
964     } else {
965          // scaling
966          // for now only if half
967          assert (!fullMatrix_);
968          const double * columnScale = model->columnScale();
969          double direction = model->objectiveScale();
970          // direction is actually scale out not scale in
971          if (direction)
972               direction = 1.0 / direction;
973          if (!columnScale) {
974               for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
975                    double valueI = solution[iColumn];
976                    CoinBigIndex j;
977                    for (j = columnQuadraticStart[iColumn];
978                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
979                         int jColumn = columnQuadratic[j];
980                         double valueJ = solution[jColumn];
981                         double elementValue = quadraticElement[j];
982                         elementValue *= direction;
983                         if (iColumn != jColumn) {
984                              c += valueI * valueJ * elementValue;
985                         } else {
986                              c += 0.5 * valueI * valueI * elementValue;
987                         }
988                    }
989               }
990          } else {
991               // scaling
992               for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
993                    double valueI = solution[iColumn];
994                    double scaleI = columnScale[iColumn] * direction;
995                    CoinBigIndex j;
996                    for (j = columnQuadraticStart[iColumn];
997                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
998                         int jColumn = columnQuadratic[j];
999                         double valueJ = solution[jColumn];
1000                         double elementValue = quadraticElement[j];
1001                         elementValue *= scaleI * columnScale[jColumn];
1002                         if (iColumn != jColumn) {
1003                              c += valueI * valueJ * elementValue;
1004                         } else {
1005                              c += 0.5 * valueI * valueI * elementValue;
1006                         }
1007                    }
1008               }
1009          }
1010     }
1011     currentObj = c + linearCost;
1012     return currentObj;
1013}
1014// Scale objective
1015void
1016ClpQuadraticObjective::reallyScale(const double * columnScale)
1017{
1018     const int * columnQuadratic = quadraticObjective_->getIndices();
1019     const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
1020     const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
1021     double * quadraticElement = quadraticObjective_->getMutableElements();
1022     for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1023          double scaleI = columnScale[iColumn];
1024          objective_[iColumn] *= scaleI;
1025          CoinBigIndex j;
1026          for (j = columnQuadraticStart[iColumn];
1027                    j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
1028               int jColumn = columnQuadratic[j];
1029               quadraticElement[j] *= scaleI * columnScale[jColumn];
1030          }
1031     }
1032}
1033/* Given a zeroed array sets nonlinear columns to 1.
1034   Returns number of nonlinear columns
1035*/
1036int
1037ClpQuadraticObjective::markNonlinear(char * which)
1038{
1039     int iColumn;
1040     const int * columnQuadratic = quadraticObjective_->getIndices();
1041     const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
1042     const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
1043     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1044          CoinBigIndex j;
1045          for (j = columnQuadraticStart[iColumn];
1046                    j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
1047               int jColumn = columnQuadratic[j];
1048               which[jColumn] = 1;
1049               which[iColumn] = 1;
1050          }
1051     }
1052     int numberNonLinearColumns = 0;
1053     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1054          if(which[iColumn])
1055               numberNonLinearColumns++;
1056     }
1057     return numberNonLinearColumns;
1058}
Note: See TracBrowser for help on using the repository browser.