/*
 * SamplingGrid.hpp
 *
 *  Created on: 25.07.2012
 *      Author: heber
 */

#ifndef SAMPLINGGRID_HPP_
#define SAMPLINGGRID_HPP_

// include config.h
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include <boost/function.hpp>
#include <iosfwd>
#include <vector>

#include "boost/serialization/export.hpp"
#include "boost/serialization/vector.hpp"

#include "LinearAlgebra/defs.hpp"

#include "Fragmentation/Summation/SetValues/SamplingGridProperties.hpp"
#include "Fragmentation/Summation/ZeroInstance.hpp"

class MPQCData;
class SamplingGridTest;

/** This class stores a sample function on a three-dimensional grid.
 *
 *  \note We do not use boost::multi_array because it is not trivial to serialize.
 *
 */
class SamplingGrid : public SamplingGridProperties
{
  //!> grant unit test access to private parts
  friend class SamplingGridTest;
  //!> grant output operator access
  friend std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other);
public:
  //!> typedef for sampled values
  typedef std::vector< double > sampledvalues_t;

  /** Constructor for class SamplingGrid for full window.
   *
   * Here, the window of sampled values spans the given domain.
   *
   * \param _begin offset for grid per axis
   * \param _end edge length of grid per axis
   * \param _level number of grid points in \f$2^{\mathrm{level}}\f$
   * \param _sampled_grid sample points
   */
  SamplingGrid(const double _begin[NDIM],
      const double _end[NDIM],
      const int _level,
      const sampledvalues_t &_sampled_grid);

  /** Constructor for class SamplingGrid for empty window.
   *
   * Here, the window is initially of size zero.
   *
   * \param _begin offset for grid per axis
   * \param _end edge length of grid per axis
   * \param _level number of grid points in \f$2^{\mathrm{level}}\f$
   */
  SamplingGrid(const double _begin[NDIM],
      const double _end[NDIM],
      const int _level);

  /** Copy constructor for class SamplingGrid with full window from SamplingGridProperties.
   *
   * Here, the window is initially empty.
   *
   * \param _props properties to copy
   */
  SamplingGrid(const SamplingGridProperties &_props);

  /** Copy constructor for class SamplingGrid with empty window from SamplingGridProperties.
   *
   * Here, the window must span the whole domain
   *
   * \param _props properties to copy
   * \param _sampled_grid sample points
   */
  SamplingGrid(
      const SamplingGridProperties &_props,
      const sampledvalues_t &_sampled_grid);

  /** Copy constructor for class SamplingGrid.
   *
   * The window of sampled values corresponds to the one on \a _grid.
   *
   * \param _grid grid to copy
   */
  SamplingGrid(const SamplingGrid &_grid);

  /** default cstor.
   */
  SamplingGrid();

  virtual ~SamplingGrid();

  /** Checks whether another instance is consistent with this one.
   *
   *  \note Conistency is stronger as grids must have the same window.
   *
   * \param _props other properties to check against
   * \return true - are consistent, false - else
   */
  bool isCongruent(const SamplingGrid &_props) const;

  /** Assignment operator.
   *
   * \param other other instance to assign ourselves to
   */
  SamplingGrid& operator=(const SamplingGrid& other);

  /** Addition operator with another SamplingGrid instance \a other.
   *
   * \param other other instance to sum onto this one.
   * \return ref to this instance
   */
  SamplingGrid& operator+=(const SamplingGrid& other)
  {
    superposeOtherGrids(other, +1.);
    return *this;
  }

  /** Element-wise multiplication operator with another SamplingGrid instance \a other.
   *
   * With non-zero windows we have to pay some more attention here.
   * Now, the windows may not be congruent but we have to find the intersection
   * of the two windows and then construct the new window only of this size and
   * multiply. The trick then is not to copy&change the other grid but to
   * access it properly.
   *
   * \param other other instance to sum onto this one.
   * \return ref to this instance
   */
  SamplingGrid& operator*=(const SamplingGrid& other);

  /** Scaling values in SamplingGrid instance by \a _value.
   *
   * With non-zero windows we have to pay some more attention here.
   * Now, the windows may not be congruent but we have to find the intersection
   * of the two windows and then construct the new window only of this size and
   * multiply. The trick then is not to copy&change the other grid but to
   * access it properly.
   *
   * \param other other instance to sum onto this one.
   * \return ref to this instance
   */
  SamplingGrid& operator*=(const double _value);

  /** Subtraction operator with another SamplingGrid instance \a other.
   *
   * \param other other instance to subtract from this one.
   * \return ref to this instance
   */
  SamplingGrid& operator-=(const SamplingGrid& other)
  {
    superposeOtherGrids(other, -1.);
    return *this;
  }

  /** Sample a given grid \a other down to grid level \a _level and store
   * in \a instance.
   *
   * \param instance given instance to store downsampled grid in
   * \param other instance to get grid from
   * \param _level level to sample down to
   */
  static void downsample(SamplingGrid& instance, const SamplingGrid& other, const int _level);

  /** Returns the numeric integral over the grid.
   *
   * @return sum of grid values times volume element
   */
  double integral() const;

  /** Returns the numeric integral over the grid where the grid is element-wise multiplied with \a weight.
   *
   * @param weight grid of weights
   * @return sum of grid values weighted by respective element in weight times volume element
   */
  double integral(const SamplingGrid &weight) const;

  /** Returns the total number of gridpoints of the discrete mesh covering the (window) volume.
   *
   * @return number of gridpoints sampled_values should have
   */
  const size_t getWindowGridPoints() const;

  /** Returns the number of gridpoints of the discrete mesh for the current
   * window size for given axis \axis.
   *
   * \param axis axis to calculate number of gridpoints for
   * \return number of gridpoints along this axis
   */
  const size_t getWindowGridPointsPerAxis(const size_t axis) const;

  /** Returns the length of the window for the given \a axis.
   *
   * \param axis axis for which to get step length
   * \return window length for the given axis, i.e. end - begin
   */
  const double getWindowLengthPerAxis(const size_t axis) const;

  /** Returns the discrete length in grid cells of the window for the given \a axis.
   *
   * \param axis axis for which to get step length
   * \return window length in grid cells for the given axis
   */
  const size_t getDiscreteWindowLengthPerAxis(const size_t axis) const;

  /** Returns the volume of the domain covered by the current window.
   *
   * @return volume of window
   */
  const double getWindowVolume() const;

  /** Sets the size of the window.
   *
   * \note also resets the sampled points so far.
   *
   * \param _begin_window start of new window
   * \param _end_window end of window
   */
  void setWindow(const double _begin_window[NDIM], const double _end_window[NDIM]);

  /** Helper function to convert begin_window and end_window that are w.r.t.
   * to domain [begin:end] to indices that can be used when traversing the grid.
   *
   * \param larger_wbegin begin of domain
   * \param larger_wend end of domain
   * \param smaller_wbegin begin of window
   * \param smaller_wend end of window
   * \param pre_offset discrete length from 0 to start of window
   * \param post_offset discrete length from end of window to end
   * \param length discrete length of window
   * \param total total number of points for checking, should be sum of other three
   */
  void getDiscreteWindowCopyIndices(
      const double *larger_wbegin,
      const double *larger_wend,
      const double *smaller_wbegin,
      const double *smaller_wend,
      size_t *pre_offset,
      size_t *post_offset,
      size_t *length,
      size_t *total) const;

  /** Returns begin, length and end of window relative to the full domain in discrete
   * grid points, i.e. begin gives the first grid points with the window and end its
   * last plus 1.
   */
  void getDiscreteWindowIndices(
      size_t _wbegin[NDIM],
      size_t _wlength[NDIM],
      size_t _wend[NDIM]) const;

  /** Returns number of grid points before the window, during the window, and
   * after the window including the total length of the domain for check.
   *
   * \param _pre_offset grid points before start of window
   * \param _post_offset grid points after end of window
   * \param _length grid points in window
   * \param _total grid points in domain
   */
  void getDiscreteWindowOffsets(
      size_t _pre_offset[NDIM],
      size_t _post_offset[NDIM],
      size_t _length[NDIM],
      size_t _total[NDIM]) const;

  /** Equality operator.
   *
   * @param other other instance to check against
   * @return true - both are equal, false - grids differ
   */
  bool operator==(const SamplingGrid& other) const;

  bool operator!=(const SamplingGrid& other) const
  {
    return (!(*this == other));
  }

