// Copyright (2019-2023) Paul Scherrer Institute

#include <catch2/catch.hpp>

#include "../image_analysis/RadialIntegrationProfile.h"

TEST_CASE("RadialIntegrationMapping_Constructor","[RadialIntegration]") {
    DiffractionExperiment x;

    std::unique_ptr<RadialIntegrationMapping> radial;

    x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(5);
    REQUIRE_NOTHROW(radial = std::make_unique<RadialIntegrationMapping>(x));
}

TEST_CASE("RadialIntegrationMapping_GetBinNumber","[RadialIntegration]") {
    DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36));
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(4);
    RadialIntegrationMapping mapping(x);
    REQUIRE(mapping.GetBinNumber() == 39);
}

TEST_CASE("RadialIntegrationMapping_GetBinNumber_DetectorLimit","[RadialIntegration]") {
    DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36));
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);

    x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(9.9);
    RadialIntegrationMapping mapping(x);
    REQUIRE(mapping.GetBinNumber() < 80);
}

TEST_CASE("RadialIntegrationMapping_GetBinToQ","[RadialIntegration]") {
    DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36));
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(4);

    RadialIntegrationMapping mapping(x);
    auto bin_to_q = mapping.GetBinToQ();

    REQUIRE(bin_to_q[0] == Approx(0.15));
    REQUIRE(bin_to_q[1] == Approx(0.25));
    REQUIRE(bin_to_q[15] == Approx(1.65));
    REQUIRE(bin_to_q[38] == Approx(3.95));
}

TEST_CASE("RadialIntegrationMapping_QToBin","[RadialIntegration]") {
    DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36));
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(4);

    RadialIntegrationMapping mapping(x);

    REQUIRE(mapping.QToBin(0.0) == 0);
    REQUIRE(std::floor(mapping.QToBin(0.200001)) == 1);
    REQUIRE(mapping.QToBin(0.6) == Approx(5));
    REQUIRE(mapping.QToBin(50.0) == Approx(38));
}

TEST_CASE("RadialIntegrationProfile","[RadialIntegration]") {
    DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36));
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(4);

    RadialIntegrationMapping mapping(x);

    RadialIntegrationProfile profile(mapping, x);

    std::vector<int64_t> sum(mapping.GetBinNumber());
    std::vector<uint64_t> count(mapping.GetBinNumber());

    for (int i = 0; i < mapping.GetBinNumber(); i++) {
        sum[i] = (i * i * 4) * (1LU<<24);
        count[i] = i;
    }
    REQUIRE_NOTHROW(profile.Add(sum, count));
    REQUIRE_NOTHROW(profile.Add(sum, count));

    std::vector<int64_t> sum_wr(mapping.GetBinNumber() - 1);
    REQUIRE_THROWS(profile.Add(sum_wr, count));

    Plot plot = profile.GetPlot();

    REQUIRE(plot.x.size() == mapping.GetBinNumber());
    REQUIRE(plot.y.size() == mapping.GetBinNumber());
    for (int i = 0; i < mapping.GetBinNumber(); i++) {
        REQUIRE(plot.x[i] == Approx(mapping.GetBinToQ()[i]));
        REQUIRE(plot.y[i] == Approx(i * 4));
    }
}

TEST_CASE("RadialIntegrationProfile_operatorAdd","[RadialIntegration]") {
    DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36));
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(4);

    RadialIntegrationMapping mapping(x);

    RadialIntegrationProfile profile0(mapping, x), profile1(mapping, x);

    std::vector<int64_t> sum(mapping.GetBinNumber());
    std::vector<uint64_t> count(mapping.GetBinNumber());

    for (int i = 0; i < mapping.GetBinNumber(); i++) {
        sum[i] = (i * i * 4) * (1LU<<24);
        count[i] = i;
    }
    REQUIRE_NOTHROW(profile0.Add(sum, count));
    REQUIRE_NOTHROW(profile1 += profile0);

    auto plot = profile1.GetPlot();

    REQUIRE(plot.x.size() == mapping.GetBinNumber());
    REQUIRE(plot.y.size() == mapping.GetBinNumber());
    for (int i = 0; i < mapping.GetBinNumber(); i++) {
        REQUIRE(plot.x[i] == Approx(mapping.GetBinToQ()[i]));
        REQUIRE(plot.y[i] == Approx(i * 4));
    }
}

TEST_CASE("RadialIntegrationProfile_GetMeanValueOfBins","[RadialIntegration]") {
    DiffractionExperiment x(DetectorGeometry(8, 2, 8, 36));
    x.DetectorDistance_mm(50).BeamX_pxl(1000).BeamY_pxl(1000);
    x.QSpacingForRadialInt_recipA(0.1).LowQForRadialInt_recipA(0.1).HighQForRadialInt_recipA(4);

    RadialIntegrationMapping mapping(x);

    RadialIntegrationProfile profile(mapping, x);

    std::vector<int64_t> sum(mapping.GetBinNumber());
    std::vector<uint64_t> count(mapping.GetBinNumber());

    for (int i = 0; i < mapping.GetBinNumber(); i++) {
        sum[i] = (i * i * 4) * (1LU<<24);
        count[i] = i;
    }
    REQUIRE_NOTHROW(profile.Add(sum, count));

    REQUIRE(profile.GetMeanValueOfBins(0,2) == Approx ((sum[0] + sum[1] + sum[2]) /(1LU<<24) / double(count[0] + count[1] + count[2])));
    REQUIRE(profile.GetMeanValueOfBins(5,7) == Approx ((sum[5] + sum[6] + sum[7]) /(1LU<<24) / double (count[5] + count[6] + count[7])));
}