#ifndef __FRAME_FIELDS__
#define __FRAME_FIELDS__

// ========================================================================== //
// Local libraries

#include "Matrix/Matrix.hpp"
#include "Utilities/BBoxMap.h"

// External libraries
#include <GEL/CGLA/CGLA.h>

// Standard libraries
#include <vector>

// ========================================================================== //
// Start of the parent class

// Frame Field Parent Class
//
// This class define the necessary framework for handling frame fields on a 3D
// domain. Specializations of this class can introduce underlying data
// structures, such as grid aligned fields.
//
class FrameField {
  public:
    // ====================================================================== //
    // Data Specific operations. These should be defined in the inherited
    // classes.

    /** File reading function
     *
     * @param data_location folder containing the input files.
     */
    virtual void read(std::string data_location);

    /** Get the field index.
     *
     * @param position
     * @return element index at position.
     */
    virtual int get_field_index(const CGLA::Vec3d &position);

    /** Function designed to detect singularities in the frames.
     *
     * This is a costly operation. The list of singular elements are saved to
     * the specified filename to minimize later computations.
     *
     * @param filename
     */
    virtual void detect_singular(std::string filename);

    /** Exporting the detected singularities.
     *
     * @param export_location
     */
    virtual void export_singularity(std::string export_location);

    double element_size = 0.0; // Size of elements
  protected:
    // ====================================================================== //
    // Design domain variables

    int    dimF  = 3; // Number of frame vectors.
    size_t numel = 0; // Total number of elements
    BBox   bbox;      // Bounding box

    std::vector<CGLA::Vec3d>              pos_list;   // Frame locations
    std::vector<std::vector<CGLA::Vec3d>> frame_list; // Frame direction field
    std::vector<std::vector<double>>      layer_list; // Layer thicknesses
    std::vector<double>                   densi_list; // Densities
    std::vector<size_t>                   singular_elements; // Singularity

    virtual std::vector<size_t> neighbours(const CGLA::Vec3d &position);
    virtual std::vector<double> get_local(const CGLA::Vec3d &position);

  public:
    // ====================================================================== //
    // Public API

    FrameField();

    /** Read the list of singular elements from the file.
     *
     * We assume the default naming is used.
     *
     * @param file_location
     */
    void read_singular(std::string file_location);
    /** Export function for density list.
     *
     * @param export_location
     */
    void export_densities(std::string export_location);

    /** Get the bbox object */
    BBox get_bbox() { return this->bbox; }

    // ====================================================================== //
    // Indexing functions

    int lookup_frame_index(const CGLA::Vec3d &position,
                           const CGLA::Vec3d &direction);
    int interpolate_frame_index(const CGLA::Vec3d &position,
                                const CGLA::Vec3d &direction);
    int nearest_index(const std::vector<CGLA::Vec3d> &field,
                      const CGLA::Vec3d &             d);

    // ====================================================================== //
    // Direct lookup functions

    /** Direct lookup of frame at position.
     *
     * @param position
     * @return List of frame directions
     */
    std::vector<CGLA::Vec3d> lookup_frame(const CGLA::Vec3d &position);

    /** Direct lookup of layer at position.
     *
     * @param position
     * @return List of layer thicknesses
     */
    std::vector<double> lookup_layer(const CGLA::Vec3d &position);

    /** Direct lookup of density at position.
     *
     * @param position
     * @return Density value
     */
    double lookup_density(const CGLA::Vec3d &position);

    /** Direct lookup of nearest frame location.
     *
     * @param position
     * @return Location
     */
    CGLA::Vec3d lookup_position(const CGLA::Vec3d &position);

    // ====================================================================== //
    // Directional lookup functions

    /** Directional lookup function for frames.
     *
     * This function will do a lookup of the frame system at a given position.
     * In addition to the lookup, we determine which frame is nearest a
     * specified direction.
     *
     * @param position of the lookup.
     * @param direction desired direction.
     * @return nearest frame direction, sign corrected.
     */
    CGLA::Vec3d lookup_frame_nearest(const CGLA::Vec3d &position,
                                     const CGLA::Vec3d &direction);

