text:
@testset "NGramIterator and friends" begin
    for s in [randstring(100) for _ in 1:10]
        it = NGramIterator(s)
        @test length(it) == 100 + 3 - 1
        c = collect(it)
        @test c == ngrams(s)
        @test c == NGramIterator(codeunits(s)) |> collect
    end
end

@testset "NGramMatrix construction types and basics" begin
    M1 = NGramMatrix(["hello", "world"])
    M2 = NGramMatrix(["a part", "is", missing] |> PooledArray)
    M3 = NGramMatrix([missing, missing])
    @test M1 isa NGramMatrix{String}
    @test M1 isa AbstractMatrix{Int64}
    @test M2 isa NGramMatrix{Union{String, Missing}}
    @test M2 isa AbstractMatrix{Union{Int64, Missing}}
    @test M3 isa NGramMatrix{Missing}
    @test M3 isa AbstractMatrix{Missing}

    @test size(M1) == (2053, 2)
    @test size(M2) == (2053, 3)
    @test size(M3) == (2053, 2)

    M1 = NGramMatrix([codeunits("hello"), codeunits("world")])
    M2 = NGramMatrix([codeunits("a part"), codeunits("is"), missing])
    M3 = NGramMatrix([missing, missing] |> PooledArray)
    @test M1 isa NGramMatrix{<:CodeUnits}
    @test M1 isa AbstractMatrix{Int64}
    @test M2 isa NGramMatrix{<:Union{CodeUnits, Missing}}
    @test M2 isa AbstractMatrix{Union{Int64, Missing}}
    @test M3 isa NGramMatrix{Missing}
    @test M3 isa AbstractMatrix{Missing}

    @test size(M1) == (2053, 2)
    @test size(M2) == (2053, 3)
    @test size(M3) == (2053, 2)

    M1 = NGramMatrix([[1,2,3], [4,5,6]] |> PooledArray)
    M2 = NGramMatrix([Int[], [0], missing])
    M3 = NGramMatrix([missing, missing])
    @test M1 isa NGramMatrix{Vector{Int64}}
    @test M1 isa AbstractMatrix{Int64}
    @test M2 isa NGramMatrix{Union{Vector{Int64}, Missing}}
    @test M2 isa AbstractMatrix{Union{Int64, Missing}}
    @test M3 isa NGramMatrix{Missing}
    @test M3 isa AbstractMatrix{Missing}

    @test size(M1) == (2053, 2)
    @test size(M2) == (2053, 3)
    @test size(M3) == (2053, 2)
end

@testset "Indexing" begin
    S1 = ["hello", "world"]
    S2 = ["a part", "is", missing] |> PooledArray
    S3 = [missing, missing]
    M1 = NGramMatrix(S1)
    M2 = NGramMatrix(S2)
    M3 = NGramMatrix(S3)

    for i in eachindex(S1)
        @test M1[:, i] == NGramMatrix(S1[i])
    end
    for i in eachindex(S2)
        @test isequal(M2[:, i], NGramMatrix([S2[i]] |> PooledArray))
    end
    for i in eachindex(S3)
        @test isequal(M3[:, i], NGramMatrix(S3[i]))
    end

    @test M1[:, [1,2]] == NGramMatrix(S1[[1,2]])
    @test isequal(M2[:, [3,2]], NGramMatrix(S2[[3,2]]))
    @test isequal(M3[:, [1]], NGramMatrix(S3[[1]]))

    @test M1[:, :] == M1
    @test isequal(M2[:, :], M2)
    @test isequal(M3[:, :], M3)
    @test M1[:, :] !== M1
    @test M2[:, :] !== M2
    @test M3[:, :] !== M3
    @test hash(M1[:, :]) == hash(M1)
    @test hash(M2[:, :]) == hash(M2)
    @test hash(M3[:, :]) == hash(M3)

    @test_throws MethodError M1[1, 2]
    @test_throws MethodError M1[1, :]
    @test_throws MethodError M1[[1,2], 2]
    @test_throws MethodError M1[[2,1], :]

    @test_throws BoundsError M1[:, 0]
    @test_throws BoundsError M1[:, 3]
    @test_throws BoundsError M1[:, [1,0]]
    @test_throws BoundsError M1[:, [1,2,3]]
end

@testset "hcat" begin
    S1 = ["hello", "world"]
    S2 = ["a part", "is", missing]
    S3 = [missing, missing]
    M1 = NGramMatrix(S1)
    M2 = NGramMatrix(S2)
    M3 = NGramMatrix(S3)

    @test areequal(hcat(M1, M2), reduce(hcat, [M1, M2]), NGramMatrix(vcat(S1, S2)))
    @test areequal(hcat(M1, M3), reduce(hcat, [M1, M3]), NGramMatrix(vcat(S1, S3)))
    @test areequal(hcat(M2, M3), reduce(hcat, [M2, M3]), NGramMatrix(vcat(S2, S3)))
    @test areequal(hcat(M1, M2, M3), reduce(hcat, [M1, M2, M3]), NGramMatrix(vcat(S1, S2, S3)))

    @inferred hcat(M1, M1)
    @inferred reduce(hcat, [M1, M1])
    @inferred hcat(M2, M2)
    @inferred reduce(hcat, [M2, M2])
    @inferred hcat(M3, M3)
    @inferred reduce(hcat, [M3, M3])
    @inferred hcat(M1, M2)
    @inferred reduce(hcat, [M1, M2])
    @inferred hcat(M2, M3)
    @inferred reduce(hcat, [M2, M3])
    @inferred hcat(M3, M1)
    @inferred reduce(hcat, [M3, M1])
    @inferred hcat(M1, M2, M3)
    @inferred reduce(hcat, [M1, M2, M3])

    @test_throws DimensionMismatch hcat(NGramMatrix(S1, 1), NGramMatrix(S1, 2))
    @test_throws DimensionMismatch hcat(NGramMatrix(S1, 2, 256), NGramMatrix(S1, 2, 128))
    @test_throws DimensionMismatch hcat(NGramMatrix(S1, 2, 256, 7), NGramMatrix(S1, 2, 256, 2053))
    @test_throws DimensionMismatch reduce(hcat, [NGramMatrix(S1, 1), NGramMatrix(S1, 2)])
    @test_throws DimensionMismatch reduce(hcat, [NGramMatrix(S1, 2, 256), NGramMatrix(S1, 2, 128)])
    @test_throws DimensionMismatch reduce(hcat, [NGramMatrix(S1, 2, 256, 7), NGramMatrix(S1, 2, 256, 2053)])
end

@testset "ngram computation" begin
    x = [1,3,5,2,6,8,3]
    b = 8 + 1

    slicer = (x,n) -> map(i -> vcat(fill(Mill.string_start_code(), max(0, n-i)),
                                    x[max(i-n+1,1):min(i,length(x))],
                                    fill(Mill.string_end_code(), max(0, i-length(x)))),
                                    1:length(x)+n-1)
    indexes = (x,b) -> mapreduce(i -> i[2] * b^(i[1]-1), +, enumerate(reverse(x)))
    function idx2vec(i,n)
        o = zeros(Int,n)
        for v in i
            o[mod(v,n)+1] +=1
        end
        o
    end

    @testset "testing ngrams on vector of Ints" begin
        @test ngrams(x,3,b) == map(x -> indexes(x,b),slicer(x,3))
        @test ngrams(x,2,b) == map(x -> indexes(x,b),slicer(x,2))
        @test ngrams(x,1,b) == map(x -> indexes(x,b),slicer(x,1))
    end

    @testset "testing frequency of ngrams on vector of Ints and on Strings" begin
        @test countngrams(x,3,b,10) == idx2vec(map(x -> indexes(x,b), slicer(x,3)), 10)
        for s in split("Lorem ipsum dolor sit amet, consectetur adipiscing elit")
            @test countngrams(s,3,256,10) == idx2vec(ngrams(s,3,256),10)
        end

        s = split("Lorem ipsum dolor sit amet, consectetur adipiscing elit")
        @test countngrams(s,3,256,10) == hcat(map(ss -> idx2vec(ngrams(ss,3,256), 10),s)...)
    end
end

@testset "NGramMatrix to SparseMatrix" begin
    for (n, m) in product([2,3,5], [10, 100, 1000])
        b = 256
        S = ["hello", "world", "!!!"]
        Si = map(codeunits, S)
        Sc = map(i -> Int.(i), Si)
        B = NGramMatrix(S, n, b, m)
        Bi = NGramMatrix(Si, n, b, m)
        Bc = NGramMatrix(Sc, n, b, m)
        @test SparseMatrixCSC(B) == SparseMatrixCSC(Bi) == SparseMatrixCSC(Bc)
        @test SparseMatrixCSC(countngrams(S, n, b, m)) == SparseMatrixCSC(B)
        @test SparseMatrixCSC(countngrams(Si, n, b, m)) == SparseMatrixCSC(Bi)
        @test SparseMatrixCSC(countngrams(Sc, n, b, m)) == SparseMatrixCSC(Bc)
    end
end

@testset "NGramMatrix multiplication" begin
    for (n, m) in product([2,3,5], [10, 100, 1000])
        b = 256
        A = randn(10, m)

        S = [randstring(100) for _ in 1:10]
        Si = map(codeunits, S)
        Sc = map(i -> Int.(i), Si)
        B = NGramMatrix(S, n, b, m)
        Bi = NGramMatrix(Si, n, b, m)
        Bc = NGramMatrix(Sc, n, b, m)

        @test all(A * B ≈ A * countngrams(S, n, b, m))
        @test all(A * Bi ≈ A * B)
        @test all(A * Bc ≈ A * B)

        @test all(A * B ≈ A * SparseMatrixCSC(B))
        @test all(A * B ≈ A * SparseMatrixCSC(Bi))
        @test all(A * B ≈ A * SparseMatrixCSC(Bc))

        @inferred A * B
        @inferred A * Bi
        @inferred A * Bc

        @test eltype(A * B) == eltype(A)
        @test eltype(A * Bi) == eltype(A)
        @test eltype(A * Bc) == eltype(A)
    end
end

@testset "NGramMatrix multiplication gradtest" begin
    for (n, m) in product([2,3,5], [10, 20])
        b = 256
        A = randn(10, m)

        s = [randstring(100) for _ in 1:10] |> PooledArray
        si = map(codeunits, s)
        sc = map(i -> Int.(i), si)
        Ns = countngrams(s, n, b, m)
        Nsi = countngrams(si, n, b, m)
        Nsc = countngrams(sc, n, b, m)
        B = NGramMatrix(s, n, b, m)
        Bi = NGramMatrix(si, n, b, m)
        Bc = NGramMatrix(sc, n, b, m)
        Ss = SparseMatrixCSC(B)
        Ssi = SparseMatrixCSC(Bi)
        Ssc = SparseMatrixCSC(Bc)

        f = gradf(*, A, B)
        dNs = gradient(f, A, Ns)[1]
        dSs = gradient(f, A, Ss)[1]

        dA, dB = gradient(f, A, B)
        @test dA ≈ dNs
        @test dA ≈ dSs
        @test isnothing(dB)
        @test gradtest(A -> A * B, A)

        dA, dBi = gradient(f, A, Bi)
        @test dA ≈ dNs
        @test dA ≈ dSs
        @test isnothing(dBi)
        @test gradtest(A -> A * Bi, A)

        dA, dBc = gradient(f, A, Bc)
        @test dA ≈ dNs
        @test dA ≈ dSs
        @test isnothing(dBi)
        @test gradtest(A -> A * Bc, A)
    end
end

@testset "NGramMatrix multiplication second derivative" begin
    for (n, m) in product([2, 3], [10, 20])
        S = [randstring(10) for _ in 1:10]
        B = NGramMatrix(S, n, 256, m)
        A = randn(10, m)
        f = gradf(A -> sin.(A * B), A)
        df = gradf(A -> gradient(f, A)[1], A)
        @test gradtest(df, A)
    end
end

@testset "integration with MiLL & Flux" begin
    S = ["hello", "world", "!!!"]
    Si = map(i -> Int.(codeunits(i)), S)
    for (A, S) in [(NGramMatrix(S, 3, 256, 2057), S), (NGramMatrix(Si, 3, 256, 2057), Si)]
        @test reduce(catobs, [A, A]).S == vcat(S, S)
        @test hcat(A, A).S == vcat(S, S)

        W = randn(40, 2057)
        @test gradtest(W -> W * A, W)

        n = ArrayNode(A, nothing)
        @test reduce(catobs, [n, n]).data.S == vcat(S, S)
        n = BagNode(n, [1:3], nothing)
        @test reduce(catobs, [n, n]).data.data.S == vcat(S, S)
    end
end

@testset "equals with missings" begin
    M1 = NGramMatrix(["hello", "world"])
    M2 = NGramMatrix(["a part", "is", missing])
    M3 = NGramMatrix([missing, missing])
    @test M1 == M1
    @test isequal(M1, M1)
    @test M2 != M2
    @test isequal(M2, M2)
    @test M3 != M3
    @test isequal(M3, M3)
    @test M1 != M2
    @test !isequal(M1, M2)
    @test M1 != M3
    @test !isequal(M1, M3)
    @test M2 != M3
    @test !isequal(M2, M3)
