//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      Base/Axis/Scale.cpp
//! @brief     Implements interface Scale.
//!
//! @homepage  http://www.bornagainproject.org
//! @license   GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
//  ************************************************************************************************

#include "Base/Axis/Scale.h"
#include "Base/Const/Units.h"
#include "Base/Math/Numeric.h"
#include "Base/Util/Assert.h"
#include "Base/Util/StringUtil.h"
#include <cmath>
#include <iomanip>
#include <iostream>
#include <numbers>
#include <stdexcept>

using Numeric::almostEqual;
using std::numbers::pi;

const int ulp = 7; // required precision for isEqui... tests, in units of epsilon

Scale::Scale(const Coordinate& coord, const std::vector<Bin1D>& bins)
    : m_bins(bins)
    , m_coord(std::make_unique<Coordinate>(coord))
{
    if (size() == 0)
        throw std::runtime_error("Axis \"" + coord.name() + "\" has no bins");
    for (size_t i = 0; i < size() - 1; ++i) {
        if (bin(i).max() > bin(i + 1).min() && bin(i + 1).max() > bin(i).min())
            throw std::runtime_error("Axis \"" + coord.name() + "\" has overlapping bins");
        if (bin(i) == bin(i + 1))
            throw std::runtime_error("Axis \"" + coord.name() + "\" has repeating bins");
    }
    if (isScan()) {
        for (const Bin1D& b : m_bins)
            if (b.binSize() != 0)
                throw std::runtime_error("Finite bin(s) in scan");
    } else {
        for (const Bin1D& b : m_bins)
            if (b.binSize() == 0)
                throw std::runtime_error("Empty bin(s) in sweep");
    }
}

Scale::Scale(const Scale& other)
    : m_bins(other.m_bins)
{
    ASSERT(other.m_coord);
    m_coord = std::make_unique<Coordinate>(*other.m_coord);
}

Scale* Scale::clone() const
{
    return new Scale(*this);
}

std::string Scale::axisLabel() const
{
    ASSERT(m_coord);
    return m_coord->label();
}

std::string Scale::coordName() const
{
    return Coordinate(axisLabel()).name();
}

std::string Scale::unit() const
{
    return Coordinate(axisLabel()).unit();
}

size_t Scale::size() const
{
    return m_bins.size();
}

double Scale::min() const
{
    return m_bins.front().min();
}

double Scale::max() const
{
    return m_bins.back().max();
}

std::pair<double, double> Scale::bounds() const
{
    return {min(), max()};
}

bool Scale::rangeComprises(double value) const
{
    return value >= min() && value < max();
}

double Scale::span() const
{
    return max() - min();
}

const Bin1D& Scale::bin(size_t i) const
{
    return m_bins.at(i);
}

double Scale::binCenter(size_t i) const
{
    return bin(i).center();
}

std::vector<double> Scale::binCenters() const
{
    std::vector<double> result;
    result.reserve(m_bins.size());
    for (const Bin1D& b : m_bins)
        result.emplace_back(b.center());
    return result;
}

size_t Scale::closestIndex(double value) const
{
    for (size_t i = 0; i < size() - 1; ++i)
        if (value < (bin(i).max() + bin(i + 1).min()) / 2)
            return i;
    return size() - 1;
}

bool Scale::isEquiDivision() const
{
    const size_t N = size();
    for (size_t i = 0; i < N; ++i) {
        const Bin1D& b = bin(i);
        // exactly replicate the computation of bin bounds in the EquiDivision factory function
        if (!almostEqual(b.min(), (N - i) * (min() / N) + i * (max() / N), ulp)
            || !almostEqual(b.max(), (N - i - 1) * (min() / N) + (i + 1) * (max() / N), ulp))
            return false;
    }
    return true;
}

bool Scale::isEquiScan() const
{
    const size_t N = size();
    ASSERT(N);
    if (N == 1)
        return !bin(0).binSize();
    for (size_t i = 0; i < N; ++i) {
        const Bin1D& b = bin(i);
        if (b.binSize())
            return false;
        // exactly replicate the computation of bin bounds in the EquiScan factory function
        if (!almostEqual(b.min(), (N - 1 - i) * (min() / (N - 1)) + i * (max() / (N - 1)), ulp))
            return false;
    }
    return true;
}

