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