2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2009,2010,2011,2012,2013, by the GROMACS development team, led by
5 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
6 * others, as listed in the AUTHORS file in the top-level source
7 * directory and at http://www.gromacs.org.
9 * GROMACS is free software; you can redistribute it and/or
10 * modify it under the terms of the GNU Lesser General Public License
11 * as published by the Free Software Foundation; either version 2.1
12 * of the License, or (at your option) any later version.
14 * GROMACS is distributed in the hope that it will be useful,
15 * but WITHOUT ANY WARRANTY; without even the implied warranty of
16 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
17 * Lesser General Public License for more details.
19 * You should have received a copy of the GNU Lesser General Public
20 * License along with GROMACS; if not, see
21 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
22 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
24 * If you want to redistribute modifications to GROMACS, please
25 * consider that scientific software is very special. Version
26 * control is crucial - bugs must be traceable. We will be happy to
27 * consider code for inclusion in the official distribution, but
28 * derived work must not be called official GROMACS. Details are found
29 * in the README & COPYING files - if they are missing, get the
30 * official version at http://www.gromacs.org.
32 * To help us fund GROMACS development, we humbly ask that you cite
33 * the research papers on the package. Check out http://www.gromacs.org.
37 * Implements internal selection method for comparison expressions.
39 * \author Teemu Murtola <teemu.murtola@gmail.com>
40 * \ingroup module_selection
42 #include "gromacs/legacyheaders/maths.h"
43 #include "gromacs/legacyheaders/macros.h"
44 #include "gromacs/legacyheaders/smalloc.h"
46 #include "gromacs/selection/selmethod.h"
47 #include "gromacs/utility/common.h"
48 #include "gromacs/utility/exceptions.h"
50 /** Defines the comparison operator for comparison expressions. */
53 CMP_INVALID, /**< Indicates an error */
58 CMP_EQUAL, /**< '==' */
62 /** The operand has a single value. */
63 #define CMP_SINGLEVAL 1
64 /** The operand value is dynamic. */
65 #define CMP_DYNAMICVAL 2
66 /** The value is real. */
68 /** The integer array is allocated. */
69 #define CMP_ALLOCINT 16
70 /** The real array is allocated. */
71 #define CMP_ALLOCREAL 32
74 * Data structure for comparison expression operand values.
78 /** Flags that describe the type of the operand. */
80 /** (Array of) integer value(s). */
82 /** (Array of) real value(s). */
87 * Data structure for comparison expression evaluation.
91 /** Comparison operator as a string. */
93 /** Comparison operator type. */
98 t_compare_value right;
99 } t_methoddata_compare;
102 * Allocates data for comparison expression evaluation.
104 * \param[in] npar Not used (should be 5).
105 * \param[in,out] param Method parameters (should point to a copy of
106 * \ref smparams_compare).
107 * \returns Pointer to the allocated data (\c t_methoddata_compare).
109 * Allocates memory for a \c t_methoddata_compare structure.
112 init_data_compare(int npar, gmx_ana_selparam_t *param);
114 * Initializes data for comparison expression evaluation.
116 * \param[in] top Not used.
117 * \param[in] npar Not used (should be 5).
118 * \param[in] param Method parameters (should point to \ref smparams_compare).
119 * \param[in] data Should point to a \c t_methoddata_compare.
120 * \returns 0 if the input data is valid, -1 on error.
123 init_compare(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
124 /** Frees the memory allocated for comparison expression evaluation. */
126 free_data_compare(void *data);
127 /** Evaluates comparison expressions. */
129 evaluate_compare(t_topology *top, t_trxframe *fr, t_pbc *pbc,
130 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
132 /** Parameters for comparison expression evaluation. */
133 static gmx_ana_selparam_t smparams_compare[] = {
134 {"int1", {INT_VALUE, -1, {NULL}}, NULL,
135 SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
136 {"real1", {REAL_VALUE, -1, {NULL}}, NULL,
137 SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
138 {"op", {STR_VALUE, 1, {NULL}}, NULL, 0},
139 {"int2", {INT_VALUE, -1, {NULL}}, NULL,
140 SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
141 {"real2", {REAL_VALUE, -1, {NULL}}, NULL,
142 SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
145 /** \internal Selection method data for comparison expression evaluation. */
146 gmx_ana_selmethod_t sm_compare = {
147 "cmp", GROUP_VALUE, SMETH_SINGLEVAL,
148 asize(smparams_compare), smparams_compare,
161 * Returns a \c e_comparison_t value corresponding to an operator.
163 * \param[in] str String to process.
164 * \returns The comparison type corresponding to the first one or two
165 * characters of \p str.
167 * \p str can contain any number of characters; only the first two
169 * If the beginning of \p str does not match any of the recognized types,
170 * \ref CMP_INVALID is returned.
172 static e_comparison_t
173 comparison_type(char *str)
177 case '<': return (str[1] == '=') ? CMP_LEQ : CMP_LESS;
178 case '>': return (str[1] == '=') ? CMP_GEQ : CMP_GTR;
179 case '=': return (str[1] == '=') ? CMP_EQUAL : CMP_INVALID;
180 case '!': return (str[1] == '=') ? CMP_NEQ : CMP_INVALID;
186 * Returns a string corresponding to a \c e_comparison_t value.
188 * \param[in] cmpt Comparison type to convert.
189 * \returns Pointer to a string that corresponds to \p cmpt.
191 * The return value points to a string constant and should not be \p free'd.
193 * The function returns NULL if \p cmpt is not one of the valid values.
196 comparison_type_str(e_comparison_t cmpt)
200 case CMP_INVALID: return "INVALID"; break;
201 case CMP_LESS: return "<"; break;
202 case CMP_LEQ: return "<="; break;
203 case CMP_GTR: return ">"; break;
204 case CMP_GEQ: return ">="; break;
205 case CMP_EQUAL: return "=="; break;
206 case CMP_NEQ: return "!="; break;
212 * \param[in] fp File to receive the output.
213 * \param[in] data Should point to a \c t_methoddata_compare.
216 _gmx_selelem_print_compare_info(FILE *fp, void *data)
218 t_methoddata_compare *d = (t_methoddata_compare *)data;
221 /* Print the left value */
222 if ((d->left.flags & CMP_SINGLEVAL) && !(d->left.flags & CMP_DYNAMICVAL))
224 if (d->left.flags & CMP_REALVAL)
226 fprintf(fp, "%f ", d->left.r[0]);
230 fprintf(fp, "%d ", d->left.i[0]);
233 /* Print the operator */
234 if (d->cmpt != CMP_INVALID)
236 fprintf(fp, "%s", comparison_type_str(d->cmpt));
240 fprintf(fp, "%s", d->cmpop);
242 /* Print the right value */
243 if ((d->right.flags & CMP_SINGLEVAL) && !(d->right.flags & CMP_DYNAMICVAL))
245 if (d->right.flags & CMP_REALVAL)
247 fprintf(fp, " %f", d->right.r[0]);
251 fprintf(fp, " %d", d->right.i[0]);
258 init_data_compare(int /* npar */, gmx_ana_selparam_t *param)
260 t_methoddata_compare *data;
263 param[2].val.u.s = &data->cmpop;
268 * Reverses a comparison operator.
270 * \param[in] type Comparison operator to reverse.
271 * \returns The correct comparison operator that equals \p type when the
272 * left and right sides are interchanged.
274 static e_comparison_t
275 reverse_comparison_type(e_comparison_t type)
279 case CMP_LESS: return CMP_GTR;
280 case CMP_LEQ: return CMP_GEQ;
281 case CMP_GTR: return CMP_LESS;
282 case CMP_GEQ: return CMP_LEQ;
289 * Initializes the value storage for comparison expression.
291 * \param[out] val Value structure to initialize.
292 * \param[in] param Parameters to use for initialization.
293 * \returns The number of values provided for the value, 0 on error.
296 init_comparison_value(t_compare_value *val, gmx_ana_selparam_t param[2])
301 if (param[0].flags & SPAR_SET)
303 val->flags |= (param[0].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
304 val->flags |= !(param[0].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL : 0;
306 val->i = param[0].val.u.i;
308 else if (param[1].flags & SPAR_SET)
310 val->flags |= (param[1].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
311 val->flags |= !(param[1].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL : 0;
312 val->flags |= CMP_REALVAL;
314 val->r = param[1].val.u.r;
326 * Converts an integer value to floating point.
328 * \param[in] n Number of values in the \p val->u array.
329 * \param[in,out] val Value to convert.
332 convert_int_real(int n, t_compare_value *val)
338 for (i = 0; i < n; ++i)
340 rv[i] = (real)val->i[i];
342 /* Free the previous value if one is present. */
345 val->flags |= CMP_REALVAL | CMP_ALLOCREAL;
349 * Converts a floating point value to integer.
351 * \param[in] n Number of values in the \p val->u array.
352 * \param[in,out] val Value to convert.
353 * \param[in] cmpt Comparison operator type.
354 * \param[in] bRight true if \p val appears on the right hand size of
356 * \returns 0 on success, EINVAL on error.
358 * The values are rounded such that the same comparison operator can be used.
361 convert_real_int(int n, t_compare_value *val, e_comparison_t cmpt, bool bRight)
368 cmpt = reverse_comparison_type(cmpt);
371 /* Round according to the comparison type */
372 for (i = 0; i < n; ++i)
378 iv[i] = (int)ceil(val->r[i]);
382 iv[i] = (int)floor(val->r[i]);
387 /* TODO: Implement, although it isn't typically very useful.
388 * Implementation is only a matter or proper initialization,
389 * the evaluation function can already handle this case with
390 * proper preparations. */
391 GMX_THROW(gmx::NotImplementedError("Equality comparison between dynamic integer and static real expressions not implemented"));
392 case CMP_INVALID: /* Should not be reached */
394 GMX_THROW(gmx::InternalError("Invalid comparison type"));
397 /* Free the previous value if one is present. */
400 val->flags &= ~CMP_REALVAL;
401 val->flags |= CMP_ALLOCINT;
405 init_compare(t_topology * /* top */, int /* npar */, gmx_ana_selparam_t *param, void *data)
407 t_methoddata_compare *d = (t_methoddata_compare *)data;
410 /* Store the values */
411 n1 = init_comparison_value(&d->left, ¶m[0]);
412 n2 = init_comparison_value(&d->right, ¶m[3]);
413 /* Store the comparison type */
414 d->cmpt = comparison_type(d->cmpop);
415 if (d->cmpt == CMP_INVALID)
417 GMX_THROW(gmx::InternalError("Invalid comparison type"));
419 /* Convert the values to the same type */
420 /* TODO: Currently, there are no dynamic integer-valued selection methods,
421 * which means that only the branches with convert_int_real() will ever be
422 * taken. It should be considered whether it is necessary to support these
423 * other cases at all.
425 if ((d->left.flags & CMP_REALVAL) && !(d->right.flags & CMP_REALVAL))
427 if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
429 /* Nothing can be done */
431 else if (!(d->right.flags & CMP_DYNAMICVAL))
433 convert_int_real(n2, &d->right);
435 else /* d->left is static */
437 convert_real_int(n1, &d->left, d->cmpt, false);
440 else if (!(d->left.flags & CMP_REALVAL) && (d->right.flags & CMP_REALVAL))
442 if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
444 /* Reverse the sides to place the integer on the right */
446 d->left.r = d->right.r;
448 d->right.i = d->left.i;
450 flags = d->left.flags;
451 d->left.flags = d->right.flags;
452 d->right.flags = flags;
453 d->cmpt = reverse_comparison_type(d->cmpt);
455 else if (!(d->left.flags & CMP_DYNAMICVAL))
457 convert_int_real(n1, &d->left);
459 else /* d->right is static */
461 convert_real_int(n2, &d->right, d->cmpt, true);
467 * \param data Data to free (should point to a \c t_methoddata_compare).
469 * Frees the memory allocated for \c t_methoddata_compare.
472 free_data_compare(void *data)
474 t_methoddata_compare *d = (t_methoddata_compare *)data;
477 if (d->left.flags & CMP_ALLOCINT)
481 if (d->left.flags & CMP_ALLOCREAL)
485 if (d->right.flags & CMP_ALLOCINT)
489 if (d->right.flags & CMP_ALLOCREAL)
497 * Implementation for evaluate_compare() for integer values.
499 * \param[in] top Not used.
500 * \param[in] fr Not used.
501 * \param[in] pbc Not used.
502 * \param[in] g Evaluation index group.
503 * \param[out] out Output data structure (\p out->u.g is used).
504 * \param[in] data Should point to a \c t_methoddata_compare.
507 evaluate_compare_int(t_topology *top, t_trxframe *fr, t_pbc *pbc,
508 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
510 t_methoddata_compare *d = (t_methoddata_compare *)data;
515 GMX_UNUSED_VALUE(top);
516 GMX_UNUSED_VALUE(fr);
517 GMX_UNUSED_VALUE(pbc);
518 for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
525 case CMP_INVALID: break;
526 case CMP_LESS: bAccept = a < b; break;
527 case CMP_LEQ: bAccept = a <= b; break;
528 case CMP_GTR: bAccept = a > b; break;
529 case CMP_GEQ: bAccept = a >= b; break;
530 case CMP_EQUAL: bAccept = a == b; break;
531 case CMP_NEQ: bAccept = a != b; break;
535 out->u.g->index[ig++] = g->index[i];
537 if (!(d->left.flags & CMP_SINGLEVAL))
541 if (!(d->right.flags & CMP_SINGLEVAL))
546 out->u.g->isize = ig;
550 * Implementation for evaluate_compare() if either value is non-integer.
552 * \param[in] top Not used.
553 * \param[in] fr Not used.
554 * \param[in] pbc Not used.
555 * \param[in] g Evaluation index group.
556 * \param[out] out Output data structure (\p out->u.g is used).
557 * \param[in] data Should point to a \c t_methoddata_compare.
559 * Left value is assumed to be real-valued; right value can be either.
560 * This is ensured by the initialization method.
563 evaluate_compare_real(t_topology *top, t_trxframe *fr, t_pbc *pbc,
564 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
566 t_methoddata_compare *d = (t_methoddata_compare *)data;
571 GMX_UNUSED_VALUE(top);
572 GMX_UNUSED_VALUE(fr);
573 GMX_UNUSED_VALUE(pbc);
574 for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
577 b = (d->right.flags & CMP_REALVAL) ? d->right.r[i2] : d->right.i[i2];
581 case CMP_INVALID: break;
582 case CMP_LESS: bAccept = a < b; break;
583 case CMP_LEQ: bAccept = a <= b; break;
584 case CMP_GTR: bAccept = a > b; break;
585 case CMP_GEQ: bAccept = a >= b; break;
586 case CMP_EQUAL: bAccept = gmx_within_tol(a, b, GMX_REAL_EPS); break;
587 case CMP_NEQ: bAccept = !gmx_within_tol(a, b, GMX_REAL_EPS); break;
591 out->u.g->index[ig++] = g->index[i];
593 if (!(d->left.flags & CMP_SINGLEVAL))
597 if (!(d->right.flags & CMP_SINGLEVAL))
602 out->u.g->isize = ig;
606 * \param[in] top Not used.
607 * \param[in] fr Not used.
608 * \param[in] pbc Not used.
609 * \param[in] g Evaluation index group.
610 * \param[out] out Output data structure (\p out->u.g is used).
611 * \param[in] data Should point to a \c t_methoddata_compare.
614 evaluate_compare(t_topology *top, t_trxframe *fr, t_pbc *pbc,
615 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
617 t_methoddata_compare *d = (t_methoddata_compare *)data;
619 if (!((d->left.flags | d->right.flags) & CMP_REALVAL))
621 evaluate_compare_int(top, fr, pbc, g, out, data);
625 evaluate_compare_real(top, fr, pbc, g, out, data);