src/gromacs/selection/sm_merge.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2011,2012,2013,2014, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements the \p merge and \p plus selection modifiers.
  38  *
  39  * \author Teemu Murtola <teemu.murtola@gmail.com>
  40  * \ingroup module_selection
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "gromacs/legacyheaders/macros.h"
  45 #include "gromacs/math/vec.h"
  46 #include "gromacs/selection/position.h"
  47 #include "gromacs/utility/common.h"
  48 #include "gromacs/utility/exceptions.h"
  49 #include "gromacs/utility/smalloc.h"
  50
  51 #include "selmethod.h"
  52
  53 /*! \internal \brief
  54  * Data structure for the merging selection modifiers.
  55  */
  56 typedef struct
  57 {
  58     /** Input positions. */
  59     gmx_ana_pos_t    p1;
  60     /** Other input positions. */
  61     gmx_ana_pos_t    p2;
  62     /** Stride for merging (\c stride values from \c p1 for each in \c p2). */
  63     int              stride;
  64 } t_methoddata_merge;
  65
  66 /** Allocates data for the merging selection modifiers. */
  67 static void *
  68 init_data_merge(int npar, gmx_ana_selparam_t *param);
  69 /*! \brief
  70  * Initializes data for the merging selection modifiers.
  71  *
  72  * \param[in] top   Not used.
  73  * \param[in] npar  Not used (should be 2 or 3).
  74  * \param[in] param Method parameters (should point to \ref smparams_merge).
  75  * \param[in] data  Should point to a \p t_methoddata_merge.
  76  * \returns   0 if everything is successful, -1 on error.
  77  */
  78 static void
  79 init_merge(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
  80 /** Initializes output for the \p merge selection modifier. */
  81 static void
  82 init_output_merge(t_topology *top, gmx_ana_selvalue_t *out, void *data);
  83 /** Initializes output for the \p plus selection modifier. */
  84 static void
  85 init_output_plus(t_topology *top, gmx_ana_selvalue_t *out, void *data);
  86 /** Frees the memory allocated for the merging selection modifiers. */
  87 static void
  88 free_data_merge(void *data);
  89 /*! \brief
  90  * Evaluates the \p merge selection modifier.
  91  *
  92  * \param[in]  top   Not used.
  93  * \param[in]  fr    Not used.
  94  * \param[in]  pbc   Not used.
  95  * \param[in]  p     Positions to merge (should point to \p data->p1).
  96  * \param[out] out   Output data structure (\p out->u.p is used).
  97  * \param[in]  data  Should point to a \p t_methoddata_merge.
  98  */
  99 static void
 100 evaluate_merge(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 101                gmx_ana_pos_t * p, gmx_ana_selvalue_t *out, void *data);
 102 /*! \brief
 103  * Evaluates the \p plus selection modifier.
 104  *
 105  * \param[in]  top   Not used.
 106  * \param[in]  fr    Not used.
 107  * \param[in]  pbc   Not used.
 108  * \param[in]  p     Positions to merge (should point to \p data->p1).
 109  * \param[out] out   Output data structure (\p out->u.p is used).
 110  * \param[in]  data  Should point to a \p t_methoddata_merge.
 111  */
 112 static void
 113 evaluate_plus(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 114               gmx_ana_pos_t * p, gmx_ana_selvalue_t *out, void *data);
 115
 116 /** Parameters for the merging selection modifiers. */
 117 static gmx_ana_selparam_t smparams_merge[] = {
 118     {NULL,       {POS_VALUE, -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_VARNUM},
 119     {NULL,       {POS_VALUE, -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_VARNUM},
 120     {"stride",   {INT_VALUE,  1, {NULL}}, NULL, SPAR_OPTIONAL},
 121 };
 122
 123 /** Help text for the merging selection modifiers. */
 124 static const char *help_merge[] = {
 125     "MERGING SELECTIONS[PAR]",
 126
 127     "[TT]POSEXPR merge POSEXPR [stride INT][tt][BR]",
 128     "[TT]POSEXPR merge POSEXPR [merge POSEXPR ...][tt][BR]",
 129     "[TT]POSEXPR plus POSEXPR [plus POSEXPR ...][tt][PAR]",
 130
 131     "Basic selection keywords can only create selections where each atom",
 132     "occurs at most once. The [TT]merge[tt] and [TT]plus[tt] selection",
 133     "keywords can be used to work around this limitation. Both create",
 134     "a selection that contains the positions from all the given position",
 135     "expressions, even if they contain duplicates.",
 136     "The difference between the two is that [TT]merge[tt] expects two or more",
 137     "selections with the same number of positions, and the output contains",
 138     "the input positions selected from each expression in turn, i.e.,",
 139     "the output is like A1 B1 A2 B2 and so on. It is also possible to merge",
 140     "selections of unequal size as long as the size of the first is a",
 141     "multiple of the second one. The [TT]stride[tt] parameter can be used",
 142     "to explicitly provide this multiplicity.",
 143     "[TT]plus[tt] simply concatenates the positions after each other, and",
 144     "can work also with selections of different sizes.",
 145     "These keywords are valid only at the selection level, not in any",
 146     "subexpressions.",
 147 };
 148
 149 /** Selection method data for the \p plus modifier. */
 150 gmx_ana_selmethod_t sm_merge = {
 151     "merge", POS_VALUE, SMETH_MODIFIER,
 152     asize(smparams_merge), smparams_merge,
 153     &init_data_merge,
 154     NULL,
 155     &init_merge,
 156     &init_output_merge,
 157     &free_data_merge,
 158     NULL,
 159     NULL,
 160     &evaluate_merge,
 161     {"merge POSEXPR", asize(help_merge), help_merge},
 162 };
 163
 164 /** Selection method data for the \p plus modifier. */
 165 gmx_ana_selmethod_t sm_plus = {
 166     "plus", POS_VALUE, SMETH_MODIFIER,
 167     asize(smparams_merge)-1, smparams_merge,
 168     &init_data_merge,
 169     NULL,
 170     &init_merge,
 171     &init_output_plus,
 172     &free_data_merge,
 173     NULL,
 174     NULL,
 175     &evaluate_plus,
 176     {"plus POSEXPR", asize(help_merge), help_merge},
 177 };
 178
 179 /*!
 180  * \param[in]     npar  Should be 2 for \c plus and 3 for \c merge.
 181  * \param[in,out] param Method parameters (should point to a copy of
 182  *   \ref smparams_merge).
 183  * \returns Pointer to the allocated data (\p t_methoddata_merge).
 184  *
 185  * Allocates memory for a \p t_methoddata_merge structure.
 186  */
 187 static void *
 188 init_data_merge(int npar, gmx_ana_selparam_t *param)
 189 {
 190     t_methoddata_merge *data = new t_methoddata_merge();
 191     data->stride     = 0;
 192     param[0].val.u.p = &data->p1;
 193     param[1].val.u.p = &data->p2;
 194     if (npar > 2)
 195     {
 196         param[2].val.u.i = &data->stride;
 197     }
 198     return data;
 199 }
 200
 201 static void
 202 init_merge(t_topology * /* top */, int /* npar */, gmx_ana_selparam_t * /* param */, void *data)
 203 {
 204     t_methoddata_merge *d = (t_methoddata_merge *)data;
 205
 206     if (d->stride < 0)
 207     {
 208         GMX_THROW(gmx::InvalidInputError("Stride for merging should be positive"));
 209     }
 210     /* If no stride given, deduce it from the input sizes */
 211     if (d->stride == 0)
 212     {
 213         d->stride = d->p1.count() / d->p2.count();
 214     }
 215     if (d->p1.count() != d->stride*d->p2.count())
 216     {
 217         GMX_THROW(gmx::InconsistentInputError("The number of positions to be merged are not compatible"));
 218     }
 219 }
 220
 221 /*! \brief
 222  * Does common initialization to all merging modifiers.
 223  *
 224  * \param[in]     top   Topology data structure.
 225  * \param[in,out] out   Pointer to output data structure.
 226  * \param[in,out] data  Should point to \c t_methoddata_merge.
 227  */
 228 static void
 229 init_output_common(t_topology *top, gmx_ana_selvalue_t *out, void *data)
 230 {
 231     t_methoddata_merge *d = (t_methoddata_merge *)data;
 232
 233     GMX_UNUSED_VALUE(top);
 234     if (d->p1.m.type != d->p2.m.type)
 235     {
 236         /* TODO: Maybe we could pick something else here? */
 237         out->u.p->m.type = INDEX_UNKNOWN;
 238     }
 239     else
 240     {
 241         out->u.p->m.type = d->p1.m.type;
 242     }
 243     gmx_ana_pos_reserve_for_append(out->u.p, d->p1.count() + d->p2.count(),
 244                                    d->p1.m.b.nra + d->p2.m.b.nra,
 245                                    d->p1.v != NULL, d->p1.f != NULL);
 246     gmx_ana_pos_empty_init(out->u.p);
 247 }
 248
 249 /*!
 250  * \param[in]     top   Topology data structure.
 251  * \param[in,out] out   Pointer to output data structure.
 252  * \param[in,out] data  Should point to \c t_methoddata_merge.
 253  */
 254 static void
 255 init_output_merge(t_topology *top, gmx_ana_selvalue_t *out, void *data)
 256 {
 257     t_methoddata_merge *d = (t_methoddata_merge *)data;
 258     int                 i, j;
 259
 260     init_output_common(top, out, data);
 261     for (i = 0; i < d->p2.count(); ++i)
 262     {
 263         for (j = 0; j < d->stride; ++j)
 264         {
 265             gmx_ana_pos_append_init(out->u.p, &d->p1, d->stride * i + j);
 266         }
 267         gmx_ana_pos_append_init(out->u.p, &d->p2, i);
 268     }
 269 }
 270
 271 /*!
 272  * \param[in]     top   Topology data structure.
 273  * \param[in,out] out   Pointer to output data structure.
 274  * \param[in,out] data  Should point to \c t_methoddata_merge.
 275  */
 276 static void
 277 init_output_plus(t_topology *top, gmx_ana_selvalue_t *out, void *data)
 278 {
 279     t_methoddata_merge *d = (t_methoddata_merge *)data;
 280     int                 i;
 281
 282     init_output_common(top, out, data);
 283     for (i = 0; i < d->p1.count(); ++i)
 284     {
 285         gmx_ana_pos_append_init(out->u.p, &d->p1, i);
 286     }
 287     for (i = 0; i < d->p2.count(); ++i)
 288     {
 289         gmx_ana_pos_append_init(out->u.p, &d->p2, i);
 290     }
 291 }
 292
 293 /*!
 294  * \param data Data to free (should point to a \p t_methoddata_merge).
 295  *
 296  * Frees the memory allocated for \c t_methoddata_merge.
 297  */
 298 static void
 299 free_data_merge(void *data)
 300 {
 301     t_methoddata_merge *d = (t_methoddata_merge *)data;
 302     delete d;
 303 }
 304
 305 static void
 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.count() != d->stride*d->p2.count())
 314     {
 315         GMX_THROW(gmx::InconsistentInputError("The number of positions to be merged are not compatible"));
 316     }
 317     gmx_ana_pos_empty(out->u.p);
 318     for (i = 0; i < d->p2.count(); ++i)
 319     {
 320         for (j = 0; j < d->stride; ++j)
 321         {
 322             refid = d->p1.m.refid[d->stride*i+j];
 323             if (refid != -1)
 324             {
 325                 refid = (d->stride+1) * (refid / d->stride) + (refid % d->stride);
 326             }
 327             gmx_ana_pos_append(out->u.p, &d->p1, d->stride*i+j, refid);
 328         }
 329         refid = (d->stride+1)*d->p2.m.refid[i] + d->stride;
 330         gmx_ana_pos_append(out->u.p, &d->p2, i, refid);
 331     }
 332     gmx_ana_pos_append_finish(out->u.p);
 333 }
 334
 335 static void
 336 evaluate_plus(t_topology * /* top */, t_trxframe * /* fr */, t_pbc * /* pbc */,
 337               gmx_ana_pos_t * /* p */, gmx_ana_selvalue_t *out, void *data)
 338 {
 339     t_methoddata_merge *d = (t_methoddata_merge *)data;
 340     int                 i;
 341     int                 refid;
 342
 343     gmx_ana_pos_empty(out->u.p);
 344     for (i = 0; i < d->p1.count(); ++i)
 345     {
 346         refid = d->p1.m.refid[i];
 347         gmx_ana_pos_append(out->u.p, &d->p1, i, refid);
 348     }
 349     for (i = 0; i < d->p2.count(); ++i)
 350     {
 351         refid = d->p2.m.refid[i];
 352         if (refid != -1)
 353         {
 354             refid += d->p1.m.b.nr;
 355         }
 356         gmx_ana_pos_append(out->u.p, &d->p2, i, refid);
 357     }
 358     gmx_ana_pos_append_finish(out->u.p);
 359 }