    /** Directional lookup function for layer thicknesses.
     *
     * This function will do a lookup and return the layer thickness of the
     * frame nearest the specified direction.
     *
     * @param position of the lookup.
     * @param direction desired direction.
     * @return layer thickness associated with the nearest frame.
     */
    double lookup_layer_nearest(const CGLA::Vec3d &position,
                                const CGLA::Vec3d &direction);

    // ====================================================================== //
    // Interpolation functions

    virtual double      interpolate(const std::vector<double> &local,
                                    const std::vector<double> &value);
    virtual CGLA::Vec3d interpolate(const std::vector<double> &     local,
                                    const std::vector<CGLA::Vec3d> &value);

    /** Frame interpolation.
     *
     * This function return the interpolated frames at a given point.
     *
     * @param position
     * @return list of frame vectors
     */
    std::vector<CGLA::Vec3d> interpolate_frame(const CGLA::Vec3d &position);
    /** Interpolate layer thickness at given position.
     *
     *
     * @param position
     * @return list of layer thicknesses
     */
    std::vector<double> interpolate_layer(const CGLA::Vec3d &position);
    /** Interpolate the density at given position.
     *
     * @param position
     * @return density value
     */
    double interpolate_density(const CGLA::Vec3d &position);

    // ====================================================================== //
    // Directional interpolation functions

    /** Directional lookup function for frames.
     *
     * This function will do a lookup of the frame system at a given position.
     * In addition to the lookup, we determine which frame is nearest a
     * specified direction.
     *
     * @param position of the lookup.
     * @param direction desired direction.
     * @return nearest frame direction, sign corrected.
     */
    CGLA::Vec3d interpolate_frame_nearest(const CGLA::Vec3d &position,
                                          const CGLA::Vec3d &direction);

    /** Directional lookup function for layer thicknesses.
     *
     * This function will do a lookup and return the layer thickness of the
     * frame nearest the specified direction.
     *
     * @param position of the lookup.
     * @param direction desired direction.
     * @return layer thickness associated with the nearest frame.
     */
    double interpolate_layer_nearest(const CGLA::Vec3d &position,
                                     const CGLA::Vec3d &direction);

  protected:
    std::vector<std::vector<int>>
    align_index(const std::vector<size_t> &idx_list);

    // Index based lookup of elements
    inline const std::vector<CGLA::Vec3d> lookup_frame(const int &idx) {
        return this->frame_list[idx];
    }
    inline const std::vector<double> lookup_layer(const int &idx) {
        return this->layer_list[idx];
    }
    inline const double lookup_density(const int &idx) {
        return this->densi_list[idx];
    }
    inline const CGLA::Vec3d lookup_position(const int &idx) {
        return this->pos_list[idx];
    }
};

// ========================================================================== //
// Grid based frame fields

/** Specialized class for axis aligned grid of frames.
 *
 * This class handle an axis aligned grid of frames. These are assumed to be
 * located input cubes. The number of cubes along each dimension should be
 * specified in the information.csv.
 * Each line in the main file 'field.csv' should contain first the index of each
 * point in base 1, then the frame vectors.
 */
class FrameGrid : public FrameField {

  private:
    // Data Fields
    int            dimX, dimY, dimZ; // Number of elements in each dimension.
    Matrix<size_t> frame_grid;       // Lookup Grid

  public:
    // Constructor and IO
    FrameGrid();
    FrameGrid(int X, int Y, int Z, int F);

    /** Reading of frame field files.
     *
     * Each line in the main file 'field.csv' should contain first the index of
     * each point in base 1, then the frame vectors.
     *
     * @param data_location
     */
    void read(std::string data_location) override;

    // Interpolators

    int  get_field_index(const CGLA::Vec3d &position) override;
    void detect_singular(std::string filename) override;
    void export_singularity(std::string export_location) override;

  private:
    std::vector<size_t> neighbours(const CGLA::Vec3d &position) override;
    std::vector<double> get_local(const CGLA::Vec3d &position) override;

    double      interpolate(const std::vector<double> &local,
                            const std::vector<double> &value) override;
    CGLA::Vec3d interpolate(const std::vector<double> &     local,
                            const std::vector<CGLA::Vec3d> &value) override;
};

#endif /* __FRAME_FIELDS__ */