1/* $Id: ClpConstraintQuadratic.cpp 1665 2011-01-04 17:55:54Z lou $ */
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}
290