OOPSE/libmdtools/Constraint.cpp

#include "Constraint.hpp"
#include "SymMatrix.hpp"

ConstraintBase::ConstraintBase(){
  init_ndim = false;
  ndim = 0;
}

ConstraintBase::ConstraintBase(int dim){
  ndim = dim;
  init_ndim = true;
}

void ConstraintBase::setDim(int dim){

  if(isDimSet() && dim != ndim){
    cout << "ConstraintBase Warning : About to change ndim which is already set" << endl;
  }

  ndim = dim;
  init_dim = true;
}

SimpleBoundCons::SimpleBoundCons(int index, double bound, bool flag)
                        : ConstraintBase(){

  bound = b;

  consType = ConsType::simpleBound;

  if (flag){
    boundType = BoundType::upper;
  }
  else{
    boundType = BoundType::lower;
  }
}

SimpleBoundCons::SimpleBoundCons(int ndim, int index, double bound, bool flag)
                        : ConstraintBase(ndim){

  bound = b;

  consType = ConsType::simpleBound;

  if (flag){
    boundType = BoundType::upper;
  }
  else{
    boundType = BoundType::lower;
  }
}
double SimpleBoundCons::calcResidual(vector<double>& x){
  double residual;

  residual = x[index] - bound;

  if(boundType == BoundType::lower)
    return -residual;
  else if (boundType == BoundType::upper)
    return residual;
   else{
    cout << "SimpleBoundCons Error: BoundType of SimpleBoundCons can not be BoundType::equ" << endl;
    exit(ERROR_CONSTRAINT);
   }
}

vector<double> SimpleBoundCons::calcConsGrad(vector<double>& x){
  vector<double> result(ndim);
  result  = 0;

  if(boundType == BoundType::lower)
    result[index] = -1.0;
  else if (boundType == BoundType::upper)
    result[index] = 1.0;
   else{
    cout << "SimpleBoundCons Error: BoundType of SimpleBoundCons can not be BoundType::equ" << endl;
    exit(ERROR_CONSTRAINT);
   }
  
  result result;
}
SymMatrix SimpleBoundCons::calcConsHessian(vector<double>& x){
    SymMatrix H(ndim);
    H = 0;
    return H;
}

LinearCons::LinearCons(vector<int>&  theIndex, vector<double>&  theCoeff, double b, BoundType bType)
              :ConstraintBase(){
              
  bound = b;
  boundType = bType;
  if(bType == BoundType::equ){
    consType = ConsType::linearEqu;
  }
  else{
    consType = ConsType::linearInequ;
  }

  index = theIndex;
  coeff = theCoeff;
}

LinearCons::LinearCons(int dim, vector<int>&  theIndex, vector<double>&  theCoeff, double b, BoundType bType)
              :ConstraintBase(dim) {

  if (dim != theCoeff.size()){
    cout << "LinearCons Error: the dimension of index and coeff does not match" << endl;
    exit(ERROR_CONSTRAINT);
  }
  
  bound = b;
  boundType = bType;
  if(bType == BoundType::equ){
    consType = ConsType::linearEqu;
  }
  else{
    consType = ConsType::linearInequ;
  }

  index = theIndex;
  coeff = theCoeff;
}

double LinearCons::calcResidual(vector<double>& x){
  double residue;
  double valueOfLinearCons;

  valueOfLinearCons = 0;
  
  for (int i = 0; i < coeff.siz(); i++)
    valueOfLinearCons += coeff[i] * x[index[i]];
  
  residue = valueOfLinearCons - bound;

  if(boundType == BoundType::lower)
    residue = -residue;
  
  return residue;  
}

vector<double> LinearCons::calcConsGrad(vector<double>& x){
  vector<double> consGrad(ndim);
  double sign;

  consGrad = 0;

  if(boundType == BoundType::lower)
    sign = -1.0;
  else
    sign = 1.0;
  
  for(int i = 0; i < coeff.size(); i++)      
    result[index[i]] = coeff[i] * sign;
  return consGrad;
}

SymMatrix LinearCons::calcConsHessian(vector<double>& x){
    SymMatrix H(ndim);
    H = 0;
    return H;  
}

NonlinearCons::NonlinearCons(vector<int>& theIndex, NLModel* theModel , double b, BoundType bType)
                  : ConstraintBase(){

  if(theIndex.size() != theModel->getDim()){
    cout << "NonlinearCons Error: the dimension of index and the model does not match" << endl;
    exit(ERROR_CONSTRAINT);   
  }
  
  bound = b;
  boundType = bType;
  
  if(bType == BoundType::equ){
    consType = ConsType::nonlinearEqu;
  }
  else{
    consType = ConsType::nonlinearInequ;
  }

  index = theIndex;
  model = theModel;
}

