/*
+ * This file is part of the GROMACS molecular simulation package.
*
- * This source code is part of
+ * Copyright (c) 2009,2010,2011,2012,2013,2014, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
*
- * G R O M A C S
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
*
- * GROningen MAchine for Chemical Simulations
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
*
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2009, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
*
* To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
+ * the research papers on the package. Check out http://www.gromacs.org.
*/
/*! \internal \file
* \brief
* Implements functions in indexutil.h.
*
- * \author Teemu Murtola <teemu.murtola@cbr.su.se>
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
* \ingroup module_selection
*/
-#include "gromacs/legacyheaders/index.h"
-#include "gromacs/legacyheaders/gmx_fatal.h"
-#include "gromacs/legacyheaders/smalloc.h"
-#include "gromacs/legacyheaders/string2.h"
-#include "gromacs/legacyheaders/typedefs.h"
+#include "gmxpre.h"
+
+#include "indexutil.h"
-#include "gromacs/selection/indexutil.h"
+#include <cstdlib>
+#include <cstring>
+
+#include <algorithm>
+#include <string>
+#include <vector>
+
+#include "gromacs/topology/block.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/smalloc.h"
/********************************************************************
* gmx_ana_indexgrps_t functions
*/
struct gmx_ana_indexgrps_t
{
+ //! Initializes an empty set of groups.
+ explicit gmx_ana_indexgrps_t(int nr) : nr(nr), g(NULL)
+ {
+ names.reserve(nr);
+ snew(g, nr);
+ }
+ ~gmx_ana_indexgrps_t()
+ {
+ for (int i = 0; i < nr; ++i)
+ {
+ gmx_ana_index_deinit(&g[i]);
+ }
+ sfree(g);
+ }
+
/** Number of index groups. */
- int nr;
+ int nr;
/** Array of index groups. */
- gmx_ana_index_t *g;
+ gmx_ana_index_t *g;
+ /** Group names. */
+ std::vector<std::string> names;
};
-/*!
- * \param[out] g Index group structure.
- * \param[in] ngrps Number of groups for which memory is allocated.
- */
-void
-gmx_ana_indexgrps_alloc(gmx_ana_indexgrps_t **g, int ngrps)
-{
- snew(*g, 1);
- (*g)->nr = ngrps;
- snew((*g)->g, ngrps);
-}
-
/*!
* \param[out] g Index group structure.
* \param[in] top Topology structure.
{
t_blocka *block = NULL;
char **names = NULL;
- int i, j;
if (fnm)
{
}
else
{
- snew(*g, 1);
- (*g)->nr = 0;
- (*g)->g = NULL;
+ *g = new gmx_ana_indexgrps_t(0);
return;
}
- gmx_ana_indexgrps_alloc(g, block->nr);
- for (i = 0; i < block->nr; ++i)
+ try
{
- gmx_ana_index_t *grp = &(*g)->g[i];
+ *g = new gmx_ana_indexgrps_t(block->nr);
+ for (int i = 0; i < block->nr; ++i)
+ {
+ gmx_ana_index_t *grp = &(*g)->g[i];
- grp->isize = block->index[i+1] - block->index[i];
- snew(grp->index, grp->isize);
- for (j = 0; j < grp->isize; ++j)
+ grp->isize = block->index[i+1] - block->index[i];
+ snew(grp->index, grp->isize);
+ for (int j = 0; j < grp->isize; ++j)
+ {
+ grp->index[j] = block->a[block->index[i]+j];
+ }
+ grp->nalloc_index = grp->isize;
+ (*g)->names.push_back(names[i]);
+ }
+ }
+ catch (...)
+ {
+ for (int i = 0; i < block->nr; ++i)
{
- grp->index[j] = block->a[block->index[i]+j];
+ sfree(names[i]);
}
- grp->name = names[i];
- grp->nalloc_index = grp->isize;
+ sfree(names);
+ done_blocka(block);
+ sfree(block);
+ throw;
}
-
+ for (int i = 0; i < block->nr; ++i)
+ {
+ sfree(names[i]);
+ }
+ sfree(names);
done_blocka(block);
sfree(block);
- sfree(names);
}
/*!
void
gmx_ana_indexgrps_free(gmx_ana_indexgrps_t *g)
{
- int i;
-
- if (g->nr == 0)
- {
- sfree(g);
- return;
- }
- for (i = 0; i < g->nr; ++i)
- {
- gmx_ana_index_deinit(&g->g[i]);
- }
- sfree(g->g);
- g->nr = 0;
- g->g = NULL;
- sfree(g);
-}
-
-/*!
- * \param[out] dest Destination index groups.
- * \param[in] src Source index groups.
- *
- * A deep copy is made for all fields, including the group names.
- */
-void
-gmx_ana_indexgrps_clone(gmx_ana_indexgrps_t **dest, gmx_ana_indexgrps_t *src)
-{
- int g;
-
- gmx_ana_indexgrps_alloc(dest, src->nr);
- for (g = 0; g < src->nr; ++g)
- {
- gmx_ana_index_copy(&(*dest)->g[g], &src->g[g], true);
- }
+ delete g;
}
/*!
}
/*!
- * \param[out] dest Output structure.
- * \param[in] src Input index groups.
- * \param[in] n Number of the group to extract.
+ * \param[out] dest Output structure.
+ * \param[out] destName Receives the name of the group if found.
+ * \param[in] src Input index groups.
+ * \param[in] n Number of the group to extract.
* \returns true if \p n is a valid group in \p src, false otherwise.
*/
bool
-gmx_ana_indexgrps_extract(gmx_ana_index_t *dest, gmx_ana_indexgrps_t *src, int n)
+gmx_ana_indexgrps_extract(gmx_ana_index_t *dest, std::string *destName,
+ gmx_ana_indexgrps_t *src, int n)
{
+ destName->clear();
if (n < 0 || n >= src->nr)
{
dest->isize = 0;
return false;
}
+ if (destName != NULL)
+ {
+ *destName = src->names[n];
+ }
gmx_ana_index_copy(dest, &src->g[n], true);
return true;
}
/*!
- * \param[out] dest Output structure.
- * \param[in] src Input index groups.
- * \param[in] name Name (or part of the name) of the group to extract.
+ * \param[out] dest Output structure.
+ * \param[out] destName Receives the name of the group if found.
+ * \param[in] src Input index groups.
+ * \param[in] name Name (or part of the name) of the group to extract.
* \returns true if \p name is a valid group in \p src, false otherwise.
*
* Uses the Gromacs routine find_group() to find the actual group;
* the comparison is case-insensitive.
*/
bool
-gmx_ana_indexgrps_find(gmx_ana_index_t *dest, gmx_ana_indexgrps_t *src, char *name)
+gmx_ana_indexgrps_find(gmx_ana_index_t *dest, std::string *destName,
+ gmx_ana_indexgrps_t *src,
+ const char *name)
{
- int i;
- char **names;
+ const char **names;
+ destName->clear();
snew(names, src->nr);
- for (i = 0; i < src->nr; ++i)
+ for (int i = 0; i < src->nr; ++i)
{
- names[i] = src->g[i].name;
+ names[i] = src->names[i].c_str();
}
- i = find_group(name, src->nr, names);
+ int n = find_group(const_cast<char *>(name), src->nr,
+ const_cast<char **>(names));
sfree(names);
- if (i == NOTSET)
+ if (n < 0)
{
dest->isize = 0;
return false;
}
- return gmx_ana_indexgrps_extract(dest, src, i);
+ return gmx_ana_indexgrps_extract(dest, destName, src, n);
}
/*!
void
gmx_ana_indexgrps_print(FILE *fp, gmx_ana_indexgrps_t *g, int maxn)
{
- int i;
-
- for (i = 0; i < g->nr; ++i)
+ for (int i = 0; i < g->nr; ++i)
{
- fprintf(fp, " %2d: ", i);
- gmx_ana_index_dump(fp, &g->g[i], i, maxn);
+ fprintf(fp, " Group %2d \"%s\" ", i, g->names[i].c_str());
+ gmx_ana_index_dump(fp, &g->g[i], maxn);
}
}
{
g->isize = 0;
g->index = NULL;
- g->name = NULL;
g->nalloc_index = 0;
}
* \param[out] g Output structure.
* \param[in] isize Number of atoms in the new group.
* \param[in] index Array of \p isize atoms (can be NULL if \p isize is 0).
- * \param[in] name Name for the new group (can be NULL).
* \param[in] nalloc Number of elements allocated for \p index
* (if 0, \p index is not freed in gmx_ana_index_deinit())
*
* No copy if \p index is made.
*/
void
-gmx_ana_index_set(gmx_ana_index_t *g, int isize, atom_id *index, char *name,
- int nalloc)
+gmx_ana_index_set(gmx_ana_index_t *g, int isize, atom_id *index, int nalloc)
{
g->isize = isize;
g->index = index;
- g->name = name;
g->nalloc_index = nalloc;
}
/*!
* \param[out] g Output structure.
* \param[in] natoms Number of atoms.
- * \param[in] name Name for the new group (can be NULL).
*/
void
-gmx_ana_index_init_simple(gmx_ana_index_t *g, int natoms, char *name)
+gmx_ana_index_init_simple(gmx_ana_index_t *g, int natoms)
{
int i;
{
g->index[i] = i;
}
- g->name = name;
g->nalloc_index = natoms;
}
{
sfree(g->index);
}
- sfree(g->name);
gmx_ana_index_clear(g);
}
snew(dest->index, dest->isize);
dest->nalloc_index = dest->isize;
}
- memcpy(dest->index, src->index, dest->isize*sizeof(*dest->index));
- }
- if (bAlloc && src->name)
- {
- dest->name = strdup(src->name);
- }
- else if (bAlloc || src->name)
- {
- dest->name = src->name;
+ std::memcpy(dest->index, src->index, dest->isize*sizeof(*dest->index));
}
}
/*!
* \param[in] fp Where to print the output.
* \param[in] g Index group to print.
- * \param[in] i Group number to use if the name is NULL.
* \param[in] maxn Maximum number of indices to print (-1 = print all).
*/
void
-gmx_ana_index_dump(FILE *fp, gmx_ana_index_t *g, int i, int maxn)
+gmx_ana_index_dump(FILE *fp, gmx_ana_index_t *g, int maxn)
{
int j, n;
- if (g->name)
- {
- fprintf(fp, "\"%s\"", g->name);
- }
- else
- {
- fprintf(fp, "Group %d", i+1);
- }
- fprintf(fp, " (%d atoms)", g->isize);
+ fprintf(fp, "(%d atoms)", g->isize);
if (maxn != 0)
{
fprintf(fp, ":");
fprintf(fp, "\n");
}
-/*!
- * \param[in] g Input index group.
- * \param[in] natoms Number of atoms to check against.
- *
- * If any atom index in the index group is less than zero or >= \p natoms,
- * gmx_fatal() is called.
- */
-void
-gmx_ana_index_check(gmx_ana_index_t *g, int natoms)
+int
+gmx_ana_index_get_max_index(gmx_ana_index_t *g)
{
- int j;
-
- for (j = 0; j < g->isize; ++j)
+ if (g->isize == 0)
{
- if (g->index[j] >= natoms)
- {
- gmx_fatal(FARGS, "Atom index (%d) in index group %s (%d atoms) "
- "larger than number of atoms in trajectory (%d atoms)",
- g->index[j], g->name, g->isize, natoms);
- }
- else if (g->index[j] < 0)
- {
- gmx_fatal(FARGS, "Atom index (%d) in index group %s (%d atoms) "
- "is less than zero",
- g->index[j], g->name, g->isize);
- }
+ return 0;
+ }
+ else
+ {
+ return *std::max_element(g->index, g->index + g->isize);
}
}
return true;
}
+bool
+gmx_ana_index_check_range(gmx_ana_index_t *g, int natoms)
+{
+ for (int i = 0; i < g->isize; ++i)
+ {
+ if (g->index[i] < 0 || g->index[i] >= natoms)
+ {
+ return false;
+ }
+ }
+ return true;
+}
+
/********************************************************************
* Set operations
********************************************************************/
void
gmx_ana_index_sort(gmx_ana_index_t *g)
{
- qsort(g->index, g->isize, sizeof(*g->index), cmp_atomid);
+ std::qsort(g->index, g->isize, sizeof(*g->index), cmp_atomid);
}
/*!
srenew(t->a, g->isize);
t->nalloc_a = g->isize;
}
- memcpy(t->a, g->index, g->isize*sizeof(*(t->a)));
+ std::memcpy(t->a, g->index, g->isize*sizeof(*(t->a)));
}
/* Allocate memory for the block index. We don't know in advance
break;
default: /* Should not be reached */
- gmx_bug("internal error");
+ GMX_RELEASE_ASSERT(false, "Unreachable code was reached");
break;
}
}
* \returns true if \p g consists of one or more complete elements of type
* \p type, false otherwise.
*
+ * \p g is assumed to be sorted, otherwise may return false negatives.
+ *
* If \p type is \ref INDEX_ATOM, the return value is always true.
* If \p type is \ref INDEX_UNKNOWN or \ref INDEX_ALL, the return value is
* always false.
gmx_ana_index_has_complete_elems(gmx_ana_index_t *g, e_index_t type,
t_topology *top)
{
+ // TODO: Consider whether unsorted groups need to be supported better.
switch (type)
{
case INDEX_UNKNOWN:
gmx_ana_indexmap_clear(gmx_ana_indexmap_t *m)
{
m->type = INDEX_UNKNOWN;
- m->nr = 0;
m->refid = NULL;
m->mapid = NULL;
m->mapb.nr = 0;
m->mapb.index = NULL;
m->mapb.nalloc_index = 0;
+ m->mapb.nra = 0;
+ m->mapb.a = NULL;
+ m->mapb.nalloc_a = 0;
m->orgid = NULL;
m->b.nr = 0;
m->b.index = NULL;
m->b.nalloc_index = 0;
m->b.nalloc_a = 0;
m->bStatic = true;
- m->bMapStatic = true;
}
/*!
m->type = type;
gmx_ana_index_make_block(&m->b, top, g, type, false);
gmx_ana_indexmap_reserve(m, m->b.nr, m->b.nra);
- m->nr = m->b.nr;
- for (i = mi = 0; i < m->nr; ++i)
+ for (i = mi = 0; i < m->b.nr; ++i)
{
ii = (type == INDEX_UNKNOWN ? 0 : m->b.a[m->b.index[i]]);
switch (type)
break;
}
}
- for (i = 0; i < m->nr; ++i)
+ for (i = 0; i < m->b.nr; ++i)
{
m->refid[i] = i;
m->mapid[i] = m->orgid[i];
}
- m->mapb.nr = m->nr;
- memcpy(m->mapb.index, m->b.index, (m->nr+1)*sizeof(*(m->mapb.index)));
- m->bStatic = true;
- m->bMapStatic = true;
+ m->mapb.nr = m->b.nr;
+ m->mapb.nra = m->b.nra;
+ m->mapb.a = m->b.a;
+ std::memcpy(m->mapb.index, m->b.index, (m->b.nr+1)*sizeof(*(m->mapb.index)));
+ m->bStatic = true;
}
/*!
gmx_ana_indexmap_set_static(gmx_ana_indexmap_t *m, t_blocka *b)
{
sfree(m->mapid);
- m->mapid = m->orgid;
- sfree(m->b.index);
- m->b.nalloc_index = 0;
- m->b.index = b->index;
sfree(m->mapb.index);
- m->mapb.nalloc_index = 0;
- m->mapb.index = m->b.index;
+ sfree(m->b.index);
sfree(m->b.a);
- m->b.nalloc_a = 0;
- m->b.a = b->a;
+ m->mapb.nalloc_index = 0;
+ m->mapb.nalloc_a = 0;
+ m->b.nalloc_index = 0;
+ m->b.nalloc_a = 0;
+ m->mapid = m->orgid;
+ m->mapb.index = b->index;
+ m->mapb.a = b->a;
+ m->b.index = b->index;
+ m->b.a = b->a;
}
/*!
dest->type = src->type;
dest->b.nr = src->b.nr;
dest->b.nra = src->b.nra;
- memcpy(dest->orgid, src->orgid, dest->b.nr*sizeof(*dest->orgid));
- memcpy(dest->b.index, src->b.index, (dest->b.nr+1)*sizeof(*dest->b.index));
- memcpy(dest->b.a, src->b.a, dest->b.nra*sizeof(*dest->b.a));
+ std::memcpy(dest->orgid, src->orgid, dest->b.nr*sizeof(*dest->orgid));
+ std::memcpy(dest->b.index, src->b.index, (dest->b.nr+1)*sizeof(*dest->b.index));
+ std::memcpy(dest->b.a, src->b.a, dest->b.nra*sizeof(*dest->b.a));
}
- dest->nr = src->nr;
dest->mapb.nr = src->mapb.nr;
- memcpy(dest->refid, src->refid, dest->nr*sizeof(*dest->refid));
- memcpy(dest->mapid, src->mapid, dest->nr*sizeof(*dest->mapid));
- memcpy(dest->mapb.index, src->mapb.index, (dest->mapb.nr+1)*sizeof(*dest->mapb.index));
- dest->bStatic = src->bStatic;
- dest->bMapStatic = src->bMapStatic;
+ dest->mapb.nra = src->mapb.nra;
+ if (src->mapb.nalloc_a > 0)
+ {
+ if (bFirst)
+ {
+ snew(dest->mapb.a, src->mapb.nalloc_a);
+ dest->mapb.nalloc_a = src->mapb.nalloc_a;
+ }
+ std::memcpy(dest->mapb.a, src->mapb.a, dest->mapb.nra*sizeof(*dest->mapb.a));
+ }
+ else
+ {
+ dest->mapb.a = src->mapb.a;
+ }
+ std::memcpy(dest->refid, src->refid, dest->mapb.nr*sizeof(*dest->refid));
+ std::memcpy(dest->mapid, src->mapid, dest->mapb.nr*sizeof(*dest->mapid));
+ std::memcpy(dest->mapb.index, src->mapb.index, (dest->mapb.nr+1)*sizeof(*dest->mapb.index));
+ dest->bStatic = src->bStatic;
+}
+
+/*! \brief
+ * Helper function to set the source atoms in an index map.
+ *
+ * \param[in,out] m Mapping structure.
+ * \param[in] isize Number of atoms in the \p index array.
+ * \param[in] index List of atoms.
+ */
+static void
+set_atoms(gmx_ana_indexmap_t *m, int isize, int *index)
+{
+ m->mapb.nra = isize;
+ if (m->mapb.nalloc_a == 0)
+ {
+ m->mapb.a = index;
+ }
+ else
+ {
+ for (int i = 0; i < isize; ++i)
+ {
+ m->mapb.a[i] = index[i];
+ }
+ }
}
/*!
bool bMaskOnly)
{
int i, j, bi, bj;
- bool bStatic;
/* Process the simple cases first */
if (m->type == INDEX_UNKNOWN && m->b.nra == 0)
}
if (m->type == INDEX_ALL)
{
+ set_atoms(m, g->isize, g->index);
if (m->b.nr > 0)
{
m->mapb.index[1] = g->isize;
return;
}
/* Reset the reference IDs and mapping if necessary */
- bStatic = (g->isize == m->b.nra && m->nr == m->b.nr);
- if (bStatic || bMaskOnly)
+ const bool bToFull = (g->isize == m->b.nra);
+ const bool bWasFull = (m->mapb.nra == m->b.nra);
+ if (bToFull || bMaskOnly)
{
if (!m->bStatic)
{
m->refid[bj] = bj;
}
}
- if (!m->bMapStatic)
+ if (!bWasFull)
{
for (bj = 0; bj < m->b.nr; ++bj)
{
{
m->mapb.index[bj] = m->b.index[bj];
}
- m->bMapStatic = true;
}
+ set_atoms(m, m->b.nra, m->b.a);
+ m->mapb.nr = m->b.nr;
}
/* Exit immediately if the group is static */
- if (bStatic)
+ if (bToFull)
{
m->bStatic = true;
return;
if (bMaskOnly)
{
- m->nr = m->b.nr;
for (i = j = bj = 0; i < g->isize; ++i, ++j)
{
/* Find the next atom in the block */
}
else
{
+ set_atoms(m, g->isize, g->index);
for (i = j = bi = 0, bj = -1; i < g->isize; ++i)
{
/* Find the next atom in the block */
}
/* Update the number of blocks */
m->mapb.index[bi] = g->isize;
- m->nr = bi;
- m->bMapStatic = false;
+ m->mapb.nr = bi;
}
- m->mapb.nr = m->nr;
m->bStatic = false;
}
{
sfree(m->mapb.index);
}
+ if (m->mapb.nalloc_a > 0)
+ {
+ sfree(m->mapb.a);
+ }
sfree(m->orgid);
if (m->b.nalloc_index > 0)
{