Use Case OpenMOLE : Effective management strategies for an invasive plant, the Japanese knotweed

Preamble

This file describes a simulation experiment. It includes both the
explanations and the source code used to generate the results. All the
files, scripts and data necessary to reproduce the experiment are
available in a git repository available here. Their execution requires the use of OpenMOLE and R. Detailed instructions for reproducing the results are available at the end of this document. The OpenMOLE or R executable script files are generated from the code blocks visible in this document with the script ./make.hs. Running it requires to have Stack installed. On linux or OS X, Stack can be installed with the command:

curl -sSL https://get.haskellstack.org/ | sh

This document can be read at two levels:

  • A scientific level that focuses on the numerical experiment process,
    for readers interested in modeling and who want to discover the type of scientific question that modeling, and particularly the OpenMOLE platform, can address.

  • A technical level focused on the implementation of the approach with OpenMOLE. This level of speech is contained in the blocks of code. It is intended for readers who have already been exposed to the OpenMOLE core concepts (workflow, task, hook, if needed, refer to the OpenMOLE documentation), wish to reproduce the study presented here or follow in detail the progress of a complete numerical experiment with OpenMOLE.

Introduction

We propose to present a study of an invasive plant growth model, the
Japanese knotweed. Our objective is to find effective management
strategies to control the evolution of this plant. The model was
developed in scala and its source code is available
here. Numerical experiments were performed with the OpenMOLE platform. The source code written to achieve them is detailed in this document.

Japanese knotweed stand

Japanese knotweed, Fallopia japonica, is an invasive plant. It is
native of Southeast Asia and was introduced in the British Isles during
the 19th century. It is present in North America, Northern and Central
Europe, Australia and New Zealand. Japanese knotweed grows in a wide
variety of sandy, swampy, rocky habitats. It mainly invades
human-modified habitats: roadsides, wasteland, as well as rivers banks.
It is perennial (lives several years) and geophyte: its rhizome allows
it to spend the winter season buried in the ground.

The length of its aerial stems varies between 1 and 3 meters. Its
rhizome can reach up to 8 cm in diameter, and can extend over 15 to 20
meters. The rhizome is also capable of burrowing deeply, from 2 to 3
meters, and represents 2/3 of the total biomass of the plant.

A commonly used management technique is mowing of aerial parts. Managers can vary its intensity and frequency.

The model describes the clonal growth of a knotweed, taking into account
mowing as a management strategy. It is a stochastic individual based
model, that allows individual variability in order to evaluate
management strategies as a function of the stage of proliferation
(deterministic models are often valid in large populations, that is,
when the invasive plant has already been present for a long time).

The individuals in the model are the crowns (place on the rhizome, on
the surface, from which sprout one or more stems). Each crown is
characterized by two traits (information): its position (in the plane)
and its biomass (the biomass of the underground rhizome, which connects
a crown to his or her daughter crowns). The set of crowns forms a
stand.

The simulation model takes place in two stages: first initialization,
then evolution. Both stages are based on the same plant dynamics but are
subject to different management methods. The initialization step is used
to generate a stand with a desired size, and the second to study the
effect of the given management strategy on the task obtained. We will
study the influence of the duration and frequency of mowing during this
second phase, on the area and the number of crowns in a knotweed stand.

Young sprout of Japanese knotweed

This document is organized into three parts. The first one presents the
simulation of the model, its input parameters and the output variables
being studied. The second one describes the model calibration, that is
the method for setting the model input parameters that controls the
plant growth dynamics. The last one meets the objective of the study by
addressing the influence of the management parameters on the plant.

Model Simulation

Inputs and outputs

The model has two types of parameters:

  • parameters concerning plant dynamics, for example the maximum
    biomass of a crown, or the intraspecific competition distance:
file: src/workflow.oms

val  K = Val[Double] // maximum biomass of a crown
val  L = Val[Double] // biomass growth speed 
val  distanceCompetition =  Val[Double] // intra-specific competition distance
val  distanceParent =  Val[Double] // distance around a crown that takes part in the birth rate
val  shape = Val[Double] // parameter of the Gamma law for the distance between a mother crown and its daughter
val  scale =  Val[Double] // parameter of the Gamma law for the distance between a mother crown and its daughter
val  deathParameterDecrease =  Val[Double] // decrease in the mortality rate as a function of biomass 
val  deathParameterScaling = Val[Double] // maximum mortality rate (for a null biomass)
val  mowingParameter = Val[Double] // effect of mowing
val  bbar = Val[Double] // birth rate in ideal condiations
val  a0 = Val[Double] // biomass of a crown at birth
  • management parameters for the evolution stage: the duration T of
    the management project, the average number of mowing event per year
    tau and the size of the initial population (number of crowns),
