/// Generate LDA data - returns an array of dictionaries mapping unique word index
		/// to word count per document.
		/// <param name="trueTheta">Known Theta</param>
		/// <param name="truePhi">Known Phi</param>
		/// <param name="averageNumWords">Average number of words to sample per doc</param>
		/// <returns></returns>
		public static Dictionary<int, int>[] GenerateLDAData(Dirichlet[] trueTheta, Dirichlet[] truePhi, int averageNumWords)
		{
			int numVocab = truePhi[0].Dimension;
			int numTopics = truePhi.Length;
			int numDocs = trueTheta.Length;

			// Sample from the model
			Vector[] topicDist = new Vector[numDocs];
			Vector[] wordDist = new Vector[numTopics];
			for (int i = 0; i < numDocs; i++)
				topicDist[i] = trueTheta[i].Sample();
			for (int i = 0; i < numTopics; i++)
				wordDist[i] = truePhi[i].Sample();

			var wordCounts = new Dictionary<int, int>[numDocs];
			for (int i=0; i < numDocs; i++)
			{
				int LengthOfDoc = Poisson.Sample((double)averageNumWords);

				var counts = new Dictionary<int, int>();
				for (int j=0; j < LengthOfDoc; j++)
				{
					int topic = Discrete.Sample(topicDist[i]);
					int w = Discrete.Sample(wordDist[topic]);
					if (!counts.ContainsKey(w))
						counts.Add(w, 1);
					else
						counts[w] = counts[w] + 1;
				}
				wordCounts[i] = counts;
			}
			return wordCounts;
		}

 public void CanGetAlpha()
 {
     var d = new Dirichlet(0.3, 10);
     for (var i = 0; i < 10; i++)
     {
         Assert.AreEqual(0.3, d.Alpha[i]);
     }
 }

        public void CanCreateSymmetricDirichlet()
        {
            var d = new Dirichlet(0.3, 5);

            for (var i = 0; i < 5; i++)
            {
                Assert.AreEqual(0.3, d.Alpha[i]);
            }
        }

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="sample">Incoming message from 'sample'.</param>
		/// <param name="trialCount">Constant value for 'trialCount'.</param>
		/// <param name="p">Incoming message from 'p'.</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_(sample,p) p(sample,p) factor(sample,trialCount,p))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(IList<int> sample, int trialCount, Dirichlet p)
		{
			double result = MMath.GammaLn(trialCount+1);
			for (int i = 0; i < sample.Count; i++) {
				result += MMath.GammaLn(sample[i]+p.PseudoCount[i]) + MMath.GammaLn(p.PseudoCount[i])
					-MMath.GammaLn(sample[i]+1);
			}
			result += MMath.GammaLn(p.TotalCount) - MMath.GammaLn(p.TotalCount + trialCount);
			return result;
		}

        public void CanGetAlpha()
        {
            Dirichlet d = new Dirichlet(0.3, 10);

            double[] alpha = new double[10];
            for (int i = 0; i < 10; i++)
            {
                Assert.AreEqual(0.3, d.Alpha[i]);
            }
        }

        public void CanCreateDirichlet()
        {
            var alpha = new double[10];
            for (var i = 0; i < 10; i++)
            {
                alpha[i] = i;
            }

            var d = new Dirichlet(alpha);

            for (var i = 0; i < 5; i++)
            {
                Assert.AreEqual(i, d.Alpha[i]);
            }
        }

		/// <summary>
		/// Randomly create true theta and phi arrays
		/// </summary>
		/// <param name="numVocab">Vocabulary size</param>
		/// <param name="numTopics">Number of topics</param>
		/// <param name="numDocs">Number of documents</param>
		/// <param name="averageDocLength">Avarage document length</param>
		/// <param name="trueTheta">Theta array (output)</param>
		/// <param name="truePhi">Phi array (output)</param>
		public static void CreateTrueThetaAndPhi(
			int numVocab, int numTopics, int numDocs, int averageDocLength, int averageWordsPerTopic,
			out Dirichlet[] trueTheta, out Dirichlet[] truePhi)
		{
			truePhi = new Dirichlet[numTopics];
			for (int i=0; i < numTopics; i++)
			{
				truePhi[i] = Dirichlet.Uniform(numVocab);
				truePhi[i].PseudoCount.SetAllElementsTo(0.0);
				// Draw the number of unique words in the topic.
				int numUniqueWordsPerTopic = Poisson.Sample((double)averageWordsPerTopic);
				if (numUniqueWordsPerTopic >= numVocab) numUniqueWordsPerTopic = numVocab;
				if (numUniqueWordsPerTopic < 1) numUniqueWordsPerTopic = 1;
				double expectedRepeatOfWordInTopic = 
					((double)numDocs) * averageDocLength / numUniqueWordsPerTopic;
				int[] shuffledWordIndices = Rand.Perm(numVocab);
				for (int j = 0; j < numUniqueWordsPerTopic; j++)
				{
					int wordIndex = shuffledWordIndices[j];
					// Draw the count for that word
					int cnt = Poisson.Sample(expectedRepeatOfWordInTopic);
					truePhi[i].PseudoCount[wordIndex] = cnt + 1.0;
				}
			}

			trueTheta = new Dirichlet[numDocs];
			for (int i=0; i < numDocs; i++)
			{
				trueTheta[i] = Dirichlet.Uniform(numTopics);
				trueTheta[i].PseudoCount.SetAllElementsTo(0.0);
				// Draw the number of unique topics in the doc.
				int numUniqueTopicsPerDoc = Math.Min(1 + Poisson.Sample(1.0), numTopics);
				double expectedRepeatOfTopicInDoc = 
					averageDocLength / numUniqueTopicsPerDoc;
				int[] shuffledTopicIndices = Rand.Perm(numTopics);
				for (int j = 0; j < numUniqueTopicsPerDoc; j++)
				{
					int topicIndex = shuffledTopicIndices[j];
					// Draw the count for that topic
					int cnt = Poisson.Sample(expectedRepeatOfTopicInDoc);
					trueTheta[i].PseudoCount[topicIndex] = cnt + 1.0;
				}
			}
		}

		/// <summary>
		/// Run a single test for a single model
		/// </summary>
		/// <param name="sizeVocab">Size of the vocabulary</param>
		/// <param name="numTopics">Number of topics</param>
		/// <param name="trainWordsInTrainDoc">Lists of words in training documents used for training</param>
		/// <param name="testWordsInTrainDoc">Lists of words in training documents used for testing</param>
		/// <param name="alpha">Background pseudo-counts for distributions over topics</param>
		/// <param name="beta">Background pseudo-counts for distributions over words</param>
		/// <param name="shared">If true, uses shared variable version of the model</param>
		/// <param name="trueThetaTest">The true topic distributions for the documents in the test set</param>
		/// <param name="wordsInTestDoc">Lists of words in test documents</param>
		/// <param name="vocabulary">Vocabulary</param>
		static void RunTest(
			int sizeVocab, 
			int numTopics,
			Dictionary<int, int>[] trainWordsInTrainDoc, 
			Dictionary<int, int>[] testWordsInTrainDoc,
			double alpha, 
			double beta, 
			bool shared,
			Dirichlet[] trueThetaTest, 
			Dictionary<int, int>[] wordsInTestDoc,
			Dictionary<int, string> vocabulary = null
			)
		{
			Stopwatch stopWatch = new Stopwatch();
			// Square root of number of documents is the optimal for memory
			int batchCount = (int)Math.Sqrt((double)trainWordsInTrainDoc.Length);
			Rand.Restart(5);
			ILDA model;
			LDAPredictionModel predictionModel;
			LDATopicInferenceModel topicInfModel;
			if (shared)
			{
				model = new LDAShared(batchCount, sizeVocab, numTopics);
				((LDAShared)model).IterationsPerPass = Enumerable.Repeat(10, 5).ToArray();
			}
			else
			{
				model = new LDAModel(sizeVocab, numTopics);
				model.Engine.NumberOfIterations = 50;
			}

			Console.WriteLine("\n\n************************************");
			Console.WriteLine(
				String.Format("\nTraining {0}LDA model...\n",
				shared ? "batched " : "non-batched "));

			// Train the model - we will also get rough estimates of execution time and memory
			Dirichlet[] postTheta, postPhi;
			GC.Collect();
			PerformanceCounter memCounter = new PerformanceCounter("Memory", "Available MBytes");
			float preMem = memCounter.NextValue();
			stopWatch.Reset();
			stopWatch.Start();
			double logEvidence = model.Infer(trainWordsInTrainDoc, alpha, beta, out postTheta, out postPhi);
			stopWatch.Stop();
			float postMem = memCounter.NextValue();
			double approxMB = preMem - postMem;
			GC.KeepAlive(model); // Keep the model alive to this point (for the	memory counter)
			Console.WriteLine(String.Format("Approximate memory usage: {0:F2} MB", approxMB));
			Console.WriteLine(String.Format("Approximate execution time (including model compilation): {0} seconds", stopWatch.ElapsedMilliseconds/1000));

			// Calculate average log evidence over total training words
			int totalWords = trainWordsInTrainDoc.Sum(doc => doc.Sum(w => w.Value));
			Console.WriteLine("\nTotal number of training words = {0}", totalWords);
			Console.WriteLine(String.Format("Average log evidence of model: {0:F2}", logEvidence / (double)totalWords));

			if (vocabulary != null)
			{
				int numWordsToPrint = 20;
				// Print out the top n words for each topic
				for (int i = 0; i < postPhi.Length; i++)
				{
					double[] pc = postPhi[i].PseudoCount.ToArray();
					int[] wordIndices = new int[pc.Length];
					for (int j=0; j < wordIndices.Length; j++)
						wordIndices[j] = j;
					Array.Sort(pc, wordIndices);
					Console.WriteLine("Top {0} words in topic {1}:", numWordsToPrint, i);
					int idx = wordIndices.Length;
					for (int j = 0; j < numWordsToPrint; j++)
						Console.Write("\t{0}", vocabulary[wordIndices[--idx]]);
					Console.WriteLine();
				}
			}

			if (testWordsInTrainDoc != null)
			{
				// Test on unseen words in training documents
				Console.WriteLine("\n\nCalculating perplexity on test words in training documents...");
				predictionModel = new LDAPredictionModel(sizeVocab, numTopics);
				predictionModel.Engine.NumberOfIterations = 5;
				var predDist = predictionModel.Predict(postTheta, postPhi);
				var perplexity = Utilities.Perplexity(predDist, testWordsInTrainDoc);
				Console.WriteLine(String.Format("\nPerplexity = {0:F3}", perplexity));
			}

			if (wordsInTestDoc != null)
//.........这里部分代码省略.........

        /// <summary>
        /// Attachs the data to the workers labels with and sets the workers' confusion matrix priors.
        /// </summary>
        /// <param name="taskIndices">The matrix of the task indices (columns) of each worker (rows).</param>
        /// <param name="workerLabels">The matrix of the labels (columns) of each worker (rows).</param>
        /// <param name="confusionMatrixPrior">The workers' confusion matrix priors.</param>
        protected virtual void AttachData(int[][] taskIndices, int[][] workerLabels, Dirichlet[][] confusionMatrixPrior)
        {
            int numClasses = c.SizeAsInt;
            WorkerCount.ObservedValue = taskIndices.Length;
            WorkerTaskCount.ObservedValue = taskIndices.Select(tasks => tasks.Length).ToArray();
            WorkerTaskIndex.ObservedValue = taskIndices;
            // Prediction mode is indicated by none of the workers having a label.
            // We can just look at the first one
            if (workerLabels[0] != null)
            {
                WorkerLabel.ObservedValue = workerLabels;
            }
            else
            {
                WorkerLabel.ClearObservedValue();
            }

            if (confusionMatrixPrior != null)
            {
                ConfusionMatrixPrior.ObservedValue = Util.ArrayInit(confusionMatrixPrior.Length, worker => Util.ArrayInit(numClasses, lab => confusionMatrixPrior[worker][lab]));
            }
        }

		/// <summary>
		/// VMP message to 'mean'
		/// </summary>
		/// <param name="mean">Incoming message from 'mean'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <param name="totalCount">Constant value for 'totalCount'.</param>
		/// <param name="prob">Incoming message from 'prob'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <param name="to_mean">Previous outgoing message to 'mean'.</param>
		/// <returns>The outgoing VMP message to the 'mean' argument</returns>
		/// <remarks><para>
		/// The outgoing message is the exponential of the average log-factor value, where the average is over all arguments except 'mean'.
		/// The formula is <c>exp(sum_(prob) p(prob) log(factor(prob,mean,totalCount)))</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="mean"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="prob"/> is not a proper distribution</exception>
		public static Dirichlet MeanAverageLogarithm([Proper] Dirichlet mean, double totalCount, [SkipIfUniform] Dirichlet prob, Dirichlet to_mean)
		{
			return MeanAverageLogarithm(mean, Gamma.PointMass(totalCount), prob, to_mean);
		}

		public static Dirichlet ProbAverageConditional([SkipIfUniform] Gamma alpha, Dirichlet result)
		{
			throw new NotSupportedException(NotSupportedMessage);
		}

		public static double LogEvidenceRatio(Dirichlet prob, double alpha) { return 0.0; }

		/// <summary>
		/// Evidence message for VMP
		/// </summary>
		/// <param name="sample">Incoming message from 'sampleFromPseudoCounts'.</param>
		/// <param name="pseudoCounts">Constant value for 'pseudoCount'.</param>
		/// <param name="to_sample">Outgoing message to 'sample'.</param>
		/// <returns>Average of the factor's log-value across the given argument distributions</returns>
		/// <remarks><para>
		/// The formula for the result is <c>sum_(sampleFromPseudoCounts) p(sampleFromPseudoCounts) log(factor(sampleFromPseudoCounts,pseudoCount))</c>.
		/// Adding up these values across all factors and variables gives the log-evidence estimate for VMP.
		/// </para></remarks>
		public static double AverageLogFactor(Dirichlet sample, Vector pseudoCounts, [Fresh] Dirichlet to_sample)
		{
			return to_sample.GetAverageLog(sample);
		}

        public void ValidateVariance()
        {
            var alpha = new double[10];
            var sum = 0.0;
            for (var i = 0; i < 10; i++)
            {
                alpha[i] = i;
                sum += i;
            }

            var d = new Dirichlet(alpha);
            for (var i = 0; i < 10; i++)
            {
                AssertHelpers.AlmostEqual(i * (sum - i) / (sum * sum * (sum + 1.0)), d.Variance[i], 15);
            }
        }

        public void ValidateMean()
        {
            var d = new Dirichlet(0.3, 5);

            for (var i = 0; i < 5; i++)
            {
                AssertHelpers.AlmostEqual(0.3 / 1.5, d.Mean[i], 15);
            }
        }

		public static Dirichlet ProbAverageLogarithm([SkipIfUniform] Gamma alpha, Dirichlet result)
		{
			result.PseudoCount.SetAllElementsTo(alpha.GetMean());
            result.TotalCount = result.PseudoCount.Sum(); 
			return result;
		}

 public void CanGetDimension()
 {
     var d = new Dirichlet(0.3, 10);
     Assert.AreEqual(10, d.Dimension);
 }

		/// <summary>
		/// Evidence message for EP
		/// </summary>
		/// <param name="prob">Incoming message from 'prob'.</param>
		/// <param name="alpha">Constant value for 'alpha'.</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 int_prob p(prob) factor(prob,alpha,K))</c>.
		/// </para></remarks>
		public static double LogAverageFactor(Dirichlet prob, double alpha)
		{
			int dim = prob.Dimension;
			double sum = MMath.GammaLn(dim*alpha) - dim*MMath.GammaLn(alpha) - prob.GetLogNormalizer() - MMath.GammaLn(prob.TotalCount+dim*alpha);
			for (int i = 0; i < dim; i++) {
				sum += MMath.GammaLn(prob.PseudoCount[i]+alpha);
			}
			return sum;
		}

		/// <summary>
		/// Perform the quadrature required for the VMP evidence message
		/// </summary>
		/// <param name="meanQ">Incoming message from m='mean'.</param>
		/// <param name="totalCountQ">Incoming message from s='totalCount'.</param>
		/// <returns>Vector of E[ LogGamma(s*m_k)].</returns>
		/// <remarks><para>
		/// The quadrature over 'totalCount' (which is Gamma-distributed) is 
		/// peformed by a change of variable x=log(s) followed by Gauss-Hermite 
		/// quadrature. The quadrature over m is performed using Gauss-Legendre. 
		/// </para></remarks>
		public static Vector EvidenceMessageExpectations(
				Dirichlet meanQ,
				Gamma totalCountQ)
		{
			// Get shape and scale of the distribution
			double at, bt;
			totalCountQ.GetShapeAndScale(out at, out bt);
			bt = 1 / bt; // want rate not scale

			// Mean in the transformed domain
			double proposalMean = totalCountQ.GetMeanLog();
			// Laplace approximation of variance in transformed domain 
			double proposalVariance = 1 / at;

			// Quadrature coefficient
			int nt = 32;
			Vector nodes = Vector.Zero(nt);
			Vector weights = Vector.Zero(nt);
			Vector expx = Vector.Zero(nt);
			if (!totalCountQ.IsPointMass) {
				Quadrature.GaussianNodesAndWeights(proposalMean, proposalVariance, nodes, weights);
				// Precompute weights for each m slice
				for (int i = 0; i < nt; i++) {
					double x = nodes[i];
					expx[i] = Math.Exp(x);
					double p = at * x - bt * expx[i] - Gaussian.GetLogProb(x, proposalMean, proposalVariance);
					weights[i] *= Math.Exp(p);
				}
			}

			int nm = 20;
			Vector mnodes = Vector.Zero(nm);
			Vector mweight = Vector.Zero(nm);
			Quadrature.UniformNodesAndWeights(0, 1, mnodes, mweight);
			int K = meanQ.Dimension;
			Vector[] mweights = new Vector[K];
			Beta[] mkDist = new Beta[K];
			double[] EELogGamma = new double[K];
			for (int i = 0; i < K; i++) {
				mweights[i] = Vector.Copy(mweight);
				mkDist[i] = new Beta(meanQ.PseudoCount[i], meanQ.TotalCount - meanQ.PseudoCount[i]);
				EELogGamma[i] = 0;
			}

			for (int j = 0; j < nm; j++) {
				double m = mnodes[j];
				double ELogGamma = 0;
				if (totalCountQ.IsPointMass)
					ELogGamma = MMath.GammaLn(m * totalCountQ.Point);
				else {
					// Calculate expectations in x=log(s) space using Gauss-Hermite quadrature
					for (int i = 0; i < nt; i++) {
						double x = nodes[i];
						ELogGamma += weights[i] * (MMath.GammaLn(m * expx[i]) + x);
					}
					// Normalise and add removed components
					double normalisation = Math.Pow(bt, at) / MMath.Gamma(at);
					ELogGamma = normalisation * ELogGamma - proposalMean;
				}
				for (int i = 0; i < K; i++) {
					mweights[i][j] *= Math.Exp(mkDist[i].GetLogProb(m));
					EELogGamma[i] += mweights[i][j] * ELogGamma;
				}
			}
			return Vector.FromArray(EELogGamma);
		}

        /// <summary>
        /// Returns the confusion matrix prior of each worker.
        /// </summary>
        /// <returns>The confusion matrix prior of each worker.</returns>
        public Dirichlet[] GetConfusionMatrixPrior()
        {
            var confusionMatrixPrior = new Dirichlet[LabelCount];
            for (int d = 0; d < LabelCount; d++)
            {
                //confusionMatrixPrior[d] = new Dirichlet(Util.ArrayInit(LabelCount, i => i == d ? (InitialWorkerBelief / (1 - InitialWorkerBelief)) * (LabelCount - 1) : 1.0));
                confusionMatrixPrior[d] = new Dirichlet(Util.ArrayInit(LabelCount, i => i == d ? 0.8 : 0.2));
            }

            return confusionMatrixPrior;
        }