Wavelet solver

build.sbt

@@ -30,7 +30,9 @@ libraryDependencies ++= Seq(
   "com.fasterxml.jackson.core" % "jackson-databind" % "2.5.3",
   "com.github.tototoshi" %% "scala-csv" % "1.3.10",
   "org.apache.xmlgraphics" % "batik-all" % "1.16",
-  "com.fasterxml.jackson.module" % "jackson-module-scala_2.12" % "2.8.8")
+  "com.fasterxml.jackson.module" % "jackson-module-scala_2.12" % "2.8.8",
+  "de.sciss" % "jwave" % "1.0.3",
+)
 
 lazy val originalCBackend = (project in file("src") / "native" / "Original")
   .settings(nativeCompile / sourceDirectory := baseDirectory.value)

src/main/scala/core/solver/wavelet/CoefficientAveragingWaveletSolver.scala

@@ -0,0 +1,117 @@
+package core.solver.wavelet
+
+import util.Profiler
+
+class CoefficientAveragingWaveletSolver(override val querySize: Int, override val debug: Boolean = false)
+  extends
+    WaveletSolver[Double]("CoefficientAveragingWaveletSolver", querySize, debug) {
+
+
+  override def solve(): Array[Double] = {
+    val finalWavelet: Array[Double] = new Array(N)
+
+    // target order is ascending variables
+    val finalCuboidVariables = (0 until querySize).sorted(Ordering[Int])
+
+    // multiple cuboids can have the same variable,
+    // each corresponding wavelet will have coefficients for that variable
+    // we need to combine or select the coefficients from all the wavelets
+
+    // here we will average the coefficients for each variable,
+    // after extrapolating and permuting the marginal wavelets into the target order
+
+    // for each variable, find the number of cuboids that it appears in
+    val variableCounts: Map[Int, Int] = Profiler(s"$solverName Coefficient Averaging Preparation") {
+      cuboids.keys.flatten.groupBy(identity).mapValues(_.size)
+    }
+
+    for (((variables, cuboid), index) <- cuboids.zipWithIndex) {
+      val wavelet: Array[Double] = Profiler(s"$solverName Forward Transform") {
+        // transform cuboid to wavelet
+        transformer.forward(cuboid)
+      }
+
+      val extrapolatedWavelet: Array[Double] = Profiler(s"$solverName Wavelet Extrapolation") {
+        // extrapolate wavelet to the correct size, by filling the missing values with 0
+        val extrapolatedWavelet = Array.fill(N)(0.0)
+        System.arraycopy(wavelet, 0, extrapolatedWavelet, 0, wavelet.length)
+
+        extrapolatedWavelet
+      }
+
+      // copy the coefficients to the final wavelet
+      Profiler(s"$solverName Wavelet Permutation & Combination") {
+        // copy the sum to the final wavelet
+        if (index == 0) {
+          finalWavelet(0) = extrapolatedWavelet(0)
+        }
+
+        // fill the wavelet values into the right positions in the final wavelet
+
+        /* eg:
+        querySize = 3
+        extrapolatedWavelet (0, 2) :: (a, b, c, d)
+
+        for variable = 2 (index 0): variablePos = 2^0 = 1, numCoefficients = 2^0 = 1, variableCoefficients = (b)
+        for variable = 0 (index 1): variablePos = 2^1 = 2, numCoefficients = 2^1 = 2, variableCoefficients = (c, d)
+
+        in final wavelet (0, 1, 2) :: (t, 1, 2, 3, 4, 5, 6, 7)
+
+        for variable = 2: variablePos = 2^0 = 1, numCoefficients = 2^0 = 1, finalWavelet(1) = (b)
+        for variable = 0: variablePos = 2^2 = 4, numCoefficients = 2^2 = 4, finalWavelet([4...7]) = (c/2, c/2, d/2, d/2)
+        */
+        // for each variable
+        variables.reverse.zipWithIndex.foreach { case (variable, i) =>
+          val variablePos = Math.pow(2, i).toInt
+          val numCoefficients = variablePos
+
+          // get all the coefficients of the variable
+          val variableCoefficients = extrapolatedWavelet.slice(variablePos, variablePos + numCoefficients)
+
+          // set the coefficients of the variable in the final wavelet
+          val finalVariablePos = Math.pow(2, (querySize - 1) - variable).toInt
+          val finalNumCoefficients = finalVariablePos
+
+          // average the coefficients for each variable
+          var adjustedCoefficients = variableCoefficients.map { c => c / variableCounts(variable) }
+
+          if (numCoefficients < finalNumCoefficients) {
+            // we need to equally split the coefficients, giving us double the number of coefficients
+            val divisor = finalNumCoefficients / numCoefficients
+            adjustedCoefficients = variableCoefficients.flatMap(c => Array.fill(divisor)(c / divisor))
+
+          } /* else if (numCoefficient > finalNumCoefficients)
+            we need to sum the coefficients pairwise, giving us half the number of coefficients
+
+            eg. given source cuboid (2, 4, 5) :: [total, 5, 4_0, 4_1, 2_0, 2_1, 2_2, 2_3]
+            in the final cuboid (..., 2, ..., 4, 5, ...)
+            there will always be equal or more variables succeeding variable 2
+            hence, 2 can never acquire an earlier position that requires less coefficients than it has
+           */
+
+          System.arraycopy(adjustedCoefficients, 0, finalWavelet, finalVariablePos, finalNumCoefficients)
+        }
+      }
+
+      if (debug) {
+        println(s"Source Cuboid: $variables :: ${cuboid.mkString(", ")}")
+        println(s"Source Wavelet:  $variables :: ${extrapolatedWavelet.mkString(", ")}")
+        println(s"Combined Wavelet: $finalCuboidVariables :: ${finalWavelet.mkString(", ")}")
+        println()
+      }
+    }
+
+    val finalCuboid: Array[Double] = Profiler(s"$solverName Reverse Transform") {
+      transformer.reverse(finalWavelet)
+    }
+
+    if (debug) {
+      println(s"Final Wavelet: $finalCuboidVariables :: ${finalWavelet.mkString(", ")}")
+      println(s"Final Cuboid: $finalCuboidVariables :: ${finalCuboid.mkString(", ")}")
+      println()
+    }
+
+    solution = finalCuboid
+    finalCuboid
+  }
+}