file: src/workflow.oms

val  T = Val[Double]
val  tau = Val[List[Double]]
val  proportionMowing = Val[List[Double]]
val  SeqInitialPopSize = Val[List[Double]] 
/* The model can evolve several stands at the same time,
 * and each can be applied different management settings.
 * This explains that these variables are of type
 * `List[Double]`: each list has one element for each stand.
 */

The model is designed to simulate the evolution of several knotweed
stands simultaneously and independently of each other. All stands share
the same dynamics, the same management for model initialization, and the
same duration of the management project. However, each stand can be
given a different proportion of mowing (proportionMowing) and average
number of annual mowing (tau).

The number of simulated years is controlled by the variable T. We may
want to stop it prematurely when it exceeds a certain population size,
and we must thus give as input of the model a population size not to be
exceeded in the simulation.

file: src/workflow.oms

val  SeqFinalSizeMax = Val[List[Int]]

The model is stochastic. Its execution depends on a seed to initialize
the random number generator.

file: src/workflow.oms

val seed = Val[Int]

The outputs of the models are the positions of the crowns of the initial
and final population.

Simulated stands. In black, the positions of the crowns of the initialpopulation (1200 crowns), in blue those of the final population (480crowns), after 8 years of mowing.

file: src/workflow.oms

/* The model writes the output about the initial and final populations 
 * of each stand in as many files as there are stands (remember that 
 * the model simulates the evolution of several stands in parallel) 
 * in two separate directories. 
 * We need variables to represent these directories.
 */
val dirSimOutputPopInit = Val[File]
val dirSimOutputPopFinal = Val[File]

Simulation

A simulation consists in executing the model with given parameter values
and the getting its output back. Here, we cover how a simulation is set
up and run.

file: src/workflow.oms

