colvargrid.h

// -*- c++ -*-
 
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/Colvars/colvars
// Please update all Colvars source files before making any changes.
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
 
#ifndef COLVARGRID_H
#define COLVARGRID_H
 
#include <iosfwd>
#include <memory>
 
#include "colvar.h"
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
 
 
class colvar_grid_params  {
 
public:
  size_t nd = 0;
 
  std::vector<int> nx;
 
  std::vector<int> nxc;
 
  std::vector<colvarvalue>  lower_boundaries;
 
  std::vector<colvarvalue>  upper_boundaries;
 
  std::vector<cvm::real>    widths;
};
 
 
template <class T> class colvar_grid : public colvar_grid_params, public colvarparse {
 
  //protected:
public: // TODO create accessors for these after all instantiations work
 
  size_t mult = 1;
 
  size_t nt;
 
  std::vector<T> data;
 
  std::vector<size_t> new_data;
 
  std::vector<colvar *> cv;
 
  std::vector<bool> use_actual_value;
 
  inline size_t address(std::vector<int> const &ix) const
  {
    size_t addr = 0;
    for (size_t i = 0; i < nd; i++) {
      addr += ix[i]*static_cast<size_t>(nxc[i]);
      if (cvm::debug()) {
        if (ix[i] >= nx[i] || ix[i] < 0) {
          cvmodule->error("Error: index along dimension " + cvm::to_str(i) + " (" + cvm::to_str(ix[i])
            + ") exceeds bounds in colvar_grid.\n", COLVARS_BUG_ERROR);
          return 0;
        }
      }
    }
    return addr;
  }
 
  inline void index(size_t address, std::vector<int> &ix) {
    for (int dim = nd - 1; dim >= 0; --dim) {
      ix[dim] = address % nx[dim];
      address /= nx[dim];
    }
  }
 
public:
  std::vector<bool>        periodic;
 
  std::vector<bool>        hard_lower_boundaries;
 
  std::vector<bool>        hard_upper_boundaries;
 
  bool has_parent_data;
 
  bool has_data;
 
  inline size_t num_variables() const
  {
    return nd;
  }
 
  inline std::vector<int> const &number_of_points_vec() const
  {
    return nx;
  }
 
  inline size_t number_of_points(int const icv = -1) const
  {
    if (icv < 0) {
      return nt;
    } else {
      return nx[icv];
    }
  }
 
  inline std::vector<int> const & sizes() const
  {
    return nx;
  }
 
  inline void set_sizes(std::vector<int> const &new_sizes)
  {
    nx = new_sizes;
  }
 
  inline size_t multiplicity() const
  {
    return mult;
  }
 
  inline void request_actual_value(bool b = true)
  {
    size_t i;
    for (i = 0; i < use_actual_value.size(); i++) {
      use_actual_value[i] = b;
    }
  }
 
  int setup(std::vector<int> const &nx_i,
            T const &t = T(),
            size_t const &mult_i = 1)
  {
    if (cvm::debug()) {
      cvmodule->log("Allocating grid: multiplicity = "+cvm::to_str(mult_i)+
               ", dimensions = "+cvm::to_str(nx_i)+".\n");
    }
 
    mult = mult_i;
 
    data.clear();
 
    nx = nx_i;
    nd = nx.size();
 
    nxc.resize(nd);
 
    // setup dimensions
    nt = mult;
    for (int i = nd-1; i >= 0; i--) {
      if (nx[i] <= 0) {
        cvm::error_static("Error: providing an invalid number of grid points, "+
                   cvm::to_str(nx[i])+".\n", COLVARS_BUG_ERROR);
        return COLVARS_ERROR;
      }
      nxc[i] = nt;
      nt *= nx[i];
    }
 
    if (cvm::debug()) {
      cvmodule->log("Total number of grid elements = "+cvm::to_str(nt)+".\n");
    }
 
    data.reserve(nt);
    data.assign(nt, t);
 
    return COLVARS_OK;
  }
 
  int setup()
  {
    return setup(this->nx, T(), this->mult);
  }
 
  void reset(T const &t = T())
  {
    data.assign(nt, t);
  }
 
 
  // This constructor depends on a static cvm pointer and is deprecated
  colvar_grid() : colvarparse(cvm::main()), has_data(false)
  {
    nd = nt = 0;
    mult = 1;
    has_parent_data = false;
    this->setup();
  }
 
  virtual ~colvar_grid()
  {}
 
  // This constructor depends on a static cvm pointer and is deprecated
  colvar_grid(colvar_grid<T> const &g) : colvar_grid_params(colvar_grid_params(g)),
                                         colvarparse(cvm::main()),
                                         mult(g.mult),
                                         data(),
                                         cv(g.cv),
                                         use_actual_value(g.use_actual_value),
                                         periodic(g.periodic),
                                         hard_lower_boundaries(g.hard_lower_boundaries),
                                         hard_upper_boundaries(g.hard_upper_boundaries),
                                         has_parent_data(false),
                                         has_data(false)
  {}
 
  colvar_grid(std::vector<int> const &nx_i,
              T const &t = T(),
              size_t mult_i = 1)
    : colvarparse(cvm::main()), has_parent_data(false), has_data(false)
  {
    this->setup(nx_i, t, mult_i);
  }
 
  colvar_grid(std::vector<colvar *> const &colvars,
              T const &t = T(),
              size_t mult_i = 1,
              bool add_extra_bin = false,
              std::shared_ptr<const colvar_grid_params> params = nullptr,
              std::string config = std::string())
    : colvarparse(cvm::main()), has_parent_data(false), has_data(false)
  {
    (void) t;
    this->init_from_colvars(colvars, mult_i, add_extra_bin, params, config);
  }
 
  colvar_grid(std::string const &filename, size_t mult_i = 1);
 
  int init_from_colvars(std::vector<colvar *> const &colvars,
                        size_t mult_i = 1,
                        bool add_extra_bin = false,
                        std::shared_ptr<const colvar_grid_params> params = nullptr,
                        std::string config = std::string())
  {
    if (cvm::debug()) {
      cvmodule->log("Reading grid configuration from collective variables.\n");
    }
 
    cv = colvars;
    nd = colvars.size();
    mult = mult_i;
 
    size_t i;
 
    if (cvm::debug()) {
      cvmodule->log("Allocating a grid for "+cvm::to_str(colvars.size())+
               " collective variables, multiplicity = "+cvm::to_str(mult_i)+".\n");
    }
 
    for (i =  0; i < nd; i++) {
      if (cv[i]->value().type() != colvarvalue::type_scalar) {
        cvm::error_static("Colvar grids can only be automatically "
                   "constructed for scalar variables.  "
                   "ABF and histogram can not be used; metadynamics "
                   "can be used with useGrids disabled.\n", COLVARS_INPUT_ERROR);
        return COLVARS_ERROR;
      }
 
      if (cv[i]->width <= 0.0) {
        cvm::error_static("Tried to initialize a grid on a "
                   "variable with negative or zero width.\n", COLVARS_INPUT_ERROR);
        return COLVARS_ERROR;
      }
 
      widths.push_back(cv[i]->width);
      hard_lower_boundaries.push_back(cv[i]->is_enabled(colvardeps::f_cv_hard_lower_boundary));
      hard_upper_boundaries.push_back(cv[i]->is_enabled(colvardeps::f_cv_hard_upper_boundary));
 
      // By default, get reported colvar value (for extended Lagrangian colvars)
      use_actual_value.push_back(false);
 
      // except if a colvar is specified twice in a row
      // then the first instance is the actual value
      // For histograms of extended-system coordinates
      if (i > 0 && cv[i-1] == cv[i]) {
        use_actual_value[i-1] = true;
      }
 
      // This needs to work if the boundaries are undefined in the colvars
      lower_boundaries.push_back(cv[i]->lower_boundary);
      upper_boundaries.push_back(cv[i]->upper_boundary);
    }
 
    // Replace widths and boundaries with optional custom configuration
    if (!config.empty()) {
      this->parse_params(config);
      this->check_keywords(config, "grid");
 
      if (params) {
        cvm::error_static("Error: init_from_colvars was passed both a grid config and a template grid.", COLVARS_BUG_ERROR);
        return COLVARS_BUG_ERROR;
      }
    } else if (params) {
      // Match grid sizes with template
 
      if (params->nd != nd) {
        cvm::error_static("Trying to initialize grid from template with wrong dimension (" +
                    cvm::to_str(params->nd) + " instead of " +
                    cvm::to_str(this->nd) + ").");
        return COLVARS_ERROR;
      }
 
      widths =params->widths;
      lower_boundaries =params->lower_boundaries;
      upper_boundaries =params->upper_boundaries;
    }
 
    // Only now can we determine periodicity
    for (i =  0; i < nd; i++) {
      periodic.push_back(cv[i]->periodic_boundaries(lower_boundaries[i].real_value,
                                                    upper_boundaries[i].real_value));
 
      if (add_extra_bin) {
        // Shift the grid by half the bin width (values at edges instead of center of bins)
        lower_boundaries[i] -= 0.5 * widths[i];
 
        if (periodic[i]) {
          // Just shift
          upper_boundaries[i] -= 0.5 * widths[i];
        } else {
          // Widen grid by one bin width
          upper_boundaries[i] += 0.5 * widths[i];
        }
      }
    }
 
    // Reset grid sizes based on widths and boundaries
    this->init_from_boundaries();
    return this->setup();
  }
 
  int init_from_boundaries()
  {
    if (cvm::debug()) {
      cvmodule->log("Configuring grid dimensions from colvars boundaries.\n");
    }
 
    // these will have to be updated
    nx.clear();
    nxc.clear();
    nt = 0;
 
    for (size_t i =  0; i < lower_boundaries.size(); i++) {
      // Re-compute periodicity using current grid boundaries
      periodic[i] = cv[i]->periodic_boundaries(lower_boundaries[i].real_value,
                                               upper_boundaries[i].real_value);
 
      cvm::real nbins = ( upper_boundaries[i].real_value -
                          lower_boundaries[i].real_value ) / widths[i];
      int nbins_round = (int)(nbins+0.5);
 
      if (cvm::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) {
        cvmodule->log("Warning: grid interval("+
                 cvm::to_str(lower_boundaries[i], cvmodule->cv_width, cvmodule->cv_prec)+" - "+
                 cvm::to_str(upper_boundaries[i], cvmodule->cv_width, cvmodule->cv_prec)+
                 ") is not commensurate to its bin width("+
                 cvm::to_str(widths[i], cvmodule->cv_width, cvmodule->cv_prec)+").\n");
        upper_boundaries[i].real_value = lower_boundaries[i].real_value +
          (nbins_round * widths[i]);
      }
 
      if (cvm::debug())
        cvmodule->log("Number of points is "+cvm::to_str((int) nbins_round)+
                 " for the colvar no. "+cvm::to_str(i+1)+".\n");
 
      nx.push_back(nbins_round);
    }
 
    return COLVARS_OK;
  }
 
  inline void wrap(std::vector<int> & ix) const
  {
    for (size_t i = 0; i < nd; i++) {
      if (periodic[i]) {
        ix[i] = (ix[i] + nx[i]) % nx[i]; // Avoid modulo with negative operands (implementation-defined)
      } else {
        if (ix[i] < 0 || ix[i] >= nx[i]) {
          cvm::error_static("Trying to wrap illegal index vector (non-PBC) for a grid point: "
                     + cvm::to_str(ix), COLVARS_BUG_ERROR);
          return;
        }
      }
    }
  }
 
  inline bool wrap_detect_edge(std::vector<int> & ix) const
  {
    bool edge = false;
    for (size_t i = 0; i < nd; i++) {
      if (periodic[i]) {
        ix[i] = (ix[i] + nx[i]) % nx[i]; // Avoid modulo with negative operands (implementation-defined)
      } else if (ix[i] < 0 || ix[i] >= nx[i]) {
        edge = true;
      }
    }
    return edge;
  }
 
  inline bool wrap_to_edge(std::vector<int> & ix, std::vector<int> & edge_bin) const
  {
    bool edge = false;
    edge_bin = ix;
    for (size_t i = 0; i < nd; i++) {
      if (periodic[i]) {
        ix[i] = (ix[i] + nx[i]) % nx[i]; // Avoid modulo with negative operands (implementation-defined)
        edge_bin[i] = ix[i];
      } else if (ix[i] < 0) {
        edge = true;
        edge_bin[i] = 0;
      } else if (ix[i] >= nx[i]) {
        edge = true;
        edge_bin[i] = nx[i] - 1;
      }
    }
    return edge;
  }
 
 
  inline int current_bin_scalar(int const i) const
  {
    return value_to_bin_scalar(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
  }
 
  inline int current_bin_flat_bound() const
  {
    std::vector<int> index = new_index();
    for (size_t i = 0; i < nd; i++) {
      index[i] = current_bin_scalar_bound(i);
    }
    return address(index);
  }
 
  inline int current_bin_scalar_bound(int const i) const
  {
    return value_to_bin_scalar_bound(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
  }
 
  inline int current_bin_scalar(int const i, int const iv) const
  {
    return value_to_bin_scalar(use_actual_value[i] ?
                               cv[i]->actual_value().vector1d_value[iv] :
                               cv[i]->value().vector1d_value[iv], i);
  }
 
  inline int value_to_bin_scalar(colvarvalue const &value, const int i) const
  {
    return (int) cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
  }
 
  inline cvm::real current_bin_scalar_fraction(int const i) const
  {
    return value_to_bin_scalar_fraction(use_actual_value[i] ? cv[i]->actual_value() : cv[i]->value(), i);
  }
 
  inline cvm::real value_to_bin_scalar_fraction(colvarvalue const &value, const int i) const
  {
    cvm::real x = (value.real_value - lower_boundaries[i].real_value) / widths[i];
    return x - cvm::floor(x);
  }
 
  inline int value_to_bin_scalar_bound(colvarvalue const &value, const int i) const
  {
    int bin_index = cvm::floor( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
 
    // Wrap bins for periodic dimensions before truncating
    if (periodic[i]) bin_index %= nx[i];
    if (bin_index < 0) bin_index=0;
    if (bin_index >=int(nx[i])) bin_index=int(nx[i])-1;
    return (int) bin_index;
  }
 
  inline int value_to_bin_scalar(colvarvalue const &value,
                                 colvarvalue const &new_offset,
                                 cvm::real   const &new_width) const
  {
    return (int) cvm::floor( (value.real_value - new_offset.real_value) / new_width );
  }
 
  inline colvarvalue bin_to_value_scalar(int const &i_bin, int const i) const
  {
    return lower_boundaries[i].real_value + widths[i] * (0.5 + i_bin);
  }
 
  inline colvarvalue bin_to_value_scalar(int const &i_bin,
                                         colvarvalue const &new_offset,
                                         cvm::real const &new_width) const
  {
    return new_offset.real_value + new_width * (0.5 + i_bin);
  }
 
  inline void set_value(std::vector<int> const &ix,
                        T const &t,
                        size_t const &imult = 0)
  {
    data[this->address(ix)+imult] = t;
    has_data = true;
  }
 
  inline void set_value(size_t i, T const &t)
  {
    data[i] = t;
  }
 
  inline T get_value(size_t i) const
  {
    return data[i];
  }
 
 
  void delta_grid(colvar_grid<T> const &other_grid)
  {
 
    if (other_grid.multiplicity() != this->multiplicity()) {
      cvm::error_static("Error: trying to subtract two grids with "
                 "different multiplicity.\n");
      return;
    }
 
    if (other_grid.data.size() != this->data.size()) {
      cvm::error_static("Error: trying to subtract two grids with "
                 "different size.\n");
      return;
    }
 
    for (size_t i = 0; i < data.size(); i++) {
      data[i] = other_grid.data[i] - data[i];
    }
    has_data = true;
  }
 
 
  void copy_grid(colvar_grid<T> const &other_grid)
  {
    if (other_grid.multiplicity() != this->multiplicity()) {
      cvm::error_static("Error: trying to copy two grids with "
                 "different multiplicity.\n");
      return;
    }
 
    if (other_grid.data.size() != this->data.size()) {
      cvm::error_static("Error: trying to copy two grids with "
                 "different size.\n");
      return;
    }
 
 
    for (size_t i = 0; i < data.size(); i++) {
      data[i] = other_grid.data[i];
    }
    has_data = true;
  }
 
  void raw_data_out(T* out_data) const
  {
    for (size_t i = 0; i < data.size(); i++) out_data[i] = data[i];
  }
  void raw_data_out(std::vector<T>& out_data) const
  {
    out_data = data;
  }
  void raw_data_in(const T* in_data)
  {
    for (size_t i = 0; i < data.size(); i++) data[i] = in_data[i];
    has_data = true;
  }
  void raw_data_in(const std::vector<T>& in_data)
  {
    data = in_data;
    has_data = true;
  }
  size_t raw_data_num() const { return data.size(); }
 
 
  inline T const & value(std::vector<int> const &ix,
                         size_t const &imult = 0) const
  {
    return data[this->address(ix) + imult];
  }
 
  inline T const & value(size_t i) const
  {
    return data[i];
  }
 
  inline void add_constant(T const &t)
  {
    for (size_t i = 0; i < nt; i++)
      data[i] += t;
    has_data = true;
  }
 
  inline void multiply_constant(cvm::real const &a)
  {
    for (size_t i = 0; i < nt; i++)
      data[i] *= a;
  }
 
  inline void remove_small_values(cvm::real const &a)
  {
    for (size_t i = 0; i < nt; i++)
      if(data[i]<a) data[i] = a;
  }
 
 
  inline std::vector<int> const get_colvars_index(std::vector<colvarvalue> const &values) const
  {
    std::vector<int> index = new_index();
    for (size_t i = 0; i < nd; i++) {
      index[i] = value_to_bin_scalar(values[i], i);
    }
    return index;
  }
 
  inline std::vector<int> const get_colvars_index() const
  {
    std::vector<int> index = new_index();
    for (size_t i = 0; i < nd; i++) {
      index[i] = current_bin_scalar(i);
    }
    return index;
  }
 
  inline std::vector<int> const get_colvars_index_bound() const
  {
    std::vector<int> index = new_index();
    for (size_t i = 0; i < nd; i++) {
      index[i] = current_bin_scalar_bound(i);
    }
    return index;
  }
 
  inline cvm::real bin_distance_from_boundaries(std::vector<colvarvalue> const &values,
                                                bool skip_hard_boundaries = false)
  {
    cvm::real minimum = 1.0E+16;
    for (size_t i = 0; i < nd; i++) {
 
      if (periodic[i]) continue;
 
      cvm::real dl = cvm::sqrt(cv[i]->dist2(values[i], lower_boundaries[i])) / widths[i];
      cvm::real du = cvm::sqrt(cv[i]->dist2(values[i], upper_boundaries[i])) / widths[i];
 
      if (values[i].real_value < lower_boundaries[i])
        dl *= -1.0;
      if (values[i].real_value > upper_boundaries[i])
        du *= -1.0;
 
      if ( ((!skip_hard_boundaries) || (!hard_lower_boundaries[i])) && (dl < minimum))
        minimum = dl;
      if ( ((!skip_hard_boundaries) || (!hard_upper_boundaries[i])) && (du < minimum))
        minimum = du;
    }
 
    return minimum;
  }
 
 
  void map_grid(colvar_grid<T> const &other_grid)
  {
    if (other_grid.multiplicity() != this->multiplicity()) {
      cvm::error_static("Error: trying to merge two grids with values of "
                 "different multiplicity.\n");
      return;
    }
 
    std::vector<colvarvalue> const &gb  = this->lower_boundaries;
    std::vector<cvm::real> const &gw    = this->widths;
    std::vector<colvarvalue> const &ogb = other_grid.lower_boundaries;
    std::vector<cvm::real> const &ogw   = other_grid.widths;
 
    std::vector<int> ix = this->new_index();
    std::vector<int> oix = other_grid.new_index();
 
    if (cvm::debug())
      cvmodule->log("Remapping grid...\n");
    for ( ; this->index_ok(ix); this->incr(ix)) {
 
      for (size_t i = 0; i < nd; i++) {
        oix[i] =
          value_to_bin_scalar(bin_to_value_scalar(ix[i], gb[i], gw[i]),
                              ogb[i],
                              ogw[i]);
      }
 
      if (! other_grid.index_ok(oix)) {
        continue;
      }
 
      for (size_t im = 0; im < mult; im++) {
        this->set_value(ix, other_grid.value(oix, im), im);
      }
    }
 
    has_data = true;
    if (cvm::debug())
      cvmodule->log("Remapping done.\n");
  }
 
  void add_grid(colvar_grid<T> const &other_grid,
                cvm::real scale_factor = 1.0)
  {
    if (other_grid.multiplicity() != this->multiplicity()) {
      cvm::error_static("Error: trying to sum togetehr two grids with values of "
                 "different multiplicity.\n");
      return;
    }
    if (scale_factor != 1.0)
      for (size_t i = 0; i < data.size(); i++) {
        data[i] += static_cast<T>(scale_factor * other_grid.data[i]);
      }
    else
      // skip multiplication if possible
      for (size_t i = 0; i < data.size(); i++) {
        data[i] += other_grid.data[i];
      }
    has_data = true;
  }
 
  virtual T value_output(std::vector<int> const &ix,
                         size_t const &imult = 0) const
  {
    return value(ix, imult);
  }
 
  virtual void value_input(std::vector<int> const &ix,
                           T const &t,
                           size_t const &imult = 0,
                           bool add = false)
  {
    if ( add )
      data[address(ix) + imult] += t;
    else
      data[address(ix) + imult] = t;
    has_data = true;
  }
 
 
  //   /// Get the pointer to the binned value indexed by ix
  //   inline T const *value_p (std::vector<int> const &ix)
  //   {
  //     return &(data[address (ix)]);
  //   }
 
  inline std::vector<int> const new_index() const
  {
    return std::vector<int> (nd, 0);
  }
 
  inline bool index_ok(std::vector<int> const &ix) const
  {
    for (size_t i = 0; i < nd; i++) {
      if ( (ix[i] < 0) || (ix[i] >= int(nx[i])) )
        return false;
    }
    return true;
  }
 
  inline void incr(std::vector<int> &ix) const
  {
    for (int i = ix.size()-1; i >= 0; i--) {
 
      ix[i]++;
 
      if (ix[i] >= nx[i]) {
 
        if (i > 0) {
          ix[i] = 0;
          continue;
        } else {
          // this is the last iteration, a non-valid index is being
          // set for the outer index, which will be caught by
          // index_ok()
          ix[0] = nx[0];
          return;
        }
      } else {
        return;
      }
    }
  }
 
  std::string get_state_params() const;
 
  int parse_params(std::string const &conf,
                   colvarparse::Parse_Mode const parse_mode = colvarparse::parse_normal);
 
  void check_consistency()
  {
    for (size_t i = 0; i < nd; i++) {
      if ( (cvm::sqrt(cv[i]->dist2(cv[i]->lower_boundary,
                                   lower_boundaries[i])) > 1.0E-10) ||
           (cvm::sqrt(cv[i]->dist2(cv[i]->upper_boundary,
                                   upper_boundaries[i])) > 1.0E-10) ||
           (cvm::sqrt(cv[i]->dist2(cv[i]->width,
                                   widths[i])) > 1.0E-10) ) {
        cvm::error_static("Error: restart information for a grid is "
                   "inconsistent with that of its colvars.\n");
        return;
      }
    }
  }
 
 
  void check_consistency(colvar_grid<T> const &other_grid)
  {
    for (size_t i = 0; i < nd; i++) {
      // we skip dist2(), because periodicities and the like should
      // matter: boundaries should be EXACTLY the same (otherwise,
      // map_grid() should be used)
      if ( (cvm::fabs(other_grid.lower_boundaries[i] -
                      lower_boundaries[i]) > 1.0E-10) ||
           (cvm::fabs(other_grid.upper_boundaries[i] -
                      upper_boundaries[i]) > 1.0E-10) ||
           (cvm::fabs(other_grid.widths[i] -
                      widths[i]) > 1.0E-10) ||
           (data.size() != other_grid.data.size()) ) {
        cvm::error_static("Error: inconsistency between "
                   "two grids that are supposed to be equal, "
                   "aside from the data stored.\n");
        return;
      }
    }
  }
 
  std::istream & read_restart(std::istream &is);
 
  cvm::memory_stream & read_restart(cvm::memory_stream &is);
 
  std::ostream & write_restart(std::ostream &os);
 
  cvm::memory_stream & write_restart(cvm::memory_stream &os);
 
  std::istream &read_raw(std::istream &is);
 
  cvm::memory_stream &read_raw(cvm::memory_stream &is);
 
private:
 
  template <class IST> IST& read_restart_template_(IST& is);
 
  template <class IST> IST& read_raw_template_(IST& is);
 
public:
 
  std::ostream &write_raw(std::ostream &os, size_t const buf_size = 3) const;
 
  cvm::memory_stream &write_raw(cvm::memory_stream &os, size_t const buf_size = 3) const;
 
  std::istream & read_multicol(std::istream &is, bool add = false);
 
  int read_multicol(std::string const &filename,
                    std::string description = "grid file",
                    bool add = false);
 
  std::ostream & write_multicol(std::ostream &os) const;
 
  int write_multicol(std::string const &filename,
                     std::string description = "grid file") const;
 
  std::ostream & write_opendx(std::ostream &os) const;
 
  int write_opendx(std::string const &filename,
                   std::string description = "grid file") const;
};
 
 
 
class colvar_grid_count : public colvar_grid<size_t>
{
public:
 
  colvar_grid_count();
 
  virtual ~colvar_grid_count()
  {}
 
  colvar_grid_count(std::vector<colvar *>  &colvars,
                    std::shared_ptr<const colvar_grid_params> params = nullptr);
 
  colvar_grid_count(std::vector<colvar *>  &colvars,
                    std::string            config);
 
  colvar_grid_count(std::string &filename);
 
  inline void incr_count(std::vector<int> const &ix)
  {
    ++(data[this->address(ix)]);
  }
 
  inline size_t const & new_value(std::vector<int> const &ix)
  {
    return new_data[address(ix)];
  }
 
  std::string get_state_params() const;
 
  int parse_params(std::string const &conf,
                   colvarparse::Parse_Mode const parse_mode = colvarparse::parse_normal);
 
  std::istream & read_restart(std::istream &is);
 
  cvm::memory_stream & read_restart(cvm::memory_stream &is);
 
  std::ostream & write_restart(std::ostream &os);
 
  cvm::memory_stream & write_restart(cvm::memory_stream &os);
 
  std::istream &read_raw(std::istream &is);
 
  cvm::memory_stream &read_raw(cvm::memory_stream &is);
 
  std::ostream &write_raw(std::ostream &os, size_t const buf_size = 3) const;
 
  cvm::memory_stream &write_raw(cvm::memory_stream &os, size_t const buf_size = 3) const;
 
  std::istream & read_multicol(std::istream &is, bool add = false);
 
  int read_multicol(std::string const &filename,
                            std::string description = "grid file",
                            bool add = false);
 
  std::ostream & write_multicol(std::ostream &os) const;
 
  int write_multicol(std::string const &filename,
                     std::string description = "grid file") const;
 
  std::ostream & write_opendx(std::ostream &os) const;
 
  int write_opendx(std::string const &filename,
                   std::string description = "grid file") const;
 
  virtual void value_input(std::vector<int> const &ix,
                           size_t const &t,
                           size_t const &imult = 0,
                           bool add = false)
  {
    (void) imult;
    if (add) {
      data[address(ix)] += t;
      if (this->has_parent_data) {
        // save newly read data for inputting parent grid
        new_data[address(ix)] = t;
      }
    } else {
      data[address(ix)] = t;
    }
    has_data = true;
  }
 
  inline int local_sample_count(int radius)
  {
    std::vector<int> ix0 = new_index();
    std::vector<int> ix = new_index();
 
    for (size_t i = 0; i < nd; i++) {
      ix0[i] = current_bin_scalar_bound(i);
    }
    if (radius < 1) {
      // Simple case: no averaging
      if (index_ok(ix0))
        return value(ix0);
      else
        return 0;
    }
    size_t count = 0;
    size_t nbins = 0;
    int i, j, k;
    bool edge;
    ix = ix0;
    // Treat each dimension separately to simplify code
    switch (nd)
    {
    case 1:
      for (i = -radius; i <= radius; i++) {
        ix[0] = ix0[0] + i;
        edge = wrap_detect_edge(ix);
        if (!edge) {
          nbins++;
          count += value(ix);
        }
      }
      break;
    case 2:
      for (i = -radius; i <= radius; i++) {
        ix[0] = ix0[0] + i;
        for (j = -radius; j <= radius; j++) {
          ix[1] = ix0[1] + j;
          edge = wrap_detect_edge(ix);
          if (!edge) {
            nbins++;
            count += value(ix);
          }
        }
      }
      break;
    case 3:
      for (i = -radius; i <= radius; i++) {
        ix[0] = ix0[0] + i;
        for (j = -radius; j <= radius; j++) {
          ix[1] = ix0[1] + j;
          for (k = -radius; k <= radius; k++) {
            ix[2] = ix0[2] + k;
            edge = wrap_detect_edge(ix);
            if (!edge) {
              nbins++;
              count += value(ix);
            }
          }
        }
      }
      break;
    default:
      cvm::error_static("Error: local_sample_count is not implemented for grids of dimension > 3", COLVARS_NOT_IMPLEMENTED);
      break;
    }
 
    if (nbins)
      // Integer division - an error on the order of 1 doesn't matter
      return count / nbins;
    else
      return 0.0;
  }
 
 
  inline cvm::real log_gradient_finite_diff(const std::vector<int> &ix0,
                                            int n = 0, int offset = 0)
  {
    cvm::real Am1, A0, A1, A2; // values at locations -1, 0, +1, +2
    std::vector<int> ix = ix0;
 
    // TODO this can be rewritten more concisely with wrap_edge()
    if (periodic[n]) {
      ix[n]--; wrap(ix);
      Am1 = value(ix) + offset;
      ix = ix0; wrap(ix);
      A0 = value(ix) + offset; // Used just to detect gaps
      ix[n]++; wrap(ix);
      A1 = value(ix) + offset;
      if (Am1 * A0 * A1 == 0) {
        return 0.; // can't handle 0, and don't compute derivative over a gap
      } else {
        return (cvmodule->logn(A1) - cvmodule->logn(Am1))
          / (widths[n] * 2.);
      }
    } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
      ix[n]--;
      Am1 = value(ix) + offset;
      ix = ix0;
      A0 = value(ix) + offset;
      ix[n]++;
      A1 = value(ix) + offset;
      if (Am1 * A0 * A1 == 0) {
        return 0.; // can't handle 0, and don't compute derivative over a gap
      } else {
        return (cvmodule->logn(A1) - cvmodule->logn(Am1))
          / (widths[n] * 2.);
      }
    } else {
      // edge: use 2nd order derivative
      int increment = (ix[n] == 0 ? 1 : -1);
      // move right from left edge, or the other way around
      A0 = value(ix) + offset;
      ix[n] += increment; A1 = value(ix) + offset;
      ix[n] += increment; A2 = value(ix) + offset;
      if (A0 * A1 * A2 == 0) {
        return 0.; // can't handle 0
      } else {
        return (-1.5 * cvmodule->logn(A0) + 2. * cvmodule->logn(A1)
          - 0.5 * cvmodule->logn(A2)) * increment / widths[n];
      }
    }
  }
 
 
  inline cvm::real gradient_finite_diff(const std::vector<int> &ix0,
                                        int n = 0)
  {
    cvm::real Am1, A0, A1, A2;
    std::vector<int> ix = ix0;
 
    // FIXME this can be rewritten more concisely with wrap_edge()
    if (periodic[n]) {
      ix[n]--; wrap(ix);
      Am1 = value(ix);
      ix = ix0; wrap(ix);
      A0 = value(ix); // Used just to detect gaps
      ix[n]++; wrap(ix);
      A1 = value(ix);
      if (Am1 * A0 * A1 == 0) {
        return 0.; // can't handle empty bins
      } else {
        return (A1 - Am1) / (widths[n] * 2.);
      }
    } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
      ix[n]--;
      Am1 = value(ix);
      ix = ix0;
      A0 = value(ix);
      ix[n]++;
      A1 = value(ix);
      if (Am1 * A0 * A1 == 0) {
        return 0.; // can't handle empty bins
      } else {
        return (A1 - Am1) / (widths[n] * 2.);
      }
    } else {
      // edge: use 2nd order derivative
      int increment = (ix[n] == 0 ? 1 : -1);
      // move right from left edge, or the other way around
      A0 = value(ix);
      ix[n] += increment; A1 = value(ix);
      ix[n] += increment; A2 = value(ix);
      return (-1.5 * A0 + 2. * A1
          - 0.5 * A2) * increment / widths[n];
    }
  }
};
 
 
class colvar_grid_scalar : public colvar_grid<cvm::real>
{
public:
 
  std::shared_ptr<colvar_grid_count> samples;
 
  colvar_grid_scalar();
 
  colvar_grid_scalar(colvar_grid_scalar const &g);
 
  virtual ~colvar_grid_scalar();
 
  colvar_grid_scalar(std::vector<colvar *> &colvars,
                     std::shared_ptr<const colvar_grid_params> params = nullptr,
                     bool add_extra_bin = false,
                     std::string config = std::string());
 
  colvar_grid_scalar(std::string const &filename);
 
  inline void acc_value(std::vector<int> const &ix,
                        cvm::real const &new_value,
                        size_t const &imult = 0)
  {
    (void) imult;
    // only legal value of imult here is 0
    data[address(ix)] += new_value;
    if (samples)
      samples->incr_count(ix);
    has_data = true;
  }
 
  std::string get_state_params() const;
 
  int parse_params(std::string const &conf,
                   colvarparse::Parse_Mode const parse_mode = colvarparse::parse_normal);
 
  std::istream & read_restart(std::istream &is);
 
  cvm::memory_stream & read_restart(cvm::memory_stream &is);
 
  std::ostream & write_restart(std::ostream &os);
 
  cvm::memory_stream & write_restart(cvm::memory_stream &os);
 
  std::istream &read_raw(std::istream &is);
 
  cvm::memory_stream &read_raw(cvm::memory_stream &is);
 
  std::ostream &write_raw(std::ostream &os, size_t const buf_size = 3) const;
 
  cvm::memory_stream &write_raw(cvm::memory_stream &os, size_t const buf_size = 3) const;
 
  std::istream & read_multicol(std::istream &is, bool add = false);
 
  int read_multicol(std::string const &filename,
                    std::string description = "grid file",
                    bool add = false);
 
  std::ostream & write_multicol(std::ostream &os) const;
 
  int write_multicol(std::string const &filename,
                     std::string description = "grid file") const;
 
  std::ostream & write_opendx(std::ostream &os) const;
 
  int write_opendx(std::string const &filename,
                   std::string description = "grid file") const;
 
  inline void vector_gradient_finite_diff( const std::vector<int> &ix0, std::vector<cvm::real> &grad)
  {
    cvm::real A0, A1;
    std::vector<int> ix;
    size_t i, j, k, n;
 
    if (nd == 2) {
      for (n = 0; n < 2; n++) {
        ix = ix0;
        A0 = value(ix);
        ix[n]++; wrap(ix);
        A1 = value(ix);
        ix[1-n]++; wrap(ix);
        A1 += value(ix);
        ix[n]--; wrap(ix);
        A0 += value(ix);
        grad[n] = 0.5 * (A1 - A0) / widths[n];
      }
    } else if (nd == 3) {
 
      cvm::real p[8]; // potential values within cube, indexed in binary (4 i + 2 j + k)
      ix = ix0;
      int index = 0;
      for (i = 0; i<2; i++) {
        ix[1] = ix0[1];
        for (j = 0; j<2; j++) {
          ix[2] = ix0[2];
          for (k = 0; k<2; k++) {
            wrap(ix);
            p[index++] = value(ix);
            ix[2]++;
          }
          ix[1]++;
        }
        ix[0]++;
      }
 
      // The following would be easier to read using binary literals
      //                  100    101    110    111      000    001    010   011
      grad[0] = 0.25 * ((p[4] + p[5] + p[6] + p[7]) - (p[0] + p[1] + p[2] + p[3])) / widths[0];
      //                  010     011    110   111      000    001    100   101
      grad[1] = 0.25 * ((p[2] + p[3] + p[6] + p[7]) - (p[0] + p[1] + p[4] + p[5])) / widths[1];
      //                  001    011     101   111      000    010   100    110
      grad[2] = 0.25 * ((p[1] + p[3] + p[5] + p[7]) - (p[0] + p[2] + p[4] + p[6])) / widths[2];
    } else {
      cvm::error_static("Finite differences available in dimension 2 and 3 only.");
    }
  }
 
 
  inline cvm::real log_gradient_finite_diff(const std::vector<int> &ix0,
                                            int n = 0, int offset = 0)
  {
    cvm::real A0, A1, A2;
    std::vector<int> ix = ix0;
 
    // TODO this can be rewritten more concisely with wrap_edge()
    if (periodic[n]) {
      ix[n]--; wrap(ix);
      A0 = value(ix) + offset;
      ix = ix0;
      ix[n]++; wrap(ix);
      A1 = value(ix) + offset;
      if (A0 * A1 == 0) {
        return 0.; // can't handle empty bins
      } else {
        return (cvmodule->logn(A1) - cvmodule->logn(A0))
          / (widths[n] * 2.);
      }
    } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
      ix[n]--;
      A0 = value(ix) + offset;
      ix = ix0;
      ix[n]++;
      A1 = value(ix) + offset;
      if (A0 * A1 == 0) {
        return 0.; // can't handle empty bins
      } else {
        return (cvmodule->logn(A1) - cvmodule->logn(A0))
          / (widths[n] * 2.);
      }
    } else {
      // edge: use 2nd order derivative
      int increment = (ix[n] == 0 ? 1 : -1);
      // move right from left edge, or the other way around
      A0 = value(ix) + offset;
      ix[n] += increment; A1 = value(ix) + offset;
      ix[n] += increment; A2 = value(ix) + offset;
      if (A0 * A1 * A2 == 0) {
        return 0.; // can't handle empty bins
      } else {
        return (-1.5 * cvmodule->logn(A0) + 2. * cvmodule->logn(A1)
          - 0.5 * cvmodule->logn(A2)) * increment / widths[n];
      }
    }
  }
 
 
  inline cvm::real gradient_finite_diff(const std::vector<int> &ix0,
                                        int n = 0)
  {
    cvm::real Am1, A0, A1, A2;
    std::vector<int> ix = ix0;
 
    // FIXME this can be rewritten more concisely with wrap_edge()
    if (periodic[n]) {
      ix[n]--; wrap(ix);
      Am1 = value(ix);
      ix = ix0; wrap(ix);
      A0 = value(ix); // Used just to detect gaps
      ix[n]++; wrap(ix);
      A1 = value(ix);
      if (Am1 * A0 * A1 == 0) {
        return 0.; // can't handle empty bins
      } else {
        return (A1 - Am1) / (widths[n] * 2.);
      }
    } else if (ix[n] > 0 && ix[n] < nx[n]-1) { // not an edge
      ix[n]--;
      Am1 = value(ix);
      ix = ix0;
      A0 = value(ix);
      ix[n]++;
      A1 = value(ix);
      if (Am1 * A0 * A1 == 0) {
        return 0.; // can't handle empty bins
      } else {
        return (A1 - Am1) / (widths[n] * 2.);
      }
    } else {
      // edge: use 2nd order derivative
      int increment = (ix[n] == 0 ? 1 : -1);
      // move right from left edge, or the other way around
      A0 = value(ix);
      ix[n] += increment; A1 = value(ix);
      ix[n] += increment; A2 = value(ix);
      return (-1.5 * A0 + 2. * A1
          - 0.5 * A2) * increment / widths[n];
    }
  }
 
 
  virtual inline cvm::real value_output(std::vector<int> const &ix,
                                        size_t const &imult = 0) const override
  {
    int s;
    if (imult > 0) {
      cvm::error_static("Error: trying to access a component "
                 "larger than 1 in a scalar data grid.\n");
      return 0.;
    }
    if (samples) {
      return ( (s = samples->value(ix)) > 0) ?
        (data[address(ix) + imult] / cvm::real(s)) :
        0.0;
    } else {
      return data[address(ix) + imult];
    }
  }
 
  virtual void value_input(std::vector<int> const &ix,
                           cvm::real const &new_value,
                           size_t const &imult = 0,
                           bool add = false) override
  {
    if (imult > 0) {
      cvm::error_static("Error: trying to access a component "
                 "larger than 1 in a scalar data grid.\n");
      return;
    }
    if (add) {
      if (samples)
        // Special case if samples == 0 and value != 0
        data[address(ix)] += new_value * samples->new_value(ix);
      else
        data[address(ix)] += new_value;
    } else {
      if (samples)
        data[address(ix)] = new_value * samples->value(ix);
      else
        data[address(ix)] = new_value;
    }
    has_data = true;
  }
 
  cvm::real maximum_value() const;
 
  cvm::real minimum_value() const;
 
  cvm::real minimum_pos_value() const;
 
  cvm::real integral() const;
 
  cvm::real entropy() const;
 
  cvm::real grid_rmsd(colvar_grid_scalar const &other_grid) const;
};
 
 
 
