text:
using sem, Arrow, ModelingToolkit, Statistics, LinearAlgebra,
    SparseArrays, BenchmarkTools, Optim, LineSearches, DataFrames

## Observed Data
dat = DataFrame(Arrow.Table("test/comparisons/reg_1.arrow"))
par = DataFrame(Arrow.Table("test/comparisons/reg_1_par.arrow"))
start_lav = DataFrame(Arrow.Table("test/comparisons/reg_1_start.arrow"))

semobserved = SemObsCommon(data = Matrix{Float64}(dat))
rel_tol = 3.1956e-13

diff_fin = SemFiniteDiff(LBFGS(
    m = 50,
    alphaguess = InitialHagerZhang(),
    linesearch = HagerZhang()),
    Optim.Options(;
    f_tol = 1e-10,
    x_tol = 1.5e-8))

## Model definition
@ModelingToolkit.variables x[1:243]

#F
Ind = collect(1:25)
Jnd = collect(1:25)
V = fill(1,25)
F = sparse(Ind, Jnd, V, 25, 46)

#A
Ind = [collect(2:25); 5; collect(6:25); fill(27, 19); 1; 4]
Jnd = [fill(26, 24); 27; collect(27:46); collect(28:46); 26; 27]
V =[x[1:64]; 1.0; 1.0]
A = sparse(Ind, Jnd, V, 46, 46)

#S
Ind = [collect(1:6); 26; 27; collect(28:46)]
Jnd = [collect(1:6); 26; 27; collect(28:46)]
V = [x[65:72]; fill(1.0, 19)]
S = sparse(Ind, Jnd, V, 46, 46)
xind = 73
for i = 28:46
    for j = (i+1):46
       S[i,j] = x[xind]
       S[j,i] = x[xind]
       xind += 1
    end
end

par_order = [
    collect(2:6);
    collect(11:29);
    8;9;
    collect(49:67);
    collect(106:124);
    collect(125:130);
    131; 132;
    collect(133:303)]

start_val = Vector{Float64}(start_lav.est[par_order])

#= start_val = vcat(
    fill(1, 2),
    fill(0.5, 3),
    fill(0.1, 19),
    fill(1,2),
    fill(1,19),
    fill(0,19),
    fill(0.5, 6),
    fill(0.05, 2),
    fill(0, 171)
    ) =#

loss = Loss(
    [SemML(semobserved, [0.0], similar(start_val))])

imply = ImplySymbolic(A, S, F, x, start_val)

model_fin = Sem(semobserved, imply, loss, diff_fin)

solution_fin = sem_fit(model_fin)
@btime model_fin(start_val)
grad = similar(start_val)

using FiniteDiff
@btime FiniteDiff.finite_difference_gradient!(grad, model_fin, start_val)

FiniteDiff.GradientCache(
    c1         :: Union{Nothing,AbstractArray{<:Number}},
    c2         :: Union{Nothing,AbstractArray{<:Number}},
    fx         :: Union{Nothing,<:Number,AbstractArray{<:Number}} = nothing,
    fdtype     :: Type{T1} = Val{:central},
    returntype :: Type{T2} = eltype(df),
    inplace    :: Type{Val{T3}} = Val{true})

cache = FiniteDiff.GradientCache(grad, start_val)
@btime FiniteDiff.finite_difference_gradient!(grad, model_fin, start_val, cache)

FiniteDiff.finite_difference_gradient!(
    df,
    f,
    x,
    fdtype::Type{T1}=Val{:central},
    returntype::Type{T2}=eltype(df),
    inplace::Type{Val{T3}}=Val{true};
    [epsilon_factor])

# Cached
FiniteDiff.finite_difference_gradient!(
    df::AbstractArray{<:Number},
    f,
    x::AbstractArray{<:Number},
    cache::GradientCache;
    [epsilon_factor])

using FiniteDifferences

@btime FiniteDifferences.grad(central_fdm(2, 1), model_fin, start_val)

function (model::Sem{A, B, C, D} where {A, B, C, D <: SemFiniteDiff})(par::Vector, grad::Vector)
    if length(grad) > 0
        grad .= FiniteDifferences.grad(central_fdm(2, 1), model_fin, par)[1]
    end
    return model(par)
end

@code_warntype model_fin.loss(start_val, model_fin)

using ProfileView

function profile_test(n)
    for i = 1:n
        sem_fit(model_fin)
    end
end

ProfileView.@profview profile_test(1)

@btime model_fin(start_val)

grad = similar(start_val)

@btime model_fin(start_val, grad)

using NLopt
diff_nlopt = SemFiniteDiff(:LD_LBFGS, nothing)
model_nlopt = Sem(semobserved, imply, loss, diff_nlopt)
@btime solution_nlopt = sem.sem_fit_nlopt(model_nlopt, 3.1956e-13)

function profile_test_nlopt(n)
    for i = 1:n
        sem.sem_fit_nlopt(model_nlopt, 3.1956e-13)
    end
end

ProfileView.@profview profile_test_nlopt(10)

all(
    abs.(solution_fin.minimizer .- par.est[par_order]
        ) .< 0.05*abs.(par.est[par_order]))

##### regularized
solutions = []

ml = SemML(semobserved, [0.0], similar(start_val))
con = SemConstant(-(logdet(semobserved.obs_cov) + semobserved.n_man))

for i = 1:50
    ridge = sem.SemRidge(0.01*i, 46:64)
    loss = Loss([ml, con, ridge])
    model_fin = Sem(semobserved, imply, loss, diff_fin)
    push!(solutions, sem_fit(model_fin))
end

getindex.(getfield.(solutions, :minimizer), 46)

maximum(abs.(getfield.(solutions, :minimizer)[50] .- solution_fin.minimizer))

fieldnames(typeof(solutions[1]))

pars = getfield.(solutions, :minimizer)
pars = getindex.(pars, [46:64])

pars = permutedims(hcat(pars...))

#using DataFrames

#pars = DataFrame(pars)

#using Plots

plot(pars)

savefig("pars.pdf")


text:
function [ zu,zw,xu,xs ] = griddef(nz,nx,dz,dx,zu,zw,xu,xs)
#UNTITLED2 Summary of this function goes here
#   Detailed explanation goes here

# Assign height values to u-scalar heights

for k = 1 : nz
  zu[k] = (k - 1.5)*dz;
end

# Assign height values to w heights

for k = 2 : nz
  zw[k] = dz*(k-2);
end

# Assign x-dist values to u on the x grid-pt

for i = 2 : nx

    xu[i] = dx*(i-2);

end

# Assign x-dist values to scalars on the x grid-pt

# xs(1) = dx/2;

for i = 2 : nx

#    xs(i) = dx/2 + dx*(i-1);
     xs[i] = (i-1.5)*dx;

end

end

ext: jl

max_stars_repo_name: jeffreycunn/toy-atmospheric-model

max_issues_repo_name: jeffreycunn/toy-atmospheric-model

max_forks_repo_name: jeffreycunn/toy-atmospheric-model


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "DecFP"
version = v"2.2.0" # 2.0 Update 2

# Collection of sources required to build DecFP
sources = [
    "https://www.netlib.org/misc/intel/IntelRDFPMathLib20U2.tar.gz" =>
    "93c0c78e0989df88f8540bf38d6743734804cef1e40706fd8fe5c6a03f79e173",
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd IntelRDFPMathLib20U2
sed -i -e 's/^#if \(!defined _MSC_VER || defined __INTEL_COMPILER\)/#if !defined __ENABLE_BINARY80__ \&\& (\1)/' LIBRARY/src/bid_functions.h TESTS/test_bid_functions.h
cd LIBRARY
if [ "$nbits" == 64 ]; then
     _HOST_ARCH=x86_64
else
     _HOST_ARCH=x86
fi

if [[ $target == *"-w64-"* ]]; then
    make CC_NAME=cc _HOST_OS=Windows_NT AR_CMD="ar rv" _HOST_ARCH=$_HOST_ARCH CALL_BY_REF=0 GLOBAL_RND=1 GLOBAL_FLAGS=1 UNCHANGED_BINARY_FLAGS=1
    $CC $LDFLAGS -shared -o libbid.$dlext *.obj
    mkdir -p $prefix/bin
    cp libbid.$dlext $prefix/bin/
else
    CFLAGS_OPT="-fPIC -fsigned-char -D__ENABLE_BINARY80__=0"
    if [[ $target == *"-musl"* ]]; then
        CFLAGS_OPT+=" -D__QNX__"
    elif [[ $target == *"freebsd"* ]]; then
        CFLAGS_OPT+=" -D__linux"
    fi
    make CC_NAME=cc CFLAGS_OPT="$CFLAGS_OPT" CFLAGS="$CFLAGS_OPT" _HOST_ARCH=$_HOST_ARCH CALL_BY_REF=0 GLOBAL_RND=1 GLOBAL_FLAGS=1 UNCHANGED_BINARY_FLAGS=1
    $CC $LDFLAGS -shared -o libbid.$dlext *.o
    mkdir -p $prefix/lib
    cp libbid.$dlext $prefix/lib/
fi

"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Linux(:i686, libc=:glibc),
    Linux(:x86_64, libc=:glibc),
    Linux(:aarch64, libc=:glibc),
    Linux(:armv7l, libc=:glibc, call_abi=:eabihf),
    Linux(:powerpc64le, libc=:glibc),
    Linux(:i686, libc=:musl),
    Linux(:x86_64, libc=:musl),
    Linux(:aarch64, libc=:musl),
    Linux(:armv7l, libc=:musl, call_abi=:eabihf),
    MacOS(:x86_64),
    FreeBSD(:x86_64),
    Windows(:i686),
    Windows(:x86_64)
]

# The products that we will ensure are always built
products(prefix) = [
    LibraryProduct(prefix, "libbid", :libbid)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jmkuhn/DecFPBuilder

max_issues_repo_name: jmkuhn/DecFPBuilder

max_forks_repo_name: jmkuhn/DecFPBuilder


text:
"""
    change_sense(objective_sense)

Change the objective sense of optimization.
Possible arguments are `MOI.MAX_SENSE` and `MOI.MIN_SENSE`.

If you want to change the objective and sense at the same time, use
[`change_objective`](@ref) instead to do both at once.
"""
function change_sense(objective_sense)
    (model, opt_model) -> COBREXA.JuMP.set_objective_sense(opt_model, objective_sense)
end

"""
    change_optimizer(optimizer)

Change the JuMP optimizer used to run the optimization.

This may be used to try different approaches for reaching the optimum, and in
problems that may require different optimizers for different parts, such as the
[`parsimonious_flux_balance_analysis`](@ref).
"""
function change_optimizer(optimizer)
    (model, opt_model) -> COBREXA.JuMP.set_optimizer(opt_model, optimizer)
end

"""
    change_optimizer_attribute(attribute_key, value)

Change a JuMP optimizer attribute. The attributes are optimizer-specific, refer
to the JuMP documentation and the documentation of the specific optimizer for
usable keys and values.
"""
function change_optimizer_attribute(attribute_key, value)
    (model, opt_model) ->
        COBREXA.JuMP.set_optimizer_attribute(opt_model, attribute_key, value)
end

"""
    constrain_objective_value(tolerance)

Limit the objective value to `tolerance`-times the current objective value, as
with [`objective_bounds`](@ref).
"""
function constrain_objective_value(tolerance)
    return (model, opt_model) -> begin
        λmin, λmax = objective_bounds(tolerance)(objective_value(opt_model))
        old_objective = objective_function(opt_model)
        @constraint(opt_model, λmin <= sum(old_objective) <= λmax)
    end
end

"""
    change_constraint(id::String, lb, ub)

Change the lower and upper bounds (`lb` and `ub` respectively) of reaction `id`.
"""
function change_constraint(id::String, lb, ub)
    (model, opt_model) -> begin
        ind = first(indexin([id], reactions(model)))
        set_bound(ind, opt_model, lb = lb, ub = ub)
    end
end

"""
    change_objective(new_objective::Union{String,Vector{String}}; weights=[], sense=MOI.MAX_SENSE)

Modification that changes the objective function used in a constraint based
analysis function.  `new_objective` can be a single reaction identifier, or an
array of reactions identifiers.

Optionally, the objective can be weighted by a vector of `weights`, and a
optimization `sense` can be set.
"""
function change_objective(
    new_objective::Union{String,Vector{String}};
    weights = [],
    sense = MOI.MAX_SENSE,
)
    (model, opt_model) -> begin

        # Construct objective_indices array
        if typeof(new_objective) == String
            objective_indices = [first(indexin([new_objective], reactions(model)))]
        else
            objective_indices =
                [first(indexin([rxnid], reactions(model))) for rxnid in new_objective]
        end

        # Initialize weights
        opt_weights = spzeros(n_reactions(model))

        isempty(weights) && (weights = ones(length(objective_indices))) # equal weights

        for (j, i) in enumerate(objective_indices)
            opt_weights[i] = weights[j]
        end

        v = opt_model[:x]
        @objective(opt_model, sense, sum(opt_weights[i] * v[i] for i in objective_indices))
    end
end


text:
immutable PokeRequest
    rid::Int64                 # request id
    data::Vector{PokeRecord}   # data
end

immutable PokeResponse
    rid::Int64                 # id of corresponding request
    status::Int16              # request status code
    message::ASCIIString       # human-readable error message (if any)
end


text:
@testset "Types" begin
    # Test the classifier
    svm = LSSVC()
    @test svm.kernel == :rbf
    @test svm.γ == 1.0
    @test svm.σ == 1.0
    @test svm.degree == 0.0

    # Test the regressor
    svr = LSSVR()
    @test svm.kernel == :rbf
    @test svm.γ == 1.0
    @test svm.σ == 1.0
    @test svm.degree == 0.0
end

max_stars_repo_name: gmagannaDevelop/Elysivm

max_issues_repo_name: gmagannaDevelop/Elysivm

max_forks_repo_name: gmagannaDevelop/Elysivm


text:
function correct_grid_size(ft::DataType)
    g = RegularCartesianGrid(ft, (4, 6, 8), (2π, 4π, 9π))

    # Checking ≈ as the grid could be storing Float32 values.
    (g.Nx ≈ 4  && g.Ny ≈ 6  && g.Nz ≈ 8 &&
     g.Lx ≈ 2π && g.Ly ≈ 4π && g.Lz ≈ 9π)
end

function correct_cell_volume(ft::DataType)
    Nx, Ny, Nz = 19, 13, 7
    Δx, Δy, Δz = 0.1, 0.2, 0.3
    Lx, Ly, Lz = Nx*Δx, Ny*Δy, Nz*Δz
    g = RegularCartesianGrid(ft, (Nx, Ny, Nz), (Lx, Ly, Lz))

    # Checking ≈ as the grid could be storing Float32 values.
    g.V ≈ Δx*Δy*Δz
end

function faces_start_at_zero(ft::DataType)
    g = RegularCartesianGrid(ft, (10, 10, 10), (2π, 2π, 2π))
    g.xF[1] == 0 && g.yF[1] == 0 && g.zF[1] == 0
end

function end_faces_match_grid_length(ft::DataType)
    g = RegularCartesianGrid(ft, (12, 13, 14), (π, π^2, π^3))
    (g.xF[end] - g.xF[1] ≈ π   &&
     g.yF[end] - g.yF[1] ≈ π^2 &&
     g.zF[1] - g.zF[end] ≈ π^3)
end

@testset "Grids" begin
    println("Testing grids...")

    @testset "Grid initialization" begin
        println("  Testing grid initialization...")
        for FT in float_types
            @test correct_grid_size(FT)
            @test correct_cell_volume(FT)
            @test faces_start_at_zero(FT)
            @test end_faces_match_grid_length(FT)
        end
    end

    @testset "Grid dimensions" begin
        println("  Testing grid dimensions...")
        L = (100, 100, 100)
        for FT in float_types
            @test isbitstype(typeof(RegularCartesianGrid(FT, (16, 16, 16), (1, 1, 1))))

            @test_throws ArgumentError RegularCartesianGrid(FT, (32,), L)
            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 64), L)
            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32, 16), L)
            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32), (100,))
            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32), (100, 100))
            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32), (100, 100, 1, 1))
            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32), (100, 100, -100))

            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32.0), (1, 1, 1))
            @test_throws ArgumentError RegularCartesianGrid(FT, (20.1, 32, 32), (1, 1, 1))
            @test_throws ArgumentError RegularCartesianGrid(FT, (32, nothing, 32), (1, 1, 1))
            @test_throws ArgumentError RegularCartesianGrid(FT, (32, "32", 32), (1, 1, 1))
            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32), (1, nothing, 1))
            @test_throws ArgumentError RegularCartesianGrid(FT, (32, 32, 32), (1, "1", 1))
        end
    end