NonlinearCons::NonlinearCons(int dim, vector<int>& theIndex, NLModel* theModel , double b, BoundType bType)
                   :ConstraintBase(dim){

  if(theIndex.size() != theModel->getDim()){
    cout << "NonlinearCons Error: the dimension of index and the model does not match" << endl;
    exit(ERROR_CONSTRAINT);   
  }
  
  bound = b;
  boundType = bType;
  
  if(bType == BoundType::equ){
    consType = ConsType::nonlinearEqu;
  }
  else{
    consType = ConsType::nonlinearInequ;
  }

  index = theIndex;
  model = theModel;
}


void NonlinearCons::setDim(int dim){
  
  if(theIndex.size() != theModel->getDim()){
    cout << "NonlinearCons Error: the dimension of index and the model does not match" << endl;
    exit(ERROR_CONSTRAINT);   
  }

  ConstraintBase::setDm(dim);
  
}

double NonlinearCons::::calcResidual(vector<double>& x){
  double fVal;

  fVal = model->calcF(x);

  fVal -= bound;

  if(boundType == BoundType::lower)
    fVal = -fVal;
  
  return fVal;  
}

vector<double> NonlinearCons::calcConsGrad(vector<double>& x){
  vector<double> consGrad(ndim);

  consGrad = model->calcGrad(x);

  if(boundType == BoundType::lower)
    consGrad = -consGrad;

  return consGrad;
}

SymMatrix NonlinearCons::calcConsHessian(vector<double>& x){
  SymMatrix H(ndim);
  
  H = model->calcHessian(x);

  if(boundType == BoundType::lower)
    H = -H;

  return H;  
}

ConstraintList::ConstraintList(){
  consType = 0;
}

