/***************************************************************************
 *   Copyright (C) 2005 by Nicolas PASCAL                                  *
 *   nicolas.pascal@icare.univ-lille1.fr                                   *
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 *   This program is distributed in the hope that it will be useful,       *
 *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
 *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
 *   GNU General Public License for more details.                          *
 *                                                                         *
 *   You should have received a copy of the GNU General Public License     *
 *   along with this program; if not, write to the                         *
 *   Free Software Foundation, Inc.,                                       *
 *   59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.             *
 ***************************************************************************/

/***************************************************************************
 *  History :                                                              *
 *  16/08/05 : start                                                        *
 ***************************************************************************/

#ifndef TOOLS_H
#define TOOLS_H

#include <iostream>
#include <string>
#include <list>
#include <vector>
#include <sstream>
#include <ctype.h>
#include <cmath>
#include <cfloat>
#include <typeinfo>
#include <algorithm>
#include <cassert>
#include <limits>

#include "tools_exceptions.h"
#include "geometry_common.h"

// #define isnan(x) ((x) != (x))

namespace MyTools {
    /**
     * @brief decompose an integer in tens digits, by increasing tens power (v[0] : units, v[1] : tens, v[2] : hundreds...)
     * For example, 1977 will return a vector of [7, 7, 9, 1]
     * @param n an unsigned integer
     * @return a vector of decimal digits, sorted by increasing tens power
     */
    vector <unsigned short> get_digits (unsigned long long n);

    /**
     * @brief compute the point that is the 2D intersection of the 2 lines (each define by 2 different points)
     * If the lines are parallel (or superimposed), the method mill return false, and the intersection point will be {inf,inf}
     * source of the algorithm : http://astronomy.swin.edu.au/~pbourke/geometry/lineline2d/
     * @param l1 the 2 points of the first line
     * @param l2 the 2 points of the 2nd line
     * @param intersection the intersection point. Must be a 1D vector
     * @return true if an intersection has been found (false if the lines are parallel)
     */
    template <typename NumericType>
    const bool get_intersection(const std::vector < std::vector <NumericType> > &l1, const std::vector < std::vector <NumericType> > &l2, std::vector< NumericType > &intersection) {
          // extract the coordinates of the 2 points -> use double to be as precise as possible
        double x1=l1[0][0];double y1=l1[0][1];
        double x2=l1[1][0];double y2=l1[1][1];
        double x3=l2[0][0];double y3=l2[0][1];
        double x4=l2[1][0];double y4=l2[1][1];
        // some denominators
        double denom=((y4-y3)*(x2-x1)-(x4-x3)*(y2-y1));
        // inetrsection coordinates
        double x=DBL_MAX;
        double y=DBL_MAX;

        if (denom==0.)
            // the 2 lines are parallels
            return false;
        double div_denom=1./denom;

        double ua=((x4-x3)*(y1-y3)-(y4-y3)*(x1-x3))*div_denom;

        // set intersection
        x=x1+ua*(x2-x1);
        y=y1+ua*(y2-y1);

        // fill the output vector
        intersection.clear();
        intersection.push_back(static_cast<NumericType>(x));
        intersection.push_back(static_cast<NumericType>(y));

        return true;
    }

    /**
     * @brief search for a value in a range
     * @param first1 start of search
     * @param last1 end of search
     * @param value value to search
     * @return the position of the value if found
     */
//     template<class ForwardIterator, class T_Value>
//     ForwardIterator find_first_of ( ForwardIterator first1, ForwardIterator last1, const T_Value &value) {
//       return find_first_of(first1,last1,(T_Value*)(&value),(T_Value*)(&value+1));
//     }

    /**
    * @brief test the system byte order
    */
    inline static bool is_big_endian(){
        short int word = 0x0001;
        char *byte = (char *) &word;
        return(byte[0] ? true : false);
    };
    /**
     * @brief return the endianness of the system
     * @return 0 if little endian, 1 if big endian
     */
    inline static int get_endianness() {
        int i = 1;
        char *p = (char *)&i;
        return (p[0] == 1 ? 0 : 1);
    }
    /**
     * @brief swap the bytes of a
     * @param b the number to swap
     * @param n the number of bytes of the real type of @a b
     */
    inline void byte_swapper(unsigned char* b, int n) {
      register int i = 0;
      register int j = n-1;
      while (i<j) {
        std::swap(b[i], b[j]);
        i++, j--;
      }
    };

    /**
    * @brief Swap data of all C types
    * @param nb the number to swap
    * @return the swapped number
    */
    template<class T>
    inline static void byte_swap (T &nb) {
      byte_swapper((unsigned char*)(&nb),sizeof(nb));
    };

