Integrate relevant changes into PSE. Check pre initialisation vectors are already handled. Extend a couple of tests to convert from modellist to mapping. TODO : multiscale testing

Author: Samuel-amap

Project.toml

@@ -12,6 +12,7 @@ FLoops = "cc61a311-1640-44b5-9fba-1b764f453329"
 Markdown = "d6f4376e-aef5-505a-96c1-9c027394607a"
 MultiScaleTreeGraph = "dd4a991b-8a45-4075-bede-262ee62d5583"
 PlantMeteo = "4630fe09-e0fb-4da5-a846-781cb73437b6"
+SHA = "ea8e919c-243c-51af-8825-aaa63cd721ce"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 Tables = "bd369af6-aec1-5ad0-b16a-f7cc5008161c"
 Term = "22787eb5-b846-44ae-b979-8e399b8463ab"
@@ -25,6 +26,7 @@ FLoops = "0.2"
 Markdown = "1.9"
 MultiScaleTreeGraph = "0.13, 0.14"
 PlantMeteo = "0.6"
+SHA = "0.7.0"
 Statistics = "1.9"
 Tables = "1"
 Term = "1, 2"

src/PlantSimEngine.jl

@@ -22,6 +22,9 @@ import MultiScaleTreeGraph: symbol, node_id
 # To compute mean:
 import Statistics
 
+# For avoiding name conflicts when generating models from status vectors
+import SHA: sha1
+
 using PlantMeteo
 
 # UninitializedVar + PreviousTimeStep:
@@ -89,6 +92,9 @@ include("run.jl")
 # Fitting
 include("evaluation/fit.jl")
 
+# Utilities for mapping initialisation
+include("mtg/mapping/model_generation_from_status_vectors.jl")
+
 # Examples
 include("examples_import.jl")
 

src/mtg/initialisation.jl

