src/gromacs/selection/sm_same.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,2021, 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 the \p same selection method.
  40  *
  41  * \author Teemu Murtola <teemu.murtola@gmail.com>
  42  * \ingroup module_selection
  43  */
  44 #include "gmxpre.h"
  45
  46 #include <cstdlib>
  47 #include <cstring>
  48
  49 #include "gromacs/utility/arraysize.h"
  50 #include "gromacs/utility/exceptions.h"
  51 #include "gromacs/utility/smalloc.h"
  52
  53 #include "keywords.h"
  54 #include "parsetree.h"
  55 #include "selelem.h"
  56 #include "selmethod.h"
  57 #include "selmethod_impl.h"
  58
  59 /*! \internal
  60  * \brief
  61  * Data structure for the \p same selection method.
  62  *
  63  * To avoid duplicate initialization code, the same data structure is used
  64  * for matching both integer and string keywords; hence the unions.
  65  *
  66  * \ingroup module_selection
  67  */
  68 typedef struct
  69 {
  70     /** Value for each atom to match. */
  71     union
  72     {
  73         int*   i;
  74         char** s;
  75         void*  ptr;
  76     } val;
  77     /*! \brief
  78      * Number of values in the \p as array.
  79      *
  80      * For string values, this is actually the number of values in the
  81      * \p as_s_sorted array.
  82      */
  83     int nas;
  84     /** Values to match against. */
  85     union
  86     {
  87         int*   i;
  88         char** s;
  89         void*  ptr;
  90     } as;
  91     /*! \brief
  92      * Separate array for sorted \p as.s array.
  93      *
  94      * The array of strings returned as the output value of a parameter should
  95      * not be messed with to avoid memory corruption (the pointers in the array
  96      * may be reused for several evaluations), so we keep our own copy for
  97      * modifications.
  98      */
  99     char** as_s_sorted;
 100     /** Whether simple matching can be used. */
 101     bool bSorted;
 102 } t_methoddata_same;
 103
 104 /*! \brief
 105  * Allocates data for the \p same selection method.
 106  *
 107  * \param[in]     npar  Not used (should be 2).
 108  * \param[in,out] param Method parameters (should point to a copy of
 109  *      ::smparams_same_int or ::smparams_same_str).
 110  * \returns Pointer to the allocated data (\ref t_methoddata_same).
 111  */
 112 static void* init_data_same(int npar, gmx_ana_selparam_t* param);
 113 /*! \brief
 114  * Initializes the \p same selection method.
 115  *
 116  * \param   top   Not used.
 117  * \param   npar  Not used (should be 2).
 118  * \param   param Initialized method parameters (should point to a copy of
 119  *      ::smparams_same_int or ::smparams_same_str).
 120  * \param   data  Pointer to \ref t_methoddata_same to initialize.
 121  * \returns 0 on success, -1 on failure.
 122  */
 123 static void init_same(const gmx_mtop_t* top, int npar, gmx_ana_selparam_t* param, void* data);
 124 /** Frees the data allocated for the \p same selection method. */
 125 static void free_data_same(void* data);
 126 /*! \brief
 127  * Initializes the evaluation of the \p same selection method for a frame.
 128  *
 129  * \param[in]  context  Not used.
 130  * \param      data Should point to a \ref t_methoddata_same.
 131  *
 132  * Sorts the \c data->as.i array and removes identical values for faster and
 133  * simpler lookup.
 134  */
 135 static void init_frame_same_int(const gmx::SelMethodEvalContext& context, void* data);
 136 /** Evaluates the \p same selection method. */
 137 static void evaluate_same_int(const gmx::SelMethodEvalContext& context,
 138                               gmx_ana_index_t*                 g,
 139                               gmx_ana_selvalue_t*              out,
 140                               void*                            data);
 141 /*! \brief
 142  * Initializes the evaluation of the \p same selection method for a frame.
 143  *
 144  * \param[in]  context  Not used.
 145  * \param      data Should point to a \ref t_methoddata_same.
 146  *
 147  * Sorts the \c data->as.s array and removes identical values for faster and
 148  * simpler lookup.
 149  */
 150 static void init_frame_same_str(const gmx::SelMethodEvalContext& context, void* data);
 151 /** Evaluates the \p same selection method. */
 152 static void evaluate_same_str(const gmx::SelMethodEvalContext& context,
 153                               gmx_ana_index_t*                 g,
 154                               gmx_ana_selvalue_t*              out,
 155                               void*                            data);
 156
 157 /** Parameters for the \p same selection method. */
 158 static gmx_ana_selparam_t smparams_same_int[] = {
 159     { nullptr, { INT_VALUE, -1, { nullptr } }, nullptr, SPAR_DYNAMIC | SPAR_ATOMVAL },
 160     { "as", { INT_VALUE, -1, { nullptr } }, nullptr, SPAR_DYNAMIC | SPAR_VARNUM },
 161 };
 162
 163 /** Parameters for the \p same selection method. */
 164 static gmx_ana_selparam_t smparams_same_str[] = {
 165     { nullptr, { STR_VALUE, -1, { nullptr } }, nullptr, SPAR_DYNAMIC | SPAR_ATOMVAL },
 166     { "as", { STR_VALUE, -1, { nullptr } }, nullptr, SPAR_DYNAMIC | SPAR_VARNUM },
 167 };
 168
 169 /** Help text for the \p same selection method. */
 170 static const char* const help_same[] = {
 171     "::",
 172     "",
 173     "  same KEYWORD as ATOM_EXPR",
 174     "",
 175
 176     "The keyword [TT]same[tt] can be used to select all atoms for which",
 177     "the given [TT]KEYWORD[tt] matches any of the atoms in [TT]ATOM_EXPR[tt].",
 178     "Keywords that evaluate to integer or string values are supported.",
 179 };
 180
 181 /*! \internal \brief Selection method data for the \p same method. */
 182 gmx_ana_selmethod_t sm_same = {
 183     "same",
 184     GROUP_VALUE,
 185     0,
 186     asize(smparams_same_int),
 187     smparams_same_int,
 188     &init_data_same,
 189     nullptr,
 190     &init_same,
 191     nullptr,
 192     &free_data_same,
 193     &init_frame_same_int,
 194     &evaluate_same_int,
 195     nullptr,
 196     { "same KEYWORD as ATOM_EXPR", "Extending selections", asize(help_same), help_same },
 197 };
 198
 199 /*! \brief
 200  * Selection method data for the \p same method.
 201  *
 202  * This selection method is used for matching string keywords. The parser
 203  * never sees this method; _gmx_selelem_custom_init_same() replaces sm_same
 204  * with this method in cases where it is required.
 205  */
 206 static gmx_ana_selmethod_t sm_same_str = {
 207     "same",
 208     GROUP_VALUE,
 209     SMETH_SINGLEVAL,
 210     asize(smparams_same_str),
 211     smparams_same_str,
 212     &init_data_same,
 213     nullptr,
 214     &init_same,
 215     nullptr,
 216     &free_data_same,
 217     &init_frame_same_str,
 218     &evaluate_same_str,
 219     nullptr,
 220     { nullptr, nullptr, 0, nullptr },
 221 };
 222
 223 static void* init_data_same(int /* npar */, gmx_ana_selparam_t* param)
 224 {
 225     t_methoddata_same* data;
 226
 227     snew(data, 1);
 228     data->as_s_sorted = nullptr;
 229     param[1].nvalptr  = &data->nas;
 230     return data;
 231 }
 232
 233 /*!
 234  * \param[in,out] method  The method to initialize.
 235  * \param[in,out] params  Pointer to the first parameter.
 236  * \param[in]     scanner Scanner data structure.
 237  *
 238  * If \p *method is not a \c same method, this function returns
 239  * immediately.
 240  */
 241 void _gmx_selelem_custom_init_same(gmx_ana_selmethod_t**                           method,
 242                                    const gmx::SelectionParserParameterListPointer& params,
 243                                    void*                                           scanner)
 244 {
 245
 246     /* Do nothing if this is not a same method. */
 247     if (!*method || (*method)->name != sm_same.name || params->empty())
 248     {
 249         return;
 250     }
 251
 252     const gmx::SelectionParserValueList& kwvalues = params->front().values();
 253     if (kwvalues.size() != 1 || !kwvalues.front().hasExpressionValue()
 254         || kwvalues.front().expr->type != SEL_EXPRESSION)
 255     {
 256         GMX_THROW(gmx::InvalidInputError("'same' should be followed by a single keyword"));
 257     }
 258     gmx_ana_selmethod_t* kwmethod = kwvalues.front().expr->u.expr.method;
 259     if (kwmethod->type == STR_VALUE)
 260     {
 261         *method = &sm_same_str;
 262     }
 263
 264     /* We do custom processing for the "as" parameter. */
 265     gmx::SelectionParserParameterList::iterator asparam = ++params->begin();
 266     if (asparam != params->end() && asparam->name() == sm_same.param[1].name)
 267     {
 268         const gmx::SelectionParserValueList& asvalues = asparam->values();
 269         if (asvalues.size() != 1 || !asvalues.front().hasExpressionValue())
 270         {
 271             // TODO: Think about providing more informative context.
 272             GMX_THROW(gmx::InvalidInputError(
 273                     "'same ... as' should be followed by a single expression"));
 274         }
 275         const gmx::SelectionTreeElementPointer& child = asvalues.front().expr;
 276         /* Create a second keyword evaluation element for the keyword given as
 277          * the first parameter, evaluating the keyword in the group given by the
 278          * second parameter. */
 279         gmx::SelectionTreeElementPointer kwelem =
 280                 _gmx_sel_init_keyword_evaluator(kwmethod, child, scanner);
 281         /* Replace the second parameter with one with a value from \p kwelem. */
 282         std::string pname = asparam->name();
 283         *asparam          = gmx::SelectionParserParameter::createFromExpression(pname, kwelem);
 284     }
 285 }
 286
 287 static void init_same(const gmx_mtop_t* /* top */, int /* npar */, gmx_ana_selparam_t* param, void* data)
 288 {
 289     t_methoddata_same* d = static_cast<t_methoddata_same*>(data);
 290
 291     d->val.ptr = param[0].val.u.ptr;
 292     d->as.ptr  = param[1].val.u.ptr;
 293     if (param[1].val.type == STR_VALUE)
 294     {
 295         snew(d->as_s_sorted, d->nas);
 296     }
 297     if (!(param[0].flags & SPAR_ATOMVAL))
 298     {
 299         GMX_THROW(
 300                 gmx::InvalidInputError("The 'same' selection keyword combined with a "
 301                                        "non-keyword does not make sense"));
 302     }
 303 }
 304
 305 /*!
 306  * \param data Data to free (should point to a \ref t_methoddata_same).
 307  */
 308 static void free_data_same(void* data)
 309 {
 310     t_methoddata_same* d = static_cast<t_methoddata_same*>(data);
 311
 312     sfree(d->as_s_sorted);
 313     sfree(d);
 314 }
 315
 316 /*! \brief
 317  * Helper function for comparison of two integers.
 318  */
 319 static int cmp_int(const void* a, const void* b)
 320 {
 321     if (*reinterpret_cast<const int*>(a) < *reinterpret_cast<const int*>(b))
 322     {
 323         return -1;
 324     }
 325     if (*reinterpret_cast<const int*>(a) > *reinterpret_cast<const int*>(b))
 326     {
 327         return 1;
 328     }
 329     return 0;
 330 }
 331
 332 static void init_frame_same_int(const gmx::SelMethodEvalContext& /*context*/, void* data)
 333 {
 334     t_methoddata_same* d = static_cast<t_methoddata_same*>(data);
 335     int                i, j;
 336
 337     /* Collapse adjacent values, and check whether the array is sorted. */
 338     d->bSorted = true;
 339     if (d->nas == 0)
 340     {
 341         return;
 342     }
 343     for (i = 1, j = 0; i < d->nas; ++i)
 344     {
 345         if (d->as.i[i] != d->as.i[j])
 346         {
 347             if (d->as.i[i] < d->as.i[j])
 348             {
 349                 d->bSorted = false;
 350             }
 351             ++j;
 352             d->as.i[j] = d->as.i[i];
 353         }
 354     }
 355     d->nas = j + 1;
 356
 357     if (!d->bSorted)
 358     {
 359         qsort(d->as.i, d->nas, sizeof(d->as.i[0]), &cmp_int);
 360         /* More identical values may become adjacent after sorting. */
 361         for (i = 1, j = 0; i < d->nas; ++i)
 362         {
 363             if (d->as.i[i] != d->as.i[j])
 364             {
 365                 ++j;
 366                 d->as.i[j] = d->as.i[i];
 367             }
 368         }
 369         d->nas = j + 1;
 370     }
 371 }
 372
 373 /*!
 374  * See sel_updatefunc() for description of the parameters.
 375  * \p data should point to a \c t_methoddata_same.
 376  *
 377  * Calculates which values in \c data->val.i can be found in \c data->as.i
 378  * (assumed sorted), and writes the corresponding atoms to output.
 379  * If \c data->val is sorted, uses a linear scan of both arrays, otherwise a
 380  * binary search of \c data->as is performed for each block of values in
 381  * \c data->val.
 382  */
 383 static void evaluate_same_int(const gmx::SelMethodEvalContext& /*context*/,
 384                               gmx_ana_index_t*    g,
 385                               gmx_ana_selvalue_t* out,
 386                               void*               data)
 387 {
 388     t_methoddata_same* d = static_cast<t_methoddata_same*>(data);
 389     int                i, j;
 390
 391     out->u.g->isize = 0;
 392     i = j = 0;
 393     while (j < g->isize)
 394     {
 395         if (d->bSorted)
 396         {
 397             /* If we are sorted, we can do a simple linear scan. */
 398             while (i < d->nas && d->as.i[i] < d->val.i[j])
 399             {
 400                 ++i;
 401             }
 402         }
 403         else
 404         {
 405             /* If not, we must do a binary search of all the values. */
 406             int i1, i2;
 407
 408             i1 = 0;
 409             i2 = d->nas;
 410             while (i2 - i1 > 1)
 411             {
 412                 int itry = (i1 + i2) / 2;
 413                 if (d->as.i[itry] <= d->val.i[j])
 414                 {
 415                     i1 = itry;
 416                 }
 417                 else
 418                 {
 419                     i2 = itry;
 420                 }
 421             }
 422             i = (d->as.i[i1] == d->val.i[j] ? i1 : d->nas);
 423         }
 424         /* Check whether the value was found in the as list. */
 425         if (i == d->nas || d->as.i[i] != d->val.i[j])
 426         {
 427             /* If not, skip all atoms with the same value. */
 428             int tmpval = d->val.i[j];
 429             ++j;
 430             while (j < g->isize && d->val.i[j] == tmpval)
 431             {
 432                 ++j;
 433             }
 434         }
 435         else
 436         {
 437             /* Copy all the atoms with this value to the output. */
 438             while (j < g->isize && d->val.i[j] == d->as.i[i])
 439             {
 440                 out->u.g->index[out->u.g->isize++] = g->index[j];
 441                 ++j;
 442             }
 443         }
 444         if (j < g->isize && d->val.i[j] < d->val.i[j - 1])
 445         {
 446             d->bSorted = false;
 447         }
 448     }
 449 }
 450
 451 /*! \brief
 452  * Helper function for comparison of two strings.
 453  */
 454 static int cmp_str(const void* a, const void* b)
 455 {
 456     return strcmp(*static_cast<char* const*>(a), *static_cast<char* const*>(b));
 457 }
 458
 459 static void init_frame_same_str(const gmx::SelMethodEvalContext& /*context*/, void* data)
 460 {
 461     t_methoddata_same* d = static_cast<t_methoddata_same*>(data);
 462     int                i, j;
 463
 464     /* Collapse adjacent values.
 465      * For strings, it's unlikely that the values would be sorted originally,
 466      * so set bSorted always to false. */
 467     d->bSorted = false;
 468     if (d->nas == 0)
 469     {
 470         return;
 471     }
 472     d->as_s_sorted[0] = d->as.s[0];
 473     for (i = 1, j = 0; i < d->nas; ++i)
 474     {
 475         if (strcmp(d->as.s[i], d->as_s_sorted[j]) != 0)
 476         {
 477             ++j;
 478             d->as_s_sorted[j] = d->as.s[i];
 479         }
 480     }
 481     d->nas = j + 1;
 482
 483     qsort(d->as_s_sorted, d->nas, sizeof(d->as_s_sorted[0]), &cmp_str);
 484     /* More identical values may become adjacent after sorting. */
 485     for (i = 1, j = 0; i < d->nas; ++i)
 486     {
 487         if (strcmp(d->as_s_sorted[i], d->as_s_sorted[j]) != 0)
 488         {
 489             ++j;
 490             d->as_s_sorted[j] = d->as_s_sorted[i];
 491         }
 492     }
 493     d->nas = j + 1;
 494 }
 495
 496 /*!
 497  * See sel_updatefunc() for description of the parameters.
 498  * \p data should point to a \c t_methoddata_same.
 499  *
 500  * Calculates which strings in \c data->val.s can be found in \c data->as.s
 501  * (assumed sorted), and writes the corresponding atoms to output.
 502  * A binary search of \c data->as is performed for each block of values in
 503  * \c data->val.
 504  */
 505 static void evaluate_same_str(const gmx::SelMethodEvalContext& /*context*/,
 506                               gmx_ana_index_t*    g,
 507                               gmx_ana_selvalue_t* out,
 508                               void*               data)
 509 {
 510     t_methoddata_same* d = static_cast<t_methoddata_same*>(data);
 511     int                j;
 512
 513     out->u.g->isize = 0;
 514     j               = 0;
 515     while (j < g->isize)
 516     {
 517         /* Do a binary search of the strings. */
 518         void* ptr = nullptr;
 519         if (d->nas > 0)
 520         {
 521             ptr = bsearch(&d->val.s[j], d->as_s_sorted, d->nas, sizeof(d->as_s_sorted[0]), &cmp_str);
 522         }
 523         /* Check whether the value was found in the as list. */
 524         if (ptr == nullptr)
 525         {
 526             /* If not, skip all atoms with the same value. */
 527             const char* tmpval = d->val.s[j];
 528             ++j;
 529             while (j < g->isize && strcmp(d->val.s[j], tmpval) == 0)
 530             {
 531                 ++j;
 532             }
 533         }
 534         else
 535         {
 536             const char* tmpval = d->val.s[j];
 537             /* Copy all the atoms with this value to the output. */
 538             while (j < g->isize && strcmp(d->val.s[j], tmpval) == 0)
 539             {
 540                 out->u.g->index[out->u.g->isize++] = g->index[j];
 541                 ++j;
 542             }
 543         }
 544     }
 545 }