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