Simulation update funcs

src/core/defs.jl

@@ -38,7 +38,8 @@ end
 
 Return the name of `def`.  Possible `NamedDef`s include `DatumDef`, and `ComponentDef`.
 """
-name(def::NamedDef) = def.name
+name(def::NamedDef) = def.name          # old definition; should deprecate this...
+Base.nameof(def::NamedDef) = def.name   # 'nameof' is the more julian name
 
 number_type(md::ModelDef) = md.number_type
 

src/mcs/defmcs.jl

@@ -169,4 +169,67 @@ macro defsim(expr)
                     Tuple{Symbol, Symbol}[$(_saves...)],
                     $data))
     end
-end
+end
+
+#
+# Simulation update methods
+#
+function _update_nt_type!(sim::Simulation{T}) where T <: AbstractSimulationData
+    names = (keys(sim.rvdict)...,)
+    types = [eltype(fld) for fld in values(sim.rvdict)]
+    sim.nt_type = NamedTuple{names, Tuple{types...}}
+    nothing
+end
+
+function deleteRV!(sim::Simulation, name::Symbol)
+    deleteTransform!(sim, name)
+    delete!(sim.rvdict, name)
+    _update_nt_type!(sim)
+end
+
+function addRV!(sim::Simulation, rv::RandomVariable)
+    name = rv.name
+    haskey(sim.rvdict, name) && error("Simulation already has RV :$name. Use replaceRV! to replace it.")
+    sim.rvdict[name] = rv
+    _update_nt_type!(sim)
+end
+
+addRV!(sim::Simulation, name::Symbol, dist::Distribution) = addRV!(sim, RandomVariable(name, dist))
+
+function replaceRV!(sim::Simulation, rv::RandomVariable)
+    sim.rvdict[rv.name] = rv
+    _update_nt_type!(sim)
+end
+
+replaceRV!(sim::Simulation, name::Symbol, dist::Distribution) = replaceRV!(sim, RandomVariable(name, dist))
+
+function deleteTransform!(sim::Simulation, name::Symbol)
+    pos = findall(t -> t.rvname == name, sim.translist)
+    deleteat!(sim.translist, pos)    
+    _update_nt_type!(sim)
+end
+
+function addTransform!(sim::Simulation, t::TransformSpec)
+    push!(sim.translist, t)
+    _update_nt_type!(sim)
+end
+
+function addTransform!(sim::Simulation, paramname::Symbol, op::Symbol, rvname::Symbol, dims::Vector{T}=[]) where T
+    addTransform!(sim, TransformSpec(paramname, op, rvname, dims))
+end
+
+
+function deleteSave!(sim::Simulation, key::Tuple{Symbol, Symbol})
+    pos = findall(isequal(key), sim.savelist)
+    deleteat!(sim.savelist, pos)
+end
+
+deleteSave!(sim::Simulation, comp_name::Symbol, datum_name::Symbol) = deleteSave!(sim, (comp_name, datum_name))
+
+function addSave!(sim::Simulation, key::Tuple{Symbol, Symbol})
+    deleteSave!(sim, key) 
+    push!(sim.savelist, key)
+    nothing
+end
+
+addSave!(sim::Simulation, comp_name::Symbol, datum_name::Symbol) = addSave!(sim, (comp_name, datum_name))

src/mcs/lhs.jl

@@ -163,18 +163,21 @@ function lhs(rvlist::Vector{RandomVariable}, trials::Int;
     return asDataFrame ? DataFrame(samples, map(rv->rv.name, rvlist)) : samples
 end
 
+function lhs(sim::LatinHypercubeSimulation; 
+             corrmatrix::Union{Matrix{Float64},Nothing}=nothing,
+             asDataFrame::Bool=true)
+    rvlist = collect(values(sim.rvdict))
+    lhs(rvlist, sim.trials, corrmatrix=corrmatrix, asDataFrame=asDataFrame)
+end
+
 function lhs!(sim::LatinHypercubeSimulation; corrmatrix::Union{Matrix{Float64},Nothing}=nothing)
     # TBD: verify that any correlated values are actual distributions, not stored vectors?
 
-    trials = sim.trials
     rvdict = sim.rvdict
-    num_rvs = length(rvdict)
-    rvlist = sim.dist_rvs
 
-    samples = lhs(rvlist, sim.trials, corrmatrix=corrmatrix, asDataFrame=false)
+    samples = lhs(sim, corrmatrix=corrmatrix, asDataFrame=false)
 
-    for (i, rv) in enumerate(rvlist)
-        dist = rv.dist
+    for (i, rv) in enumerate(values(rvdict))
         name = rv.name
         values = samples[:, i]
         rvdict[name] = RandomVariable(name, SampleStore(values))

src/mcs/mcs_types.jl

@@ -1,6 +1,8 @@
 using Distributions
 using Statistics
 
+@enum ScenarioLoopPlacement OUTER INNER
+
 """
     RandomVariable{T}
 