end


text:
#module LazyJSONBenchmark

using LazyJSON
using DataFrames
using JSON
using JSON2

using Mmap

#j = String(read("ec2-2016-11-15.normal.json"))
f = open("ec2-2016-11-15.normal.json", "r")
j = String(Mmap.mmap(f));

j = string("[", join(fill(j, 10), ", "), "]")

j2 = String(read("geo.json"))

struct HTTP
    method::String
    requestUri::String
end

struct IOType
    shape::String
end

Base.convert(::Type{IOType}, ::Nothing) = IOType("")

struct Operation
    name::String
    http::HTTP
    input::IOType
    output::IOType
    documentation::String
end

JSON2.@format Operation begin
    output => (default=IOType(""),)
end

struct TestLazyA end
struct TestLazyB end
struct TestLazyC end
struct TestJSON end
struct TestJSON2 end

name(::Type{TestLazyA}) = "Lazy"
name(::Type{TestLazyB}) = "Lazy (B)"
name(::Type{TestLazyC}) = "Lazy (C)"
name(::Type{TestJSON}) = "JSON"
name(::Type{TestJSON2}) = "JSON2"

println("""
    For these tests, the content of `ec2-2016-11-15.normal.json` has been
    duplicated 10 times into a top level JSON array "[ , , , ...]" this
    results in an overall input data size of ~10MB.

    ## test1
    Reads `ec2-2016-11-15.normal.json` and extracts a single value:
    `operations.AcceptReservedInstancesExchangeQuote.input.shape`.
    This value is close to the start of the input data.

    Variants:
     - Lazy: LazyJSON.jl `AbstractDict` interface.
     - Lazy (B): LazyJSON.jl `getproperty` interface.
     - Lazy (C): LazyJSON.jl `lazy=false` (parse whole input to Dicts etc)
     - JSON: JSON.jl `parse` interface.
     - JSON2: JSON2.jl `read -> NamedTuple` interface.

    """)

test1(::Type{TestLazyA}, j) =
    LazyJSON.parse(j)[1]["operations"
                    ]["AcceptReservedInstancesExchangeQuote"
                    ]["input"
                    ]["shape"]

test1(::Type{TestLazyB}, j) =
    LazyJSON.parse(j; getproperty=true)[1].
                      operations.
                      AcceptReservedInstancesExchangeQuote.
                      input.
                      shape

test1(::Type{TestLazyC}, j) =
    LazyJSON.parse(j; lazy=false)[1]["operations"
                ]["AcceptReservedInstancesExchangeQuote"
                ]["input"
                ]["shape"]

test1(::Type{TestJSON}, j) =
    JSON.parse(j)[1]["operations"
                ]["AcceptReservedInstancesExchangeQuote"
                ]["input"
                ]["shape"]

test1(::Type{TestJSON2}, j) =
    JSON2.read(j)[1].operations.
                  AcceptReservedInstancesExchangeQuote.
                  input.
                  shape

const n = 10

results = DataFrame(Test=String[],
                    Variant=String[],
                    µs=Int[],
                    bytes=Int[],
                    poolalloc=Int[],
                    bigalloc=Int[])

function run_test1(T, j; warmup=false)

    count = 0
    val, t, bytes, gctime, memallocs = @timed for i in 1:n
        v = test1(T, j)
        if v == "AcceptReservedInstancesExchangeQuoteRequest"
            count += 1
        end
    end
    @assert count == n

    if !warmup
        t = round(Int, t/n * 1000000)
        push!(results, ["test1", name(T), t, bytes/n,
                        round(Int, memallocs.poolalloc/n),
                        round(Int, memallocs.bigalloc/n)])
    end
end

tests = (TestLazyA,
         TestLazyB,
         TestLazyC,
         TestJSON,
         TestJSON2)

for warmup in [true, false]
    for t in tests
        run_test1(t, j; warmup=warmup)
    end
end

println("```")
@show results
println("```\n")
println("""
    Note: LazyJSON.jl is similar to JSON.jl in speed and memory use in non-lazy mode.

    """)
results = results[1:0,:]

println("""
    ## test2
    Read `ec2-2016-11-15.normal.json` and extracts an array value:
    `shapes.scope.enum`
    This value is close to the end of the input data.

    Variants:
     - Lazy: LazyJSON.jl `AbstractDict` interface.
     - Lazy (B): LazyJSON.jl `getproperty` interface.
     - Lazy (C): LazyJSON.jl `lazy=false` (parse whole input to Dicts etc)
     - JSON: JSON.jl `parse` interface.
     - JSON2: JSON2.jl `read -> NamedTuple` interface.
    """)

test2(::Type{TestLazyA}, j) =
    LazyJSON.parse(j)[10]["shapes"]["scope"]["enum"]

test2(::Type{TestLazyB}, j) =
        LazyJSON.parse(j; getproperty=true)[10].shapes.scope.enum

test2(::Type{TestLazyC}, j) =
        LazyJSON.parse(j; lazy=false)[10]["shapes"]["scope"]["enum"]

test2(::Type{TestJSON}, j) =
    JSON.parse(j)[10]["shapes"]["scope"]["enum"]

test2(::Type{TestJSON2}, j) =
    JSON2.read(j)[10].shapes.scope.enum

function run_test2(T, j; warmup=false)

    count = 0
    val, t, bytes, gctime, memallocs = @timed for i in 1:n
        v = test2(T, j)
        if v[1] == "Availability Zone" && v[2] == "Region"
            count += 1
        end
    end
    @assert count == n

    if !warmup
        t = round(Int, t/n * 1000000)
        push!(results, ["test2", name(T), t, bytes/n,
                        round(Int, memallocs.poolalloc/n),
                        round(Int, memallocs.bigalloc/n)])
    end
end

for warmup in [true, false]
    for t in tests
        run_test2(t, j; warmup=warmup)
    end
end

println("```")
@show results
println("```\n")
println("""
    Note: It takes LazyJSON.jl a bit longer to access values near the end of
    the input.

    """)

results = results[1:0,:]

println("""
    ## test3
    Modifes `ec2-2016-11-15.normal.json` by replacing a value near the
    start of the file and two values near the end.

    Variants:
     - Lazy: LazyJSON.jl `getproperty` interface finds values and
             `LazyJSON.splice` modifies the JSON data in-place.
     - JSON: JSON.jl `parse` to `Dict`, modify, then write new JSON text.
     - JSON2: _Parses to immutable `NamedTuples`. Modificaiton not supported._
    """)

function test3(::Type{TestLazyA}, j)
    r = LazyJSON.value(j; getproperty=true)
    r = LazyJSON.splice(r, r[10].shapes.ZoneNameStringList.member.shape, "foo")
    r = LazyJSON.splice(r, r[10].shapes.scope.enum[1], "bar")
    r = LazyJSON.splice(r, r[10].shapes.scope.enum[2], "fum")
    r = LazyJSON.value(string(r); getproperty=true)
    return r[10].shapes.scope.enum[1]
end

function test3(::Type{TestJSON}, j)
    r = JSON.parse(j)
    r[10]["shapes"]["ZoneNameStringList"]["member"]["shape"] = "foo"
    r[10]["shapes"]["scope"]["enum"][1] = "bar"
    r[10]["shapes"]["scope"]["enum"][2] = "fum"
    s = JSON.json(r)
    r = JSON.parse(s)
    return r[10]["shapes"]["scope"]["enum"][1]
end

tests = (TestLazyA,
         TestJSON)

function run_test3(T, j; warmup=false)

    count = 0
    val, t, bytes, gctime, memallocs = @timed for i in 1:n
        v = test3(T, j)
        if v == "bar"
            count += 1
        end
    end
    @assert count == n

    if !warmup
        t = round(Int, t/n * 1000000)
        push!(results, ["test3", name(T), t, bytes/n,
                        round(Int, memallocs.poolalloc/n),
                        round(Int, memallocs.bigalloc/n)])
    end
end

for warmup in [true, false]
    for t in tests
        run_test3(t, j; warmup=warmup)
    end
end

println("```")
@show results
println("```\n")
results = results[1:0,:]

println("""
    ## test4
    Reads a 1.2MB GeoJSON file an extracts a country name near the middle
    of the file.

    Variants:
     - Lazy: LazyJSON.jl `AbstractDict` interface.
       `LazyJSON.parse(j)["features"][15]["properties"]["formal_en"]`
     - Lazy (B): LazyJSON.jl `getproperty` interface.
       `LazyJSON.parse(j; getproperty=true).features[15].properties.formal_en`
     - Lazy (C): LazyJSON.jl `lazy=false` (parse whole input to Dicts etc)
       `LazyJSON.parse(j; lazy=false)["features"][15]["properties"]["formal_en"]`
     - JSON: JSON.jl `parse` interface.
       `JSON.parse(j)["features"][15]["properties"]["formal_en"]`
     - JSON2: JSON2.jl `read -> NamedTuple` interface.
       `JSON2.read(j).features[15].properties.formal_en`
    """)

test4(::Type{TestLazyA}, j) =
    LazyJSON.parse(j)["features"][15]["properties"]["formal_en"]

test4(::Type{TestLazyB}, j) =
    LazyJSON.parse(j; getproperty=true).features[15].properties.formal_en

test4(::Type{TestLazyC}, j) =
    LazyJSON.parse(j; lazy=false)["features"][15]["properties"]["formal_en"]

test4(::Type{TestJSON}, j) =
    JSON.parse(j)["features"][15]["properties"]["formal_en"]

test4(::Type{TestJSON2}, j) =
    JSON2.read(j).features[15].properties.formal_en

function run_test4(T, j; warmup=false)

    count = 0
    val, t, bytes, gctime, memallocs = @timed for i in 1:n
        v = test4(T, j)
        if v == "Republic of Palau"
            count += 1
        end
    end
    @assert count == n

    if !warmup
        t = round(Int, t/n * 1000000)
        push!(results, ["test4", name(T), t, bytes/n,
                        round(Int, memallocs.poolalloc/n),
                        round(Int, memallocs.bigalloc/n)])
    end
end

tests = (TestLazyA,
         TestLazyB,
         TestLazyC,
         TestJSON,
         TestJSON2)

for warmup in [true, false]
    for t in tests
        run_test4(t, j2; warmup=warmup)
    end
end

println("```")
@show results
println("```\n")

println("""
    Note: LazyJSON.jl in non-lazy mode is a bit faster than JSON.jl for this
    input.

    """)

results = results[1:0,:]

println("""
    ## test5
    Reads a 1.2MB GeoJSON file and checks that the outline polygon for
    a single country is within an expected lat/lon range.
    ```
    r = r["features"][15]["geometry"]["coordinates"][6][1]
    @assert r[1][1] == 134.41651451900023
    for (x, y) in r
       @assert 134.2 < x < 134.5
       @assert 7.21 < y < 7.32
    end
    ```
    """)

