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 (20 years, 5 months ago) by tim
File size:	5983 byte(s)
Log Message:	revision of constraint for Nonlinear Optimization Model
#	Content
1	#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	}