src/gromacs/selection/sm_compare.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   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.
   8  *
   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.
  13  *
  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.
  18  *
  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.
  23  *
  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.
  31  *
  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.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements internal selection method for comparison expressions.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_selection
  41  */
  42 #include "gromacs/legacyheaders/maths.h"
  43 #include "gromacs/legacyheaders/macros.h"
  44 #include "gromacs/legacyheaders/smalloc.h"
  45
  46 #include "gromacs/selection/selmethod.h"
  47 #include "gromacs/utility/common.h"
  48 #include "gromacs/utility/exceptions.h"
  49
  50 /** Defines the comparison operator for comparison expressions. */
  51 typedef enum
  52 {
  53     CMP_INVALID,        /**< Indicates an error */
  54     CMP_LESS,           /**< '<' */
  55     CMP_LEQ,            /**< '<=' */
  56     CMP_GTR,            /**< '>' */
  57     CMP_GEQ,            /**< '>=' */
  58     CMP_EQUAL,          /**< '==' */
  59     CMP_NEQ             /**< '!=' */
  60 } e_comparison_t;
  61
  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. */
  67 #define CMP_REALVAL    4
  68 /** The integer array is allocated. */
  69 #define CMP_ALLOCINT   16
  70 /** The real array is allocated. */
  71 #define CMP_ALLOCREAL  32
  72
  73 /*! \internal \brief
  74  * Data structure for comparison expression operand values.
  75  */
  76 typedef struct
  77 {
  78     /** Flags that describe the type of the operand. */
  79     int             flags;
  80     /** (Array of) integer value(s). */
  81     int            *i;
  82     /** (Array of) real value(s). */
  83     real           *r;
  84 } t_compare_value;
  85
  86 /*! \internal \brief
  87  * Data structure for comparison expression evaluation.
  88  */
  89 typedef struct
  90 {
  91     /** Comparison operator as a string. */
  92     char            *cmpop;
  93     /** Comparison operator type. */
  94     e_comparison_t   cmpt;
  95     /** Left value. */
  96     t_compare_value  left;
  97     /** Right value. */
  98     t_compare_value  right;
  99 } t_methoddata_compare;
 100
 101 /*! \brief
 102  * Allocates data for comparison expression evaluation.
 103  *
 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).
 108  *
 109  * Allocates memory for a \c t_methoddata_compare structure.
 110  */
 111 static void *
 112 init_data_compare(int npar, gmx_ana_selparam_t *param);
 113 /*! \brief
 114  * Initializes data for comparison expression evaluation.
 115  *
 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.
 121  */
 122 static void
 123 init_compare(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
 124 /** Frees the memory allocated for comparison expression evaluation. */
 125 static void
 126 free_data_compare(void *data);
 127 /** Evaluates comparison expressions. */
 128 static void
 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);
 131
 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},
 143 };
 144
 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,
 149     &init_data_compare,
 150     NULL,
 151     &init_compare,
 152     NULL,
 153     &free_data_compare,
 154     NULL,
 155     &evaluate_compare,
 156     NULL,
 157     {NULL, 0, NULL},
 158 };
 159
 160 /*! \brief
 161  * Returns a \c e_comparison_t value corresponding to an operator.
 162  *
 163  * \param[in] str  String to process.
 164  * \returns   The comparison type corresponding to the first one or two
 165  *   characters of \p str.
 166  *
 167  * \p str can contain any number of characters; only the first two
 168  * are used.
 169  * If the beginning of \p str does not match any of the recognized types,
 170  * \ref CMP_INVALID is returned.
 171  */
 172 static e_comparison_t
 173 comparison_type(char *str)
 174 {
 175     switch (str[0])
 176     {
 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;
 181     }
 182     return CMP_INVALID;
 183 }
 184
 185 /*! \brief
 186  * Returns a string corresponding to a \c e_comparison_t value.
 187  *
 188  * \param[in] cmpt  Comparison type to convert.
 189  * \returns   Pointer to a string that corresponds to \p cmpt.
 190  *
 191  * The return value points to a string constant and should not be \p free'd.
 192  *
 193  * The function returns NULL if \p cmpt is not one of the valid values.
 194  */
 195 static const char *
 196 comparison_type_str(e_comparison_t cmpt)
 197 {
 198     switch (cmpt)
 199     {
 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;
 207     }
 208     return NULL;
 209 }
 210
 211 /*!
 212  * \param[in] fp    File to receive the output.
 213  * \param[in] data  Should point to a \c t_methoddata_compare.
 214  */
 215 void
 216 _gmx_selelem_print_compare_info(FILE *fp, void *data)
 217 {
 218     t_methoddata_compare *d = (t_methoddata_compare *)data;
 219
 220     fprintf(fp, " \"");
 221     /* Print the left value */
 222     if ((d->left.flags & CMP_SINGLEVAL) && !(d->left.flags & CMP_DYNAMICVAL))
 223     {
 224         if (d->left.flags & CMP_REALVAL)
 225         {
 226             fprintf(fp, "%f ", d->left.r[0]);
 227         }
 228         else
 229         {
 230             fprintf(fp, "%d ", d->left.i[0]);
 231         }
 232     }
 233     /* Print the operator */
 234     if (d->cmpt != CMP_INVALID)
 235     {
 236         fprintf(fp, "%s", comparison_type_str(d->cmpt));
 237     }
 238     else
 239     {
 240         fprintf(fp, "%s", d->cmpop);
 241     }
 242     /* Print the right value */
 243     if ((d->right.flags & CMP_SINGLEVAL) && !(d->right.flags & CMP_DYNAMICVAL))
 244     {
 245         if (d->right.flags & CMP_REALVAL)
 246         {
 247             fprintf(fp, " %f", d->right.r[0]);
 248         }
 249         else
 250         {
 251             fprintf(fp, " %d", d->right.i[0]);
 252         }
 253     }
 254     fprintf(fp, "\"");
 255 }
 256
 257 static void *
 258 init_data_compare(int /* npar */, gmx_ana_selparam_t *param)
 259 {
 260     t_methoddata_compare *data;
 261
 262     snew(data, 1);
 263     param[2].val.u.s = &data->cmpop;
 264     return data;
 265 }
 266
 267 /*! \brief
 268  * Reverses a comparison operator.
 269  *
 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.
 273  */
 274 static e_comparison_t
 275 reverse_comparison_type(e_comparison_t type)
 276 {
 277     switch (type)
 278     {
 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;
 283         default:       break;
 284     }
 285     return type;
 286 }
 287
 288 /*! \brief
 289  * Initializes the value storage for comparison expression.
 290  *
 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.
 294  */
 295 static int
 296 init_comparison_value(t_compare_value *val, gmx_ana_selparam_t param[2])
 297 {
 298     int  n;
 299
 300     val->flags = 0;
 301     if (param[0].flags & SPAR_SET)
 302     {
 303         val->flags |=  (param[0].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
 304         val->flags |= !(param[0].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL  : 0;
 305         n           = param[0].val.nr;
 306         val->i      = param[0].val.u.i;
 307     }
 308     else if (param[1].flags & SPAR_SET)
 309     {
 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;
 313         n           = param[1].val.nr;
 314         val->r      = param[1].val.u.r;
 315     }
 316     else
 317     {
 318         n           = 0;
 319         val->i      = NULL;
 320         val->r      = NULL;
 321     }
 322     return n;
 323 }
 324
 325 /*! \brief
 326  * Converts an integer value to floating point.
 327  *
 328  * \param[in]     n   Number of values in the \p val->u array.
 329  * \param[in,out] val Value to convert.
 330  */
 331 static void
 332 convert_int_real(int n, t_compare_value *val)
 333 {
 334     int   i;
 335     real *rv;
 336
 337     snew(rv, n);
 338     for (i = 0; i < n; ++i)
 339     {
 340         rv[i] = (real)val->i[i];
 341     }
 342     /* Free the previous value if one is present. */
 343     sfree(val->r);
 344     val->r      = rv;
 345     val->flags |= CMP_REALVAL | CMP_ALLOCREAL;
 346 }
 347
 348 /*! \brief
 349  * Converts a floating point value to integer.
 350  *
 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
 355  *   \p cmpt.
 356  * \returns       0 on success, EINVAL on error.
 357  *
 358  * The values are rounded such that the same comparison operator can be used.
 359  */
 360 static void
 361 convert_real_int(int n, t_compare_value *val, e_comparison_t cmpt, bool bRight)
 362 {
 363     int   i;
 364     int  *iv;
 365
 366     if (!bRight)
 367     {
 368         cmpt = reverse_comparison_type(cmpt);
 369     }
 370     snew(iv, n);
 371     /* Round according to the comparison type */
 372     for (i = 0; i < n; ++i)
 373     {
 374         switch (cmpt)
 375         {
 376             case CMP_LESS:
 377             case CMP_GEQ:
 378                 iv[i] = (int)ceil(val->r[i]);
 379                 break;
 380             case CMP_GTR:
 381             case CMP_LEQ:
 382                 iv[i] = (int)floor(val->r[i]);
 383                 break;
 384             case CMP_EQUAL:
 385             case CMP_NEQ:
 386                 sfree(iv);
 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 */
 393                 sfree(iv);
 394                 GMX_THROW(gmx::InternalError("Invalid comparison type"));
 395         }
 396     }
 397     /* Free the previous value if one is present. */
 398     sfree(val->i);
 399     val->i      = iv;
 400     val->flags &= ~CMP_REALVAL;
 401     val->flags |= CMP_ALLOCINT;
 402 }
 403
 404 static void
 405 init_compare(t_topology * /* top */, int /* npar */, gmx_ana_selparam_t *param, void *data)
 406 {
 407     t_methoddata_compare *d = (t_methoddata_compare *)data;
 408     int                   n1, n2;
 409
 410     /* Store the values */
 411     n1 = init_comparison_value(&d->left, &param[0]);
 412     n2 = init_comparison_value(&d->right, &param[3]);
 413     /* Store the comparison type */
 414     d->cmpt = comparison_type(d->cmpop);
 415     if (d->cmpt == CMP_INVALID)
 416     {
 417         GMX_THROW(gmx::InternalError("Invalid comparison type"));
 418     }
 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.
 424      */
 425     if ((d->left.flags & CMP_REALVAL) && !(d->right.flags & CMP_REALVAL))
 426     {
 427         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 428         {
 429             /* Nothing can be done */
 430         }
 431         else if (!(d->right.flags & CMP_DYNAMICVAL))
 432         {
 433             convert_int_real(n2, &d->right);
 434         }
 435         else /* d->left is static */
 436         {
 437             convert_real_int(n1, &d->left, d->cmpt, false);
 438         }
 439     }
 440     else if (!(d->left.flags & CMP_REALVAL) && (d->right.flags & CMP_REALVAL))
 441     {
 442         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 443         {
 444             /* Reverse the sides to place the integer on the right */
 445             int    flags;
 446             d->left.r      = d->right.r;
 447             d->right.r     = NULL;
 448             d->right.i     = d->left.i;
 449             d->left.i      = NULL;
 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);
 454         }
 455         else if (!(d->left.flags & CMP_DYNAMICVAL))
 456         {
 457             convert_int_real(n1, &d->left);
 458         }
 459         else /* d->right is static */
 460         {
 461             convert_real_int(n2, &d->right, d->cmpt, true);
 462         }
 463     }
 464 }
 465
 466 /*!
 467  * \param data Data to free (should point to a \c t_methoddata_compare).
 468  *
 469  * Frees the memory allocated for \c t_methoddata_compare.
 470  */
 471 static void
 472 free_data_compare(void *data)
 473 {
 474     t_methoddata_compare *d = (t_methoddata_compare *)data;
 475
 476     sfree(d->cmpop);
 477     if (d->left.flags & CMP_ALLOCINT)
 478     {
 479         sfree(d->left.i);
 480     }
 481     if (d->left.flags & CMP_ALLOCREAL)
 482     {
 483         sfree(d->left.r);
 484     }
 485     if (d->right.flags & CMP_ALLOCINT)
 486     {
 487         sfree(d->right.i);
 488     }
 489     if (d->right.flags & CMP_ALLOCREAL)
 490     {
 491         sfree(d->right.r);
 492     }
 493     sfree(d);
 494 }
 495
 496 /*! \brief
 497  * Implementation for evaluate_compare() for integer values.
 498  *
 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.
 505  */
 506 static void
 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)
 509 {
 510     t_methoddata_compare *d = (t_methoddata_compare *)data;
 511     int                   i, i1, i2, ig;
 512     int                   a, b;
 513     bool                  bAccept;
 514
 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)
 519     {
 520         a       = d->left.i[i1];
 521         b       = d->right.i[i2];
 522         bAccept = false;
 523         switch (d->cmpt)
 524         {
 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;
 532         }
 533         if (bAccept)
 534         {
 535             out->u.g->index[ig++] = g->index[i];
 536         }
 537         if (!(d->left.flags & CMP_SINGLEVAL))
 538         {
 539             ++i1;
 540         }
 541         if (!(d->right.flags & CMP_SINGLEVAL))
 542         {
 543             ++i2;
 544         }
 545     }
 546     out->u.g->isize = ig;
 547 }
 548
 549 /*! \brief
 550  * Implementation for evaluate_compare() if either value is non-integer.
 551  *
 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.
 558  *
 559  * Left value is assumed to be real-valued; right value can be either.
 560  * This is ensured by the initialization method.
 561  */
 562 static void
 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)
 565 {
 566     t_methoddata_compare *d = (t_methoddata_compare *)data;
 567     int                   i, i1, i2, ig;
 568     real                  a, b;
 569     bool                  bAccept;
 570
 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)
 575     {
 576         a       = d->left.r[i1];
 577         b       = (d->right.flags & CMP_REALVAL) ? d->right.r[i2] : d->right.i[i2];
 578         bAccept = false;
 579         switch (d->cmpt)
 580         {
 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;
 588         }
 589         if (bAccept)
 590         {
 591             out->u.g->index[ig++] = g->index[i];
 592         }
 593         if (!(d->left.flags & CMP_SINGLEVAL))
 594         {
 595             ++i1;
 596         }
 597         if (!(d->right.flags & CMP_SINGLEVAL))
 598         {
 599             ++i2;
 600         }
 601     }
 602     out->u.g->isize = ig;
 603 }
 604
 605 /*!
 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.
 612  */
 613 static void
 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)
 616 {
 617     t_methoddata_compare *d = (t_methoddata_compare *)data;
 618
 619     if (!((d->left.flags | d->right.flags) & CMP_REALVAL))
 620     {
 621         evaluate_compare_int(top, fr, pbc, g, out, data);
 622     }
 623     else
 624     {
 625         evaluate_compare_real(top, fr, pbc, g, out, data);
 626     }
 627 }