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 #include "gromacs/legacyheaders/maths.h"
  39 #include "gromacs/legacyheaders/macros.h"
  40 #include "gromacs/legacyheaders/smalloc.h"
  41
  42 #include "gromacs/selection/selmethod.h"
  43 #include "gromacs/utility/exceptions.h"
  44
  45 /** Defines the comparison operator for comparison expressions. */
  46 typedef enum
  47 {
  48     CMP_INVALID,        /**< Indicates an error */
  49     CMP_LESS,           /**< '<' */
  50     CMP_LEQ,            /**< '<=' */
  51     CMP_GTR,            /**< '>' */
  52     CMP_GEQ,            /**< '>=' */
  53     CMP_EQUAL,          /**< '==' */
  54     CMP_NEQ             /**< '!=' */
  55 } e_comparison_t;
  56
  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. */
  62 #define CMP_REALVAL    4
  63 /** The integer array is allocated. */
  64 #define CMP_ALLOCINT   16
  65 /** The real array is allocated. */
  66 #define CMP_ALLOCREAL  32
  67
  68 /*! \internal \brief
  69  * Data structure for comparison expression operand values.
  70  */
  71 typedef struct
  72 {
  73     /** Flags that describe the type of the operand. */
  74     int             flags;
  75     /** (Array of) integer value(s). */
  76     int            *i;
  77     /** (Array of) real value(s). */
  78     real           *r;
  79 } t_compare_value;
  80
  81 /*! \internal \brief
  82  * Data structure for comparison expression evaluation.
  83  */
  84 typedef struct
  85 {
  86     /** Comparison operator as a string. */
  87     char            *cmpop;
  88     /** Comparison operator type. */
  89     e_comparison_t   cmpt;
  90     /** Left value. */
  91     t_compare_value  left;
  92     /** Right value. */
  93     t_compare_value  right;
  94 } t_methoddata_compare;
  95
  96 /** Allocates data for comparison expression evaluation. */
  97 static void *
  98 init_data_compare(int npar, gmx_ana_selparam_t *param);
  99 /** Initializes data for comparison expression evaluation. */
 100 static void
 101 init_compare(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
 102 /** Frees the memory allocated for comparison expression evaluation. */
 103 static void
 104 free_data_compare(void *data);
 105 /** Evaluates comparison expressions. */
 106 static void
 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);
 109
 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},
 121 };
 122
 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,
 127     &init_data_compare,
 128     NULL,
 129     &init_compare,
 130     NULL,
 131     &free_data_compare,
 132     NULL,
 133     &evaluate_compare,
 134     NULL,
 135     {NULL, 0, NULL},
 136 };
 137
 138 /*! \brief
 139  * Returns a \c e_comparison_t value corresponding to an operator.
 140  *
 141  * \param[in] str  String to process.
 142  * \returns   The comparison type corresponding to the first one or two
 143  *   characters of \p str.
 144  *
 145  * \p str can contain any number of characters; only the first two
 146  * are used.
 147  * If the beginning of \p str does not match any of the recognized types,
 148  * \ref CMP_INVALID is returned.
 149  */
 150 static e_comparison_t
 151 comparison_type(char *str)
 152 {
 153     switch (str[0])
 154     {
 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;
 159     }
 160     return CMP_INVALID;
 161 }
 162
 163 /*! \brief
 164  * Returns a string corresponding to a \c e_comparison_t value.
 165  *
 166  * \param[in] cmpt  Comparison type to convert.
 167  * \returns   Pointer to a string that corresponds to \p cmpt.
 168  *
 169  * The return value points to a string constant and should not be \p free'd.
 170  *
 171  * The function returns NULL if \p cmpt is not one of the valid values.
 172  */
 173 static const char *
 174 comparison_type_str(e_comparison_t cmpt)
 175 {
 176     switch (cmpt)
 177     {
 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;
 185     }
 186     return NULL;
 187 }
 188
 189 /*!
 190  * \param[in] fp    File to receive the output.
 191  * \param[in] data  Should point to a \c t_methoddata_compare.
 192  */
 193 void
 194 _gmx_selelem_print_compare_info(FILE *fp, void *data)
 195 {
 196     t_methoddata_compare *d = (t_methoddata_compare *)data;
 197
 198     fprintf(fp, " \"");
 199     /* Print the left value */
 200     if ((d->left.flags & CMP_SINGLEVAL) && !(d->left.flags & CMP_DYNAMICVAL))
 201     {
 202         if (d->left.flags & CMP_REALVAL)
 203         {
 204             fprintf(fp, "%f ", d->left.r[0]);
 205         }
 206         else
 207         {
 208             fprintf(fp, "%d ", d->left.i[0]);
 209         }
 210     }
 211     /* Print the operator */
 212     if (d->cmpt != CMP_INVALID)
 213     {
 214         fprintf(fp, "%s", comparison_type_str(d->cmpt));
 215     }
 216     else
 217     {
 218         fprintf(fp, "%s", d->cmpop);
 219     }
 220     /* Print the right value */
 221     if ((d->right.flags & CMP_SINGLEVAL) && !(d->right.flags & CMP_DYNAMICVAL))
 222     {
 223         if (d->right.flags & CMP_REALVAL)
 224         {
 225             fprintf(fp, " %f", d->right.r[0]);
 226         }
 227         else
 228         {
 229             fprintf(fp, " %d", d->right.i[0]);
 230         }
 231     }
 232     fprintf(fp, "\"");
 233 }
 234
 235 /*!
 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).
 240  *
 241  * Allocates memory for a \c t_methoddata_compare structure.
 242  */
 243 static void *
 244 init_data_compare(int npar, gmx_ana_selparam_t *param)
 245 {
 246     t_methoddata_compare *data;
 247
 248     snew(data, 1);
 249     param[2].val.u.s = &data->cmpop;
 250     return data;
 251 }
 252
 253 /*! \brief
 254  * Reverses a comparison operator.
 255  *
 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.
 259  */
 260 static e_comparison_t
 261 reverse_comparison_type(e_comparison_t type)
 262 {
 263     switch (type)
 264     {
 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;
 269         default:       break;
 270     }
 271     return type;
 272 }
 273
 274 /*! \brief
 275  * Initializes the value storage for comparison expression.
 276  *
 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.
 280  */
 281 static int
 282 init_comparison_value(t_compare_value *val, gmx_ana_selparam_t param[2])
 283 {
 284     int  n;
 285
 286     val->flags = 0;
 287     if (param[0].flags & SPAR_SET)
 288     {
 289         val->flags |=  (param[0].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
 290         val->flags |= !(param[0].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL  : 0;
 291         n           = param[0].val.nr;
 292         val->i      = param[0].val.u.i;
 293     }
 294     else if (param[1].flags & SPAR_SET)
 295     {
 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;
 299         n           = param[1].val.nr;
 300         val->r      = param[1].val.u.r;
 301     }
 302     else
 303     {
 304         n           = 0;
 305         val->i      = NULL;
 306         val->r      = NULL;
 307     }
 308     return n;
 309 }
 310
 311 /*! \brief
 312  * Converts an integer value to floating point.
 313  *
 314  * \param[in]     n   Number of values in the \p val->u array.
 315  * \param[in,out] val Value to convert.
 316  */
 317 static void
 318 convert_int_real(int n, t_compare_value *val)
 319 {
 320     int   i;
 321     real *rv;
 322
 323     snew(rv, n);
 324     for (i = 0; i < n; ++i)
 325     {
 326         rv[i] = (real)val->i[i];
 327     }
 328     /* Free the previous value if one is present. */
 329     sfree(val->r);
 330     val->r      = rv;
 331     val->flags |= CMP_REALVAL | CMP_ALLOCREAL;
 332 }
 333
 334 /*! \brief
 335  * Converts a floating point value to integer.
 336  *
 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
 341  *   \p cmpt.
 342  * \returns       0 on success, EINVAL on error.
 343  *
 344  * The values are rounded such that the same comparison operator can be used.
 345  */
 346 static void
 347 convert_real_int(int n, t_compare_value *val, e_comparison_t cmpt, bool bRight)
 348 {
 349     int   i;
 350     int  *iv;
 351
 352     if (!bRight)
 353     {
 354         cmpt = reverse_comparison_type(cmpt);
 355     }
 356     snew(iv, n);
 357     /* Round according to the comparison type */
 358     for (i = 0; i < n; ++i)
 359     {
 360         switch (cmpt)
 361         {
 362             case CMP_LESS:
 363             case CMP_GEQ:
 364                 iv[i] = (int)ceil(val->r[i]);
 365                 break;
 366             case CMP_GTR:
 367             case CMP_LEQ:
 368                 iv[i] = (int)floor(val->r[i]);
 369                 break;
 370             case CMP_EQUAL:
 371             case CMP_NEQ:
 372                 sfree(iv);
 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 */
 379                 sfree(iv);
 380                 GMX_THROW(gmx::InternalError("Invalid comparison type"));
 381         }
 382     }
 383     /* Free the previous value if one is present. */
 384     sfree(val->i);
 385     val->i      = iv;
 386     val->flags &= ~CMP_REALVAL;
 387     val->flags |= CMP_ALLOCINT;
 388 }
 389
 390 /*!
 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.
 396  */
 397 static void
 398 init_compare(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
 399 {
 400     t_methoddata_compare *d = (t_methoddata_compare *)data;
 401     int                   n1, n2;
 402
 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)
 407     {
 408         GMX_THROW(gmx::InternalError("One of the values for comparison missing"));
 409     }
 410     /* Store the comparison type */
 411     d->cmpt = comparison_type(d->cmpop);
 412     if (d->cmpt == CMP_INVALID)
 413     {
 414         GMX_THROW(gmx::InternalError("Invalid comparison type"));
 415     }
 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.
 421      */
 422     if ((d->left.flags & CMP_REALVAL) && !(d->right.flags & CMP_REALVAL))
 423     {
 424         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 425         {
 426             /* Nothing can be done */
 427         }
 428         else if (!(d->right.flags & CMP_DYNAMICVAL))
 429         {
 430             convert_int_real(n2, &d->right);
 431         }
 432         else /* d->left is static */
 433         {
 434             convert_real_int(n1, &d->left, d->cmpt, false);
 435         }
 436     }
 437     else if (!(d->left.flags & CMP_REALVAL) && (d->right.flags & CMP_REALVAL))
 438     {
 439         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 440         {
 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);
 451         }
 452         else if (!(d->left.flags & CMP_DYNAMICVAL))
 453         {
 454             convert_int_real(n1, &d->left);
 455         }
 456         else /* d->right is static */
 457         {
 458             convert_real_int(n2, &d->right, d->cmpt, true);
 459         }
 460     }
 461 }
 462
 463 /*!
 464  * \param data Data to free (should point to a \c t_methoddata_compare).
 465  *
 466  * Frees the memory allocated for \c t_methoddata_compare.
 467  */
 468 static void
 469 free_data_compare(void *data)
 470 {
 471     t_methoddata_compare *d = (t_methoddata_compare *)data;
 472
 473     sfree(d->cmpop);
 474     if (d->left.flags & CMP_ALLOCINT)
 475     {
 476         sfree(d->left.i);
 477     }
 478     if (d->left.flags & CMP_ALLOCREAL)
 479     {
 480         sfree(d->left.r);
 481     }
 482     if (d->right.flags & CMP_ALLOCINT)
 483     {
 484         sfree(d->right.i);
 485     }
 486     if (d->right.flags & CMP_ALLOCREAL)
 487     {
 488         sfree(d->right.r);
 489     }
 490     sfree(d);
 491 }
 492
 493 /*!
 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.
 500  */
 501 static void
 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)
 504 {
 505     t_methoddata_compare *d = (t_methoddata_compare *)data;
 506     int                   i, i1, i2, ig;
 507     int                   a, b;
 508     bool                  bAccept;
 509
 510     for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
 511     {
 512         a       = d->left.i[i1];
 513         b       = d->right.i[i2];
 514         bAccept = false;
 515         switch (d->cmpt)
 516         {
 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;
 524         }
 525         if (bAccept)
 526         {
 527             out->u.g->index[ig++] = g->index[i];
 528         }
 529         if (!(d->left.flags & CMP_SINGLEVAL))
 530         {
 531             ++i1;
 532         }
 533         if (!(d->right.flags & CMP_SINGLEVAL))
 534         {
 535             ++i2;
 536         }
 537     }
 538     out->u.g->isize = ig;
 539 }
 540
 541 /*!
 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.
 548  */
 549 static void
 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)
 552 {
 553     t_methoddata_compare *d = (t_methoddata_compare *)data;
 554     int                   i, i1, i2, ig;
 555     real                  a, b;
 556     bool                  bAccept;
 557
 558     for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
 559     {
 560         a       = d->left.r[i1];
 561         b       = (d->right.flags & CMP_REALVAL) ? d->right.r[i2] : d->right.i[i2];
 562         bAccept = false;
 563         switch (d->cmpt)
 564         {
 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;
 572         }
 573         if (bAccept)
 574         {
 575             out->u.g->index[ig++] = g->index[i];
 576         }
 577         if (!(d->left.flags & CMP_SINGLEVAL))
 578         {
 579             ++i1;
 580         }
 581         if (!(d->right.flags & CMP_SINGLEVAL))
 582         {
 583             ++i2;
 584         }
 585     }
 586     out->u.g->isize = ig;
 587 }
 588
 589 /*!
 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.
 596  */
 597 static void
 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)
 600 {
 601     t_methoddata_compare *d = (t_methoddata_compare *)data;
 602
 603     if (!((d->left.flags | d->right.flags) & CMP_REALVAL))
 604     {
 605         evaluate_compare_int(top, fr, pbc, g, out, data);
 606     }
 607     else
 608     {
 609         evaluate_compare_real(top, fr, pbc, g, out, data);
 610     }
 611 }