src/main/scala/core/solver/wavelet/CoefficientSelectionWaveletSolver.scala

@@ -0,0 +1,144 @@
+package core.solver.wavelet
+
+import util.Profiler
+
+import scala.math.Ordering.Implicits.seqDerivedOrdering
+
+class CoefficientSelectionWaveletSolver(override val querySize: Int,
+                                        override val debug: Boolean = false)
+  extends
+    WaveletSolver[Double]("CoefficientSelectionWaveletSolver", querySize, debug) {
+
+
+  override def solve(): Array[Double] = {
+    val finalWavelet: Array[Double] = new Array(N)
+
+    // target order is ascending variables
+    val finalCuboidVariables = (0 until querySize).sorted(Ordering[Int])
+
+    // multiple cuboids can have the same variable,
+    // each corresponding wavelet will have coefficients for that variable
+    // we need to combine or select the coefficients from all the wavelets
+
+    // here we sort the stored cuboids by the variable seq in ascending order
+    // eg: if we have [[0, 1, 4], [0, 1, 2], [1], [2, 3, 4], [2, 4], [1, 2]]
+    // we will sort them to [[0, 1, 2], [0, 1, 4], [1], [1, 2], [2, 3, 4], [2, 4]]
+    // then do a stable sort on the length of the element seq
+    // we will finally get [[1], [1, 2], [2, 4], [0, 1, 2], [0, 1, 4], [2, 3, 4]]
+    // thus the lower dimensional marginals will be overwritten by the higher dimensional marginals
+    val sortedCuboids = Profiler(s"$solverName Coefficient Priority Selection") {
+      cuboids.toSeq
+        .sortBy(_._1)(seqDerivedOrdering)
+        .sortBy(_._1.length)(Ordering[Int])
+    }
+
+    if (debug) {
+      println(s"Coefficient priority: ")
+      println(s"\t${
+        sortedCuboids.reverse.zipWithIndex.map {
+          case ((variables, _), index) => s"$index : (${variables.mkString(",")})"
+        }.mkString("\n\t")
+      }")
+      println()
+    }
+
+    for (((variables, cuboid), index) <- sortedCuboids.zipWithIndex) {
+      val wavelet: Array[Double] = Profiler(s"$solverName Forward Transform") {
+        // transform cuboid to wavelet
+        transformer.forward(cuboid)
+      }
+
+      val extrapolatedWavelet: Array[Double] = Profiler(s"$solverName Wavelet Extrapolation") {
+        // extrapolate wavelet to the correct size, by filling the missing values with 0
+        val extrapolatedWavelet = Array.fill(N)(0.0)
+        System.arraycopy(wavelet, 0, extrapolatedWavelet, 0, wavelet.length)
+
+        extrapolatedWavelet
+      }
+
+
+
+      // copy the coefficients to the final wavelet
+      Profiler(s"$solverName Wavelet Permutation & Combination") {
+        // copy the sum to the final wavelet
+        if (index == 0) {
+          finalWavelet(0) = extrapolatedWavelet(0)
+        }
+
+        // fill the wavelet values into the right positions in the final wavelet
+
+        /* eg:
+        querySize = 3
+        extrapolatedWavelet (0, 2) :: (a, b, c, d)
+
+        for variable = 2 (index 0): variablePos = 2^0 = 1, numCoefficients = 2^0 = 1, variableCoefficients = (b)
+        for variable = 0 (index 1): variablePos = 2^1 = 2, numCoefficients = 2^1 = 2, variableCoefficients = (c, d)
+
+        in final wavelet (0, 1, 2) :: (t, 1, 2, 3, 4, 5, 6, 7)
+
+        for variable = 2: variablePos = 2^0 = 1, numCoefficients = 2^0 = 1, finalWavelet(1) = (b)
+        for variable = 0: variablePos = 2^2 = 4, numCoefficients = 2^2 = 4, finalWavelet([4...7]) = (c/2, c/2, d/2, d/2)
+        */
+        // for each variable
+        variables.reverse.zipWithIndex.foreach { case (variable, i) =>
+          val variablePos = Math.pow(2, i).toInt
+          val numCoefficients = variablePos
+
+          // get all the coefficients of the variable
+          val variableCoefficients = extrapolatedWavelet.slice(variablePos, variablePos + numCoefficients)
+
+          // set the coefficients of the variable in the final wavelet
+          val finalVariablePos = Math.pow(2, (querySize - 1) - variable).toInt
+          val finalNumCoefficients = finalVariablePos
+
+          var adjustedCoefficients = Array.empty[Double]
+
+          if (numCoefficients == finalNumCoefficients) {
+            // we can just copy the coefficients
+            adjustedCoefficients = variableCoefficients
+
+          } else if (numCoefficients < finalNumCoefficients) {
+            // we need to equally split the coefficients, giving us double the number of coefficients
+            val divisor = finalNumCoefficients / numCoefficients
+            adjustedCoefficients = variableCoefficients.flatMap(c => Array.fill(divisor)(c / divisor))
+
+          } else /* if numCoefficient > finalNumCoefficients  */ {
+            // we need to sum the coefficients pairwise, giving us half the number of coefficients
+
+            /*
+            eg. given marginal source cuboid (2, 4, 5) :: [total, 5, 4_0, 4_1, 2_0, 2_1, 2_2, 2_3]
+            the final cuboid will be (..., 2, ..., 4, 5, ...)
+            In the final cuboid, 2 cannot have less number of variables succeeding it than it has in the marginal cuboid.
+            Hence, 2 will never require less number of coefficients than it has in the marginal cuboid.
+            */
+            assert(assertion = false,
+              "Unexpected condition: numCoefficients > finalNumCoefficients for \n" +
+                s"MarginalCuboid: ${variables.mkString(", ")}")
+          }
+
+          System.arraycopy(adjustedCoefficients, 0, finalWavelet, finalVariablePos, finalNumCoefficients)
+        }
+      }
+
+      if (debug) {
+        println(s"Source Cuboid: $variables :: ${cuboid.mkString(", ")}")
+        println(s"Source Wavelet:  $variables :: ${extrapolatedWavelet.mkString(", ")}")
+        println(s"Combined Wavelet: $finalCuboidVariables :: ${finalWavelet.mkString(", ")}")
+        println()
+      }
+    }
+
+    val finalCuboid: Array[Double] = Profiler(s"$solverName Reverse Transform") {
+      transformer.reverse(finalWavelet)
+    }
+
+    if (debug) {
+      println(s"Final Wavelet: $finalCuboidVariables :: ${finalWavelet.mkString(", ")}")
+      println(s"Final Cuboid: $finalCuboidVariables :: ${finalCuboid.mkString(", ")}")
+      println()
+    }
+
+    solution = finalCuboid
+    finalCuboid
+  }
+}