class colvar_grid_gradient : public colvar_grid<cvm::real>
{
public:
 
  std::shared_ptr<colvar_grid_count> samples;
 
  colvar_grid_gradient();
 
  virtual ~colvar_grid_gradient()
  {}
 
  // /// Constructor from specific sizes arrays
  // colvar_grid_gradient(std::vector<int> const &nx_i);
 
  // /// Constructor from a vector of colvars
  // colvar_grid_gradient(std::vector<colvar *>  &colvars,
  //                      std::string config = std::string());
 
  colvar_grid_gradient(std::string const &filename, std::shared_ptr<colvar_grid_count> samples_in = nullptr);
 
  colvar_grid_gradient(std::vector<colvar *> &colvars,
                       std::shared_ptr<colvar_grid_count> samples_in = nullptr,
                       std::shared_ptr<const colvar_grid_params> params = nullptr,
                       std::string config = std::string());
 
  int full_samples;
  int min_samples;
 
  std::string get_state_params() const;
 
  int parse_params(std::string const &conf,
                   colvarparse::Parse_Mode const parse_mode = colvarparse::parse_normal);
 
  std::istream & read_restart(std::istream &is);
 
  cvm::memory_stream & read_restart(cvm::memory_stream &is);
 
  std::ostream & write_restart(std::ostream &os);
 
