Merge 07566c3 into db3cfd6

REQUIRE

@@ -10,6 +10,7 @@ IterTools
 IterableTables
 FileIO
 GraphPlot 0.3.0
+GlobalSensitivityAnalysis
 JSON
 LightGraphs
 Logging

docs/src/internals/montecarlo.md

@@ -249,7 +249,7 @@ By default, all defined models are run. In some cases, you may want to run some
 
 **Example**
 
-The following example is available in `"Mimi.jl/test/mcs/test_defmcs.jl"`.
+The following example is available in `"Mimi.jl/test/mcs/test_defsim.jl"`.
 
 ```julia
 using Mimi

src/Mimi.jl

@@ -10,7 +10,7 @@ using StringBuilders
 
 export
     @defcomp,
-    @defmcs,
+    @defsim,
     MarginalModel,
     Model,
     add_comp!,  
@@ -39,7 +39,7 @@ export
     parameters, 
     plot_comp_graph,
     replace_comp!, 
-    run_mcs,
+    run_sim,
     set_dimension!, 
     set_leftover_params!, 
     set_models!,

src/core/defs.jl

@@ -761,7 +761,7 @@ end
     copy_external_params(md::ModelDef)
 
 Make copies of ModelParameter subtypes representing external parameters of model `md`. 
-This is used both in the copy() function below, and in the MCS subsystem 
+This is used both in the copy() function below, and in the simulation subsystem 
 to restore values between trials.
 
 """

src/mcs/defmcs.jl

@@ -41,7 +41,7 @@ function _make_dims(args)
         elseif isa(arg, Expr) && arg.head == :tuple  # tuple of Strings/Symbols (@capture didn't work...)
             argtype = typeof(arg.args[1])            # ensure all are same type as first element
             if ! isempty(filter(s -> typeof(s) != argtype, arg.args))
-                error("A parameter dimension tuple must all String or all Symbol (got $arg)")
+                error("A parameter dimension tuple must be all Strings or all Symbols (got $arg)")
             end
             dim = :(convert(Vector{$argtype}, $(arg.args)))
 
@@ -54,12 +54,14 @@ function _make_dims(args)
     return dims
 end
 
-macro defmcs(expr)
+macro defsim(expr)
     let # to make vars local to each macro invocation
         local _rvs        = []
-        local _corrs      = []
         local _transforms = []
         local _saves      = []
+        local _sim_args  = Dict{Symbol, Any}()
+
+        simdatatype = nothing
 
         # distilled into a function since it's called from two branches below
         function saverv(rvname, distname, distargs)
@@ -102,11 +104,6 @@ macro defmcs(expr)
                     end
                 end
 
-            # e.g., name1:name2 = 0.7
-            elseif @capture(elt, name1_:name2_ = value_)
-                expr = :(CorrelationSpec($(QuoteNode(name1)), $(QuoteNode(name2)), $value))
-                push!(_corrs, esc(expr))
-
             # e.g., ext_var5[2010:2050, :] *= name2
             # A bug in Macrotools prevents this shorter expression from working:
             # elseif @capture(elt, ((extvar_  = rvname_Symbol) | 
@@ -152,13 +149,31 @@ macro defmcs(expr)
                     expr = :(TransformSpec($(QuoteNode(extvar)), $(QuoteNode(op)), $(QuoteNode(rvname))))
                 end
                 push!(_transforms, esc(expr))
+
+            elseif @capture(elt, sampling(simdatatype_, simargs__))
+                # General method for setting an SA method's required parameters
+                for arg in simargs
+                    if !@capture(arg, name_=value_)
+                        error("simargs must be given as keyword arguments (got $simargs)")
+                    end
+                    _sim_args[name] = __module__.eval(value)
+                end
             else
-                error("Unrecognized expression '$elt' in @defmcs")
+                error("Unrecognized expression '$elt' in @defsim")
             end
         end
-        return :(MonteCarloSimulation([$(_rvs...)], 
-                                      [$(_transforms...)], 
-                                      CorrelationSpec[$(_corrs...)], 
-                                      Tuple{Symbol, Symbol}[$(_saves...)]))
+
+        # set default
+        simdatatype = (simdatatype == nothing) ? nameof(MCSData) : simdatatype
+
+        # call constructor on given args
+        data = :($simdatatype(; $(_sim_args)...))
+
+        # TBD: need to generalize this to support other methods
+        return :(Simulation{$simdatatype}(
+                    [$(_rvs...)],
+                    [$(_transforms...)],  
+                    Tuple{Symbol, Symbol}[$(_saves...)],
+                    $data))
     end
 end

src/mcs/lhs.jl

@@ -15,6 +15,29 @@ import StatsBase
 using Statistics
 using LinearAlgebra
 
+# Defines a target rank correlation to establish between the two named random vars.
+const CorrelationSpec = Tuple{Symbol, Symbol, Float64}
+
+mutable struct LHSData <: AbstractSimulationData
+    corrlist::Vector{CorrelationSpec}
+
+    function LHSData(;corrlist::Vector{CorrelationSpec}=CorrelationSpec[])
+        new(corrlist)
+    end
+end
+
+const LatinHypercubeSimulation = Simulation{LHSData}
+
+function Base.show(data::LHSData)
+    Mimi.print_nonempty("corrlist", data.corrlist)
+end
+
+function sample!(sim::LatinHypercubeSimulation, samplesize::Int)
+    sim.trials = samplesize
+    corrmatrix = correlation_matrix(sim)
+    lhs!(sim, corrmatrix=corrmatrix)
+end
+
 """
     rank_corr_coef(m::Matrix{Float64})
 
@@ -140,15 +163,15 @@ function lhs(rvlist::Vector{RandomVariable}, trials::Int;
     return asDataFrame ? DataFrame(samples, map(rv->rv.name, rvlist)) : samples
 end
 
-function lhs!(mcs::MonteCarloSimulation; corrmatrix::Union{Matrix{Float64},Nothing}=nothing)
+function lhs!(sim::LatinHypercubeSimulation; corrmatrix::Union{Matrix{Float64},Nothing}=nothing)
     # TBD: verify that any correlated values are actual distributions, not stored vectors?
 
-    trials = mcs.trials
-    rvdict = mcs.rvdict
+    trials = sim.trials
+    rvdict = sim.rvdict
     num_rvs = length(rvdict)
-    rvlist = mcs.dist_rvs
+    rvlist = sim.dist_rvs
 
-    samples = lhs(rvlist, mcs.trials, corrmatrix=corrmatrix, asDataFrame=false)
+    samples = lhs(rvlist, sim.trials, corrmatrix=corrmatrix, asDataFrame=false)
 
     for (i, rv) in enumerate(rvlist)
         dist = rv.dist
@@ -170,13 +193,14 @@ df: Generated by prior call to LHS or something similar.
 rvlist: The random variables to add to the df
 trials: the number of trials to generate for each parameter
 """
-function lhs_amend!(df::DataFrame, rvlist::Vector{RandomVariable}, trials::Int, sampling::SamplingOptions=LHS)
+function lhs_amend!(df::DataFrame, rvlist::Vector{RandomVariable}, trials::Int, 
+                    simtype::T) where T <: AbstractSimulationData
     for rv in rvlist
-        if sampling == LHS
+        if T == LatinHypercubeSimulation
             values = quantile.(rv.dist, _get_percentiles(trials))  # extract values from the RV for these percentiles
             shuffle!(values)                                       # randomize the stratified samples
 
-        else # sampling == Random
+        elseif T == MonteCarloSimulation
             values = rv.dist.rand(trials)
         end
 
@@ -186,31 +210,31 @@ function lhs_amend!(df::DataFrame, rvlist::Vector{RandomVariable}, trials::Int,
 end
 
 """
-    correlation_matrix(mcs::MonteCarloSimulation)
+    correlation_matrix(sim::LatinHypercubeSimulation)
 
 Return a Matrix holding the correlations between random variables
-as indicated in the MonteCarloSimulation, or nothing if no correlations
+as indicated in the Simulation, or nothing if no correlations
 have been defined.
 
 TBD: if needed, compute correlation matrix only for correlated
      RVs, leaving all uncorrelated RVs alone.
 """
-function correlation_matrix(mcs::MonteCarloSimulation)
-    if length(mcs.corrlist) == 0
+function correlation_matrix(sim::LatinHypercubeSimulation)
+    if length(sim.data.corrlist) == 0
         return nothing
     end
 
-    rvdict = mcs.rvdict
+    rvdict = sim.rvdict
 
     # create a mapping of names to RV position in list
     names = Dict([(rv.name, i) for (i, rv) in enumerate(values(rvdict))])
 
     count = length(rvdict)
     corrmatrix = Matrix(1.0I, count, count)
 
-    for corr in mcs.corrlist
-        n1 = corr.name1
-        n2 = corr.name2
+    for corr in sim.data.corrlist
+        n1 = corr[1]
+        n2 = corr[2]
 
         # We don't support correlation between stored samples
         # if rvdict[n1].dist isa SampleStore || rvdict[n2].dist isa SampleStore
@@ -220,7 +244,7 @@ function correlation_matrix(mcs::MonteCarloSimulation)
         i = names[n1]
         j = names[n2]
 
-        corrmatrix[i, j] = corrmatrix[j, i] = corr.value
+        corrmatrix[i, j] = corrmatrix[j, i] = corr[3]
     end
 
     return corrmatrix

src/mcs/mcs.jl

@@ -7,12 +7,14 @@ using StatsBase
 using IterTools
 
 include("mcs_types.jl")
-include("defmcs.jl")
 include("EmpiricalDistribution.jl")
-include("lhs.jl")
 include("montecarlo.jl")
+include("lhs.jl")
+include("sobol.jl")
+include("defmcs.jl")
 
 export 
-    @defmcs, generate_trials!, run_mcs, save_trial_inputs, save_trial_results, set_models!,
-    EmpiricalDistribution, RandomVariable, TransformSpec, CorrelationSpec, MonteCarloSimulation,
-    RANDOM, LHS, INNER, OUTER
+    @defsim, generate_trials!, run_sim, save_trial_inputs, save_trial_results, set_models!,
+    EmpiricalDistribution, RandomVariable, TransformSpec, CorrelationSpec, Simulation, AbstractSimulationData,
+    LHSData, LatinHypercubeSimulation, MCSData, MonteCarloSimulation, SobolData, SobolSimulation,
+    INNER, OUTER, sample!, analyze

src/mcs/mcs_types.jl

@@ -27,7 +27,6 @@ distribution(rv::RandomVariable) = rv.dist
 
 
 @enum ScenarioLoopPlacement OUTER INNER
-@enum SamplingOptions LHS RANDOM
 
 # SampleStore is a faux Distribution that implements base.rand() 
 # to yield stored values.
@@ -77,18 +76,6 @@ struct TransformSpec
     end 
 end
 
-"""
-    CorrelationSpec
-
-Defines a target rank correlation to establish between the two named random vars.
-"""
-struct CorrelationSpec
-    name1::Symbol
-    name2::Symbol
-    value::Float64
-end
-
-
 """
     Base.rand(s::SampleStore{T}, n::Int=1) where T
 
@@ -107,38 +94,36 @@ function Base.reset(s::SampleStore{T}) where T
     return nothing
 end
 
-# debugging tool
-# global save_dict = Dict()
+abstract type AbstractSimulationData end
 
 """
-    MonteCarloSimulation
+    Simulation
     
-Holds all the data that defines a Monte Carlo simulation.
+Holds all the data that defines a simulation.
 """
-mutable struct MonteCarloSimulation
+mutable struct Simulation{T}
     trials::Int
     current_trial::Int
     current_data::Any               # holds data for current_trial when current_trial > 0
     rvdict::OrderedDict{Symbol, RandomVariable}
     translist::Vector{TransformSpec}
-    corrlist::Vector{CorrelationSpec}
     savelist::Vector{Tuple{Symbol, Symbol}}
     dist_rvs::Vector{RandomVariable}
     nt_type::Any                    # a generated NamedTuple type to hold data for a single trial
     models::Vector{Model}
     results::Vector{Dict{Tuple, DataFrame}}
+    data::T
 
-    function MonteCarloSimulation(rvlist::Vector, 
-                                  translist::Vector{TransformSpec}, 
-                                  corrlist::Vector{CorrelationSpec},
-                                  savelist::Vector{Tuple{Symbol, Symbol}})
+    function Simulation{T}(rvlist::Vector, 
+                           translist::Vector{TransformSpec}, 
+                           savelist::Vector{Tuple{Symbol, Symbol}},
+                           data::T) where T <: AbstractSimulationData
         self = new()
         self.trials = 0
         self.current_trial = 0
         self.current_data = nothing
         self.rvdict = OrderedDict([rv.name => rv for rv in rvlist])
         self.translist = translist
-        self.corrlist = corrlist
         self.savelist = savelist
         self.dist_rvs = [rv for rv in rvlist]
 
@@ -149,16 +134,22 @@ mutable struct MonteCarloSimulation
         # These are parallel arrays; each model has a corresponding results dict
         self.models = Vector{Model}(undef, 0)
         self.results = [Dict{Tuple, DataFrame}()]
-        # save_dict[:mcs] = self
+
+        # data specific to a given sensitivity analysis method
+        self.data = data
+
         return self
     end
 end
 
+struct MCSData <: AbstractSimulationData end
+
+const MonteCarloSimulation = Simulation{MCSData}
 
-struct MCSIterator{NT}
-    mcs::MonteCarloSimulation
+struct SimIterator{NT, T}
+    sim::Simulation{T}
 
-    function MCSIterator{NT}(mcs::MonteCarloSimulation) where NT
-        return new{NT}(mcs)
+    function SimIterator{NT}(sim::Simulation{T}) where {NT <: NamedTuple, T <: AbstractSimulationData}
+        return new{NT, T}(sim)
     end
 end