src/gromacs/selection/sm_permute.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009,2010,2011,2012,2013,2014,2015, 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 permute selection modifier.
  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/exceptions.h"
  48 #include "gromacs/utility/smalloc.h"
  49 #include "gromacs/utility/stringutil.h"
  50
  51 #include "selmethod.h"
  52
  53 /*! \internal \brief
  54  * Data structure for the \p permute selection modifier.
  55  */
  56 typedef struct
  57 {
  58     /** Positions to permute. */
  59     gmx_ana_pos_t    p;
  60     /** Number of elements in the permutation. */
  61     int              n;
  62     /** Array describing the permutation. */
  63     int             *perm;
  64     /** Array that has the permutation reversed. */
  65     int             *rperm;
  66 } t_methoddata_permute;
  67
  68 /*! \brief
  69  * Allocates data for the \p permute selection modifier.
  70  *
  71  * \param[in]     npar  Not used (should be 2).
  72  * \param[in,out] param Method parameters (should point to a copy of
  73  *   \ref smparams_permute).
  74  * \returns Pointer to the allocated data (\p t_methoddata_permute).
  75  *
  76  * Allocates memory for a \p t_methoddata_permute structure.
  77  */
  78 static void *
  79 init_data_permute(int npar, gmx_ana_selparam_t *param);
  80 /*! \brief
  81  * Initializes data for the \p permute selection modifier.
  82  *
  83  * \param[in] top   Not used.
  84  * \param[in] npar  Not used (should be 2).
  85  * \param[in] param Method parameters (should point to \ref smparams_permute).
  86  * \param[in] data  Should point to a \p t_methoddata_permute.
  87  * \returns   0 if the input permutation is valid, -1 on error.
  88  */
  89 static void
  90 init_permute(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
  91 /*! \brief
  92  * Initializes output for the \p permute selection modifier.
  93  *
  94  * \param[in]     top   Topology data structure.
  95  * \param[in,out] out   Pointer to output data structure.
  96  * \param[in,out] data  Should point to \c t_methoddata_permute.
  97  */
  98 static void
  99 init_output_permute(t_topology *top, gmx_ana_selvalue_t *out, void *data);
 100 /** Frees the memory allocated for the \p permute selection modifier. */
 101 static void
 102 free_data_permute(void *data);
 103 static void
 104 /*! \brief
 105  * Evaluates the \p permute selection modifier.
 106  *
 107  * \param[in]  top   Not used.
 108  * \param[in]  fr    Not used.
 109  * \param[in]  pbc   Not used.
 110  * \param[in]  p     Positions to permute (should point to \p data->p).
 111  * \param[out] out   Output data structure (\p out->u.p is used).
 112  * \param[in]  data  Should point to a \p t_methoddata_permute.
 113  * \returns    0 if \p p could be permuted, -1 on error.
 114  *
 115  * Returns -1 if the size of \p p is not divisible by the number of
 116  * elements in the permutation.
 117  */
 118 evaluate_permute(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 119                  gmx_ana_pos_t *p, gmx_ana_selvalue_t *out, void *data);
 120
 121 /** Parameters for the \p permute selection modifier. */
 122 static gmx_ana_selparam_t smparams_permute[] = {
 123     {NULL,       {POS_VALUE, -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_VARNUM},
 124     {NULL,       {INT_VALUE, -1, {NULL}}, NULL, SPAR_VARNUM},
 125 };
 126
 127 /** Help text for the \p permute selection modifier. */
 128 static const char *const help_permute[] = {
 129     "::",
 130     "",
 131     "  permute P1 ... PN",
 132     "",
 133     "By default, all selections are evaluated such that the atom indices are",
 134     "returned in ascending order. This can be changed by appending",
 135     "[TT]permute P1 P2 ... PN[tt] to an expression.",
 136     "The [TT]Pi[tt] should form a permutation of the numbers 1 to N.",
 137     "This keyword permutes each N-position block in the selection such that",
 138     "the i'th position in the block becomes Pi'th.",
 139     "Note that it is the positions that are permuted, not individual atoms.",
 140     "A fatal error occurs if the size of the selection is not a multiple of n.",
 141     "It is only possible to permute the whole selection expression, not any",
 142     "subexpressions, i.e., the [TT]permute[tt] keyword should appear last in",
 143     "a selection.",
 144 };
 145
 146 /** Selection method data for the \p permute modifier. */
 147 gmx_ana_selmethod_t sm_permute = {
 148     "permute", POS_VALUE, SMETH_MODIFIER,
 149     asize(smparams_permute), smparams_permute,
 150     &init_data_permute,
 151     NULL,
 152     &init_permute,
 153     &init_output_permute,
 154     &free_data_permute,
 155     NULL,
 156     NULL,
 157     &evaluate_permute,
 158     {"POSEXPR permute P1 ... PN",
 159      "Permuting selections", asize(help_permute), help_permute},
 160 };
 161
 162 static void *
 163 init_data_permute(int /* npar */, gmx_ana_selparam_t *param)
 164 {
 165     t_methoddata_permute *data = new t_methoddata_permute();
 166     data->n          = 0;
 167     data->perm       = NULL;
 168     data->rperm      = NULL;
 169     param[0].val.u.p = &data->p;
 170     return data;
 171 }
 172
 173 static void
 174 init_permute(t_topology * /* top */, int /* npar */, gmx_ana_selparam_t *param, void *data)
 175 {
 176     t_methoddata_permute *d = (t_methoddata_permute *)data;
 177     int                   i;
 178
 179     d->n    = param[1].val.nr;
 180     d->perm = param[1].val.u.i;
 181     if (d->p.count() % d->n != 0)
 182     {
 183         GMX_THROW(gmx::InconsistentInputError(
 184                           gmx::formatString("The number of positions to be permuted is not divisible by %d", d->n)));
 185     }
 186     snew(d->rperm, d->n);
 187     for (i = 0; i < d->n; ++i)
 188     {
 189         d->rperm[i] = -1;
 190     }
 191     for (i = 0; i < d->n; ++i)
 192     {
 193         d->perm[i]--;
 194         if (d->perm[i] < 0 || d->perm[i] >= d->n)
 195         {
 196             GMX_THROW(gmx::InvalidInputError("Invalid permutation"));
 197         }
 198         if (d->rperm[d->perm[i]] >= 0)
 199         {
 200             GMX_THROW(gmx::InvalidInputError("Invalid permutation"));
 201         }
 202         d->rperm[d->perm[i]] = i;
 203     }
 204 }
 205
 206 static void
 207 init_output_permute(t_topology * /* top */, gmx_ana_selvalue_t *out, void *data)
 208 {
 209     t_methoddata_permute *d = (t_methoddata_permute *)data;
 210     int                   i, j, b;
 211
 212     out->u.p->m.type = d->p.m.type;
 213     gmx_ana_pos_reserve_for_append(out->u.p, d->p.count(), d->p.m.b.nra,
 214                                    d->p.v != NULL, d->p.f != NULL);
 215     gmx_ana_pos_empty_init(out->u.p);
 216     for (i = 0; i < d->p.count(); i += d->n)
 217     {
 218         for (j = 0; j < d->n; ++j)
 219         {
 220             b = i + d->rperm[j];
 221             gmx_ana_pos_append_init(out->u.p, &d->p, b);
 222         }
 223     }
 224 }
 225
 226 /*!
 227  * \param data Data to free (should point to a \p t_methoddata_permute).
 228  *
 229  * Frees the memory allocated for \c t_methoddata_permute.
 230  */
 231 static void
 232 free_data_permute(void *data)
 233 {
 234     t_methoddata_permute *d = (t_methoddata_permute *)data;
 235
 236     sfree(d->rperm);
 237     delete d;
 238 }
 239
 240 static void
 241 evaluate_permute(t_topology * /* top */, t_trxframe * /* fr */, t_pbc * /* pbc */,
 242                  gmx_ana_pos_t * /*p*/, gmx_ana_selvalue_t *out, void *data)
 243 {
 244     t_methoddata_permute *d = (t_methoddata_permute *)data;
 245     int                   i, j, b;
 246     int                   refid;
 247
 248     if (d->p.count() % d->n != 0)
 249     {
 250         GMX_THROW(gmx::InconsistentInputError(
 251                           gmx::formatString("The number of positions to be permuted is not divisible by %d", d->n)));
 252     }
 253     gmx_ana_pos_empty(out->u.p);
 254     for (i = 0; i < d->p.count(); i += d->n)
 255     {
 256         for (j = 0; j < d->n; ++j)
 257         {
 258             b     = i + d->rperm[j];
 259             refid = d->p.m.refid[b];
 260             if (refid != -1)
 261             {
 262                 /* De-permute the reference ID */
 263                 refid = refid - (refid % d->n) + d->perm[refid % d->n];
 264             }
 265             gmx_ana_pos_append(out->u.p, &d->p, b, refid);
 266         }
 267     }
 268     gmx_ana_pos_append_finish(out->u.p);
 269 }