val model = ScalaTask("""
/* Before running the model, the directories in which the model will write 
 * the csv files that contain the positions of the crowns are created
 */
  val dirSimOutputPopInit = mkDir() 
  val dirSimOutputPopFinal = mkDir()

  /* The random number generator is initialized with a seed.
   */
  val rng = newRNG(seed)

The simulation is stopped prematurely in the following cases:

  • the iteration number needed for the initialization step exceeds a
    threshold which depends on the requested size,
  • the population size during the evolution stage exceeds a threshold
    that depends on initial or final reference population sizes
    (typically derived from field data),
  • the number of iterations of the model exceeds 4000000,
  • the desired population size isn’t reached at the end of
    initialization step after 4 trials.
file: src/workflow.oms

  val SeqNmaxPopIni = SeqInitialPopSize.map(s => math.round(2*s / 0.005).toInt).toSeq
  val taillePopFinaleMaxVect = ( SeqInitialPopSize zip SeqFinalSizeMax).map(t => math.max(10*t._1,10*t._2).toInt).toSeq
  val NmaxEvol = 4000000
  val compteurMax = 4

For the initialization step of the model, the proportion of mowed crown
at each mowing event is fixed at 0.9, the average number of mowing event
per year during the initialization is fixed at 1, and the first crown
created has a biomass equal to 3.

file: src/workflow.oms

  val proportionMowingInitial = 0.9
  val tauInitial = 1.0
  val biomassFirstIndiv = 3.0

Running the model in simulation takes several sets of parameters.

  • the initial population,
file: src/workflow.oms

  renouee.RunCalibration.evolutionWriteABC ( SeqInitialPopSize )(
  • the file paths where to write the simulation outputs,
file: src/workflow.oms

    dirSimOutputPopInit, dirSimOutputPopFinal )(
  • the stopping criteria for the simulation,
file: src/workflow.oms

    SeqNmaxPopIni, compteurMax ,NmaxEvol, SeqFinalSizeMax )(
  • the management parameters of respectively the initial and evolution
    stages ,
file: src/workflow.oms

    renouee.Management(
        tau = tauInitial, 
        proportionMowing = proportionMowingInitial ),
      renouee.ManagementSeveralEvolution(
        T = T,
        tau = tau, 
        proportionMowing = proportionMowing, 
        xAxisMowLimit = List.fill(tau.length)(0)),
        // Rq: the xAxisMowLimit is useless with the random mowing technique, 
	//but we have to set a value
  • the parameters of the plant dynamics,
file: src/workflow.oms

/* Note: In the instantiation of the PlantGrowth class, the names 
 * on the left refer to the parameters of the class, and on the right 
 * to the OpenMole variables defined earlier.
 */
      renouee.PlantGrowth(
        distanceCompetition = distanceCompetition, 
        distanceParent = distanceParent, 
        shape = shape, 
        scale = scale, 
        K = K, 
        L = L,
        mowingParameter = mowingParameter,
        deathParameterDecrease = deathParameterDecrease, 
        deathParameterScaling = deathParameterScaling,
        a0 = a0, 
        bbar = bbar,
        biomassFirstIndiv = biomassFirstIndiv),
      renouee.ManagementTechnique.Alea )(
  • and the random number generator..
file: src/workflow.oms

    rng )

The simulation is controled by variables left free which concern the
plant dynamics,

file: src/workflow.oms

""") set(
  inputs += (K, L, distanceCompetition, distanceParent, shape, scale, deathParameterDecrease, deathParameterScaling, mowingParameter, bbar, a0),

management,

file: src/workflow.oms

  inputs += (tau, proportionMowing, T),

the initialization seed of the random number generator,

file: src/workflow.oms

  inputs += (seed), 

and the initial and maximum desired population sizes.

file: src/workflow.oms

  inputs += ( SeqInitialPopSize, SeqFinalSizeMax), 

It produces as result the initial and final populations written in csv
files.

file: src/workflow.oms

/* Directories containing csv output simulation files 
 * for initial and final populations
 */
  outputs += (dirSimOutputPopInit, dirSimOutputPopFinal), 

file: src/workflow.oms

/* For technical reasons related to OpenMOLE, the variables used after the simulation
 * must be set as outputs of the task, even if they are not modified by the 
 * simulation and are not strictly speaking part of its results.
 */
  outputs += (K, L, distanceCompetition, distanceParent, shape, scale, deathParameterDecrease, deathParameterScaling, mowingParameter, bbar, a0),
  outputs += (tau, proportionMowing, T), 
  outputs += (seed), 
  outputs += ( SeqInitialPopSize, SeqFinalSizeMax),
  plugins += pluginsOf(renouee.RunCalibration)
)

Computation of the area

Let us recall that we want to study the impact of a management strategy
on the area occupied by a knotweed stand. As we have just seen, given
the parameters of plant dynamics and management, the model gives us the
positions and biomasses of the crowns of the initial and final
population. For each simulation, we need to calculate the area occupied
by the initial and final population and their sizes.

file: src/workflow.oms

val simu_areasIni = Val[Array[Double]]
val simu_areasFinale = Val[Array[Double]]
val simu_taillePopIni = Val[Array[Double]]
val simu_taillePopFinale = Val[Array[Double]]

The calculation of these quantities from the outputs of the model is
performed with the software R. The R GeoRange library provides a
function that calculates the convex envelope of a point cloud and its
area.

file: src/OMcalculAireEtPopSize.R

library(GeoRange)

 # We encapsulate these calculations in a function R which takes as argument
 # a name of the directory in which are located the output files of
 # simulation and a prefix of these files. 
 
CalculAireEtTaillePop = function(dirSimOutputPop,filePrefix){
  # We calculate the quantities of interest for each file 
  # in the simulation output directory with the given prefix 
 
  fileNames = list.files(dirSimOutputPop,pattern="*.txt")
  nbFiles=length(fileNames)

  # We initialize to 0 the vectors of quantities that interest us. 
   
  areas = numeric(nbFiles)
  popSize = numeric(nbFiles)

  # For each knotweed stand represented by a file, we calculate 
  # the area and population size.
    
  for (i in 1:nbFiles){
    pop=read.table(paste(dirSimOutputPop,"/",filePrefix,i-1,".txt",sep=""),
                   header=TRUE, sep= ",")
    x=pop$xPos
    y=pop$yPos
    
    if (length(x)>2){
      
      liste = chull(x,y)
      hull <- cbind(x,y) [liste,]
      
      areas[i]=CHullArea(hull[,1],hull[,2])
      popSize[i]=length(x)
      
    }
  }
  return(list(areas,popSize))
}

