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