@@ -307,6 +307,10 @@ automatically passed as is.
 """
 function init_simulation(mtg, mapping; nsteps=1, outputs=nothing, type_promotion=nothing, check=true, verbose=false)
 
+    # Ensure the user called the model generation function to handle vectors passed into a status
+    # before we keep going
+    check_statuses_contain_no_remaining_vectors(mapping)
+
     # Get the status of each node by node type, pre-initialised considering multi-scale variables:
     statuses, status_templates, reverse_multiscale_mapping, vars_need_init =
         init_statuses(mtg, mapping, first(hard_dependencies(mapping; verbose=false)); type_promotion=type_promotion, verbose=verbose, check=check)

src/mtg/mapping/model_generation_from_status_vectors.jl

@@ -0,0 +1,240 @@
+# Utilities to handle vectors passed into a mapping's statuses
+# This is a convenience when prototyping, not recommended for production or proper fitting
+# (Write the actual finalised model explicitely instead)
+
+
+# The way we generate models from status vectors is to eval() code at runtime.
+# A simple custom timestep model provides the correct index to the generated models
+# This approach feels a little brittle but works. A (possible ?) improvement would be to directly fiddle with the AST, but it's a little more involved
+# Another approach might be to generate a string to be included with include_string, that might avoid awkward global variables and world age problems
+
+# There will still be brittleness given that it's not trivial to handle user/modeler errors : 
+# For instance, providing a vector that is called in a scale mapping is likely to cause things to go badly
+
+# May need some more complex timestep models in the future
+# TODO : unhandled case : what if the timestep models are already in the provided modellist ?
+
+# These models might be worth exposing in the future ?
+PlantSimEngine.@process "basic_current_timestep" verbose = false
+
+struct HelperCurrentTimestepModel <: AbstractBasic_Current_TimestepModel
+end
+
+PlantSimEngine.inputs_(::HelperCurrentTimestepModel) = (next_timestep=1,)
+PlantSimEngine.outputs_(m::HelperCurrentTimestepModel) = (current_timestep=1,)
+
+function PlantSimEngine.run!(m::HelperCurrentTimestepModel, models, status, meteo, constants=nothing, extra=nothing)
+    status.current_timestep = status.next_timestep
+ end 
+
+ PlantSimEngine.ObjectDependencyTrait(::Type{<:HelperCurrentTimestepModel}) = PlantSimEngine.IsObjectDependent()
+ PlantSimEngine.TimeStepDependencyTrait(::Type{<:HelperCurrentTimestepModel}) = PlantSimEngine.IsTimeStepDependent()
+
+ PlantSimEngine.@process "basic_next_timestep" verbose = false
+ struct HelperNextTimestepModel <: AbstractBasic_Next_TimestepModel
+ end
+ 
+ PlantSimEngine.inputs_(::HelperNextTimestepModel) = (current_timestep=1,)
+ PlantSimEngine.outputs_(m::HelperNextTimestepModel) = (next_timestep=1,)
+ 
+ function PlantSimEngine.run!(m::HelperNextTimestepModel, models, status, meteo, constants=nothing, extra=nothing)
+     status.next_timestep = status.current_timestep + 1
+  end 
+
+PlantSimEngine.ObjectDependencyTrait(::Type{<:HelperNextTimestepModel}) = PlantSimEngine.IsObjectDependent()
+PlantSimEngine.TimeStepDependencyTrait(::Type{<:HelperNextTimestepModel}) = PlantSimEngine.IsTimeStepDependent()
+
+
+# TODO should the new_status be copied ?
+# Note : User specifies at which level they want the basic timestep model to be inserted at, as well as the meteo length
+function replace_mapping_status_vectors_with_generated_models(mapping_with_vectors_in_status, timestep_model_organ_level, nsteps)
+    
+    mapping = Dict(organ => models for (organ, models) in mapping_with_vectors_in_status)
+    for (organ,models) in mapping
+        for status in models
+            if isa(status, Status)
+                new_status, generated_models = generate_model_from_status_vector_variable(mapping, timestep_scale, status, organ, nsteps)
+
+                if length(generated_models) > 0
+                    if organ == timestep_model_organ_level
+                        mapping[organ] = (
+                            models..., 
+                            generated_models..., 
+                            HelperNextTimestepModel(),
+                            MultiScaleModel(
+                            model=HelperCurrentTimestepModel(),
+                            mapping=[PreviousTimeStep(:next_timestep),],
+                            ),
+                            new_status,)
+                    else
+                        mapping[organ] = (
+                            models..., 
+                            generated_models...,                          
+                            new_status,)
+                    end
+                end
+            end
+        end
+    end
+
+    return mapping
+end
+
+# Note : eval works in global scope, and state synchronisation doesn't occur until one returns to top-level
+# This is to enable optimisations. See 'world-age problem'. The doc for eval currently isn't detailed enough.
+# Essentially, generating a struct with a process_ method and then immediately creating a simulation graph
+# that calls process_ will fail as it won't yet be defined since state hasn't synchronised. 
+# Returning a new mapping to top-level and *then* creating the graph will work.
+# The fact that eval works in global scope is also why we make use of some global variables here
+function generate_model_from_status_vector_variable(mapping, timestep_scale, status, organ, nsteps)
+    
+    # Note : 534f1c161f91bb346feba1a84a55e8251f5ad446 is a prefix to reduce likelihood of global variable name conflicts
+    # it is the hash generated by bytes2hex(sha1("PlantSimEngine_prototype"))
+    # If this function is hard to read, copy it into a temporary file and remove the hash suffix
+
+    # Ah, another point that remains to be seen is that those CSV.SentinelArrays.ChainedVector obtained from the meteo file isn't an AbstractVector
+    # meaning currently we won't generate models from them unless the conversion is made before that
+    # So another minor potential improvement would be to return a warning to the user and do the conversion when generating the model
+    # See the test code in test-mapping.jl : cumsum(meteo_day.TT) returns such a data structure
+
+    global generated_models_534f1c161f91bb346feba1a84a55e8251f5ad446 = ()
+    global new_status_534f1c161f91bb346feba1a84a55e8251f5ad446 = NamedTuple()
+   
+    for symbol in keys(status)
+        global value_534f1c161f91bb346feba1a84a55e8251f5ad446 = getproperty(status, symbol)
+        if isa(value_534f1c161f91bb346feba1a84a55e8251f5ad446, AbstractVector)
+            @assert length(value_534f1c161f91bb346feba1a84a55e8251f5ad446) > 0 "Error during generation of models from vector values provided at the $organ-level status : provided $symbol vector is empty"
+            # TODO : Might need to fiddle with timesteps here in the future in case of varying timestep models
+            @assert nsteps == length(value_534f1c161f91bb346feba1a84a55e8251f5ad446) "Error during generation of models from vector values provided at the $organ-level status : provided $symbol vector length doesn't match the expected # of timesteps"
+            var_type = eltype(value_534f1c161f91bb346feba1a84a55e8251f5ad446)
+            base_name = string(symbol) * bytes2hex(sha1(join(value_534f1c161f91bb346feba1a84a55e8251f5ad446)))
+            process_name = lowercase(base_name)
+  
+            var_titlecase::String = titlecase(base_name)
+            model_name = "My$(var_titlecase)Model"
+            process_abstract_name = "Abstract$(var_titlecase)Model"
+            var_vector = "$(symbol)_vector"
+
+            abstract_process_decl = "abstract type $process_abstract_name <: PlantSimEngine.AbstractModel end"
+            eval(Meta.parse(abstract_process_decl))
+            
+            process_name_decl = "PlantSimEngine.process_(::Type{$process_abstract_name}) = :$process_name"
+            eval(Meta.parse(process_name_decl))
+    
+            struct_decl::String = "struct $model_name <: $process_abstract_name \n$var_vector::Vector{$var_type} \nend\n\n"
+            eval(Meta.parse(struct_decl))
+            
+            inputs_decl::String = "function PlantSimEngine.inputs_(::$model_name)\n(current_timestep=1,)\nend\n\n"
+            eval(Meta.parse(inputs_decl))
+    
+            default_value = value_534f1c161f91bb346feba1a84a55e8251f5ad446[1]
+            outputs_decl::String = "function PlantSimEngine.outputs_(::$model_name)\n($symbol=$default_value,)\nend\n\n"
+            eval(Meta.parse(outputs_decl))
+    
+            constructor_decl =  "$model_name(; $var_vector = Vector{$var_type}()) = $model_name($var_vector)\n\n"
+            eval(Meta.parse(constructor_decl))
+    
+            run_decl = "function PlantSimEngine.run!(m::$model_name, models, status, meteo, constants=nothing, extra_args=nothing)\nstatus.$symbol = m.$var_vector[status.current_timestep]\nend\n\n"
+            eval(Meta.parse(run_decl))
+    
+            # add name to vector of models
+            if timestep_scale != organ
+                mapping_decl = "mapping[\"($organ)\"] = MultiScaleModel($process_name($var_vector), mapping=\"($timestep_scale)\" => (:current_timestep,))"
+                eval(Meta.parse(mapping_decl))
+            else
+            end 
+       
+        model_add_decl = "generated_models_534f1c161f91bb346feba1a84a55e8251f5ad446 = (generated_models_534f1c161f91bb346feba1a84a55e8251f5ad446..., $model_name($var_vector=$value_534f1c161f91bb346feba1a84a55e8251f5ad446),)"
+        eval(Meta.parse(model_add_decl))
+        else
+            new_status_decl = "new_status_534f1c161f91bb346feba1a84a55e8251f5ad446 = Status(; NamedTuple(new_status_534f1c161f91bb346feba1a84a55e8251f5ad446)..., $symbol=$value_534f1c161f91bb346feba1a84a55e8251f5ad446)"
+            eval(Meta.parse(new_status_decl))
+        end
+    end
+    
+    @assert length(status) == length(new_status_534f1c161f91bb346feba1a84a55e8251f5ad446) + length(generated_models_534f1c161f91bb346feba1a84a55e8251f5ad446) "Error during generation of models from vector values provided at the $organ-level status"
+    return new_status_534f1c161f91bb346feba1a84a55e8251f5ad446, generated_models_534f1c161f91bb346feba1a84a55e8251f5ad446
+end
+
+
+# This is a helper function only for testing purposes, but it makes sense to include it here since it calls 
+# generate_model_from_status_vector_variable, which has those awkward global variables
+function modellist_to_mapping(modellist_original::ModelList, modellist_status; nsteps=nothing, outputs=nothing)
+    
+    modellist = Base.copy(modellist_original, modellist_original.status)
+
+    default_scale = "Default"
+    mtg = MultiScaleTreeGraph.Node(MultiScaleTreeGraph.NodeMTG("/", default_scale, 0, 0),)
+
+    models = modellist.models
+
+    mapping_incomplete = Dict(
+        default_scale => (
+        models..., 
+        MultiScaleModel(
+        model=HelperCurrentTimestepModel(),
+        mapping=[PreviousTimeStep(:next_timestep),],
+        ),
+        Status((modellist_status..., current_timestep=1,next_timestep=1,))
+        ),
+    )
+    
+    timestep_scale = "Default"
+    organ = "Default"
+ 
+    # todo improve on this
+    st = (last(mapping_incomplete["Default"]))
+    new_status, generated_models = generate_model_from_status_vector_variable(mapping_incomplete, timestep_scale, st, organ, nsteps)
+
+    mapping = Dict(default_scale => (
+        models..., generated_models..., 
+        HelperNextTimestepModel(),
+        MultiScaleModel(
+            model=HelperCurrentTimestepModel(),
+            mapping=[PreviousTimeStep(:next_timestep),],
+            ),
+            new_status,
+    ),
+    )
+
+    if isnothing(outputs)
+        f = []
+        for i in 1:length(modellist.models)
+            aa = init_variables(modellist.models[i])
+            bb = keys(aa)
+            for j in 1:length(bb)
+                push!(f, bb[j])
+            end
+            #f = (f..., bb...)
+        end
+
+        f = unique!(f)
+        all_vars = (f...,)
+        #all_vars = merge((keys(init_variables(object.models[i])) for i in 1:length(object.models))...)
+    else 
+        all_vars = outputs
+        # TODO sanity check
+    end
+
+    return mtg, mapping, Dict(default_scale => all_vars)
+end
+
+function check_statuses_contain_no_remaining_vectors(mapping)
+    for (organ,models) in mapping
+
+        # Special case (scales that map to a single-model don't need to be declared as a tuple for user-convenience)
+        if isa(models, AbstractModel) || isa(models, MultiScaleModel)
+            continue
+        end
+
+        for status in models
+            if isa(status, Status)
+                for symbol in keys(status)
+                    value = getproperty(status, symbol)
+                    @assert !isa(value, AbstractVector) "Error : Mapping status at $organ level contains a vector. If this was intentional, call the function generate_models_from_status_vectors on your mapping before calling run!. And bear in mind this is not meant for production. If this wasn't intentional, then it's likely an issue on the mapping definition, or an unusual model."
+                end
+            end
+        end
+    end
+    return true
+end

test/helper-functions.jl

@@ -0,0 +1,45 @@
+# Simple helper functions that can be used in various tests here and there
+
+
+function compare_outputs_modellist_mapping(models, graphsim)
+    graphsim_df = outputs(graphsim, DataFrame)
+
+    graphsim_df_outputs_only = select(graphsim_df, Not([:timestep, :organ, :node]))
+    models_df = DataFrame(status(models))
+    
+    models_df_sorted = models_df[:, sortperm(names(models_df))]
+    graphsim_df_outputs_only_sorted = graphsim_df_outputs_only[:, sortperm(names(graphsim_df_outputs_only))]
+    return graphsim_df_outputs_only_sorted == models_df_sorted
+end
+
+# doesn't check for mtg equality
+function compare_outputs_graphsim(graphsim, graphsim2)
+    graphsim_df = outputs(graphsim, DataFrame)
+    graphsim_df_sorted = graphsim_df[:, sortperm(names(graphsim_df))]
+    
+    graphsim2_df = outputs(graphsim2, DataFrame)
+    graphsim2_df_sorted = graphsim2_df[:, sortperm(names(graphsim2_df))]
+    return graphsim_df_sorted == graphsim2_df_sorted
+end
+
+# Breaking this function into two to ensure eval() state synchronisation happens (see comments around the modellist_to_mapping definition)
+# Naming could be better
+function check_multiscale_simulation_is_equivalent_begin(models::ModelList, status, meteo)
+
+    mtg, mapping, out = PlantSimEngine.modellist_to_mapping(models, status; nsteps=length(meteo), outputs=nothing)
+    return mtg, mapping, out
+end
+
+function check_multiscale_simulation_is_equivalent_end(models::ModelList, mtg, mapping, out, meteo)
+    graph_sim = PlantSimEngine.GraphSimulation(mtg, mapping, nsteps=length(meteo), check=true, outputs=out)
+
+    sim = run!(graph_sim,
+        meteo,
+        PlantMeteo.Constants(),
+        nothing;
+        check=true,
+        executor=SequentialEx()
+    );
+
+    return compare_outputs_modellist_mapping(models, graph_sim)
+end

test/runtests.jl

@@ -7,6 +7,8 @@ using MultiScaleTreeGraph
 using PlantMeteo, Statistics
 using Documenter # for doctests
 
+include("helper-functions.jl")
+
 @testset "Testing PlantSimEngine" begin
     Aqua.test_all(PlantSimEngine, ambiguities=false)
     Aqua.test_ambiguities([PlantSimEngine])