The calculation of the area and the population size is performed for all
the initial and final knotweed stands at the end of a simulation.

file: src/workflow.oms

val rTask =
  RTask("""
     source("OMcalculAireEtPopSize.R")
     resIni = CalculAireEtTaillePop("dirSimOutputPopInit","pos_popIni_")
     resFinal = CalculAireEtTaillePop("dirSimOutputPopFinal","pos_popFinal_")
     areasIni = resIni[[1]]
     popSizeIni = resIni[[2]]
     areasFinal = resFinal[[1]]
     popSizeFinal = resFinal[[2]]
    """,
  /* Running the R function requires installation of 
   * the libgdal and libproj libraries
   */
  install = Seq("apt update", "apt install -y libgdal-dev libproj-dev "), 
    libraries = Seq("GeoRange")
  ) set(
    resources += workDirectory / "OMcalculAireEtPopSize.R",
    
    // Directories where to find csv files in simulation output
    inputFiles += (dirSimOutputPopInit, "dirSimOutputPopInit"),
    inputFiles += (dirSimOutputPopFinal, "dirSimOutputPopFinal"),
    
    rOutputs += ("areasIni", simu_areasIni),
    rOutputs += ("areasFinal", simu_areasFinale), 
    rOutputs += ("popSizeIni", simu_taillePopIni),
    rOutputs += ("popSizeFinal", simu_taillePopFinale),
     
    /* For the same technical reasons as those mentioned above, 
     * it is necessary to declare in the inputs and outputs of the task 
     * all the variables that we want to use after it in the workflow.
     */
    inputs += (K, L, distanceCompetition, distanceParent, shape, scale, deathParameterDecrease, deathParameterScaling, mowingParameter, bbar, a0),
    inputs += (tau, proportionMowing, T, SeqInitialPopSize), 
    inputs += (seed),
    outputs += (K, L, distanceCompetition, distanceParent, shape, scale, deathParameterDecrease, deathParameterScaling, mowingParameter, bbar, a0),
    outputs += (tau, proportionMowing, T, SeqInitialPopSize), 
    outputs += (seed)
  )

The results of the calculations are saved in a csv file, one line for
each stand.

file: src/workflow.oms

def simResultToCSV(f: String) = CSVHook(workDirectory / f, seed,  tau, proportionMowing,  T, SeqInitialPopSize, simu_taillePopIni, simu_areasIni , simu_taillePopFinale, simu_areasFinale) set (
  csvHeader := "index, tau, proportionMowing,  T, initialPopSize, TailleIni, AireIni, TailleFinale, AireFinale" )

Example: A single simulation run

Running a simulation requires giving a value to its free input
variables. The following values are for illustrative purposes.

file: src/workflow.oms

val arguments = ScalaTask("""

  val seed = 3
  
  val proportionMowing = List(0.75, 1.0) 
  val T = 5.0
  val tau = List(4.0, 0.0)
  
  val SeqInitialPopSize = List(150.0, 500.0) 
  val SeqFinalSizeMax = List(1000, 1000) 

  val K = 5.0 
  val L = 0.33 
  val distanceCompetition  =  0.15
  val distanceParent =  0.15
  val shape = 10.0 
  val scale =  1.0 
  val deathParameterDecrease =  1.5
  val deathParameterScaling = 1.0 
  val mowingParameter = 0.6 
  val bbar = 1.0 
  val a0 = 0.5 
  
  """ )  set (
  outputs += (K, L, distanceCompetition, distanceParent, shape, scale, deathParameterDecrease, deathParameterScaling, mowingParameter, bbar, a0),
  outputs += (tau, proportionMowing, T), 
  outputs += (seed),
  outputs += (SeqInitialPopSize, SeqFinalSizeMax)
)

We can then transmit these values to the simulation, then calculate the
area of the initial and final knots, and write them in a csv file.

Workflow diagram for a modelsimulation

file: src/oneSimulation.oms

