[WIP] Adding Warped Gaussian Process model

Author: mandar2812

Diff for: dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/analysis/PushforwardMap.scala

@@ -18,6 +18,7 @@ under the License.
 * */
 package io.github.mandar2812.dynaml.analysis
 
+import io.github.mandar2812.dynaml.models.ContinuousProcess
 import io.github.mandar2812.dynaml.pipes.{DataPipe, Encoder}
 import io.github.mandar2812.dynaml.probability.{ContinuousDistrRV, MeasurableDistrRV}
 import spire.algebra.Field
@@ -43,6 +44,8 @@ abstract class PushforwardMap[
   def ->[R <: ContinuousDistrRV[Source]](r: R): MeasurableDistrRV[Source, Destination, Jacobian] =
     new MeasurableDistrRV[Source, Destination, Jacobian](r)(self)
 
+  //def ->[P <: ContinuousProcess[_, _, ]]
+
   /**
     * Return the implementation used for calculating
     * determinant of the inverse function's Jacobian

Diff for: dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/models/Model.scala

@@ -45,6 +45,8 @@ trait Model[T, Q, R] {
     * */
   protected val g: T
 
+  def data = g
+
   /**
     * Predict the value of the
     * target variable given a

Diff for: dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/models/gp/WarpedGP.scala

@@ -1,15 +1,118 @@
 package io.github.mandar2812.dynaml.models.gp
 
-import io.github.mandar2812.dynaml.kernels.LocalScalarKernel
+import breeze.linalg.{DenseMatrix, det, diag}
+import io.github.mandar2812.dynaml.algebra.{PartitionedMatrix, PartitionedVector}
+import io.github.mandar2812.dynaml.analysis.{DifferentiableMap, PartitionedVectorField, PushforwardMap}
+import io.github.mandar2812.dynaml.models.{ContinuousProcess, SecondOrderProcess}
+import io.github.mandar2812.dynaml.optimization.GloballyOptWithGrad
+import io.github.mandar2812.dynaml.pipes.DataPipe
+import io.github.mandar2812.dynaml.probability.{E, MeasurableDistrRV}
+import io.github.mandar2812.dynaml.utils
 
 import scala.reflect.ClassTag
 
 /**
   * Created by mandar on 02/01/2017.
   */
-abstract class WarpedGP[T, I](
-  cov: LocalScalarKernel[I], n: LocalScalarKernel[I],
-  data: T, num: Int)(implicit ev: ClassTag[I]) extends
-  AbstractGPRegressionModel[T, I](cov, n, data, num) {
+abstract class WarpedGP[T, I](p: AbstractGPRegressionModel[T, I])(
+  wFuncT: PushforwardMap[Double, Double, Double])(
+  implicit ev: ClassTag[I], pf: PartitionedVectorField)
+  extends ContinuousProcess[
+    T, I, Double,
+    MeasurableDistrRV[PartitionedVector, PartitionedVector, PartitionedMatrix]]
+  with SecondOrderProcess[
+    T, I, Double, Double, DenseMatrix[Double],
+    MeasurableDistrRV[PartitionedVector, PartitionedVector, PartitionedMatrix]]
+  with GloballyOptWithGrad {
 
+
+  //Define the default determinant implementation
+  implicit val detImpl = DataPipe(
+    (m: PartitionedMatrix) => m.filterBlocks(c => c._1 == c._2).map(c => det(c._2)).product)
+
+  //Define the push forward map for the multivariate case
+  val wFuncPredDistr: PushforwardMap[PartitionedVector, PartitionedVector, PartitionedMatrix] =
+    PushforwardMap(
+      DataPipe((v: PartitionedVector) => v.map(c => (c._1, c._2.map(wFuncT.run)))),
+      DifferentiableMap(
+        (v: PartitionedVector) => v.map(c => (c._1, c._2.map(wFuncT.i.run))),
+        (v: PartitionedVector) => new PartitionedMatrix(
+          v._data.map(l => ((l._1, l._1), diag(l._2.map(wFuncT.i.J)))) ++
+            utils.combine(Seq((0 until v.rowBlocks.toInt).toList, (0 until v.rowBlocks.toInt).toList))
+              .map(c =>
+                (c.head.toLong, c.last.toLong))
+              .filter(c => c._2 != c._1)
+              .map(c => (c, DenseMatrix.zeros[Double](v.rows.toInt/v.rowBlocks.toInt, v.rows.toInt/v.rowBlocks.toInt)))
+              .toStream, num_cols = v.rows, num_rows = v.rows))
+      )
+
+  /**
+    * Draw three predictions from the posterior predictive distribution
+    * 1) Mean or MAP estimate Y
+    * 2) Y- : The lower error bar estimate (mean - sigma*stdDeviation)
+    * 3) Y+ : The upper error bar. (mean + sigma*stdDeviation)
+    **/
+  override def predictionWithErrorBars[U <: Seq[I]](testData: U, sigma: Int) = ???
+
+  /**
+    * Mean Function: Takes a member of the index set (input)
+    * and returns the corresponding mean of the distribution
+    * corresponding to input.
+    **/
+  override val mean = p.mean
+  /**
+    * Underlying covariance function of the
+    * Gaussian Processes.
+    **/
+  override val covariance = p.covariance
+  /**
+    * Stores the names of the hyper-parameters
+    **/
+  override protected var hyper_parameters: List[String] = p._hyper_parameters
+  /**
+    * A Map which stores the current state of
+    * the system.
+    **/
+  override protected var current_state: Map[String, Double] = p._current_state
+
+  /**
+    * Calculates the energy of the configuration,
+    * in most global optimization algorithms
+    * we aim to find an approximate value of
+    * the hyper-parameters such that this function
+    * is minimized.
+    *
+    * @param h       The value of the hyper-parameters in the configuration space
+    * @param options Optional parameters about configuration
+    * @return Configuration Energy E(h)
+    **/
+  override def energy(h: Map[String, Double], options: Map[String, String]) = p.energy(h, options)
+
+  /** Calculates posterior predictive distribution for
+    * a particular set of test data points.
+    *
+    * @param test A Sequence or Sequence like data structure
+    *             storing the values of the input patters.
+    **/
+  override def predictiveDistribution[U <: Seq[I]](test: U) = wFuncPredDistr -> p.predictiveDistribution(test)
+
+  /**
+    * Convert from the underlying data structure to
+    * Seq[(I, Y)] where I is the index set of the GP
+    * and Y is the value/label type.
+    **/
+  override def dataAsSeq(data: T) = p.dataAsSeq(data)
+
+  /**
+    * The training data
+    **/
+  override protected val g: T = p.data
+
+  /**
+    * Predict the value of the
+    * target variable given a
+    * point.
+    *
+    **/
+  override def predict(point: I) = wFuncT(p.predictionWithErrorBars(Seq(point), 1).head._2)
 }

Diff for: dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/optimization/GloballyOptimizable.scala

@@ -39,6 +39,10 @@ trait GloballyOptimizable {
    * */
   protected var current_state: Map[String, Double]
 
+  def _current_state = current_state
+
+  def _hyper_parameters = hyper_parameters
+
   /**
    * Calculates the energy of the configuration,
    * in most global optimization algorithms

Diff for: dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/probability/GaussianRV.scala

@@ -8,17 +8,19 @@ import io.github.mandar2812.dynaml.probability.distributions.BlockedMultiVariate
 import spire.implicits._
 import spire.algebra.Field
 
+abstract class AbstractGaussianRV[T, V] extends ContinuousDistrRV[T]
+
 /**
   * Created by mandar on 26/7/16.
   */
-case class GaussianRV(mu: Double, sigma: Double) extends ContinuousDistrRV[Double] {
+case class GaussianRV(mu: Double, sigma: Double) extends AbstractGaussianRV[Double, Double] {
   override val underlyingDist = new Gaussian(mu, sigma)
 }
 
 case class MultGaussianRV(
   mu: DenseVector[Double], covariance: DenseMatrix[Double])(
   implicit ev: Field[DenseVector[Double]])
-  extends ContinuousDistrRV[DenseVector[Double]] {
+  extends AbstractGaussianRV[DenseVector[Double], DenseMatrix[Double]] {
 
   override val underlyingDist = MultivariateGaussian(mu, covariance)
 
@@ -39,7 +41,7 @@ object MultGaussianRV {
 
 case class MultGaussianPRV(mu: PartitionedVector, covariance: PartitionedPSDMatrix)(
   implicit ev: Field[PartitionedVector])
-  extends ContinuousDistrRV[PartitionedVector] {
+  extends AbstractGaussianRV[PartitionedVector, PartitionedPSDMatrix] {
 
   override val underlyingDist: ContinuousDistr[PartitionedVector] = BlockedMultiVariateGaussian(mu, covariance)
 

Diff for: dynaml-core/src/main/scala-2.11/io/github/mandar2812/dynaml/probability/MeasurableFunction.scala

@@ -77,7 +77,7 @@ object RealValuedMeasurableFunction {
 class MeasurableDistrRV[Domain1, Domain2, Jacobian](
   baseRV: ContinuousDistrRV[Domain1])(p: PushforwardMap[Domain1, Domain2, Jacobian])
   extends MeasurableFunction[Domain1, Domain2](baseRV)(p)
-    with RandomVarWithDistr[Domain2, ContinuousDistr[Domain2]] {
+    with ContinuousDistrRV[Domain2] {
 
   override val underlyingDist = new ContinuousDistr[Domain2] {
     override def unnormalizedLogPdf(x: Domain2) =