function test5(::Type{TestLazyA}, j)
    r = LazyJSON.parse(j)
    r = r["features"][15]["geometry"]["coordinates"][6][1]
    @assert r[1][1] == 134.41651451900023
    for (x, y) in r
       @assert 134.2 < x < 134.5
       @assert 7.21 < y < 7.32
    end
end

function test5(::Type{TestJSON}, j)
    r = JSON.parse(j)
    r = r["features"][15]["geometry"]["coordinates"][6][1]
    @assert r[1][1] == 134.41651451900023
    for (x, y) in r
       @assert 134.2 < x < 134.5
       @assert 7.21 < y < 7.32
    end
end

function test5(::Type{TestJSON2}, j)
    r = JSON2.read(j)
    r = r.features[15].geometry.coordinates[6][1]
    @assert r[1][1] == 134.41651451900023
    for (x, y) in r
       @assert 134.2 < x < 134.5
       @assert 7.21 < y < 7.32
    end
end

function run_test5(T, j; warmup=false)

    val, t, bytes, gctime, memallocs = @timed for i in 1:n
        test5(T, j)
    end

    if !warmup
        t = round(Int, t/n * 1000000)
        push!(results, ["test5", name(T), t, bytes/n,
                        round(Int, memallocs.poolalloc/n),
                        round(Int, memallocs.bigalloc/n)])
    end
end

tests = (TestLazyA,
         TestJSON,
         TestJSON2)

for warmup in [true, false]
    for t in tests
        run_test5(t, j2; warmup=warmup)
    end
end

println("```")
@show results
println("```\n")
results = results[1:0,:]

println("""
    ## test6
    Defines `struct Operation`, `struct IOType` and `struct HTTP` with
    fields that match the API operations data in `ec2-2016-11-15.normal.json`.
    It then does JSON2-style direct-to-struct parsing to read the JSON data
    into a Julia object `Dict{String,Operation}`
    (LazyJSON provides `@generated` `Base.convert` methods for this).

    Variants:
     - Lazy: LazyJSON.jl `AbstractDict` interface.
        `convert(Dict{String,Operation}, LazyJSON.parse(j))`
     - JSON2: JSON2.jl `read -> NamedTuple` interface.
        `JSON2.read(j, Dict{String,Operation})`
    """)

function test6(::Type{TestLazyA}, j)
    v = convert(Dict{String,Operation}, LazyJSON.parse(j))
    return v["AcceptReservedInstancesExchangeQuote"].input.shape
end

function test6(::Type{TestJSON2}, j)
    v = JSON2.read(j, Dict{String,Operation})
    return v["AcceptReservedInstancesExchangeQuote"].input.shape
end

function run_test6(T, j; warmup=false)

    count = 0
    val, t, bytes, gctime, memallocs = @timed for i in 1:n
        v = test6(T, j)
        if v == "AcceptReservedInstancesExchangeQuoteRequest"
            count += 1
        end
    end

    if !warmup
        t = round(Int, t/n * 1000000)
        push!(results, ["test6", name(T), t, bytes/n,
                        round(Int, memallocs.poolalloc/n),
                        round(Int, memallocs.bigalloc/n)])
    end
end

tests = (TestLazyA,
         TestJSON2)

jops = string(LazyJSON.parse(j)[1]["operations"])

for warmup in [true, false]
    for t in tests
        run_test6(t, jops; warmup=warmup)
    end
end
display(results)
println("```")
@show results
println("```\n")
#end


text:
module TestDeprecated

using Test, DataFrames

const ≅ = isequal

@testset "by and aggregate" begin
    @test_throws ArgumentError by()
    @test_throws ArgumentError aggregate()
end

@testset "deprecated broadcasting assignment" begin
    df = DataFrame(a=1:4, b=1, c=2)
    df.a .= 'a':'d'
    @test df == DataFrame(a=97:100, b=1, c=2)
    dfv = view(df, 2:3, 2:3)
    dfv.b .= 0
    @test df.b == [1, 0, 0, 1]
end

@testset "indicator in joins" begin
    name = DataFrame(ID = [1, 2, 3], Name = ["John Doe", "Jane Doe", "Joe Blogs"])
    job = DataFrame(ID = [1, 2, 4], Job = ["Lawyer", "Doctor", "Farmer"])

    @test outerjoin(name, job, on = :ID, indicator=:source) ≅
          outerjoin(name, job, on = :ID, source=:source)
    @test leftjoin(name, job, on = :ID, indicator=:source) ≅
          leftjoin(name, job, on = :ID, source=:source)
    @test rightjoin(name, job, on = :ID, indicator=:source) ≅
          rightjoin(name, job, on = :ID, source=:source)

    @test_throws ArgumentError outerjoin(name, job, on = :ID,
                                         indicator=:source, source=:source)
    @test_throws ArgumentError leftjoin(name, job, on = :ID,
                                       indicator=:source, source=:source)
    @test_throws ArgumentError rightjoin(name, job, on = :ID,
                                         indicator=:source, source=:source)
end

end # module


