Remove MathNet.Numerics & NumSharp dependencies

OnnxStack.StableDiffusion/Helpers/ArrayHelpers.cs

@@ -14,5 +14,19 @@ public static float[] Linspace(float start, float end, int partitions, bool roun
                 : result.Select(x => MathF.Round(x)).ToArray();
         }
 
+
+        public static float[] Range(int start, int end)
+        {
+            return Enumerable.Range(start, end)
+                .Select(x => (float)x)
+                .ToArray();
+        }
+
+        public static float[] Log(float[] array)
+        {
+            return array
+                .Select(x => MathF.Log(x))
+                .ToArray();
+        }
     }
 }

OnnxStack.StableDiffusion/Helpers/MathHelpers.cs

@@ -0,0 +1,71 @@
+using System;
+
+namespace OnnxStack.StableDiffusion.Helpers
+{
+    internal class MathHelpers
+    {
+        /// <summary>
+        /// Approximation of the definite integral of an analytic smooth function on a closed interval.
+        /// </summary>
+        /// <param name="function">The analytic smooth function to integrate.</param>
+        /// <param name="start">The start.</param>
+        /// <param name="end">The end.</param>
+        /// <param name="epsilon">The expected relative accuracy.</param>
+        /// <returns></returns>
+        public static float IntegrateOnClosedInterval(Func<double, double> function, double start, double end, double epsilon = 1e-4)
+        {
+            return (float)AdaptiveSimpson(function, start, end, epsilon);
+        }
+
+
+        /// <summary>
+        /// Initializes the adaptive Simpson's rule by calculating initial values and calling the auxiliary method.
+        /// </summary>
+        /// <param name="f">The f.</param>
+        /// <param name="a">a.</param>
+        /// <param name="b">The b.</param>
+        /// <param name="epsilon">The epsilon.</param>
+        /// <returns></returns>
+        private static double AdaptiveSimpson(Func<double, double> f, double a, double b, double epsilon)
+        {
+            double c = (a + b) / 2.0;
+            double h = b - a;
+            double fa = f(a);
+            double fb = f(b);
+            double fc = f(c);
+            double s = (h / 6) * (fa + 4 * fc + fb);
+            return AdaptiveSimpsonAux(f, a, b, epsilon, s, fa, fb, fc);
+        }
+
+
+        /// <summary>
+        /// Recursively applies the Simpson's rule and adapts the interval size based on the estimated error.
+        /// </summary>
+        /// <param name="f">The f.</param>
+        /// <param name="a">a.</param>
+        /// <param name="b">The b.</param>
+        /// <param name="epsilon">The epsilon.</param>
+        /// <param name="s">The s.</param>
+        /// <param name="fa">The fa.</param>
+        /// <param name="fb">The fb.</param>
+        /// <param name="fc">The fc.</param>
+        /// <returns></returns>
+        private static double AdaptiveSimpsonAux(Func<double, double> f, double a, double b, double epsilon, double s, double fa, double fb, double fc)
+        {
+            double c = (a + b) / 2.0;
+            double h = b - a;
+            double d = (a + c) / 2.0;
+            double e = (c + b) / 2.0;
+            double fd = f(d);
+            double fe = f(e);
+            double s1 = (h / 12) * (fa + 4 * fd + fc);
+            double s2 = (h / 12) * (fc + 4 * fe + fb);
+            double s_ = s1 + s2;
+            if (Math.Abs(s_ - s) <= 15 * epsilon)
+            {
+                return s_ + (s_ - s) / 15.0;
+            }
+            return AdaptiveSimpsonAux(f, a, c, epsilon / 2, s1, fa, fc, fd) + AdaptiveSimpsonAux(f, c, b, epsilon / 2, s2, fc, fb, fe);
+        }
+    }
+}

OnnxStack.StableDiffusion/OnnxStack.StableDiffusion.csproj

