sensitivityALUMSS.oms

// to execute from the command line: openmole -s sensitivityALUMSS.oms

// output files
val agreOut = Val[File]
val f = Val[File]

// random number generator seed
val mySeed = Val[Int]

// model parameters
val ksi = Val[Double]
val y0 = Val[Double]
val sar = Val[Double]
val a = Val[Double]
val w = Val[Double]
val Tag = Val[Double]
val Tab = Val[Double]
val Tr = Val[Double]
val Td = Val[Double]
val d = Val[Double]

// simulation parameters
val T = Val[Double]
val dtp = Val[Double]
val L = Val[Double]
val d0 = Val[Double]
val a0 = Val[Double]
val dtsave = Val[Double]

// variable outputs
val t = Val[Double]
val nFrag = Val[Double]
val meanSize = Val[Double]
val stdSize = Val[Double]
val maxSize = Val[Double]
val meanES = Val[Double]
val stdES = Val[Double]
val P = Val[Double]
val N = Val[Double]
val D = Val[Double]
val A0 = Val[Double]
val A1 = Val[Double]
val Y = Val[Double]
val C = Val[Double]
val nMax = Val[Double]
val nMin = Val[Double]
val pMax = Val[Double]
val pMin = Val[Double]
val ripleyN = Val[Double]
val ripleyD = Val[Double]
val ripleyA0 = Val[Double]
val ripleyA1 = Val[Double]
val corrLen = Val[Double]


///////////////////////////////////////////////////////////////////////////////
// 1- model plugin with all the parameter inputs as well as the output files
///////////////////////////////////////////////////////////////////////////////

val scriptTask =
  SystemExecTask(
    command = Seq(
      "./alumss-exec ${T} ${dtp} ${L} ${a0} ${d0} ${ksi} ${y0} ${sar} ${a} ${w} ${Tag} ${Tab} ${Tr} ${Td} ${d} ${dtsave} ${mySeed} "),
    )set(
      inputs += (T,dtp,L,a0,d0,ksi,y0,sar,a,w,Tag,Tab,Tr,Td,d,dtsave,mySeed),
      // outputFiles += ("DATA_AGRE_T_${T}_dtp_${dtp}_n_${L}_a0_${a0}_d0_${d0}_ksi_${ksi}_y0_${y0}_sar_${sar}_a_${a}_w_${w}_Tag_${Tag}_Tab_${Tab}_Tr_${Tr}_Td_${Td}_d_${d}_dtsave_${dtsave}_expid_${mySeed}.dat",agreOut),
      outputFiles += ("DATA_AGRE",agreOut),
      resources+="/home/AD.LSM.CNRS.FR/diego.bengochea-paz/gillespie-land-use/alumss-exec",

      // default parameter and initialization values
      T:=5000.0,
      dtp:=0.1,
      L:=40,
      a0:=0.1,
      d0:=0.0,
      ksi:=3.0,
      y0:=0.2,
      sar:=0.25,
      a:=0.2,
      w:=0.0,
      Tag:=0.1,
      Tab:=50.0,
      Tr:=5.0,
      Td:=50.0,
      d:=1.0,
      dtsave:=7000.0,
      mySeed:=1111
    )