text:
@testset "Testing kalman.jl" begin

    # set up
    A = [.95 0; 0. .95]
    Q = Matrix(I, 2, 2) .* 0.5
    G = Matrix(I, 2, 2) .* 0.5
    R = Matrix(I, 2, 2) .* 0.2
    kf = Kalman(A, G, Q, R)

    rough_kwargs = Dict(:atol => 1e-2, :rtol => 1e-4)

    sig_inf, kal_gain = stationary_values(kf)

    # Compute the Kalman gain and sigma infinity according to the
    # recursive equations and compare
    mat_inv = inv(G * sig_inf * G' + R)

    kal_recursion = A * sig_inf * (G') * mat_inv
    sig_recursion = A * sig_inf * A' - kal_recursion * G * sig_inf * A' + Q
    # test stationary
    @test isapprox(sig_inf, sig_recursion; rough_kwargs...)
    @test isapprox(kal_gain, kal_recursion; rough_kwargs...)

    # test_update_using_stationary
    set_state!(kf, [0.0 0.0]', sig_inf)
    update!(kf, [0.0 0.0]')
    @test isapprox(kf.cur_sigma, sig_inf; rough_kwargs...)
    @test isapprox(kf.cur_x_hat, [0.0 0.0]'; rough_kwargs...)

    # test update nonstationary
    curr_x, curr_sigma = fill(one(Float64), 2, 1), Matrix(I, 2, 2) .* .75
    y_observed = fill(0.75, 2, 1)
    set_state!(kf, curr_x, curr_sigma)
    update!(kf, y_observed)
    mat_inv = inv(G * curr_sigma * G' + R)
    curr_k = A * curr_sigma * (G') * mat_inv
    new_sigma = A * curr_sigma * A' - curr_k * G * curr_sigma * A' + Q
    new_xhat = A * curr_x + curr_k * (y_observed - G * curr_x)
    @test isapprox(kf.cur_sigma, new_sigma; rough_kwargs...)
    @test isapprox(kf.cur_x_hat, new_xhat; rough_kwargs...)

    # test smooth
    A = [.5 .3;
         .1 .6]
    Q = [1 .1;
        .1 .8]
    G = A
    R = Q/10
    kn = Kalman(A, G, Q, R)
    cov_init = [1.76433188153014    0.657255445961981
                0.657255445961981    1.40080308176678]
    set_state!(kn, zeros(2,1), cov_init)
    y =
    [-0.852059587559099 0.366646911007266   1.16654837064550    -0.183684880920978  -0.429545082282421
     -0.365266818777927  -1.16089364984681   0.0441137903162171  -0.641364287040669  0.289104799234833]
    x_smoothed, logL, P_smoothed = smooth(kn, y)
    x =
    [-0.842404354244332 1.09920664072649    1.72823043243333    0.180454617610392   -0.711307666172619
     -0.616562264406391  -1.46832223756820   -0.275840062339632  -0.802470118677828  0.282423885770331]
    @test isapprox(x, x_smoothed)
    @test isapprox(-12.7076583170714, logL)

    set_state!(kn, zeros(2,1), cov_init)
    @test isapprox(-12.7076583170714, compute_loglikelihood(kn, y))
end  # @testset


text:
using SubHunt
using Base.Test

using POMDPs
using POMDPToolbox
using DiscreteValueIteration
using QMDP
using ParticleFilters
using Plots

pomdp = SubHuntPOMDP()

# show(STDOUT, MIME("text/plain"), SubVis(pomdp))

rng = MersenneTwister(6000)
# policy = RandomPolicy(pomdp, rng=rng)

solver = QMDPSolver()
if !isdefined(:policy) || policy.pomdp != pomdp
    policy = solve(solver, pomdp, verbose=true)
end

s = initial_state(pomdp, rng)
@show value(policy, s)
@show value(policy, SubState(s.own, s.target, s.goal, true))

pfp = PingFirst(policy)

filter = SIRParticleFilter(pomdp, 100000, rng=rng)

for (s,a,r,sp,o,bp) in stepthrough(pomdp, pfp, filter, "sarspobp", max_steps=20, rng=rng)
    v = SubVis(pomdp, s=s, a=a, r=r)
    show(STDOUT, MIME("text/plain"), v)
    vp = SubVis(pomdp, s=sp, b=bp)
    plot(vp)
    gui()
    sleep(0.5)
end


text:
module TimeSeriesA6

#Use this to print
#=
using Weave
codeName = "TimeSeriesA6"
weave(Pkg.dir("$(pwd())\\$(codeName).jl"),
  informat="script",
  out_path = "$(pwd())\\output\\$(codeName)_Appendix.html",
  doctype = "md2html")
=#

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

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

importall CTIO, CTReg, CTStat, CTVAR, CTNumerical

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 NIPA_DATE_FORMAT = "m/d/yyyy"
const SP500_NAME = "GSPC"
const DBV_NAME = "DBV"
const FRED_NAME = "CPI"
const PPI_NAME = "PPI"
const A6_DATA_NAME = "MktRet_DP_TermSpread"

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, 3, Ints=true)

#preprocess rate data
function preProcessA6()::Void
  rateDF::DataFrame = CSV.read("$DATA_PATH\\$A6_DATA_NAME.csv")

  #showcols(rateDF)
  rename!(rateDF, [:MktExRet, :Mkt_DP,	:y10minFedFunds], [:re, :dp, :y10e])

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

  return nothing
end

#Calculates the summary statistics
function summaryA6(rateDF, cols::Vector{Symbol} = [:re, :dp, :y10e],
      rowNames::Vector{String} = ["mean", "\$\\sigma\$", "\$\\rho_{t-1}\$", "Half Life", "N"]
    )
  #label the table
  colNames::Vector{Vector{String}} =
    [((s::Symbol)->replace(String(s),"_","")).(cols)]
  numCols = length(colNames[end])

  numRows::Int = length(rowNames)
  content::Vector{Vector{String}} = [Vector{String}(numCols) for i::Int ∈ 1:numRows]

  #get the statistics and record the results
  for c::Int ∈ 1:numCols
    ρ₁::Float64 = autoρ(rateDF[:, cols[c]], 1)

    content[1][c] = n2s(mean(rateDF[:, cols[c]]))
    content[2][c] = n2s(std(rateDF[:, cols[c]]))
    content[3][c] = n2s(ρ₁)
    content[4][c] = n2s(log(0.5)/log(ρ₁))
    content[5][c] = n2s(length(rateDF[:, cols[c]]))
  end

  #write the table
  summaryTable::String = texTable( "Summary of Data Series",
    """See Tex File""", #caption
    colNames, #colNames
    Vector{String}(),#contentRowNames
    Vector{Matrix{String}}(), #content
    rowNames, #descRowNames
    content, #descContent
    Vector{String}(), #notes
  )

  return summaryTable
end

#VAR table
function VARTableA6(rateDF::DataFrame, LHS::Vector{Symbol}, RHS::Vector{Symbol}, est::VAREstimate;
    title="Estimate")
  #label the table

  K::Int = length(LHS)
  colNames::Vector{Vector{String}} =
    [["\$\\phi_0\$"; (String).(RHS); "N"; "\$R^2\$";"\$R^2 (Adj)\$"; "F"]]
  numCols = length(colNames[end])

  rowNames::Vector{String} = ((i::Int)->isodd(i)?(String(LHS[(i+1)÷2])):"").(1:(2K))
  numRows::Int = length(rowNames)
  content::Vector{Vector{String}} = [Vector{String}(numCols) for i::Int ∈ 1:numRows]

  #run the tests and record the results
  for k::Int ∈ 1:K
    r::Int = 2(k-1) + 1 #pointer to the beginning of each double row

    #need these for the adjusted R2 and F statistics
    MSE::Float64 = est.SSE[k] / (est.N[k]-K-1.)
    MSR::Float64 = est.SSR[k] / K
    MST::Float64 = est.SST[k] /(est.N[k] - 1.)

    content[r][1] = n2s(est.var.ϕ₀[k])
    content[r+1][1] = est.coeffErrors[k] ≠ nothing ? "($(n2s(est.coeffErrors[k][1])))" : ""
    for c::Int ∈ 2:(numCols-4)
      content[r][c] = n2s(est.var.ϕ₁[k, c-1])
      content[r+1][c] = est.coeffErrors[k] ≠ nothing ? "($(n2s(est.coeffErrors[k][c])))" : ""
    end

    content[r][end-3] = est.N[k] ≠ nothing ? "$(n2s(est.N[k]))" : ""
    content[r+1][end-3] = ""
    content[r][end-2] = est.SSE[k] ≠ nothing ? "$(n2s(est.SSR[k]/est.SST[k]))" : ""
    content[r+1][end-2] = ""
    content[r][end-1] = est.SSE[k] ≠ nothing ? "$(n2s(1.-MSE/MST))" : ""
    content[r+1][end-1] = ""
    content[r][end] = est.SSE[k] ≠ nothing ? "$(n2s(MSR/MSE))" : ""
    content[r+1][end] = ""

  end

  #write the table
  summaryTable::String = texTable( title,
    """See Tex File""", #caption
    colNames, #colNames
    Vector{String}(),#contentRowNames
    Vector{Matrix{String}}(), #content
    rowNames, #descRowNames
    content, #descContent
    Vector{String}(), #notes
  )

  return summaryTable
end

#make the problem-specific graphs
function visualizeA6(rateDF::DataFrame, predictions::Vector{VARPrediction},
      responses::Vector{Matrix{Float64}})::Void

  plots::Vector{PlotContainer} = Vector{PlotContainer}()
  plotNames::Vector{String} = Vector{String}()

  Gadfly.push_theme(:default)
  outDF = deepcopy(rateDF)
  dropmissing!(outDF)

  push!(plotNames, "summaryRatesAIO")
  push!(plots, plot(
    layer(outDF, x=:DATE, y=:y10e,  Geom.line, Theme(default_color=COLORS[1])),
    layer(outDF, x=:DATE, y=:dp,  Geom.line, Theme(default_color=COLORS[2])),
    layer(outDF, x=:DATE, y=:re,  Geom.line, Theme(default_color=COLORS[3])),
    Guide.manual_color_key("Legend", ["Excess Yield", "dp",
      "Excess Return"], COLORS[1:3]),
    Guide.xlabel("Date"),Guide.ylabel("Rate"),
    Coord.Cartesian(xmin=1950., xmax=2018.),
    Guide.title("Excess 10-year yield, dividend yield, and excess market return")))

  push!(plotNames, "10Y Excess Yield Predictions")
  push!(plots, plot(
    layer(x=predictions[1].dates, y=predictions[1].Y[1,:],  Geom.line,
      Theme(default_color=COLORS[1])),
    layer(x=predictions[1].dates, y=predictions[1].YHat[1,:],  Geom.line,
      Theme(default_color=COLORS[2])),
    layer(x=predictions[2].dates, y=predictions[2].YHat[1,:],  Geom.line,
      Theme(default_color=COLORS[3])),
    layer(x=predictions[3].dates, y=predictions[3].YHat[1,:],  Geom.line,
      Theme(default_color=COLORS[4])),
    Guide.manual_color_key("Legend", ["Actual", "t+1 Prediction",
      "t+4 Prediction", "t+20 Prediction"], COLORS[1:4]),
    Coord.Cartesian(xmin=1950., xmax=2018.),
    Guide.xlabel("Date"),Guide.ylabel("Excess Yield"),
    Guide.title("Excess yield predictions")))

  push!(plotNames, "Dividend yield predictions")
  push!(plots, plot(
    layer(x=predictions[1].dates, y=predictions[1].Y[2,:],  Geom.line,
      Theme(default_color=COLORS[1])),
    layer(x=predictions[1].dates, y=predictions[1].YHat[2,:],  Geom.line,
      Theme(default_color=COLORS[2])),
    layer(x=predictions[2].dates, y=predictions[2].YHat[2,:],  Geom.line,
      Theme(default_color=COLORS[3])),
    layer(x=predictions[3].dates, y=predictions[3].YHat[2,:],  Geom.line,
      Theme(default_color=COLORS[4])),
    Guide.manual_color_key("Legend", ["Actual", "t+1 Prediction",
      "t+4 Prediction", "t+20 Prediction"], COLORS[1:4]),
    Guide.xlabel("Date"),Guide.ylabel("Dividend Yield"),
    Coord.Cartesian(xmin=1950., xmax=2018.),
    Guide.title("Dividend yield predictions")))

  push!(plotNames, "Excess Market Predictions")
  push!(plots, plot(
    layer(x=predictions[1].dates, y=predictions[1].Y[3,:],  Geom.line,
      Theme(default_color=COLORS[1])),
    layer(x=predictions[1].dates, y=predictions[1].YHat[3,:],  Geom.line,
      Theme(default_color=COLORS[2])),
    layer(x=predictions[2].dates, y=predictions[2].YHat[3,:],  Geom.line,
      Theme(default_color=COLORS[3])),
    layer(x=predictions[3].dates, y=predictions[3].YHat[3,:],  Geom.line,
      Theme(default_color=COLORS[4])),
    Guide.manual_color_key("Legend", ["Actual", "t+1 Prediction",
      "t+4 Prediction", "t+20 Prediction"], COLORS[1:4]),
    Guide.xlabel("Date"),Guide.ylabel("Rate of Return"),
    Coord.Cartesian(xmin=1950., xmax=2018.),
    Guide.title("Excess market return predictions")))

    #Impulse responses
    xLags::Vector{Int} = collect(0:(size(responses[1],2)-1))

    push!(plotNames, "Yield Impulse Response")
    push!(plots, plot(
      layer(x=xLags, y=responses[1][1,:],  Geom.line,
        Theme(default_color=COLORS[1])),
      layer(x=xLags, y=responses[1][2,:],  Geom.line,
        Theme(default_color=COLORS[2])),
      layer(x=xLags, y=responses[1][3,:],  Geom.line,
        Theme(default_color=COLORS[3])),
      Guide.manual_color_key("Legend", ["Excess Yield", "dp",
        "Excess Return"], COLORS[1:3]),
      Guide.xlabel("Lag"),Guide.ylabel("Rate"), Guide.title("Yield Impulse Response")))

    push!(plotNames, "Dividend Yield Impulse Response")
    push!(plots, plot(
      layer(x=xLags, y=responses[2][1,:],  Geom.line,
        Theme(default_color=COLORS[1])),
      layer(x=xLags, y=responses[2][2,:],  Geom.line,
        Theme(default_color=COLORS[2])),
      layer(x=xLags, y=responses[2][3,:],  Geom.line,
        Theme(default_color=COLORS[3])),
      Guide.manual_color_key("Legend", ["Excess Yield", "dp",
        "Excess Return"], COLORS[1:3]),
      Guide.xlabel("Lag"),Guide.ylabel("Rate"), Guide.title("Dividend Yield Impulse Response")))

    push!(plotNames, "Market Return Impulse Response")
    push!(plots, plot(
      layer(x=xLags, y=responses[3][1,:],  Geom.line,
        Theme(default_color=COLORS[1])),
      layer(x=xLags, y=responses[3][2,:],  Geom.line,
        Theme(default_color=COLORS[2])),
      layer(x=xLags, y=responses[3][3,:],  Geom.line,
        Theme(default_color=COLORS[3])),
      Guide.manual_color_key("Legend", ["Excess Yield", "dp",
        "Excess Return"], COLORS[1:3]),
      Guide.xlabel("Lag"),Guide.ylabel("Rate"), Guide.title("Excess Return Impulse")))

    push!(plotNames, "DRₜ-EDRₜ and dpₜ-Edpₜ")
    push!(plots, plot(
      layer(outDF, x=:DATE, y=:rDREDR,  Geom.line,
        Theme(default_color=COLORS[1])),
      layer(outDF, x=:DATE, y=:ldpEdp,  Geom.line,
        Theme(default_color=COLORS[2])),
      Guide.manual_color_key("Legend", ["DRₜ-EDRₜ", "dpₜ"], COLORS[1:2]),
      Guide.xlabel("Date"),Guide.ylabel("Rate"), Guide.title("DRₜ-EDRₜ and dpₜ-Edpₜ")))

  #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

#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 A6StationarityTest(var::VAR; verbose::Bool = true)::Void

  #get the eigenvectors
  λFact::Base.LinAlg.Eigen = eigfact(var.ϕ₁)
  λ::Vector{T where T<:Number} = λFact[:values]

  stationaryTests::BitArray = (abs2).(λ) .< 1.

  if verbose
    for i::Int ∈ 1:length(stationaryTests)
      println("$(stationaryTests[i]?"Passed":"Failed"): λ=$(n2s(λ[i])) |λ|=$(n2s(abs(λ[i])))")
    end
  end

  println("The VAR is $(prod(stationaryTests)?"Stationary":"Not Stationary")")

  return nothing
end

function CSA6(var::VAR, Y::Matrix{Float64}; ρ::Float64=0.96)
  K::Int = size(Y,1)
  T::Int = size(Y,2)
  I::Matrix{Float64} = eye(K)

  pre::Matrix{Float64} = ρ*(I .- ρ*var.ϕ₁)\I
  μ::Vector{Float64} = (I .- var.ϕ₁)\I*var.ϕ₀

  M::Matrix{Float64} = pre * (Y .- μ)

  return M
end

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

  lags::Int = 20
  LHS::Vector{Symbol} = [:y10e, :dp, :re]
  RHS::Vector{Symbol} = [:y10eL1, :dpL1, :reL1]

  #get the dataframe
  rateDF::DataFrame = getDF("$DATA_PATH\\$A6_DATA_NAME.jls", refreshData, preProcessA6)
  T::Int = size(rateDF, 1)

  #make the summary tables
  tables::Vector{String} = Vector{String}()
  push!(tables, summaryA6(rateDF))

  lagDF!(rateDF, LHS) #make the lags

  #esitmate the VAR
  est = VAREstimate(rateDF, LHS, RHS, errors=getNeweyWestFunc(40))
  push!(tables, VARTableA6(rateDF, LHS, RHS, est, title="VAR 1 for Rates"))

  #make the predictions
  validXY::Vector{Bool} = completecases(rateDF[RHS]) .& completecases(rateDF[LHS])
  X::Matrix{Float64} = Matrix{Float64}(rateDF[validXY, RHS])'
  Y::Matrix{Float64} = Matrix{Float64}(rateDF[validXY, LHS])'
  dates::Vector{Float64} = rateDF[validXY, :DATE]

  #container for holding all predictions
  predictions::Vector{VARPrediction} = Vector{VARPrediction}()
  push!(predictions, VARPrediction(est.var, dates, X, Y, τ = 1))
  push!(predictions, VARPrediction(est.var, dates, X, Y, τ = 4))
  push!(predictions, VARPrediction(est.var, dates, X, Y, τ = 20))

  responses::Vector{Matrix{Float64}} = Vector{Matrix{Float64}}()
  push!(responses, propagateImpulse(est.var, [1.,0.,0.]))
  push!(responses, propagateImpulse(est.var, [0.,1.,0.]))
  push!(responses, propagateImpulse(est.var, [0.,0.,1.]))

  #get the DR-EDR quantity and store in the dataframe
  rateDF[:, :rDREDR] = Vector{Union{Float64,Missing}}(T)
  rateDF[:, :rDREDR] = missing
  rateDF[validXY, :rDREDR] = CSA6(est.var, Y)[3,:]
  rateDF[:,:ldpEdp] = (log).(rateDF[:,:dp])

  rateDF[:,:ldpEdp] -= mean(skipmissing(rateDF[:,:ldpEdp]))

  #make the graphs

  visualizeA6(rateDF, predictions, responses)

  A6StationarityTest(est.var)

  println("Predicted return volatility: ", ((diag(varΣ(est.var)).*4).^0.5)[3])
  println("Historical Vol of Return: $(std(Y[3,:])*2)")
  println("Conditional volatility of return: $((est.var.L*est.var.L')[3,3]^0.5*2)")

  #write out the tables
  writeTables2File(tables,
    HEADER_NAME, FOOTER_NAME, path=OUTPUT_PATH,
    outName = "TimeSeries A6.tex")

  return nothing
end

@time begin
  A6Script(refreshData = true)
end

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
module geometry

export chain, triangle, rectangle, square, hexagonal, box, cube, ring, excitation_phases, lhc1, lhc1_p

"""
    geometry.chain(a, N)

Positions of spins on a chain in x-direction.

The chain starts at the origin and continues into positive x-direction.

# Arguments
* `a`: Spin-spin distance.
* `N`: Number of spins
"""
chain(a::T, N::Int) where T<:Real = [[i*a, zero(T), zero(T)] for i=0:N-1]

"""
    geometry.triangle(a)

Positions of spins on a equilateral triangle in the xy-plane with edge length `a`.

The positions are: (0,0,0), (a,0,0), (a/2, h, 0)
"""
triangle(a::T) where T<:Real = Vector{float(T)}[[0, 0, 0], [a, 0, 0], [a/2, a*sqrt(3)/2, 0]]

"""
    geometry.rectangle(a, b; Nx=2, Ny=2)

Positions of spins on a rectangular lattice in the xy-plane.

The lattice starts at the origin and continues into positive x and y direction.

# Arguments
* `a`: Spin-spin distance in x-direction.
* `b`: Spin-spin distance in y-direction.
* `Nx=2`: Number of spins into x direction.
* `Ny=2`: Number of spins into y direction.
"""
rectangle(a::S, b::T; Nx::Int=2, Ny::Int=2) where {S<:Real,T<:Real} = vec([[i*a,zero(S),zero(S)]+[zero(T),j*b,zero(T)] for i=0:Nx-1, j=0:Ny-1])

"""
    geometry.square(a; Nx=2, Ny=2)

Positions of spins on a square lattice in the xy-plane with distance `a`.

The lattice starts at the origin and continues into positive x and y direction.
"""
square(a::Real; Nx::Int=2, Ny::Int=2) = rectangle(a, a; Nx=Nx, Ny=Ny)

"""
    geometry.hexagonal(a; Nr=1)

Positions of spins on a hexagonal lattice in the xy-plane.

The hexagonal lattice consists of `Nr` rings.
"""
function hexagonal(a::T; Nr::Int=1) where T<:Real
    positions = Vector{float(T)}[]
    ax = sqrt(3.0/4)*a
    for ix=-Nr:Nr
        if isodd(ix)
            Ny = div(2*Nr+1-abs(ix),2)
            for iy=-Ny:Ny-1
                push!(positions, [ax*ix, (0.5+iy)*a, 0])
            end
        else
            Ny = div(2*Nr-abs(ix),2)
            for iy=-Ny:Ny
                push!(positions, [ax*ix, iy*a, 0])
            end
        end
    end
    return positions
end

"""
    geometry.box(a, b, c; Nx=2, Ny=2, Nz=2)

Positions of spins on a orthorhombic lattice.

The lattice starts at the origin and continues into
positive x, y and z direction.
"""
box(a::S, b::T, c::U; Nx::Int=2, Ny::Int=2, Nz::Int=2) where {S<:Real,T<:Real,U<:Real}= vec([[i*a,zero(S),zero(S)] + [zero(T),j*b,zero(T)] + [zero(U),zero(U),k*c] for i=0:Nx-1, j=0:Ny-1, k=0:Nz-1])

"""
    geometry.cube(a; Nx=2, Ny=2, Nz=2)

Positions of spins on a cubic lattice with edge length `a`.

The lattice starts at the origin and continues into
positive x, y and z direction.
"""
cube(a::Real; Nx::Int=2, Ny::Int=2, Nz::Int=2) = box(a, a, a, Nx=Nx, Ny=Ny, Nz=Nz)

"""
    geometry.ring(a, N; distance = false)

Ring of N particles with radius a.
If distance is set to true, then `a` gives the distance between particles
and the radius is determined accordingly.
"""
function ring(a::Real, N::Int; distance = false)
    x = distance ? 0.5a/sin(pi/N) : a
    return [[x*cos(2*pi*j/N), x*sin(2*pi*j/N), 0] for j=0:(N-1)]
end

"""

    geometry.ring_rd(a, N)

    Calculate the radius of a ring when the distance between particles is a.

"""
ring_rd(a::Real, N::Int) = 0.5*a/sin(pi/N)

"""
    geometry.ring_p_tangential(N)

    Return tangential polarization vectors for a ring of N atoms in the xy-plane created by ring(a, N).

    Arguments:
    * `N`: Number of emitters.
"""
ring_p_tangential(N::Int) = [[-sin(2*pi*j/N), cos(2*pi*j/N), 0.] for j=0:(N-1)]

"""
    geometry.excitation_phases(k, pos)

    Calculate the excitation phases created in a geometry by an incident plane wave with wave vector k.

    Arguments:
    * `k`: Wave vector of the incident plane wave.
    * `pos`: List of atomic positions.
"""
function excitation_phases(k::Vector{T}, pos::Vector) where T
    @assert length(k) == 3
    k_ = normalize(k)

    PHI = T[]

    for p=pos
        @assert length(p) == 3
        phi = (dot(k_, p) * 2*pi) % (2*pi)
        push!(PHI, phi)
    end

    return PHI
end

"""

    geometry.lhc1()

    Geometry of the LHC-I light harvesting complex in the xy-plane centered around `z=0` with a mean interatomic distance of `1.0`. Positions `1` to `32` constitute the outer ring, while the inner ring comprises positions `33` to `38`.
"""
function lhc1()
    R = [[4.84223, -1.36377, 0.0],[5.10879, -0.398266, 0.0214459],[4.99561, 0.593013, 0.0],[4.87235, 1.58711, 0.0214459],[4.3883, 2.45968, 0.0],[3.89406, 3.33084, 0.0214459],[3.11302, 3.95179, 0.0],[2.32309, 4.56744, 0.0214459],[1.36377, 4.84223, 0.0],[0.398266, 5.10879, 0.0214459],[-0.593013, 4.99561, 0.0],[-1.58711, 4.87235, 0.0214459],[-2.45968, 4.3883, 0.0],[-3.33084, 3.89406, 0.0214459],[-3.95179, 3.11302, 0.0],[-4.56744, 2.32309, 0.0214459],[-4.84223, 1.36377, 0.0],[-5.10879, 0.398266, 0.0214459],[-4.99561, -0.593013, 0.0],[-4.87235, -1.58711, 0.0214459],[-4.3883, -2.45968, 0.0],[-3.89406, -3.33084, 0.0214459],[-3.11302, -3.95179, 0.0],[-2.32309, -4.56744, 0.0214459],[-1.36377, -4.84223, 0.0],[-0.398266, -5.10879, 0.0214459],[0.593013, -4.99561, 0.0],[1.58711, -4.87235, 0.0214459],[2.45968, -4.3883, 0.0],[3.33084, -3.89406, 0.0214459],[3.95179, -3.11302, 0.0],[4.56744, -2.32309, 0.0214459],[0.359274, 0.249228, -0.403357],[-0.36848, -0.221933, -0.376279],[-0.308908, 1.14151, -0.119144],[0.364906, -1.07327, 0.0079068],[0.108506, 1.35747, 1.00961],[0.0429294, -1.1805, 1.16267]]

    return R
end

"""

    geometry.lhc1_p()

    Normalized polarization vectors in an LHC-I light harvesting complex. Entries `1` to `32` pertain to the outer ring, while entries `33` to `38` give the inner ring's dipole orientations.

    Note: This function does not return a geometry, butrather te polarizations of dipoles in the LHC-I complex.

"""
function lhc1_p()
    p = [[0.634, 0.76, 0.147],[0.452, 0.886, 0.098],[0.295, 0.944, 0.147],[0.079, 0.992, 0.098],[-0.089, 0.985, 0.147],[-0.307, 0.947, 0.098],[-0.459, 0.876, 0.147],[-0.646, 0.757, 0.098],[-0.76, 0.634, 0.147],[-0.886, 0.452, 0.098],[-0.944, 0.295, 0.147],[-0.992, 0.079, 0.098],[-0.985, -0.089, 0.147],[-0.947, -0.307, 0.098],[-0.876, -0.459, 0.147],[-0.757, -0.646, 0.098],[-0.634, -0.76, 0.147],[-0.452, -0.886, 0.098],[-0.295, -0.944, 0.147],[-0.079, -0.992, 0.098],[0.089, -0.985, 0.147],[0.307, -0.947, 0.098],[0.459, -0.876, 0.147],[0.646, -0.757, 0.098],[0.76, -0.634, 0.147],[0.886, -0.452, 0.098],[0.944, -0.295, 0.147],[0.992, -0.079, 0.098],[0.985, 0.089, 0.147],[0.947, 0.307, 0.098],[0.876, 0.459, 0.147],[0.757, 0.646, 0.098],[0.371, 0.871, -0.321],[-0.379, -0.86, -0.342],[-0.3, -0.87, -0.391],[0.208, 0.855, -0.476],[0.53, 0.019, 0.848],[-0.367, -0.053, 0.929]]

    return p
end
end # module


text:
FLOAT = Float64
INT = Int
FLOATCUDA = Float32
INTCUDA = Int32

VALID_TYPES = [
    :Identity,
    :Local,
    :Global,
    :GlobalArray,
    :Medium
]

VALID_UPDATES = [
    :UpdateGlobal,
    :UpdateLocal,
    :UpdateLocalInteraction,
    :UpdateInteraction,
    :UpdateMedium,
    :UpdateMediumInteraction,
    :Equation,
]

POSITIONSYMBOLS = [:x,:y,:z]

BASICARGS=[:t,:N,:dt,:simulationBox]

PREDECLAREDPARAMETERS = Dict("Local"=>[:x,:y,:z],"Identity"=>:id,"Integration"=>[:dt,:dW],"Community"=>[:N],"Evolve"=>[:nMax])

UPDATINGOPERATORS = [:+= ,:-=,:*=,:/=,:\=,:÷=,:%=,:^=,:&=,:|=,:⊻=,:>>>=,:>>=,:<<=]

RESERVEDSYMBOLS = [:x,:y,:z,:id,:t,:N,:dt,:dW,:nMax,
                    :idMediumX,:idMediumY,:idMediumZ,
                    :ic1_,:ic2_,:nnic2_,:pos_,
                    :localV,:identityV,:globalV,:mediumV,
                    :localVCopy,:identityVCopy,:globalVCopy,:mediumVCopy,
                    :simulationBox,:radiusInteraction,:n_,
                    :ALGORITHM_,:ARGS_,:AUX1_,:AUX2_,
                    :index_,:stride_,:lockadd_,
                    :∇,:∇x,:∇y,:∇z,:Δ,:Δx,:Δy,:Δz,:δx,:δy,:δz,:xₘ,:yₘ,:zₘ,:∑ₐ,:∑ₙ];

GLOBALARRAYCOPY = "_Copy"

DIFFSYMBOL = :d

INTERACTIONSYMBOLS = [:i,:j]
INTERACTIONSYMBOLSDICT = Dict([:i=>"[ic1_,VAR]",:j=>"[nnic2_,VAR]"])

DIFFMEDIUMSYMBOL = :∂t
MEDIUMSYMBOLS = [:∂t,:∇,:∇x,:∇y,:∇z,:Δ,:Δx,:Δy,:Δz,:δr]

#Adaptation functions
GPUINFUNCTION =
"index_ = (threadIdx().x) + (blockIdx().x - 1) * blockDim().x
stride_ = blockDim().x * gridDim().x
"
CPUINFUNCTION =
"Threads.@threads"

GPUOUTFUNCTION =
"CUDA.@cuda threads=threads_ blocks=nBlocks_ "
CPUOUTFUNCTION =
""

GPUARRAY =
"CUDA."
CPUARRAY =
""

GPUARRAYEMPTY =
"CUDA.CuArray{Float32}"
CPUARRAYEMPTY =
"Array{Float64}"

GPUARRAYEMPTYINT =
"CUDA.CuArray{Int32}"
CPUARRAYEMPTYINT =
"Array{Int64}"


text:
const strip_patterns                = UInt32(0)
const strip_corrupt_utf8            = UInt32(0x1) << 0
const strip_case                    = UInt32(0x1) << 1
const stem_words                    = UInt32(0x1) << 2
const tag_part_of_speech            = UInt32(0x1) << 3

const strip_whitespace              = UInt32(0x1) << 5
const strip_punctuation             = UInt32(0x1) << 6
const strip_numbers                 = UInt32(0x1) << 7
const strip_non_letters             = UInt32(0x1) << 8

const strip_indefinite_articles     = UInt32(0x1) << 9
const strip_definite_articles       = UInt32(0x1) << 10
const strip_articles                = (strip_indefinite_articles |
                                       strip_definite_articles)

const strip_prepositions            = UInt32(0x1) << 13
const strip_pronouns                = UInt32(0x1) << 14

const strip_stopwords               = UInt32(0x1) << 16
const strip_sparse_terms            = UInt32(0x1) << 17
const strip_frequent_terms          = UInt32(0x1) << 18

const strip_html_tags               = UInt32(0x1) << 20

const alpha_sparse      = 0.05
const alpha_frequent    = 0.95

const regex_cache = Dict{AbstractString, Regex}()
function mk_regex(regex_string)
    d = haskey(regex_cache, regex_string) ?
            regex_cache[regex_string] :
            (regex_cache[regex_string] = Regex(regex_string))
    (length(regex_cache) > 50) && empty!(regex_cache)
    d
end

"""
    remove_corrupt_utf8(str)
Remove corrupt UTF8 characters in `str`.
See also: [`remove_corrupt_utf8!`](@ref)
"""
function remove_corrupt_utf8(s::AbstractString)
    return map(x->isvalid(x) ? x : ' ', s)
end

remove_corrupt_utf8!(d::FileDocument) = error("FileDocument cannot be modified")

"""
    remove_corrupt_utf8!(doc)
    remove_corrupt_utf8!(crps)
Remove corrupt UTF8 characters for `doc` or documents in `crps`.
Does not support `FileDocument` or Corpus containing `FileDocument`.
See also: [`remove_corrupt_utf8`](@ref)
"""
function remove_corrupt_utf8!(d::StringDocument)
    d.text = remove_corrupt_utf8(d.text)
    nothing
end

function remove_corrupt_utf8!(d::TokenDocument)
    for i in 1:length(d.tokens)
        d.tokens[i] = remove_corrupt_utf8(d.tokens[i])
    end
end

function remove_corrupt_utf8!(d::NGramDocument)
    new_ngrams = Dict{AbstractString, Int}()
    for token in keys(d.ngrams)
        new_token = remove_corrupt_utf8(token)
        if haskey(new_ngrams, new_token)
            new_ngrams[new_token] = new_ngrams[new_token] + 1
        else
            new_ngrams[new_token] = 1
        end
    end
    d.ngrams = new_ngrams
end

function remove_corrupt_utf8!(crps::Corpus)
    for doc in crps
        remove_corrupt_utf8!(doc)
    end
end

"""
    remove_case(str)
Convert `str` to lowercase.
See also: [`remove_case!`](@ref)
"""
remove_case(s::T) where {T <: AbstractString} = lowercase(s)

"""
    remove_case!(doc)
    remove_case!(crps)
Convert the text of `doc` or `crps` to lowercase.
Does not support `FileDocument` or `crps` containing `FileDocument`.
# Example
```julia-repl
julia> str = "The quick brown fox jumps over the lazy dog"
julia> sd = StringDocument(str)
A StringDocument{String}
 * Language: Languages.English()
 * Title: Untitled Document
 * Author: Unknown Author
 * Timestamp: Unknown Time
 * Snippet: The quick brown fox jumps over the lazy dog
julia> remove_case!(sd)
julia> sd.text
"the quick brown fox jumps over the lazy dog"
```
See also: [`remove_case`](@ref)
"""
remove_case!(d::FileDocument) = error("FileDocument cannot be modified")

function remove_case!(d::StringDocument)
    d.text = remove_case(d.text)
    nothing
end

function remove_case!(d::TokenDocument)
    for i in 1:length(d.tokens)
        d.tokens[i] = remove_case(d.tokens[i])
    end
end

function remove_case!(d::NGramDocument)
    new_ngrams = Dict{AbstractString, Int}()
    for token in keys(d.ngrams)
        new_token = remove_case(token)
        if haskey(new_ngrams, new_token)
            new_ngrams[new_token] = new_ngrams[new_token] + 1
        else
            new_ngrams[new_token] = 1
        end
    end
    d.ngrams = new_ngrams
end

function remove_case!(crps::Corpus)
    for doc in crps
        remove_case!(doc)
    end
end

const script_tags = Regex("<script\\b[^>]*>([\\s\\S]*?)</script>")
const style_tags = Regex("<style\\b[^>]*>([\\s\\S]*?)</style>")
const html_tags = Regex("<[^>]*>")

"""
    remove_html_tags(str)
Remove html tags from `str`, including the style and script tags.
See also: [`remove_html_tags!`](@ref)
"""
function remove_html_tags(s::AbstractString)
    s = remove_patterns(s, script_tags)
    s = remove_patterns(s, style_tags)
    remove_patterns(s, html_tags)
end

"""
    remove_html_tags!(doc::StringDocument)
    remove_html_tags!(crps)
Remove html tags from the `StringDocument` or documents `crps`.
Does not work for documents other than `StringDocument`.
# Example
```julia-repl
julia> html_doc = StringDocument(
             "
               <html>
                   <head><script language=\"javascript\">x = 20;</script></head>
                   <body>
                       <h1>Hello</h1><a href=\"world\">world</a>
                   </body>
               </html>
             "
            )
A StringDocument{String}
 * Language: Languages.English()
 * Title: Untitled Document
 * Author: Unknown Author
 * Timestamp: Unknown Time
 * Snippet:  <html> <head><s
julia> remove_html_tags!(html_doc)
julia> strip(text(html_doc))
"Hello world"
```
See also: [`remove_html_tags`](@ref)
"""
function remove_html_tags!(d::AbstractDocument)
    error("HTML tags can be removed only from a StringDocument")
end

function remove_html_tags!(d::StringDocument)
    d.text = remove_html_tags(d.text)
    nothing
end

function remove_html_tags!(crps::Corpus)
    for doc in crps
        remove_html_tags!(doc)
    end
end

"""
    remove_words!(doc, words::Vector{AbstractString})
    remove_words!(crps, words::Vector{AbstractString})
Remove the occurences of words from `doc` or `crps`.
# Example
```julia-repl
julia> str="The quick brown fox jumps over the lazy dog"
julia> sd=StringDocument(str);
julia> remove_words = ["fox", "over"]
julia> remove_words!(sd, remove_words)
julia> sd.text
"the quick brown   jumps   the lazy dog"
```
"""
function remove_words!(entity::(Union{AbstractDocument,Corpus}),
               words::Vector{T}) where T <: AbstractString
    skipwords = Set{AbstractString}()
    union!(skipwords, words)
    prepare!(entity, strip_patterns, skip_words = skipwords)
end

##############################################################################
#
# Part-of-Speech tagging
#
##############################################################################

function tag_pos!(entity::Union{Corpus, TokenDocument, StringDocument})
    @warn "tag_pos! is deprecated, Use Perceptrontagger instead"
    tagger = PerceptronTagger(true)
    if typeof(entity) == Corpus
        entity.documents = broadcast(x -> predict(tagger, x), entity.documents)
    elseif typeof(entity) == StringDocument
        entity = predict(tagger, tokenize(entity.text))
    end
end

"""
    sparse_terms(crps, alpha=0.05])
Find the sparse terms from Corpus, occuring in less than `alpha` percentage of the documents.
# Example
```
julia> crps = Corpus([StringDocument("This is Document 1"),
                      StringDocument("This is Document 2")])
A Corpus with 2 documents:
* 2 StringDocument's
* 0 FileDocument's
* 0 TokenDocument's
* 0 NGramDocument's
Corpus's lexicon contains 0 tokens
Corpus's index contains 0 tokens
julia> sparse_terms(crps, 0.5)
2-element Array{String,1}:
 "1"
 "2"
```
See also: [`remove_sparse_terms!`](@ref), [`frequent_terms`](@ref)
"""
function sparse_terms(crps::Corpus, alpha::Real = alpha_sparse)
    update_lexicon!(crps)
    update_inverse_index!(crps)
    res = Array{String}(undef, 0)
    ndocs = length(crps.documents)
    for term in keys(crps.lexicon)
        f = length(crps.inverse_index[term]) / ndocs
        if f <= alpha
            push!(res, String(term))
        end
    end
    return res
end

"""
    frequent_terms(crps, alpha=0.95)
Find the frequent terms from Corpus, occuring more than `alpha` percentage of the documents.
# Example
```
julia> crps = Corpus([StringDocument("This is Document 1"),
                      StringDocument("This is Document 2")])
A Corpus with 2 documents:
 * 2 StringDocument's
 * 0 FileDocument's
 * 0 TokenDocument's
 * 0 NGramDocument's
Corpus's lexicon contains 0 tokens
Corpus's index contains 0 tokens
julia> frequent_terms(crps)
3-element Array{String,1}:
 "is"
 "This"
 "Document"
```
See also: [`remove_frequent_terms!`](@ref), [`sparse_terms`](@ref)
"""
function frequent_terms(crps::Corpus, alpha::Real = alpha_frequent)
    update_lexicon!(crps)
    update_inverse_index!(crps)
    res = Array{String}(undef, 0)
    ndocs = length(crps.documents)
    for term in keys(crps.lexicon)
        f = length(crps.inverse_index[term]) / ndocs
        if f >= alpha
            push!(res, String(term))
        end
    end
    return res
end

"""
    remove_sparse_terms!(crps, alpha=0.05)
Remove sparse terms in crps, occuring less than `alpha` percent of documents.
# Example
```julia-repl
julia> crps = Corpus([StringDocument("This is Document 1"),
                      StringDocument("This is Document 2")])
A Corpus with 2 documents:
 * 2 StringDocument's
 * 0 FileDocument's
 * 0 TokenDocument's
 * 0 NGramDocument's
Corpus's lexicon contains 0 tokens
Corpus's index contains 0 tokens
julia> remove_sparse_terms!(crps, 0.5)
julia> crps[1].text
"This is Document "
julia> crps[2].text
"This is Document "
```
See also: [`remove_frequent_terms!`](@ref), [`sparse_terms`](@ref)
"""
remove_sparse_terms!(crps::Corpus, alpha::Real = alpha_sparse) = remove_words!(crps, sparse_terms(crps, alpha))

"""
    remove_frequent_terms!(crps, alpha=0.95)
Remove terms in `crps`, occuring more than `alpha` percent of documents.
# Example
```julia-repl
julia> crps = Corpus([StringDocument("This is Document 1"),
                      StringDocument("This is Document 2")])
A Corpus with 2 documents:
* 2 StringDocument's
* 0 FileDocument's
* 0 TokenDocument's
* 0 NGramDocument's
Corpus's lexicon contains 0 tokens
Corpus's index contains 0 tokens
julia> remove_frequent_terms!(crps)
julia> text(crps[1])
"     1"
julia> text(crps[2])
"     2"
```
See also: [`remove_sparse_terms!`](@ref), [`frequent_terms`](@ref)
"""
remove_frequent_terms!(crps::Corpus, alpha::Real = alpha_frequent) = remove_words!(crps, frequent_terms(crps, alpha))

"""
    prepare!(doc, flags)
    prepare!(crps, flags)
Preprocess document or corpus based on the input flags.
# List of Flags
* strip_patterns
* strip_corrupt_utf8
* strip_case
* stem_words
* tag_part_of_speech
* strip_whitespace
* strip_punctuation
* strip_numbers
* strip_non_letters
* strip_indefinite_articles
* strip_definite_articles
* strip_articles
* strip_prepositions
* strip_pronouns
* strip_stopwords
* strip_sparse_terms
* strip_frequent_terms
* strip_html_tags
# Example
```julia-repl
julia> doc = StringDocument("This is a document of mine")
A StringDocument{String}
 * Language: Languages.English()
 * Title: Untitled Document
 * Author: Unknown Author
 * Timestamp: Unknown Time
 * Snippet: This is a document of mine
julia> prepare!(doc, strip_pronouns | strip_articles)
julia> text(doc)
"This is   document of "
```
"""
function prepare!(crps::Corpus, flags::UInt32; skip_patterns = Set{AbstractString}(), skip_words = Set{AbstractString}())
    ((flags & strip_sparse_terms) > 0) && union!(skip_words, sparse_terms(crps))
    ((flags & strip_frequent_terms) > 0) && union!(skip_words, frequent_terms(crps))

    ((flags & strip_corrupt_utf8) > 0) && remove_corrupt_utf8!(crps)
    ((flags & strip_case) > 0) && remove_case!(crps)
    ((flags & strip_html_tags) > 0) && remove_html_tags!(crps)

    lang = language(crps.documents[1])   # assuming all documents are of the same language - practically true
    r = _build_regex(lang, flags, skip_patterns, skip_words)
    !isempty(r.pattern) && remove_patterns!(crps, r)

    ((flags & strip_whitespace) > 0) && remove_whitespace!(crps)

    ((flags & stem_words) > 0) && stem!(crps)
    ((flags & tag_part_of_speech) > 0) && tag_pos!(crps)
    nothing
end

function prepare!(d::AbstractDocument, flags::UInt32; skip_patterns = Set{AbstractString}(), skip_words = Set{AbstractString}())
    ((flags & strip_corrupt_utf8) > 0) && remove_corrupt_utf8!(d)
    ((flags & strip_case) > 0) && remove_case!(d)
    ((flags & strip_html_tags) > 0) && remove_html_tags!(d)

    r = _build_regex(language(d), flags, skip_patterns, skip_words)
    !isempty(r.pattern) && remove_patterns!(d, r)
    ((flags & strip_whitespace) > 0) && remove_whitespace!(d)

    ((flags & stem_words) > 0) && stem!(d)
    ((flags & tag_part_of_speech) > 0) && tag_pos!(d)
    nothing
end

"""
    remove_whitespace(str)
Squash multiple whitespaces to a single one.
And remove all leading and trailing whitespaces.
See also: [`remove_whitespace!`](@ref)
"""
remove_whitespace(str::AbstractString) = replace(strip(str), r"\s+"=>" ")

"""
    remove_whitespace!(doc)
    remove_whitespace!(crps)
Squash multiple whitespaces to a single space and remove all leading and trailing whitespaces in document or crps.
Does no-op for `FileDocument`, `TokenDocument` or `NGramDocument`.
See also: [`remove_whitespace`](@ref)
"""
function remove_whitespace!(d::StringDocument)
    d.text = remove_whitespace(d.text)
end

function remove_whitespace!(crps::Corpus)
    for doc in crps
        remove_whitespace!(doc)
    end
end

function remove_whitespace!(d::AbstractDocument)
    nothing
end

"""
    remove_patterns(str, rex::Regex)
Remove the part of str matched by rex.
See also: [`remove_patterns!`](@ref)
"""
function remove_patterns(s::AbstractString, rex::Regex)
    return replace(s, rex => "")
end

function remove_patterns(s::SubString{T}, rex::Regex) where T <: String
    iob = IOBuffer()
    ioffset = s.offset
    data = codeunits(s.string)
    ibegin = 1
    for m in eachmatch(rex, s)
        len = m.match.offset-ibegin
	next = nextind(s, lastindex(m.match)+m.match.offset)
        if len > 0
            write(iob, SubString(s, ibegin, ibegin+len))
            if  next != length(s)+1
            	write(iob, ' ')
	    end
        end
        ibegin = next
    end
    len = lastindex(s) - ibegin
    (len > 0) && write(iob, SubString(s, ibegin, ibegin+len))
    String(take!(iob))
end

"""
    remove_patterns!(doc, rex::Regex)
    remove_patterns!(crps, rex::Regex)
Remove patterns matched by `rex` in document or Corpus.
Does not modify `FileDocument` or Corpus containing `FileDocument`.
See also: [`remove_patterns`](@ref)
"""
remove_patterns!(d::FileDocument, rex::Regex) = error("FileDocument cannot be modified")

function remove_patterns!(d::StringDocument, rex::Regex)
    d.text = remove_patterns(d.text, rex)
    nothing
end

function remove_patterns!(d::TokenDocument, rex::Regex)
    for i in 1:length(d.tokens)
        d.tokens[i] = remove_patterns(d.tokens[i], rex)
    end
end

function remove_patterns!(d::NGramDocument, rex::Regex)
    new_ngrams = Dict{AbstractString, Int}()
    for token in keys(d.ngrams)
        new_token = remove_patterns(token, rex)
        if haskey(new_ngrams, new_token)
            new_ngrams[new_token] = new_ngrams[new_token] + 1
        else
            new_ngrams[new_token] = 1
        end
    end
    d.ngrams = new_ngrams
    nothing
end

function remove_patterns!(crps::Corpus, rex::Regex)
    for doc in crps
        remove_patterns!(doc, rex)
    end
end

##
# internal helper methods

_build_regex(lang, flags::UInt32) = _build_regex(lang, flags, Set{AbstractString}(), Set{AbstractString}())
_build_regex(lang, flags::UInt32, patterns::Set{T}, words::Set{T}) where {T <: AbstractString} = _combine_regex(_build_regex_patterns(lang, flags, patterns, words))

function _combine_regex(regex_parts::Set{T}) where T <: AbstractString
    l = length(regex_parts)
    (0 == l) && return r""
    (1 == l) && return mk_regex(pop!(regex_parts))

    iob = IOBuffer()
    write(iob, "($(pop!(regex_parts)))")
    for part in regex_parts
        write(iob, "|($part)")
    end
    mk_regex(String(take!(iob)))
end

function _build_regex_patterns(lang, flags::UInt32, patterns::Set{T}, words::Set{T}) where T <: AbstractString
    #((flags & strip_whitespace) > 0) && push!(patterns, "\\s+")
    if (flags & strip_non_letters) > 0
        push!(patterns, "[^a-zA-Z\\s]")
    else
        ((flags & strip_punctuation) > 0) && push!(patterns, "[-.,:;,!?'\"\\[\\]\\(\\)\\{\\}|\\`#\$%@^&*_+<>“”—’‘/]+")
        ((flags & strip_numbers) > 0) && push!(patterns, "\\d+")
    end
    if (flags & strip_articles) > 0
        union!(words, articles(lang))
    else
        ((flags & strip_indefinite_articles) > 0) && union!(words, indefinite_articles(lang))
        ((flags & strip_definite_articles) > 0) && union!(words, definite_articles(lang))
    end
    ((flags & strip_prepositions) > 0) && union!(words, prepositions(lang))
    ((flags & strip_pronouns) > 0) && union!(words, pronouns(lang))
    ((flags & strip_stopwords) > 0) && union!(words, stopwords(lang))

    words_pattern = _build_words_pattern(collect(words))
    !isempty(words_pattern) && push!(patterns, words_pattern)
    patterns
end

function _build_words_pattern(words::Vector{T}) where T <: AbstractString
    isempty(words) && return ""

    iob = IOBuffer()
    write(iob, "\\b(")
    write(iob, words[1])
    l = length(words)
    for idx in 2:l
        write(iob, '|')
        write(iob, words[idx])
    end
    write(iob, ")\\b")
    String(take!(iob))
end


text:
"""
    Env

A WSTP library environment.

# External links
 - [WSENV](https://reference.wolfram.com/language/ref/c/WSENV.html)
"""
mutable struct Env
    ptr::Ptr{Cvoid}
end

const CEnv = Ptr{Cvoid}
Base.cconvert(::Type{CEnv}, env::Env) = env.ptr

# global environment object
const env = Env(C_NULL)
# reference counter to ensure all links are closed before deinitialize
const REFCOUNT = Threads.Atomic{Int}(1)

# this should be called when opening a link
function refcount_inc()
    Threads.atomic_add!(REFCOUNT, 1)
end

# this should be called in the finalizer of a link
function refcount_dec()
    # refcount zero, all objects finalized, now finalize MathLink
    if Threads.atomic_sub!(REFCOUNT, 1) == 1
        # void WSDeinitialize(WSENV env)
        ccall((:MLDeinitialize, mlib), Cvoid, (CEnv,), env)
        env.ptr = C_NULL
    end
end

function __init__()
    # WSENV WSInitialize(WSEnvironmentParameter p)
    env.ptr = ccall((:MLInitialize, mlib), CEnv, (Ptr{Cvoid},), C_NULL)
    if env.ptr == C_NULL
        error("Could not initialize MathLink library")
    end
    atexit(refcount_dec)
end


text:
module TestStepSize
using LinearAlgebra
using Test
using OptimizationAlgorithms
using OptimizationAlgorithms: fixedpoint!, StoppingCriterion
using ForwardDiff: gradient, derivative

n = 2
A = randn(n, n)
A = A'A + I
μ = randn(n)
f(x) = (x-μ)'A*(x-μ)

ε = 1e-6
maxiter = 1024
@testset "descent" begin
    x = randn(n)
    G = OptimizationAlgorithms.Gradient(f, x)
    U = OptimizationAlgorithms.DecreasingStep(G, x)
    fixedpoint!(U, x, StoppingCriterion(x, dx = 1e-2ε, maxiter = maxiter))
    @test norm(gradient(f, x)) < ε

    # also test scalar case
    x = 2rand()-1
    g(x) = -exp(-(x-1)^2)
    G = OptimizationAlgorithms.Gradient(g, x)
    U = OptimizationAlgorithms.DecreasingStep(G, x)
    xn, t = fixedpoint!(U, x, StoppingCriterion(x, dx = 1e-2ε, maxiter = maxiter))
    @test abs(derivative(g, xn)) < ε
end

@testset "armijo" begin
    x = randn(n)
    G = OptimizationAlgorithms.Gradient(f, x)
    U = OptimizationAlgorithms.ArmijoStep(G, x)
    fixedpoint!(U, x, StoppingCriterion(x, dx = 1e-2ε, maxiter = maxiter))
    @test norm(gradient(f, x)) < ε
end

end # TestStepSize


text:
module OperationsResearchApplications

# export key functions
export search_ora, run_ora

"""
    OperationsResearchApplications

A collection of different applications (apps) of Operations Research (OR).

These apps are tagged by both the subject area of the application (ex: manufacturing, supply chain, etc.) and the model(s) that they employ (LP, MILP, Simulation, etc.)
"""

struct App
    title::String
    code_path::String
    documentation::String
end

include("search_ora.jl")  # search for apps based on title, area tags, or model types

function run_ora(app::App)
    include(app.code_path)
end

end


text:
using Compat
using ProtoBuf
import ProtoBuf.meta
import Base: hash, isequal, ==

type __enum_EventType <: ProtoEnum
    EVENT_CREATE::Int32
    EVENT_CLOSE::Int32
    EVENT_APPEND::Int32
    EVENT_RENAME::Int32
    EVENT_METADATA::Int32
    EVENT_UNLINK::Int32
    __enum_EventType() = new(0,1,2,3,4,5)
end #type __enum_EventType
const EventType = __enum_EventType()

type __enum_INodeType <: ProtoEnum
    I_TYPE_FILE::Int32
    I_TYPE_DIRECTORY::Int32
    I_TYPE_SYMLINK::Int32
    __enum_INodeType() = new(0,1,2)
end #type __enum_INodeType
const INodeType = __enum_INodeType()

type __enum_MetadataUpdateType <: ProtoEnum
    META_TYPE_TIMES::Int32
    META_TYPE_REPLICATION::Int32
    META_TYPE_OWNER::Int32
    META_TYPE_PERMS::Int32
    META_TYPE_ACLS::Int32
    META_TYPE_XATTRS::Int32
    __enum_MetadataUpdateType() = new(0,1,2,3,4,5)
end #type __enum_MetadataUpdateType
const MetadataUpdateType = __enum_MetadataUpdateType()

type EventProto
    _type::Int32
    contents::Array{UInt8,1}
    EventProto() = (o=new(); fillunset(o); o)
end #type EventProto
const __req_EventProto = Symbol[:_type,:contents]
meta(t::Type{EventProto}) = meta(t, __req_EventProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::EventProto) = ProtoBuf.protohash(v)
isequal(v1::EventProto, v2::EventProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::EventProto, v2::EventProto) = ProtoBuf.protoeq(v1, v2)

type EventBatchProto
    txid::Int64
    events::Array{EventProto,1}
    EventBatchProto() = (o=new(); fillunset(o); o)
end #type EventBatchProto
const __req_EventBatchProto = Symbol[:txid]
meta(t::Type{EventBatchProto}) = meta(t, __req_EventBatchProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::EventBatchProto) = ProtoBuf.protohash(v)
isequal(v1::EventBatchProto, v2::EventBatchProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::EventBatchProto, v2::EventBatchProto) = ProtoBuf.protoeq(v1, v2)

type CreateEventProto
    _type::Int32
    path::AbstractString
    ctime::Int64
    ownerName::AbstractString
    groupName::AbstractString
    perms::FsPermissionProto
    replication::Int32
    symlinkTarget::AbstractString
    overwrite::Bool
    defaultBlockSize::Int64
    CreateEventProto() = (o=new(); fillunset(o); o)
end #type CreateEventProto
const __req_CreateEventProto = Symbol[:_type,:path,:ctime,:ownerName,:groupName,:perms]
const __val_CreateEventProto = @compat Dict(:defaultBlockSize => 0)
meta(t::Type{CreateEventProto}) = meta(t, __req_CreateEventProto, ProtoBuf.DEF_FNUM, __val_CreateEventProto, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::CreateEventProto) = ProtoBuf.protohash(v)
isequal(v1::CreateEventProto, v2::CreateEventProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::CreateEventProto, v2::CreateEventProto) = ProtoBuf.protoeq(v1, v2)

type CloseEventProto
    path::AbstractString
    fileSize::Int64
    timestamp::Int64
    CloseEventProto() = (o=new(); fillunset(o); o)
end #type CloseEventProto
const __req_CloseEventProto = Symbol[:path,:fileSize,:timestamp]
meta(t::Type{CloseEventProto}) = meta(t, __req_CloseEventProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::CloseEventProto) = ProtoBuf.protohash(v)
isequal(v1::CloseEventProto, v2::CloseEventProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::CloseEventProto, v2::CloseEventProto) = ProtoBuf.protoeq(v1, v2)

type AppendEventProto
    path::AbstractString
    newBlock::Bool
    AppendEventProto() = (o=new(); fillunset(o); o)
end #type AppendEventProto
const __req_AppendEventProto = Symbol[:path]
const __val_AppendEventProto = @compat Dict(:newBlock => false)
meta(t::Type{AppendEventProto}) = meta(t, __req_AppendEventProto, ProtoBuf.DEF_FNUM, __val_AppendEventProto, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::AppendEventProto) = ProtoBuf.protohash(v)
isequal(v1::AppendEventProto, v2::AppendEventProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::AppendEventProto, v2::AppendEventProto) = ProtoBuf.protoeq(v1, v2)

type RenameEventProto
    srcPath::AbstractString
    destPath::AbstractString
    timestamp::Int64
    RenameEventProto() = (o=new(); fillunset(o); o)
end #type RenameEventProto
const __req_RenameEventProto = Symbol[:srcPath,:destPath,:timestamp]
meta(t::Type{RenameEventProto}) = meta(t, __req_RenameEventProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::RenameEventProto) = ProtoBuf.protohash(v)
isequal(v1::RenameEventProto, v2::RenameEventProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::RenameEventProto, v2::RenameEventProto) = ProtoBuf.protoeq(v1, v2)

type MetadataUpdateEventProto
    path::AbstractString
    _type::Int32
    mtime::Int64
    atime::Int64
    replication::Int32
    ownerName::AbstractString
    groupName::AbstractString
    perms::FsPermissionProto
    acls::Array{AclEntryProto,1}
    xAttrs::Array{XAttrProto,1}
    xAttrsRemoved::Bool
    MetadataUpdateEventProto() = (o=new(); fillunset(o); o)
end #type MetadataUpdateEventProto
const __req_MetadataUpdateEventProto = Symbol[:path,:_type]
meta(t::Type{MetadataUpdateEventProto}) = meta(t, __req_MetadataUpdateEventProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::MetadataUpdateEventProto) = ProtoBuf.protohash(v)
isequal(v1::MetadataUpdateEventProto, v2::MetadataUpdateEventProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::MetadataUpdateEventProto, v2::MetadataUpdateEventProto) = ProtoBuf.protoeq(v1, v2)

type UnlinkEventProto
    path::AbstractString
    timestamp::Int64
    UnlinkEventProto() = (o=new(); fillunset(o); o)
end #type UnlinkEventProto
const __req_UnlinkEventProto = Symbol[:path,:timestamp]
meta(t::Type{UnlinkEventProto}) = meta(t, __req_UnlinkEventProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::UnlinkEventProto) = ProtoBuf.protohash(v)
isequal(v1::UnlinkEventProto, v2::UnlinkEventProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::UnlinkEventProto, v2::UnlinkEventProto) = ProtoBuf.protoeq(v1, v2)

type EventsListProto
    events::Array{EventProto,1}
    firstTxid::Int64
    lastTxid::Int64
    syncTxid::Int64
    batch::Array{EventBatchProto,1}
    EventsListProto() = (o=new(); fillunset(o); o)
end #type EventsListProto
const __req_EventsListProto = Symbol[:firstTxid,:lastTxid,:syncTxid]
meta(t::Type{EventsListProto}) = meta(t, __req_EventsListProto, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, true, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES)
hash(v::EventsListProto) = ProtoBuf.protohash(v)
isequal(v1::EventsListProto, v2::EventsListProto) = ProtoBuf.protoisequal(v1, v2)
==(v1::EventsListProto, v2::EventsListProto) = ProtoBuf.protoeq(v1, v2)

export EventType, INodeType, MetadataUpdateType, EventProto, EventBatchProto, CreateEventProto, CloseEventProto, AppendEventProto, RenameEventProto, MetadataUpdateEventProto, UnlinkEventProto, EventsListProto


text:
##########################################################
# Functions related to plotting case series from SIR model
##########################################################

using PyCall, PyPlot, Plots
pygui(:tk)

"""
plotEvolution(N::Float64,states::Matrix{Float64},
    d::Dynamics;log10::Bool=false)

Generate a plot of states evolved for the SIR model.

- N: the population size
- states: a matrix of states evoloved over time. See the function, 'evolve'.
- d: a type of Dynamics of the SIR model
- log10: a keyword argument (default is false) of scaling state values by log10

"""

function plotEvolution(N::Float64,states::Matrix{Float64},
    d::Dynamics;log10::Bool=false)
    if (log10==false)
        Plots.plot(N*Matrix(states'),label=stateNames(d))
    else
        Plots.plot(N*Matrix(states'),label=stateNames(d),ylab=:log10)
    end
end

"""
plotfit(fit::CaseModelFitResult,nt::Int64=0)

- fit: output from fitCaseModel
- nt: how long we want to extend predictions (use zero for not extrapolating)
"""
function plotFit(fit::CaseModelFitResult,nt::Int64=0)
    if (nt==0)
        nt = fit.inputs.nt
    end
    Plots.plot(fit.inputs.C,yaxis=:log,seriestype=:scatter,
                color=:black,label="actual")
    Plots.plot!(estimatedStates(fit,nt)[!,:I],
                yaxis=:log,label="infected")
    Plots.plot!(fitted(fit,nt),yaxis=:log, label="fitted",color=:blue)
end

"""
pyplotEvolution(N::Float64, states::Matrix{Float64},
              d::Dynamics;  log10 = "Null", figTitle = "",
              grid = "off")

Generate a plot of states evolved for the SIR model.

# Arguments
- N: the population size
- states: a matrix of states evoloved over time. See the function, `evolve`.
- d: a type of Dynamics of the SIR model
- log10: a keyword argument (default is "Null") where "semilog" makes a plot
          with log scaling on the y axis, "loglog" makes a plot with log scaling
          on the x axis and y axis, and "Null"makes a plot without log scaling.
- figTitle: a keyword argument (default is empty) that add a title to the figure.
- grid: a keyword argument (default is `false`) that adds a grid on the plotting if
         it is `true`.
- fsize: a keyword argument (default is [5, 5]]) that set the size of the figure.
- savefigure: a keyword argument (default is empty "") that save the figure in
               a given path.(e.g., savefigure = "../images/Plot.png") The extension
               of the filename indicates the format that is used to save the figure
               such as .pdf, .png, .jpg,...
"""

function pyPlotEvolution(N::Float64, states::Matrix{Float64},
                         d::Dynamics; log10 = "Null", figTitle = "",
                         grid = false, fsize= [10, 7], savefigure = "")

    fig, (ax) = PyPlot.subplots(nrows=1, figsize=(fsize[1], fsize[2]))

    ###############
    #  PLOT DATA  #
    ###############

    if log10 == "semilog"
        ax.semilogy(N*Matrix(states'))
    elseif log10 == "loglog"
        ax.loglog(N*Matrix(states'))
    elseif log10 == "Null"
        ax.plot(N*Matrix(states'))
    else
        println("Incorrect log10 option, no log is applied.")
        ax.plot(N*Matrix(states'))
    end

    #####################
    #  AXES AND LABELS  #
    #####################

    ax.spines["top"].set_visible(false) # Hide the top edge of the axis
    ax.spines["right"].set_visible(false) # Hide the right edge of the axis
    ax.xaxis.set_ticks_position("bottom")
    # ax1.spines["left"].set_position(("axes",-0.03)) # Offset the left scale from the axis
    # ax1.spines["bottom"].set_position(("axes",-0.05)) # Offset the bottom scale from the axis

    ax.set_xlabel("Time")
    ax.set_ylabel("Confirmed")

    # axis("tight")

    PyPlot.grid(grid)
    PyPlot.legend(stateNames(d)[:], loc="center right", fancybox="true")
    PyPlot.title(figTitle)

    if !isempty(savefigure)
        PyPlot.savefig(savefigure)
    end
    PyPlot.display_figs()

end

"""
pyplotFit(fit::CaseModelFitResult, nt::Int64 = 0;
        log10 = "Null", figTitle = "", grid = false,
        fsize= [5, 5], savefigure = "", )

Generate a plot of fitted state evolution for the SIR model.

# Arguments
- pop: a type SIRXPopulation that contains population size
- states: a matrix of states evoloved over time. See the function, `evolve`.
- d: a type of Dynamics of the SIR model
- actualState: actual case numbers
- log10: a keyword argument (default is "Null") where "semilog" makes a plot
          with log scaling on the y axis, "loglog" makes a plot with log scaling
          on the x axis and y axis, and "Null"makes a plot without log scaling.
- figTitle: a keyword argument (default is empty) that add a title to the figure.
- grid: a keyword argument (default is `false`) that adds a grid on the plotting if
         it is `true`.
- fsize: a keyword argument (default is [5, 5]]) that set the size of the figure.
- savefigure: a keyword argument (default is empty "") that save the figure in
               a given path.(e.g., savefigure = "../images/Plot.png") The extension
               of the filename indicates the format that is used to save the figure
               such as .pdf, .png, .jpg,...
"""

function pyplotFit(fit::CaseModelFitResult, nt::Int64 = 0;
                 log10 = "Null", figTitle = "", grid = false,
                 fsize= [10, 7], savefigure = "", ref = [0])

    fig, (ax) = PyPlot.subplots(nrows=1, figsize=(fsize[1], fsize[2]))

    ###############
    #  PLOT DATA  #
    ###############

    if (nt==0)
        nt = fit.inputs.nt
    end

    if log10 == "semilog"
        # ax.semilogy(fit.inputs.C, linestyle ="-", color = "#004477",
        #           marker=".", markersize = 15, markerfacecolor= "none",
        #           label= "actual")
        ax.semilogy(fit.inputs.C, linestyle ="", color = "#004477",
                marker=".", markersize = 15, markerfacecolor= "#004477",
                label= "actual")
        ax.semilogy(estimatedStates(fit, nt)[!, :I], linestyle ="--", color = "red",
                  label= "infected")
        ax.semilogy(fitted(fit, nt), label="fitted")
        if (length(ref) > 1)
            ax.semilogy(ref, color = "green",  label="paper")
        end
        xlim(0, nt)

    elseif log10 == "loglog"
        ax.loglog(fit.inputs.C, linestyle ="", color = "#004477",
                  marker=".", markersize = 15, markerfacecolor= "#004477",
                  label= "actual")
        ax.loglog(estimatedStates(fit, nt)[!, :I], linestyle ="--", color = "red",
                  label= "infected")
        ax.loglog(fitted(fit, nt), label="fitted")
        if (length(ref) > 1)
            ax.semilogy(ref, color = "green",  label="paper")
        end
        # xlim(0, log10(nt+2))

    elseif log10 == "null"
        ax.plot(fit.inputs.C, linestyle ="", color = "#004477",
                  marker=".", markersize = 15, markerfacecolor= "#004477",
                  label= "actual")
        ax.plot(estimatedStates(fit, nt)[!, :I], linestyle ="--", color = "red",
                  label= "infected")
        ax.plot(fitted(fit, nt), label="fitted")
        if (length(ref) > 1)
            ax.semilogy(ref, color = "green",  label="paper")
        end
        xlim(0, nt)
    else
        println("Incorrect log10 option, no log is applied.")
        ax.plot(fit.inputs.C, linestyle ="", color = "#004477",
                  marker=".", markersize = 15, markerfacecolor= "#004477",
                  label= "actual")
        ax.plot(estimatedStates(fit, nt)[!, :I], linestyle ="--", color = "red",
                  label= "infected")
        ax.plot(fitted(fit, nt), label="fitted")
        if (length(ref) > 1)
            ax.semilogy(ref, color = "green",  label="paper")
        end
        xlim(0, nt)
    end

    #####################
    #  AXES AND LABELS  #
    #####################

    ax.spines["top"].set_visible(false) # Hide the top edge of the axis
    ax.spines["right"].set_visible(false) # Hide the right edge of the axis
    ax.xaxis.set_ticks_position("bottom")

    # ax.spines["left"].set_position(("axes",-0.03)) # Offset the left scale from the axis
    # ax.spines["bottom"].set_position(("axes",-0.05)) # Offset the bottom scale from the axis

    ax.autoscale(false)

    ax.set_xlabel("Time")
    ax.set_ylabel("Confirmed")

    # axis("tight")

    PyPlot.grid(grid)
    PyPlot.legend( loc="lower center")
    PyPlot.title(figTitle)

    if !isempty(savefigure)
        PyPlot.savefig(savefigure)
    end
    PyPlot.display_figs()

end


text:
@testset "258.add-digits.jl" begin
    @test all(add_digits(i) == i for i in 0:9)
    @test add_digits(18) == 9
    @test add_digits(10) == 1
    @test add_digits(38) == 2
end


text:
module MvNormalMeanCovarianceNodeTest

using Test
using ReactiveMP
using Random

import ReactiveMP: make_node

@testset "MvNormalMeanCovarianceNode" begin

    @testset "Creation" begin
        node = make_node(MvNormalMeanCovariance)

        @test functionalform(node)          === MvNormalMeanCovariance
        @test sdtype(node)                  === Stochastic()
        @test name.(interfaces(node))       === (:out, :μ, :Σ)
        @test factorisation(node)           === ((1, 2, 3), )
        @test localmarginalnames(node)      === (:out_μ_Σ, )
        @test metadata(node)                === nothing

        node = make_node(MvNormalMeanCovariance, meta = 1)

        @test functionalform(node)          === MvNormalMeanCovariance
        @test sdtype(node)                  === Stochastic()
        @test name.(interfaces(node))       === (:out, :μ, :Σ)
        @test factorisation(node)           === ((1, 2, 3), )
        @test localmarginalnames(node)      === (:out_μ_Σ, )
        @test metadata(node)                === 1

        node = make_node(MvNormalMeanCovariance, factorisation = ((1, ), (2, 3)))

        @test functionalform(node)          === MvNormalMeanCovariance
        @test sdtype(node)                  === Stochastic()
        @test name.(interfaces(node))       === (:out, :μ, :Σ)
        @test factorisation(node)           === ((1, ), (2, 3, ))
        @test localmarginalnames(node)      === (:out, :μ_Σ, )
        @test metadata(node)                === nothing
    end

    @testset "AverageEnergy" begin

        begin
            q_out = PointMass([ 1.0, 1.0 ])
            q_μ   = MvNormalMeanCovariance([ 1.0, 1.0 ], [ 1.0 0.0; 0.0 1.0 ])
            q_Σ   = PointMass([ 2.0 0.0; 0.0 2.0 ])

            for N in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision)
                marginals = (Marginal(q_out, false, false), Marginal(convert(N, q_μ), false, false), Marginal(q_Σ, false, false))
                @test score(AverageEnergy(), MvNormalMeanCovariance, Val{(:out, :μ, :Σ)}, marginals, nothing) ≈ 3.0310242469692907
            end
        end

        begin
            q_out = MvNormalMeanCovariance([0.34777017741128646, 0.478300212208703], [0.36248838736696753 0.42718825177834396; 0.42718825177834396 0.6613345138178703])
            q_μ   = MvNormalMeanCovariance([0.28114649146383175, 0.7773140083754762], [0.5753219837395999 0.57742814389364; 0.57742814389364 0.7164979237943061])
            q_Σ   = PointMass([0.9575163948991356 0.952247049670764; 0.952247049670764 1.001998675194526])

            for N1 in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), N2 in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision)
                marginals = (Marginal(convert(N1, q_out), false, false), Marginal(convert(N2, q_μ), false, false), Marginal(q_Σ, false, false))
                @test score(AverageEnergy(), MvNormalMeanCovariance, Val{(:out, :μ, :Σ)}, marginals, nothing) ≈ 4.863921806562866
            end
        end

        begin
            q_out = MvNormalMeanCovariance([0.3685447210709054, 0.7059290804186025], [0.47091887707487257 0.1715872501382935; 0.1715872501382935 0.07937045688177925])
            q_μ   = MvNormalMeanCovariance([0.9821121935372943, 0.8740184196009864], [1.2740987247403437 0.376849163380472; 0.376849163380472 0.17831073369063927])
            q_Σ   = PointMass([0.761731832494188 0.21222845781825198; 0.21222845781825198 0.5989429144303828])

            for N1 in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision), N2 in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision)
                marginals = (Marginal(convert(N1, q_out), false, false), Marginal(convert(N2, q_μ), false, false), Marginal(q_Σ, false, false))
                @test score(AverageEnergy(), MvNormalMeanCovariance, Val{(:out, :μ, :Σ)}, marginals, nothing) ≈ 2.867156770180465
            end
        end

        begin
            q_out_μ = MvNormalMeanCovariance([0.7110546689402384, 0.2550803989358359, 0.7284714504680616, 0.4523569391382931], [0.9354883429919909 0.979541773915805 0.11235208550731414 0.720928227797172; 0.979541773915805 1.5518339076149474 0.27953313117565737 0.9410472233142694; 0.11235208550731414 0.27953313117565737 0.09215639567249542 0.15210811270662963; 0.720928227797172 0.9410472233142694 0.15210811270662963 0.7114840285830272])
            q_Σ   = PointMass([0.6588719526947034 0.3784788140113857; 0.3784788140113857 0.2611292983123214])

            for N in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision)
                marginals = (Marginal(convert(N, q_out_μ), false, false), Marginal(q_Σ, false, false))
                @test score(AverageEnergy(), MvNormalMeanCovariance, Val{(:out_μ, :Σ)}, marginals, nothing) ≈ 6.741420950973408
            end
        end

        begin
            q_out_μ = MvNormalMeanCovariance([0.6273836373763129, 0.2611455272882084, 0.15015693447418932, 0.7907870618820343], [2.4624395254748563 1.3049173253859347 2.2581875619221057 1.8418052081652991; 1.3049173253859347 1.0230726719658598 1.0344117375112738 0.8684398658701927; 2.2581875619221057 1.0344117375112738 2.164616392026885 1.7042534061619456; 1.8418052081652991 0.8684398658701927 1.7042534061619456 1.5840784999702135])
            q_Σ   = PointMass([0.9918745791319346 0.7102392994361006; 0.7102392994361006 0.5209829748135837])

            for N in (MvNormalMeanPrecision, MvNormalMeanCovariance, MvNormalWeightedMeanPrecision)
                marginals = (Marginal(convert(N, q_out_μ), false, false), Marginal(q_Σ, false, false))
                @test score(AverageEnergy(), MvNormalMeanCovariance, Val{(:out_μ, :Σ)}, marginals, nothing) ≈ 60.075730792149585
            end
        end

    end

end
end