Dev/working (#59)

* Improve special function performance and increase resiliance, especially to NaNs. Improve Poisson and Beta random deviate performance.

* UInt128, Int128 working. DoubleDouble improved. Complete Elliptic Inegral of 3rd kind. Improved Gamma speed and accuracy using series aound 1 and 2. Explicit ScaledModifiedBessel and SphericalBessel.

* Mostly tests and documentation. RegressionResult layer.

Numerics/Analysis/MultiFunctionMath_Extrema_ModelTrust.cs

@@ -141,20 +141,22 @@ private static MultiExtremum FindMinimum_ModelTrust (MultiFunctor f, IReadOnlyLi
                 double value = f.Evaluate(point);
 
                 double delta = model.MinimumValue - value;
-                double tol = settings.ComputePrecision(value);
+                double tol = settings.ComputePrecision(Math.Min(model.MinimumValue, value));
+                // Note value can be way off, so use better of old best and new value to compute tol.
+                // When we didn't do this before, we got value = infinity, so tol = infinity, and thus terminated!
 
-                if (delta > 0 && settings.Listener != null) {
+                if (delta > 0.0 && settings.Listener != null) {
                     MultiExtremum report = new MultiExtremum(f.EvaluationCount, settings, point, value, Math.Max(Math.Abs(delta), 0.75 * tol), model.GetHessian());
                     settings.Listener(report);
                 }
 
                 // To terminate, we demand: a reduction, that the reduction be small, that the reduction be in line with
-                // its expected value, that we have not run up against trust boundary, and that the gradient is small.
+                // its expected value, that we have not run up against the trust boundary, and that the gradient is small.
                 // I had wanted to demand delta > 0, but we run into some cases where delta keeps being very slightly
                 // negative, typically orders of magnitude less than tol, causing the trust radius to shrink in an
                 // endless cycle that causes our approximation to ultimately go sour, even though terminating on the original
                 // very slightly negative delta would have produced an accurate estimate. So we tolerate this case for now.
-                if ((-tol / 4.0 <= delta) && (delta <= tol)) {
+                if ((delta <= tol) && (-0.25 * tol <= delta)) {
                     // We demand that the model be decent, i.e. that the expected delta was within tol of the measured delta.
                     if (Math.Abs(delta - deltaExpected) <= tol) {
                         // We demand that the step not just be small because it ran up against the trust radius.
@@ -181,7 +183,7 @@ private static MultiExtremum FindMinimum_ModelTrust (MultiFunctor f, IReadOnlyLi
                 // If the actual change was very far from the expected change, reduce the trust radius.
                 // If the expected change did a good job of predicting the actual change, increase the trust radius.
                 if ((delta < 0.25 * deltaExpected) /*|| (8.0 * deltaExpected < delta)*/) {
-                    trustRadius = trustRadius / 2.0;
+                    trustRadius = 0.5 * trustRadius;
                 } else if ((0.75 * deltaExpected <= delta) /*&& (delta <= 2.0 * deltaExpected)*/) {
                     trustRadius = 2.0 * trustRadius;
                 }
@@ -195,8 +197,9 @@ private static MultiExtremum FindMinimum_ModelTrust (MultiFunctor f, IReadOnlyLi
                     double bad = model.ComputeBadness(i, z, point, value);
                     if (bad > fBad) { iBad = i; fBad = bad; }
                 }
+                // Use the new point as long as it is better than our worst existing point.
                 if (value < fMax) {
-                    //Debug.WriteLine("iMax={0}, iBad={1}", iMax, iBad);
+                    Debug.Assert(!Double.IsPositiveInfinity(value) && !Double.IsNaN(value));
                     model.ReplacePoint(iBad, point, z, value);
                 }
                 // There is some question about how best to choose which point to replace.
@@ -710,21 +713,28 @@ public double[] ConvertPoint (double[] z) {
             return (x);
         }
 
+        // Badness is used as a criteria for which point to throw out of collection.
+        // Basic idea is to compute a number which (i) increases with value above
+        // minimum and also (ii) increases with distance from minimum.
+
         public double ComputeBadness (int index, double[] z, double[] point, double value) {
             double s = 0.0;
             if (value < MinimumValue) {
+                // If new candidate is best, compute badness relative to it.
                 for (int i = 0; i < point.Length; i++) {
                     s += MoreMath.Sqr(points[index][i] - point[i]);
                 }
                 s = Math.Pow(s, 3.0 / 2.0) * (values[index] - value);
             } else {
+                // Otherwise compute badness relative to best existing point.
                 if (index == minValueIndex) return (0.0);
                 for (int i = 0; i < point.Length; i++) {
                     s += MoreMath.Sqr(points[index][i] - origin[i]);
                 }
                 s = Math.Pow(s, 3.0 / 2.0) * (values[index] - values[minValueIndex]);
             }
-            return (s);
+            Debug.Assert(s >= 0.0);
+            return s;
         }
 
         public double FindMinimumSeperation () {

Numerics/Analysis/NamespaceDoc.cs

@@ -3,7 +3,7 @@
 namespace Meta.Numerics.Analysis {
 
     /// <summary>
-    /// Contains types used to analyze user-supplied functions.
+    /// Contains types used to solve, maximize, integrate, and otherwise perform analysis on user-supplied functions.
     /// </summary>
     /// <remarks>
     /// <para>This namespace contains types for the numerical analysis of user-supplied functions. Examples

Numerics/Core/DiscreteInterval.cs

@@ -1,8 +1,6 @@
 using System;
-using System.Collections.Generic;
-using System.Linq;
-using System.Text;
-using System.Threading.Tasks;
+using System.Diagnostics;
+using System.Globalization;
 
 namespace Meta.Numerics {
 
@@ -12,6 +10,7 @@ namespace Meta.Numerics {
     public struct DiscreteInterval : IEquatable<DiscreteInterval> {
 
         internal DiscreteInterval (int min, int max) {
+            Debug.Assert(max >= min);
             this.min = min;
             this.max = max;
         }
@@ -63,31 +62,68 @@ internal int Width {
 
         // Equality
 
-        public bool Equals (DiscreteInterval other) {
-            return ((this.min == other.min) && (this.max == other.max));
+        private static bool Equals (DiscreteInterval u, DiscreteInterval v) {
+            return (u.min == v.min) && (u.max == v.max);
         }
 
+
+        /// <summary>
+        /// Tests whether two discrete intervals are equal.
+        /// </summary>
+        /// <param name="a">The first discrete interval.</param>
+        /// <param name="b">The second discrete interval.</param>
+        /// <returns><see langword="true"/> if <paramref name="a"/> and <paramref name="b"/> are equal, otherwise <see langword="false"/>.</returns>
         public static bool operator == (DiscreteInterval a, DiscreteInterval b) {
-            return (a.Equals(b));
+            return Equals(a, b);
         }
 
+        /// <summary>
+        /// Tests whether two intervals are not equal.
+        /// </summary>
+        /// <param name="a">The first interval.</param>
+        /// <param name="b">The second interval.</param>
+        /// <returns><see langword="false"/> if <paramref name="a"/> and <paramref name="b"/> are equal, otherwise <see langword="true"/>.</returns>
         public static bool operator != (DiscreteInterval a, DiscreteInterval b) {
-            return (!a.Equals(b));
+            return !Equals(a, b);
         }
 
+        /// <summary>
+        /// Tests wither the current instance is equal to another discrete interval.
+        /// </summary>
+        /// <param name="other">Another discrete interval.</param>
+        /// <returns><see langword="true"/> if the current instance equals <paramref name="other"/>, otherwise <see langword="false"/>.</returns>
+        public bool Equals (DiscreteInterval other) {
+            return Equals(this, other);
+        }
+
+        /// <summary>
+        /// Tests whether a given object is equal to the current interval.
+        /// </summary>
+        /// <param name="obj">An object.</param>
+        /// <returns><see langword="true"/> if <paramref name="obj"/> is an equal <see cref="DiscreteInterval"/>, otherwise <see langword="false"/>.</returns>
         public override bool Equals (object obj) {
             if (obj is DiscreteInterval other) {
-                return (this.Equals(other));
+                return Equals(this, other);
             } else {
-                return (false);
+                return false;
             }
         }
 
+        /// <inheritdoc />
         public override int GetHashCode () {
             unchecked {
                 return (min.GetHashCode() + 31 * max.GetHashCode());
             }
         }
 
+        /// <inheritdoc />
+        public override string ToString () {
+            return ToString(CultureInfo.CurrentCulture);
+        }
+
+        private string ToString (IFormatProvider format) {
+            return String.Format(format, "{{{0}..{1}}}", min, max);
+        }
+
     }
 }

Numerics/Core/Interval.cs

@@ -13,7 +13,7 @@ namespace Meta.Numerics {
     /// </remarks>
     public struct Interval : IEquatable<Interval> {
 
-        private double a, b, w;
+        private readonly double a, b, w;
 
         private Interval (double a, double b, double w) {
             this.a = a;
@@ -120,55 +120,67 @@ public static Interval FromMidpointAndWidth (double midpoint, double width) {
 
         // equality
 
+        private static bool Equals (Interval u, Interval v) {
+            return ((u.a == v.a) && (u.w == v.w));
+        }
+
+        /// <summary>
+        /// Tests wither the current instance is equal to another interval.
+        /// </summary>
+        /// <param name="other">Another iteterval.</param>
+        /// <returns><see langword="true"/> if the current instance equals <paramref name="other"/>, otherwise <see langword="false"/>.</returns>
         public bool Equals (Interval other) {
-            return ((this.a == other.a) && (this.w == other.w));
+            return Equals(this, other);
         }
 
         /// <summary>
-        /// Determines whether two intervals are equal.
+        /// Tests whether two intervals are equal.
         /// </summary>
-        /// <param name="i1">The first interval.</param>
-        /// <param name="i2">The second interval.</param>
-        /// <returns>True if <paramref name="i1"/> and <paramref name="i2"/> are equal, otherwise false.</returns>
-        public static bool operator == (Interval i1, Interval i2) {
-            return (i1.Equals(i2));
+        /// <param name="u">The first interval.</param>
+        /// <param name="v">The second interval.</param>
+        /// <returns><see langword="true"/> if <paramref name="u"/> and <paramref name="v"/> are equal, otherwise <see langword="false"/>.</returns>
+        public static bool operator == (Interval u, Interval v) {
+            return Equals(u, v);
         }
 
         /// <summary>
-        /// Determines whether two intervals are not equal.
+        /// Tests whether two intervals are not equal.
         /// </summary>
-        /// <param name="i1">The first interval.</param>
-        /// <param name="i2">The second interval.</param>
-        /// <returns>True if <paramref name="i1"/> and <paramref name="i2"/> are not equal, otherwise false.</returns>
-        public static bool operator != (Interval i1, Interval i2) {
-            return (!i1.Equals(i2));
+        /// <param name="u">The first interval.</param>
+        /// <param name="v">The second interval.</param>
+        /// <returns><see langword="false"/> if <paramref name="u"/> and <paramref name="v"/> are equal, otherwise <see langword="true"/>.</returns>
+        public static bool operator != (Interval u, Interval v) {
+            return !Equals(u, v);
         }
 
         /// <summary>
-        /// Determines whether a given object is an equal interval.
+        /// Tests whether a given object is equal to the current interval.
         /// </summary>
         /// <param name="obj">An object.</param>
-        /// <returns>True if <paramref name="obj"/> is an equal <see cref="Interval"/>, otherwise false.</returns>
+        /// <returns><see langword="true"/> if <paramref name="obj"/> is an equal <see cref="Interval"/>, otherwise <see langword="false"/>.</returns>
         public override bool Equals (object obj) {
             if (obj is Interval other) {
-                return(this.Equals(other));
+                return Equals(this, other);
             } else {
-                return (false);
+                return false;
             }
         }
 
         /// <inheritdoc />
         public override int GetHashCode () {
-            return (a.GetHashCode() ^ w.GetHashCode());
+            return a.GetHashCode() ^ w.GetHashCode();
         }
 
-        // text representation
         /// <summary>
         /// Produces a string representation of the interval.
         /// </summary>
         /// <returns>A string representation of the interval.</returns>
         public override string ToString () {
-            return (String.Format(CultureInfo.CurrentCulture, "[{0},{1}]", a, b));
+            return ToString(CultureInfo.CurrentCulture);
+        }
+
+        private string ToString (IFormatProvider format) {
+            return String.Format(CultureInfo.CurrentCulture, "[{0},{1}]", a, b);
         }
 
 #if SHO

Numerics/Core/Polynomial.cs

@@ -17,7 +17,7 @@ public class Polynomial {
         /// <returns>The specified polynomial.</returns>
         /// <remarks>
         /// <para>Coefficients should be arranged from low to high order, so that the kth entry is the coefficient of x<sup>k</sup>. For example,
-        /// to specify the polynomial 5 - 6 x + 7 x<sup>2</sup>, give the values 5, -6, 7.</para>
+        /// to specify the polynomial 5 - 6 x + 7 x<sup>3</sup>, give the values 5, -6, , 0, 7.</para>
         /// </remarks>
         public static Polynomial FromCoefficients (params double[] coefficients) {
             if (coefficients == null) throw new ArgumentNullException(nameof(coefficients));

Numerics/Data/FrameColumn.cs

@@ -17,6 +17,8 @@ namespace Meta.Numerics.Data
     public sealed class FrameColumn : IEnumerable {
 
         internal FrameColumn(FrameView frame, int c) {
+            Debug.Assert(frame != null);
+            Debug.Assert(c >= 0 && c < frame.Columns.Count);
             this.column = frame.columns[c];
             this.map = frame.map;
         }
@@ -77,8 +79,7 @@ public object this[int r] {
         /// </remarks>
         public IReadOnlyList<T> As<T> () {
             // If the requested column is of the requested type, expose it directly.
-            IReadOnlyList<T> typedColumn = column as IReadOnlyList<T>;
-            if (typedColumn != null) {
+            if (column is IReadOnlyList<T> typedColumn) {
                 return (new TypedFrameColumn<T>(column.Name, typedColumn, map));
             } else {
                 // If the requested column is cast-able to the requested type, cast it.

Numerics/Data/FrameColumnCollection.cs

@@ -33,10 +33,11 @@ public int Count {
         /// </summary>
         /// <param name="index">The index of the column.</param>
         /// <returns>The column with index <paramref name="index"/>.</returns>
+        /// <exception cref="IndexOutOfRangeException"><paramref name="index"/> is less than zero or greater than or equal to the number of rows <see cref="Count"/>.</exception>
         public FrameColumn this [int index] {
             get {
-                if ((index < 0) || (index >= frame.columns.Count)) throw new ArgumentOutOfRangeException(nameof(index));
-                return (new FrameColumn(frame, index));
+                if ((index < 0) || (index >= frame.columns.Count)) throw new IndexOutOfRangeException();
+                return new FrameColumn(frame, index);
             }
         }
 
@@ -45,10 +46,12 @@ public FrameColumn this [int index] {
         /// </summary>
         /// <param name="name">The name of the column.</param>
         /// <returns>The column with the given name.</returns>
+        /// <exception cref="IndexOutOfRangeException">There is no column with the name <paramref name="name"/> in the collection.</exception>
         public FrameColumn this [string name] {
             get {
                 int index = frame.GetColumnIndex(name);
-                return (new FrameColumn(frame, index));
+                if (index < 0) throw new IndexOutOfRangeException();
+                return new FrameColumn(frame, index);
             }
         }
 

Numerics/Data/FrameRow.cs

@@ -52,9 +52,10 @@ IEnumerable<object> IReadOnlyDictionary<string, object>.Values {
         /// </summary>
         /// <param name="c">The index of the column.</param>
         /// <returns>The value of the <paramref name="c"/>th column.</returns>
+        /// <exception cref="IndexOutOfRangeException"><paramref name="c"/> is less than zero or greater than or equal to the number of columns <see cref="Count"/>.</exception>
         public object this[int c] {
             get {
-                if ((c < 0) || (c >= Count)) throw new ArgumentOutOfRangeException(nameof(c));
+                if ((c < 0) || (c >= Count)) throw new IndexOutOfRangeException();
                 return (frame.columns[c].GetItem(frame.map[r]));
             }
         }
@@ -64,6 +65,7 @@ public object this[int c] {
         /// </summary>
         /// <param name="name">The name of the column.</param>
         /// <returns>The value.</returns>
+        /// <exception cref="IndexOutOfRangeException">There is no column in the row with the given <paramref name="name"/>.</exception>
         public object this[string name] {
             get {
                 int c = frame.GetColumnIndex(name);

Numerics/Extended/AdvancedDoubleDoubleMath_Gamma.cs

@@ -132,7 +132,7 @@ private static DoubleDouble LogGamma_ShiftToZetaSeries (DoubleDouble x) {
             -((DoubleDouble) 74611) / 330, ((DoubleDouble) 854513) / 138, -((DoubleDouble) 236364091) / 2730, ((DoubleDouble) 8553103) / 6, -((DoubleDouble) 23749461029) / 870
         };
 
-        public static DoubleDouble BernoulliSum (DoubleDouble x) {
+        private static DoubleDouble BernoulliSum (DoubleDouble x) {
             DoubleDouble rxPower = 1.0 / x;
             DoubleDouble rxSquared = rxPower * rxPower;
             DoubleDouble f = 0.5 * Bernoulli[1] * rxPower;