kevinwolz / hisafer

An R toolbox for the Hi-sAFe biophysical agroforestry model
6 stars 4 forks source link

Create function to access crop varietal parameters through define_hisafe() #146

Open Guillaume-Blanchet opened 4 years ago

Guillaume-Blanchet commented 4 years ago

Hi Kevin, I dug more into the parametrization of the Pea plant file according to field data. You can have a look to the draft of a modelling report I am preparing here: https://guillaume-blanchet.github.io/hisafe-winter-pea-calibration-restinclieres/calibration-WinterPea.html https://guillaume-blanchet.github.io/hisafe-winter-pea-calibration-restinclieres/calibration-WinterPea-sim-outputs_visualization.html

I did not find any way to access to the varietal parameter through define_hisafe().

The way I worked so far was to prepare different "varieties" with a range of values for a given parameter. However, it is not very efficient as I had to prepare files by hand. Do you think it'd be possible to access to varietal parameters with a function similar to layer_params() ; root_init_params() ; tree_init_params() ?

I would like such function to perform optimization of parameters with an algorithm, like a Nelder-Mead Simplex algorithm, generally used by the STICS team to optimize varietal parameters.

Looking at the code of those functions, they seems pretty similar, so I guess it should not be to complicated, but you have probably a better overview on this. I can work on this, but if you had few tips for guiding me, it'd be appreciated !

Cheers, Guillaume

kevinwolz commented 4 years ago

@Guillaume-Blanchet this is definitely a need for hisafer. I think I was lazy in not including a function like this along with the other table definition functions because I just didn't have a need to modify the variety parameters like this.

I've created a new function, variety_params(), that works exactly the same as the other table definition functions that you mention. Install the latest version (1.4.24) from GitHub to get this improvement.

When using these functions, there are two things that are critical to remember:

(1) In order to define multiple simulations using these table functions, you must call the the table function multiple times. That is, if, for example, you want to define three simulations in an experiment with different values of croirac, this call does NOT get the desired effect:

hip <- define_hisafe(path = ".", template = "monocrop", varieties = variety_params(template = "monocrop", croirac = seq(0.05, 0.15, 0.05)))

Each call to variety_params() creates a single variety table. So, instead you need to do something like this:

VARIETIES <- list()
croirac.values <- seq(0.05, 0.15, 0.05)
for(i in croirac.values) VARIETIES <- c(VARIETIES, variety_params(template = "monocrop", croirac = i))

hip <- define_hisafe(path = ".", template = "monocrop", varieties = VARIETIES)
hip$exp.plan$varieties

Of course, rather than a for loop, I would recommend the use of purrr::map():

croirac.values <- seq(0.05, 0.15, 0.05)
VARIETIES <- purrr::map(croirac.values, function(x) variety_params(template = "monocrop", croirac = x)[[1]])

hip <- define_hisafe(path = ".", template = "monocrop", varieties = VARIETIES)
hip$exp.plan$varieties

(2) Any parameters passed to define_hisafe() that are located in the .TREE, .PLT, and .TEC files will be applied to all .TREE, .PLT, and .TEC files in each simulation. Consequently, there is really no good way right now to use define_hisafe(), even with the new variety_params(), to parameterize different variety values to different crops or in different years. This shouldn't impact your intended use, as I imagine you are running one-year simulations with winter pea and just modifying the .PLT file used in each simluation?

Guillaume-Blanchet commented 4 years ago

Wow ! Thanks for your reactivity Kevin ! :+1:

kevinwolz commented 4 years ago

Just let me know if you need help utilizing this within an optimization algorithm.

Guillaume-Blanchet commented 4 years ago

Hey @kevinwolz , just tried to use the new function. However :

  1. by adding the argument "varieties", there is a risk of confusion with the other argument "variete", formely used to select the crop variety
  2. when I use different values of a parameter as you suggested (croirac for instance), then the other values taken by default are those from durum wheat ("Allur" according to the new list created) and not the parameters of the crop defined in the .sim file.

Not sure about the best solution regarding this, but one idea would be to add one additional argument to explicitly name the plant file we want to modify.

kevinwolz commented 4 years ago

Please share your code.

Let me think about if/how renaming the table might work...

Guillaume-Blanchet commented 4 years ago

Yep. Here is the code I tried.

stlevamf.values <- seq(800, 860, 10)
VARIETIES <- purrr::map(stlevamf.values, function(x) variety_params(template = "./template-v1_winter_pea_monocrop_A2_2018-2019",
                                                                   stlevamf = x)[[1]])

hip_s2 <- define_hisafe(path = simulation_folder_path,
                        exp.name = "calibration_plt_s2",
                        profiles = c("cells", "plot"), # reduced output -> interest only for phenology
                        template = "template-v1_winter_pea_monocrop_A2_2018-2019",
                        weatherFile = paste0("./weather_files/restinclieres_A2_weather_1994-2019",c("","_EXC"),".wth"),
                        simulationDayStart = c(335),
                        layer.initialization = layer_init_params(template         = "./template-v1_winter_pea_monocrop_A2_2018-2019",
                                                                 no3Concentration = c(51.1,21.1, 5, 2, 0),
                                                                 nh4concentration = c(2.71,10.7,0,0,0),
                                                                 waterContent = c(0.185,0.280,0.308,0.335,0.354)),
                        #variete = c(3),
                        varieties = VARIETIES,
                        factorial = T,
                        densitesem = c(190),
                        profsem = c(2))

Template could changed to "restinclieres_monocrop_A2", and performed over a year with winter pea (2014-2015 or 2018-2019). Tell me if you want a full reprex folder.

kevinwolz commented 4 years ago

Yes, please share the folder. You can just zip it and attach it here.

Guillaume-Blanchet commented 4 years ago

Hi @kevinwolz , here is a folder. I ran this example simulation this morning. I confirm that the plant file generated with variety_params() is based on durum wheat only. reprex_folder.zip

kevinwolz commented 4 years ago

I took one look at your template folder and realized, in horror, what's going on. I will explain more later, but for right now to fix it all you need to do is remove all the .plt files from your template folder except the ONE that you want to use. The function will then pull from that one file.