source: branches/BSP/trunk/Clp/src/ClpConstraintQuadratic.cpp @ 1069

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

add defaulthandler and stuff for nonlinear

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