end

begin
    println("Benchmarking multiplication")
    # begin block body
    A = randn(80,2053);
    S = [randstring(10) for i in 1:1000];
    B = NGramMatrix(S, 3, 256, 2053)
    C = sparse(countngrams(S, 3, 256, size(A, 2)));
    println("A * B::NGramMatrix (This should be the fastest)");
    @btime A*B;                                                 # 526.456 μs (2002 allocations: 671.95 KiB)
    println("A * countngrams(S, 3, 256, size(A, 2))")
    @btime A * countngrams(S, 3, 256, size(A, 2));                   # 154.646 ms (3013 allocations: 16.38 MiB)
    println("A * sparse(countngrams(S, 3, 256, size(A, 2)))")
    @btime A*sparse(countngrams(S, 3, 256, size(A, 2)));           # 7.525 ms (3013 allocations: 16.57 MiB)
    print("A * C where C = sparse(countngrams(S, 3, 256, size(A, 2)));");
    @btime A*C;                                                 # 1.527 ms (2 allocations: 625.08 KiB)
end


text:
# if !@isdefined(AbstractModel)
#     include(Pkg.dir("DifferentialDynamicProgramming","src","interfaces.jl"))
# end
#
# include("gmmutils.jl")
#
#
# """
# `M, dynamics, ass, T = GMMModel(model::GMMModel, x, u, y, K [,d1])`
#
# Fits a Gaussian mixture model to a set of data. The objective is to predict `y∈ ℜ N × ny` given `x,u`
#
# `K` is the number of clusters to use
#
# `d1` is the number of dimensions to project the data down to when fitting the model.
# """
# function GMMModel(model::GMMModel, x,u, xnew, K;
#     d1     = 2size(x,2)+size(u,2),
#     doplot = false,
#     nTries = 5)::GMMModel
#
#     # Put all data in a large matrix
#     Xorig = [x u xnew]
#     X, Xmin, Xmax = normalize_data(Xorig, 0.1)
#     X     = X'
#
#     D,N   = size(X)
#     U,S,V = svd(X)
#     X̄     = diagm(S)*V'
#     X̄1    = X̄[1:d1,:]
#     U1    = U[:,1:d1]
#     # X   = [x xnew]
#     T     = AffineTransform(U1,vec(Xmin),vec(Xmax))
#
#     # Fit a GMM to determine cluster assignments
#
#     M̄         = fit_good_model(X̄1,K,nTries)
#     Q,ass     = get_assignments(M̄, X̄1)
#
#     # Fit locally linear dynamics to data within each cluster
#     dynamics = Vector{LinearDynamics}(undef,K)
#     for k = 1:K
#         dynamics[k] = fit_dynamics(xnew[ass .== k,:], x[ass .== k,:], u[ass .== k,:], λ = 1e-3)
#     end
#
#     if doplot
#         # Get predictions of the next state using the fitted models
#         xp = similar(x)
#         # dynamics = Vector{LinearDynamics}(N)
#         for n = 1:N
#             #     dynamics[n] = conditionalDynamics(M, xp_indices, vec(Xorig[n,1:14]'), dof, T)
#             xp[n,:] = predict_dynamics(dynamics[ass[n]], x[n,:], u[n,:])'
#         end
#         plot(layout=size(x,2))
#         plot!(x[2:end,:], c=:blue, lab="State")
#         plot!(xp[1:end-1,:], c=:red, lab="Predictions")
#         plot!(title="Recorded states and model predictions")
#
#         plot(u,title="Recorded control signals")
#
#         # colors = ["c","b","g","r","m","k","y","w"]
#         # matshow(Q);  title("Q - soft assignment matrix")
#
#         scatter3d(X̄1[1,:], X̄1[2,:],X̄1[3,:], zcolor=ass, markersize=10,  title="Points with color coded cluster assignment")
#
#     end
#     model.M        = M̄
#     model.dynamics = dynamics
#     model.T        = T
#     yhat = predict(model, x,u)
#     fit = nrmse(xnew,yhat)
#     println("Modelfit: ", round(fit, digits=3))
#     return model
# end
#
# function GMMModel(x, u, y, args...; kwargs...)::GMMModel
#     model = GMMModel(0,0,0)
#     GMMModel(model, x, u, y, args...; kwargs...)
#     model
# end
#
#
# function predict(model::GMMModel,x,u)
#     u[isnan(u)] = 0
#     X    = [x u]'
#     Tr   = strip(model.T,1:size(X,1))
#     ass  = get_assignments_robust(model.M, X, Tr)
#     T,n  = size(x)
#     xnew = zeros(T,n)
#     for t = 1:T
#         dyn       = model.dynamics[ass[t]]
#         xnew[t,:] = dyn.fx*x[t,:] + dyn.fu*u[t,:]
#     end
#     return xnew
# end
#
# function predict(model::GMMModel,x::Vector,u,i)
#     n    = length(x)
#     X    = [x; u]
#     Tr   = strip(model.T,1:size(X,1))
#     ass  = get_assignments_robust(model.M, X, Tr)
#     xnew = zeros(T,n)
#     dyn  = model.dynamics[ass[1]]
#     xnew = dyn.fx*x + dyn.fu*u
# end
#
# """
#     df(model::GMMModel,x,u)
# x ∈ ℜ n × N
# u ∈ ℜ m × N
# """
# function df(model::GMMModel,x,u)
#     u[isnan(u)] = 0
#     n,N = size(x)
#     m   = size(u,1)
#     fx  = Array{Float64}(undef,n,n,N)
#     fu  = Array{Float64}(undef,n,m,N)
#     X   = [x; u]
#     Tr  = strip(model.T,1:size(X,1))
#     ass = get_assignments_robust(model.M, X, Tr)
#     for i = 1:N
#         dyn = model.dynamics[ass[i]]
#         fx[:,:,i] = dyn.fx
#         fu[:,:,i] = dyn.fu
#     end
#     fxx=fxu=fuu = []
#     return fx,fu,fxx,fxu,fuu
# end

ext: jl


text:
# # Entries
#
# The key ingredient for data representations are `AxisEntries`.
#
# ## The `AxisEntries` type
#
# An `AxisEntries` object is
# made of four components:
# - axis,
# - entries,
# - scales,
# - labels.

using AlgebraOfGraphics, CairoMakie
using AlgebraOfGraphics: CategoricalScale, ContinuousScale
resolution = (600, 600)
fig = Figure(; resolution)
N = 11
rg = range(1, 2, length=N)
ae = AxisEntries(
    Axis(fig[1, 1]),
    [
        Entry(
            plottype=Scatter,
            positional=(rg, cosh.(rg)),
            named=(color=1:N, marker=fill("b", N));
            attributes=Dict(:markersize => 15)
        ),
        Entry(
            plottype=Scatter,
            positional=(rg, sinh.(rg)),
            named=(color=1:N, marker=fill("c", N));
            attributes=Dict(:markersize => 15)
        ),
    ],
    Dict(
        1 => ContinuousScale(identity, (0, 4)),
        2 => ContinuousScale(identity, (0, 4)),
        :color => ContinuousScale(identity, (1, N)),
        :marker => CategoricalScale(["a", "b", "c"], [:circle, :utriangle, :dtriangle]),
    ), # scales
    Dict(
        1 => "x",
        2 => "y",
        :color => "identity",
        :marker => "function"
    ), # labels
)
plot!(ae)
fig


text:
using Transpiler
using Base.Test
import Transpiler: gli, GLMethod
using Sugar: getsource!, dependencies!
# empty caches
Transpiler.empty_caches!()

cl_mapkernel = GLMethod((test, (Float32, Float32)))
source = getsource!(cl_mapkernel)
testsource = """float test(float a, float b)
{
    float y;
    float x;
    x = sqrt(sin(a) * b) / float(10.0);
    y = float(33.0) * x + cos(b);
    return y * float(10.0);
}"""
@testset "test kernel" begin
    @test source == testsource
end

decl = GLMethod((fortest, (Float32,)))
source = Sugar.getsource!(decl)
#TODO remove xxtempx4, which is unused now...

testsource = """float fortest(float x)
{
    float acc;
    int xxtempx4;
    int i;
    acc = x;
    for(i = 1; i <= 5; i++){
        if(i == 1){
            acc = acc + x;
            continue;
        };
        if(i == 2){
            acc = acc - x;
            continue;
        };
        acc = acc + x * x;
    };
    return acc;
}"""
@testset "for loops + if elseif " begin
    @test testsource == source
end

using Colors, Transpiler, StaticArrays
const Vec = SVector
const Vec3f0 = Vec{3, Float32}

immutable Light{T}
    position::Vec{3,T}
    ambient::Vec{3,T}
    diffuse::Vec{3,T}
    diffuse_power::T
    specular::Vec{3,T}
    specular_power::T
end

immutable Shading{T}
    ambient::Vec{3, T}
    specular::Vec{3, T}
    shininess::T
end

function blinnphong{NV, T}(V::Vec{NV, T}, N, L, color, shading, light)
    lambertian = max(dot(L, N), 0f0)
    half_direction = normalize(L .+ V)
    specular_angle = max(dot(half_direction, N), 0.0)
    specular = specular_angle ^ shading.shininess
    surf_color = (lambertian * color) .+ (specular * shading.specular)
    return light.ambient .+ surf_color
end

decl = Transpiler.GLMethod((blinnphong, (
    Vec3f0, Vec3f0, Vec3f0, Vec3f0, Shading{Float32}, Light{Float32}
)))

source = Sugar.getsource!(decl);
testsource = """vec3 blinnphong(vec3 V, vec3 N, vec3 L, vec3 color, Shading_float shading, Light_float light)
{
    vec3 surf_color;
    float specular;
    float specular_angle;
    vec3 half_direction;
    float lambertian;
    lambertian = max(dot(L, N), 0.0);
    half_direction = normalize(L + V);
    specular_angle = max(dot(half_direction, N), 0.0);
    specular = pow(specular_angle, shading.shininess);
    surf_color = lambertian * color + specular * shading.specular;
    return light.ambient + surf_color;
}"""

@testset "blinnphong" begin
    @test source == testsource
end


text:
push!(LOAD_PATH, joinpath(pwd(), "src"))

using JuML
using Pkg

#fannie headers:
colnames = ["loan_id", "monthly_rpt_prd", "servicer_name", "last_rt", "last_upb", "loan_age",
    "months_to_legal_mat" , "adj_month_to_mat", "maturity_date", "msa", "delq_status",
    "mod_flag", "zero_bal_code", "zb_dte", "lpi_dte", "fcc_dte","disp_dt", "fcc_cost",
    "pp_cost", "ar_cost", "ie_cost", "tax_cost", "ns_procs", "ce_procs", "rmw_procs",
    "o_procs", "non_int_upb", "prin_forg_upb_fhfa", "repch_flag", "prin_forg_upb_oth",
    "transfer_flg"]

importcsv("E:\\FannieMae\\Performance_All\\Performance_2017Q4.txt"; sep = "|", colnames = colnames,
           isinteger = ((colname, freq) -> colname == "loan_id"),
           isdatetime = ((colname, freq) ->
                            if colname == "monthly_rpt_prd" || colname == "lpi_dte" || colname == "fcc_dte" || colname == "disp_dt"
                                true, "mm/dd/yyyy"
                            elseif colname == "maturity_date" || colname == "zb_dte"
                                true, "mm/yyyy"
                            else
                                false, ""
                            end
                        ));

fannie_df = DataFrame("D:\\FannieMae\\Performance_All\\Performance_2017Q4"; preload = true)

@time s = summary(fannie_df[:servicer_name])

@time s = summary(factor(fannie_df[:loan_age]))
@time f = factor(fannie_df[:loan_age])
@time summary(f)

max_stars_repo_name: Statfactory/JuML

max_issues_repo_name: Statfactory/JuML

max_forks_repo_name: Statfactory/JuML


text:
#=
BINARYBUILDER_RUNNER=privileged julia build_tarballs.jl --verbose --debug
=#

using BinaryBuilder
using Pkg

version = v"0.0.1"

sources = [
	DirectorySource("./iglwrap"),
	# 2.3.0 stable release:
	GitSource("https://github.com/libigl/libigl.git",
		"e60423e28c86b6aa2a3f6eb0112e8fd881f96777")
]

script = raw"""
cd $WORKSPACE/srcdir
echo a
mkdir build||true; cd build
cmake \
  -DCMAKE_INSTALL_PREFIX=$prefix \
	-DCMAKE_MODULE_PATH=${WORKSPACE}/srcdir/libigl/cmake\; \
	..
make && make install
install_license ${WORKSPACE}/srcdir/LICENSE
"""

products = [
	LibraryProduct("libiglwrap", :libiglwrap;
		dlopen_flags=[:RTLD_NOW,:RTLD_LOCAL])
]

platforms = [
	Platform("x86_64", "linux"),
]
# platforms = supported_platforms()

platforms = expand_cxxstring_abis(platforms)

dependencies = [
  Dependency("boost_jll", compat="~1.71"),
	Dependency("GMP_jll"),
	Dependency("MPFR_jll"),
	Dependency("OpenBLAS_jll"),
 ]

