text:
const fluxTypes = [:in, :out, :within]

function generateNodeFluxReport( mpSim::MPsim, timeGrid::Vector{Float64},
    fluxType::Symbol, node::String, simNr::Int )

    if lowercase(node) ∈ ["active", ""]
        return generatePopFluxReport( mpSim, timeGrid, fluxType, simNr )
    elseif haskey( mpSim.baseNodeList, node )
        return generateBaseNodeFluxReport( mpSim, timeGrid, fluxType, node,
            simNr )
    else
        return generateCompoundNodeFluxReport( mpSim, timeGrid, fluxType, node,
            simNr )
    end  # if lowercase( node ) ∈ ["active", ""]

end  # generateNodeFluxReport( mpSim, timeGrid, fluxType, node )

function generatePopFluxReport( mpSim::MPsim, timeGrid::Vector{Float64},
    fluxType::Symbol, simNr::Int )

    queryPartCmd = string( "SELECT *, count( `", mpSim.idKey,
        "` ) counts FROM `", mpSim.transDBname, simNr == 0 ? "" : simNr,
        "` WHERE" )

    # Construct (part of) the query.
    if fluxType === :in
        queryPartCmd = string( queryPartCmd,
            "\n    ", mpSim.sNode, " IS NULL AND",
            "\n    ", mpSim.tNode, " IS NOT NULL AND",
            "\n    ", mpSim.tNode, " IS NOT 'active'" )
    elseif fluxType === :out
        queryPartCmd = string( queryPartCmd,
            "\n    ", mpSim.sNode, " IS NOT NULL AND",
            "\n    ", mpSim.tNode, " IS NULL AND",
            "\n    ", mpSim.sNode, " IS NOT 'active'" )
    else
        queryPartCmd = string( queryPartCmd,
            "\n    ", mpSim.sNode, " IS NOT NULL AND",
            "\n    ", mpSim.tNode, " IS NOT NULL AND ",
            "\n    ", mpSim.sNode, " IS NOT ", mpSim.tNode )
    end  # if fluxType === :in

    countDict = determineTransitionTypes( mpSim, queryPartCmd,
        length( timeGrid ) )
    performFluxCounts( mpSim, queryPartCmd, timeGrid, countDict )
    results, names = putFluxResultsInArray( mpSim, countDict,
        length( timeGrid ) )

    # Final touches.
    if fluxType === :in
        names = vcat( "external => active", names )
    elseif fluxType === :out
        names = vcat( "active => external", names )
    else
        names = vcat( "within active", names )
    end  # if fluxType === :in

    names = vcat( "timeStart", "timePoint", names )

    return DataFrame( hcat( vcat( 0.0, timeGrid[1:(end - 1)] ), timeGrid,
        results ), Symbol.( names ) )

end  # generatePopFluxReport( mpSim, timeGrid, fluxType, simNr )

function generateBaseNodeFluxReport( mpSim::MPsim, timeGrid::Vector{Float64},
    fluxType::Symbol, node::String, simNr::Int )

    queryPartCmd = string( "SELECT *, count( `", mpSim.idKey,
        "` ) counts FROM `", mpSim.transDBname, simNr == 0 ? "" : simNr,
        "` WHERE" )

    # Construct (part of) the query.
    if fluxType === :in
        queryPartCmd = string( queryPartCmd,
            "\n    ", mpSim.sNode, " IS NOT '", node, "' AND",
            "\n    ", mpSim.tNode, " IS '", node, "'" )
    else
        queryPartCmd = string( queryPartCmd,
            "\n    ", mpSim.sNode, " IS '", node, "' AND",
            "\n    ", mpSim.tNode, " IS NOT '", node, "'" )
    end  # if fluxType === :in

    countDict = determineTransitionTypes( mpSim, queryPartCmd,
        length( timeGrid ) )
    performFluxCounts( mpSim, queryPartCmd, timeGrid, countDict )
    results, names = putFluxResultsInArray( mpSim, countDict,
        length( timeGrid ) )

    # Final touches.
    if fluxType === :in
        names = vcat( string( "other => ", node ), names )
    else
        names = vcat( string( node, " => other" ), names )
    end  # if fluxType === :in

    names = vcat( "timeStart", "timePoint", names )

    return DataFrame( hcat( vcat( 0.0, timeGrid[1:(end - 1)] ), timeGrid,
        results ), Symbol.( names ) )

end  # generateBaseNodeFluxReport( mpSim, timeGrid, fluxType, node, simNr )

function generateCompoundNodeFluxReport( mpSim::MPsim,
    timeGrid::Vector{Float64}, fluxType::Symbol, node::String, simNr::Int )

    queryPartCmd = string( "SELECT *, count( `", mpSim.idKey,
        "` ) counts FROM `", mpSim.transDBname, simNr == 0 ? "" : simNr,
        "` WHERE" )

    compoundNode = mpSim.compoundNodeList[node]
    nodeList = compoundNode.baseNodeList
    names = Vector{String}()
    results = zeros( Int, length( timeGrid ) )

    if !isempty( nodeList )
        if fluxType === :in
            # The clause "OR ", mpSim.sNode, " IS NULL" is needed because the SQL statement "field NOT IN (values)" gives FALSE for missing values.
            queryPartCmd = string( queryPartCmd,
                "\n    (", mpSim.sNode, " NOT IN ('", join( nodeList, "', '" ),
                "') OR ", mpSim.sNode, " IS NULL) AND",
                "\n    ", mpSim.tNode, " IN ('", join( nodeList, "', '" ),
                "')" )
        elseif fluxType === :out
            queryPartCmd = string( queryPartCmd,
                "\n    ", mpSim.sNode, " IN ('", join( nodeList, "', '" ),
                "') AND",
                "\n    (", mpSim.tNode, " NOT IN ('", join( nodeList, "', '" ),
                "') OR ", mpSim.tNode, " IS NULL)" )
        else
            queryPartCmd = string( queryPartCmd,
                "\n    ", mpSim.sNode, " IN ('", join( nodeList, "', '" ),
                "') AND",
                "\n    ", mpSim.tNode, " IN ('", join( nodeList, "', '" ),
                "') AND",
                "\n    ", mpSim.sNode, " IS NOT ", mpSim.tNode )
        end  # if fluxType === :in

        countDict = determineTransitionTypes( mpSim, queryPartCmd,
            length( timeGrid ) )
        performFluxCounts( mpSim, queryPartCmd, timeGrid, countDict )
        results, names = putFluxResultsInArray( mpSim, countDict,
            length( timeGrid ) )
    end  # if !isempty( nodeList )

    # Final touches.
    if fluxType === :in
        names = vcat( string( "other => ", node ), names )
    elseif fluxType === :out
        names = vcat( string( node, " => other" ), names )
    else
        names = vcat( string( "within ", node ), names )
    end  # if fluxType === :in

    names = vcat( "timeStart", "timePoint", names )

    return DataFrame( hcat( vcat( 0.0, timeGrid[1:(end - 1)] ), timeGrid,
        results ), Symbol.( names ) )

end  # generateCompoundNodeFluxReport( mpSim, timeGrid, fluxType, node, simNr )

function determineTransitionTypes( mpSim::MPsim, queryCmd::String,
    nTimePoints::Int )

    # Determine transition types.
    result = DataFrame( DBInterface.execute( mpSim.simDB, string( queryCmd,
        "\n   GROUP BY transition, ", mpSim.sNode, ", ", mpSim.tNode ) ) )
    countDict = Dict{Tuple, Vector{Int}}()

    if isempty( result )
        return countDict
    end  # if isempty( result )

    for ii in eachindex( result[:, :transition] )
        countDict[(result[ii, :transition],
            result[ii, Symbol( mpSim.sNode )],
            result[ii, Symbol( mpSim.tNode )])] = zeros( Int, nTimePoints )
    end  # for ii in eachindex( result[:transition] )

    return countDict

end  # determineTransitionTypes( mpSim, queryCmd, nTimePoints )

function performFluxCounts( mpSim::MPsim, queryPartCmd::String,
    timeGrid::Vector{Float64}, countDict::Dict{Tuple,Vector{Int}} )

    queryPartCmd = string( queryPartCmd, " AND\n    " )

    for ii in eachindex( timeGrid )
        queryCmd = string( queryPartCmd, generateTimeFork( ii == 1 ?
            timeGrid[1] : timeGrid[ii - 1], timeGrid[ii] ) )

        queryCmd = string( queryCmd,
            "\n   GROUP BY transition, ", mpSim.sNode, ", ", mpSim.tNode )
        result = nothing
        try
            result = DataFrame( DBInterface.execute( mpSim.simDB, queryCmd ) )
        catch err
            println(queryCmd)
            rethrow(err)
        end

        nrecords = size( result, 1 )

        for jj in 1:nrecords
            countDict[(result[jj, :transition],
                result[jj, Symbol( mpSim.sNode )],
                result[jj, Symbol( mpSim.tNode )])][ii] =
                result[jj, :counts]
        end  # for jj in eachindex( result[:transition] )
    end  # for ii in eachindex( timeGrid )

end  # performFluxCounts( mpSim, queryPartCmd, timeGrid, countDict )

function putFluxResultsInArray( mpSim::MPsim,
    countDict::Dict{Tuple,Vector{Int}}, nTimePoints::Int )

    dictNames = collect( keys( countDict ) )
    results = zeros( Int, nTimePoints, length( dictNames ) )
    names = Vector{String}( undef, length( dictNames ) )

    for ii in eachindex( dictNames )
        results[:, ii] = countDict[dictNames[ii]]
        transitionName = dictNames[ii][1]
        sourceNode = dictNames[ii][2] isa Missing ? "external" :
            dictNames[ii][2]
        targetNode = dictNames[ii][3] isa Missing ? "external" :
            dictNames[ii][3]
        names[ii] = string( transitionName, ": ", sourceNode, " => ",
            targetNode )
    end  # for ii in eachindex( dictNames )

    results = hcat( sum( results, dims = 2 ), results )
    return results, names

end  # putFluxResultsInArray( mpSim, countDict, nTimePoints )


text:
__precompile__(true)

using LinearAlgebra

module ExtendedKronigPennyMatrix
using Unitful

include("basic.jl")

export
   Alternates,
   get_E10,
   get_potential,
   update!

include("KronigPenny.jl")

export
   KronigPennyPotential, KronigPennyModel,
   constructMatrix

end

ext: jl


