diff --git a/README.md b/README.md index a8784d6..3a6175d 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,15 @@ +# Forked version of SMM.jl + +This is a modified version of SMM.jl with two changes for the STBNews project: + +1. Change to parallel dispatch to improve memory usage (in `algoBGP.jl`). Note STB code uses `algoSTB.jl`, in `stb-model` project. +2. Allow passing options to objective function (in `mprob.jl::evaluateObjective`) + +To install, in your julia REPL, type + +```julia +] add https://github.com/gregobad/SMM.jl#better_parallel +``` # SMM.jl: Simulated Method of Moments for Julia @@ -11,15 +23,6 @@ This package provides a `Julia` infrastructure for *[Simulated Method of Moments](http://en.wikipedia.org/wiki/Method_of_simulated_moments)* estimation, or other problems where we want to optimize a non-differentiable objective function. The setup is suitable for all kinds of **likelihood-free estimators** - in general, those require evaluating the objective at many regions. The user can supply their own algorithms for generating successive new parameter guesses. We provide a set of MCMC template algorithms. The code can be run in serial or on a cluster. - -## Installation - -In your julia REPL, type - -```julia -] add SMM -``` - ## Documentation [![][docs-stable-img]][docs-stable-url] @@ -105,7 +108,7 @@ New algorithms: ## History -This package grew out of the [R package `mopt`](https://github.com/tlamadon/mopt). +This package grew out of the [R package `mopt`](https://github.com/tlamadon/mopt). ## Thanks to all Contributors! diff --git a/src/SMM.jl b/src/SMM.jl index 052caec..43e09a8 100644 --- a/src/SMM.jl +++ b/src/SMM.jl @@ -5,7 +5,7 @@ module SMM # Dependencies # ############ -using Distributions +using Distributions using DataFrames, DataFramesMeta using OrderedCollections using PDMats @@ -28,16 +28,17 @@ import Base.get, Base.write, Base.== import Base.show, Statistics.mean, Statistics.median # exports: Types -export MProb, Eval,MAlgo, MAlgoBGP +export MProb, Eval, MAlgo, MAlgoBGP # exports: methods -export addParam!, - addSampledParam!, - addMoment!, - addEvalFunc!, - setMoments!, +export addParam!, + addSampledParam!, + addMoment!, + addEvalFunc!, + addEvalFuncOpts!, + setMoments!, setValue!, - readEvalArray , + readEvalArray , dataMoment, dataMomentd, dataMomentW, @@ -75,7 +76,3 @@ include("mopt/econometrics.jl") end # module - - - - diff --git a/src/mopt/AlgoBGP.jl b/src/mopt/AlgoBGP.jl index 76ad779..5b6d65b 100644 --- a/src/mopt/AlgoBGP.jl +++ b/src/mopt/AlgoBGP.jl @@ -56,7 +56,7 @@ mutable struct BGPChain <: AbstractChain sigma_update_steps :: Int64 # update sampling vars every sigma_update_steps iterations sigma_adjust_by :: Float64 # adjust sampling vars by sigma_adjust_by percent up or down smpl_iters :: Int64 # max number of trials to get a new parameter from MvNormal that lies within support - min_improve :: Float64 + min_improve :: Float64 batches :: Vector{UnitRange{Int}} # vector of indices to update together. """ @@ -93,7 +93,7 @@ mutable struct BGPChain <: AbstractChain this.m = m # how many bundles + rest nb, rest = divrem(np,batch_size) - this.sigma = sig + this.sigma = sig this.batches = UnitRange{Int}[] i = 1 for ib in 1:nb @@ -262,7 +262,7 @@ end Computes the next `Eval` for chain `c`: -1. Get last accepted param +1. Get last accepted param 2. get a new param via [`proposal`](@ref) 3. [`evaluateObjective`](@ref) 4. Accept or Reject the new value via [`doAcceptReject!`](@ref) @@ -293,6 +293,28 @@ function next_eval(c::BGPChain) end +""" + next_acceptreject(c::BGPChain, ev::Eval) + +Stores `ev` to chain `c`, given already computed obj fn value: + +1. Accept or Reject the new value via [`doAcceptReject!`](@ref) +2. Store `Eval` on chain `c`. + +""" +function next_acceptreject(c::BGPChain, ev::Eval) + @debug "iteration = $(c.iter)" + + # accept reject + doAcceptReject!(c,ev) + + # save eval on STBChain + set_eval!(c,ev) + + return c + +end + """ doAcceptReject!(c::BGPChain,eval_new::Eval) @@ -392,7 +414,7 @@ end """ proposal(c::BGPChain) -Gaussian Transition Kernel centered on current parameter value. +Gaussian Transition Kernel centered on current parameter value. 1. Map all ``k`` parameters into ``\\mu \\in [0,1]^k``. 2. update all parameters by sampling from `MvNormal`, ``N(\\mu,\\sigma)``, where ``sigma`` is `c.sigma` until all params are in ``[0,1]^k`` @@ -411,8 +433,8 @@ function proposal(c::BGPChain) # map into [0,1] # (x-a)/(b-a) = z \in [0,1] - mu01 = mapto_01(mu,lb,ub) - + mu01 = mapto_01(mu,lb,ub) + # Transition Kernel is q(.|theta(t-1)) ~ TruncatedN(theta(t-1), Sigma,lb,ub) # if there is only one batch of params @@ -571,21 +593,26 @@ function computeNextIteration!( algo::MAlgoBGP ) # incrementBGPChainIter!(algo.chains) - # TODO - # this is probably inefficeint - # ideally, would only pmap evaluateObjective, to avoid large data transfers to each worker (now we're transferring each chain back and forth to each worker.) + if get(algo.opts, "parallel", false) + for i in 1:length(algo.chains) + algo.chains[i].iter +=1 #increment iteration on master + end - # if get(algo.opts, "parallel", false) - cs = pmap( x->next_eval(x), algo.chains ) # this does proposal, evaluateObjective, doAcceptRecject - # else + pps = map(proposal, algo.chains) # proposals on master + + evs = pmap(x -> evaluateObjective(algo.m, x), get(algo.opts, "worker_pool", default_worker_pool()), pps) # pmap only the objective function evaluation step + + cs = map(next_acceptreject, algo.chains, evs) # doAcceptRecject, set_eval + + else # for i in algo.chains # @debug(logger," ") # @debug(logger," ") # @debug(logger,"debugging chain id $(i.id)") # next_eval!(i) # end - # cs = map( x->next_eval(x), algo.chains ) # this does proposal, evaluateObjective, doAcceptRecject - # end + cs = map( x->next_eval(x), algo.chains ) # this does proposal, evaluateObjective, doAcceptRecject + end # reorder and insert into algo for i in 1:algo.opts["N"] algo.chains[i] = cs[map(x->getfield(x,:id) == i,cs)][1] @@ -624,7 +651,7 @@ function exchangeMoves!(algo::MAlgoBGP) # algo["N"] exchange moves are proposed props = [(i,j) for i in 1:algo["N"], j in 1:algo["N"] if (i 0 - r[:MSE_SD] = (r[:distance] ./ r[:data_sd] ).^2 + r[:MSE_SD] = (r[:distance] ./ r[:data_sd] ).^2 r[:abs_percent_SD_weighted] = abs.(r[:distance]./ r[:data_sd] ) r[:MSE_SD_1000] = r[:MSE_SD] ./ 1000 r[:abs_percent_SD_weighted_1000] = r[:abs_percent_SD_weighted] ./ 1000 @@ -282,3 +282,8 @@ function show(io::IO,e::Eval) print(io,collect(keys(e.dataMoments))) print(io,"\n===========================\n") end + + +function reInit!(m::MProb, ev::Eval) + reInit!(m, ev.params) +end diff --git a/src/mopt/mprob.jl b/src/mopt/mprob.jl index fca11c3..3c5744e 100644 --- a/src/mopt/mprob.jl +++ b/src/mopt/mprob.jl @@ -19,12 +19,12 @@ of datamoments `moments`. The key idea here is the one of [simulated method of m ```julia-repl pb = Dict(p1" => [0.2,-2,2] , "p2" => [-0.2,-2,2] ) moms = DataFrame(name=["mu2","mu1"],value=[0.0,0.0],weight=rand(2)) -m = MProb() -addSampledParam!(m,pb) -addMoment!(m,moms) +m = MProb() +addSampledParam!(m,pb) +addMoment!(m,moms) SMM.addEvalFunc!(m,SMM.objfunc_norm) ``` - + """ mutable struct MProb @@ -80,7 +80,7 @@ Add parameters to be sampled to an `MProb`. """ function addSampledParam!(m::MProb,name::Any,init::Any, lb::Any, ub::Any) @assert ub>lb - m.initial_value[Symbol(name)] = init + m.initial_value[Symbol(name)] = init m.params_to_sample[Symbol(name)] = Dict( :lb => lb , :ub => ub) return m end @@ -97,6 +97,18 @@ function addSampledParam!(m::MProb,d=OrderedDict{Any,Array{Any,1}}) return m end + +# --------------------- CHANGING STARTING VALUE ---------- +function reInit!(m::MProb, name::Any, init::Any) + m.initial_value[Symbol(name)] = init +end + +function reInit!(m::MProb, d = OrderedDict{Any, Any}) + for (k,v) in d + reInit!(m, Symbol(k), v) + end +end + # -------------------- ADDING MOMENTS -------------------- """ @@ -110,7 +122,7 @@ Add Moments to an `MProb`. Adds a single moment to the mprob. """ function addMoment!(m::MProb,name::Symbol,value,weight) m.moments[Symbol(name)] = Dict( :value => value , :weight => weight ) - return m + return m end """ @@ -132,14 +144,14 @@ function addMoment!(m::MProb,d::Dict) for k in keys(d) addMoment!(m,Symbol(k),d[k][:value],d[k][:weight]) end - return m + return m end function addMoment!(m::MProb,d::DataFrame) for i in 1:nrow(d) m = addMoment!(m,Symbol(d[i,:name]),d[i,:value],d[i,:weight]) end - return m + return m end @@ -148,10 +160,15 @@ function addEvalFunc!(m::MProb,f::Function) m.objfunc = f end +# -------------------- ADDING OBJ FUNCTION OPTIONS -------------------- +function addEvalFuncOpts!(m::MProb,d::Dict) + m.objfunc_opts = d +end + """ evaluateObjective(m::MProb,p::Union{Dict,OrderedDict};noseed=false) -Evaluate the objective function of an [`MProb`](@ref) at a given parameter vector `p`. +Evaluate the objective function of an [`MProb`](@ref) at a given parameter vector `p`. Set `noseed` to `true` if you want to generate new random draws for shocks in the objective function (necessary for estimation of standard errors in [`get_stdErrors`](@ref), for example) """ @@ -162,7 +179,7 @@ function evaluateObjective(m::MProb,p::Union{Dict,OrderedDict};noseed=false) end try # ev = eval(Expr(:call,m.objfunc,ev)) - ev = m.objfunc(ev) + ev = m.objfunc(ev;m.objfunc_opts...) catch ex @warn "caught exception $ex at param $p" ev.status = -2 @@ -179,7 +196,7 @@ function evaluateObjective(m::MProb,ev) # catching errors try # ev = eval(Expr(:call,m.objfunc,ev)) - ev = m.objfunc(ev) + ev = m.objfunc(ev;m.objfunc_opts...) catch ex @warn "caught exception $ex at param $(ev.params)" ev.status = -2 @@ -231,15 +248,39 @@ function mapto_01(p::OrderedDict,lb::Vector{Float64},ub::Vector{Float64}) return (mu .- lb) ./ (ub .- lb) end +function mapto_01(p::Vector{Float64},lb::Vector{Float64},ub::Vector{Float64}) + return (p .- lb) ./ (ub .- lb) +end + +function mapto_01!(p::Vector{Float64},lb::Vector{Float64},ub::Vector{Float64}) + p .-= lb + p ./= (ub .- lb) +end + +function mapto_01!(p::Float64,lb::Float64,ub::Float64) + p -= lb + p /= (ub - lb) +end + """ mapto_ab(p::Vector{Float64},lb::Vector{Float64},ub::Vector{Float64}) - + map param from [0,1] to [a,b] """ function mapto_ab(p::Vector{Float64},lb::Vector{Float64},ub::Vector{Float64}) p .* (ub .- lb) .+ lb end +function mapto_ab!(p::Vector{Float64},lb::Vector{Float64},ub::Vector{Float64}) + p .*= (ub .- lb) + p .+= lb +end + +function mapto_ab!(p::Float64,lb::Float64,ub::Float64) + p *= (ub - lb) + p += lb +end + function show(io::IO,m::MProb) print(io,"MProb Object:\n") @@ -253,7 +294,7 @@ function show(io::IO,m::MProb) # print(io,"\nobjective function call: $(Expr(:call,m.objfunc,m.current_param,m.moments,m.moments_to_use))\n") # if !m.prepared # print(io,"\ncall MoptPrepare(m) to setup cluster\n") - # else + # else # print(io,"\nNumber of chains: $(m.N)\n") # end print(io,"\nEND SHOW MProb\n")