ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/group/trunk/OOPSE-4/src/selection/SelectionEvaluator.cpp
Revision: 3533
Committed: Wed Oct 7 20:49:50 2009 UTC (14 years, 10 months ago) by cli2
File size: 13843 byte(s)
Log Message:
Bug fixes in the SelectionEvaluator and SelectionCompiler
Added print option in Restraints

File Contents

# Content
1 /*
2 * Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3 *
4 * The University of Notre Dame grants you ("Licensee") a
5 * non-exclusive, royalty free, license to use, modify and
6 * redistribute this software in source and binary code form, provided
7 * that the following conditions are met:
8 *
9 * 1. Acknowledgement of the program authors must be made in any
10 * publication of scientific results based in part on use of the
11 * program. An acceptable form of acknowledgement is citation of
12 * the article in which the program was described (Matthew
13 * A. Meineke, Charles F. Vardeman II, Teng Lin, Christopher
14 * J. Fennell and J. Daniel Gezelter, "OOPSE: An Object-Oriented
15 * Parallel Simulation Engine for Molecular Dynamics,"
16 * J. Comput. Chem. 26, pp. 252-271 (2005))
17 *
18 * 2. Redistributions of source code must retain the above copyright
19 * notice, this list of conditions and the following disclaimer.
20 *
21 * 3. Redistributions in binary form must reproduce the above copyright
22 * notice, this list of conditions and the following disclaimer in the
23 * documentation and/or other materials provided with the
24 * distribution.
25 *
26 * This software is provided "AS IS," without a warranty of any
27 * kind. All express or implied conditions, representations and
28 * warranties, including any implied warranty of merchantability,
29 * fitness for a particular purpose or non-infringement, are hereby
30 * excluded. The University of Notre Dame and its licensors shall not
31 * be liable for any damages suffered by licensee as a result of
32 * using, modifying or distributing the software or its
33 * derivatives. In no event will the University of Notre Dame or its
34 * licensors be liable for any lost revenue, profit or data, or for
35 * direct, indirect, special, consequential, incidental or punitive
36 * damages, however caused and regardless of the theory of liability,
37 * arising out of the use of or inability to use software, even if the
38 * University of Notre Dame has been advised of the possibility of
39 * such damages.
40 */
41
42 #include <stack>
43 #include "selection/SelectionEvaluator.hpp"
44 #include "primitives/Atom.hpp"
45 #include "primitives/DirectionalAtom.hpp"
46 #include "primitives/RigidBody.hpp"
47 #include "primitives/Molecule.hpp"
48 #include "io/basic_ifstrstream.hpp"
49
50 namespace oopse {
51
52
53 SelectionEvaluator::SelectionEvaluator(SimInfo* si)
54 : info(si), nameFinder(info), distanceFinder(info), indexFinder(info),
55 isLoaded_(false){
56 nStuntDouble = info->getNGlobalAtoms() + info->getNGlobalRigidBodies();
57 }
58
59 bool SelectionEvaluator::loadScript(const std::string& filename,
60 const std::string& script) {
61 clearDefinitionsAndLoadPredefined();
62 this->filename = filename;
63 this->script = script;
64 if (! compiler.compile(filename, script)) {
65 error = true;
66 errorMessage = compiler.getErrorMessage();
67
68 sprintf( painCave.errMsg,
69 "SelectionCompiler Error: %s\n", errorMessage.c_str());
70 painCave.severity = OOPSE_ERROR;
71 painCave.isFatal = 1;
72 simError();
73 return false;
74 }
75
76 pc = 0;
77 aatoken = compiler.getAatokenCompiled();
78 linenumbers = compiler.getLineNumbers();
79 lineIndices = compiler.getLineIndices();
80
81 std::vector<std::vector<Token> >::const_iterator i;
82
83 isDynamic_ = false;
84 for (i = aatoken.begin(); i != aatoken.end(); ++i) {
85 if (containDynamicToken(*i)) {
86 isDynamic_ = true;
87 break;
88 }
89 }
90
91 isLoaded_ = true;
92 return true;
93 }
94
95 void SelectionEvaluator::clearState() {
96 error = false;
97 errorMessage = "";
98 }
99
100 bool SelectionEvaluator::loadScriptString(const std::string& script) {
101 clearState();
102 return loadScript("", script);
103 }
104
105 bool SelectionEvaluator::loadScriptFile(const std::string& filename) {
106 clearState();
107 return loadScriptFileInternal(filename);
108 }
109
110 bool SelectionEvaluator::loadScriptFileInternal(const std::string & filename) {
111 ifstrstream ifs(filename.c_str());
112 if (!ifs.is_open()) {
113 return false;
114 }
115
116 const int bufferSize = 65535;
117 char buffer[bufferSize];
118 std::string script;
119 while(ifs.getline(buffer, bufferSize)) {
120 script += buffer;
121 }
122 return loadScript(filename, script);
123 }
124
125 void SelectionEvaluator::instructionDispatchLoop(OOPSEBitSet& bs){
126
127 while ( pc < aatoken.size()) {
128 statement = aatoken[pc++];
129 statementLength = statement.size();
130 Token token = statement[0];
131 switch (token.tok) {
132 case Token::define:
133 define();
134 break;
135 case Token::select:
136 select(bs);
137 break;
138 default:
139 unrecognizedCommand(token);
140 return;
141 }
142 }
143
144 }
145
146 OOPSEBitSet SelectionEvaluator::expression(const std::vector<Token>& code,
147 int pcStart) {
148 OOPSEBitSet bs;
149 std::stack<OOPSEBitSet> stack;
150
151 for (int pc = pcStart; pc < code.size(); ++pc) {
152 Token instruction = code[pc];
153
154 switch (instruction.tok) {
155 case Token::expressionBegin:
156 break;
157 case Token::expressionEnd:
158 break;
159 case Token::all:
160 bs = OOPSEBitSet(nStuntDouble);
161 bs.setAll();
162 stack.push(bs);
163 break;
164 case Token::none:
165 bs = OOPSEBitSet(nStuntDouble);
166 stack.push(bs);
167 break;
168 case Token::opOr:
169 bs = stack.top();
170 stack.pop();
171 stack.top() |= bs;
172 break;
173 case Token::opAnd:
174 bs = stack.top();
175 stack.pop();
176 stack.top() &= bs;
177 break;
178 case Token::opNot:
179 stack.top().flip();
180 break;
181 case Token::within:
182 withinInstruction(instruction, stack.top());
183 break;
184 //case Token::selected:
185 // stack.push(getSelectionSet());
186 // break;
187 case Token::name:
188 stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value)));
189 break;
190 case Token::index:
191 stack.push(indexInstruction(instruction.value));
192 break;
193 case Token::identifier:
194 stack.push(lookupValue(boost::any_cast<std::string>(instruction.value)));
195 break;
196 case Token::opLT:
197 case Token::opLE:
198 case Token::opGE:
199 case Token::opGT:
200 case Token::opEQ:
201 case Token::opNE:
202 stack.push(comparatorInstruction(instruction));
203 break;
204 default:
205 unrecognizedExpression();
206 }
207 }
208 if (stack.size() != 1)
209 evalError("atom expression compiler error - stack over/underflow");
210
211 return stack.top();
212 }
213
214
215
216 OOPSEBitSet SelectionEvaluator::comparatorInstruction(const Token& instruction) {
217 int comparator = instruction.tok;
218 int property = instruction.intValue;
219 float comparisonValue = boost::any_cast<float>(instruction.value);
220 float propertyValue;
221 OOPSEBitSet bs(nStuntDouble);
222 bs.clearAll();
223
224 SimInfo::MoleculeIterator mi;
225 Molecule* mol;
226 Molecule::AtomIterator ai;
227 Atom* atom;
228 Molecule::RigidBodyIterator rbIter;
229 RigidBody* rb;
230
231 for (mol = info->beginMolecule(mi); mol != NULL;
232 mol = info->nextMolecule(mi)) {
233
234 for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
235 compareProperty(atom, bs, property, comparator, comparisonValue);
236 }
237
238 for (rb = mol->beginRigidBody(rbIter); rb != NULL;
239 rb = mol->nextRigidBody(rbIter)) {
240 compareProperty(rb, bs, property, comparator, comparisonValue);
241 }
242 }
243
244 return bs;
245 }
246
247 void SelectionEvaluator::compareProperty(StuntDouble* sd, OOPSEBitSet& bs,
248 int property, int comparator,
249 float comparisonValue) {
250 RealType propertyValue = 0.0;
251 switch (property) {
252 case Token::mass:
253 propertyValue = sd->getMass();
254 break;
255 case Token::charge:
256 if (sd->isAtom()){
257 Atom* atom = static_cast<Atom*>(sd);
258 propertyValue = getCharge(atom);
259 } else if (sd->isRigidBody()) {
260 RigidBody* rb = static_cast<RigidBody*>(sd);
261 RigidBody::AtomIterator ai;
262 Atom* atom;
263 for (atom = rb->beginAtom(ai); atom != NULL; atom = rb->nextAtom(ai)) {
264 propertyValue+= getCharge(atom);
265 }
266 }
267 break;
268 case Token::x:
269 propertyValue = sd->getPos().x();
270 break;
271 case Token::y:
272 propertyValue = sd->getPos().y();
273 break;
274 case Token::z:
275 propertyValue = sd->getPos().z();
276 break;
277 default:
278 unrecognizedAtomProperty(property);
279 }
280
281 bool match = false;
282 switch (comparator) {
283 case Token::opLT:
284 match = propertyValue < comparisonValue;
285 break;
286 case Token::opLE:
287 match = propertyValue <= comparisonValue;
288 break;
289 case Token::opGE:
290 match = propertyValue >= comparisonValue;
291 break;
292 case Token::opGT:
293 match = propertyValue > comparisonValue;
294 break;
295 case Token::opEQ:
296 match = propertyValue == comparisonValue;
297 break;
298 case Token::opNE:
299 match = propertyValue != comparisonValue;
300 break;
301 }
302 if (match)
303 bs.setBitOn(sd->getGlobalIndex());
304
305 }
306
307 void SelectionEvaluator::withinInstruction(const Token& instruction,
308 OOPSEBitSet& bs){
309
310 boost::any withinSpec = instruction.value;
311 float distance;
312 if (withinSpec.type() == typeid(float)){
313 distance = boost::any_cast<float>(withinSpec);
314 } else if (withinSpec.type() == typeid(int)) {
315 distance = boost::any_cast<int>(withinSpec);
316 } else {
317 evalError("casting error in withinInstruction");
318 bs.clearAll();
319 }
320
321 bs = distanceFinder.find(bs, distance);
322 }
323
324 void SelectionEvaluator::define() {
325 assert(statement.size() >= 3);
326
327 std::string variable = boost::any_cast<std::string>(statement[1].value);
328
329 variables.insert(VariablesType::value_type(variable,
330 expression(statement, 2)));
331 }
332
333
334 /** @todo */
335 void SelectionEvaluator::predefine(const std::string& script) {
336
337 if (compiler.compile("#predefine", script)) {
338 std::vector<std::vector<Token> > aatoken = compiler.getAatokenCompiled();
339 if (aatoken.size() != 1) {
340 evalError("predefinition does not have exactly 1 command:"
341 + script);
342 return;
343 }
344 std::vector<Token> statement = aatoken[0];
345 if (statement.size() > 2) {
346 int tok = statement[1].tok;
347 if (tok == Token::identifier ||
348 (tok & Token::predefinedset) == Token::predefinedset) {
349 std::string variable = boost::any_cast<std::string>(statement[1].value);
350 variables.insert(VariablesType::value_type(variable, statement));
351
352 } else {
353 evalError("invalid variable name:" + script);
354 }
355 }else {
356 evalError("bad predefinition length:" + script);
357 }
358
359 } else {
360 evalError("predefined set compile error:" + script +
361 "\ncompile error:" + compiler.getErrorMessage());
362 }
363 }
364
365 void SelectionEvaluator::select(OOPSEBitSet& bs){
366 bs = expression(statement, 1);
367 }
368
369 OOPSEBitSet SelectionEvaluator::lookupValue(const std::string& variable){
370
371 OOPSEBitSet bs(nStuntDouble);
372 std::map<std::string, boost::any>::iterator i = variables.find(variable);
373
374 if (i != variables.end()) {
375 if (i->second.type() == typeid(OOPSEBitSet)) {
376 return boost::any_cast<OOPSEBitSet>(i->second);
377 } else if (i->second.type() == typeid(std::vector<Token>)){
378 bs = expression(boost::any_cast<std::vector<Token> >(i->second), 2);
379 i->second = bs; /**@todo fixme */
380 return bs;
381 }
382 } else {
383 unrecognizedIdentifier(variable);
384 }
385
386 return bs;
387 }
388
389 OOPSEBitSet SelectionEvaluator::nameInstruction(const std::string& name){
390 return nameFinder.match(name);
391 }
392
393 bool SelectionEvaluator::containDynamicToken(const std::vector<Token>& tokens){
394 std::vector<Token>::const_iterator i;
395 for (i = tokens.begin(); i != tokens.end(); ++i) {
396 if (i->tok & Token::dynamic) {
397 return true;
398 }
399 }
400
401 return false;
402 }
403
404 void SelectionEvaluator::clearDefinitionsAndLoadPredefined() {
405 variables.clear();
406 //load predefine
407 //predefine();
408 }
409
410 OOPSEBitSet SelectionEvaluator::evaluate() {
411 OOPSEBitSet bs(nStuntDouble);
412 if (isLoaded_) {
413 pc = 0;
414 instructionDispatchLoop(bs);
415 }
416
417 return bs;
418 }
419
420 OOPSEBitSet SelectionEvaluator::indexInstruction(const boost::any& value) {
421 OOPSEBitSet bs(nStuntDouble);
422
423 if (value.type() == typeid(int)) {
424 int index = boost::any_cast<int>(value);
425 if (index < 0 || index >= bs.size()) {
426 invalidIndex(index);
427 } else {
428 bs = indexFinder.find(index);
429 }
430 } else if (value.type() == typeid(std::pair<int, int>)) {
431 std::pair<int, int> indexRange= boost::any_cast<std::pair<int, int> >(value);
432 assert(indexRange.first <= indexRange.second);
433 if (indexRange.first < 0 || indexRange.second >= bs.size()) {
434 invalidIndexRange(indexRange);
435 }else {
436 bs = indexFinder.find(indexRange.first, indexRange.second);
437 }
438 }
439
440 return bs;
441 }
442
443
444 RealType SelectionEvaluator::getCharge(Atom* atom) {
445 RealType charge =0.0;
446 AtomType* atomType = atom->getAtomType();
447 if (atomType->isCharge()) {
448 GenericData* data = atomType->getPropertyByName("Charge");
449 if (data != NULL) {
450 DoubleGenericData* doubleData= dynamic_cast<DoubleGenericData*>(data);
451
452 if (doubleData != NULL) {
453 charge = doubleData->getData();
454
455 } else {
456 sprintf( painCave.errMsg,
457 "Can not cast GenericData to DoubleGenericData\n");
458 painCave.severity = OOPSE_ERROR;
459 painCave.isFatal = 1;
460 simError();
461 }
462 }
463 }
464
465 return charge;
466 }
467
468 }