val readOutput =
    ScalaTask("""

    // read file as an array of lines and get last line
    val lastLineAgre = scala.io.Source.fromFile(agreOut).getLines.toList.last.mkString
    // split the string and store it in an array
    val valuesAgre = lastLineAgre.split(" ")

    val P = valuesAgre(1).toDouble
    val N = valuesAgre(2).toDouble
    val D = valuesAgre(3).toDouble
    val A0 = valuesAgre(4).toDouble
    val A1 = valuesAgre(5).toDouble
    // val Y = valuesAgre(6).toDouble
    val nFrag = valuesAgre(7).toDouble
    val meanSize = valuesAgre(8).toDouble
    val stdSize = valuesAgre(9).toDouble
    val maxSize = valuesAgre(10).toDouble
    val meanES = valuesAgre(11).toDouble
    val stdES = valuesAgre(12).toDouble
    // val C = valuesAgre(13).toDouble
    // val nMax = valuesAgre(14).toDouble
    // val nMin = valuesAgre(15).toDouble
    // val pMax = valuesAgre(16).toDouble
    // val pMin = valuesAgre(17).toDouble
    // val ripleyN = valuesAgre(18).toDouble
    // val ripleyD = valuesAgre(19).toDouble
    // val ripleyA0 = valuesAgre(20).toDouble
    // val ripleyA1 = valuesAgre(21).toDouble
    val corrLen = valuesAgre(24).toDouble

    """
    )set(
        inputs+=agreOut,
        outputs+=(P,N,D,A0,A1,nFrag,meanSize,stdSize,maxSize,meanES,stdES,corrLen)
    )

val moleTask = MoleTask(scriptTask -- readOutput)

val replications = Replication(
  evaluation = moleTask,
  seed = mySeed,
  sample = 10,
  aggregation = Seq (
    P aggregate average,
    N aggregate average,
    A0 aggregate average,
    A1 aggregate average,
    D aggregate average,
    nFrag aggregate average,
    meanSize aggregate average,
    stdSize aggregate average,
    maxSize aggregate average,
    meanES aggregate average,
    stdES aggregate average,
    corrLen aggregate average
  )
)

val env = LocalEnvironment(40)

// SensitivityMorris(
//     evaluation = moleTask on env,
//     inputs = Seq(
//         ksi in (0.5,2.0),
//         sar in (0.1,1.0),
//         a in (0.0,1.0),
//         w in (0.0,8.0),
//         Tag in (0.01,1.0),
//         Tab in (0.1,100.0),
//         Tr in (0.1,10.0),
//         Td in (1.0,100.0),
//         a0 in (0.1,0.9),
//         y0 in (0.05,1.0)
//         ),
//     outputs = Seq(
//         P,
//         N,
//         D,
//         A0,
//         A1,
//         nFrag,
//         meanSize,
//         stdSize,
//         maxSize,
//         meanES,
//         stdES,
//         ripleyN,
//         ripleyD,
//         ripleyA0,
//         ripleyA1
//         ),
//     sample = 100,
//     level = 100
// ) hook (workDirectory / "gillespieSensitivityResults6")

SensitivityMorris(
    evaluation = replications on env,
    inputs = Seq(
        sar in (0.05,1.0),
        a in (0.0,1.0),
        w in (0.0,10.0),
        Tag in (0.01,1.0),
        y0 in (0.01,2.0)
        ),
    outputs = Seq(
        P,
        N,
        A0,
        A1,
        D,
        nFrag,
        meanSize,
        stdSize,
        maxSize,
        meanES,
        stdES,
        corrLen
        ),
    sample = 1000,
    level = 20
) hook (workDirectory / "morrisSensitivityResults")

// SensitivitySaltelli(
//     evaluation = moleTask on env,
//     inputs = Seq(
//         ksi in (0.5,2.0),
//         sar in (0.1,1.0),
//         a in (0.0,1.0),
//         w in (0.0,8.0),
//         Tag in (0.01,1.0),
//         Tab in (0.1,100.0),
//         Tr in (0.1,10.0),
//         Td in (1.0,100.0),
//         a0 in (0.1,0.9),
//         y0 in (0.05,1.0)
//         ),
//     outputs = Seq(
//         P,
//         N,
//         D,
//         A0,
//         A1,
//         nFrag,
//         meanSize,
//         stdSize,
//         maxSize,
//         meanES,
//         stdES,
//         ripleyN,
//         ripleyD,
//         ripleyA0,
//         ripleyA1
//         ),
//     sample = 2000,
// ) hook (workDirectory / "gillespieSensitivityResultsSaltelli")