build_tarballs(ARGS, "IGLWrap", version,
    sources,
    script,
    platforms,
    products,
    dependencies,
)

max_stars_repo_name: plut/iglwrap

max_issues_repo_name: plut/iglwrap

max_forks_repo_name: plut/iglwrap


text:
export broadsum, broadsum!, bsum, bsum!

"""
	broadsum(*, A, B, C)       = sum(A .* B .* C)
	broadsum(f,*, A,B)         = sum(f, A .* B)
These aim to work exactly like `sum(broadcast(...))`, but without materialising the broadcast array.
Simplest case works & is fast. Version with `f` is slow.

	broadsum(*, A, B; dims=1)  = sum(A .* B; dims=1)
    broadsum!(Z, *,A,B)        = sum!(Z, A .* B)
Similarly immitates `sum!(Z, broadcast(...))` without materialising.
Now uses `LazyArrays.BroadcastArray`. The in-place form is actually a little slower.
"""
bsum(op, As...; dims=:) =
    _bsum(dims, identity, Broadcast.broadcasted(op, As...))

bsum(f::Function, op::Function, As...) =
    _bsum((:), f, Broadcast.broadcasted(op, As...))

@inline function _bsum(::Colon, fun::Function, bc)
    @assert length(bc.args) >= 1

    T = Broadcast.combine_eltypes(fun∘bc.f, bc.args)
    tot = zero(T)
    @simd for I in eachindex(bc)
        @inbounds tot += fun(bc[I])
    end
    tot
end

using LazyArrays: BroadcastArray

_bsum(dims::Union{Int,NTuple}, ::typeof(identity), bc) = sum(BroadcastArray(bc); dims=dims)

bsum!(Z::AbstractArray, op, As...) = sum!(Z, BroadcastArray(op, As...))

const broadsum = bsum
const broadsum! = bsum!

# also dispatch complete sum to method which has better reverse:
# _bsum(::Colon, ::typeof(identity), bc) = sum_(BroadcastArray(bc))
# or not, crashes?

# Base._sum(A::BroadcastArray, ::Colon) = _bsum((:), identity, A.broadcasted)


text:
using Pkg

function general_registry_location()
    general_registry_dir = joinpath(DEPOT_PATH[1], "registries", "General")
    registry_toml_file = joinpath(general_registry_dir, "Registry.toml")
    return general_registry_dir, registry_toml_file
end

function general_registry_exists()
    general_registry_dir, registry_toml_file = general_registry_location()
    if !isdir(general_registry_dir)
        return false
    elseif !isfile(registry_toml_file)
        return false
    else
        return true
    end
end

function add_general_registry()
    @info("Attempting to clone the General registry")
    general_registry_dir, registry_toml_file = general_registry_location()
    rm(general_registry_dir; force = true, recursive = true)
    Pkg.Registry.add("General")
    isfile(registry_toml_file) || throw(ErrorException("the Registry.toml file does not exist"))
    return nothing
end

function main(; n = 10, max_delay = 120)
    VERSION >= v"1.5-" || return

    if general_registry_exists()
        @info("The General registry already exists locally")
        return
    end

    delays = ExponentialBackOff(; n = n, max_delay = max_delay)
    try
        retry(add_general_registry; delays = delays)()
        @info("Successfully add the General registry")
    catch ex
        msg = "I was unable to added the General registry. However, the build will continue."
        @error(msg, exception=(ex,catch_backtrace()))
    end

    return
end

main()

max_stars_repo_name: anandijain/julia-runtest

max_issues_repo_name: anandijain/julia-runtest

max_forks_repo_name: anandijain/julia-runtest


text:
#=
  From
  "Probabilistic logic programming and its applications"
  Luc De Raedt, Leuven
  https://www.youtube.com/watch?v=3lnVBqxjC88
  @ 3:44

  """
  Mike has a bag of marbles with 4 white, 8 blue, and
  6 red marbles. He pulls out one marble from the bag
  and it is red. What is the probability that the
  second marble he pulls out of the bag is white?

  The answer is 0.234931.
  """

  Cf ~/cplint/bag_of_marbles.pl
     ~/blog/bag_of_marbles.blog
     ~/psi/bag_of_marbles.psi
     ~/webppl/bag_of_marbles.wppl

  Cf bag_of_marbles.jl. This current model use a for loop.

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

#=
draw1: The probability of drawing white is the probability for 1.0

Distributions of variable draw[2] (num:0)
2.00000 =>  188281  (0.470703)
3.00000 =>  117185  (0.292963)
1.00000 =>   94534  (0.236335)
=#

@model function bag_of_marbles()

    white = 1
    blue  = 2
    red   = 3

    start = [4,8,6]

    draw = tzeros(2)
    for t in 1:2
        # We have to use sum(Int8[...]) since sum([]) throws an error
        draw[t] ~ Categorical(
                     simplex([
                            start[white] - sum(Int8[draw[i]==white for i in 1:t-1]),
                            start[blue]  - sum(Int8[draw[i]==blue  for i in 1:t-1]),
                            start[red]   - sum(Int8[draw[i]==red   for i in 1:t-1])
                            ])
                            )
    end

    true ~ Dirac(draw[1] == red)

    return draw[2]==white

end

model = bag_of_marbles()
num_chains = 4

# HH has problem with this!
chains = sample(model, MH(), MCMCThreads(), 100_000, num_chains)
# chains = sample(model, MH(), 10_000)

# chains = sample(model, PG(15), MCMCThreads(), 1_000, num_chains)
# chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, SMC(1000), 10_000)

# Note: IS don't generate chains the same way as MH, PG, and SMC!
# chains = sample(model, IS(), MCMCThreads(), 1000, num_chains)

display(chains)

show_var_dist_pct(chains,Symbol("draw[1]"))
println("\ndraw1: The probability of drawing white is the probability for 1.0")
show_var_dist_pct(chains,Symbol("draw[2]"))

# println("\ndraw1==white:")
# genq = generated_quantities(model,chains)
# show_var_dist_pct(genq,20)
# println("quantile:")
# println(quantile(vec(genq), [0.0, 0.025, 0.25, 0.5, 0.75, 0.975, 1.0]))

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
density(f::PolyElem) = length(findall(x->!iszero(x), coefficients(f)))/length(f)

#move elsewhere? Not used in here
function Hecke.representation_matrix(a::ResElem{<:PolyElem})
  R = parent(a)
  S = base_ring(base_ring(R))
  n = degree(modulus(R))
  m = zero_matrix(S, n, n)
  for i=0:n-1
    m[1, i+1] = coeff(a.data, i)
  end
  for j=2:n
    a *= gen(R)
    for i=0:n-1
      m[j, i+1] = coeff(a.data, i)
    end
  end
  return m
end


text:
# Plot sparse matrices

using TrajOptSOCPs
using Plots, SparseArrays

function plotHessianAL(al::augLag, traj)
    h = evalHessAl(al, traj)
    plth = spy(sparse(h))
    return plth
end

max_stars_repo_name: danineamati/TrajOptSOCPs


text:
using LinearAlgebra, Statistics
using Distributions, Plots, Printf, Random
using QuantEcon

## Simulation
d = Categorical([0.5, 0.3, 0.1, 0.1]);
@show rand(d, 5);
@show supertype(typeof(d));
@show pdf(d, 4);
@show support(d);
@show pdf.(d, support(d));

function mc_sample_path(P; init = 1, sample_size = 1000)
    @assert size(P)[1] == size(P)[2] # square required
    N = size(P)[1]; # should be square

    # create vector of discrete RVs for each row of stochastic matrix P
    dists = [Categorical(P[i, :]) for i in 1:N]

    # set up Simulation
    X = fill(0, sample_size); # allocate memory, or zeros(Int64, sample_size)
    X[1] = init # set the initial state

    for t in 2:sample_size
        dist = dists[X[t - 1]];
        X[t] = rand(dist);
    end

    return X
end

P = [0.4 0.6; 0.2 0.8];

X = mc_sample_path(P, sample_size = 100_000); # note 100_000 = 100000
μ_1 = count(X .== 1) / length(X);  # .== broadcasts test for equality. Could use mean(X .== 1)


text:
export fgh_FEOL_2, FixedEffectsOrderedLogit_2, DiD!

"""
    FixedEffectsOrderedLogit_2(formula,data,isymbol,tsymbol; b0, solver)

computes the composite conditional MLE for the FEOL model with two time periods.

The formula, data, isymbol, tsymbol tell the algorithm which model and data to use.
If you want to supply a starting value, assign `b0` (or it will start from zeros).
The default solver is `LBFGS` (using `Optim.jl`). Set `solver=Newton()` to use NewtonRaphson
(analytical Hessian will be used), or set it to...

- `NewtonTrustRegion`
- `AcceleratedGradientDescent`
- `NelderMead`
- ... (see the documentation for `Optim.jl`),

but in that case you must issue `using Optim` where you are issuing the estimation command.
"""
function FixedEffectsOrderedLogit_2(y,X,cuts1,cuts2;b0=false,solver = LBFGS(), verbose = true)

    # set starting value.
    if b0 == false
        n,K = size(X)
        # Count how many additional columns there will be.
        n_1 = length(cuts1) - 1
        n_2 = length(cuts2)
        b_start = zeros(K+n_1+n_2)

        println()
    else
        b_start = b0
    end

    # Define the objective, gradient, and Hessian
    objective, gradient!, Hessian! = fgh_FEOL_2(b_start,y,X,cuts1,cuts2)

    # Perform the optimization
    println("========================================")
    println("Starting the first minimization attempt.")
    println("========================================")
    opt_res = optimize(objective, gradient!, Hessian!, b_start, solver)

    # Report optimization diagnostics to user.
    if verbose
        print(opt_res)
    end

    if opt_res.g_converged
        println("")
        println("===============================")
        println("Converged on the first attempt.")
        println("===============================")
        println("")
    else
        opt_res_1 = opt_res
        println("")
        println("")
        println("====================================================")
        println("Starting a second round, using Newton Trust Region..")
        println("====================================================")
        println("")
        opt_res = optimize(objective, gradient!, b_start, NewtonTrustRegion())
        if verbose
            print(opt_res)
        end
        println("")
        println("")
        println("======================================================")
        println("Further decrease of $(opt_res_1.minimum - opt_res.minimum).")
        println("======================================================")
    end

    return opt_res.minimizer
end

"""
    fgh_FEOL_2(theta,y,X)

computes the (negative of the) log-likelihood, score, and
Hessian for the fixed effects ordered logit model.
"""
function fgh_FEOL_2(theta::Array,y::Array,X::Array, cuts1, cuts2)

    # The data (y,X) comes in cross-sectionalized: n high, 2 cq K wide.
    # Settle some dimensions first.
    n, K = size(X)
    n_1 = length(cuts1) - 1
    n_2 = length(cuts2)

    #
    # Probably should do some checks on the number of observations for each set of cut points.
    #

    # theta is a vector of coefficients. First K are beta.
    function min_ll(theta::Array{Float64})

#         println(cuts1)
#         println(cuts2)
#         println(n_1)
#         println(n_2)
#         println(K)

        # Initialize
        Xb = X*theta[1:K]
        ll = 0.

        for i in 1:(n_1+1)

            d1 = y[:,1].>=cuts1[i]

            for j in 1:n_2

                d2 = y[:,2].>=cuts2[j]
                dbar = (d1.+d2).==1

                if i == 1
                    h = logistic.(Xb .- theta[K+n_1+j])
                else
                    h = logistic.(Xb .+ theta[K+i-1] .- theta[K+n_1+j])
                end

                ll = ll - mean(dbar.*(d2.*log.(h) .+ (1.-d2).*log.(1.-h)))

            end
        end

        return ll

    end

    function min_score!(J::Array,theta::Array)

       # Initialize the Jacobian to zero.
        J[:] .= 0.
        # Set the linear index.
        Xb = X*theta[1:K]

        for i in 1:length(cuts1)

            d1 = y[:,1].>=cuts1[i]

            for j in 1:length(cuts2)

                Jnew = zeros(K+n_1+n_2)

                d2 = y[:,2].>=cuts2[j]
                dbar = (d1 .+ d2) .== 1

                if i == 1
                    scalar_in = dbar.*(d2 .- logistic.(Xb .- theta[K+n_1+j]))
                else
                    scalar_in = dbar.*(d2 .- logistic.(Xb .+ theta[K+i-1] .- theta[K+n_1+j]))
                end

                Jnew[1:K] = X'*scalar_in/n

                cut_H = mean(scalar_in)
                if i>1
                    Jnew[K+i-1] = cut_H
                end
                Jnew[K+n_1+j] = -cut_H

                # Assemble it.
                J[:] = J - Jnew
            end
        end
    end

    function min_Hess!(H::Array,theta::AbstractArray)

        # Initialize the Jacobian to zero.
        H[:] .= 0.

        # Predefine the linear index.
        Xb = X*theta[1:K]

        for i in 1:length(cuts1)

            d1 = y[:,1].>=cuts1[i]

            for j in 1:length(cuts2)

                # Initialize the Hessian for this case
                Hnew = zeros(K+n_1+n_2,K+n_1+n_2)

                d2 = y[:,2].>=cuts2[j]
                dbar = (d1 + d2) .== 1

                # Compute the top-left block.
                if i == 1
                    Lambda = dbar.*(logistic.(Xb - theta[K+n_1+j]))
                else
                    Lambda = dbar.*(logistic.(Xb + theta[K+i-1] - theta[K+n_1+j]))
                end
                Hnew[1:K,1:K] = -(1/n)*(Lambda.*(1-Lambda) .* X)'*X

                # Compute the bottom-right block.
                cut_H = -mean(Lambda.*(1-Lambda))
                if i>1
                    Hnew[K+i-1,K+i-1] = cut_H
                end
                Hnew[K+n_1+j,K+n_1+j] = cut_H

                # Compute the off-diagonal blocks.
                H_off = [ -mean( Lambda.*(1-Lambda) .* X[:,k]) for k in 1:K ]

                if i>1
                    Hnew[1:K,K+i-1] = H_off
                    Hnew[K+i-1,1:K] = H_off'
                end
                Hnew[1:K,K+n_1+j] = -H_off
                Hnew[K+n_1+j,1:K] = -H_off'

                H[:] = H - Hnew

            end
        end
    end

    return min_ll, min_score!, min_Hess!
