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