private:
  /** Extend the window such that the number of sample points stored is an
   * even number per axis. This is used by downsample()
   */
  void padWithZerosForEvenNumberedSamples();

  /** Sets the size of the domain.
   *
   * \note also resets the sampled points so far and the window.
   *
   * \param _begin start of new window
   * \param _end end of window
   */
  void setDomain(const double _begin[NDIM], const double _end[NDIM]);

  /** Sets the size of the domain.
   *
   * \note this is just internally used for easing the array setting.
   *
   * \param _begin start of domain
   * \param _end end of domain
   */
  void setDomainSize(const double _begin[NDIM], const double _end[NDIM]);

  /** Extends the window while keeping the values.
   *
   * \param _begin_window new start of window
   * \param _end_window new end of window
   */
  void extendWindow(const double _begin_window[NDIM], const double _end_window[NDIM]);

  /** Shrinks the window while keeping the values.
   *
   * \param _begin_window new start of window
   * \param _end_window new end of window
   */
  void shrinkWindow(const double _begin_window[NDIM], const double _end_window[NDIM]);

  /** Adds another (smaller) window onto the one in this instance.
   *
   * \note We assume here that the given window fits on the this one.
   *
   * \param _begin_window start of other window
   * \param _end_window end of other window
   * \param _sampled_grid other set of sampled values
   * @param prefactor +1. is then addition, -1. is subtraction.
   */
  void addOntoWindow(
      const double _begin_window[NDIM],
      const double _end_window[NDIM],
      const sampledvalues_t &_sampled_grid,
      const double prefactor);

  /** Adds another (larger) window into the one in this instance.
   *
   * \note We assume here that the given window is larger than this one.
   *
   * \param _begin_window start of other window
   * \param _end_window end of other window
   * \param _sampled_grid other set of sampled values
   * @param prefactor +1. is then addition, -1. is subtraction.
   */
  void addIntoWindow(
      const double _begin_window[NDIM],
      const double _end_window[NDIM],
      const sampledvalues_t &_sampled_grid,
      const double prefactor);

  /** Enum to help in addWindowOntoWindow() decide which iterator needs to be
   * advanced.
   */
  enum eLargerWindow {
    destwindow,
    sourcewindow
  };

  /** Helper function to copy one (larger) window into a (smaller) window.
   *
   * \note Why do we need the extra \a choice? We need to know which window
   * tuples is associated with which sampled values that are constrained by
   * one of them being constant, hence the source values
   *
   * \param larger_wbegin start of larger window
   * \param larger_wend end of larger window
   * \param smaller_wbegin start of smaller window
   * \param smaller_wend end of smaller window
   * \param dest_sampled_grid larger set of sampled values
   * \param source_sampled_grid smaller set of sampled values
   * \param op operation to perform with the two elements
   * \param larger_window indicates which is the larger window
   */
  void addWindowOntoWindow(
      const double larger_wbegin[NDIM],
      const double larger_wend[NDIM],
      const double smaller_wbegin[NDIM],
      const double smaller_wend[NDIM],
      sampledvalues_t &dest_sampled_grid,
      const sampledvalues_t &source_sampled_grid,
      boost::function<void (double &, const double &)> op,
      enum eLargerWindow larger_window);

  /** Helper function that contains all the logic of how to superpose two
   * grids.
   *
   * Is called by SamplingGrid::operator+=() and SamplingGrid::operator-=()
   *
   * @param other other histogram
   * @param prefactor +1. is then addition, -1. is subtraction.
   */
  void superposeOtherGrids(const SamplingGrid &other, const double prefactor);

  /** Sets the size of the window.
   *
   * \note also resets the sampled points so far.
   *
   * \param _begin_window start of new window
   * \param _end_window end of window
   */
  void setWindowSize(const double _begin_window[NDIM], const double _end_window[NDIM]);