@@ -1,7 +1,7 @@
 <Project Sdk="Microsoft.NET.Sdk">
 
   <PropertyGroup>
-    <Version>0.31.0</Version>
+    <Version>0.31.16</Version>
     <TargetFramework>net7.0</TargetFramework>
     <ImplicitUsings>disable</ImplicitUsings>
     <Nullable>disable</Nullable>
@@ -50,12 +50,6 @@
     <ProjectReference Include="..\OnnxStack.Core\OnnxStack.Core.csproj" Condition=" '$(Configuration)' == 'Debug' OR '$(Configuration)'=='Debug-Nvidia'" />
   </ItemGroup>
 
-
-  <ItemGroup>
-    <PackageReference Include="MathNet.Numerics" Version="5.0.0" />
-    <PackageReference Include="NumSharp" Version="0.30.0" />
-  </ItemGroup>
-
   <ItemGroup>
     <None Update="README.md">
       <Pack>True</Pack>

OnnxStack.StableDiffusion/Schedulers/StableDiffusion/DDIMScheduler.cs

@@ -114,27 +114,26 @@ public override SchedulerStepResult Step(DenseTensor<float> modelOutput, int tim
             DenseTensor<float> predOriginalSample = null;
             if (Options.PredictionType == PredictionType.Epsilon)
             {
-                var sampleBeta = sample.SubtractTensors(modelOutput.MultiplyTensorByFloat((float)Math.Sqrt(betaProdT)));
-                predOriginalSample = sampleBeta.DivideTensorByFloat((float)Math.Sqrt(alphaProdT));
+                var sampleBeta = sample.SubtractTensors(modelOutput.MultiplyTensorByFloat(MathF.Sqrt(betaProdT)));
+                predOriginalSample = sampleBeta.DivideTensorByFloat(MathF.Sqrt(alphaProdT));
                 predEpsilon = modelOutput;
             }
             else if (Options.PredictionType == PredictionType.Sample)
             {
                 predOriginalSample = modelOutput;
                 predEpsilon = sample.SubtractTensors(predOriginalSample
-                    .MultiplyTensorByFloat((float)Math.Sqrt(alphaProdT)))
-                    .DivideTensorByFloat((float)Math.Sqrt(betaProdT));
+                    .MultiplyTensorByFloat(MathF.Sqrt(alphaProdT)))
+                    .DivideTensorByFloat(MathF.Sqrt(betaProdT));
             }
             else if (Options.PredictionType == PredictionType.VariablePrediction)
             {
-                var alphaSqrt = (float)Math.Sqrt(alphaProdT);
-                var betaSqrt = (float)Math.Sqrt(betaProdT);
-                predOriginalSample = sample
-                    .MultiplyTensorByFloat(alphaSqrt)
-                    .SubtractTensors(modelOutput.MultiplyTensorByFloat(betaSqrt));
-                predEpsilon = modelOutput
-                    .MultiplyTensorByFloat(alphaSqrt)
-                    .AddTensors(sample.MultiplyTensorByFloat(betaSqrt));
+                var tmp = sample.MultiplyTensorByFloat(MathF.Pow(alphaProdT, 0.5f));
+                var tmp2 = modelOutput.MultiplyTensorByFloat(MathF.Pow(betaProdT, 0.5f));
+                predOriginalSample = tmp.Subtract(tmp2);
+
+                var tmp3 = modelOutput.MultiplyTensorByFloat(MathF.Pow(alphaProdT, 0.5f));
+                var tmp4 = sample.MultiplyTensorByFloat(MathF.Pow(betaProdT, 0.5f));
+                predEpsilon = tmp3.Add(tmp4);
             }
 
 