@@ -24,17 +26,34 @@ Base.eltype(rv::RandomVariable) = eltype(rv.dist)
 
 distribution(rv::RandomVariable) = rv.dist
 
+abstract type PseudoDistribution end
+
+"""     ReshapedDistribution
+A pseudo-distribution that returns a reshaped array of values from the
+stored distribution and dimensions.
+
+Example:
+    rd = ReshapedDistribution([5, 5], Dirichlet(25,1))
+"""
+struct ReshapedDistribution <: PseudoDistribution
+    dims::Vector{Int}
+    dist::Distribution
+end
+
+function Base.rand(rd::ReshapedDistribution, draws::Int=1)
+    return [reshape(rand(rd.dist), rd.dims...) for i in 1:draws]
+end
 
-@enum ScenarioLoopPlacement OUTER INNER
 
 # SampleStore is a faux Distribution that implements base.rand() 
 # to yield stored values.
-mutable struct SampleStore{T}
-    values::Vector{T}   # generally Int or Float64
-    idx::Int            # index of next value to return
+mutable struct SampleStore{T} <: PseudoDistribution
+    values::Vector{T}       # generally Int or Float64
+    idx::Int                # index of next value to return
+    dist::Union{Nothing, Distribution, PseudoDistribution}  # original distribution, if any
 
-    function SampleStore(values::Vector{T}) where T
-        return new{T}(values, 1)
+    function SampleStore(values::Vector{T}, dist::Union{Nothing, Distribution, PseudoDistribution}=nothing) where T
+        return new{T}(values, 1, dist)
     end
 end
 
@@ -54,31 +73,13 @@ end
 # Probably shouldn't use correlation on values loaded from a file rather than 
 # from a proper distribution.
 function Statistics.quantile(ss::SampleStore{T}, probs::AbstractArray) where T
-    return quantile.(ss, probs)
-end
-
-
-"""     ReshapedDistribution
-A pseudo-distribution that returns a reshaped array of values from the
-stored distribution and dimensions.
-
-Example:
-    rd = ReshapedDistribution([5, 5], Dirichlet(25,1))
-"""
-struct ReshapedDistribution
-    dims::Vector{Int}
-    dist::Distribution
+    return quantile(sort(ss.values), probs)
 end
 
-# function Base.rand(rd::ReshapedDistribution, draws::Int=1)
-#     values = rand(rd.dist, draws)
-#     dims = (draws == 1 ? rd.dims : [rd.dims..., draws])
-#     return reshape(values, dims...)
-# end
+Base.length(ss::SampleStore{T}) where T = length(ss.values)
 
-function Base.rand(rd::ReshapedDistribution, draws::Int=1)
-    return [reshape(rand(rd.dist), rd.dims...) for i in 1:draws]
-end
+Base.iterate(ss::SampleStore{T}) where T = iterate(ss.values)
+Base.iterate(ss::SampleStore{T}, idx) where T = iterate(ss.values, idx)
 
 struct TransformSpec
     paramname::Symbol
@@ -127,7 +128,6 @@ mutable struct Simulation{T}
     rvdict::OrderedDict{Symbol, RandomVariable}
     translist::Vector{TransformSpec}
     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}}
@@ -144,7 +144,6 @@ mutable struct Simulation{T}
         self.rvdict = OrderedDict([rv.name => rv for rv in rvlist])
         self.translist = translist
         self.savelist = savelist
-        self.dist_rvs = [rv for rv in rvlist]
 
         names = (keys(self.rvdict)...,)
         types = [eltype(fld) for fld in values(self.rvdict)]

src/mcs/montecarlo.jl

@@ -174,9 +174,12 @@ function Random.rand!(sim::Simulation{T}) where T <: AbstractSimulationData
     rvdict = sim.rvdict
     trials = sim.trials
 
-    for rv in sim.dist_rvs
-        values = rand(rv.dist, trials)
-        rvdict[rv.name] = RandomVariable(rv.name, SampleStore(values))
+    for rv in values(sim.rvdict)
+        # use underlying distribution, if known
+        orig_dist = (rv.dist isa SampleStore ? rv.dist.dist : rv.dist)
+        dist = (orig_dist === nothing ? rv.dist : orig_dist)
+        values = rand(dist, trials)
+        rvdict[rv.name] = RandomVariable(rv.name, SampleStore(values, orig_dist))
     end
 end
 

