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