    /**
    convert a std::string into integer, assuming s represent an integer
    */
    int string2int(const std::string &s) ;
    float string2float(const std::string &s) ;
    double string2double(const std::string &s) ;
    int char2int(const char &c);
    std::string char2string(const char &c);
    /**
    * @brief convert an integer to a std::string.
    * If nb_digit is given, some 0 will be added on the left of the integer to reach nb_digit.
    * For example, int2string(12,4) will return "0012" : 2 zeros will be added to represent the input integer with 4 digits.
    * If nb_digit is less than the needed digits to represent the input integer, it is ignored.
    * @param i the integer to be converted
    * @param nb_digit the number of digits of the result string. If -1 (default), ignored
    * @return the string representation of the input integer
    */
    std::string int2string(const int &i, const int &nb_digit = -1) ;
    /**
     * @brief convert a numeric to a string
     * @param nb the number to convert
     * @return @a nb as a string representation
     */
    template<typename NumericType>
    std::string to_string( const NumericType & nb ) {
        std::ostringstream oss;
        oss << nb;
        return oss.str();
    }
    /**
     * @brief convert a string to a numeric.
     * Depending of the compiler, the output numeric type can be required to be explicitally specified as a template argument
     * @param s_nb the number as a string representation
     * @return @a s_nb as a numeric type
     */
    template<typename NumericType>
    NumericType to_num( const string &s_nb ) {
      // generate a flux with the string to be converted
      std::istringstream iss(s_nb);
      NumericType nb ;
      iss >> nb;
      return nb;
    }
    /**
     * @brief print out a bitfield coded in a char type
     * @param bit_field the char or numeric to convert
     * @param size the size of @a bit_field (in bytes)
     * @return the bits string
     */
    std::string to_bit( void* bit_field, const int size=1 );

    /**
     * @brief round a number to the given type
     * @param x the number to round
     * @return the rounded number as a long
     */
    template<typename T_out>
    T_out round(const double &x) {
        return static_cast<T_out>(x + 0.5);
    }
    /**
    Compute the angle from the origin between 2 points, positive anticlockwise.
    @param p1 a 2 value vector containing the carthesian coordinates [x,y] of the first point
    @param p2 a 2 value vector containing the carthesian coordinates [x,y] of the second point
    @return the angle between -pi->pi
    */
    double angle_2D(double* p1, double* p2);

    const bool is_int(const std::string &s);
    const bool is_float( const std::string & s );
    const bool is_int(const char * start, const int sz);
    /**
     * @brief check if the given string represents a data of the type given as template argument
     * This source code is taken from http://c.developpez.com/faq/cpp/?page=strings#STRINGS_is_type
     * @param str the string to test
     * @return true if @a str has the type T
     */
    template<typename T>
    bool has_type( const std::string & str ) {
      std::istringstream iss( str );
      T val;
      // try to fully convert it. If successfull, the input string has the good type
      return ( iss >> val ) && ( iss.eof() );
    }

    const int count( const std::string & elt, const std::string & s );
    const int count( const char elt, const std::string & s );

    void remove_space(std::string &);
    /**
     * @brief remove starting and trailing characters
     * @param str a string
     * @param c the character to strip
     * @return @a str stripped from charcter @a c
     */
    string strip (const std::string &str, const char c = ' ');

    /**
    * average the values of an array. It skips the fill values. If no good values have been found, it returns -LDBL_MAX
    * @param array the array to be averaged
    * @param nb_val the number of values of the array
    * @param fill_value the fill value that represents bad ones
    * @return the average, or -LDBL_MAX
    */
    template <class T>
    const long double avg_array(const T* array, const int& nb_val, const T& fill_value);

    template <class T>
    T abs(const T& val=T(0)) {
        if (val<T(0)) return -val;
        else return val;
    }

    template <class T>
    const long double avg_array(const T* array, const int& nb_val, const T& fill_value) {
        bool good_val_found=false;
        long double acc=0.;
        int nb_good_val=0;
        for (int i = 0 ; i<nb_val ; ++i)
            if (array[i]!=fill_value) {
                if (!good_val_found)
                    good_val_found=true;
                acc+=static_cast<long double>(array[i]);
                ++nb_good_val;
            }
        if (good_val_found)
            return acc/static_cast<long double>(nb_good_val);
        else
            return -LDBL_MAX;
    }

    /**
     * @brief stricly greater than predicate
     * @param i first value to compare
     * @param j second value to compare
     * @return true if ( i > j )
     */
    template<typename T>
    bool gt_comparator ( T i, T j) { return ( i > j ); };
    /**
     * @brief operator function to use with the STL algorithm
     * @param i operand 1
     * @param j operand 2
     * @return the sum of the 2 operands
     */
    template<typename T>
    T sum (const T &i, const T &j) { return i+j; }
    /**
     * @brief operator function to use with the STL algorithm
     * @param i operand 1
     * @param j operand 2
     * @return the product of the 2 operands
     */
    template<typename T>
    T prod (const T &i, const T &j) { return i * j; }
    /**
     * @brief operator function to use with the STL algorithm
     * @param i operand
     * @return the natural logarithm of the operand
     */
    template<typename T>
    T log_op (const T &i) { return log((double)(i)); }

