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