end

"""
    DID!(f)

computes the linear DiD estimate using a simple comparison of means.
Should be updated to allow for covariates.
"""
function DiD!(f::NLDID)
    dif_treat = mean(f.y_treat[f.t.==1]) - mean(f.y_treat[f.t.==0])
    dif_control = mean(f.felt.y[:,2]) - mean(f.felt.y[:,1])
    f.linearATT = dif_treat - dif_control
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: healthlake
using AWS.Compat
using AWS.UUIDs

"""
    create_fhirdatastore(datastore_type_version)
    create_fhirdatastore(datastore_type_version, params::Dict{String,<:Any})

Creates a Data Store that can ingest and export FHIR formatted data.

# Arguments
- `datastore_type_version`: The FHIR version of the Data Store. The only supported version
  is R4.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"ClientToken"`: Optional user provided token used for ensuring idempotency.
- `"DatastoreName"`: The user generated name for the Data Store.
- `"PreloadDataConfig"`: Optional parameter to preload data upon creation of the Data
  Store. Currently, the only supported preloaded data is synthetic data generated from
  Synthea.
"""
create_fhirdatastore(DatastoreTypeVersion; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("CreateFHIRDatastore", Dict{String, Any}("DatastoreTypeVersion"=>DatastoreTypeVersion, "ClientToken"=>string(uuid4())); aws_config=aws_config)
create_fhirdatastore(DatastoreTypeVersion, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("CreateFHIRDatastore", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("DatastoreTypeVersion"=>DatastoreTypeVersion, "ClientToken"=>string(uuid4())), params)); aws_config=aws_config)