text:
function cleanup()
    try
        if isdir("PlutoNotebooks/180--Working_With_Genie_Apps/MyGenieApp")
            rm("PlutoNotebooks/180--Working_With_Genie_Apps/MyGenieApp", recursive=true)
        end

        if isdir("PlutoNotebooks/41--Developing_MVC_Web_Apps/Watchtonight/")
            rm("PlutoNotebooks/41--Developing_MVC_Web_Apps/Watchtonight/", recursive=true)
        end

        if isdir("PlutoNotebooks/9--Publishing_Your_Julia_Code_Online_With_Genie_Apps/MyGenieApp")
            rm("PlutoNotebooks/9--Publishing_Your_Julia_Code_Online_With_Genie_Apps/MyGenieApp", recursive=true)
        end

        if isdir("PlutoNotebooks/4--Developing_Web_Services/MyGenieApp")
            rm("PlutoNotebooks/4--Developing_Web_Services/MyGenieApp", recursive=true)
        end
    catch
        print("cleanup failed")
    end
end

cleanup()


text:
include("context.jl")
include("attributes.jl")
include("misc.jl")


text:
# ---
# title: 313. Super Ugly Number
# id: problem313
# author: Indigo
# date: 2021-06-10
# difficulty: Medium
# categories: Math, Heap
# link: <https://leetcode.com/problems/super-ugly-number/description/>
# hidden: true
# ---
#
# Write a program to find the `nth` super ugly number.
#
# Super ugly numbers are positive numbers whose all prime factors are in the
# given prime list `primes` of size `k`.
#
# **Example:**
#
#
#
#     Input: n = 12, primes = [2,7,13,19]
#     Output: 32
#     Explanation:[1,2,4,7,8,13,14,16,19,26,28,32] is the sequence of the first 12
#                  super ugly numbers given primes = [2,7,13,19] of size 4.
#
# **Note:**
#
#   * `1` is a super ugly number for any given `primes`.
#   * The given numbers in `primes` are in ascending order.
#   * 0 < `k` ≤ 100, 0 < `n` ≤ 106, 0 < `primes[i]` < 1000.
#   * The nth super ugly number is guaranteed to fit in a 32-bit signed integer.
#
#
## @lc code=start
using LeetCode

function nth_super_ugly_number(n::Int, primes::Vector{Int})
    dp = fill(1, n)
    len = length(primes)
    cur_min = primes[:]
    ptrs = fill(1, len)

    @inbounds for i in 2:n
        minn = minimum(cur_min)
        dp[i] = minn
        @simd for pidx in 1:len
            if minn == cur_min[pidx]
                ptrs[pidx] += 1
                cur_min[pidx] = dp[ptrs[pidx]] * primes[pidx]
            end
        end
    end
    return dp[end]
end
## @lc code=end


text:
using Random, Econometrics
"""
nonparametric regression and quantile regression
* uses pre-generated kernel weights
* can be local constant, local linear, or local quadratic

y is n X p

x is n X dimx

xeval is neval X dimx

weights is n X neval (different weights for each eval. point)

weights should sum to one by columns, they are typical
nonparametric weights from a kernel.

keyword order=0,1 or 2 for local consant, local linear, or local quadratic
(default local linear)

keyword do_median (default false)
keyword do_ci to compute .05 and 0.95 quantiles

execute npreg() for an example
"""
function npreg()
    println("npreg(), with no arguments, runs a simple example")
    println("execute edit(npreg,()) to see the code")
    #using Plots
    k = 3 # number of regressors
    Random.seed!(1) # set seed to enable testing
    n = 10000
    bandwidth = 0.25*n^(-1.0/(4 + k))
    neval = 100
    x = rand(n,k)*pi*2.0
    xeval = [pi*ones(neval,k-1) range(pi/2., stop=pi*1.5, length=neval)]
    y = cos.(sum(x,dims=2))
    y = y + 0.1*randn(size(y))
    ytrue = cos.(sum(xeval,dims=2))
    weights = kernelweights(x, xeval, bandwidth, true, "knngaussian", 200)
    yhat, y50, y05, y95 = npreg(y, x, xeval, weights, order=1, do_median=true, do_ci=true)
    println("true, mean, median, q05, q95")
    prettyprint([ytrue yhat y50 y05 y95])
    return yhat[1,1] # for testing
end

#=
# to do testing, should use nonrandom data
function npreg(x::Int64)
    #println("npreg(), with a single integer argument, runs a test, using the npreg() example")
    #using Plots
    k = 3 # number of regressors
    Random.seed!(1) # set seed to enable testing
    n = 10000
    bandwidth = 0.25*n^(-1.0/(4 + k))
    neval = 100
    x = rand(n,k)*pi*2.0
    xeval = [pi*ones(neval,k-1) range(pi/2., stop=pi*1.5, length=neval)]
    y = cos.(sum(x,dims=2))
    y = y + 0.1*randn(size(y))
    ytrue = cos.(sum(xeval,dims=2))
    weights = kernelweights(x, xeval, bandwidth, true, "knngaussian", 200)
    yhat, y50, y05, y95 = npreg(y, x, xeval, weights, order=1, do_median=true, do_ci=true)
    return yhat[1,1] # for testing
end
=#

function npreg(y, x, xeval, weights; order=1, do_median=false, do_ci=false)
    weights = sqrt.(weights)
    neval, dimx = size(xeval)
    n, dimy = size(y)
    # local constant
    if order==0
        X = ones(n,1)
        Xeval = ones(neval,1)
    elseif order==1
    # local linear
        X = [ones(n,1) x]
        Xeval = [ones(neval,1) xeval]
    else
    # local quadratic
        stacked = [x; xeval]
        # cross products
        CP = zeros(n+neval, Int((dimx-1)*dimx/2))
        cpind = 0
        for i = 1:dimx-1
            for j = (i+1):dimx
                cpind += 1
                CP[:,cpind] = stacked[:,[i]].*stacked[:,[j]]
            end
        end
        ZZ = [ones(n+neval,1) stacked CP]
        X = view(ZZ,1:n,:)
        Xeval = view(ZZ,(n+1):n+neval,:)
    end
    # do the fit
    yhat = zeros(neval, dimy)
    y50 = nothing
    y05 = nothing
    y95 = nothing
    if do_median y50 = zeros(neval, dimy) end
    if do_ci
        y05 = zeros(neval, dimy)
        y95 = zeros(neval, dimy)
    end
    for i = 1:neval
        WX = weights[:,i] .* X
        Wy = weights[:,i] .* y
        b = WX\Wy
        yhat[i,:] = Xeval[i:i,:]*b
        if do_median
            for j = 1:dimy
                y50[i,j] = (Xeval[i:i,:]*qreg_coef(Wy[:,j], WX, 0.5))[1,1]
            end
        end
        if do_ci
            for j = 1:dimy
                y05[i,j] = (Xeval[i:i,:]*qreg_coef(Wy[:,j], WX, 0.05))[1,1]
                y95[i,j] = (Xeval[i:i,:]*qreg_coef(Wy[:,j], WX, 0.95))[1,1]
            end
        end
    end
    if do_median==false
        return yhat
    elseif do_ci==false
        return yhat, y50
    else
        return yhat, y50, y05, y95
    end
end

max_issues_repo_name: DawievLill/Econometrics

max_forks_repo_name: DawievLill/Econometrics


text:
module TimeSeriesA4P2

#Use this to print
#=
weave(Pkg.dir("$(pwd())\\FamaMacbeth.jl"),
  informat="script",
  out_path = "$(pwd())\\FamaMacbethOUT.html",
  doctype = "md2html")
=#

if pwd() ∉ LOAD_PATH
  push!(LOAD_PATH,pwd())
end

using DataFrames, Distributions, StatsBase, GZip, JLD, ForwardDiff,
  Gadfly, JuMP, NLopt, HCubature, StaticArrays, CTNumerical,
  Measures, Formatting

importall CTIO, CTReg, CTStat

const DATA_PATH = pwd() * "\\data" #path where we store the data
const OUTPUT_PATH = pwd() * "\\output"
const YAHOO_DATE_FORMAT = "m/d/yyyy"
const FRED_DATE_FORMAT = "m/d/yyyy"
const PPI_DATE_FORMAT = "m/d/yyyy"
const SP500_NAME = "GSPC"
const DBV_NAME = "DBV"
const FRED_NAME = "CPI"
const PPI_NAME = "PPI"

const FOOTER_NAME = "footer.tex" #these are just headers for the tex tables
const HEADER_NAME = "header.tex"
const DECIMALS = 2
const NUM_FORMAT = "{:.$(DECIMALS)f}"
const PlotContainer = Union{Plot,Gadfly.Compose.Context,Vector{Plot}} #holds graphs
const OVERRIDE_STAR_STRINGS = ["\\ostar","\\dstar","\\tstar"]
const COLORS = ["DodgerBlue", "OrangeRed", "Green",
  "Brown", "Black", "BlueViolet"]

n2s(x::T where T<:Real) = num2Str(x, 4, Ints=true)

#preprocess PPI data
function preProcessPPI()::Void

  ppiDF::DataFrame = readtable("$DATA_PATH\\$PPI_NAME.csv")

  #read and pre-process the data
  dateFormat::String = PPI_DATE_FORMAT
  if :Date ∈ names(ppiDF)
    rename!(ppiDF, :Date, :DATE)
  elseif :_Date ∈ names(ppiDF)
    rename!(ppiDF, :_Date, :DATE)
  elseif :_DATE ∈ names(ppiDF)
    rename!(ppiDF, :_DATE, :DATE)
  end

  ppiDF[:,:DATE] = ((s::String)->Date(s, dateFormat)::Date).(ppiDF[:,:DATE])::DataVector{Date}
  sort!(ppiDF, cols = [:DATE])

  #difference

  ppiDF[:,:Δppi] = similar(ppiDF[:,:ppi])
  ppiDF[:,:Δppi] .= NA
  ppiDF[2:end,:Δppi] = ppiDF[2:end,:ppi] .- ppiDF[1:(end-1),:ppi]

  #log and log difference
  ppiDF[:,:lppi] = log.(ppiDF[:,:ppi])
  ppiDF[:,:Δlppi] = similar(ppiDF[:,:lppi])
  ppiDF[:,:Δlppi] .= NA
  ppiDF[2:end,:Δlppi] = ppiDF[2:end,:lppi] .- ppiDF[1:(end-1),:lppi]

  stream::IOStream = open("$DATA_PATH\\$PPI_NAME.jls", "w")
  serialize(stream, ppiDF)
  close(stream)

  return nothing
end

#to get the df
function getDF(pathStr::String, refreshData::Bool, preProcessFunc::Function)::DataFrame

  #re-process the data if needed or desired
  if !isfile(pathStr) || refreshData
    preProcessFunc()
  end

  #read the binary file
  stream::IOStream = open(pathStr)
  DF::DataFrame = deserialize(stream)
  close(stream)

  return DF
end

function A4P2RMS(Θ::ARMA, Y::Vector{Float64}, YSample::Vector{Float64};
    YHat::Vector{Float64}=similar(Y))

  TSample::Int = length(YSample)
  TOOS::Int = length(Y) - TSample
  ARMAFillYHat!(Θ, Y, YHat)

  return mean((abs2).(Y[(TSample+1):end] .- YHat[(TSample+1):end]))
