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