@@ -154,24 +153,23 @@ public override SchedulerStepResult Step(DenseTensor<float> modelOutput, int tim
             //# σ_t = sqrt((1 − α_t−1)/(1 − α_t)) * sqrt(1 − α_t/α_t−1)
             var eta = 0f;
             var variance = GetVariance(currentTimestep, previousTimestep);
-            var stdDevT = eta * (float)Math.Sqrt(variance);
+            var stdDevT = eta * MathF.Pow(variance, 0.5f);
 
             var useClippedModelOutput = false;
             if (useClippedModelOutput)
             {
                 //# the pred_epsilon is always re-derived from the clipped x_0 in Glide
                 predEpsilon = sample
-                    .SubtractTensors(predOriginalSample.MultiplyTensorByFloat((float)Math.Sqrt(alphaProdT)))
-                    .DivideTensorByFloat((float)Math.Sqrt(betaProdT));
+                    .SubtractTensors(predOriginalSample.MultiplyTensorByFloat(MathF.Pow(alphaProdT, 0.5f)))
+                    .DivideTensorByFloat(MathF.Pow(betaProdT, 0.5f));
             }
 
 
             //# 5. compute "direction pointing to x_t" of formula (12) from https://arxiv.org/pdf/2010.02502.pdf
-            var predSampleDirection = predEpsilon.MultiplyTensorByFloat((float)Math.Sqrt(1f - alphaProdTPrev - Math.Pow(stdDevT, 2f)));
-
+            var predSampleDirection = predEpsilon.MultiplyTensorByFloat(MathF.Pow(1.0f - alphaProdTPrev - MathF.Pow(stdDevT, 2f), 0.5f));
 
             //# 6. compute x_t without "random noise" of formula (12) from https://arxiv.org/pdf/2010.02502.pdf
-            var prevSample = predSampleDirection.AddTensors(predOriginalSample.MultiplyTensorByFloat((float)Math.Sqrt(alphaProdTPrev)));
+            var prevSample = predSampleDirection.AddTensors(predOriginalSample.MultiplyTensorByFloat(MathF.Pow(alphaProdTPrev, 0.5f)));
 
             if (eta > 0)
                 prevSample = prevSample.AddTensors(CreateRandomSample(modelOutput.Dimensions).MultiplyTensorByFloat(stdDevT));
@@ -192,8 +190,8 @@ public override DenseTensor<float> AddNoise(DenseTensor<float> originalSamples,
             // Ref: https://github.com/huggingface/diffusers/blob/main/src/diffusers/schedulers/scheduling_ddpm.py#L456
             int timestep = timesteps[0];
             float alphaProd = _alphasCumProd[timestep];
-            float sqrtAlpha = (float)Math.Sqrt(alphaProd);
-            float sqrtOneMinusAlpha = (float)Math.Sqrt(1.0f - alphaProd);
+            float sqrtAlpha = MathF.Sqrt(alphaProd);
+            float sqrtOneMinusAlpha = MathF.Sqrt(1.0f - alphaProd);
 
             return noise
                 .MultiplyTensorByFloat(sqrtOneMinusAlpha)

OnnxStack.StableDiffusion/Schedulers/StableDiffusion/DDPMScheduler.cs

@@ -1,5 +1,4 @@
 using Microsoft.ML.OnnxRuntime.Tensors;
-using NumSharp;
 using OnnxStack.Core;
 using OnnxStack.StableDiffusion.Config;
 using OnnxStack.StableDiffusion.Enums;
@@ -115,7 +114,7 @@ public override SchedulerStepResult Step(DenseTensor<float> modelOutput, int tim
             if (Options.Thresholding)
             {
                 // TODO: https://github.com/huggingface/diffusers/blob/main/src/diffusers/schedulers/scheduling_ddpm.py#L322
-                predOriginalSample = ThresholdSample(predOriginalSample);
+               // predOriginalSample = ThresholdSample(predOriginalSample);
             }
             else if (Options.ClipSample)
             {
@@ -255,87 +254,6 @@ private float GetVariance(int timestep, float predictedVariance = 0f)
         }
 
 
-        /// <summary>
-        /// Thresholds the sample.
-        /// </summary>
-        /// <param name="input">The input.</param>
-        /// <param name="dynamicThresholdingRatio">The dynamic thresholding ratio.</param>
-        /// <param name="sampleMaxValue">The sample maximum value.</param>
-        /// <returns></returns>
-        private DenseTensor<float> ThresholdSample(DenseTensor<float> input, float dynamicThresholdingRatio = 0.995f, float sampleMaxValue = 1f)
-        {
-            var sample = new NDArray(input.ToArray(), new Shape(input.Dimensions.ToArray()));
-            var batch_size = sample.shape[0];
-            var channels = sample.shape[1];
-            var height = sample.shape[2];
-            var width = sample.shape[3];
-
-            // Flatten sample for doing quantile calculation along each image
-            var flatSample = sample.reshape(batch_size, channels * height * width);
-
-            // Calculate the absolute values of the sample
-            var absSample = np.abs(flatSample);
-
-            // Calculate the quantile for each row
-            var quantiles = new List<float>();
-            for (int i = 0; i < batch_size; i++)
-            {
-                var row = absSample[$"{i},:"].MakeGeneric<float>();
-                var percentileValue = CalculatePercentile(row, dynamicThresholdingRatio);
-                percentileValue = Math.Clamp(percentileValue, 1f, sampleMaxValue);
-                quantiles.Add(percentileValue);
-            }
-
-            // Create an NDArray from quantiles
-            var quantileArray = np.array(quantiles.ToArray());
-
-            // Calculate the thresholded sample
-            var sExpanded = np.expand_dims(quantileArray, 1); // Expand to match the sample shape
-            var negSExpanded = np.negative(sExpanded); // Get the negation of sExpanded
-            var thresholdedSample = sample - negSExpanded; // Element-wise subtraction
-            thresholdedSample = np.maximum(thresholdedSample, negSExpanded); // Ensure values are not less than -sExpanded
-            thresholdedSample = np.minimum(thresholdedSample, sExpanded); // Ensure values are not greater than sExpanded
-            thresholdedSample = thresholdedSample / sExpanded;
-
-            // Reshape to the original shape
-            thresholdedSample = thresholdedSample.reshape(batch_size, channels, height, width);
-
-            return new DenseTensor<float>(thresholdedSample.ToArray<float>(), thresholdedSample.shape);
-        }
-
-
-        /// <summary>
-        /// Calculates the percentile.
-        /// </summary>
-        /// <param name="data">The data.</param>
-        /// <param name="percentile">The percentile.</param>
-        /// <returns></returns>
-        private float CalculatePercentile(NDArray data, float percentile)
-        {
-            // Sort the data indices in ascending order
-            var sortedIndices = np.argsort<float>(data);
-
-            // Calculate the index corresponding to the percentile
-            var index = (int)Math.Ceiling(percentile / 100f * (data.Shape[0] - 1));
-
-            // Retrieve the value at the calculated index
-            var percentileValue = data[sortedIndices[index]];
-
-            return percentileValue.GetSingle();
-        }
-
-
-        /// <summary>
-        /// Determines whether the VarianceType is learned.
-        /// </summary>
-        /// <returns>
-        ///   <c>true</c> if the VarianceType is learned; otherwise, <c>false</c>.
-        /// </returns>
-        private bool IsVarianceTypeLearned()
-        {
-            return Options.VarianceType == VarianceType.Learned || Options.VarianceType == VarianceType.LearnedRange;
-        }
-
         protected override void Dispose(bool disposing)
         {
             _alphasCumProd = null;

OnnxStack.StableDiffusion/Schedulers/StableDiffusion/EulerAncestralScheduler.cs

@@ -1,8 +1,6 @@
 using Microsoft.ML.OnnxRuntime.Tensors;
-using NumSharp;
 using OnnxStack.Core;
 using OnnxStack.StableDiffusion.Config;
-using OnnxStack.StableDiffusion.Enums;
 using OnnxStack.StableDiffusion.Helpers;
 using System;
 using System.Collections.Generic;
@@ -55,14 +53,14 @@ protected override int[] SetTimesteps()
         {
             var sigmas = _sigmas.ToArray();
             var timesteps = GetTimesteps();
-            var log_sigmas = np.log(sigmas).ToArray<float>();
-            var range = np.arange(0, (float)_sigmas.Length).ToArray<float>();
+            var logSigmas = ArrayHelpers.Log(sigmas);
+            var range = ArrayHelpers.Range(0, _sigmas.Length);
             sigmas = Interpolate(timesteps, range, _sigmas);
 
             if (Options.UseKarrasSigmas)
             {
                 sigmas = ConvertToKarras(sigmas);
-                timesteps = SigmaToTimestep(sigmas, log_sigmas);
+                timesteps = SigmaToTimestep(sigmas, logSigmas);
             }
 
             _sigmas = sigmas

OnnxStack.StableDiffusion/Schedulers/StableDiffusion/EulerScheduler.cs

@@ -1,5 +1,4 @@
 using Microsoft.ML.OnnxRuntime.Tensors;
-using NumSharp;
 using OnnxStack.Core;
 using OnnxStack.StableDiffusion.Config;
 using OnnxStack.StableDiffusion.Enums;
@@ -55,19 +54,20 @@ protected override int[] SetTimesteps()
         {
             var sigmas = _sigmas.ToArray();
             var timesteps = GetTimesteps();
-            var log_sigmas = np.log(sigmas).ToArray<float>();
-            var range = np.arange(0, (float)_sigmas.Length).ToArray<float>();
+            var logSigmas = ArrayHelpers.Log(sigmas);
+            var range = ArrayHelpers.Range(0, sigmas.Length);
 
             // TODO: Implement "interpolation_type"
-            var interpolation_type = "linear";
-            sigmas = interpolation_type == "log_linear"
-                ? np.exp(np.linspace(np.log(sigmas.Last()), np.log(sigmas.First()), timesteps.Length + 1)).ToArray<float>()
-                : Interpolate(timesteps, range, _sigmas);
+            //var interpolation_type = "linear";
+            //sigmas = interpolation_type == "log_linear"
+            //    ? np.exp(np.linspace(np.log(sigmas.Last()), np.log(sigmas.First()), timesteps.Length + 1)).ToArray<float>()
+            //    : Interpolate(timesteps, range, _sigmas);
 
+            sigmas = Interpolate(timesteps, range, _sigmas);
             if (Options.UseKarrasSigmas)
             {
                 sigmas = ConvertToKarras(sigmas);
-                timesteps = SigmaToTimestep(sigmas, log_sigmas);
+                timesteps = SigmaToTimestep(sigmas, logSigmas);
             }
 
             _sigmas = sigmas

OnnxStack.StableDiffusion/Schedulers/StableDiffusion/KDPM2Scheduler.cs

@@ -1,8 +1,8 @@
 using Microsoft.ML.OnnxRuntime.Tensors;
-using NumSharp;
 using OnnxStack.Core;
 using OnnxStack.StableDiffusion.Config;
 using OnnxStack.StableDiffusion.Enums;
+using OnnxStack.StableDiffusion.Helpers;
 using System;
 using System.Collections.Generic;
 using System.Linq;
@@ -63,8 +63,8 @@ protected override int[] SetTimesteps()
             // Create timesteps based on the specified strategy
             var sigmas = _sigmas.ToArray();
             var timesteps = GetTimesteps();
-            var logSigmas = np.log(sigmas).ToArray<float>();
-            var range = np.arange(0, (float)_sigmas.Length).ToArray<float>();
+            var logSigmas = ArrayHelpers.Log(sigmas);
+            var range = ArrayHelpers.Range(0, _sigmas.Length);
             sigmas = Interpolate(timesteps, range, _sigmas);
 
             if (Options.UseKarrasSigmas)