    /**
     * @class KeyComparator pseudo functionnal class used to search a given message key
     */
    class KeyComparator {
    public:
        /** the key to search for */
        string key;
        /**
         * @brief constructor
         * @param key the key to search for
         */
        KeyComparator (const string & key) : key(key){};
        /**
         * @brief comparator function
         * @param key_val a (key, value) pair
         * @return true if the key of @a key_val is equal to the searched key
         */
        bool operator ()(pair <string, string> const& key_val) const {
            return key_val.first == key;
        }
    };
    /**
     * @class NotEqualComparator pseudo functionnal class used for container search
     */
    template <class T>
    class NotEqualComparator {
    public:
        /** the value to ignore */
        T val;
        /**
         * @brief constructor
         * @param key the key to ignore
         */
        NotEqualComparator (const T &val) : val(val){};
        /**
         * @brief comparator function
         * @param ignored_val the value to skip
         * @return true if the value is not equal to ignored_val
         */
        bool operator ()(const T ignored_val) const {
            return val != ignored_val;
        }
    };
    /**
     * @brief computes the abscissa of the value @a v in the vector @a v_abs. STL vgector version
     * If the value @a v is strictly less than the minimum of @a v_abs or striclty greater than the maximum, will return NaN
     * @param val the searched value
     * @param vec the vector of abscissa. Must be in a strict monotonic ordering
     * @return the abscissa as a floating point number
     */
    template<typename T>
    double get_abscissa ( const T val, const vector < T > & vec ) {
        return get_abscissa ( val, vec.begin(), vec.size() );
    };

    /**
     * @brief computes the abscissa of the value @a v in the vector @a v_abs, using a binary search
     * If the value @a v is strictly less than the minimum of @a v_abs or striclty greater than the maximum, will return NaN
     * @param v the searched value
     * @param v_abs the vector of abscissa. Must be in a strict monotonic ordering
     * @param n_abs the number of values of the vector
     * @return the abscissa as a floating point number
     */
    template<typename T>
    double get_abscissa ( const T val, const T * v_abs, const int n_abs ) {
        double invalid = numeric_limits<double>::signaling_NaN();
        if ( n_abs <= 1 )
            return invalid;

        // find insertion position of v in the abscissa vector
        T * v_abs_begin = const_cast < T * > ( v_abs );
        T * v_abs_end   = v_abs_begin + n_abs;

        bool is_ascending_order = ( v_abs [1] > v_abs [0] );
        T * it = v_abs_begin;

        if ( is_ascending_order ) { // ascending order
            it = lower_bound ( v_abs_begin, v_abs_end, val );
        } else { // descending order -> use ">" as comparator
            it = lower_bound ( v_abs_begin, v_abs_end, val, gt_comparator < T > );
        }

        if ( it == v_abs_begin || it == v_abs_end )
            return invalid;

        size_t ilow = ( it - 1 ) - v_abs;
        T lower_bound = * ( it - 1 );
        T upper_bound = * ( it );
// printf ( "bound = ( %f, %f )\tratio = %f\n", lower_bound, upper_bound, ( val - lower_bound ) / ( upper_bound - lower_bound ) );
        return ilow + ( val - lower_bound ) / ( upper_bound - lower_bound );
    };
    /**
     * @brief computes the abscissa of the logarithm of the output coordinates @a x_out in the logarithm of the reference coordinates @a x_in using a sequential search.
     * If the value @a v is strictly less than the minimum of @a v_abs or striclty greater than the maximum, will return NaN
     * @param y_out output abscissa
     * @param x_out searched coordinates
     * @param sz_out number of searched coordinates
     * @param x_in reference coordinates system. Must be strictly monotonic
     * @param sz_in number of reference coordinates
     */
    template <typename Txout, typename Tyout, typename Txin>
    void get_abscissa_sequential (Tyout * y_out, const Txout * x_out, const int sz_out, const Txin * x_in, const int sz_in) {
        // check if x_in are sorted in ascending or descending order
        bool is_asc = ( x_in [1] > x_in [0] );

        // compute the linear abscissa of x_out in x_in coordinates
        Tyout nan = numeric_limits<Tyout>::signaling_NaN();
        int j;
        for (int i = 0 ; i < sz_out ; ++i) {
            j = 0;
            while ((j < (sz_in - 1)) &&
                   ((  is_asc && (x_out[i] > x_in [j+1])) ||
                    (! is_asc && (x_out[i] < x_in [j+1]))
                   ))
                ++j;
            if (j >= (sz_in - 1) || isnan (x_out[i]))
                y_out [i] = nan;
            else {
                y_out [i] = j + (x_out[i] - x_in [j]) / (x_in [j+1] - x_in [j]);
            }
        }
    }

