#include "core/algo/GeoVolumeBuilder.hpp"

#include <gtest/gtest.h>

#include <cmath>
#include <cstdint>
#include <filesystem>
#include <fstream>
#include <iomanip>
#include <string>
#include <vector>

#include "data/store/ChunkedVolumeStore.hpp"

namespace fs = std::filesystem;
using namespace geopro::core;

// ---- principalAxisAngle 纯函数 ----
TEST(GeoVolumeBuilder, PcaFindsRoadDirection) {
  // 点沿 +Y（北向）排布 → 主轴角约 ±90°(±pi/2)。
  std::vector<double> xs, ys;
  for (int i = 0; i < 20; ++i) {
    xs.push_back(0.01 * i);  // 微小横向噪声
    ys.push_back(static_cast<double>(i));
  }
  const double ang = principalAxisAngle(xs, ys);
  EXPECT_NEAR(std::abs(ang), 3.14159265358979323846 / 2.0, 0.05);
}

TEST(GeoVolumeBuilder, PcaAlongEastIsZero) {
  std::vector<double> xs, ys;
  for (int i = 0; i < 20; ++i) {
    xs.push_back(static_cast<double>(i));
    ys.push_back(0.01 * i);
  }
  EXPECT_NEAR(principalAxisAngle(xs, ys), 0.0, 0.05);
}

TEST(GeoVolumeBuilder, PcaDegenerateReturnsZero) {
  EXPECT_DOUBLE_EQ(principalAxisAngle({}, {}), 0.0);
  EXPECT_DOUBLE_EQ(principalAxisAngle({1.0}, {1.0}), 0.0);
}

// ---- 合成小线建体 ----
namespace {

// 写一通道 .iprb + .iprh：值恒为 fixedVal（便于校验重叠均值）。
void writeChannel(const fs::path& iprb, int samples, int traces,
                  std::int16_t fixedVal) {
  fs::path iprh = fs::path(iprb).replace_extension(".iprh");
  std::ofstream h(iprh);
  h << "SAMPLES: " << samples << "\n";
  h << "LAST TRACE: " << (traces - 1) << "\n";
  h << "CHANNELS: 2\n";
  h << "TIMEWINDOW: 100.0\n";
  h << "SOIL VELOCITY: 100.0\n";  // m/µs → 1e8 m/s
  h << "DISTANCE INTERVAL: 0.05\n";
  h.close();

  std::ofstream b(iprb, std::ios::binary);
  for (int t = 0; t < traces; ++t)
    for (int s = 0; s < samples; ++s)
      b.write(reinterpret_cast<const char*>(&fixedVal), sizeof(fixedVal));
}

// 写一条南北直线的 .gps（lat 从 lat0 递增到 lat1，lon 固定）。
void writeGps(const fs::path& path, double lat0, double lat1, double lon,
              int pts) {
  std::ofstream f(path);
  for (int i = 0; i < pts; ++i) {
    const double frac = pts > 1 ? static_cast<double>(i) / (pts - 1) : 0.0;
    const double lat = lat0 + (lat1 - lat0) * frac;
    f << "2023-06-03\t00:00:00:000\t" << std::fixed
      << std::setprecision(10) << lat << "\tN\t" << lon << "\tE\t9.0\tM\t4\n";
  }
}

// 写两通道 .ord（横偏 -0.5 / +0.5，末列=1 有效）。
void writeOrd(const fs::path& path) {
  std::ofstream f(path);
  f << "0 -0.500000 -1.5 1\n";
  f << "1 0.500000 -1.5 1\n";
}

GeoLineInput makeLine(const fs::path& dir, const std::string& tag, double lat0,
                      double lat1, double lon, int traces, std::int16_t val) {
  const int samples = 8;
  writeChannel(dir / (tag + "_A01.iprb"), samples, traces, val);
  writeChannel(dir / (tag + "_A02.iprb"), samples, traces, val);
  writeOrd(dir / (tag + ".ord"));
  writeGps(dir / (tag + ".gps"), lat0, lat1, lon, traces);
  GeoLineInput in;
  in.iprb = {(dir / (tag + "_A01.iprb")).string(),
             (dir / (tag + "_A02.iprb")).string()};
  in.ord = (dir / (tag + ".ord")).string();
  in.gps = (dir / (tag + ".gps")).string();
  return in;
}

// 写一条 L 形弯轨迹的 .gps：先沿北(纬增,lon 固定)走半程，再沿东(经增,lat 固定)走半程。
void writeLGps(const fs::path& path, double lat0, double lon0, double dLatHalf,
               double dLonHalf, int pts) {
  std::ofstream f(path);
  const int half = pts / 2;
  for (int i = 0; i < pts; ++i) {
    double lat, lon;
    if (i <= half) {
      const double frac = half > 0 ? static_cast<double>(i) / half : 0.0;
      lat = lat0 + dLatHalf * frac;  // 北段：纬增
      lon = lon0;
    } else {
      const double frac =
          (pts - 1 - half) > 0
              ? static_cast<double>(i - half) / (pts - 1 - half)
              : 0.0;
      lat = lat0 + dLatHalf;            // 拐角后纬固定
      lon = lon0 + dLonHalf * frac;     // 东段：经增
    }
    f << "2023-06-03\t00:00:00:000\t" << std::fixed << std::setprecision(10)
      << lat << "\tN\t" << lon << "\tE\t9.0\tM\t4\n";
  }
}

GeoLineInput makeLLine(const fs::path& dir, const std::string& tag, double lat0,
                       double lon0, double dLatHalf, double dLonHalf, int traces,
                       std::int16_t val) {
  const int samples = 8;
  writeChannel(dir / (tag + "_A01.iprb"), samples, traces, val);
  writeChannel(dir / (tag + "_A02.iprb"), samples, traces, val);
  writeOrd(dir / (tag + ".ord"));
  writeLGps(dir / (tag + ".gps"), lat0, lon0, dLatHalf, dLonHalf, traces);
  GeoLineInput in;
  in.iprb = {(dir / (tag + "_A01.iprb")).string(),
             (dir / (tag + "_A02.iprb")).string()};
  in.ord = (dir / (tag + ".ord")).string();
  in.gps = (dir / (tag + ".gps")).string();
  return in;
}

}  // namespace

