src/gmxlib/selection/sm_merge.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 merging selection modifier.
  33  */
  34 #ifdef HAVE_CONFIG_H
  35 #include <config.h>
  36 #endif
  37
  38 #include <macros.h>
  39 #include <smalloc.h>
  40 #include <vec.h>
  41
  42 #include <position.h>
  43 #include <selmethod.h>
  44
  45 /*! \internal \brief
  46  * Data structure for the merging selection modifiers.
  47  */
  48 typedef struct
  49 {
  50     /** Input positions. */
  51     gmx_ana_pos_t    p1;
  52     /** Other input positions. */
  53     gmx_ana_pos_t    p2;
  54     /** Group to store the output atom indices. */
  55     gmx_ana_index_t  g;
  56     /** Stride for merging (\c stride values from \c p1 for each in \c p2). */
  57     int              stride;
  58 } t_methoddata_merge;
  59
  60 /** Allocates data for the merging selection modifiers. */
  61 static void *
  62 init_data_merge(int npar, gmx_ana_selparam_t *param);
  63 /** Initializes data for the merging selection modifiers. */
  64 static int
  65 init_merge(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
  66 /** Initializes output for the \p merge selection modifier. */
  67 static int
  68 init_output_merge(t_topology *top, gmx_ana_selvalue_t *out, void *data);
  69 /** Initializes output for the \p plus selection modifier. */
  70 static int
  71 init_output_plus(t_topology *top, gmx_ana_selvalue_t *out, void *data);
  72 /** Frees the memory allocated for the merging selection modifiers. */
  73 static void
  74 free_data_merge(void *data);
  75 /** Evaluates the \p merge selection modifier. */
  76 static int
  77 evaluate_merge(t_topology *top, t_trxframe *fr, t_pbc *pbc,
  78                gmx_ana_pos_t *p, gmx_ana_selvalue_t *out, void *data);
  79 /** Evaluates the \p plus selection modifier. */
  80 static int
  81 evaluate_plus(t_topology *top, t_trxframe *fr, t_pbc *pbc,
  82               gmx_ana_pos_t *p, gmx_ana_selvalue_t *out, void *data);
  83
  84 /** Parameters for the merging selection modifiers. */
  85 static gmx_ana_selparam_t smparams_merge[] = {
  86     {NULL,       {POS_VALUE, -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_VARNUM},
  87     {NULL,       {POS_VALUE, -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_VARNUM},
  88     {"stride",   {INT_VALUE,  1, {NULL}}, NULL, SPAR_OPTIONAL},
  89 };
  90
  91 /** Help text for the merging selection modifiers. */
  92 static const char *help_merge[] = {
  93     "MERGING SELECTIONS[PAR]",
  94
  95     "[TT]POSEXPR merge POSEXPR [stride INT][tt][BR]",
  96     "[TT]POSEXPR merge POSEXPR [merge POSEXPR ...][tt][BR]",
  97     "[TT]POSEXPR plus POSEXPR [plus POSEXPR ...][tt][PAR]",
  98
  99     "Basic selection keywords can only create selections where each atom",
 100     "occurs at most once. The [TT]merge[tt] and [TT]plus[tt] selection",
 101     "keywords can be used to work around this limitation. Both create",
 102     "a selection that contains the positions from all the given position",
 103     "expressions, even if they contain duplicates.",
 104     "The difference between the two is that [TT]merge[tt] expects two or more",
 105     "selections with the same number of positions, and the output contains",
 106     "the input positions selected from each expression in turn, i.e.,",
 107     "the output is like A1 B1 A2 B2 and so on. It is also possible to merge",
 108     "selections of unequal size as long as the size of the first is a",
 109     "multiple of the second one. The [TT]stride[tt] parameter can be used",
 110     "to explicitly provide this multiplicity.",
 111     "[TT]plus[tt] simply concatenates the positions after each other, and",
 112     "can work also with selections of different sizes.",
 113     "These keywords are valid only at the selection level, not in any",
 114     "subexpressions.[PAR]",
 115 };
 116
 117 /** \internal Selection method data for the \p plus modifier. */
 118 gmx_ana_selmethod_t sm_merge = {
 119     "merge", POS_VALUE, SMETH_MODIFIER,
 120     asize(smparams_merge), smparams_merge,
 121     &init_data_merge,
 122     NULL,
 123     &init_merge,
 124     &init_output_merge,
 125     &free_data_merge,
 126     NULL,
 127     NULL,
 128     &evaluate_merge,
 129     {"merge POSEXPR", asize(help_merge), help_merge},
 130 };
 131
 132 /** \internal Selection method data for the \p plus modifier. */
 133 gmx_ana_selmethod_t sm_plus = {
 134     "plus", POS_VALUE, SMETH_MODIFIER,
 135     asize(smparams_merge)-1, smparams_merge,
 136     &init_data_merge,
 137     NULL,
 138     &init_merge,
 139     &init_output_plus,
 140     &free_data_merge,
 141     NULL,
 142     NULL,
 143     &evaluate_plus,
 144     {"plus POSEXPR", asize(help_merge), help_merge},
 145 };
 146
 147 /*!
 148  * \param[in]     npar  Should be 2 for \c plus and 3 for \c merge.
 149  * \param[in,out] param Method parameters (should point to a copy of
 150  *   \ref smparams_merge).
 151  * \returns Pointer to the allocated data (\p t_methoddata_merge).
 152  *
 153  * Allocates memory for a \p t_methoddata_merge structure.
 154  */
 155 static void *
 156 init_data_merge(int npar, gmx_ana_selparam_t *param)
 157 {
 158     t_methoddata_merge *data;
 159
 160     snew(data, 1);
 161     data->stride = 0;
 162     param[0].val.u.p = &data->p1;
 163     param[1].val.u.p = &data->p2;
 164     if (npar > 2)
 165     {
 166         param[2].val.u.i = &data->stride;
 167     }
 168     return data;
 169 }
 170
 171 /*!
 172  * \param[in] top   Not used.
 173  * \param[in] npar  Not used (should be 2 or 3).
 174  * \param[in] param Method parameters (should point to \ref smparams_merge).
 175  * \param[in] data  Should point to a \p t_methoddata_merge.
 176  * \returns   0 if everything is successful, -1 on error.
 177  */
 178 static int
 179 init_merge(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
 180 {
 181     t_methoddata_merge *d = (t_methoddata_merge *)data;
 182     int                 i;
 183
 184     if (d->stride < 0)
 185     {
 186         fprintf(stderr, "error: stride for merging should be positive\n");
 187         return -1;
 188     }
 189     /* If no stride given, deduce it from the input sizes */
 190     if (d->stride == 0)
 191     {
 192         d->stride = d->p1.nr / d->p2.nr;
 193     }
 194     if (d->p1.nr != d->stride*d->p2.nr)
 195     {
 196         fprintf(stderr, "error: the number of positions to be merged are not compatible\n");
 197         return -1;
 198     }
 199     gmx_ana_index_reserve(&d->g, d->p1.g->isize + d->p2.g->isize);
 200     d->g.isize = d->p1.g->isize + d->p2.g->isize;
 201     return 0;
 202 }
 203
 204 /*! \brief
 205  * Does common initialization to all merging modifiers.
 206  *
 207  * \param[in]     top   Topology data structure.
 208  * \param[in,out] out   Pointer to output data structure.
 209  * \param[in,out] data  Should point to \c t_methoddata_merge.
 210  * \returns       0 for success.
 211  */
 212 static int
 213 init_output_common(t_topology *top, gmx_ana_selvalue_t *out, void *data)
 214 {
 215     t_methoddata_merge *d = (t_methoddata_merge *)data;
 216
 217     if (d->p1.m.type != d->p2.m.type)
 218     {
 219         /* TODO: Maybe we could pick something else here? */
 220         out->u.p->m.type = INDEX_UNKNOWN;
 221     }
 222     else
 223     {
 224         out->u.p->m.type = d->p1.m.type;
 225     }
 226     gmx_ana_pos_reserve(out->u.p, d->p1.nr + d->p2.nr, d->g.isize);
 227     gmx_ana_pos_set_evalgrp(out->u.p, &d->g);
 228     gmx_ana_pos_empty_init(out->u.p);
 229     d->g.isize = 0;
 230     return 0;
 231 }
 232
 233 /*!
 234  * \param[in]     top   Topology data structure.
 235  * \param[in,out] out   Pointer to output data structure.
 236  * \param[in,out] data  Should point to \c t_methoddata_merge.
 237  * \returns       0 for success.
 238  */
 239 static int
 240 init_output_merge(t_topology *top, gmx_ana_selvalue_t *out, void *data)
 241 {
 242     t_methoddata_merge *d = (t_methoddata_merge *)data;
 243     int                 i, j;
 244
 245     init_output_common(top, out, data);
 246     for (i = 0; i < d->p2.nr; ++i)
 247     {
 248         for (j = 0; j < d->stride; ++j)
 249         {
 250             gmx_ana_pos_append_init(out->u.p, &d->g, &d->p1, d->stride*i+j);
 251         }
 252         gmx_ana_pos_append_init(out->u.p, &d->g, &d->p2, i);
 253     }
 254     return 0;
 255 }
 256
 257 /*!
 258  * \param[in]     top   Topology data structure.
 259  * \param[in,out] out   Pointer to output data structure.
 260  * \param[in,out] data  Should point to \c t_methoddata_merge.
 261  * \returns       0 for success.
 262  */
 263 static int
 264 init_output_plus(t_topology *top, gmx_ana_selvalue_t *out, void *data)
 265 {
 266     t_methoddata_merge *d = (t_methoddata_merge *)data;
 267     int                 i;
 268
 269     init_output_common(top, out, data);
 270     for (i = 0; i < d->p1.nr; ++i)
 271     {
 272         gmx_ana_pos_append_init(out->u.p, &d->g, &d->p1, i);
 273     }
 274     for (i = 0; i < d->p2.nr; ++i)
 275     {
 276         gmx_ana_pos_append_init(out->u.p, &d->g, &d->p2, i);
 277     }
 278     return 0;
 279 }
 280
 281 /*!
 282  * \param data Data to free (should point to a \p t_methoddata_merge).
 283  *
 284  * Frees the memory allocated for \c t_methoddata_merge.
 285  */
 286 static void
 287 free_data_merge(void *data)
 288 {
 289     t_methoddata_merge *d = (t_methoddata_merge *)data;
 290
 291     gmx_ana_pos_deinit(&d->p1);
 292     gmx_ana_pos_deinit(&d->p2);
 293     gmx_ana_index_deinit(&d->g);
 294 }
 295
 296 /*!
 297  * \param[in]  top   Not used.
 298  * \param[in]  fr    Not used.
 299  * \param[in]  pbc   Not used.
 300  * \param[in]  p     Positions to merge (should point to \p data->p1).
 301  * \param[out] out   Output data structure (\p out->u.p is used).
 302  * \param[in]  data  Should point to a \p t_methoddata_merge.
 303  * \returns    0 on success.
 304  */
 305 static int
 306 evaluate_merge(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 307                gmx_ana_pos_t *p, gmx_ana_selvalue_t *out, void *data)
 308 {
 309     t_methoddata_merge *d = (t_methoddata_merge *)data;
 310     int                 i, j;
 311     int                 refid;
 312
 313     if (d->p1.nr != d->stride*d->p2.nr)
 314     {
 315         fprintf(stderr, "error: the number of positions to be merged are not compatible\n");
 316         return -1;
 317     }
 318     d->g.isize = 0;
 319     gmx_ana_pos_empty(out->u.p);
 320     for (i = 0; i < d->p2.nr; ++i)
 321     {
 322         for (j = 0; j < d->stride; ++j)
 323         {
 324             refid = d->p1.m.refid[d->stride*i+j];
 325             if (refid != -1)
 326             {
 327                 refid = (d->stride+1) * (refid / d->stride) + (refid % d->stride);
 328             }
 329             gmx_ana_pos_append(out->u.p, &d->g, &d->p1, d->stride*i+j, refid);
 330         }
 331         refid = (d->stride+1)*d->p2.m.refid[i]+d->stride;
 332         gmx_ana_pos_append(out->u.p, &d->g, &d->p2, i, refid);
 333     }
 334     gmx_ana_pos_append_finish(out->u.p);
 335     return 0;
 336 }
 337
 338 /*!
 339  * \param[in]  top   Not used.
 340  * \param[in]  fr    Not used.
 341  * \param[in]  pbc   Not used.
 342  * \param[in]  p     Positions to merge (should point to \p data->p1).
 343  * \param[out] out   Output data structure (\p out->u.p is used).
 344  * \param[in]  data  Should point to a \p t_methoddata_merge.
 345  * \returns    0 on success.
 346  */
 347 static int
 348 evaluate_plus(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 349               gmx_ana_pos_t *p, gmx_ana_selvalue_t *out, void *data)
 350 {
 351     t_methoddata_merge *d = (t_methoddata_merge *)data;
 352     int                 i;
 353     int                 refid;
 354
 355     d->g.isize = 0;
 356     gmx_ana_pos_empty(out->u.p);
 357     for (i = 0; i < d->p1.nr; ++i)
 358     {
 359         refid = d->p1.m.refid[i];
 360         gmx_ana_pos_append(out->u.p, &d->g, &d->p1, i, refid);
 361     }
 362     for (i = 0; i < d->p2.nr; ++i)
 363     {
 364         refid = d->p2.m.refid[i];
 365         if (refid != -1)
 366         {
 367             refid += d->p1.m.b.nr;
 368         }
 369         gmx_ana_pos_append(out->u.p, &d->g, &d->p2, i, refid);
 370     }
 371     gmx_ana_pos_append_finish(out->u.p);
 372     return 0;
 373 }