Merge branch 'fixsaltelli' into 'master'

Fixsaltelli

See merge request !23

openmole/build.sbt

@@ -467,7 +467,7 @@ lazy val abc = OsgiProject(pluginDir, "org.openmole.plugin.method.abc", imports
 
 lazy val directSampling = OsgiProject(pluginDir, "org.openmole.plugin.method.directsampling", imports = Seq("*")) dependsOn(openmoleDSL, distributionDomain, pattern, modifierDomain, fileHook) settings (pluginSettings: _*)
 
-lazy val sensitivity = OsgiProject(pluginDir, "org.openmole.plugin.method.sensitivity", imports = Seq("*")) dependsOn(exception, workflow, workspace, openmoleDSL, lhsSampling, directSampling) settings (pluginSettings: _*)
+lazy val sensitivity = OsgiProject(pluginDir, "org.openmole.plugin.method.sensitivity", imports = Seq("*")) dependsOn(exception, workflow, workspace, openmoleDSL, lhsSampling, directSampling, collectionDomain % "test", boundsDomain % "test") settings (pluginSettings: _*)
 
 
 /* Sampling */

openmole/plugins/org.openmole.plugin.method.sensitivity/src/main/scala/org/openmole/plugin/method/sensitivity/MorrisAggregation.scala

@@ -112,9 +112,17 @@ object MorrisAggregation {
     modelInputs:  Seq[ScalarOrSequenceOfDouble[_]],
     modelOutputs: Seq[Val[Double]])(implicit name: sourcecode.Name, definitionScope: DefinitionScope) = {
 
-    val muOutputs = Sensitivity.outputs(modelInputs, modelOutputs).map { case (i, o) ⇒ Morris.mu(i, o) }
-    val muStarOutputs = Sensitivity.outputs(modelInputs, modelOutputs).map { case (i, o) ⇒ Morris.muStar(i, o) }
-    val sigmaOutputs = Sensitivity.outputs(modelInputs, modelOutputs).map { case (i, o) ⇒ Morris.sigma(i, o) }
+    def morrisOutputs(
+      modelInputs:  Seq[ScalarOrSequenceOfDouble[_]],
+      modelOutputs: Seq[Val[Double]]) =
+      for {
+        i ← ScalarOrSequenceOfDouble.prototypes(modelInputs)
+        o ← modelOutputs
+      } yield (i, o)
+
+    val muOutputs = morrisOutputs(modelInputs, modelOutputs).map { case (i, o) ⇒ Morris.mu(i, o) }
+    val muStarOutputs = morrisOutputs(modelInputs, modelOutputs).map { case (i, o) ⇒ Morris.muStar(i, o) }
+    val sigmaOutputs = morrisOutputs(modelInputs, modelOutputs).map { case (i, o) ⇒ Morris.sigma(i, o) }
 
     FromContextTask("MorrisAggregation") { p ⇒
       import p._
@@ -126,7 +134,7 @@ object MorrisAggregation {
       // for each part of the space we were asked to explore, compute the elementary effects and returns them
       // into the variables passed by the user
       val List(mu, muStar, sigma) =
-        Sensitivity.outputs(modelInputs, modelOutputs).map {
+        morrisOutputs(modelInputs, modelOutputs).map {
           case (input, output) ⇒
             val outputValues: Array[Double] = context(output.toArray)
             MorrisSampling.Log.logger.fine("Processing the elementary change for input " + input + " on " + output)

openmole/plugins/org.openmole.plugin.method.sensitivity/src/main/scala/org/openmole/plugin/method/sensitivity/SaltelliAggregation.scala

@@ -28,7 +28,9 @@ object SaltelliAggregation {
     val NB = fB.size
     val k = fC.size //Number of parameters
     val f02 = math.pow(fB.sum / NB.toDouble, 2)
-    val varY = fB.map(fBj ⇒ math.pow(fBj, 2)).sum / NB.toDouble - f02
+    // val varY = fB.map(fBj ⇒ math.pow(fBj, 2)).sum / NB.toDouble - f02
+    val f0 = fB.sum / NB.toDouble
+    val varY = fB.map(fBj ⇒ math.pow(fBj - f0, 2)).sum / NB.toDouble
     def avgProduct(u: Array[Double], v: Array[Double]): Double = {
       val prods = (u zip v).map({ case (uj, vj) ⇒ uj * vj })
       prods.sum / prods.size.toDouble
@@ -63,8 +65,16 @@ object SaltelliAggregation {
     firstOrderSI: Val[Array[Array[Double]]]        = Val[Array[Array[Double]]]("firstOrderSI"),
     totalOrderSI: Val[Array[Array[Double]]]        = Val[Array[Array[Double]]]("totalOrderSI"))(implicit name: sourcecode.Name, definitionScope: DefinitionScope) = {
 
-    val fOOutputs = Sensitivity.outputs(modelInputs, modelOutputs).map { case (i, o) ⇒ Saltelli.firstOrder(i, o) }
-    val tOOutputs = Sensitivity.outputs(modelInputs, modelOutputs).map { case (i, o) ⇒ Saltelli.totalOrder(i, o) }
+    def saltelliOutputs(
+      modelInputs:  Seq[ScalarOrSequenceOfDouble[_]],
+      modelOutputs: Seq[Val[Double]]) =
+      for {
+        o ← modelOutputs
+        i ← ScalarOrSequenceOfDouble.prototypes(modelInputs)
+      } yield (i, o)
+
+    val fOOutputs = saltelliOutputs(modelInputs, modelOutputs).map { case (i, o) ⇒ Saltelli.firstOrder(i, o) }
+    val tOOutputs = saltelliOutputs(modelInputs, modelOutputs).map { case (i, o) ⇒ Saltelli.totalOrder(i, o) }
 
     FromContextTask("SaltelliAggregation") { p ⇒
       import p._
@@ -93,11 +103,17 @@ object SaltelliAggregation {
       val fC: Array[Array[Array[Double]]] =
         modelInputs.map { i ⇒ reindex("c$" ++ i.prototype.name) }.toArray.transpose
 
+      // ftoi(o)._1(i) contains first order index for input i on output o.
+      // ftoi(o)._2(i) contains total order index for input i on output o.
       val ftoi: Array[(Array[Double], Array[Double])] =
         (fA zip fB zip fC).map { case ((fAo, fBo), fCo) ⇒ firstAndTotalOrderIndices(fAo, fBo, fCo) }
 
       // first order indices
+      // fosi(o)(i) contains first order index for input i on output o.
       val fosi = ftoi.map { _._1.toArray }.toArray
+
+      // total order indices
+      // tosi(o)(i) contains total order index for input i on output o.
       val tosi = ftoi.map { _._2.toArray }.toArray
 
       val fosiv =

openmole/plugins/org.openmole.plugin.method.sensitivity/src/main/scala/org/openmole/plugin/method/sensitivity/SaltelliHook.scala

@@ -19,7 +19,7 @@ object SaltelliHook {
             Sensitivity.writeResults(ps, inputs, dsl.data.outputs, Saltelli.firstOrder(_, _)).from(context)
           }
 
-          (dir / "firstOrderTotalOrder.csv").withPrintStream(overwrite = true) { ps ⇒
+          (dir / "totalOrderIndices.csv").withPrintStream(overwrite = true) { ps ⇒
             Sensitivity.writeResults(ps, inputs, dsl.data.outputs, Saltelli.totalOrder(_, _)).from(context)
           }
         case WritableOutput.PrintStreamValue(ps) ⇒
@@ -34,4 +34,4 @@ object SaltelliHook {
       context
     }
 
-}
\ No newline at end of file
+}

openmole/plugins/org.openmole.plugin.method.sensitivity/src/main/scala/org/openmole/plugin/method/sensitivity/SaltelliSampling.scala

@@ -19,7 +19,7 @@ object SaltelliSampling {
 
   def buildLineOfC(i: Int, lineOfA: Array[Double], lineOfB: Array[Double]) =
     (lineOfA zip lineOfB zipWithIndex) map {
-      case ((a, b), index) ⇒ if (index == i) a else b
+      case ((a, b), index) ⇒ if (index == i) b else a
     }
 
   val namespace = Namespace("saltelli")
@@ -79,5 +79,4 @@ class SaltelliSampling(val samples: FromContext[Int], val factors: ScalarOrSeque
 
     (aVariables ++ bVariables ++ cVariables).toIterator
   }
-
 }

openmole/plugins/org.openmole.plugin.method.sensitivity/src/main/scala/org/openmole/plugin/method/sensitivity/package.scala

@@ -63,14 +63,6 @@ package object sensitivity {
       }
     }
 
-    def outputs(
-      modelInputs: Seq[ScalarOrSequenceOfDouble[_]],
-      modelOutputs: Seq[Val[Double]]) =
-      for {
-        i ← ScalarOrSequenceOfDouble.prototypes(modelInputs)
-        o ← modelOutputs
-      } yield (i, o)
-
 
     case class SaltelliParams(inputs: Seq[ScalarOrSequenceOfDouble[_]], outputs: Seq[Val[_]])
     case class MorrisParams(inputs: Seq[ScalarOrSequenceOfDouble[_]], outputs: Seq[Val[_]])

openmole/plugins/org.openmole.plugin.method.sensitivity/src/test/scala/org/openmole/plugin/method/sensitivity/SaltelliSpec.scala

+package org.openmole.plugin.method.sensitivity
+
+import org.scalatest.{ FlatSpec, Matchers }
+import org.openmole.core.dsl._
+import org.openmole.core.workflow.test._
+import scala.util.Random
+
+import org.openmole.core.dsl._
+import org.openmole.core.workflow.mole.Mole
+import org.openmole.core.workflow.task.FromContextTask
+import org.openmole.core.workflow.validation._
+import org.openmole.plugin.domain.collection._
+import org.scalatest._
+import org.openmole.plugin.domain.bounds._
+
+class Saltelli extends FlatSpec with Matchers {
+
+  import org.openmole.core.workflow.test.Stubs._
+
+  val rng = new Random()
+
+  val x1 = Val[Double]
+  val x2 = Val[Double]
+  val y1 = Val[Double]
+  val y2 = Val[Double]
+  val y3 = Val[Double]
+
+  /* Expected values of first order (SI) and total order (STI) sensitivity indices.
+   *
+   * - y1, x1: SI1 = 4/5, STI1 = 4/5
+   * - y1, x2: SI2 = 1/5, STI2 = 1/5
+   * - y2, x1: SI1 = (9 / 4) * (12 / 42) ~= 0.643,
+   *           STI1 = (7.0 / 36.0) / (40.0 / 144.0) = 0.7
+   * - y2, x2: SI2 = 12 / 42 ~= 0.286,
+   *           STI2 = (13.0 / 144.0) / (40.0 / 144.0) ~= 0.325
+   */
+
+  val model = TestTask { context ⇒
+    context +
+      (y1 -> (context(x1) + 0.5 * context(x2))) +
+      (y2 -> (context(x1) + 0.5 * context(x2) + context(x1) * context(x2))) +
+      (y3 -> (context(x1)))
+  } set (
+    inputs += (x1, x2),
+    outputs += (y1, y2, y3)
+  )
+
+  val sen = SensitivitySaltelli(
+    evaluation = model,
+    samples = 10000,
+    inputs = Seq(x1 in (0.0, 1.0), x2 in (0.0, 1.0)),
+    outputs = Seq(y1, y2, y3)
+  )
+
+  "Saltelli" should "run" in {
+    sen run ()
+  }
+
+}
+