import _file_.workflow._
arguments -- model -- rTask hook(simResultToCSV("Resultat_Aire_Taille.csv"))

Each line of the csv file describes a simulated knotweed stand, giving
the management parameters that were applied to it and its initial and
final population sizes and areas.

index tau proportionMowing T initialPopSize TailleIni AireIni (m²) TailleFinale AireFinale
3 4.0 0.75 5.0 150.0 150.0 3.963 108.0 3.83
3 0.0 1.0 5.0 500.0 500.0 13.93 611.0 6.06

Model calibration

Our goal is to find the influence of management parameters on the stand
dynamics. We must therefore set the parameters of the plant dynamics.
For some, we do not find values in the literature on the Japanese
knotweed. We thus perform a calibration which consists in finding the
values of the parameters of the plant dynamics with which the model will
best reproduce field data. Comparisons between simulation outputs and
field data are made on the size and area of the initial and final
population of each stand (see below). For the calibration, the
management parameters can be read from field data.

Management parameters

We have information relating to the management of the stands in the
field. We can thus directly set model management parameters (T,tau,
proportionMowing, etc.) for each field stand without resorting to the
more complex calibration procedure.

file: src/workflow.oms

val argumentsCalibration = ScalaTask("""
  val SeqInitialPopSize = (workDirectory / "Martin_taillePop2008_allegee.txt").content.split(" ").map( _.trim.toDouble).toList
  val taillePop2015 = (workDirectory / "Martin_taillePop2015_allegee.txt").content.split(" ").map( _.trim.toDouble).toList 
  val proportionMowingSingleData = 0.9   			// on met la même proportion de fauche pour toutes les taches qui ne sont pas fauchées entièrement
  val proportionMowing = (workDirectory / "p_allegee.txt").content.split(" ").map( _.trim.toDouble).toList.map(t=> if(t==1)(1) else (proportionMowingSingleData)).toList
  val tau = (workDirectory / "tau_allegee.txt").content.split(" ").map( _.trim.toDouble).toList
  val T = 8.0
  val SeqFinalSizeMax =  (SeqInitialPopSize zip taillePop2015.toList).map(t => math.floor(math.max(10*t._1,10*t._2)).toInt).toSeq
""" )   set (
  resources += workDirectory / "data_allegee2" / "Martin_taillePop2008_allegee.txt",
  resources += workDirectory / "data_allegee2" / "Martin_taillePop2015_allegee.txt",
  resources += workDirectory / "data_allegee2" / "p_allegee.txt",
  resources += workDirectory / "data_allegee2" / "tau_allegee.txt",
  outputs += (K, L, distanceCompetition, distanceParent, shape, scale, deathParameterDecrease, deathParameterScaling, mowingParameter, bbar, a0),
  outputs += (seed),
  outputs += (tau, proportionMowing, T), 
  outputs += (SeqInitialPopSize, SeqFinalSizeMax),
  inputs += (K, L, distanceCompetition, distanceParent, shape, scale, deathParameterDecrease, deathParameterScaling, mowingParameter, bbar, a0),
  inputs += (seed)
)

Reading data

To compare the field data to the outputs of the model during
calibration, we will read the areas and population sizes of each
knotweed stand of the field.

file: src/workflow.oms

val data_taillePop2015 = Val[List[Double]]
val data_area2008 = Val[Array[Double]]
val data_area2015 = Val[Array[Double]] 
/* Note: we do not create an OpenMOLE variable 
 * `val data_taillePop2008 = Val [List [Double]]` 
 * because we already have a `SeqInitialPopSize` variable 
 * for the initial population size to make evolve in the model
 */

Here is an example of a dataset of 3 knotweed stands located in the
French Alps. (The result presented below, is realized with 18 stands
measured in 2008 and in 2015.)

Stand Size 2008 Area 2008 Size 2015 Area 2015 τ FullMow
1 261 14.525 112 10.815 0 0
2 1878 52.177 872 42.187 1 0
3 1063 75.899 1493 72.829 1 1

Data to be read are contained in text files.

file: src/workflow.oms