    /**
     * @brief computes the abscissa of the logarithm of the output coordinates @a x_out in the logarithm of the reference coordinates @a x_in using a sequential search.
     * If the value @a v is strictly less than the minimum of @a v_abs or striclty greater than the maximum, will return NaN
     * @param y_out output abscissa
     * @param x_out searched coordinates
     * @param sz_out number of searched coordinates
     * @param x_in reference coordinates system. Must be strictly monotonic
     * @param sz_in number of reference coordinates
     */
    template <typename Txout, typename Tyout, typename Txin>
    void get_log_abscissa_sequential (Tyout * y_out, const Txout * x_out, const int sz_out, const Txin * x_in, const int sz_in) {
        // build log of searched x
        double x_out_log[sz_out];
        for (int i = 0 ; i < sz_out ; ++i )
            x_out_log [i] = log(double(x_out[i]));

        // build log(P) of reference x coordinates
        double x_in_log [sz_in];
        for (int i = 0 ; i < sz_in ; ++i )
            x_in_log [i] = log((double)(x_in[i]));

        get_abscissa_sequential (y_out, x_out_log, sz_out, x_in_log, sz_in);
    }

    /**
     * @brief generate the vector of values between val_min ( included ) and val_max ( excluded ), spaced by inc
     * @param val_min start of range generation
     * @param val_max end of range generation
     * @param inc increment between 2 successive elements
     * @return the vector of values between val_min and val_max spaced by inc, in increasing order
     */
    template<typename T>
    vector <T> arange (const T &val_min, const T &val_max, const T &inc = (T)(1));
    /**
     * @brief generate the vector of values between val_min ( included ) and val_max ( excluded ), spaced by inc
     * @param v [OUT] vector of values between val_min and val_max spaced by inc, in increasing order
     * @param val_min [IN] start of range generation
     * @param val_max [IN] end of range generation
     * @param inc [IN] increment between 2 successive elements. +1 by default
     */
    template<typename T>
    void arange (vector <T> & v, const T &val_min, const T &val_max, const T &inc = (T)(1));
    /**
     * @brief push the content of @a v_in at the end of @a v_out
     * @param v_in vector to be pushed
     * @param v_in destination vector
     */
    template<typename T>
    void push_back (const vector <T> & v_in, vector <T> & v_out);

    /**
     * @brief return a string representation of the given vector. C style version
     * @param vec a vector
     * @return vector content as a string
     */
    template <typename T>
    std::string vec2str (const T * vec, const size_t sz) {
        std::ostringstream oss;
        oss << "[";
        for (size_t i = 0 ; i < sz ; ++i) {
            oss << vec[i];
            if (i < (sz - 1))
                oss << ", ";
        }
        oss << "]";
        return oss.str();
    };
    /**
     * @brief return a string representation of the given vector. STL version
     * @param vec a vector
     * @return vector content as a string
     */
    template <typename T>
    std::string vec2str ( const vector <T> & vec ) {
        return vec2str ((T*)(&vec[0]), vec.size());
    };
    /**
     * @brief print vector content
     * @param vec a vector as a C array
     * @param sz number of elements of the vector
     */
    template <typename T>
    void print_vec ( const T * vec, const size_t sz ) {
        cout << vec2str(vec, sz) << endl;
    };
    /**
     * @brief print vector content
     * @param vec a vector
     */
    template <typename T>
    void print_vec ( const vector <T> & vec ) {
        print_vec ((T*)(&vec[0]), vec.size());
    };
};

template<typename T>
void MyTools::push_back(const std::vector<T> & v_in, std::vector<T> & v_out) {
    if (! v_in.empty()) {
        v_out.reserve (v_out.size() + v_in.size());
        for (typename std::vector<T>::const_iterator it = v_in.begin() ; it != v_in.end() ; ++it)
            v_out.push_back(*it);
    }
};

template<typename T>
vector< T > MyTools::arange(const T & val_min, const T & val_max, const T & inc) {
    vector <T> v(0);
    arange (v, val_min, val_max, inc);
    return v;
}

template<typename T>
void MyTools::arange(vector< T > & v, const T & val_min, const T & val_max, const T & inc) {
    v.clear();
    size_t sz = (size_t)((val_max - val_min) / inc);
    v.reserve (sz);
    for (T val = val_min ; val < val_max ; val = (T)(val + inc))
        v.push_back (val);
}

#endif