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.
   5  * Copyright (c) 2014,2015,2016,2017,2018 by the GROMACS development team.
   6  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   7  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   8  * and including many others, as listed in the AUTHORS file in the
   9  * top-level source directory and at http://www.gromacs.org.
  10  *
  11  * GROMACS is free software; you can redistribute it and/or
  12  * modify it under the terms of the GNU Lesser General Public License
  13  * as published by the Free Software Foundation; either version 2.1
  14  * of the License, or (at your option) any later version.
  15  *
  16  * GROMACS is distributed in the hope that it will be useful,
  17  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  18  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  19  * Lesser General Public License for more details.
  20  *
  21  * You should have received a copy of the GNU Lesser General Public
  22  * License along with GROMACS; if not, see
  23  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  24  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  25  *
  26  * If you want to redistribute modifications to GROMACS, please
  27  * consider that scientific software is very special. Version
  28  * control is crucial - bugs must be traceable. We will be happy to
  29  * consider code for inclusion in the official distribution, but
  30  * derived work must not be called official GROMACS. Details are found
  31  * in the README & COPYING files - if they are missing, get the
  32  * official version at http://www.gromacs.org.
  33  *
  34  * To help us fund GROMACS development, we humbly ask that you cite
  35  * the research papers on the package. Check out http://www.gromacs.org.
  36  */
  37 /*! \internal \file
  38  * \brief
  39  * Implements internal selection method for comparison expressions.
  40  *
  41  * \author Teemu Murtola <teemu.murtola@gmail.com>
  42  * \ingroup module_selection
  43  */
  44 #include "gmxpre.h"
  45
  46 #include <cmath>
  47
  48 #include "gromacs/math/utilities.h"
  49 #include "gromacs/utility/arraysize.h"
  50 #include "gromacs/utility/basedefinitions.h"
  51 #include "gromacs/utility/exceptions.h"
  52 #include "gromacs/utility/smalloc.h"
  53
  54 #include "keywords.h"
  55 #include "selmethod.h"
  56
  57 /** Defines the comparison operator for comparison expressions. */
  58 typedef enum
  59 {
  60     CMP_INVALID, /**< Indicates an error */
  61     CMP_LESS,    /**< '<' */
  62     CMP_LEQ,     /**< '<=' */
  63     CMP_GTR,     /**< '>' */
  64     CMP_GEQ,     /**< '>=' */
  65     CMP_EQUAL,   /**< '==' */
  66     CMP_NEQ      /**< '!=' */
  67 } e_comparison_t;
  68
  69 /** The operand has a single value. */
  70 #define CMP_SINGLEVAL 1
  71 /** The operand value is dynamic. */
  72 #define CMP_DYNAMICVAL 2
  73 /** The value is real. */
  74 #define CMP_REALVAL 4
  75 /** The integer array is allocated. */
  76 #define CMP_ALLOCINT 16
  77 /** The real array is allocated. */
  78 #define CMP_ALLOCREAL 32
  79
  80 /*! \internal \brief
  81  * Data structure for comparison expression operand values.
  82  */
  83 typedef struct
  84 {
  85     /** Flags that describe the type of the operand. */
  86     int flags;
  87     /** (Array of) integer value(s). */
  88     int* i;
  89     /** (Array of) real value(s). */
  90     real* r;
  91 } t_compare_value;
  92
  93 /*! \internal \brief
  94  * Data structure for comparison expression evaluation.
  95  */
  96 typedef struct
  97 {
  98     /** Comparison operator as a string. */
  99     char* cmpop;
 100     /** Comparison operator type. */
 101     e_comparison_t cmpt;
 102     /** Left value. */
 103     t_compare_value left;
 104     /** Right value. */
 105     t_compare_value right;
 106 } t_methoddata_compare;
 107
 108 /*! \brief
 109  * Allocates data for comparison expression evaluation.
 110  *
 111  * \param[in]     npar  Not used (should be 5).
 112  * \param[in,out] param Method parameters (should point to a copy of
 113  *   \ref smparams_compare).
 114  * \returns       Pointer to the allocated data (\c t_methoddata_compare).
 115  *
 116  * Allocates memory for a \c t_methoddata_compare structure.
 117  */
 118 static void* init_data_compare(int npar, gmx_ana_selparam_t* param);
 119 /*! \brief
 120  * Initializes data for comparison expression evaluation.
 121  *
 122  * \param[in] top   Not used.
 123  * \param[in] npar  Not used (should be 5).
 124  * \param[in] param Method parameters (should point to \ref smparams_compare).
 125  * \param[in] data  Should point to a \c t_methoddata_compare.
 126  * \returns   0 if the input data is valid, -1 on error.
 127  */
 128 static void 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 free_data_compare(void* data);
 131 /*! \brief
 132  * Evaluates comparison expressions.
 133  *
 134  * \param[in]  context  Not used.
 135  * \param[in]  g        Evaluation index group.
 136  * \param[out] out      Output data structure (\p out->u.g is used).
 137  * \param[in]  data     Should point to a \c t_methoddata_compare.
 138  */
 139 static void evaluate_compare(const gmx::SelMethodEvalContext& context,
 140                              gmx_ana_index_t*                 g,
 141                              gmx_ana_selvalue_t*              out,
 142                              void*                            data);
 143
 144 /** Parameters for comparison expression evaluation. */
 145 static gmx_ana_selparam_t smparams_compare[] = {
 146     { "int1", { INT_VALUE, -1, { nullptr } }, nullptr, SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL },
 147     { "real1", { REAL_VALUE, -1, { nullptr } }, nullptr, SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL },
 148     { "op", { STR_VALUE, 1, { nullptr } }, nullptr, 0 },
 149     { "int2", { INT_VALUE, -1, { nullptr } }, nullptr, SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL },
 150     { "real2", { REAL_VALUE, -1, { nullptr } }, nullptr, SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL },
 151 };
 152
 153 /** \internal Selection method data for comparison expression evaluation. */
 154 gmx_ana_selmethod_t sm_compare = {
 155     "cmp",
 156     GROUP_VALUE,
 157     SMETH_SINGLEVAL,
 158     asize(smparams_compare),
 159     smparams_compare,
 160     &init_data_compare,
 161     nullptr,
 162     &init_compare,
 163     nullptr,
 164     &free_data_compare,
 165     nullptr,
 166     &evaluate_compare,
 167     nullptr,
 168     { nullptr, nullptr, 0, nullptr },
 169 };
 170
 171 /*! \brief
 172  * Returns a \c e_comparison_t value corresponding to an operator.
 173  *
 174  * \param[in] str  String to process.
 175  * \returns   The comparison type corresponding to the first one or two
 176  *   characters of \p str.
 177  *
 178  * \p str can contain any number of characters; only the first two
 179  * are used.
 180  * If the beginning of \p str does not match any of the recognized types,
 181  * \ref CMP_INVALID is returned.
 182  */
 183 static e_comparison_t comparison_type(const char* str)
 184 {
 185     switch (str[0])
 186     {
 187         case '<': return (str[1] == '=') ? CMP_LEQ : CMP_LESS;
 188         case '>': return (str[1] == '=') ? CMP_GEQ : CMP_GTR;
 189         case '=': return (str[1] == '=') ? CMP_EQUAL : CMP_INVALID;
 190         case '!': return (str[1] == '=') ? CMP_NEQ : CMP_INVALID;
 191     }
 192     return CMP_INVALID;
 193 }
 194
 195 /*! \brief
 196  * Returns a string corresponding to a \c e_comparison_t value.
 197  *
 198  * \param[in] cmpt  Comparison type to convert.
 199  * \returns   Pointer to a string that corresponds to \p cmpt.
 200  *
 201  * The return value points to a string constant and should not be \p free'd.
 202  *
 203  * The function returns NULL if \p cmpt is not one of the valid values.
 204  */
 205 static const char* comparison_type_str(e_comparison_t cmpt)
 206 {
 207     const char* p = nullptr;
 208     switch (cmpt)
 209     {
 210         case CMP_INVALID: p = "INVALID"; break;
 211         case CMP_LESS: p = "<"; break;
 212         case CMP_LEQ: p = "<="; break;
 213         case CMP_GTR: p = ">"; break;
 214         case CMP_GEQ: p = ">="; break;
 215         case CMP_EQUAL: p = "=="; break;
 216         case CMP_NEQ:
 217             p = "!=";
 218             break;
 219             // No default clause so we intentionally get compiler errors
 220             // if new selection choices are added later.
 221     }
 222     return p;
 223 }
 224
 225 /*!
 226  * \param[in] fp    File to receive the output.
 227  * \param[in] data  Should point to a \c t_methoddata_compare.
 228  */
 229 void _gmx_selelem_print_compare_info(FILE* fp, void* data)
 230 {
 231     t_methoddata_compare* d = static_cast<t_methoddata_compare*>(data);
 232
 233     fprintf(fp, " \"");
 234     /* Print the left value */
 235     if ((d->left.flags & CMP_SINGLEVAL) && !(d->left.flags & CMP_DYNAMICVAL))
 236     {
 237         if (d->left.flags & CMP_REALVAL)
 238         {
 239             fprintf(fp, "%f ", d->left.r[0]);
 240         }
 241         else
 242         {
 243             fprintf(fp, "%d ", d->left.i[0]);
 244         }
 245     }
 246     /* Print the operator */
 247     if (d->cmpt != CMP_INVALID)
 248     {
 249         fprintf(fp, "%s", comparison_type_str(d->cmpt));
 250     }
 251     else
 252     {
 253         fprintf(fp, "%s", d->cmpop);
 254     }
 255     /* Print the right value */
 256     if ((d->right.flags & CMP_SINGLEVAL) && !(d->right.flags & CMP_DYNAMICVAL))
 257     {
 258         if (d->right.flags & CMP_REALVAL)
 259         {
 260             fprintf(fp, " %f", d->right.r[0]);
 261         }
 262         else
 263         {
 264             fprintf(fp, " %d", d->right.i[0]);
 265         }
 266     }
 267     fprintf(fp, "\"");
 268 }
 269
 270 static void* init_data_compare(int /* npar */, gmx_ana_selparam_t* param)
 271 {
 272     t_methoddata_compare* data;
 273
 274     snew(data, 1);
 275     param[2].val.u.s = &data->cmpop;
 276     return data;
 277 }
 278
 279 /*! \brief
 280  * Reverses a comparison operator.
 281  *
 282  * \param[in] type  Comparison operator to reverse.
 283  * \returns   The correct comparison operator that equals \p type when the
 284  *   left and right sides are interchanged.
 285  */
 286 static e_comparison_t reverse_comparison_type(e_comparison_t type)
 287 {
 288     switch (type)
 289     {
 290         case CMP_LESS: return CMP_GTR;
 291         case CMP_LEQ: return CMP_GEQ;
 292         case CMP_GTR: return CMP_LESS;
 293         case CMP_GEQ: return CMP_LEQ;
 294         default: break;
 295     }
 296     return type;
 297 }
 298
 299 /*! \brief
 300  * Initializes the value storage for comparison expression.
 301  *
 302  * \param[out] val   Value structure to initialize.
 303  * \param[in]  param Parameters to use for initialization.
 304  * \returns    The number of values provided for the value, 0 on error.
 305  */
 306 static int init_comparison_value(t_compare_value* val, gmx_ana_selparam_t param[2])
 307 {
 308     int n;
 309
 310     val->flags = 0;
 311     if (param[0].flags & SPAR_SET)
 312     {
 313         val->flags |= (param[0].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
 314         val->flags |= !(param[0].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL : 0;
 315         n      = param[0].val.nr;
 316         val->i = param[0].val.u.i;
 317     }
 318     else if (param[1].flags & SPAR_SET)
 319     {
 320         val->flags |= (param[1].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
 321         val->flags |= !(param[1].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL : 0;
 322         val->flags |= CMP_REALVAL;
 323         n      = param[1].val.nr;
 324         val->r = param[1].val.u.r;
 325     }
 326     else
 327     {
 328         n      = 0;
 329         val->i = nullptr;
 330         val->r = nullptr;
 331     }
 332     return n;
 333 }
 334
 335 /*! \brief
 336  * Converts an integer value to floating point.
 337  *
 338  * \param[in]     n   Number of values in the \p val->u array.
 339  * \param[in,out] val Value to convert.
 340  */
 341 static void convert_int_real(int n, t_compare_value* val)
 342 {
 343     int   i;
 344     real* rv;
 345
 346     snew(rv, n);
 347     for (i = 0; i < n; ++i)
 348     {
 349         rv[i] = static_cast<real>(val->i[i]);
 350     }
 351     /* Free the previous value if one is present. */
 352     sfree(val->r);
 353     val->r = rv;
 354     val->flags |= CMP_REALVAL | CMP_ALLOCREAL;
 355 }
 356
 357 /*! \brief
 358  * Converts a floating point value to integer.
 359  *
 360  * \param[in]     n      Number of values in the \p val->u array.
 361  * \param[in,out] val    Value to convert.
 362  * \param[in]     cmpt   Comparison operator type.
 363  * \param[in]     bRight true if \p val appears on the right hand size of
 364  *   \p cmpt.
 365  * \returns       0 on success, EINVAL on error.
 366  *
 367  * The values are rounded such that the same comparison operator can be used.
 368  */
 369 static void convert_real_int(int n, t_compare_value* val, e_comparison_t cmpt, bool bRight)
 370 {
 371     int  i;
 372     int* iv;
 373
 374     if (!bRight)
 375     {
 376         cmpt = reverse_comparison_type(cmpt);
 377     }
 378     snew(iv, n);
 379     /* Round according to the comparison type */
 380     for (i = 0; i < n; ++i)
 381     {
 382         switch (cmpt)
 383         {
 384             case CMP_LESS:
 385             case CMP_GEQ: iv[i] = static_cast<int>(std::ceil(val->r[i])); break;
 386             case CMP_GTR:
 387             case CMP_LEQ: iv[i] = static_cast<int>(std::floor(val->r[i])); break;
 388             case CMP_EQUAL:
 389             case CMP_NEQ:
 390                 sfree(iv);
 391                 /* TODO: Implement, although it isn't typically very useful.
 392                  * Implementation is only a matter of proper initialization,
 393                  * the evaluation function can already handle this case with
 394                  * proper preparations. */
 395                 GMX_THROW(
 396                         gmx::NotImplementedError("Equality comparison between dynamic integer and "
 397                                                  "static real expressions not implemented"));
 398             case CMP_INVALID: /* Should not be reached */
 399                 sfree(iv);
 400                 GMX_THROW(gmx::InternalError("Invalid comparison type"));
 401         }
 402     }
 403     /* Free the previous value if one is present. */
 404     sfree(val->i);
 405     val->i = iv;
 406     val->flags &= ~CMP_REALVAL;
 407     val->flags |= CMP_ALLOCINT;
 408 }
 409
 410 static void init_compare(const gmx_mtop_t* /* top */, int /* npar */, gmx_ana_selparam_t* param, void* data)
 411 {
 412     t_methoddata_compare* d = static_cast<t_methoddata_compare*>(data);
 413     int                   n1, n2;
 414
 415     /* Store the values */
 416     n1 = init_comparison_value(&d->left, &param[0]);
 417     n2 = init_comparison_value(&d->right, &param[3]);
 418     /* Store the comparison type */
 419     d->cmpt = comparison_type(d->cmpop);
 420     if (d->cmpt == CMP_INVALID)
 421     {
 422         GMX_THROW(gmx::InternalError("Invalid comparison type"));
 423     }
 424     /* Convert the values to the same type */
 425     /* TODO: Currently, there are no dynamic integer-valued selection methods,
 426      * which means that only the branches with convert_int_real() will ever be
 427      * taken. It should be considered whether it is necessary to support these
 428      * other cases at all.
 429      */
 430     if ((d->left.flags & CMP_REALVAL) && !(d->right.flags & CMP_REALVAL))
 431     {
 432         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 433         {
 434             /* Nothing can be done */
 435         }
 436         else if (!(d->right.flags & CMP_DYNAMICVAL))
 437         {
 438             convert_int_real(n2, &d->right);
 439         }
 440         else /* d->left is static */
 441         {
 442             convert_real_int(n1, &d->left, d->cmpt, false);
 443         }
 444     }
 445     else if (!(d->left.flags & CMP_REALVAL) && (d->right.flags & CMP_REALVAL))
 446     {
 447         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 448         {
 449             /* Reverse the sides to place the integer on the right */
 450             int flags;
 451             d->left.r      = d->right.r;
 452             d->right.r     = nullptr;
 453             d->right.i     = d->left.i;
 454             d->left.i      = nullptr;
 455             flags          = d->left.flags;
 456             d->left.flags  = d->right.flags;
 457             d->right.flags = flags;
 458             d->cmpt        = reverse_comparison_type(d->cmpt);
 459         }
 460         else if (!(d->left.flags & CMP_DYNAMICVAL))
 461         {
 462             convert_int_real(n1, &d->left);
 463         }
 464         else /* d->right is static */
 465         {
 466             convert_real_int(n2, &d->right, d->cmpt, true);
 467         }
 468     }
 469 }
 470
 471 /*!
 472  * \param data Data to free (should point to a \c t_methoddata_compare).
 473  *
 474  * Frees the memory allocated for \c t_methoddata_compare.
 475  */
 476 static void free_data_compare(void* data)
 477 {
 478     t_methoddata_compare* d = static_cast<t_methoddata_compare*>(data);
 479
 480     sfree(d->cmpop);
 481     if (d->left.flags & CMP_ALLOCINT)
 482     {
 483         sfree(d->left.i);
 484     }
 485     if (d->left.flags & CMP_ALLOCREAL)
 486     {
 487         sfree(d->left.r);
 488     }
 489     if (d->right.flags & CMP_ALLOCINT)
 490     {
 491         sfree(d->right.i);
 492     }
 493     if (d->right.flags & CMP_ALLOCREAL)
 494     {
 495         sfree(d->right.r);
 496     }
 497     sfree(d);
 498 }
 499
 500 /*! \brief
 501  * Implementation for evaluate_compare() for integer values.
 502  *
 503  * \param[in]  g     Evaluation index group.
 504  * \param[out] out   Output data structure (\p out->u.g is used).
 505  * \param[in]  data  Should point to a \c t_methoddata_compare.
 506  */
 507 static void evaluate_compare_int(gmx_ana_index_t* g, gmx_ana_selvalue_t* out, void* data)
 508 {
 509     t_methoddata_compare* d = static_cast<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 /*! \brief
 546  * Implementation for evaluate_compare() if either value is non-integer.
 547  *
 548  * \param[in]  g     Evaluation index group.
 549  * \param[out] out   Output data structure (\p out->u.g is used).
 550  * \param[in]  data  Should point to a \c t_methoddata_compare.
 551  *
 552  * Left value is assumed to be real-valued; right value can be either.
 553  * This is ensured by the initialization method.
 554  */
 555 static void evaluate_compare_real(gmx_ana_index_t* g, gmx_ana_selvalue_t* out, void* data)
 556 {
 557     t_methoddata_compare* d = static_cast<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 static void evaluate_compare(const gmx::SelMethodEvalContext& /*context*/,
 594                              gmx_ana_index_t*    g,
 595                              gmx_ana_selvalue_t* out,
 596                              void*               data)
 597 {
 598     t_methoddata_compare* d = static_cast<t_methoddata_compare*>(data);
 599
 600     if (!((d->left.flags | d->right.flags) & CMP_REALVAL))
 601     {
 602         evaluate_compare_int(g, out, data);
 603     }
 604     else
 605     {
 606         evaluate_compare_real(g, out, data);
 607     }
 608 }