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