val data =
  ScalaTask("""
               val data_area2008 = (workDirectory / "Martin_area2008_allegee.txt").content.split(" ").map( _.trim.toDouble).toArray
               val data_area2015 = (workDirectory / "Martin_area2015_allegee.txt").content.split(" ").map( _.trim.toDouble).toArray
               val data_taillePop2015 = (workDirectory / "Martin_taillePop2015_allegee.txt").content.split(" ").map( _.trim.toDouble).toList
  """) set (
    resources += workDirectory / "data_allegee2" / "Martin_area2015_allegee.txt",
    resources += workDirectory / "data_allegee2" / "Martin_taillePop2015_allegee.txt",
    resources += workDirectory / "data_allegee2" / "Martin_area2008_allegee.txt",
    outputs += (data_area2008,data_area2015,data_taillePop2015) 
  )

Comparison of model outputs and field data

To find the parameter values with which the model best reproduces the
data, we compare the data to the outputs of a simulation by defining a
distance delta as the sum of the absolute differences between each
output variables of the model and its corresponding variable in the
data. Since the areas of the stands and the population sizes do not have
the same order of magnitude, the differences between simulated data and
field data are renormalized by the field data.

file: src/workflow.oms

val delta = Val[Double]

val objective = 
  ScalaTask("""
  val delta = 
      (SeqInitialPopSize zip simu_taillePopIni  ).map(z => ( math.abs( z._1 - z._2)/z._1 ).toDouble ).sum +
      (data_area2008 zip simu_areasIni  ).map(z => ( math.abs( z._1 - z._2)/z._1 ).toDouble ).sum +
      (data_area2015 zip simu_areasFinale  ).map(z => ( math.abs( z._1 - z._2)/z._1 ).toDouble ).sum +
      (data_taillePop2015 zip simu_taillePopFinale  ).map(z => ( math.abs( z._1 - z._2)/z._1 ).toDouble ).sum 
    """) set (
    inputs += (SeqInitialPopSize , simu_taillePopIni, data_area2008,simu_areasIni,data_area2015,simu_areasFinale,data_taillePop2015,simu_taillePopFinale),
    outputs += (delta)
    )

Calibration

Calibration consists of finding values for the plant dynamics parameters
that minimize the delta distance previously defined. We use the NSGA2
optimization algorithm which will look for values for free parameters in
predefined intervals. The intervals are chosen wide enough so that the
optimal values are not on the border.

file: src/workflow.oms

/* Technical note: To ensure the consistency of the workflow, we must
 * create a `Slot` on the `objective` task. The latter will only get executed
 * when the two tasks on which it depends are completed: 
 * `rTask` and `data`. Likewise, we need to create a common `evalFirst`
 * starting task for these both parallel tasks.
 */
val evalFirst = Strain(EmptyTask())
val objSlot = Slot(objective)  // we create a slot over the task objective

val nsga2 = 
  NSGA2Evolution(
    mu = 200,
    genome = Seq(
            K in (7.0 , 20.0) , 
            L in (0.1 , 1.0 )   ,
            distanceCompetition in (0.05 , 0.3),
            distanceParent  in (0.1 , 0.3) ,
            shape  in (1.0 , 10.0) ,
            scale  in (0.5 , 3.0)  ,
            deathParameterDecrease   in (0.2,3.0 ) ,
            deathParameterScaling  in (0.8, 1.5) , 
            mowingParameter in  (0.01 , 0.2 ),
            bbar  in (0.01 , 0.6 ) ,
            a0  in (0.2 , 4.0 ) 
            ) ,
    objectives = Seq(delta),
    stochastic = Stochastic(seed = seed, reevaluate = 0.2, replications = 100),
    evaluation = (evalFirst -- argumentsCalibration -- model -- rTask -- objSlot) & (evalFirst -- data -- objSlot),
    parallelism = 2000,
    termination = 300000
  )

The calibration algorithm is an iterative algorithm, which at each step
gives us a set of solutions. As the steps progress, the solutions
improve. The solutions are stored at each step in a csv file. Each
solution is a set of values for each free parameter.

file: src/workflow.oms

val savePopulationHook = SavePopulationHook(nsga2, workDirectory / "resultsCalibration")

Calibration execution

The calibration of the model required many simulations and mobilized the
resources of the European grid for one week.

file: src/calibrate.oms

import _file_.workflow._

// Requires the appropriate rights.
val env = EGIEnvironment("vo.complex-systems.eu") 
(nsga2 on env hook savePopulationHook) 

Workflow diagram for calibration