TEST(GeoVolumeBuilder, BuildsSyntheticLinesOverlapAveraged) {
  const fs::path tmp = fs::temp_directory_path() / "geopro_geovol_test";
  std::error_code ec;
  fs::remove_all(tmp, ec);
  fs::create_directories(tmp);

  // 两条同位置南北线（lat 30.200→30.201，~111m），值不同 → 重叠 cell 取均值。
  const double lat0 = 30.200, lat1 = 30.201, lon = 120.244;
  std::vector<GeoLineInput> lines = {
      makeLine(tmp, "synA_001", lat0, lat1, lon, /*traces=*/40, /*val=*/100),
      makeLine(tmp, "synB_002", lat0, lat1, lon, /*traces=*/40, /*val=*/300),
  };

  const std::string store = (tmp / "store").string();
  GeoGridSpec spec{/*cellXY=*/0.5, /*cellZ=*/0.1};
  GeoBuildResult r = buildGeoVolume(lines, spec, store, /*pyramidLevels=*/1);

  // 维度合理：约 111m 长 → nx 在数百量级；横路窄(~1m 阵列) → ny 小；nz>1。
  EXPECT_GT(r.nx, 50);
  EXPECT_GE(r.ny, 1);
  EXPECT_GT(r.nz, 1);
  EXPECT_GT(r.filled, 0);
  EXPECT_LE(r.filled, r.total);

  // store 可读，维度一致。
  geopro::data::ChunkedVolumeStore s(store);
  EXPECT_EQ(s.meta().nx, r.nx);
  EXPECT_EQ(s.meta().ny, r.ny);
  EXPECT_EQ(s.meta().nz, r.nz);
  EXPECT_EQ(s.levels(), 2);  // level0 + 1

  // 重叠均值：两线值 100/300，命中同 cell → 均值 200。扫所有块找非 blank 体素，
  // 其物理值应接近 200（量化误差内）。
  const auto& m = s.meta();
  bool foundNonBlank = false;
  double sampleVal = 0.0;
  for (int bz = 0; bz < s.bricksZ() && !foundNonBlank; ++bz)
    for (int by = 0; by < s.bricksY() && !foundNonBlank; ++by)
      for (int bx = 0; bx < s.bricksX() && !foundNonBlank; ++bx) {
        auto vox = s.readBrick(bx, by, bz);
        for (std::int16_t q : vox) {
          if (q != geopro::core::ScalarVolumeI16::kBlank) {
            sampleVal = m.quant.toPhys(q);
            foundNonBlank = true;
            break;
          }
        }
      }
  ASSERT_TRUE(foundNonBlank);
  EXPECT_NEAR(sampleVal, 200.0, 2.0);

  fs::remove_all(tmp, ec);
}

