#include <iostream>
#include <vector>
#include <cmath>
#include <algorithm>
#include <iomanip>
#include <fstream>
#include <string>
#include <sstream>
#include <Eigen/Dense> 
#include <omp.h> 

#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

namespace NItem
{
    // --- 基础结构 ---
    struct TrendParam {
        double angle; double ratio;
        double cx, cy; double weight;
    };

    struct GridResult {
        int nx, ny;
        double xMin, xMax, yMin, yMax, step;
        std::vector<double> gx, gy, zi;
    };

    class GeoEngine
    {
    public:
        void initGridStrict(const std::vector<double>& sx, const std::vector<double>& sy, double step, GridResult& grid) {
            auto x_mm = std::minmax_element(sx.begin(), sx.end());
            auto y_mm = std::minmax_element(sy.begin(), sy.end());
            double dx = *x_mm.second - *x_mm.first;
            double dy = *y_mm.second - *y_mm.first;

            // 10% 边界扩展
            grid.xMin = *x_mm.first - dx * 0.1;
            grid.yMin = *y_mm.first - dy * 0.1;
            double targetMaxX = *x_mm.second + dx * 0.1;
            double targetMaxY = *y_mm.second + dy * 0.1;

            grid.step = step;
            grid.nx = 0;
            for (double x = grid.xMin; x < targetMaxX + step - 1e-8; x += step) grid.nx++;
            grid.ny = 0;
            for (double y = grid.yMin; y < targetMaxY + step - 1e-8; y += step) grid.ny++;

            // 计算实际导出边界
            grid.xMax = grid.xMin + (grid.nx - 1) * step;
            grid.yMax = grid.yMin + (grid.ny - 1) * step;

            size_t total = (size_t)grid.nx * grid.ny;
            grid.gx.assign(total, 0.0);
            grid.gy.assign(total, 0.0);
            grid.zi.assign(total, 0.0);

            // 生成扁平化网格：Y 为行 (从小到大)，X 为列 (从小到大)
            for (int r = 0; r < grid.ny; ++r) {
                double cur_y = grid.yMin + r * step;
                for (int c = 0; c < grid.nx; ++c) {
                    size_t idx = (size_t)r * grid.nx + c;
                    grid.gx[idx] = grid.xMin + c * step;
                    grid.gy[idx] = cur_y;
                }
            }
        }

        /**
         * 2. 最小曲率插值 (TPS) + 并行处理
         */
        void calculateSpline(const std::vector<double>& sx, const std::vector<double>& sy, const std::vector<double>& sz,
            const std::vector<TrendParam>& trends, int fusionMode, GridResult& grid) {
            int n = sx.size();
            if (n == 0 || trends.empty()) return;

            std::vector<std::vector<double>> trendResults(trends.size());

            for (size_t t_idx = 0; t_idx < trends.size(); ++t_idx) {
                const auto& t = trends[t_idx];
                double rad = (90.0 - t.angle) * M_PI / 180.0;
                double cosA = std::cos(rad), sinA = std::sin(rad);

                std::vector<double> tx(n), ty(n);
                for (int i = 0; i < n; ++i) {
                    double dx = sx[i] - t.cx;
                    double dy = sy[i] - t.cy;
                    tx[i] = (dx * cosA + dy * sinA) / t.ratio;
                    ty[i] = -dx * sinA + dy * cosA;
                }

                Eigen::MatrixXd A(n, n);
                Eigen::VectorXd b(n);
                for (int i = 0; i < n; ++i) {
                    b(i) = sz[i];
                    for (int j = 0; j < n; ++j) {
                        double r = std::sqrt(std::pow(tx[i] - tx[j], 2) + std::pow(ty[i] - ty[j], 2));
                        A(i, j) = (r < 1e-12) ? 0 : r * r * std::log(r);
                    }
                }
                A.diagonal().array() += 1e-12; // 稳定性因子
                Eigen::VectorXd weights = A.fullPivLu().solve(b);

                std::vector<double> cur_zi(grid.gx.size());
#pragma omp parallel for schedule(dynamic)
                for (long long i = 0; i < (long long)grid.gx.size(); ++i) {
                    double dx = grid.gx[i] - t.cx;
                    double dy = grid.gy[i] - t.cy;
                    double rx = (dx * cosA + dy * sinA) / t.ratio;
                    double ry = -dx * sinA + dy * cosA;
                    double val = 0;
                    for (int j = 0; j < n; ++j) {
                        double r = std::sqrt(std::pow(rx - tx[j], 2) + std::pow(ry - ty[j], 2));
                        val += weights(j) * ((r < 1e-12) ? 0 : r * r * std::log(r));
                    }
                    cur_zi[i] = val;
                }
                trendResults[t_idx] = cur_zi;
            }

            grid.zi.assign(grid.gx.size(), 0.0);
#pragma omp parallel for schedule(static)
            for (long long i = 0; i < (long long)grid.gx.size(); ++i) {
                if (fusionMode == 1) {
                    double maxVal = -1e18;
                    for (const auto& res : trendResults) if (res[i] > maxVal) maxVal = res[i];
                    grid.zi[i] = maxVal;
                }
                else {
                    double sumV = 0, sumW = 0;
                    for (size_t t = 0; t < trends.size(); ++t) {
                        sumV += trendResults[t][i] * trends[t].weight;
                        sumW += trends[t].weight;
                    }
                    grid.zi[i] = (sumW > 1e-12) ? sumV / sumW : 0;
                }
            }
        }