The results below correspond to the best solution (the one for which the
delta distance is the smallest) after 150 000 steps.

Variable Valeur Calibration
K 12.72
L 0.26
distanceCompetition 0.15
distanceParent 0.20
shape 4.34
scale 2.36
deathParameterDecrease 2.32
deathParameterScaling 1.12
mowingParameter 0.11
bbar 0.18
a0 1.73
delta 26.06
evolution.samples 79

Influence of the management parameters.

In the previous section, we were particularly interested in the
parameters of the plant dynamics. These parameters are now set to the
values found at the end of the previous section. We are now interested
in the influence of management parameters on the area of a simulated
knotweed stand, the objective of our study. In particular, the variables

  • tau : mean number of mowing event a year,
  • T: the duration of the management (year),

will be sampled according to a regular grid (cf below). The proportion
of crowns mown at each mowing event is set at 0.9, and the initial
population size at 1500. 50 replications are performed for each
parameter value by changing the seed of the random number generator. We
will study the average area on the replications of the same management
strategy.

Set the plant dynamics parameters

The parameters of the plant dynamics are fixed at the values previously
found by calibration.

file: src/workflow.oms

val argumentsExploration = ScalaTask("""

  // Parameters are read from the file "ParamMin.csv"
  val lines = (workDirectory / "ParamMin.csv").content.split("\n").toVector
  def doubleQuoteFilter(c: Char) = c != '"'
  val tempNames = lines.map(p => p.split(",").toList(0) )
  val tempVal = lines.map(p => p.split(",").toList(1).toDouble )

  // Plant parameters
  val K = tempVal(0)
  val L = tempVal(1)
  val distanceCompetition = tempVal(2)
  val distanceParent = tempVal(3)
  val shape = tempVal(4)
  val scale = tempVal(5)
  val deathParameterDecrease = tempVal(6)
  val deathParameterScaling = tempVal(7)
  val mowingParameter = tempVal(8)
  val bbar = tempVal(9)
  val a0 = tempVal(10)

  // Population size not to be exceeded
  val SeqFinalSizeMax = SeqInitialPopSize.map(p => math.max(p*10,20000).toInt) 

""" )  set (
    resources += workDirectory / "ParamMin.csv",
    outputs += (K, L, distanceCompetition, distanceParent, shape, scale, deathParameterDecrease, deathParameterScaling, mowingParameter, bbar, a0),
    outputs += (SeqFinalSizeMax),

    // Transmission of management parameters (controlled by the model's exploration task)
    inputs += (tau, proportionMowing, T), 
    outputs += (tau, proportionMowing, T), 
    inputs += (seed),
    outputs += (seed),
    inputs += (SeqInitialPopSize),
    outputs += (SeqInitialPopSize)
)

Sampling of management parameters

The average number of mowing per year, tau, varies from 0 to 15, the
management duration in years, T, varies from 0 to 20, each with a step
of 0.5. This results in a sample of 1271 pairs of values.

file: src/workflow.oms

def completeSampling =
                (tau in (0.0 to 15.0 by  0.5).map(p=>List(p))) x     
                (T in  (0.0 to 20.0 by 0.5))  x 
                (proportionMowing in  List( List(0.9))) x       
                (SeqInitialPopSize in List(List(1500.0)))

Running simulations

For each pair of values sampled for the management parameters, we run 50
simulations of the model by varying the seed of the random number
generator.

file: src/workflow.oms

// The simulations are to be run on the European grid.
val env = EGIEnvironment("vo.complex-systems.eu")

val exploration = DirectSampling(
  evaluation =  ( MoleTask(argumentsExploration -- model -- rTask) on env by 5 ) hook(simResultToCSV("resultatExploration.csv")) ,
  sampling = (completeSampling x (seed in (UniformDistribution[Int]()) take 50 )).shuffle
/* Mixing the points sampled by the call to the shuffle method 
 * avoids the longest simulations being all together, which can 
 * increase the total computation time when combining the simulations 
 * together for the parallel calculation.
 */
)

file: src/paramSampling.oms

import _file_.workflow._
exploration

Results

We have already seen that the results of each simulation are written in
a csv file, one line per knotweed stand, where each line gives the
management parameters used for the stand, its initial and final area.

