src/gromacs/selection/sm_permute.cpp

   1 /*
   2  *
   3  *                This source code is part of
   4  *
   5  *                 G   R   O   M   A   C   S
   6  *
   7  *          GROningen MAchine for Chemical Simulations
   8  *
   9  * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
  10  * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
  11  * Copyright (c) 2001-2009, The GROMACS development team,
  12  * check out http://www.gromacs.org for more information.
  13
  14  * This program is free software; you can redistribute it and/or
  15  * modify it under the terms of the GNU General Public License
  16  * as published by the Free Software Foundation; either version 2
  17  * of the License, or (at your option) any later version.
  18  *
  19  * If you want to redistribute modifications, please consider that
  20  * scientific software is very special. Version control is crucial -
  21  * bugs must be traceable. We will be happy to consider code for
  22  * inclusion in the official distribution, but derived work must not
  23  * be called official GROMACS. Details are found in the README & COPYING
  24  * files - if they are missing, get the official version at www.gromacs.org.
  25  *
  26  * To help us fund GROMACS development, we humbly ask that you cite
  27  * the papers on the package - you can find them in the top README file.
  28  *
  29  * For more info, check our website at http://www.gromacs.org
  30  */
  31 /*! \internal \file
  32  * \brief
  33  * Implements the \p permute selection modifier.
  34  *
  35  * \author Teemu Murtola <teemu.murtola@cbr.su.se>
  36  * \ingroup module_selection
  37  */
  38 #include "gromacs/legacyheaders/macros.h"
  39 #include "gromacs/legacyheaders/smalloc.h"
  40 #include "gromacs/legacyheaders/vec.h"
  41
  42 #include "gromacs/selection/position.h"
  43 #include "gromacs/selection/selmethod.h"
  44 #include "gromacs/utility/exceptions.h"
  45 #include "gromacs/utility/stringutil.h"
  46
  47 /*! \internal \brief
  48  * Data structure for the \p permute selection modifier.
  49  */
  50 typedef struct
  51 {
  52     /** Positions to permute. */
  53     gmx_ana_pos_t    p;
  54     /** Group to receive the output permutation. */
  55     gmx_ana_index_t  g;
  56     /** Number of elements in the permutation. */
  57     int              n;
  58     /** Array describing the permutation. */
  59     int             *perm;
  60     /** Array that has the permutation reversed. */
  61     int             *rperm;
  62 } t_methoddata_permute;
  63
  64 /** Allocates data for the \p permute selection modifier. */
  65 static void *
  66 init_data_permute(int npar, gmx_ana_selparam_t *param);
  67 /** Initializes data for the \p permute selection modifier. */
  68 static void
  69 init_permute(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
  70 /** Initializes output for the \p permute selection modifier. */
  71 static void
  72 init_output_permute(t_topology *top, gmx_ana_selvalue_t *out, void *data);
  73 /** Frees the memory allocated for the \p permute selection modifier. */
  74 static void
  75 free_data_permute(void *data);
  76 /** Evaluates the \p permute selection modifier. */
  77 static void
  78 evaluate_permute(t_topology *top, t_trxframe *fr, t_pbc *pbc,
  79                  gmx_ana_pos_t *p, gmx_ana_selvalue_t *out, void *data);
  80
  81 /** Parameters for the \p permute selection modifier. */
  82 static gmx_ana_selparam_t smparams_permute[] = {
  83     {NULL,       {POS_VALUE, -1, {NULL}}, NULL, SPAR_DYNAMIC | SPAR_VARNUM},
  84     {NULL,       {INT_VALUE, -1, {NULL}}, NULL, SPAR_VARNUM},
  85 };
  86
  87 /** Help text for the \p permute selection modifier. */
  88 static const char *help_permute[] = {
  89     "PERMUTING SELECTIONS[PAR]",
  90
  91     "[TT]permute P1 ... PN[tt][PAR]",
  92
  93     "By default, all selections are evaluated such that the atom indices are",
  94     "returned in ascending order. This can be changed by appending",
  95     "[TT]permute P1 P2 ... PN[tt] to an expression.",
  96     "The [TT]Pi[tt] should form a permutation of the numbers 1 to N.",
  97     "This keyword permutes each N-position block in the selection such that",
  98     "the i'th position in the block becomes Pi'th.",
  99     "Note that it is the positions that are permuted, not individual atoms.",
 100     "A fatal error occurs if the size of the selection is not a multiple of n.",
 101     "It is only possible to permute the whole selection expression, not any",
 102     "subexpressions, i.e., the [TT]permute[tt] keyword should appear last in",
 103     "a selection.",
 104 };
 105
 106 /** \internal Selection method data for the \p permute modifier. */
 107 gmx_ana_selmethod_t sm_permute = {
 108     "permute", POS_VALUE, SMETH_MODIFIER,
 109     asize(smparams_permute), smparams_permute,
 110     &init_data_permute,
 111     NULL,
 112     &init_permute,
 113     &init_output_permute,
 114     &free_data_permute,
 115     NULL,
 116     NULL,
 117     &evaluate_permute,
 118     {"permute P1 ... PN", asize(help_permute), help_permute},
 119 };
 120
 121 /*!
 122  * \param[in]     npar  Not used (should be 2).
 123  * \param[in,out] param Method parameters (should point to a copy of
 124  *   \ref smparams_permute).
 125  * \returns Pointer to the allocated data (\p t_methoddata_permute).
 126  *
 127  * Allocates memory for a \p t_methoddata_permute structure.
 128  */
 129 static void *
 130 init_data_permute(int npar, gmx_ana_selparam_t *param)
 131 {
 132     t_methoddata_permute *data;
 133
 134     snew(data, 1);
 135     param[0].val.u.p = &data->p;
 136     return data;
 137 }
 138
 139 /*!
 140  * \param[in] top   Not used.
 141  * \param[in] npar  Not used (should be 2).
 142  * \param[in] param Method parameters (should point to \ref smparams_permute).
 143  * \param[in] data  Should point to a \p t_methoddata_permute.
 144  * \returns   0 if the input permutation is valid, -1 on error.
 145  */
 146 static void
 147 init_permute(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data)
 148 {
 149     t_methoddata_permute *d = (t_methoddata_permute *)data;
 150     int                   i;
 151
 152     gmx_ana_index_reserve(&d->g, d->p.g->isize);
 153     d->n    = param[1].val.nr;
 154     d->perm = param[1].val.u.i;
 155     if (d->p.nr % d->n != 0)
 156     {
 157         GMX_THROW(gmx::InconsistentInputError(
 158                     gmx::formatString("The number of positions to be permuted is not divisible by %d", d->n)));
 159     }
 160     snew(d->rperm, d->n);
 161     for (i = 0; i < d->n; ++i)
 162     {
 163         d->rperm[i] = -1;
 164     }
 165     for (i = 0; i < d->n; ++i)
 166     {
 167         d->perm[i]--;
 168         if (d->perm[i] < 0 || d->perm[i] >= d->n)
 169         {
 170             GMX_THROW(gmx::InvalidInputError("Invalid permutation"));
 171         }
 172         if (d->rperm[d->perm[i]] >= 0)
 173         {
 174             GMX_THROW(gmx::InvalidInputError("Invalid permutation"));
 175         }
 176         d->rperm[d->perm[i]] = i;
 177     }
 178 }
 179
 180 /*!
 181  * \param[in]     top   Topology data structure.
 182  * \param[in,out] out   Pointer to output data structure.
 183  * \param[in,out] data  Should point to \c t_methoddata_permute.
 184  */
 185 static void
 186 init_output_permute(t_topology *top, gmx_ana_selvalue_t *out, void *data)
 187 {
 188     t_methoddata_permute *d = (t_methoddata_permute *)data;
 189     int                   i, j, b;
 190
 191     gmx_ana_pos_copy(out->u.p, &d->p, true);
 192     gmx_ana_pos_set_evalgrp(out->u.p, &d->g);
 193     d->g.isize = 0;
 194     gmx_ana_pos_empty_init(out->u.p);
 195     for (i = 0; i < d->p.nr; i += d->n)
 196     {
 197         for (j = 0; j < d->n; ++j)
 198         {
 199             b = i + d->rperm[j];
 200             gmx_ana_pos_append_init(out->u.p, &d->g, &d->p, b);
 201         }
 202     }
 203 }
 204
 205 /*!
 206  * \param data Data to free (should point to a \p t_methoddata_permute).
 207  *
 208  * Frees the memory allocated for \c t_methoddata_permute.
 209  */
 210 static void
 211 free_data_permute(void *data)
 212 {
 213     t_methoddata_permute *d = (t_methoddata_permute *)data;
 214
 215     gmx_ana_index_deinit(&d->g);
 216     sfree(d->rperm);
 217     sfree(d);
 218 }
 219
 220 /*!
 221  * \param[in]  top   Not used.
 222  * \param[in]  fr    Not used.
 223  * \param[in]  pbc   Not used.
 224  * \param[in]  p     Positions to permute (should point to \p data->p).
 225  * \param[out] out   Output data structure (\p out->u.p is used).
 226  * \param[in]  data  Should point to a \p t_methoddata_permute.
 227  * \returns    0 if \p p could be permuted, -1 on error.
 228  *
 229  * Returns -1 if the size of \p p is not divisible by the number of
 230  * elements in the permutation.
 231  */
 232 static void
 233 evaluate_permute(t_topology *top, t_trxframe *fr, t_pbc *pbc,
 234                  gmx_ana_pos_t *p, gmx_ana_selvalue_t *out, void *data)
 235 {
 236     t_methoddata_permute *d = (t_methoddata_permute *)data;
 237     int                   i, j, b;
 238     int                   refid;
 239
 240     if (d->p.nr % d->n != 0)
 241     {
 242         GMX_THROW(gmx::InconsistentInputError(
 243                     gmx::formatString("The number of positions to be permuted is not divisible by %d", d->n)));
 244     }
 245     d->g.isize = 0;
 246     gmx_ana_pos_empty(out->u.p);
 247     for (i = 0; i < d->p.nr; 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, &d->g, p, b, refid);
 259         }
 260     }
 261     gmx_ana_pos_append_finish(out->u.p);
 262 }