"""
    delete_fhirdatastore()
    delete_fhirdatastore(params::Dict{String,<:Any})

Deletes a Data Store.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"DatastoreId"`:  The AWS-generated ID for the Data Store to be deleted.
"""
delete_fhirdatastore(; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("DeleteFHIRDatastore"; aws_config=aws_config)
delete_fhirdatastore(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("DeleteFHIRDatastore", params; aws_config=aws_config)

"""
    describe_fhirdatastore()
    describe_fhirdatastore(params::Dict{String,<:Any})

Gets the properties associated with the FHIR Data Store, including the Data Store ID, Data
Store ARN, Data Store name, Data Store status, created at, Data Store type version, and
Data Store endpoint.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"DatastoreId"`: The AWS-generated Data Store id. This is part of the
  ‘CreateFHIRDatastore’ output.
"""
describe_fhirdatastore(; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("DescribeFHIRDatastore"; aws_config=aws_config)
describe_fhirdatastore(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("DescribeFHIRDatastore", params; aws_config=aws_config)

"""
    describe_fhirexport_job(datastore_id, job_id)
    describe_fhirexport_job(datastore_id, job_id, params::Dict{String,<:Any})

Displays the properties of a FHIR export job, including the ID, ARN, name, and the status
of the job.

# Arguments
- `datastore_id`: The AWS generated ID for the Data Store from which files are being
  exported from for an export job.
- `job_id`: The AWS generated ID for an export job.

"""
describe_fhirexport_job(DatastoreId, JobId; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("DescribeFHIRExportJob", Dict{String, Any}("DatastoreId"=>DatastoreId, "JobId"=>JobId); aws_config=aws_config)
describe_fhirexport_job(DatastoreId, JobId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("DescribeFHIRExportJob", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("DatastoreId"=>DatastoreId, "JobId"=>JobId), params)); aws_config=aws_config)

"""
    describe_fhirimport_job(datastore_id, job_id)
    describe_fhirimport_job(datastore_id, job_id, params::Dict{String,<:Any})

Displays the properties of a FHIR import job, including the ID, ARN, name, and the status
of the job.

# Arguments
- `datastore_id`: The AWS-generated ID of the Data Store.
- `job_id`: The AWS-generated job ID.

"""
describe_fhirimport_job(DatastoreId, JobId; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("DescribeFHIRImportJob", Dict{String, Any}("DatastoreId"=>DatastoreId, "JobId"=>JobId); aws_config=aws_config)
describe_fhirimport_job(DatastoreId, JobId, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("DescribeFHIRImportJob", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("DatastoreId"=>DatastoreId, "JobId"=>JobId), params)); aws_config=aws_config)

"""
    list_fhirdatastores()
    list_fhirdatastores(params::Dict{String,<:Any})

Lists all FHIR Data Stores that are in the user’s account, regardless of Data Store
status.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"Filter"`: Lists all filters associated with a FHIR Data Store request.
- `"MaxResults"`: The maximum number of Data Stores returned in a single page of a
  ListFHIRDatastoresRequest call.
- `"NextToken"`: Fetches the next page of Data Stores when results are paginated.
"""
list_fhirdatastores(; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("ListFHIRDatastores"; aws_config=aws_config)
list_fhirdatastores(params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("ListFHIRDatastores", params; aws_config=aws_config)

"""
    start_fhirexport_job(client_token, data_access_role_arn, datastore_id, output_data_config)
    start_fhirexport_job(client_token, data_access_role_arn, datastore_id, output_data_config, params::Dict{String,<:Any})

Begins a FHIR export job.

# Arguments
- `client_token`: An optional user provided token used for ensuring idempotency.
- `data_access_role_arn`: The Amazon Resource Name used during the initiation of the job.
- `datastore_id`: The AWS generated ID for the Data Store from which files are being
  exported for an export job.
- `output_data_config`: The output data configuration that was supplied when the export job
  was created.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"JobName"`: The user generated name for an export job.
"""
start_fhirexport_job(ClientToken, DataAccessRoleArn, DatastoreId, OutputDataConfig; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("StartFHIRExportJob", Dict{String, Any}("ClientToken"=>ClientToken, "DataAccessRoleArn"=>DataAccessRoleArn, "DatastoreId"=>DatastoreId, "OutputDataConfig"=>OutputDataConfig); aws_config=aws_config)
start_fhirexport_job(ClientToken, DataAccessRoleArn, DatastoreId, OutputDataConfig, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("StartFHIRExportJob", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClientToken"=>ClientToken, "DataAccessRoleArn"=>DataAccessRoleArn, "DatastoreId"=>DatastoreId, "OutputDataConfig"=>OutputDataConfig), params)); aws_config=aws_config)

"""
    start_fhirimport_job(client_token, data_access_role_arn, datastore_id, input_data_config)
    start_fhirimport_job(client_token, data_access_role_arn, datastore_id, input_data_config, params::Dict{String,<:Any})

Begins a FHIR Import job.

# Arguments
- `client_token`: Optional user provided token used for ensuring idempotency.
- `data_access_role_arn`: The Amazon Resource Name (ARN) that gives Amazon HealthLake
  access permission.
- `datastore_id`: The AWS-generated Data Store ID.
- `input_data_config`: The input properties of the FHIR Import job in the StartFHIRImport
  job request.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"JobName"`: The name of the FHIR Import job in the StartFHIRImport job request.
"""
start_fhirimport_job(ClientToken, DataAccessRoleArn, DatastoreId, InputDataConfig; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("StartFHIRImportJob", Dict{String, Any}("ClientToken"=>ClientToken, "DataAccessRoleArn"=>DataAccessRoleArn, "DatastoreId"=>DatastoreId, "InputDataConfig"=>InputDataConfig); aws_config=aws_config)
start_fhirimport_job(ClientToken, DataAccessRoleArn, DatastoreId, InputDataConfig, params::AbstractDict{String}; aws_config::AbstractAWSConfig=global_aws_config()) = healthlake("StartFHIRImportJob", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClientToken"=>ClientToken, "DataAccessRoleArn"=>DataAccessRoleArn, "DatastoreId"=>DatastoreId, "InputDataConfig"=>InputDataConfig), params)); aws_config=aws_config)


text:
using DrWatson
@quickactivate "StatReth"

# %%
using Distributions
using Plots
using StatsPlots

# %% 4.1
pos = [sum(rand(Uniform(-1, 1), 16)) for _ in 1:1000]

# %% 4.2
prod(1 .+ rand(Uniform(0.0, 0.1), 12))

# %% 4.3
growth = [prod(1 .+ rand(Uniform(0.0, 0.1), 12)) for _ in 1:10_000]
density(growth)

# %% 4.4
big = [prod(1 .+ rand(Uniform(0.0, 0.5), 12)) for _ in 1:10_000]
small = [prod(1 .+ rand(Uniform(0.0, 0.01), 12)) for _ in 1:10_000]

# %% 4.5
log_big = [prod(1 .+ log.(rand(Uniform(0.0, 0.5), 12))) for _ in 1:10_000]
density(log_big)
density(log_big, xscale = :log10)


text:
function chisq_goodness(x::AbstractVector{T}) where {T <: Real}
    # figure out number of bins, which depends on length(x)
    nbins = round(Int, ceil(2 * (length(x)^(2/5))))

    # create a `range` object
    bin_range = range(-3, 3, length = nbins + 1)

    # use StatsBase `fit` function to count number of observations per bin
    bincounts = fit(Histogram, x, bin_range).weights

    # use `pdf` function to find expected number of observations per bin
    bin_centers = bin_range[1:end-1] .+ bin_range.step.hi / 2
    expected_pdf = pdf.(Normal(0, 1), bin_centers)
    bin_expected = floor.(Int, expected_pdf .* sum(bincounts) / sum(expected_pdf))

    # calcuate pearson stats (X2)
    𝜒² = pearson_chi(bincounts, bin_expected)

    # calculate X2/DF and return it
    return 𝜒² / nbins
end

function pearson_chi(
    observed::AbstractVector{T},
    expected::AbstractVector{T}
) where {T <: Integer}
    χ² = 0
    for i = 1:length(observed)
        χ² += (observed[i] - expected[i])^2 / expected[i]
    end
    return χ²
end


text:
#---------------------------------------------------------------------#
#This code computes the Legendre cardinal functions at the Quadrature
#points
#Written by F.X. Giraldo on 4/2000
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey; CA 93943-5216
#---------------------------------------------------------------------#
function lagrange_basis(P,Q,xgl,xnq)

    #Initialize arrays
    psi=zeros(P,Q)
    dpsi=zeros(P,Q)

    #Perform Quadrature/Sample the Function
    for l=1:Q
        xl=xnq[l]

        #Construct Basis
        for i=1:P
            xi=xgl[i]
            psi[i,l]=1
            dpsi[i,l]=0
            for j=1:P
                xj=xgl[j]
                #Basis
                if (j != i)
                    psi[i,l]=psi[i,l]*(xl-xj)/(xi-xj)
                end
                ddpsi=1
                if (j != i)
                    for k=1:P
                        xk=xgl[k]
                        #Derivative of Basis
                        if (k !=i && k !=j)
                            ddpsi=ddpsi*(xl-xk)/(xi-xk)
                        end
                    end
                    dpsi[i,l]=dpsi[i,l] + ddpsi/(xi-xj)
                end
            end
        end
    end

    return (psi,dpsi)
end


text:
function parse_command(command_text::String)::Tuple{String, Dict{String, String}}
    command_name = try string(match(r"/([^@\s]+)", command_text)[1]) catch; "" end
    params = try match(r"\s(.*)$", command_text)[1] |> String |> split catch; "" end

    keys = map(
        p -> replace(p, r"[-—]" => "") |> lowercase,
        params[1:2:end]
    )
    values = params[2:2:end]

    command_params = Dict{String, String}(zip(keys, values))

    return command_name, command_params
end

function execute_command(command_name::String, command_params::Dict)::String
    try
        return HTTP.request("POST", API_URL * command_name; query = command_params).body |>
            String |>
            JSON.parse |>
            response -> response["message"]
    catch
        return "Something is wrong with the params or the API is unavailable, master."
    end
end

max_stars_repo_name: iwsylit/spotify-utils-bot

max_issues_repo_name: iwsylit/spotify-utils-bot

max_forks_repo_name: iwsylit/spotify-utils-bot


text:
@testset "Test interpolators" begin

    using jFEMTools
    import Tensors
    jF = jFEMTools

    # Monomials
    ###############################################
    P1_basis = jF.Monomials(2,1,Tensors.Vec{2}((0.0,0.0)),1.0)
    @test jF.value(P1_basis,1,Tensors.Vec{2}((1.0,2.0))) == 1.0
    @test jF.value(P1_basis,2,Tensors.Vec{2}((1.0,2.0))) == 1.0
    @test jF.value(P1_basis,3,Tensors.Vec{2}((1.0,2.0))) == 2.0
    @test jF.getnbasefunctions(P1_basis) == 3

end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools


text:
module MinimumVolumeEllipsoids

using Distributions
using LinearAlgebra
using PDMats
using SparseArrays
using SpecialFunctions

export Ellipsoid

export minimum_volume_ellipsoid
export rand
export volume

struct Ellipsoid
    H::AbstractPDMat
    x::AbstractVector
end

function minimum_volume_ellipsoid(
    X::AbstractMatrix,
    tol::Real=1e-7,
    KKY::Integer=0,
    maxit::Integer=100000;
    centered::Bool=false,
)
    n, m = size(X)

    if centered
        _, R, ϕ = _minvol(X, tol, KKY, maxit)

        H = inv(R.L * R.L' / ϕ)
        x = zeros(n)
    else
        Y = [X; ones(1, m)]

        u, R, ϕ = _minvol(Y, tol, KKY, maxit)

        H = inv(R.L * R.L' / ϕ)[1:n, 1:n]
        x = X * u
    end

    H = (H + H') / 2
    H = PDMat(H)

    return Ellipsoid(H, x)
end

function _minvol(X::AbstractMatrix, tol::Real=1e-7, KKY::Integer=0, maxit::Integer=100000)
    n, m = size(X)

    n100 = max(n, 100)
    n50000 = max(n, 50000)

    iter = 0

    if KKY >= 1
        u = (1 / m) * ones(m)
    else
        u = initwt(X)
    end

    # Initialize Cholesky Factor
    upos = findall(u .> 0)
    R, var = _compute_R_and_var(u, X, upos)
    ϕ = 1

    ω₊ = maximum(var)

    # Use the Harman-Pronzato test to see if columns of X can be eliminated.
    essential_indices = _harman_pronzato_elimination(ω₊, n, var, u)
    XX = X[:, essential_indices]

    # If only n columns remain, recompute u and R
    if length(essential_indices) == n
        u = (1 / n) * ones(n)
        upos = findall(u .> 1e-8)
        R, var = _compute_R_and_var(u, XX, upos)
    else
        var = var[essential_indices]
        u = u[essential_indices] / sum(u[essential_indices])
        upos = findall(u .> 1e-8)
    end

    # Find "furthest" and "closest" points
    ω₊, j = findmax(var)
    ω₋, i = findmin(var[upos])
    i = upos[i]

    while ((ω₊ > (1 + tol) * n) || (ω₋ < (1 - tol) * n)) && iter < maxit
        iter += 1

        (j, ω) = ω₊ + ω₋ > 2 * n ? (j, ω₊) : (i, ω₋)

        # compute new λ
        λ = (ω - n) / ((n - 1) * ω)
        λ = max(λ, -u[j])

        # Update u and make sure it stays nonnegative
        uold = u
        u[j] = max(u[j] + λ, 0)
        u = (1 / (1 + λ)) * u
        upos = findall(u .> 1e-8)

        Rxj = R.U' \ XX[:, j]
        x = ϕ * (R.U \ Rxj)
        ωⱼ = ϕ * (Rxj' * Rxj)

        recompute_R = abs(ωⱼ - ω) / max(1, ω) > 1e-8

        ω = ωⱼ

        # Every 50.000 iterations: Check if R should be recomputed
        if mod(iter, n50000) == 0
            upos = findall(uold .> 0)
            M = XX[:, upos] * Diagonal(uold[upos]) * XX[:, upos]'
            if norm(ϕ * M - R.L * R.L') / (ϕ * norm(M)) > 1e-8
                recompute_R = true
            end
        end

        if recompute_R
            upos = findall(u .> 0)
            R, var = _compute_R_and_var(u, XX, upos)
            ϕ = 1
        else
            xx = sqrt(abs(λ) * ϕ) * XX[:, j]
            if λ > 0
                lowrankupdate!(R, xx)
            else
                lowrankdowndate!(R, xx)
            end
            ϕ *= (1 + λ)
            mult = λ / (1 + λ * ω)
            var = (1 + λ) * (var - mult * transpose((transpose(x) * XX)) .^ 2)
        end

        ω₊, j = findmax(var)

        # Every 100 iterations: Check if more points can be eliminated
        if mod(iter, n100) == 0
            essential = _harman_pronzato_elimination(ω₊, n, var, u)
            if length(essential) < length(essential_indices)
                essential_indices = essential_indices[essential]
                XX = XX[:, essential]
                if length(essential_indices) == n
                    u = (1 / n) * ones(n, 1)
                    R, var = _compute_R_and_var(u, XX, upos)
                    ϕ = 1
                else
                    var = var[essential]
                    u = u[essential] / sum(u[essential])
                end
                ω₊, j = findmax(var)
            end
        end

        upos = findall(u .> 0)
        ω₋, i = findmin(var[upos])
        i = upos[i]
    end

    uu = zeros(m)
    uu[essential_indices] = u
    u = uu

    return u, R, ϕ
end

function _harman_pronzato_elimination(
    ω₊::Float64, n::Integer, var::AbstractVector, u::AbstractVector
)
    δn = ω₊ - n
    # TODO: Figure out when and why this happens
    if δn < 0
        δn = 0
    end
    threshold = n * (1 + δn / 2 - √(δn - δn / n + ((δn / n)^2 * n^2) / 4))
    return findall((var .> threshold) .| (u .> 1e-8))
end

function _compute_R_and_var(u::AbstractVector, X::AbstractMatrix, upos::AbstractVector)
    A = spdiagm(sqrt.(u[upos])) * X[:, upos]'
    F = qr(A)
    R = Cholesky(F.R, :U, 0)
    RX = R.U' \ X
    var = vec(sum(RX .* RX; dims=1))
    return R, var
end

"""
    initwt(X)

Obtain the initial weights `u` using the Kumar-Yildirim algorithm, taking into account that
`X` represents [X, -X].
"""
function initwt(X::AbstractMatrix)
    n, m = size(X)
    u = zeros(m)
    Q = 1.0 * I(n)
    d = Q[:, 1]

    for j in 1:n
        # compute the maximizer of | d'*x | over the columns of X.
        dX = vec(abs.(d' * X))
        _, ind = findmax(dX)
        u[ind] = 1

        j == n && break

        # update Q
        y = X[:, ind]
        z = Q' * y

        if j > 1
            z[1:(j - 1)] = zeros(j - 1, 1)
        end

        ζ = norm(z)
        zj = z[j]
        if zj < 0
            ζ = -ζ
        end
        zj += ζ
        z[j] = zj
        Q = Q - (Q * z) * ((1 / (ζ * zj)) * z')
        d = Q[:, j + 1]
    end

    u /= n
    return u
end

function Base.rand(ϵ::Ellipsoid, m::Integer)
    n = size(ϵ.H, 1)
    X = randn(n, m)
    X = X ./ kron(ones(n, 1), sqrt.(sum(X .^ 2; dims=1)))
    R = ones(n, 1) * rand(1, m) .^ (1 / n)
    sphere = R .* X
    ellipsoid = sqrt(n) * cholesky(inv(ϵ.H)).L * sphere + ϵ.x .* ones(1, m)
    return ellipsoid
end

include("volume.jl")

end # module


text:
Pkg.activate(".")
import SINDYPI_Julia

using SINDYPI_Julia

using Test

#Test 1: Glycolysis
using MAT
vars = matread("C:/Users/Raj/Desktop/SINDY PI Brunton/SINDy-PI-master/SINDy-PI-master/Comparison/DataLength/YeastGlycolysis/SINDy_PI/Datas/TrainingData.mat")
X = vars["xt"]
DX = vars["dxt"]

Polyorder =  Int64[6 6 2 2 2 6 2]
Trigorder = Int64[0 0 0 0 0 0 0]

SINDYPI(X, DX, Polyorder, Trigorder)

#Test 2: Pendulum on Cart
using DifferentialEquations, Plots
# Setup
initial = [0.3, 0, 1.0, 0]
tspan = (0,16.0)
dt = 0.001;

 m = 1
 M = 1
 L = 1
 g = 9.81

# Define the function
function Cart_Pendulum!(du, u,p,t)
    du[1] = u[3]
    du[2] = u[4]
    du[3] = -((M +m)*g*sin(u[1]) +m*(L^2)*sin(u[1])*cos(u[1])*u[3]*u[3] )/(L*L*(M + m -m*cos(u[1])*cos(u[1])))
    du[4] = (m*L*L*sin(u[1])*u[3]*u[3] + m*g*sin(u[1])*cos(u[1]))/(L*(M + m -m*cos(u[1])*cos(u[1])))
    return [du[1] du[2] du[3] du[4]]
end

#Pass to solvers
prob = ODEProblem(Cart_Pendulum!, initial, tspan)

sol = solve(prob, Tsit5(), saveat = dt)

plot(sol,linewidth=2,
     title ="Pendulum on Cart Problem",
     xaxis = "Time",
     yaxis = "Height",
     label = ["Theta","S", "dTheta", "dS"])

X = Array(sol)
DX = similar(X)
du = zeros(4,1)
for (i, xi) in enumerate(eachcol(X))
    DX[:,i] = Cart_Pendulum!(du, xi, [], 0.0)
end
X = X'
DX = DX'

Polyorder =  Int64[1 1 0 0]
Trigorder = Int64[0 0 3 4]

SINDYPI(X, DX, Polyorder, Trigorder)

#TEST 3: Michaelis - Menten
jx = 0.6
Vm = 1.5
Km = 0.3

f(u,p,t) = jx - (Vm*u/(Km + u))

u0 = 1/2
tspan = (0.0,3.0)
dt = 0.01
prob = ODEProblem(f,u0,tspan)

sol = DifferentialEquations.solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8, saveat = dt)

X = sol.u
DX = jx .- (Vm*X./(Km .+ X))

Polyorder =  Int64[1]
Trigorder = Int64[0]

SINDYPI(X, DX, Polyorder, Trigorder)

@testset "SINDYPI_Julia.jl" begin
    # Write your own tests here.
end


text:
module Nbnhhsh

using HTTP, JSON3

const POSTURL = "https://lab.magiconch.com/api/nbnhhsh/guess"
const BASEHEADER = Dict(
    "User-Agent" => "Mozilla/5.0 (X11; Linux x86_64; rv:89.0) Gecko/20100101 Firefox/89.0",
    "content-type" => "application/json",
    "Origin" => "https://lab.magiconch.com",
)

function hhsh(sx::AbstractString)
    post_data = """ {"text":"$(HTTP.escapeuri(sx))"} """
    return only(JSON3.read(HTTP.post(POSTURL, BASEHEADER, post_data).body))
end

end


text:
println("Loading packages...")
using PyPlot
using DataFrames
using Dates
using QuandlAccess
using Statistics
using MarketData
println("Done.\n")

# REFERENCES
# http://www.financialsense.com/contributors/cris-sheridan/financial-stress-measures

# FOR MORE COMPLETE US STATISTICS
# https://www.quandl.com/c/usa
# http://www.multpl.com/sitemap

#########################################
# number of years to obtain data
macro_yrs_ultralong = 20
future_yrs_long = 8
future_yrs_short = 1
date_plotend = today()
#########################################

MonthNameToNumber = Dict("JAN"=>1, "FEB"=>2, "MAR"=>3, "APR"=>4, "MAY"=>5, "JUN"=>6, "JUL"=>7, "AUG"=>8, "SEP"=>9, "OCT"=>10, "NOV"=>11, "DEC"=>12)

legend_fontsize = 10
maxattempts = 2
dpi = 109
figsize = (1920/dpi, 1080/dpi)
fmt = "yyyy-mm-dd" # date format string
isverbose = true
nrows_verbose = 5

function get_daily_data_from_fred(series_name::String, date_start::Date, date_end::Date, isverbose::Bool=true, nrows_verbose::Int=5)::DataFrame
    data = DataFrame(MarketData.fred(series_name))
    if isverbose
        println(data[1:nrows_verbose, :], "\n...")
    end
    DataFrames.rename!(data, [:date, :valuestr])
    data = data[date_start .<= data[!, :date] .<= date_end, :]
    data[!, :value] = zeros(Float64, size(data, 1))
    if typeof(data[1, :valuestr]) <: String
        for n = 1:size(data, 1)
            if data[n, :valuestr] == "."
                data[n, :value] = NaN
            else
                data[n, :value] = parse.(Float64, data[n, :valuestr])
            end
        end
    else
        data[!, :value] = data[!, :valuestr]
    end
    data = data[!, ["date", "value"]]
    data = data[completecases(data), :]
    data = filter(x->!isnan(x.value), data)
    sort!(data, [:date])
    return data
end

function get_daily_data_from_quandl(quandl::Quandl, series_name::String, date_start::Date, date_end::Date, cols::Vector{String}, isverbose::Bool=true, nrows_verbose::Int=5)
    data = quandl(QuandlAccess.TimeSeries(series_name); start_date = date_start, end_date = date_end)
    if isverbose
        println(data[1:nrows_verbose, :], "\n...")
    end
    data = data[!, cols]
    data = data[completecases(data), :]
    if length(cols) == 2
        DataFrames.rename!(data, [:date, :value])
        sort!(data, [:date])
        return data
    else
        datas = []
        for n = 1:(length(cols)-1)
            datatmp = data[!, [cols[1], cols[1+n]]]
            DataFrames.rename!(datatmp, [:date, :value])
            sort!(datatmp, [:date])
            push!(datas, datatmp)
        end
        return Tuple(datas)
    end
end

function get_daily_data_from_yahoo(symbol::String, date_start::Date, date_end::Date, isverbose::Bool=true, nrows_verbose::Int=5)::DataFrame
    data = DataFrame(MarketData.yahoo(symbol, MarketData.YahooOpt(period1 = DateTime(date_start), period2 = DateTime(date_end))))
    if isverbose
        println(data[1:nrows_verbose, :], "\n...")
    end
    data = data[!, ["timestamp", "Close"]]
    DataFrames.rename!(data, [:date, :value])
    data = data[completecases(data), :]
    data = filter(x->!isnan(x.value), data)
    sort!(data, [:date])
    return data
end

function calc_two_dataframes(data1::DataFrame, operator::String, data2::DataFrame)::DataFrame
    data = innerjoin(data1, data2, on=:date, makeunique=true)
    data = data[completecases(data), :]
    if operator == "/"
        data[!, :value] = data[!, :value] ./ data[!, :value_1]
    elseif operator == "*"
        data[!, :value] = data[!, :value] .* data[!, :value_1]
    elseif operator == "+"
        data[!, :value] = data[!, :value] .+ data[!, :value_1]
    elseif operator == "-"
        data[!, :value] = data[!, :value] .- data[!, :value_1]
    else
        error("Unknown dataframe operator $operator")
    end
    data = data[!, ["date", "value"]]

    return data
end

isdownloaded = false
quandl = Quandl(ENV["QUANDL_API_KEY"])

date_plotstart = date_plotend - Day(Int(round(macro_yrs_ultralong*365.25, digits=0)))
date_future_plotstart_long = date_plotend -  Day(Int(round(future_yrs_long*365.25, digits=0)))
date_future_plotstart_short = date_plotend - Day(Int(round(future_yrs_short*365.25, digits=0)))

cities_of_interest = ["National", "Chicago", "SanFrancisco", "LosAngeles", "SanDiego", "Portland", "Seattle", "Phoenix", "Dallas"]
colors = ["black", "green", "blue", "cyan", "magenta", "red", "yellow", "darkblue", "pink", "purple", "orange",
            "brown"]

close("all")

if !isdownloaded
    println("Downloading data from Fred, Yahoo, and Quandl...\n")

    # Permalink of any Quandl data set: https://www.quandl.com/data/xxx/yyy
    # where xxx/yyy is the Quandl Code

    # S&P500 index
    println("************** S&P500 (Fred) **************")
    SP500 = get_daily_data_from_yahoo("^GSPC", date_plotstart, date_plotend, isverbose, nrows_verbose)

    # Gold price
    println("\n************** Gold price (Quandl) **************")
    gold = get_daily_data_from_quandl(quandl, "LBMA/GOLD", date_plotstart, date_plotend, ["Date", "USD (PM)"], isverbose, nrows_verbose)

    # S&P500 / Gold
    SP500_gold = calc_two_dataframes(SP500, "/", gold)

    # Shiller P/E 10 ratio
    println("\n************** Shiller P/E 10 (Quandl) **************")
    ShillerPE10 = get_daily_data_from_quandl(quandl, "MULTPL/SHILLER_PE_RATIO_MONTH", date_plotstart, date_plotend, ["Date", "Value"], isverbose, nrows_verbose)

    # Tobin"s Q ratio
    # Nonfinancial Corporate Business; Corporate Equities; Liability, Level
    println("\n************** Nonfinancial Corporate Business; Corporate Equities; Liability, Level (Fred) **************")
    equity = get_daily_data_from_fred("MVEONWMVBSNNCB", date_plotstart, date_plotend, isverbose, nrows_verbose)
    # Nonfinancial Corporate Business; Net Worth, Level
    println("\n************** Nonfinancial Corporate Business; Net Worth, Level (Fred) **************")
    networth = get_daily_data_from_fred("TNWMVBSNNCB", date_plotstart, date_plotend, isverbose, nrows_verbose)
    TobinQ = calc_two_dataframes(equity, "/", networth)

    # Inflation - CPI
    println("\n************** CPI (Fred) **************")
    cpi = get_daily_data_from_fred("CPIAUCSL", date_plotstart, date_plotend, isverbose, nrows_verbose)
    # CPI By Group
    println("\n************** CPI: Food (Fred) **************")
    cpi_food = get_daily_data_from_fred("CPIUFDSL", date_plotstart, date_plotend, isverbose, nrows_verbose)
    println("\n************** CPI: Housing (Fred) **************")
    cpi_housing = get_daily_data_from_fred("CPIHOSSL", date_plotstart, date_plotend, isverbose, nrows_verbose)
    println("\n************** CPI: Medical (Fred) **************")
    cpi_medical = get_daily_data_from_fred("CPIMEDSL", date_plotstart, date_plotend, isverbose, nrows_verbose)
    println("\n************** CPI: Education (Fred) **************")
    cpi_education = get_daily_data_from_fred("CUSR0000SAE1", date_plotstart, date_plotend, isverbose, nrows_verbose)

    # Inflation - GDP Deflator
    println("\n************** GDP Deflator (Fred) **************")
    gdpdef = get_daily_data_from_fred("GDPDEF", date_plotstart, date_plotend, isverbose, nrows_verbose)

    # S&P500 / GDP Deflator
    SP500_gdpdef = calc_two_dataframes(SP500, "/", gdpdef)

    # Money Supply - Monetary Base, millions of dollars
    println("\n************** Monetary Base (Fred) **************")
    MB = get_daily_data_from_fred("BOGMBASEW", date_plotstart, date_plotend, isverbose, nrows_verbose)
    MB[!, :value] /= 1000  # convert millions of dollars to billions of dollars

    # Money Supply - M2, billions of dollars
    println("\n************** M2 (Fred) **************")
    M2 = get_daily_data_from_fred("M2", date_plotstart, date_plotend, isverbose, nrows_verbose)

    # S&P500 / M2
    SP500_M2 = calc_two_dataframes(SP500, "/", M2)

    # US Treasury Zero-Coupon Yield Curve
    println("\n************** US Treasury Zero-Coupon Yield Curve (Quandl) **************")
    treasury_yield1, treasury_yield2, treasury_yield5, treasury_yield10, treasury_yield15, treasury_yield20 =
        get_daily_data_from_quandl(quandl, "FED/SVENY", date_plotstart, date_plotend, ["Date", "SVENY01", "SVENY02", "SVENY05", "SVENY10", "SVENY15", "SVENY20"], isverbose, nrows_verbose)
    # short term interest rate (1 Yr) / long term interest rate (15 Yr)
    treasury_yield_spread = calc_two_dataframes(treasury_yield1, "/", treasury_yield15)

    # US GDP, billions of dollars
    println("\n************** GDP (Fred) **************")
    GDP = get_daily_data_from_fred("GDP", date_plotstart, date_plotend, isverbose, nrows_verbose)
    # Real GDP, Billions of Chained 2012 Dollars,
    println("\n************** Real GDP (Fred) **************")
    RealGDP = get_daily_data_from_fred("GDPC1", date_plotstart, date_plotend, isverbose, nrows_verbose)

    # S&P500 / GDP
    SP500_gdp = calc_two_dataframes(SP500, "/", GDP)

    # Deflated GDP
    GDP_deflated = calc_two_dataframes(GDP, "/", gdpdef)
    GDP_deflated[!, :value] *= 100

    # S&P500 / Deflated GDP
    SP500_deflgdp = calc_two_dataframes(SP500, "/", GDP_deflated)

    # Excess Monetary Base Explansion: MB / GDP
    MB_GDP = calc_two_dataframes(MB, "/", GDP)
    M2_GDP = calc_two_dataframes(M2, "/", GDP)
    # Normalized to pre-2008 era (1982 to May 2008) which was pretty flat
    MB_GDP_norm = copy(MB_GDP)
    MB_GDP_norm[!, :value] = MB_GDP[!, :value] / mean(filter(x -> Date(1982, 1, 1) < x.date < Date(2008, 5, 1), MB_GDP)[!, :value])

    # Adjust treasury yield spread using excess monetary base expansion
    treasury_yield_spread_adj = calc_two_dataframes(treasury_yield_spread, "*", MB_GDP_norm)

    # TED spread
    println("\n************** TED Spread (Fred) **************")
    ted_spread = get_daily_data_from_fred("TEDRATE", date_plotstart, date_plotend, isverbose, nrows_verbose)

    # VIX
    println("\n************** VIX (Yahoo) **************")
    vix = get_daily_data_from_yahoo("^VIX", date_plotstart, date_plotend, isverbose, nrows_verbose)

    # Financial Stress Indices
    # St. Louis Fed Financial Stress Index
    println("\n************** St. Louis Fed Financial Stress Index **************")
    stl_fsi = get_daily_data_from_fred("STLFSI", date_plotstart, date_plotend, isverbose, nrows_verbose)
    # Kansas City Financial Stress Index
    println("\n************** Kansas City Financial Stress Index **************")
    kc_fsi = get_daily_data_from_fred("KCFSI", date_plotstart, date_plotend, isverbose, nrows_verbose)
    # Cleveland Financial Stress Index
    println("\n************** Cleveland Financial Stress Index **************")
    c_fsi = get_daily_data_from_fred("CFSI", date_plotstart, date_plotend, isverbose, nrows_verbose)
    # Chicago Fed National Financial Conditions Index
    println("\n************** Chicago Fed National Financial Conditions Index **************")
    n_fci = get_daily_data_from_fred("NFCI", date_plotstart, date_plotend, isverbose, nrows_verbose)

    # Population
    println("\n************** US Population (Fred) **************")
    population = get_daily_data_from_fred("POP", date_plotstart, date_plotend, isverbose, nrows_verbose)  # unit: thousands of persons
    # Working age population (age 15-64)
    println("\n************** Working age population (age 15-64) (Fred) **************")
    wa_population = get_daily_data_from_fred("LFWA64TTUSM647N", date_plotstart, date_plotend, isverbose, nrows_verbose)  # unit: persons
    wa_population[!, :value] /= 1000  # convert unit to thousands of persons
    # population values below are all in thousands of persons
    println("\n************** White population (Fred) **************")
    ratio_white = calc_two_dataframes(get_daily_data_from_fred("LNU00000003", date_plotstart, date_plotend, isverbose, nrows_verbose), "/", population)
    println("\n************** Black population (Fred) **************")
    ratio_black = calc_two_dataframes(get_daily_data_from_fred("LNU00000006", date_plotstart, date_plotend, isverbose, nrows_verbose), "/", population)
    println("\n************** Hispanic population (Fred) **************")
    ratio_hispanic = calc_two_dataframes(get_daily_data_from_fred("LNU00000009", date_plotstart, date_plotend, isverbose, nrows_verbose), "/", population)
    println("\n************** Asian population (Fred) **************")
    ratio_asian = calc_two_dataframes(get_daily_data_from_fred("LNU00032183", date_plotstart, date_plotend, isverbose, nrows_verbose), "/", population)
    gdp_per_capita = calc_two_dataframes(GDP, "/", population)  # billions of dollars per thousands of persons
    gdp_per_capita[!, :value] *= (1e6 / 1e3)   # thousand $ per person
    realgdp_per_capita = calc_two_dataframes(RealGDP, "/", population)
    realgdp_per_capita[!, :value] *= (1e6 / 1e3)   # thousand $ per person

    # Labor Market
    # Civilian Employment-Population ratio
    println("\n************** Civilian Employment-Population ratio (Fred) **************")
    epr = get_daily_data_from_fred("EMRATIO", date_plotstart, date_plotend, isverbose, nrows_verbose)  # thousands
    # Civilian Unemployment Rate
    println("\n************** Civilian Unemployment Rate (Fred) **************")
    uer = get_daily_data_from_fred("UNRATE", date_plotstart, date_plotend, isverbose, nrows_verbose)  # thousands
    # Civilian Labor Force Participation Rate
    println("\n************** Civilian Labor Force Participation Rate (Fred) **************")
    lfpr = get_daily_data_from_fred("CIVPART", date_plotstart, date_plotend, isverbose, nrows_verbose)  # thousands

    # S&P/Case-Shiller Housing Indexes
    CaseShillerIndexID = Dict("City20"=>"SPCS20RSA", "Chicago"=>"CHXRSA", "SanFrancisco"=>"SFXRSA", "LosAngeles"=>"LXXRSA",
    "SanDiego"=>"SDXRSA", "NewYork"=>"NYXRSA", "Portland"=>"POXRSA", "Seattle"=>"SEXRSA",
    "Atlanta"=>"ATXRSA", "Boston"=>"BOXRSA", "Charlotte"=>"CRXRSA", "Cleveland"=>"CEXRSA",
    "Dallas"=>"DAXRSA", "Denver"=>"DNXRSA", "Detroit"=>"DEXRSA", "LasVegas"=>"LVXRSA", "Miami"=>"MIXRSA",
    "Minneapolis"=>"MNXRSA", "Phoenix"=>"PHXRSA", "Tampa"=>"TPXRSA", "WashingtonDC"=>"WDXRSA",
    "City10"=>"SPCS10RSA", "National"=>"CSUSHPISA")
    caseshiller = Dict()
    for city in cities_of_interest
        println("\n************** Case Shiller Index: " * city * " (Fred) **************")
        caseshiller[city] = get_daily_data_from_fred(CaseShillerIndexID[city], date_plotstart, date_plotend, isverbose, nrows_verbose)  # thousands
    end

    # Commodity and Energy Prices represented in Continuous Front-Month Future Contracts
    futures_underlying = ["USDIndex", "EURIndex", "JPYIndex", "5YrYield", "10YrYield", "Gold", "Silver", "Copper",
                            "CrudeOil", "BrentCrudeOil", "Gasoline", "NaturalGas", "Wheat", "Corn", "LiveCattle",
                            "Cotton", "Sugar", "Coffee", "Cocoa", "OrangeJuice"]

    futures_symbols = ["DX=F", "6E=F", "6J=F", "^FVX", "^TNX", "GC=F", "SI=F", "HG=F",
                        "CL=F", "BZ=F", "RB=F", "NG=F", "ZW=F", "ZC=F", "LE=F",
                        "CT=F", "SB=F", "KC=F", "CC=F", "OJ=F"]

    futures_contracts = Dict(zip(futures_underlying, futures_symbols))
    futures_prices = Dict()
    for comdty in futures_underlying
        println("\n************** Futures: " * comdty * " (Yahoo) **************")
        futures_prices[comdty] = get_daily_data_from_yahoo(futures_contracts[comdty], date_plotstart, date_plotend, isverbose, nrows_verbose)
    end

    println("\nAll downloads finished!")
    isdownloaded = true
end

if isdownloaded
    println("\nPlotting...")
    nfig = 0
    #########################################################################
    # Plots

    fig = figure(facecolor="w", figsize=figsize, dpi=dpi)

    nrows = 3
    ncols = 2

    #...........................
    todaystr = Dates.format(today(), fmt)

    nplot = 0

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    ax.plot(SP500[!, :date], SP500[!, :value], "b-", label="S&P500")
    ax.set_xlim([date_plotstart, date_plotend])
    ax.legend(prop=Dict("size" => legend_fontsize), loc="upper left")
    ax2 = ax.twinx()
    ax2.plot(gold[!, :date], gold[!, :value], "r-", label="Gold")
    ax2.set_ylabel("USD/OZ")
    ax2.set_xlim([date_plotstart, date_plotend])
    ax2.legend(prop=Dict("size" => legend_fontsize), loc="lower right")
    grid(true, linestyle=":")
    PyPlot.title("Stock and Gold as of " * todaystr)

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    ax.plot(SP500_gold[!, :date], SP500_gold[!, :value], "b-")
    ax.set_xlim([date_plotstart, date_plotend])
    # ax.legend(prop=Dict("size" => legend_fontsize), loc="upper left")
    grid(true, linestyle=":")
    PyPlot.title("S&P500 / Gold as of " * todaystr)

    #...........................

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    baseline_yearsago = 10
    baseline_year = year(Dates.today()) - baseline_yearsago
    baseline_date = Date(baseline_year, month(Dates.today()), day(Dates.today()))
    baseline_cpi = filter(x->x.date >= baseline_date, cpi)[1, :value]
    ax.plot(cpi[!, :date], cpi[!, :value] / baseline_cpi * 100, ":", color="blue", linewidth=4, label="CPI")
    baseline_cpi_food = filter(x->x.date >= baseline_date, cpi_food)[1, :value]
    ax.plot(cpi_food[!, :date], cpi_food[!, :value] / baseline_cpi_food * 100, "-", color="green", label="CPI:Food")
    baseline_cpi_housing = filter(x->x.date >= baseline_date, cpi_housing)[1, :value]
    ax.plot(cpi_housing[!, :date], cpi_housing[!, :value] / baseline_cpi_housing * 100, "-", color="aqua", label="CPI:Housing")
    baseline_cpi_medical = filter(x->x.date >= baseline_date, cpi_medical)[1, :value]
    ax.plot(cpi_medical[!, :date], cpi_medical[!, :value] / baseline_cpi_medical * 100, "-", color="magenta", label="CPI:Medical")
    baseline_cpi_education = filter(x->x.date >= baseline_date, cpi_education)[1, :value]
    ax.plot(cpi_education[!, :date], cpi_education[!, :value] / baseline_cpi_education * 100, "-", color="gold", label="CPI:Education")
    baseline_gdpdef = filter(x->x.date >= baseline_date, gdpdef)[1, :value]
    ax.plot(gdpdef[!, :date], gdpdef[!, :value] / baseline_gdpdef * 100, ":", color="red", linewidth=4, label="GDP Deflator")
    ax.set_xlim([date_plotstart, date_plotend])
    ax.set_ylabel("Index")
    ax.legend(prop=Dict("size" => legend_fontsize), loc="upper left")
    grid(true, linestyle=":")
    plt.title("Inflation Index ($baseline_yearsago years ago = 100) as of " * todaystr)

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    ax.plot(SP500_gdp[!, :date], SP500_gdp[!, :value], "b-", label="S&P500 / GDP")
    ax.set_xlim([date_plotstart, date_plotend])
    ax.legend(prop=Dict("size" => legend_fontsize), loc="upper left")
    ax2 = ax.twinx()
    ax2.plot(SP500_M2[!, :date], SP500_M2[!, :value], "r-", label="S&P500 / M2")
    ax2.set_xlim([date_plotstart, date_plotend])
    ax2.legend(prop=Dict("size" => legend_fontsize), loc="upper right")
    plt.title("S&P500 / M2 as of " * todaystr)
    grid(true, linestyle=":")
    plt.title("S&P500 vs. GDP and M2 as of " * todaystr)

    #...........................

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    ax.plot(MB[!, :date], MB[!, :value], "b-", label="MB")
    ax.plot(M2[!, :date], M2[!, :value], "r-", label="M2")
    ax.plot(GDP[!, :date], GDP[!, :value], "-", label="GDP")
    ax.plot(RealGDP[!, :date], RealGDP[!, :value], "-", color="darkgreen", label="Real GDP (2009 USD)")
    ax.set_xlim([date_plotstart, date_plotend])
    ax.set_ylabel("USD Bln")
    ax.legend(prop=Dict("size" => legend_fontsize), loc="upper left")
    ax2 = ax.twinx()
    ax2.plot(MB_GDP[!, :date], MB_GDP[!, :value], "m--", label="MB/GDP")
    ax2.plot(M2_GDP[!, :date], M2_GDP[!, :value], "c--", label="M2/GDP")
    ax2.set_xlim([date_plotstart, date_plotend])
    ax2.legend(prop=Dict("size" => legend_fontsize), loc="upper center")
    grid(true, linestyle=":")
    plt.title("Money Supply and GDP as of " * todaystr)

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    ax.plot(ShillerPE10[!, :date], ShillerPE10[!, :value], "b-", label="Shiller P/E 10")
    ax.set_xlim([date_plotstart, date_plotend])
    ax.legend(prop=Dict("size" => legend_fontsize), loc="upper left")
    ax2 = ax.twinx()
    ax2.plot(TobinQ[!, :date], TobinQ[!, :value], "r-", label="Tobin's Q")
    ax2.set_xlim([date_plotstart, date_plotend])
    ax2.legend(prop=Dict("size" => legend_fontsize), loc="upper right")
    grid(true, linestyle=":")
    plt.title("Stock Market Valuation - Ratios as of " * todaystr)

    tight_layout()

    nfig += 1
    savefig(joinpath(@__DIR__, "pictures", "BigPicture$nfig.png"))

    #-------------------------------------------------------------------------

    figure(facecolor="w", figsize=figsize, dpi=dpi)

    nrows = 3
    ncols = 2
    nplot = 0

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    ax.plot(treasury_yield20[!, :date], treasury_yield20[!, :value], "-",  color="blue", linewidth=0.5, label="20 Yr")
    ax.plot(treasury_yield10[!, :date], treasury_yield10[!, :value], "-",  color="green", linewidth=0.5, label="10 Yr")
    ax.plot(treasury_yield5[!, :date], treasury_yield5[!, :value], "-",  color="yellow", linewidth=0.5, label="5 Yr")
    ax.plot(treasury_yield2[!, :date], treasury_yield2[!, :value], "-", color="orange", linewidth=0.5, label="2 Yr")
    ax.plot(treasury_yield1[!, :date], treasury_yield1[!, :value], "-", color="red", linewidth=0.5, label="1 Yr")
    ax.set_xlim([date_plotstart, date_plotend])
    ax.set_ylabel("%")
    ax.legend(prop=Dict("size" => legend_fontsize), loc="upper left")
    grid(true, linestyle=":")
    plt.title("Treasury Zero-Coupon Yield as of " * todaystr)

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    ax.plot(SP500[!, :date], SP500[!, :value], "-", color="blue", label="S&P500")
    ax.set_xlim([date_plotstart, date_plotend])
    # ax.set_yscale("log")
    # # ax.yaxis.set_major_formatter(fmt)
    ax.legend(prop=Dict("size" => legend_fontsize), loc="upper left")
    ax2 = ax.twinx()
    ax2.plot(treasury_yield_spread[!, :date], treasury_yield_spread[!, :value], "-",
                color="magenta", linewidth=1, label="1Yr/15Yr")
    ax2.plot(treasury_yield_spread_adj[!, :date], treasury_yield_spread_adj[!, :value], "-",
                color="cyan", linewidth=1, label="1Yr/15Yr_MBAdj")
    ax2.plot([treasury_yield_spread[!, :date][1], treasury_yield_spread[!, :date][end]], [1, 1], "-",
                color="red")
    ax2.set_xlim([date_plotstart, date_plotend])
    ax2.legend(prop=Dict("size" => legend_fontsize), loc="upper center")
    grid(true, linestyle=":")
    plt.title("Stock Market and Interest Rate Structure as of " * todaystr)

    #...........................

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    # ax.plot(SP500[!, :date], SP500[!, :value], "-", color="blue", label="S&P500")
    ax.plot(vix[!, :date], vix[!, :value], "-", color="blue", label="VIX")
    ax.set_xlim([date_plotstart, date_plotend])
    # ax.set_yscale("log")
    # # ax.yaxis.set_major_formatter(fmt)
    ax.legend(prop=Dict("size" => legend_fontsize), loc="upper left")
    ax2 = ax.twinx()
    ax2.plot(ted_spread[!, :date], ted_spread[!, :value], "-",
                color="magenta", linewidth=1, label="TED Spread")
    ax2.set_xlim([date_plotstart, date_plotend])
    ax2.set_ylabel("%")
    # ax.set_yscale("log")
    ax2.legend(prop=Dict("size" => legend_fontsize), loc="upper right")
    grid(true, linestyle=":")
    plt.title("Financial Stress Indicators (1) as of " * todaystr)

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    ax.plot(SP500[!, :date], SP500[!, :value], "-", color="blue", label="S&P500")
    ax.set_xlim([date_plotstart, date_plotend])
    # ax.set_yscale("log")
    # # ax.yaxis.set_major_formatter(fmt)
    ax.legend(prop=Dict("size" => legend_fontsize), loc="upper left")
    ax2 = ax.twinx()
    ax2.plot(stl_fsi[!, :date], stl_fsi[!, :value], "-",
                color="red", linewidth=1, label="St. Louis Fed FSI")
    ax2.plot(kc_fsi[!, :date], kc_fsi[!, :value], "-",
                color="green", linewidth=1, label="Kansas City FSI")
    ax2.set_xlim([date_plotstart, date_plotend])
    # ax.set_yscale("log")
    ax2.legend(prop=Dict("size" => legend_fontsize), loc="upper center")
    grid(true, linestyle=":")
    plt.title("Financial Stress Indicators (2) as of " * todaystr)

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    tedspreadvix = innerjoin(ted_spread, vix, on=:date, makeunique=true)
    ax.plot(tedspreadvix[!, :value_1], tedspreadvix[!, :value], "b.")
    ax.set_xlabel("VIX")
    ax.set_ylabel("TED Spread")
    corr = cor(tedspreadvix[!, :value_1], tedspreadvix[!, :value])
    grid(true, linestyle=":")
    plt.title("VIX ~ TED Spread, corr = $(round(corr*100, digits=2))%")

    #...........................

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    ax.plot(SP500[!, :date], SP500[!, :value], "-", color="blue", label="S&P500")
    ax.set_xlim([date_plotstart, date_plotend])
    # ax.set_yscale("log")
    # ax.yaxis.set_major_formatter(fmt)
    ax.legend(prop=Dict("size" => legend_fontsize), loc="upper left")
    ax2 = ax.twinx()
    ax2.plot(c_fsi[!, :date], c_fsi[!, :value], "-",
                color="red", linewidth=1, label="Cleveland FSI")
    ax2.plot(n_fci[!, :date], n_fci[!, :value], "-",
                color="green", linewidth=1, label="Chicago Fed National FCI")
    ax2.set_xlim([date_plotstart, date_plotend])
    # ax.set_yscale("log")
    ax2.legend(prop=Dict("size" => legend_fontsize), loc="upper center")
    grid(true, linestyle=":")
    plt.title("Financial Stress Indicators (3) as of " * todaystr)

    tight_layout()

    nfig += 1
    savefig(joinpath(@__DIR__, "pictures", "BigPicture$nfig.png"))
    #-----------------------------------------------------------------------

    figure(facecolor="w", figsize=figsize, dpi=dpi)

    nrows = 2
    ncols = 2
    nplot = 0

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    ax.plot(population[!, :date], population[!, :value]/1e6, "--",
            color="blue", linewidth=1, label="Population")
    ax.plot(wa_population[!, :date], wa_population[!, :value]/1e6, "--",
            color="magenta", linewidth=1, label="Working Age (15-64) Population")
    ax.set_ylabel("Bln Persons")
    ax.set_xlim([date_plotstart, date_plotend])
    # ax.set_yscale("log")
    # ax.yaxis.set_major_formatter(fmt)
    ax.legend(prop=Dict("size" => legend_fontsize), loc="upper left")
    ax2 = ax.twinx()
    ax2.plot(gdp_per_capita[!, :date], gdp_per_capita[!, :value], "-",
                color="red", linewidth=1, label="GDP Per Capita")
    ax2.plot(realgdp_per_capita[!, :date], realgdp_per_capita[!, :value], "-",
                color="green", linewidth=1, label="Real GDP Per Capita (2009 USD)")
    ax2.set_ylabel("K USD")
    ax2.set_xlim([date_plotstart, date_plotend])
    # ax.set_yscale("log")
    ax2.legend(prop=Dict("size" => legend_fontsize), loc="lower right")
    grid(true, linestyle=":")
    plt.title("Population as of " * todaystr)

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    ax.plot(epr[!, :date], epr[!, :value], "-",
            color="green", linewidth=1, label="Employment-Population Ratio")
    ax.plot(lfpr[!, :date], lfpr[!, :value], "-",
            color="blue", linewidth=1, label="Labor Force Participation Rate")
    ax.set_xlim([date_plotstart, date_plotend])
    ax.set_ylabel("%")
    # ax.set_yscale("log")
    # ax.yaxis.set_major_formatter(fmt)
    ax.legend(prop=Dict("size" => legend_fontsize), loc="lower left")
    ax2 = ax.twinx()
    ax2.plot(uer[!, :date], uer[!, :value], "-",
                color="red", linewidth=1, label="Unemployment Rate")
    ax2.set_xlim([date_plotstart, date_plotend])
    ax2.set_ylabel("%")
    # ax.set_yscale("log")
    ax2.legend(prop=Dict("size" => legend_fontsize), loc="upper center")
    grid(true, linestyle=":")
    plt.title("Labor Market Condition as of " * todaystr)

    nplot += 1
    ax = subplot(nrows, ncols, nplot)
    for (n, city) in enumerate(cities_of_interest)
        ax.plot(caseshiller[city][!, :date], caseshiller[city][!, :value], "-", color=colors[n],
                linewidth =  city in ["National", "Chicago", "SanFrancisco"] ? 3 : 1, label=city)
    end
    ax.set_xlim([date_plotstart, date_plotend])
    # ax.set_yscale("log")
    # ax.yaxis.set_major_formatter(fmt)
    ax.legend(prop=Dict("size" => legend_fontsize), loc="upper left")
    grid(true, linestyle=":")
    plt.title("S&P/Case-Shiller Home Price Indices as of " * todaystr)

    tight_layout()

    nfig += 1
    savefig(joinpath(@__DIR__, "pictures", "BigPicture$nfig.png"))
    #----------------------------------------------------------------------

    figure(facecolor="w", figsize=figsize, dpi=dpi)

    nrows = 4
    ncols = 5
    nplot = 0

    for comdty in futures_underlying
        global nplot
        local ax
        nplot += 1
        ax = subplot(nrows, ncols, nplot)
        ax.plot(futures_prices[comdty][!, :date], futures_prices[comdty][!, :value], "-",
                color="blue", linewidth=1, label=comdty)
        ax.set_xlim([date_future_plotstart_long, date_plotend])
        # ax.yaxis.set_major_formatter(fmt)
        ax.legend(prop=Dict("size" => legend_fontsize), loc=0)
        grid(true, linestyle=":")
        if nplot == Int(round(ncols/2+0.01, digits=0))
            ax.set_title("Futures - Long Term ($future_yrs_long-year) as of $todaystr")
        end
        tick_params(axis="both", which="major", labelsize=8)
        tick_params(axis="both", which="minor", labelsize=8)
    end

    tight_layout()

    nfig += 1
    savefig(joinpath(@__DIR__, "pictures", "BigPicture$nfig.png"))
    #--------------------------------------------------------------------------
    figure(facecolor="w", figsize=figsize, dpi=dpi)

    nrows = 4
    ncols = 5
    nplot = 0

    for comdty in futures_underlying
        global nplot
        local ax
        nplot += 1
        ax = subplot(nrows, ncols, nplot)
        ind = (date_future_plotstart_short .<= futures_prices[comdty][!, :date] .<= date_plotend)
        ax.plot(futures_prices[comdty][ind, :date], futures_prices[comdty][ind, :value], "-",
                color="red", linewidth=1, label=comdty)
        ax.set_xlim([date_future_plotstart_short, date_plotend])
        # ax.yaxis.set_major_formatter(fmt)
        ax.legend(prop=Dict("size" => legend_fontsize), loc=0)
        grid(true, linestyle=":")
        if nplot == Int(round(ncols/2+0.01, digits=0))
            ax.set_title("Futures - Short Term ($future_yrs_short-year) as of $todaystr")
        end
        tick_params(axis="both", which="major", labelsize=7)
        tick_params(axis="both", which="minor", labelsize=7)
    end

    tight_layout()

    nfig += 1
    savefig(joinpath(@__DIR__, "pictures", "BigPicture$nfig.png"))
    #------------------------------------------------------------------------

    println("Plots opened in new windows.\nThe script finished.")
else
    println("Downloading failed!\nThe script aborted.")
end


text:
struct Vec3{T} <: FieldVector{3, T}
    x::T
    y::T
    z::T
end

StaticArrays.similar_type(::Type{<:Vec3}, ::Type{T}, ::Size{(3,)}) where {T} = Vec3{T}

function gaussian_product_center(
        a::Float64,
        RA::Vec3{Float64},
        b::Float64,
        RB::Vec3{Float64},
    )
    p = a + b
    return (a * RA + b * RB) / p
end

function norm2(R::Vec3{Float64})
    s = 0.0
    for r in R
        s += abs2(r)
    end
    return s
end

gamma_inc_p(a, x) = gamma_inc(a, x, 0)[1]

function boys(n::Int, x::Float64)
    if x < 1e-7
        return 1 / (2 * n + 1) - x / (2 * n + 3) + x^2 / (2 * n + 5) / 2
    else
        return 1 / (2 * x^(n + 0.5)) * gamma(n + 0.5) * gamma_inc_p(n + 0.5, x)
    end
end


text:
# Solves the Lorenz problem over different periods, and writes time average
module LorenzMod

using DelimitedFiles, DifferentialEquations, GeometricIntegratorsDiffEq

include("gregory.jl")
using .Gregory

function calcZAverage(sol, ibeg::Int, iend::Int,
    num_steps::Int, order::Int)::Float64
    @assert( iend - ibeg == num_steps)
    #data = zeros(num_steps+1)
    #for i = 1:num_steps+1
    #    data[i] = u[ibeg+i-1][3]
    #end
    numbnd = order-1
    if 2*numbnd > num_steps
        numbnd = div(num_steps, 2)
        order = numbnd + 1
        println("WARNING: had to adopt order ", order, " quadrature.")
    end
    val = integrate(sol[3,ibeg:iend], numbnd, order, h = 1/num_steps)
    return val
end

"""
Print3DArray(file, data)

Print a 3D array to file using column-major ordering.
"""
function Print3DArray(file, arr::AbstractArray{Float64,3})
    for k = 1:size(arr,3)
        for j = 1:size(arr,2)
            for i = 1:size(arr,1)
                print(file, arr[i,j,k], " ")
            end
            println(file)
        end
    end
end

# define the Lorenz problem
function lorenz!(du,u,p,t)
    du[1] = 10.0*(u[2]-u[1])
    du[2] = u[1]*(28.0-u[3]) - u[2]
    du[3] = u[1]*u[2] - (8/3)*u[3]
end

"""
gatherStatsLorenz()

Runs the Lorenz problem with different discretization orders, steps sizes and
integration periods, and writes results to file.
"""
function gatherStatsLorenz()

    u0 = [-5.276131687990789; -8.274832825254473; 16.36301350713599]
    time = 1000.0
    probLorenz = ODEProblem(lorenz!, u0, (0, time))

    time_fac = 10
    #num_steps = [6250; 12500; 25000; 50000; 100000]
    base_time = time/time_fac^3
    base_steps = [10; 20; 40; 80; 160]
    order = [2, 4, 8]
    #file_prefix = ["Heun_u_", "RK4_u_", "GLRK(4)_u_"]
    method = [Heun(), RK4(), Vern8() ]

    # loop over the different order methods
    for p = 1:size(order,1)
        println("Gathering statistics for order ",order[p])
        # each method gets its own file

        dt = zeros(Float64, (time_fac^3, size(base_steps,1), 3) )
        z_avg = zeros(size(dt))
        cputime = zeros(size(dt))

        println("Initial compile...")
        cpu_full = @elapsed sol = solve(probLorenz, method[p],
                                        dt=base_time/base_steps[1],
                                        adaptive=false)

        for n = 1:size(base_steps,1)
            dt[:,n,:] .= base_time/base_steps[n]
            println("\tSolve using num steps = ",time_fac^3*base_steps[n])
            cpu_full = @elapsed sol = solve(probLorenz, method[p],
                                            dt=base_time/base_steps[n],
                                            adaptive=false)
            #if any(isnan, sol)
            if sol.retcode != :Success
                println("\t\tLorenz failed to solve due to instability.")
                continue
            end

            # We partition the simulation into smaller periods for statistics
            # gathering.  For example, if the simulation time is 1000 units
            # then we have 10 samples for 100, 100 samples for 10 units, and
            # 1000 samples for 1 unit.
            # The statistics for the medium and long periods are duplicated to
            # simplify the file I/O (NumPy loadtxt is easier)
            ptr = 1
            for s1 = 1:time_fac
                ibeg = base_steps[n]*(s1-1)*time_fac*time_fac + 1
                iend = base_steps[n]*s1*time_fac*time_fac + 1
                z_avg_long = calcZAverage(sol, ibeg, iend,
                                          time_fac^2*base_steps[n],
                order[p])
                for s2 = 1:time_fac
                    ibeg = base_steps[n]*((s1-1)*time_fac + s2 - 1)*time_fac + 1
                    iend = base_steps[n]*((s1-1)*time_fac + s2)*time_fac + 1
                    z_avg_med = calcZAverage(sol, ibeg, iend,
                    time_fac*base_steps[n], order[p])
                    for s3 = 1:time_fac
                        # get statistics for the shortest period
                        ibeg = base_steps[n]*(
                        ( (s1-1)*time_fac + s2-1)*time_fac + s3 - 1) + 1
                        iend = base_steps[n]*(
                        ( (s1-1)*time_fac + s2-1)*time_fac + s3) + 1
                        z_avg[ptr,n,1] = calcZAverage(sol, ibeg, iend,
                        base_steps[n], order[p])
                        cputime[ptr,n,1] = cpu_full/(time_fac^3)

                        # get statistics for the medium period
                        ibeg = base_steps[n]*((s1-1)*time_fac+s2-1)*time_fac + 1
                        iend = base_steps[n]*((s1-1)*time_fac + s2)*time_fac + 1
                        z_avg[ptr,n,2] = z_avg_med
                        cputime[ptr,n,2] = cpu_full/(time_fac^2)

                        # get statistics for the long period
                        z_avg[ptr,n,3] = z_avg_long
                        cputime[ptr,n,3] = cpu_full/time_fac

                        ptr += 1
                    end
                end
            end
        end

        # write the data
        data = open("raw_statistics_lorenz_order$(order[p]).dat", "w")
        Print3DArray(data, dt)
        Print3DArray(data, z_avg)
        Print3DArray(data, cputime)
        close(data)
    end
end

end # LorenzMod

max_stars_repo_name: jehicken/high-order-and-chaos

max_issues_repo_name: jehicken/high-order-and-chaos

max_forks_repo_name: jehicken/high-order-and-chaos