// L 形弯线：PCA 版把拐后段甩进大 Y'(虚高 ny、稀疏)，曲线版沿 X 拉直展开
// (ny 接近真实横向、填充更密)。断言曲线版 ny < PCA 版 ny 且填充率更高。
TEST(GeoVolumeBuilder, CurvilinearStraightensBendVsPca) {
  const fs::path tmp = fs::temp_directory_path() / "geopro_geovol_curv_test";
  std::error_code ec;
  fs::remove_all(tmp, ec);
  fs::create_directories(tmp);

  // L 形：北 ~111m（0.001° 纬）+ 东 ~96m（0.001° 经）。
  const double lat0 = 30.200, lon0 = 120.244;
  std::vector<GeoLineInput> lines = {makeLLine(
      tmp, "synL_001", lat0, lon0, /*dLatHalf=*/0.001, /*dLonHalf=*/0.001,
      /*traces=*/120, /*val=*/150)};

  GeoGridSpec spec{/*cellXY=*/0.5, /*cellZ=*/0.1};

  const std::string storePca = (tmp / "pca").string();
  GeoBuildResult rPca =
      buildGeoVolume(lines, spec, storePca, /*pyramidLevels=*/0,
                     /*curvilinear=*/false);

  const std::string storeCurv = (tmp / "curv").string();
  GeoBuildResult rCurv =
      buildGeoVolume(lines, spec, storeCurv, /*pyramidLevels=*/0,
                     /*curvilinear=*/true);

  const double fillPca =
      rPca.total > 0 ? static_cast<double>(rPca.filled) / rPca.total : 0.0;
  const double fillCurv =
      rCurv.total > 0 ? static_cast<double>(rCurv.filled) / rCurv.total : 0.0;

  // 曲线版把弯路拉直：横向维度 ny 明显小于 PCA 版。
  EXPECT_LT(rCurv.ny, rPca.ny);
  // 填充更密。
  EXPECT_GT(fillCurv, fillPca);
  // 曲线版仍建出有效体。
  EXPECT_GT(rCurv.filled, 0);
  EXPECT_GT(rCurv.nx, 50);

  fs::remove_all(tmp, ec);
}

// 距离权纯函数：距 cell 中心越近权越大，单调递减，中心处取最大 1。
TEST(GeoVolumeBuilder, DistanceWeightDecreasesWithDistance) {
  EXPECT_DOUBLE_EQ(distanceWeight(0.0), 1.0);          // 中心权最大
  EXPECT_LT(distanceWeight(1.0), distanceWeight(0.5)); // 单调递减
  EXPECT_LT(distanceWeight(2.0), distanceWeight(1.0));
  EXPECT_GT(distanceWeight(0.5), 0.0);
}

// 加权均值偏向近点：同 cell 两样本，近点(dist=0.1)值 100、远点(dist=0.9)值 300。
// 等权均值 = 200；距离加权均值应明显 < 200（偏向近点的 100）。
TEST(GeoVolumeBuilder, WeightedMeanFavorsNearerPoint) {
  const double vNear = 100.0, vFar = 300.0;
  const double wNear = distanceWeight(0.1);
  const double wFar = distanceWeight(0.9);
  const double weighted =
      (wNear * vNear + wFar * vFar) / (wNear + wFar);
  const double equal = 0.5 * (vNear + vFar);  // 200
  EXPECT_LT(weighted, equal);       // 偏离等权
  EXPECT_LT(weighted - vNear, vFar - weighted);  // 更靠近近点值
}