minor fixes in ditribution files
[gromacs/qmmm-gamess-us.git] / src / gmxlib / selection / sm_compare.c
blob54c932b457564062dcffb9b36552a91d5dfd9ecd
1 /*
3 * This source code is part of
5 * G R O M A C S
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
31 /*! \internal \file
32 * \brief
33 * Implementation of internal selection method for comparison expressions.
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
39 #include <maths.h>
40 #include <macros.h>
41 #include <smalloc.h>
42 #include <gmx_fatal.h>
44 #include <selmethod.h>
46 /** Defines the comparison operator for comparison expressions. */
47 typedef enum
49 CMP_INVALID, /**< Indicates an error */
50 CMP_LESS, /**< '<' */
51 CMP_LEQ, /**< '<=' */
52 CMP_GTR, /**< '>' */
53 CMP_GEQ, /**< '>=' */
54 CMP_EQUAL, /**< '==' */
55 CMP_NEQ, /**< '!=' */
56 } e_comparison_t;
58 /** The operand has a single value. */
59 #define CMP_SINGLEVAL 1
60 /** The operand value is dynamic. */
61 #define CMP_DYNAMICVAL 2
62 /** The value is real. */
63 #define CMP_REALVAL 4
64 /** The integer array is allocated. */
65 #define CMP_ALLOCINT 16
66 /** The real array is allocated. */
67 #define CMP_ALLOCREAL 32
69 /*! \internal \brief
70 * Data structure for comparison expression operand values.
72 typedef struct
74 /** Flags that describe the type of the operand. */
75 int flags;
76 /** (Array of) integer value(s). */
77 int *i;
78 /** (Array of) real value(s). */
79 real *r;
80 } t_compare_value;
82 /*! \internal \brief
83 * Data structure for comparison expression evaluation.
85 typedef struct
87 /** Comparison operator as a string. */
88 char *cmpop;
89 /** Comparison operator type. */
90 e_comparison_t cmpt;
91 /** Left value. */
92 t_compare_value left;
93 /** Right value. */
94 t_compare_value right;
95 } t_methoddata_compare;
97 /** Allocates data for comparison expression evaluation. */
98 static void *
99 init_data_compare(int npar, gmx_ana_selparam_t *param);
100 /** Initializes data for comparison expression evaluation. */
101 static int
102 init_compare(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
103 /** Frees the memory allocated for comparison expression evaluation. */
104 static void
105 free_data_compare(void *data);
106 /** Evaluates comparison expressions. */
107 static int
108 evaluate_compare(t_topology *top, t_trxframe *fr, t_pbc *pbc,
109 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
111 /** Parameters for comparison expression evaluation. */
112 static gmx_ana_selparam_t smparams_compare[] = {
113 {"int1", {INT_VALUE, -1, {NULL}}, NULL,
114 SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
115 {"real1", {REAL_VALUE, -1, {NULL}}, NULL,
116 SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
117 {"op", {STR_VALUE, 1, {NULL}}, NULL, 0},
118 {"int2", {INT_VALUE, -1, {NULL}}, NULL,
119 SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
120 {"real2", {REAL_VALUE, -1, {NULL}}, NULL,
121 SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
124 /** \internal Selection method data for comparison expression evaluation. */
125 gmx_ana_selmethod_t sm_compare = {
126 "cmp", GROUP_VALUE, SMETH_SINGLEVAL,
127 asize(smparams_compare), smparams_compare,
128 &init_data_compare,
129 NULL,
130 &init_compare,
131 NULL,
132 &free_data_compare,
133 NULL,
134 &evaluate_compare,
135 NULL,
136 {NULL, 0, NULL},
139 /*! \brief
140 * Returns a \c e_comparison_t value corresponding to an operator.
142 * \param[in] str String to process.
143 * \returns The comparison type corresponding to the first one or two
144 * characters of \p str.
146 * \p str can contain any number of characters; only the first two
147 * are used.
148 * If the beginning of \p str does not match any of the recognized types,
149 * \ref CMP_INVALID is returned.
151 static e_comparison_t
152 comparison_type(char *str)
154 switch (str[0])
156 case '<': return (str[1] == '=') ? CMP_LEQ : CMP_LESS;
157 case '>': return (str[1] == '=') ? CMP_GEQ : CMP_GTR;
158 case '=': return (str[1] == '=') ? CMP_EQUAL : CMP_INVALID;
159 case '!': return (str[1] == '=') ? CMP_NEQ : CMP_INVALID;
161 return CMP_INVALID;
164 /*! \brief
165 * Returns a string corresponding to a \c e_comparison_t value.
167 * \param[in] cmpt Comparison type to convert.
168 * \returns Pointer to a string that corresponds to \p cmpt.
170 * The return value points to a string constant and should not be \p free'd.
172 * The function returns NULL if \p cmpt is not one of the valid values.
174 static const char *
175 comparison_type_str(e_comparison_t cmpt)
177 switch (cmpt)
179 case CMP_INVALID: return "INVALID"; break;
180 case CMP_LESS: return "<"; break;
181 case CMP_LEQ: return "<="; break;
182 case CMP_GTR: return ">"; break;
183 case CMP_GEQ: return ">="; break;
184 case CMP_EQUAL: return "=="; break;
185 case CMP_NEQ: return "!="; break;
187 return NULL;
191 * \param[in] fp File to receive the output.
192 * \param[in] data Should point to a \c t_methoddata_compare.
194 void
195 _gmx_selelem_print_compare_info(FILE *fp, void *data)
197 t_methoddata_compare *d = (t_methoddata_compare *)data;
199 fprintf(fp, " \"");
200 /* Print the left value */
201 if ((d->left.flags & CMP_SINGLEVAL) && !(d->left.flags & CMP_DYNAMICVAL))
203 if (d->left.flags & CMP_REALVAL)
205 fprintf(fp, "%f ", d->left.r[0]);
207 else
209 fprintf(fp, "%d ", d->left.i[0]);
212 /* Print the operator */
213 if (d->cmpt != CMP_INVALID)
215 fprintf(fp, "%s", comparison_type_str(d->cmpt));
217 else
219 fprintf(fp, "%s", d->cmpop);
221 /* Print the right value */
222 if ((d->right.flags & CMP_SINGLEVAL) && !(d->right.flags & CMP_DYNAMICVAL))
224 if (d->right.flags & CMP_REALVAL)
226 fprintf(fp, " %f", d->right.r[0]);
228 else
230 fprintf(fp, " %d", d->right.i[0]);
233 fprintf(fp, "\"");
237 * \param[in] npar Not used (should be 5).
238 * \param[in,out] param Method parameters (should point to a copy of
239 * \ref smparams_compare).
240 * \returns Pointer to the allocated data (\c t_methoddata_compare).
242 * Allocates memory for a \c t_methoddata_compare structure.
244 static void *
245 init_data_compare(int npar, gmx_ana_selparam_t *param)
247 t_methoddata_compare *data;
249 snew(data, 1);
250 param[2].val.u.s = &data->cmpop;
251 return data;
254 /* \brief
255 * Reverses a comparison operator.
257 * \param[in] type Comparison operator to reverse.
258 * \returns The correct comparison operator that equals \p type when the
259 * left and right sides are interchanged.
261 static e_comparison_t
262 reverse_comparison_type(e_comparison_t type)
264 switch (type)
266 case CMP_LESS: return CMP_GTR;
267 case CMP_LEQ: return CMP_GEQ;
268 case CMP_GTR: return CMP_LESS;
269 case CMP_GEQ: return CMP_LEQ;
270 default: break;
272 return type;
275 /*! \brief
276 * Initializes the value storage for comparison expression.
278 * \param[out] val Value structure to initialize.
279 * \param[in] param Parameters to use for initialization.
280 * \returns The number of values provided for the value, 0 on error.
282 static int
283 init_comparison_value(t_compare_value *val, gmx_ana_selparam_t param[2])
285 int n;
287 val->flags = 0;
288 if (param[0].flags & SPAR_SET)
290 val->flags |= (param[0].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
291 val->flags |= !(param[0].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL : 0;
292 n = param[0].val.nr;
293 val->i = param[0].val.u.i;
295 else if (param[1].flags & SPAR_SET)
297 val->flags |= (param[1].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
298 val->flags |= !(param[1].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL : 0;
299 val->flags |= CMP_REALVAL;
300 n = param[1].val.nr;
301 val->r = param[1].val.u.r;
303 else
305 n = 0;
306 val->i = NULL;
307 val->r = NULL;
309 return n;
312 /* \brief
313 * Converts an integer value to floating point.
315 * \param[in] n Number of values in the \p val->u array.
316 * \param[in,out] val Value to convert.
318 static void
319 convert_int_real(int n, t_compare_value *val)
321 int i;
322 real *rv;
324 snew(rv, n);
325 for (i = 0; i < n; ++i)
327 rv[i] = (real)val->i[i];
329 /* Free the previous value if one is present. */
330 sfree(val->r);
331 val->r = rv;
332 val->flags |= CMP_REALVAL | CMP_ALLOCREAL;
335 /* \brief
336 * Converts a floating point value to integer.
338 * \param[in] n Number of values in the \p val->u array.
339 * \param[in,out] val Value to convert.
340 * \param[in] cmpt Comparison operator type.
341 * \param[in] bRight TRUE if \p val appears on the right hand size of
342 * \p cmpt.
343 * \returns 0 on success, EINVAL on error.
345 * The values are rounded such that the same comparison operator can be used.
347 static int
348 convert_real_int(int n, t_compare_value *val, e_comparison_t cmpt, bool bRight)
350 int i;
351 int *iv;
353 if (!bRight)
355 cmpt = reverse_comparison_type(cmpt);
357 snew(iv, n);
358 /* Round according to the comparison type */
359 for (i = 0; i < n; ++i)
361 switch (cmpt)
363 case CMP_LESS:
364 case CMP_GEQ:
365 iv[i] = (int)ceil(val->r[i]);
366 break;
367 case CMP_GTR:
368 case CMP_LEQ:
369 iv[i] = (int)floor(val->r[i]);
370 break;
371 case CMP_EQUAL:
372 case CMP_NEQ:
373 fprintf(stderr, "comparing equality an integer expression and a real value\n");
374 sfree(iv);
375 return EINVAL;
376 case CMP_INVALID: /* Should not be reached */
377 gmx_bug("internal error");
378 sfree(iv);
379 return EINVAL;
382 /* Free the previous value if one is present. */
383 sfree(val->i);
384 val->i = iv;
385 val->flags &= ~CMP_REALVAL;
386 val->flags |= CMP_ALLOCINT;
387 return 0;
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.
397 static int
398 init_compare(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
400 t_methoddata_compare *d = (t_methoddata_compare *)data;
401 int n1, n2;
403 /* Store the values */
404 n1 = init_comparison_value(&d->left, &param[0]);
405 n2 = init_comparison_value(&d->right, &param[3]);
406 if (n1 == 0 || n2 == 0)
408 gmx_bug("one of the values for comparison missing");
409 return -1;
411 /* Store the comparison type */
412 d->cmpt = comparison_type(d->cmpop);
413 if (d->cmpt == CMP_INVALID)
415 gmx_bug("invalid comparison type");
416 return -1;
418 /* Convert the values to the same type */
419 if ((d->left.flags & CMP_REALVAL) && !(d->right.flags & CMP_REALVAL))
421 if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
423 /* Nothing can be done */
425 else if (!(d->right.flags & CMP_DYNAMICVAL))
427 convert_int_real(n2, &d->right);
429 else /* d->left is static */
431 if (convert_real_int(n1, &d->left, d->cmpt, FALSE))
433 return -1;
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 */
442 int flags;
443 d->left.r = d->right.r;
444 d->right.r = NULL;
445 d->right.i = d->left.i;
446 d->left.i = NULL;
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 if (convert_real_int(n2, &d->right, d->cmpt, TRUE))
460 return -1;
464 return 0;
468 * \param data Data to free (should point to a \c t_methoddata_compare).
470 * Frees the memory allocated for \c t_methoddata_compare.
472 static void
473 free_data_compare(void *data)
475 t_methoddata_compare *d = (t_methoddata_compare *)data;
477 sfree(d->cmpop);
478 if (d->left.flags & CMP_ALLOCINT)
480 sfree(d->left.i);
482 if (d->left.flags & CMP_ALLOCREAL)
484 sfree(d->left.r);
486 if (d->right.flags & CMP_ALLOCINT)
488 sfree(d->right.i);
490 if (d->right.flags & CMP_ALLOCREAL)
492 sfree(d->right.r);
497 * \param[in] top Not used.
498 * \param[in] fr Not used.
499 * \param[in] pbc Not used.
500 * \param[in] g Evaluation index group.
501 * \param[out] out Output data structure (\p out->u.g is used).
502 * \param[in] data Should point to a \c t_methoddata_compare.
503 * \returns 0 for success.
505 static int
506 evaluate_compare_int(t_topology *top, t_trxframe *fr, t_pbc *pbc,
507 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
509 t_methoddata_compare *d = (t_methoddata_compare *)data;
510 int i, i1, i2, ig;
511 int a, b;
512 bool bAccept;
514 for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
516 a = d->left.i[i1];
517 b = d->right.i[i2];
518 bAccept = FALSE;
519 switch (d->cmpt)
521 case CMP_INVALID: break;
522 case CMP_LESS: bAccept = a < b; break;
523 case CMP_LEQ: bAccept = a <= b; break;
524 case CMP_GTR: bAccept = a > b; break;
525 case CMP_GEQ: bAccept = a >= b; break;
526 case CMP_EQUAL: bAccept = a == b; break;
527 case CMP_NEQ: bAccept = a != b; break;
529 if (bAccept)
531 out->u.g->index[ig++] = g->index[i];
533 if (!(d->left.flags & CMP_SINGLEVAL))
535 ++i1;
537 if (!(d->right.flags & CMP_SINGLEVAL))
539 ++i2;
542 out->u.g->isize = ig;
543 return 0;
547 * \param[in] top Not used.
548 * \param[in] fr Not used.
549 * \param[in] pbc Not used.
550 * \param[in] g Evaluation index group.
551 * \param[out] out Output data structure (\p out->u.g is used).
552 * \param[in] data Should point to a \c t_methoddata_compare.
553 * \returns 0 for success.
555 static int
556 evaluate_compare_real(t_topology *top, t_trxframe *fr, t_pbc *pbc,
557 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
559 t_methoddata_compare *d = (t_methoddata_compare *)data;
560 int i, i1, i2, ig;
561 real a, b;
562 bool bAccept;
564 for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
566 a = d->left.r[i1];
567 b = (d->right.flags & CMP_REALVAL) ? d->right.r[i2] : d->right.i[i2];
568 bAccept = FALSE;
569 switch (d->cmpt)
571 case CMP_INVALID: break;
572 case CMP_LESS: bAccept = a < b; break;
573 case CMP_LEQ: bAccept = a <= b; break;
574 case CMP_GTR: bAccept = a > b; break;
575 case CMP_GEQ: bAccept = a >= b; break;
576 case CMP_EQUAL: bAccept = gmx_within_tol(a, b, GMX_REAL_EPS); break;
577 case CMP_NEQ: bAccept = !gmx_within_tol(a, b, GMX_REAL_EPS); break;
579 if (bAccept)
581 out->u.g->index[ig++] = g->index[i];
583 if (!(d->left.flags & CMP_SINGLEVAL))
585 ++i1;
587 if (!(d->right.flags & CMP_SINGLEVAL))
589 ++i2;
592 out->u.g->isize = ig;
593 return 0;
597 * \param[in] top Not used.
598 * \param[in] fr Not used.
599 * \param[in] pbc Not used.
600 * \param[in] g Evaluation index group.
601 * \param[out] out Output data structure (\p out->u.g is used).
602 * \param[in] data Should point to a \c t_methoddata_compare.
603 * \returns 0 for success.
605 static int
606 evaluate_compare(t_topology *top, t_trxframe *fr, t_pbc *pbc,
607 gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
609 t_methoddata_compare *d = (t_methoddata_compare *)data;
611 if (!((d->left.flags | d->right.flags) & CMP_REALVAL))
613 return evaluate_compare_int(top, fr, pbc, g, out, data);
615 else
617 return evaluate_compare_real(top, fr, pbc, g, out, data);
619 return 0;