/*
 * BondVectors.hpp
 *
 *  Created on: Jun 13, 2017
 *      Author: heber
 */


#ifndef DYNAMICS_BONDVECTORS_HPP_
#define DYNAMICS_BONDVECTORS_HPP_

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

#include <map>
#include <vector>

#include <boost/function.hpp>

#include "CodePatterns/Assert.hpp"

#include "LinearAlgebra/Vector.hpp"

#include "Bond/bond.hpp"

/** This class represents all bond vectors, i.e. the normalized direction
 * along a list of bonds, and provides means to extract them from a set of
 * atoms such that for an arbitrary bond the vector can be quickly retrieved.
 */
class BondVectors
{
public:
  //!> typedef for the internal container of the bonds
  typedef std::vector<bond::ptr> container_t;

  //!> typedef for the association of bonds to bond vectors
  typedef std::map<bond::ptr, Vector> mapped_t;

  /** Default cstor for class BondVectors.
   *
   */
  BondVectors();

  /** Prepares the internal container from the bonds of a range of atoms.
   *
   * \param _start start of range
   * \param _end end of range
   * \param _step time step to request bonds for
   */
  template <class T>
  void setFromAtomRange(
      typename T::iterator _start,
      typename T::iterator _end,
      const size_t &_step);

  /** Returns the number of bonds whose bond vectors have been registered.
   *
   * \param number of bonds
   */
  size_t size() const
  { return container.size(); }

  /** Getter for the sorted bonds.
   *
   * \return const ref to internal container
   */
  const container_t& getSorted() const;

  /** Getter for the Bondvectors.
   *
   * \param _step time step for which the bond vector is request
   * \return a map from bond to bond vector
   */
  const mapped_t& getBondVectorsAtStep(const size_t &_step) const;

  /** Get the position in the internal container for a specific bond.
   *
   * \param _bond given bond
   * \return position in the vector, -1 if not present
   */
  size_t getIndexForBond(const bond::ptr &_bond) const;

  /** Gather the subset of BondVectors for the given atom.
   *
   * \param _walker atom to get BondVectors for
   * \param _step time step for which the bond vector is request
   */
  std::vector<Vector> getAtomsBondVectorsAtStep(
      const atom &_walker,
      const size_t &_step) const;

  //!> typedef for the weights for the Bondvectors of a single atom
  typedef std::deque<double> weights_t;

  /** Calculates the weights for a frame where each Bondvector of the
   * given atom is a vector of the frame.
   *
   * The idea is that we can represent any vector by appropriate weights such
   * that is still sums up to one.
   *
   * \param _walker atom to get BondVectors for
   * \param _bondvectors precalculated bond vectors for given \a _walker
   * \param _step time step for which the bond vector is request
   */
  weights_t getWeightsForAtomAtStep(
      const atom &_walker,
      const std::vector<Vector> &_bondvectors,
      const size_t &_step) const;

  /** Calculates the weights for a frame where each Bondvector of the
   * given atom is a vector of the frame.
   *
   * The idea is that we can represent any vector by appropriate weights such
   * that is still sums up to one.
   *
   * \param _walker atom to get BondVectors for
   * \param _step time step for which the bond vector is request
   */
  weights_t getWeightsForAtomAtStep(
      const atom &_walker,
      const size_t &_step) const;

  /** Function typedef to store the bond gradient into a specific container
   * depending on the atom, its current bond and the time step.
   */
  typedef boost::function<void (
      const atom &,
      const bond::ptr &,
      const size_t &,
      const double)> forcestore_t;

  /** Function calculates the remaining part of the atomic gradient that is
   * not captured by the sum of the force along the Bond Vectors.
   *
   * \param _walker atom to get BondVectors for
   * \param _walkerGradient gradient of atom to get BondVectors for
   * \param _BondVectors precalculated bond vectors for given \a _walker
   * \param _weights weight per bond vector (as it is a frame, not a basis)
   * \param _step time step for which the bond vector is request
   * \param _forcestore additional function which may be used to store each
   *        calculated bond force in a bound container
   */
  Vector getRemnantGradientForAtomAtStep(
      const atom &_walker,
      const Vector &_walkerGradient,
      const std::vector<Vector> _BondVectors,
      const BondVectors::weights_t &_weights,
      const size_t &_step,
      forcestore_t _forcestore) const;

  /** Calculates the \a _walkkerGradient projected onto the bond vector for
   * every of \a _walker's bonds.
   *
   * \param _walker atom
   * \param _walkerGradient atom's gradient
   * \param _timestep time step
   * \param _projected_forces list of forces for every \a whichatom_t and every bond
   */
  void getProjectedGradientsForAtomAtStep(
      const atom &_walker,
      const Vector &_walkerGradient,
      const size_t _timestep,
      std::vector< std::vector<double> > &_projected_forces) const;

  // knowing the number of bonds in total, we can setup the storage for the
  // projected forces
  enum whichatom_t {
    leftside=0,
    rightside=1,
    MAX_sides
  };

private:
  /** Calculates the bond vector for each bond in the internal container.
   *
   * \param _step time step for which the bond vector is request
   */
  void recalculateBondVectorsAtStep(const size_t &_step) const;

  /** Helper function to check whether weights sum up to one for each
   * Bond Vector.
   *
   * \param _walker atom to get BondVectors for
   * \param _BondVectors precalculated bond vectors for given \a _walker
   * \param _weights weight per bond vector (as it is a frame, not a basis)
   * \param _step time step for which the bond vector is request
   */
  bool getCheckWeightSumForAtomAtStep(
      const atom &_walker,
      const std::vector<Vector> _BondVectors,
      const BondVectors::weights_t &_weights,
      const size_t &_step) const;

private:
  //!> internal container for sorted bonds
  container_t container;

  //!> states whether map needs update or not
  mutable bool map_is_dirty;

  //!> contains the step for which the map was calculated
  mutable size_t current_step_for_map;

  //!> internal map for bond Bondvector association
  mutable mapped_t current_mapped_vectors;

};

#include "BondVectors_impl.hpp"

#endif /* DYNAMICS_BONDVECTORS_HPP_ */