//     Copyright (c) 2012 Vadym Kliuchnikov sqct(dot)software(at)gmail(dot)com, Dmitri Maslov, Michele Mosca
//
//     This file is part of SQCT.
// 
//     SQCT is free software: you can redistribute it and/or modify
//     it under the terms of the GNU Lesser General Public License as published by
//     the Free Software Foundation, either version 3 of the License, or
//     (at your option) any later version.
// 
//     SQCT 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 Lesser General Public License for more details.
// 
//     You should have received a copy of the GNU Lesser General Public License
//     along with SQCT.  If not, see <http://www.gnu.org/licenses/>.
// 

#include <fstream>
#include <algorithm>
#include <cassert>
#include <iostream>
#include <iterator>
#include "epsilonnet.h"
#include "output.h"

#include <boost/range.hpp>
using boost::make_iterator_range;

using namespace std;

/// \brief Header of the file with epsilon net
struct epsilonnetHeader
{
    size_t nodes_count;
};

std::ostream& operator<<(std::ostream& out, const enetNode& node )
{
    out << "enetNode[" << node.ipxx << "," << node.ipQxx << ","
        << node.num_offset << "," << node.compl_offset << "]";
    return out;
}

bool epsilonnet::loadFromFile(const char* filename)
{
    epsilonnetHeader eh;
    ifstream ifs( filename, ios_base::binary );
    if( !ifs )
        return false;
    ifs.read( (char*) &eh,sizeof(eh));
    nodes.clear();
    nodes.resize( eh.nodes_count );
    ifs.read( (char*) &nodes[0], nodes.size() * sizeof(enetNode) );
    numbers.clear();
    numbers.resize( nodes.back().num_offset );
    ifs.read( (char*) &numbers[0], numbers.size() * sizeof(ri) );
    ifs.close();
    denominator_exponent = denominatorExponent2();
    return true;
}

epsilonnet::epsilonnet()
{
    enetNode nd={0,0,0,0};
    nodes.push_back(nd);
}

size_t epsilonnet::nodesCount(const char* filename) const
{
    epsilonnetHeader eh;
    ifstream ifs( filename, ios_base::binary );
    ifs.read( (char*) &eh,sizeof(eh));
    ifs.close();
    return eh.nodes_count;
}

void epsilonnet::saveToFile  (const char* filename) const
{
    if( numbers.size() == 0 )
        return;

    epsilonnetHeader eh;
    eh.nodes_count = nodes.size();
    ofstream ofs( filename, ios_base::binary );
    ofs.write( (const char*) &eh,sizeof(eh));

    assert( nodes.back().ipxx == 0 );
    assert( nodes.back().ipQxx == 0 );
    assert( nodes.back().num_offset == numbers.size() );

    ofs.write( (const char*) &nodes[0], nodes.size() * sizeof(enetNode) );
    ofs.write( (const char*) &numbers[0], numbers.size() * sizeof(ri) );

    ofs.close();
}

/// \brief Comparator based on pointers data
template< class T>
struct eq_ptr
{
    bool operator() ( const T* a , const T* b )
    {
        return *a == *b;
    }
};

/// \brief Modified stub of back_insert_iterator from http://www.cplusplus.com/reference/std/iterator/back_insert_iterator/
template <class Container>
  class back_insert_iterator_ptr :
    public iterator<output_iterator_tag,void,void,void,void>
{
protected:
  Container* container;

public:
  typedef Container container_type;
  explicit back_insert_iterator_ptr (Container& x) : container(&x) {}
  back_insert_iterator_ptr<Container>& operator= (const typename Container::value_type* value)
    { container->push_back(*value); return *this; }
  back_insert_iterator_ptr<Container>& operator* ()
    { return *this; }
  back_insert_iterator_ptr<Container>& operator++ ()
    { return *this; }
  back_insert_iterator_ptr<Container> operator++ (int)
    { return *this; }
};

void epsilonnet::addNode( const pair< ip_type,ip_type>& ip, const nodeRangesPtr &ranges)
{
    nodes.back().ipxx = ip.first;
    nodes.back().ipQxx = ip.second;

    eq_ptr< ri > eq_ri;
    back_insert_iterator_ptr< vector<ri> > biit( numbers );

    int size = numbers.size();
    unique_copy( ranges.nums_begin, ranges.nums_end , biit , eq_ri );
    nodes.back().compl_offset = numbers.size() - size;
    unique_copy( ranges.nums_compl_begin, ranges.nums_compl_end , biit , eq_ri );

    enetNode en = {0,0,numbers.size(),0};
    nodes.push_back(en);
}

