parallel_mc_iid_iv_2.jl

n = 31

ttt()=(randn(100,1), randn(100,1), randn(100,5)) 

function rvproducer()
    for i = 1:n
        produce(ttt())
    end
end

p = Task(rvproducer)

pmap(identity, p)


@everywhere using MinimumDivergence
@everywhere using KNITRO


@everywhere hstderr(mdp::MinimumDivergence.MinimumDivergenceProblems) = vcov!(mdp, :hessian).chol.UL[1]
@everywhere wstderr(mdp::MinimumDivergence.MinimumDivergenceProblems) = vcov!(mdp, :weighted).chol.UL[1]

function dgp(n::Int64         = 100,
             m::Int64         = 5,
             k::Int64         = 1,
             theta0::Float64  = 0.0,
             rho::Float64     = 0.9,
             CP::Int64        = 20)
    ## Generate IV Model - Return a IV object
    tau     = fill(sqrt(CP/(m*n)), m)
    z       = randn(n, m)
    vi      = randn(n, 1)
    eta     = randn(n, 1)
    epsilon = rho*eta+sqrt(1-rho^2)*vi
    x      = z*tau .+ eta    
    y      = x[:,1]*theta0 + epsilon
    ## BLAS.gemm!('N', 'N', 1.0, z, tau, 1.0, eta)
    ## BLAS.gemm!('N', 'N', 1.0, eta, [theta0], 1.0, epsilon)
    #IV(epsilon, eta, z)
    (y, x, z)
end

@everywhere lb = [-20.0]
@everywhere ub = [20.0]
@everywhere θ₀  = [0.0]

@everywhere HellingerDistance()  = CressieRead(-0.5)
@everywhere ContinuousUpdating() = CressieRead(1.0)
@everywhere ModifiedHellingerDistance(ϑ::Real)  = ModifiedCressieRead(-0.5, ϑ)

#function test(args, nsim, seed)
## @everywhere stddiv  = [:ReverseKullbackLeibler, :HellingerDistance, :ContinuousUpdating]
## @everywhere moddiv  = [:ModifiedKullbackLeibler, :ModifiedReverseKullbackLeibler, :ModifiedHellingerDistance]

## @everywhere stdprob = [:RKL, :HD, :CUE]
## @everywhere modprob = [:MKL, :MRKL, :MHD]
## @everywhere probnames = [stdprob, modprob]
## @everywhere pnames = [:KL, probnames]

@everywhere solver = IpoptSolver(print_level=0)

@everywhere function sim(xx)

    iv = IV(xx...)

    #KL = MinDivProb(iv, KullbackLeibler(), solver = KnitroSolver(ms_enable = 1, ms_maxsolves = 15, outlev=0))
    KL = MinDivProb(iv, KullbackLeibler(), solver = solver)
    solve(KL)
    
    coefkl = coef(KL)
    pikl   = getmdweights(KL)
    
    RKL = MinDivProb(iv, ReverseKullbackLeibler(), coefkl, lb, ub, pikl, solver = solver)
    solve(RKL)
    HD = MinDivProb(iv, HellingerDistance(), coefkl,  lb, ub, pikl, solver = solver)
    solve(HD)
    CUE = MinDivProb(iv, ContinuousUpdating(), coefkl, lb, ub, pikl, solver = solver)
    solve(CUE)
    MKL = MinDivProb(iv, ModifiedKullbackLeibler(1.0), coefkl, lb, ub, pikl, solver = solver)
    solve(MKL)
    MRKL = MinDivProb(iv, ModifiedReverseKullbackLeibler(1.0), coefkl, lb, ub, pikl,solver = solver)
    solve(MRKL)
    MHD = MinDivProb(iv, ModifiedHellingerDistance(1.0), coefkl, lb, ub, pikl, solver = solver)
    solve(MHD)

    [coef(KL) coef(RKL) coef(HD) coef(CUE) coef(MKL) coef(MRKL) coef(MHD) wstderr(KL) wstderr(RKL) wstderr(HD) wstderr(CUE) wstderr(MKL) wstderr(MRKL) wstderr(MHD) hstderr(KL) hstderr(RKL) hstderr(HD) hstderr(CUE) hstderr(MKL) hstderr(MRKL) hstderr(MHD) getobjval(KL) getobjval(RKL) getobjval(HD) getobjval(CUE) getobjval(MKL) getobjval(MRKL) getobjval(MHD), MinimumDivergence.status_plain(KL) MinimumDivergence.status_plain(RKL) MinimumDivergence.status_plain(HD) MinimumDivergence.status_plain(CUE) MinimumDivergence.status_plain(MKL) MinimumDivergence.status_plain(MRKL) MinimumDivergence.status_plain(MHD)]
    
end

function replicate(sim::Function, dgp::Function, n::Integer, args;
                   err_retry=true, err_stop=true)
    function rvproducer()
        for i = 1:n
            produce(dgp())
        end
    end
    p = Task(rvproducer)
    pmap(sim, p, err_retry=err_retry, err_stop=err_stop)
end

replicate(dgp::Function, n::Integer, args; err_retry=true, err_stop=true) =
    replicate(identity, dgp, n, args, err_retry=err_retry, err_stop=err_stop)


function sendto(p::Int; args...)
    for (nm, val) in args
        @spawnat(p, eval(Main, Expr(:(=), nm, val)))
    end
end


function sendto(ps::Vector{Int}; args...)
    for p in ps
        sendto(p; args...)
    end
end



## The results are saved in an array (nsim, , ngdp)
## The arguments are in an array (6,

nsim = 6


@everywhere list_of_args = [(n, m, 1, 0.0, rho, CP) for m = (3, 10, 15), rho = (0.3, 0.5, 0.9), CP = (10, 20, 35), n = (200)]

@everywhere list_of_args = reshape(list_of_args, *(size(list_of_args)...), 1)

mc_results = Array(Float64, nsim, 35, length(list_of_args));

for j = 1:length(list_of_args)
    args = list_of_args[j]
    sendto(workers(), args = args) 
    mc_results[:, :, j] = vcat(replicate(sim, dgp, nsim, args, err_stop=false, err_retry=false)...)
end



## ttt()=(randn(100,1), randn(100,1), randn(100,5)) 

## function rvproducer()
##     for i = 1:n
##         produce(ttt())
##     end
## end

## p = Task(rvproducer)

## pmap(sim, p)