src/gromacs/selection/sm_same.cpp

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2009, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \internal \file
  32  * \brief
  33  * Implements the \p same selection method.
  34  *
  35  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  36  * \ingroup module_selection
  37  */
  38 #include <stdlib.h>
  39
  40 #include "gromacs/legacyheaders/macros.h"
  41 #include "gromacs/legacyheaders/smalloc.h"
  42 #include "gromacs/legacyheaders/string2.h"
  43
  44 #include "gromacs/selection/selmethod.h"
  45 #include "gromacs/utility/exceptions.h"
  46
  47 #include "keywords.h"
  48 #include "parsetree.h"
  49 #include "selelem.h"
  50
  51 /*! \internal \brief
  52  * Data structure for the \p same selection method.
  53  *
  54  * To avoid duplicate initialization code, the same data structure is used
  55  * for matching both integer and string keywords; hence the unions.
  56  */
  57 typedef struct
  58 {
  59     /** Value for each atom to match. */
  60     union
  61     {
  62         int                 *i;
  63         char               **s;
  64         void                *ptr;
  65     }                        val;
  66     /*! \brief
  67      * Number of values in the \p as array.
  68      *
  69      * For string values, this is actually the number of values in the
  70      * \p as_s_sorted array.
  71      */
  72     int                      nas;
  73     /** Values to match against. */
  74     union
  75     {
  76         int                 *i;
  77         char               **s;
  78         void                *ptr;
  79     }                        as;
  80     /*! \brief
  81      * Separate array for sorted \p as.s array.
  82      *
  83      * The array of strings returned as the output value of a parameter should
  84      * not be messed with to avoid memory corruption (the pointers in the array
  85      * may be reused for several evaluations), so we keep our own copy for
  86      * modifications.
  87      */
  88     char                   **as_s_sorted;
  89     /** Whether simple matching can be used. */
  90     bool                     bSorted;
  91 } t_methoddata_same;
  92
  93 /** Allocates data for the \p same selection method. */
  94 static void *
  95 init_data_same(int npar, gmx_ana_selparam_t *param);
  96 /** Initializes the \p same selection method. */
  97 static void
  98 init_same(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
  99 /** Frees the data allocated for the \p same selection method. */
 100 static void
 101 free_data_same(void *data);
 102 /** Initializes the evaluation of the \p same selection method for a frame. */
 103 static void
 104 init_frame_same_int(t_topology *top, t_trxframe *fr, t_pbc *pbc, void *data);
 105 /** Evaluates the \p same selection method. */
 106 static void
 107 evaluate_same_int(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 108                   gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
 109 /** Initializes the evaluation of the \p same selection method for a frame. */
 110 static void
 111 init_frame_same_str(t_topology *top, t_trxframe *fr, t_pbc *pbc, void *data);
 112 /** Evaluates the \p same selection method. */
 113 static void
 114 evaluate_same_str(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 115                  gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data);
 116
 117 /** Parameters for the \p same selection method. */
 118 static gmx_ana_selparam_t smparams_same_int[] = {
 119     {NULL, {INT_VALUE, -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_ATOMVAL},
 120     {"as", {INT_VALUE, -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_VARNUM},
 121 };
 122
 123 /** Parameters for the \p same selection method. */
 124 static gmx_ana_selparam_t smparams_same_str[] = {
 125     {NULL, {STR_VALUE, -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_ATOMVAL},
 126     {"as", {STR_VALUE, -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_VARNUM},
 127 };
 128
 129 /** Help text for the \p same selection method. */
 130 static const char *help_same[] = {
 131     "EXTENDING SELECTIONS[PAR]",
 132
 133     "[TT]same KEYWORD as ATOM_EXPR[tt][PAR]",
 134
 135     "The keyword [TT]same[tt] can be used to select all atoms for which",
 136     "the given [TT]KEYWORD[tt] matches any of the atoms in [TT]ATOM_EXPR[tt].",
 137     "Keywords that evaluate to integer or string values are supported.",
 138 };
 139
 140 /*! \internal \brief Selection method data for the \p same method. */
 141 gmx_ana_selmethod_t sm_same = {
 142     "same", GROUP_VALUE, 0,
 143     asize(smparams_same_int), smparams_same_int,
 144     &init_data_same,
 145     NULL,
 146     &init_same,
 147     NULL,
 148     &free_data_same,
 149     &init_frame_same_int,
 150     &evaluate_same_int,
 151     NULL,
 152     {"same KEYWORD as ATOM_EXPR", asize(help_same), help_same},
 153 };
 154
 155 /*! \brief
 156  * Selection method data for the \p same method.
 157  *
 158  * This selection method is used for matching string keywords. The parser
 159  * never sees this method; _gmx_selelem_custom_init_same() replaces sm_same
 160  * with this method in cases where it is required.
 161  */
 162 static gmx_ana_selmethod_t sm_same_str = {
 163     "same", GROUP_VALUE, SMETH_SINGLEVAL,
 164     asize(smparams_same_str), smparams_same_str,
 165     &init_data_same,
 166     NULL,
 167     &init_same,
 168     NULL,
 169     &free_data_same,
 170     &init_frame_same_str,
 171     &evaluate_same_str,
 172     NULL,
 173     {"same KEYWORD as ATOM_EXPR", asize(help_same), help_same},
 174 };
 175
 176 /*!
 177  * \param[in]     npar  Not used (should be 2).
 178  * \param[in,out] param Method parameters (should point to a copy of
 179  *      ::smparams_same_int or ::smparams_same_str).
 180  * \returns Pointer to the allocated data (\ref t_methoddata_same).
 181  */
 182 static void *
 183 init_data_same(int npar, gmx_ana_selparam_t *param)
 184 {
 185     t_methoddata_same *data;
 186
 187     snew(data, 1);
 188     data->as_s_sorted = NULL;
 189     param[1].nvalptr = &data->nas;
 190     return data;
 191 }
 192
 193 /*!
 194  * \param[in,out] method  The method to initialize.
 195  * \param[in,out] params  Pointer to the first parameter.
 196  * \param[in]     scanner Scanner data structure.
 197  * \returns       0 on success, a non-zero error code on error.
 198  *
 199  * If \p *method is not a \c same method, this function returns zero
 200  * immediately.
 201  */
 202 int
 203 _gmx_selelem_custom_init_same(gmx_ana_selmethod_t **method,
 204                               const gmx::SelectionParserParameterListPointer &params,
 205                               void *scanner)
 206 {
 207
 208     /* Do nothing if this is not a same method. */
 209     if (!*method || (*method)->name != sm_same.name || params->empty())
 210     {
 211         return 0;
 212     }
 213
 214     const gmx::SelectionParserValueList &kwvalues = params->front().values();
 215     if (kwvalues.size() != 1 || !kwvalues.front().hasExpressionValue()
 216         || kwvalues.front().expr->type != SEL_EXPRESSION)
 217     {
 218         _gmx_selparser_error(scanner, "'same' should be followed by a single keyword");
 219         return -1;
 220     }
 221     gmx_ana_selmethod_t *kwmethod = kwvalues.front().expr->u.expr.method;
 222     if (kwmethod->type == STR_VALUE)
 223     {
 224         *method = &sm_same_str;
 225     }
 226
 227     /* We do custom processing for the "as" parameter. */
 228     gmx::SelectionParserParameterList::iterator asparam = ++params->begin();
 229     if (asparam != params->end() && asparam->name() == sm_same.param[1].name)
 230     {
 231         gmx::SelectionParserParameterList kwparams;
 232         gmx::SelectionParserValueListPointer values(
 233                 new gmx::SelectionParserValueList(asparam->values()));
 234         kwparams.push_back(
 235                 gmx::SelectionParserParameter::create(NULL, move(values)));
 236
 237         /* Create a second keyword evaluation element for the keyword given as
 238          * the first parameter, evaluating the keyword in the group given by the
 239          * second parameter. */
 240         gmx::SelectionTreeElementPointer kwelem
 241             = _gmx_sel_init_keyword_evaluator(kwmethod, kwparams, scanner);
 242         // FIXME: Use exceptions.
 243         if (!kwelem)
 244         {
 245             return -1;
 246         }
 247         /* Replace the second parameter with one with a value from \p kwelem. */
 248         std::string pname = asparam->name();
 249         *asparam = gmx::SelectionParserParameter::createFromExpression(pname, kwelem);
 250     }
 251     return 0;
 252 }
 253
 254 /*!
 255  * \param   top   Not used.
 256  * \param   npar  Not used (should be 2).
 257  * \param   param Initialized method parameters (should point to a copy of
 258  *      ::smparams_same_int or ::smparams_same_str).
 259  * \param   data  Pointer to \ref t_methoddata_same to initialize.
 260  * \returns 0 on success, -1 on failure.
 261  */
 262 static void
 263 init_same(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
 264 {
 265     t_methoddata_same *d = (t_methoddata_same *)data;
 266
 267     d->val.ptr = param[0].val.u.ptr;
 268     d->as.ptr  = param[1].val.u.ptr;
 269     if (param[1].val.type == STR_VALUE)
 270     {
 271         snew(d->as_s_sorted, d->nas);
 272     }
 273     if (!(param[0].flags & SPAR_ATOMVAL))
 274     {
 275         GMX_THROW(gmx::InvalidInputError(
 276                     "The 'same' selection keyword combined with a "
 277                     "non-keyword does not make sense"));
 278     }
 279 }
 280
 281 /*!
 282  * \param data Data to free (should point to a \ref t_methoddata_same).
 283  */
 284 static void
 285 free_data_same(void *data)
 286 {
 287     t_methoddata_same *d = (t_methoddata_same *)data;
 288
 289     sfree(d->as_s_sorted);
 290 }
 291
 292 /*! \brief
 293  * Helper function for comparison of two integers.
 294  */
 295 static int
 296 cmp_int(const void *a, const void *b)
 297 {
 298     if (*(int *)a < *(int *)b)
 299     {
 300         return -1;
 301     }
 302     if (*(int *)a > *(int *)b)
 303     {
 304         return 1;
 305     }
 306     return 0;
 307 }
 308
 309 /*!
 310  * \param[in]  top  Not used.
 311  * \param[in]  fr   Current frame.
 312  * \param[in]  pbc  PBC structure.
 313  * \param      data Should point to a \ref t_methoddata_same.
 314  *
 315  * Sorts the \c data->as.i array and removes identical values for faster and
 316  * simpler lookup.
 317  */
 318 static void
 319 init_frame_same_int(t_topology *top, t_trxframe *fr, t_pbc *pbc, void *data)
 320 {
 321     t_methoddata_same *d = (t_methoddata_same *)data;
 322     int                i, j;
 323
 324     /* Collapse adjacent values, and check whether the array is sorted. */
 325     d->bSorted = true;
 326     for (i = 1, j = 0; i < d->nas; ++i)
 327     {
 328         if (d->as.i[i] != d->as.i[j])
 329         {
 330             if (d->as.i[i] < d->as.i[j])
 331             {
 332                 d->bSorted = false;
 333             }
 334             ++j;
 335             d->as.i[j] = d->as.i[i];
 336         }
 337     }
 338     d->nas = j + 1;
 339
 340     if (!d->bSorted)
 341     {
 342         qsort(d->as.i, d->nas, sizeof(d->as.i[0]), &cmp_int);
 343         /* More identical values may become adjacent after sorting. */
 344         for (i = 1, j = 0; i < d->nas; ++i)
 345         {
 346             if (d->as.i[i] != d->as.i[j])
 347             {
 348                 ++j;
 349                 d->as.i[j] = d->as.i[i];
 350             }
 351         }
 352         d->nas = j + 1;
 353     }
 354 }
 355
 356 /*!
 357  * See sel_updatefunc() for description of the parameters.
 358  * \p data should point to a \c t_methoddata_same.
 359  *
 360  * Calculates which values in \c data->val.i can be found in \c data->as.i
 361  * (assumed sorted), and writes the corresponding atoms to output.
 362  * If \c data->val is sorted, uses a linear scan of both arrays, otherwise a
 363  * binary search of \c data->as is performed for each block of values in
 364  * \c data->val.
 365  */
 366 static void
 367 evaluate_same_int(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 368                   gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
 369 {
 370     t_methoddata_same *d = (t_methoddata_same *)data;
 371     int                    i, j;
 372
 373     out->u.g->isize = 0;
 374     i = j = 0;
 375     while (j < g->isize)
 376     {
 377         if (d->bSorted)
 378         {
 379             /* If we are sorted, we can do a simple linear scan. */
 380             while (i < d->nas && d->as.i[i] < d->val.i[j]) ++i;
 381         }
 382         else
 383         {
 384             /* If not, we must do a binary search of all the values. */
 385             int i1, i2;
 386
 387             i1 = 0;
 388             i2 = d->nas;
 389             while (i2 - i1 > 1)
 390             {
 391                 int itry = (i1 + i2) / 2;
 392                 if (d->as.i[itry] <= d->val.i[j])
 393                 {
 394                     i1 = itry;
 395                 }
 396                 else
 397                 {
 398                     i2 = itry;
 399                 }
 400             }
 401             i = (d->as.i[i1] == d->val.i[j] ? i1 : d->nas);
 402         }
 403         /* Check whether the value was found in the as list. */
 404         if (i == d->nas || d->as.i[i] != d->val.i[j])
 405         {
 406             /* If not, skip all atoms with the same value. */
 407             int tmpval = d->val.i[j];
 408             ++j;
 409             while (j < g->isize && d->val.i[j] == tmpval) ++j;
 410         }
 411         else
 412         {
 413             /* Copy all the atoms with this value to the output. */
 414             while (j < g->isize && d->val.i[j] == d->as.i[i])
 415             {
 416                 out->u.g->index[out->u.g->isize++] = g->index[j];
 417                 ++j;
 418             }
 419         }
 420         if (j < g->isize && d->val.i[j] < d->val.i[j - 1])
 421         {
 422             d->bSorted = false;
 423         }
 424     }
 425 }
 426
 427 /*! \brief
 428  * Helper function for comparison of two strings.
 429  */
 430 static int
 431 cmp_str(const void *a, const void *b)
 432 {
 433     return strcmp(*(char **)a, *(char **)b);
 434 }
 435
 436 /*!
 437  * \param[in]  top  Not used.
 438  * \param[in]  fr   Current frame.
 439  * \param[in]  pbc  PBC structure.
 440  * \param      data Should point to a \ref t_methoddata_same.
 441  *
 442  * Sorts the \c data->as.s array and removes identical values for faster and
 443  * simpler lookup.
 444  */
 445 static void
 446 init_frame_same_str(t_topology *top, t_trxframe *fr, t_pbc *pbc, void *data)
 447 {
 448     t_methoddata_same *d = (t_methoddata_same *)data;
 449     int                i, j;
 450
 451     /* Collapse adjacent values.
 452      * For strings, it's unlikely that the values would be sorted originally,
 453      * so set bSorted always to false. */
 454     d->bSorted = false;
 455     d->as_s_sorted[0] = d->as.s[0];
 456     for (i = 1, j = 0; i < d->nas; ++i)
 457     {
 458         if (strcmp(d->as.s[i], d->as_s_sorted[j]) != 0)
 459         {
 460             ++j;
 461             d->as_s_sorted[j] = d->as.s[i];
 462         }
 463     }
 464     d->nas = j + 1;
 465
 466     qsort(d->as_s_sorted, d->nas, sizeof(d->as_s_sorted[0]), &cmp_str);
 467     /* More identical values may become adjacent after sorting. */
 468     for (i = 1, j = 0; i < d->nas; ++i)
 469     {
 470         if (strcmp(d->as_s_sorted[i], d->as_s_sorted[j]) != 0)
 471         {
 472             ++j;
 473             d->as_s_sorted[j] = d->as_s_sorted[i];
 474         }
 475     }
 476     d->nas = j + 1;
 477 }
 478
 479 /*!
 480  * See sel_updatefunc() for description of the parameters.
 481  * \p data should point to a \c t_methoddata_same.
 482  *
 483  * Calculates which strings in \c data->val.s can be found in \c data->as.s
 484  * (assumed sorted), and writes the corresponding atoms to output.
 485  * A binary search of \c data->as is performed for each block of values in
 486  * \c data->val.
 487  */
 488 static void
 489 evaluate_same_str(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 490                   gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
 491 {
 492     t_methoddata_same *d = (t_methoddata_same *)data;
 493     int                    j;
 494
 495     out->u.g->isize = 0;
 496     j = 0;
 497     while (j < g->isize)
 498     {
 499         /* Do a binary search of the strings. */
 500         void *ptr;
 501         ptr = bsearch(&d->val.s[j], d->as_s_sorted, d->nas,
 502                       sizeof(d->as_s_sorted[0]), &cmp_str);
 503         /* Check whether the value was found in the as list. */
 504         if (ptr == NULL)
 505         {
 506             /* If not, skip all atoms with the same value. */
 507             const char *tmpval = d->val.s[j];
 508             ++j;
 509             while (j < g->isize && strcmp(d->val.s[j], tmpval) == 0) ++j;
 510         }
 511         else
 512         {
 513             const char *tmpval = d->val.s[j];
 514             /* Copy all the atoms with this value to the output. */
 515             while (j < g->isize && strcmp(d->val.s[j], tmpval) == 0)
 516             {
 517                 out->u.g->index[out->u.g->isize++] = g->index[j];
 518                 ++j;
 519             }
 520         }
 521     }
 522 }