bool Scale::isScan() const
{
    for (const Bin1D& b : bins())
        if (b.binSize())
            return false;
    return true;
}

Scale Scale::clipped(double lower, double upper) const
{
    if (lower > upper)
        throw std::runtime_error("Scale::clipped called with invalid bounds (lower > upper)");
    std::vector<Bin1D> out_bins;
    const bool is_scan = isScan();
    for (const Bin1D& b : m_bins)
        if (auto bc = b.clipped_or_nil(lower, upper))
            if (is_scan || bc.value().binSize() > 0)
                out_bins.emplace_back(bc.value());
    return {m_coord->label(), out_bins};
}

Scale Scale::clipped(std::pair<double, double> bounds) const
{
    return clipped(bounds.first, bounds.second);
}

bool Scale::operator==(const Scale& other) const
{
    return axisLabel() == other.axisLabel() && m_bins == other.m_bins;
}

std::ostream& operator<<(std::ostream& ostr, const Scale& ax)
{
    size_t N = ax.size();
    ASSERT(N > 0);

    ostr << std::setprecision(15);

    if (ax.isScan()) {
        ostr << "ListScan(\"" << ax.axisLabel() << "\", [";
        for (double v : ax.binCenters())
            ostr << v << ",";
        ostr << "])";
        return ostr;
    }

    if (ax.isEquiDivision()) {
        ostr << "EquiDivision(\"" << ax.axisLabel() << "\", " << ax.size() << ", " << ax.min()
             << ", " << ax.max() << ")";
        return ostr;
    }

    ostr << "GenericScale(\"" << ax.axisLabel() << "\", [";
    for (const Bin1D& b : ax.bins())
        ostr << b.min() << "," << b.max() << ",";
    ostr << "])";
    return ostr;
}

Scale Scale::plottableScale() const
{
    ASSERT(m_coord);
    if (m_coord->unit() == "rad")
        return transformedScale(Coordinate(m_coord->name(), "deg"), Units::rad2deg);
    return {m_coord->label(), m_bins};
}

Scale Scale::transformedScale(const Coordinate& coord, const trafo_t& axTrafo) const
{
    std::vector<Bin1D> outvector;
    for (const Bin1D& b : m_bins) {
        double bmi = axTrafo(b.min());
        double bma = axTrafo(b.max());
        outvector.emplace_back(Bin1D::FromTo(bmi, bma));
    }
    return {coord, outvector};
}

Scale Scale::reversedScale() const
{
    std::vector<Bin1D> outvector;
    for (const Bin1D& b : m_bins)
        outvector.emplace(outvector.begin(), Bin1D::FromTo(b.max(), b.min()));
    return {*m_coord, outvector};
}

Scale Scale::alpha_f_Scale(double lambda, double alpha_i) const
{
    if (m_coord->unit() == "1/nm")
        return transformedScale("alpha_f (rad)", [lambda, alpha_i](double qz) {
            return std::asin(qz * lambda / 2 / pi - std::sin(alpha_i));
        });
    return {m_coord->label(), m_bins};
}

Scale Scale::phi_f_Scale(double lambda) const
{
    if (m_coord->unit() == "1/nm")
        return transformedScale("phi_f (rad)",
                                [lambda](double qy) { return std::asin(qy * lambda / 2 / pi); });
    return {m_coord->label(), m_bins};
}

Scale Scale::qz_Scale(double lambda, double alpha_i) const
{
    if (m_coord->unit() == "rad")
        return transformedScale("q_z (1/nm)", [lambda, alpha_i](double alpha_f) {
            return 2 * pi / lambda * (std::sin(alpha_i) + std::sin(alpha_f));
        });
    return {m_coord->label(), m_bins};
}

Scale Scale::qy_Scale(double lambda) const
{
    if (m_coord->unit() == "rad")
        return transformedScale(
            "q_y (1/nm)", [lambda](double phi_f) { return 2 * pi / lambda * std::sin(phi_f); });
    return {m_coord->label(), m_bins};
}