src/gromacs/selection/sm_compare.cpp

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   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.
  13
  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.
  18  *
  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.
  25  *
  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.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \internal \file
  32  * \brief
  33  * Implements internal selection method for comparison expressions.
  34  *
  35  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  36  * \ingroup module_selection
  37  */
  38 #ifdef HAVE_CONFIG_H
  39 #include <config.h>
  40 #endif
  41
  42 #include "maths.h"
  43 #include "macros.h"
  44 #include "smalloc.h"
  45
  46 #include "gromacs/selection/selmethod.h"
  47 #include "gromacs/utility/exceptions.h"
  48
  49 /** Defines the comparison operator for comparison expressions. */
  50 typedef enum
  51 {
  52     CMP_INVALID,        /**< Indicates an error */
  53     CMP_LESS,           /**< '<' */
  54     CMP_LEQ,            /**< '<=' */
  55     CMP_GTR,            /**< '>' */
  56     CMP_GEQ,            /**< '>=' */
  57     CMP_EQUAL,          /**< '==' */
  58     CMP_NEQ             /**< '!=' */
  59 } e_comparison_t;
  60
  61 /** The operand has a single value. */
  62 #define CMP_SINGLEVAL  1
  63 /** The operand value is dynamic. */
  64 #define CMP_DYNAMICVAL 2
  65 /** The value is real. */
  66 #define CMP_REALVAL    4
  67 /** The integer array is allocated. */
  68 #define CMP_ALLOCINT   16
  69 /** The real array is allocated. */
  70 #define CMP_ALLOCREAL  32
  71
  72 /*! \internal \brief
  73  * Data structure for comparison expression operand values.
  74  */
  75 typedef struct
  76 {
  77     /** Flags that describe the type of the operand. */
  78     int             flags;
  79     /** (Array of) integer value(s). */
  80     int        *i;
  81     /** (Array of) real value(s). */
  82     real       *r;
  83 } t_compare_value;
  84
  85 /*! \internal \brief
  86  * Data structure for comparison expression evaluation.
  87  */
  88 typedef struct
  89 {
  90     /** Comparison operator as a string. */
  91     char            *cmpop;
  92     /** Comparison operator type. */
  93     e_comparison_t   cmpt;
  94     /** Left value. */
  95     t_compare_value  left;
  96     /** Right value. */
  97     t_compare_value  right;
  98 } t_methoddata_compare;
  99
 100 /** Allocates data for comparison expression evaluation. */
 101 static void *
 102 init_data_compare(int npar, gmx_ana_selparam_t *param);
 103 /** Initializes data for comparison expression evaluation. */
 104 static void
 105 init_compare(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
 106 /** Frees the memory allocated for comparison expression evaluation. */
 107 static void
 108 free_data_compare(void *data);
 109 /** Evaluates comparison expressions. */
 110 static void
 111 evaluate_compare(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 112                  gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
 113
 114 /** Parameters for comparison expression evaluation. */
 115 static gmx_ana_selparam_t smparams_compare[] = {
 116     {"int1",  {INT_VALUE,  -1, {NULL}}, NULL,
 117      SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
 118     {"real1", {REAL_VALUE, -1, {NULL}}, NULL,
 119      SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
 120     {"op",    {STR_VALUE,   1, {NULL}}, NULL, 0},
 121     {"int2",  {INT_VALUE,  -1, {NULL}}, NULL,
 122      SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
 123     {"real2", {REAL_VALUE, -1, {NULL}}, NULL,
 124      SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL},
 125 };
 126
 127 /** \internal Selection method data for comparison expression evaluation. */
 128 gmx_ana_selmethod_t sm_compare = {
 129     "cmp", GROUP_VALUE, SMETH_SINGLEVAL,
 130     asize(smparams_compare), smparams_compare,
 131     &init_data_compare,
 132     NULL,
 133     &init_compare,
 134     NULL,
 135     &free_data_compare,
 136     NULL,
 137     &evaluate_compare,
 138     NULL,
 139     {NULL, 0, NULL},
 140 };
 141
 142 /*! \brief
 143  * Returns a \c e_comparison_t value corresponding to an operator.
 144  *
 145  * \param[in] str  String to process.
 146  * \returns   The comparison type corresponding to the first one or two
 147  *   characters of \p str.
 148  *
 149  * \p str can contain any number of characters; only the first two
 150  * are used.
 151  * If the beginning of \p str does not match any of the recognized types,
 152  * \ref CMP_INVALID is returned.
 153  */
 154 static e_comparison_t
 155 comparison_type(char *str)
 156 {
 157     switch (str[0])
 158     {
 159         case '<': return (str[1] == '=') ? CMP_LEQ   : CMP_LESS;
 160         case '>': return (str[1] == '=') ? CMP_GEQ   : CMP_GTR;
 161         case '=': return (str[1] == '=') ? CMP_EQUAL : CMP_INVALID;
 162         case '!': return (str[1] == '=') ? CMP_NEQ   : CMP_INVALID;
 163     }
 164     return CMP_INVALID;
 165 }
 166
 167 /*! \brief
 168  * Returns a string corresponding to a \c e_comparison_t value.
 169  *
 170  * \param[in] cmpt  Comparison type to convert.
 171  * \returns   Pointer to a string that corresponds to \p cmpt.
 172  *
 173  * The return value points to a string constant and should not be \p free'd.
 174  *
 175  * The function returns NULL if \p cmpt is not one of the valid values.
 176  */
 177 static const char *
 178 comparison_type_str(e_comparison_t cmpt)
 179 {
 180     switch (cmpt)
 181     {
 182         case CMP_INVALID: return "INVALID"; break;
 183         case CMP_LESS:    return "<";  break;
 184         case CMP_LEQ:     return "<="; break;
 185         case CMP_GTR:     return ">";  break;
 186         case CMP_GEQ:     return ">="; break;
 187         case CMP_EQUAL:   return "=="; break;
 188         case CMP_NEQ:     return "!="; break;
 189     }
 190     return NULL;
 191 }
 192
 193 /*!
 194  * \param[in] fp    File to receive the output.
 195  * \param[in] data  Should point to a \c t_methoddata_compare.
 196  */
 197 void
 198 _gmx_selelem_print_compare_info(FILE *fp, void *data)
 199 {
 200     t_methoddata_compare *d = (t_methoddata_compare *)data;
 201
 202     fprintf(fp, " \"");
 203     /* Print the left value */
 204     if ((d->left.flags & CMP_SINGLEVAL) && !(d->left.flags & CMP_DYNAMICVAL))
 205     {
 206         if (d->left.flags & CMP_REALVAL)
 207         {
 208             fprintf(fp, "%f ", d->left.r[0]);
 209         }
 210         else
 211         {
 212             fprintf(fp, "%d ", d->left.i[0]);
 213         }
 214     }
 215     /* Print the operator */
 216     if (d->cmpt != CMP_INVALID)
 217     {
 218         fprintf(fp, "%s", comparison_type_str(d->cmpt));
 219     }
 220     else
 221     {
 222         fprintf(fp, "%s", d->cmpop);
 223     }
 224     /* Print the right value */
 225     if ((d->right.flags & CMP_SINGLEVAL) && !(d->right.flags & CMP_DYNAMICVAL))
 226     {
 227         if (d->right.flags & CMP_REALVAL)
 228         {
 229             fprintf(fp, " %f", d->right.r[0]);
 230         }
 231         else
 232         {
 233             fprintf(fp, " %d", d->right.i[0]);
 234         }
 235     }
 236     fprintf(fp, "\"");
 237 }
 238
 239 /*!
 240  * \param[in]     npar  Not used (should be 5).
 241  * \param[in,out] param Method parameters (should point to a copy of
 242  *   \ref smparams_compare).
 243  * \returns       Pointer to the allocated data (\c t_methoddata_compare).
 244  *
 245  * Allocates memory for a \c t_methoddata_compare structure.
 246  */
 247 static void *
 248 init_data_compare(int npar, gmx_ana_selparam_t *param)
 249 {
 250     t_methoddata_compare *data;
 251
 252     snew(data, 1);
 253     param[2].val.u.s = &data->cmpop;
 254     return data;
 255 }
 256
 257 /* \brief
 258  * Reverses a comparison operator.
 259  *
 260  * \param[in] type  Comparison operator to reverse.
 261  * \returns   The correct comparison operator that equals \p type when the
 262  *   left and right sides are interchanged.
 263  */
 264 static e_comparison_t
 265 reverse_comparison_type(e_comparison_t type)
 266 {
 267     switch (type)
 268     {
 269         case CMP_LESS: return CMP_GTR;
 270         case CMP_LEQ:  return CMP_GEQ;
 271         case CMP_GTR:  return CMP_LESS;
 272         case CMP_GEQ:  return CMP_LEQ;
 273         default:       break;
 274     }
 275     return type;
 276 }
 277
 278 /*! \brief
 279  * Initializes the value storage for comparison expression.
 280  *
 281  * \param[out] val   Value structure to initialize.
 282  * \param[in]  param Parameters to use for initialization.
 283  * \returns    The number of values provided for the value, 0 on error.
 284  */
 285 static int
 286 init_comparison_value(t_compare_value *val, gmx_ana_selparam_t param[2])
 287 {
 288     int  n;
 289
 290     val->flags = 0;
 291     if (param[0].flags & SPAR_SET)
 292     {
 293         val->flags |=  (param[0].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
 294         val->flags |= !(param[0].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL  : 0;
 295         n           = param[0].val.nr;
 296         val->i      = param[0].val.u.i;
 297     }
 298     else if (param[1].flags & SPAR_SET)
 299     {
 300         val->flags |=  (param[1].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
 301         val->flags |= !(param[1].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL  : 0;
 302         val->flags |= CMP_REALVAL;
 303         n           = param[1].val.nr;
 304         val->r      = param[1].val.u.r;
 305     }
 306     else
 307     {
 308         n           = 0;
 309         val->i      = NULL;
 310         val->r      = NULL;
 311     }
 312     return n;
 313 }
 314
 315 /* \brief
 316  * Converts an integer value to floating point.
 317  *
 318  * \param[in]     n   Number of values in the \p val->u array.
 319  * \param[in,out] val Value to convert.
 320  */
 321 static void
 322 convert_int_real(int n, t_compare_value *val)
 323 {
 324     int   i;
 325     real *rv;
 326
 327     snew(rv, n);
 328     for (i = 0; i < n; ++i)
 329     {
 330         rv[i] = (real)val->i[i];
 331     }
 332     /* Free the previous value if one is present. */
 333     sfree(val->r);
 334     val->r      = rv;
 335     val->flags |= CMP_REALVAL | CMP_ALLOCREAL;
 336 }
 337
 338 /* \brief
 339  * Converts a floating point value to integer.
 340  *
 341  * \param[in]     n      Number of values in the \p val->u array.
 342  * \param[in,out] val    Value to convert.
 343  * \param[in]     cmpt   Comparison operator type.
 344  * \param[in]     bRight true if \p val appears on the right hand size of
 345  *   \p cmpt.
 346  * \returns       0 on success, EINVAL on error.
 347  *
 348  * The values are rounded such that the same comparison operator can be used.
 349  */
 350 static void
 351 convert_real_int(int n, t_compare_value *val, e_comparison_t cmpt, bool bRight)
 352 {
 353     int   i;
 354     int  *iv;
 355
 356     if (!bRight)
 357     {
 358         cmpt = reverse_comparison_type(cmpt);
 359     }
 360     snew(iv, n);
 361     /* Round according to the comparison type */
 362     for (i = 0; i < n; ++i)
 363     {
 364         switch (cmpt)
 365         {
 366             case CMP_LESS:
 367             case CMP_GEQ:
 368                 iv[i] = (int)ceil(val->r[i]);
 369                 break;
 370             case CMP_GTR:
 371             case CMP_LEQ:
 372                 iv[i] = (int)floor(val->r[i]);
 373                 break;
 374             case CMP_EQUAL:
 375             case CMP_NEQ:
 376                 sfree(iv);
 377                 /* TODO: Implement, although it isn't typically very useful.
 378                  * Implementation is only a matter or proper initialization,
 379                  * the evaluation function can already handle this case with
 380                  * proper preparations. */
 381                 GMX_THROW(gmx::NotImplementedError("Equality comparison between dynamic integer and static real expressions not implemented"));
 382             case CMP_INVALID: /* Should not be reached */
 383                 sfree(iv);
 384                 GMX_THROW(gmx::InternalError("Invalid comparison type"));
 385         }
 386     }
 387     /* Free the previous value if one is present. */
 388     sfree(val->i);
 389     val->i      = iv;
 390     val->flags &= ~CMP_REALVAL;
 391     val->flags |= CMP_ALLOCINT;
 392 }
 393
 394 /*!
 395  * \param[in] top   Not used.
 396  * \param[in] npar  Not used (should be 5).
 397  * \param[in] param Method parameters (should point to \ref smparams_compare).
 398  * \param[in] data  Should point to a \c t_methoddata_compare.
 399  * \returns   0 if the input data is valid, -1 on error.
 400  */
 401 static void
 402 init_compare(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
 403 {
 404     t_methoddata_compare *d = (t_methoddata_compare *)data;
 405     int                   n1, n2;
 406
 407     /* Store the values */
 408     n1 = init_comparison_value(&d->left, &param[0]);
 409     n2 = init_comparison_value(&d->right, &param[3]);
 410     if (n1 == 0 || n2 == 0)
 411     {
 412         GMX_THROW(gmx::InternalError("One of the values for comparison missing"));
 413     }
 414     /* Store the comparison type */
 415     d->cmpt = comparison_type(d->cmpop);
 416     if (d->cmpt == CMP_INVALID)
 417     {
 418         GMX_THROW(gmx::InternalError("Invalid comparison type"));
 419     }
 420     /* Convert the values to the same type */
 421     if ((d->left.flags & CMP_REALVAL) && !(d->right.flags & CMP_REALVAL))
 422     {
 423         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 424         {
 425             /* Nothing can be done */
 426         }
 427         else if (!(d->right.flags & CMP_DYNAMICVAL))
 428         {
 429             convert_int_real(n2, &d->right);
 430         }
 431         else /* d->left is static */
 432         {
 433             convert_real_int(n1, &d->left, d->cmpt, false);
 434         }
 435     }
 436     else if (!(d->left.flags & CMP_REALVAL) && (d->right.flags & CMP_REALVAL))
 437     {
 438         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 439         {
 440             /* Reverse the sides to place the integer on the right */
 441             int    flags;
 442             d->left.r  = d->right.r;
 443             d->right.r = NULL;
 444             d->right.i = d->left.i;
 445             d->left.i  = NULL;
 446             flags          = d->left.flags;
 447             d->left.flags  = d->right.flags;
 448             d->right.flags = flags;
 449             d->cmpt = reverse_comparison_type(d->cmpt);
 450         }
 451         else if (!(d->left.flags & CMP_DYNAMICVAL))
 452         {
 453             convert_int_real(n1, &d->left);
 454         }
 455         else /* d->right is static */
 456         {
 457             convert_real_int(n2, &d->right, d->cmpt, true);
 458         }
 459     }
 460 }
 461
 462 /*!
 463  * \param data Data to free (should point to a \c t_methoddata_compare).
 464  *
 465  * Frees the memory allocated for \c t_methoddata_compare.
 466  */
 467 static void
 468 free_data_compare(void *data)
 469 {
 470     t_methoddata_compare *d = (t_methoddata_compare *)data;
 471
 472     sfree(d->cmpop);
 473     if (d->left.flags & CMP_ALLOCINT)
 474     {
 475         sfree(d->left.i);
 476     }
 477     if (d->left.flags & CMP_ALLOCREAL)
 478     {
 479         sfree(d->left.r);
 480     }
 481     if (d->right.flags & CMP_ALLOCINT)
 482     {
 483         sfree(d->right.i);
 484     }
 485     if (d->right.flags & CMP_ALLOCREAL)
 486     {
 487         sfree(d->right.r);
 488     }
 489 }
 490
 491 /*!
 492  * \param[in]  top   Not used.
 493  * \param[in]  fr    Not used.
 494  * \param[in]  pbc   Not used.
 495  * \param[in]  g     Evaluation index group.
 496  * \param[out] out   Output data structure (\p out->u.g is used).
 497  * \param[in]  data  Should point to a \c t_methoddata_compare.
 498  */
 499 static void
 500 evaluate_compare_int(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 501                      gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
 502 {
 503     t_methoddata_compare *d = (t_methoddata_compare *)data;
 504     int                   i, i1, i2, ig;
 505     int                   a, b;
 506     bool                  bAccept;
 507
 508     for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
 509     {
 510         a = d->left.i[i1];
 511         b = d->right.i[i2];
 512         bAccept = false;
 513         switch (d->cmpt)
 514         {
 515             case CMP_INVALID: break;
 516             case CMP_LESS:    bAccept = a <  b; break;
 517             case CMP_LEQ:     bAccept = a <= b; break;
 518             case CMP_GTR:     bAccept = a >  b; break;
 519             case CMP_GEQ:     bAccept = a >= b; break;
 520             case CMP_EQUAL:   bAccept = a == b; break;
 521             case CMP_NEQ:     bAccept = a != b; break;
 522         }
 523         if (bAccept)
 524         {
 525             out->u.g->index[ig++] = g->index[i];
 526         }
 527         if (!(d->left.flags & CMP_SINGLEVAL))
 528         {
 529             ++i1;
 530         }
 531         if (!(d->right.flags & CMP_SINGLEVAL))
 532         {
 533             ++i2;
 534         }
 535     }
 536     out->u.g->isize = ig;
 537 }
 538
 539 /*!
 540  * \param[in]  top   Not used.
 541  * \param[in]  fr    Not used.
 542  * \param[in]  pbc   Not used.
 543  * \param[in]  g     Evaluation index group.
 544  * \param[out] out   Output data structure (\p out->u.g is used).
 545  * \param[in]  data  Should point to a \c t_methoddata_compare.
 546  */
 547 static void
 548 evaluate_compare_real(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 549                       gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
 550 {
 551     t_methoddata_compare *d = (t_methoddata_compare *)data;
 552     int                   i, i1, i2, ig;
 553     real                  a, b;
 554     bool                  bAccept;
 555
 556     for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
 557     {
 558         a = d->left.r[i1];
 559         b = (d->right.flags & CMP_REALVAL) ? d->right.r[i2] : d->right.i[i2];
 560         bAccept = false;
 561         switch (d->cmpt)
 562         {
 563             case CMP_INVALID: break;
 564             case CMP_LESS:    bAccept = a <  b; break;
 565             case CMP_LEQ:     bAccept = a <= b; break;
 566             case CMP_GTR:     bAccept = a >  b; break;
 567             case CMP_GEQ:     bAccept = a >= b; break;
 568             case CMP_EQUAL:   bAccept =  gmx_within_tol(a, b, GMX_REAL_EPS); break;
 569             case CMP_NEQ:     bAccept = !gmx_within_tol(a, b, GMX_REAL_EPS); break;
 570         }
 571         if (bAccept)
 572         {
 573             out->u.g->index[ig++] = g->index[i];
 574         }
 575         if (!(d->left.flags & CMP_SINGLEVAL))
 576         {
 577             ++i1;
 578         }
 579         if (!(d->right.flags & CMP_SINGLEVAL))
 580         {
 581             ++i2;
 582         }
 583     }
 584     out->u.g->isize = ig;
 585 }
 586
 587 /*!
 588  * \param[in]  top   Not used.
 589  * \param[in]  fr    Not used.
 590  * \param[in]  pbc   Not used.
 591  * \param[in]  g     Evaluation index group.
 592  * \param[out] out   Output data structure (\p out->u.g is used).
 593  * \param[in]  data  Should point to a \c t_methoddata_compare.
 594  */
 595 static void
 596 evaluate_compare(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 597                  gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
 598 {
 599     t_methoddata_compare *d = (t_methoddata_compare *)data;
 600
 601     if (!((d->left.flags | d->right.flags) & CMP_REALVAL))
 602     {
 603         evaluate_compare_int(top, fr, pbc, g, out, data);
 604     }
 605     else
 606     {
 607         evaluate_compare_real(top, fr, pbc, g, out, data);
 608     }
 609 }