Gaussian Process Mixtures

 - Added `ProbGPMixtureMachine` to create GP mixtures from a single GP
model

Author: mandar2812

dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/models/gp/AbstractGPRegressionModel.scala

@@ -74,8 +74,6 @@ abstract class AbstractGPRegressionModel[T, I: ClassTag](
 
   override protected val g: T = data
 
-  def _data = g
-
   val npoints = num
 
   protected var blockSize = 1000

dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/optimization/AbstractGridSearch.scala

@@ -18,9 +18,7 @@ under the License.
 * */
 package io.github.mandar2812.dynaml.optimization
 
-import breeze.linalg.DenseVector
 import org.apache.log4j.Logger
-import io.github.mandar2812.dynaml.utils
 
 /**
   * @author mandar2812 datum 01/12/15.

dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/optimization/CoupledSimulatedAnnealing.scala

@@ -18,13 +18,6 @@ under the License.
 * */
 package io.github.mandar2812.dynaml.optimization
 
-import breeze.linalg.DenseVector
-import breeze.stats.distributions.CauchyDistribution
-import io.github.mandar2812.dynaml.probability.RandomVariable
-import io.github.mandar2812.dynaml.utils
-
-import scala.util.Random
-
 /**
   * Implementation of the Coupled Simulated Annealing algorithm
   * for global optimization.

dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/optimization/GridSearch.scala

@@ -18,7 +18,6 @@ under the License.
 * */
 package io.github.mandar2812.dynaml.optimization
 
-import org.apache.log4j.Logger
 
 /**
  * @author mandar2812 datum 24/6/15.

dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/optimization/ProbGPCommMachine.scala

@@ -222,10 +222,10 @@ object ProbGPCommMachine {
 
 }
 
-/*
+
 class ProbGPMixtureMachine[T, I: ClassTag](
   model: AbstractGPRegressionModel[T, I]) extends
-  ModelTuner[AbstractGPRegressionModel[T, I], GaussianProcessMixture[I]] {
+  AbstractCSA[AbstractGPRegressionModel[T, I], GaussianProcessMixture[I]](model) {
 
   private var policy: String = "CSA"
 
@@ -256,14 +256,67 @@ class ProbGPMixtureMachine[T, I: ClassTag](
     this
   }
 
-  override val system = model
+  private def calculateEnergyLandscape(initialConfig: Map[String, Double], options: Map[String, String]) =
+    if(policy == "CSA") performCSA(initialConfig, options)
+    else getEnergyLandscape(initialConfig, options, meanFieldPrior)
+
+  private def modelProbabilities = DataPipe(ProbGPCommMachine.calculateModelWeightsSigmoid(baselinePolicy) _)
 
   override def optimize(
     initialConfig: Map[String, Double],
     options: Map[String, String]) = {
 
+    //Find out the blocked hyper parameters and their values
+    val blockedHypParams = system.covariance.blocked_hyper_parameters ++ system.noiseModel.blocked_hyper_parameters
+
+    val (kernelPipe, noisePipe) = (system.covariance.asPipe, system.noiseModel.asPipe)
+
+    val (kernelParams, noiseParams) = (
+      system.covariance.hyper_parameters,
+      system.noiseModel.hyper_parameters)
+
+    val blockedState = system._current_state.filterKeys(blockedHypParams.contains)
+
+    val energyLandscape = calculateEnergyLandscape(initialConfig, options)
+
+    val data = system.data
+
+    //Calculate the weights of each configuration
+    val (weights, models) = modelProbabilities(energyLandscape).map(c => {
+
+      val model_state = c._2 ++ blockedState
+
+      implicit val transform = DataPipe(system.dataAsSeq _)
+
+      val model = AbstractGPRegressionModel(
+        kernelPipe(model_state), noisePipe(model_state),
+        system.mean)(
+        data, system.npoints)
+
+
+      (c._1, model)
+    }).unzip
+
+
+    val configsAndWeights = modelProbabilities(energyLandscape).map(c => (c._1, c._2 ++ blockedState))
+
+    logger.info("===============================================")
+    logger.info("Constructing Gaussian Process Mixture")
+
+    logger.info("Number of model instances = "+weights.length)
+    logger.info("--------------------------------------")
+    logger.info(
+      "Calculated model probabilities/weights are \n"+
+        configsAndWeights.map(wc =>
+          "\nConfiguration: \n"+
+            GlobalOptimizer.prettyPrint(wc._2)+
+            "\nProbability = "+wc._1+"\n"
+        ).reduceLeft((a, b) => a++b)
+    )
+    logger.info("--------------------------------------")
+
 
 
+    (new GaussianProcessMixture[I](models, DenseVector(weights.toArray)), Map())
   }
 }
-*/

scripts/stochasticPriors.sc

@@ -5,6 +5,7 @@ import io.github.mandar2812.dynaml.models.bayes.{LinearTrendESGPrior, LinearTren
 import io.github.mandar2812.dynaml.probability._
 import com.quantifind.charts.Highcharts._
 import io.github.mandar2812.dynaml.analysis.implicits._
+import io.github.mandar2812.dynaml.optimization.ProbGPMixtureMachine
 import io.github.mandar2812.dynaml.pipes.Encoder
 
 val rbfc = new RBFCovFunc(1.5)
@@ -55,8 +56,16 @@ sgp_prior.globalOptConfig_(Map("gridStep" -> "0.15", "gridSize" -> "40"))
 val gpModel = gp_prior.posteriorModel(dataset)
 val sgpModel = sgp_prior.posteriorModel(dataset)
 
+
+gp_prior.globalOptConfig_(Map("gridStep" -> "0.0", "gridSize" -> "1", "globalOpt" -> "GS", "policy" -> "GS"))
+val gpModel1 = gp_prior.posteriorModel(dataset)
+val mixt_machine = new ProbGPMixtureMachine(gpModel1)
+val (mix_model, _) = mixt_machine.optimize(gp_prior.covariance.state ++ gp_prior.noiseCovariance.state)
+
+
 val zs: MultGaussianPRV = gpModel.predictiveDistribution(xs)
 val sgp_zs: BlockedMESNRV = sgpModel.predictiveDistribution(xs)
+val mix_zs = mix_model.predictiveDistribution(xs)
 
 val MultGaussianPRV(m, c) = zs
 val eigD = eig(c.toBreezeMatrix)
@@ -75,14 +84,23 @@ if(eValuesPositive) {
   println("Predictive Covariance Ill-Posed!")
 }
 
-val samplesSGPPost = sgp_zs.iid(8).sample().map(s => s.toBreezeVector.toArray.toSeq)
+val samplesSGPPost = sgp_zs.iid(8).sample().map(_.toBreezeVector.toArray.toSeq)
 
 spline(xs, samplesSGPPost.head)
 hold()
 samplesSGPPost.tail.foreach((s: Seq[Double]) => spline(xs, s))
 unhold()
 title("Ext. Skew Gaussian Process posterior samples")
 
+val samplesMixPost = mix_zs.iid(8).sample().map(_.toBreezeVector.toArray.toSeq)
+
+spline(xs, samplesMixPost.head)
+hold()
+samplesMixPost.tail.foreach((s: Seq[Double]) => spline(xs, s))
+unhold()
+title("Gaussian Process Mixture posterior samples")
+
+
 val (dx, dy) = dataset.sorted.unzip
 
 spline(dx, dy)