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, 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/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     /* TODO: Currently, there are no dynamic integer-valued selection methods,
 422      * which means that only the branches with convert_int_real() will ever be
 423      * taken. It should be considered whether it is necessary to support these
 424      * other cases at all.
 425      */
 426     if ((d->left.flags & CMP_REALVAL) && !(d->right.flags & CMP_REALVAL))
 427     {
 428         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 429         {
 430             /* Nothing can be done */
 431         }
 432         else if (!(d->right.flags & CMP_DYNAMICVAL))
 433         {
 434             convert_int_real(n2, &d->right);
 435         }
 436         else /* d->left is static */
 437         {
 438             convert_real_int(n1, &d->left, d->cmpt, false);
 439         }
 440     }
 441     else if (!(d->left.flags & CMP_REALVAL) && (d->right.flags & CMP_REALVAL))
 442     {
 443         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 444         {
 445             /* Reverse the sides to place the integer on the right */
 446             int    flags;
 447             d->left.r      = d->right.r;
 448             d->right.r     = NULL;
 449             d->right.i     = d->left.i;
 450             d->left.i      = NULL;
 451             flags          = d->left.flags;
 452             d->left.flags  = d->right.flags;
 453             d->right.flags = flags;
 454             d->cmpt        = reverse_comparison_type(d->cmpt);
 455         }
 456         else if (!(d->left.flags & CMP_DYNAMICVAL))
 457         {
 458             convert_int_real(n1, &d->left);
 459         }
 460         else /* d->right is static */
 461         {
 462             convert_real_int(n2, &d->right, d->cmpt, true);
 463         }
 464     }
 465 }
 466
 467 /*!
 468  * \param data Data to free (should point to a \c t_methoddata_compare).
 469  *
 470  * Frees the memory allocated for \c t_methoddata_compare.
 471  */
 472 static void
 473 free_data_compare(void *data)
 474 {
 475     t_methoddata_compare *d = (t_methoddata_compare *)data;
 476
 477     sfree(d->cmpop);
 478     if (d->left.flags & CMP_ALLOCINT)
 479     {
 480         sfree(d->left.i);
 481     }
 482     if (d->left.flags & CMP_ALLOCREAL)
 483     {
 484         sfree(d->left.r);
 485     }
 486     if (d->right.flags & CMP_ALLOCINT)
 487     {
 488         sfree(d->right.i);
 489     }
 490     if (d->right.flags & CMP_ALLOCREAL)
 491     {
 492         sfree(d->right.r);
 493     }
 494     sfree(d);
 495 }
 496
 497 /*!
 498  * \param[in]  top   Not used.
 499  * \param[in]  fr    Not used.
 500  * \param[in]  pbc   Not used.
 501  * \param[in]  g     Evaluation index group.
 502  * \param[out] out   Output data structure (\p out->u.g is used).
 503  * \param[in]  data  Should point to a \c t_methoddata_compare.
 504  */
 505 static void
 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)
 508 {
 509     t_methoddata_compare *d = (t_methoddata_compare *)data;
 510     int                   i, i1, i2, ig;
 511     int                   a, b;
 512     bool                  bAccept;
 513
 514     for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
 515     {
 516         a       = d->left.i[i1];
 517         b       = d->right.i[i2];
 518         bAccept = false;
 519         switch (d->cmpt)
 520         {
 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;
 528         }
 529         if (bAccept)
 530         {
 531             out->u.g->index[ig++] = g->index[i];
 532         }
 533         if (!(d->left.flags & CMP_SINGLEVAL))
 534         {
 535             ++i1;
 536         }
 537         if (!(d->right.flags & CMP_SINGLEVAL))
 538         {
 539             ++i2;
 540         }
 541     }
 542     out->u.g->isize = ig;
 543 }
 544
 545 /*!
 546  * \param[in]  top   Not used.
 547  * \param[in]  fr    Not used.
 548  * \param[in]  pbc   Not used.
 549  * \param[in]  g     Evaluation index group.
 550  * \param[out] out   Output data structure (\p out->u.g is used).
 551  * \param[in]  data  Should point to a \c t_methoddata_compare.
 552  */
 553 static void
 554 evaluate_compare_real(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 555                       gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
 556 {
 557     t_methoddata_compare *d = (t_methoddata_compare *)data;
 558     int                   i, i1, i2, ig;
 559     real                  a, b;
 560     bool                  bAccept;
 561
 562     for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
 563     {
 564         a       = d->left.r[i1];
 565         b       = (d->right.flags & CMP_REALVAL) ? d->right.r[i2] : d->right.i[i2];
 566         bAccept = false;
 567         switch (d->cmpt)
 568         {
 569             case CMP_INVALID: break;
 570             case CMP_LESS:    bAccept = a <  b; break;
 571             case CMP_LEQ:     bAccept = a <= b; break;
 572             case CMP_GTR:     bAccept = a >  b; break;
 573             case CMP_GEQ:     bAccept = a >= b; break;
 574             case CMP_EQUAL:   bAccept =  gmx_within_tol(a, b, GMX_REAL_EPS); break;
 575             case CMP_NEQ:     bAccept = !gmx_within_tol(a, b, GMX_REAL_EPS); break;
 576         }
 577         if (bAccept)
 578         {
 579             out->u.g->index[ig++] = g->index[i];
 580         }
 581         if (!(d->left.flags & CMP_SINGLEVAL))
 582         {
 583             ++i1;
 584         }
 585         if (!(d->right.flags & CMP_SINGLEVAL))
 586         {
 587             ++i2;
 588         }
 589     }
 590     out->u.g->isize = ig;
 591 }
 592
 593 /*!
 594  * \param[in]  top   Not used.
 595  * \param[in]  fr    Not used.
 596  * \param[in]  pbc   Not used.
 597  * \param[in]  g     Evaluation index group.
 598  * \param[out] out   Output data structure (\p out->u.g is used).
 599  * \param[in]  data  Should point to a \c t_methoddata_compare.
 600  */
 601 static void
 602 evaluate_compare(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 603                  gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
 604 {
 605     t_methoddata_compare *d = (t_methoddata_compare *)data;
 606
 607     if (!((d->left.flags | d->right.flags) & CMP_REALVAL))
 608     {
 609         evaluate_compare_int(top, fr, pbc, g, out, data);
 610     }
 611     else
 612     {
 613         evaluate_compare_real(top, fr, pbc, g, out, data);
 614     }
 615 }