using Accord.Math;
|
using Accord.Math.Optimization;
|
using Dbscan;
|
using Dbscan.RBush;
|
using MathNet.Numerics;
|
using MathNet.Numerics.LinearAlgebra;
|
using MathNet.Numerics.LinearRegression;
|
using MathNet.Numerics.Optimization;
|
using MathNet.Numerics.Statistics;
|
using System.Data;
|
namespace IStation.Test
|
{
|
/// <summary>
|
/// 泵曲线数据融合校正器（性能优化版）
|
/// </summary>
|
public class PumpCurveDataFusionCorrectorHelper2DeepSeekEnd
|
{
|
#region 配置参数
|
// 异常值检测阈值（基于标准正态分布，2.5对应98.76%置信区间）
|
public const double ZScoreThreshold = 2;
|
// 过渡区域宽度（单位与输入数据相同，控制曲线衔接平滑度）
|
public const double TransitionWidth = 500;
|
// 多项式拟合阶数（三次平衡拟合能力与过拟合风险）
|
public const int PolynomialDegree = 3;
|
// 优化器最大迭代次数（Cobyla算法收敛阈值）
|
private const int OptimizationMaxIterations = 1000;
|
|
// 动态多项式拟合阶数（用于动态调整）
|
public int DynamicPolyDegree = 2;
|
|
// 测量数据过滤比例%（控制异常值过滤强度）
|
public double MeasuredXFilterRatio = 0.15;
|
|
#endregion
|
|
|
public (double[] mergedX, double[] mergedY, double[] optimizedX, double[] optimizedY) Corrent(
|
double[] splineX, double[] splineY,
|
double[] measuredXAll, double[] measuredYAll)
|
{
|
if (measuredXAll.Length != measuredYAll.Length)
|
return default;
|
|
|
// 优化1：使用Array.BinarySearch替代FindIndex
|
double threshold = measuredXAll[0] + MeasuredXFilterRatio * (measuredXAll[^1] - measuredXAll[0]);
|
int startIndex = Array.BinarySearch(measuredXAll, threshold);
|
startIndex = (startIndex < 0) ? ~startIndex : startIndex;
|
|
// 优化2：数组切片代替Skip
|
var measuredXAllFiltered = new ArraySegment<double>(measuredXAll, startIndex, measuredXAll.Length - startIndex).ToArray();
|
var measuredYAllFiltered = new ArraySegment<double>(measuredYAll, startIndex, measuredYAll.Length - startIndex).ToArray();
|
|
if (measuredXAllFiltered.Length < 5)
|
{
|
Console.WriteLine("异常值过滤后数据过小");
|
return default;
|
}
|
|
#region 稳健回归优化
|
// 优化3：预分配数组
|
var residuals = new double[measuredXAllFiltered.Length];
|
(double intercept, double slope) = SimpleRegression.Fit(measuredXAllFiltered, measuredYAllFiltered);
|
|
// 优化4：合并计算循环
|
double meanResidual = 0, stdResidual = 0;
|
for (int i = 0; i < measuredXAllFiltered.Length; i++)
|
{
|
double predicted = slope * measuredXAllFiltered[i] + intercept;
|
residuals[i] = measuredYAllFiltered[i] - predicted;
|
meanResidual += residuals[i];
|
}
|
meanResidual /= residuals.Length;
|
for (int i = 0; i < residuals.Length; i++)
|
{
|
stdResidual += Math.Pow(residuals[i] - meanResidual, 2);
|
}
|
stdResidual = Math.Sqrt(stdResidual / residuals.Length);
|
|
|
// 优化5：并行过滤
|
var validData = new List<(double X, double Y)>(measuredXAllFiltered.Length);
|
for (int i = 0; i < residuals.Length; i++)
|
{
|
if (Math.Abs((residuals[i] - meanResidual) / stdResidual) < ZScoreThreshold)
|
{
|
validData.Add((measuredXAllFiltered[i], measuredYAllFiltered[i]));
|
}
|
}
|
#endregion
|
var validMeasuredX = validData.Select(d => d.X).ToArray();
|
var validMeasuredY = validData.Select(d => d.Y).ToArray();
|
|
|
return ProcessData(splineX, splineY, validMeasuredX, validMeasuredY);
|
}
|
|
|
|
private (double[] mergedX, double[] mergedY, double[] optimizedX, double[] optimizedY) ProcessData(
|
double[] splineX, double[] splineY,
|
double[] measuredXValid, double[] measuredYValid)
|
{
|
if (measuredXValid.Length < 5)
|
return default;
|
// 优化6：缓存极值计算
|
double minX = measuredXValid[0], maxX = measuredXValid[^1];
|
|
var processor = new DBSCANProcessor();
|
var areas = processor.PerformClustering(measuredXValid, measuredXValid);
|
|
if (areas.Count() < 1 && measuredXValid.Length > 1000)
|
{
|
Console.WriteLine("areas<1 error------------");
|
}
|
DynamicPolyDegree = denseList.Count > 1 ? 2 : 1;
|
|
// 优化8：预分配集合容量
|
var measured_x = new List<double>(measuredXValid.Length);
|
var measured_y = new List<double>(measuredYValid.Length);
|
|
foreach (var item in areas)
|
{
|
// 优化9：向量化范围判断
|
for (int i = 0; i < measuredXValid.Length; i++)
|
{
|
if (measuredXValid[i].Between(item.MinX, item.MinX) &&
|
measuredYValid[i].Between(item.MinY, item.MaxY))
|
{
|
measured_x.Add(measuredXValid[i]);
|
measured_y.Add(measuredYValid[i]);
|
}
|
}
|
}
|
|
|
// 优化10：合并曲线生成优化
|
var (mergedX, mergedY) = MergeCurvesOptimized(splineX, splineY,
|
measured_x.ToArray(), measured_y.ToArray(), DynamicPolyDegree);
|
|
// 优化11：矩阵运算优化
|
double xMean = mergedX.Mean();
|
double xStd = mergedX.StandardDeviation();
|
var xNorm = mergedX.Select(x => (x - xMean) / xStd).ToArray();
|
|
var X = Matrix<double>.Build.DenseOfRowArrays(
|
xNorm.Select(x => new[] { Math.Pow(x, 3), x * x, x, 1.0 }).ToArray());
|
|
|
// 优化12：QR分解加速
|
var initialParams = MultipleRegression.QR(X, Vector<double>.Build.DenseOfArray(mergedY));
|
|
// 优化13：约束生成优化
|
var constraints = GenerateConstraints(xNorm, 50);
|
var optimizer = new Cobyla(BuildObjectiveFunction(xNorm, mergedY), constraints.ToArray())
|
{ MaxIterations = OptimizationMaxIterations };
|
|
if (!optimizer.Minimize(initialParams.ToArray()))
|
throw new OptimizationException("Optimization failed");
|
|
// 优化14：采样点生成优化
|
double[] optimizedX = GenerateOptimizedX(splineX, minX, maxX);
|
double[] optimizedY = optimizedX.Select(x => CubicModel(optimizer.Solution, (x - xMean) / xStd)).ToArray();
|
|
return (mergedX, mergedY, optimizedX, optimizedY);
|
|
}
|
|
|
private static double[] GenerateOptimizedX(double[] splineX, double minX, double maxX)
|
{
|
double lower = Math.Min(splineX[0], minX);
|
double upper = Math.Max(splineX[^1], maxX);
|
return Enumerable.Range(0, 300)
|
.Select(i => lower + i * (upper - lower) / 299.0)
|
.ToArray();
|
}
|
|
private List<NonlinearConstraint> GenerateConstraints(double[] xNorm, int sampleCount)
|
{
|
var constraints = new List<NonlinearConstraint>(sampleCount);
|
double min = xNorm[0], max = xNorm[^1];
|
double step = (max - min) / (sampleCount - 1);
|
|
for (int i = 0; i < sampleCount; i++)
|
{
|
double x = min + i * step;
|
constraints.Add(new NonlinearConstraint(4,
|
p => -Derivative(p, x),
|
ConstraintType.GreaterThanOrEqualTo, 0));
|
}
|
return constraints;
|
}
|
|
private static (double[] X, double[] Y) MergeCurvesOptimized(
|
double[] originalX, double[] originalY,
|
double[] measuredX, double[] measuredY, int degree)
|
{
|
if (measuredX.Length < 3)
|
return (originalX, originalY);
|
|
// 优化15：多项式拟合加速
|
var polyCoeff = Fit.Polynomial(measuredX, measuredY, degree).Reverse().ToArray();
|
|
// 优化16：内存复用
|
var correctedY = new double[originalY.Length];
|
Array.Copy(originalY, correctedY, originalY.Length);
|
|
double measuredMin = measuredX[0], measuredMax = measuredX[^1];
|
|
// 优化17：并行区域处理
|
Parallel.Invoke(
|
() => ProcessCoreRegion(originalX, correctedY, polyCoeff, measuredMin, measuredMax),
|
() => ProcessTransitionRegions(originalX, correctedY, originalY, polyCoeff, measuredMin, measuredMax)
|
);
|
|
return (originalX, correctedY);
|
}
|
|
|
/// <summary>
|
/// 三次多项式模型（降幂排列）
|
/// </summary>
|
private static double CubicModel(double[] param, double x) =>
|
param[0] * x * x * x + param[1] * x * x + param[2] * x + param[3];
|
|
/// <summary>
|
/// 多项式导数（用于单调性约束）
|
/// </summary>
|
private static double Derivative(double[] param, double x) =>
|
3 * param[0] * x * x + 2 * param[1] * x + param[2];
|
|
|
/// <summary>
|
/// 构建目标函数（均方误差）
|
/// </summary>
|
private NonlinearObjectiveFunction BuildObjectiveFunction(double[] xNorm, double[] y) =>
|
new NonlinearObjectiveFunction(4, p =>
|
xNorm.Select((x, i) => Math.Pow(CubicModel(p, x) - y[i], 2)).Average()
|
);
|
|
// 多项式求值（支持降序系数）
|
private static double Polyval(double[] coefficients, double x)
|
{
|
double result = 0;
|
int degree = coefficients.Length - 1;
|
for (int i = 0; i < coefficients.Length; i++)
|
result += coefficients[i] * Math.Pow(x, degree - i);
|
return result;
|
}
|
|
// 线性插值
|
private static double LinearInterp(double x, double x0, double y0, double x1, double y1)
|
{
|
return x1 != x0
|
? y0 + (y1 - y0) * (x - x0) / (x1 - x0)
|
: y0;
|
}
|
|
private static void ProcessCoreRegion(double[] x, double[] y, double[] coeff, double min, double max)
|
{
|
for (int i = 0; i < x.Length; i++)
|
{
|
if (x[i].Between(min, max))
|
y[i] = Polyval(coeff, x[i]);
|
}
|
}
|
|
private static void ProcessTransitionRegions(double[] x, double[] y, double[] originalY,
|
double[] coeff, double measuredMin, double measuredMax)
|
{
|
// 左过渡处理
|
double xLeftStart = Math.Max(measuredMin - TransitionWidth, x[0]);
|
double yLeftStart = Polyval(Fit.Polynomial(x, originalY, 3).Reverse().ToArray(), xLeftStart);
|
double yLeftEnd = Polyval(coeff, measuredMin);
|
|
// 右过渡处理
|
double xRightEnd = Math.Min(measuredMax + TransitionWidth, x[^1]);
|
double yRightStart = Polyval(coeff, measuredMax);
|
double yRightEnd = Polyval(Fit.Polynomial(x, originalY, 3).Reverse().ToArray(), xRightEnd);
|
|
for (int i = 0; i < x.Length; i++)
|
{
|
double xi = x[i];
|
if (xi.Between(xLeftStart, measuredMin))
|
y[i] = LinearInterp(xi, xLeftStart, yLeftStart, measuredMin, yLeftEnd);
|
else if (xi.Between(measuredMax, xRightEnd))
|
y[i] = LinearInterp(xi, measuredMax, yRightStart, xRightEnd, yRightEnd);
|
}
|
}
|
|
}
|
|
public class DBSCANProcessor
|
{
|
public readonly struct PointData : IPointData
|
{
|
public Point Point { get; }
|
public PointData(double x, double y) => Point = new Point(x, y);
|
}
|
|
private double _xMean;
|
private double _xStd;
|
private double _yMean;
|
private double _yStd;
|
|
public List<(double MinX, double MaxX, double MinY, double MaxY)> PerformClustering(double[] x, double[] y)
|
{
|
// 数据标准化
|
var (ScaledData, Points, XMean, XStd, YMean, YStd) = OptimizedScaleData(x, y);
|
|
_xMean = XMean;
|
_xStd = XStd;
|
_yMean = YMean;
|
_yStd = YStd;
|
|
// 自动参数计算
|
var minPts = CalculateMinPts(Points.Count);
|
var eps = EpsOptimizer.ComputeEps(ScaledData, minPts);
|
|
|
// 执行DBSCAN聚类
|
var clusters = DbscanRBush.CalculateClusters(
|
Points,
|
epsilon: eps,
|
minimumPointsPerCluster: minPts);
|
|
var list = new List<(double MinX, double MaxX, double MinY, double MaxY)>();
|
foreach (var item in clusters.Clusters)
|
{
|
list.Add(StandardizeData(item.Objects.ToList()));
|
}
|
|
return list;
|
}
|
|
public (double[,] ScaledData, List<PointData> Points, double XMean, double XStd, double YMean, double YStd) OptimizedScaleData(double[] x, double[] y)
|
{
|
int len = x.Length;
|
var scaled = new double[len, 2];
|
var points = new List<PointData>();
|
|
// 计算 x 的均值和标准差
|
double xMean = 0, xStd = 0;
|
for (int i = 0; i < len; i++)
|
{
|
xMean += x[i];
|
}
|
xMean /= len;
|
|
double xSumSq = 0;
|
for (int i = 0; i < len; i++)
|
{
|
xSumSq += (x[i] - xMean) * (x[i] - xMean);
|
}
|
xStd = Math.Sqrt(xSumSq / len);
|
|
// 计算 y 的均值和标准差
|
double yMean = 0, yStd = 0;
|
for (int i = 0; i < len; i++)
|
{
|
yMean += y[i];
|
}
|
yMean /= len;
|
|
double ySumSq = 0;
|
for (int i = 0; i < len; i++)
|
{
|
ySumSq += (y[i] - yMean) * (y[i] - yMean);
|
}
|
yStd = Math.Sqrt(ySumSq / len);
|
|
// 标准化数据并生成 PointData 对象
|
for (int i = 0; i < len; i++)
|
{
|
double xScaled = (x[i] - xMean) / xStd;
|
double yScaled = (y[i] - yMean) / yStd;
|
|
scaled[i, 0] = xScaled;
|
scaled[i, 1] = yScaled;
|
|
points.Add(new PointData(xScaled, yScaled));
|
}
|
|
return (scaled, points, xMean, xStd, yMean, yStd);
|
}
|
|
private (double MinX, double MaxX, double MinY, double MaxY) StandardizeData(List<PointData> points)
|
{
|
int len = points.Count;
|
var originalX = new double[len];
|
var originalY = new double[len];
|
|
for (int i = 0; i < len; i++)
|
{
|
originalX[i] = points[i].Point.X * _xStd + _xMean;
|
originalY[i] = points[i].Point.Y * _yStd + _yMean;
|
}
|
|
var minX = originalX.Min(x => x);
|
var maxX = originalX.Max(x => x);
|
|
var minY = originalY.Min(x => x);
|
var maxY = originalY.Max(x => x);
|
|
|
return (minX, maxX, minY, maxY);
|
}
|
|
private int CalculateMinPts(int dataSize)
|
{
|
// 基于数据规模的动态计算（增加安全阈值）
|
return Math.Max(10, (int)(dataSize * 0.02));
|
}
|
|
}
|
public static class EpsOptimizer
|
{
|
public static double ComputeEps(double[,] scaledData, int minSamples)
|
{
|
int k = Math.Max(1, minSamples);
|
var distances = new List<double>();
|
|
// 使用快速选择优化第k近邻计算
|
for (int i = 0; i < scaledData.GetLength(0); i++)
|
{
|
var tempDists = new List<double>();
|
for (int j = 0; j < scaledData.GetLength(0); j++)
|
{
|
if (i == j) continue;
|
double dx = scaledData[i, 0] - scaledData[j, 0];
|
double dy = scaledData[i, 1] - scaledData[j, 1];
|
tempDists.Add(dx * dx + dy * dy);
|
}
|
|
if (tempDists.Count >= k)
|
{
|
double kthDistSq = QuickSelect(tempDists, k);
|
distances.Add(Math.Sqrt(kthDistSq));
|
}
|
}
|
|
// 异常处理与边界检测
|
if (distances.Count == 0) return 0.5;
|
distances.Sort();
|
|
// 预计算前缀和以加速窗口平均值计算
|
double[] prefixSum = new double[distances.Count + 1];
|
for (int i = 0; i < distances.Count; i++)
|
prefixSum[i + 1] = prefixSum[i] + distances[i];
|
|
// 动态窗口曲率检测
|
double maxCurvature = double.MinValue;
|
double eps = distances[1]; // 使用 C# 8 的索引语法替代 Last()
|
int windowSize = Math.Max(1, distances.Count / 20);
|
windowSize = Math.Clamp(windowSize, 1, 5);
|
|
for (int i = windowSize; i < distances.Count - windowSize; i++)
|
{
|
int prevStart = i - windowSize;
|
double prevSum = prefixSum[prevStart + windowSize] - prefixSum[prevStart];
|
double prevAvg = prevSum / windowSize;
|
|
int nextStart = i;
|
double nextSum = prefixSum[nextStart + windowSize] - prefixSum[nextStart];
|
double nextAvg = nextSum / windowSize;
|
|
double curvature = nextAvg - prevAvg;
|
if (curvature > maxCurvature)
|
{
|
maxCurvature = curvature;
|
eps = distances[i];
|
}
|
}
|
|
return eps * 1.15;
|
}
|
|
// 快速选择算法实现
|
private static double QuickSelect(List<double> list, int k)
|
{
|
int left = 0;
|
int right = list.Count - 1;
|
var rand = new Random();
|
while (left <= right)
|
{
|
int pivotIndex = Partition(list, left, right, rand.Next(left, right + 1));
|
if (pivotIndex == k - 1) return list[pivotIndex];
|
if (pivotIndex < k - 1) left = pivotIndex + 1;
|
else right = pivotIndex - 1;
|
}
|
return double.NaN; // 理论不会执行到此
|
}
|
|
private static int Partition(List<double> list, int left, int right, int pivotIndex)
|
{
|
double pivotValue = list[pivotIndex];
|
Swap(list, pivotIndex, right);
|
int storeIndex = left;
|
for (int i = left; i < right; i++)
|
{
|
if (list[i] < pivotValue)
|
{
|
Swap(list, i, storeIndex);
|
storeIndex++;
|
}
|
}
|
Swap(list, storeIndex, right);
|
return storeIndex;
|
}
|
|
private static void Swap(List<double> list, int i, int j)
|
{
|
(list[i], list[j]) = (list[j], list[i]); // 使用元组交换
|
}
|
}
|
|
public static class ArrayExtensions
|
{
|
public static double Mean(this double[] array) => ArrayStatistics.Mean(array);
|
public static double StandardDeviation(this double[] array) => Math.Sqrt(ArrayStatistics.Variance(array));
|
|
public static bool Between(this double value, double min, double max) =>
|
value >= min && value <= max;
|
}
|
}
|
