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

Last change on this file since 2030 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

  • Property svn:keywords set to Id
File size: 10.4 KB
Line 
1/* $Id: ClpConstraintQuadratic.cpp 1665 2011-01-04 17:55:54Z forrest $ */
2// Copyright (C) 2007, 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 "ClpSimplex.hpp"
10#include "ClpConstraintQuadratic.hpp"
11#include "CoinSort.hpp"
12//#############################################################################
13// Constructors / Destructor / Assignment
14//#############################################################################
15
16//-------------------------------------------------------------------
17// Default Constructor
18//-------------------------------------------------------------------
19ClpConstraintQuadratic::ClpConstraintQuadratic ()
20     : ClpConstraint()
21{
22     type_ = 0;
23     start_ = NULL;
24     column_ = NULL;
25     coefficient_ = NULL;
26     numberColumns_ = 0;
27     numberCoefficients_ = 0;
28     numberQuadraticColumns_ = 0;
29}
30
31//-------------------------------------------------------------------
32// Useful Constructor
33//-------------------------------------------------------------------
34ClpConstraintQuadratic::ClpConstraintQuadratic (int row, int numberQuadraticColumns ,
35          int numberColumns, const CoinBigIndex * start,
36          const int * column, const double * coefficient)
37     : ClpConstraint()
38{
39     type_ = 0;
40     rowNumber_ = row;
41     numberColumns_ = numberColumns;
42     numberQuadraticColumns_ = numberQuadraticColumns;
43     start_ = CoinCopyOfArray(start, numberQuadraticColumns + 1);
44     int numberElements = start_[numberQuadraticColumns_];
45     column_ = CoinCopyOfArray(column, numberElements);
46     coefficient_ = CoinCopyOfArray(coefficient, numberElements);
47     char * mark = new char [numberQuadraticColumns_];
48     memset(mark, 0, numberQuadraticColumns_);
49     int iColumn;
50     for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
51          CoinBigIndex j;
52          for (j = start_[iColumn]; j < start_[iColumn+1]; j++) {
53               int jColumn = column_[j];
54               if (jColumn >= 0) {
55                    assert (jColumn < numberQuadraticColumns_);
56                    mark[jColumn] = 1;
57               }
58               mark[iColumn] = 1;
59          }
60     }
61     numberCoefficients_ = 0;
62     for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
63          if (mark[iColumn])
64               numberCoefficients_++;
65     }
66     delete [] mark;
67}
68
69//-------------------------------------------------------------------
70// Copy constructor
71//-------------------------------------------------------------------
72ClpConstraintQuadratic::ClpConstraintQuadratic (const ClpConstraintQuadratic & rhs)
73     : ClpConstraint(rhs)
74{
75     numberColumns_ = rhs.numberColumns_;
76     numberCoefficients_ = rhs.numberCoefficients_;
77     numberQuadraticColumns_ = rhs.numberQuadraticColumns_;
78     start_ = CoinCopyOfArray(rhs.start_, numberQuadraticColumns_ + 1);
79     int numberElements = start_[numberQuadraticColumns_];
80     column_ = CoinCopyOfArray(rhs.column_, numberElements);
81     coefficient_ = CoinCopyOfArray(rhs.coefficient_, numberElements);
82}
83
84
85//-------------------------------------------------------------------
86// Destructor
87//-------------------------------------------------------------------
88ClpConstraintQuadratic::~ClpConstraintQuadratic ()
89{
90     delete [] start_;
91     delete [] column_;
92     delete [] coefficient_;
93}
94
95//----------------------------------------------------------------
96// Assignment operator
97//-------------------------------------------------------------------
98ClpConstraintQuadratic &
99ClpConstraintQuadratic::operator=(const ClpConstraintQuadratic& rhs)
100{
101     if (this != &rhs) {
102          delete [] start_;
103          delete [] column_;
104          delete [] coefficient_;
105          numberColumns_ = rhs.numberColumns_;
106          numberCoefficients_ = rhs.numberCoefficients_;
107          numberQuadraticColumns_ = rhs.numberQuadraticColumns_;
108          start_ = CoinCopyOfArray(rhs.start_, numberQuadraticColumns_ + 1);
109          int numberElements = start_[numberQuadraticColumns_];
110          column_ = CoinCopyOfArray(rhs.column_, numberElements);
111          coefficient_ = CoinCopyOfArray(rhs.coefficient_, numberElements);
112     }
113     return *this;
114}
115//-------------------------------------------------------------------
116// Clone
117//-------------------------------------------------------------------
118ClpConstraint * ClpConstraintQuadratic::clone() const
119{
120     return new ClpConstraintQuadratic(*this);
121}
122
123// Returns gradient
124int
125ClpConstraintQuadratic::gradient(const ClpSimplex * model,
126                                 const double * solution,
127                                 double * gradient,
128                                 double & functionValue,
129                                 double & offset,
130                                 bool useScaling,
131                                 bool refresh) const
132{
133     if (refresh || !lastGradient_) {
134          offset_ = 0.0;
135          functionValue_ = 0.0;
136          if (!lastGradient_)
137               lastGradient_ = new double[numberColumns_];
138          CoinZeroN(lastGradient_, numberColumns_);
139          bool scaling = (model && model->rowScale() && useScaling);
140          if (!scaling) {
141               int iColumn;
142               for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
143                    double valueI = solution[iColumn];
144                    CoinBigIndex j;
145                    for (j = start_[iColumn]; j < start_[iColumn+1]; j++) {
146                         int jColumn = column_[j];
147                         if (jColumn >= 0) {
148                              double valueJ = solution[jColumn];
149                              double elementValue = coefficient_[j];
150                              if (iColumn != jColumn) {
151                                   offset_ -= valueI * valueJ * elementValue;
152                                   double gradientI = valueJ * elementValue;
153                                   double gradientJ = valueI * elementValue;
154                                   lastGradient_[iColumn] += gradientI;
155                                   lastGradient_[jColumn] += gradientJ;
156                              } else {
157                                   offset_ -= 0.5 * valueI * valueI * elementValue;
158                                   double gradientI = valueI * elementValue;
159                                   lastGradient_[iColumn] += gradientI;
160                              }
161                         } else {
162                              // linear part
163                              lastGradient_[iColumn] += coefficient_[j];
164                              functionValue_ += valueI * coefficient_[j];
165                         }
166                    }
167               }
168               functionValue_ -= offset_;
169          } else {
170               abort();
171               // do scaling
172               const double * columnScale = model->columnScale();
173               for (int i = 0; i < numberCoefficients_; i++) {
174                    int iColumn = column_[i];
175                    double value = solution[iColumn]; // already scaled
176                    double coefficient = coefficient_[i] * columnScale[iColumn];
177                    functionValue_ += value * coefficient;
178                    lastGradient_[iColumn] = coefficient;
179               }
180          }
181     }
182     functionValue = functionValue_;
183     offset = offset_;
184     CoinMemcpyN(lastGradient_, numberColumns_, gradient);
185     return 0;
186}
187// Resize constraint
188void
189ClpConstraintQuadratic::resize(int newNumberColumns)
190{
191     if (numberColumns_ != newNumberColumns) {
192          abort();
193#ifndef NDEBUG
194          int lastColumn = column_[numberCoefficients_-1];
195#endif
196          assert (newNumberColumns > lastColumn);
197          delete [] lastGradient_;
198          lastGradient_ = NULL;
199          numberColumns_ = newNumberColumns;
200     }
201}
202// Delete columns in  constraint
203void
204ClpConstraintQuadratic::deleteSome(int numberToDelete, const int * which)
205{
206     if (numberToDelete) {
207          abort();
208          int i ;
209          char * deleted = new char[numberColumns_];
210          memset(deleted, 0, numberColumns_ * sizeof(char));
211          for (i = 0; i < numberToDelete; i++) {
212               int j = which[i];
213               if (j >= 0 && j < numberColumns_ && !deleted[j]) {
214                    deleted[j] = 1;
215               }
216          }
217          int n = 0;
218          for (i = 0; i < numberCoefficients_; i++) {
219               int iColumn = column_[i];
220               if (!deleted[iColumn]) {
221                    column_[n] = iColumn;
222                    coefficient_[n++] = coefficient_[i];
223               }
224          }
225          numberCoefficients_ = n;
226     }
227}
228// Scale constraint
229void
230ClpConstraintQuadratic::reallyScale(const double * )
231{
232     abort();
233}
234/* Given a zeroed array sets nonquadratic columns to 1.
235   Returns number of nonlinear columns
236*/
237int
238ClpConstraintQuadratic::markNonlinear(char * which) const
239{
240     int iColumn;
241     for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
242          CoinBigIndex j;
243          for (j = start_[iColumn]; j < start_[iColumn+1]; j++) {
244               int jColumn = column_[j];
245               if (jColumn >= 0) {
246                    assert (jColumn < numberQuadraticColumns_);
247                    which[jColumn] = 1;
248                    which[iColumn] = 1;
249               }
250          }
251     }
252     int numberCoefficients = 0;
253     for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
254          if (which[iColumn])
255               numberCoefficients++;
256     }
257     return numberCoefficients;
258}
259/* Given a zeroed array sets possible nonzero coefficients to 1.
260   Returns number of nonzeros
261*/
262int
263ClpConstraintQuadratic::markNonzero(char * which) const
264{
265     int iColumn;
266     for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
267          CoinBigIndex j;
268          for (j = start_[iColumn]; j < start_[iColumn+1]; j++) {
269               int jColumn = column_[j];
270               if (jColumn >= 0) {
271                    assert (jColumn < numberQuadraticColumns_);
272                    which[jColumn] = 1;
273               }
274               which[iColumn] = 1;
275          }
276     }
277     int numberCoefficients = 0;
278     for (iColumn = 0; iColumn < numberQuadraticColumns_; iColumn++) {
279          if (which[iColumn])
280               numberCoefficients++;
281     }
282     return numberCoefficients;
283}
284// Number of coefficients
285int
286ClpConstraintQuadratic::numberCoefficients() const
287{
288     return numberCoefficients_;
289}
Note: See TracBrowser for help on using the repository browser.