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,2019, 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* 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 init_compare(const gmx_mtop_t* top, int npar, gmx_ana_selparam_t* param, void* data);
 127 /** Frees the memory allocated for comparison expression evaluation. */
 128 static void free_data_compare(void* data);
 129 /*! \brief
 130  * Evaluates comparison expressions.
 131  *
 132  * \param[in]  context  Not used.
 133  * \param[in]  g        Evaluation index group.
 134  * \param[out] out      Output data structure (\p out->u.g is used).
 135  * \param[in]  data     Should point to a \c t_methoddata_compare.
 136  */
 137 static void evaluate_compare(const gmx::SelMethodEvalContext& context,
 138                              gmx_ana_index_t*                 g,
 139                              gmx_ana_selvalue_t*              out,
 140                              void*                            data);
 141
 142 /** Parameters for comparison expression evaluation. */
 143 static gmx_ana_selparam_t smparams_compare[] = {
 144     { "int1", { INT_VALUE, -1, { nullptr } }, nullptr, SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL },
 145     { "real1", { REAL_VALUE, -1, { nullptr } }, nullptr, SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL },
 146     { "op", { STR_VALUE, 1, { nullptr } }, nullptr, 0 },
 147     { "int2", { INT_VALUE, -1, { nullptr } }, nullptr, SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL },
 148     { "real2", { REAL_VALUE, -1, { nullptr } }, nullptr, SPAR_OPTIONAL | SPAR_DYNAMIC | SPAR_ATOMVAL },
 149 };
 150
 151 /** \internal Selection method data for comparison expression evaluation. */
 152 gmx_ana_selmethod_t sm_compare = {
 153     "cmp",
 154     GROUP_VALUE,
 155     SMETH_SINGLEVAL,
 156     asize(smparams_compare),
 157     smparams_compare,
 158     &init_data_compare,
 159     nullptr,
 160     &init_compare,
 161     nullptr,
 162     &free_data_compare,
 163     nullptr,
 164     &evaluate_compare,
 165     nullptr,
 166     { nullptr, nullptr, 0, nullptr },
 167 };
 168
 169 /*! \brief
 170  * Returns a \c e_comparison_t value corresponding to an operator.
 171  *
 172  * \param[in] str  String to process.
 173  * \returns   The comparison type corresponding to the first one or two
 174  *   characters of \p str.
 175  *
 176  * \p str can contain any number of characters; only the first two
 177  * are used.
 178  * If the beginning of \p str does not match any of the recognized types,
 179  * \ref CMP_INVALID is returned.
 180  */
 181 static e_comparison_t comparison_type(const char* str)
 182 {
 183     switch (str[0])
 184     {
 185         case '<': return (str[1] == '=') ? CMP_LEQ : CMP_LESS;
 186         case '>': return (str[1] == '=') ? CMP_GEQ : CMP_GTR;
 187         case '=': return (str[1] == '=') ? CMP_EQUAL : CMP_INVALID;
 188         case '!': return (str[1] == '=') ? CMP_NEQ : CMP_INVALID;
 189     }
 190     return CMP_INVALID;
 191 }
 192
 193 /*! \brief
 194  * Returns a string corresponding to a \c e_comparison_t value.
 195  *
 196  * \param[in] cmpt  Comparison type to convert.
 197  * \returns   Pointer to a string that corresponds to \p cmpt.
 198  *
 199  * The return value points to a string constant and should not be \p free'd.
 200  *
 201  * The function returns NULL if \p cmpt is not one of the valid values.
 202  */
 203 static const char* comparison_type_str(e_comparison_t cmpt)
 204 {
 205     const char* p = nullptr;
 206     switch (cmpt)
 207     {
 208         case CMP_INVALID: p = "INVALID"; break;
 209         case CMP_LESS: p = "<"; break;
 210         case CMP_LEQ: p = "<="; break;
 211         case CMP_GTR: p = ">"; break;
 212         case CMP_GEQ: p = ">="; break;
 213         case CMP_EQUAL: p = "=="; break;
 214         case CMP_NEQ:
 215             p = "!=";
 216             break;
 217             // No default clause so we intentionally get compiler errors
 218             // if new selection choices are added later.
 219     }
 220     return p;
 221 }
 222
 223 /*!
 224  * \param[in] fp    File to receive the output.
 225  * \param[in] data  Should point to a \c t_methoddata_compare.
 226  */
 227 void _gmx_selelem_print_compare_info(FILE* fp, void* data)
 228 {
 229     t_methoddata_compare* d = static_cast<t_methoddata_compare*>(data);
 230
 231     fprintf(fp, " \"");
 232     /* Print the left value */
 233     if ((d->left.flags & CMP_SINGLEVAL) && !(d->left.flags & CMP_DYNAMICVAL))
 234     {
 235         if (d->left.flags & CMP_REALVAL)
 236         {
 237             fprintf(fp, "%f ", d->left.r[0]);
 238         }
 239         else
 240         {
 241             fprintf(fp, "%d ", d->left.i[0]);
 242         }
 243     }
 244     /* Print the operator */
 245     if (d->cmpt != CMP_INVALID)
 246     {
 247         fprintf(fp, "%s", comparison_type_str(d->cmpt));
 248     }
 249     else
 250     {
 251         fprintf(fp, "%s", d->cmpop);
 252     }
 253     /* Print the right value */
 254     if ((d->right.flags & CMP_SINGLEVAL) && !(d->right.flags & CMP_DYNAMICVAL))
 255     {
 256         if (d->right.flags & CMP_REALVAL)
 257         {
 258             fprintf(fp, " %f", d->right.r[0]);
 259         }
 260         else
 261         {
 262             fprintf(fp, " %d", d->right.i[0]);
 263         }
 264     }
 265     fprintf(fp, "\"");
 266 }
 267
 268 static void* init_data_compare(int /* npar */, gmx_ana_selparam_t* param)
 269 {
 270     t_methoddata_compare* data;
 271
 272     snew(data, 1);
 273     param[2].val.u.s = &data->cmpop;
 274     return data;
 275 }
 276
 277 /*! \brief
 278  * Reverses a comparison operator.
 279  *
 280  * \param[in] type  Comparison operator to reverse.
 281  * \returns   The correct comparison operator that equals \p type when the
 282  *   left and right sides are interchanged.
 283  */
 284 static e_comparison_t reverse_comparison_type(e_comparison_t type)
 285 {
 286     switch (type)
 287     {
 288         case CMP_LESS: return CMP_GTR;
 289         case CMP_LEQ: return CMP_GEQ;
 290         case CMP_GTR: return CMP_LESS;
 291         case CMP_GEQ: return CMP_LEQ;
 292         default: break;
 293     }
 294     return type;
 295 }
 296
 297 /*! \brief
 298  * Initializes the value storage for comparison expression.
 299  *
 300  * \param[out] val   Value structure to initialize.
 301  * \param[in]  param Parameters to use for initialization.
 302  * \returns    The number of values provided for the value, 0 on error.
 303  */
 304 static int init_comparison_value(t_compare_value* val, gmx_ana_selparam_t param[2])
 305 {
 306     int n;
 307
 308     val->flags = 0;
 309     if (param[0].flags & SPAR_SET)
 310     {
 311         val->flags |= (param[0].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
 312         val->flags |= !(param[0].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL : 0;
 313         n      = param[0].val.nr;
 314         val->i = param[0].val.u.i;
 315     }
 316     else if (param[1].flags & SPAR_SET)
 317     {
 318         val->flags |= (param[1].flags & SPAR_DYNAMIC) ? CMP_DYNAMICVAL : 0;
 319         val->flags |= !(param[1].flags & SPAR_ATOMVAL) ? CMP_SINGLEVAL : 0;
 320         val->flags |= CMP_REALVAL;
 321         n      = param[1].val.nr;
 322         val->r = param[1].val.u.r;
 323     }
 324     else
 325     {
 326         n      = 0;
 327         val->i = nullptr;
 328         val->r = nullptr;
 329     }
 330     return n;
 331 }
 332
 333 /*! \brief
 334  * Converts an integer value to floating point.
 335  *
 336  * \param[in]     n   Number of values in the \p val->u array.
 337  * \param[in,out] val Value to convert.
 338  */
 339 static void 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] = static_cast<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 convert_real_int(int n, t_compare_value* val, e_comparison_t cmpt, bool bRight)
 368 {
 369     int  i;
 370     int* iv;
 371
 372     if (!bRight)
 373     {
 374         cmpt = reverse_comparison_type(cmpt);
 375     }
 376     snew(iv, n);
 377     /* Round according to the comparison type */
 378     for (i = 0; i < n; ++i)
 379     {
 380         switch (cmpt)
 381         {
 382             case CMP_LESS:
 383             case CMP_GEQ: iv[i] = static_cast<int>(std::ceil(val->r[i])); break;
 384             case CMP_GTR:
 385             case CMP_LEQ: iv[i] = static_cast<int>(std::floor(val->r[i])); 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(
 394                         gmx::NotImplementedError("Equality comparison between dynamic integer and "
 395                                                  "static real expressions not implemented"));
 396             case CMP_INVALID: /* Should not be reached */
 397                 sfree(iv);
 398                 GMX_THROW(gmx::InternalError("Invalid comparison type"));
 399         }
 400     }
 401     /* Free the previous value if one is present. */
 402     sfree(val->i);
 403     val->i = iv;
 404     val->flags &= ~CMP_REALVAL;
 405     val->flags |= CMP_ALLOCINT;
 406 }
 407
 408 static void init_compare(const gmx_mtop_t* /* top */, int /* npar */, gmx_ana_selparam_t* param, void* data)
 409 {
 410     t_methoddata_compare* d = static_cast<t_methoddata_compare*>(data);
 411     int                   n1, n2;
 412
 413     /* Store the values */
 414     n1 = init_comparison_value(&d->left, &param[0]);
 415     n2 = init_comparison_value(&d->right, &param[3]);
 416     /* Store the comparison type */
 417     d->cmpt = comparison_type(d->cmpop);
 418     if (d->cmpt == CMP_INVALID)
 419     {
 420         GMX_THROW(gmx::InternalError("Invalid comparison type"));
 421     }
 422     /* Convert the values to the same type */
 423     /* TODO: Currently, there are no dynamic integer-valued selection methods,
 424      * which means that only the branches with convert_int_real() will ever be
 425      * taken. It should be considered whether it is necessary to support these
 426      * other cases at all.
 427      */
 428     if ((d->left.flags & CMP_REALVAL) && !(d->right.flags & CMP_REALVAL))
 429     {
 430         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 431         {
 432             /* Nothing can be done */
 433         }
 434         else if (!(d->right.flags & CMP_DYNAMICVAL))
 435         {
 436             convert_int_real(n2, &d->right);
 437         }
 438         else /* d->left is static */
 439         {
 440             convert_real_int(n1, &d->left, d->cmpt, false);
 441         }
 442     }
 443     else if (!(d->left.flags & CMP_REALVAL) && (d->right.flags & CMP_REALVAL))
 444     {
 445         if (d->left.flags & d->right.flags & CMP_DYNAMICVAL)
 446         {
 447             /* Reverse the sides to place the integer on the right */
 448             int flags;
 449             d->left.r      = d->right.r;
 450             d->right.r     = nullptr;
 451             d->right.i     = d->left.i;
 452             d->left.i      = nullptr;
 453             flags          = d->left.flags;
 454             d->left.flags  = d->right.flags;
 455             d->right.flags = flags;
 456             d->cmpt        = reverse_comparison_type(d->cmpt);
 457         }
 458         else if (!(d->left.flags & CMP_DYNAMICVAL))
 459         {
 460             convert_int_real(n1, &d->left);
 461         }
 462         else /* d->right is static */
 463         {
 464             convert_real_int(n2, &d->right, d->cmpt, true);
 465         }
 466     }
 467 }
 468
 469 /*!
 470  * \param data Data to free (should point to a \c t_methoddata_compare).
 471  *
 472  * Frees the memory allocated for \c t_methoddata_compare.
 473  */
 474 static void free_data_compare(void* data)
 475 {
 476     t_methoddata_compare* d = static_cast<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]  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 evaluate_compare_int(gmx_ana_index_t* g, gmx_ana_selvalue_t* out, void* data)
 506 {
 507     t_methoddata_compare* d = static_cast<t_methoddata_compare*>(data);
 508     int                   i, i1, i2, ig;
 509     int                   a, b;
 510     bool                  bAccept;
 511
 512     for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
 513     {
 514         a       = d->left.i[i1];
 515         b       = d->right.i[i2];
 516         bAccept = false;
 517         switch (d->cmpt)
 518         {
 519             case CMP_INVALID: break;
 520             case CMP_LESS: bAccept = a < b; break;
 521             case CMP_LEQ: bAccept = a <= b; break;
 522             case CMP_GTR: bAccept = a > b; break;
 523             case CMP_GEQ: bAccept = a >= b; break;
 524             case CMP_EQUAL: bAccept = a == b; break;
 525             case CMP_NEQ: bAccept = a != b; break;
 526         }
 527         if (bAccept)
 528         {
 529             out->u.g->index[ig++] = g->index[i];
 530         }
 531         if (!(d->left.flags & CMP_SINGLEVAL))
 532         {
 533             ++i1;
 534         }
 535         if (!(d->right.flags & CMP_SINGLEVAL))
 536         {
 537             ++i2;
 538         }
 539     }
 540     out->u.g->isize = ig;
 541 }
 542
 543 /*! \brief
 544  * Implementation for evaluate_compare() if either value is non-integer.
 545  *
 546  * \param[in]  g     Evaluation index group.
 547  * \param[out] out   Output data structure (\p out->u.g is used).
 548  * \param[in]  data  Should point to a \c t_methoddata_compare.
 549  *
 550  * Left value is assumed to be real-valued; right value can be either.
 551  * This is ensured by the initialization method.
 552  */
 553 static void evaluate_compare_real(gmx_ana_index_t* g, gmx_ana_selvalue_t* out, void* data)
 554 {
 555     t_methoddata_compare* d = static_cast<t_methoddata_compare*>(data);
 556     int                   i, i1, i2, ig;
 557     real                  a, b;
 558     bool                  bAccept;
 559
 560     for (i = i1 = i2 = ig = 0; i < g->isize; ++i)
 561     {
 562         a       = d->left.r[i1];
 563         b       = (d->right.flags & CMP_REALVAL) ? d->right.r[i2] : d->right.i[i2];
 564         bAccept = false;
 565         switch (d->cmpt)
 566         {
 567             case CMP_INVALID: break;
 568             case CMP_LESS: bAccept = a < b; break;
 569             case CMP_LEQ: bAccept = a <= b; break;
 570             case CMP_GTR: bAccept = a > b; break;
 571             case CMP_GEQ: bAccept = a >= b; break;
 572             case CMP_EQUAL: bAccept = gmx_within_tol(a, b, GMX_REAL_EPS); break;
 573             case CMP_NEQ: bAccept = !gmx_within_tol(a, b, GMX_REAL_EPS); break;
 574         }
 575         if (bAccept)
 576         {
 577             out->u.g->index[ig++] = g->index[i];
 578         }
 579         if (!(d->left.flags & CMP_SINGLEVAL))
 580         {
 581             ++i1;
 582         }
 583         if (!(d->right.flags & CMP_SINGLEVAL))
 584         {
 585             ++i2;
 586         }
 587     }
 588     out->u.g->isize = ig;
 589 }
 590
 591 static void evaluate_compare(const gmx::SelMethodEvalContext& /*context*/,
 592                              gmx_ana_index_t*    g,
 593                              gmx_ana_selvalue_t* out,
 594                              void*               data)
 595 {
 596     t_methoddata_compare* d = static_cast<t_methoddata_compare*>(data);
 597
 598     if (!((d->left.flags | d->right.flags) & CMP_REALVAL))
 599     {
 600         evaluate_compare_int(g, out, data);
 601     }
 602     else
 603     {
 604         evaluate_compare_real(g, out, data);
 605     }
 606 }