  cvm::memory_stream & write_restart(cvm::memory_stream &os);
 
  std::istream &read_raw(std::istream &is);
 
  cvm::memory_stream &read_raw(cvm::memory_stream &is);
 
  std::ostream &write_raw(std::ostream &os, size_t const buf_size = 3) const;
 
  cvm::memory_stream &write_raw(cvm::memory_stream &os, size_t const buf_size = 3) const;
 
  std::istream & read_multicol(std::istream &is, bool add = false);
 
  int read_multicol(std::string const &filename,
                            std::string description = "grid file",
                            bool add = false);
 
  std::ostream & write_multicol(std::ostream &os) const;
 
  int write_multicol(std::string const &filename,
                             std::string description = "grid file") const;
 
  std::ostream & write_opendx(std::ostream &os) const;
 
  int write_opendx(std::string const &filename,
                           std::string description = "grid file") const;
 
  inline void vector_value(std::vector<int> const &ix, std::vector<cvm::real> &v) const
  {
    cvm::real const * p = &value(ix);
    if (samples) {
      int count = samples->value(ix);
      if (count) {
        cvm::real invcount = 1.0 / count;
        for (size_t i = 0; i < mult; i++) {
          v[i] = invcount * p[i];
        }
      } else {
        for (size_t i = 0; i < mult; i++) {
          v[i] = 0.0;
        }
      }
    } else {
      for (size_t i = 0; i < mult; i++) {
        v[i] = p[i];
      }
    }
  }
 
 
  inline void acc_value(std::vector<int> const &ix, std::vector<colvarvalue> const &values) {
    for (size_t imult = 0; imult < mult; imult++) {
      data[address(ix) + imult] += values[imult].real_value;
    }
    if (samples)
      samples->incr_count(ix);
  }
 
