ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/OpenMD/branches/development/src/selection/SelectionEvaluator.cpp
Revision: 1665
Committed: Tue Nov 22 20:38:56 2011 UTC (13 years, 5 months ago) by gezelter
File size: 14278 byte(s)
Log Message:
updated copyright notices

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

Properties

Name Value
svn:keywords Author Id Revision Date