src/gmxlib/selection/sm_permute.c

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