src/main/scala/core/solver/wavelet/Demo.scala

@@ -0,0 +1,120 @@
+package core.solver.wavelet
+
+import core.solver.SolverTools
+import core.solver.iterativeProportionalFittingSolver.IPFUtils
+import util.BitUtils
+
+import scala.util.Random
+
+object Demo {
+
+  /**
+   * Run the solver with the given parameters, and print the actual data, calculated result and error
+   *
+   * @param cuboidIndices list of cuboids to provide to the solver written as Seq(marginal_variables)
+   * @param solver        solver to use
+   */
+  def run(cuboidIndices: Seq[Seq[Int]], solver: WaveletSolver[Double]): Unit = {
+    val rng = new Random(42)
+
+    val querySize = solver.querySize
+
+    // Generate random result
+    val actual: Array[Double] = Array.fill(1 << querySize)(rng.nextInt(100))
+
+    implicit def seqToInt: Seq[Int] => Int = BitUtils.SetToInt
+
+    // List of cuboids represented as (column indices -> flattened array of values)
+    val cuboids: Map[Seq[Int], Array[Double]] =
+      cuboidIndices.map { marginalVariables =>
+        marginalVariables -> IPFUtils.getMarginalDistribution(querySize, actual, marginalVariables.size, marginalVariables)
+      }.toMap
+
+    // Add cuboids to solver
+    cuboids.foreach { case (columnIndices, values) => solver.addCuboid(columnIndices, values) }
+
+    // Solve
+    val result: Array[Double] = solver.solve()
+
+    // print list of doubles with 2 decimal places
+    println(s"Actual: ${PrintUtils.toString(actual)}")
+    println(s"Result: ${PrintUtils.toString(result)}")
+    println(s"Error: ${SolverTools.error(actual, result)}")
+  }
+}
+
+object TransformDemo {
+
+  def main(args: Array[String]): Unit = {
+    val t = new HaarTransformer[Double]()
+
+    val a = Array(1.0, 5.0, 11.0, 2.0)
+
+    val w = t.forward(a)
+
+    val a_ = t.reverse(w)
+
+    println(s"Original: ${PrintUtils.toString(a)}")
+    println(s"Wavelet: ${PrintUtils.toString(w)}")
+    println(s"Reversed: ${PrintUtils.toString(a_)}")
+  }
+}
+
+
+object Demo1 {
+  def main(args: Array[String]): Unit = {
+    Demo.run(
+      Seq(Seq(0, 1, 4), Seq(1, 5), Seq(2, 3, 5), Seq(0, 2), Seq(3, 4), Seq(1)),
+      new CoefficientSelectionWaveletSolver(6, debug = true))
+  }
+}
+
+object Demo2 {
+  def main(args: Array[String]): Unit = {
+    Demo.run(
+      Seq(Seq(0, 1, 4), Seq(1, 5), Seq(2, 3, 5), Seq(0, 2), Seq(3, 4), Seq(1)),
+      new CoefficientAveragingWaveletSolver(6, debug = true))
+  }
+}
+
+object Demo3 {
+
+  def main(args: Array[String]): Unit = {
+
+    val runs: Seq[Run] = Seq(
+      //      Run(3, Seq((Seq(0, 1, 2), Array(1.0, 19.0, 13.0, 1.0, 4.0, 7.0, 1.0, 13.0))), new SingleCuboidWaveletSolver(_, _)),
+      //      Run(3, Seq((Seq(0, 1), Array(5.0, 26.0, 14.0, 14.0))), new SingleCuboidWaveletSolver(_, _)),
+      //      Run(3, Seq((Seq(0, 2), Array(14.0, 20.0, 5.0, 20.0))), new SingleCuboidWaveletSolver(_, _)),
+      //      Run(3, Seq((Seq(1), Array(31.0, 28.0))), new SingleCuboidWaveletSolver(_, _)),
+      //      Run(3, Seq((Seq(0), Array(19.0, 40.0))), new SingleCuboidWaveletSolver(_, _)),
+      //      Run(3,
+      //        Seq(
+      //          (Seq(1), Array(31.0, 28.0)),
+      //          (Seq(0, 2), Array(14.0, 20.0, 5.0, 20.0))
+      //        ),
+      //        new MutuallyExclusiveWaveletSolver(_, _)
+      //      ),
+      //      Run(3,
+      //        Seq(
+      //          (Seq(0, 1, 2), Array(1.0, 19.0, 13.0, 1.0, 4.0, 7.0, 1.0, 13.0)),
+      //        ),
+      //        new MutuallyExclusiveWaveletSolver(_, _)
+      //      ),
+      Run(3, Seq((Seq(0, 1, 2), Array(1.0, 5.0, 11.0, 7.0, 12.0, 4.0, 6.0, 3.0))), new SingleCuboidWaveletSolver(_, _)),
+
+
+    )
+
+    runs.foreach { case Run(querySize, cuboids, createSolver) =>
+      val solver = createSolver(querySize, true)
+      cuboids.foreach { case (columnIndices, values) =>
+        solver.addCuboid(columnIndices, values)
+      }
+      solver.solve()
+
+      println("================================")
+    }
+  }
+
+  private case class Run(querySize: Int, cuboids: Seq[(Seq[Int], Array[Double])], solver: (Int, Boolean) => WaveletSolver[Double])
+}