src/gmxlib/selection/sm_same.c

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