diff --git a/src/GModel.jl b/src/GModel.jl index e42aa48eb..01ef69357 100644 --- a/src/GModel.jl +++ b/src/GModel.jl @@ -1,32 +1,30 @@ module GModel using DocStringExtensions -using Pkg; Pkg.add(PackageSpec(name="Cloudy", version="0.1.0")) -using Cloudy.ParticleDistributions -using Cloudy.Sources -# packages +# TODO: Remove build (which currently prevents segfault): +using Pkg; Pkg.build() + using Random using Sundials # CVODE_BDF() solver for ODE using Distributions using LinearAlgebra using DifferentialEquations -# exports export run_G export run_G_ensemble # TODO: It would be nice to have run_G_ensemble take a pointer to the # function G (called run_G below), which maps the parameters u to G(u), -# as an input argument. In the example below, run_G runs Cloudy, but in -# general run_G will be supplied by the user and run the specific model -# the user wants to do UQ with. +# as an input argument. In the example below, run_G runs the model, but +# in general run_G will be supplied by the user and run the specific +# model the user wants to do UQ with. # So, the interface would ideally be something like: # run_G_ensemble(params, G) where G is user-defined """ GSettings{FT<:AbstractFloat, KT, D} -Structure to hold all information to run the forward model G (Cloudy) +Structure to hold all information to run the forward model G # Fields $(DocStringExtensions.FIELDS) @@ -44,23 +42,28 @@ end """ - run_G_ensemble(params::Array{FT, 2}, settings::GSettings; - rng_seed=42) - - Run the forward model G for an array of parameters by iteratively - calling run_G for each of the N_ensemble parameter values. - Return g_ens, an array of size N_ensemble x N_data, where - g_ens[j,:] = G(params[j,:]) - - - `params` - array of size N_ensemble x N_parameters containing the - parameters for which G will be run - - `settings` - a GSetttings struct + run_G_ensemble(params::Array{FT, 2}, + settings::GSettings{FT}, + update_params, + moment, + get_src; + rng_seed=42) where {FT<:AbstractFloat} + +Run the forward model G for an array of parameters by iteratively +calling run_G for each of the N_ensemble parameter values. +Return g_ens, an array of size N_ensemble x N_data, where +g_ens[j,:] = G(params[j,:]) + + - `params` - array of size N_ensemble x N_parameters containing the + parameters for which G will be run + - `settings` - a GSetttings struct """ function run_G_ensemble(params::Array{FT, 2}, settings::GSettings{FT}, update_params, - moment; + moment, + get_src; rng_seed=42) where {FT<:AbstractFloat} N_ens = size(params, 1) # params is N_ens x N_params @@ -69,8 +72,8 @@ function run_G_ensemble(params::Array{FT, 2}, Random.seed!(rng_seed) for i in 1:N_ens - # run Cloudy with the current parameters, i.e., map θ to G(θ) - g_ens[i, :] = run_G(params[i, :], settings, update_params, moment) + # run the model with the current parameters, i.e., map θ to G(θ) + g_ens[i, :] = run_G(params[i, :], settings, update_params, moment, get_src) end return g_ens @@ -78,19 +81,23 @@ end """ - run_G(u::Array{FT, 1}, settings::GSettings) + run_G(u::Array{FT, 1}, + settings::GSettings{FT}, + update_params, + moment, + get_src) where {FT<:AbstractFloat} - Return g_u = G(u), a vector of length N_data (=N_moments in the case - of Cloudy) +Return g_u = G(u), a vector of length N_data. - - `u` - parameter vector of length N_parameters - - `settings` - a GSettings struct + - `u` - parameter vector of length N_parameters + - `settings` - a GSettings struct """ function run_G(u::Array{FT, 1}, settings::GSettings{FT}, update_params, - moment) where {FT<:AbstractFloat} + moment, + get_src) where {FT<:AbstractFloat} # generate the initial distribution dist = update_params(settings.dist, u) @@ -107,7 +114,7 @@ function run_G(u::Array{FT, 1}, end # Set up ODE problem: dM/dt = f(M,p,t) - rhs(M, p, t) = get_int_coalescence(M, dist, settings.kernel) + rhs(M, p, t) = get_src(M, dist, settings.kernel) prob = ODEProblem(rhs, moments_init, settings.tspan) # Solve the ODE sol = solve(prob, CVODE_BDF(), alg_hints=[:stiff], reltol=tol, abstol=tol) @@ -117,4 +124,4 @@ function run_G(u::Array{FT, 1}, return moments_final end -end # module GModel +end # module diff --git a/test/Cloudy/runtests.jl b/test/Cloudy/runtests.jl index b405dd66d..aa8a97ed0 100644 --- a/test/Cloudy/runtests.jl +++ b/test/Cloudy/runtests.jl @@ -2,8 +2,6 @@ # Import Cloudy modules using Pkg; Pkg.add(PackageSpec(name="Cloudy", version="0.1.0")) using Cloudy -using Cloudy.KernelTensors -using Cloudy.ParticleDistributions const PDistributions = Cloudy.ParticleDistributions Pkg.add("Plots") @@ -49,7 +47,7 @@ dist_true = PDistributions.Gamma(N0_true, θ_true, k_true) # Collision-coalescence kernel to be used in Cloudy coalescence_coeff = 1/3.14/4 kernel_func = x -> coalescence_coeff -kernel = CoalescenceTensor(kernel_func, 0, 100.0) +kernel = Cloudy.KernelTensors.CoalescenceTensor(kernel_func, 0, 100.0) # Time period over which to run Cloudy tspan = (0., 0.5) @@ -62,7 +60,8 @@ tspan = (0., 0.5) g_settings_true = GModel.GSettings(kernel, dist_true, moments, tspan) yt = GModel.run_G(u_true, g_settings_true, PDistributions.update_params, - PDistributions.moment) + PDistributions.moment, + Cloudy.Sources.get_int_coalescence) n_samples = 100 samples = zeros(n_samples, length(yt)) noise_level = 0.05 @@ -106,7 +105,8 @@ for i in 1:N_iter g_ens = GModel.run_G_ensemble(exp_transform(ekiobj.u[end]), g_settings, PDistributions.update_params, - PDistributions.moment + PDistributions.moment, + Cloudy.Sources.get_int_coalescence ) EKI.update_ensemble!(ekiobj, g_ens) end