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