00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #define NOMINMAX 1 // prevent MSVS from defining min/max macros
00024
00025 #include <stdio.h>
00026 #include <stdlib.h>
00027 #include <string.h>
00028 #include <math.h>
00029
00030 #include "AtomParser.h"
00031 #include "y.tab.h"
00032 #include "ParseTree.h"
00033 #include "Inform.h"
00034 #include "JRegex.h"
00035 #include "AtomSel.h"
00036 #include "Timestep.h"
00037 #include "DrawMolecule.h"
00038 #include "SpatialSearch.h"
00039
00040 #include <vector>
00041 #include <algorithm>
00042
00043
00044
00045 #define case_compare_numeric_macro(switchcase, symbol) \
00046 case switchcase: \
00047 l->convert(SymbolTableElement::IS_FLOAT); \
00048 r->convert(SymbolTableElement::IS_FLOAT); \
00049 ldval = l->dval; \
00050 rdval = r->dval; \
00051 flg = flgs; \
00052 for (i=num-1; i>=0; i--) { \
00053 *flg &= (*ldval symbol *rdval); \
00054 ldval += lincr; rdval += rincr; flg++; \
00055 } \
00056 break;
00057
00058 #define case_compare_string_macro(switchcase, symbol) \
00059 case switchcase: \
00060 l->convert(SymbolTableElement::IS_STRING); \
00061 r->convert(SymbolTableElement::IS_STRING); \
00062 lsptr = l->sval; \
00063 rsptr = r->sval; \
00064 flg = flgs; \
00065 for (i=num-1; i>=0; i--) { \
00066 if (*flg) \
00067 *flg &= (strcmp(*lsptr, *rsptr) symbol 0); \
00068 lsptr += lincr; rsptr += rincr; flg++; \
00069 } \
00070 break;
00071
00072
00074 ParseTree::ParseTree(SymbolTable *parser, atomparser_node *parse_tree)
00075 {
00076 tree = parse_tree;
00077 table = parser;
00078 selected_array = NULL;
00079 num_selected = 0;
00080 context = NULL;
00081 }
00082
00083 ParseTree::~ParseTree(void) {
00084 if (selected_array != NULL)
00085 delete [] selected_array;
00086 delete tree;
00087 }
00088
00089 void ParseTree::eval_compare(atomparser_node *node, int num, int *flgs) {
00090 int i;
00091 double *ldval, *rdval;
00092 char **lsptr, **rsptr;
00093 int lincr, rincr;
00094 int *flg;
00095
00096
00097 symbol_data *l = eval(node->left, num, flgs);
00098 symbol_data *r = eval(node->right, num, flgs);
00099
00100
00101
00102
00103 lincr = l->num == num ? 1 : 0;
00104 rincr = r->num == num ? 1 : 0;
00105
00106 switch (node->ival) {
00107 case_compare_numeric_macro(NLT, < )
00108 case_compare_numeric_macro(NLE, <= )
00109 case_compare_numeric_macro(NEQ, == )
00110 case_compare_numeric_macro(NGE, >= )
00111 case_compare_numeric_macro(NGT, > )
00112 case_compare_numeric_macro(NNE, != )
00113
00114 case_compare_string_macro(SLT, < )
00115 case_compare_string_macro(SLE, <= )
00116 case_compare_string_macro(SEQ, == )
00117 case_compare_string_macro(SGE, >= )
00118 case_compare_string_macro(SGT, > )
00119 case_compare_string_macro(SNE, != )
00120
00121 case MATCH: {
00122 l->convert(SymbolTableElement::IS_STRING);
00123 r->convert(SymbolTableElement::IS_STRING);
00124 lsptr = l->sval;
00125 rsptr = r->sval;
00126 flg = flgs;
00127 JRegex *rgx = NULL;
00128 const char *first = *rsptr;
00129
00130 for (i=0; i<num; i++) {
00131 if (i==0 || strcmp(*rsptr, first)) {
00132 if (rgx)
00133 delete rgx;
00134 rgx = new JRegex(*rsptr);
00135 first = *rsptr;
00136 }
00137 if (rgx) {
00138 if (*flg)
00139 *flg &= (rgx->match(*lsptr, strlen(*lsptr)) != -1);
00140 } else {
00141 *flg = 0;
00142 }
00143 lsptr += lincr; rsptr += rincr; flg++;
00144 }
00145 if (rgx) {
00146 delete rgx;
00147 }
00148
00149 break;
00150 }
00151
00152 default:
00153 msgWarn << "ParseTree::eval_compare() missing operator!" << sendmsg;
00154 }
00155
00156 delete l;
00157 delete r;
00158 }
00159
00160
00161
00162 symbol_data * ParseTree::eval_mathop(atomparser_node *node, int num, int *flgs)
00163 {
00164 symbol_data *l = eval(node->left, num, flgs);
00165 symbol_data *r = eval(node->right, num, flgs);
00166
00167
00168
00169 int lincr = l->num == num ? (~0) : 0;
00170 int rincr = r->num == num ? (~0) : 0;
00171 l->convert(SymbolTableElement::IS_FLOAT);
00172 r->convert(SymbolTableElement::IS_FLOAT);
00173 symbol_data *tmp = new symbol_data(SymbolTableElement::IS_FLOAT, num);
00174 int i;
00175 const double *lval = l->dval;
00176 const double *rval = r->dval;
00177 double *tmpval = tmp->dval;
00178
00179
00180
00181
00182
00183
00184 int firstsel = 0;
00185 int lastsel = -1;
00186 if (!find_first_selection_aligned(num, flgs, &firstsel) &&
00187 !find_last_selection_aligned(num, flgs, &lastsel)) {
00188
00189
00190 switch (node->node_type) {
00191 case ADD:
00192 for (i=firstsel; i<=lastsel; i++) {
00193 if (flgs[i]) tmpval[i] = lval[lincr & i] + rval[rincr & i];
00194 }
00195 break;
00196 case SUB:
00197 for (i=firstsel; i<=lastsel; i++) {
00198 if (flgs[i]) tmpval[i] = lval[lincr & i] - rval[rincr & i];
00199 }
00200 break;
00201 case MULT:
00202 for (i=firstsel; i<=lastsel; i++) {
00203 if (flgs[i]) tmpval[i] = lval[lincr & i] * rval[rincr & i];
00204 }
00205 break;
00206 case DIV:
00207 for (i=firstsel; i<=lastsel; i++) {
00208 if (flgs[i]) tmpval[i] = lval[lincr & i] / rval[rincr & i];
00209 }
00210 break;
00211 case MOD:
00212 for (i=firstsel; i<=lastsel; i++) {
00213 if (flgs[i]) tmpval[i] = fmod(lval[lincr & i], rval[rincr & i]);
00214 }
00215 break;
00216 case EXP:
00217 for (i=firstsel; i<=lastsel; i++) {
00218 if (flgs[i]) tmpval[i] = pow(lval[lincr & i], rval[rincr & i]);
00219 }
00220 break;
00221 }
00222 }
00223
00224 delete l;
00225 delete r;
00226 return tmp;
00227 }
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238 void ParseTree::eval_stringfctn(atomparser_node *node, int num, int *flgs) {
00239 int count = 0;
00240 atomparser_node *left;
00241
00242
00243 for (left = node->left; left != NULL; left = left->left) {
00244 count++;
00245 }
00246 if (count == 0)
00247 return;
00248
00249
00250 char **argv= (char **) calloc(1, count * sizeof(char *));
00251 int *types = new int[count];
00252 int i=0;
00253 for (left = node->left; left != NULL; left = left -> left, i++) {
00254
00255 switch (left->sele.st) {
00256 case RAW_STRING:
00257 types[i] = 0;
00258 argv[i] = (char *) ((const char *) left->sele.s);
00259 break;
00260
00261 case SQ_STRING:
00262 types[i] = 1;
00263 argv[i] = (char *) ((const char *) left->sele.s);
00264 break;
00265
00266 case DQ_STRING:
00267 types[i] = 2;
00268 argv[i] = (char *) ((const char *) left->sele.s);
00269 break;
00270 }
00271
00272 if (left->extra_type != -1) {
00273 types[i] += 3;
00274 }
00275 }
00276
00277
00278 int *tmp_flgs = new int[num];
00279 memcpy(tmp_flgs, flgs, num * sizeof(int));
00280 SymbolTableElement *elem = table->fctns.data(node->extra_type);
00281 elem->keyword_stringfctn(context, count, (const char **)argv, types, num, tmp_flgs);
00282
00283
00284 for (i = num-1; i>=0; i--) {
00285 if (flgs[i]) flgs[i] = tmp_flgs[i];
00286 }
00287 delete [] tmp_flgs;
00288 delete [] types;
00289 free(argv);
00290 }
00291
00292 static void same_string(symbol_data *tmp, symbol_data *tmp2, int num,
00293 int *subselect, int *flgs) {
00294 hash_t hash;
00295 hash_init(&hash, num);
00296
00297
00298 int i;
00299 for (i=0; i<num; i++)
00300 if (subselect[i])
00301 hash_insert(&hash, tmp2->sval[i], 0);
00302
00303
00304
00305
00306 for (i=0; i<num; i++)
00307 if (flgs[i])
00308 flgs[i] = (hash_lookup(&hash, tmp->sval[i]) != HASH_FAIL);
00309
00310 hash_destroy(&hash);
00311 }
00312
00313 static void same_int(symbol_data *tmp, symbol_data *tmp2, int num,
00314 int *subselect, int *flgs) {
00315 int firstsubsel = -1;
00316 if (find_first_selection_aligned(num, subselect, &firstsubsel)) {
00317
00318 memset(flgs, 0, num*sizeof(int));
00319 return;
00320 }
00321
00322 int lastsubsel = -1;
00323 if (find_last_selection_aligned(num, subselect, &lastsubsel)) {
00324
00325
00326 msgErr << "internal inconsistency in first/last selection search results"
00327 << sendmsg;
00328
00329 memset(flgs, 0, num*sizeof(int));
00330 return;
00331 }
00332
00333
00334
00335
00336
00337
00338 int *int_table = NULL, int_min, int_max;
00339
00340
00341 int_min = int_max = tmp2->ival[firstsubsel];
00342 int i;
00343 for (i=firstsubsel; i<=lastsubsel; i++) {
00344 if (subselect[i]) {
00345 int ival = tmp2->ival[i];
00346 if (ival > int_max)
00347 int_max = ival;
00348 if (ival < int_min)
00349 int_min = ival;
00350 }
00351 }
00352
00353 int_table = (int *) calloc(1+int_max-int_min, sizeof(int));
00354 for (i=firstsubsel; i<=lastsubsel; i++) {
00355 if (subselect[i]) {
00356 int_table[tmp2->ival[i]-int_min] = 1;
00357 }
00358 }
00359
00360
00361 for (i=0; i<num; i++) {
00362 if (flgs[i]) {
00363 int ival = tmp->ival[i];
00364 if (ival >= int_min && ival <= int_max)
00365 flgs[i] = int_table[ival-int_min];
00366 else
00367 flgs[i] = 0;
00368 }
00369 }
00370 delete [] int_table;
00371 }
00372
00373 static void same_double(symbol_data *tmp, symbol_data *tmp2, int num,
00374 int *subselect, int *flgs) {
00375 int firstsubsel, lastsubsel, subselselected;
00376 if (analyze_selection_aligned(num, subselect,
00377 &firstsubsel, &lastsubsel, &subselselected)) {
00378 return;
00379 }
00380
00381
00382
00383 hash_t hash;
00384
00385 char *doublestring = new char[25*subselselected];
00386 char *istring = doublestring;
00387 hash_init(&hash, subselselected);
00388 int i;
00389 for (i=firstsubsel; i<=lastsubsel; i++) {
00390 if (subselect[i]) {
00391 sprintf(istring,"%f", (double) tmp2->dval[i]);
00392 hash_insert(&hash, istring, 0);
00393 istring += 25;
00394 }
00395 }
00396
00397 char tmpstring[25];
00398 for (i=0; i<num; i++) {
00399 sprintf(tmpstring,"%f", (double) tmp->dval[i]);
00400 flgs[i] &= (hash_lookup(&hash, tmpstring) != HASH_FAIL);
00401 }
00402 hash_destroy(&hash);
00403
00404 delete [] doublestring;
00405 }
00406
00407
00408
00409
00410
00411
00412 void ParseTree::eval_same(atomparser_node *node, int num, int *flgs) {
00413 int i;
00414 int *subselect = new int[num];
00415 for (i=0; i<num; i++)
00416 subselect[i]=1;
00417
00418
00419 if (eval(node->left, num, subselect)) {
00420 delete [] subselect;
00421 msgErr << "eval of a 'same' returned data when it shouldn't have"
00422 << sendmsg;
00423 return;
00424 }
00425
00426
00427
00428
00429 SymbolTableElement *elem = table->fctns.data(node->extra_type);
00430 SymbolTableElement::symtype has_type = elem->returns_a;
00431 symbol_data *tmp, *tmp2;
00432 tmp = new symbol_data(has_type, num);
00433 tmp2 = new symbol_data(has_type, num);
00434
00435
00436
00437 switch (has_type) {
00438 case SymbolTableElement::IS_INT:
00439 elem->keyword_int(context, num, tmp->ival, flgs);
00440 elem->keyword_int(context, num, tmp2->ival, subselect);
00441 same_int(tmp, tmp2, num, subselect, flgs);
00442 break;
00443
00444 case SymbolTableElement::IS_FLOAT:
00445 elem->keyword_double(context, num, tmp->dval, flgs);
00446 elem->keyword_double(context, num, tmp2->dval, subselect);
00447 same_double(tmp, tmp2, num, subselect, flgs);
00448 break;
00449
00450 case SymbolTableElement::IS_STRING:
00451 elem->keyword_string(context, num, (const char **)tmp->sval, flgs);
00452 elem->keyword_string(context, num, (const char **)tmp2->sval, subselect);
00453 same_string(tmp, tmp2, num, subselect, flgs);
00454 break;
00455 }
00456
00457 delete tmp;
00458 delete tmp2;
00459 delete [] subselect;
00460 }
00461
00462
00463
00464
00465 symbol_data *ParseTree::eval_key(atomparser_node *node, int num, int *flgs) {
00466
00467 SymbolTableElement *elem = table->fctns.data(node->extra_type);
00468 SymbolTableElement::symtype has_type = elem->returns_a;
00469 symbol_data *tmp;
00470 tmp = new symbol_data(has_type, num);
00471
00472 switch (has_type) {
00473 case SymbolTableElement::IS_INT:
00474 elem->keyword_int(context, num, tmp->ival, flgs);
00475 break;
00476 case SymbolTableElement::IS_FLOAT:
00477 elem->keyword_double(context, num, tmp->dval, flgs);
00478 break;
00479 case SymbolTableElement::IS_STRING:
00480 elem->keyword_string(context, num, (const char **)tmp->sval, flgs);
00481 break;
00482 }
00483
00484
00485
00486 int *int_table = NULL;
00487
00488
00489
00490
00491 int firstsel = 0;
00492 int lastsel = -1;
00493 if (find_first_selection_aligned(num, flgs, &firstsel)) {
00494
00495
00496 firstsel=0;
00497 } else {
00498 if (find_last_selection_aligned(num, flgs, &lastsel)) {
00499 msgErr << "internal inconsistency in first/last selection search results"
00500 << sendmsg;
00501 }
00502 }
00503
00504 int int_min=0, int_max=0;
00505 if (has_type == SymbolTableElement::IS_INT) {
00506 if (lastsel != -1) {
00507 int_min = int_max = tmp->ival[firstsel];
00508 }
00509
00510
00511
00512
00513 for (int i=firstsel; i<=lastsel; i++) {
00514 if (flgs[i]) {
00515 const int ival = tmp->ival[i];
00516 if (ival > int_max)
00517 int_max = ival;
00518 if (ival < int_min)
00519 int_min = ival;
00520 }
00521 }
00522 int_table = (int *) calloc(1+int_max-int_min, sizeof(int));
00523 }
00524
00525
00526
00527
00528
00529
00530
00531
00532
00533
00534 int *newflgs = new int[num];
00535
00536 memset(newflgs, 0, num*sizeof(int));
00537
00538 if (node->left) {
00539 atomparser_node *left = node->left;
00540 while (left) {
00541 if (left->extra_type == -1) {
00542
00543 switch(has_type) {
00544 case SymbolTableElement::IS_INT:
00545 {
00546 int ival = atoi(left->sele.s);
00547 if (ival >= int_min && ival <= int_max)
00548 int_table[ival-int_min] = 1;
00549 }
00550 break;
00551
00552 case SymbolTableElement::IS_FLOAT:
00553 {
00554
00555 double dval = atof(left->sele.s);
00556 double delta = fabs(dval / 1000);
00557 double maxval = dval+delta;
00558 double minval = dval-delta;
00559 for (int i=firstsel; i<=lastsel; i++) {
00560 if (flgs[i])
00561 newflgs[i] |= (minval <= tmp->dval[i] && maxval >= tmp->dval[i]);
00562 }
00563 }
00564 break;
00565
00566 case SymbolTableElement::IS_STRING:
00567 {
00568 switch (left->sele.st) {
00569 case SQ_STRING:
00570 case RAW_STRING:
00571 {
00572 for (int i=firstsel; i<=lastsel; i++) {
00573
00574
00575
00576 if (flgs[i] && (tmp->sval[i] != NULL)) {
00577 newflgs[i] |= !strcmp(left->sele.s, tmp->sval[i]);
00578 }
00579 }
00580 }
00581 break;
00582
00583 case DQ_STRING:
00584 default:
00585 {
00586
00587
00588
00589
00590
00591
00592 JString temps = "^("+left->sele.s+")$";
00593 JRegex r(temps, 1);
00594 for (int i=firstsel; i<=lastsel; i++) {
00595 if (flgs[i]) {
00596 newflgs[i] |= (r.match(tmp->sval[i], strlen(tmp->sval[i])) != -1);
00597 }
00598 }
00599 }
00600 break;
00601 }
00602 }
00603 }
00604 } else {
00605 switch(has_type) {
00606 case SymbolTableElement::IS_INT:
00607 {
00608 int ltval = atoi(left->sele.s);
00609 int gtval = atoi(left->left->sele.s);
00610 if (ltval < int_min) ltval = int_min;
00611 if (gtval > int_max) gtval = int_max;
00612 for (int i=ltval-int_min; i<= gtval-int_min; i++)
00613 int_table[i] = 1;
00614 }
00615 break;
00616 case SymbolTableElement::IS_FLOAT:
00617 {
00618 double ltval = atof(left->sele.s);
00619 double gtval = atof(left->left->sele.s);
00620 for (int i=firstsel; i<=lastsel; i++) {
00621 if (flgs[i])
00622 newflgs[i] |= ((ltval <= tmp->dval[i]) && (gtval >= tmp->dval[i]));
00623 }
00624 }
00625 break;
00626 default:
00627 {
00628
00629 for (int i=firstsel; i<=lastsel; i++) {
00630 if (flgs[i])
00631 newflgs[i] |= (flgs[i] && strcmp(left->sele.s, tmp->sval[i]) <= 0
00632 && strcmp(left->left->sele.s, tmp->sval[i]) >= 0);
00633
00634 }
00635 }
00636 }
00637 left = left->left;
00638 }
00639 left = left->left;
00640 }
00641
00642
00643 if (has_type == SymbolTableElement::IS_INT) {
00644 for (int i=firstsel; i<=lastsel; i++) {
00645 if (flgs[i])
00646 flgs[i] = int_table[tmp->ival[i]-int_min];
00647 }
00648 free(int_table);
00649 } else {
00650 for (int i=firstsel; i<=lastsel; i++) {
00651 if (flgs[i])
00652 flgs[i] = newflgs[i];
00653 }
00654 }
00655
00656 firstsel=0;
00657 lastsel=num-1;
00658
00659 delete [] newflgs;
00660 delete tmp;
00661 return NULL;
00662 } else {
00663
00664
00665
00666 delete [] newflgs;
00667 if (int_table) free(int_table);
00668 return tmp;
00669 }
00670 }
00671
00672
00673 void ParseTree::eval_single(atomparser_node *node, int num, int *flgs) {
00674
00675
00676 atomsel_ctxt *ctxt = (atomsel_ctxt *)context;
00677 ctxt->singleword = table->fctns.name(node->extra_type);
00678 table->fctns.data(node->extra_type)->keyword_single(context, num, flgs);
00679 }
00680
00681
00682 void ParseTree::eval_exwithin(atomparser_node *node, int num, int *flgs) {
00683 eval_within(node, num, flgs);
00684
00685
00686 int *others = new int[num];
00687 int i;
00688 for (i=0; i<num; others[i++] = 1);
00689
00690
00691 if (eval(node->left, num, others)) {
00692 delete [] others;
00693 msgErr << "eval of a 'within' returned data when it shouldn't have." << sendmsg;
00694 return;
00695 }
00696 for (i=0; i<num; i++) {
00697 if (others[i]) flgs[i] = 0;
00698 }
00699 delete [] others;
00700 }
00701
00702
00703
00704 static Timestep *selframe(DrawMolecule *atom_sel_mol, int which_frame) {
00705 switch (which_frame) {
00706 case AtomSel::TS_LAST: return atom_sel_mol->get_last_frame();
00707 case AtomSel::TS_NOW : return atom_sel_mol->current();
00708 default: {
00709 if (!atom_sel_mol->get_frame(which_frame)) {
00710 return atom_sel_mol->get_last_frame();
00711
00712 } else {
00713 return atom_sel_mol->get_frame(which_frame);
00714 }
00715 }
00716 }
00717 return NULL;
00718 }
00719
00720 void ParseTree::eval_pbwithin(atomparser_node *node, int num, int *flgs) {
00721
00722
00723 if ((float) node->dval <= 0.0f) {
00724 eval(node->left, num, flgs);
00725 return;
00726 }
00727
00728
00729
00730
00731 int i;
00732
00733
00734 ResizeArray<float> coords(3*2*num);
00735
00736
00737
00738 ResizeArray<int> others(2*num);
00739
00740
00741
00742
00743 ResizeArray<int> repflgs(2*num);
00744
00745
00746 ResizeArray<int> repindexes(num);
00747
00748
00749 atomsel_ctxt *ctxt = (atomsel_ctxt *)context;
00750 const Timestep *ts = selframe(ctxt->atom_sel_mol, ctxt->which_frame);
00751 if (!ts) {
00752 msgErr << "No timestep available for 'within' search!" << sendmsg;
00753 return;
00754 }
00755 for (i=0; i<num; ++i) others.append(1);
00756 if (eval(node->left, num, &others[0])) {
00757 msgErr << "eval of a 'within' returned data when it shouldn't have." << sendmsg;
00758 return;
00759 }
00760
00761
00762 const float * pos=ts->pos;
00763 for (i=0; i<num; i++) {
00764 coords.append(pos[0]);
00765 coords.append(pos[1]);
00766 coords.append(pos[2]);
00767 pos += 3;
00768 repflgs.append(flgs[i]);
00769 }
00770
00771
00772 float min[3], max[3];
00773 if (!find_minmax_selected(num, &others[0], ts->pos,
00774 min[0], min[1], min[2], max[0], max[1], max[2])) {
00775 memset(flgs, 0, num*sizeof(int));
00776 return;
00777 }
00778
00779
00780 const float cutoff = (float)node->dval;
00781 for (i=0; i<3; i++) {
00782 min[i] -= cutoff;
00783 max[i] += cutoff;
00784 }
00785
00786
00787 float A[3], B[3], C[3];
00788 ts->get_transform_vectors(A, B, C);
00789 for (i=-1; i<=1; i++) {
00790 float v1[3];
00791 vec_scale(v1, (float) i, A);
00792 for (int j=-1; j<=1; j++) {
00793 float v2[3];
00794 vec_scale(v2, (float) j, B);
00795 for (int k=-1; k<=1; k++) {
00796
00797 if (!i && !j && !k) continue;
00798 float v3[3];
00799 vec_scale(v3, (float) k, C);
00800 float vx = v1[0] + v2[0] + v3[0];
00801 float vy = v1[1] + v2[1] + v3[1];
00802 float vz = v1[2] + v2[2] + v3[2];
00803 pos = ts->pos;
00804 for (int ind=0; ind<num; ind++) {
00805 if (flgs[ind]) {
00806 const float x = pos[0] + vx;
00807 const float y = pos[1] + vy;
00808 const float z = pos[2] + vz;
00809 if (x>min[0] && x<=max[0] &&
00810 y>min[1] && y<=max[1] &&
00811 z>min[2] && z<=max[2]) {
00812 repindexes.append(ind);
00813 coords.append(x);
00814 coords.append(y);
00815 coords.append(z);
00816 repflgs.append(1);
00817 others.append(0);
00818 }
00819 }
00820 pos += 3;
00821 }
00822 }
00823 }
00824 }
00825
00826 find_within(&coords[0], &repflgs[0], &others[0], others.num(), cutoff);
00827
00828
00829 memcpy(flgs, &repflgs[0], num*sizeof(int));
00830
00831
00832 for (i=0; i<repindexes.num(); i++) {
00833 if (repflgs[num+i]) {
00834 flgs[repindexes[i]] = 1;
00835 }
00836 }
00837 }
00838
00839
00840 void ParseTree::eval_within(atomparser_node *node, int num, int *flgs) {
00841
00842 if ((float) node->dval > 0.0f) {
00843
00844 int *others = new int[num];
00845 int i;
00846 for (i=0; i<num; ++i)
00847 others[i] = 1;
00848
00849 if (eval(node->left, num, others)) {
00850 delete [] others;
00851 msgErr << "eval of a 'within' returned data when it shouldn't have." << sendmsg;
00852 return;
00853 }
00854
00855
00856 atomsel_ctxt *ctxt = (atomsel_ctxt *)context;
00857 Timestep *ts = selframe(ctxt->atom_sel_mol, ctxt->which_frame);
00858 if (!ts) {
00859 msgErr << "No timestep available for 'within' search!" << sendmsg;
00860 return;
00861 }
00862
00863 find_within(ts->pos, flgs, others, num, (float) node->dval);
00864
00865 delete [] others;
00866 } else {
00867
00868
00869 eval(node->left, num, flgs);
00870 }
00871 }
00872
00873
00874 void ParseTree::eval_within_bonds(atomparser_node *node, int num, int *flgs) {
00875 atomsel_ctxt *ctxt = (atomsel_ctxt *)context;
00876 int *others = new int[num];
00877
00878 int i;
00879 for (i=0; i<num; ++i)
00880 others[i] = 1;
00881
00882 if (eval(node->left, num, others)) {
00883 delete [] others;
00884 msgErr << "eval of a 'within' returned data when it shouldn't have." << sendmsg;
00885 return;
00886 }
00887
00888
00889 int *bondedsel = new int[num];
00890 memcpy(bondedsel, others, num*sizeof(int));
00891
00892
00893 int bonddepth;
00894 for (bonddepth=0; bonddepth<node->ival; bonddepth++) {
00895
00896 for (i=0; i<num; i++) {
00897 if (bondedsel[i]) {
00898 MolAtom *atom = ctxt->atom_sel_mol->atom(i);
00899 int j;
00900 for (j=0; j<atom->bonds; j++) {
00901 others[atom->bondTo[j]] = 1;
00902 }
00903 }
00904 }
00905
00906
00907 memcpy(bondedsel, others, num*sizeof(int));
00908 }
00909
00910 for (i=0; i<num; ++i)
00911 if (bondedsel[i] && flgs[i])
00912 flgs[i] = 1;
00913 else
00914 flgs[i] = 0;
00915
00916 delete [] bondedsel;
00917 delete [] others;
00918 }
00919
00920
00921
00922
00923
00924
00925
00926
00927
00928
00929 namespace {
00930 struct PointDistance {
00931 float o;
00932 int i;
00933 PointDistance() {}
00934 PointDistance(float o_, int i_) : o(o_), i(i_) {}
00935 bool operator<(const PointDistance& p) const {
00936 return o<p.o;
00937 }
00938 };
00939 }
00940
00941 void ParseTree::eval_k_nearest(atomparser_node *node, int num, int *flgs) {
00942 atomsel_ctxt *ctxt = (atomsel_ctxt *)context;
00943 const Timestep *ts = selframe(ctxt->atom_sel_mol, ctxt->which_frame);
00944 if (!ts) {
00945 msgErr << "No timestep available for 'nearest' search!" << sendmsg;
00946 return;
00947 }
00948 const int N = node->ival;
00949 int i;
00950
00951
00952 std::vector<int> others(num);
00953 for (i=0; i<num; ++i) others[i] = 1;
00954 if (eval(node->left, num, &others[0])) {
00955 msgErr << "eval of a 'within' returned data when it shouldn't have." << sendmsg;
00956 return;
00957 }
00958
00959
00960 for (i=0; i<num; i++) if (others[i]) break;
00961 if (i==num) {
00962 memset(flgs,0,num*sizeof(*flgs));
00963 return;
00964 }
00965
00966 std::vector<PointDistance> distances;
00967 for (i=0; i<num; i++) {
00968 if (others[i] || !flgs[i]) continue;
00969 #if 1
00970 float d2=1e37f;
00971 #else
00972 float d2=std::numeric_limits<float>::max();
00973 #endif
00974 for (int j=0; j<num; j++) {
00975 if (!others[j]) continue;
00976 float d2_j=distance2(ts->pos+3*i, ts->pos+3*j);
00977 if (d2_j<d2) d2=d2_j;
00978 }
00979 distances.push_back(PointDistance(d2,i));
00980 }
00981 std::sort(distances.begin(), distances.end());
00982 int n=N;
00983 if (n>distances.size()) n=distances.size();
00984 memset(flgs,0,num*sizeof(*flgs));
00985 for (i=0; i<n; i++) {
00986 flgs[distances[i].i]=1;
00987 }
00988 }
00989
00990
00991
00992
00993
00994 void ParseTree::find_rings(int num, int *flgs, int *others,
00995 int minringsize, int maxringsize) {
00996 #ifdef VMDWITHCARBS
00997 int i;
00998
00999
01000
01001
01002 atomsel_ctxt *ctxt = (atomsel_ctxt *)context;
01003 ctxt->atom_sel_mol->find_small_rings_and_links(5, maxringsize);
01004 SmallRing *ring;
01005 memset(flgs, 0, num*sizeof(int));
01006 for (i=0; i < ctxt->atom_sel_mol->smallringList.num(); i++) {
01007 ring = ctxt->atom_sel_mol->smallringList[i];
01008 int N = ring->num();
01009 if (N >= minringsize && N <= maxringsize) {
01010 int j;
01011 for (j=0; j<N; j++) {
01012 int ind = (*ring)[j];
01013 flgs[ind] = others[ind];
01014 }
01015 }
01016 }
01017 #else
01018 memset(flgs, 0, num*sizeof(int));
01019 #endif
01020 }
01021
01022
01023 void ParseTree::eval_maxringsize(atomparser_node *node, int num, int *flgs) {
01024
01025 int *others = new int[num];
01026 int i;
01027 for (i=0; i<num; ++i)
01028 others[i] = 1;
01029
01030 if (eval(node->left, num, others)) {
01031 delete [] others;
01032 msgErr << "eval of a 'maxringsize' returned data when it shouldn't have." << sendmsg;
01033 return;
01034 }
01035
01036 find_rings(num, flgs, others, 1, node->ival);
01037
01038 delete [] others;
01039 }
01040
01041
01042 void ParseTree::eval_ringsize(atomparser_node *node, int num, int *flgs) {
01043
01044 int *others = new int[num];
01045 int i;
01046 for (i=0; i<num; ++i)
01047 others[i] = 1;
01048
01049 if (eval(node->left, num, others)) {
01050 delete [] others;
01051 msgErr << "eval of a 'ringsize' returned data when it shouldn't have." << sendmsg;
01052 return;
01053 }
01054
01055 find_rings(num, flgs, others, node->ival, node->ival);
01056
01057 delete [] others;
01058 }
01059
01060
01061
01062
01063 symbol_data *ParseTree::eval(atomparser_node *node, int num, int *flgs) {
01064 int i;
01065 int *flg1, *flg2;
01066 symbol_data *tmp;
01067 switch(node->node_type) {
01068 case AND:
01069 eval(node->left, num, flgs);
01070 eval(node->right, num, flgs);
01071 return NULL;
01072
01073 case NOT:
01074 flg1 = new int[num];
01075 memcpy(flg1, flgs, num*sizeof(int));
01076
01077 eval(node->left, num, flg1);
01078
01079 for (i=num-1; i>=0; i--) {
01080 if (flgs[i])
01081 flgs[i] = !flg1[i];
01082 }
01083 delete [] flg1;
01084 break;
01085
01086 case OR:
01087 flg1 = new int[num];
01088 memcpy(flg1, flgs, num*sizeof(int));
01089 eval(node->left, num, flg1);
01090 flg2 = new int[num];
01091 memcpy(flg2, flgs, num*sizeof(int));
01092 eval(node->right, num, flg2);
01093 for (i=num-1; i>=0; i--) {
01094 flgs[i] = flgs[i] && (flg1[i] || flg2[i]);
01095 }
01096 delete [] flg1;
01097 delete [] flg2;
01098 break;
01099
01100 case FLOATVAL:
01101 tmp = new symbol_data(SymbolTableElement::IS_FLOAT, 1);
01102 tmp->dval[0] = node->dval;
01103 return tmp;
01104
01105 case INTVAL:
01106 tmp = new symbol_data(SymbolTableElement::IS_INT, 1);
01107 tmp->ival[0] = node->ival;
01108 return tmp;
01109
01110 case STRWORD:
01111 tmp = new symbol_data(SymbolTableElement::IS_STRING, 1);
01112 tmp->sval[0] = (char *)(const char *)node->sele.s;
01113 return tmp;
01114
01115 case KEY:
01116 return eval_key(node, num, flgs);
01117
01118 case STRFCTN:
01119 eval_stringfctn(node, num, flgs);
01120 break;
01121
01122 case FUNC:
01123 {
01124
01125
01126
01127 symbol_data *inp = eval(node->left, num, flgs);
01128 inp->convert(SymbolTableElement::IS_FLOAT);
01129
01130
01131 symbol_data *ret = new symbol_data(SymbolTableElement::IS_FLOAT, num);
01132 SymbolTableElement *elem = table->fctns.data(node->extra_type);
01133
01134
01135
01136
01137 if (inp->num == num) {
01138 for (i=0; i<num; i++)
01139 ret->dval[i] = elem->fctn(inp->dval[i]);
01140 } else {
01141
01142
01143 double d = elem->fctn(inp->dval[0]);
01144 for (i=0; i<num; i++)
01145 ret->dval[i] = d;
01146 }
01147 delete inp;
01148 return ret;
01149 }
01150
01151 case ADD:
01152 case SUB:
01153 case MULT:
01154 case MOD:
01155 case EXP:
01156 case DIV:
01157 return eval_mathop(node, num, flgs);
01158
01159 case UMINUS:
01160 tmp = eval(node->left, num, flgs);
01161 tmp->convert(SymbolTableElement::IS_FLOAT);
01162 for (i=0; i<tmp->num; i++) {
01163 tmp->dval[i] = -tmp->dval[i];
01164 }
01165 return tmp;
01166
01167 case COMPARE:
01168 eval_compare(node, num, flgs);
01169 break;
01170
01171 case WITHIN:
01172 eval_within(node, num, flgs);
01173 break;
01174
01175 case EXWITHIN:
01176 eval_exwithin(node, num, flgs);
01177 break;
01178
01179 case PBWITHIN:
01180 eval_pbwithin(node, num, flgs);
01181 break;
01182
01183 #if defined(NEAREST)
01184 case NEAREST:
01185 eval_k_nearest(node, num, flgs);
01186 break;
01187 #endif
01188
01189 #if defined(WITHINBONDS)
01190 case WITHINBONDS:
01191 eval_within_bonds(node, num, flgs);
01192 break;
01193 #endif
01194
01195 #if defined(MAXRINGSIZE)
01196 case MAXRINGSIZE:
01197 eval_maxringsize(node, num, flgs);
01198 break;
01199 #endif
01200
01201 #if defined(RINGSIZE)
01202 case RINGSIZE:
01203 eval_ringsize(node, num, flgs);
01204 break;
01205 #endif
01206
01207 case SAME:
01208 eval_same(node, num, flgs);
01209 break;
01210
01211 case SINGLE:
01212 eval_single(node, num, flgs);
01213 break;
01214
01215 default:
01216 msgWarn << "ParseTree::eval() unknown node type: " << node->node_type << sendmsg;
01217 break;
01218 }
01219
01220 return NULL;
01221 }
01222
01223
01224
01225 void ParseTree::eval_find_recursion(atomparser_node *node, int *found,
01226 hash_t *hash) {
01227
01228
01229
01230 switch (node->node_type) {
01231 case AND:
01232 case OR:
01233 eval_find_recursion(node->left, found, hash);
01234 eval_find_recursion(node->right, found, hash);
01235
01236
01237 break;
01238
01239 case NOT:
01240 case UMINUS:
01241 case WITHIN:
01242 case EXWITHIN:
01243 case SAME:
01244 case FUNC:
01245 eval_find_recursion(node->left, found, hash);
01246 break;
01247
01248 case SINGLE:
01249 {
01250 const char *thisword = table->fctns.name(node->extra_type);
01251 const char *macro = table->get_custom_singleword(thisword);
01252 if (macro) {
01253 if (hash_insert(hash, thisword, 0) != HASH_FAIL) {
01254 *found = 1;
01255 } else {
01256 ParseTree *subtree = table->parse(macro);
01257 if (subtree != NULL) {
01258 eval_find_recursion(subtree->tree, found, hash);
01259 delete subtree;
01260 } else {
01261
01262
01263
01264
01265 msgErr << "ParseTree) internal processing error, NULL "
01266 << "subtree value while checking recursion" << sendmsg;
01267 }
01268 hash_delete(hash, thisword);
01269 }
01270 }
01271 }
01272 break;
01273 }
01274 }
01275
01276
01277
01278 int ParseTree::find_recursion(const char *head) {
01279 hash_t hash;
01280 hash_init(&hash, 10);
01281 hash_insert(&hash, head, 0);
01282 int found = 0;
01283 eval_find_recursion(tree, &found, &hash);
01284 hash_destroy(&hash);
01285 return found;
01286 }
01287
01288
01289
01290
01291
01292 int ParseTree::evaluate(int num_atoms, int *flgs) {
01293 int num = num_atoms;
01294 if (!tree || num < 0 ) {
01295 return 0;
01296 }
01297
01298
01299 for (int i=0; i<num; i++) {
01300 flgs[i] = 1;
01301 }
01302
01303
01304 symbol_data *retdat = eval(tree, num, flgs);
01305 if (retdat) {
01306 msgErr << "Atom selection returned data when it shouldn't\n" << sendmsg;
01307 delete retdat;
01308 }
01309
01310 return 1;
01311 }
01312
01313
01314
01315 void symbol_data::make_space(void) {
01316 free_space();
01317 switch(type) {
01318 case SymbolTableElement::IS_FLOAT:
01319 dval = new double[num];
01320 break;
01321
01322 case SymbolTableElement::IS_INT:
01323 ival = new int[num];
01324 break;
01325
01326 case SymbolTableElement::IS_STRING:
01327 sval = new char *[num];
01328 memset(sval, 0, num*sizeof(char *));
01329 break;
01330 }
01331 }
01332
01333
01334
01335 void symbol_data::free_space(void) {
01336 switch (type) {
01337 case SymbolTableElement::IS_FLOAT:
01338 if (dval)
01339 delete [] dval;
01340 dval = NULL;
01341 break;
01342
01343 case SymbolTableElement::IS_INT:
01344 if (ival)
01345 delete [] ival;
01346 ival = NULL;
01347 break;
01348
01349 case SymbolTableElement::IS_STRING:
01350 if (sval) {
01351
01352 if (free_sval)
01353 for (int i=0; i<num; i++) free(sval[i]);
01354
01355 delete [] sval;
01356 sval = NULL;
01357 }
01358 free_sval = 0;
01359 break;
01360
01361 default:
01362 msgErr << "Unknown data type " << (int)type
01363 << " in symbol_data::free_space" << sendmsg;
01364 }
01365 }
01366
01367
01368
01369 symbol_data::symbol_data(SymbolTableElement::symtype new_type, int new_num) {
01370 type = new_type;
01371 num = new_num;
01372 dval = NULL;
01373 ival = NULL;
01374 sval = NULL;
01375 free_sval = 0;
01376 make_space();
01377 }
01378
01379
01380 symbol_data::~symbol_data(void) {
01381 free_space();
01382 }
01383
01384
01385 void symbol_data::convert(SymbolTableElement::symtype totype) {
01386
01387 if (totype == type)
01388 return;
01389
01390
01391 if (totype == SymbolTableElement::IS_FLOAT) {
01392 double *tmp = new double[num];
01393 if (type == SymbolTableElement::IS_INT) {
01394 for (int i=num-1; i>=0; i--) {
01395 tmp[i] = (double) ival[i];
01396 }
01397 } else {
01398 for (int i=num-1; i>=0; i--) {
01399
01400
01401
01402
01403
01404
01405
01406
01407 if (sval[i] != NULL) {
01408 tmp[i] = atof(sval[i]);
01409 } else {
01410 for (int j=num-1; j>=0; j--) {
01411 tmp[i] = 0.0f;
01412 }
01413 msgErr << "ParseTree) internal processing error, NULL string value "
01414 << "while converting to floating point" << sendmsg;
01415 break;
01416 }
01417 }
01418 }
01419 free_space();
01420 type = totype;
01421 dval = tmp;
01422 return;
01423 }
01424
01425
01426 if (totype == SymbolTableElement::IS_STRING) {
01427 char **tmp = new char*[num];
01428 memset(tmp, 0, num*sizeof(char *));
01429 char s[100];
01430 if (type == SymbolTableElement::IS_INT) {
01431 for (int i=num-1; i>=0; i--) {
01432 sprintf(s, "%ld", (long) ival[i]);
01433 tmp[i] = strdup(s);
01434 }
01435 } else {
01436 for (int i=num-1; i>=0; i--) {
01437 sprintf(s, "%f", (double) dval[i]);
01438 tmp[i] = strdup(s);
01439 }
01440 }
01441 free_space();
01442 type = totype;
01443 sval = tmp;
01444 free_sval = TRUE;
01445 return;
01446 }
01447
01448
01449 if (totype == SymbolTableElement::IS_INT) {
01450 int *tmp = new int[num];
01451 if (type == SymbolTableElement::IS_FLOAT) {
01452 for (int i=num-1; i>=0; i--) {
01453 tmp[i] = (int) dval[i];
01454 }
01455 } else {
01456 for (int i=num-1; i>=0; i--) {
01457
01458
01459
01460
01461
01462
01463
01464
01465 if (sval[i] != NULL) {
01466 tmp[i] = atoi(sval[i]);
01467 } else {
01468 for (int j=num-1; j>=0; j--) {
01469 tmp[i] = 0;
01470 }
01471 msgErr << "ParseTree) internal processing error, NULL string value "
01472 << "while converting to integer" << sendmsg;
01473 break;
01474 }
01475 }
01476 }
01477 free_space();
01478 type = totype;
01479 ival = tmp;
01480 return;
01481 }
01482 }
01483