src/gmxlib/selection/sm_same.c

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