Merge pull request #409 from mimiframework/mcs-reshape-dist

- Added support for distributions that produce a matrix of values

Author: lrennels

.travis.yml

@@ -3,7 +3,6 @@ os:
   - linux
   - osx
 julia:
-  - 1.0
   - 1.1
   - nightly
 matrix:

REQUIRE

@@ -1,4 +1,4 @@
-julia 1.0
+julia 1.1
 Compose
 CSVFiles
 DataFrames

appveyor.yml

@@ -1,6 +1,5 @@
 environment:
   matrix:
-  - julia_version: 1.0
   - julia_version: 1.1
   - julia_version: nightly
 

src/mcs/defmcs.jl

@@ -108,13 +108,6 @@ macro defmcs(expr)
                 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) | 
-            #                       (extvar_ += rvname_Symbol) |
-            #                       (extvar_ *= rvname_Symbol) |
-            #                       (extvar_  = distname_(distargs__)) | 
-            #                       (extvar_ += distname_(distargs__)) |
-            #                       (extvar_ *= distname_(distargs__))))
             elseif (@capture(elt, extvar_  = rvname_Symbol) ||
                     @capture(elt, extvar_ += rvname_Symbol) ||
                     @capture(elt, extvar_ *= rvname_Symbol) ||

src/mcs/mcs.jl

@@ -15,4 +15,4 @@ include("montecarlo.jl")
 export 
     @defmcs, generate_trials!, run_mcs, save_trial_inputs, save_trial_results, set_models!,
     EmpiricalDistribution, RandomVariable, TransformSpec, CorrelationSpec, MonteCarloSimulation,
-    RANDOM, LHS, INNER, OUTER
+    RANDOM, LHS, INNER, OUTER, ReshapedDistribution

src/mcs/mcs_types.jl

@@ -42,8 +42,6 @@ end
 
 Base.eltype(ss::SampleStore) = eltype(ss.values)
 
-distribution(ss::SampleStore) = ss.dist
-
 #
 # TBD: maybe have a different SampleStore subtype for values drawn from a dist
 # versus those loaded from a file, which would be treated as immutable?
@@ -62,6 +60,28 @@ function Statistics.quantile(ss::SampleStore{T}, probs::AbstractArray) where T
 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
+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
+
+function Base.rand(rd::ReshapedDistribution, draws::Int=1)
+    return [reshape(rand(rd.dist), rd.dims...) for i in 1:draws]
+end
+
 struct TransformSpec
     paramname::Symbol
     op::Symbol

src/mcs/montecarlo.jl

@@ -3,6 +3,7 @@ using IterableTables
 using TableTraits
 using Random
 using ProgressMeter
+using Serialization
 
 function Base.show(io::IO, mcs::MonteCarloSimulation)
     println("MonteCarloSimulation")
@@ -83,11 +84,11 @@ function _store_trial_results(mcs::MonteCarloSimulation, trialnum::Int)
 end
 
 """
-    save_trial_results(mcs::MonteCarloSimulation, output_dir::String)
+    save_trial_results(mcs::MonteCarloSimulation, output_dir)
 
 Save the stored MCS results to files in the directory `output_dir`
 """
-function save_trial_results(mcs::MonteCarloSimulation, output_dir::AbstractString)
+function save_trial_results(mcs::MonteCarloSimulation, output_dir)
     multiple_results = (length(mcs.results) > 1)
 
     for (i, results) in enumerate(mcs.results)
@@ -106,8 +107,44 @@ function save_trial_results(mcs::MonteCarloSimulation, output_dir::AbstractStrin
     end
 end
 
-function save_trial_inputs(mcs::MonteCarloSimulation, filename::String)
-    mkpath(dirname(filename), mode=0o770)   # ensure that the specified path exists
+function _dummy_writer(mcs, pname, dir)
+    @info "_dummy_writer(pname=$pname, dir=$dir)"
+end
+
+
+"""
+    _serialize(mcs::MonteCarloSimulation, pname, dir)
+
+This is the default function used for writing non-bit-type (array-valued)
+random vars, e.g., those produced by ReshapedDistribution.
+
+This default method uses Serialization.serialize to write out the array.
+No application-defined types are written, so the file should be robust.
+"""
+function _serialize(mcs::MonteCarloSimulation, pname, dir)
+    rv = mcs.rvdict[pname]
+    if rv isa RandomVariable{Mimi.SampleStore{T}} where T
+        name = first(split(string(pname), "!"))     # e.g., :alpha!1 --> "alpha"
+        path = joinpath(dir, "$name.dat")
+        serialize(path, mcs.rvdict[pname].dist.values)
+    else
+        error("Tried to _serialize $pname, which isn't a RandomVariable{SampleStore{T}}")
+    end
+end
+
+function save_trial_inputs(mcs::MonteCarloSimulation, filename, non_bit_writer=_serialize)
+    dir = dirname(filename)
+    mkpath(dir, mode=0o770)   # ensure that the specified path exists
+
+    # If ReshapedDistribution was used, a single trial value for a field
+    # will be an Array of values. CSV format doesn't handle these well,
+    # so we serialize the arrays into a file using the field's name.
+    non_bits = [name for (name, T) in zip(column_names(mcs), column_types(mcs)) if ! isbitstype(T)]
+    for pname in non_bits
+        non_bit_writer(mcs, pname, dir)
+    end
+
+    # TBD: avoid writing array values to CSV files...
     save(filename, mcs)
     return nothing
 end
@@ -190,7 +227,7 @@ function _copy_mcs_params(mcs::MonteCarloSimulation)
 
     for (i, m) in enumerate(mcs.models)
         md = m.mi.md
-        param_vec[i] = Dict{Symbol, ModelParameter}(trans.paramname => copy(external_param(md, trans.paramname)) for trans in mcs.translist)
+        param_vec[i] = Dict{Symbol, ModelParameter}(trans.paramname => deepcopy(external_param(md, trans.paramname)) for trans in mcs.translist)
     end
 
     return param_vec
@@ -225,7 +262,13 @@ function _param_indices(param::ArrayModelParameter{T}, md::ModelDef, trans::Tran
     num_pdims = length(pdims)
 
     tdims  = trans.dims
-    num_dims = length(tdims)
+    num_dims = length(tdims) 
+
+    # special case for handling reshaped data where a single draw returns a matrix of values
+    if num_dims == 0
+        indices = repeat([Colon()], num_pdims)
+        return indices
+    end
 
     if num_pdims != num_dims
         pname = trans.paramname
@@ -256,7 +299,8 @@ function _perturb_param!(param::ScalarModelParameter{T}, md::ModelDef, trans::Tr
     end
 end
 
-function _perturb_param!(param::ArrayModelParameter{T}, md::ModelDef, trans::TransformSpec, rvalue::Number) where T
+function _perturb_param!(param::ArrayModelParameter{T}, md::ModelDef, 
+                         trans::TransformSpec, rvalue::Union{Number, Array{<: Number, N}}) where {T, N}
     op = trans.op
     pvalue = value(param)
     indices = _param_indices(param, md, trans)
@@ -539,6 +583,9 @@ IterableTables.getiterator(mcs::MonteCarloSimulation) = MCSIterator{mcs.nt_type}
 column_names(mcs::MonteCarloSimulation) = fieldnames(mcs.nt_type)
 column_types(mcs::MonteCarloSimulation) = [eltype(fld) for fld in values(mcs.rvdict)]
 
+# TBD: strip the "!1" off the end of each field name?
+# column_names(iter::MCSIterator) = Tuple([first(split(string(name), "!")) for name in fieldnames(iter.mcs.nt_type)])
+
 column_names(iter::MCSIterator) = column_names(iter.mcs)
 column_types(iter::MCSIterator) = IterableTables.column_types(iter.mcs)
 

test/mcs/runtests.jl

@@ -9,4 +9,6 @@ using Test
     @info("test_defmcs.jl")
     include("test_defmcs.jl")
 
+    @info("test_reshaping.jl")
+    include("test_reshaping.jl")    
 end

test/mcs/test-model/emissions.jl

@@ -0,0 +1,24 @@
+@defcomp emissions begin
+    regions     = Index()                           # The regions index must be specified for each component
+
+    E           = Variable(index=[time, regions])   # Total greenhouse gas emissions
+    E_Global    = Variable(index=[time])            # Global emissions (sum of regional emissions)
+    sigma       = Parameter(index=[time, regions])  # Emissions output ratio
+    YGROSS      = Parameter(index=[time, regions])  # Gross output - Note that YGROSS is now a parameter
+
+    # function init(p, v, d)
+    # end
+    
+    function run_timestep(p, v, d, t)
+        # Define an equation for E
+        for r in d.regions
+            v.E[t,r] = p.YGROSS[t,r] * p.sigma[t,r]
+        end
+
+        # Define an equation for E_Global
+        for r in d.regions
+            v.E_Global[t] = sum(v.E[t,:])
+        end
+    end
+
+end

test/mcs/test-model/gross_economy.jl

@@ -0,0 +1,32 @@
+@defcomp grosseconomy begin
+    regions = Index()                           #Note that a regional index is defined here
+
+    YGROSS  = Variable(index=[time, regions])   #Gross output
+    K       = Variable(index=[time, regions])   #Capital
+    l       = Parameter(index=[time, regions])  #Labor
+    tfp     = Parameter(index=[time, regions])  #Total factor productivity
+    s       = Parameter(index=[time, regions])  #Savings rate
+    depk    = Parameter(index=[regions])        #Depreciation rate on capital - Note that it only has a region index
+    k0      = Parameter(index=[regions])        #Initial level of capital
+    share   = Parameter()                       #Capital share
+
+    tester  = Parameter(index=[time, regions])  # test for reshaped Dirichlet distribution
+
+    function run_timestep(p, v, d, t)
+    # Note that the regional dimension is defined in d and parameters and variables are indexed by 'r'
+
+        # Define an equation for K
+        for r in d.regions
+            if is_first(t)
+                v.K[t,r] = p.k0[r]
+            else
+                v.K[t,r] = (1 - p.depk[r])^5 * v.K[t-1,r] + v.YGROSS[t-1,r] * p.s[t-1,r] * 5
+            end
+        end
+
+        # Define an equation for YGROSS
+        for r in d.regions
+            v.YGROSS[t,r] = p.tfp[t,r] * v.K[t,r]^p.share * p.l[t,r]^(1-p.share)
+        end
+    end
+end