void epsilonnet::addNode(epsilonnet::ip_type ipxx, epsilonnet::ip_type ipQxx, const nodeRanges &ranges)
{
    nodes.back().ipxx = ipxx;
    nodes.back().ipQxx = ipQxx;

    auto bi = back_inserter( numbers );
    int size = numbers.size();
    unique_copy( ranges.nums_begin, ranges.nums_end ,bi );
    nodes.back().compl_offset = numbers.size() - size;
    unique_copy( ranges.nums_compl_begin, ranges.nums_compl_end , bi );

    enetNode en = {0,0,numbers.size(),0};
    nodes.push_back(en);
}

void epsilonnet::getNode(size_t node_id, nodeRanges &ranges) const
{
    const enetNode& node = nodes[ node_id ];
    ranges.nums_begin = numbers.begin() + node.num_offset;
    ranges.nums_end = numbers.begin() + complOffset( node_id );
    ranges.nums_compl_begin = ranges.nums_end;
    ranges.nums_compl_end = numbers.begin() + nodes[ node_id + 1 ].num_offset;
}

size_t epsilonnet::complOffset(size_t nodeId) const
{
    const enetNode& node = nodes[ nodeId ];
    return node.num_offset + node.compl_offset;
}

int epsilonnet::denominatorExponent2() const
{
    size_t offset = complOffset(0);
    auto denom = numbers[offset].ipxx() + numbers[0].ipxx();
    int val = ri::gde2( denom );
    // overflow check
    decltype( denom ) d = 1;
    d <<= val;
    assert( d == denom && "Sum of ipxx for number and complementary number is not a power of 2." );
    //denominator_exponent = val;
    return val;
}

int epsilonnet::sde() const
{
    int res = ::sde( 2 * denominatorExponent2(), numbers[0].abs2().gde() );
    if( res == std::numeric_limits<int>::max() ) res = 0;
    return res;
}

double epsilonnet::findExhaustiveApproximation(const epsilonnet::vector2double &vec, epsilonnet::vi &result) const
{
    //denominator_exponent = denominatorExponent2();
    vi current;
    double current_dist = 1.0;
    double best_dist = 1.0;
    vector2double canonical_vec = vec;
    bool conj_1 = false, conj_2 = false;
    int w_pow1 = 0, w_pow2 = 0;
    canonical_vec.first = canonical( vec.first, w_pow1, conj_1 );
    canonical_vec.second = canonical( vec.second, w_pow2, conj_2 );

    for( int i = 0; i < nodes.size(); ++i )
    {
        current_dist = findExhaustiveApproximation( canonical_vec, current, i );
        if( best_dist > current_dist )
        {
            best_dist = current_dist;
            result  = current;
        }
    }

//    cout << result.de << ":" << endl;
//    cout <<  canonical_vec.first << result.d[0].toComplex( result.de ) << endl;
//    cout <<  canonical_vec.second << result.d[1].toComplex( result.de ) << endl;

    if( conj_1 ) result.d[0].conjugate_eq();
    if( conj_2 ) result.d[1].conjugate_eq();
    result.d[0].mul_eq_w( w_pow1 );
    result.d[1].mul_eq_w( w_pow2 );
    result.de = denominator_exponent;

    return best_dist;
}

typedef  epsilonnet::vector2double vd;

/// \brief Inline Euclidean distance computation
inline double dist_squared( const vd& a, const vd& b )
{
    double d1 = a.first.real() - b.first.real();
    double d2 = a.second.real() - b.second.real();
    double d3 = a.first.imag() - b.first.imag();
    double d4 = a.second.imag() - b.second.imag();
    return d1*d1 +d2*d2 +d3*d3 + d4*d4;
}

/// \brief Inline Euclidean distance computation
inline double dist_squared( const vd& a, double fre,  double fim , double sre ,  double sim )
{
    double d1 = a.first.real() - fre;
    double d2 = a.second.real() - sre;
    double d3 = a.first.imag() - fim;
    double d4 = a.second.imag() - sim;
    return d1*d1 +d2*d2 +d3*d3 + d4*d4;
}

/// \todo Remove code duplication by changing one of epsilonnet::findExhaustiveApproximation
double epsilonnet::findExhaustiveApproximation(const epsilonnet::vector2double &vec, epsilonnet::vi &result, int node_id) const
{
    double best = 3.0;
    nodeRanges nr;
    getNode( node_id, nr );
    auto ic = nr.nums_begin;
    auto jc = nr.nums_compl_begin;
    bool tw = true;
    for( auto i = nr.nums_begin; i != nr.nums_end ;++i )
    {
        double iim, ire;
        i->toComplex( denominator_exponent, ire, iim );
        for( auto j = nr.nums_compl_begin; j != nr.nums_compl_end; ++j )
        {
            double jim, jre;
            j->toComplex( denominator_exponent, jre, jim  );
            double d1 = dist_squared( vec, ire, iim, jre, jim );
            double d2 = dist_squared( vec, jre, jim, ire, iim );

            if( d1 < best )
            {
                best = d1;
                ic = i; jc = j;
                tw = false;
            }

            if( d2 < best )
            {
                best = d2;
                ic = i; jc = j;
                tw = true;
            }
        }
    }

    if( tw )
    {
        result.d[0] = *jc;
        result.d[1] = *ic;
        //result.de = denominator_exponent; // -avoid extra operations
    }
    else
    {
        result.d[0] = *ic;
        result.d[1] = *jc;
        //result.de = denominator_exponent; // -avoid extra operations
    }

    return best;
}

