gunshi-project-ss/src/main/java/com/gunshi/project/hsz/util/ProjectCalculateUtil.java

package com.gunshi.project.hsz.util;

import com.gunshi.project.hsz.common.model.vo.OsmoticPressDetailVo;
import com.gunshi.project.hsz.model.RegressionEquation;

import java.math.BigDecimal;
import java.math.RoundingMode;
import java.util.List;
import java.util.stream.Collectors;

/**
 * 工程安全分析-线性回归方程计算
 */
public class ProjectCalculateUtil {

    private static List<OsmoticPressDetailVo> filterValidData(List<OsmoticPressDetailVo> data) {
        return data.stream()
                .filter(vo -> vo.getRz() != null && vo.getValue() != null)
                .collect(Collectors.toList());
    }

    private static BigDecimal calculateRSquared(List<OsmoticPressDetailVo> data,
                                                List<BigDecimal> coefficients,
                                                BigDecimal meanY) {
        BigDecimal ssTot = BigDecimal.ZERO;
        BigDecimal ssRes = BigDecimal.ZERO;

        for (OsmoticPressDetailVo vo : data) {
            BigDecimal x = vo.getRz();
            BigDecimal y = vo.getValue();
            BigDecimal yPredicted = predictWithCoefficients(x, coefficients);

            BigDecimal diffActual = y.subtract(meanY);
            ssTot = ssTot.add(diffActual.multiply(diffActual));

            BigDecimal diffResidual = y.subtract(yPredicted);
            ssRes = ssRes.add(diffResidual.multiply(diffResidual));
        }

        if (ssTot.compareTo(BigDecimal.ZERO) == 0) {
            return BigDecimal.ONE;
        }

        return BigDecimal.ONE.subtract(ssRes.divide(ssTot, 30, RoundingMode.HALF_UP));
    }

    private static BigDecimal predictWithCoefficients(BigDecimal x, List<BigDecimal> coefficients) {
        BigDecimal result = BigDecimal.ZERO;
        BigDecimal xPower = BigDecimal.ONE;

        for (BigDecimal coeff : coefficients) {
            result = result.add(coeff.multiply(xPower));
            xPower = xPower.multiply(x);
        }

        return result;
    }

    private static BigDecimal calculateDeterminant3x3(BigDecimal a11, BigDecimal a12, BigDecimal a13,
                                                      BigDecimal a21, BigDecimal a22, BigDecimal a23,
                                                      BigDecimal a31, BigDecimal a32, BigDecimal a33) {
        return a11.multiply(a22.multiply(a33).subtract(a23.multiply(a32)))
                .subtract(a12.multiply(a21.multiply(a33).subtract(a23.multiply(a31))))
                .add(a13.multiply(a21.multiply(a32).subtract(a22.multiply(a31))));
    }

    /**
     * 计算线性回归方程（1阶） 库水位x轴，检测值y轴
     */
    public static RegressionEquation calculateLinear(List<OsmoticPressDetailVo> data) {
        List<OsmoticPressDetailVo> validData = filterValidData(data);
        if (validData.isEmpty()) {
            throw new IllegalArgumentException("没有有效的数据点进行计算");
        }

        int n = validData.size();

        BigDecimal sumX = BigDecimal.ZERO;
        BigDecimal sumY = BigDecimal.ZERO;
        BigDecimal sumXY = BigDecimal.ZERO;
        BigDecimal sumX2 = BigDecimal.ZERO;
        BigDecimal sumY2 = BigDecimal.ZERO;

        for (OsmoticPressDetailVo vo : validData) {
            BigDecimal x = vo.getRz();
            BigDecimal y = vo.getValue();

            sumX = sumX.add(x);
            sumY = sumY.add(y);
            sumXY = sumXY.add(x.multiply(y));
            sumX2 = sumX2.add(x.multiply(x));
            sumY2 = sumY2.add(y.multiply(y));
        }

        BigDecimal nBig = new BigDecimal(n);
        BigDecimal meanX = sumX.divide(nBig, 30, RoundingMode.HALF_UP);
        BigDecimal meanY = sumY.divide(nBig, 30, RoundingMode.HALF_UP);

        BigDecimal numerator = sumXY.subtract(nBig.multiply(meanX).multiply(meanY));
        BigDecimal denominator = sumX2.subtract(nBig.multiply(meanX).multiply(meanX));

        if (denominator.compareTo(BigDecimal.ZERO) == 0) {
            throw new IllegalArgumentException("数据点x值相同，无法计算斜率");
        }

        BigDecimal slope = numerator.divide(denominator, 30, RoundingMode.HALF_UP);
        BigDecimal intercept = meanY.subtract(slope.multiply(meanX));

        // 计算R²
        BigDecimal rSquared = calculateRSquared(validData,
                List.of(intercept, slope), meanY);

        // 系数列表：[常数项, 一次项系数]
        List<BigDecimal> coefficients = List.of(
                intercept.setScale(10, RoundingMode.HALF_UP),
                slope.setScale(10, RoundingMode.HALF_UP)
        );

        return new RegressionEquation(1, coefficients, rSquared, n);
    }