  inline void acc_force(std::vector<int> const &ix, cvm::real const *forces) {
    for (size_t imult = 0; imult < mult; imult++) {
      data[address(ix) + imult] -= forces[imult];
    }
    if (samples)
      samples->incr_count(ix);
  }
 
  virtual cvm::real value_output(std::vector<int> const &ix,
                                 size_t const &imult = 0) const override
  {
    int s;
    if (samples) {
      return ( (s = samples->value(ix)) > 0) ?
        (data[address(ix) + imult] / cvm::real(s)) :
        0.0;
    } else {
      return data[address(ix) + imult];
    }
  }
 
  inline cvm::real smooth_inverse_weight(cvm::real weight)
  {
    cvm::real fact;
    if ( weight <= min_samples ) {
      fact = 0.0;
    } else if ( weight < full_samples ) {
      fact = (weight - min_samples) / (weight * cvm::real(full_samples - min_samples));
    } else {
      fact = 1.0 / weight;
    }
    return fact;
  }
 
 
  virtual inline cvm::real value_output_smoothed(std::vector<int> const &ix, bool smoothed = true)
  {
    cvm::real weight, fact;
 
    if (samples) {
      weight = cvm::real(samples->value(ix));
    } else {
      weight = 1.;
    }
 
    if (smoothed) {
      fact = smooth_inverse_weight(weight);
    } else {
      fact = weight > 0. ? 1. / weight : 0.;
    }
 
    return fact * data[address(ix)];
  }
 
 
  inline void vector_value_smoothed(std::vector<int> const &ix, cvm::real *grad, bool smoothed = true)
  {
    cvm::real weight, fact;
 
    if (samples) {
      weight = cvm::real(samples->value(ix));
    } else {
      weight = 1.;
    }
 
    if (smoothed) {
      fact = smooth_inverse_weight(weight);
    } else {
      fact = weight > 0. ? 1. / weight : 0.;
    }
    cvm::real *p = &(data[address(ix)]);
 
    for (size_t imult = 0; imult < mult; imult++) {
      grad[imult] = fact * p[imult];
    }
  }
 
  virtual void value_input(std::vector<int> const &ix,
                           cvm::real const &new_value,
                           size_t const &imult = 0,
                           bool add = false) override
  {
    if (add) {
      if (samples)
        data[address(ix) + imult] += new_value * samples->new_value(ix);
      else
        data[address(ix) + imult] += new_value;
    } else {
      if (samples)
        data[address(ix) + imult] = new_value * samples->value(ix);
      else
        data[address(ix) + imult] = new_value;
    }
    has_data = true;
  }
 
 
  inline cvm::real average(bool smoothed = false)
  {
    if (nd != 1 || nx[0] == 0) {
      return 0.0;
    }
 
    cvm::real sum = 0.0;
    for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix)) {
      sum += value_output_smoothed(ix, smoothed);
    }
 
    return (sum / cvm::real(nx[0]));
  }
 
  cvm::real grid_rmsd(colvar_grid_gradient const &other_grid) const;
 
  void write_1D_integral(std::ostream &os);
 
};
 
 #endif