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     sfree(d);
 291 }
 292
 293 /*! \brief
 294  * Helper function for comparison of two integers.
 295  */
 296 static int
 297 cmp_int(const void *a, const void *b)
 298 {
 299     if (*(int *)a < *(int *)b)
 300     {
 301         return -1;
 302     }
 303     if (*(int *)a > *(int *)b)
 304     {
 305         return 1;
 306     }
 307     return 0;
 308 }
 309
 310 /*!
 311  * \param[in]  top  Not used.
 312  * \param[in]  fr   Current frame.
 313  * \param[in]  pbc  PBC structure.
 314  * \param      data Should point to a \ref t_methoddata_same.
 315  *
 316  * Sorts the \c data->as.i array and removes identical values for faster and
 317  * simpler lookup.
 318  */
 319 static void
 320 init_frame_same_int(t_topology *top, t_trxframe *fr, t_pbc *pbc, void *data)
 321 {
 322     t_methoddata_same *d = (t_methoddata_same *)data;
 323     int                i, j;
 324
 325     /* Collapse adjacent values, and check whether the array is sorted. */
 326     d->bSorted = true;
 327     for (i = 1, j = 0; i < d->nas; ++i)
 328     {
 329         if (d->as.i[i] != d->as.i[j])
 330         {
 331             if (d->as.i[i] < d->as.i[j])
 332             {
 333                 d->bSorted = false;
 334             }
 335             ++j;
 336             d->as.i[j] = d->as.i[i];
 337         }
 338     }
 339     d->nas = j + 1;
 340
 341     if (!d->bSorted)
 342     {
 343         qsort(d->as.i, d->nas, sizeof(d->as.i[0]), &cmp_int);
 344         /* More identical values may become adjacent after sorting. */
 345         for (i = 1, j = 0; i < d->nas; ++i)
 346         {
 347             if (d->as.i[i] != d->as.i[j])
 348             {
 349                 ++j;
 350                 d->as.i[j] = d->as.i[i];
 351             }
 352         }
 353         d->nas = j + 1;
 354     }
 355 }
 356
 357 /*!
 358  * See sel_updatefunc() for description of the parameters.
 359  * \p data should point to a \c t_methoddata_same.
 360  *
 361  * Calculates which values in \c data->val.i can be found in \c data->as.i
 362  * (assumed sorted), and writes the corresponding atoms to output.
 363  * If \c data->val is sorted, uses a linear scan of both arrays, otherwise a
 364  * binary search of \c data->as is performed for each block of values in
 365  * \c data->val.
 366  */
 367 static void
 368 evaluate_same_int(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 369                   gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
 370 {
 371     t_methoddata_same *d = (t_methoddata_same *)data;
 372     int                    i, j;
 373
 374     out->u.g->isize = 0;
 375     i = j = 0;
 376     while (j < g->isize)
 377     {
 378         if (d->bSorted)
 379         {
 380             /* If we are sorted, we can do a simple linear scan. */
 381             while (i < d->nas && d->as.i[i] < d->val.i[j]) ++i;
 382         }
 383         else
 384         {
 385             /* If not, we must do a binary search of all the values. */
 386             int i1, i2;
 387
 388             i1 = 0;
 389             i2 = d->nas;
 390             while (i2 - i1 > 1)
 391             {
 392                 int itry = (i1 + i2) / 2;
 393                 if (d->as.i[itry] <= d->val.i[j])
 394                 {
 395                     i1 = itry;
 396                 }
 397                 else
 398                 {
 399                     i2 = itry;
 400                 }
 401             }
 402             i = (d->as.i[i1] == d->val.i[j] ? i1 : d->nas);
 403         }
 404         /* Check whether the value was found in the as list. */
 405         if (i == d->nas || d->as.i[i] != d->val.i[j])
 406         {
 407             /* If not, skip all atoms with the same value. */
 408             int tmpval = d->val.i[j];
 409             ++j;
 410             while (j < g->isize && d->val.i[j] == tmpval) ++j;
 411         }
 412         else
 413         {
 414             /* Copy all the atoms with this value to the output. */
 415             while (j < g->isize && d->val.i[j] == d->as.i[i])
 416             {
 417                 out->u.g->index[out->u.g->isize++] = g->index[j];
 418                 ++j;
 419             }
 420         }
 421         if (j < g->isize && d->val.i[j] < d->val.i[j - 1])
 422         {
 423             d->bSorted = false;
 424         }
 425     }
 426 }
 427
 428 /*! \brief
 429  * Helper function for comparison of two strings.
 430  */
 431 static int
 432 cmp_str(const void *a, const void *b)
 433 {
 434     return strcmp(*(char **)a, *(char **)b);
 435 }
 436
 437 /*!
 438  * \param[in]  top  Not used.
 439  * \param[in]  fr   Current frame.
 440  * \param[in]  pbc  PBC structure.
 441  * \param      data Should point to a \ref t_methoddata_same.
 442  *
 443  * Sorts the \c data->as.s array and removes identical values for faster and
 444  * simpler lookup.
 445  */
 446 static void
 447 init_frame_same_str(t_topology *top, t_trxframe *fr, t_pbc *pbc, void *data)
 448 {
 449     t_methoddata_same *d = (t_methoddata_same *)data;
 450     int                i, j;
 451
 452     /* Collapse adjacent values.
 453      * For strings, it's unlikely that the values would be sorted originally,
 454      * so set bSorted always to false. */
 455     d->bSorted = false;
 456     d->as_s_sorted[0] = d->as.s[0];
 457     for (i = 1, j = 0; i < d->nas; ++i)
 458     {
 459         if (strcmp(d->as.s[i], d->as_s_sorted[j]) != 0)
 460         {
 461             ++j;
 462             d->as_s_sorted[j] = d->as.s[i];
 463         }
 464     }
 465     d->nas = j + 1;
 466
 467     qsort(d->as_s_sorted, d->nas, sizeof(d->as_s_sorted[0]), &cmp_str);
 468     /* More identical values may become adjacent after sorting. */
 469     for (i = 1, j = 0; i < d->nas; ++i)
 470     {
 471         if (strcmp(d->as_s_sorted[i], d->as_s_sorted[j]) != 0)
 472         {
 473             ++j;
 474             d->as_s_sorted[j] = d->as_s_sorted[i];
 475         }
 476     }
 477     d->nas = j + 1;
 478 }
 479
 480 /*!
 481  * See sel_updatefunc() for description of the parameters.
 482  * \p data should point to a \c t_methoddata_same.
 483  *
 484  * Calculates which strings in \c data->val.s can be found in \c data->as.s
 485  * (assumed sorted), and writes the corresponding atoms to output.
 486  * A binary search of \c data->as is performed for each block of values in
 487  * \c data->val.
 488  */
 489 static void
 490 evaluate_same_str(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 491                   gmx_ana_index_t *g, gmx_ana_selvalue_t *out, void *data)
 492 {
 493     t_methoddata_same *d = (t_methoddata_same *)data;
 494     int                    j;
 495
 496     out->u.g->isize = 0;
 497     j = 0;
 498     while (j < g->isize)
 499     {
 500         /* Do a binary search of the strings. */
 501         void *ptr;
 502         ptr = bsearch(&d->val.s[j], d->as_s_sorted, d->nas,
 503                       sizeof(d->as_s_sorted[0]), &cmp_str);
 504         /* Check whether the value was found in the as list. */
 505         if (ptr == NULL)
 506         {
 507             /* If not, skip all atoms with the same value. */
 508             const char *tmpval = d->val.s[j];
 509             ++j;
 510             while (j < g->isize && strcmp(d->val.s[j], tmpval) == 0) ++j;
 511         }
 512         else
 513         {
 514             const char *tmpval = d->val.s[j];
 515             /* Copy all the atoms with this value to the output. */
 516             while (j < g->isize && strcmp(d->val.s[j], tmpval) == 0)
 517             {
 518                 out->u.g->index[out->u.g->isize++] = g->index[j];
 519                 ++j;
 520             }
 521         }
 522     }
 523 }