diff --git a/pipeline/Project.toml b/pipeline/Project.toml index 7163245c4..e34070f2b 100644 --- a/pipeline/Project.toml +++ b/pipeline/Project.toml @@ -20,7 +20,6 @@ JLD2 = "033835bb-8acc-5ee8-8aae-3f567f8a3819" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" LogExpFunctions = "2ab3a3ac-af41-5b50-aa03-7779005ae688" MCMCChains = "c7f686f2-ff18-58e9-bc7b-31028e88f75d" -Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" RCall = "6f49c342-dc21-5d91-9882-a32aef131414" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" ReverseDiff = "37e2e3b7-166d-5795-8a7a-e32c996b4267" diff --git a/pipeline/data/truth_sims/placeholder.md b/pipeline/data/truth_data/placeholder.md similarity index 100% rename from pipeline/data/truth_sims/placeholder.md rename to pipeline/data/truth_data/placeholder.md diff --git a/pipeline/scripts/run_testmode_pipeline.jl b/pipeline/scripts/run_analysis_pipeline.jl similarity index 100% rename from pipeline/scripts/run_testmode_pipeline.jl rename to pipeline/scripts/run_analysis_pipeline.jl diff --git a/pipeline/scripts/run_pipeline.jl b/pipeline/scripts/run_pipeline.jl deleted file mode 100644 index 2c59c5756..000000000 --- a/pipeline/scripts/run_pipeline.jl +++ /dev/null @@ -1,25 +0,0 @@ -# Local environment script to run the analysis pipeline -using Pkg -Pkg.activate(joinpath(@__DIR__(), "..")) -using Dagger - -@info(""" - Running the analysis pipeline. - -------------------------------------------- - """) - -# Define the backend resources to use for the pipeline -# in this case we are using distributed local workers with loaded modules -using Distributed -pids = addprocs(; exeflags = ["--project=$(Base.active_project())"]) - -@everywhere using EpiAwarePipeline - -# Create an instance of the pipeline behaviour -pipeline = RtwithoutRenewalPipeline() - -# Run the pipeline -do_pipeline(pipeline) - -# Remove the workers -rmprocs(pids) diff --git a/pipeline/scripts/run_priorpred_pipeline.jl b/pipeline/scripts/run_priorpred_pipeline.jl index 504dbe245..d22091f7a 100644 --- a/pipeline/scripts/run_priorpred_pipeline.jl +++ b/pipeline/scripts/run_priorpred_pipeline.jl @@ -3,8 +3,11 @@ using Pkg Pkg.activate(joinpath(@__DIR__(), "..")) using Dagger +@assert !isempty(ARGS) "Test mode script requires the number of draws as an argument." +ndraws = parse(Int64, ARGS[1]) + @info(""" - Running the analysis pipeline. + Running the prior predictive pipeline in test mode with $(ndraws) draws per model. -------------------------------------------- """) @@ -15,11 +18,17 @@ pids = addprocs(; exeflags = ["--project=$(Base.active_project())"]) @everywhere using EpiAwarePipeline -# Create an instance of the pipeline behaviour -pipeline = RtwithoutRenewalPriorPipeline() +# Create instances of the pipeline behaviour + +pipelines = [ + SmoothOutbreakPipeline(ndraws = ndraws, nchains = 1, priorpredictive = true), + MeasuresOutbreakPipeline(ndraws = ndraws, nchains = 1, priorpredictive = true), + SmoothEndemicPipeline(ndraws = ndraws, nchains = 1, priorpredictive = true), + RoughEndemicPipeline(ndraws = ndraws, nchains = 1, priorpredictive = true) +] # Run the pipeline -do_pipeline(pipeline) +do_pipeline(pipelines) # Remove the workers rmprocs(pids) diff --git a/pipeline/src/EpiAwarePipeline.jl b/pipeline/src/EpiAwarePipeline.jl index d5349606b..501987799 100644 --- a/pipeline/src/EpiAwarePipeline.jl +++ b/pipeline/src/EpiAwarePipeline.jl @@ -11,9 +11,11 @@ with execution determined by available computational resources. module EpiAwarePipeline using CSV, Dagger, DataFramesMeta, Dates, Distributions, DocStringExtensions, DrWatson, - EpiAware, Plots, Statistics, ADTypes, AbstractMCMC, Plots, JLD2, MCMCChains, Turing, - DynamicPPL, LogExpFunctions, RCall, LinearAlgebra, Random, AlgebraOfGraphics, - CairoMakie, ReverseDiff + EpiAware, Statistics, ADTypes, AbstractMCMC, JLD2, MCMCChains, Turing, DynamicPPL, + LogExpFunctions, RCall, LinearAlgebra, Random, AlgebraOfGraphics, CairoMakie, + ReverseDiff + +using EpiAware.EpiInfModels: oneexpy # Exported pipeline types export AbstractEpiAwarePipeline, EpiAwarePipeline, AbstractRtwithoutRenewalPipeline, @@ -56,7 +58,7 @@ export make_prediction_dataframe_from_output, make_truthdata_dataframe, export figureone, figuretwo # Exported functions: plot functions -export plot_truth_data, plot_Rt +export plot_truth_data, plot_Rt, prior_predictive_plot include("docstrings.jl") include("pipeline/pipeline.jl") @@ -67,6 +69,5 @@ include("infer/infer.jl") include("forecast/forecast.jl") include("scoring/scoring.jl") include("analysis/analysis.jl") -include("mainplots/mainplots.jl") -include("plot_functions.jl") +include("plotting/plotting.jl") end diff --git a/pipeline/src/constructors/make_inference_method.jl b/pipeline/src/constructors/make_inference_method.jl index 548f3d664..497f87be0 100644 --- a/pipeline/src/constructors/make_inference_method.jl +++ b/pipeline/src/constructors/make_inference_method.jl @@ -8,7 +8,7 @@ Constructs an inference method for the given pipeline. This is a default method. - An inference method. """ -function make_inference_method(pipeline::AbstractEpiAwarePipeline; ndraws::Integer = 2000, +function make_inference_method(ndraws::Integer, pipeline::AbstractEpiAwarePipeline; mcmc_ensemble::AbstractMCMC.AbstractMCMCEnsemble = MCMCSerial(), nruns_pthf::Integer = 4, maxiters_pthf::Integer = 100, nchains::Integer = 4) return EpiMethod( @@ -19,27 +19,28 @@ function make_inference_method(pipeline::AbstractEpiAwarePipeline; ndraws::Integ end """ -Method for sampling from prior predictive distribution of the model. +Example pipeline. """ -function make_inference_method(pipeline::RtwithoutRenewalPriorPipeline; n_samples = 2_000) +function make_inference_method( + pipeline::EpiAwareExamplePipeline; ndraws::Integer = 20, + mcmc_ensemble::AbstractMCMC.AbstractMCMCEnsemble = MCMCThreads(), + nruns_pthf::Integer = 4, maxiters_pthf::Integer = 100, nchains::Integer = 4) return EpiMethod( - pre_sampler_steps = AbstractEpiOptMethod[], - sampler = DirectSample(n_samples = n_samples) + pre_sampler_steps = [ManyPathfinder(nruns = nruns_pthf, maxiters = maxiters_pthf)], + sampler = NUTSampler( + target_acceptance = 0.9, adtype = AutoReverseDiff(; compile = true), + ndraws = ndraws, nchains = nchains, mcmc_parallel = mcmc_ensemble) ) end """ -Pipeline test mode method for sampling from prior predictive distribution of the model. +Method for sampling from prior predictive distribution of the model. """ function make_inference_method( - pipeline::EpiAwareExamplePipeline; ndraws::Integer = 20, - mcmc_ensemble::AbstractMCMC.AbstractMCMCEnsemble = MCMCThreads(), - nruns_pthf::Integer = 4, maxiters_pthf::Integer = 100, nchains::Integer = 4) + pipeline::AbstractRtwithoutRenewalPipeline, ::Val{:priorpredictive}) return EpiMethod( - pre_sampler_steps = [ManyPathfinder(nruns = nruns_pthf, maxiters = maxiters_pthf)], - sampler = NUTSampler( - target_acceptance = 0.9, adtype = AutoReverseDiff(; compile = true), ndraws = ndraws, - nchains = nchains, mcmc_parallel = mcmc_ensemble) + pre_sampler_steps = AbstractEpiOptMethod[], + sampler = DirectSample(n_samples = pipeline.ndraws) ) end @@ -55,7 +56,12 @@ Constructs an inference method for the Rt-without-renewal pipeline. # Examples """ function make_inference_method(pipeline::AbstractRtwithoutRenewalPipeline) - return make_inference_method(pipeline; ndraws = pipeline.ndraws, - mcmc_ensemble = pipeline.mcmc_ensemble, nruns_pthf = pipeline.nruns_pthf, - maxiters_pthf = pipeline.maxiters_pthf, nchains = pipeline.nchains) + if pipeline.priorpredictive + return make_inference_method(pipeline, Val(:priorpredictive)) + else + return make_inference_method( + pipeline.ndraws, pipeline; mcmc_ensemble = pipeline.mcmc_ensemble, + nruns_pthf = pipeline.nruns_pthf, + maxiters_pthf = pipeline.maxiters_pthf, nchains = pipeline.nchains) + end end diff --git a/pipeline/src/constructors/selector.jl b/pipeline/src/constructors/selector.jl index 65cf4d7b8..edcbd7f9b 100644 --- a/pipeline/src/constructors/selector.jl +++ b/pipeline/src/constructors/selector.jl @@ -13,3 +13,11 @@ Example/test mode is to return a randomly selected item from the list. function _selector(list, pipeline::EpiAwareExamplePipeline) return [rand(list)] end + +""" +Internal method for selecting from a list of items based on the pipeline type. +Example/test mode is to return a randomly selected item from the list. +""" +function _selector(list, pipeline::AbstractRtwithoutRenewalPipeline) + return pipeline.testmode ? [rand(list)] : list +end diff --git a/pipeline/src/infer/InferenceConfig.jl b/pipeline/src/infer/InferenceConfig.jl index 61cdd3b38..a7b53f4dd 100644 --- a/pipeline/src/infer/InferenceConfig.jl +++ b/pipeline/src/infer/InferenceConfig.jl @@ -12,13 +12,18 @@ Inference configuration struct for specifying the parameters and models used in - `epimethod::E`: Inference method. - `transformation::F`: Transformation function. - `log_I0_prior::Distribution`: Prior for log initial infections. Default is `Normal(log(100.0), 1e-5)`. +- `lookahead::X`: Number of days to forecast ahead. +- `latent_model_name::String`: Name of the latent model. +- `pipeline`: Pipeline type defining the inference scenario. # Constructors - `InferenceConfig(igp, latent_model, observation_model; gi_mean, gi_std, case_data, tspan, epimethod, transformation = exp)`: Constructs an `InferenceConfig` object with the specified parameters. - `InferenceConfig(inference_config::Dict; case_data, tspan, epimethod)`: Constructs an `InferenceConfig` object from a dictionary of configuration values. """ -struct InferenceConfig{T, F, IGP, L, O, E, D <: Distribution, X <: Integer} +struct InferenceConfig{ + T, F, IGP, L, O, E, D <: Distribution, X <: Integer, + P <: AbstractRtwithoutRenewalPipeline} gi_mean::T gi_std::T igp::IGP @@ -32,19 +37,23 @@ struct InferenceConfig{T, F, IGP, L, O, E, D <: Distribution, X <: Integer} transformation::F log_I0_prior::D lookahead::X + latent_model_name::String + pipeline::P - function InferenceConfig(igp, latent_model, observation_model; gi_mean, gi_std, - case_data, truth_I_t, truth_I0, tspan, epimethod, - transformation = exp, log_I0_prior, lookahead) - new{typeof(gi_mean), typeof(transformation), - typeof(igp), typeof(latent_model), typeof(observation_model), - typeof(epimethod), typeof(log_I0_prior), typeof(lookahead)}( + function InferenceConfig( + igp, latent_model, observation_model; gi_mean, gi_std, case_data, + truth_I_t, truth_I0, tspan, epimethod, transformation = oneexpy, + log_I0_prior, lookahead, latent_model_name, pipeline) + new{typeof(gi_mean), typeof(transformation), typeof(igp), + typeof(latent_model), typeof(observation_model), typeof(epimethod), + typeof(log_I0_prior), typeof(lookahead), typeof(pipeline)}( gi_mean, gi_std, igp, latent_model, observation_model, - case_data, truth_I_t, truth_I0, tspan, epimethod, transformation, log_I0_prior, lookahead) + case_data, truth_I_t, truth_I0, tspan, epimethod, + transformation, log_I0_prior, lookahead, latent_model_name, pipeline) end function InferenceConfig( - inference_config::Dict; case_data, truth_I_t, truth_I0, tspan, epimethod) + inference_config::Dict; case_data, truth_I_t, truth_I0, tspan, epimethod, pipeline) InferenceConfig( inference_config["igp"], inference_config["latent_namemodels"].second, @@ -57,11 +66,71 @@ struct InferenceConfig{T, F, IGP, L, O, E, D <: Distribution, X <: Integer} tspan = tspan, epimethod = epimethod, log_I0_prior = inference_config["log_I0_prior"], - lookahead = inference_config["lookahead"] + lookahead = inference_config["lookahead"], + latent_model_name = inference_config["latent_namemodels"].first, + pipeline ) end end +""" +Internal function that returns a dictionary containing key configuration fields from the given `InferenceConfig` object. +The dictionary includes the following keys: + +- `"igp"`: The string representation of the `igp` field. +- `"latent_model"`: The name of the latent model. +- `"gi_mean"`: The mean of the generation interval. +- `"gi_std"`: The standard deviation of the generation interval. +- `"tspan"`: The time span for the inference. +- `"priorpredictive"`: The prior predictive setting. + +# Arguments +- `config::InferenceConfig`: An instance of `InferenceConfig` containing the configuration details. + +# Returns +- `Dict{String, Any}`: A dictionary with the key configuration fields. +""" +function _saved_config_fields(config::InferenceConfig) + return Dict( + "igp" => string(config.igp), + "latent_model" => config.latent_model_name, + "gi_mean" => config.gi_mean, + "gi_std" => config.gi_std, + "tspan" => string(config.tspan[1]) * "_" * string(config.tspan[2]), + "priorpredictive" => config.pipeline.priorpredictive, + "scenario" => config.pipeline |> typeof |> string + ) +end + +""" +This function makes inference on the underlying parameters of the model specified +in the `InferenceConfig` object `config`. + +# Arguments +- `config::InferenceConfig`: The configuration object containing the case data +to make inference on and model configuration. +- `epiprob::EpiProblem`: The EpiProblem object containing the model to make inference on. + +# Returns +- `inference_results`: The results of the simulation or inference. + +""" +function create_inference_results(config, epiprob) + #Return the sampled infections and observations + idxs = config.tspan[1]:config.tspan[2] + y_t = ismissing(config.case_data) ? missing : + Vector{Union{Missing, Int64}}(config.case_data[idxs]) + inference_results = apply_method(epiprob, + config.epimethod, + (y_t = y_t,) + ) + inference_results = apply_method(epiprob, + config.epimethod, + (y_t = y_t,); + ) + return inference_results +end + """ This method makes inference on the underlying parameters of the model specified in the `InferenceConfig` object `config`. @@ -77,22 +146,46 @@ to make inference on and model configuration. function infer(config::InferenceConfig) #Define the EpiProblem epiprob = define_epiprob(config) - idxs = config.tspan[1]:config.tspan[2] + + #Define savable configuration + save_config = _saved_config_fields(config) #Return the sampled infections and observations - y_t = ismissing(config.case_data) ? missing : config.case_data[idxs] - inference_results = apply_method(epiprob, - config.epimethod, - (y_t = y_t,); - ) + inference_results = try + create_inference_results(config, epiprob) + catch e + string(e) + end - forecast_results = generate_forecasts( - inference_results.samples, inference_results.data, epiprob, config.lookahead) + if config.pipeline.priorpredictive + if inference_results isa String + return Dict("priorpredictive" => inference_results, + "inference_config" => save_config) + else + fig = prior_predictive_plot( + config, inference_results, epiprob; ps = [0.025, 0.1, 0.25]) + figdir = config.pipeline.testmode ? mktempdir() : plotsdir("priorpredictive") + figpath = joinpath(figdir, "priorpred_" * savename(save_config) * ".png") + CairoMakie.save(figpath, fig) + return Dict("priorpredictive" => "Pass", "inference_config" => save_config) + end + else + forecast_results = try + generate_forecasts( + inference_results.samples, inference_results.data, epiprob, config.lookahead) + catch e + string(e) + end - epidata = epiprob.epi_model.data - score_results = summarise_crps(config, inference_results, forecast_results, epidata) + epidata = epiprob.epi_model.data + score_results = try + summarise_crps(config, inference_results, forecast_results, epidata) + catch e + string(e) + end - return Dict("inference_results" => inference_results, - "epiprob" => epiprob, "inference_config" => config, - "forecast_results" => forecast_results, "score_results" => score_results) + return Dict( + "inference_results" => inference_results, "inference_config" => save_config, + "forecast_results" => forecast_results, "score_results" => score_results) + end end diff --git a/pipeline/src/infer/define_epiprob.jl b/pipeline/src/infer/define_epiprob.jl index 047afd9b8..4290556db 100644 --- a/pipeline/src/infer/define_epiprob.jl +++ b/pipeline/src/infer/define_epiprob.jl @@ -16,7 +16,7 @@ function define_epiprob(config::InferenceConfig) model_data = EpiData( gen_distribution = gen_distribution, transformation = config.transformation) #Build the epidemiological model - epi = config.igp(model_data, config.log_I0_prior) + epi = config.igp(data = model_data, initialisation_prior = config.log_I0_prior) epi_prob = EpiProblem(epi_model = epi, latent_model = config.latent_model, diff --git a/pipeline/src/infer/generate_inference_results.jl b/pipeline/src/infer/generate_inference_results.jl index 32199be0b..03c63d61b 100644 --- a/pipeline/src/infer/generate_inference_results.jl +++ b/pipeline/src/infer/generate_inference_results.jl @@ -14,19 +14,20 @@ Generate inference results based on the given configuration of inference model o - `inference_results`: The generated inference results. """ function generate_inference_results( - truthdata, inference_config, pipeline::AbstractEpiAwarePipeline; - inference_method, datadir_name = "epiaware_observables") + truthdata, inference_config, pipeline::AbstractEpiAwarePipeline) tspan = make_tspan( pipeline; T = inference_config["T"], lookback = inference_config["lookback"]) + inference_method = make_inference_method(pipeline) config = InferenceConfig( inference_config; case_data = truthdata["y_t"], truth_I_t = truthdata["I_t"], - truth_I0 = truthdata["truth_I0"], tspan, epimethod = inference_method) + truth_I0 = truthdata["truth_I0"], tspan, epimethod = inference_method, pipeline = pipeline) # produce or load inference results prfx = _inference_prefix(truthdata, inference_config, pipeline) + _datadir_str = _get_inferencedatadir_str(pipeline) inference_results, inferencefile = produce_or_load( - infer, config, datadir(datadir_name); prefix = prfx) + infer, config, _datadir_str; prefix = prfx) return inference_results end @@ -49,12 +50,13 @@ which is deleted after the function call. - `inference_results`: The generated inference results. """ function generate_inference_results( - truthdata, inference_config, pipeline::EpiAwareExamplePipeline; inference_method) + truthdata, inference_config, pipeline::EpiAwareExamplePipeline) tspan = make_tspan( pipeline; T = inference_config["T"], lookback = inference_config["lookback"]) + inference_method = make_inference_method(pipeline) config = InferenceConfig( inference_config; case_data = truthdata["y_t"], truth_I_t = truthdata["I_t"], - truth_I0 = truthdata["truth_I0"], tspan = tspan, epimethod = inference_method) + truth_I0 = truthdata["truth_I0"], tspan = tspan, epimethod = inference_method, pipeline = pipeline) # produce or load inference results prfx = _inference_prefix(truthdata, inference_config, pipeline) @@ -65,23 +67,3 @@ function generate_inference_results( infer, config, datadir_name; prefix = prfx) return inference_results end - -""" -Method for prior predictive modelling. -""" -function generate_inference_results( - inference_config, pipeline::RtwithoutRenewalPriorPipeline) - tspan = make_tspan( - pipeline; T = inference_config["T"], lookback = inference_config["lookback"]) - config = InferenceConfig( - inference_config; case_data = missing, tspan, epimethod = DirectSample()) - - # produce or load inference results - prfx = _inference_prefix(truthdata, inference_config, pipeline) - - datadir_name = mktempdir() - - inference_results, inferencefile = produce_or_load( - infer, config, datadir(datadir_name); prefix = prfx) - return inference_results -end diff --git a/pipeline/src/infer/infer.jl b/pipeline/src/infer/infer.jl index 3054f01ce..56724123a 100644 --- a/pipeline/src/infer/infer.jl +++ b/pipeline/src/infer/infer.jl @@ -1,5 +1,5 @@ include("InferenceConfig.jl") -include("inference_prefix.jl") +include("inference_helpers.jl") include("generate_inference_results.jl") include("map_inference_results.jl") include("define_epiprob.jl") diff --git a/pipeline/src/infer/inference_prefix.jl b/pipeline/src/infer/inference_helpers.jl similarity index 79% rename from pipeline/src/infer/inference_prefix.jl rename to pipeline/src/infer/inference_helpers.jl index 4cf735b6a..057afbaa3 100644 --- a/pipeline/src/infer/inference_prefix.jl +++ b/pipeline/src/infer/inference_helpers.jl @@ -16,7 +16,11 @@ This is an internal method that generates the part of the prefix for the inferen results file name from the pipeline. """ _prefix_from_pipeline(pipeline::AbstractEpiAwarePipeline) = "observables" -_prefix_from_pipeline(pipeline::AbstractRtwithoutRenewalPipeline) = pipeline.prefix +function _prefix_from_pipeline(pipeline::AbstractRtwithoutRenewalPipeline) + pipeline.priorpredictive ? + "priorpredictive_" * "" * pipeline.prefix : + "inference_" * pipeline.prefix +end """ This is an internal method that generates the prefix for the inference results file name. @@ -53,3 +57,12 @@ function _inference_prefix( return _prefix_from_pipeline(pipeline) * _prefix_from_config(truthdata, inference_config) end + +""" +Internal method for setting the data directory path for the inference data. +""" +_get_inferencedatadir_str(pipeline::AbstractEpiAwarePipeline) = datadir("epiaware_observables") +function _get_inferencedatadir_str(pipeline::AbstractRtwithoutRenewalPipeline) + pipeline.testmode ? mktempdir() : + pipeline.priorpredictive ? datadir("priorpredictive") : datadir("epiaware_observables") +end diff --git a/pipeline/src/infer/map_inference_results.jl b/pipeline/src/infer/map_inference_results.jl index c33b95bee..36d2c21e1 100644 --- a/pipeline/src/infer/map_inference_results.jl +++ b/pipeline/src/infer/map_inference_results.jl @@ -15,9 +15,9 @@ tasks from `Dagger.@spawn`. """ function map_inference_results( - truthdata, inference_configs, pipeline::AbstractEpiAwarePipeline; inference_method) + truthdata, inference_configs, pipeline::AbstractEpiAwarePipeline) map(inference_configs) do inference_config Dagger.@spawn generate_inference_results( - truthdata, inference_config, pipeline; inference_method) + truthdata, inference_config, pipeline) end end diff --git a/pipeline/src/pipeline/do_inference.jl b/pipeline/src/pipeline/do_inference.jl index f9cfccb34..846b671f5 100644 --- a/pipeline/src/pipeline/do_inference.jl +++ b/pipeline/src/pipeline/do_inference.jl @@ -11,8 +11,7 @@ An array of inference results. """ function do_inference(truthdata, pipeline::AbstractEpiAwarePipeline) inference_configs = make_inference_configs(pipeline) - inference_method = make_inference_method(pipeline) inference_results = map_inference_results( - truthdata, inference_configs, pipeline; inference_method) + truthdata, inference_configs, pipeline) return inference_results end diff --git a/pipeline/src/pipeline/do_truthdata.jl b/pipeline/src/pipeline/do_truthdata.jl index 3781bc3f5..e3db49afe 100644 --- a/pipeline/src/pipeline/do_truthdata.jl +++ b/pipeline/src/pipeline/do_truthdata.jl @@ -16,3 +16,32 @@ function do_truthdata(pipeline::AbstractEpiAwarePipeline) end return truthdata_from_configs end + +""" +Generate truth data for the given pipeline. + +# Arguments +- `pipeline::AbstractRtwithoutRenewalPipeline`: The pipeline object for which to generate truth data. + +# Returns +- If `pipeline.priorpredictive` is `true`, returns a list containing a dictionary with missing data and initial values. +- Otherwise, returns a list of spawned tasks that generate truth data based on configurations. + +# Details +- When `pipeline.priorpredictive` is `true`, the function returns a dictionary with keys `"y_t"`, `"I_t"`, `"truth_I0"`, and `"truth_gi_mean"`, where `"y_t"` is set to `missing`, `"I_t"` is a vector of 100 elements all set to `1.0`, and both `"truth_I0"` and `"truth_gi_mean"` are set to `1.0`. +- When `pipeline.priorpredictive` is `false`, the function generates truth data configurations using `make_truth_data_configs(pipeline)` and spawns tasks to generate truth data for each configuration using `Dagger.@spawn`. +""" +function do_truthdata(pipeline::AbstractRtwithoutRenewalPipeline) + if pipeline.priorpredictive + missingdata = Dict("y_t" => missing, "I_t" => fill(1.0, 100), "truth_I0" => 1.0, + "truth_gi_mean" => 1.0) + return [missingdata] + else + truth_data_configs = make_truth_data_configs(pipeline) + truthdata_from_configs = map(truth_data_configs) do truth_data_config + return Dagger.@spawn cache=true generate_truthdata( + truth_data_config, pipeline; plot = false) + end + return truthdata_from_configs + end +end diff --git a/pipeline/src/pipeline/pipelinetypes.jl b/pipeline/src/pipeline/pipelinetypes.jl index a2146ff42..3dc7d017a 100644 --- a/pipeline/src/pipeline/pipelinetypes.jl +++ b/pipeline/src/pipeline/pipelinetypes.jl @@ -44,6 +44,8 @@ Rt = make_Rt(pipeline) |> Rt -> plot(Rt, maxiters_pthf::Integer = 100 nchains::Integer = 4 prefix::String = "smooth_outbreak" + testmode::Bool = false + priorpredictive::Bool = false end """ @@ -57,6 +59,8 @@ The pipeline type for the Rt pipeline for an outbreak scenario where Rt has maxiters_pthf::Integer = 100 nchains::Integer = 4 prefix::String = "measures_outbreak" + testmode::Bool = false + priorpredictive::Bool = false end """ @@ -70,6 +74,8 @@ The pipeline type for the Rt pipeline for an endemic scenario where Rt changes i maxiters_pthf::Integer = 100 nchains::Integer = 4 prefix::String = "smooth_endemic" + testmode::Bool = false + priorpredictive::Bool = false end """ @@ -83,4 +89,6 @@ The pipeline type for the Rt pipeline for an endemic scenario where Rt changes i maxiters_pthf::Integer = 100 nchains::Integer = 4 prefix::String = "rough_endemic" + testmode::Bool = false + priorpredictive::Bool = false end diff --git a/pipeline/src/plot_functions.jl b/pipeline/src/plot_functions.jl deleted file mode 100644 index aa8b28600..000000000 --- a/pipeline/src/plot_functions.jl +++ /dev/null @@ -1,50 +0,0 @@ -""" -Plot the true cases and latent infections. This is the default method for plotting. - -# Arguments -- `data`: A dictionary containing the data for plotting. -- `config`: The configuration for the truth data scenario. -- `pipeline::AbstractEpiAwarePipeline`: The pipeline object which sets pipeline - behavior. - -# Returns -- `plt_cases`: The plot object representing the cases and latent infections. -""" -function plot_truth_data( - data, config, pipeline::AbstractEpiAwarePipeline; plotsname = "truth_data") - plt_cases = Plots.scatter( - data["y_t"], label = "Cases", xlabel = "Time", ylabel = "Daily cases", - title = "Cases and latent infections", legend = :bottomright) - Plots.plot!(plt_cases, data["I_t"], label = "True latent infections") - - if !isdir(plotsdir(plotsname)) - mkdir(plotsdir(plotsname)) - end - _plotsname = _simulate_prefix(pipeline) * plotsname - savefig(plt_cases, plotsdir(plotsname, savename(_plotsname, config, "png"))) - return plt_cases -end - -""" -Plot and save the plot of the true Rt values over time. - -# Arguments -- `true_Rt`: An array of true Rt values. -- `pipeline::AbstractEpiAwarePipeline`: The pipeline object which sets pipeline - behavior. - -# Returns -- `plt_Rt`: The plot object. - -""" -function plot_Rt(true_Rt, pipeline::AbstractEpiAwarePipeline) - plt_Rt = plot(true_Rt, label = "True Rt", xlabel = "Time", ylabel = "Rt", - title = "True Rt", legend = :topright) - - if !isdir(plotsdir("truth_data")) - mkdir(plotsdir("truth_data")) - end - savefig(plt_Rt, plotsdir("truth_data", "true_Rt")) - - return plt_Rt -end diff --git a/pipeline/src/plotting/basicplots.jl b/pipeline/src/plotting/basicplots.jl new file mode 100644 index 000000000..142ace5cc --- /dev/null +++ b/pipeline/src/plotting/basicplots.jl @@ -0,0 +1,74 @@ +""" +Internal method for creating a path/name for a plot. This also creates the directory +if it does not exist. +""" +function _mkplotpath(pipeline::AbstractEpiAwarePipeline, config, plotsubdir) + !isdir(plotsdir(plotsubdir)) ? mkdir(plotsdir(plotsubdir)) : nothing + _plotsname = _simulate_prefix(pipeline) + return plotsdir(plotsubdir, savename(_plotsname, config, "png")) +end + +""" +Plot the true cases and latent infections. This is the default method for plotting. + +# Arguments +- `data`: A dictionary containing the data for plotting. +- `config`: The configuration for the truth data scenario. +- `pipeline::AbstractEpiAwarePipeline`: The pipeline object which sets pipeline + behavior. +- `plotsubdir`: The subdirectory to save the plot. +- `saveplot`: A boolean indicating whether to save the plot. + +# Returns +- `f`: The plot object representing the cases and latent infections. +- `plotpath`: The path of the plot. +""" +function plot_truth_data(data, config, pipeline::AbstractEpiAwarePipeline; + plotsubdir = "truth_data", saveplot = true, kwargs...) + f = Figure(; backgroundcolor = :white) + ax = Axis(f[1, 1], + title = "Cases and latent infections", + xlabel = "Time", + ylabel = "Daily cases", + kwargs... + ) + scatter!(ax, data["y_t"], color = :black, label = "Daily cases") + lines!(ax, data["I_t"], label = "True latent infections", color = :red, + linestyle = :dash) + axislegend(position = :rt, framevisible = false, backgroundcolor = (:white, 0.0)) + + plotpath = _mkplotpath(pipeline, config, plotsubdir) + saveplot && CairoMakie.save(plotpath, f) + return f, plotpath +end + +""" +Plot and save the plot of the true Rt values over time. + +# Arguments +- `true_Rt`: An array of true Rt values. +- `config`: A scenario configuration. +- `pipeline::AbstractEpiAwarePipeline`: The pipeline object which sets pipeline + behavior. +- `plotsubdir`: The subdirectory to save the plot. +- `saveplot`: A boolean indicating whether to save the plot. + +# Returns +- `f`: The plot object representing the cases and latent infections. +- `plotpath`: The path of the plot. + +""" +function plot_Rt(true_Rt, config, pipeline::AbstractEpiAwarePipeline; + plotsubdir = "truth_data", saveplot = true) + f = Figure(; backgroundcolor = :white) + ax = Axis(f[1, 1], + title = "True Reproduction number", + xlabel = "Time", + ylabel = "Rt" + ) + lines!(ax, true_Rt) + + plotpath = EpiAwarePipeline._mkplotpath(pipeline, config, plotsubdir) + saveplot && CairoMakie.save(plotpath, f) + return f, plotpath +end diff --git a/pipeline/src/mainplots/df_checking.jl b/pipeline/src/plotting/df_checking.jl similarity index 100% rename from pipeline/src/mainplots/df_checking.jl rename to pipeline/src/plotting/df_checking.jl diff --git a/pipeline/src/mainplots/figureone.jl b/pipeline/src/plotting/figureone.jl similarity index 100% rename from pipeline/src/mainplots/figureone.jl rename to pipeline/src/plotting/figureone.jl diff --git a/pipeline/src/mainplots/figuretwo.jl b/pipeline/src/plotting/figuretwo.jl similarity index 100% rename from pipeline/src/mainplots/figuretwo.jl rename to pipeline/src/plotting/figuretwo.jl diff --git a/pipeline/src/mainplots/mainplots.jl b/pipeline/src/plotting/plotting.jl similarity index 54% rename from pipeline/src/mainplots/mainplots.jl rename to pipeline/src/plotting/plotting.jl index ddcc65ba1..5e594171f 100644 --- a/pipeline/src/mainplots/mainplots.jl +++ b/pipeline/src/plotting/plotting.jl @@ -1,3 +1,5 @@ +include("basicplots.jl") include("df_checking.jl") include("figureone.jl") include("figuretwo.jl") +include("prior_predictive_plot.jl") diff --git a/pipeline/src/plotting/prior_predictive_plot.jl b/pipeline/src/plotting/prior_predictive_plot.jl new file mode 100644 index 000000000..c1427381c --- /dev/null +++ b/pipeline/src/plotting/prior_predictive_plot.jl @@ -0,0 +1,101 @@ +""" +Internal method for generating a title string for a prior predictive plot based on the provided + configuration. + +# Arguments +- `config::InferenceConfig`: `InferenceConfig` object containing the configuration for the + prior predictive plot. +# Returns +- `String`: A formatted title string for the prior predictive plot. + +""" +function _make_prior_plot_title(config::InferenceConfig) + igp_str = string(config.igp) + latent_model_str = config.latent_model_name |> uppercase + gi_mean_str = config.gi_mean |> string + T_str = config.tspan[2] |> string + return "Prior pred. IGP: $(igp_str), latent model: $(latent_model_str), truth gi mean: $(gi_mean_str), T: $(T_str)" +end + +""" +Generate a prior predictive plot for the given configuration, output, and epidemiological probabilities. + +# Arguments +- `config`: Configuration object containing settings for the plot. +- `output`: Output object containing generated data for plotting. +- `epiprob`: Epidemiological probabilities object. +- `ps`: Array of percentiles for quantile calculations (default: [0.05, 0.1, 0.25]). +- `bottom_alpha`: Opacity for the lowest percentile band (default: 0.1). +- `top_alpha`: Opacity for the highest percentile band (default: 0.5). +- `case_color`: Color for the cases plot (default: :black). +- `logI_color`: Color for the log(Incidence) plot (default: :purple). +- `rt_color`: Color for the exponential growth rate plot (default: :blue). +- `Rt_color`: Color for the reproduction number plot (default: :green). +- `figsize`: Tuple specifying the size of the figure (default: (750, 600)). + +# Returns +- `fig`: A Figure object containing the prior predictive plots. + +# Notes +- The function asserts that all percentiles in `ps` are in the range [0, 0.5). +- The function creates a 2x2 grid of subplots with linked x-axes for the top and bottom rows. +- The function plots the median and percentile bands for cases, log(Incidence), exponential growth rate, and reproduction number. +""" +function prior_predictive_plot(config, output, epiprob; + ps = [0.05, 0.1, 0.25], + bottom_alpha = 0.1, + top_alpha = 0.5, + case_color = :black, + logI_color = :purple, + rt_color = :blue, + Rt_color = :green, + figsize = (750, 600)) + @assert all(0 .<= ps .< 0.5) "Percentiles must be in the range [0, 0.5)" + prior_pred_plot_title = _make_prior_plot_title(config) + qs, n_levels = _setup_levels(sort(ps)) + opacity_scale = range(bottom_alpha, top_alpha, length = n_levels) |> collect + + # Create the figure and axes + fig = Figure(size = figsize) + ax11 = Axis(fig[1, 1]; xlabel = "t", ylabel = "Cases") + ax12 = Axis(fig[1, 2]; xlabel = "t", ylabel = "log(Incidence)") + ax21 = Axis(fig[2, 1]; xlabel = "t", ylabel = "Exp. growth rate") + ax22 = Axis(fig[2, 2]; xlabel = "t", ylabel = "Reproduction number") + linkxaxes!(ax11, ax21) + linkxaxes!(ax12, ax22) + Label(fig[0, :]; text = prior_pred_plot_title, fontsize = 16) + + # Quantile calculations + gen_y_t = mapreduce(hcat, output.generated) do gen + gen.generated_y_t + end |> X -> timeseries_samples_into_quantiles(X, qs) + gen_quantities = generate_quantiles_for_targets(output, epiprob.epi_model.data, qs) + + # Plot the prior predictive samples + # Cases + f = findfirst(!ismissing, gen_y_t[:, 1]) + lines!(ax11, 1:size(gen_y_t, 1), gen_y_t[:, 1], + color = case_color, linewidth = 3, label = "Median") + for i in 1:n_levels + band!(ax11, f:size(gen_y_t, 1), gen_y_t[f:size(gen_y_t, 1), (2 * i)], + gen_y_t[f:size(gen_y_t, 1), (2 * i) + 1], + color = (case_color, opacity_scale[i]), + label = "($(ps[i]*100)-$((1 - ps[i])*100))%") + end + vlines!(ax11, [f], color = case_color, linestyle = :dash, label = "Obs. window") + axislegend(ax11; position = :lt, framevisible = false) + + # Other quantities + for (ax, target, c) in zip( + [ax12, ax21, ax22], [gen_quantities.log_I_t, gen_quantities.rt, gen_quantities.Rt], + [logI_color, rt_color, Rt_color]) + lines!(ax, 1:size(target, 1), target[:, 1], + color = logI_color, linewidth = 3, label = "Median") + for i in 1:n_levels + band!(ax, 1:size(target, 1), target[:, (2 * i)], target[:, (2 * i) + 1], + color = (c, opacity_scale[i]), label = "") + end + end + + fig +end diff --git a/pipeline/src/simulate/TruthSimulationConfig.jl b/pipeline/src/simulate/TruthSimulationConfig.jl index 09dc2203e..da0f0dc45 100644 --- a/pipeline/src/simulate/TruthSimulationConfig.jl +++ b/pipeline/src/simulate/TruthSimulationConfig.jl @@ -46,7 +46,7 @@ function simulate(config::TruthSimulationConfig) transformation = config.transformation) #Sample infections NB: prior is arbitrary because I0 is fixed. - epi = config.igp(model_data, Normal(log(config.I0), 0.01)) + epi = config.igp(model_data; initialisation_prior = Normal(log(config.I0), 0.01)) inf_mdl = generate_latent_infs(epi, log.(config.truth_process)) |> mdl -> fix(mdl, (init_incidence = log(config.I0),)) @@ -67,7 +67,8 @@ function simulate(config::TruthSimulationConfig) (std = 1.0, cluster_factor = config.cluster_factor, ϵ_t = config.logit_daily_ascertainment)) - y_t = obs_model() .|> y -> ismissing(y) ? missing : Int(y) + y_t = obs_model() .|> (y -> ismissing(y) ? missing : Int64(y)) |> + Vector{Union{Missing, Int64}} #Return the sampled infections and observations diff --git a/pipeline/src/simulate/generate_truthdata.jl b/pipeline/src/simulate/generate_truthdata.jl index e27666bf0..373d5cb3d 100644 --- a/pipeline/src/simulate/generate_truthdata.jl +++ b/pipeline/src/simulate/generate_truthdata.jl @@ -20,8 +20,7 @@ the `simulate` function to generate the truth data. This is the default method. - `truthfile`: The file path where the truth data is saved. """ function generate_truthdata( - truth_data_config, pipeline::AbstractEpiAwarePipeline; plot = true, - datadir_str = "truth_data") + truth_data_config, pipeline::AbstractEpiAwarePipeline; plot = true) default_params = make_default_params(pipeline) config = TruthSimulationConfig( truth_process = default_params["Rt"], gi_mean = truth_data_config["gi_mean"], @@ -29,9 +28,10 @@ function generate_truthdata( cluster_factor = default_params["cluster_factor"], I0 = default_params["I0"]) prefix = _simulate_prefix(pipeline) + _datadir_str = _get_truthdatadir_str(pipeline) truthdata, truthfile = produce_or_load( - simulate, config, datadir(datadir_str); prefix = prefix) + simulate, config, _datadir_str; prefix = prefix) if plot plot_truth_data(truthdata, config, pipeline) end diff --git a/pipeline/src/simulate/simulate.jl b/pipeline/src/simulate/simulate.jl index 7d81fa1ba..f7fbccd77 100644 --- a/pipeline/src/simulate/simulate.jl +++ b/pipeline/src/simulate/simulate.jl @@ -1,3 +1,3 @@ include("TruthSimulationConfig.jl") -include("simulate_prefix.jl") +include("simulate_helpers.jl") include("generate_truthdata.jl") diff --git a/pipeline/src/simulate/simulate_helpers.jl b/pipeline/src/simulate/simulate_helpers.jl new file mode 100644 index 000000000..59fedfbb7 --- /dev/null +++ b/pipeline/src/simulate/simulate_helpers.jl @@ -0,0 +1,15 @@ +""" +Internal method for setting the prefix for the truth data file name. +""" +_simulate_prefix(pipeline::AbstractEpiAwarePipeline) = "truth_data" +function _simulate_prefix(pipeline::AbstractRtwithoutRenewalPipeline) + "truth_data_" * pipeline.prefix +end + +""" +Internal method for setting the data directory path for the truth data. +""" +_get_truthdatadir_str(pipeline::AbstractEpiAwarePipeline) = datadir("truth_data") +function _get_truthdatadir_str(pipeline::AbstractRtwithoutRenewalPipeline) + pipeline.testmode ? mktempdir() : datadir("truth_data") +end diff --git a/pipeline/src/simulate/simulate_prefix.jl b/pipeline/src/simulate/simulate_prefix.jl deleted file mode 100644 index 15d0d6d99..000000000 --- a/pipeline/src/simulate/simulate_prefix.jl +++ /dev/null @@ -1,7 +0,0 @@ -""" -Internal method for setting the prefix for the truth data file name. -""" -_simulate_prefix(pipeline::AbstractEpiAwarePipeline) = "truth_data" -function _simulate_prefix(pipeline::AbstractRtwithoutRenewalPipeline) - "truth_data_" * pipeline.prefix -end diff --git a/pipeline/src/utils/timeseries_utils.jl b/pipeline/src/utils/timeseries_utils.jl index a24e85628..ce94e49c2 100644 --- a/pipeline/src/utils/timeseries_utils.jl +++ b/pipeline/src/utils/timeseries_utils.jl @@ -1,4 +1,3 @@ - """ Transforms a matrix of time series samples into quantiles. @@ -14,8 +13,12 @@ A matrix where each row represents the quantiles of a time series. """ function timeseries_samples_into_quantiles(X, qs) mapreduce(vcat, eachrow(X)) do row - _row = filter(x -> !isnan(x), row) - quantile(_row, qs)' + if any(ismissing, row) + return fill(missing, length(qs))' + else + _row = filter(x -> !isnan(x), row) + return quantile(_row, qs)' + end end end @@ -39,10 +42,29 @@ An array of quantiles for each target. """ function generate_quantiles_for_targets(output, D::EpiData, qs) - mapreduce(_process_reduction, output["forecast_results"].generated, - output["inference_results"].samples[:init_incidence]) do gen, logI0 - calculate_processes(gen.I_t, exp(logI0), D, EpiAwareExamplePipeline()) + mapreduce(_process_reduction, output.generated, + output.samples[:init_incidence]) do gen, logI0 + calculate_processes(gen.I_t, exp(logI0), D) end |> res -> map(res) do X timeseries_samples_into_quantiles(X, qs) end end + +""" +Internal method that given a list of percentiles `ps`, this function calculates the corresponding quantile levels. +This function returns both the `p/2` and `1 - p/2` quantiles for each percentile `p` in `ps`. +# Arguments +- `ps::Vector{Float64}`: A vector of percentiles. + +# Returns +- `qs::Vector{Float64}`: A vector of quantile levels, including the median (0.5) and the calculated quantiles based on the input percentiles. +- `n_levels::Int`: The number of input percentiles. + +""" +function _setup_levels(ps) + n_levels = length(ps) + qs = mapreduce(vcat, ps) do percentile + [percentile / 2, 1 - percentile / 2] + end |> x -> [0.5; x] + return qs, n_levels +end diff --git a/pipeline/test/constructors/test_constructors.jl b/pipeline/test/constructors/test_constructors.jl index eba8600c1..c608729bd 100644 --- a/pipeline/test/constructors/test_constructors.jl +++ b/pipeline/test/constructors/test_constructors.jl @@ -1,5 +1,4 @@ @testset "make_gi_params: returns a dictionary with correct keys" begin - using EpiAwarePipeline pipeline = EpiAwareExamplePipeline() params = make_gi_params(pipeline) @@ -9,15 +8,12 @@ end @testset "make_inf_generating_processes" begin - using EpiAwarePipeline, EpiAware pipeline = EpiAwareExamplePipeline() igps = make_inf_generating_processes(pipeline) @test igps == [DirectInfections, ExpGrowthRate, Renewal] end @testset "make_Rt: returns an array" begin - using EpiAwarePipeline - map([EpiAwareExamplePipeline(), SmoothOutbreakPipeline(), MeasuresOutbreakPipeline(), SmoothEndemicPipeline(), RoughEndemicPipeline()]) do pipeline Rt = make_Rt(pipeline) @@ -26,7 +22,6 @@ end end @testset "default_tspan: returns an Tuple{Integer, Integer}" begin - using EpiAwarePipeline pipeline = EpiAwareExamplePipeline() tspan = make_tspan(pipeline) @@ -34,7 +29,7 @@ end end @testset "make_model_priors: generates a dict with correct keys and distributions" begin - using EpiAwarePipeline, Distributions + using Distributions pipeline = EpiAwareExamplePipeline() priors_dict = make_model_priors(pipeline) @@ -49,7 +44,6 @@ end end @testset "make_epiaware_name_latentmodel_pairs: generates a vector of Pairs with correct keys and latent models" begin - using EpiAwarePipeline, EpiAware pipeline = EpiAwareExamplePipeline() namemodel_vect = make_epiaware_name_latentmodel_pairs(pipeline) @@ -59,7 +53,7 @@ end end @testset "make_inference_method: constructor and defaults" begin - using EpiAwarePipeline, EpiAware, ADTypes, AbstractMCMC + using ADTypes, AbstractMCMC pipeline = EpiAwareExamplePipeline() method = make_inference_method(pipeline) @@ -77,7 +71,9 @@ end @testset "make_inference_method: for prior predictive checking" begin using EpiAwarePipeline, EpiAware, ADTypes, AbstractMCMC - pipeline = RtwithoutRenewalPriorPipeline() + pipetype = [SmoothOutbreakPipeline, MeasuresOutbreakPipeline, + SmoothEndemicPipeline, RoughEndemicPipeline] |> rand + pipeline = pipetype(; ndraws = 100, testmode = true, priorpredictive = true) method = make_inference_method(pipeline) @@ -86,7 +82,6 @@ end end @testset "make_truth_data_configs" begin - using EpiAwarePipeline pipeline = SmoothOutbreakPipeline() example_pipeline = EpiAwareExamplePipeline() @testset "make_truth_data_configs should return a dictionary" begin @@ -107,8 +102,6 @@ end end @testset "default inference configurations" begin - using EpiAwarePipeline - pipeline = SmoothOutbreakPipeline() example_pipeline = EpiAwareExamplePipeline() @@ -134,7 +127,6 @@ end end @testset "make_default_params" begin - using EpiAwarePipeline pipeline = SmoothOutbreakPipeline() # Expected default parameters @@ -155,7 +147,7 @@ end end @testset "make_delay_distribution" begin - using EpiAwarePipeline, Distributions + using Distributions pipeline = SmoothOutbreakPipeline() delay_distribution = make_delay_distribution(pipeline) @test delay_distribution isa Distribution @@ -165,7 +157,6 @@ end end @testset "make_observation_model" begin - using EpiAware # Mock pipeline object pipeline = SmoothOutbreakPipeline() default_params = make_default_params(pipeline) diff --git a/pipeline/test/end-to-end/test_full_inference.jl b/pipeline/test/end-to-end/test_full_inference.jl deleted file mode 100644 index ed05b7271..000000000 --- a/pipeline/test/end-to-end/test_full_inference.jl +++ /dev/null @@ -1,20 +0,0 @@ -using Test -@testset "run inference for random scenario with short toy data" begin - using DrWatson - quickactivate(@__DIR__(), "EpiAwarePipeline") - - using EpiAwarePipeline, EpiAware - pipeline = EpiAwareExamplePipeline() - - tspan = (1, 28) - inference_method = make_inference_method(pipeline) - truth_data_config = make_truth_data_configs(pipeline)[1] - inference_configs = make_inference_configs(pipeline) - inference_config = rand(inference_configs) - truthdata = Dict("y_t" => fill(100, 28), "truth_gi_mean" => 1.5) - - results = generate_inference_results( - truthdata, inference_config, pipeline; tspan, inference_method) - @test results["inference_results"] isa EpiAwareObservables - @test results["forecast_results"] isa EpiAwareObservables -end diff --git a/pipeline/test/end-to-end/test_prior_predictive.jl b/pipeline/test/end-to-end/test_prior_predictive.jl deleted file mode 100644 index 0169ace8f..000000000 --- a/pipeline/test/end-to-end/test_prior_predictive.jl +++ /dev/null @@ -1,20 +0,0 @@ -using Test -@testset "run prior predictive modelling for random scenario" begin - using DrWatson, EpiAware - quickactivate(@__DIR__(), "EpiAwarePipeline") - - using EpiAwarePipeline - pipeline = RtwithoutRenewalPriorPipeline() - - tspan = (1, 28) - inference_method = make_inference_method(pipeline) - truth_data_config = make_truth_data_configs(pipeline)[1] - inference_configs = make_inference_configs(pipeline) - inference_config = rand(inference_configs) - truthdata = Dict("y_t" => fill(100, 28), "truth_gi_mean" => 1.5) - - inference_results = generate_inference_results( - truthdata, inference_config, pipeline; tspan, inference_method) - - @test inference_results["inference_results"] isa EpiAwareObservables -end diff --git a/pipeline/test/forecast/test_forecast.jl b/pipeline/test/forecast/test_forecast.jl index c0b854985..010674ca1 100644 --- a/pipeline/test/forecast/test_forecast.jl +++ b/pipeline/test/forecast/test_forecast.jl @@ -1,5 +1,4 @@ @testset "define_forecast_epiprob" begin - using EpiAwarePipeline pipeline = SmoothOutbreakPipeline() inference_configs = make_inference_configs(pipeline) @@ -11,7 +10,7 @@ epimethod = make_inference_method(pipeline) epiprob = InferenceConfig(rand(inference_configs); case_data, truth_I_t = I_t, - truth_I0 = I0, tspan, epimethod) |> + truth_I0 = I0, tspan, epimethod, pipeline = pipeline) |> define_epiprob @test_throws AssertionError define_forecast_epiprob(epiprob, -1) diff --git a/pipeline/test/infer/test_InferenceConfig.jl b/pipeline/test/infer/test_InferenceConfig.jl index e0d378584..5509fa25b 100644 --- a/pipeline/test/infer/test_InferenceConfig.jl +++ b/pipeline/test/infer/test_InferenceConfig.jl @@ -1,5 +1,5 @@ @testset "InferenceConfig: constructor function" begin - using Distributions, EpiAwarePipeline, EpiAware + using Distributions struct TestLatentModel <: AbstractLatentModel end @@ -16,6 +16,7 @@ latent_model = TestLatentModel() observation_model = TestObsModel() epimethod = TestMethod() + latent_model_name = "TestLatentModel" case_data = [10, 20, 30, 40, 50] I_t = [10, 20, 30, 40, 50] ./ 2 I0 = 1.0 @@ -31,7 +32,9 @@ tspan = tspan, epimethod = epimethod, log_I0_prior = Normal(log(100.0), 1e-5), - lookahead = lookahead + lookahead = lookahead, + latent_model_name = latent_model_name, + pipeline = SmoothOutbreakPipeline() ) @test config.gi_mean == gi_mean @@ -42,13 +45,15 @@ @test config.case_data == case_data @test config.tspan == tspan @test config.epimethod == epimethod + @test config.latent_model_name == latent_model_name end @testset "construct from config dictionary" begin pipeline = SmoothOutbreakPipeline() inference_configs = make_inference_configs(pipeline) @test [InferenceConfig(ic; case_data, truth_I_t = I_t, - truth_I0 = I0, tspan, epimethod) isa InferenceConfig + truth_I0 = I0, tspan, epimethod, pipeline = pipeline) isa + InferenceConfig for ic in inference_configs] |> all end end diff --git a/pipeline/test/infer/test_define_epiprob.jl b/pipeline/test/infer/test_define_epiprob.jl index 0f097aaf2..a59175b96 100644 --- a/pipeline/test/infer/test_define_epiprob.jl +++ b/pipeline/test/infer/test_define_epiprob.jl @@ -1,5 +1,4 @@ @testset "test define_epiprob" begin - using EpiAwarePipeline pipeline = SmoothOutbreakPipeline() inference_configs = make_inference_configs(pipeline) @@ -11,7 +10,7 @@ epimethod = make_inference_method(pipeline) epiprob = InferenceConfig(rand(inference_configs); case_data, tspan, - epimethod, truth_I_t = I_t, truth_I0 = I0) |> + epimethod, truth_I_t = I_t, truth_I0 = I0, pipeline = pipeline) |> define_epiprob @test epiprob isa EpiProblem diff --git a/pipeline/test/infer/test_infer.jl b/pipeline/test/infer/test_infer.jl new file mode 100644 index 000000000..e236ced3b --- /dev/null +++ b/pipeline/test/infer/test_infer.jl @@ -0,0 +1,2 @@ +include("test_define_epiprob.jl") +include("test_InferenceConfig.jl") diff --git a/pipeline/test/mainplots/test_utils.jl b/pipeline/test/mainplots/test_utils.jl deleted file mode 100644 index f2c694f2f..000000000 --- a/pipeline/test/mainplots/test_utils.jl +++ /dev/null @@ -1,21 +0,0 @@ - -@testset "timeseries_samples_into_quantiles tests" begin - using EpiAwarePipeline - X = [1 2 3; 4 5 6; 7 8 9] - qs = [0.25, 0.5, 0.75] - expected_output = [1.5 2.0 2.5 - 4.5 5.0 5.5 - 7.5 8.0 8.5] - @test timeseries_samples_into_quantiles(X, qs) == expected_output -end - -@testset "generate_quantiles_for_targets tests" begin - output = Dict("inference_results" => Dict( - "generated" => [1, 2, 3], "samples" => Dict(:init_incidence => [0.1, 0.2, 0.3]))) - D = EpiData() - qs = [0.25, 0.5, 0.75] - expected_output = [ - [0.1 0.2 0.3; 0.1 0.2 0.3; 0.1 0.2 0.3], [0.1 0.2 0.3; 0.1 0.2 0.3; 0.1 0.2 0.3], - [0.1 0.2 0.3; 0.1 0.2 0.3; 0.1 0.2 0.3]] - @test generate_quantiles_for_targets(output, D, qs) == expected_output -end diff --git a/pipeline/test/pipeline/test_pipeline.jl b/pipeline/test/pipeline/test_pipeline.jl new file mode 100644 index 000000000..ad5f71815 --- /dev/null +++ b/pipeline/test/pipeline/test_pipeline.jl @@ -0,0 +1,2 @@ +include("test_pipelinetypes.jl") +include("test_pipelinefunctions.jl") diff --git a/pipeline/test/pipeline/test_pipelinefunctions.jl b/pipeline/test/pipeline/test_pipelinefunctions.jl index 1a2b87d15..b9836a6bf 100644 --- a/pipeline/test/pipeline/test_pipelinefunctions.jl +++ b/pipeline/test/pipeline/test_pipelinefunctions.jl @@ -1,41 +1,56 @@ @testset "do_truthdata tests" begin - using EpiAwarePipeline, Dagger - pipeline = EpiAwareExamplePipeline() - truthdata_dg_task = do_truthdata(pipeline) - truthdata = fetch.(truthdata_dg_task) - - @test length(truthdata) == 1 - @test all([data["y_t"] isa Vector{Union{Missing, T}} where {T <: Integer} - for data in truthdata]) + using Dagger + for pipetype in [SmoothOutbreakPipeline, MeasuresOutbreakPipeline, + SmoothEndemicPipeline, RoughEndemicPipeline] + pipeline = pipetype(; testmode = true) + truthdata_dg_task = do_truthdata(pipeline) + truthdata = fetch.(truthdata_dg_task) + + @test length(truthdata) == 1 + @test all([data["y_t"] isa Vector{Union{Missing, T}} where {T <: Integer} + for data in truthdata]) + end end @testset "do_inference tests" begin - using EpiAwarePipeline, Dagger, EpiAware - pipeline = EpiAwareExamplePipeline() + using Dagger - function make_inference() - truthdata = do_truthdata(pipeline) + function make_inference(pipeline) + truthdata_dg_task = do_truthdata(pipeline) + truthdata = fetch.(truthdata_dg_task) do_inference(truthdata[1], pipeline) end - inference_results_tsk = make_inference() - inference_results = fetch.(inference_results_tsk) - @test length(inference_results) == 1 - @test all([result["inference_results"] isa EpiAwareObservables - for result in inference_results]) + for pipetype in [SmoothOutbreakPipeline, MeasuresOutbreakPipeline, + SmoothEndemicPipeline, RoughEndemicPipeline] + pipeline = pipetype(; ndraws = 20, nchains = 1, testmode = true) + inference_results_tsk = make_inference(pipeline) + inference_results = fetch.(inference_results_tsk) + @test length(inference_results) == 1 + @test all([result["inference_results"] isa EpiAwareObservables + for result in inference_results]) + end end -@testset "do_pipeline test: just run" begin - using EpiAwarePipeline - pipeline = EpiAwareExamplePipeline() - res = do_pipeline(pipeline) +@testset "do_pipeline test: just run all pipeline objects" begin + using Dagger + pipelines = map([SmoothOutbreakPipeline, MeasuresOutbreakPipeline, + SmoothEndemicPipeline, RoughEndemicPipeline]) do pipetype + pipetype(; ndraws = 10, nchains = 1, testmode = true) + end + + res = do_pipeline(pipelines) fetch(res) @test isnothing(res) end -@testset "do_pipeline test: just run as a vector" begin - using EpiAwarePipeline - pipelines = fill(EpiAwareExamplePipeline(), 2) +@testset "do_pipeline test: prior predictive" begin + using Dagger + pipelines = map([SmoothOutbreakPipeline, MeasuresOutbreakPipeline, + SmoothEndemicPipeline, RoughEndemicPipeline]) do pipetype + pipetype(; ndraws = 10, nchains = 1, testmode = true, priorpredictive = true) + end + res = do_pipeline(pipelines) fetch(res) @test isnothing(res) diff --git a/pipeline/test/pipeline/test_pipelinetypes.jl b/pipeline/test/pipeline/test_pipelinetypes.jl index d36f8332d..f92360de0 100644 --- a/pipeline/test/pipeline/test_pipelinetypes.jl +++ b/pipeline/test/pipeline/test_pipelinetypes.jl @@ -1,5 +1,4 @@ @testset "EpiAwarePipeline Tests" begin - using EpiAwarePipeline @testset "AbstractEpiAwarePipeline" begin @test_throws MethodError AbstractEpiAwarePipeline() end diff --git a/pipeline/test/plotting/plotting_tests.jl b/pipeline/test/plotting/plotting_tests.jl new file mode 100644 index 000000000..bbe9bdd61 --- /dev/null +++ b/pipeline/test/plotting/plotting_tests.jl @@ -0,0 +1,3 @@ +include("test_utils.jl") +include("test_plot_funcs.jl") +include("prior_pred_plot.jl") diff --git a/pipeline/test/plotting/prior_pred_plot.jl b/pipeline/test/plotting/prior_pred_plot.jl new file mode 100644 index 000000000..e34d74dcd --- /dev/null +++ b/pipeline/test/plotting/prior_pred_plot.jl @@ -0,0 +1,17 @@ +@testset "Prior pred plotting" begin + using CairoMakie + #pick a random scenario + pipetype = [SmoothOutbreakPipeline, MeasuresOutbreakPipeline, + SmoothEndemicPipeline, RoughEndemicPipeline] |> rand + P = pipetype(; testmode = true, nchains = 1, ndraws = 2000, priorpredictive = true) + inference_config = make_inference_configs(P) |> rand + inference_config["T"] = 21 + #Add missing data + missingdata = Dict("y_t" => missing, "I_t" => fill(1.0, 100), "truth_I0" => 1.0, + "truth_gi_mean" => inference_config["gi_mean"]) + results = generate_inference_results(missingdata, inference_config, P) + + fig = results["priorpredictive"] + + @test fig isa String #Figure object no longer returned +end diff --git a/pipeline/test/plotting/test_plot_funcs.jl b/pipeline/test/plotting/test_plot_funcs.jl new file mode 100644 index 000000000..5dd2aa8c0 --- /dev/null +++ b/pipeline/test/plotting/test_plot_funcs.jl @@ -0,0 +1,30 @@ +@testset "test plot_truth_data" begin + using CairoMakie + #Toy data + y_t = rand(1:100, 28) + I_t = randn(28) .+ 50 + data = Dict("y_t" => y_t, "I_t" => I_t) + config = Dict("truth_gi_mean" => 1.5) + subdirname = "test" + testpipeline = RtwithoutRenewalPriorPipeline() + + f, path = plot_truth_data( + data, config, testpipeline; plotsubdir = subdirname, saveplot = false) + @test f isa Figure + @test path isa String + @test (splitdir ∘ dirname)(path)[end] == subdirname +end + +@testset "test plot_Rt" begin + using CairoMakie + #Toy data + R_t = randn(100) |> cumsum .|> exp + config = Dict("truth_gi_mean" => 1.5) + subdirname = "test" + testpipeline = RtwithoutRenewalPriorPipeline() + + f, path = plot_Rt(R_t, config, testpipeline; plotsubdir = subdirname, saveplot = false) + @test f isa Figure + @test path isa String + @test (splitdir ∘ dirname)(path)[end] == subdirname +end diff --git a/pipeline/test/plotting/test_utils.jl b/pipeline/test/plotting/test_utils.jl new file mode 100644 index 000000000..bf1952b48 --- /dev/null +++ b/pipeline/test/plotting/test_utils.jl @@ -0,0 +1,9 @@ + +@testset "timeseries_samples_into_quantiles tests" begin + X = [1 2 3; 4 5 6; 7 8 9] + qs = [0.25, 0.5, 0.75] + expected_output = [1.5 2.0 2.5 + 4.5 5.0 5.5 + 7.5 8.0 8.5] + @test timeseries_samples_into_quantiles(X, qs) == expected_output +end diff --git a/pipeline/test/runtests.jl b/pipeline/test/runtests.jl index cd848d7ef..2e111aed5 100644 --- a/pipeline/test/runtests.jl +++ b/pipeline/test/runtests.jl @@ -1,15 +1,14 @@ using DrWatson, Test quickactivate(@__DIR__(), "EpiAwarePipeline") - +using EpiAwarePipeline, EpiAware +import Random +Random.seed!(123) # Run tests -include("pipeline/test_pipelinetypes.jl"); -include("pipeline/test_pipelinefunctions.jl"); -include("utils/test_calculate_processes.jl"); -include("utils/test_simple_crps.jl"); +include("utils/test_utils.jl"); include("constructors/test_constructors.jl"); -include("simulate/test_TruthSimulationConfig.jl"); -include("simulate/test_SimulationConfig.jl"); -include("infer/test_InferenceConfig.jl"); -include("infer/test_define_epiprob.jl"); +include("simulate/test_simulate.jl"); +include("infer/test_infer.jl"); include("forecast/test_forecast.jl"); include("scoring/test_score_parameters.jl"); +include("plotting/plotting_tests.jl"); +include("pipeline/test_pipeline.jl"); diff --git a/pipeline/test/scoring/test_score_parameters.jl b/pipeline/test/scoring/test_score_parameters.jl index e7d25c67d..4de0ea2b1 100644 --- a/pipeline/test/scoring/test_score_parameters.jl +++ b/pipeline/test/scoring/test_score_parameters.jl @@ -1,5 +1,5 @@ @testset "score_parameter tests" begin - using MCMCChains, EpiAwarePipeline + using MCMCChains samples = MCMCChains.Chains(0.5 .+ randn(1000, 2, 1), [:a, :b]) truths = fill(0.5, 2) diff --git a/pipeline/test/simulate/test_SimulationConfig.jl b/pipeline/test/simulate/test_SimulationConfig.jl index 8f8e5bd4b..9cdc7e417 100644 --- a/pipeline/test/simulate/test_SimulationConfig.jl +++ b/pipeline/test/simulate/test_SimulationConfig.jl @@ -1,7 +1,7 @@ # Test the TruthSimulationConfig struct constructor @testset "TruthSimulationConfig" begin - using Distributions, EpiAwarePipeline, EpiAware, LogExpFunctions + using Distributions, LogExpFunctions logit_daily_ascertainment = [fill(1.0, 5); fill(0.5, 2)] normed_daily_ascertainment = logit_daily_ascertainment |> x -> 7 * softmax(x) config = TruthSimulationConfig( diff --git a/pipeline/test/simulate/test_TruthSimulationConfig.jl b/pipeline/test/simulate/test_TruthSimulationConfig.jl index c92732807..c309aed58 100644 --- a/pipeline/test/simulate/test_TruthSimulationConfig.jl +++ b/pipeline/test/simulate/test_TruthSimulationConfig.jl @@ -1,5 +1,5 @@ @testset "simulate runs" begin - using Distributions, EpiAwarePipeline, EpiAware, LogExpFunctions + using Distributions, LogExpFunctions # Define a mock TruthSimulationConfig object for testing logit_daily_ascertainment = [zeros(5); -0.5 * ones(2)] @@ -26,8 +26,6 @@ end @testset "generate_truthdata" begin - using EpiAwarePipeline - pipeline = EpiAwareExamplePipeline() truth_data_config = Dict("gi_mean" => 2.5, "gi_std" => 1.5) truthdata = generate_truthdata(truth_data_config, pipeline; plot = false) diff --git a/pipeline/test/simulate/test_simulate.jl b/pipeline/test/simulate/test_simulate.jl new file mode 100644 index 000000000..f7f9de503 --- /dev/null +++ b/pipeline/test/simulate/test_simulate.jl @@ -0,0 +1,2 @@ +include("test_SimulationConfig.jl") +include("test_TruthSimulationConfig.jl") diff --git a/pipeline/test/utils/test_calculate_processes.jl b/pipeline/test/utils/test_calculate_processes.jl index 150089dd8..f00eb34da 100644 --- a/pipeline/test/utils/test_calculate_processes.jl +++ b/pipeline/test/utils/test_calculate_processes.jl @@ -1,5 +1,4 @@ @testset "calculate_processes" begin - using EpiAware, EpiAwarePipeline using Random rng = MersenneTwister(1234) I0 = 10.0 diff --git a/pipeline/test/utils/test_utils.jl b/pipeline/test/utils/test_utils.jl new file mode 100644 index 000000000..9191e4577 --- /dev/null +++ b/pipeline/test/utils/test_utils.jl @@ -0,0 +1,2 @@ +include("test_calculate_processes.jl") +include("test_simple_crps.jl")