end

#displays information about the estimated model
function A4P2EstimateInfo(ΘEst::Vector{ARMAEstimate}, ΘEstSample::Vector{ARMAEstimate},
    Y::Vector{Float64}, YSample::Vector{Float64})

  n::Int = length(Y)

  for i::Int ∈ 1:length(ΘEst)
    k::Int = ΘEst[i].Θ.p+ΘEst[i].Θ.q+2
    l::Float64 = ARMALogLike(ΘEst[i].Θ,Y)

    println("""\n*************************************************
      Info for ARMA($(ΘEst[i].Θ.p),$(ΘEst[i].Θ.q))""")
    println("Coefficients: ")
    print("c=$(n2s(ΘEst[i].Θ.c)) ($(n2s(ΘEst[i].cError)))\n")
    for j::Int ∈ 1:ΘEst[i].Θ.p
      print("ϕ$j=$(n2s(ΘEst[i].Θ.ϕ[j])) ($(n2s(ΘEst[i].ϕErrors[j]))), ")
    end
    print("\n")
    for j::Int ∈ 1:ΘEst[i].Θ.q
      print("ϑ$j=$(n2s(ΘEst[i].Θ.ϑ[j])) ($(n2s(ΘEst[i].ϑErrors[j]))), ")
    end
    print("\nσ²x1000=$(n2s(ΘEst[i].Θ.σ²*1000.)) ($(n2s(ΘEst[i].σ²Error*1000.)))\n")
    println("MSPE (x1000): $(n2s(1000.*A4P2RMS(ΘEstSample[i].Θ, Y, YSample)))")
    println("AIC: $(n2s(2*k-2*l))")
    println("BIC: $(n2s(log(n)*k-2*l))")

    for j::Int ∈ 8:12
      println("Q$j: $(n2s(ljungBox(Y .- ΘEst[i].YHat, j:j)))")
    end
  end

  return nothing

end

#problem specific graphs
function visualizeA4P2(ΘEstimate::Vector{ARMAEstimate}, ppiDF::DataFrame)::Void

  ACFRange::Range = 1:24
  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  Gadfly.push_theme(:default)
  outDF = deepcopy(ppiDF)
  completecases!(outDF)

  push!(plotNames, "Summary Graphs")
  push!(plots, Vector{Plot}())
  for series::Symbol ∈ [:ppi, :Δppi, :lppi, :Δlppi]
    push!(plots[end],
      plot(outDF, x=:DATE, y=series,  Geom.line, Geom.point,
        Guide.xlabel("Date"),Guide.ylabel("$series"),
        Guide.title("Summary $series")))
  end
  plots[end] = gridstack(reshape(plots[end], :, 2))

  push!(plotNames, "ACF Graphs")
  push!(plots, Vector{Plot}())
  for series::Symbol ∈ [:ppi, :Δppi, :lppi, :Δlppi]
    push!(plots[end],
      plot(x=collect(ACFRange), y=autoρ(Vector{Float64}(outDF[:,series]),ACFRange),
        Geom.line, Geom.point,
        Guide.xlabel("Date"),Guide.ylabel("$series ρ"),
        Guide.title("ACF $series")))
  end
  plots[end] = gridstack(reshape(plots[end], :, 2))

  push!(plotNames, "PACF Graphs")
  push!(plots, Vector{Plot}())
  for series::Symbol ∈ [:ppi, :Δppi, :lppi, :Δlppi]
    push!(plots[end],
      plot(x=collect(ACFRange), y=pAutoρ(Vector{Float64}(outDF[:,series]),ACFRange),
        Geom.line, Geom.point,
        Guide.xlabel("Date"),Guide.ylabel("$series ∂ρ"),
        Guide.title("PACF $series")))
  end
  plots[end] = gridstack(reshape(plots[end], :, 2))

  #write the graphs
  for i::Int ∈ 1:length(plots)
    draw(SVG("$OUTPUT_PATH\\$(plotNames[i]).svg", 9inch, 7inch), plots[i])
    #display(plots[i])
  end

  return nothing
end