double epsilonnet::findExhaustiveApproximation(const epsilonnet::vector2double &vec, epsilonnet::vi &result, int node_id, double &bdist) const
{
    nodeRanges nr;
    getNode( node_id, nr );
    bool found = false;
    auto ic = nr.nums_begin;
    auto jc = nr.nums_compl_begin;
    bool tw = true;
    for( auto i = nr.nums_begin; i != nr.nums_end ;++i )
    {
        double iim, ire;
        i->toComplex( denominator_exponent, ire, iim );
        for( auto j = nr.nums_compl_begin; j != nr.nums_compl_end; ++j )
        {
            double jim, jre;
            j->toComplex( denominator_exponent, jre, jim  );
            double d1 = dist_squared( vec, ire, iim, jre, jim );
            double d2 = dist_squared( vec, jre, jim, ire, iim );

            if( d1 < bdist )
            {
                found = true;
                bdist = d1;
                ic = i; jc = j;
                tw = false;
            }

            if( d2 < bdist )
            {
                found = true;
                bdist = d2;
                ic = i; jc = j;
                tw = true;
            }
        }
    }

    if( found )
    {
        if( tw )
        {
            result.d[0] = *jc;
            result.d[1] = *ic;
            result.de = denominator_exponent;
        }
        else
        {
            result.d[0] = *ic;
            result.d[1] = *jc;
            result.de = denominator_exponent;
        }
    }

    return bdist;
}

//// Code that was used to iterate through epsilon net. Does not required now

//private:
//    /// \brief Executes functor for all complex values in epsilon net
//    template< class T>
//    void for_all_complex( T& functor ) const;
//    /// \brief Executes functor for all complex values in epsilon net within node
//    template< class T>
//    void for_all_complex( T& functor, int node_id ) const;
//    /// \brief Executes functor for all complex values in epsilon net within node with fixed first element
//    template< class T>
//    void for_all_complex( T& functor, nodeRanges& nr, ri& a ) const;
//    /// \brief Executes functor for a, b and elements that are of the form \f$ a\omega^k, a\omega^k \f$
//    template< class T>
//    void for_all_complex( T& functor, std::complex<double>& a, ri& b ) const;

//template< class T>
//void epsilonnet::for_all_complex( T& functor ) const
//{
//    denominator_exponent = denominatorExponent2();
//    for( int i = 0; i < nodes.size() - 1 ; ++i )
//        for_all_complex( functor, i );
//}

//template< class T>
//void epsilonnet::for_all_complex( T& functor, int node_id ) const
//{
//    nodeRanges nr;
//    getNode( node_id, nr );
//    for( auto i : make_iterator_range( nr.nums_begin, nr.nums_end) )
//    {
//        ri a(i);
//        if( i.is_im_eq0() )
//        {
//            for_all_complex( functor, nr, a );
//        }
//        else
//        {
//            for_all_complex( functor, nr, a );
//            a.conjugate_eq();
//            for_all_complex( functor, nr, a );
//        }
//    }
//}

//template< class T>
//void epsilonnet::for_all_complex( T& functor, nodeRanges& nr, ri& a ) const
//{
//    for( int i = 0; i < 8 ; ++i )
//    {
//        a.mul_eq_w();
//        auto first = a.toComplex( denominator_exponent );
//        for( auto j : make_iterator_range(nr.nums_compl_begin, nr.nums_compl_end) )
//        {
//            ri b(j);
//            if( j.is_im_eq0() )
//            {
//                for_all_complex( functor, first, b );
//            }
//            else
//            {
//                for_all_complex( functor, first, b );
//                b.conjugate_eq();
//                for_all_complex( functor, first, b );
//            }
//        }
//    }
//}

//template< class T>
//void epsilonnet::for_all_complex( T& functor, std::complex<double>& first, ri& b ) const
//{
//    int max_ph = 1;
//    if( all_phases ) max_ph = 8;

//    for( int i = 0; i < max_ph ; ++i )
//    {
//        b.mul_eq_w();
//        auto second = b.toComplex( denominator_exponent );
//        functor( first, second );
//        functor( second, first );
//    }
//}