Remove GModel dependence on Cloudy

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,28 +72,32 @@ 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
 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

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