    /**
     * 计算二次回归方程（2阶）
     */
    public static RegressionEquation calculateQuadratic(List<OsmoticPressDetailVo> data) {
        List<OsmoticPressDetailVo> validData = filterValidData(data);
        if (validData.size() < 3) {
            throw new IllegalArgumentException("二次回归至少需要3个数据点");
        }

        int n = validData.size();

        BigDecimal sumX = BigDecimal.ZERO, sumY = BigDecimal.ZERO;
        BigDecimal sumX2 = BigDecimal.ZERO, sumX3 = BigDecimal.ZERO, sumX4 = BigDecimal.ZERO;
        BigDecimal sumXY = BigDecimal.ZERO, sumX2Y = BigDecimal.ZERO;

        for (OsmoticPressDetailVo vo : validData) {
            BigDecimal x = vo.getRz();
            BigDecimal y = vo.getValue();

            BigDecimal x2 = x.multiply(x);
            BigDecimal x3 = x2.multiply(x);
            BigDecimal x4 = x3.multiply(x);

            sumX = sumX.add(x);
            sumY = sumY.add(y);
            sumX2 = sumX2.add(x2);
            sumX3 = sumX3.add(x3);
            sumX4 = sumX4.add(x4);
            sumXY = sumXY.add(x.multiply(y));
            sumX2Y = sumX2Y.add(x2.multiply(y));
        }

        BigDecimal nBig = new BigDecimal(n);

        // 使用克莱姆法则求解
        BigDecimal det = calculateDeterminant3x3(
                nBig, sumX, sumX2,
                sumX, sumX2, sumX3,
                sumX2, sumX3, sumX4
        );

        if (det.compareTo(BigDecimal.ZERO) == 0) {
            throw new IllegalArgumentException("矩阵奇异，无法求解二次回归方程");
        }

        BigDecimal detA0 = calculateDeterminant3x3(
                sumY, sumX, sumX2,
                sumXY, sumX2, sumX3,
                sumX2Y, sumX3, sumX4
        );

        BigDecimal detA1 = calculateDeterminant3x3(
                nBig, sumY, sumX2,
                sumX, sumXY, sumX3,
                sumX2, sumX2Y, sumX4
        );

        BigDecimal detA2 = calculateDeterminant3x3(
                nBig, sumX, sumY,
                sumX, sumX2, sumXY,
                sumX2, sumX3, sumX2Y
        );

        BigDecimal a0 = detA0.divide(det, 30, RoundingMode.HALF_UP);
        BigDecimal a1 = detA1.divide(det, 30, RoundingMode.HALF_UP);
        BigDecimal a2 = detA2.divide(det, 30, RoundingMode.HALF_UP);

        // 计算R²
        BigDecimal meanY = sumY.divide(nBig, 30, RoundingMode.HALF_UP);
        BigDecimal rSquared = calculateRSquared(validData,
                List.of(a0, a1, a2), meanY);

        // 系数列表：[常数项, 一次项系数, 二次项系数]
        List<BigDecimal> coefficients = List.of(
                a0.setScale(10, RoundingMode.HALF_UP),
                a1.setScale(10, RoundingMode.HALF_UP),
                a2.setScale(10, RoundingMode.HALF_UP)
        );

        return new RegressionEquation(2, coefficients, rSquared, n);
    }

}