diff --git a/src/Gemstone.Numeric/Analysis/CurveFit.cs b/src/Gemstone.Numeric/Analysis/CurveFit.cs index 00fd6be05..a5108aa81 100644 --- a/src/Gemstone.Numeric/Analysis/CurveFit.cs +++ b/src/Gemstone.Numeric/Analysis/CurveFit.cs @@ -233,4 +233,55 @@ public static void LeastSquares(double[] zValues, double[] xValues, double[] yVa b = (coeff10 - c * coeff13) / coeff12; a = (coeff00 - b * coeff02 - c * coeff03) / coeff01; } + + /// + /// Uses least squares linear regression to estimate the coefficients a and b + /// from the given (x,y) data points for the equation y = a + bx. + /// + /// x-value array + /// y-value array + /// the out a coefficient + /// the out b coefficient + public static void LeastSquares(double[] xValues, double[] yValues, out double a, out double b) + { + double n = xValues.Length; + + double xSum = 0; + double ySum = 0; + + double xySum = 0; + + double xxSum = 0; + double yySum = 0; + + for (int i = 0; i < n; i++) + { + double x = xValues[i]; + double y = yValues[i]; + + xSum += x; + ySum += y; + + xySum += x * y; + + xxSum += x * x; + yySum += y * y; + } + + // XtX = [n, x; x, xx] + + // inv(XtX) = 1/div * [xx, -x; -x, n] + + // XtY = [y; xy] + + // veta = inv(XtX) * XtY = 1/d [xx*y - x*xy; -x*y + n*xy] + double coeff00 = xxSum*ySum - xSum*xySum; + double coeff01 = n * xySum - xSum * ySum; + + double divisor = 1.0D/(-xSum*xSum + xxSum*n); + + b = coeff00 * divisor; + a = coeff01 * divisor; + } + } diff --git a/src/Gemstone.Numeric/Analysis/DigitalFilter.cs b/src/Gemstone.Numeric/Analysis/DigitalFilter.cs new file mode 100644 index 000000000..e60ee6069 --- /dev/null +++ b/src/Gemstone.Numeric/Analysis/DigitalFilter.cs @@ -0,0 +1,258 @@ +//****************************************************************************************************** +// DigitalFilter.cs - Gbtc +// +// Copyright © 2022, Grid Protection Alliance. All Rights Reserved. +// +// Licensed to the Grid Protection Alliance (GPA) under one or more contributor license agreements. See +// the NOTICE file distributed with this work for additional information regarding copyright ownership. +// The GPA licenses this file to you under the MIT License (MIT), the "License"; you may not use this +// file except in compliance with the License. You may obtain a copy of the License at: +// +// http://opensource.org/licenses/MIT +// +// Unless agreed to in writing, the subject software distributed under the License is distributed on an +// "AS-IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. Refer to the +// License for the specific language governing permissions and limitations. +// +// Code Modification History: +// ---------------------------------------------------------------------------------------------------- +// 05/10/2025 - C. Lackner +// Generated original version of source code. +// +//****************************************************************************************************** + +using System; +using System.Collections.Generic; +using System.Linq; + +namespace Gemstone.Numeric.Analysis; + +/// +/// A class that implements a digital filter. +/// +public class DigitalFilter +{ + private List m_b; // Numerator coefficients + private List m_a; // Denominator coefficients + + private DigitalFilter(double[] b, double[] a, int size) + { + m_b = new List(size); + m_b.AddRange(b); + + m_a = new List(size); + m_a.AddRange(a); + + Resize(m_b, size); + Resize(m_a, size); + } + + // Inspired by https://github.com/pgii/FiltfiltSharp + private double[] FiltFilt(double[] x) + { + if (x == null) + throw new ArgumentNullException(nameof(x)); + + int nx = x.Length; + int na = m_a.Count; + int nb = m_b.Count; + int order = Math.Max(nb, na); + int factor = 3 * (order - 1); + + if (nx <= factor) + throw new ArgumentOutOfRangeException(nameof(x), "FiltFilt signal 'x' size is too small"); + + //Resize(m_b, order); + //Resize(m_a, order); + + List rows = new(order); + List cols = new(order); + + AddIndexRange(rows, 0, order - 2); + + if (order > 2) + { + AddIndexRange(rows, 1, order - 2); + AddIndexRange(rows, 0, order - 3); + } + + AddIndexCount(cols, 0, order - 1); + + if (order > 2) + { + AddIndexRange(cols, 1, order - 2); + AddIndexRange(cols, 1, order - 2); + } + + int count = rows.Count; + List data = new(count); + + Resize(data, count); + + data[0] = 1 + m_a[1]; + int j = 1; + + if (order > 2) + { + for (int i = 2; i < order; i++) + data[j++] = m_a[i]; + + for (int i = 0; i < order - 2; i++) + data[j++] = 1.0; + + for (int i = 0; i < order - 2; i++) + data[j++] = -1.0; + } + + List leftPad = SubvectorReverse(x, factor, 1); + leftPad = leftPad.Select(q => 2 * x[0] - q).ToList(); + + List rightPad = SubvectorReverse(x, nx - 2, nx - factor - 1); + rightPad = rightPad.Select(q => 2 * x[nx - 1] - q).ToList(); + + int maxRowCount = rows.Max() + 1; + int maxColCount = cols.Max() + 1; + + List signal1 = new(nx + factor * 2); + List signal2 = new(signal1.Capacity); + List zi = new(maxRowCount); + + signal1.AddRange(leftPad); + signal1.AddRange(x); + signal1.AddRange(rightPad); + + Matrix sp = new(maxRowCount, maxColCount, 0); + + for (int k = 0; k < count; ++k) + sp[rows[k]][cols[k]] = data[k]; + + double[] segment1 = Segment(m_b, 1, order - 1); + double[] segment2 = Segment(m_a, 1, order - 1); + Matrix mzi = sp.Inverse() * Calc(segment1, m_b.ToArray()[0], segment2); + + Resize(zi, mzi.NRows, 1); + ScaleZi(mzi, zi, signal1[0]); + Filter(signal1, signal2, zi); + + signal2.Reverse(); + + ScaleZi(mzi, zi, signal2[0]); + Filter(signal2, signal1, zi); + + return SubvectorReverse(signal1, signal1.Count - factor - 1, factor).ToArray(); + } + + private void Filter(List x, List y, List zi) + { + if (m_a.Count == 0) + throw new Exception("FiltFilt coefficient 'a' is empty"); + + if (m_a.All(val => val == 0.0D)) + throw new Exception("FiltFilt coefficient 'a' must have at least one non-zero number"); + + if (m_a[0] == 0) + throw new Exception("FiltFilt coefficient 'a' first element cannot be zero"); + + m_a = m_a.Select(q => q / m_a[0]).ToList(); + m_b = m_b.Select(q => q / m_a[0]).ToList(); + + int nx = x.Count; + int order = Math.Max(m_a.Count, m_b.Count); + + Resize(m_b, order); + Resize(m_a, order); + Resize(zi, order); + Resize(y, nx); + + for (int i = 0; i < nx; i++) + { + int index = order - 1; + + while (index != 0) + { + if (i >= index) + zi[index - 1] = m_b[index] * x[i - index] - m_a[index] * y[i - index] + zi[index]; + + index--; + } + + y[i] = m_b[0] * x[i] + zi[0]; + } + + zi.RemoveAt(zi.Count - 1); + } + + private static void Resize(List vector, int length, double empty = 0.0D) + { + if (vector.Count >= length) + return; + + vector.AddRange(Enumerable.Repeat(empty, length - vector.Count)); + } + + private static void AddIndexRange(List vector, int start, int stop, int increment = 1) + { + for (int i = start; i <= stop; i += increment) + vector.Add(i); + } + + private static void AddIndexCount(List vector, int value, int count) + { + while (count-- != 0) + vector.Add(value); + } + + private static List SubvectorReverse(IReadOnlyList vector, int stop, int start) + { + int length = stop - start + 1; + List result = new(Enumerable.Repeat(0.0D, length)); + int endIndex = length - 1; + + for (int i = start; i <= stop; i++) + result[endIndex--] = vector[i]; + + return result; + } + + private static void ScaleZi(Matrix mzi, List zi, double factor) + { + for (int i = 0; i < mzi.NRows; i++) + zi[i] = mzi[i][0] * factor; + } + + private static Matrix Calc(double[] segment1, double d, double[] segment2) + { + Matrix result = new (segment1.Length, 1, 0.0D); + + for (int i = 0; i < segment1.Length; i++) + result[i][0] = segment1[i] - d * segment2[i]; + + return result; + } + + private static double[] Segment(List vector, int start, int stop) + { + int length = stop - start + 1; + double[] result = new double[length]; + + for (int i = 0; i < length; i++) + result[i] = vector[start + i]; + + return result; + } + + public static double[] FiltFilt(double[] b, double[] a, double[] x) + { + if (b == null) + throw new ArgumentNullException(nameof(b)); + + if (a == null) + throw new ArgumentNullException(nameof(a)); + + if (b.Length == 0 || a.Length == 0) + throw new ArgumentException("FiltFilt coefficients b and a cannot be empty."); + + DigitalFilter filter = new(b, a, Math.Max(b.Length, a.Length)); + return filter.FiltFilt(x); + } +} diff --git a/src/Gemstone.Numeric/Analysis/NumericAnalysisExtensions.cs b/src/Gemstone.Numeric/Analysis/NumericAnalysisExtensions.cs index 1cce89683..6cacd6aa6 100644 --- a/src/Gemstone.Numeric/Analysis/NumericAnalysisExtensions.cs +++ b/src/Gemstone.Numeric/Analysis/NumericAnalysisExtensions.cs @@ -29,7 +29,10 @@ using System; using System.Collections.Generic; +using System.DirectoryServices.ActiveDirectory; using System.Linq; +using System.Numerics; +using Gemstone.Collections.CollectionExtensions; namespace Gemstone.Numeric.Analysis; @@ -166,4 +169,27 @@ public static float StandardDeviation(this IEnumerable source, Func + /// median calculates the median values of a sequence of values, while being aware of NaN values. + /// + /// Source data sample. + /// A Flag indicating whether any NaN should be skipped. + /// + /// is null + public static double NaNAwareMedian(this IEnumerable source, bool ommitNaN=true) + { + if (source is null) + throw new ArgumentNullException(nameof(source), "source is null"); + + if (!ommitNaN && source.Any(d => double.IsNaN(d))) + return double.NaN; + + IEnumerable sampleData = source.Where(d => !double.IsNaN(d)); + + if (sampleData.Count() == 0) + return double.NaN; + + return sampleData.Median()?.Average() ?? double.NaN; + } } diff --git a/src/Gemstone.Numeric/Analysis/VariableModeDecomposition.cs b/src/Gemstone.Numeric/Analysis/VariableModeDecomposition.cs new file mode 100644 index 000000000..3e8114bf3 --- /dev/null +++ b/src/Gemstone.Numeric/Analysis/VariableModeDecomposition.cs @@ -0,0 +1,419 @@ +//****************************************************************************************************** +// VariableModeDecomposition.cs - Gbtc +// +// Copyright © 2014, Grid Protection Alliance. All Rights Reserved. +// +// Licensed to the Grid Protection Alliance (GPA) under one or more contributor license agreements. See +// the NOTICE file distributed with this work for additional information regarding copyright ownership. +// The GPA licenses this file to you under the MIT License (MIT), the "License"; you may +// not use this file except in compliance with the License. You may obtain a copy of the License at: +// +// http://www.opensource.org/licenses/MIT +// +// Unless agreed to in writing, the subject software distributed under the License is distributed on an +// "AS-IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. Refer to the +// License for the specific language governing permissions and limitations. +// +// Code Modification History: +// ---------------------------------------------------------------------------------------------------- +// 04/30/2014 - Stephen E. Chisholm +// Generated original version of source code. +// +//****************************************************************************************************** + +//------------------------------------------------------------------------------------------------------ +// Code base on Wikipedia article http://en.wikipedia.org/wiki/Box–Muller_transform +//------------------------------------------------------------------------------------------------------ + +using System; +using System.Collections; +using System.Collections.Generic; +using System.ComponentModel.Design; +using System.Diagnostics; +using System.DirectoryServices.ActiveDirectory; +using System.Drawing; +using System.Linq; +using System.Reflection; +using System.Security.Cryptography; +using System.Security.Cryptography.Xml; +using System.Xml.Linq; +using Gemstone.Numeric.Random.Normal; +using Gemstone.Numeric.Random.Uniform; +using JetBrains.Annotations; +using MathNet.Numerics; +using MathNet.Numerics.Data.Matlab; +using MathNet.Numerics.LinearAlgebra; +using MathNet.Numerics.LinearAlgebra.Factorization; +using MathNet.Numerics.LinearAlgebra.Solvers; +using MathNet.Numerics.Random; +using Microsoft.Extensions.Hosting; + +namespace Gemstone.Numeric.Analysis; + +/// +/// Implements a BoxMuller method for generating statistically normal random numbers. +/// +public static class VariableModeDecomposition +{ + + //public static void Test() + //{ + // string fileName = "D:\\Matlab\\DEF Algorithm\\m-code\\Test.mat"; + // Matrix x = MatlabReader.Read(fileName, "x"); + // vmd(x.ToColumnArrays().First()); + //} + + + /// + /// returns intrinsic mode functions (IMFs) and a residual signal corresponding to the variational mode decomposition (VMD) of X, with default decomposition parameters. + /// + /// The input signal to oeprate on + /// The maximum Numebr of Itterations. + /// The number of intrinsic mode Functions to be found. + /// The penalty factor for the VMD algorithm. + public static Matrix vmd(IEnumerable x, int MaxIterations = 500, int NumIMFs = 5, double PenaltyFactor = 1000) + { + // not sure how we get the opts... + int nFFt = 2 * x.Count(); + + double tau = 0.01; + int nSignalLength = x.Count(); + int nHalfSignalLenght = (int)Math.Floor(nSignalLength * 0.5); + double relativeDiff = double.PositiveInfinity; + double relativeTolerance = 0.005; + double absoluteDiff = double.PositiveInfinity; + double absoluteTolerance = 1e-6; + int nMirroredSignalLength = nSignalLength * 2 + (nHalfSignalLenght - (int)Math.Ceiling(nSignalLength * 0.5)); + int nFFTLength = nMirroredSignalLength; + + int NumHalfFreqSamples; + if (nFFTLength % 2 == 0) + NumHalfFreqSamples = nFFTLength / 2 + 1; + else + NumHalfFreqSamples = (nFFTLength + 1) / 2; + + Matrix initIMFfd = new(NumHalfFreqSamples, NumIMFs, new ComplexNumber(0.0D, 0.0D)); + Matrix initialLM = new(NumHalfFreqSamples, 1, new(0.0D, 0)); + + Matrix xComplex = new(x.Select(v => new ComplexNumber(v, 0.0D)).ToArray(), 1); + Matrix sigFD = SignalBoundary(xComplex, nHalfSignalLenght, nMirroredSignalLength, nSignalLength, false).GetSubmatrix(0, 0, NumHalfFreqSamples, 1); + + // fft for initial IMFs and get half of bandwidth + initIMFfd.OperateByColumn((c) => + { + Complex32[] fft = c.Select((v) => new Complex32((float)v.Real, (float)v.Imaginary)).Concat(Enumerable.Repeat(new Complex32(0, 0), nFFt - c.Length)).ToArray(); + + MathNet.Numerics.IntegralTransforms.Fourier.Forward(fft, MathNet.Numerics.IntegralTransforms.FourierOptions.Matlab); + c = fft.Take(NumHalfFreqSamples).Select((c) => new ComplexNumber((double)c.Real, (double)c.Imaginary)).ToArray(); + }); + + + initIMFfd = initIMFfd + 2.2204e-16; + + //IMFfd = initIMFfd; + ComplexNumber[] sumIMF = initIMFfd.RowSums; + Matrix LM = (Matrix)initialLM.Clone(); + + // Frequency vector from[0, 0.5) for odd nfft and[0, 0.5] for even nfft + Matrix f; + if (nFFt % 2 == 0) + f = new(Enumerable.Range(0, (nFFt / 2)).Select((v) => (double)(v / (double)nFFt)).ToArray(), 1); + else + f = new(Enumerable.Range(0, (nFFt / 2) + 1 ).Select((v) => (double)(v / (double)nFFt)).ToArray(), 1); + + // Get the initial central frequencies + double[] centralFreq = InitialCentralFreqByFindPeaks(sigFD.GetColumnEnumerable(0).Select(v => v.Magnitude), f.GetColumnEnumerable(0), nFFTLength, NumIMFs); + + int iter = 0; + Matrix initIMFNorm = initIMFfd.TransformByValue(((v, i, j) => v.Magnitude * v.Magnitude)); + Matrix normIMF = new(initIMFNorm.NRows, initIMFNorm.NColumns, 0.0D); + + Matrix imffd = (Matrix)initIMFfd.Clone(); + + while (iter < MaxIterations && (relativeDiff > relativeTolerance || absoluteDiff > absoluteTolerance)) + { + for (int kk = 0; kk < NumIMFs; kk++) + { + sumIMF = sumIMF.Select((s, i) => s - imffd[i][kk]).ToArray(); + imffd.ReplaceSubmatrix(sigFD.TransformByValue((v, i, j) => (v - sumIMF[i] + LM[i][j] * 0.5D) / (PenaltyFactor * (f[i][j] - centralFreq[kk]) * (f[i][j] - centralFreq[kk]) + 1.0D)), 0, kk); + normIMF.ReplaceSubmatrix(new(imffd.GetColumnEnumerable(kk).Select((v, i) => v.Magnitude * v.Magnitude).ToArray(), 1), 0, kk); + centralFreq[kk] = f.TransposeAndMultiply(normIMF.GetColumn(kk))[0][0] / normIMF.GetColumnEnumerable(kk).Sum(); + sumIMF = sumIMF.Select((s, i) => s + imffd[i][kk]).ToArray(); + } + + + LM = LM + tau * (sigFD - sumIMF); + double[] absDiff = (imffd - initIMFfd).TransformByValue((v, i, j) => v.Magnitude * v.Magnitude / imffd.NRows).ColumnSums; + absoluteDiff = absDiff.Sum(); + relativeDiff = absDiff.Select((v, i) => v / initIMFNorm.GetColumn(i).Average()).Sum(); + + int[] sortedIndex = imffd.TransformByValue((v, i, j) => v.Magnitude * v.Magnitude).ColumnSums.Select((v, i) => new Tuple(i, v)).OrderByDescending((v) => v.Item2).Select(v => v.Item1).ToArray(); + imffd = Matrix.Combine(sortedIndex.Select((i) => new Matrix(imffd.GetColumn(i), 1)).ToArray()); + + centralFreq = sortedIndex.Take(centralFreq.Length).Select((i) => centralFreq[i]).ToArray(); + + initIMFfd = (Matrix)imffd.Clone(); + initIMFNorm = (Matrix)normIMF; + + iter = iter + 1; + } + + // Transform to time domain + Matrix IMFfdFull = new(nFFt, NumIMFs, new ComplexNumber(0.0D, 0.0D)); + IMFfdFull.ReplaceSubmatrix(imffd, 0, 0); + + + if (nFFTLength % 2 == 0) + IMFfdFull.ReplaceSubmatrix(imffd.FlipUpsideDown().GetSubmatrix(1,0,imffd.NRows - 2, imffd.NColumns).TransformByValue((c) => c.Conjugate), imffd.NRows, 0); + else + IMFfdFull.ReplaceSubmatrix(imffd.FlipUpsideDown().GetSubmatrix(0, 0, imffd.NRows - 1, imffd.NColumns), imffd.NRows, 0); + + + + IEnumerable sortIndex = centralFreq.Select((v, i) => new Tuple(i, v)).OrderByDescending((v) => v.Item2).Select(v => v.Item1); + + return Matrix.Combine(sortIndex.Select((i) => SignalBoundary(new(IMFfdFull.GetColumn(i), 1), nHalfSignalLenght, nMirroredSignalLength, nSignalLength, true).TransformByValue((v, i, j) => v.Real)).ToArray()); + } + + /// + /// signalBoundary applies mirroring to signal if ifInverse is 0 and removes mirrored signal otherwise.Mirror extension of the signal by half its + /// length on each side.Removing mirrored signal is a inverse process of the mirror extension. + /// + /// + /// + /// + /// + private static Matrix SignalBoundary(Matrix x, int nHalfSignalLenght, int nMirroredSignalLength, int nSignalLength, bool isInverse) + { + Complex32[] fft; + + if (isInverse) //removed mirrored signal + { + fft = x.GetColumn(0).Select((v) => new Complex32((float)v.Real, (float)v.Imaginary)).ToArray(); + + MathNet.Numerics.IntegralTransforms.Fourier.Inverse(fft, MathNet.Numerics.IntegralTransforms.FourierOptions.Matlab); + return new(fft.Select((c) => new ComplexNumber((double)c.Real,(double)c.Imaginary)).Skip(nHalfSignalLenght).Take(nMirroredSignalLength-2*nHalfSignalLenght).ToArray(),1); + } + + fft = x.GetColumnEnumerable(0).Take(nHalfSignalLenght).Reverse().Concat(x.GetColumnEnumerable(0)).Concat(x.GetColumnEnumerable(0).Skip((int)Math.Ceiling(nSignalLength / 2.0D)).Reverse()) + .Select((v) => new Complex32((float)v.Real, (float)v.Imaginary)).ToArray(); + + MathNet.Numerics.IntegralTransforms.Fourier.Forward(fft, MathNet.Numerics.IntegralTransforms.FourierOptions.Matlab); + return new (fft.Select((c) => new ComplexNumber((double)c.Real, (double)c.Imaginary)).ToArray(),1); + } + + /// + /// Initialize central frequencies by finding the locations of signal peaks + /// in frequency domain by using findpeaks function.The number of peaks is + /// determined by NumIMFs. + /// + /// + /// + /// + /// + private static double[] InitialCentralFreqByFindPeaks(IEnumerable x, IEnumerable f, int FFTLength, int NumIMFs) + { + double BW = 2.0D / FFTLength; + double minBWGapIndex = 2.0D * BW / f.Skip(1).First(); + double xMean = x.Average(); + IEnumerable xfilt = x.Select(v => v < xMean ? xMean : v); + bool[] TF = IsLocalMax(xfilt, (int)minBWGapIndex); + + + IEnumerable> peaks = xfilt.Zip(f, (v, freq) => new Tuple(v, freq)).Where((v,i) => TF[i]); + + int numpPeaks = peaks.Count(); + + + // Check for DC component + if (x.First() >= x.Skip(1).First()) + peaks = peaks.Prepend(new(x.First(), f.First())); + + + UniformRandomNumberGenerator random = new UniformRandomNumberGenerator(0); + IEnumerable centralFreq = random.Next(NumIMFs).Select(v => v.Value); + + if (peaks.Count() < NumIMFs) + centralFreq = peaks.Select(v => v.Item2).Concat(centralFreq).Take(NumIMFs); + else + centralFreq = peaks.OrderByDescending((v) => v.Item1).Take(NumIMFs).Select(v => v.Item2); + + return centralFreq.ToArray(); + } + + /// + /// returns a logical array whose elements are true when a local maximum is detected in the corresponding element of x. + /// + /// + /// + /// + private static bool[] IsLocalMax(IEnumerable x, int minSepperation = 0) + { + + IEnumerable s = x.Skip(1).Zip(x, (a, b) => a > b ? 1.0 : (a == b ? 0 : -1.0)); + if (!s.Any(v => v != 0)) + return Enumerable.Repeat(false, x.Count()).ToArray(); + + + s = FillMissing(s, v => v == 0).ToArray(); + + IEnumerable maxVals = s.Skip(1).Zip(s, (a, b) => a < b).Prepend(false).Append(false); + + if (minSepperation == 0) + return maxVals.ToArray(); + + + // Compute inflectionPts + IEnumerable inflectionPoints = x.Skip(1).Zip(x, (a, b) => a != b).Prepend(true); + inflectionPoints = s.Skip(1).Zip(s, (a, b) => a != b).Zip(inflectionPoints.Skip(1), (a, b) => a && b).Prepend(true).Append(true); + + // This will also restrict to the top N most prominent maxima. + int[] locMaxima = maxVals.Select((t, i) => new Tuple(i, t)).Where(v => v.Item2).Select(v => v.Item1).ToArray(); + + double[] P = ComputeProminence(x, locMaxima).ToArray(); + + Tuple[] flatIndices = locMaxima.Select((v) => { + double val = x.ElementAt(v); + var right = x.Select((val2, i) => (Value: val2, Index: i-1)).Skip(v).Where((o) => o.Value != val).FirstOrDefault(); + if (right == default) + return new Tuple(v, x.Count() - 1); + + return new Tuple(v,right.Index); + }).ToArray(); + + // Iterate through each local maxima. + int left = 0; + int right = 0; + + bool[] filteredMaxVals = maxVals.ToArray(); + for(int i = 0; i < locMaxima.Length; i++) + { + while ((flatIndices[i].Item1 - flatIndices[left].Item2) >= minSepperation) + left++; + + right = Math.Max(right, i); + + while ((right <= (locMaxima.Length - 2)) && ((flatIndices[right + 1].Item1 - flatIndices[i].Item2) < minSepperation)) + right = right + 1; + + IEnumerable leftIdx = locMaxima.Skip(left).Take(i - left); + leftIdx = leftIdx.Where((i) => filteredMaxVals[i]); + + if (leftIdx.Count() > 0) + { + double leftMax = leftIdx.Select((i) => P[i]).Max(); + if (leftMax >= P[locMaxima[i]]) + filteredMaxVals[locMaxima[i]] = false; + } + if (right - i > 0) + { + double rightMax = locMaxima.Skip(i+1).Take(right-i).Select((index) => P[index]).Max(); + if (rightMax > P[locMaxima[i]]) + filteredMaxVals[locMaxima[i]] = false; + } + } + + IEnumerable changeIndices = x.Skip(1).Zip(x, (a, b) => a - b).Select((v, i) => new Tuple(v, i)).Where((s) => s.Item1 != 0).Select((v) => v.Item2 + 1).Prepend(1); + Tuple[] flatRanges = changeIndices.Skip(1).Zip(changeIndices,(i1,i2) => new Tuple(i2,i1 -1)).Append(new Tuple(changeIndices.Last() + 1, x.Count() - 1)).ToArray(); + + foreach (Tuple range in flatRanges) + { + if (filteredMaxVals[range.Item1] && range.Item1 < range.Item2) + { + filteredMaxVals[range.Item1] = false; + int center = (int)Math.Floor((range.Item1 + range.Item2) * 0.5D); + filteredMaxVals[center] = true; + } + } + + return filteredMaxVals; + + + } + + private static IEnumerable ComputeProminence(IEnumerable data, IEnumerable localMaxIndices) + { + IEnumerator dataEnumerator = data.GetEnumerator(); + IEnumerator localMaxEnumerator = localMaxIndices.GetEnumerator(); + + localMaxEnumerator.MoveNext(); + int i = 0; + while (dataEnumerator.MoveNext()) + { + if (i != localMaxEnumerator.Current) + { + yield return 0.0D; + i++; + continue; + } + + int right = i; + int left = i; + double v = dataEnumerator.Current; + while (v <= dataEnumerator.Current && left > 0) + { + left--; + v = data.ElementAt(left); + } + v = dataEnumerator.Current; + while (v <= dataEnumerator.Current && right < (data.Count() - 1)) + { + right++; + v = data.ElementAt(right); + } + + double minLeft = data.Skip(left).Take(i - left).Min(); + double minRight = data.Skip(i + 1).Take(right - i).Min(); + + yield return dataEnumerator.Current - Math.Max(minLeft, minRight); + i++; + + if (!localMaxEnumerator.MoveNext()) + break; + } + while (dataEnumerator.MoveNext()) + yield return 0.0D; + } + + /// + /// Replace values meeting condition in with the next Value not meeting this criteria. + /// + /// + /// + /// + /// + /// + private static IEnumerable FillMissing(IEnumerable data, Func replace) + { + + + T last = data.Reverse().First((v) => !replace(v)); + + IEnumerator enumerator = data.GetEnumerator(); + IEnumerator checkEnumerator = data.GetEnumerator(); + + int i = 0, j = 0; + + if(!checkEnumerator.MoveNext()) + yield break; + + while (enumerator.MoveNext()) + { + if (!replace(enumerator.Current)) + yield return enumerator.Current; + else + { + while ((i > j || replace(checkEnumerator.Current)) && checkEnumerator.MoveNext()) + { + j++; + } + + if (replace(checkEnumerator.Current)) + yield return last; + else + yield return checkEnumerator.Current; + } + i++; + } + } +} diff --git a/src/Gemstone.Numeric/ComplexNumber.cs b/src/Gemstone.Numeric/ComplexNumber.cs index 4132f8315..92de795d8 100644 --- a/src/Gemstone.Numeric/ComplexNumber.cs +++ b/src/Gemstone.Numeric/ComplexNumber.cs @@ -72,6 +72,8 @@ PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR using Gemstone.Numeric.ComplexExtensions; using Gemstone.Units; using System.Numerics; +using System.Globalization; +using System.Diagnostics.CodeAnalysis; // ReSharper disable CompareOfFloatsByEqualityOperator // ReSharper disable PossibleInvalidOperationException @@ -80,7 +82,7 @@ namespace Gemstone.Numeric; /// /// Represents a complex number. /// -public struct ComplexNumber : IEquatable +public struct ComplexNumber : INumberBase, IComparisonOperators { #region [ Members ] @@ -442,6 +444,32 @@ public static implicit operator ComplexNumber(Complex value) => return new ComplexNumber(real, imaginary); } + public static bool operator >(ComplexNumber left, ComplexNumber right) => + left.Magnitude > right.Magnitude; + + public static bool operator >=(ComplexNumber left, ComplexNumber right) => + left.Magnitude >= right.Magnitude; + + public static bool operator <(ComplexNumber left, ComplexNumber right) => + left.Magnitude < right.Magnitude; + + public static bool operator <=(ComplexNumber left, ComplexNumber right) => + left.Magnitude <= right.Magnitude; + + public static ComplexNumber operator --(ComplexNumber value) + { + value.Magnitude--; + return value; + } + + public static ComplexNumber operator ++(ComplexNumber value) + { + value.Magnitude++; + return value; + } + + public static ComplexNumber operator +(ComplexNumber value) => value; + /// /// Returns specified raised to the specified power. /// @@ -478,5 +506,158 @@ public static ComplexNumber op_Exponent(ComplexNumber z, double y) => public static ComplexNumber Parse(string str) => str.FromComplexNotation(); + public static ComplexNumber One => new(1.0,0.0); + + public static int Radix => 10; + + public static ComplexNumber Zero => new(0.0,0.0); + + public static ComplexNumber AdditiveIdentity => Zero; + + public static ComplexNumber MultiplicativeIdentity => One; + + public static ComplexNumber Abs(ComplexNumber value) => value.Magnitude; + + + public static bool IsCanonical(ComplexNumber value) + { + throw new NotImplementedException(); + } + + public static bool IsComplexNumber(ComplexNumber value) => value.Imaginary != 0.0D && value.Real != 0.0D; + + public static bool IsEvenInteger(ComplexNumber value) => double.IsEvenInteger(value.Real) && double.IsEvenInteger(value.Imaginary); + + public static bool IsFinite(ComplexNumber value) => double.IsFinite(value.Magnitude); + + public static bool IsImaginaryNumber(ComplexNumber value) => value.Imaginary != 0.0D && value.Real == 0.0D; + + public static bool IsInfinity(ComplexNumber value) => double.IsInfinity(value.Magnitude); + + public static bool IsInteger(ComplexNumber value) => double.IsInteger(value.Magnitude); + + public static bool IsNaN(ComplexNumber value) => double.IsNaN(value.Real) || double.IsNaN(value.Imaginary) ; + + public static bool IsNegative(ComplexNumber value) => value.Real < 0.0D; + + public static bool IsNegativeInfinity(ComplexNumber value) => double.IsNegativeInfinity(value.Real) || double.IsNegativeInfinity(value.Imaginary); + + public static bool IsNormal(ComplexNumber value) + { + throw new NotImplementedException(); + } + + public static bool IsOddInteger(ComplexNumber value) => double.IsOddInteger(value.Real) && double.IsOddInteger(value.Imaginary); + + public static bool IsPositive(ComplexNumber value) => value.Real > 0.0D; + + public static bool IsPositiveInfinity(ComplexNumber value) => double.IsPositiveInfinity(value.Real) || double.IsPositiveInfinity(value.Imaginary); + + public static bool IsRealNumber(ComplexNumber value) => value.Imaginary == 0.0D && value.Real != 0.0D; + + public static bool IsSubnormal(ComplexNumber value) + { + throw new NotImplementedException(); + } + + public static bool IsZero(ComplexNumber value) => value.Magnitude == 0.0D; + + public static ComplexNumber MaxMagnitude(ComplexNumber x, ComplexNumber y) => x.Magnitude > y.Magnitude ? x : y; + + public static ComplexNumber MaxMagnitudeNumber(ComplexNumber x, ComplexNumber y) + { + throw new NotImplementedException(); + } + + public static ComplexNumber MinMagnitude(ComplexNumber x, ComplexNumber y) + { + throw new NotImplementedException(); + } + + public static ComplexNumber MinMagnitudeNumber(ComplexNumber x, ComplexNumber y) + { + throw new NotImplementedException(); + } + + public static ComplexNumber Parse(ReadOnlySpan s, NumberStyles style, IFormatProvider? provider) + { + throw new NotImplementedException(); + } + + public static ComplexNumber Parse(string s, NumberStyles style, IFormatProvider? provider) + { + throw new NotImplementedException(); + } + + public static bool TryParse(ReadOnlySpan s, NumberStyles style, IFormatProvider? provider, [MaybeNullWhen(false)] out ComplexNumber result) + { + throw new NotImplementedException(); + } + + public static bool TryParse([NotNullWhen(true)] string? s, NumberStyles style, IFormatProvider? provider, [MaybeNullWhen(false)] out ComplexNumber result) + { + throw new NotImplementedException(); + } + + public bool TryFormat(Span destination, out int charsWritten, ReadOnlySpan format, IFormatProvider? provider) + { + throw new NotImplementedException(); + } + + public string ToString(string? format, IFormatProvider? formatProvider) + { + throw new NotImplementedException(); + } + + public static ComplexNumber Parse(ReadOnlySpan s, IFormatProvider? provider) + { + throw new NotImplementedException(); + } + + public static bool TryParse(ReadOnlySpan s, IFormatProvider? provider, [MaybeNullWhen(false)] out ComplexNumber result) + { + throw new NotImplementedException(); + } + + public static ComplexNumber Parse(string s, IFormatProvider? provider) + { + throw new NotImplementedException(); + } + + public static bool TryParse([NotNullWhen(true)] string? s, IFormatProvider? provider, [MaybeNullWhen(false)] out ComplexNumber result) + { + throw new NotImplementedException(); + } + + static bool INumberBase.TryConvertFromChecked(TOther value, out ComplexNumber result) + { + throw new NotImplementedException(); + } + + static bool INumberBase.TryConvertFromSaturating(TOther value, out ComplexNumber result) + { + throw new NotImplementedException(); + } + + static bool INumberBase.TryConvertFromTruncating(TOther value, out ComplexNumber result) + { + throw new NotImplementedException(); + } + + static bool INumberBase.TryConvertToChecked(ComplexNumber value, out TOther result) + { + throw new NotImplementedException(); + } + + static bool INumberBase.TryConvertToSaturating(ComplexNumber value, out TOther result) + { + throw new NotImplementedException(); + } + + static bool INumberBase.TryConvertToTruncating(ComplexNumber value, out TOther result) + { + throw new NotImplementedException(); + } + #endregion } diff --git a/src/Gemstone.Numeric/Gemstone.Numeric.csproj b/src/Gemstone.Numeric/Gemstone.Numeric.csproj index a1054ad58..1fdb74a8c 100644 --- a/src/Gemstone.Numeric/Gemstone.Numeric.csproj +++ b/src/Gemstone.Numeric/Gemstone.Numeric.csproj @@ -59,6 +59,8 @@ + + diff --git a/src/Gemstone.Numeric/Interpolation/Pchip.cs b/src/Gemstone.Numeric/Interpolation/Pchip.cs new file mode 100644 index 000000000..7b452976c --- /dev/null +++ b/src/Gemstone.Numeric/Interpolation/Pchip.cs @@ -0,0 +1,169 @@ +//****************************************************************************************************** +// Pchip.cs - Gbtc +// +// Copyright © 2022, Grid Protection Alliance. All Rights Reserved. +// +// Licensed to the Grid Protection Alliance (GPA) under one or more contributor license agreements. See +// the NOTICE file distributed with this work for additional information regarding copyright ownership. +// The GPA licenses this file to you under the MIT License (MIT), the "License"; you may not use this +// file except in compliance with the License. You may obtain a copy of the License at: +// +// http://opensource.org/licenses/MIT +// +// Unless agreed to in writing, the subject software distributed under the License is distributed on an +// "AS-IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. Refer to the +// License for the specific language governing permissions and limitations. +// +// Code Modification History: +// ---------------------------------------------------------------------------------------------------- +// 04/04/2022 - Stephen C. Wills +// Generated original version of source code. +// +//****************************************************************************************************** + +using System; +using System.Collections; +using System.Collections.Generic; +using System.Linq; +using System.Numerics; + +namespace Gemstone.Numeric.Interpolation; + +/// +/// Represents a piecebut-three-order Hermite Interpolation similar to Matlab Pchip. +/// +public static class Pchip +{ + private static double ExteriorSlope(double d1, double d2, double h1, double h2) + { + double s = ((2.0 * h1 + h2) * d1 - h1 * d2) / (h1 + h2); + double signd1 = Math.Sign(d1); + double signs = Math.Sign(s); + + if (signs != signd1) + { + s = 0.0; + } + else + { + signs = Math.Sign(d2); + + if (signd1 != signs && Math.Abs(s) > Math.Abs(3.0 * d1)) + s = 3.0 * d1; + } + + return s; + } + + /// + /// interpolates to find Vq, the values of the underlying function y=F(x) at the query points Xq. + /// + /// The x-values provided for the estimation of F(x). + /// The y-values provided for the estimation of y = F(x) + /// The x values to be estimated + /// the estimated y-values at location x = + public static double[] Interp1(double[] x, double[] y, double[] new_x) + { + int x_len = x.Count(); + int new_x_len = new_x.Count(); + double[] new_y =new double[new_x_len]; + + int low_ip1, ix, low_i, high_i, mid_i; + double hs, hs3, w1; + double[] del = new double[x_len - 1]; + double[] slopes = new double[x_len]; + double[] h = new double[x_len - 1]; + double[] pp_coefs = new double[x_len - 1 + 3 * (x_len - 1)]; + + for (low_ip1 = 0; low_ip1 < x_len - 1; low_ip1++) + { + hs = x[low_ip1 + 1] - x[low_ip1]; + del[low_ip1] = (y[low_ip1 + 1] - y[low_ip1]) / hs; + h[low_ip1] = hs; + } + + for (low_ip1 = 0; low_ip1 < x_len - 2; low_ip1++) + { + hs = h[low_ip1] + h[low_ip1 + 1]; + hs3 = 3.0 * hs; + w1 = (h[low_ip1] + hs) / hs3; + hs = (h[low_ip1 + 1] + hs) / hs3; + hs3 = 0.0; + + if (del[low_ip1] < 0.0) + { + if (del[low_ip1 + 1] <= del[low_ip1]) + { + hs3 = del[low_ip1] / (w1 * (del[low_ip1] / del[low_ip1 + 1]) + hs); + } + else + { + if (del[low_ip1 + 1] < 0.0) + hs3 = del[low_ip1 + 1] / (w1 + hs * (del[low_ip1 + 1] / del[low_ip1])); + } + } + else + { + if (del[low_ip1] > 0.0) + { + if (del[low_ip1 + 1] >= del[low_ip1]) + { + hs3 = del[low_ip1] / (w1 * (del[low_ip1] / del[low_ip1 + 1]) + hs); + } + else + { + if (del[low_ip1 + 1] > 0.0) + hs3 = del[low_ip1 + 1] / (w1 + hs * (del[low_ip1 + 1] / del[low_ip1])); + } + } + } + + slopes[low_ip1 + 1] = hs3; + } + + slopes[0] = ExteriorSlope(del[0], del[1], h[0], h[1]); + slopes[x_len - 1] = ExteriorSlope(del[x_len - 2], del[x_len - 3], h[x_len - 2], h[x_len - 3]); + + for (low_ip1 = 0; low_ip1 < x_len - 1; low_ip1++) + { + hs = (del[low_ip1] - slopes[low_ip1]) / h[low_ip1]; + hs3 = (slopes[low_ip1 + 1] - del[low_ip1]) / h[low_ip1]; + pp_coefs[low_ip1] = (hs3 - hs) / h[low_ip1]; + pp_coefs[low_ip1 + x_len - 1] = 2.0 * hs - hs3; + pp_coefs[low_ip1 + 2 * (x_len - 1)] = slopes[low_ip1]; + pp_coefs[low_ip1 + 3 * (x_len - 1)] = y[low_ip1]; + } + + for (ix = 0; ix < new_x_len; ix++) + { + low_i = 0; + low_ip1 = 2; + high_i = x_len; + + while (high_i > low_ip1) + { + mid_i = (low_i + high_i + 1) >> 1; + + if (new_x[ix] >= x[mid_i - 1]) + { + low_i = mid_i - 1; + low_ip1 = mid_i + 1; + } + else + { + high_i = mid_i; + } + } + + hs = new_x[ix] - x[low_i]; + hs3 = pp_coefs[low_i]; + + for (low_ip1 = 0; low_ip1 < 3; low_ip1++) + hs3 = hs * hs3 + pp_coefs[low_i + (low_ip1 + 1) * (x_len - 1)]; + + new_y[ix] = hs3; + } + + return new_y; + } +} diff --git a/src/Gemstone.Numeric/Matrix/Matrix.cs b/src/Gemstone.Numeric/Matrix/Matrix.cs new file mode 100644 index 000000000..18b014707 --- /dev/null +++ b/src/Gemstone.Numeric/Matrix/Matrix.cs @@ -0,0 +1,778 @@ +//****************************************************************************************************** +// Matrix.cs - Gbtc +// +// Copyright © 2012, Grid Protection Alliance. All Rights Reserved. +// +// Licensed to the Grid Protection Alliance (GPA) under one or more contributor license agreements. See +// the NOTICE file distributed with this work for additional information regarding copyright ownership. +// The GPA licenses this file to you under the MIT License (MIT), the "License"; you may +// not use this file except in compliance with the License. You may obtain a copy of the License at: +// +// http://www.opensource.org/licenses/MIT +// +// Unless agreed to in writing, the subject software distributed under the License is distributed on an +// "AS-IS" BASIS, WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. Refer to the +// License for the specific language governing permissions and limitations. +// +// Code Modification History: +// ---------------------------------------------------------------------------------------------------- +// 04/08/2025 - C Lackner +// Initial version of source generated. +// +//****************************************************************************************************** + + +using System; +using System.ComponentModel; +using System.Runtime.CompilerServices; +using Gemstone.Numeric.ComplexExtensions; +using Gemstone.Units; +using System.Numerics; +using System.Linq; +using System.Collections.Generic; +using System.Data; + +// ReSharper disable CompareOfFloatsByEqualityOperator +// ReSharper disable PossibleInvalidOperationException +namespace Gemstone.Numeric; + +/// +/// Represents a complex number. +/// +public struct Matrix : ICloneable where T : struct, INumberBase, IComparisonOperators +{ + #region [ Members ] + + // Fields + private T[][] m_value; + + #endregion + + #region [ Constructors ] + + /// + /// Creates a of given size, prepopulated with a value. + /// + /// The number of columns of this . + /// The number of rows of this . + /// The value to populate the with. + public Matrix(int rows, int cols, T value) + : this() + { + T[] row = new T[cols]; + for (int i = 0; i < cols; i++) + { + row[i] = value; + } + m_value = new T[rows][]; + + for (int j = 0; j < rows; j++) + { + m_value[j] = (T[])row.Clone(); + } + + } + + /// + /// Creates a from the given data. + /// + /// The data of this . + public Matrix(T[][] data) + : this() + { + m_value = data; + + } + + /// + /// Creates a of repeating rows. + /// + /// The number of rows of this . + /// The values of a row of this . + public Matrix(int rows, T[] row) + : this() + { + m_value = new T[rows][]; + + for (int j = 0; j < rows; j++) + { + m_value[j] = (T[])row.Clone(); + } + } + + /// + /// Creates a of repeating columns. + /// + /// The number of columns of this . + /// The values of a collumn of this . + public Matrix(T[] column, int cols) + : this() + { + m_value = column.Select(v => Enumerable.Repeat(v, cols).ToArray()).ToArray(); + } + + /// + /// Gets the Transpose of the . + /// + public Matrix Transpose + { + get + { + Matrix transposed = new Matrix(NColumns, NRows, default(T)); + for (int j = 0; j < NColumns; j++) + { + for (int i = 0; i < NRows; i++) + { + transposed[j][i] = this[i][j]; + } + } + + return transposed; + } + } + #endregion + + #region [ Properties ] + + /// + /// Gets the value of this . + /// + public T[][] Value + { + get => m_value; + } + + /// + /// Gets the number of columns in this . + /// + public int NColumns + { + get => NRows > 0 ? m_value[0].Length : 0; + } + + /// + /// Gets the number of rows in this . + /// + public int NRows + { + get => m_value.Length; + } + + /// + /// Gets the sum of each column of the . + /// + public T[] ColumnSums + { + get + { + T[] columnSums = new T[NColumns]; + + for (int i = 0; i < NRows; i++) + { + for (int j = 0; j < NColumns; j++) + { + if (i == 0) + columnSums[j] = m_value[i][j]; + else + columnSums[j] = columnSums[j] + m_value[i][j]; + } + } + return columnSums; + } + } + + /// + /// Gets the sum of each row of the . + /// + public T[] RowSums + { + get + { + T[] rowSums = new T[NRows]; + + for (int i = 0; i < NRows; i++) + { + for (int j = 0; j < NColumns; j++) + { + if (j == 0) + rowSums[i] = m_value[i][j]; + else + rowSums[i] = rowSums[i] + m_value[i][j]; + } + } + return rowSums; + } + } + + /// + /// Accesses the element stored at (i,j) + /// + /// + /// + public T this[Tuple key] + { + get { return m_value[key.Item1][key.Item2]; } + set { m_value[key.Item1][key.Item2] = value; } + } + + /// + /// Accesses the element stored at (i,j) + /// + /// + /// + public T[] this[int key] + { + get { return m_value[key]; } + set { m_value[key] = value; } + } + + #endregion + + #region [ Methods ] + + /// + /// Returns a resulting from transposing this matrix and multiplying it with a + /// + /// + /// + /// + /// + public Matrix TransposeAndMultiply(Matrix value) where U : struct, INumberBase, IMultiplyOperators, IComparisonOperators + { + if (NRows != value.NRows) + throw new ArgumentException("Cannot multiply matrices due to dimension missmatch."); + Matrix data = new Matrix(NColumns, value.NColumns, default(T)); + for (int i = 0; i < NColumns; ++i) + { + for (int k = 0; k < NRows; ++k) + { + for (int j = 0; j < value.NColumns; ++j) + { + if (k == 0) + data[i][j] = value[k][j] * m_value[k][i]; + else + data[i][j] += value[k][j] * m_value[k][i]; + } + } + } + return data; + } + + + /// + /// Returns a resulting from transposing this matrix and multiplying it with a 1xN + /// + /// + /// + /// + /// + public Matrix TransposeAndMultiply(U[] value) where U : struct, INumberBase, IMultiplyOperators, IComparisonOperators + { + if (NColumns == value.Length) + return TransposeAndMultiply(new Matrix(1, value)); + if (NRows == value.Length) + return TransposeAndMultiply(new Matrix(value, 1)); + throw new ArgumentException("Cannot multiply matrices due to dimension missmatch."); + } + + + /// + /// Applies the given function to each row of the . + /// + /// + public void OperateByRow(Action func) + { + for (int i = 0; i < NRows; ++i) + func.Invoke(m_value[i]); + } + + /// + /// Applies the given function to each column of the . + /// + /// + public void OperateByColumn(Action func) + { + for (int i = 0; i < NColumns; ++i) + { + T[] col = GetColumn(i); + func.Invoke(col); + for (int j = 0; j < col.Length; j++) + m_value[j][i] = col[j]; + } + } + + /// + /// Applies the given function to each value of the . + /// + /// + public void OperateByValue(Func func) + { + for (int i = 0; i < NRows; ++i) + { + for (int k = 0; k < NColumns; ++k) + { + this[i][k] = func.Invoke(this[i][k], i, k); + } + } + } + + public Matrix TransformByValue(Func func) where U : struct, INumberBase, IComparisonOperators + { + Matrix data = new Matrix(NRows, NColumns, default(U)); + for (int i = 0; i < NRows; ++i) + { + for (int k = 0; k < NColumns; ++k) + { + data[i][k] = func.Invoke(this[i][k], i, k); + } + } + return data; + } + + public Matrix TransformByValue(Func func) where U : struct, INumberBase, IComparisonOperators + => TransformByValue((v, i, j) => func.Invoke(v)); + + + /// + /// returns the speciefied row of the Matrix + /// + /// + /// + public T[] GetRow(int index) => m_value[index]; + + /// + /// returns the specified column of the Matrix + /// + /// + /// + public T[] GetColumn(int index) + { + return GetColumnEnumerable(index).ToArray(); + } + + public IEnumerable GetColumnEnumerable(int index) + { + if (index < 0 || index >= NColumns) + throw new ArgumentOutOfRangeException(nameof(index), "Index out of range."); + + int i = 0; + while (i < NRows) + { + yield return m_value[i][index]; + i++; + } + } + + /// + /// Gets a submatrix of the starting at the given row and column, with the given number of rows and columns. + /// + /// the first row of the submatrix. + /// the first column of the submarix. + /// the number of rows of the submatrix. + /// the number of columns of the submatrix. + /// A new that is the submatrix speciffied. + /// + public Matrix GetSubmatrix(int startRow, int startColumn, int numRows, int numCols) + { + if (startRow < 0 || startRow >= NRows) + throw new ArgumentOutOfRangeException(nameof(startRow), "Index out of range."); + if (startColumn < 0 || startColumn >= NColumns) + throw new ArgumentOutOfRangeException(nameof(startColumn), "Index out of range."); + + if (startColumn + numCols > NColumns) + throw new ArgumentOutOfRangeException(nameof(startColumn), "Index out of range."); + if (startRow + numRows > NRows) + throw new ArgumentOutOfRangeException(nameof(startRow), "Index out of range."); + + Matrix matrix = new Matrix(numRows, numCols, default(T)); + for (int i = 0; i < numRows; ++i) + { + for (int j = 0; j < numCols; ++j) + { + matrix[i][j] = m_value[startRow + i][startColumn + j]; + } + } + + return matrix; + } + + public void ReplaceSubmatrix(Matrix subMatrix, int startRow, int startColumn) + { + if (startRow < 0 || startRow >= NRows) + throw new ArgumentOutOfRangeException(nameof(startRow), "Index out of range."); + if (startColumn < 0 || startColumn >= NColumns) + throw new ArgumentOutOfRangeException(nameof(startColumn), "Index out of range."); + + if (startColumn + subMatrix.NColumns > NColumns) + throw new ArgumentOutOfRangeException(nameof(startColumn), "Index out of range."); + if (startRow + subMatrix.NRows > NRows) + throw new ArgumentOutOfRangeException(nameof(startRow), "Index out of range."); + + for (int i = 0; i < subMatrix.NRows; ++i) + { + for (int j = 0; j < subMatrix.NColumns; ++j) + { + this[i + startRow][j + startColumn] = subMatrix[i][j]; + } + } + } + + public Matrix PointWhiseMultiply(Matrix matrix) + { + if (this.NColumns != matrix.NColumns || matrix.NRows != this.NRows) + throw new ArgumentException("Matrix sizes do not match."); + + Matrix data = new Matrix(NRows, NColumns, default(T)); + for (int i = 0; i < NRows; ++i) + { + for (int k = 0; k < NColumns; ++k) + { + data[i][k] = this[i][k] * matrix[i][k]; + } + } + return data; + } + + public object Clone() + { + return new Matrix(m_value.Select(v => (T[])v.Clone()).ToArray()); + } + + public Matrix FlipUpsideDown() + { + Matrix data = new Matrix(NRows, NColumns, default(T)); + for (int i = 0; i < NRows; ++i) + { + data[NRows - i - 1] = (T[])m_value[i].Clone(); + } + return data; + } + + + /// + /// Returns the inverse of the matrix. + /// + public Matrix Inverse() + { + Matrix result = (Matrix)this.Clone(); + + int[] Permutation; + Matrix Lower; + Matrix Upper; + + this.LUDecomposition(out Lower, out Upper, out Permutation); + + T[] b = Enumerable.Repeat(T.Zero, NRows).ToArray(); + + for (int i = 0; i < NRows; ++i) + { + int j = Permutation.TakeWhile(x => x != i).Count() - 1; + b[j] = T.One; + + T[] x = HelperSolve(Lower, Upper, b); // + b[j] = T.Zero; + for (int k = 0; k < NRows; k++) + result[k][i] = x[k]; + } + return result; + } + + /// + /// Does an LUP Decomposition. + /// + /// The Lower Triangular Matrix + /// The upper Triangular Matrix + /// The Permutation Matrix P + /// + /// + public void LUDecomposition(out Matrix Lower, out Matrix Upper, out int[] Permutation) + { + // Doolittle LUP decomposition with partial pivoting. + // rerturns: result is L (with 1s on diagonal) and U; + // perm holds row permutations; toggle is +1 or -1 (even or odd) + if (NRows != NColumns) + throw new Exception("Attempt to decompose a non-square m"); + + Lower = new Matrix(NRows, NColumns, default(T)); + Upper = new Matrix(this); + + Permutation = Enumerable.Range(0, NRows).ToArray(); + + + + for (int j = 0; j < NColumns; ++j) // each column + { + T colMax = this[j][j]; + int pRow = j; // row of largest value in column j + for (int i = j + 1; i < NRows; ++i) + { + if (T.Abs(this[i][j]) > T.Abs(colMax)) + { + colMax = T.Abs(this[i][j]); + pRow = i; + } + } + + if (pRow != j) // if largest value not on pivot, swap rows + { + T[] row = Upper[pRow]; + Upper[j] = Upper[j].Select((v, i) => (i <= j ? v : row[i])).ToArray(); + + row = Lower[pRow]; + Lower[j] = Lower[j].Select((v, i) => (i > j ? v : row[i])).ToArray(); + + int tmp = Permutation[pRow]; // and swap perm info + Permutation[pRow] = Permutation[j]; + Permutation[j] = tmp; + } + + // if (result[j][j] == 0.0) + // { + // // find a good row to swap + // int goodRow = -1; + // for (int row = j + 1; row less - than n; + // ++row) + //{ + // if (result[row][j] != 0.0) + // goodRow = row; + // } + + // if (goodRow == -1) + // throw new Exception("Cannot use Doolittle's method"); + + // // swap rows so 0.0 no longer on diagonal + // double[] rowPtr = result[goodRow]; + // result[goodRow] = result[j]; + // result[j] = rowPtr; + + // int tmp = perm[goodRow]; // and swap perm info + // perm[goodRow] = perm[j]; + // perm[j] = tmp; + + // toggle = -toggle; // adjust the row-swap toggle + // } + // // -------------------------------------------------- + // // if diagonal after swap is zero . . + // //if (Math.Abs(result[j][j]) less-than 1.0E-20) + // // return null; // consider a throw + + for (int i = j + 1; i < NRows; ++i) + { + T lowerVal = Upper[i][j] / Upper[j][j]; + Lower[i][j] = lowerVal; + T[] row = Upper[j]; + Upper[i] = Upper[i].Select((v, k) => (k < j ? v : v - lowerVal * row[k])).ToArray(); + } + } + } + + private T[] HelperSolve(Matrix Lower, Matrix Upper, T[] b) + { + // before calling this helper, permute b using the perm array + // from MatrixDecompose that generated luMatrix + + T[] x = new T[b.Length]; + x[0] = b[0]; + + for (int i = 1; i < Lower.NRows; ++i) + { + x[i] = b[i]; + for (int j = 0; j < i; ++j) + x[i] -= Lower[i][j] * x[j]; + } + + x[Lower.NRows - 1] /= Upper[Lower.NRows - 1][Lower.NRows - 1]; + + for (int i = Lower.NRows - 2; i >= 0; --i) + { + for (int j = i + 1; j < Lower.NRows; ++j) + x[i] -= Lower[i][j] * x[j]; + x[i] = x[i] / Upper[i][i]; + } + + return x; + } + + #endregion + + #region [ Operators ] + + /// + /// Implicitly converts a to a . + /// + /// Operand. + /// ComplexNumber representing the result of the operation. + public static implicit operator Matrix(T[][] value) => + new(value); + + /// + /// Implicitly converts a to a .NET value. + /// + /// Operand. + /// Complex representing the result of the operation. + public static implicit operator T[][](Matrix value) => value.Value; + + + /// + /// Returns the negated value. + /// + /// Left hand operand. + /// ComplexNumber representing the result of the unary negation operation. + public static Matrix operator -(Matrix z) => + new(z.Value.Select(v => v.Select((x) => -x).ToArray()).ToArray()); + + /// + /// Returns computed sum of values. + /// + /// Left hand operand. + /// Right hand operand. + /// representing the result of the addition operation. + public static Matrix operator +(Matrix value1, Matrix value2) + { + if (value1.NColumns != value2.NColumns || value1.NRows != value2.NRows) + throw new ArgumentException("Matrix sizes do not match."); + return new(value1.Value.Select((v, i) => v.Select((x, j) => x + value2.Value[i][j]).ToArray()).ToArray()); + } + + /// + /// Returns computed sum of values. + /// + /// Left hand operand. + /// Right hand operand. + /// representing the result of the addition operation. + public static Matrix operator +(Matrix value1, T value2) + { + return new(value1.Value.Select((v, i) => v.Select((x, j) => x + value2).ToArray()).ToArray()); + } + + public static Matrix operator +(T value1, Matrix value2) => value2 + value1; + + + /// + /// Returns computed difference of values. + /// + /// Left hand operand. + /// Right hand operand. + /// representing the result of the subtraction operation. + public static Matrix operator -(Matrix value1, Matrix value2) + { + if (value1.NColumns != value2.NColumns || value1.NRows != value2.NRows) + throw new ArgumentException("Matrix sizes do not match."); + return new(value1.Value.Select((v, i) => v.Select((x, j) => x - value2.Value[i][j]).ToArray()).ToArray()); + } + + public static Matrix operator -(Matrix value1, T value2) + { + return new(value1.Value.Select((v, i) => v.Select((x, j) => x - value2).ToArray()).ToArray()); + } + + public static Matrix operator -(T value1, Matrix value2) => (-value2) + value1; + + public static Matrix operator -(Matrix value1, T[] value2) + { + if (value1.NColumns == value2.Length) + return new(value1.Value.Select((v, i) => v.Select((x, j) => x - value2[j]).ToArray()).ToArray()); + if (value1.NRows == value2.Length) + return new(value1.Value.Select((v, i) => v.Select((x, j) => x - value2[i]).ToArray()).ToArray()); + + throw new ArgumentException("Matrix sizes do not match."); + + + } + + /// + /// Returns computed product of values. + /// + /// Left hand operand. + /// Right hand operand. + /// representing the result of the multiplication operation. + public static Matrix operator *(Matrix value1, Matrix value2) + { + if (value1.NColumns != value2.NRows) + throw new ArgumentException("Cannot multiply matrices due to dimension missmatch."); + + Matrix data = new Matrix(value2.NRows, value2.NColumns, default(T)); + for (int i = 0; i < value1.NRows; ++i) + { + for (int k = 0; k < value2.NColumns; ++k) + { + data[i][k] = value1.GetRow(i).Zip(value2.GetColumnEnumerable(k), (v1, v2) => v1 * v2).Aggregate(default(T), (s, v) => s + v); + } + } + return data; + } + + /// + /// Returns computed product of values. + /// + /// Left hand operand. + /// Right hand operand. + /// representing the result of the multiplication operation. + public static Matrix operator *(Matrix value1, T value2) + { + return new(value1.Value.Select((v, i) => v.Select((x, j) => x * value2).ToArray()).ToArray()); + } + + /// + /// Returns computed product of values. + /// + /// Left hand operand. + /// Right hand operand. + /// representing the result of the multiplication operation. + public static Matrix operator *(T value1, Matrix value2) => value2 * value1; + + public static Matrix operator /(Matrix value1, T value2) + { + return new(value1.Value.Select((v, i) => v.Select((x, j) => x / value2).ToArray()).ToArray()); + } + + public static Matrix operator /(Matrix value1, Matrix value2) + { + if (value1.NColumns != value2.NColumns || value1.NRows != value2.NRows) + throw new ArgumentException("Matrix sizes do not match."); + + Matrix data = new Matrix(value1.NRows, value1.NColumns, default(T)); + for (int i = 0; i < value1.NRows; ++i) + { + for (int k = 0; k < value1.NColumns; ++k) + { + data[i][k] = value1[i][k] * value2[i][k]; + } + } + return data; + } + + #endregion + + #region [ Static ] + + public static Matrix Combine(Matrix[] matrices) + { + if (matrices.Length == 0) + throw new ArgumentException("List of matrices cannot be empty."); + int nCols = matrices.Sum(m => m.NColumns); + int nRows = matrices[0].NRows; + + if (matrices.Any(m => m.NRows != nRows)) + throw new ArgumentException("All matrices must have the same number of rows."); + + Matrix matrix = new Matrix(nRows, nCols, default(T)); + + int count = 0; + for (int j = 0; j < nRows; j++) + { + count = 0; + for (int i = 0; i < matrices.Length; i++) + { + for (int k = 0; k < matrices[i].NColumns; k++) + { + matrix[j][count] = matrices[i][j][k]; + count++; + } + } + } + + return matrix; + } + #endregion +}