public:
  /// We do not store the whole grid if many entries are actually zero
  /// but only a window wherein the sampled function is non-zero.

  //!> sample points of the window
  sampledvalues_t sampled_grid;

  //!> start of the window relative to SamplingGridProperties::begin and SamplingGridProperties::size
  double begin_window[NDIM];
  //!> end of the window relative to SamplingGridProperties::begin and SamplingGridProperties::size
  double end_window[NDIM];

private:
  friend class MPQCData;

  friend class boost::serialization::access;
  // serialization
  template <typename Archive>
  void serialize(Archive& ar, const unsigned int version)
  {
    ar & boost::serialization::base_object<SamplingGridProperties>(*this);
    ar & const_cast< sampledvalues_t &>(sampled_grid);
    for(size_t i=0;i<NDIM;++i) {
      ar & begin_window[i];
      ar & end_window[i];
    }
  }

  //!> static typedef to use in cstor when no initial values are given
  static const double zeroOffset[NDIM];
};

/** Output operator for class SamplingGrid.
 *
 * \param ost output stream to print to
 * \param other instance to print
 * \return ref to stream for concatenation
 */
std::ostream & operator<<(std::ostream &ost, const SamplingGrid& other);

template<> SamplingGrid ZeroInstance<SamplingGrid>();

// we need to give this class a unique key for serialization
// its is only serialized through its base class FragmentJob
BOOST_CLASS_EXPORT_KEY(SamplingGrid)

// define inline functions
#include "SamplingGrid_inline.hpp"

#endif /* SAMPLINGGRID_HPP_ */