3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
9 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
10 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
11 * Copyright (c) 2001-2009, The GROMACS development team,
12 * check out http://www.gromacs.org for more information.
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
19 * If you want to redistribute modifications, please consider that
20 * scientific software is very special. Version control is crucial -
21 * bugs must be traceable. We will be happy to consider code for
22 * inclusion in the official distribution, but derived work must not
23 * be called official GROMACS. Details are found in the README & COPYING
24 * files - if they are missing, get the official version at www.gromacs.org.
26 * To help us fund GROMACS development, we humbly ask that you cite
27 * the papers on the package - you can find them in the top README file.
29 * For more info, check our website at http://www.gromacs.org
33 * Implements internal selection method for comparison expressions.
35 * \author Teemu Murtola <teemu.murtola@cbr.su.se>
36 * \ingroup module_selection
38 #include "gromacs/legacyheaders/maths.h"
39 #include "gromacs/legacyheaders/macros.h"
40 #include "gromacs/legacyheaders/smalloc.h"
42 #include "gromacs/selection/selmethod.h"
43 #include "gromacs/utility/exceptions.h"
45 /** Defines the comparison operator for comparison expressions. */
48 CMP_INVALID, /**< Indicates an error */
53 CMP_EQUAL, /**< '==' */
57 /** The operand has a single value. */
58 #define CMP_SINGLEVAL 1
59 /** The operand value is dynamic. */
60 #define CMP_DYNAMICVAL 2
61 /** The value is real. */
63 /** The integer array is allocated. */
64 #define CMP_ALLOCINT 16
65 /** The real array is allocated. */
66 #define CMP_ALLOCREAL 32
69 * Data structure for comparison expression operand values.
73 /** Flags that describe the type of the operand. */
75 /** (Array of) integer value(s). */
77 /** (Array of) real value(s). */
82 * Data structure for comparison expression evaluation.
86 /** Comparison operator as a string. */
88 /** Comparison operator type. */
93 t_compare_value right;
94 } t_methoddata_compare;
96 /** Allocates data for comparison expression evaluation. */
98 init_data_compare(int npar, gmx_ana_selparam_t *param);
99 /** Initializes data for comparison expression evaluation. */
101 init_compare(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
102 /** Frees the memory allocated for comparison expression evaluation. */
104 free_data_compare(void *data);
105 /** Evaluates comparison expressions. */
107 evaluate_compare(t_topology *top, t_trxframe *fr, t_pbc *pbc,
108 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
110 /** Parameters for comparison expression evaluation. */
111 static gmx_ana_selparam_t smparams_compare[] = {
112 {"int1", {INT_VALUE, -1, {NULL}}, NULL,
113 SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
114 {"real1", {REAL_VALUE, -1, {NULL}}, NULL,
115 SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
116 {"op", {STR_VALUE, 1, {NULL}}, NULL, 0},
117 {"int2", {INT_VALUE, -1, {NULL}}, NULL,
118 SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
119 {"real2", {REAL_VALUE, -1, {NULL}}, NULL,
120 SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
123 /** \internal Selection method data for comparison expression evaluation. */
124 gmx_ana_selmethod_t sm_compare = {
125 "cmp", GROUP_VALUE, SMETH_SINGLEVAL,
126 asize(smparams_compare), smparams_compare,
139 * Returns a \c e_comparison_t value corresponding to an operator.
141 * \param[in] str String to process.
142 * \returns The comparison type corresponding to the first one or two
143 * characters of \p str.
145 * \p str can contain any number of characters; only the first two
147 * If the beginning of \p str does not match any of the recognized types,
148 * \ref CMP_INVALID is returned.
150 static e_comparison_t
151 comparison_type(char *str)
155 case '<': return (str[1] == '=') ? CMP_LEQ : CMP_LESS;
156 case '>': return (str[1] == '=') ? CMP_GEQ : CMP_GTR;
157 case '=': return (str[1] == '=') ? CMP_EQUAL : CMP_INVALID;
158 case '!': return (str[1] == '=') ? CMP_NEQ : CMP_INVALID;
164 * Returns a string corresponding to a \c e_comparison_t value.
166 * \param[in] cmpt Comparison type to convert.
167 * \returns Pointer to a string that corresponds to \p cmpt.
169 * The return value points to a string constant and should not be \p free'd.
171 * The function returns NULL if \p cmpt is not one of the valid values.
174 comparison_type_str(e_comparison_t cmpt)
178 case CMP_INVALID: return "INVALID"; break;
179 case CMP_LESS: return "<"; break;
180 case CMP_LEQ: return "<="; break;
181 case CMP_GTR: return ">"; break;
182 case CMP_GEQ: return ">="; break;
183 case CMP_EQUAL: return "=="; break;
184 case CMP_NEQ: return "!="; break;
190 * \param[in] fp File to receive the output.
191 * \param[in] data Should point to a \c t_methoddata_compare.
194 _gmx_selelem_print_compare_info(FILE *fp, void *data)
196 t_methoddata_compare *d = (t_methoddata_compare *)data;
199 /* Print the left value */
200 if ((d->left.flags & CMP_SINGLEVAL) && !(d->left.flags & CMP_DYNAMICVAL))
202 if (d->left.flags & CMP_REALVAL)
204 fprintf(fp, "%f ", d->left.r[0]);
208 fprintf(fp, "%d ", d->left.i[0]);
211 /* Print the operator */
212 if (d->cmpt != CMP_INVALID)
214 fprintf(fp, "%s", comparison_type_str(d->cmpt));
218 fprintf(fp, "%s", d->cmpop);
220 /* Print the right value */
221 if ((d->right.flags & CMP_SINGLEVAL) && !(d->right.flags & CMP_DYNAMICVAL))
223 if (d->right.flags & CMP_REALVAL)
225 fprintf(fp, " %f", d->right.r[0]);
229 fprintf(fp, " %d", d->right.i[0]);
236 * \param[in] npar Not used (should be 5).
237 * \param[in,out] param Method parameters (should point to a copy of
238 * \ref smparams_compare).
239 * \returns Pointer to the allocated data (\c t_methoddata_compare).
241 * Allocates memory for a \c t_methoddata_compare structure.
244 init_data_compare(int npar, gmx_ana_selparam_t *param)
246 t_methoddata_compare *data;
249 param[2].val.u.s = &data->cmpop;
254 * Reverses a comparison operator.
256 * \param[in] type Comparison operator to reverse.
257 * \returns The correct comparison operator that equals \p type when the
258 * left and right sides are interchanged.
260 static e_comparison_t
261 reverse_comparison_type(e_comparison_t type)
265 case CMP_LESS: return CMP_GTR;
266 case CMP_LEQ: return CMP_GEQ;
267 case CMP_GTR: return CMP_LESS;
268 case CMP_GEQ: return CMP_LEQ;
275 * Initializes the value storage for comparison expression.
277 * \param[out] val Value structure to initialize.
278 * \param[in] param Parameters to use for initialization.
279 * \returns The number of values provided for the value, 0 on error.
282 init_comparison_value(t_compare_value *val, gmx_ana_selparam_t param[2])
287 if (param[0].flags & SPAR_SET)
289 val->flags |= (param[0].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
290 val->flags |= !(param[0].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL : 0;
292 val->i = param[0].val.u.i;
294 else if (param[1].flags & SPAR_SET)
296 val->flags |= (param[1].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
297 val->flags |= !(param[1].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL : 0;
298 val->flags |= CMP_REALVAL;
300 val->r = param[1].val.u.r;
312 * Converts an integer value to floating point.
314 * \param[in] n Number of values in the \p val->u array.
315 * \param[in,out] val Value to convert.
318 convert_int_real(int n, t_compare_value *val)
324 for (i = 0; i < n; ++i)
326 rv[i] = (real)val->i[i];
328 /* Free the previous value if one is present. */
331 val->flags |= CMP_REALVAL | CMP_ALLOCREAL;
335 * Converts a floating point value to integer.
337 * \param[in] n Number of values in the \p val->u array.
338 * \param[in,out] val Value to convert.
339 * \param[in] cmpt Comparison operator type.
340 * \param[in] bRight true if \p val appears on the right hand size of
342 * \returns 0 on success, EINVAL on error.
344 * The values are rounded such that the same comparison operator can be used.
347 convert_real_int(int n, t_compare_value *val, e_comparison_t cmpt, bool bRight)
354 cmpt = reverse_comparison_type(cmpt);
357 /* Round according to the comparison type */
358 for (i = 0; i < n; ++i)
364 iv[i] = (int)ceil(val->r[i]);
368 iv[i] = (int)floor(val->r[i]);
373 /* TODO: Implement, although it isn't typically very useful.
374 * Implementation is only a matter or proper initialization,
375 * the evaluation function can already handle this case with
376 * proper preparations. */
377 GMX_THROW(gmx::NotImplementedError("Equality comparison between dynamic integer and static real expressions not implemented"));
378 case CMP_INVALID: /* Should not be reached */
380 GMX_THROW(gmx::InternalError("Invalid comparison type"));
383 /* Free the previous value if one is present. */
386 val->flags &= ~CMP_REALVAL;
387 val->flags |= CMP_ALLOCINT;
391 * \param[in] top Not used.
392 * \param[in] npar Not used (should be 5).
393 * \param[in] param Method parameters (should point to \ref smparams_compare).
394 * \param[in] data Should point to a \c t_methoddata_compare.
395 * \returns 0 if the input data is valid, -1 on error.
398 init_compare(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
400 t_methoddata_compare *d = (t_methoddata_compare *)data;
403 /* Store the values */
404 n1 = init_comparison_value(&d->left, ¶m[0]);
405 n2 = init_comparison_value(&d->right, ¶m[3]);
406 if (n1 == 0 || n2 == 0)
408 GMX_THROW(gmx::InternalError("One of the values for comparison missing"));
410 /* Store the comparison type */
411 d->cmpt = comparison_type(d->cmpop);
412 if (d->cmpt == CMP_INVALID)
414 GMX_THROW(gmx::InternalError("Invalid comparison type"));
416 /* Convert the values to the same type */
417 /* TODO: Currently, there are no dynamic integer-valued selection methods,
418 * which means that only the branches with convert_int_real() will ever be
419 * taken. It should be considered whether it is necessary to support these
420 * other cases at all.
422 if ((d->left.flags & CMP_REALVAL) && !(d->right.flags & CMP_REALVAL))
424 if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
426 /* Nothing can be done */
428 else if (!(d->right.flags & CMP_DYNAMICVAL))
430 convert_int_real(n2, &d->right);
432 else /* d->left is static */
434 convert_real_int(n1, &d->left, d->cmpt, false);
437 else if (!(d->left.flags & CMP_REALVAL) && (d->right.flags & CMP_REALVAL))
439 if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
441 /* Reverse the sides to place the integer on the right */
443 d->left.r = d->right.r;
445 d->right.i = d->left.i;
447 flags = d->left.flags;
448 d->left.flags = d->right.flags;
449 d->right.flags = flags;
450 d->cmpt = reverse_comparison_type(d->cmpt);
452 else if (!(d->left.flags & CMP_DYNAMICVAL))
454 convert_int_real(n1, &d->left);
456 else /* d->right is static */
458 convert_real_int(n2, &d->right, d->cmpt, true);
464 * \param data Data to free (should point to a \c t_methoddata_compare).
466 * Frees the memory allocated for \c t_methoddata_compare.
469 free_data_compare(void *data)
471 t_methoddata_compare *d = (t_methoddata_compare *)data;
474 if (d->left.flags & CMP_ALLOCINT)
478 if (d->left.flags & CMP_ALLOCREAL)
482 if (d->right.flags & CMP_ALLOCINT)
486 if (d->right.flags & CMP_ALLOCREAL)
494 * \param[in] top Not used.
495 * \param[in] fr Not used.
496 * \param[in] pbc Not used.
497 * \param[in] g Evaluation index group.
498 * \param[out] out Output data structure (\p out->u.g is used).
499 * \param[in] data Should point to a \c t_methoddata_compare.
502 evaluate_compare_int(t_topology *top, t_trxframe *fr, t_pbc *pbc,
503 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
505 t_methoddata_compare *d = (t_methoddata_compare *)data;
510 for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
517 case CMP_INVALID: break;
518 case CMP_LESS: bAccept = a < b; break;
519 case CMP_LEQ: bAccept = a <= b; break;
520 case CMP_GTR: bAccept = a > b; break;
521 case CMP_GEQ: bAccept = a >= b; break;
522 case CMP_EQUAL: bAccept = a == b; break;
523 case CMP_NEQ: bAccept = a != b; break;
527 out->u.g->index[ig++] = g->index[i];
529 if (!(d->left.flags & CMP_SINGLEVAL))
533 if (!(d->right.flags & CMP_SINGLEVAL))
538 out->u.g->isize = ig;
542 * \param[in] top Not used.
543 * \param[in] fr Not used.
544 * \param[in] pbc Not used.
545 * \param[in] g Evaluation index group.
546 * \param[out] out Output data structure (\p out->u.g is used).
547 * \param[in] data Should point to a \c t_methoddata_compare.
550 evaluate_compare_real(t_topology *top, t_trxframe *fr, t_pbc *pbc,
551 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
553 t_methoddata_compare *d = (t_methoddata_compare *)data;
558 for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
561 b = (d->right.flags & CMP_REALVAL) ? d->right.r[i2] : d->right.i[i2];
565 case CMP_INVALID: break;
566 case CMP_LESS: bAccept = a < b; break;
567 case CMP_LEQ: bAccept = a <= b; break;
568 case CMP_GTR: bAccept = a > b; break;
569 case CMP_GEQ: bAccept = a >= b; break;
570 case CMP_EQUAL: bAccept = gmx_within_tol(a, b, GMX_REAL_EPS); break;
571 case CMP_NEQ: bAccept = !gmx_within_tol(a, b, GMX_REAL_EPS); break;
575 out->u.g->index[ig++] = g->index[i];
577 if (!(d->left.flags & CMP_SINGLEVAL))
581 if (!(d->right.flags & CMP_SINGLEVAL))
586 out->u.g->isize = ig;
590 * \param[in] top Not used.
591 * \param[in] fr Not used.
592 * \param[in] pbc Not used.
593 * \param[in] g Evaluation index group.
594 * \param[out] out Output data structure (\p out->u.g is used).
595 * \param[in] data Should point to a \c t_methoddata_compare.
598 evaluate_compare(t_topology *top, t_trxframe *fr, t_pbc *pbc,
599 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
601 t_methoddata_compare *d = (t_methoddata_compare *)data;
603 if (!((d->left.flags | d->right.flags) & CMP_REALVAL))
605 evaluate_compare_int(top, fr, pbc, g, out, data);
609 evaluate_compare_real(top, fr, pbc, g, out, data);