src/mcs/sobol.jl

@@ -15,34 +15,43 @@ end
 
 const SobolSimulation = Simulation{SobolData}
 
-function sample!(sim::SobolSimulation, samplesize::Int)
+function _compute_N(sim::SobolSimulation)
+    num_rvs = length(sim.rvdict)
+    factor = (sim.data.calc_second_order ? 2 : 1)
+    N = sim.trials / (factor * num_rvs + 2)
+    return N
+end
 
-    rvdict = sim.rvdict
-    rvlist = sim.dist_rvs
-    num_rvs = length(rvdict)
+function _compute_trials(sim::SobolSimulation, N::Int)
+    num_rvs = length(sim.rvdict)
+    factor = (sim.data.calc_second_order ? 2 : 1)
+    sim.trials = N * (factor * num_rvs + 2)
+end
 
-    if sim.data.calc_second_order
-        sim.trials = samplesize * (2 * num_rvs + 2)
-    else
-        sim.trials = samplesize * (num_rvs + 2)
-    end
+# Use original distribution when resampling from SampleStore
+_get_dist(rv::RandomVariable) = (rv.dist isa SampleStore ? rv.dist.dist : rv.dist)
+
+function sample!(sim::SobolSimulation, samplesize::Int)
+    rvdict = sim.rvdict
+    sim.trials = _compute_trials(sim, samplesize)
 
     # get the samples to plug in to trials
     payload = create_GSA_payload(sim)
     samples = GlobalSensitivityAnalysis.sample(payload)
 
-    for (i, rv) in enumerate(rvlist)
-        dist = rv.dist
+    for (i, rv) in enumerate(values(rvdict))
+        # use underlying distribution, if known
+        orig_dist = _get_dist(rv)
+
         name = rv.name
         values = samples[:, i]
-        rvdict[name] = RandomVariable(name, SampleStore(values))
+        rvdict[name] = RandomVariable(name, SampleStore(values, orig_dist))
     end
 end
 
 function analyze(sim::SobolSimulation, model_output::AbstractArray{<:Number, N1}; N::Union{Nothing, Int}=nothing) where N1
-    if N != nothing
-        sim.data.calc_second_order ? factor = 2 * length(sim.rvdict) : factor = length(sim.rvdict)
-        sim.trials = N * (factor + 2)
+    if N !== nothing
+        sim.trials = _compute_trials(sim, N)
     end
 
     if sim.trials == 0
@@ -54,25 +63,11 @@ function analyze(sim::SobolSimulation, model_output::AbstractArray{<:Number, N1}
 end
 
 function create_GSA_payload(sim::SobolSimulation)
-
-    rvlist = sim.dist_rvs
-
-    # add all distinct rvs to the rv_info dictionary to be passed to GSA's 
-    # SobolData payload
-    rv_info = OrderedDict{Symbol, Any}(rvlist[1].name => rvlist[1].dist)
-    for rv in rvlist[2:end] 
-        rv_info[rv.name] = rv.dist
-    end
+    rv_info = OrderedDict{Symbol, Any}([name => _get_dist(rv) for (name, rv) in sim.rvdict])
 
     # back out N
-    num_rvs = length(sim.rvdict)
-    if sim.data.calc_second_order
-        samples = sim.trials / (2 * num_rvs + 2)
-    else
-        samples = sim.trials / (num_rvs + 2)
-    end
+    N = _compute_N(sim)
 
     # create payload
-    return GlobalSensitivityAnalysis.SobolData(params = rv_info, calc_second_order = sim.data.calc_second_order, N = samples)
-
+    return GlobalSensitivityAnalysis.SobolData(params = rv_info, calc_second_order = sim.data.calc_second_order, N = N)
 end

test/mcs/test_defmcs.jl

@@ -11,7 +11,6 @@ using Test
 using Mimi: reset_compdefs, modelinstance, compinstance, 
             get_var_value, OUTER, INNER, ReshapedDistribution
 
-reset_compdefs()
 include("../../examples/tutorial/02-two-region-model/two-region-model.jl")
 using .MyModel
 m = construct_MyModel()
@@ -46,7 +45,6 @@ sim = @defsim begin
 end
 
 
-
 # Optionally, user functions can be called just before or after a trial is run
 function print_result(m::Model, sim::Simulation, trialnum::Int)
     ci = Mimi.compinstance(m.mi, :emissions)
@@ -227,4 +225,3 @@ trial2 = copy(sim2.rvdict[:name1].dist.values)
 
 @test length(trial1) == length(trial2)
 @test trial1 != trial2
-