src/gromacs/selection/indexutil.h

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2009-2018, The GROMACS development team.
   5  * Copyright (c) 2019, by the GROMACS development team, led by
   6  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   7  * and including many others, as listed in the AUTHORS file in the
   8  * top-level source directory and at http://www.gromacs.org.
   9  *
  10  * GROMACS is free software; you can redistribute it and/or
  11  * modify it under the terms of the GNU Lesser General Public License
  12  * as published by the Free Software Foundation; either version 2.1
  13  * of the License, or (at your option) any later version.
  14  *
  15  * GROMACS is distributed in the hope that it will be useful,
  16  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  17  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  18  * Lesser General Public License for more details.
  19  *
  20  * You should have received a copy of the GNU Lesser General Public
  21  * License along with GROMACS; if not, see
  22  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  23  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  24  *
  25  * If you want to redistribute modifications to GROMACS, please
  26  * consider that scientific software is very special. Version
  27  * control is crucial - bugs must be traceable. We will be happy to
  28  * consider code for inclusion in the official distribution, but
  29  * derived work must not be called official GROMACS. Details are found
  30  * in the README & COPYING files - if they are missing, get the
  31  * official version at http://www.gromacs.org.
  32  *
  33  * To help us fund GROMACS development, we humbly ask that you cite
  34  * the research papers on the package. Check out http://www.gromacs.org.
  35  */
  36 /*! \file
  37  * \brief API for handling index files and index groups.
  38  *
  39  * The API contains functions and data structures for handling index
  40  * files more conveniently than as several separate variables.
  41  * In addition to basic functions for initializing the data structures and
  42  * making copies, functions are provided for performing (most) set operations
  43  * on sorted index groups.
  44  * There is also a function for partitioning a index group based on
  45  * topology information such as residues or molecules.
  46  * Finally, there is a set of functions for constructing mappings between
  47  * an index group and its subgroups such.
  48  * These can be used with dynamic index group in calculations if one
  49  * needs to have a unique ID for each possible atom/residue/molecule in the
  50  * selection, e.g., for analysis of dynamics or for look-up tables.
  51  *
  52  * Mostly, these functions are used internally by the selection engine, but
  53  * it is necessary to use some of these functions in order to provide external
  54  * index groups to a gmx::SelectionCollection.
  55  * Some of the checking functions can be useful outside the selection engine to
  56  * check the validity of input groups.
  57  *
  58  * \author Teemu Murtola <teemu.murtola@gmail.com>
  59  * \inlibraryapi
  60  * \ingroup module_selection
  61  */
  62 #ifndef GMX_SELECTION_INDEXUTIL_H
  63 #define GMX_SELECTION_INDEXUTIL_H
  64
  65 #include <cstdio>
  66
  67 #include <string>
  68 #include <vector>
  69
  70 #include "gromacs/topology/block.h"
  71 #include "gromacs/utility/arrayref.h"
  72
  73 namespace gmx
  74 {
  75 class TextWriter;
  76
  77 /*!\internal \brief Bundle index groups with their names.
  78  *
  79  * \note This class has to outlive the t_blocka &indexGroup it is given.
  80  *       We want to avoid a deep-copy due to the potentially very large data
  81  *       stored in the t_blocka.
  82  *
  83  * \todo update this class, once the two legacy data structures, t_blocka and
  84  *       the char ** of group names are refactored.
  85  */
  86 class IndexGroupsAndNames
  87 {
  88 public:
  89     /*!\brief Construct from index group and group names
  90      * \param[in] indexGroup
  91      * \param[in] groupNames names of the index groups
  92      */
  93     IndexGroupsAndNames(const t_blocka& indexGroup, ArrayRef<char const* const> groupNames);
  94
  95     /*!\brief Return if a group name is contained in the groups.
  96      *
  97      * String comparison is case insensitive
  98      *
  99      * \param[in] groupName the group name to be queried
 100      * \returns true if index group name is contained
 101      */
 102     bool containsGroupName(const std::string& groupName) const;
 103
 104     /*!\brief Return the integer indices of a group.
 105      *
 106      * If two index groups share a name, return the one found first.
 107      *
 108      * Indices may be empty.
 109      *
 110      * \param[in] groupName the name of the group whose indices shall be returned
 111      * \returns atom indices of the selected index group
 112      * \throws if groupName is not present as index group
 113      */
 114     std::vector<index> indices(const std::string& groupName) const;
 115
 116 private:
 117     const t_blocka&          indexGroup_;
 118     std::vector<std::string> groupNames_;
 119 };
 120
 121 } // namespace gmx
 122
 123 struct gmx_mtop_t;
 124
 125 /** Stores a set of index groups. */
 126 struct gmx_ana_indexgrps_t;
 127
 128 /*! \brief
 129  * Specifies the type of index partition or index mapping in several contexts.
 130  *
 131  * \see gmx_ana_index_make_block(), gmx_ana_indexmap_init()
 132  */
 133 typedef enum
 134 {
 135     INDEX_UNKNOWN, /**< Unknown index type.*/
 136     INDEX_ATOM,    /**< Each atom in a separate block.*/
 137     INDEX_RES,     /**< Each residue in a separate block.*/
 138     INDEX_MOL,     /**< Each molecule in a separate block.*/
 139     INDEX_ALL      /**< All atoms in a single block.*/
 140 } e_index_t;
 141
 142 /*! \brief
 143  * Stores a single index group.
 144  */
 145 struct gmx_ana_index_t
 146 {
 147     /** Number of atoms. */
 148     int isize;
 149     /** List of atoms. */
 150     int* index;
 151     /** Number of items allocated for \p index. */
 152     int nalloc_index;
 153 };
 154
 155 /*! \brief
 156  * Data structure for calculating index group mappings.
 157  */
 158 struct gmx_ana_indexmap_t
 159 {
 160     /** Type of the mapping. */
 161     e_index_t type;
 162     /*! \brief
 163      * Current reference IDs.
 164      *
 165      * This array provides a mapping from the current index group (last given
 166      * to gmx_ana_indexmap_update()) to the blocks in \p b, i.e., the
 167      * original index group used in gmx_ana_indexmap_init().
 168      * The mapping is zero-based.
 169      * If \p bMaskOnly is provided to gmx_ana_indexmap_update(), the indices
 170      * for blocks not present in the current group are set to -1, otherwise
 171      * they are removed completely and the \p nr field updated.
 172      */
 173     int* refid;
 174     /*! \brief
 175      * Current mapped IDs.
 176      *
 177      * This array provides IDs for the current index group.  Instead of a
 178      * zero-based mapping that \p refid provides, the values from the \p orgid
 179      * array are used, thus allowing the mapping to be customized.
 180      * In other words, `mapid[i] = orgid[refid[i]]`.
 181      * If \p bMaskOnly is provided to gmx_ana_indexmap_update(), this array
 182      * equals \p orgid.
 183      */
 184     int* mapid;
 185     /*! \brief
 186      * Mapped block structure.
 187      *
 188      * A block structure that corresponds to the current index group.
 189      * \c mapb.nra and \c mapb.a correspond to the last mapped index group.
 190      */
 191     t_blocka mapb;
 192
 193     /*! \brief
 194      * Customizable ID numbers for the original blocks.
 195      *
 196      * This array has \p b.nr elements, each defining an original ID number for
 197      * a block in \p b (i.e., in the original group passed to
 198      * gmx_ana_indexmap_init()).
 199      * These are initialized in gmx_ana_indexmap_init() based on the type:
 200      *  - \ref INDEX_ATOM : the atom indices
 201      *  - \ref INDEX_RES :  the residue indices
 202      *  - \ref INDEX_MOL :  the molecule indices
 203      *
 204      * All the above numbers are zero-based.
 205      * After gmx_ana_indexmap_init(), the caller is free to change these values
 206      * if the above are not appropriate.
 207      * The mapped values can be read through \p mapid.
 208      */
 209     int* orgid;
 210
 211     /*! \brief
 212      * Block data that defines the mapping (internal use only).
 213      *
 214      * The data is initialized by gmx_ana_indexmap_init() and is not changed
 215      * after that.
 216      * Hence, it cannot be directly applied to the index group passed to
 217      * gmx_ana_indexmap_update() unless \p bMaskOnly was specified or the
 218      * index group is identical to the one provided to gmx_ana_indexmap_init().
 219      */
 220     t_blocka b;
 221     /*! \brief
 222      * true if the current reference IDs are for the whole group (internal use only).
 223      *
 224      * This is used internally to optimize the evaluation such that
 225      * gmx_ana_indexmap_update() does not take any time if the group is
 226      * actually static.
 227      */
 228     bool bStatic;
 229 };
 230
 231
 232 /*! \name Functions for handling gmx_ana_indexgrps_t
 233  */
 234 /*@{*/
 235 /** Reads index groups from a file or constructs them from topology. */
 236 void gmx_ana_indexgrps_init(gmx_ana_indexgrps_t** g, gmx_mtop_t* top, const char* fnm);
 237 /** Frees memory allocated for index groups. */
 238 void gmx_ana_indexgrps_free(gmx_ana_indexgrps_t* g);
 239 /** Returns true if the index group structure is emtpy. */
 240 bool gmx_ana_indexgrps_is_empty(gmx_ana_indexgrps_t* g);
 241
 242 /** Returns a pointer to an index group. */
 243 gmx_ana_index_t* gmx_ana_indexgrps_get_grp(gmx_ana_indexgrps_t* g, int n);
 244 /** Extracts a single index group. */
 245 bool gmx_ana_indexgrps_extract(gmx_ana_index_t* dest, std::string* destName, gmx_ana_indexgrps_t* src, int n);
 246 /** Finds and extracts a single index group by name. */
 247 bool gmx_ana_indexgrps_find(gmx_ana_index_t*     dest,
 248                             std::string*         destName,
 249                             gmx_ana_indexgrps_t* src,
 250                             const char*          name);
 251
 252 /** Writes out a list of index groups. */
 253 void gmx_ana_indexgrps_print(gmx::TextWriter* writer, gmx_ana_indexgrps_t* g, int maxn);
 254 /*@}*/
 255
 256 /*! \name Functions for handling gmx_ana_index_t
 257  */
 258 /*@{*/
 259 /** Reserves memory to store an index group of size \p isize. */
 260 void gmx_ana_index_reserve(gmx_ana_index_t* g, int isize);
 261 /** Frees any memory not necessary to hold the current contents. */
 262 void gmx_ana_index_squeeze(gmx_ana_index_t* g);
 263 /** Initializes an empty index group. */
 264 void gmx_ana_index_clear(gmx_ana_index_t* g);
 265 /** Constructs a \c gmx_ana_index_t from given values. */
 266 void gmx_ana_index_set(gmx_ana_index_t* g, int isize, int* index, int nalloc);
 267 /** Creates a simple index group from the first to the \p natoms'th atom. */
 268 void gmx_ana_index_init_simple(gmx_ana_index_t* g, int natoms);
 269 /** Frees memory allocated for an index group. */
 270 void gmx_ana_index_deinit(gmx_ana_index_t* g);
 271 /** Copies a \c gmx_ana_index_t. */
 272 void gmx_ana_index_copy(gmx_ana_index_t* dest, gmx_ana_index_t* src, bool bAlloc);
 273
 274 /** Writes out the contents of a index group. */
 275 void gmx_ana_index_dump(gmx::TextWriter* writer, gmx_ana_index_t* g, int maxn);
 276
 277 /*! \brief
 278  * Returns maximum atom index that appears in an index group.
 279  *
 280  * \param[in]  g      Index group to query.
 281  * \returns    Largest atom index that appears in \p g, or zero if \p g is empty.
 282  */
 283 int gmx_ana_index_get_max_index(gmx_ana_index_t* g);
 284 /** Checks whether an index group is sorted. */
 285 bool gmx_ana_index_check_sorted(gmx_ana_index_t* g);
 286 /*! \brief
 287  * Checks whether an index group has atoms from a defined range.
 288  *
 289  * \param[in]  g      Index group to check.
 290  * \param[in]  natoms Largest atom number allowed.
 291  * \returns    true if all atoms in the index group are in the
 292  *     range 0 to \p natoms (i.e., no atoms over \p natoms are referenced).
 293  */
 294 bool gmx_ana_index_check_range(gmx_ana_index_t* g, int natoms);
 295 /*@}*/
 296
 297 /*! \name Functions for set operations on gmx_ana_index_t
 298  */
 299 /*@{*/
 300 /** Sorts the indices within an index group. */
 301 void gmx_ana_index_sort(gmx_ana_index_t* g);
 302 /*! \brief
 303  * Removes duplicates from a sorted index group.
 304  *
 305  * \param[in,out] g  Index group to be processed.
 306  */
 307 void gmx_ana_index_remove_duplicates(gmx_ana_index_t* g);
 308 /** Checks whether two index groups are equal. */
 309 bool gmx_ana_index_equals(gmx_ana_index_t* a, gmx_ana_index_t* b);
 310 /** Checks whether a sorted index group contains another sorted index group. */
 311 bool gmx_ana_index_contains(gmx_ana_index_t* a, gmx_ana_index_t* b);
 312
 313 /** Calculates the intersection between two sorted index groups. */
 314 void gmx_ana_index_intersection(gmx_ana_index_t* dest, gmx_ana_index_t* a, gmx_ana_index_t* b);
 315 /** Calculates the set difference between two sorted index groups. */
 316 void gmx_ana_index_difference(gmx_ana_index_t* dest, gmx_ana_index_t* a, gmx_ana_index_t* b);
 317 /** Calculates the size of the difference between two sorted index groups. */
 318 int gmx_ana_index_difference_size(gmx_ana_index_t* a, gmx_ana_index_t* b);
 319 /** Calculates the union of two sorted index groups. */
 320 void gmx_ana_index_union(gmx_ana_index_t* dest, gmx_ana_index_t* a, gmx_ana_index_t* b);
 321 /*! \brief
 322  * Calculates the union of two index groups, where the second group may not be sorted.
 323  *
 324  * \param[out] dest Output index group (the union of \p a and \p b).
 325  * \param[in]  a    First index group (must be sorted).
 326  * \param[in]  b    Second index group.
 327  *
 328  * \p a and \p b can have common items.
 329  * \p dest can equal \p a or \p b.
 330  */
 331 void gmx_ana_index_union_unsorted(gmx_ana_index_t* dest, gmx_ana_index_t* a, gmx_ana_index_t* b);
 332 /** Merges two distinct sorted index groups. */
 333 void gmx_ana_index_merge(gmx_ana_index_t* dest, gmx_ana_index_t* a, gmx_ana_index_t* b);
 334 /** Calculates the intersection and the difference in one call. */
 335 void gmx_ana_index_partition(gmx_ana_index_t* dest1,
 336                              gmx_ana_index_t* dest2,
 337                              gmx_ana_index_t* src,
 338                              gmx_ana_index_t* g);
 339 /*@}*/
 340
 341 /*! \name Functions for handling gmx_ana_indexmap_t and related things
 342  */
 343 /*@{*/
 344 /** Partition a group based on topology information. */
 345 void gmx_ana_index_make_block(t_blocka* t, const gmx_mtop_t* top, gmx_ana_index_t* g, e_index_t type, bool bComplete);
 346 /** Checks whether a group consists of full blocks. */
 347 bool gmx_ana_index_has_full_blocks(const gmx_ana_index_t* g, const gmx::RangePartitioning* b);
 348 /** Checks whether a group consists of full blocks. */
 349 bool gmx_ana_index_has_full_ablocks(gmx_ana_index_t* g, t_blocka* b);
 350 /** Checks whether a group consists of full residues/molecules. */
 351 bool gmx_ana_index_has_complete_elems(gmx_ana_index_t* g, e_index_t type, const gmx_mtop_t* top);
 352
 353 /** Initializes an empty index group mapping. */
 354 void gmx_ana_indexmap_clear(gmx_ana_indexmap_t* m);
 355 /** Reserves memory for an index group mapping. */
 356 void gmx_ana_indexmap_reserve(gmx_ana_indexmap_t* m, int nr, int isize);
 357 /** Initializes an index group mapping. */
 358 void gmx_ana_indexmap_init(gmx_ana_indexmap_t* m, gmx_ana_index_t* g, const gmx_mtop_t* top, e_index_t type);
 359 /*! \brief
 360  * Initializes `orgid` entries based on topology grouping.
 361  *
 362  * \param[in,out] m     Mapping structure to use/initialize.
 363  * \param[in]     top   Topology structure
 364  *     (can be NULL if not required for \p type).
 365  * \param[in]     type  Type of groups to use.
 366  * \returns  The number of groups of type \p type that were present in \p m.
 367  * \throws   InconsistentInputError if the blocks in \p m do not have a unique
 368  *     group (e.g., contain atoms from multiple residues with `type == INDEX_RES`).
 369  *
 370  * By default, the gmx_ana_indexmap_t::orgid fields are initialized to
 371  * atom/residue/molecule indices from the topology (see documentation for the
 372  * struct for more details).
 373  * This function can be used to set the field to a zero-based group index
 374  * instead.  The first block will always get `orgid[0] = 0`, and all following
 375  * blocks that belong to the same residue/molecule (\p type) will get the same
 376  * index.  Each time a block does not belong to the same group, it will get the
 377  * next available number.
 378  * If `type == INDEX_ATOM`, the `orgid` field is initialized as 0, 1, ...,
 379  * independent of whether the blocks are single atoms or not.
 380  *
 381  * Strong exception safety guarantee.
 382  */
 383 int gmx_ana_indexmap_init_orgid_group(gmx_ana_indexmap_t* m, const gmx_mtop_t* top, e_index_t type);
 384 /** Sets an index group mapping to be static. */
 385 void gmx_ana_indexmap_set_static(gmx_ana_indexmap_t* m, t_blocka* b);
 386 /** Frees memory allocated for index group mapping. */
 387 void gmx_ana_indexmap_deinit(gmx_ana_indexmap_t* m);
 388 /** Makes a deep copy of an index group mapping. */
 389 void gmx_ana_indexmap_copy(gmx_ana_indexmap_t* dest, gmx_ana_indexmap_t* src, bool bFirst);
 390 /** Updates an index group mapping. */
 391 void gmx_ana_indexmap_update(gmx_ana_indexmap_t* m, gmx_ana_index_t* g, bool bMaskOnly);
 392 /*@}*/
 393
 394 #endif