        /**
         * 3. 导出 GRD (标准 DSAA)
         * 修复点：输出 yMin yMax，确保步长为正值 (+20)
         */
        void exportToGRD(const std::string& filename, const GridResult& grid) {
            std::ofstream f(filename);
            if (!f.is_open()) return;

            auto mm = std::minmax_element(grid.zi.begin(), grid.zi.end());

            f << "DSAA\n";
            f << grid.nx << " " << grid.ny << "\n";
            f << std::fixed << std::setprecision(6) << grid.xMin << " " << grid.xMax << "\n";
            f << grid.yMin << " " << grid.yMax << "\n"; // 修改为 yMin yMax
            f << std::scientific << std::setprecision(6) << *mm.first << " " << *mm.second << "\n";

            // 按 Surfer 10个数值一行的标准写入
            for (size_t i = 0; i < grid.zi.size(); ++i) {
                f << grid.zi[i];
                if ((i + 1) % 10 == 0 || (i + 1) == grid.zi.size()) f << "\n";
                else f << " ";
            }
        }
    };

    bool loadSandstoneCSV(const std::string& fname, std::vector<double>& x, std::vector<double>& y, std::vector<double>& z) {
        std::ifstream f(fname);
        if (!f.is_open()) return false;
        std::string line; std::getline(f, line);
        while (std::getline(f, line)) {
            if (line.empty()) continue;
            std::stringstream ss(line);
            std::string item; std::vector<std::string> row;
            while (std::getline(ss, item, ',')) row.push_back(item);
            if (row.size() >= 4) {
                try {
                    x.push_back(std::stod(row[1])); y.push_back(std::stod(row[2])); z.push_back(std::stod(row[3]));
                }
                catch (...) { continue; }
            }
        }
        return !x.empty();
    }
}
//
//int main() {
//    GeoEngine engine;
//    GridResult grid;
//    std::vector<double> sx, sy, sz;
//
//    // 1. 读取数据
//    if (!loadSandstoneCSV("G:/Desktop/333.csv", sx, sy, sz)) return -1;
//
//    // 2. 初始化网格 (步长 20)
//    engine.initGridStrict(sx, sy, 20.0, grid);
//
//    // 3. 确定旋转中心 (与 Python fixed 版对齐)
//    auto x_mm = std::minmax_element(sx.begin(), sx.end());
//    auto y_mm = std::minmax_element(sy.begin(), sy.end());
//    double midX = (*x_mm.first + *x_mm.second) / 2.0;
//    double midY = (*y_mm.first + *y_mm.second) / 2.0;
//
//    // 4. 计算与导出
//    TrendParam p1 = { 0.0, 0.0, midX, midY, 1.0 };
//    engine.calculateSpline(sx, sy, sz, { p1 }, 0, grid);
//    engine.exportToGRD("G:/Desktop/result_cpp_fixed.grd", grid);
//}