My generic method fails because bool isOrderedField = typeof(T).IsAssignableTo(typeof(INumber<>)); evaluates (according to my intent, erroneously) to false for a generic extension method of signature public static T[][] MovingStatistics<T>(this IReadOnlyList<T> source, int radius, T[][] target, bool calculateStandardDeviation = false, bool calculateExtrema = false, int startIndex = 0, int length = -1) where T : INumberBase<T> when the latter is called on a double array. I have tried to investigate and come up with the following minimum example (I was suspicious of the fact that the culprit is an extension method, that is why I tested it both ways):
using System.Numerics;
namespace InferredGenericMethodTypeBehaviour
{
internal class Program
{
static void Main(string[] args)
{
double[] array = new double[] { 11.0 };
Console.WriteLine($"The method [{nameof(Program.IsItAListOfNumbers)}] says the statement that [{nameof(array)}] comprises numbers is {IsItAListOfNumbers(array)}.");
Console.WriteLine($"The method [{nameof(Extensions.AmIAListOfNumbers)}] says the statement that [{nameof(array)}] comprises numbers is {array.AmIAListOfNumbers()}.");
}
static bool IsItAListOfNumbers<T>(IReadOnlyList<T> list) where T : INumberBase<T> => typeof(T).IsAssignableTo(typeof(INumber<>));
}
public static class Extensions
{
public static bool AmIAListOfNumbers<T>(this IReadOnlyList<T> array) where T : INumberBase<T> => typeof(T).IsAssignableTo(typeof(INumber<>));
}
}
This is NOT an MWE as it cannot reproduce the behaviour that got me investigating in the first place. The output will be
The method [IsItAListOfNumbers] says the statement that [array] comprises numbers is True.
The method [AmIAListOfNumbers] says the statement that [array] comprises numbers is True.
In my desperation, I inserted the following line into the offending code:
bool amIEvil = new double[] { 11.0 }.ConsistsOfNumbers();
where the generic extension method ConsistsOfNumbers() is identical to the minimum example:
public static bool ConsistsOfNumbers<T>(this IReadOnlyList<T> source) where T : INumberBase<T> => typeof(T).IsAssignableTo(typeof(INumber<>));
There, in the program I am actually developing, amIEvil will evaluate to false. Both my actual program and the minimum (non-)working example run in .NET 7.0, on the same machine, under the same settings.
What can cause this difference? I can of course circumvent it in a number of ways, I have found several answers here on stack overflow, such as this one, and after all since the whole thing is to check whether the type of the number belongs to an ordered field and not a non-ordered one like complex numbers, I can simply replace it with typeof(T) != typeof(Complex), but I cannot wrap my head round this glaring inconsistency.
Your help is much appreciated.
EDIT:
In the comments, I was instructed to publish the full code and the project files, so I shall do so:
Calling code:
private async Task TestCrossSpectrumOnAcqRecording()
{
OpenFileDialog dialogue = new OpenFileDialog();
SettingsWindow settingsWindow;
double[][] data;
double[][] spectrumMovingStatistics;
double fullDuration;
double samplingInterval;
double overlap = 0.8;
double windowWidth = 1.1;
double startFrequency = 0.1;
double endFrequency = 10.0;
double frequencyStep;
int startFrequencyIndex = 0;
int frequencyCount = -1;
int xIndex = 0;
int yIndex = 1;
int radius = 33;
IReadOnlyList<AcqChannelProperties> channels;
string[] channelNames;
TetheredSun.Spectral.Window window;
WindowPropagator windowPropagator;
WindowedCrossSpectral crossSpectral;
Recording recording;
FixedStepDataSequence t;
FixedStepDataSequence averagedF;
FixedStepDataSequence totalF;
ComplexDataSequence xSpectrum;
ComplexDataSequence ySpectrum;
RecordingViewModel recordingViewModel;
dialogue.Title = "Select AcqKnowledge file to open";
dialogue.Filter = AcqFile.Filter;
if (dialogue.ShowDialog() is not true) return;
using (Stream source = new FileStream(dialogue.FileName, FileMode.Open)) {
using (AcqFile acqFile = new AcqFile(source)) {
data = await Task.Run(() => acqFile.ReadAllChannels()).ConfigureAwait(true);
samplingInterval = acqFile.SamplingInterval / 1000.0;
fullDuration = (acqFile.Count - 1) * samplingInterval;
channels = acqFile.ChannelProperties;
channelNames = channels.Select(channel => channel.Name).ToArray();
}
}
settingsWindow = new SettingsWindow();
settingsWindow.Title = "Set parameters for cross-spectral investigations";
settingsWindow.Settings.Add("Window width [s]", windowWidth);
settingsWindow.Settings.Add("Window overlap [%]", overlap * 100.0);
settingsWindow.Settings.Add("Window type", new ListMenu<Type>(windowTypes, windowType => windowType.Name.Decamelise()));
settingsWindow.Settings.Add("First term of cross spectrum (x)", new ListMenu<string>(channelNames, channelName => channelName, xIndex));
settingsWindow.Settings.Add("Second term of cross spectrum (y)", new ListMenu<string>(channelNames, channelName => channelName, yIndex));
if (settingsWindow.ShowDialog() is not true) return;
windowWidth = settingsWindow.Settings[0].GetValue<double>();
overlap = settingsWindow.Settings[1].GetValue<double>() / 100.0;
window = settingsWindow.Settings[2].GetValue<ListMenu<Type>>().Selected.New() as TetheredSun.Spectral.Window;
xIndex = Array.IndexOf(channelNames, settingsWindow.Settings[3].GetValue<ListMenu<string>>().Selected);
yIndex = Array.IndexOf(channelNames, settingsWindow.Settings[4].GetValue<ListMenu<string>>().Selected);
frequencyStep = 1.0 / windowWidth;
startFrequencyIndex = (int)(startFrequency / frequencyStep);
frequencyCount = (int)(endFrequency / frequencyStep) - startFrequencyIndex + 1;
windowPropagator = new WindowPropagator() {
Width = windowWidth,
Step = (1.0 - overlap) * windowWidth,
Window = window
};
windowPropagator.SamplingInterval = samplingInterval;
crossSpectral = new WindowedCrossSpectral() {
StartFrequencyIndex = startFrequencyIndex,
FrequencyCount = frequencyCount
};
using (WindowedFourierAnalyser windowedFourierAnalyser = new WindowedFourierAnalyser(windowPropagator)) {
windowedFourierAnalyser.Targets.Add(crossSpectral);
await Task.Run(() => windowedFourierAnalyser.Analyse(data[xIndex].Standardise().ToArray(), data[yIndex].Standardise().ToArray())).ConfigureAwait(true);
}
t = new FixedStepDataSequence() { Count = data[xIndex].Length, Step = samplingInterval, Name = "time", Unit = "s" };
averagedF = new FixedStepDataSequence() { Count = crossSpectral.CrossSpectrum.Count, Minimum = startFrequencyIndex * frequencyStep, Step = frequencyStep, Name = "frequency", Unit = "Hz" };
frequencyStep = 1.0 / (data[xIndex].Length * samplingInterval);
startFrequencyIndex = (int)(startFrequency / frequencyStep);
frequencyCount = (int)(endFrequency / frequencyStep) - startFrequencyIndex + 1;
using (FourierTransformer fourierTransformer = new FourierTransformer(data[xIndex].Length, (TetheredSun.Spectral.Window)window.GetType().New())) {
xSpectrum = new ComplexDataSequence(Spectrum.Amplitude(fourierTransformer.Transform(data[xIndex]), false, startFrequencyIndex, frequencyCount));
xSpectrum.Name = $"Spectrum of {channels[xIndex].Name}";
xSpectrum.Unit = $"{channels[xIndex].Unit} rms";
ySpectrum = new ComplexDataSequence(Spectrum.Amplitude(fourierTransformer.Transform(data[yIndex]), false, startFrequencyIndex, frequencyCount));
ySpectrum.Name = $"Spectrum of {channels[yIndex].Name}";
ySpectrum.Unit = $"{channels[yIndex].Unit} rms";
}
totalF = new FixedStepDataSequence() { Count = frequencyCount, Minimum = startFrequencyIndex * frequencyStep, Step = frequencyStep, Name = "frequency", Unit = "Hz" };
// xSpectrum is a custom Complex sequence that iterates over an alternating real-imaginary double array, it has a custom iterator that enumerates the magnitude as a double sequence, but it is converted to double[] with Enumerable.ToArray<T>(), so this should not be a problem:
spectrumMovingStatistics = xSpectrum.Magnitude.ToArray().MovingStatistics(radius, null, true, true); // This is the offending line. Will throw an exception as there will be too few items in the array.
/* I inserted the following lines for debugging:
double[] debug = new double[] { 11.0 };
bool amIEvil = debug.ConsistsOfNumbers(); xSpectrum.Magnitude.ToArray().ConsistsOfNumbers();*/
for (int i = 0; i < spectrumMovingStatistics[0].Length; i++) {
spectrumMovingStatistics[1][i] += spectrumMovingStatistics[0][i];
}
recording = new Recording() { Name = $"Cross spectrum: {channelNames[xIndex]} v {channelNames[yIndex]}" };
recording.Add(new Signal<double, double>(t, new DataSequence<double>(data[xIndex], channels[xIndex].Name, channels[xIndex].Unit)) { Name = channels[xIndex].Name });
recording.Add(new Signal<double, double>(t, new DataSequence<double>(data[yIndex], channels[yIndex].Name, channels[yIndex].Unit)) { Name = channels[yIndex].Name });
recording.Add(new Signal<double, Complex>(totalF, xSpectrum) { Name = "Spectrum 1" });
recording.Add(new Signal<double, Complex>(totalF, ySpectrum) { Name = "Spectrum 2" });
recording.Add(new Signal<double, Complex>(averagedF, new DataSequence<Complex>(crossSpectral.CrossSpectrum, $"Cross-spectrum of {channels[xIndex].Name} and {channels[yIndex].Name}", $"{channels[xIndex].Unit} ∙ {channels[yIndex].Unit} rms")) { Name = "Cross-spectrum" });
recording.Add(new Signal<double, double>(averagedF, new DataSequence<double>(crossSpectral.CrossSpectrum.Phase.Select(phase => phase / Math.PI).MovingStatistics(10).Select(stats => stats[0]).ToArray(), $"Phase lag between {channels[xIndex].Name} and {channels[yIndex].Name}", $"∙ π rad")) { Name = "Phase lag" });
recording.Add(new Signal<double, double>(averagedF, new DataSequence<double>(crossSpectral.Coherence, $"Coherence between {channels[xIndex].Name} and {channels[yIndex].Name}", "relative units")) { Name = "Coherence" });
ViewObject = (recordingViewModel = new RecordingViewModel() { Recording = recording });
recordingViewModel[2].Series.Add(new SeriesViewModel(totalF, spectrumMovingStatistics[0]) { Name = "Moving average" });
recordingViewModel[2].Series.Add(new SeriesViewModel(totalF, spectrumMovingStatistics[1]) { Name = "Upper Bollinger band" });
recordingViewModel[2].Series.Add(new SeriesViewModel(totalF, spectrumMovingStatistics[2]) { Name = "Local minimum" });
recordingViewModel[2].Series.Add(new SeriesViewModel(totalF, spectrumMovingStatistics[3]) { Name = "Local maximum" });
await Task.Run(() => recording.Save(@$"E:\Temp\{channels[xIndex].Name}-{channels[yIndex].Name} - crossSpectrum.mba")).ConfigureAwait(false);
}
The offending extension method is
public static T[][] MovingStatistics<T>(this IReadOnlyList<T> source, int radius, T[][] target, bool calculateStandardDeviation = false, bool calculateExtrema = false, int startIndex = 0, int length = -1) where T : INumberBase<T>
{
if (startIndex < 0 || startIndex >= source.Count) startIndex = 0;
int lower = startIndex - radius;
int upper = startIndex + radius;
int count;
int windowSize = 2 * radius + 1;
int targetIndex = 0;
int statisticsCount = 1;
int minimumIndex = -1;
int maximumIndex = -1;
// This is where the trouble is. This evaluates to false and thus the target array will only have 2 elements instead of 4.
bool isOrderedField = typeof(T).IsAssignableTo(typeof(INumber<>));
IComparer<T> comparer = Comparer<T>.Default;
T incoming;
T outgoing;
T nextSize;
T increment;
T incrementPerSize;
T oldMean;
T mean = T.Zero;
T q = T.Zero;
T size = T.Zero;
T minimum = T.Zero;
T maximum = T.Zero;
ReadOnlySpan<T> sourceSpan = source.ToReadOnlySpan();
if (length <= 0) length = source.Count;
if (lower < 0) lower = 0;
if (upper >= source.Count) upper = source.Count - 1;
// Here is where what went wrong above propagates its wrongness:
calculateExtrema &= isOrderedField;
if (calculateStandardDeviation) ++statisticsCount;
if (calculateExtrema) {
minimumIndex = statisticsCount;
maximumIndex = statisticsCount + 1;
minimum = T.CreateSaturating(Double.PositiveInfinity);
maximum = T.CreateSaturating(Double.NegativeInfinity);
statisticsCount += 2;
}
count = 1;
InitialiseTarget(length);
// The same iteration as below but with Span<T> objects and Unsafe.Add(ref reference, int offset) instead of array indexing -- actually worth it for 3× speed.
if (!sourceSpan.IsEmpty) return MovingStatisticsSpanIterator<T>(sourceSpan, radius, target, calculateStandardDeviation, calculateExtrema, startIndex, length);
// Initial stage: filling up the first half of the window.
for (int i = lower; i <= upper; i++) {
incoming = source[i];
increment = incoming - mean; // (x_k - A_(k-1))
size = T.CreateChecked(count++); // k in https://en.wikipedia.org/wiki/Standard_deviation#Rapid_calculation_methods
incrementPerSize = increment / size; // (x_k - A_(k-1)) / k
mean += incrementPerSize; // A_k = A_(k-1) + (x_k - A_(k-1)) /
if (calculateStandardDeviation) {
q += increment * (increment - incrementPerSize); // Q_k = Q_(k-1) + (x_k - A_(k-1)) * (x_k - A_k), => dQ_k := Q_k - Q_(k-1) = (x_k - A_(k-1)) * (x_k - A_k) =
// = increment * (x_k - [A_(k-1) + (x_k - A_(k-1)) / k]) = increment * (x_k - A_(k-1) - increment / k) = increment * (increment - incrementPerSize)
}
if (calculateExtrema) {
if (comparer.Compare(incoming, minimum) < 0) {
minimum = incoming;
}
if (comparer.Compare(incoming, maximum) > 0) {
maximum = incoming;
}
}
}
target[0][targetIndex] = mean;
if (calculateStandardDeviation) {
target[1][targetIndex] = count > 1 ? (q / (size - T.One)).SquareRoot() : T.Zero;
}
if (calculateExtrema) {
target[minimumIndex][targetIndex] = minimum;
target[maximumIndex][targetIndex] = maximum;
}
++targetIndex;
if (targetIndex >= length) return target;
// Next stage: filling up the second half of the window:
upper++;
while (upper < source.Count && count <= windowSize) {
if (targetIndex >= length) return target;
incoming = source[upper++];
increment = incoming - mean; // (x_k - A_(k-1))
size = T.CreateChecked(count++); // k in https://en.wikipedia.org/wiki/Standard_deviation#Rapid_calculation_methods
incrementPerSize = increment / size; // (x_k - A_(k-1)) / k
mean += incrementPerSize; // A_k = A_(k-1) + (x_k - A_(k-1)) / k
target[0][targetIndex] = mean;
if (calculateStandardDeviation) {
q += increment * (increment - incrementPerSize); // Q_k = Q_(k-1) + (x_k - A_(k-1)) * (x_k - A_k), => dQ_k := Q_k - Q_(k-1) = (x_k - A_(k-1)) * (x_k - A_k) =
// = increment * (x_k - [A_(k-1) + (x_k - A_(k-1)) / k]) = increment * (x_k - A_(k-1) - increment / k) = increment * (increment - incrementPerSize)
target[1][targetIndex] = (q / (size - T.One)).SquareRoot();
}
if (calculateExtrema) {
if (comparer.Compare(incoming, minimum) < 0) {
minimum = incoming;
}
if (comparer.Compare(incoming, maximum) > 0) {
maximum = incoming;
}
target[minimumIndex][targetIndex] = minimum;
target[maximumIndex][targetIndex] = maximum;
}
++targetIndex;
}
// Next stage: advance the window and correct with outgoing and incoming elements.
if (count > windowSize) {
count = windowSize;
}
size = T.CreateChecked(count);
while (upper < source.Count) {
if (targetIndex >= length) return target;
incoming = source[upper++];
outgoing = source[lower++];
increment = incoming - outgoing;
incrementPerSize = increment / size;
oldMean = mean;
mean += incrementPerSize; // A_k = A_(k-1) + (x_k - A_(k-1)) / k
target[0][targetIndex] = mean;
if (calculateStandardDeviation) {
q += increment * (incoming + outgoing - mean - oldMean);
target[1][targetIndex] = (q / (size - T.One)).SquareRoot();
}
if (calculateExtrema) {
if (outgoing == minimum || outgoing == maximum) {
(minimum, maximum) = source.Extrema(lower, count);
}
if (comparer.Compare(incoming, minimum) < 0) {
minimum = incoming;
}
if (comparer.Compare(incoming, maximum) > 0) {
maximum = incoming;
}
(target[minimumIndex][targetIndex], target[maximumIndex][targetIndex]) = (minimum, maximum);
}
++targetIndex;
}
// Last stage, wherein there are only outgoing elements:
while (count > radius) {
if (targetIndex >= length) return target;
size = T.CreateChecked(count--);
nextSize = size - T.One;
outgoing = source[lower++];
oldMean = mean;
mean = mean * (size / nextSize) - outgoing / nextSize;
target[0][targetIndex] = mean;
if (calculateStandardDeviation) {
q += size * oldMean * oldMean - outgoing * outgoing - nextSize * mean * mean;
target[1][targetIndex] = (q / (nextSize - T.One)).SquareRoot();
}
if (calculateExtrema) {
if (outgoing == minimum || outgoing == maximum) {
(minimum, maximum) = source.Extrema(lower, count);
}
(target[minimumIndex][targetIndex], target[maximumIndex][targetIndex]) = (minimum, maximum);
}
++targetIndex;
}
return target;
// Local functions:
void InitialiseTarget(int targetSize)
{
if (target is null) {
target = new T[statisticsCount][];
for (int i = 0; i < statisticsCount; i++) {
target[i] = new T[targetSize];
}
} else if (target.Length < statisticsCount) {
throw new ArgumentException($"An array of row size {target.Length} cannot accommodate moving statistics of {statisticsCount} elements.", nameof(target));
} else if (target[0].Length < targetSize) {
throw new ArgumentException($"An array of column size {target[0].Length} cannot accommodate the moving statistics of a sequence with a count of {targetSize}.", nameof(target));
}
}
}
The project file of the calling project:
<Project Sdk="Microsoft.NET.Sdk">
<PropertyGroup>
<OutputType>WinExe</OutputType>
<TargetFramework>net7.0-windows</TargetFramework>
<UseWPF>true</UseWPF>
<Platforms>AnyCPU;x64</Platforms>
<AllowUnsafeBlocks>True</AllowUnsafeBlocks>
</PropertyGroup>
<ItemGroup>
<PackageReference Include="MahApps.Metro" Version="2.4.9" />
<PackageReference Include="MahApps.Metro.IconPacks" Version="4.11.0" />
</ItemGroup>
<ItemGroup>
<ProjectReference Include="..\TetheredSun.2.0\TetheredSun.Acq\TetheredSun.Acq.csproj" />
<ProjectReference Include="..\TetheredSun.2.0\TetheredSun.Charting\TetheredSun.Charting.csproj" />
<ProjectReference Include="..\TetheredSun.2.0\TetheredSun.MahAppsExtensions\TetheredSun.MahAppsExtensions.csproj" />
<ProjectReference Include="..\TetheredSun.2.0\TetheredSun.Office\TetheredSun.Office.csproj" />
<ProjectReference Include="..\TetheredSun.2.0\TetheredSun.Spectral\TetheredSun.Spectral.csproj" />
<ProjectReference Include="..\TetheredSun.2.0\TetheredSun.Wpf\TetheredSun.Wpf.csproj" />
<ProjectReference Include="..\TetheredSun.2.0\TetheredSun\TetheredSun.csproj" />
</ItemGroup>
</Project>
and that of the project containing the called extension method
<Project Sdk="Microsoft.NET.Sdk">
<PropertyGroup>
<TargetFramework>net7.0</TargetFramework>
<Platforms>AnyCPU;x64;x86</Platforms>
</PropertyGroup>
</Project>
I could not find any global.json file in any of the projects concerned.
