diff --git a/src/ChebyshevSharp/ChebyshevApproximation.cs b/src/ChebyshevSharp/ChebyshevApproximation.cs index be1e812..b519896 100644 --- a/src/ChebyshevSharp/ChebyshevApproximation.cs +++ b/src/ChebyshevSharp/ChebyshevApproximation.cs @@ -11,6 +11,14 @@ namespace ChebyshevSharp; /// public class ChebyshevApproximation { + private double[][] _domain = Array.Empty(); + private int[] _nNodes = Array.Empty(); + private double[][] _nodeArrays = Array.Empty(); + private double[]? _tensorValues; + private double[][]? _weights; + private double[][,]? _diffMatrices; + private double[][]? _diffMatricesTFlat; + /// The function to approximate. Null after load or from_values. public Func? Function { get; internal set; } @@ -18,28 +26,92 @@ public class ChebyshevApproximation public int NumDimensions { get; internal set; } /// Domain bounds for each dimension, as list of [lo, hi]. - public double[][] Domain { get; internal set; } = Array.Empty(); + public double[][] Domain + { + get => CloneHelpers.DeepCopy(_domain)!; + internal set => _domain = value ?? Array.Empty(); + } /// Number of Chebyshev nodes per dimension. - public int[] NNodes { get; internal set; } = Array.Empty(); + public int[] NNodes + { + get => CloneHelpers.DeepCopy(_nNodes)!; + internal set => _nNodes = value ?? Array.Empty(); + } /// Maximum supported derivative order. public int MaxDerivativeOrder { get; internal set; } = 2; /// Chebyshev nodes per dimension, each sorted ascending. - public double[][] NodeArrays { get; internal set; } = Array.Empty(); + public double[][] NodeArrays + { + get => CloneHelpers.DeepCopy(_nodeArrays)!; + internal set => _nodeArrays = value ?? Array.Empty(); + } /// Flat tensor of function values at all node combinations (C-order). - public double[]? TensorValues { get; internal set; } + public double[]? TensorValues + { + get => CloneHelpers.DeepCopy(_tensorValues); + internal set => _tensorValues = value; + } /// Barycentric weights per dimension. - public double[][]? Weights { get; internal set; } + public double[][]? Weights + { + get => CloneHelpers.DeepCopy(_weights); + internal set => _weights = value; + } /// Spectral differentiation matrices per dimension. - public double[][,]? DiffMatrices { get; internal set; } + public double[][,]? DiffMatrices + { + get => CloneHelpers.DeepCopy(_diffMatrices); + internal set => _diffMatrices = value; + } /// Pre-transposed diff matrices flattened to double[] for BLAS GEMM (row-major). - internal double[][]? DiffMatricesTFlat { get; set; } + internal double[][]? DiffMatricesTFlat + { + get => _diffMatricesTFlat; + set => _diffMatricesTFlat = value; + } + + internal double[][] DomainStorage + { + get => _domain; + set => _domain = value; + } + + internal int[] NNodesStorage + { + get => _nNodes; + set => _nNodes = value; + } + + internal double[][] NodeArraysStorage + { + get => _nodeArrays; + set => _nodeArrays = value; + } + + internal double[]? TensorValuesStorage + { + get => _tensorValues; + set => _tensorValues = value; + } + + internal double[][]? WeightsStorage + { + get => _weights; + set => _weights = value; + } + + internal double[][,]? DiffMatricesStorage + { + get => _diffMatrices; + set => _diffMatrices = value; + } /// Time taken by Build() in seconds. public double BuildTime { get; internal set; } @@ -109,8 +181,8 @@ public ChebyshevApproximation( Function = function; NumDimensions = numDimensions; - Domain = domain.Select(d => (double[])d.Clone()).ToArray(); - NNodes = (int[])nNodes.Clone(); + _domain = domain.Select(d => (double[])d.Clone()).ToArray(); + _nNodes = (int[])nNodes.Clone(); MaxDerivativeOrder = maxDerivativeOrder; _additionalData = additionalData; _nWorkers = Internal.ParallelBuild.NormalizeNWorkers(nWorkers); @@ -119,15 +191,15 @@ public ChebyshevApproximation( if (!deferBuild) { // Generate Chebyshev nodes for each dimension - NodeArrays = new double[numDimensions][]; + _nodeArrays = new double[numDimensions][]; for (int d = 0; d < numDimensions; d++) { - NodeArrays[d] = BarycentricKernel.MakeNodesForDim(domain[d][0], domain[d][1], nNodes[d]); + _nodeArrays[d] = BarycentricKernel.MakeNodesForDim(domain[d][0], domain[d][1], nNodes[d]); } } else { - NodeArrays = Array.Empty(); + _nodeArrays = Array.Empty(); } } @@ -195,7 +267,7 @@ public ChebyshevApproximation( Function = function; NumDimensions = numDimensions; - Domain = domain.Select(d => (double[])d.Clone()).ToArray(); + _domain = domain.Select(d => (double[])d.Clone()).ToArray(); ErrorThreshold = errorThreshold; MaxN = maxN; MaxDerivativeOrder = maxDerivativeOrder; @@ -204,20 +276,20 @@ public ChebyshevApproximation( _nWorkers = Internal.ParallelBuild.NormalizeNWorkers(nWorkers); _progress = progress; - // If all entries are non-null, populate NNodes + nodes immediately (matches existing fixed-N behavior). + // If all entries are non-null, populate resolved node counts immediately. if (resolved.All(n => n != null)) { - NNodes = resolved.Select(n => n!.Value).ToArray(); - ValidateNodeCounts(nameof(ChebyshevApproximation), NNodes); - NodeArrays = new double[numDimensions][]; + _nNodes = resolved.Select(n => n!.Value).ToArray(); + ValidateNodeCounts(nameof(ChebyshevApproximation), _nNodes); + _nodeArrays = new double[numDimensions][]; for (int d = 0; d < numDimensions; d++) - NodeArrays[d] = BarycentricKernel.MakeNodesForDim(domain[d][0], domain[d][1], NNodes[d]); + _nodeArrays[d] = BarycentricKernel.MakeNodesForDim(domain[d][0], domain[d][1], _nNodes[d]); } else { - // Auto-N path: NNodes left empty until Build() resolves. - NNodes = Array.Empty(); - NodeArrays = Array.Empty(); + // Auto-N path: resolved node counts are empty until Build(). + _nNodes = Array.Empty(); + _nodeArrays = Array.Empty(); } } @@ -250,9 +322,9 @@ public void Build(bool verbose = true) internal void BuildFixedGrid(bool verbose = true) { int total = TensorShape.RequireArrayLength( - TensorShape.CheckedProduct(NNodes, nameof(BuildFixedGrid)), + TensorShape.CheckedProduct(_nNodes, nameof(BuildFixedGrid)), nameof(BuildFixedGrid), - NNodes); + _nNodes); if (verbose) Console.WriteLine($"Building {NumDimensions}D Chebyshev approximation ({total:N0} evaluations)..."); @@ -269,29 +341,29 @@ internal void BuildFixedGrid(bool verbose = true) int rem = flat; for (int d = NumDimensions - 1; d >= 0; d--) { - indices[d] = rem % NNodes[d]; - rem /= NNodes[d]; + indices[d] = rem % _nNodes[d]; + rem /= _nNodes[d]; } var pt = new double[NumDimensions]; for (int d = 0; d < NumDimensions; d++) - pt[d] = NodeArrays[d][indices[d]]; + pt[d] = _nodeArrays[d][indices[d]]; points[flat] = pt; } var tensorValues = Internal.ParallelBuild.EvaluateInParallel( Function!, points, _additionalData, _nWorkers, _progress); ValidateFiniteBuildValues(tensorValues); - TensorValues = tensorValues; + _tensorValues = tensorValues; NEvaluations = total; // Step 2: Pre-compute barycentric weights - Weights = new double[NumDimensions][]; + _weights = new double[NumDimensions][]; for (int d = 0; d < NumDimensions; d++) - Weights[d] = BarycentricKernel.ComputeBarycentricWeights(NodeArrays[d]); + _weights[d] = BarycentricKernel.ComputeBarycentricWeights(_nodeArrays[d]); // Step 3: Pre-compute differentiation matrices - DiffMatrices = new double[NumDimensions][,]; + _diffMatrices = new double[NumDimensions][,]; for (int d = 0; d < NumDimensions; d++) - DiffMatrices[d] = BarycentricKernel.ComputeDifferentiationMatrix(NodeArrays[d], Weights[d]); + _diffMatrices[d] = BarycentricKernel.ComputeDifferentiationMatrix(_nodeArrays[d], _weights[d]); // Step 4: Pre-transpose diff matrices for VectorizedEval PrecomputeTransposedDiffMatrices(); @@ -301,7 +373,7 @@ internal void BuildFixedGrid(bool verbose = true) if (verbose) { - int totalWeights = Weights.Sum(w => w.Length); + int totalWeights = _weights.Sum(w => w.Length); Console.WriteLine($" Built in {BuildTime:F3}s ({totalWeights} weights, {totalWeights * 8} bytes)"); } @@ -332,23 +404,23 @@ private static void ValidateFiniteBuildValues(double[] values) /// Interpolated value or derivative at the query point. public double Eval(double[] point, int[] derivativeOrder) { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException("Call Build() first"); - EvaluationArguments.ValidatePointInDomain(point, NumDimensions, Domain); + EvaluationArguments.ValidatePointInDomain(point, NumDimensions, _domain); EvaluationArguments.ValidateDerivativeOrder(derivativeOrder, NumDimensions); // Current working data and its shape - double[] current = TensorValues; - int[] currentShape = (int[])NNodes.Clone(); + double[] current = _tensorValues; + int[] currentShape = (int[])_nNodes.Clone(); for (int d = NumDimensions - 1; d >= 0; d--) { double x = point[d]; int deriv = derivativeOrder[d]; - double[] nodes = NodeArrays[d]; - double[] weights = Weights![d]; - double[,] diffMatrix = DiffMatrices![d]; - int nNodesD = NNodes[d]; + double[] nodes = _nodeArrays[d]; + double[] weights = _weights![d]; + double[,] diffMatrix = _diffMatrices![d]; + int nNodesD = _nNodes[d]; if (d == 0) { @@ -406,10 +478,10 @@ public double Eval(double[] point, int[] derivativeOrder) /// Interpolated value at the query point. public double Eval(double[] point) { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException( "Cannot evaluate an unbuilt interpolant. Call Build() or SetOriginalFunctionValues() first."); - EvaluationArguments.ValidatePointInDomain(point, NumDimensions, Domain); + EvaluationArguments.ValidatePointInDomain(point, NumDimensions, _domain); return Eval(point, new int[NumDimensions]); } @@ -422,12 +494,12 @@ public double Eval(double[] point) /// Interpolated value or derivative. public double VectorizedEval(double[] point, int[] derivativeOrder) { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException("Call Build() first"); - EvaluationArguments.ValidatePointInDomain(point, NumDimensions, Domain); + EvaluationArguments.ValidatePointInDomain(point, NumDimensions, _domain); EvaluationArguments.ValidateDerivativeOrder(derivativeOrder, NumDimensions); - double[] current = TensorValues; + double[] current = _tensorValues; // Track tensor dimensions without shape array allocations. // leadSize = product of all dims before current last dim. @@ -438,13 +510,13 @@ public double VectorizedEval(double[] point, int[] derivativeOrder) { double x = point[d]; int deriv = derivativeOrder[d]; - int lastDim = NNodes[d]; + int lastDim = _nNodes[d]; int leadSize = totalSize / lastDim; // Apply differentiation matrix if derivative order > 0 if (deriv > 0) { - double[] dTFlat = DiffMatricesTFlat![d]; + double[] dTFlat = _diffMatricesTFlat![d]; for (int o = 0; o < deriv; o++) current = BarycentricKernel.MatmulLastAxisMatrixFlat(current, leadSize, lastDim, dTFlat, lastDim); } @@ -453,7 +525,7 @@ public double VectorizedEval(double[] point, int[] derivativeOrder) int exactIdx = -1; for (int i = 0; i < lastDim; i++) { - if (Math.Abs(x - NodeArrays[d][i]) < 1e-14) + if (Math.Abs(x - _nodeArrays[d][i]) < 1e-14) { exactIdx = i; break; @@ -475,7 +547,7 @@ public double VectorizedEval(double[] point, int[] derivativeOrder) double sumW = 0.0; for (int i = 0; i < lastDim; i++) { - double wod = Weights![d][i] / (x - NodeArrays[d][i]); + double wod = _weights![d][i] / (x - _nodeArrays[d][i]); wNorm[i] = wod; sumW += wod; } @@ -502,17 +574,17 @@ public double VectorizedEval(double[] point, int[] derivativeOrder) /// Results array of length N. public double[] VectorizedEvalBatch(double[][] points, int[] derivativeOrder) { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException("Call Build() first"); - EvaluationArguments.ValidatePointsInDomain(points, NumDimensions, Domain); + EvaluationArguments.ValidatePointsInDomain(points, NumDimensions, _domain); EvaluationArguments.ValidateDerivativeOrder(derivativeOrder, NumDimensions); // Hoist: apply all derivative-matrix matmuls once — they are point-independent. // Process from last dimension to first to match VectorizedEval ordering. - double[] tensorWithDerivs = ApplyDerivativePasses(TensorValues, NNodes, derivativeOrder); + double[] tensorWithDerivs = ApplyDerivativePasses(_tensorValues, _nNodes, derivativeOrder); double[] results = new double[points.Length]; - int totalSize = TensorValues.Length; + int totalSize = _tensorValues.Length; for (int i = 0; i < points.Length; i++) { @@ -523,14 +595,14 @@ public double[] VectorizedEvalBatch(double[][] points, int[] derivativeOrder) for (int d = NumDimensions - 1; d >= 0; d--) { double x = points[i][d]; - int lastDim = NNodes[d]; + int lastDim = _nNodes[d]; int leadSize = curSize / lastDim; // Barycentric contraction along last axis (no diff-matrix here — already hoisted) int exactIdx = -1; for (int j = 0; j < lastDim; j++) { - if (Math.Abs(x - NodeArrays[d][j]) < 1e-14) + if (Math.Abs(x - _nodeArrays[d][j]) < 1e-14) { exactIdx = j; break; @@ -550,7 +622,7 @@ public double[] VectorizedEvalBatch(double[][] points, int[] derivativeOrder) double sumW = 0.0; for (int j = 0; j < lastDim; j++) { - double wod = Weights![d][j] / (x - NodeArrays[d][j]); + double wod = _weights![d][j] / (x - _nodeArrays[d][j]); wNorm[j] = wod; sumW += wod; } @@ -585,7 +657,7 @@ private double[] ApplyDerivativePasses(double[] tensor, int[] shape, int[] deriv int deriv = derivativeOrder[d]; if (deriv > 0) { - double[,] dm = DiffMatrices![d]; + double[,] dm = _diffMatrices![d]; for (int o = 0; o < deriv; o++) result = BarycentricKernel.MatmulAlongAxis(result, shape, d, dm); } @@ -601,9 +673,9 @@ private double[] ApplyDerivativePasses(double[] tensor, int[] shape, int[] deriv /// One result per derivative order. public double[] VectorizedEvalMulti(double[] point, int[][] derivativeOrders) { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException("Call Build() first"); - EvaluationArguments.ValidatePointInDomain(point, NumDimensions, Domain); + EvaluationArguments.ValidatePointInDomain(point, NumDimensions, _domain); EvaluationArguments.ValidateDerivativeOrders(derivativeOrders, NumDimensions); // Pre-compute dimension info (shared across all derivative orders) @@ -612,9 +684,9 @@ public double[] VectorizedEvalMulti(double[] point, int[][] derivativeOrders) { double x = point[d]; int exactIdx = -1; - for (int i = 0; i < NNodes[d]; i++) + for (int i = 0; i < _nNodes[d]; i++) { - if (Math.Abs(x - NodeArrays[d][i]) < 1e-14) + if (Math.Abs(x - _nodeArrays[d][i]) < 1e-14) { exactIdx = i; break; @@ -627,19 +699,19 @@ public double[] VectorizedEvalMulti(double[] point, int[][] derivativeOrders) } else { - double[] diff = new double[NNodes[d]]; - for (int i = 0; i < NNodes[d]; i++) - diff[i] = x - NodeArrays[d][i]; + double[] diff = new double[_nNodes[d]]; + for (int i = 0; i < _nNodes[d]; i++) + diff[i] = x - _nodeArrays[d][i]; - double[] wOverDiff = new double[NNodes[d]]; + double[] wOverDiff = new double[_nNodes[d]]; double sumW = 0.0; - for (int i = 0; i < NNodes[d]; i++) + for (int i = 0; i < _nNodes[d]; i++) { - wOverDiff[i] = Weights![d][i] / diff[i]; + wOverDiff[i] = _weights![d][i] / diff[i]; sumW += wOverDiff[i]; } - double[] wNorm = new double[NNodes[d]]; - for (int i = 0; i < NNodes[d]; i++) + double[] wNorm = new double[_nNodes[d]]; + for (int i = 0; i < _nNodes[d]; i++) wNorm[i] = wOverDiff[i] / sumW; dimInfo[d] = (false, -1, wNorm); @@ -647,23 +719,23 @@ public double[] VectorizedEvalMulti(double[] point, int[][] derivativeOrders) } double[] results = new double[derivativeOrders.Length]; - int tensorSize = TensorValues.Length; + int tensorSize = _tensorValues.Length; for (int q = 0; q < derivativeOrders.Length; q++) { int[] derivOrder = derivativeOrders[q]; - double[] current = TensorValues; + double[] current = _tensorValues; int totalSize = tensorSize; for (int d = NumDimensions - 1; d >= 0; d--) { int deriv = derivOrder[d]; - int lastDim = NNodes[d]; + int lastDim = _nNodes[d]; int leadSize = totalSize / lastDim; if (deriv > 0) { - double[] dTFlat = DiffMatricesTFlat![d]; + double[] dTFlat = _diffMatricesTFlat![d]; for (int o = 0; o < deriv; o++) current = BarycentricKernel.MatmulLastAxisMatrixFlat(current, leadSize, lastDim, dTFlat, lastDim); } @@ -702,14 +774,14 @@ public double[] VectorizedEvalMulti(double[] point, int[][] derivativeOrders) /// Per-dimension last-coefficient magnitudes, one entry per dim. public double[] ErrorEstimatePerDim() { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException("Call Build() first"); var perDim = new double[NumDimensions]; for (int d = 0; d < NumDimensions; d++) { double maxErrThisDim = 0.0; - int[] otherShape = NNodes.Where((_, i) => i != d).ToArray(); + int[] otherShape = _nNodes.Where((_, i) => i != d).ToArray(); int otherTotal = TensorShape.RequireArrayLength( TensorShape.CheckedProduct(otherShape, nameof(ErrorEstimatePerDim)), nameof(ErrorEstimatePerDim), @@ -717,7 +789,7 @@ public double[] ErrorEstimatePerDim() for (int otherFlat = 0; otherFlat < otherTotal; otherFlat++) { - double[] values1d = Extract1DSlice(TensorValues, NNodes, d, otherFlat, otherShape); + double[] values1d = Extract1DSlice(_tensorValues, _nNodes, d, otherFlat, otherShape); double[] coeffs = BarycentricKernel.ChebyshevCoefficients1D(values1d); double lastCoeff = Math.Abs(coeffs[^1]); if (lastCoeff > maxErrThisDim) @@ -765,7 +837,7 @@ public static int GetOptimalN1( function, 1, new[] { new[] { domain.lo, domain.hi } }, nNodes: null, errorThreshold: errorThreshold, maxN: maxN); cheb.Build(verbose: false); - return cheb.NNodes[0]; + return cheb._nNodes[0]; } /// @@ -789,7 +861,7 @@ public static double[] ChebyshevCoefficients1D(double[] values) /// portable .pcb format readable by C/Rust/Julia consumers. public void Save(string path, string format = "json") { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException( "Cannot save an unbuilt interpolant. Call Build() first."); @@ -813,13 +885,13 @@ private void SaveJson(string path) var state = new SerializationState { NumDimensions = NumDimensions, - Domain = Domain, - NNodes = NNodes, + Domain = _domain, + NNodes = _nNodes, MaxDerivativeOrder = MaxDerivativeOrder, - NodeArrays = NodeArrays, - TensorValues = TensorValues!, - Weights = Weights!, - DiffMatrices = DiffMatrices!.Select(Flatten2D).ToArray(), + NodeArrays = _nodeArrays, + TensorValues = _tensorValues!, + Weights = _weights!, + DiffMatrices = _diffMatrices!.Select(Flatten2D).ToArray(), BuildTime = BuildTime, NEvaluations = NEvaluations, OriginalNNodes = OriginalNNodes, @@ -843,7 +915,7 @@ private void SaveBinary(string path) using var fs = File.Create(path); using var w = new BinaryWriter(fs); Internal.PcbFormat.WriteHeader(w, Internal.PcbFormat.ClassTagApproximation); - Internal.PcbFormat.WriteApproximationBody(w, Domain, NNodes, TensorValues!); + Internal.PcbFormat.WriteApproximationBody(w, _domain, _nNodes, _tensorValues!); } /// @@ -900,23 +972,23 @@ private static ChebyshevApproximation LoadJson(string path) { Function = null, NumDimensions = state.NumDimensions, - Domain = state.Domain, - NNodes = state.NNodes, + _domain = state.Domain, + _nNodes = state.NNodes, MaxDerivativeOrder = state.MaxDerivativeOrder ?? 2, - NodeArrays = state.NodeArrays, - TensorValues = state.TensorValues, - Weights = state.Weights, + _nodeArrays = state.NodeArrays, + _tensorValues = state.TensorValues, + _weights = state.Weights, BuildTime = state.BuildTime, NEvaluations = state.NEvaluations, _cachedErrorEstimate = null, }; // Reconstruct diff matrices from flat arrays - obj.DiffMatrices = new double[state.NumDimensions][,]; + obj._diffMatrices = new double[state.NumDimensions][,]; for (int d = 0; d < state.NumDimensions; d++) { int n = state.NNodes[d]; - obj.DiffMatrices[d] = Unflatten2D(state.DiffMatrices[d], n, n); + obj._diffMatrices[d] = Unflatten2D(state.DiffMatrices[d], n, n); } obj.PrecomputeTransposedDiffMatrices(); @@ -925,7 +997,7 @@ private static ChebyshevApproximation LoadJson(string path) if (state.OriginalNNodes != null) obj.OriginalNNodes = state.OriginalNNodes; else - obj.OriginalNNodes = obj.NNodes.Select(n => (int?)n).ToArray(); + obj.OriginalNNodes = obj._nNodes.Select(n => (int?)n).ToArray(); obj.ErrorThreshold = state.ErrorThreshold; obj.MaxN = state.MaxN ?? 64; @@ -1135,8 +1207,8 @@ public static ChebyshevApproximation FromValues( { Function = null, NumDimensions = numDimensions, - Domain = domain.Select(d => (double[])d.Clone()).ToArray(), - NNodes = (int[])nNodes.Clone(), + _domain = domain.Select(d => (double[])d.Clone()).ToArray(), + _nNodes = (int[])nNodes.Clone(), MaxDerivativeOrder = maxDerivativeOrder, BuildTime = 0.0, NEvaluations = 0, @@ -1144,21 +1216,21 @@ public static ChebyshevApproximation FromValues( }; // Chebyshev nodes - obj.NodeArrays = new double[numDimensions][]; + obj._nodeArrays = new double[numDimensions][]; for (int d = 0; d < numDimensions; d++) - obj.NodeArrays[d] = BarycentricKernel.MakeNodesForDim(domain[d][0], domain[d][1], nNodes[d]); + obj._nodeArrays[d] = BarycentricKernel.MakeNodesForDim(domain[d][0], domain[d][1], nNodes[d]); - obj.TensorValues = (double[])tensorValues.Clone(); + obj._tensorValues = (double[])tensorValues.Clone(); // Barycentric weights - obj.Weights = new double[numDimensions][]; + obj._weights = new double[numDimensions][]; for (int d = 0; d < numDimensions; d++) - obj.Weights[d] = BarycentricKernel.ComputeBarycentricWeights(obj.NodeArrays[d]); + obj._weights[d] = BarycentricKernel.ComputeBarycentricWeights(obj._nodeArrays[d]); // Differentiation matrices - obj.DiffMatrices = new double[numDimensions][,]; + obj._diffMatrices = new double[numDimensions][,]; for (int d = 0; d < numDimensions; d++) - obj.DiffMatrices[d] = BarycentricKernel.ComputeDifferentiationMatrix(obj.NodeArrays[d], obj.Weights[d]); + obj._diffMatrices[d] = BarycentricKernel.ComputeDifferentiationMatrix(obj._nodeArrays[d], obj._weights[d]); obj.PrecomputeTransposedDiffMatrices(); @@ -1178,14 +1250,14 @@ internal static ChebyshevApproximation FromGrid(ChebyshevApproximation source, d { Function = null, NumDimensions = source.NumDimensions, - Domain = source.Domain.Select(d => (double[])d.Clone()).ToArray(), - NNodes = (int[])source.NNodes.Clone(), + _domain = source._domain.Select(d => (double[])d.Clone()).ToArray(), + _nNodes = (int[])source._nNodes.Clone(), MaxDerivativeOrder = source.MaxDerivativeOrder, - NodeArrays = Internal.CloneHelpers.DeepCopy(source.NodeArrays)!, - Weights = Internal.CloneHelpers.DeepCopy(source.Weights), - DiffMatrices = Internal.CloneHelpers.DeepCopy(source.DiffMatrices), - DiffMatricesTFlat = Internal.CloneHelpers.DeepCopy(source.DiffMatricesTFlat), - TensorValues = tensorValues, + _nodeArrays = Internal.CloneHelpers.DeepCopy(source._nodeArrays)!, + _weights = Internal.CloneHelpers.DeepCopy(source._weights), + _diffMatrices = Internal.CloneHelpers.DeepCopy(source._diffMatrices), + _diffMatricesTFlat = Internal.CloneHelpers.DeepCopy(source._diffMatricesTFlat), + _tensorValues = tensorValues, BuildTime = 0.0, NEvaluations = 0, _cachedErrorEstimate = null, @@ -1202,18 +1274,18 @@ internal static ChebyshevApproximation FromGrid(ChebyshevApproximation source, d /// public ChebyshevApproximation Extrude(params (int dimIndex, double[] bounds, int nNodes)[] extrudeParams) { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException("Call Build() first"); var sorted = ExtrudeSlice.NormalizeExtrusionParams(extrudeParams, NumDimensions); - double[] tensor = (double[])TensorValues.Clone(); - int[] shape = (int[])NNodes.Clone(); - var nodes = NodeArrays.ToList(); - var weights = Weights!.ToList(); - var diffMats = DiffMatrices!.ToList(); - var domain = Domain.Select(d => (double[])d.Clone()).ToList(); - var nNodes = NNodes.ToList(); + double[] tensor = (double[])_tensorValues.Clone(); + int[] shape = (int[])_nNodes.Clone(); + var nodes = _nodeArrays.ToList(); + var weights = _weights!.ToList(); + var diffMats = _diffMatrices!.ToList(); + var domain = _domain.Select(d => (double[])d.Clone()).ToList(); + var nNodes = _nNodes.ToList(); foreach (var (dimIdx, bounds, n) in sorted) { @@ -1238,13 +1310,13 @@ public ChebyshevApproximation Extrude(params (int dimIndex, double[] bounds, int { Function = null, NumDimensions = newNdim, - Domain = domain.ToArray(), - NNodes = nNodes.ToArray(), + _domain = domain.ToArray(), + _nNodes = nNodes.ToArray(), MaxDerivativeOrder = MaxDerivativeOrder, - NodeArrays = nodes.ToArray(), - Weights = weights.ToArray(), - DiffMatrices = diffMats.ToArray(), - TensorValues = tensor, + _nodeArrays = nodes.ToArray(), + _weights = weights.ToArray(), + _diffMatrices = diffMats.ToArray(), + _tensorValues = tensor, BuildTime = 0.0, NEvaluations = 0, _cachedErrorEstimate = null, @@ -1259,7 +1331,7 @@ public ChebyshevApproximation Extrude(params (int dimIndex, double[] bounds, int /// public ChebyshevApproximation Slice(params (int dimIndex, double value)[] sliceParams) { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException("Call Build() first"); var sorted = ExtrudeSlice.NormalizeSlicingParams(sliceParams, NumDimensions); @@ -1267,20 +1339,20 @@ public ChebyshevApproximation Slice(params (int dimIndex, double value)[] sliceP // Validate values within domain foreach (var (dimIdx, value) in sorted) { - double lo = Domain[dimIdx][0]; - double hi = Domain[dimIdx][1]; + double lo = _domain[dimIdx][0]; + double hi = _domain[dimIdx][1]; if (value < lo || value > hi) throw new ArgumentException( $"Slice value {value} for dim {dimIdx} is outside domain [{lo}, {hi}]"); } - double[] tensor = (double[])TensorValues.Clone(); - int[] shape = (int[])NNodes.Clone(); - var nodes = NodeArrays.ToList(); - var weights = Weights!.ToList(); - var diffMats = DiffMatrices!.ToList(); - var domain = Domain.Select(d => (double[])d.Clone()).ToList(); - var nNodes = NNodes.ToList(); + double[] tensor = (double[])_tensorValues.Clone(); + int[] shape = (int[])_nNodes.Clone(); + var nodes = _nodeArrays.ToList(); + var weights = _weights!.ToList(); + var diffMats = _diffMatrices!.ToList(); + var domain = _domain.Select(d => (double[])d.Clone()).ToList(); + var nNodes = _nNodes.ToList(); foreach (var (dimIdx, value) in sorted) { @@ -1303,13 +1375,13 @@ public ChebyshevApproximation Slice(params (int dimIndex, double value)[] sliceP { Function = null, NumDimensions = newNdim, - Domain = domain.ToArray(), - NNodes = nNodes.ToArray(), + _domain = domain.ToArray(), + _nNodes = nNodes.ToArray(), MaxDerivativeOrder = MaxDerivativeOrder, - NodeArrays = nodes.ToArray(), - Weights = weights.ToArray(), - DiffMatrices = diffMats.ToArray(), - TensorValues = tensor, + _nodeArrays = nodes.ToArray(), + _weights = weights.ToArray(), + _diffMatrices = diffMats.ToArray(), + _tensorValues = tensor, BuildTime = 0.0, NEvaluations = 0, _cachedErrorEstimate = null, @@ -1331,7 +1403,7 @@ public ChebyshevApproximation Slice(params (int dimIndex, double value)[] sliceP /// Scalar if all dims integrated, otherwise a lower-dimensional interpolant. public object Integrate(int[]? dims = null, (double lo, double hi)[]? bounds = null) { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException("Call Build() first"); if (dims == null) @@ -1345,18 +1417,18 @@ public object Integrate(int[]? dims = null, (double lo, double hi)[]? bounds = n throw new ArgumentException($"dim {d} out of range [0, {NumDimensions - 1}]"); } - var perDimBounds = Calculus.NormalizeBounds(sortedDims, bounds, Domain); + var perDimBounds = Calculus.NormalizeBounds(sortedDims, bounds, _domain); var dimToIdx = new Dictionary(); for (int i = 0; i < sortedDims.Length; i++) dimToIdx[sortedDims[i]] = i; - double[] tensor = (double[])TensorValues.Clone(); - int[] shape = (int[])NNodes.Clone(); - var nodes = NodeArrays.ToList(); - var wts = Weights!.ToList(); - var diffMats = DiffMatrices!.ToList(); - var domain = Domain.Select(d => (double[])d.Clone()).ToList(); - var nNodes = NNodes.ToList(); + double[] tensor = (double[])_tensorValues.Clone(); + int[] shape = (int[])_nNodes.Clone(); + var nodes = _nodeArrays.ToList(); + var wts = _weights!.ToList(); + var diffMats = _diffMatrices!.ToList(); + var domain = _domain.Select(d => (double[])d.Clone()).ToList(); + var nNodes = _nNodes.ToList(); // Process dimensions in descending order foreach (int d in sortedDims.OrderByDescending(x => x)) @@ -1405,13 +1477,13 @@ public object Integrate(int[]? dims = null, (double lo, double hi)[]? bounds = n { Function = null, NumDimensions = newNdim, - Domain = domain.ToArray(), - NNodes = nNodes.ToArray(), + _domain = domain.ToArray(), + _nNodes = nNodes.ToArray(), MaxDerivativeOrder = MaxDerivativeOrder, - NodeArrays = nodes.ToArray(), - Weights = wts.ToArray(), - DiffMatrices = diffMats.ToArray(), - TensorValues = tensor, + _nodeArrays = nodes.ToArray(), + _weights = wts.ToArray(), + _diffMatrices = diffMats.ToArray(), + _tensorValues = tensor, BuildTime = 0.0, NEvaluations = 0, _cachedErrorEstimate = null, @@ -1426,13 +1498,13 @@ public object Integrate(int[]? dims = null, (double lo, double hi)[]? bounds = n /// public double[] Roots(int? dim = null, Dictionary? fixedDims = null) { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException("Call Build() first"); - var (targetDim, sliceParams) = Calculus.ValidateCalculusArgs(NumDimensions, dim, fixedDims, Domain); + var (targetDim, sliceParams) = Calculus.ValidateCalculusArgs(NumDimensions, dim, fixedDims, _domain); ChebyshevApproximation sliced = sliceParams.Length > 0 ? Slice(sliceParams) : this; - return Calculus.Roots1D(sliced.TensorValues!, sliced.Domain[0]); + return Calculus.Roots1D(sliced._tensorValues!, sliced._domain[0]); } /// @@ -1440,15 +1512,15 @@ public double[] Roots(int? dim = null, Dictionary? fixedDims = null /// public (double value, double location) Minimize(int? dim = null, Dictionary? fixedDims = null) { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException("Call Build() first"); - var (targetDim, sliceParams) = Calculus.ValidateCalculusArgs(NumDimensions, dim, fixedDims, Domain); + var (targetDim, sliceParams) = Calculus.ValidateCalculusArgs(NumDimensions, dim, fixedDims, _domain); ChebyshevApproximation sliced = sliceParams.Length > 0 ? Slice(sliceParams) : this; return Calculus.Optimize1D( - sliced.TensorValues!, sliced.NodeArrays[0], sliced.Weights![0], - sliced.DiffMatrices![0], sliced.Domain[0], "min"); + sliced._tensorValues!, sliced._nodeArrays[0], sliced._weights![0], + sliced._diffMatrices![0], sliced._domain[0], "min"); } /// @@ -1456,15 +1528,15 @@ public double[] Roots(int? dim = null, Dictionary? fixedDims = null /// public (double value, double location) Maximize(int? dim = null, Dictionary? fixedDims = null) { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException("Call Build() first"); - var (targetDim, sliceParams) = Calculus.ValidateCalculusArgs(NumDimensions, dim, fixedDims, Domain); + var (targetDim, sliceParams) = Calculus.ValidateCalculusArgs(NumDimensions, dim, fixedDims, _domain); ChebyshevApproximation sliced = sliceParams.Length > 0 ? Slice(sliceParams) : this; return Calculus.Optimize1D( - sliced.TensorValues!, sliced.NodeArrays[0], sliced.Weights![0], - sliced.DiffMatrices![0], sliced.Domain[0], "max"); + sliced._tensorValues!, sliced._nodeArrays[0], sliced._weights![0], + sliced._diffMatrices![0], sliced._domain[0], "max"); } // ------------------------------------------------------------------ @@ -1477,9 +1549,9 @@ public double[] Roots(int? dim = null, Dictionary? fixedDims = null if (a.GetType() != b.GetType()) throw new InvalidOperationException("Cannot combine different types"); Algebra.CheckCompatible(a, b); - double[] newValues = new double[a.TensorValues!.Length]; + double[] newValues = new double[a._tensorValues!.Length]; for (int i = 0; i < newValues.Length; i++) - newValues[i] = a.TensorValues[i] + b.TensorValues![i]; + newValues[i] = a._tensorValues[i] + b._tensorValues![i]; return FromGrid(a, newValues); } @@ -1489,21 +1561,21 @@ public double[] Roots(int? dim = null, Dictionary? fixedDims = null if (a.GetType() != b.GetType()) throw new InvalidOperationException("Cannot combine different types"); Algebra.CheckCompatible(a, b); - double[] newValues = new double[a.TensorValues!.Length]; + double[] newValues = new double[a._tensorValues!.Length]; for (int i = 0; i < newValues.Length; i++) - newValues[i] = a.TensorValues[i] - b.TensorValues![i]; + newValues[i] = a._tensorValues[i] - b._tensorValues![i]; return FromGrid(a, newValues); } /// Multiply interpolant by a scalar. public static ChebyshevApproximation operator *(ChebyshevApproximation a, double scalar) { - if (a.TensorValues == null) + if (a._tensorValues == null) throw new InvalidOperationException("Operand is not built. Call Build() first."); - double[] newValues = new double[a.TensorValues!.Length]; + double[] newValues = new double[a._tensorValues!.Length]; for (int i = 0; i < newValues.Length; i++) - newValues[i] = a.TensorValues[i] * scalar; + newValues[i] = a._tensorValues[i] * scalar; return FromGrid(a, newValues); } @@ -1532,21 +1604,21 @@ public double[] Roots(int? dim = null, Dictionary? fixedDims = null /// public override string ToString() { - bool built = TensorValues != null; - long totalNodes = TensorShape.CheckedProduct(NNodes, nameof(ToString)); + bool built = _tensorValues != null; + long totalNodes = TensorShape.CheckedProduct(_nNodes, nameof(ToString)); string status = built ? "built" : "not built"; int maxDisplay = 6; string nodesStr, domainStr; if (NumDimensions > maxDisplay) { - nodesStr = "[" + string.Join(", ", NNodes.Take(maxDisplay)) + ", ...]"; - domainStr = string.Join(" x ", Domain.Take(maxDisplay).Select(d => $"[{d[0]}, {d[1]}]")) + " x ..."; + nodesStr = "[" + string.Join(", ", _nNodes.Take(maxDisplay)) + ", ...]"; + domainStr = string.Join(" x ", _domain.Take(maxDisplay).Select(d => $"[{d[0]}, {d[1]}]")) + " x ..."; } else { - nodesStr = "[" + string.Join(", ", NNodes) + "]"; - domainStr = string.Join(" x ", Domain.Select(d => $"[{d[0]}, {d[1]}]")); + nodesStr = "[" + string.Join(", ", _nNodes) + "]"; + domainStr = string.Join(" x ", _domain.Select(d => $"[{d[0]}, {d[1]}]")); } var sb = new StringBuilder(); @@ -1569,8 +1641,8 @@ public override string ToString() /// public string ToReprString() { - bool built = TensorValues != null; - return $"ChebyshevApproximation(dims={NumDimensions}, nodes=[{string.Join(", ", NNodes)}], built={built})"; + bool built = _tensorValues != null; + return $"ChebyshevApproximation(dims={NumDimensions}, nodes=[{string.Join(", ", _nNodes)}], built={built})"; } // ------------------------------------------------------------------ @@ -1630,16 +1702,16 @@ private static void ValidateNodeCounts(string caller, int[] nNodes) /// /// Pre-compute transposed diff matrices as flat arrays for BLAS GEMM. - /// Called after DiffMatrices is set in Build, FromValues, Load, Extrude, Slice, and Integrate. + /// Called after differentiation matrices are set in Build, FromValues, Load, Extrude, Slice, and Integrate. /// internal void PrecomputeTransposedDiffMatrices() { - if (DiffMatrices == null) return; - DiffMatricesTFlat = new double[DiffMatrices.Length][]; - for (int d = 0; d < DiffMatrices.Length; d++) + if (_diffMatrices == null) return; + _diffMatricesTFlat = new double[_diffMatrices.Length][]; + for (int d = 0; d < _diffMatrices.Length; d++) { - int rows = DiffMatrices[d].GetLength(0); - int cols = DiffMatrices[d].GetLength(1); + int rows = _diffMatrices[d].GetLength(0); + int cols = _diffMatrices[d].GetLength(1); // Transpose and flatten in one pass (row-major) int flatLength = TensorShape.RequireArrayLength( TensorShape.CheckedProduct(new[] { rows, cols }, nameof(PrecomputeTransposedDiffMatrices)), @@ -1648,8 +1720,8 @@ internal void PrecomputeTransposedDiffMatrices() var flat = new double[flatLength]; for (int i = 0; i < rows; i++) for (int j = 0; j < cols; j++) - flat[j * rows + i] = DiffMatrices[d][i, j]; - DiffMatricesTFlat[d] = flat; + flat[j * rows + i] = _diffMatrices[d][i, j]; + _diffMatricesTFlat[d] = flat; } } @@ -1738,7 +1810,7 @@ internal static double[] Flatten2D(double[,] matrix) public string GetConstructorType() => _constructorType; /// Per-dimension Chebyshev node counts actually used. After auto-N construction, these are the resolved values. - public int[] GetUsedNs() => (int[])NNodes.Clone(); + public int[] GetUsedNs() => (int[])_nNodes.Clone(); /// Maximum derivative order this approximation supports. public int GetMaxDerivativeOrder() => MaxDerivativeOrder; @@ -1760,9 +1832,9 @@ internal static double[] Flatten2D(double[,] matrix) public int GetNumEvaluationPoints() { return TensorShape.RequireArrayLength( - TensorShape.CheckedProduct(NNodes, nameof(GetNumEvaluationPoints)), + TensorShape.CheckedProduct(_nNodes, nameof(GetNumEvaluationPoints)), nameof(GetNumEvaluationPoints), - NNodes); + _nNodes); } /// @@ -1788,12 +1860,12 @@ public double[] GetEvaluationPoints() int rem = flat; for (int d = ndim - 1; d >= 0; d--) { - indices[d] = rem % NNodes[d]; - rem /= NNodes[d]; + indices[d] = rem % _nNodes[d]; + rem /= _nNodes[d]; } for (int d = 0; d < ndim; d++) { - points[flat * ndim + d] = NodeArrays[d][indices[d]]; + points[flat * ndim + d] = _nodeArrays[d][indices[d]]; } } @@ -1819,41 +1891,41 @@ public double[] GetEvaluationPoints() public void SetOriginalFunctionValues(double[] values) { ArgumentNullException.ThrowIfNull(values); - if (_isConstructionFinished || TensorValues != null) + if (_isConstructionFinished || _tensorValues != null) throw new InvalidOperationException( "interpolant is already constructed; SetOriginalFunctionValues is for unconstructed deferred objects"); int expected = TensorShape.RequireArrayLength( - TensorShape.CheckedProduct(NNodes, nameof(SetOriginalFunctionValues)), + TensorShape.CheckedProduct(_nNodes, nameof(SetOriginalFunctionValues)), nameof(SetOriginalFunctionValues), - NNodes); + _nNodes); if (values.Length != expected) throw new ArgumentException( - $"values has {values.Length} entries, expected {expected} for nNodes=[{string.Join(",", NNodes)}]"); + $"values has {values.Length} entries, expected {expected} for nNodes=[{string.Join(",", _nNodes)}]"); for (int i = 0; i < values.Length; i++) { if (!double.IsFinite(values[i])) throw new ArgumentException("values contains NaN or Inf", nameof(values)); } - // Materialize NodeArrays now if deferred (NodeArrays is empty when deferBuild was true). - if (NodeArrays.Length == 0) + // Materialize node arrays now if deferred. + if (_nodeArrays.Length == 0) { - NodeArrays = new double[NumDimensions][]; + _nodeArrays = new double[NumDimensions][]; for (int d = 0; d < NumDimensions; d++) - NodeArrays[d] = BarycentricKernel.MakeNodesForDim(Domain[d][0], Domain[d][1], NNodes[d]); + _nodeArrays[d] = BarycentricKernel.MakeNodesForDim(_domain[d][0], _domain[d][1], _nNodes[d]); } // Mirror FromValues precomputation (bit-identical). - TensorValues = (double[])values.Clone(); + _tensorValues = (double[])values.Clone(); - Weights = new double[NumDimensions][]; + _weights = new double[NumDimensions][]; for (int d = 0; d < NumDimensions; d++) - Weights[d] = BarycentricKernel.ComputeBarycentricWeights(NodeArrays[d]); + _weights[d] = BarycentricKernel.ComputeBarycentricWeights(_nodeArrays[d]); - DiffMatrices = new double[NumDimensions][,]; + _diffMatrices = new double[NumDimensions][,]; for (int d = 0; d < NumDimensions; d++) - DiffMatrices[d] = BarycentricKernel.ComputeDifferentiationMatrix(NodeArrays[d], Weights[d]); + _diffMatrices[d] = BarycentricKernel.ComputeDifferentiationMatrix(_nodeArrays[d], _weights[d]); PrecomputeTransposedDiffMatrices(); @@ -1913,13 +1985,13 @@ public ChebyshevApproximation Clone() { var copy = new ChebyshevApproximation(); copy.NumDimensions = NumDimensions; - copy.NNodes = Internal.CloneHelpers.DeepCopy(NNodes)!; - copy.Domain = Internal.CloneHelpers.DeepCopy(Domain)!; - copy.NodeArrays = Internal.CloneHelpers.DeepCopy(NodeArrays)!; - copy.TensorValues = Internal.CloneHelpers.DeepCopy(TensorValues); - copy.Weights = Internal.CloneHelpers.DeepCopy(Weights); - copy.DiffMatrices = Internal.CloneHelpers.DeepCopy(DiffMatrices); - copy.DiffMatricesTFlat = Internal.CloneHelpers.DeepCopy(DiffMatricesTFlat); + copy.NNodes = Internal.CloneHelpers.DeepCopy(_nNodes)!; + copy.Domain = Internal.CloneHelpers.DeepCopy(_domain)!; + copy.NodeArrays = Internal.CloneHelpers.DeepCopy(_nodeArrays)!; + copy.TensorValues = Internal.CloneHelpers.DeepCopy(_tensorValues); + copy.Weights = Internal.CloneHelpers.DeepCopy(_weights); + copy.DiffMatrices = Internal.CloneHelpers.DeepCopy(_diffMatrices); + copy.DiffMatricesTFlat = Internal.CloneHelpers.DeepCopy(_diffMatricesTFlat); copy.MaxDerivativeOrder = MaxDerivativeOrder; copy.MaxN = MaxN; copy.ErrorThreshold = ErrorThreshold; @@ -1954,11 +2026,11 @@ public ChebyshevApproximation Clone() /// If has not been called. public SobolResult SobolIndices() { - if (TensorValues == null) + if (_tensorValues == null) throw new InvalidOperationException( "SobolIndices requires a built ChebyshevApproximation. Call Build() first."); - var coeffs = Internal.Sensitivity.ChebyshevCoefficientsND(TensorValues, NNodes); - return Internal.Sensitivity.ComputeSobolFromCoeffs(coeffs, NNodes); + var coeffs = Internal.Sensitivity.ChebyshevCoefficientsND(_tensorValues, _nNodes); + return Internal.Sensitivity.ComputeSobolFromCoeffs(coeffs, _nNodes); } // ------------------------------------------------------------------ diff --git a/src/ChebyshevSharp/ChebyshevSpline.cs b/src/ChebyshevSharp/ChebyshevSpline.cs index a19d59d..22c9d94 100644 --- a/src/ChebyshevSharp/ChebyshevSpline.cs +++ b/src/ChebyshevSharp/ChebyshevSpline.cs @@ -1020,22 +1020,22 @@ private static ChebyshevSpline LoadJson(string path) { Function = null, NumDimensions = ps.NumDimensions, - Domain = ps.Domain, - NNodes = ps.NNodes, + DomainStorage = ps.Domain, + NNodesStorage = ps.NNodes, MaxDerivativeOrder = ps.MaxDerivativeOrder ?? 2, - NodeArrays = ps.NodeArrays, - TensorValues = ps.TensorValues, - Weights = ps.Weights, + NodeArraysStorage = ps.NodeArrays, + TensorValuesStorage = ps.TensorValues, + WeightsStorage = ps.Weights, BuildTime = ps.BuildTime, NEvaluations = ps.NEvaluations, }; // Reconstruct diff matrices - piece.DiffMatrices = new double[ps.NumDimensions][,]; + piece.DiffMatricesStorage = new double[ps.NumDimensions][,]; for (int d = 0; d < ps.NumDimensions; d++) { int n = ps.NNodes[d]; - piece.DiffMatrices[d] = ChebyshevApproximation.Unflatten2D(ps.DiffMatrices[d], n, n); + piece.DiffMatricesStorage[d] = ChebyshevApproximation.Unflatten2D(ps.DiffMatrices[d], n, n); } piece.PrecomputeTransposedDiffMatrices(); return piece; diff --git a/src/ChebyshevSharp/Internal/AdaptiveBuild.cs b/src/ChebyshevSharp/Internal/AdaptiveBuild.cs index 7765d50..a13750f 100644 --- a/src/ChebyshevSharp/Internal/AdaptiveBuild.cs +++ b/src/ChebyshevSharp/Internal/AdaptiveBuild.cs @@ -48,11 +48,11 @@ public static void RunDoublingLoop(ChebyshevApproximation approx, bool verbose) while (true) { // Apply current grid - approx.NNodes = (int[])current.Clone(); - approx.NodeArrays = new double[numDim][]; + approx.NNodesStorage = (int[])current.Clone(); + approx.NodeArraysStorage = new double[numDim][]; for (int d = 0; d < numDim; d++) - approx.NodeArrays[d] = BarycentricKernel.MakeNodesForDim( - approx.Domain[d][0], approx.Domain[d][1], current[d]); + approx.NodeArraysStorage[d] = BarycentricKernel.MakeNodesForDim( + approx.DomainStorage[d][0], approx.DomainStorage[d][1], current[d]); approx.BuildFixedGrid(verbose: false); totalEvals += approx.NEvaluations; @@ -149,7 +149,7 @@ private static (double[] PerDim, int Evaluations) ValidationErrorPerAutoDim( validationShape); double[] probeNodes = BarycentricKernel.MakeNodesForDim( - approx.Domain[dim][0], approx.Domain[dim][1], probeN); + approx.DomainStorage[dim][0], approx.DomainStorage[dim][1], probeN); double maxErr = 0.0; for (int flat = 0; flat < total; flat++) @@ -160,7 +160,7 @@ private static (double[] PerDim, int Evaluations) ValidationErrorPerAutoDim( { int idx = rem % validationShape[d]; rem /= validationShape[d]; - point[d] = d == dim ? probeNodes[idx] : approx.NodeArrays[d][idx]; + point[d] = d == dim ? probeNodes[idx] : approx.NodeArraysStorage[d][idx]; } double expected = function(point, approx.AdditionalData); diff --git a/src/ChebyshevSharp/Internal/Algebra.cs b/src/ChebyshevSharp/Internal/Algebra.cs index ed81659..f7e1914 100644 --- a/src/ChebyshevSharp/Internal/Algebra.cs +++ b/src/ChebyshevSharp/Internal/Algebra.cs @@ -39,27 +39,31 @@ internal static void CheckCompatible(ChebyshevApproximation a, ChebyshevApproxim $"Cannot combine {a.GetType().Name} with {b.GetType().Name}; " + "operands must be the same type."); - if (a.TensorValues == null) + if (a.TensorValuesStorage == null) throw new InvalidOperationException("Left operand is not built. Call Build() first."); - if (b.TensorValues == null) + if (b.TensorValuesStorage == null) throw new InvalidOperationException("Right operand is not built. Call Build() first."); if (a.NumDimensions != b.NumDimensions) throw new ArgumentException( $"Dimension mismatch: {a.NumDimensions} vs {b.NumDimensions}"); - if (!a.NNodes.SequenceEqual(b.NNodes)) + int[] aNNodes = a.NNodesStorage; + int[] bNNodes = b.NNodesStorage; + if (!aNNodes.SequenceEqual(bNNodes)) throw new ArgumentException( - $"Node count mismatch: [{string.Join(", ", a.NNodes)}] vs [{string.Join(", ", b.NNodes)}]"); + $"Node count mismatch: [{string.Join(", ", aNNodes)}] vs [{string.Join(", ", bNNodes)}]"); // v0.21.1: numerical comparison on Domain[d] (was SequenceEqual = exact). // Tolerates sub-ULP drift between equivalent allocations. + double[][] aDomain = a.DomainStorage; + double[][] bDomain = b.DomainStorage; for (int d = 0; d < a.NumDimensions; d++) { - if (!DoublesAllClose(a.Domain[d], b.Domain[d])) + if (!DoublesAllClose(aDomain[d], bDomain[d])) throw new ArgumentException( $"Domain mismatch at dim {d}: " + - $"[{a.Domain[d][0]}, {a.Domain[d][1]}] vs [{b.Domain[d][0]}, {b.Domain[d][1]}]"); + $"[{aDomain[d][0]}, {aDomain[d][1]}] vs [{bDomain[d][0]}, {bDomain[d][1]}]"); } if (a.MaxDerivativeOrder != b.MaxDerivativeOrder) diff --git a/tests/ChebyshevSharp.Tests/ApproxPublicStateOwnershipTests.cs b/tests/ChebyshevSharp.Tests/ApproxPublicStateOwnershipTests.cs new file mode 100644 index 0000000..dc24451 --- /dev/null +++ b/tests/ChebyshevSharp.Tests/ApproxPublicStateOwnershipTests.cs @@ -0,0 +1,182 @@ +using Xunit; + +namespace ChebyshevSharp.Tests; + +public class ApproxPublicStateOwnershipTests +{ + private static ChebyshevApproximation BuildApproximation() + { + var approx = new ChebyshevApproximation( + (p, _) => p[0] * p[0] + 2.0 * p[1], + numDimensions: 2, + domain: new[] { new[] { -1.0, 1.0 }, new[] { -1.0, 1.0 } }, + nNodes: new[] { 7, 7 }); + approx.Build(verbose: false); + return approx; + } + + [Fact] + public void Domain_and_NNodes_properties_return_snapshots() + { + var approx = BuildApproximation(); + + double[][] domain = approx.Domain; + int[] nNodes = approx.NNodes; + domain[0][0] = 10.0; + nNodes[0] = 999; + + Assert.Equal(-1.0, approx.Domain[0][0]); + Assert.Equal(new[] { 7, 7 }, approx.NNodes); + Assert.Throws(() => approx.Eval(new[] { 2.0, 0.0 })); + } + + [Fact] + public void NodeArrays_and_TensorValues_properties_return_snapshots() + { + var approx = BuildApproximation(); + double before = approx.Eval(new[] { 0.25, -0.5 }); + + double[][] nodes = approx.NodeArrays; + double[] tensor = approx.TensorValues!; + nodes[0][0] = 123.0; + tensor[0] = 123.0; + + double after = approx.Eval(new[] { 0.25, -0.5 }); + Assert.Equal(before, after, precision: 12); + Assert.NotEqual(123.0, approx.NodeArrays[0][0]); + Assert.NotEqual(123.0, approx.TensorValues![0]); + } + + [Fact] + public void Weight_and_diff_matrix_properties_return_snapshots() + { + var approx = BuildApproximation(); + double valueBefore = approx.Eval(new[] { 0.25, -0.5 }); + double derivativeBefore = approx.Eval(new[] { 0.25, -0.5 }, new[] { 1, 0 }); + + double[][] weights = approx.Weights!; + double[][,] diffMatrices = approx.DiffMatrices!; + weights[0][0] = 123.0; + diffMatrices[0][0, 0] = 123.0; + + double valueAfter = approx.Eval(new[] { 0.25, -0.5 }); + double derivativeAfter = approx.Eval(new[] { 0.25, -0.5 }, new[] { 1, 0 }); + Assert.Equal(valueBefore, valueAfter, precision: 12); + Assert.Equal(derivativeBefore, derivativeAfter, precision: 12); + Assert.NotEqual(123.0, approx.Weights![0][0]); + Assert.NotEqual(123.0, approx.DiffMatrices![0][0, 0]); + } + + [Fact] + public void Internal_storage_accessors_remain_live_while_public_properties_snapshot() + { + var approx = new ChebyshevApproximation(); + var domain = new[] { new[] { -1.0, 1.0 } }; + var nNodes = new[] { 3 }; + var nodes = new[] { new[] { -0.5, 0.0, 0.5 } }; + var tensor = new[] { 1.0, 2.0, 3.0 }; + var weights = new[] { new[] { -1.0, 2.0, -1.0 } }; + var diffMatrices = new[] { new double[,] { { 0.0, 1.0 }, { -1.0, 0.0 } } }; + var diffMatricesFlat = new[] { new[] { 0.0, -1.0, 1.0, 0.0 } }; + + approx.Domain = domain; + approx.NNodes = nNodes; + approx.NodeArrays = nodes; + approx.TensorValues = tensor; + approx.Weights = weights; + approx.DiffMatrices = diffMatrices; + approx.DiffMatricesTFlat = diffMatricesFlat; + + Assert.Same(domain, approx.DomainStorage); + Assert.Same(nNodes, approx.NNodesStorage); + Assert.Same(nodes, approx.NodeArraysStorage); + Assert.Same(tensor, approx.TensorValuesStorage); + Assert.Same(weights, approx.WeightsStorage); + Assert.Same(diffMatrices, approx.DiffMatricesStorage); + Assert.Same(diffMatricesFlat, approx.DiffMatricesTFlat); + Assert.NotSame(domain, approx.Domain); + Assert.NotSame(nNodes, approx.NNodes); + Assert.NotSame(nodes, approx.NodeArrays); + Assert.NotSame(tensor, approx.TensorValues); + Assert.NotSame(weights, approx.Weights); + Assert.NotSame(diffMatrices, approx.DiffMatrices); + + var replacementDomain = new[] { new[] { 0.0, 2.0 } }; + var replacementNNodes = new[] { 5 }; + var replacementNodes = new[] { new[] { -1.0, -0.5, 0.0, 0.5, 1.0 } }; + var replacementTensor = new[] { 1.0 }; + var replacementWeights = new[] { new[] { 1.0 } }; + var replacementDiffMatrices = new[] { new double[,] { { 1.0 } } }; + + approx.DomainStorage = replacementDomain; + approx.NNodesStorage = replacementNNodes; + approx.NodeArraysStorage = replacementNodes; + approx.TensorValuesStorage = replacementTensor; + approx.WeightsStorage = replacementWeights; + approx.DiffMatricesStorage = replacementDiffMatrices; + + Assert.Same(replacementDomain, approx.DomainStorage); + Assert.Same(replacementNNodes, approx.NNodesStorage); + Assert.Same(replacementNodes, approx.NodeArraysStorage); + Assert.Same(replacementTensor, approx.TensorValuesStorage); + Assert.Same(replacementWeights, approx.WeightsStorage); + Assert.Same(replacementDiffMatrices, approx.DiffMatricesStorage); + } + + [Fact] + public void Internal_setters_accept_null_without_exposing_mutable_empty_state() + { + var approx = new ChebyshevApproximation(); + + approx.Domain = null!; + approx.NNodes = null!; + approx.NodeArrays = null!; + approx.TensorValues = null; + approx.Weights = null; + approx.DiffMatrices = null; + approx.DiffMatricesTFlat = null; + + Assert.Empty(approx.Domain); + Assert.Empty(approx.NNodes); + Assert.Empty(approx.NodeArrays); + Assert.Null(approx.TensorValues); + Assert.Null(approx.Weights); + Assert.Null(approx.DiffMatrices); + Assert.Null(approx.DiffMatricesTFlat); + + approx.PrecomputeTransposedDiffMatrices(); + + Assert.Null(approx.DiffMatricesTFlat); + } + + [Fact] + public void Nullable_nnodes_constructor_with_explicit_counts_initializes_private_storage() + { + var approx = new ChebyshevApproximation( + (p, _) => p[0] + p[1], + numDimensions: 2, + domain: new[] { new[] { -2.0, 2.0 }, new[] { 0.0, 1.0 } }, + nNodes: new int?[] { 4, 5 }); + + Assert.Equal(new[] { 4, 5 }, approx.NNodes); + Assert.Equal(2, approx.NodeArrays.Length); + Assert.Equal(4, approx.NodeArrays[0].Length); + Assert.Equal(5, approx.NodeArrays[1].Length); + Assert.Same(approx.NNodesStorage, approx.NNodesStorage); + Assert.Same(approx.NodeArraysStorage, approx.NodeArraysStorage); + } + + [Fact] + public void VectorizedEvalBatch_requires_built_interpolant() + { + var approx = new ChebyshevApproximation( + (p, _) => p[0], + numDimensions: 1, + domain: new[] { new[] { -1.0, 1.0 } }, + nNodes: new[] { 3 }, + deferBuild: true); + + Assert.Throws( + () => approx.VectorizedEvalBatch(new[] { new[] { 0.0 } }, new[] { 0 })); + } +}