private static void Test_BPMIncremental(
			int nClass, int totalFeatures, double noisePrec, int maxItemsInBatch,
			int nChunks, string trainingFile, Vector[] testData)
		{
			Console.WriteLine("\n------- BPM Train Incremental -------");
			VectorGaussian[] wInfer = new VectorGaussian[nClass];
			BPM bpmIncremental = new BPM(nClass, totalFeatures, noisePrec);
			bpmIncremental.TrainingEngine.ShowProgress = false;
			bpmIncremental.TestEngine.ShowProgress = false;
			int LocToStart = 0;
			for (int c = 0; c < nChunks; c++)
			{
				List<Vector>[] dataChunks = DataFromFile.Read(trainingFile, nClass, maxItemsInBatch, ref LocToStart);
				wInfer = bpmIncremental.TrainIncremental(dataChunks);
			}
#if ShowWeights
			for (int i = 0; i < wInfer.GetLength(0); i++)
			{
				Console.WriteLine(wInfer[i].ToString());
			}
#endif
			Console.WriteLine("\nPredictions:");
			Discrete[] predictions = bpmIncremental.Test(testData);
			foreach (Discrete pred in predictions)
				Console.WriteLine(pred);
			Console.WriteLine();
		}

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="concat">Constant value for 'concat'.</param>
		/// <param name="first">Constant value for 'first'.</param>
		/// <param name="second">Incoming message from 'second'.</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_(second) p(second) factor(concat,first,second))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(Vector concat, Vector first, VectorGaussian second)
		{
			for (int i = 0; i < first.Count; i++) {
				if (concat[i] != first[i]) return Double.NegativeInfinity;
			}
			Vector concat2 = Vector.Subvector(concat, first.Count, second.Dimension);
			return second.GetLogProb(concat2);
		}

		/// <summary>
		/// EP message to 'B'
		/// </summary>
		/// <param name="sum">Incoming message from 'Sum'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="A">Constant value for 'A'.</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 'B' as the random arguments are varied.
		/// The formula is <c>proj[p(B) sum_(Sum) p(Sum) factor(Sum,A,B)]/p(B)</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="sum"/> is not a proper distribution</exception>
		public static VectorGaussian BAverageConditional([SkipIfUniform] Gaussian sum, bool[] A, VectorGaussian result)
		{
			if (result == default(VectorGaussian)) result = new VectorGaussian(A.Length);
			// (m - a'b)^2/v = (a'bb'a - 2a'bm + m^2)/v
			var ma = Vector.FromArray(A.Select(x => x ? 1.0 : 0.0).ToArray());
			result.Precision.SetToOuter(ma, ma);
			result.Precision.Scale(sum.Precision);
			result.MeanTimesPrecision.SetToProduct(ma, sum.MeanTimesPrecision);
			return result;
		}

        /// <summary>
        /// Test the test data.
        /// </summary>
        /// <param name="priorValue">Prior distribution from training.</param>
        /// <param name="testData">The test data.</param>
        /// <returns>The prediction.</returns>
        public Discrete[] Test(VectorGaussian[] priorValue, Vector[] testData)
        {
            // Set the prior of all classes
            for (int i = 0; i < numOfClasses; i++)
            {
                this.classes[i].Prior.ObservedValue = priorValue[i];
            }

            // Set the observed test data
            this.numOfVectors.ObservedValue = testData.Length;
            this.featureVectors.ObservedValue = testData;

            // Infer the test model outputs
            return Distribution.ToArray<Discrete[]>(this.Engine.Infer(this.modelOutput));
        }

        public static VectorGaussianWishart Combine(VectorGaussian position, Wishart orientation, VectorGaussianWishart result)
        {
            if (orientation.IsUniform())
            {
                result.SetToUniform();
            }
            else if (position.IsUniform())
            {
                result.SetTo(orientation.Shape, orientation.Rate, Vector.Zero(2), 0);
            }
            else
            {
                PositiveDefiniteMatrix rateTimesPrecision = new PositiveDefiniteMatrix(2, 2);
                rateTimesPrecision.SetToProduct(orientation.Rate, position.Precision);
                double trace = MathHelpers.Invert(rateTimesPrecision).Trace();
                Vector positionMean = position.MeanTimesPrecision * MathHelpers.Invert(position.Precision);
                result.SetTo(orientation.Shape, orientation.Rate, positionMean, orientation.Dimension / (orientation.Shape * trace));
            }

            return result;
        }

		public static VectorGaussian RotateAverageLogarithm([SkipIfUniform] Gaussian x, [SkipIfUniform] Gaussian y, [Proper] WrappedGaussian angle, VectorGaussian result)
		{
			// for x ~ N(m,v):
			// E[cos(x)] = cos(m)*exp(-v/2)
			// E[sin(x)] = sin(m)*exp(-v/2)
			if (angle.Period != 2*Math.PI) throw new ArgumentException("angle.Period ("+angle.Period+") != 2*PI ("+2*Math.PI+")");
			double angleMean, angleVar;
			angle.Gaussian.GetMeanAndVariance(out angleMean, out angleVar);
			double expHalfVar = Math.Exp(-0.5*angleVar);
			double mCos = Math.Cos(angleMean)*expHalfVar;
			double mSin = Math.Sin(angleMean)*expHalfVar;
			double mCos2 = mCos*mCos;
			double mSin2 = mSin*mSin;
			//  E[cos(x)^2] = 0.5 E[1+cos(2x)] = 0.5 (1 + cos(2m) exp(-2v))
			//  E[sin(x)^2] = E[1 - cos(x)^2] = 0.5 (1 - cos(2m) exp(-2v))
			double expVar = expHalfVar*expHalfVar;
			// cos2m = cos(2m)*exp(-v)
			double cos2m = 2*mCos2 - expVar;
			double mCosSqr = 0.5*(1 + cos2m*expVar);
			double mSinSqr = 1 - mCosSqr;
			double mSinCos = mSin*mCos*expVar;
			if (result.Dimension != 2) throw new ArgumentException("result.Dimension ("+result.Dimension+") != 2");
			double mx, vx, my, vy;
			x.GetMeanAndVariance(out mx, out vx);
			y.GetMeanAndVariance(out my, out vy);
			Vector mean = Vector.Zero(2);
			mean[0] = mCos*mx - mSin*my;
			mean[1] = mSin*mx + mCos*my;
			double mx2 = mx*mx + vx;
			double my2 = my*my + vy;
			double mxy = mx*my;
			PositiveDefiniteMatrix variance = new PositiveDefiniteMatrix(2, 2);
			variance[0, 0] = mx2*mCosSqr - 2*mxy*mSinCos + my2*mSinSqr - mean[0]*mean[0];
			variance[0, 1] = (mx2 - my2)*mSinCos + mxy*(mCosSqr - mSinSqr) - mean[0]*mean[1];
			variance[1, 0] = variance[0, 1];
			variance[1, 1] = mx2*mSinSqr + 2*mxy*mSinCos + my2*mCosSqr - mean[1]*mean[1];
			result.SetMeanAndVariance(mean, variance);
			return result;
		}

        private static IEnumerable<Vector> CreateTestData(int dimensions, int clusters)
        {
            var dimRange = new Range(dimensions);
            var clusterRange = new Range(clusters);
            var trueMeans = Variable.Array<Vector>(clusterRange);
            var trueCovariance = Variable.Array<PositiveDefiniteMatrix>(clusterRange);
            var trueVectorGaussian = new VectorGaussian[clusters];
            for (var i = 0; i < clusters; i++)
            {
                var trueMeansVector = Vector.FromArray(CreateRandomArray(dimensions));
                Console.WriteLine("True means {0}", trueMeansVector);
                trueVectorGaussian[i] = VectorGaussian.FromMeanAndPrecision(trueMeansVector, PositiveDefiniteMatrix.Identity(dimensions));
            }

            var vectorGaussians = Variable.Array<VectorGaussian>(clusterRange);
            var dirichletMixture = Dirichlet.Uniform(clusters).Sample();
            Console.WriteLine("True mixture {0}", dirichletMixture);

            while (true)
            {
                var index = Rand.Sample(dirichletMixture);
                yield return trueVectorGaussian[index].Sample();
            }
        }

		/// <summary>
		/// Evidence message for VMP
		/// </summary>
		/// <param name="product">Constant value for 'product'.</param>
		/// <param name="A">Constant value for 'a'.</param>
		/// <param name="B">Incoming message from 'b'.</param>
		/// <param name="BMean">Buffer 'BMean'.</param>
		/// <param name="BVariance">Buffer 'BVariance'.</param>
		/// <returns>Zero</returns>
		/// <remarks><para>
		/// The formula for the result is <c>log(factor(product,a,b))</c>.
		/// Adding up these values across all factors and variables gives the log-evidence estimate for VMP.
		/// </para></remarks>
		public static double AverageLogFactor(Vector product, Matrix A, VectorGaussian B, [Fresh] Vector BMean, [Fresh] PositiveDefiniteMatrix BVariance) { return LogAverageFactor(product, A, B, BMean, BVariance); }

		public static double LogEvidenceRatio(VectorGaussian product, Matrix A, VectorGaussian B) { return 0.0; }

		/// <summary>
		/// VMP message to 'B'
		/// </summary>
		/// <param name="sum">Incoming message from 'Sum'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="A">Incoming message from 'A'.</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 'B'.
		/// Because the factor is deterministic, 'Sum' is integrated out before taking the logarithm.
		/// The formula is <c>exp(sum_(A) p(A) log(sum_Sum p(Sum) factor(Sum,A,B)))</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="sum"/> is not a proper distribution</exception>
		public static VectorGaussian BAverageLogarithm([SkipIfUniform] Gaussian sum, DistributionStructArray<Bernoulli, bool> A, VectorGaussian result)
		{
			if (result == default(VectorGaussian)) result = new VectorGaussian(A.Count);
			// E[log N(x; ab, 0)] = -0.5 E[(x-ab)^2]/0 = -0.5 (E[x^2] - 2 E[x] a' E[b] + trace(aa' E[bb']))/0
			// message to a = N(a; E[x]*inv(var(b)+E[b]E[b]')*E[b], var(x)*inv(var(b)+E[b]E[b]'))
			// result.Precision = (var(b)+E[b]*E[b]')/var(x)
			// result.MeanTimesPrecision = E[x]/var(x)*E[b] = E[b]*X.MeanTimesPrecision
			Vector ma = Vector.FromArray(A.Select(x => x.GetMean()).ToArray());
			Vector va = Vector.FromArray(A.Select(x => x.GetVariance()).ToArray());
			result.Precision.SetToDiagonal(va);
			result.Precision.SetToSumWithOuter(result.Precision, 1, ma, ma);
			result.Precision.SetToProduct(result.Precision, sum.Precision);
			result.MeanTimesPrecision.SetToProduct(ma, sum.MeanTimesPrecision);
			return result;
		}

		public static VectorGaussian BAverageLogarithm(Vector product, Matrix A, VectorGaussian result)
		{
			throw new NotSupportedException(MatrixVectorProductOp.LowRankNotSupportedMessage);
		}

		/// <summary>
		/// EP message to 'b'
		/// </summary>
		/// <param name="product">Incoming message from 'product'.</param>
		/// <param name="A">Constant value for 'a'.</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 'b' as the random arguments are varied.
		/// The formula is <c>proj[p(b) sum_(product) p(product) factor(product,a,b)]/p(b)</c>.
		/// </para></remarks>
		public static VectorGaussian BAverageConditional(VectorGaussian product, Matrix A, VectorGaussian result)
		{
			if (product.IsPointMass) return BAverageConditional(product.Point, A, result);
			//   (p.mean - A*B)'*p.prec*(p.mean - A*B)
			// = B'*(A'*p.prec*A)*B - B'*A'*p.prec*p.mean - ...
			// B.prec = A'*p.prec*A
			// B.precTimesMean = A'*p.precTimesMean
			Matrix temp = (product.Precision*A).Transpose();
			result.Precision.SetToProduct(temp, A);
			result.MeanTimesPrecision.SetToProduct(product.MeanTimesPrecision, A);
			return result;
		}

		/// <summary>
		/// Performs inference on this Bayes point machine
		/// </summary>
		/// <param name="xValuesData">Lists of vectors for each component</param>
		/// <returns></returns>
		private VectorGaussian[] InferW(List<Vector>[] xValuesData)
		{
			// Set the observed data
			for (int c = 0; c < nClass; c++)
			{
				trainModel.nItems[c].ObservedValue = xValuesData[c].Count;
				trainModel.xValues[c].ObservedValue = xValuesData[c].ToArray();
			}
			// Infer the weights
			VectorGaussian[] wInfer = new VectorGaussian[nClass];
			// Reset the priors to support incremental training
			for (int c = 0; c < nClass; c++)
			{
				wInfer[c] = (trainModel.ie).Infer<VectorGaussian>(trainModel.w[c]);
				trainModel.wInit[c].ObservedValue = wInfer[c];
			}
			return wInfer;
		}

        /// <summary>
        /// Train the machine with training data.
        /// </summary>
        /// <param name="trainingData">Training data</param>
        private void TrainModel(IList<DataVector> trainingData)
        {
            // Create the array to store classification results
            VariableArray<bool> trainResults =
                Variable.Observed(trainingData.Select(r => r.ClassId == 1).ToArray()).Named("trainResults");

            InitWeight();

            // Create an array to represent the observed training data
            Range resultRange = trainResults.Range.Named("relevance");
            VariableArray<Vector> observedData =
                Variable.Observed(trainingData.Select(r => r.FeatureVector).ToArray(), resultRange).Named("observedData");

            // Create Bayes Point Machine
            trainResults[resultRange] =
                Variable.GaussianFromMeanAndVariance(
                    Variable.InnerProduct(weight, observedData[resultRange]).Named("innerProduct"), this.Noise) > 0;

            // Using expectation propagation to infer the posterior over test data
            this.Posterior = Engine.Infer<VectorGaussian>(this.weight);
        }

		public static double AverageLogFactor(VectorGaussian rotate) { return 0.0; }

		/// <summary>
		/// EP message to 'product'
		/// </summary>
		/// <param name="A">Constant value for 'a'.</param>
		/// <param name="BMean">Buffer 'BMean'.</param>
		/// <param name="BVariance">Buffer 'BVariance'.</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 'product' conditioned on the given values.
		/// </para></remarks>
		public static VectorGaussian ProductAverageConditional(Matrix A, [Fresh] Vector BMean, [Fresh] PositiveDefiniteMatrix BVariance, VectorGaussian result)
		{
			// P.mean = A*B.mean
			// P.var = A*B.var*A'
			// if A is invertible, then
			// P.prec = inv(A)'*inv(B.var)*inv(A)
			// P.precTimesMean = inv(A)'*B.precTimesMean
			Vector rmean = A * BMean;
			PositiveDefiniteMatrix rvariance = new PositiveDefiniteMatrix(result.Dimension, result.Dimension);
			Matrix temp = (A*BVariance).Transpose();
			rvariance.SetToProduct(A, temp);
			result.SetMeanAndVariance(rmean, rvariance);
			return result;
		}

		/// <summary>
		/// VMP message to 'product'
		/// </summary>
		/// <param name="A">Constant value for 'a'.</param>
		/// <param name="BMean">Buffer 'BMean'.</param>
		/// <param name="BVariance">Buffer 'BVariance'.</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 'product' conditioned on the given values.
		/// </para></remarks>
		public static VectorGaussian ProductAverageLogarithm(Matrix A, [Fresh] Vector BMean, [Fresh] PositiveDefiniteMatrix BVariance, VectorGaussian result)
		{
			return ProductAverageConditional(A, BMean, BVariance, result);
		}

		/// <summary>
		/// Runs the BPM using dense features and incremental training.
		/// </summary>
		/// <param name="numClasses">The number of classes</param>
		/// <param name="noisePrecision">The noise precision</param>
		/// <param name="trainingSetFile">The file containing the training set</param>
		/// <param name="testSet">The test set</param>
		/// <param name="maxItemsPerChunk">The maximum number of items per chunk</param>
		/// <param name="numChunks">The total number of chunks</param>
		private static void RunIncrementalBPM(int numClasses, double noisePrecision, string trainingSetFile, bool labelAtEnd, bool addBias, Vector[] testSet, int maxItemsPerChunk, int numChunks)
		{
			Console.WriteLine("\n------- Incremental BPM -------");
			BPM bpm = new BPM(numClasses, noisePrecision);

			int[] labels;
			Vector[] featureVectors;
			VectorGaussian[] posteriorWeights = new VectorGaussian[numClasses];

			int locationToStart = 0;
			for (int c = 0; c < numChunks; c++)
			{
				featureVectors = DataUtils.Read(trainingSetFile, maxItemsPerChunk, labelAtEnd, addBias, out labels, ref locationToStart);
				posteriorWeights = bpm.Train(featureVectors, labels);
			}

#if ShowWeights
			Console.WriteLine("Weights=" + StringUtil.ArrayToString(posteriorWeights));
#endif

			Console.WriteLine("\nPredictions:");
			Discrete[] predictions = bpm.Test(testSet);
			foreach (Discrete pred in predictions)
				Console.WriteLine(pred);
			Console.WriteLine();
		}

		/// <summary>
		/// EP message to 'A'
		/// </summary>
		/// <param name="matrixMultiply">Incoming message from 'matrixMultiply'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <param name="A">Incoming message from 'A'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <param name="B">Constant value for 'B'.</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 'A' as the random arguments are varied.
		/// The formula is <c>proj[p(A) sum_(matrixMultiply) p(matrixMultiply) factor(matrixMultiply,A,B)]/p(A)</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="matrixMultiply"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="A"/> is not a proper distribution</exception>
		public static GaussianArray2D AAverageConditional([SkipIfUniform] GaussianArray2D matrixMultiply, [SkipIfUniform] GaussianArray2D A, double[,] B, GaussianArray2D result)
		{
			int rows = matrixMultiply.GetLength(0);
			int cols = matrixMultiply.GetLength(1);
			int inner = B.GetLength(0);
			if (result == null) result = new GaussianArray2D(rows, inner);
			// sum_{i,j} (m[i,j] - a[i,:]*b[:,j])^2/v[i,j] = 
			// sum_{i,j} (m[i,j]^2 - 2m[i,j]a[i,:]*b[:,j] + a[i,:]*(b[:,j] b[:,j]')*a[i,:]')/v[i,j]
			// meanTimesPrec(a[i,:]) = sum_j (m[i,j]/v[i,j]) b[:,j]
			// prec(a[i,:]) = sum_j b[:,j]*b[:,j]'/v[i,j]
			Vector bj = Vector.Zero(inner);
			Vector mean = Vector.Zero(inner);
			VectorGaussian ai = new VectorGaussian(inner);
			PositiveDefiniteMatrix variance = new PositiveDefiniteMatrix(inner, inner);
			for (int i = 0; i < rows; i++) {
				ai.Precision.SetAllElementsTo(0.0);
				ai.MeanTimesPrecision.SetAllElementsTo(0.0);
				// we are projecting from family of full covariance Gaussians to diagonal
				// covariance, so we should include the context
				for (int c = 0; c < inner; c++) {
					ai.Precision[c, c] = A[i, c].Precision;
					ai.MeanTimesPrecision[c] = A[i, c].MeanTimesPrecision;
				}
				for (int j = 0; j < cols; j++) {
					Gaussian xij = matrixMultiply[i, j];
					for (int k = 0; k < inner; k++) {
						bj[k] = B[k, j];
					}
					if (xij.IsPointMass) throw new NotImplementedException(LowRankNotSupportedMessage);
					ai.Precision.SetToSumWithOuter(ai.Precision, xij.Precision, bj, bj);
					ai.MeanTimesPrecision.SetToSum(1.0, ai.MeanTimesPrecision, xij.MeanTimesPrecision, bj);
				}
				ai.GetMeanAndVariance(mean, variance);
				for (int k = 0; k < inner; k++) {
					Gaussian rik = result[i, k];
					rik.SetMeanAndVariance(mean[k], variance[k, k]);
					result[i, k] = rik / A[i, k];
				}
			}
			return result;
		}

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="product">Constant value for 'product'.</param>
		/// <param name="A">Constant value for 'a'.</param>
		/// <param name="BMean">Buffer 'BMean'.</param>
		/// <param name="BVariance">Buffer 'BVariance'.</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(factor(product,a,b))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(Vector product, Matrix A, VectorGaussian B, [Fresh] Vector BMean, [Fresh] PositiveDefiniteMatrix BVariance)
		{
			VectorGaussian toProduct = ProductAverageConditional(A, BMean, BVariance, new VectorGaussian(A.Rows));
			return toProduct.GetLogProb(product);
		}