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

Last change on this file since 2030 was 1723, checked in by forrest, 9 years ago

out some printf statements

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 46.7 KB
Line 
1/* $Id: ClpQuadraticObjective.cpp 1723 2011-04-17 15:07:10Z forrest $ */
2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include "CoinPragma.hpp"
7#include "CoinHelperFunctions.hpp"
8#include "CoinIndexedVector.hpp"
9#include "ClpFactorization.hpp"
10#include "ClpSimplex.hpp"
11#include "ClpQuadraticObjective.hpp"
12//#############################################################################
13// Constructors / Destructor / Assignment
14//#############################################################################
15//-------------------------------------------------------------------
16// Default Constructor
17//-------------------------------------------------------------------
18ClpQuadraticObjective::ClpQuadraticObjective ()
19     : ClpObjective()
20{
21     type_ = 2;
22     objective_ = NULL;
23     quadraticObjective_ = NULL;
24     gradient_ = NULL;
25     numberColumns_ = 0;
26     numberExtendedColumns_ = 0;
27     activated_ = 0;
28     fullMatrix_ = false;
29}
30
31//-------------------------------------------------------------------
32// Useful Constructor
33//-------------------------------------------------------------------
34ClpQuadraticObjective::ClpQuadraticObjective (const double * objective ,
35          int numberColumns,
36          const CoinBigIndex * start,
37          const int * column, const double * element,
38          int numberExtendedColumns)
39     : ClpObjective()
40{
41     type_ = 2;
42     numberColumns_ = numberColumns;
43     if (numberExtendedColumns >= 0)
44          numberExtendedColumns_ = CoinMax(numberColumns_, numberExtendedColumns);
45     else
46          numberExtendedColumns_ = numberColumns_;
47     if (objective) {
48          objective_ = new double [numberExtendedColumns_];
49          CoinMemcpyN(objective, numberColumns_, objective_);
50          memset(objective_ + numberColumns_, 0, (numberExtendedColumns_ - numberColumns_)*sizeof(double));
51     } else {
52          objective_ = new double [numberExtendedColumns_];
53          memset(objective_, 0, numberExtendedColumns_ * sizeof(double));
54     }
55     if (start)
56          quadraticObjective_ = new CoinPackedMatrix(true, numberColumns, numberColumns,
57                    start[numberColumns], element, column, start, NULL);
58     else
59          quadraticObjective_ = NULL;
60     gradient_ = NULL;
61     activated_ = 1;
62     fullMatrix_ = false;
63}
64
65//-------------------------------------------------------------------
66// Copy constructor
67//-------------------------------------------------------------------
68ClpQuadraticObjective::ClpQuadraticObjective (const ClpQuadraticObjective & rhs,
69          int type)
70     : ClpObjective(rhs)
71{
72     numberColumns_ = rhs.numberColumns_;
73     numberExtendedColumns_ = rhs.numberExtendedColumns_;
74     fullMatrix_ = rhs.fullMatrix_;
75     if (rhs.objective_) {
76          objective_ = new double [numberExtendedColumns_];
77          CoinMemcpyN(rhs.objective_, numberExtendedColumns_, objective_);
78     } else {
79          objective_ = NULL;
80     }
81     if (rhs.gradient_) {
82          gradient_ = new double [numberExtendedColumns_];
83          CoinMemcpyN(rhs.gradient_, numberExtendedColumns_, gradient_);
84     } else {
85          gradient_ = NULL;
86     }
87     if (rhs.quadraticObjective_) {
88          // see what type of matrix wanted
89          if (type == 0) {
90               // just copy
91               quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
92          } else if (type == 1) {
93               // expand to full symmetric
94               fullMatrix_ = true;
95               const int * columnQuadratic1 = rhs.quadraticObjective_->getIndices();
96               const CoinBigIndex * columnQuadraticStart1 = rhs.quadraticObjective_->getVectorStarts();
97               const int * columnQuadraticLength1 = rhs.quadraticObjective_->getVectorLengths();
98               const double * quadraticElement1 = rhs.quadraticObjective_->getElements();
99               CoinBigIndex * columnQuadraticStart2 = new CoinBigIndex [numberExtendedColumns_+1];
100               int * columnQuadraticLength2 = new int [numberExtendedColumns_];
101               int iColumn;
102               int numberColumns = rhs.quadraticObjective_->getNumCols();
103               int numberBelow = 0;
104               int numberAbove = 0;
105               int numberDiagonal = 0;
106               CoinZeroN(columnQuadraticLength2, numberExtendedColumns_);
107               for (iColumn = 0; iColumn < numberColumns; iColumn++) {
108                    for (CoinBigIndex j = columnQuadraticStart1[iColumn];
109                              j < columnQuadraticStart1[iColumn] + columnQuadraticLength1[iColumn]; j++) {
110                         int jColumn = columnQuadratic1[j];
111                         if (jColumn > iColumn) {
112                              numberBelow++;
113                              columnQuadraticLength2[jColumn]++;
114                              columnQuadraticLength2[iColumn]++;
115                         } else if (jColumn == iColumn) {
116                              numberDiagonal++;
117                              columnQuadraticLength2[iColumn]++;
118                         } else {
119                              numberAbove++;
120                         }
121                    }
122               }
123               if (numberAbove > 0) {
124                    if (numberAbove == numberBelow) {
125                         // already done
126                         quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
127                         delete [] columnQuadraticStart2;
128                         delete [] columnQuadraticLength2;
129                    } else {
130                         printf("number above = %d, number below = %d, error\n",
131                                numberAbove, numberBelow);
132                         abort();
133                    }
134               } else {
135                    int numberElements = numberDiagonal + 2 * numberBelow;
136                    int * columnQuadratic2 = new int [numberElements];
137                    double * quadraticElement2 = new double [numberElements];
138                    columnQuadraticStart2[0] = 0;
139                    numberElements = 0;
140                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
141                         int n = columnQuadraticLength2[iColumn];
142                         columnQuadraticLength2[iColumn] = 0;
143                         numberElements += n;
144                         columnQuadraticStart2[iColumn+1] = numberElements;
145                    }
146                    for (iColumn = 0; iColumn < numberColumns; iColumn++) {
147                         for (CoinBigIndex j = columnQuadraticStart1[iColumn];
148                                   j < columnQuadraticStart1[iColumn] + columnQuadraticLength1[iColumn]; j++) {
149                              int jColumn = columnQuadratic1[j];
150                              if (jColumn > iColumn) {
151                                   // put in two places
152                                   CoinBigIndex put = columnQuadraticLength2[jColumn] + columnQuadraticStart2[jColumn];
153                                   columnQuadraticLength2[jColumn]++;
154                                   quadraticElement2[put] = quadraticElement1[j];
155                                   columnQuadratic2[put] = iColumn;
156                                   put = columnQuadraticLength2[iColumn] + columnQuadraticStart2[iColumn];
157                                   columnQuadraticLength2[iColumn]++;
158                                   quadraticElement2[put] = quadraticElement1[j];
159                                   columnQuadratic2[put] = jColumn;
160                              } else if (jColumn == iColumn) {
161                                   CoinBigIndex put = columnQuadraticLength2[iColumn] + columnQuadraticStart2[iColumn];
162                                   columnQuadraticLength2[iColumn]++;
163                                   quadraticElement2[put] = quadraticElement1[j];
164                                   columnQuadratic2[put] = iColumn;
165                              } else {
166                                   abort();
167                              }
168                         }
169                    }
170                    // Now create
171                    quadraticObjective_ =
172                         new CoinPackedMatrix (true,
173                                               rhs.numberExtendedColumns_,
174                                               rhs.numberExtendedColumns_,
175                                               numberElements,
176                                               quadraticElement2,
177                                               columnQuadratic2,
178                                               columnQuadraticStart2,
179                                               columnQuadraticLength2, 0.0, 0.0);
180                    delete [] columnQuadraticStart2;
181                    delete [] columnQuadraticLength2;
182                    delete [] columnQuadratic2;
183                    delete [] quadraticElement2;
184               }
185          } else {
186               fullMatrix_ = false;
187               abort(); // code when needed
188          }
189
190     } else {
191          quadraticObjective_ = NULL;
192     }
193}
194/* Subset constructor.  Duplicates are allowed
195   and order is as given.
196*/
197ClpQuadraticObjective::ClpQuadraticObjective (const ClpQuadraticObjective &rhs,
198          int numberColumns,
199          const int * whichColumn)
200     : ClpObjective(rhs)
201{
202     fullMatrix_ = rhs.fullMatrix_;
203     objective_ = NULL;
204     int extra = rhs.numberExtendedColumns_ - rhs.numberColumns_;
205     numberColumns_ = 0;
206     numberExtendedColumns_ = numberColumns + extra;
207     if (numberColumns > 0) {
208          // check valid lists
209          int numberBad = 0;
210          int i;
211          for (i = 0; i < numberColumns; i++)
212               if (whichColumn[i] < 0 || whichColumn[i] >= rhs.numberColumns_)
213                    numberBad++;
214          if (numberBad)
215               throw CoinError("bad column list", "subset constructor",
216                               "ClpQuadraticObjective");
217          numberColumns_ = numberColumns;
218          objective_ = new double[numberExtendedColumns_];
219          for (i = 0; i < numberColumns_; i++)
220               objective_[i] = rhs.objective_[whichColumn[i]];
221          CoinMemcpyN(rhs.objective_ + rhs.numberColumns_,      (numberExtendedColumns_ - numberColumns_),
222                      objective_ + numberColumns_);
223          if (rhs.gradient_) {
224               gradient_ = new double[numberExtendedColumns_];
225               for (i = 0; i < numberColumns_; i++)
226                    gradient_[i] = rhs.gradient_[whichColumn[i]];
227               CoinMemcpyN(rhs.gradient_ + rhs.numberColumns_,  (numberExtendedColumns_ - numberColumns_),
228                           gradient_ + numberColumns_);
229          } else {
230               gradient_ = NULL;
231          }
232     } else {
233          gradient_ = NULL;
234          objective_ = NULL;
235     }
236     if (rhs.quadraticObjective_) {
237          quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_,
238                    numberColumns, whichColumn,
239                    numberColumns, whichColumn);
240     } else {
241          quadraticObjective_ = NULL;
242     }
243}
244
245
246//-------------------------------------------------------------------
247// Destructor
248//-------------------------------------------------------------------
249ClpQuadraticObjective::~ClpQuadraticObjective ()
250{
251     delete [] objective_;
252     delete [] gradient_;
253     delete quadraticObjective_;
254}
255
256//----------------------------------------------------------------
257// Assignment operator
258//-------------------------------------------------------------------
259ClpQuadraticObjective &
260ClpQuadraticObjective::operator=(const ClpQuadraticObjective& rhs)
261{
262     if (this != &rhs) {
263          fullMatrix_ = rhs.fullMatrix_;
264          delete quadraticObjective_;
265          quadraticObjective_ = NULL;
266          delete [] objective_;
267          delete [] gradient_;
268          ClpObjective::operator=(rhs);
269          numberColumns_ = rhs.numberColumns_;
270          numberExtendedColumns_ = rhs.numberExtendedColumns_;
271          if (rhs.objective_) {
272               objective_ = new double [numberExtendedColumns_];
273               CoinMemcpyN(rhs.objective_, numberExtendedColumns_, objective_);
274          } else {
275               objective_ = NULL;
276          }
277          if (rhs.gradient_) {
278               gradient_ = new double [numberExtendedColumns_];
279               CoinMemcpyN(rhs.gradient_, numberExtendedColumns_, gradient_);
280          } else {
281               gradient_ = NULL;
282          }
283          if (rhs.quadraticObjective_) {
284               quadraticObjective_ = new CoinPackedMatrix(*rhs.quadraticObjective_);
285          } else {
286               quadraticObjective_ = NULL;
287          }
288     }
289     return *this;
290}
291
292// Returns gradient
293double *
294ClpQuadraticObjective::gradient(const ClpSimplex * model,
295                                const double * solution, double & offset, bool refresh,
296                                int includeLinear)
297{
298     offset = 0.0;
299     bool scaling = false;
300     if (model && (model->rowScale() ||
301                   model->objectiveScale() != 1.0 || model->optimizationDirection() != 1.0))
302          scaling = true;
303     const double * cost = NULL;
304     if (model)
305          cost = model->costRegion();
306     if (!cost) {
307          // not in solve
308          cost = objective_;
309          scaling = false;
310     }
311     if (!scaling) {
312          if (!quadraticObjective_ || !solution || !activated_) {
313               return objective_;
314          } else {
315               if (refresh || !gradient_) {
316                    if (!gradient_)
317                         gradient_ = new double[numberExtendedColumns_];
318                    const int * columnQuadratic = quadraticObjective_->getIndices();
319                    const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
320                    const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
321                    const double * quadraticElement = quadraticObjective_->getElements();
322                    offset = 0.0;
323                    // use current linear cost region
324                    if (includeLinear == 1)
325                         CoinMemcpyN(cost, numberExtendedColumns_, gradient_);
326                    else if (includeLinear == 2)
327                         CoinMemcpyN(objective_, numberExtendedColumns_, gradient_);
328                    else
329                         memset(gradient_, 0, numberExtendedColumns_ * sizeof(double));
330                    if (activated_) {
331                         if (!fullMatrix_) {
332                              int iColumn;
333                              for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
334                                   double valueI = solution[iColumn];
335                                   CoinBigIndex j;
336                                   for (j = columnQuadraticStart[iColumn];
337                                             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
338                                        int jColumn = columnQuadratic[j];
339                                        double valueJ = solution[jColumn];
340                                        double elementValue = quadraticElement[j];
341                                        if (iColumn != jColumn) {
342                                             offset += valueI * valueJ * elementValue;
343                                             //if (fabs(valueI*valueJ*elementValue)>1.0e-12)
344                                             //printf("%d %d %g %g %g -> %g\n",
345                                             //       iColumn,jColumn,valueI,valueJ,elementValue,
346                                             //       valueI*valueJ*elementValue);
347                                             double gradientI = valueJ * elementValue;
348                                             double gradientJ = valueI * elementValue;
349                                             gradient_[iColumn] += gradientI;
350                                             gradient_[jColumn] += gradientJ;
351                                        } else {
352                                             offset += 0.5 * valueI * valueI * elementValue;
353                                             //if (fabs(valueI*valueI*elementValue)>1.0e-12)
354                                             //printf("XX %d %g %g -> %g\n",
355                                             //       iColumn,valueI,elementValue,
356                                             //       0.5*valueI*valueI*elementValue);
357                                             double gradientI = valueI * elementValue;
358                                             gradient_[iColumn] += gradientI;
359                                        }
360                                   }
361                              }
362                         } else {
363                              // full matrix
364                              int iColumn;
365                              offset *= 2.0;
366                              for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
367                                   CoinBigIndex j;
368                                   double value = 0.0;
369                                   double current = gradient_[iColumn];
370                                   for (j = columnQuadraticStart[iColumn];
371                                             j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
372                                        int jColumn = columnQuadratic[j];
373                                        double valueJ = solution[jColumn] * quadraticElement[j];
374                                        value += valueJ;
375                                   }
376                                   offset += value * solution[iColumn];
377                                   gradient_[iColumn] = current + value;
378                              }
379                              offset *= 0.5;
380                         }
381                    }
382               }
383               if (model)
384                    offset *= model->optimizationDirection() * model->objectiveScale();
385               return gradient_;
386          }
387     } else {
388          // do scaling
389          assert (solution);
390          // for now only if half
391          assert (!fullMatrix_);
392          if (refresh || !gradient_) {
393               if (!gradient_)
394                    gradient_ = new double[numberExtendedColumns_];
395               double direction = model->optimizationDirection() * model->objectiveScale();
396               // direction is actually scale out not scale in
397               //if (direction)
398               //direction = 1.0/direction;
399               const int * columnQuadratic = quadraticObjective_->getIndices();
400               const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
401               const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
402               const double * quadraticElement = quadraticObjective_->getElements();
403               int iColumn;
404               const double * columnScale = model->columnScale();
405               // use current linear cost region (already scaled)
406               if (includeLinear == 1) {
407                    CoinMemcpyN(model->costRegion(), numberExtendedColumns_, gradient_);
408               }        else if (includeLinear == 2) {
409                    memset(gradient_ + numberColumns_, 0, (numberExtendedColumns_ - numberColumns_)*sizeof(double));
410                    if (!columnScale) {
411                         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
412                              gradient_[iColumn] = objective_[iColumn] * direction;
413                         }
414                    } else {
415                         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
416                              gradient_[iColumn] = objective_[iColumn] * direction * columnScale[iColumn];
417                         }
418                    }
419               } else {
420                    memset(gradient_, 0, numberExtendedColumns_ * sizeof(double));
421               }
422               if (!columnScale) {
423                    if (activated_) {
424                         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
425                              double valueI = solution[iColumn];
426                              CoinBigIndex j;
427                              for (j = columnQuadraticStart[iColumn];
428                                        j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
429                                   int jColumn = columnQuadratic[j];
430                                   double valueJ = solution[jColumn];
431                                   double elementValue = quadraticElement[j];
432                                   elementValue *= direction;
433                                   if (iColumn != jColumn) {
434                                        offset += valueI * valueJ * elementValue;
435                                        double gradientI = valueJ * elementValue;
436                                        double gradientJ = valueI * elementValue;
437                                        gradient_[iColumn] += gradientI;
438                                        gradient_[jColumn] += gradientJ;
439                                   } else {
440                                        offset += 0.5 * valueI * valueI * elementValue;
441                                        double gradientI = valueI * elementValue;
442                                        gradient_[iColumn] += gradientI;
443                                   }
444                              }
445                         }
446                    }
447               } else {
448                    // scaling
449                    if (activated_) {
450                         for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
451                              double valueI = solution[iColumn];
452                              double scaleI = columnScale[iColumn] * direction;
453                              CoinBigIndex j;
454                              for (j = columnQuadraticStart[iColumn];
455                                        j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
456                                   int jColumn = columnQuadratic[j];
457                                   double valueJ = solution[jColumn];
458                                   double elementValue = quadraticElement[j];
459                                   double scaleJ = columnScale[jColumn];
460                                   elementValue *= scaleI * scaleJ;
461                                   if (iColumn != jColumn) {
462                                        offset += valueI * valueJ * elementValue;
463                                        double gradientI = valueJ * elementValue;
464                                        double gradientJ = valueI * elementValue;
465                                        gradient_[iColumn] += gradientI;
466                                        gradient_[jColumn] += gradientJ;
467                                   } else {
468                                        offset += 0.5 * valueI * valueI * elementValue;
469                                        double gradientI = valueI * elementValue;
470                                        gradient_[iColumn] += gradientI;
471                                   }
472                              }
473                         }
474                    }
475               }
476          }
477          if (model)
478               offset *= model->optimizationDirection();
479          return gradient_;
480     }
481}
482
483//-------------------------------------------------------------------
484// Clone
485//-------------------------------------------------------------------
486ClpObjective * ClpQuadraticObjective::clone() const
487{
488     return new ClpQuadraticObjective(*this);
489}
490/* Subset clone.  Duplicates are allowed
491   and order is as given.
492*/
493ClpObjective *
494ClpQuadraticObjective::subsetClone (int numberColumns,
495                                    const int * whichColumns) const
496{
497     return new ClpQuadraticObjective(*this, numberColumns, whichColumns);
498}
499// Resize objective
500void
501ClpQuadraticObjective::resize(int newNumberColumns)
502{
503     if (numberColumns_ != newNumberColumns) {
504          int newExtended = newNumberColumns + (numberExtendedColumns_ - numberColumns_);
505          int i;
506          double * newArray = new double[newExtended];
507          if (objective_)
508               CoinMemcpyN(objective_,  CoinMin(newExtended, numberExtendedColumns_), newArray);
509          delete [] objective_;
510          objective_ = newArray;
511          for (i = numberColumns_; i < newNumberColumns; i++)
512               objective_[i] = 0.0;
513          if (gradient_) {
514               newArray = new double[newExtended];
515               if (gradient_)
516                    CoinMemcpyN(gradient_,      CoinMin(newExtended, numberExtendedColumns_), newArray);
517               delete [] gradient_;
518               gradient_ = newArray;
519               for (i = numberColumns_; i < newNumberColumns; i++)
520                    gradient_[i] = 0.0;
521          }
522          if (quadraticObjective_) {
523               if (newNumberColumns < numberColumns_) {
524                    int * which = new int[numberColumns_-newNumberColumns];
525                    int i;
526                    for (i = newNumberColumns; i < numberColumns_; i++)
527                         which[i-newNumberColumns] = i;
528                    quadraticObjective_->deleteRows(numberColumns_ - newNumberColumns, which);
529                    quadraticObjective_->deleteCols(numberColumns_ - newNumberColumns, which);
530                    delete [] which;
531               } else {
532                    quadraticObjective_->setDimensions(newNumberColumns, newNumberColumns);
533               }
534          }
535          numberColumns_ = newNumberColumns;
536          numberExtendedColumns_ = newExtended;
537     }
538
539}
540// Delete columns in  objective
541void
542ClpQuadraticObjective::deleteSome(int numberToDelete, const int * which)
543{
544     int newNumberColumns = numberColumns_ - numberToDelete;
545     int newExtended = numberExtendedColumns_ - numberToDelete;
546     if (objective_) {
547          int i ;
548          char * deleted = new char[numberColumns_];
549          int numberDeleted = 0;
550          memset(deleted, 0, numberColumns_ * sizeof(char));
551          for (i = 0; i < numberToDelete; i++) {
552               int j = which[i];
553               if (j >= 0 && j < numberColumns_ && !deleted[j]) {
554                    numberDeleted++;
555                    deleted[j] = 1;
556               }
557          }
558          newNumberColumns = numberColumns_ - numberDeleted;
559          newExtended = numberExtendedColumns_ - numberDeleted;
560          double * newArray = new double[newExtended];
561          int put = 0;
562          for (i = 0; i < numberColumns_; i++) {
563               if (!deleted[i]) {
564                    newArray[put++] = objective_[i];
565               }
566          }
567          delete [] objective_;
568          objective_ = newArray;
569          delete [] deleted;
570          CoinMemcpyN(objective_ + numberColumns_,      (numberExtendedColumns_ - numberColumns_),
571                      objective_ + newNumberColumns);
572     }
573     if (gradient_) {
574          int i ;
575          char * deleted = new char[numberColumns_];
576          int numberDeleted = 0;
577          memset(deleted, 0, numberColumns_ * sizeof(char));
578          for (i = 0; i < numberToDelete; i++) {
579               int j = which[i];
580               if (j >= 0 && j < numberColumns_ && !deleted[j]) {
581                    numberDeleted++;
582                    deleted[j] = 1;
583               }
584          }
585          newNumberColumns = numberColumns_ - numberDeleted;
586          newExtended = numberExtendedColumns_ - numberDeleted;
587          double * newArray = new double[newExtended];
588          int put = 0;
589          for (i = 0; i < numberColumns_; i++) {
590               if (!deleted[i]) {
591                    newArray[put++] = gradient_[i];
592               }
593          }
594          delete [] gradient_;
595          gradient_ = newArray;
596          delete [] deleted;
597          CoinMemcpyN(gradient_ + numberColumns_,       (numberExtendedColumns_ - numberColumns_),
598                      gradient_ + newNumberColumns);
599     }
600     numberColumns_ = newNumberColumns;
601     numberExtendedColumns_ = newExtended;
602     if (quadraticObjective_) {
603          quadraticObjective_->deleteCols(numberToDelete, which);
604          quadraticObjective_->deleteRows(numberToDelete, which);
605     }
606}
607
608// Load up quadratic objective
609void
610ClpQuadraticObjective::loadQuadraticObjective(const int numberColumns, const CoinBigIndex * start,
611          const int * column, const double * element, int numberExtended)
612{
613     fullMatrix_ = false;
614     delete quadraticObjective_;
615     quadraticObjective_ = new CoinPackedMatrix(true, numberColumns, numberColumns,
616               start[numberColumns], element, column, start, NULL);
617     numberColumns_ = numberColumns;
618     if (numberExtended > numberExtendedColumns_) {
619          if (objective_) {
620               // make correct size
621               double * newArray = new double[numberExtended];
622               CoinMemcpyN(objective_, numberColumns_, newArray);
623               delete [] objective_;
624               objective_ = newArray;
625               memset(objective_ + numberColumns_, 0, (numberExtended - numberColumns_)*sizeof(double));
626          }
627          if (gradient_) {
628               // make correct size
629               double * newArray = new double[numberExtended];
630               CoinMemcpyN(gradient_, numberColumns_, newArray);
631               delete [] gradient_;
632               gradient_ = newArray;
633               memset(gradient_ + numberColumns_, 0, (numberExtended - numberColumns_)*sizeof(double));
634          }
635          numberExtendedColumns_ = numberExtended;
636     } else {
637          numberExtendedColumns_ = numberColumns_;
638     }
639}
640void
641ClpQuadraticObjective::loadQuadraticObjective (  const CoinPackedMatrix& matrix)
642{
643     delete quadraticObjective_;
644     quadraticObjective_ = new CoinPackedMatrix(matrix);
645}
646// Get rid of quadratic objective
647void
648ClpQuadraticObjective::deleteQuadraticObjective()
649{
650     delete quadraticObjective_;
651     quadraticObjective_ = NULL;
652}
653/* Returns reduced gradient.Returns an offset (to be added to current one).
654 */
655double
656ClpQuadraticObjective::reducedGradient(ClpSimplex * model, double * region,
657                                       bool useFeasibleCosts)
658{
659     int numberRows = model->numberRows();
660     int numberColumns = model->numberColumns();
661
662     //work space
663     CoinIndexedVector  * workSpace = model->rowArray(0);
664
665     CoinIndexedVector arrayVector;
666     arrayVector.reserve(numberRows + 1);
667
668     int iRow;
669#ifdef CLP_DEBUG
670     workSpace->checkClear();
671#endif
672     double * array = arrayVector.denseVector();
673     int * index = arrayVector.getIndices();
674     int number = 0;
675     const double * costNow = gradient(model, model->solutionRegion(), offset_,
676                                       true, useFeasibleCosts ? 2 : 1);
677     double * cost = model->costRegion();
678     const int * pivotVariable = model->pivotVariable();
679     for (iRow = 0; iRow < numberRows; iRow++) {
680          int iPivot = pivotVariable[iRow];
681          double value;
682          if (iPivot < numberColumns)
683               value = costNow[iPivot];
684          else if (!useFeasibleCosts)
685               value = cost[iPivot];
686          else
687               value = 0.0;
688          if (value) {
689               array[iRow] = value;
690               index[number++] = iRow;
691          }
692     }
693     arrayVector.setNumElements(number);
694
695     // Btran basic costs
696     model->factorization()->updateColumnTranspose(workSpace, &arrayVector);
697     double * work = workSpace->denseVector();
698     ClpFillN(work, numberRows, 0.0);
699     // now look at dual solution
700     double * rowReducedCost = region + numberColumns;
701     double * dual = rowReducedCost;
702     const double * rowCost = cost + numberColumns;
703     for (iRow = 0; iRow < numberRows; iRow++) {
704          dual[iRow] = array[iRow];
705     }
706     double * dj = region;
707     ClpDisjointCopyN(costNow, numberColumns, dj);
708
709     model->transposeTimes(-1.0, dual, dj);
710     for (iRow = 0; iRow < numberRows; iRow++) {
711          // slack
712          double value = dual[iRow];
713          value += rowCost[iRow];
714          rowReducedCost[iRow] = value;
715     }
716     return offset_;
717}
718/* Returns step length which gives minimum of objective for
719   solution + theta * change vector up to maximum theta.
720
721   arrays are numberColumns+numberRows
722*/
723double
724ClpQuadraticObjective::stepLength(ClpSimplex * model,
725                                  const double * solution,
726                                  const double * change,
727                                  double maximumTheta,
728                                  double & currentObj,
729                                  double & predictedObj,
730                                  double & thetaObj)
731{
732     const double * cost = model->costRegion();
733     bool inSolve = true;
734     if (!cost) {
735          // not in solve
736          cost = objective_;
737          inSolve = false;
738     }
739     double delta = 0.0;
740     double linearCost = 0.0;
741     int numberRows = model->numberRows();
742     int numberColumns = model->numberColumns();
743     int numberTotal = numberColumns;
744     if (inSolve)
745          numberTotal += numberRows;
746     currentObj = 0.0;
747     thetaObj = 0.0;
748     for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
749          delta += cost[iColumn] * change[iColumn];
750          linearCost += cost[iColumn] * solution[iColumn];
751     }
752     if (!activated_ || !quadraticObjective_) {
753          currentObj = linearCost;
754          thetaObj = currentObj + delta * maximumTheta;
755          if (delta < 0.0) {
756               return maximumTheta;
757          } else {
758            COIN_DETAIL_PRINT(printf("odd linear direction %g\n", delta));
759               return 0.0;
760          }
761     }
762     assert (model);
763     bool scaling = false;
764     if ((model->rowScale() ||
765               model->objectiveScale() != 1.0 || model->optimizationDirection() != 1.0) && inSolve)
766          scaling = true;
767     const int * columnQuadratic = quadraticObjective_->getIndices();
768     const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
769     const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
770     const double * quadraticElement = quadraticObjective_->getElements();
771     double a = 0.0;
772     double b = delta;
773     double c = 0.0;
774     if (!scaling) {
775          if (!fullMatrix_) {
776               int iColumn;
777               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
778                    double valueI = solution[iColumn];
779                    double changeI = change[iColumn];
780                    CoinBigIndex j;
781                    for (j = columnQuadraticStart[iColumn];
782                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
783                         int jColumn = columnQuadratic[j];
784                         double valueJ = solution[jColumn];
785                         double changeJ = change[jColumn];
786                         double elementValue = quadraticElement[j];
787                         if (iColumn != jColumn) {
788                              a += changeI * changeJ * elementValue;
789                              b += (changeI * valueJ + changeJ * valueI) * elementValue;
790                              c += valueI * valueJ * elementValue;
791                         } else {
792                              a += 0.5 * changeI * changeI * elementValue;
793                              b += changeI * valueI * elementValue;
794                              c += 0.5 * valueI * valueI * elementValue;
795                         }
796                    }
797               }
798          } else {
799               // full matrix stored
800               int iColumn;
801               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
802                    double valueI = solution[iColumn];
803                    double changeI = change[iColumn];
804                    CoinBigIndex j;
805                    for (j = columnQuadraticStart[iColumn];
806                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
807                         int jColumn = columnQuadratic[j];
808                         double valueJ = solution[jColumn];
809                         double changeJ = change[jColumn];
810                         double elementValue = quadraticElement[j];
811                         valueJ *= elementValue;
812                         a += changeI * changeJ * elementValue;
813                         b += changeI * valueJ;
814                         c += valueI * valueJ;
815                    }
816               }
817               a *= 0.5;
818               c *= 0.5;
819          }
820     } else {
821          // scaling
822          // for now only if half
823          assert (!fullMatrix_);
824          const double * columnScale = model->columnScale();
825          double direction = model->optimizationDirection() * model->objectiveScale();
826          // direction is actually scale out not scale in
827          if (direction)
828               direction = 1.0 / direction;
829          if (!columnScale) {
830               for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
831                    double valueI = solution[iColumn];
832                    double changeI = change[iColumn];
833                    CoinBigIndex j;
834                    for (j = columnQuadraticStart[iColumn];
835                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
836                         int jColumn = columnQuadratic[j];
837                         double valueJ = solution[jColumn];
838                         double changeJ = change[jColumn];
839                         double elementValue = quadraticElement[j];
840                         elementValue *= direction;
841                         if (iColumn != jColumn) {
842                              a += changeI * changeJ * elementValue;
843                              b += (changeI * valueJ + changeJ * valueI) * elementValue;
844                              c += valueI * valueJ * elementValue;
845                         } else {
846                              a += 0.5 * changeI * changeI * elementValue;
847                              b += changeI * valueI * elementValue;
848                              c += 0.5 * valueI * valueI * elementValue;
849                         }
850                    }
851               }
852          } else {
853               // scaling
854               for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
855                    double valueI = solution[iColumn];
856                    double changeI = change[iColumn];
857                    double scaleI = columnScale[iColumn] * direction;
858                    CoinBigIndex j;
859                    for (j = columnQuadraticStart[iColumn];
860                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
861                         int jColumn = columnQuadratic[j];
862                         double valueJ = solution[jColumn];
863                         double changeJ = change[jColumn];
864                         double elementValue = quadraticElement[j];
865                         elementValue *= scaleI * columnScale[jColumn];
866                         if (iColumn != jColumn) {
867                              a += changeI * changeJ * elementValue;
868                              b += (changeI * valueJ + changeJ * valueI) * elementValue;
869                              c += valueI * valueJ * elementValue;
870                         } else {
871                              a += 0.5 * changeI * changeI * elementValue;
872                              b += changeI * valueI * elementValue;
873                              c += 0.5 * valueI * valueI * elementValue;
874                         }
875                    }
876               }
877          }
878     }
879     double theta;
880     //printf("Current cost %g\n",c+linearCost);
881     currentObj = c + linearCost;
882     thetaObj = currentObj + a * maximumTheta * maximumTheta + b * maximumTheta;
883     // minimize a*x*x + b*x + c
884     if (a <= 0.0) {
885          theta = maximumTheta;
886     } else {
887          theta = -0.5 * b / a;
888     }
889     predictedObj = currentObj + a * theta * theta + b * theta;
890     if (b > 0.0) {
891          if (model->messageHandler()->logLevel() & 32)
892               printf("a %g b %g c %g => %g\n", a, b, c, theta);
893          b = 0.0;
894     }
895     return CoinMin(theta, maximumTheta);
896}
897// Return objective value (without any ClpModel offset) (model may be NULL)
898double
899ClpQuadraticObjective::objectiveValue(const ClpSimplex * model, const double * solution) const
900{
901     bool scaling = false;
902     if (model && (model->rowScale() ||
903                   model->objectiveScale() != 1.0))
904          scaling = true;
905     const double * cost = NULL;
906     if (model)
907          cost = model->costRegion();
908     if (!cost) {
909          // not in solve
910          cost = objective_;
911          scaling = false;
912     }
913     double linearCost = 0.0;
914     int numberColumns = model->numberColumns();
915     int numberTotal = numberColumns;
916     double currentObj = 0.0;
917     for (int iColumn = 0; iColumn < numberTotal; iColumn++) {
918          linearCost += cost[iColumn] * solution[iColumn];
919     }
920     if (!activated_ || !quadraticObjective_) {
921          currentObj = linearCost;
922          return currentObj;
923     }
924     assert (model);
925     const int * columnQuadratic = quadraticObjective_->getIndices();
926     const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
927     const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
928     const double * quadraticElement = quadraticObjective_->getElements();
929     double c = 0.0;
930     if (!scaling) {
931          if (!fullMatrix_) {
932               int iColumn;
933               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
934                    double valueI = solution[iColumn];
935                    CoinBigIndex j;
936                    for (j = columnQuadraticStart[iColumn];
937                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
938                         int jColumn = columnQuadratic[j];
939                         double valueJ = solution[jColumn];
940                         double elementValue = quadraticElement[j];
941                         if (iColumn != jColumn) {
942                              c += valueI * valueJ * elementValue;
943                         } else {
944                              c += 0.5 * valueI * valueI * elementValue;
945                         }
946                    }
947               }
948          } else {
949               // full matrix stored
950               int iColumn;
951               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
952                    double valueI = solution[iColumn];
953                    CoinBigIndex j;
954                    for (j = columnQuadraticStart[iColumn];
955                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
956                         int jColumn = columnQuadratic[j];
957                         double valueJ = solution[jColumn];
958                         double elementValue = quadraticElement[j];
959                         valueJ *= elementValue;
960                         c += valueI * valueJ;
961                    }
962               }
963               c *= 0.5;
964          }
965     } else {
966          // scaling
967          // for now only if half
968          assert (!fullMatrix_);
969          const double * columnScale = model->columnScale();
970          double direction = model->objectiveScale();
971          // direction is actually scale out not scale in
972          if (direction)
973               direction = 1.0 / direction;
974          if (!columnScale) {
975               for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
976                    double valueI = solution[iColumn];
977                    CoinBigIndex j;
978                    for (j = columnQuadraticStart[iColumn];
979                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
980                         int jColumn = columnQuadratic[j];
981                         double valueJ = solution[jColumn];
982                         double elementValue = quadraticElement[j];
983                         elementValue *= direction;
984                         if (iColumn != jColumn) {
985                              c += valueI * valueJ * elementValue;
986                         } else {
987                              c += 0.5 * valueI * valueI * elementValue;
988                         }
989                    }
990               }
991          } else {
992               // scaling
993               for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
994                    double valueI = solution[iColumn];
995                    double scaleI = columnScale[iColumn] * direction;
996                    CoinBigIndex j;
997                    for (j = columnQuadraticStart[iColumn];
998                              j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
999                         int jColumn = columnQuadratic[j];
1000                         double valueJ = solution[jColumn];
1001                         double elementValue = quadraticElement[j];
1002                         elementValue *= scaleI * columnScale[jColumn];
1003                         if (iColumn != jColumn) {
1004                              c += valueI * valueJ * elementValue;
1005                         } else {
1006                              c += 0.5 * valueI * valueI * elementValue;
1007                         }
1008                    }
1009               }
1010          }
1011     }
1012     currentObj = c + linearCost;
1013     return currentObj;
1014}
1015// Scale objective
1016void
1017ClpQuadraticObjective::reallyScale(const double * columnScale)
1018{
1019     const int * columnQuadratic = quadraticObjective_->getIndices();
1020     const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
1021     const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
1022     double * quadraticElement = quadraticObjective_->getMutableElements();
1023     for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1024          double scaleI = columnScale[iColumn];
1025          objective_[iColumn] *= scaleI;
1026          CoinBigIndex j;
1027          for (j = columnQuadraticStart[iColumn];
1028                    j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
1029               int jColumn = columnQuadratic[j];
1030               quadraticElement[j] *= scaleI * columnScale[jColumn];
1031          }
1032     }
1033}
1034/* Given a zeroed array sets nonlinear columns to 1.
1035   Returns number of nonlinear columns
1036*/
1037int
1038ClpQuadraticObjective::markNonlinear(char * which)
1039{
1040     int iColumn;
1041     const int * columnQuadratic = quadraticObjective_->getIndices();
1042     const CoinBigIndex * columnQuadraticStart = quadraticObjective_->getVectorStarts();
1043     const int * columnQuadraticLength = quadraticObjective_->getVectorLengths();
1044     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1045          CoinBigIndex j;
1046          for (j = columnQuadraticStart[iColumn];
1047                    j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
1048               int jColumn = columnQuadratic[j];
1049               which[jColumn] = 1;
1050               which[iColumn] = 1;
1051          }
1052     }
1053     int numberNonLinearColumns = 0;
1054     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1055          if(which[iColumn])
1056               numberNonLinearColumns++;
1057     }
1058     return numberNonLinearColumns;
1059}
Note: See TracBrowser for help on using the repository browser.