ConstraintList::~ConstraintList(){
  
  for(int i = 0; i < constraint.size(); i++)
    if(!constraints[i])
      delete constraints[i];

  constraints.clear();

}
void ConstraintList::addConstraint(ConstraintBase* cons){
  constraints.push_back(cons);
  consType |= cons->getConsType();
}
Revision:	987
Committed:	Tue Jan 27 19:15:20 2004 UTC (21 years, 9 months ago) by tim
File size:	5983 byte(s)
Log Message:	revision of constraint for Nonlinear Optimization Model
#	User	Rev	Content
1	tim	987	#include "Constraint.hpp"
2			#include "SymMatrix.hpp"
3
4			ConstraintBase::ConstraintBase(){
5			init_ndim = false;
6			ndim = 0;
7			}
8
9			ConstraintBase::ConstraintBase(int dim){
10			ndim = dim;
11			init_ndim = true;
12			}
13
14			void ConstraintBase::setDim(int dim){
15
16			if(isDimSet() && dim != ndim){
17			cout << "ConstraintBase Warning : About to change ndim which is already set" << endl;
18			}
19
20			ndim = dim;
21			init_dim = true;
22			}
23
24			SimpleBoundCons::SimpleBoundCons(int index, double bound, bool flag)
25			: ConstraintBase(){
26
27			bound = b;
28
29			consType = ConsType::simpleBound;
30
31			if (flag){
32			boundType = BoundType::upper;
33			}
34			else{
35			boundType = BoundType::lower;
36			}
37			}
38
39			SimpleBoundCons::SimpleBoundCons(int ndim, int index, double bound, bool flag)
40			: ConstraintBase(ndim){
41
42			bound = b;
43
44			consType = ConsType::simpleBound;
45
46			if (flag){
47			boundType = BoundType::upper;
48			}
49			else{
50			boundType = BoundType::lower;
51			}
52			}
53			double SimpleBoundCons::calcResidual(vector<double>& x){
54			double residual;
55
56			residual = x[index] - bound;
57
58			if(boundType == BoundType::lower)
59			return -residual;
60			else if (boundType == BoundType::upper)
61			return residual;
62			else{
63			cout << "SimpleBoundCons Error: BoundType of SimpleBoundCons can not be BoundType::equ" << endl;
64			exit(ERROR_CONSTRAINT);
65			}
66			}
67
68			vector<double> SimpleBoundCons::calcConsGrad(vector<double>& x){
69			vector<double> result(ndim);
70			result = 0;
71
72			if(boundType == BoundType::lower)
73			result[index] = -1.0;
74			else if (boundType == BoundType::upper)
75			result[index] = 1.0;
76			else{
77			cout << "SimpleBoundCons Error: BoundType of SimpleBoundCons can not be BoundType::equ" << endl;
78			exit(ERROR_CONSTRAINT);
79			}
80
81			result result;
82			}
83			SymMatrix SimpleBoundCons::calcConsHessian(vector<double>& x){
84			SymMatrix H(ndim);
85			H = 0;
86			return H;
87			}
88
89			LinearCons::LinearCons(vector<int>& theIndex, vector<double>& theCoeff, double b, BoundType bType)
90			:ConstraintBase(){
91
92			bound = b;
93			boundType = bType;
94			if(bType == BoundType::equ){
95			consType = ConsType::linearEqu;
96			}
97			else{
98			consType = ConsType::linearInequ;
99			}
100
101			index = theIndex;
102			coeff = theCoeff;
103			}
104
105			LinearCons::LinearCons(int dim, vector<int>& theIndex, vector<double>& theCoeff, double b, BoundType bType)
106			:ConstraintBase(dim) {
107
108			if (dim != theCoeff.size()){
109			cout << "LinearCons Error: the dimension of index and coeff does not match" << endl;
110			exit(ERROR_CONSTRAINT);
111			}
112
113			bound = b;
114			boundType = bType;
115			if(bType == BoundType::equ){
116			consType = ConsType::linearEqu;
117			}
118			else{
119			consType = ConsType::linearInequ;
120			}
121
122			index = theIndex;
123			coeff = theCoeff;
124			}
125
126			double LinearCons::calcResidual(vector<double>& x){
127			double residue;
128			double valueOfLinearCons;
129
130			valueOfLinearCons = 0;
131
132			for (int i = 0; i < coeff.siz(); i++)
133			valueOfLinearCons += coeff[i] * x[index[i]];
134
135			residue = valueOfLinearCons - bound;
136
137			if(boundType == BoundType::lower)
138			residue = -residue;
139
140			return residue;
141			}
142
143			vector<double> LinearCons::calcConsGrad(vector<double>& x){
144			vector<double> consGrad(ndim);
145			double sign;
146
147			consGrad = 0;
148
149			if(boundType == BoundType::lower)
150			sign = -1.0;
151			else
152			sign = 1.0;
153
154			for(int i = 0; i < coeff.size(); i++)
155			result[index[i]] = coeff[i] * sign;
156			return consGrad;
157			}
158
159			SymMatrix LinearCons::calcConsHessian(vector<double>& x){
160			SymMatrix H(ndim);
161			H = 0;
162			return H;
163			}
164
165			NonlinearCons::NonlinearCons(vector<int>& theIndex, NLModel* theModel , double b, BoundType bType)
166			: ConstraintBase(){
167
168			if(theIndex.size() != theModel->getDim()){
169			cout << "NonlinearCons Error: the dimension of index and the model does not match" << endl;
170			exit(ERROR_CONSTRAINT);
171			}
172
173			bound = b;
174			boundType = bType;
175
176			if(bType == BoundType::equ){
177			consType = ConsType::nonlinearEqu;
178			}
179			else{
180			consType = ConsType::nonlinearInequ;
181			}
182
183			index = theIndex;
184			model = theModel;
185			}
186
187			NonlinearCons::NonlinearCons(int dim, vector<int>& theIndex, NLModel* theModel , double b, BoundType bType)
188			:ConstraintBase(dim){
189
190			if(theIndex.size() != theModel->getDim()){
191			cout << "NonlinearCons Error: the dimension of index and the model does not match" << endl;
192			exit(ERROR_CONSTRAINT);
193			}
194
195			bound = b;
196			boundType = bType;
197
198			if(bType == BoundType::equ){
199			consType = ConsType::nonlinearEqu;
200			}
201			else{
202			consType = ConsType::nonlinearInequ;
203			}
204
205			index = theIndex;
206			model = theModel;
207			}
208
209
210			void NonlinearCons::setDim(int dim){
211
212			if(theIndex.size() != theModel->getDim()){
213			cout << "NonlinearCons Error: the dimension of index and the model does not match" << endl;
214			exit(ERROR_CONSTRAINT);
215			}
216
217			ConstraintBase::setDm(dim);
218
219			}
220
221			double NonlinearCons::::calcResidual(vector<double>& x){
222			double fVal;
223
224			fVal = model->calcF(x);
225
226			fVal -= bound;
227
228			if(boundType == BoundType::lower)
229			fVal = -fVal;
230
231			return fVal;
232			}
233
234			vector<double> NonlinearCons::calcConsGrad(vector<double>& x){
235			vector<double> consGrad(ndim);
236
237			consGrad = model->calcGrad(x);
238
239			if(boundType == BoundType::lower)
240			consGrad = -consGrad;
241
242			return consGrad;
243			}
244
245			SymMatrix NonlinearCons::calcConsHessian(vector<double>& x){
246			SymMatrix H(ndim);
247
248			H = model->calcHessian(x);
249
250			if(boundType == BoundType::lower)
251			H = -H;
252
253			return H;
254			}
255
256			ConstraintList::ConstraintList(){
257			consType = 0;
258			}
259
260			ConstraintList::~ConstraintList(){
261
262			for(int i = 0; i < constraint.size(); i++)
263			if(!constraints[i])
264			delete constraints[i];
265
266			constraints.clear();
267
268			}
269			void ConstraintList::addConstraint(ConstraintBase* cons){
270			constraints.push_back(cons);
271			consType \|= cons->getConsType();
272			}