//-- Constant bounds --------------------------------------------------------------------------------

		/// <summary>
		/// The logarithm of the probability that L &lt;= X &lt; U.
		/// </summary>
		/// <param name="X"></param>
		/// <param name="L">Can be negative infinity.</param>
		/// <param name="U">Can be positive infinity.</param>
		/// <returns></returns>
		public static double LogProbBetween(Gaussian X, double L, double U)
		{
			if (L > U) throw new AllZeroException("low > high (" + L + " > " + U + ")");
			if (X.IsPointMass)
			{
				return Factor.IsBetween(X.Point, L, U) ? 0.0 : Double.NegativeInfinity;
			}
			else if (X.IsUniform())
			{
				if (Double.IsNegativeInfinity(L))
				{
					if (Double.IsPositiveInfinity(U)) return 0.0; // always between
					else return -MMath.Ln2;  // between half the time
				}
				else if (Double.IsPositiveInfinity(U)) return -MMath.Ln2;  // between half the time
				else return Double.NegativeInfinity;  // never between two finite numbers
			}
			else
			{
				double sqrtPrec = Math.Sqrt(X.Precision);
				double mx = X.GetMean();
				double pl = MMath.NormalCdfLn(sqrtPrec * (L - mx));  // log(p(x <= L))
				double pu = MMath.NormalCdfLn(sqrtPrec * (U - mx));  // log(p(x <= U))
				if (pl == pu) return Double.NegativeInfinity;
				if (Double.IsNegativeInfinity(pl)) return pu;
				// log(NormalCdf(yu) - NormalCdf(yl)) = NormalCdfLn(yu) + log(1 - NormalCdf(yl)/NormalCdf(yu))
				return pu + MMath.Log1MinusExp(pl - pu);
			}
		}

		// logw1 = N(mx;m1,vx+v1) phi((mx1 - m2)/sqrt(vx1+v2))
		// a1 = N(mx1;m2,vx1+v2)/phi
		internal static void ComputeStats(Gaussian max, Gaussian a, Gaussian b, out double logz,
			out double logw1, out double a1, out double vx1, out double mx1,
			out double logw2, out double a2, out double vx2, out double mx2)
		{
			double m1, v1, m2, v2;
			a.GetMeanAndVariance(out m1, out v1);
			b.GetMeanAndVariance(out m2, out v2);
			if (max.IsPointMass) {
				vx1 = 0.0;
				mx1 = max.Point;
				vx2 = 0.0;
				mx2 = max.Point;
				if (b.IsPointMass) {
					if (b.Point > max.Point) throw new AllZeroException();
					else if (b.Point == max.Point) {
						// the factor reduces to the constraint (max.Point > a)
						logw1 = Double.NegativeInfinity;
						logw2 = MMath.NormalCdfLn((max.Point - m1)/Math.Sqrt(v1));
						logz = logw2;
						a1 = 0;
						a2 = Math.Exp(Gaussian.GetLogProb(max.Point, m1, v1) - logw2);
						return;
					} else {
						// b.Point < max.Point
						// the factor reduces to the constraint (a == max.Point)
						throw new NotImplementedException();
					}
				} else if (a.IsPointMass) throw new NotImplementedException();
			} else {
				if (a.IsPointMass) {
					vx1 = 0.0;
					mx1 = a.Point;
				} else {
					vx1 = 1.0 / (max.Precision + a.Precision);
					mx1 = vx1 * (max.MeanTimesPrecision + a.MeanTimesPrecision);
				}
				if (b.IsPointMass) {
					vx2 = 0.0;
					mx2 = b.Point;
				} else {
					vx2 = 1.0 / (max.Precision + b.Precision);
					mx2 = vx2 * (max.MeanTimesPrecision + b.MeanTimesPrecision);
				}
			}
			logw1 = max.GetLogAverageOf(a);
			double logPhi1 = MMath.NormalCdfLn((mx1 - m2) / Math.Sqrt(vx1 + v2));
			logw1 += logPhi1;

			logw2 = max.GetLogAverageOf(b);
			double logPhi2 = MMath.NormalCdfLn((mx2 - m1) / Math.Sqrt(vx2 + v1));
			logw2 += logPhi2;

			logz = MMath.LogSumExp(logw1, logw2);

			double logN1 = Gaussian.GetLogProb(mx1, m2, vx1 + v2);
			a1 = Math.Exp(logN1 - logPhi1);
			double logN2 = Gaussian.GetLogProb(mx2, m1, vx2 + v1);
			a2 = Math.Exp(logN2 - logPhi2);
		}

		public Gaussian[] InferTomorrowsTime()
		{
			Gaussian[] tomorrowsTime = new Gaussian[2];

			tomorrowsTime[0] = cyclist1.InferTomorrowsTime();
			tomorrowsTime[1] = cyclist2.InferTomorrowsTime();
			return tomorrowsTime;
		}

 public void Base_ExpectedValue()
 {
     var Mu = 5.0;
     var X = new Gaussian(Mu, 2.0);
     var E = X.ExpectedValueWithConfidence().Mean;
     var Err = Math.Abs(E - Mu) / Mu;
     Assert.IsTrue(Err < eps);
 }

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="logistic">Constant value for 'logistic'.</param>
		/// <param name="x">Incoming message from 'x'.</param>
		/// <returns>Logarithm of the factor's average value across the given argument distributions</returns>
		/// <remarks><para>
		/// The formula for the result is <c>log(sum_(x) p(x) factor(logistic,x))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(double logistic, Gaussian x)
		{
			if (logistic >= 1.0 || logistic <= 0.0) return x.GetLogProb(MMath.Logit(logistic));
			// p(y,x) = delta(y - 1/(1+exp(-x))) N(x;mx,vx)
			// x = log(y/(1-y))
			// dx = 1/(y*(1-y))
			return x.GetLogProb(MMath.Logit(logistic))/(logistic*(1-logistic));
		}

		/// <summary>
		/// VMP message to 'product'
		/// </summary>
		/// <param name="A">Constant value for 'a'.</param>
		/// <param name="B">Incoming message from 'b'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="result">Modified to contain the outgoing message</param>
		/// <returns><paramref name="result"/></returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'product' as the random arguments are varied.
		/// The formula is <c>proj[sum_(b) p(b) factor(product,a,b)]</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
		public static Gaussian ProductAverageLogarithm(double A, [SkipIfUniform] Beta B)
		{
			double mb, vb;
			B.GetMeanAndVariance(out mb, out vb);	
			Gaussian result = new Gaussian();
			result.SetMeanAndVariance(A * mb, A * A * vb);
			return result;
		}

		/// <summary>
		/// VMP message to 'innerProduct'
		/// </summary>
		/// <param name="AMean">Buffer 'AMean'.</param>
		/// <param name="AVariance">Buffer 'AVariance'.</param>
		/// <param name="BMean">Buffer 'BMean'.</param>
		/// <param name="BVariance">Buffer 'BVariance'.</param>
		/// <returns>The outgoing VMP message to the 'innerProduct' argument</returns>
		/// <remarks><para>
		/// The outgoing message is the factor viewed as a function of 'innerProduct' conditioned on the given values.
		/// </para></remarks>
		public static Gaussian InnerProductAverageLogarithm(Vector AMean, PositiveDefiniteMatrix AVariance, Vector BMean, PositiveDefiniteMatrix BVariance)
		{
			Gaussian result = new Gaussian();
			// p(x|a,b) = N(E[a]'*E[b], E[b]'*var(a)*E[b] + E[a]'*var(b)*E[a] + trace(var(a)*var(b)))
			// Uses John Winn's rule for deterministic factors.
			// Strict variational inference would set the variance to 0.
			result.SetMeanAndVariance(AMean.Inner(BMean), AVariance.QuadraticForm(BMean) + BVariance.QuadraticForm(AMean) + AVariance.Inner(BVariance));
			return result;
		}

 public static Gaussian UsesAverageLogarithm(NonconjugateGaussian[] Uses, Gaussian Def, Gaussian result)
 {
     NonconjugateGaussian prod = Uses[0];
     for (int i = 1; i < Uses.Length; i++)
         prod.SetToProduct(prod, Uses[i]);
     result = prod.GetGaussian(true);
     result.SetToProduct(result, Def);
     return result;
 }

		/// <summary>
		/// VMP message to 'a'
		/// </summary>
		/// <param name="Product">Incoming message from 'product'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="B">Incoming message from 'b'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="result">Modified to contain the outgoing message</param>
		/// <returns><paramref name="result"/></returns>
		/// <remarks><para>
		/// The outgoing message is the exponential of the average log-factor value, where the average is over all arguments except 'a'.
		/// Because the factor is deterministic, 'product' is integrated out before taking the logarithm.
		/// The formula is <c>exp(sum_(b) p(b) log(sum_product p(product) factor(product,a,b)))</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="Product"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
		public static Gaussian AAverageLogarithm([SkipIfUniform] Gaussian Product, [Proper] Beta B)
		{
			if (B.IsPointMass) return GaussianProductVmpOp.AAverageLogarithm(Product, B.Point);
			if (Product.IsPointMass) return AAverageLogarithm(Product.Point, B);
			double mb, vb;
			B.GetMeanAndVariance(out mb, out vb);
			Gaussian result = new Gaussian();
			result.Precision = Product.Precision * (vb + mb*mb);
			result.MeanTimesPrecision = Product.MeanTimesPrecision * mb;
			return result;
		}

        private static void TestEnvelope()
        {
            Gaussian gaussian = new Gaussian(1, 1);
            ExponentialEnvelope env = new ExponentialEnvelope(double.NegativeInfinity, double.PositiveInfinity, gaussian, new double[] { -2, -1, 0, 1, 2 });

            for (int i = 0; i < 10000; i++)
            {
                double test = env.SampleContinuous();
                Console.WriteLine(test);
            }
            Console.WriteLine("All Worked!");
        }

		/// <summary>
		/// Gradient matching VMP message from factor to logOdds variable
		/// </summary>
		/// <param name="sample">Incoming message from 'sample'.</param>
		/// <param name="logOdds">Incoming message. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="result">Previous message sent, used for damping</param>
		/// <returns>The outgoing VMP message.</returns>
		/// <remarks><para>
		/// The outgoing message is the Gaussian approximation to the factor which results in the 
		/// same derivatives of the KL(q||p) divergence with respect to the parameters of the posterior
		/// as if the true factor had been used.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="logOdds"/> is not a proper distribution</exception>
        public static Gaussian LogRateAverageLogarithm(int sample, [Proper, SkipIfUniform] Gaussian logRate, Gaussian result)
		{
            double m,v;
            logRate.GetMeanAndVariance(out m,out v);
            Gaussian message = new Gaussian();
            message.Precision = Math.Exp(m + v / 2);
            message.MeanTimesPrecision = (m - 1) * Math.Exp(m + v / 2) + sample;
            if (damping == 0)
                return message;
            else
                return (message ^ (1 - damping)) * (result ^ damping); 
		}

		/// <summary>
		/// VMP message to 'a'
		/// </summary>
		/// <param name="Product">Incoming message from 'product'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="B">Incoming message from 'b'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <returns>The outgoing VMP message to the 'a' argument</returns>
		/// <remarks><para>
		/// The outgoing message is the exponential of the average log-factor value, where the average is over all arguments except 'a'.
		/// Because the factor is deterministic, 'product' is integrated out before taking the logarithm.
		/// The formula is <c>exp(sum_(b) p(b) log(sum_product p(product) factor(product,a,b)))</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="Product"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
		public static Gaussian AAverageLogarithm([SkipIfUniform] Gaussian Product, [Proper] Gamma B)
		{
			if (B.IsPointMass) return GaussianProductVmpOp.AAverageLogarithm(Product, B.Point);
			if (Product.IsPointMass) return AAverageLogarithm(Product.Point, B);
			if (!B.IsProper()) throw new ImproperMessageException(B);
			double mb, vb;
			B.GetMeanAndVariance(out mb, out vb);
			// catch uniform case to avoid 0*Inf
			if (Product.IsUniform()) return Product;
			Gaussian result = new Gaussian();
			result.Precision = Product.Precision * (vb + mb*mb);
			result.MeanTimesPrecision = Product.MeanTimesPrecision * mb;
			return result;
		}

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="isPositive">Incoming message from 'isPositive'.</param>
		/// <param name="x">Incoming message from 'x'.</param>
		/// <returns>Logarithm of the factor's average value across the given argument distributions</returns>
		/// <remarks><para>
		/// The formula for the result is <c>log(sum_(isPositive,x) p(isPositive,x) factor(isPositive,x))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(Bernoulli isPositive, Gaussian x)
		{
			if (isPositive.IsPointMass && x.Precision == 0) {
				double tau = x.MeanTimesPrecision;
				if (isPositive.Point && tau < 0) {
					// int I(x>0) exp(tau*x) dx = -1/tau
					return -Math.Log(-tau);
				}
				if (!isPositive.Point && tau > 0) {
					// int I(x<0) exp(tau*x) dx = 1/tau
					return -Math.Log(tau);
				}
			}
			Bernoulli to_isPositive = IsPositiveAverageConditional(x);
			return isPositive.GetLogAverageOf(to_isPositive);
		}

 public void Gaussian_Bernoulli_Conditional()
 {
     // arrange
     Uncertain<double> X = new Gaussian(1.0, 1.0);
     Uncertain<double> Y = new Gaussian(4.0, 2.0);
     // act
     if ((X > Y).Pr()) {
         Assert.Fail("X > Y evaluates true, incorrectly");
     }
     if ((Y < X).Pr()) {
         Assert.Fail("Y < X evaluates true, incorrectly");
     }
     if (!(Y > X).Pr()) {
         Assert.Fail("Y > X evaluates false, incorrectly");
     }
     if (!(X < Y).Pr()) {
         Assert.Fail("X < Y evaluates false, incorrectly");
     }
 }

        public void Gaussian_Bernoulli_Mean()
        {
            // arrange
            Uncertain<double> X = new Gaussian(1.0, 1.0);
            Uncertain<double> Y = new Gaussian(3.0, 2.0);
            var  Z = X > Y;
            var sampler = Sampler.Create(Z);
            int k = 0;
            // act

            foreach (var s in sampler.Take(100))
                if (s.Value) k += 1;
            // assert
            // Y - X is Gaussian(1, sqrt(5)), and has a 32.74% chance of being < 0
            // (i.e. 32.74% chance that X > Y)
            // The normal approximation to the binomial distribution says that the
            // 99.9997% confidence interval with n=100, p=0.3274 is +/- 0.1883.
            // So the confidence interval is roughly [0.13, 0.52].
            Assert.IsTrue(k >= 13 && k < 52); // flaky at N = 100!
        }

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="isPositive">Incoming message from 'isPositive'.</param>
		/// <param name="x">Incoming message from 'x'.</param>
		/// <returns>Logarithm of the factor's average value across the given argument distributions</returns>
		/// <remarks><para>
		/// The formula for the result is <c>log(sum_(isPositive,x) p(isPositive,x) factor(isPositive,x))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(Bernoulli isPositive, Gaussian x)
		{
			Bernoulli to_isPositive = IsPositiveAverageConditional(x);
			return isPositive.GetLogAverageOf(to_isPositive);
#if false
			// Z = p(b=T) p(x > 0) + p(b=F) p(x <= 0)
			//   = p(b=F) + (p(b=T) - p(b=F)) p(x > 0)
			if (x.IsPointMass) {
				return Factor.IsPositive(x.Point) ? isPositive.GetLogProbTrue() : isPositive.GetLogProbFalse();
			} else if(x.IsUniform()) {
				return Bernoulli.LogProbEqual(isPositive.LogOdds,0.0);
			} else {
				// m/sqrt(v) = (m/v)/sqrt(1/v)
				double z = x.MeanTimesPrecision / Math.Sqrt(x.Precision);
				if (isPositive.IsPointMass) {
					return isPositive.Point ? MMath.NormalCdfLn(z) : MMath.NormalCdfLn(-z);
				} else {
					return MMath.LogSumExp(isPositive.GetLogProbTrue() + MMath.NormalCdfLn(z), isPositive.GetLogProbFalse() + MMath.NormalCdfLn(-z));
				}
			}
#endif
		}

		/// <summary>
		/// EP message to 'sample'
		/// </summary>
		/// <param name="mean">Incoming message from 'mean'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="precision">Constant value for 'precision'.</param>
		/// <returns>The outgoing EP message to the 'sample' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'sample' as the random arguments are varied.
		/// The formula is <c>proj[p(sample) sum_(mean) p(mean) factor(sample,mean,precision)]/p(sample)</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="mean"/> is not a proper distribution</exception>
		public static Gaussian SampleAverageConditional([SkipIfUniform] Gaussian mean, double precision)
		{
			if (mean.IsPointMass) return SampleAverageConditional(mean.Point, precision);
			// if (precision < 0) throw new ArgumentException("The constant precision given to the Gaussian factor is negative", "precision");
			if (precision == 0) {
				return Gaussian.Uniform();
			} else if (double.IsPositiveInfinity(precision)) {
				return mean;
			} else {
				if (mean.Precision <= -precision) throw new ImproperMessageException(mean);
				// The formula is int_mean N(x;mean,1/prec) p(mean) = N(x; mm, mv + 1/prec)
				// sample.Precision = inv(mv + inv(prec)) = mprec*prec/(prec + mprec)
				// sample.MeanTimesPrecision = sample.Precision*mm = R*(mprec*mm)
				// R = Prec/(Prec + mean.Prec)
				// This code works for mean.IsUniform() since then mean.Precision = 0, mean.MeanTimesPrecision = 0
				Gaussian result = new Gaussian();
				double R = precision / (precision + mean.Precision);
				result.Precision = R * mean.Precision;
				result.MeanTimesPrecision = R * mean.MeanTimesPrecision;
				return result;
			}
		}

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="log">Incoming message from 'log'.</param>
		/// <param name="d">Incoming message from 'd'.</param>
		/// <param name="to_log">Previous outgoing message to 'log'.</param>
		/// <returns>Logarithm of the factor's average value across the given argument distributions</returns>
		/// <remarks><para>
		/// The formula for the result is <c>log(sum_(d,log) p(d,log) factor(log,d))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(Gaussian log, Gamma d, [Fresh] Gaussian to_log)
		{
			Gamma g = Gamma.FromShapeAndRate(d.Shape + 1, d.Rate);
			return d.Shape/d.Rate*ExpOp.LogAverageFactor(g, log, to_log);
		}

		/// <summary>
		/// Evidence message for VMP
		/// </summary>
		/// <param name="log">Incoming message from 'log'.</param>
		/// <param name="d">Incoming message from 'd'.</param>
		/// <returns>Zero</returns>
		/// <remarks><para>
		/// In Variational Message Passing, the evidence contribution of a deterministic factor is zero.
		/// Adding up these values across all factors and variables gives the log-evidence estimate for VMP.
		/// </para></remarks>
		public static double AverageLogFactor(Gaussian log, Gamma d)
		{
			double m, v;
			log.GetMeanAndVariance(out m, out v);
			double Elogd=d.GetMeanLog();
			double Elogd2;
			if (!d.IsPointMass)
				Elogd2 = MMath.Trigamma(d.Shape) + Elogd * Elogd;
			else
				Elogd2 = Math.Log(d.Point) * Math.Log(d.Point);
			return -Elogd2/(2*v)+m*Elogd/v-m*m/(2*v)-MMath.LnSqrt2PI-.5*Math.Log(v);
		}

		/// <summary>
		/// VMP message to 'log'
		/// </summary>
		/// <param name="log">Incoming message from 'log'.</param>
		/// <param name="d">Constant value for 'd'.</param>
		/// <param name="result">Modified to contain the outgoing message</param>
		/// <returns><paramref name="result"/></returns>
		/// <remarks><para>
		/// The outgoing message is the factor viewed as a function of 'log' conditioned on the given values.
		/// </para></remarks>
		public static Gaussian LogAverageLogarithm(Gaussian log, double d, Gaussian result)
		{
			result.Point = Math.Log(d);
			return result;
		}