#script for the problem
function A4P2Script(;refreshData::Bool = true)::Void

  targetSeries::Symbol = :Δlppi

  #select the models
  ARMA22Act::ARMA = ARMA(c=.007, ϕ=[-0.063, .739], ϑ=[.353,-.523],σ²=.001)
  models::Vector{ARMA} = [#=ARMA(1,1), ARMA(2,2), ARMA(2,3), =#ARMA22Act, ARMA(3,4)]
  modelsTest::Vector{ARMA} = deepcopy(models)
  ppiDF::DataFrame = getDF("$DATA_PATH\\$PPI_NAME.jls", refreshData, preProcessPPI)
  #visualizeA4P2(models, ppiDF)
  Y::Vector{Float64} = Vector{Float64}(dropna(ppiDF[:,:Δlppi]))
  YSample::Vector{Float64} =
    Vector{Float64}(dropna(ppiDF[(Dates.year).(ppiDF[:, :DATE]) .≤ 2005,:Δlppi]))
  #get the data in an array
  modelEstimates::Vector{ARMAEstimate} = ((Θ::ARMA)->ARMAEstimate(Θ,Y, verbose=true)).(models)
  modelSampleEstimates::Vector{ARMAEstimate} = ((Θ::ARMA)->ARMAEstimate(Θ,YSample)).(models)

  A4P2EstimateInfo(modelEstimates, modelSampleEstimates, Y, YSample)

  #visualizeA4P2!()

  return nothing
end

@time begin
  A4P2Script(refreshData = true)
end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiCoreV1ConfigMapList <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    items::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ConfigMap} } # spec name: items
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ListMeta } # spec name: metadata

    function IoK8sApiCoreV1ConfigMapList(;apiVersion=nothing, items=nothing, kind=nothing, metadata=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1ConfigMapList, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiCoreV1ConfigMapList, Symbol("items"), items)
        setfield!(o, Symbol("items"), items)
        validate_property(IoK8sApiCoreV1ConfigMapList, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiCoreV1ConfigMapList, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        o
    end
end # type IoK8sApiCoreV1ConfigMapList

const _property_map_IoK8sApiCoreV1ConfigMapList = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("items")=>Symbol("items"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"))
const _property_types_IoK8sApiCoreV1ConfigMapList = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("items")=>"Vector{IoK8sApiCoreV1ConfigMap}", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ListMeta")
Base.propertynames(::Type{ IoK8sApiCoreV1ConfigMapList }) = collect(keys(_property_map_IoK8sApiCoreV1ConfigMapList))
Swagger.property_type(::Type{ IoK8sApiCoreV1ConfigMapList }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1ConfigMapList[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1ConfigMapList }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1ConfigMapList[property_name]

function check_required(o::IoK8sApiCoreV1ConfigMapList)
    (getproperty(o, Symbol("items")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1ConfigMapList }, name::Symbol, val)
end


text:
using Documenter
using PublicAPI

makedocs(
    sitename = "PublicAPI",
    format = Documenter.HTML(),
    modules = [PublicAPI],
    doctest = false,  # tested via test/runtests.jl
    checkdocs = :exports,  # ignore complains about internal docs
    strict = lowercase(get(ENV, "CI", "false")) == "true",
    # See: https://juliadocs.github.io/Documenter.jl/stable/lib/public/#Documenter.makedocs
)

deploydocs(
    repo = "github.com/JuliaExperiments/PublicAPI.jl",
    devbranch = "main",  # https://github.com/JuliaDocs/Documenter.jl/issues/1443
    push_preview = true,
    # See: https://juliadocs.github.io/Documenter.jl/stable/lib/public/#Documenter.deploydocs
)


text:
using Random
using ProgressMeter
include("./utils.jl")

using Distributed
addpross(4)

@everywhere begin
    include("./mcts.jl")
end

@everywhere function playonce()
    game = Game()
    randomPlace(game, 2)
    playerX = MCTS()
    playerO = MCTS()
    chunk = Array{Array{Float32, 3}, Array{Float32, 1}, Integer}[]
    while iszero(game.winner)
        move, p = playerX(game, Clients)
        data = getAugmentation(game, p)
        append!(chunk, data)
        place(game, move)

        !iszero(game.winner) && break
        move, p = playerO(game, Clients)
        data = getAugmentation(game, p)
        append!(chunk, data)
        place(game, move)
    end
    @async closeClients(playerX.Clients)
    @async closeClients(playerO.Clients)
    for e in chunk
        @inbounds e[3] .*= game.winner
    end
    return chunk
end

function playNgames(N::Integer=64, bs::Integer=16)
    results = @showprogress pmap(1:N, batch_size=bs) do x
        playonce()
    end
    save_data(shuffle!(reduce(append!, results)))
    nothing
end


text:
module Etcd

using Requests
using Memento

import HttpCommon: Response

immutable Client
    host::String
    port::Int
    version::String
end

include("constants.jl")
include("requests.jl")
include("api.jl")
include("utils.jl")

"""
    connect(host="localhost", port=2379, version="v2")

Creates an `Etcd.Client` which can then be used for making requests to an etcd cluster.
"""
connect(host::String="localhost", port::Int=2379, version="v2") = Client(host, port, version)

export
    # methods
    machines,
    stats,
    members,
    leader,
    set,
    create,
    update,
    setdir,
    createdir,
    updatedir,
    delete,
    deletedir,
    exists,
    cas,
    cad,
    watch,
    watchloop,

    # Types
    HTTPError,
    EtcdError

end # module


text:
using Bokeh
autoopen(true)

plotfile("legend.html")

y1 = [-10:10].^2
y2 = [-10:10].^3
plot([y1 y2], legends=["x^2", "x^3"], legendsgo=:bottom_right)


text:
# This file contains a list of the currently available covariance functions

import Base.show

abstract Kernel

"""
Data to be used with a kernel object to
calculate a covariance matrix, which is independent of kernel hyperparameters.

# See also
`EmptyData`
"""
abstract KernelData

"""
Default KernelData type which is empty.
"""
type EmptyData <: KernelData
end

KernelData(k::Kernel, X::Matrix{Float64}) = EmptyData()

"""
# Description
Constructs covariance matrix from kernel and input observations

# Arguments
# `k::Kernel`: kernel for calculating covariance between pairs of points
* `X::Matrix{Float64}`: matrix of observations (each column is an observation)
# `Y::Matrix{Float64}`: another matrix of observations

# Return
`Σ::Matrix{Float64}`: covariance matrix where `Σ[i,j]` is the covariance of the Gaussian process between points `X[:,i]` and `Y[:,j]`.
"""
function cov(k::Kernel, X::Matrix{Float64}, Y::Matrix{Float64})
    d(x,y) = cov(k, x, y)
    return map_column_pairs(d, X, Y)
end

"""
# Description
Constructs covariance matrix from kernel and kernel data

# Arguments
* `k::Kernel`: kernel for calculating covariance between pairs of points
* `X::Matrix{Float64}`: matrix of input observations (each column is an observation)
* `data::KernelData`: data, constructed from input observations, used for calculating covariance matrix

# Return
`Σ::Matrix{Float64}`: covariance matrix where `Σ[i,j]` is the covariance of the Gaussian process between points `X[:,i]` and `X[:,j]`.

# See also
Kernel, KernelData
"""
function cov(k::Kernel, X::Matrix{Float64}, data::EmptyData)
    d(x,y) = cov(k, x, y)
    return map_column_pairs(d, X)
end

cov(k::Kernel, X::Matrix{Float64}) = cov(k, X, KernelData(k, X))

# Calculates the stack [dk / dθᵢ] of kernel matrix gradients
function grad_stack!(stack::AbstractArray, k::Kernel, X::Matrix{Float64}, data::EmptyData)
    d, nobsv = size(X)
    for j in 1:nobsv, i in 1:nobsv
        @inbounds stack[i,j,:] = grad_kern(k, X[:,i], X[:,j])
    end
    return stack
end

function grad_stack(k::Kernel, X::Matrix{Float64}, data::KernelData)
    n = num_params(k)
    n_obsv = size(X, 2)
    stack = Array(Float64, n_obsv, n_obsv, n)
    grad_stack!(stack, k, X, data)
    return stack
end

function grad_stack!(stack::AbstractArray, k::Kernel, X::Matrix{Float64})
    grad_stack!(stack, k, X, KernelData(k, X))
end

grad_stack(k::Kernel, X::Matrix{Float64}) = grad_stack(k, X, KernelData(k, X))

##############################
# Parameter name definitions #
##############################

# This generates names like [:ll_1, :ll_2, ...] for parameter vectors
get_param_names(n::Int, prefix::Symbol) = [symbol(prefix, :_, i) for i in 1:n]
get_param_names(v::Vector, prefix::Symbol) = get_param_names(length(v), prefix)

# Fallback. Yields names like :Matl2Iso_param_1 => 0.5
# Ideally this is never used, because the names are uninformative.
get_param_names(obj::Union{Kernel, Mean}) =
    get_param_names(num_params(obj),
                    symbol(typeof(obj).name.name, :_param_))

""" `composite_param_names(objects, prefix)`, where `objects` is a
vector of kernels/means, calls `get_param_names` on each object and prefixes the
name returned with `prefix` + object #. Eg.

    get_param_names(ProdKernel(Mat(1/2, 1/2, 1/2), SEArd([0.0, 1.0],0.0)))

yields

    :pk1_ll
    :pk1_lσ
    :pk2_ll_1
    :pk2_ll_2
    :pk2_lσ
"""
function composite_param_names(objects, prefix)
    p = Symbol[]
    for (i, obj) in enumerate(objects)
        append!(p, [symbol(prefix, i, :_, sym) for sym in get_param_names(obj)])
    end
    p
end

function show(io::IO, k::Kernel, depth::Int = 0)
    pad = repeat(" ", 2*depth)
    print(io, "$(pad)Type: $(typeof(k)), Params: ")
    # params_dict = zip(get_param_names(k), get_params(k))
    # for (k, val) in params_dict
    #     print(io, "$(k)=$(val) ")
    # end
    show(io, get_params(k))
    print(io, "\n")
end

include("stationary.jl")

include("lin.jl")               # Linear covariance function
include("se.jl")                # Squared exponential covariance function
include("rq.jl")                # Rational quadratic covariance function
include("mat.jl")               # Matern covariance function
include("periodic.jl")          # Periodic covariance function
include("poly.jl")              # Polnomial covariance function
include("noise.jl")             # White noise covariance function

# Composite kernels
include("sum_kernel.jl")        # Sum of kernels
include("prod_kernel.jl")       # Product of kernels


text:
function dispatch_callback(f::Function, argc)
    argv = map(1:argc) do i
        unbox(CLRObject(i))
    end
    result = f(argv...)
    return isnothing(result) ? 0 : box(result, 1)
end

function delegate(f::Function, delty::CLRObject)
    context = objectid(f)
    d = CLRBridge.CreateDelegate(gethandle(delty), context)
    CLRBridge.registered_callbacks[context] = argc->dispatch_callback(f, argc)
    return d
end


text:
struct AsynchronousExecution{T <: AbstractFloat,
                            A <: AbstractVector,
                            S <: LQSolver} <: AbstractExecution
    subworkers::Vector{SubWorker{T,A,S}}
    work::Vector{Work}
    progressqueue::ProgressQueue{T}
    x̄::Vector{RunningAverage{A}}
    δ::Vector{RunningAverage{T}}
    decisions::Decisions{A}
    r::IteratedValue{T}
    active_workers::Vector{Future}
    # Bookkeeping
    subobjectives::Vector{A}
    finished::Vector{Int}
    # Parameters
    κ::T

    function AsynchronousExecution(κ::T, ::Type{T}, ::Type{A}, ::Type{S}) where {T <: AbstractFloat, A <: AbstractVector, S <: LQSolver}
        return new{T,A,S}(Vector{SubWorker{T,A,S}}(undef, nworkers()),
                          Vector{Work}(undef,nworkers()),
                          RemoteChannel(() -> Channel{Progress{T}}(4*nworkers())),
                          Vector{RunningAverage{A}}(undef,nworkers()),
                          Vector{RunningAverage{T}}(undef,nworkers()),
                          RemoteChannel(() -> IterationChannel(Dict{Int,A}())),
                          RemoteChannel(() -> IterationChannel(Dict{Int,T}())),
                          Vector{Future}(undef,nworkers()),
                          Vector{A}(),
                          Vector{Int}(),
                          κ)
    end
end

function init_subproblems!(ph::AbstractProgressiveHedgingSolver, subsolver::QPSolver, execution::AsynchronousExecution{T,A,S}) where {T <: AbstractFloat, A <: AbstractVector, S <: LQSolver}
    # Create subproblems on worker processes
    m = ph.stochasticprogram
    @sync begin
        for w in workers()
            execution.subworkers[w-1] = RemoteChannel(() -> Channel{Vector{SubProblem{T,A,S}}}(1), w)
            @async load_worker!(scenarioproblems(m), m, w, execution.subworkers[w-1], subsolver)
        end
    end
    @sync begin
        # Continue preparation
        for w in workers()
            execution.work[w-1] = RemoteChannel(() -> Channel{Int}(round(Int,10/execution.κ)), w)
            @async execution.x̄[w-1] = remotecall_fetch((sw, xdim)->begin
                subproblems = fetch(sw)
                if length(subproblems) > 0
                    x̄ = sum([s.π*s.x for s in subproblems])
                    return RemoteChannel(()->RunningAverageChannel(x̄, [s.x for s in subproblems]), myid())
                else
                    return RemoteChannel(()->RunningAverageChannel(zeros(T,xdim), Vector{A}()), myid())
                end
            end, w, execution.subworkers[w-1], decision_length(m))
            @async execution.δ[w-1] = remotecall_fetch((sw, xdim)->RemoteChannel(()->RunningAverageChannel(zero(T), fill(zero(T),length(fetch(sw))))), w, execution.subworkers[w-1], decision_length(m))
            put!(execution.work[w-1], 1)
        end
        # Prepare memory
        push!(execution.subobjectives, zeros(nscenarios(ph)))
        push!(execution.finished, 0)
        log_val = ph.parameters.log
        ph.parameters.log = false
        log!(ph)
        ph.parameters.log = log_val
    end
    update_iterate!(ph)
    # Init δ₂
    @sync begin
        for w in workers()
            @async remotecall_fetch((sw,ξ,δ)->begin
                for (i,s) ∈ enumerate(fetch(sw))
                    take!(δ, i)
                    put!(δ, i, norm(s.x - ξ, 2)^2, s.π)
                end
            end, w, execution.subworkers[w-1], ph.ξ, execution.δ[w-1])
        end
    end
    return ph
end

function iterate!(ph::AbstractProgressiveHedgingSolver, execution::AsynchronousExecution{T}) where T <: AbstractFloat
    wait(execution.progressqueue)
    while isready(execution.progressqueue)
        # Add new cuts from subworkers
        t::Int, i::Int, Q::T = take!(execution.progressqueue)
        if Q == Inf
            @warn "Subproblem $(i) is infeasible, aborting procedure."
            return :Infeasible
        end
        execution.subobjectives[t][i] = Q
        execution.finished[t] += 1
        if execution.finished[t] == nscenarios(ph)
            # Update objective
            ph.Q_history[t] = current_objective_value(ph, execution.subobjectives[t])
            ph.data.Q = ph.Q_history[t]
        end
    end
    # Project and generate new iterate
    t = ph.data.iterations
    if execution.finished[t] >= execution.κ*nscenarios(ph)
        # Get dual gap
        update_dual_gap!(ph)
        # Update progress
        @unpack δ₁, δ₂ = ph.data
        ph.dual_gaps[t] = δ₂
        ph.data.δ = sqrt(δ₁ + δ₂)/(1e-10+norm(ph.ξ,2))
        # Check if optimal
        if check_optimality(ph)
            # Optimal, tell workers to stop
            map((w,aw)->!isready(aw) && put!(w,t), execution.work, execution.active_workers)
            map((w,aw)->!isready(aw) && put!(w,-1), execution.work, execution.active_workers)
            # Final log
            log!(ph)
            return :Optimal
        end
        # Update penalty (if applicable)
        update_penalty!(ph)
        # Update iterate
        update_iterate!(ph)
        # Send new work to workers
        put!(execution.decisions, t+1, ph.ξ)
        put!(execution.r, t+1, penalty(ph))
        map((w,aw)->!isready(aw) && put!(w,t+1), execution.work, execution.active_workers)
        # Prepare memory for next iteration
        push!(execution.subobjectives, zeros(nscenarios(ph)))
        push!(execution.finished, 0)
        # Log progress
        log!(ph)
    end
    # Just return a valid status for this iteration
    return :Valid
end

function init_workers!(ph::AbstractProgressiveHedgingSolver, execution::AsynchronousExecution)
    # Load initial decision
    put!(execution.decisions, 1, ph.ξ)
    put!(execution.r, 1, penalty(ph))
    for w in workers()
        execution.active_workers[w-1] = remotecall(work_on_subproblems!,
                                                   w,
                                                   execution.subworkers[w-1],
                                                   execution.work[w-1],
                                                   execution.progressqueue,
                                                   execution.x̄[w-1],
                                                   execution.δ[w-1],
                                                   execution.decisions,
                                                   execution.r)
    end
    return nothing
end

function close_workers!(ph::AbstractProgressiveHedgingSolver, execution::AsynchronousExecution)
    map(wait, execution.active_workers)
end

function resolve_subproblems!(ph::AbstractProgressiveHedgingSolver, execution::AsynchronousExecution)
    return nothing
end

function update_iterate!(ph::AbstractProgressiveHedgingSolver, execution::AsynchronousExecution)
    ξ_prev = copy(ph.ξ)
    ph.ξ[:] = sum(fetch.(execution.x̄))
    # Update δ₁
    ph.data.δ₁ = norm(ph.ξ-ξ_prev, 2)^2
    return nothing
end

function update_subproblems!(ph::AbstractProgressiveHedgingSolver, execution::AsynchronousExecution)
    return nothing
end

function update_dual_gap!(ph::AbstractProgressiveHedgingSolver, execution::AsynchronousExecution)
    ph.data.δ₂ = sum(fetch.(execution.δ))
    return nothing
end

function calculate_objective_value(ph::AbstractProgressiveHedgingSolver, execution::AsynchronousExecution)
    return calculate_objective_value(ph, execution.subworkers)
end

function
    fill_submodels!(ph::AbstractProgressiveHedgingSolver, scenarioproblems, execution::AsynchronousExecution)
    return fill_submodels!(ph, scenarioproblems, execution.subworkers)
end

# API
# ------------------------------------------------------------
"""
    Asynchronous

...
# Parameters
- `κ::Real = 0.5`: Amount of subproblems, relative to the total number of scenarios, required to generate a new iterate in master procedure.

...
"""
struct Asynchronous <: Execution
    κ::Float64

    function Asynchronous(; κ::AbstractFloat = 0.5)
        return new(κ)
    end
end

function (execution::Asynchronous)(::Type{T}, ::Type{A}, ::Type{S}) where {T <: AbstractFloat, A <: AbstractVector, S <: LQSolver}
    return AsynchronousExecution(execution.κ, T, A, S)
end

function str(::Asynchronous)
    return "asynchronous execution"
end


text:
using XCompiler, Test

@test (@enter_call sin(10); true)


text:
using IcanProj
using JLD2
using PyPlot

opt = get_options()

watershed = opt["stat_sel"][3]

# Load results

file = joinpath(opt["eval_folder"], watershed, "res_1km.jld2")

@load file res_vic

ivar = find(res_vic.var_names .== "swe")

data_swe = res_vic.data_all[:,ivar[1],:]

mean_swe = mean(data_swe,2)

plot(res_vic.time, mean_swe, label = "1km")

file = joinpath(opt["eval_folder"], watershed, "res_50km.jld2")

@load file res_vic

ivar = find(res_vic.var_names .== "swe")

data_swe = res_vic.data_all[:,ivar[1],:]

mean_swe = mean(data_swe,2)

plot(res_vic.time, mean_swe, label = "50km")

ylabel("SWE (mm)")

legend()


text:
using Documenter, SAMtools

makedocs(
    modules = [SAMtools],
    format = :html,
    checkdocs = :exports,
    sitename = "SAMtools.jl",
    pages = Any["index.md"]
)

deploydocs(
    repo = "github.com/aramirezreyes/SAMtools.jl.git",
)


text:
## make gradientinv options same for all soundings
σstart, σ0         = -2, -2
zfixed             = [-1e5]
ρfixed             = [1e12]
zstart             = 0.0
extendfrac         = 1.06
dz                 = 1.5
ρbg                = 10
nlayers            = 50
ntimesperdecade    = 10
nfreqsperdecade    = 5
modelprimary       = false
regtype            = :R1
nstepsmax          = 40
ntries             = 6
target             = nothing
lo                 = -3.
hi                 = 1.
λ²min              = -0.5
λ²max              = 8
λ²frac             = 4
ntestdivsλ²        = 50
αmin               = -4
αmax               = 0
αfrac              = 4
β²                 = 0.1
ntestdivsα         = 32
regularizeupdate   = false
knownvalue         = 0.7
firstvalue         = :last
κ                  = transD_GP.GP.Mat52()
breakonknown       = true
dobo               = false
## plot n random soundings and a background response
using Random
nplot = 2
nplot = min(nplot, length(soundings))
opt = transD_GP.SkyTEM1DInversion.makeoperatorandoptions(
                            soundings,
                            rseed = 2,
                            zfixed = zfixed,
                            ρfixed = ρfixed,
                            zstart = zstart,
                            extendfrac = extendfrac,
                            dz = dz,
                            ρbg = ρbg,
                            nlayers = nlayers,
                            ntimesperdecade = ntimesperdecade,
                            nfreqsperdecade = nfreqsperdecade)

max_stars_repo_name: GeoscienceAustralia/HiQGA

max_issues_repo_name: GeoscienceAustralia/HiQGA


text:
@testset "Escape.run(::HLAModel, data::AbstractHLAData)" begin
    ds = Escape.HLADataset("Test")
    m = Escape.HLAModel{4}()
    data = ds.data[1]

    @test @suppress Escape.run(m, data, iter = 5, warmup = 5, chains = 1) isa
        Escape.HLAModelResult
end

@testset "Escape.reduce_size!(::StanInterface.Stanfit)" begin
    ds = Escape.HLADataset("Test")
    m = Escape.HLAModel{4}()
    data = ds.data[1]

    res = @suppress Escape.run(m, data, iter = 5, warmup = 5, chains = 1,
        keep_all_parameters = true)
    previous_size = Base.summarysize(res.sf)

    @test "aux1_tau.1" in keys(res.sf.result[1])
    Escape.reduce_size!(res.sf)

    new_size = Base.summarysize(res.sf)
    @test !("aux1_tau.1" in keys(res.sf.result[1]))
    @test new_size < previous_size
end

@testset "Escape.run_model(::HLAModel, ::AbstractHLAData)" begin
    alignment_file_2_digits = joinpath(@__DIR__, "data",
        "test_large_annotated_2_digits.fasta")
    tree_file = joinpath(@__DIR__, "data", "phylogeny.tree")

    data = @suppress HLAData(
        alignment_file = alignment_file_2_digits,
        tree_file = tree_file
    )

    @test @suppress Escape.run_model(
        Escape.HLAModel{3}(),
        data, iter = 5, warmup = 5, chains = 1
    ) isa Escape.HLAModelResult

    save_file = tempname() * ".jld2"
    @test !isfile(save_file)
    @suppress Escape.run_model(
        Escape.HLAModel{3}(),
        data, iter = 5, warmup = 5, chains = 1, save_file = save_file
    )
    @test isfile(save_file)
    rm(save_file)
end

@testset "Escape.run_model(::AbstractHLAData)" begin
    alignment_file_2_digits = joinpath(@__DIR__, "data",
    "test_large_annotated_2_digits.fasta")
    tree_file = joinpath(@__DIR__, "data", "phylogeny.tree")

    data = @suppress HLAData(
        alignment_file = alignment_file_2_digits,
        tree_file = tree_file
    )

    @test @suppress Escape.run_model(
        data, iter = 5, warmup = 5, chains = 1
    ) isa Escape.HLAModelResult

    save_file = tempname() * ".jld2"
    @test !isfile(save_file)
    @suppress Escape.run_model(
        data, iter = 5, warmup = 5, chains = 1, save_file = save_file
    )
    @test isfile(save_file)
    rm(save_file)
end


text:
module Anderson
using LinearAlgebra, Parameters

include("qrtools.jl")
include("function_iteration.jl")
include("anderson_lowalloc.jl")
end # module


text:
"""
Check vtree validation, variables(parent) = variables(left) + variables(right)
"""
validate_vtree(n::Vtree)::Bool = validate_vtree(NodeType(n), n)
validate_vtree(::Leaf, n::Vtree) = nothing
function validate_vtree(::Inner, n::Vtree)
    @test num_variables(n) == num_variables(n.left) + num_variables(n.right)
    @test variables(n) == union(variables(n.left), variables(n.right))
    validate_vtree(n.right)
    validate_vtree(n.left)
end


text:
"""
```
PoolModel{T}
```

Implements several model averaging methods to pool
the predictions of different models. Currently, `PoolModel`
only works for pooling two different models, although
it can be extended to more models. Confer with the paper for details.

The available methods are

* dynamic: weights on pooled models evolve over time and are
           treated as a stochastic process.

* static:  weights on pooled models are assumed time-invariant.

* equal:   weights are set to 1/2.

* bma:     weights are updated according to Bayesian model averaging.

The default method is dynamic, which implements the Dynamic Pools method
developed by Del Negro et al. (2016). To choose another method,
use the keyword `weight_type::Symbol`, e.g.

```
pm = PoolModel(weight_type = :static) # creates a static PoolModel
```

The `weight_type` is stored as a setting, so users can retrieve
at any point by using `get_setting`.

### Fields

#### Parameters
* `parameters::Vector{AbstractParameter}`: Vector of all time-invariant hyperparameters
    for the chosen method of model averaging.

* `keys::OrderedDict{Symbol,Int}`: Maps human-readable names for predictive densities
    of pooled models.

#### Inputs to Measurement and Equilibrium Condition Equations

* `model::OrderedDict{Symbol,AbstractDSGEModel}`: Maps name to its underlying model
  object.

#### Model Specifications and Settings

* `spec::String`: The model specification identifier, \"an_schorfheide\", cached
  here for filepath computation.

* `subspec::String`: The model subspecification number, indicating that some
  parameters from the original model spec (\"ss0\") are initialized
  differently. Cached here for filepath computation.

* `settings::Dict{Symbol,Setting}`: Settings/flags that affect computation
  without changing the economic or mathematical setup of the model.

* `test_settings::Dict{Symbol,Setting}`: Settings/flags for testing mode

#### Other Fields

* `rng::MersenneTwister`: Random number generator. Can be is seeded to ensure
  reproducibility in algorithms that involve randomness (such as
  Metropolis-Hastings).

* `testing::Bool`: Indicates whether the model is in testing mode. If `true`,
  settings from `m.test_settings` are used in place of those in `m.settings`.

* `observable_mappings::OrderedDict{Symbol,Observable}`: A dictionary that
  stores data sources, series mnemonics, and transformations to/from model units.
  DSGE.jl will fetch data from the Federal Reserve Bank of
  St. Louis's FRED database; all other data must be downloaded by the
  user. See `load_data` and `Observable` for further details.
"""
mutable struct PoolModel{T} <: AbstractDSGEModel{T}
    parameters::ParameterVector{T}                         # vector of all time-invariant model parameters
    keys::OrderedDict{Symbol,Int}                          # human-readable names for all the model
    observables::OrderedDict{Symbol,Int}                   # Model names to observables (predictive densities)

    spec::String                                           # Model specification number (eg "m990")
    subspec::String                                        # Model subspecification (eg "ss0")
    settings::Dict{Symbol,Setting}                         # Settings/flags for computation
    test_settings::Dict{Symbol,Setting}                    # Settings/flags for testing mode
    rng::MersenneTwister                                   # Random number generator
    testing::Bool                                          # Whether we are in testing mode or not

    observable_mappings::OrderedDict{Symbol, Observable}
end

description(m::PoolModel) = "Julia implementation of prediction pool methods defined in 'Dynamic prediction pools: An investigation of financial frictions and forecasting performance' by Marco Del Negro, Raiden B. Hasegawa, and Frank Schorfheide: PoolModel, $(m.subspec)"

function init_model_indices!(m::PoolModel)
    # Observables
    observables = keys(m.observable_mappings)

    # Make indices
    for (i,k) in enumerate(observables);       m.observables[k]       = i end
end

function PoolModel(subspec::String="ss2";
                   custom_settings::Dict{Symbol,Setting} = Dict{Symbol,Setting}(),
                   testing = false, verbose::Symbol = :low,
                   weight_type::Symbol = :dynamic) where T<:AbstractFloat

    # Model-specific specifications
    spec               = split(basename(@__FILE__),'.')[1]
    subspec            = subspec
    settings           = Dict{Symbol,Setting}()
    test_settings      = Dict{Symbol,Setting}()
    rng                = MersenneTwister(0)

    # initialize empty model
    m = PoolModel{Float64}(
        # model parameters
        Vector{AbstractParameter{Float64}}(), OrderedDict{Symbol,Int}(),

        # Observable indices
        OrderedDict{Symbol,Int}(),

        # settings
        spec,
        subspec,
        settings,
        test_settings,
        rng,
        testing,
        OrderedDict{Symbol,Observable}())

    # Set settings
    model_settings!(m; weight_type = weight_type)
    default_test_settings!(m)
    for custom_setting in values(custom_settings)
        m <= custom_setting
    end

    # Set observable and pseudo-observable transformations
    init_observable_mappings!(m)

    # Initialize parameters
    init_parameters!(m)

    # Initialize model indices and subspec
    init_model_indices!(m)
    init_subspec!(m)

    return m
end

"""
```
init_parameters!(m::PoolModel)
```

Initializes the model's parameters, as well as empty values for the steady-state
parameters (in preparation for `steadystate!(m)` being called to initialize
those).
"""
function init_parameters!(m::PoolModel)
    # Initialize parameters
    weight_type = get_setting(m, :weight_type)
    if weight_type == :dynamic
        m <= parameter(:ρ, 0.8, (1e-5,0.999), (1e-5,0.999), SquareRoot(), Uniform(0.,1.),
                       fixed = false,
                       description="ρ: persistence of AR processing underlying λ.",
                       tex_label="\\rho")
        m <= parameter(:μ, 0., fixed = true,
                       description="μ: drift of AR processing underlying λ.",
                       tex_label="\\mu")
        m <= parameter(:σ, 1., fixed = true,
                       description="σ: volatility of AR processing underlying λ.",
                       tex_label="\\sigma")
    elseif weight_type == :equal
        m <= parameter(:λ, 0.5, (1e-5,0.999), (1e-5,0.999), SquareRoot(), Uniform(0.,1.),
                       fixed = true, description="λ: weight on model 1's predictive density",
                       tex_label="\\lambda")
    elseif weight_type == :static
        m <= parameter(:λ, 0.5, (1e-5,0.999), (1e-5,0.999), SquareRoot(), Uniform(0.,1.),
                       fixed = false, description="λ: weight on model 1's predictive density",
                       tex_label="\\lambda")
    elseif weight_type == :bma
        m <= parameter(:λ, 0.5, (1e-5,0.999), (1e-5,0.999), SquareRoot(), Uniform(0.,1.),
                       fixed = false, description="λ: weight on model 1's predictive density",
                       tex_label="\\lambda")
    end
end

function model_settings!(m::PoolModel; weight_type::Symbol = :dynamic_weight)
    default_settings!(m)

    # Weight type: dynamic, equal, or static
    if !(weight_type in [:dynamic, :equal, :static, :bma])
        error("Weight type of a PoolModel object must be :dynamic, :equal, :static, or :bma.")
    else
        m <= Setting(:weight_type, weight_type, "How to weight predictive densities")
    end

    # Data
    m <= Setting(:population_mnemonic, Nullable())
    m <= Setting(:data_id, 1922016, "Dataset identifier")
    m <= Setting(:date_presample_start, quartertodate("1992-Q1"))
    m <= Setting(:date_mainsample_start, quartertodate("1992-Q1"))
    m <= Setting(:date_forecast_start, quartertodate("2011-Q3"))

    # SMC estimation
    m <= Setting(:sampling_method, :SMC)
    m <= Setting(:n_particles, 2000)

    # Forecast
    m <= Setting(:use_population_forecast, false,
                 "Whether to use population forecasts as data")

    # Tempered particle filter
    m <= Setting(:fixed_sched, [1.],
                 "schedule for tempering in tpf; leave empty if want adaptive tempering")
    tuning = Dict(:r_star => 2., :c_init => 0.3, :target_accept_rate => 0.4,
              :resampling_method => :systematic, :n_mh_steps => 1,
              :n_particles => 1000, :n_presample_periods => 0,
              :allout => true)
    m <= Setting(:tuning, tuning, "tuning parameters for TPF")

    # MH estimation
    m <= Setting(:mh_cc, 0.15)
    m <= Setting(:mh_cc0, 0.15)
    m <= Setting(:calculate_hessian, false)
    m <= Setting(:reoptimize, false)
    m <= Setting(:n_mh_simulations, 1000)
    m <= Setting(:n_mh_blocks, 1)
    m <= Setting(:mh_thin, 1)
    m <= Setting(:n_mh_burn, 0)
end

#########################################################
# Overloading various functions for PoolModel type
#########################################################
function Base.show(io::IO, m::PoolModel)
    model_str = ""
    n_obs = n_observables(m)
    for i in 1:n_obs
        if i < n_obs
            model_str *= string(m.observable_mappings[get_key(m, :obs, i)].key) * ", "
        else
            model_str *= string(m.observable_mappings[get_key(m, :obs, i)].key) * "\n"
        end
    end
    @printf io "Dynamic Prediction Pool Method\n"
    @printf io "models: %s\n"                 model_str
    @printf io "data vintage:           %s\n" data_vintage(m)
    @printf io "description:\n %s\n"          description(m)
end

"""
```
steadystate!(m::PoolModel)
```

Computes the steady state for a `PoolModel`. Since
the current implementation does not have a useful
reason to compute the steady state, this function does nothing.

### Arguments:
- `m`: the model object
"""
function steadystate!(m::PoolModel)
end


text:
#---------------------------------------------------------
# # A Phosphorus Cycling Model
#---------------------------------------------------------

#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`p_cycle_DIP_DOP_POP.ipynb`](@__NBVIEWER_ROOT_URL__/examples/generated/p_cycle_DIP_DOP_POP.ipynb)

#---------------------------------------------------------
# ## Tracer equations
#---------------------------------------------------------

# We consider a simple model for the cycling of phosphorus with 3 state variables consisting of phosphate (PO₄) AKA dissolved inorganic phosphorus (DIP), dissolved organic phosphorus (DOP), and particulate organic phosphorus (POP).
#
# The dissolved phases are transported by advection and diffusion whereas the particulate phase sinks rapidly down the water column without any appreciable transport by the circulation.
#
# The governing equations are:
#
# $$\frac{\partial}{\partial t} DIP + \nabla \cdot \left[\boldsymbol{u} + \mathbf{K}\cdot\nabla \right] DIP = -\gamma(DIP) + \kappa_\mathsf{D} \, DOP,$$
#
# $$\frac{\partial}{\partial t} DOP + \nabla \cdot \left[\boldsymbol{u} + \mathbf{K}\cdot \nabla \right] DOP = \sigma \, \gamma(DIP) + \kappa_\mathsf{P} \, POP - \kappa_\mathsf{D} \, DOP,$$
#
# and
#
# $$\frac{\partial}{\partial t} POP + \frac{\partial}{\partial z} \left[w_\mathsf{P} \, POP\right] = (1-\sigma) \, \gamma(DIP) - \kappa_\mathsf{P} \, POP,$$
#
# where $\boldsymbol{u}$ is the fluid velocity and $\mathbf{K}$ is the eddy-diffusion tensor.
# Thus, $\nabla \cdot \left[ \boldsymbol{u} - \mathbf{K} \cdot \nabla \right]$ is a differential operator that represents the transport by the ocean circulation.
# The function $\gamma(DIP)$ represents the biological uptake of DIP by phytoplankton.

# Oxygen participates to this cycle too and satisfies its own tracer equation
#
# $$\frac{\partial}{\partial t} DO_2 + \nabla \cdot \left[\boldsymbol{u} + \mathbf{K}\cdot\nabla \right] DO_2 = -r_{\mathsf{O}_2:\mathsf{P}} \, \kappa_\mathsf{D} \, DOP + \Lambda(DO_2 - [O_2]_{\mathsf{sat}})$$
#
# where $\Lambda$ is the air-sea gas exchange operator.

# These tracer equations depend on a number of scalars, that we list below
#
# | Symbol                        | Definition                                                                  |
# |:------------------------------|:----------------------------------------------------------------------------|
# | $w_\mathsf{P}$                | depth dependent particle sinking speed                                      |
# | $\sigma$                      | fraction of the organic matter production allocated to the dissolved phase  |
# | $\kappa_\mathsf{D}$           | respiration rate for dissolved organic matter (DOP → DIP)                   |
# | $\kappa_\mathsf{P}$           | dissolution rate for particulate organic matter (POP → DOP)                 |
# | $r_{\mathsf{O}_2:\mathsf{P}}$ | number of moles of O₂ needed to respire 1 mole of DOP                       |

using AIBECS

# Load the circulation and grid

wet3D, grd, T_Circulation = OCIM0.load() ;

# Define useful constants and arrays

iwet = indices_of_wet_boxes(wet3D)
nb = length(iwet)
z = ustrip.(vector_of_depths(wet3D, grd))
ztop = vector_of_top_depths(wet3D, grd) ;

# And matrices

DIV = buildDIV(wet3D, iwet, grd)
Iabove = buildIabove(wet3D, iwet) ;

# ### Transport matrices

T_DIP(p) = T_Circulation
T_DOP(p) = T_Circulation
T_DO2(p) = T_Circulation
S₀ = buildPFD(ones(nb), DIV, Iabove)
S′ = buildPFD(ztop, DIV, Iabove)
function T_POP(p)
    w₀, w′ = p.w₀, p.w′
    return w₀ * S₀ + w′ * S′
end
T_all = (T_DIP, T_DOP, T_POP, T_DO2) ;

# Because AIBECS will solve for the steady state solution directly without time-stepping the goverining equations to equilibrium, we don't have any opportunity to specify any intial conditions.
# Initial conditions are how the total amount of conserved elements get specified in most global biogeochemical modelels.
# Thus to specify the total inventory of P in AIBECS we add a very weak resporing term to the DIP equation.
# The time-scale for this restoring term is chosen to be very long compared to the timescale with which the ocean circulation homogenizes a tracer.
# Because of this long timescale we call it the geological restoring term, but geochemists who work on geological processes don't like that name!
# In any event the long timescale allows us to prescribe the total inventory of P in a way that yields the same solution we would have gotten had we time-stepped the model to steady-state with the total inventory prescribed by the initial condition.

# ### Sources minus sinks

# ##### Geological Restoring

function geores(x, p)
    τg, xgeo = p.τg, p.xgeo
    return (xgeo .- x) / τg
end

# ##### Uptake of phosphate (DIP)

relu(x) = (x .≥ 0) .* x
function uptake(DIP, p)
    τu, ku, z₀ = p.τu, p.ku, p.z₀
    DIP⁺ = relu(DIP)
    return 1/τu * DIP⁺.^2 ./ (DIP⁺ .+ ku) .* (z .≤ z₀)
end

# ##### Remineralization DOP into DIP

function respiration(DOP, p)
    κDOP = p.κDOP
    return κDOP * DOP
end

# ##### Dissolution of POP into DOP

function dissolution(POP, p)
    κPOP = p.κPOP
    return κPOP * POP
end

# ##### Air-sea gas exchange

h = ustrip(grd.δdepth[1])             # thickness of the top layer
z = grd.depth_3D[iwet]        # depth of the gridbox centers
using WorldOceanAtlasTools
WOA = WorldOceanAtlasTools
μDO2sat , σ²DO2sat = WOA.fit_to_grid(grd,2018,"O2sat","annual","1°","an") ;
DO2sat = vec(μDO2sat)[iwet]
function airsea(DO2, p)
    κDO2 = p.κDO2
    return κDO2 * (z .< 20u"m") .* (DO2sat .- DO2) / h
end

# Add them up into sms functions (Sources Minus Sinks)

function sms_DIP(DIP, DOP, POP, DO2, p)
    return -uptake(DIP, p) + respiration(DOP, p) + geores(DIP, p)
end
function sms_DOP(DIP, DOP, POP, DO2, p)
    σ = p.σ
    return σ * uptake(DIP, p) - respiration(DOP, p) + dissolution(POP, p)
end
function sms_POP(DIP, DOP, POP, DO2, p)
    σ = p.σ
    return (1 - σ) * uptake(DIP, p) - dissolution(POP, p)
end
function sms_DO2(DIP, DOP, POP, DO2, p)
    rO2P = p.rO2P
    return airsea(DO2,p) + rO2P * (uptake(DIP,p) - respiration(DOP,p))
end
sms_all = (sms_DIP, sms_DOP, sms_POP, sms_DO2,) # bundles all the source-sink functions in a tuple

# ### Parameters
#
# Build the parameters type and p₀

t = empty_parameter_table()    # initialize table of parameters
add_parameter!(t, :xgeo, 2.17u"mmol/m^3")
add_parameter!(t, :τg, 1.0u"Myr")
add_parameter!(t, :ku, 10.0u"μmol/m^3")
add_parameter!(t, :z₀, 80.0u"m")
add_parameter!(t, :w₀, 1.0u"m/d")
add_parameter!(t, :w′, 1/4.4625u"d")
add_parameter!(t, :κDOP, 1/0.25u"yr")
add_parameter!(t, :κPOP, 1/5.25u"d")
add_parameter!(t, :σ, 0.3u"1")
add_parameter!(t, :τu, 30.0u"d")
add_parameter!(t, :κDO2, 50u"m" / 30u"d")
add_parameter!(t, :rO2P, 175.0u"mol/mol")
initialize_Parameters_type(t, "DIP_DOP_POP_O₂_Parameters")   # Generate the parameter type

# ### Generate state function and Jacobian

nt = length(T_all)    # number of tracers
n = nt * nb           # total dimension of the state vector
p = DIP_DOP_POP_O₂_Parameters() # parameters
x = p.xgeo * ones(n) # initial iterate
F, ∇ₓF = state_function_and_Jacobian(T_all, sms_all, nb)

# and solve

prob = SteadyStateProblem(F, ∇ₓF, x, p)
nothing # s = solve(prob, CTKAlg()) # Not working yet

# unpack state

DIP, DOP, POP, DO2 = state_to_tracers(x, nb, nt) # remove when line below works
nothing # DIP, DOP, POP, DO2 = state_to_tracers(s.u, nb, nt)

# We will plot the concentration of DIP at a given depth horizon

depth = grd.depth
iz = findfirst(depth .> 200u"m")
iz, depth[iz]

#-

DIP_3D = rearrange_into_3Darray(DIP, wet3D)
DIP_2D = DIP_3D[:,:,iz] * ustrip(1.0u"mol/m^3" |> u"mmol/m^3")
lat, lon = ustrip.(grd.lat), ustrip.(grd.lon)

# and plot

ENV["MPLBACKEND"]="qt5agg"
using PyPlot, PyCall
clf()
ccrs = pyimport("cartopy.crs")
ax = subplot(projection = ccrs.EqualEarth(central_longitude=-155.0))
ax.coastlines()
# making it cyclic for Cartopy
lon_cyc = [lon; 360+lon[1]]
DIP_2D_cyc = hcat(DIP_2D, DIP_2D[:,1])
# And plot
p = contourf(lon_cyc, lat, DIP_2D_cyc, levels=0:0.2:3.6, transform=ccrs.PlateCarree(), zorder=-1)
colorbar(p, orientation="horizontal");
gcf()


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

"""
    @stochastic_model(def)

Define a stochastic model capable of instantiating stochastic programs, using the syntax
```julia
@stochastic_model model_name begin
    ...
    @stage x begin
      ...
    end
    ...
end
```
where the inner blocks are [`@stage`](@ref) blocks. At least two stages must be specified in consecutive order. A stochastic model object can later be used to [`instantiate`](@ref) stochastic programs using a given set of scenarios or by using samplers. The model is referenced using `model_name` in the [`@stage`](@ref) blocks. If `model_name` is left out, the macro returns an anonymous model object, and the reserved keyword `model` must be used in the [`@stage`](@ref) blocks. Otherwise, the resulting stochastic model object is stored in a variable named `model_name`.

## Examples

The following defines a stochastic model consisitng of the first stage model given by:
```math
  minimize 100x₁ + 150x₂
    s.t  x₁ + x₂ ≤ 120
         x₁ ≥ 40
         x₂ ≥ 20
```
and the second-stage model given by:
```math
  minimize q₁(ξ)y₁ + q₂(ξ)y₂
    s.t  6y₁ + 10y₂ ≤ 60x₁
         8y₁ + 5y₂ ≤ 60x₂
         0 ≤ y₁ ≤ d₁(ξ)
         0 ≤ y₂ ≤ d₂(ξ)
```
where ``q₁(ξ), q₂(ξ), d₁(ξ), d₂(ξ)`` depend on the scenario ``ξ``.

```julia
@stochastic_model sm begin
    @stage 1 begin
        @decision(sm, x₁ >= 40)
        @decision(sm, x₂ >= 20)
        @objective(sm, Min, 100*x₁ + 150*x₂)
        @constraint(sm, x₁ + x₂ <= 120)
    end
    @stage 2 begin
        @uncertain q₁ q₂ d₁ d₂
        @variable(sm, 0 <= y₁ <= d₁)
        @variable(sm, 0 <= y₂ <= d₂)
        @objective(sm, Min, q₁*y₁ + q₂*y₂)
        @constraint(sm, 6*y₁ + 10*y₂ <= 60*x₁)
        @constraint(sm, 8*y₁ + 5*y₂ <= 80*x₂)
    end
end
```
or alternatively using anonymous syntax:
```julia
sm = @stochastic_model begin
    @stage 1 begin
        @decision(model, x₁ >= 40)
        @decision(model, x₂ >= 20)
        @objective(model, Min, 100*x₁ + 150*x₂)
        @constraint(model, x₁ + x₂ <= 120)
    end
    @stage 2 begin
        @uncertain q₁ q₂ d₁ d₂
        @variable(model, 0 <= y₁ <= d₁)
        @variable(model, 0 <= y₂ <= d₂)
        @objective(model, Min, q₁*y₁ + q₂*y₂)
        @constraint(model, 6*y₁ + 10*y₂ <= 60*x₁)
        @constraint(model, 8*y₁ + 5*y₂ <= 80*x₂)
    end
end
```
where the reserved keyword `model` is used throughout.

See also: [`@stage`](@ref), [`@parameters`](@ref), [`@decision`](@ref), [`@uncertain`](@ref)
"""
macro stochastic_model(args...)
    model_name, def = if length(args) == 1
        model_name = :model
        def = args[1]
        model_name, def
    else
        model_name = args[1]
        def = args[2]
        model_name, def
    end
    _error(x, str...) = begin
        macroname = Symbol(String(x.args[1])[2:end])
        JuMP._macro_error(macroname, prettify.(x.args[3:end]), x.args[2], str...)
    end
    stage = 0
    scenariodef = Expr(:block)
    paramdefs = Vector{Expr}()
    decisiondefs = Vector{Vector{Symbol}}()
    # Count stages
    N = 0
    prewalk(def) do x
        if @capture(x, @stage n_ arg_)
            N += 1
        end
        return x
    end
    modeldef = prewalk(def) do x
        x = if @capture(x, @stage n_ arg_)
            if @capture(arg, sp_ = stagedef_)
                x
            else
                stage += 1
                stage == n || _error(x, "Define the stages in coherent order. Expected stage $(stage) here.")
                push!(paramdefs, :(StageParameters()))
                push!(decisiondefs, Vector{Symbol}())
                if n > 1
                    pushfirst!(arg.args, :(@known $model_name $(decisiondefs[n-1]...)))
                end
                return :(@stage $n $model_name = $arg)
            end
        else
            x
        end
        paramdef = if @capture(x, @parameters arg_)
            names = Vector{Symbol}()
            default = Vector{Expr}()
            for paramdef in block(prettify(arg)).args
                if @capture(paramdef, key_Symbol = val_)
                    push!(names, key)
                    push!(default, paramdef)
                elseif @capture(paramdef, key_Symbol)
                    push!(names, key)
                else
                    _error(x, "Incorrect @parameters syntax. Specify parameter names, possibly with a default value.")
                end
            end
            code = :(StageParameters($names; $(default...)))
            Expr(:block, x.args[2], code)
        elseif @capture(x, @parameters arg__)
            names = Vector{Symbol}()
            default = Vector{Expr}()
            for paramdef in arg
                if @capture(paramdef, key_Symbol = val_)
                    push!(names, key)
                    push!(default, paramdef)
                elseif @capture(paramdef, key_Symbol)
                    push!(names, key)
                else
                    _error(x, "Incorrect @parameters syntax. Specify parameter names, possibly with a default value.")
                end
            end
            code = :(StageParameters($names; $(default...)))
            Expr(:block, x.args[2], code)
        else
            nothing
        end
        if paramdef != nothing
            paramdefs[stage] = paramdef
        end
        if @capture(x, @decision arg__)
            stage < N || _error(x, "@decision declarations cannot be used in the final stage. Consider @recourse instead.")
            if @capture(x, @decision model_ var_Symbol more__) || @capture(x, @decision model_ var_Symbol[range__] more__) ||
               @capture(x, @decision model_ var_Symbol <= ub_ more__) || @capture(x, @decision model_ var_Symbol >= lb_ more__) ||
               @capture(x, @decision model_ var_Symbol[range__] <= ub_ more__) || @capture(x, @decision model_ var_Symbol[range__] >= ln_ more__) ||
               @capture(x, @decision model_ var_Symbol in set_ more__) || @capture(x, @decision model_ var_Symbol[range__] in set_ more__) ||
               @capture(x, @decision model_ lb_ <= var_Symbol <= ub_ more__) || @capture(x, @decision model_ lb_ <= var_Symbol[range__] <= ub_ more__)
                push!(decisiondefs[stage], var)
            end
        end
        if @capture(x, @recourse arg__)
            stage == N || _error(x, "@recourse declarations can only be used in the final stage.")
        end
        if @capture(x, @uncertain arg__)
            stage == 1 && _error(x, "@uncertain declarations cannot be used in the first stage.")
        end
        scenariodef = if @capture(x, @uncertain var_Symbol::t_Symbol = scendef_)
            esc(@q begin
                @scenario $t = $scendef
            end)
        else
            scenariodef
        end
        return x
    end
    stage >= 2 || JuMP._macro_error(:stochastic_model, prettify.(def.args[2:end]), __source__, "Define at least two stages.")
    code = if length(args) == 1
        # anonymous model, return resulting expression
        code = @q begin
            $scenariodef
            StochasticModel($(esc.(paramdefs)...)) do $(esc(model_name))
                $(esc(modeldef))
            end
        end
    else
        # model name given, assign resulting expression to model object
        code = @q begin
            $scenariodef
            $(esc(model_name)) = StochasticModel($(esc.(paramdefs)...)) do $(esc(model_name))
                $(esc(modeldef))
            end
        end
    end
    return code
end


text:
mutable struct NNPolicy <: Policy
    pomdp::OCPOMDP
    nn::PyObject
    nstack::Int64
    input_shape::Tuple{Int64, Int64, Int64}
    action_map::Vector{OCAction}

    NNPolicy(pomdp::OCPOMDP, nn::PyObject) = new(pomdp, nn, 4, (84,84,1), ordered_actions(pomdp))
end

mutable struct NNBelief
    n::Int64 # number of observation to stack
    obs::Array{Float64}
end

mutable struct NNUpdater <: Updater
    problem::OCPOMDP
    n::Int64
    input_shape::Tuple{Int64, Int64, Int64}

    NNUpdater(problem::OCPOMDP) = new(problem, 4, (84,84,1))
end

function POMDPs.update(bu::NNUpdater, bold::NNBelief, a::Float64, o::OCObs)
    @assert bu.n == bold.n
    o_mat = convert_o(o, bu.problem)
    w, h = size(o_mat)
    new_obs = zeros(w, h, bu.n)
    # dequeue and append new obs
    for i=1:bu.n-1
        new_obs[:,:,i] = bold.obs[:,:,i+1]
    end
    new_obs[:,:,bu.n] = o_mat
    return NNBelief(bu.n, new_obs)
end

const OFF_KEY = -1

function POMDPs.update(bu::NNUpdater, bold::Dict{Int64,NNBelief}, a::Float64, o::Dict{Int64,OCObs})
    bnew = Dict{Int64, NNBelief}()
    for oid in keys(o)
        if haskey(bold, oid) && oid != OFF_KEY
            bnew[oid] = update(up, bold[oid], a, o[oid])
        elseif oid == OFF_KEY
            bnew[OFF_KEY] = update(up, bold[OFF_KEY], a, o[oid])
        else # cars appeared
            bnew[oid] = update(up, bold[OFF_KEY], a, o[oid])
        end
    end
    return bnew
end

function initialize_off_belief(policy::NNPolicy, ego::VehicleState)
    ped = get_off_the_grid(policy.pomdp)
    o0 = OCState(false, ego, ped)
    b0 = Dict{Int64, NNBelief}()
    o_mat = convert_o(o0, policy.pomdp)
    w, h = size(o_mat)
    o_stacked = reshape(repmat(o_mat, 1, policy.nstack), (w, h, policy.nstack))
    b0[OFF_KEY] = NNBelief(policy.nstack, o_stacked)
    return b0
end

# function normalize_o(s::OCState, pomdp::OCPOMDP)
#     o_vec = zeros(4)
#     o_vec[1], o_vec[2] = s.ego.posG.x/pomdp.x_goal, s.ego.v/pomdp.env.params.speed_limit
#     o_vec[3] = s.ped.posG.y/pomdp.y_goal
#     o_vec[4] = s.ped.v
#     if isinf(o_vec[4])
#         o_vec[4] = 0.
#     end
#     o_vec[4] /= pomdp.env.params.ped_max_speed
#     return o_vec
# end
#
# function normalize_o(x::Float64, v_ego::Float64, y::Float64, v_ped::Float64, pomdp::OCPOMDP)
#     o_vec = zeros(4)
#     o_vec[1], o_vec[2] = x/pomdp.x_goal, v_ego/pomdp.env.params.speed_limit
#     o_vec[3] = y/pomdp.y_goal
#     o_vec[4] = v_ped
#     if isinf(o_vec[4])
#         o_vec[4] = 0.
#     end
#     o_vec[4] /= pomdp.env.params.ped_max_speed
#     return o_vec
# end

function POMDPs.action(policy::NNPolicy, b::NNBelief)
    ai = policy.nn[:action](b.obs)
    return policy.action_map[ai]
end

function POMDPs.value(policy::NNPolicy, b::NNBelief)
    return policy.nn[:value](b.obs)
end

function action(policy::Policy, b::Dict{Int64, NNBelief})
    return POMDPs.action(policy, b).acc
end

const NNDriverModel = CrosswalkDriver{Float64, Dict{Int64,NNBelief},Dict{Int64, OCObs},NNPolicy,NNUpdater,POMDPSensor}

function reset_model!(model::NNDriverModel,
                      ego::Vehicle)
        model.tick = 0
        model.a = 0.
        model.b = initialize_off_belief(model.policy, ego.state)
        model.o = Dict{Int64, OCObs}()
        model.ego = ego
        reset_policy!(model.policy)
end

function reset_policy!(policy::NNPolicy)
    policy.nn[:reset]()
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
using DaemonMode
using Test
using Sockets

function init_server(port)
    task = @async serve(port, async=false)
    sleep(1)
    return task
end

function end_server(task, port)
    sendExitCode(port)
    wait(task)
end

function test_evalfile(file; port)
    task = init_server(port)
    buffer = IOBuffer()
    runfile(file, output=buffer, port=port)
    output = String(take!(buffer))
    end_server(task, port)
    return output
end

function test_evalcode_file(file; port)
    task = init_server(port)
    buffer = IOBuffer()
    code = runfile(file, output=buffer, port=port)
    output = String(take!(buffer))
    end_server(task, port)
    return output, code
end

function test_evalfiles(files; port)
    task = init_server(port)
    buffer = IOBuffer()
    outputs = String[]

    for file in files
        runfile(file, output=buffer, port=port)
        push!(outputs, String(take!(buffer)))
    end

    end_server(task, port)
    return outputs
end

@testset "Start Server" begin
    port = 3001
    @test_throws Base.IOError connect(port)
    task = @async serve(port, async=false)
    sleep(1)
    sendExitCode(port)
    wait(task)
end

@testset "runFile" begin
    files = ["hello.jl", "hello2.jl", "print.jl"]
    expected = ["Hello, World!\n", "Hello, World\n\nBye, World!\n", "Hello, World"]
    outputs = test_evalfiles(files, port=3002)

    for (expected, output) in zip(expected, outputs)
        @test output == expected
    end
end

@testset "runExpr" begin
    port = 3003
    task = @async serve(port, async=false)
    sleep(1)

    buffer = IOBuffer()
    expr = "x = 3 ; for i = 1:x ; println(i) ; end\n"
    runexpr(expr, output=buffer, port=port)
    output = String(take!(buffer))
    @test output == "1\n2\n3\n"

    buffer = IOBuffer()
    expr = "begin
        x = 2
        for i = 1:2
            println(i)
        end
    end"
    runexpr(expr, output=buffer, port=port)
    output = String(take!(buffer))
    @test output == "1\n2\n"

    expr = "begin
        x = 2
        for i = 1:2
            println(i)
        end
    end

    "

    runexpr(expr, output=buffer, port=port)
    output = String(take!(buffer))
    @test output == "1\n2\n"
    @test istaskdone(task) == false

    sendExitCode(port)
    wait(task)
end

@testset "runConflict" begin
    files = ["conflict1.jl", "conflict2.jl"]
    outputs = test_evalfiles(files, port=3004)
    expected = ["f(1) = 2\n", "f + 2 = 3\n"]

    for (expect, output) in zip(expected, outputs)
        @test expect == output
    end
end

@testset "runFileError" begin
    port = 3004
    task = init_server(port)
    files = ["bad.jl", "hello.jl"]
    expected = [1, 0]

    buffer = IOBuffer()

    for (file, return_code) in zip(files, expected)
        code = runfile(file, output=buffer, port=port)
        @test code == return_code
    end

    end_server(task, port)
end

@testset "testInclude" begin
    output = test_evalfile("include_test.jl", port=3006)
    @test output == "6\n"
end

@testset "testLogs" begin
    output = test_evalfile("test_log1.jl", port=3007)
    lines = split(output, "\n")
    # Remove colors
    lines = replace.(lines, r"\e\[.*?m"=>"")
    @test occursin("Warning: warning 1", lines[1])
    @test occursin("another line", lines[2])
    @test occursin("last one", lines[3])
    @test occursin("test_log1.jl: 4", lines[4])
    @test occursin(r"Error: error 1", lines[5])
    @test occursin("test_log1.jl: 5", lines[6])
    @test occursin("Info: info 1", lines[7])
    @test occursin("test_log1.jl: 6", lines[8])
end

@testset "testArgs" begin
    output = test_evalfile("args.jl", port=3008)
    lines = split(output, "\n")
    @test lines[1] == "String[]"
end

@testset "exit" begin
    output, code = test_evalcode_file("test_exit_0.jl", port=3009)
    @test output == "Before\n"
    @test code == 0
    output, code = test_evalcode_file("test_exit_1.jl", port=3009)
    @test output == "Before\n"
    @test code == 1
end