Sort all includes in src/gromacs
[alexxy/gromacs.git] / src / gromacs / selection / indexutil.cpp
index 819340f2a8f35a5259d091b32d15d3c4f3b0cd67..a38fe3620509088df357a6f44458da298d89c430 100644 (file)
@@ -1,10 +1,10 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2009,2010,2011,2012,2013, by the GROMACS development team, led by
- * David van der Spoel, Berk Hess, Erik Lindahl, and including many
- * others, as listed in the AUTHORS file in the top-level source
- * directory and at http://www.gromacs.org.
+ * 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.
  *
  * GROMACS is free software; you can redistribute it and/or
  * modify it under the terms of the GNU Lesser General Public License
  * \author Teemu Murtola <teemu.murtola@gmail.com>
  * \ingroup module_selection
  */
-#include "gromacs/selection/indexutil.h"
+#include "gmxpre.h"
 
+#include "indexutil.h"
+
+#include <cstdlib>
 #include <cstring>
 
+#include <algorithm>
 #include <string>
 #include <vector>
 
-#include "gromacs/legacyheaders/index.h"
-#include "gromacs/legacyheaders/gmx_fatal.h"
-#include "gromacs/legacyheaders/smalloc.h"
-#include "gromacs/legacyheaders/typedefs.h"
-
+#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
@@ -138,22 +142,22 @@ gmx_ana_indexgrps_init(gmx_ana_indexgrps_t **g, t_topology *top,
     }
     catch (...)
     {
-        done_blocka(block);
-        sfree(block);
         for (int i = 0; i < block->nr; ++i)
         {
             sfree(names[i]);
         }
         sfree(names);
+        done_blocka(block);
+        sfree(block);
         throw;
     }
-    done_blocka(block);
-    sfree(block);
     for (int i = 0; i < block->nr; ++i)
     {
         sfree(names[i]);
     }
     sfree(names);
+    done_blocka(block);
+    sfree(block);
 }
 
 /*!
@@ -246,7 +250,7 @@ gmx_ana_indexgrps_find(gmx_ana_index_t *dest, std::string *destName,
     int n = find_group(const_cast<char *>(name), src->nr,
                        const_cast<char **>(names));
     sfree(names);
-    if (n == NOTSET)
+    if (n < 0)
     {
         dest->isize = 0;
         return false;
@@ -266,7 +270,7 @@ gmx_ana_indexgrps_print(FILE *fp, gmx_ana_indexgrps_t *g, int maxn)
 {
     for (int i = 0; i < g->nr; ++i)
     {
-        fprintf(fp, " Group %2d \"%s\" ", i + 1, g->names[i].c_str());
+        fprintf(fp, " Group %2d \"%s\" ", i, g->names[i].c_str());
         gmx_ana_index_dump(fp, &g->g[i], maxn);
     }
 }
@@ -419,6 +423,19 @@ gmx_ana_index_dump(FILE *fp, gmx_ana_index_t *g, int maxn)
     fprintf(fp, "\n");
 }
 
+int
+gmx_ana_index_get_max_index(gmx_ana_index_t *g)
+{
+    if (g->isize == 0)
+    {
+        return 0;
+    }
+    else
+    {
+        return *std::max_element(g->index, g->index + g->isize);
+    }
+}
+
 /*!
  * \param[in]  g      Index group to check.
  * \returns    true if the index group is sorted and has no duplicates,
@@ -439,6 +456,19 @@ gmx_ana_index_check_sorted(gmx_ana_index_t *g)
     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
  ********************************************************************/
@@ -464,7 +494,7 @@ cmp_atomid(const void *a, const void *b)
 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);
 }
 
 /*!
@@ -844,7 +874,7 @@ gmx_ana_index_make_block(t_blocka *t, t_topology *top, gmx_ana_index_t *g,
                         break;
 
                     default: /* Should not be reached */
-                        gmx_bug("internal error");
+                        GMX_RELEASE_ASSERT(false, "Unreachable code was reached");
                         break;
                 }
             }
@@ -1024,12 +1054,14 @@ void
 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;
@@ -1038,7 +1070,6 @@ gmx_ana_indexmap_clear(gmx_ana_indexmap_t *m)
     m->b.nalloc_index    = 0;
     m->b.nalloc_a        = 0;
     m->bStatic           = true;
-    m->bMapStatic        = true;
 }
 
 /*!
@@ -1089,8 +1120,7 @@ gmx_ana_indexmap_init(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
     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)
@@ -1114,15 +1144,16 @@ gmx_ana_indexmap_init(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
                 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;
-    std::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;
 }
 
 /*!
@@ -1144,16 +1175,18 @@ void
 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;
 }
 
 /*!
@@ -1178,13 +1211,49 @@ gmx_ana_indexmap_copy(gmx_ana_indexmap_t *dest, gmx_ana_indexmap_t *src, bool bF
         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;
-    std::memcpy(dest->refid,      src->refid,      dest->nr*sizeof(*dest->refid));
-    std::memcpy(dest->mapid,      src->mapid,      dest->nr*sizeof(*dest->mapid));
+    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;
-    dest->bMapStatic = src->bMapStatic;
+    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];
+        }
+    }
 }
 
 /*!
@@ -1202,7 +1271,6 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
                         bool bMaskOnly)
 {
     int  i, j, bi, bj;
-    bool bStatic;
 
     /* Process the simple cases first */
     if (m->type == INDEX_UNKNOWN && m->b.nra == 0)
@@ -1211,6 +1279,7 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
     }
     if (m->type == INDEX_ALL)
     {
+        set_atoms(m, g->isize, g->index);
         if (m->b.nr > 0)
         {
             m->mapb.index[1] = g->isize;
@@ -1218,8 +1287,9 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
         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)
         {
@@ -1228,7 +1298,7 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
                 m->refid[bj] = bj;
             }
         }
-        if (!m->bMapStatic)
+        if (!bWasFull)
         {
             for (bj = 0; bj < m->b.nr; ++bj)
             {
@@ -1238,11 +1308,12 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
             {
                 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;
@@ -1250,7 +1321,6 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
 
     if (bMaskOnly)
     {
-        m->nr = m->b.nr;
         for (i = j = bj = 0; i < g->isize; ++i, ++j)
         {
             /* Find the next atom in the block */
@@ -1277,6 +1347,7 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
     }
     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 */
@@ -1300,10 +1371,8 @@ gmx_ana_indexmap_update(gmx_ana_indexmap_t *m, gmx_ana_index_t *g,
         }
         /* 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;
 }
 
@@ -1326,6 +1395,10 @@ gmx_ana_indexmap_deinit(gmx_ana_indexmap_t *m)
     {
         sfree(m->mapb.index);
     }
+    if (m->mapb.nalloc_a > 0)
+    {
+        sfree(m->mapb.a);
+    }
     sfree(m->orgid);
     if (m->b.nalloc_index > 0)
     {