Using these results, we draw regression curves for the average final
area (calculated with the 50 simulations for a same management
strategy), depending on the duration of the management project (T, on
the abscissa), for different τ (low in clear colours, up to 15 mows
per year in dark colours). The calculations are performed by the
resultCalibrate.R script.

Exploration: mean area as a funtion of T, for different values ofτ

In this example, we can see that with this model, we can not expect an
eradication of a knotweed stand of initial area 60m² in less than 11
years. The figure also tells us about the stand surface reduction (in
terms of final average area) when mowing once more per year. For
example, mowing 6 times a year instead of 5 for 10 years reduces the
surface of the stand by 4m² on average.

Conclusion

We presented a simulation study using the OpenMole platform to answer
questions in ecology. OpenMOLE helped us perform a calibration to infer
values for parameters for which we could find no data or information in
the literature. Then, it helped us use sampling methods and distributed
computing to simulate evolution of knotweed stands with different
management methods.

The result of this study is that the final area depends quadratically on
the duration of the management project. Thus depending on the cost of
management strategies in terms of frequency and duration, one can find
out which is most efficient according to this model. For example, one
can know whether it is more advantageous to extend the duration of
mowing for a certain period of time, or to increase the number of mows
per year.

Appendix: Reproducing the results

This section explains the practical steps to follow in order to
reproduce the results (tables and curves) presented in this document.

To reproduce the simulations, you must have downloaded the OpenMOLE software and R.

Start by cloning the git repository where the different files to be used
are located, this creates a renouee directory.

git clone https://gitlab.iscpif.fr/flavallee/renouee

Recall that it is possible to edit the source code files of this document with the command ./make.hs, in renouee (they are generated in the folder /src/reouee).

In OpenMOLE, you must import the plugin renouee_v7.jar (renouee/include, it notably contains theevolutionWriteABC function).

Then, upload all the .oms files and the .R file located in the
directory src (/renouee/src).

Then import the data files. Start importing the file ParamMin.csv
(“/renouee/code_R_Calibrate”).

Then create a new folder data_allegee2 in OpenMOLE, and import the 6
files .txt in the directory (/renouee/include/data_allegee2)

To perform a model simulation, (table of the section “Model
Simulation/Example: Running a simulation”), it is necessary to run the
oneSimulation.oms file in OpenMOLE.

It creates in your current directory a file Resultat_Aire_Taille.csv,
with the areas and sizes of the initial and final simulated poulations.

To obtain the calibration result (table in the section “Model
calibration/Calibration execution”), run the file calibrate.oms in
OpenMOLE. Note that this operation takes quite a long time (of the order
of the week using the EGI European grid).

If you do not have access to the computation grid, you can still start
the calculation by removing on env in the file calibrate.oms, line
(nsga2 on env hook savePopulationHook). It creates the firsts result
files, but it will probably not reach the convergence of the algorithm.
A directory resultsCalibration is created in the current directory. It
contains a file populationi.csv for each step i of the algorithm (i
is an integer).

We are only interested in the last file created, which gives several
sets of parameters for the plant dynamics (the lines of the file).

More precisely, we choose the line that has the best compromise between
the last two columns: delta (the quantity to be minimized) and
evolution$sample (the number of times the line has been tested).

Running the file resultCalibrate.R (“renouee/code_R_Calibrate”)
selects from a file in theresultsCalibration folder (in our case the
population167289.csv file) the line that minimizes delta and has been
tested at least 50 times, and writes it into the file ParamMin.csv.
The file produced by this script is used in OpenMOLE for the
exploration, it is generated on the command line with (to execute in
“/renouee/code_R_Calibrate”)

R -q --vanilla < resultCalibrate.R

To obtain the exploration result (curve in the section “Influence of the
management parameters./Exploitation of results”), run the file
paramSampling.oms in OpenMOLE.

Again, if you do not use the computation grid, remove “on env” in “val
Exploration”, at the end of the file workflow.oms. Running the file
paramSampling.oms creates the file resultatExploration.csv. Then,
running the file resultExploration.R (“renouee/code_R_Calibrate”) -
that require the file createMatrix.R, and the result file
resultatExploration.csv- enables to create the figure (mean area as a
funtion of T, for different values of τ).

R -q --vanilla < resultExploration.R
Use Case OpenMOLE : Effective management strategies for an invasive plant, the Japanese knotweed
Share this