ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/selection/SelectionEvaluator.cpp
Revision: 1364
Committed: Wed Oct 7 20:49:50 2009 UTC (15 years, 6 months ago) by cli2
Original Path: trunk/src/selection/SelectionEvaluator.cpp
File size: 13843 byte(s)
Log Message:
Bug fixes in the SelectionEvaluator and SelectionCompiler
Added print option in Restraints

File Contents

# User Rev Content
1 tim 277 /*
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 tim 283 #include <stack>
43 tim 277 #include "selection/SelectionEvaluator.hpp"
44 tim 283 #include "primitives/Atom.hpp"
45     #include "primitives/DirectionalAtom.hpp"
46     #include "primitives/RigidBody.hpp"
47     #include "primitives/Molecule.hpp"
48 cli2 1364 #include "io/basic_ifstrstream.hpp"
49 tim 283
50 tim 277 namespace oopse {
51    
52    
53 gezelter 507 SelectionEvaluator::SelectionEvaluator(SimInfo* si)
54 cli2 1364 : info(si), nameFinder(info), distanceFinder(info), indexFinder(info),
55     isLoaded_(false){
56     nStuntDouble = info->getNGlobalAtoms() + info->getNGlobalRigidBodies();
57 gezelter 507 }
58 tim 283
59 cli2 1364 bool SelectionEvaluator::loadScript(const std::string& filename,
60     const std::string& script) {
61 tim 288 clearDefinitionsAndLoadPredefined();
62 tim 278 this->filename = filename;
63     this->script = script;
64     if (! compiler.compile(filename, script)) {
65 gezelter 507 error = true;
66     errorMessage = compiler.getErrorMessage();
67 cli2 1364
68     sprintf( painCave.errMsg,
69     "SelectionCompiler Error: %s\n", errorMessage.c_str());
70     painCave.severity = OOPSE_ERROR;
71     painCave.isFatal = 1;
72     simError();
73 gezelter 507 return false;
74 tim 278 }
75    
76     pc = 0;
77     aatoken = compiler.getAatokenCompiled();
78     linenumbers = compiler.getLineNumbers();
79     lineIndices = compiler.getLineIndices();
80 tim 288
81     std::vector<std::vector<Token> >::const_iterator i;
82    
83     isDynamic_ = false;
84     for (i = aatoken.begin(); i != aatoken.end(); ++i) {
85 gezelter 507 if (containDynamicToken(*i)) {
86     isDynamic_ = true;
87     break;
88     }
89 tim 288 }
90    
91     isLoaded_ = true;
92 tim 278 return true;
93 gezelter 507 }
94 tim 278
95 gezelter 507 void SelectionEvaluator::clearState() {
96 tim 278 error = false;
97 tim 281 errorMessage = "";
98 gezelter 507 }
99 tim 278
100 gezelter 507 bool SelectionEvaluator::loadScriptString(const std::string& script) {
101 tim 278 clearState();
102 tim 281 return loadScript("", script);
103 gezelter 507 }
104 tim 278
105 gezelter 507 bool SelectionEvaluator::loadScriptFile(const std::string& filename) {
106 tim 278 clearState();
107     return loadScriptFileInternal(filename);
108 gezelter 507 }
109 tim 278
110 cli2 1364 bool SelectionEvaluator::loadScriptFileInternal(const std::string & filename) {
111     ifstrstream ifs(filename.c_str());
112 tim 288 if (!ifs.is_open()) {
113 gezelter 507 return false;
114 tim 288 }
115 cli2 1364
116 tim 288 const int bufferSize = 65535;
117     char buffer[bufferSize];
118     std::string script;
119     while(ifs.getline(buffer, bufferSize)) {
120 gezelter 507 script += buffer;
121 tim 288 }
122     return loadScript(filename, script);
123 gezelter 507 }
124 cli2 1364
125 tim 770 void SelectionEvaluator::instructionDispatchLoop(OOPSEBitSet& bs){
126 tim 288
127 tim 281 while ( pc < aatoken.size()) {
128 gezelter 507 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 tim 278 }
143 tim 288
144 gezelter 507 }
145 tim 278
146 cli2 1364 OOPSEBitSet SelectionEvaluator::expression(const std::vector<Token>& code,
147     int pcStart) {
148 tim 770 OOPSEBitSet bs;
149 cli2 1364 std::stack<OOPSEBitSet> stack;
150    
151 tim 288 for (int pc = pcStart; pc < code.size(); ++pc) {
152 tim 278 Token instruction = code[pc];
153 tim 282
154 tim 278 switch (instruction.tok) {
155 tim 281 case Token::expressionBegin:
156 tim 278 break;
157 tim 281 case Token::expressionEnd:
158 tim 282 break;
159 tim 281 case Token::all:
160 tim 770 bs = OOPSEBitSet(nStuntDouble);
161 tim 283 bs.setAll();
162     stack.push(bs);
163 tim 278 break;
164 tim 281 case Token::none:
165 tim 770 bs = OOPSEBitSet(nStuntDouble);
166 tim 283 stack.push(bs);
167 tim 278 break;
168 tim 281 case Token::opOr:
169 tim 283 bs = stack.top();
170     stack.pop();
171     stack.top() |= bs;
172 tim 278 break;
173 tim 281 case Token::opAnd:
174 tim 283 bs = stack.top();
175     stack.pop();
176     stack.top() &= bs;
177 tim 278 break;
178 tim 281 case Token::opNot:
179 tim 283 stack.top().flip();
180 tim 278 break;
181 tim 281 case Token::within:
182 tim 283 withinInstruction(instruction, stack.top());
183 tim 278 break;
184 gezelter 507 //case Token::selected:
185     // stack.push(getSelectionSet());
186     // break;
187 tim 281 case Token::name:
188 tim 283 stack.push(nameInstruction(boost::any_cast<std::string>(instruction.value)));
189 tim 278 break;
190 tim 295 case Token::index:
191     stack.push(indexInstruction(instruction.value));
192 tim 281 break;
193     case Token::identifier:
194 tim 283 stack.push(lookupValue(boost::any_cast<std::string>(instruction.value)));
195 tim 281 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 tim 283 stack.push(comparatorInstruction(instruction));
203 tim 278 break;
204     default:
205     unrecognizedExpression();
206     }
207     }
208 tim 283 if (stack.size() != 1)
209 tim 278 evalError("atom expression compiler error - stack over/underflow");
210 tim 283
211     return stack.top();
212 tim 278 }
213    
214    
215    
216 tim 770 OOPSEBitSet SelectionEvaluator::comparatorInstruction(const Token& instruction) {
217 tim 278 int comparator = instruction.tok;
218     int property = instruction.intValue;
219 tim 283 float comparisonValue = boost::any_cast<float>(instruction.value);
220     float propertyValue;
221 tim 770 OOPSEBitSet bs(nStuntDouble);
222 tim 295 bs.clearAll();
223 tim 283
224     SimInfo::MoleculeIterator mi;
225     Molecule* mol;
226     Molecule::AtomIterator ai;
227     Atom* atom;
228     Molecule::RigidBodyIterator rbIter;
229     RigidBody* rb;
230 tim 281
231 cli2 1364 for (mol = info->beginMolecule(mi); mol != NULL;
232     mol = info->nextMolecule(mi)) {
233    
234 gezelter 507 for(atom = mol->beginAtom(ai); atom != NULL; atom = mol->nextAtom(ai)) {
235     compareProperty(atom, bs, property, comparator, comparisonValue);
236     }
237 cli2 1364
238     for (rb = mol->beginRigidBody(rbIter); rb != NULL;
239     rb = mol->nextRigidBody(rbIter)) {
240     compareProperty(rb, bs, property, comparator, comparisonValue);
241     }
242 tim 278 }
243    
244 tim 283 return bs;
245 gezelter 507 }
246 tim 278
247 cli2 1364 void SelectionEvaluator::compareProperty(StuntDouble* sd, OOPSEBitSet& bs,
248     int property, int comparator,
249     float comparisonValue) {
250 tim 963 RealType propertyValue = 0.0;
251 gezelter 507 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 cli2 1360 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 gezelter 507 default:
278     unrecognizedAtomProperty(property);
279     }
280 tim 283
281 gezelter 507 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 cli2 1364
305 gezelter 507 }
306 tim 283
307 cli2 1364 void SelectionEvaluator::withinInstruction(const Token& instruction,
308     OOPSEBitSet& bs){
309 tim 295
310 tim 281 boost::any withinSpec = instruction.value;
311 tim 295 float distance;
312 tim 282 if (withinSpec.type() == typeid(float)){
313 gezelter 507 distance = boost::any_cast<float>(withinSpec);
314 tim 295 } else if (withinSpec.type() == typeid(int)) {
315 gezelter 507 distance = boost::any_cast<int>(withinSpec);
316 tim 295 } else {
317 gezelter 507 evalError("casting error in withinInstruction");
318     bs.clearAll();
319 tim 278 }
320 tim 282
321 tim 295 bs = distanceFinder.find(bs, distance);
322 gezelter 507 }
323 cli2 1364
324 gezelter 507 void SelectionEvaluator::define() {
325 tim 282 assert(statement.size() >= 3);
326 cli2 1364
327 tim 282 std::string variable = boost::any_cast<std::string>(statement[1].value);
328 cli2 1364
329     variables.insert(VariablesType::value_type(variable,
330     expression(statement, 2)));
331 gezelter 507 }
332 cli2 1364
333 tim 278
334 gezelter 507 /** @todo */
335     void SelectionEvaluator::predefine(const std::string& script) {
336 cli2 1364
337 tim 282 if (compiler.compile("#predefine", script)) {
338 gezelter 507 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 cli2 1364 if (tok == Token::identifier ||
348     (tok & Token::predefinedset) == Token::predefinedset) {
349 gezelter 507 std::string variable = boost::any_cast<std::string>(statement[1].value);
350     variables.insert(VariablesType::value_type(variable, statement));
351 cli2 1364
352 gezelter 507 } else {
353     evalError("invalid variable name:" + script);
354     }
355     }else {
356     evalError("bad predefinition length:" + script);
357 cli2 1364 }
358 tim 282
359     } else {
360 gezelter 507 evalError("predefined set compile error:" + script +
361     "\ncompile error:" + compiler.getErrorMessage());
362 tim 282 }
363 gezelter 507 }
364 tim 282
365 tim 770 void SelectionEvaluator::select(OOPSEBitSet& bs){
366 tim 288 bs = expression(statement, 1);
367 gezelter 507 }
368 cli2 1364
369 tim 770 OOPSEBitSet SelectionEvaluator::lookupValue(const std::string& variable){
370 cli2 1364
371 tim 770 OOPSEBitSet bs(nStuntDouble);
372 tim 282 std::map<std::string, boost::any>::iterator i = variables.find(variable);
373 tim 295
374 tim 282 if (i != variables.end()) {
375 tim 770 if (i->second.type() == typeid(OOPSEBitSet)) {
376     return boost::any_cast<OOPSEBitSet>(i->second);
377 gezelter 507 } 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 tim 283 } else {
383 gezelter 507 unrecognizedIdentifier(variable);
384 tim 282 }
385 cli2 1364
386 tim 295 return bs;
387 gezelter 507 }
388 cli2 1364
389     OOPSEBitSet SelectionEvaluator::nameInstruction(const std::string& name){
390     return nameFinder.match(name);
391 gezelter 507 }
392 tim 283
393 gezelter 507 bool SelectionEvaluator::containDynamicToken(const std::vector<Token>& tokens){
394 tim 288 std::vector<Token>::const_iterator i;
395     for (i = tokens.begin(); i != tokens.end(); ++i) {
396 gezelter 507 if (i->tok & Token::dynamic) {
397     return true;
398     }
399 tim 288 }
400 cli2 1364
401 tim 288 return false;
402 gezelter 507 }
403 tim 288
404 gezelter 507 void SelectionEvaluator::clearDefinitionsAndLoadPredefined() {
405 tim 288 variables.clear();
406     //load predefine
407     //predefine();
408 gezelter 507 }
409 tim 288
410 tim 770 OOPSEBitSet SelectionEvaluator::evaluate() {
411     OOPSEBitSet bs(nStuntDouble);
412 tim 288 if (isLoaded_) {
413 gezelter 507 pc = 0;
414     instructionDispatchLoop(bs);
415 tim 288 }
416    
417     return bs;
418 gezelter 507 }
419 tim 295
420 tim 770 OOPSEBitSet SelectionEvaluator::indexInstruction(const boost::any& value) {
421     OOPSEBitSet bs(nStuntDouble);
422 tim 295
423     if (value.type() == typeid(int)) {
424 gezelter 507 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 tim 295 } else if (value.type() == typeid(std::pair<int, int>)) {
431 gezelter 507 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 tim 295 }
439    
440     return bs;
441 gezelter 507 }
442 tim 295
443 tim 432
444 tim 963 RealType SelectionEvaluator::getCharge(Atom* atom) {
445     RealType charge =0.0;
446 tim 432 AtomType* atomType = atom->getAtomType();
447     if (atomType->isCharge()) {
448 gezelter 507 GenericData* data = atomType->getPropertyByName("Charge");
449     if (data != NULL) {
450     DoubleGenericData* doubleData= dynamic_cast<DoubleGenericData*>(data);
451 tim 432
452 gezelter 507 if (doubleData != NULL) {
453     charge = doubleData->getData();
454 tim 432
455 gezelter 507 } 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 tim 432 }
464    
465     return charge;
466 gezelter 507 }
467 tim 432
468     }