GaussianOp.PrecisionAverageConditional_Point correctly handles corner cases

src/Runtime/Distributions/Gamma.cs

@@ -319,12 +319,26 @@ namespace Microsoft.ML.Probabilistic.Distributions
         {
             if (x <= 0)
                 throw new ArgumentException("x <= 0");
-            double a = -x * x * ddLogP;
-            if (ddLogP == 0.0)
-                a = 0.0;  // in case x is infinity
+            double a;
+            if(double.IsPositiveInfinity(x))
+            {
+                if (ddLogP < 0) return Gamma.PointMass(x);
+                else if (ddLogP == 0) a = 0.0;
+                else if (forceProper)
+                {
+                    if (dLogP <= 0) return Gamma.FromShapeAndRate(1, -dLogP);
+                    else return Gamma.PointMass(x);
+                }
+                else return Gamma.FromShapeAndRate(-x, -x - dLogP);
+            }
+            else
+            {
+                a = -x * x * ddLogP;
+                if (a+1 > double.MaxValue) return Gamma.PointMass(x);
+            }
             if (forceProper)
             {
-                if (dLogP < 0)
+                if (dLogP <= 0)
                 {
                     if (a < 0)
                         a = 0;
@@ -333,7 +347,9 @@ namespace Microsoft.ML.Probabilistic.Distributions
                 {
                     double amin = x * dLogP;
                     if (a < amin)
-                        a = amin;
+                    {
+                        return Gamma.FromShapeAndRate(amin + 1, 0);
+                    }
                 }
             }
             double b = a / x - dLogP;

src/Runtime/Factors/GaussianOp.cs

@@ -502,7 +502,7 @@ namespace Microsoft.ML.Probabilistic.Factors
             // x2ddlogf = (-yv*r-0.5+r*ym*ym)/(yv*r+1)^2 - r*ym*ym/(yv*r+1)^3
             // xdlogf + x2ddlogf = (-0.5*yv*r + 0.5*r*ym*ym)/(yv*r+1)^2 - r*ym*ym/(yv*r+1)^3
             // as r->0: -0.5*r*(yv + ym*ym)
-            if (precision == 0) return Gamma.FromShapeAndRate(1.5, 0.5 * (yv + ym * ym));
+            if (precision == 0 || yv == 0) return Gamma.FromShapeAndRate(1.5, 0.5 * (yv + ym * ym));
             // point mass case
             // f(r) = N(xm;mm, xv+mv+1/r)
             // log f(r) = -0.5*log(xv+mv+1/r) - 0.5*(xm-mm)^2/(xv+mv+1/r)
@@ -512,7 +512,7 @@ namespace Microsoft.ML.Probabilistic.Factors
             // r^2 (log f)'' = -1/(yv*r + 1) + r*ym^2/(yv*r+1)^2) + 0.5/(yv*r+1)^2 - r*ym^2/(yv*r+1)^3
             double vdenom = 1 / (yv * precision + 1);
             double ymvdenom = ym * vdenom;
-            double ymvdenom2 = precision * ymvdenom * ymvdenom;
+            double ymvdenom2 = (precision > double.MaxValue) ? 0 : (precision * ymvdenom * ymvdenom);
             //dlogf = (-0.5 * denom + 0.5 * ym2denom2) * (-v2);
             //dlogf = 0.5 * (1 - ym * ymdenom) * denom * v2;
             //dlogf = 0.5 * (v - ym * ym/(yv*precision+1))/(yv*precision + 1);

test/Tests/Operators/OperatorTests.cs

@@ -103,6 +103,21 @@ namespace Microsoft.ML.Probabilistic.Tests
             Assert.True(double.IsPositiveInfinity(a) || MMath.NextDouble(a) - b > sum);
         }
 
+        [Fact]
+        public void PrecisionAverageConditional_Point_IsIncreasing()
+        {
+            foreach (var precision in DoublesAtLeastZero())
+            {
+                foreach (var yv in DoublesAtLeastZero())
+                {
+                    var result0 = GaussianOp.PrecisionAverageConditional_Point(0, yv, precision);
+                    var result = GaussianOp.PrecisionAverageConditional_Point(MMath.NextDouble(0), yv, precision);
+                    Assert.True(result.Rate >= result0.Rate);
+                    //Trace.WriteLine($"precision={precision} yv={yv}: {result0.Rate} {result.Rate}");
+                }
+            }
+        }
+
         // Test inference on a model where precision is scaled.
         internal void GammaProductTest()
         {
@@ -2063,7 +2078,7 @@ zL = (L - mx)*sqrt(prec)
             yield return MMath.NextDouble(0);
             yield return MMath.PreviousDouble(double.PositiveInfinity);
             yield return double.PositiveInfinity;
-            for (int i = 0; i <= 100; i++)
+            for (int i = 0; i <= 300; i++)
             {
                 double bigValue = System.Math.Pow(10, i);
                 yield return -bigValue;