#include "GaussPolyCoe.h"

using namespace std;

vector<double> GaussPolyCoe::GaussianElimination(vector<vector<double>> A, vector<double> b)
{
    try {
        int n = A.size();
        for (int i = 0; i < n; ++i) {
            // 寻找主元
            int maxRow = i;
            for (int j = i; j < n; ++j) {
                if (abs(A[j][i]) > abs(A[maxRow][i])) {
                    maxRow = j;
                }
            }

            // 交换行
            swap(A[i], A[maxRow]);
            swap(b[i], b[maxRow]);

            // 归一化主元行
            double pivot = A[i][i];
            if (abs(pivot) < 1e-10)
                continue; // 避免除以零
            for (int j = i; j < n; ++j) {
                A[i][j] /= pivot;
            }
            b[i] /= pivot;

            // 消去其他行
            for (int j = 0; j < n; ++j) {
                if (j != i && abs(A[j][i]) > 1e-10) {
                    double factor = A[j][i];
                    for (int k = i; k < n; ++k) {
                        A[j][k] -= factor * A[i][k];
                    }
                    b[j] -= factor * b[i];
                }
            }
        }
    } catch(const std::exception& e) {
        throw e;
    }
    return b;
}

vector<double> GaussPolyCoe::PolynomialFit(const vector<double>& x, const vector<double>& y, int degree)
{
    vector<double> coeffs;
    try {
        int n = x.size();
        int m = degree + 1; // 系数个数

        // 构造正规方程组的矩阵 A^T*A 和向量 A^T*b
        vector<vector<double>> AT_A(m, vector<double>(m, 0.0));
        vector<double> AT_b(m, 0.0);

        for (int i = 0; i < m; ++i) {
            for (int j = 0; j < m; ++j) {
                for (int k = 0; k < n; ++k) {
                    AT_A[i][j] += pow(x[k], i + j);
                }
            }

            for (int k = 0; k < n; ++k) {
                AT_b[i] += y[k] * pow(x[k], i);
            }
        }

        // 求解线性方程组得到系数
        coeffs = GaussianElimination(AT_A, AT_b);
    } catch(const std::exception& e) {
       throw e;
    }
    return coeffs;
}

double GaussPolyCoe::Predict(const vector<double>& coeffs, double x)
{
    double result = 0.0f;
    for (int i = 0; i < coeffs.size(); ++i) {
        result += coeffs[i] * pow(x, i);
    }
    result = round(result * 1000) / 1000.0f;
    return result;
}

double GaussPolyCoe::MeanSquareError(const vector<double>& x, const vector<double>& y, const vector<double>& coeffs)
{
    int n = x.size();
    double mse = 0.0;
    for (int i = 0; i < n; ++i) {
        double diff = y[i] - Predict(coeffs, x[i]);
        mse += diff * diff;
    }
    return mse / n;
}