text:
#=
plot_3d:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-08-12
=#

module Plots3D
    # packages
    using PyPlot: plt
    using PyCall
    @pyimport mpl_toolkits.mplot3d as mpl3

    # flag to switch show or not
    show_plot = true

    # methods
    function set_show_plot(flag)
        global show_plot = flag
    end

    function plot_point_3d()
        x = randn(100)
        y = randn(100)
        z = randn(100)
        fig = plt.figure()
        ax = fig.add_subplot(1, 1, 1, projection="3d")
        ax.plot(x, y, z, "o", color="r", ms=8)
        ax.set_xlabel("X")
        ax.set_ylabel("Y")
        ax.set_zlabel("Z")
        ax.set_title("3D Point Plot Sample")
        ax.grid(true)
    end

    function plot_line_3d()
        x = 1:5
        y = 1:5
        z = 1:5
        fig = plt.figure()
        ax = fig.add_subplot(1, 1, 1, projection="3d")
        ax.plot(x, y, z, "o-", color="r", ms=5)
        ax.set_xlabel("X")
        ax.set_ylabel("Y")
        ax.set_zlabel("Z")
        ax.set_title("3D Line Plot Sample")
        ax.grid(true)
    end

    function main()
        plot_point_3d()
        plot_line_3d()
        if show_plot == true
            plt.show()
        end
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .Plots3D
    Plots3D.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
using Remark
using Test

demo = joinpath(@__DIR__, "..", "demo")
isdir(demo) && rm(demo, recursive = true)
sleep(1)

mkdir(demo)
filenames = ["example.jl", "example.md"]
cp.(joinpath.(@__DIR__, "..", "examples", filenames), joinpath.(@__DIR__, "..", "demo", filenames), force = true)
sleep(1)

@testset "deps" begin
    for fn in ["remark.min.js", "katex.min.js", "katex.min.css", "auto-render.min.js"]
        @test isfile(joinpath(@__DIR__, "..", "assets", fn))
    end
end

const katexfonts = [
    "KaTeX_AMS-Regular",
    "KaTeX_Caligraphic-Bold",
    "KaTeX_Caligraphic-Regular",
    "KaTeX_Fraktur-Bold",
    "KaTeX_Fraktur-Regular",
    "KaTeX_Main-Regular",
    "KaTeX_Main-Bold",
    "KaTeX_Main-Italic",
    "KaTeX_Main-BoldItalic",
    "KaTeX_Math-Italic",
    "KaTeX_Math-BoldItalic",
    "KaTeX_SansSerif-Bold",
    "KaTeX_SansSerif-Italic",
    "KaTeX_SansSerif-Regular",
    "KaTeX_Script-Regular",
    "KaTeX_Size1-Regular",
    "KaTeX_Size2-Regular",
    "KaTeX_Size3-Regular",
    "KaTeX_Size4-Regular",
    "KaTeX_Typewriter-Regular",
]

const exts =["ttf", "woff", "woff2"]

@testset "katexfonts" begin
    for font in katexfonts
        for ext in exts
            @test isfile(joinpath(@__DIR__, "..", "assets", "fonts", join([font, ext], '.')))
        end
    end
end

isalnum(x) = isnumeric(x) || isletter(x)

@testset "slideshow" begin
    slideshowdir = Remark.slideshow(joinpath(demo, "example.jl"))
    @test slideshowdir == realpath(abspath(demo))
    @test all(isdir, joinpath.(slideshowdir, "build", "fonts", ["Lora", "Ubuntu_Mono", "Yanone_Kaffeesatz"]))
    for file in Remark.depfiles
        @test isfile(joinpath(slideshowdir, "build", splitdir(file)[2]))
    end
    v1 = filter(isalnum, read(joinpath(@__DIR__, "indexjl.html"), String))
    v2 = filter(isalnum, read(joinpath(demo, "build", "index.html"), String))
    @test split(v1, "statplotsvg")[1] == split(v2, "statplotsvg")[1]
    @test split(v1, "textarea")[3] == split(v2, "textarea")[3]
end

@testset "slideshowmd" begin
    slideshowdir = Remark.slideshow(joinpath(demo, "example.md"))
    @test slideshowdir == realpath(abspath(demo))
    @test all(isdir, joinpath.(slideshowdir, "build", "fonts", ["Lora", "Ubuntu_Mono", "Yanone_Kaffeesatz"]))
    for file in Remark.depfiles
        @test isfile(joinpath(slideshowdir, "build", splitdir(file)[2]))
    end
    v1 = filter(isalnum, read(joinpath(@__DIR__, "indexmd.html"), String))
    v2 = filter(isalnum, read(joinpath(demo, "build", "index.html"), String))
    @test split(v1, "statplotsvg")[1] == split(v2, "statplotsvg")[1]
    @test split(v1, "textarea")[3] == split(v2, "textarea")[3]
end

rm(demo, recursive = true)

ext: jl


text:
function force_dims(metnet)

    function _similar_copy(col, fill, expected_len)
        curr_len = length(col)
        if curr_len < expected_len
            scol = similar(col, expected_len)
            scol[curr_len:end] .= fill
            scol[begin:curr_len] = col
        elseif curr_len > expected_len
            scol = similar(col, expected_len)
            scol[begin:expected_len] = col[begin:expected_len]
        else
            scol = copy(col)
        end
        return scol
    end

    M, N = size(metnet)
    net = Dict()
    net[:S] = copy(metnet.S)
    net[:b] = _similar_copy(metnet.b, 0.0, M)
    net[:c] = _similar_copy(metnet.c, 0.0, N)
    net[:lb] = _similar_copy(metnet.lb, 0.0, N)
    net[:ub] = _similar_copy(metnet.ub, 0.0, N)
    net[:rxns] = copy(metnet.rxns)
    net[:mets] = copy(metnet.mets)

    for (l, k) in [(N, :rxns), (M, :mets)]
        length(net[k]) != l && error("Unsolvable miss-match $k ($(length(net[k]))) != $l")
    end

    return MetNet(metnet; reshape = false, net...)
end


text:
# Set up problem
using TrajectoryOptimization: to_array, update_constraints!, update_jacobians!, outer_loop_update_type
using LinearAlgebra
using Plots
using Test
using Logging

u_bound = 5.
model, obj = TrajectoryOptimization.Dynamics.pendulum!
opts = TrajectoryOptimization.SolverOptions()
opts.verbose = false
obj
obj.Q = 1e-3*Diagonal(I,2)
obj.R = 1e-2*Diagonal(I,1)
obj.tf = 2.
obj_c = TrajectoryOptimization.ConstrainedObjective(obj, u_min=-u_bound, u_max=u_bound) # constrained objective
obj_min = update_objective(obj_c,tf=:min,c=0.0, Q = obj.Q, R = obj.R, Qf = obj.Qf*0.)
dt = 0.1
n,m = model.n, model.m

solver = Solver(model,obj_c,integration=:rk3,dt=dt)
X_interp = line_trajectory(solver)
U = ones(m,solver.N)
results,stats = solve(solver,U)
plot(to_array(results.X)')
plot(to_array(results.U)')

solver_min = Solver(model,obj_min,integration=:rk4,N=41)
U = ones(m,solver_min.N)
solver_min.opts.verbose = true
solver_min.opts.use_static = false
solver_min.opts.max_dt = 0.2
solver_min.opts.constraint_tolerance = 0.5e-3
solver_min.opts.R_minimum_time = 1000.0
solver_min.opts.bp_reg_initial = 0
solver_min.opts.constraint_decrease_ratio = .25
solver_min.opts.penalty_scaling = 10.0
results_min,stats_min = solve(solver_min,U)
total_time(solver_min,results_min)
plot(to_array(results_min.X)[1:2,:]')
plot(to_array(results_min.U)[1:2,:]')
plot(stats_min["cost"])

TrajectoryOptimization.cost_mintime(solver_min,results_min)
max_violation(results_min)
results_min.CN
results_min.C[3]
results_min.U[1][2]
get_num_constraints(solver_min)

violations = map((x)->x.>0, results_min.Iμ) .* results_min.C
norm.(violations,Inf)

norm(to_array(violations),Inf)
max_violation(results_min)
argmax(abs.(to_array(violations)))

# Test dynamics change with dt
solver = Solver(model,obj_c,dt=dt)
solver_min = Solver(model,obj_min,N=solver.N)
@test get_num_controls(solver) == (1,1)
@test get_num_controls(solver_min) == (2,2)

x = rand(n)
u = rand(m)
S = [x;u;√dt]
xdot = zero(x)
Sdot = zero(S)

f = model.f
fd! = rk4(f,dt)
fd_aug = rk4(f_augmented!(f,n,m))
fd!(xdot,x,u)
fd_aug(Sdot,S)

fS = zeros(n+m+1,n+m+1)
ForwardDiff.jacobian!(fS,fd_aug,Sdot,S)

dft = zeros(n,1)
ft(dx,t) = fd!(dx,x,u,t[1])
ForwardDiff.jacobian!(dft,ft,xdot,[dt])

@test ~isapprox(dft,fS[1:n,end])  # make sure the augmented time jacobian is not the same as the not augmented, since S stores the square root of dt (for mintime problems)

u_dt = [u;sqrt(dt)]
xdot = zeros(n)
xdot0 = zeros(n)
solver_min.fd(xdot,x,u_dt,u_dt[end]^2)
solver.fd(xdot0,x,u)
fx,fu = solver_min.Fd(x,u_dt)
fx0,fu0 = solver.Fd(x,u)
@test xdot == xdot0
@test fx0 == fS[1:n,1:n]
@test fu0 ≈ fS[1:n,n+1]
@test fx == fx0
@test fu[:,1] ≈ fu0
@test fu ≈ fS[1:n,n+1:end]

# Change dt and make sure they are no longer equal
u_dt = [u;0.15]
solver_min.fd(xdot,x,u_dt,u_dt[end])
solver.fd(xdot0,x,u)
fx,fu = solver_min.Fd(x,u_dt)
fx0,fu0 = solver.Fd(x,u)
@test xdot ≉ xdot0
@test fx ≉ fx0
@test fu[:,1] ≉ fu0

# Run through parts of the iLQR algorithm
solver = Solver(model,obj_c,dt=dt)
solver_min = Solver(model,obj_min,N=solver.N)
n,m,N = get_sizes(solver)
p,pI,pE = TrajectoryOptimization.get_num_constraints(solver)
res_reg = ConstrainedVectorResults(n,m,p,N)
res_min = ConstrainedVectorResults(n,m+1,p+3,N)

res_reg.X[1] = obj.x0
res_min.X[1] = obj.x0

copyto!(res_reg.U, ones(m,N))
copyto!(res_min.U, [ones(m,N); ones(m,N)*sqrt(dt)])
@test to_array(res_min.U)[1,:] == vec(to_array(res_reg.U))

# Test that rollouts are the same
rollout!(res_reg, solver)
rollout!(res_min, solver_min)
@test maximum(norm.(res_reg.X .- res_min.X,Inf)) == 0

# Test constraints
res_min.U[2][2] = 0.15
res_min.U[N-1][2] = 0.2

update_constraints!(res_reg,solver)
update_constraints!(res_min,solver_min)

res_min.C[1]
res_reg.C[1]
res_reg.C[1] == res_min.C[1][[1,3]]
@test res_min.C[1][end] ≈ √dt-0.15
@test res_min.C[2][end] ≈ .15-√dt
@test res_min.C[N-2][end] == sqrt(dt)-0.2
@test res_min.C[N-1][end] == 0
res_min.Iμ[1][end,end] == 1
@test res_min.MU[1] == ones(5)
res_min.U[N-1][2] = √dt
res_min.U[2][2] = √dt

update_jacobians!(res_reg,solver)
update_jacobians!(res_min,solver_min)

@test size(res_min.fu[1]) == (n,m+1)
@test res_reg.fx[1] == res_min.fx[1]
@test res_reg.fu[1] ≈ res_min.fu[1][:,1]
@test res_reg.fu[N-1] ≈ res_min.fu[N-1][:,1]

res_min.Cu[1]
@test res_min.Cu[1] == [Matrix(I,2,2); -Matrix(I,2,2); 0 1]
@test res_min.Cu[N-1] == [Matrix(I,2,2); -Matrix(I,2,2); 0 0]

solver_min.obj
_cost(solver_min,res_min)
total_time(solver_min,res_min)

# Test cost is the same (when c=0 and R[m̄,m̄] = 0)
obj_min_c0 = update_objective(obj_min,c=0.,Q = obj_c.Q, R = obj_c.R, Qf = obj.Qf)
solver_min = Solver(model,obj_min_c0,N=solver.N)
solver_min.opts.min_time_regularization = 0
Xrand = rand(n,N)
copyto!(res_reg.X,Xrand)
copyto!(res_min.X,Xrand)

J_reg = _cost(solver,res_reg)
J_min = _cost(solver_min,res_min)
@test J_reg == J_min

# Make sure constrained cost is NOT the same
J_reg = cost(solver,res_reg)
J_min = cost(solver_min,res_min)
@test abs(J_reg-J_min) > 1e-2

# Test rollout with gains is the same
update_constraints!(res_reg,solver)
update_constraints!(res_min,solver_min)

update_jacobians!(res_reg,solver)
update_jacobians!(res_min,solver_min)

v_reg = backwardpass!(res_reg,solver)
for k = 1:N
    res_min.K[k][1,:] = res_reg.K[k]
    res_min.K[k][2,:] = zeros(2)
    res_min.d[k][1] = res_reg.d[k][1]
    res_min.d[k][2] = 0.05*0
end

rollout!(res_reg,solver,1.)
rollout!(res_min,solver_min,1.)

solver_min.opts.min_time_regularization = 0
J_reg = cost(solver, res_reg, res_reg.X_, res_reg.U_)
J_min = cost(solver_min, res_min, res_min.X_, res_min.U_)
@test J_reg == J_min
@test to_array(res_reg.X_) ≈ to_array(res_min.X_)

# Test infeasible
solver = Solver(model,obj_c,dt=dt,opts=opts); solver.state.infeasible = true;
solver_min = Solver(model,obj_min,N=21,opts=opts); solver_min.opts.infeasible = true;
m̄,mm = get_num_controls(solver)
@test (m̄,mm) == (1,3)
@test get_num_controls(solver_min) == (2,4)
p,pI,pE = get_num_constraints(solver)
@test get_num_constraints(solver_min) == (7,4,3)

res_reg = ConstrainedVectorResults(n,m+n,p,N)
res_min = ConstrainedVectorResults(n,m+1+n,p+3,solver_min.N)
res_reg.ρ[1] = 0
res_min.ρ[1] = 0

X0 = line_trajectory(solver)
U0 = ones(m,N)
U0_min = [U0; ones(m,solver_min.N)*sqrt(dt)]
ui = infeasible_controls(solver, X0, U0)
ui_min = infeasible_controls(solver_min, X0, U0_min)
@test ui == ui_min
copyto!(res_reg.X, X0)
copyto!(res_min.X, X0)
copyto!(res_reg.U, [ones(m,N); ui])
copyto!(res_min.U, [ones(m,solver_min.N); ones(m,solver_min.N)*sqrt(dt); ui_min])

update_constraints!(res_reg,solver)
update_constraints!(res_min,solver_min)
@test res_reg.C[1][1:2] == res_min.C[1][[1,3]]
@test res_min.C[1][2] == sqrt(dt) - solver_min.opts.max_dt
@test res_reg.C[1][3:4] == res_min.C[1][5:6]

rollout!(res_reg, solver)
rollout!(res_min, solver_min)
@test to_array(res_reg.X) == to_array(res_min.X)
@test to_array(res_min.X) == X0
@test to_array(res_reg.X) == X0

J = cost(solver, res_reg)
J = cost(solver_min, res_min)

update_jacobians!(res_reg,solver)
update_jacobians!(res_min,solver_min)
@test res_reg.fx[1] == res_min.fx[1]
@test res_reg.fu[N-1] == res_min.fu[N-1][:,[1,3,4]]

function init_solve(solver)
    n,m,N = get_sizes(solver)
    m̄,mm = get_num_controls(solver)
    p,pI,pE = get_num_constraints(solver)
    res = ConstrainedVectorResults(n,mm,p,N)
    dt = get_initial_dt(solver)

    U0 = ones(m,N)
    if is_min_time(solver)
        U0 = [U0; ones(m,solver.N)*sqrt(dt)]
    end
    if solver.state.infeasible
        X0 = line_trajectory(solver)
        copyto!(res.X, X0)

        ui = infeasible_controls(solver, X0, U0)
        U0 = [U0; ui]
    end
    copyto!(res.U, U0)

    update_constraints!(res,solver)
    rollout!(res, solver)
    res
end

function inner_loop(res, solver)
    update_jacobians!(res,solver)
    v = backwardpass_mintime!(res,solver)
    with_logger(logger) do
        J = forwardpass!(res,solver,v)
    end

    res.X .= deepcopy(res.X_)
    res.U .= deepcopy(res.U_)

    print_header(logger,TrajectoryOptimization.InnerLoop)
    print_row(logger,TrajectoryOptimization.InnerLoop)

    display(plot(to_array(res_min.U)'))
    total_time(solver_min, res_min)
end

obj_min = update_objective(obj_c,tf=:min,c=10., Q = obj.Q*0, R = obj.R, Qf = obj.Qf*1)
solver_min = Solver(model,obj_min,N=21,opts=opts)
solver_min.opts.infeasible = true;
solver_min.opts.min_time_regularization = 100
solver_min.opts.infeasible_regularization = 10
solver_min.opts.iterations_linesearch = 20
solver_min.opts.max_dt = 0.2
res_min = init_solve(solver_min)
[mu[4] = 100 for mu in res_min.MU]
[mu[5] = 1000 for mu in res_min.MU]
inner_loop(res_min, solver_min)
outer_loop_update_type(res_min,solver_min)

outer_loop_update_type(res_reg,solver)

res_min.U[1]

display(plot(to_array(res_min.X)'))
display(plot(to_array(res_min.U)'))

# Test inner loop
logger = TrajectoryOptimization.SolverLogger(TrajectoryOptimization.InnerIters)
TrajectoryOptimization.add_level!(logger,TrajectoryOptimization.InnerLoop,[],[],print_color=:green)
TrajectoryOptimization.add_level!(logger,TrajectoryOptimization.InnerIters,[],[],print_color=:blue)

obj_min = update_objective(obj_c,tf=:min,c=10., Q = obj.Q*0, R = obj.R, Qf = obj.Qf*0)
opts = SolverOptions()
opts.outer_loop_update_type = :individual
opts.verbose = true
solver = Solver(model,obj_c,dt=dt,opts=opts)
solver_min = Solver(model,obj_min,N=31,opts=opts)
solver_min.opts.min_time_regularization  = 0
n,m,N = get_sizes(solver)
p,pI,pE = TrajectoryOptimization.get_num_constraints(solver)
res_reg = ConstrainedVectorResults(n,m,p,N)
res_min = ConstrainedVectorResults(n,m+1,p+3,solver_min.N)
res_reg.ρ[1] = 0
res_min.ρ[1] = 0

copyto!(res_reg.U, ones(m,N))
copyto!(res_min.U, [ones(m,solver_min.N); ones(m,solver_min.N)*sqrt(dt)])
# copyto!(res_min.U, U_guess)

[mu[end] = 100 for mu in res_min.MU]
update_constraints!(res_reg,solver)
update_constraints!(res_min,solver_min)

rollout!(res_reg, solver)
rollout!(res_min, solver_min)

J = cost(solver, res_reg)
J = cost(solver_min, res_min)

iter = 1

update_jacobians!(res_reg,solver)
update_jacobians!(res_min,solver_min)

# Test that the fwp rollouts are the same given the same bwp gains
v_reg = backwardpass!(res_reg,solver)
v_min = backwardpass_mintime!(res_reg,solver)
v_reg == v_min

v_min = backwardpass_mintime!(res_min,solver_min)
# v_min = backwardpass_mintime2!(res_min,solver_min)

with_logger(logger) do

    J_reg = forwardpass!(res_reg,solver,v_reg)
    J_min = forwardpass!(res_min,solver_min,v_min)

end

print_header(logger,TrajectoryOptimization.InnerLoop)
print_row(logger,TrajectoryOptimization.InnerLoop)

res_reg.X .= deepcopy(res_reg.X_)
res_reg.U .= deepcopy(res_reg.U_)

res_min.X .= deepcopy(res_min.X_)
res_min.U .= deepcopy(res_min.U_)

display(plot(to_array(res_min.U)'))
# display(plot(to_array(res_min.X)'))
total_time(solver_min, res_min)

outer_loop_update_type(res_reg,solver)
outer_loop_update_type(res_min,solver_min)

norm(to_array(res_min.C)[end,:],Inf)
res_min.CN
res_min.IμN
res_min.μN

U_guess = to_array(res_min.U)

res_min.U[2].^2
res_reg.K[1]

res_min.K[1]

##################
#    CARTPOLE    #
##################

# Set up problem
model, obj0 = Dynamics.cartpole_analytical
n,m = model.n, model.m

obj = copy(obj0)
obj.x0 = [0;0;0;0.]
obj.xf = [0.5;pi;0;0]
obj.tf = 2.0
u_bnd = 20
x_bnd = [0.6,Inf,Inf,Inf]
obj_c = ConstrainedObjective(obj,u_min=-u_bnd, u_max=u_bnd)
obj_min = update_objective(obj_c,tf=:min,c=1.,Q = obj.Q*0., Qf = obj.Qf*1)
dt = 0.1

solver = Solver(model,obj_c,dt=dt)
U0 = ones(1,solver.N)
res0,stat0 = solve(solver,U0)
plot(to_array(res0.X)[1:2,:]')
plot(to_array(res0.U)[1:1,:]')

solver_min = Solver(model,obj_min,N=solver.N)
solver_min.opts.use_static = false
solver_min.opts.max_dt = 0.25
solver_min.opts.verbose = true
solver_min.opts.cost_tolerance = 1e-4
solver_min.opts.cost_tolerance_intermediate = 1e-2
solver_min.opts.constraint_tolerance = 0.05
solver_min.opts.outer_loop_update_type = :default
solver_min.opts.min_time_regularization = 100
solver_min.opts.μ1 = 1
U0 = ones(1,solver_min.N)
U0[:,solver.N÷2:end] *= -1
res,stats = solve(solver_min,U0)
plot(to_array(res.X)[1:2,:]')
plot(to_array(res.U)[1:2,:]',ylim=[0,0.5],markershape=:auto)

X0 = line_trajectory(solver_min)
res,stats = solve(solver_min,X0,U0)


text:
function labeldims(n)
    if n == 1
        ("\$t\$", "\$x\$")
    elseif n == 2
        ("\$x\$", "\$y\$")
    else
        ("\$x\$", "\$y\$", "\$z\$")
    end
end

function stationarycharacteristics(P::Process; n = length(getX0(P)), labels=labeldims(n), xlabel=labels[1], ylabel=labels[2], zlabel=(length(labels) == 3 ? labels[3] : nothing), size=(400, 350))
    p1 = marginaltrajectory(P; xlabel, ylabel, zlabel, size, mainalpha=0.35)
    display(p1)
    AMI = selfmutualinfo(P)
    p2 = plot(AMI.dims[1].val, Array((AMI isa AbstractVector ? AMI : AMI[:, end:-1:1]));  legend=true,   foreground_color_legend = nothing,
    background_color_legend = nothing, legendtitle=nothing, size,
    #xticks=[0, optimize_ticks(extrema(dims(AMI, 1).val)...; k_min = 3, k_max = 5)[1]...],
    yticks=optimize_ticks(extrema(AMI)...; k_min = 3, k_max = 6)[1],
    xlim=(0, Inf), xlabel="\$\\tau\$", ylabel="AMI")
    display(p2)
    return nothing
end
export stationarycharacteristics

function plotlyapunovresponse(P::Process, p, prange; size=(400, 350), N=1000, k=length(getX0(P)), kwargs...)
    λs = lyapunovresponse(P, p, prange, N, k)
    plot(prange, λs; size, xlabel="\$\\theta\$", ylabel="\$\\lambda\$", kwargs...)
end
export plotlyapunovresponse


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -19.57782063716216 <= q <= 17.16108472444359)

                     add_NL_constraint(m, :(softplus(-0.6503895050792288 + 0.5000448231799557*softplus(0.48086890691524875 + 0.8604652607112211*$(x[1]) + -0.9524910094017911*$(x[2]) + -0.8411501540974782*$(x[3])) + 0.8223864208279696*softplus(0.4084009324685236 + 0.9106618206937234*$(x[1]) + -0.8345217492048942*$(x[2]) + 0.6009855404051798*$(x[3])) + -0.39198040344556784*softplus(0.29882092390692216 + -0.9088484335369627*$(x[1]) + -0.948042223821898*$(x[2]) + 0.5765779510962643*$(x[3])) + 0.4922116220556738*softplus(0.9223388106039003 + 0.14028925335959963*$(x[1]) + -0.026802684789215103*$(x[2]) + 0.9753923898070198*$(x[3]))) + softplus(0.1539410450835783 + -0.07568393712605781*softplus(0.48086890691524875 + 0.8604652607112211*$(x[1]) + -0.9524910094017911*$(x[2]) + -0.8411501540974782*$(x[3])) + -0.902375002810281*softplus(0.4084009324685236 + 0.9106618206937234*$(x[1]) + -0.8345217492048942*$(x[2]) + 0.6009855404051798*$(x[3])) + 0.7665443117707627*softplus(0.29882092390692216 + -0.9088484335369627*$(x[1]) + -0.948042223821898*$(x[2]) + 0.5765779510962643*$(x[3])) + -0.49047645888895497*softplus(0.9223388106039003 + 0.14028925335959963*$(x[1]) + -0.026802684789215103*$(x[2]) + 0.9753923898070198*$(x[3]))) + softplus(0.6062617117178619 + -0.6503507748090844*softplus(0.48086890691524875 + 0.8604652607112211*$(x[1]) + -0.9524910094017911*$(x[2]) + -0.8411501540974782*$(x[3])) + -0.9378359170676176*softplus(0.4084009324685236 + 0.9106618206937234*$(x[1]) + -0.8345217492048942*$(x[2]) + 0.6009855404051798*$(x[3])) + -0.23896524724444168*softplus(0.29882092390692216 + -0.9088484335369627*$(x[1]) + -0.948042223821898*$(x[2]) + 0.5765779510962643*$(x[3])) + -0.3347498056547815*softplus(0.9223388106039003 + 0.14028925335959963*$(x[1]) + -0.026802684789215103*$(x[2]) + 0.9753923898070198*$(x[3]))) + softplus(-0.5413645906370745 + 0.4055343512605676*softplus(0.48086890691524875 + 0.8604652607112211*$(x[1]) + -0.9524910094017911*$(x[2]) + -0.8411501540974782*$(x[3])) + -0.13388786735965885*softplus(0.4084009324685236 + 0.9106618206937234*$(x[1]) + -0.8345217492048942*$(x[2]) + 0.6009855404051798*$(x[3])) + -0.9731450409901217*softplus(0.29882092390692216 + -0.9088484335369627*$(x[1]) + -0.948042223821898*$(x[2]) + 0.5765779510962643*$(x[3])) + 0.17849747613123057*softplus(0.9223388106039003 + 0.14028925335959963*$(x[1]) + -0.026802684789215103*$(x[2]) + 0.9753923898070198*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
using JLD2
using FileIO
using WindMixing
using OceanParameterizations
using Plots

train_files = ["wind_-5e-4_diurnal_3e-8"]

𝒟 = WindMixing.data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)

𝒟.wT.coarse

plot(1:length(𝒟.wT.coarse[end,:]), 𝒟.wT.coarse[end,:])

file = jldopen("Data/three_layer_constant_fluxes_linear_hr192_Qu0.0e+00_Qb3.0e-08_f1.0e-04_Nh256_Nz128_free_convection_8days_statistics.jld2")
# file = jldopen("Data/three_layer_constant_fluxes_linear_hr192_Qu1.0e-03_Qb0.0e+00_f1.0e-04_Nh256_Nz128_wind_mixing_8days_2_statistics.jld2")

keys(file["parameters/boundary_condition_θ_top"])
file["parameters/boundary_condition_θ_top"]
wb = file["parameters/buoyancy_flux"]
α = file["parameters/thermal_expansion_coefficient"]
g = file["parameters/gravitational_acceleration"]

wb / α / g ≈ file["parameters/boundary_condition_θ_top"]

@inline Qᶿ(x, y, t) = Qᵇ * sin(2π / (Qᵇ_period * 60 ^ 2) * t) / (α * g)
file["parameters/buoyancy_flux"]
file["parameters/momentum_flux"]

keys(file["timeseries/T"])

file["timeseries/T/$(keys(file["timeseries/T"])[5])"][1,1,:]
file["timeseries/wT/$(keys(file["timeseries/T"])[2])"][1,1,end]
close(file)

newfile = jldopen("D:\\Downloads\\test_functions.jld2", "w") do file
    file["function"] = Qᶿ
end

Qᵇ = 3
new_file = jldopen("D:\\Downloads\\test_functions.jld2")
new_file["function"](0,0,1)

train_files = ["diurnal_test"]
𝒟 = WindMixing.data(train_files, scale_type=ZeroMeanUnitVarianceScaling, enforce_surface_fluxes=false)


text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

using Printf

"""
    printmatrix(mat [, Val(F), colnames])

Display a matrix in full while controlling the formatting of each value and
optionally showing the column names.
  * `Val(F)` - display each number in the given format `F`. The format is in the
    form of a decimal point number where the whole part indicates the total
    width and the fractional part is the number of digits printed after the
    decimal point. Default is `Val(12.7)`
  * `colnames` - a list of names to display in the first row. The names are
    displayed as given, possibly with padding to match the width given in the
    `Val` argument. If any names are longer than that, they will not be
    truncated and so the display will not be aligned properly. Sorry about that!

"""
@inline printmatrix(mat::AbstractMatrix, args...) = printmatrix(stdout, mat, args...)
@inline printmatrix(io::IO, mat::AbstractMatrix, args...) = printmatrix(io, mat, Val(12.7), args...)
@generated function printmatrix(io::IO, mat::AbstractMatrix, ::Val{N}, cols = nothing) where N
    fmts = "% $(N)s "
    fmtn = "% $(N)f "
    return quote
        m, n = size(mat)
        if cols !== nothing
            s = ""
            for j ∈ 1:n
                s *= Printf.@sprintf($fmts, cols[j])
            end
            println(io, s)
        end
        for i in 1:m
            s = ""
            for j in 1:n
                s *= Printf.@sprintf($fmtn, mat[i,j])
            end
            println(io, s)
        end
        return nothing
    end
end
export printmatrix

import ModelBaseEcon: transform, inverse_transform

function transform(data::AbstractMatrix{Float64}, m::Model)
    tdata = similar(data)
    for (i, v) in enumerate(m.varshks)
        tdata[:, i] .= transform(data[:, i], v)
    end
    return tdata
end

function inverse_transform(data::AbstractMatrix{Float64}, m::Model)
    idata = similar(data)
    for (i, v) in enumerate(m.varshks)
        idata[:, i] .= inverse_transform(data[:, i], v)
    end
    return idata
end


text:
using DataFrames: DataFrames, DataFrame
import Dates
using TimeZones: TimeZones, ZonedDateTime, @tz_str
using TypedTables: TypedTables, Table
import CSV

@system StoreBase begin
    filename => "" ~ preserve::str(parameter)
    ik: indexkey => :index ~ preserve::sym(parameter)

    i(t=context.clock.tick): index => t + 1 ~ track::int
    ix: indexer ~ hold

    s: store ~ hold
end

@system DataFrameStore(StoreBase) begin
    ix(; r::Cropbox.DataFrames.DataFrameRow): indexer => Cropbox.DataFrames.row(r) ~ call::int

    df(filename, ik, ix): dataframe => begin
        df = Cropbox.CSV.File(filename) |> Cropbox.DataFrames.DataFrame |> unitfy
        df[!, ik] = map(ix, eachrow(df))
        df
    end ~ preserve::Cropbox.DataFrames.DataFrame(extern, parameter)

    gdf(df, ik): grouped_dataframe => begin
        Cropbox.DataFrames.groupby(df, ik)
    end ~ preserve::Cropbox.DataFrames.GroupedDataFrame{Cropbox.DataFrame}

    s(gdf, i): store => begin
        gdf[(i,)][1, :]
    end ~ track::Cropbox.DataFrames.DataFrameRow{Cropbox.DataFrame,Cropbox.DataFrames.Index}
end

@system DayStore(DataFrameStore) begin
    i(context.clock.time): index ~ track::int(u"d")

    daykey => :day ~ preserve::sym(parameter)
    ix(daykey; r::Cropbox.DataFrames.DataFrameRow): indexer => r[daykey] ~ call::int(u"d")
end

@system DateStore(DataFrameStore) begin
    calendar(context) ~ ::Calendar
    i(t=calendar.time): index => Cropbox.Dates.Date(t) ~ track::date

    datekey => :date ~ preserve::sym(parameter)
    ix(datekey; r::Cropbox.DataFrames.DataFrameRow): indexer => r[datekey] ~ call::date
end

@system TimeStore(DataFrameStore) begin
    calendar(context) ~ ::Calendar
    i(calendar.time): index ~ track::datetime

    datekey => :date ~ preserve::sym(parameter)
    timekey => :time ~ preserve::sym(parameter)
    tz: timezone => Cropbox.tz"UTC" ~ preserve::Cropbox.TimeZones.TimeZone(parameter)
    ix(datekey, timekey, tz; r::Cropbox.DataFrames.DataFrameRow): indexer => begin
        #HACK: handle ambiguous time conversion under DST
        occurrence = 1
        i = Cropbox.DataFrames.row(r)
        if i > 1
            r0 = parent(r)[i-1, :]
            r0[timekey] == r[timekey] && (occurrence = 2)
        end
        dt = Cropbox.Dates.DateTime(r[datekey], r[timekey])
        Cropbox.ZonedDateTime(dt, tz, occurrence)
    end ~ call::datetime
end

@system TableStore(StoreBase) begin
    #TODO: avoid dynamic dispatch on Table/NamedTuple
    ix(; i::Int, r::NamedTuple): indexer => i ~ call::int
    tb(filename, ik, ix): table => begin
        tb = Cropbox.CSV.File(filename) |> Cropbox.TypedTables.FlexTable
        setproperty!(tb, ik, map(enumerate(eachrow(tb))) do (i, r)
            ix(i, r)
        end)
        Cropbox.Table(tb)
    end ~ preserve::Cropbox.Table(extern, parameter)
    s(tb, i): store => tb[i] ~ track::NamedTuple
end

export DataFrameStore, DayStore, DateStore, TimeStore, TableStore


text:
# https://en.wikipedia.org/wiki/Day_count_convention

export yearlength, daycount, yearfraction

export Thirty360Bond, Thirty360US, ThirtyE360, ThirtyE360ISDA
export Actual365Fixed, Actual364, Actual360
export ActualActualISDA

abstract type DayCountConvention end

abstract type FixedYear{N} <: DayCountConvention end
abstract type ActualYear <: DayCountConvention end

abstract type Thirty360 <: FixedYear{360} end

struct Thirty360Bond <: Thirty360 end
struct Thirty360US <: Thirty360
    eom::Bool
end
struct ThirtyE360 <: Thirty360 end
struct ThirtyE360ISDA <: Thirty360
    maturity::Date
end

abstract type ActualFixedYear{N} <: FixedYear{N} end

struct Actual365Fixed <: ActualFixedYear{365} end
struct Actual364 <: ActualFixedYear{364} end
struct Actual360 <: ActualFixedYear{360} end

abstract type ActualActual <: ActualYear end

struct ActualActualISDA <: ActualActual end

yearlength(::FixedYear{N}, year = 0) where {N} = N

yearlength(::ActualYear, year) = isleapyear(year) ? 366 : 365

daycount(::Union{ActualFixedYear, ActualActual}, a, b) = Dates.value(b - a)

daycount(c::Thirty360, a, b) = begin
    y1, m1 = yearmonth(a)
    y2, m2 = yearmonth(b)

    d1, d2 = adjustdays(c, a, b)

    360 * (y2 - y1) + 30 * (m2 - m1) + (d2 - d1)
end

adjustdays(::Thirty360Bond, a, b) = begin
    d1, d2 = day(a), day(b)

    d1 = min(d1, 30)
    d1 == 30 && (d2 = min(d2, 30))

    d1, d2
end

adjustdays(c::Thirty360US, a, b) = begin
    m1, d1 = monthday(a)
    m2, d2 = monthday(b)

    if c.eom && m1 == 2 && d1 == day(lastdayofmonth(a))
        m2 == 2 && d2 == day(lastdayofmonth(b)) && (d2 = 30)
        d1 = 30
    end

    d1 = min(d1, 30)
    d1 == 30 && (d2 = min(d2, 30))

    d1, d2
end

adjustdays(::ThirtyE360, a, b) = begin
    d1, d2 = day(a), day(b)

    d1 = min(d1, 30)
    d2 = min(d2, 30)

    d1, d2
end

adjustdays(c::ThirtyE360ISDA, a, b) = begin
    y1, m1, d1 = yearmonthday(a)
    y2, m2, d2 = yearmonthday(b)

    d1 == day(lastdayofmonth(a)) && (d1 = 30)
    d2 == day(lastdayofmonth(b)) && (m2 != 2 || b != c.maturity) && (d2 = 30)

    d1, d2
end

yearfraction(c::Union{Thirty360, ActualFixedYear}, a, b, s::Type{S} = Float64) where {S<:Real} = convert(S, daycount(c, a, b)) / yearlength(c)

yearfraction(c::ActualActualISDA, a, b, ::Type{S} = Float64) where {S<:Real} = begin
    firstyear = 1 - yearfraction(c, a)
    betweenyears = year(b) - year(a) - 1
    lastyear = yearfraction(c, b)

    convert(S, firstyear + betweenyears + lastyear)
end

yearfraction(c::ActualActualISDA, date, ::Type{S} = Float64) where {S<:Real} = (dayofyear(date) - 1) // yearlength(c, date)


text:
import PyPlot
const plt = PyPlot

function M(p::Int64,u)
    if p == 2
        if u >= 0 && u <= 2
            return 1 - abs(u-1)
        else
            return 0.0
        end
    else
        return M(p-1,u)*u/(p-1) + (p-u)/(p-1)*M(p-1,u-1)
    end
end

function main()
    @time println( M(2,1.1) )
    @time println( M(3,1.1) )
    @time println( M(4,1.1) )

    NptsPlot = 100
    x = Array{Float64}( linspace(0.0,5.0,NptsPlot) )
    y2 = Array{Float64}(NptsPlot)
    y3 = Array{Float64}(NptsPlot)
    y4 = Array{Float64}(NptsPlot)
    y5 = Array{Float64}(NptsPlot)
    for i = 1:NptsPlot
        y2[i] = M(2,x[i])
        y3[i] = M(3,x[i])
        y4[i] = M(4,x[i])
        y5[i] = M(5,x[i])
    end
    plt.clf()
    plt.plot( x, y2, linewidth=2.0, label="M2.png" )
    plt.plot( x, y3, linewidth=2.0, label="M3.png" )
    plt.plot( x, y4, linewidth=2.0, label="M4.png" )
    plt.plot( x, y5, linewidth=2.0, label="M5.png" )
    plt.grid(true)
    plt.legend()
    plt.savefig("M2.png", dpi=300)
end

main()


text:
export BicategoryRelations, FreeBicategoryRelations,
  AbelianBicategoryRelations, FreeAbelianBicategoryRelations,
  meet, join, top, bottom, plus, zero, coplus, cozero

import Base: join, zero

# Bicategory of relations
#########################

""" Theory of *bicategories of relations*

TODO: The 2-morphisms are missing. I haven't decided how to handle them yet.

References:

- Carboni & Walters, 1987, "Cartesian bicategories I"
- Walters, 2009, blog post, "Categorical algebras of relations",
  http://rfcwalters.blogspot.com/2009/10/categorical-algebras-of-relations.html
"""
@signature MonoidalCategoryWithBidiagonals(Ob,Hom) => BicategoryRelations(Ob,Hom) begin
  # Dagger category.
  dagger(f::(A → B))::(B → A) ⊣ (A::Ob,B::Ob)

  # Self-dual compact closed category.
  dunit(A::Ob)::(munit() → (A ⊗ A))
  dcounit(A::Ob)::((A ⊗ A) → munit())

  # Logical operations.
  meet(f::(A → B), g::(A → B))::(A → B) ⊣ (A::Ob, B::Ob)
  top(A::Ob, B::Ob)::(A → B)
end

@syntax FreeBicategoryRelations(ObExpr,HomExpr) BicategoryRelations begin
  otimes(A::Ob, B::Ob) = associate_unit(new(A,B), munit)
  otimes(f::Hom, g::Hom) = associate(new(f,g))
  compose(f::Hom, g::Hom) = associate(new(f,g; strict=true))
  dagger(f::Hom) = distribute_unary(distribute_dagger(involute(new(f))),
                                    dagger, otimes)
  meet(f::Hom, g::Hom) = compose(mcopy(dom(f)), otimes(f,g), mmerge(codom(f)))
  top(A::Ob, B::Ob) = compose(delete(A), create(B))
end

""" Theory of *abelian bicategories of relations*

Unlike Carboni & Walters, we use additive notation and nomenclature.

References:

- Carboni & Walters, 1987, "Cartesian bicategories I", Sec. 5
- Baez & Erbele, 2015, "Categories in control"
"""
@signature BicategoryRelations(Ob,Hom) => AbelianBicategoryRelations(Ob,Hom) begin
  # Second diagonal and codiagonal.
  plus(A::Ob)::((A ⊗ A) → A)
  zero(A::Ob)::(munit() → A)
  coplus(A::Ob)::(A → (A ⊗ A))
  cozero(A::Ob)::(A → munit())

  # Logical operations.
  join(f::(A → B), g::(A → B))::(A → B) ⊣ (A::Ob, B::Ob)
  bottom(A::Ob, B::Ob)::(A → B)
end

@syntax FreeAbelianBicategoryRelations(ObExpr,HomExpr) AbelianBicategoryRelations begin
  otimes(A::Ob, B::Ob) = associate_unit(new(A,B), munit)
  otimes(f::Hom, g::Hom) = associate(new(f,g))
  compose(f::Hom, g::Hom) = associate(new(f,g; strict=true))
  dagger(f::Hom) = distribute_unary(distribute_dagger(involute(new(f))),
                                    dagger, otimes)
  meet(f::Hom, g::Hom) = compose(mcopy(dom(f)), otimes(f,g), mmerge(codom(f)))
  join(f::Hom, g::Hom) = compose(coplus(dom(f)), otimes(f,g), plus(codom(f)))
  top(A::Ob, B::Ob) = compose(delete(A), create(B))
  bottom(A::Ob, B::Ob) = compose(cozero(A), zero(B))
end

function show_latex(io::IO, expr::HomExpr{:plus}; kw...)
  Syntax.show_latex_script(io, expr, "\\blacktriangledown")
end
function show_latex(io::IO, expr::HomExpr{:coplus}; kw...)
  Syntax.show_latex_script(io, expr, "\\blacktriangle")
end
function show_latex(io::IO, expr::HomExpr{:zero}; kw...)
  Syntax.show_latex_script(io, expr, "\\blacksquare")
end
function show_latex(io::IO, expr::HomExpr{:cozero}; kw...)
  Syntax.show_latex_script(io, expr, "\\blacklozenge")
end


text:
function fT(x,Param)
str = lower(Param.ProfTheta)
if str == "baldaufsphere"
    (lon,lat,r)=cart2sphere(x[1],x[2],x[3]);
    r=r-Param.RadEarth;
    p=Param.p0*exp(-Param.Grav*r/(Param.Rd*Param.T0));
    Rho=p/(Param.Rd*Param.T0);
    T=p/(Rho*Param.Rd);
    d=acos(sin(Param.lat0)*sin(lat)+cos(Param.lat0)*cos(lat)*cos(lon-Param.lon0));
    T=T+Param.DeltaT*exp(-Param.ExpDist*d)*sin(pi*r/Param.H);
elseif str == "baldaufcart"
    delta=Param.Grav/(Param.Rd*Param.T0);

    dT=Param.DeltaT*exp(-(x[1]-Param.xc)^2/Param.d^2)*sin(pi*x[3]/Param.H);

    T=Param.T0+exp(delta/2*x[3])*dT;
end
return T
end


text:
### A Pluto.jl notebook ###
# v0.12.15

using Markdown
using InteractiveUtils

# ╔═╡ 2a872aae-33e9-11eb-09b2-5b2f2f54d419
md"""
## Individual Economics Formalism
"""

# ╔═╡ 847be240-33e9-11eb-3a0c-05989d6aa48e
md"""
This is how text rendering works.
"""

# ╔═╡ a470e230-33e9-11eb-1184-11506d5c302e
md"""This is how math rendering works.
	\begin{equation}
	i\hbar\frac{\partial \psi}{\partial t} = \hat H \psi
	\end{equation}
	"""

# ╔═╡ c40b3f00-33e9-11eb-3916-dff368cf3fbc
md"""
Let's take some notes.

>1. Initial conditions: Population, resources.
>2. Evolutionary traits: Population, resources, trade, demand and supply, currency.
"""

# ╔═╡ da9645d0-33e9-11eb-3e56-b58c1973201a
md"""
Sample codeblock for $Demand :: Person \times Person \to \mathbb R$.
```
demand(alice, bob) = ???
```
"""

# ╔═╡ 98695ad0-33e9-11eb-2dec-65f3e590b654
hint(text) = Markdown.MD(Markdown.Admonition("hint", "Hint", [text]));

# ╔═╡ 9419fde0-33e9-11eb-0127-c9fd0b722bb5
hint(md"""
	Click on the eye on the left to show the "code" that renders this box.
	""")

# ╔═╡ Cell order:
# ╟─2a872aae-33e9-11eb-09b2-5b2f2f54d419
# ╠═847be240-33e9-11eb-3a0c-05989d6aa48e
# ╠═a470e230-33e9-11eb-1184-11506d5c302e
# ╠═c40b3f00-33e9-11eb-3916-dff368cf3fbc
# ╠═da9645d0-33e9-11eb-3e56-b58c1973201a
# ╟─9419fde0-33e9-11eb-0127-c9fd0b722bb5
# ╟─98695ad0-33e9-11eb-2dec-65f3e590b654

max_stars_repo_name: GodotMisogi/IndEco

max_issues_repo_name: GodotMisogi/IndEco

max_forks_repo_name: GodotMisogi/IndEco


text:
G     = 6.67384e-11     # m/(kg*s^2)
dt    = 60*60*24*365.25 # Years in seconds
r_ly  = 9.4607e15       # Lightyear in m
m_sol = 1.9891e30       # Solar mass in kg
b     = 0.0

n = 10

box_size = 1
grid_spacing = 1.0 * box_size / n

iterations = 10

v_x = zeros((n,n))
v_y = zeros((n,n))
eta = ones((n,n)) # pressure
for i=1:n
	eta[i] = eta[i] * 0.1 * i
end

#temps = [ones((n,n)),
#		 ones((n,n)),
#		 ones((n,n)),
#		 ones((n,n))]

tmp = ones((n,n))
du_dt = ones((n,n))
dv_dt = ones((n,n))
tmp1 = ones((n,n))

for i in 1:iterations
	# Compute derivatives with respect to x and y by using difference in z (by shifting) / difference in x or y (going both sides)
	# Derivative of eta with respect to x (second dimension?????!???!!!??)
	roll1 = circshift(eta, (0,-1))
	roll2 = circshift(eta, (0,1))
	tmp1 = (roll1 .- roll2) ./ (grid_spacing * 2.0)

	# Derivative of eta with respect to y (first dimension?????!???!!!??)
	roll1 = circshift(eta, (-1,0))
	roll2 = circshift(eta, (1,0))
	tmp = (roll1 .- roll2) ./ (grid_spacing * 2.0)

	du_dt = -g .* tmp1 - b * u
	dv_dt = -g .* tmp - b .* v

	H = 0
	tmp1 = ((eta .+ H) .* u)
	tmp2 = tmp 1 .* v

	roll1 = circshift(tmp1, (0,-1))
	roll2 = circshift(tmp1, (0,1))
	tmp1 = -1 * (roll1 .- roll2) ./ (grid_spacing * 2.0)

	roll1 = circshift(tmp, (0,-1))
	roll2 = circshift(tmp, (0,1))
	tmp = tmp1 .- ((roll1 .- roll2) ./ (grid_spacing * 2.0))

	# On to line 106
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
# Pkg.build("QuantEcon")
import Dolo
path = path = Dolo.pkg_path

using AxisArrays

filename = joinpath(path,"examples","models","rbc_dtcc_iid_ar1.yaml")
model = Dolo.yaml_import(filename)
@time dr = Dolo.time_iteration(model, verbose=true, maxit=10000)

driving_process = Dolo.simulate(model2.exogenous, 10, 40, m0)
sim1 = Dolo.simulate(model, dr, driving_process)

sim = Dolo.simulate(model, dr; N= 1000, T=50)
# Ac= cat(model.symbols[:exogenous], model.symbols[:states], model.symbols[:controls]; dims=1)
# ll=[Symbol(i) for i in Ac]
# AA= AxisArray(sim, Axis{:N}(1:1000), Axis{:V}(ll), Axis{:T}(1:50))

Tvalues = linspace(1, sim[Axis{:T}][end], sim[Axis{:T}][end])

a_sort = sort(collect(sim[Axis{:V}(:k)]),1)
a_median = a_sort[Int(size(sim,1)*0.5),:]

import PyPlot
plt = PyPlot;

for i in 1:sim[Axis{:N}][end]
    plt.plot(Tvalues, sim[Axis{:N}(i), Axis{:V}(:k)], color="red", alpha = .05)
end
plt.title("Simulations");
plt.xlabel("Horizon");
plt.ylabel("Capital");

plt.plot(Tvalues, a_median, marker="o")

#########################################################################

##################  IRF         ###########################

s0 = model.calibration[:states]
e0 = model.calibration[:exogenous]
irf = Dolo.response(model, dr, s0, :e_z)
Dolo.response(model, dr, s0, [0.3])

############## 2 shocks
filename2 = joinpath(path,"examples","models","rbc_dtcc_iid_2ar1.yaml")
model2 = Dolo.yaml_import(filename2)
@time dr2 = Dolo.time_iteration(model2, verbose=true, maxit=10000)
s0 = model2.calibration[:states]
e0 = model2.calibration[:exogenous]

irf=Dolo.response(model2, dr2, s0, [0, 0.3]; T=40)

irf[:k]
irf2=Dolo.response(model2, dr2, s0, :e_z)

Dolo.response(model2, dr2, s0, :e_d)

#########################################################################

# filename = joinpath(path,"examples","models","rbc.yaml")
# model = Dolo.yaml_import(filename)
# # model2 = Dolo.yaml_import(filename2)
# N = 1
# T=40
# @time dr = Dolo.time_iteration(model, verbose=true, maxit=10000)
# s0 = model.calibration[:states]
# m0 = model.calibration[:exogenous]
# Dolo.response(model, dr, s0, [0.2])
# Dolo.response(model, dr, s0, [0.4])

############
# Ploting

# irf = Dolo.response(model, dr, e0; horizon = 100)

index = findfirst(isequal(:k), model.symbols[:states])

kirf = irf[:,3]
iirf = irf[:,2]
nirf = irf[:,4]
zirf = irf[:,5]
horizon=100
time = linspace(0,horizon-1,horizon)
using Gadfly

plot(x=time, y=kirf, Geom.point, Geom.line,
     Guide.xlabel("horizon"), Guide.ylabel("Capital"), Guide.title("IRF"))
plot(x=time, y=nirf, Geom.point, Geom.line,Guide.xlabel("horizon"),
     Guide.ylabel("Hours"), Guide.title("IRF"))
plot(x=time, y=iirf, Geom.point, Geom.line, Guide.xlabel("horizon"),
   Guide.ylabel("Investments"), Guide.title("IRF"))
plot(x=time, y=zirf, Geom.point, Geom.line, Guide.xlabel("horizon"),
      Guide.ylabel("AR1"), Guide.title("IRF"))

# test simulation
res = Dolo.simulation(model, sigma, dr,s0, n_exp, horizon, seed, zeros(0, 0))
res_long = Dolo.simulation(model, sigma, dr,s0, n_exp=0, horizon=1000, seed=42)

res = Dolo.simulation(model, sigma, dr,s0)
res = Dolo.simulation(model, sigma)

kvec = res_long[:,2,:]
ivec = res_long[:,4,:]
nvec = res_long[:,3,:]
horizon=1000
time = linspace(0,horizon-1,horizon)
using Gadfly

plot(x=time, y=kvec, Geom.point, Geom.line,
     Guide.xlabel("horizon"), Guide.ylabel("Capital"), Guide.title("Simulations"))
plot(x=time, y=nvec, Geom.point, Geom.line,Guide.xlabel("horizon"), Guide.ylabel("Hours"), Guide.title("Simulations"))
plot(x=time, y=ivec, Geom.point, Geom.line, Guide.xlabel("horizon"), Guide.ylabel("Investments"), Guide.title("Simulations"))

O = Theme(
    default_color=colorant"orange")

plot(layer(x=time, y=nvec, Geom.point, Geom.line, O),
     layer(x=time, y=ivec, Geom.point, Geom.line),
     Guide.xlabel("horizon"),
     Guide.ylabel("Hours (orange), Investment"),
     Guide.title("Simulations"))

# verbose=true
# verbose && @printf "%-8s%-10s%-10s%-10s%-5s\n" "t" model.symbols[:states][1] model.symbols[:states][2] model.symbols[:controls][1] model.symbols[:controls][2]
# verbose && println(repeat("-", 35))

#verbose && @printf "%-8s%-10s%-10s%-10s%-5s\n"  t round(s[1],2) round(s[2],2) round(x[1],2) round(x[2],2)


text:
import RvSpectMLBase
using Test

@testset "Theoretical Instrument" begin
    using RvSpectMLBase.TheoreticalInstrument
    num_pixels = 100
    num_orders = 10
    resolution = 1e5

    n1 = TheoreticalInstrument1D()
    lambda2d = collect(reshape(range(5000.0,stop=6000.0,length=num_pixels*num_orders),num_pixels,num_orders))
    @test_nowarn TheoreticalInstrument.TheoreticalInstrument2D(lambda2d)
    n2 = TheoreticalInstrument.TheoreticalInstrument2D(lambda2d)
    @testset "Calc λs" begin
            @test_nowarn TheoreticalInstrument.calc_λs(n1)
            @test_nowarn TheoreticalInstrument.calc_λs(n2)
    end
end

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
module LightGBM

using Dates

import Base
import Libdl
import StatsBase

# we build it so we can assert it is present...
const LGBM_library = abspath(joinpath(@__DIR__, "lib_lightgbm.$(Libdl.dlext)"))

include("wrapper.jl")
include("estimators.jl")
include("utils.jl")
include("fit.jl")
include("predict.jl")
include("cv.jl")
include("search_cv.jl")
include(joinpath(@__DIR__, "MLJInterface.jl"))

export fit!, predict, predict_classes, cv, search_cv, savemodel, loadmodel!
export LGBMEstimator, LGBMRegression, LGBMClassification

end # module LightGBM


text:
using Base.Math: @horner

#Using piecewise approximation polynomial as given in
#'Fast Computation of Complete Elliptic Integrals and Jacobian Elliptic Functions'
#                          Fukushima, Toshio. (2014). F09-FastEI. Celest Mech Dyn Astr
#                                                            DOI 10.1007/s10569-009-9228-z
#Link : https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf
#
#For m<0 , followed  Fukushima, Toshio. (2014). Precise, compact, and fast computation of complete elliptic integrals by piecewise minimax rational function approximation. Journal of Computational and Applied Mathematics. 282. 10.13140/2.1.1946.6245.
# Link:  https://www.researchgate.net/profile/Toshio_Fukushima/publication/267330394_Precise_compact_and_fast_computation_of_complete_elliptic_integrals_by_piecewise_minimax_rational_function_approximation/links/544b81a40cf2d6347f43074f/Precise-compact-and-fast-computation-of-complete-elliptic-integrals-by-piecewise-minimax-rational-function-approximation.pdf?origin=publication_detail
#Also suggested in this paper that we should consider domain only from (-inf,1].

"""
   ellipk(x)

 DLMF : https://dlmf.nist.gov/19.2#E4  , https://dlmf.nist.gov/19.2#E8
 Wiki : https://en.wikipedia.org/wiki/Elliptic_integral
Computes Complete Elliptic Integral of 1st kind at `x`-> K(x)--- given by: ``\\int_{0}^{\\pi/2} \\frac{1}{\\sqrt{1 - x(\\sin \\theta )^{2}}} d\\theta``
"""
function ellipk(a::Float64)
    flag = false
    if a<0.0
        x=a/(a-1)   #dealing with negative args
        flag=true
    elseif a>=0.0
        x=a
        flag=false
    end

    if x==0.0
        return pi/2

    elseif x==1.0
        return Inf

    elseif x>1.0
        throw(DomainError(x,"`x` must lie between -Inf and 1 ---- Domain: (-Inf,1.0]"))

    elseif 0.0 <= x < 0.1   #Table 2 from paper
        t= x-0.05
        t= @horner(t, 1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 , 0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711 , 0.108938811574293531 , 0.098853409871592910 , 0.091439629201749751 , 0.085842591595413900 , 0.081541118718303215)

    elseif 0.1 <= x <0.2  #Table 3
        t=x-0.15
        t= @horner(t , 1.635256732264579992 , 0.471190626148732291 , 0.309728410831499587 , 0.252208311773135699 , 0.226725623219684650 , 0.215774446729585976 , 0.213108771877348910 , 0.216029124605188282 , 0.223255831633057896 , 0.234180501294209925 , 0.248557682972264071 , 0.266363809892617521)

    elseif 0.2 <= x < 0.3 #Table 4
        t=x-0.25
        t= @horner(t , 1.685750354812596043 , 0.541731848613280329 , 0.401524438390690257 , 0.369642473420889090 , 0.376060715354583645 , 0.405235887085125919 , 0.453294381753999079 , 0.520518947651184205 , 0.609426039204995055 , 0.724263522282908870 , 0.871013847709812357 , 1.057652872753547036)

    elseif 0.3 <= x < 0.4 #Table 5
        t= x-0.35
        t= @horner(t , 1.744350597225613243 , 0.634864275371935304 , 0.539842564164445538 , 0.571892705193787391 , 0.670295136265406100 , 0.832586590010977199 , 1.073857448247933265 , 1.422091460675497751 , 1.920387183402304829 , 2.632552548331654201 , 3.652109747319039160 , 5.115867135558865806 , 7.224080007363877411)

    elseif 0.4 <= x < 0.5 #Table 6
        t=x-0.45
        t= @horner(t, 1.813883936816982644 , 0.763163245700557246 , 0.761928605321595831 , 0.951074653668427927 , 1.315180671703161215 , 1.928560693477410941 , 2.937509342531378755 , 4.594894405442878062 , 7.330071221881720772 , 11.87151259742530180 , 19.45851374822937738 , 32.20638657246426863 , 53.73749198700554656 , 90.27388602940998849)

    elseif 0.5 <= x < 0.6 #Table 7
        t= x-0.55
        t= @horner(t , 1.898924910271553526 , 0.950521794618244435 , 1.151077589959015808 , 1.750239106986300540 , 2.952676812636875180 , 5.285800396121450889 , 9.832485716659979747 , 18.78714868327559562 , 36.61468615273698145 , 72.45292395127771801 , 145.1079577347069102 , 293.4786396308497026 , 598.3851815055010179 , 1228.420013075863451 , 2536.529755382764488)
    elseif 0.6 <= x < 0.7 #Table 8
        t=x-0.65
        t =@horner(t , 2.007598398424376302 , 1.248457231212347337 , 1.926234657076479729 , 3.751289640087587680 , 8.119944554932045802 , 18.66572130873555361 , 44.60392484291437063 , 109.5092054309498377 , 274.2779548232413480,697.5598008606326163 , 1795.716014500247129 , 4668.381716790389910 , 12235.76246813664335 , 32290.17809718320818 , 85713.07608195964685 , 228672.1890493117096 , 612757.2711915852774)

    elseif 0.7 <= x <0.8 #Table 9
        t=x-0.75
        t= @horner(t, 2.156515647499643235 , 1.791805641849463243 , 3.826751287465713147 , 10.38672468363797208 , 31.40331405468070290 , 100.9237039498695416 , 337.3268282632272897 , 1158.707930567827917 , 4060.990742193632092 , 14454.00184034344795 , 52076.66107599404803 , 189493.6591462156887 , 695184.5762413896145 , 2567994.048255284686 , 9541921.966748386322 , 35634927.44218076174 , 133669298.4612040871 , 503352186.6866284541 , 1901975729.538660119 , 7208915015.330103756)

    elseif 0.8 <= x <0.85 #Table 10
        t = x-0.825
        t =@horner(t , 2.318122621712510589 , 2.616920150291232841 , 7.897935075731355823 , 30.50239715446672327 , 131.4869365523528456 , 602.9847637356491617 , 2877.024617809972641 , 14110.51991915180325 , 70621.44088156540229 , 358977.2665825309926 , 1847238.263723971684 , 9600515.416049214109 , 50307677.08502366879 , 265444188.6527127967 , 1408862325.028702687 , 7515687935.373774627)

    elseif 0.85 <= x <0.9 #Table 11
        t= x-0.875
        t =@horner(t, 2.473596173751343912 , 3.727624244118099310 , 15.60739303554930496 , 84.12850842805887747 , 506.9818197040613935 , 3252.277058145123644 , 21713.24241957434256 , 149037.0451890932766 , 1043999.331089990839 , 7427974.817042038995 , 53503839.67558661151 , 389249886.9948708474 , 2855288351.100810619 , 21090077038.76684053 , 156699833947.7902014 , 1170222242422.439893 , 8777948323668.937971 , 66101242752484.95041 , 499488053713388.7989 , 37859743397240299.20)

    elseif x>=0.9
        td= 1-x
        td1=td-0.05
        qd = @horner(td, 0.0, (1.0/16.0), (1.0/32.0), (21.0/1024.0), (31.0/2048.0), (6257.0/524288.0), (10293.0/1048576.0), (279025.0/33554432.0), (483127.0/67108864.0),  (435506703.0/68719476736.0), (776957575.0/137438953472.0) , (22417045555.0/4398046511104.0) , (40784671953.0/8796093022208.0) , (9569130097211.0/2251799813685248.0) , (17652604545791.0/4503599627370496.0))

        kdm = @horner(td1 , 1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415 , 0.179481482914906162 , 0.144556057087555150 , 0.123200993312427711 , 0.108938811574293531 , 0.098853409871592910 , 0.091439629201749751 , 0.085842591595413900 , 0.081541118718303215)
        km = -Base.log(qd) * (kdm/pi)
        t = km
    end
    if flag == false
        return t
    elseif flag == true
        ans = t / sqrt(1.0-a)
        return ans
    end
end

#Using piecewise approximation polynomial as given in
#'Fast Computation of Complete Elliptic Integrals and Jacobian Elliptic Functions'
#                          Fukushima, Toshio. (2014). F09-FastEI. Celest Mech Dyn Astr
#                                                            DOI 10.1007/s10569-009-9228-z
#Link : https://pdfs.semanticscholar.org/8112/c1f56e833476b61fc54d41e194c962fbe647.pdf
#
#For m<0 , followed  Fukushima, Toshio. (2014). Precise, compact, and fast computation of complete elliptic integrals by piecewise minimax rational function approximation. Journal of Computational and Applied Mathematics. 282. 10.13140/2.1.1946.6245.
# Link:  https://www.researchgate.net/profile/Toshio_Fukushima/publication/267330394_Precise_compact_and_fast_computation_of_complete_elliptic_integrals_by_piecewise_minimax_rational_function_approximation/links/544b81a40cf2d6347f43074f/Precise-compact-and-fast-computation-of-complete-elliptic-integrals-by-piecewise-minimax-rational-function-approximation.pdf?origin=publication_detail
#Also suggested in this paper that we should consider domain only from (-inf,1].

"""
   ellipe(x)

 DLMF : https://dlmf.nist.gov/19.2#E5  , https://dlmf.nist.gov/19.2#E8
 Wiki : https://en.wikipedia.org/wiki/Elliptic_integral
Computes the Complete Elliptic Integral of 2nd kind at `x` ->E(x)---gven by: ``\\int_{0}^{\\pi/2} \\sqrt{1 - x(\\sin \\theta )^{2}} d\\theta``
"""
function ellipe(a::Float64)
    flag=false
    if a<0.0
        x=a/(a-1)    #for dealing with negative args
        flag=true
    elseif a>=0.0
        x=a
        flag=false

    end

    if x==0.0
        return pi/2
    elseif x==1.0
        return 1.0

    elseif x>1.0
        throw(DomainError(x,"`x` must lie between -Inf and 1 ---- Domain : (-inf,1.0]"))

    elseif 0.0 <= x < 0.1   #Table 2 from paper
        t= x-0.05
        t= @horner(t  ,  +1.550973351780472328  ,-0.400301020103198524  ,-0.078498619442941939  ,-0.034318853117591992  ,-0.019718043317365499  ,-0.013059507731993309  ,-0.009442372874146547  ,-0.007246728512402157  ,-0.005807424012956090  ,-0.004809187786009338)

    elseif 0.1 <= x < 0.2  #Table 3
        t=x-0.15
        t= @horner(t , +1.510121832092819728 , -0.417116333905867549 , -0.090123820404774569 , -0.043729944019084312, -0.027965493064761785 , -0.020644781177568105 , -0.016650786739707238 , -0.014261960828842520 , -0.012759847429264803 , -0.011799303775587354 , -0.011197445703074968)

    elseif 0.2 <= x < 0.3 #Table 4
        t=x-0.25
        t= @horner(t , 1.467462209339427155 , -0.436576290946337775 , -0.105155557666942554 , -0.057371843593241730 , -0.041391627727340220 , -0.034527728505280841 , -0.031495443512532783 , -0.030527000890325277 , -0.030916984019238900 , -0.032371395314758122 , -0.034789960386404158)

    elseif 0.3 <= x < 0.4 #Table 5
        t= x-0.35
        t= @horner(t,   +1.422691133490879171,  -0.459513519621048674,  -0.125250539822061878,  -0.078138545094409477,  -0.064714278472050002,  -0.062084339131730311,  -0.065197032815572477,-0.072793895362578779,  -0.084959075171781003,  -0.102539850131045997,  -0.127053585157696036,  -0.160791120691274606)

    elseif 0.4 <= x < 0.5 #Table 6
        t=x-0.45
        t= @horner(t  ,+1.375401971871116291  ,-0.487202183273184837  ,-0.153311701348540228  ,-0.111849444917027833  ,-0.108840952523135768  ,-0.122954223120269076  ,-0.152217163962035047  ,-0.200495323642697339  ,-0.276174333067751758  ,-0.393513114304375851  ,-0.575754406027879147  ,-0.860523235727239756  ,-1.308833205758540162)

    elseif 0.5 <= x < 0.6 #Table 7
        t= x-0.55
        t= @horner(t  ,+1.325024497958230082  ,-0.521727647557566767  ,-0.194906430482126213  ,-0.171623726822011264  ,-0.202754652926419141  ,-0.278798953118534762  ,-0.420698457281005762  ,-0.675948400853106021  ,-1.136343121839229244  ,-1.976721143954398261  ,-3.531696773095722506  ,-6.446753640156048150  ,-11.97703130208884026)
    elseif 0.6 <= x < 0.7 #Table 8
        t=x-0.65
        t= @horner(t,  +1.270707479650149744,  -0.566839168287866583,  -0.262160793432492598,  -0.292244173533077419,  -0.440397840850423189,  -0.774947641381397458,  -1.498870837987561088,  -3.089708310445186667,  -6.667595903381001064,  -14.89436036517319078,  -34.18120574251449024,  -80.15895841905397306,  -191.3489480762984920,  -463.5938853480342030,  -1137.380822169360061)

    elseif 0.7 <= x < 0.8 #Table 9
        t=x-0.75
        t= @horner(t,  +1.211056027568459525,  -0.630306413287455807,  -0.387166409520669145,  -0.592278235311934603,  -1.237555584513049844,  -3.032056661745247199,  -8.181688221573590762,  -23.55507217389693250,  -71.04099935893064956,  -221.8796853192349888,  -712.1364793277635425,  -2336.125331440396407,  -7801.945954775964673,  -26448.19586059191933,  -90799.48341621365251,  -315126.0406449163424,  -1104011.344311591159)

    elseif 0.8 <= x < 0.85 #Table 10
        t = x-0.825
        t= @horner(t,  +1.161307152196282836,  -0.701100284555289548,  -0.580551474465437362,  -1.243693061077786614,  -3.679383613496634879,  -12.81590924337895775,  -49.25672530759985272,  -202.1818735434090269,  -869.8602699308701437,  -3877.005847313289571,  -17761.70710170939814,  -83182.69029154232061,  -396650.4505013548170,  -1920033.413682634405)

    elseif 0.85 <= x < 0.9 #Table 11
        t= x-0.875
        t =  @horner(t,  +1.124617325119752213,  -0.770845056360909542,  -0.844794053644911362,  -2.490097309450394453,  -10.23971741154384360,  -49.74900546551479866,  -267.0986675195705196,  -1532.665883825229947,  -9222.313478526091951,  -57502.51612140314030,  -368596.1167416106063,  -2415611.088701091428,  -16120097.81581656797,  -109209938.5203089915,  -749380758.1942496220,  -5198725846.725541393,  -36409256888.12139973)

    elseif x>=0.9
        td= 1-x
        td1=td-0.05

        kdm = @horner(td1 , 1.591003453790792180 , 0.416000743991786912 , 0.245791514264103415  , 0.179481482914906162  , 0.144556057087555150 , 0.123200993312427711,0.108938811574293531 , 0.098853409871592910  ,0.091439629201749751  ,0.085842591595413900  ,0.081541118718303215)
        edm =  @horner(td1  ,+1.550973351780472328  ,-0.400301020103198524  ,-0.078498619442941939  ,-0.034318853117591992,-0.019718043317365499  ,-0.013059507731993309  ,-0.009442372874146547  ,-0.007246728512402157  ,-0.005807424012956090  ,-0.004809187786009338)
        hdm = kdm -edm
        km = ellipk(Float64(x))
        #em =  km + (pi/2 - km*edm)/kdm
        em = (pi/2 + km*hdm)/kdm   #to avoid precision loss near 1
        t= em
    end
    if flag == false
        return t
    elseif flag == true
        return t * sqrt(1.0-a)

    end
end

for f in (:ellipk,:ellipe)
    @eval begin
        ($f)(x::Float16) = Float16(($f)(Float64(x)))
	($f)(x::Float32) = Float32(($f)(Float64(x)))
        ($f)(x::Real) = ($f)(float(x))
	($f)(x::AbstractFloat) = throw(MethodError($f,(x,"")))
    end
end


text:
using Test
using LegendrePolynomials
using OffsetArrays
using HyperDualNumbers
using Aqua

@testset "Project quality assurance" begin
    Aqua.test_all(LegendrePolynomials)
end

import LegendrePolynomials: LegendrePolynomialIterator

tohyper(x) = Hyper(x, one(x), one(x), zero(x))

@testset "LegendrePolynomialIterator" begin
    lmax = 5
    iter = LegendrePolynomialIterator(0.5, lmax);
    @test length(iter) == lmax + 1
    @test eltype(iter) == Float64
    @test eltype(typeof(iter)) == Float64
    @test size(iter) == (lmax+1,)
    @test axes(iter) == (0:lmax,)
    @test keys(iter) == 0:lmax

    iter2 = copy(iter)
    @test typeof(iter2) == typeof(iter)
    @test iter2.x == iter.x
    @test iter2.lmax == iter.lmax

    iter = LegendrePolynomialIterator(0.5);
    @test Base.IteratorSize(typeof(iter)) == Base.IsInfinite()

    @test_throws DomainError LegendrePolynomialIterator(-1.1, lmax)
    @test_throws DomainError LegendrePolynomialIterator(-1.1)
end

@testset "Pl" begin
	x = 2rand() - 1
    lmax = 5
    P = collectPl(x, lmax = lmax)
    P2 = collectdnPl(x, lmax = lmax, n = 0)
    @test P[0] == P2[0] == 1
    @test P[1] == P2[1] == x
    @test P[2] ≈ P2[2] ≈ (3x^2 - 1)/2
    @test P[3] ≈ P2[3] ≈ (5x^3 - 3x)/2
    @test P[4] ≈ P2[4] ≈ (35x^4 - 30x^2 + 3)/8
    @test P[5] ≈ P2[5] ≈ (63x^5 - 70x^3 + 15x)/8

    @testset "x = 0" begin
        for l = 1:2:101
            @test Pl(0, l) == dnPl(0, l, 0) == 0
        end
    end

    @testset "x = ±1" begin
        lmax = 10

        P = collectPl(1, lmax = lmax)
        @test all(==(1), P)
        P2 = collectdnPl(1, lmax = lmax, n = 0)
        @test all(==(1), P2)

        P = collectPl(-1, lmax = lmax)
        P2 = collectdnPl(-1, lmax = lmax, n = 0)
        for l in axes(P, 1)
            @test P[l] == P2[l] == (-1)^l
        end
    end

    @test_throws DomainError collectPl(-2, lmax = lmax)
    @test_throws ArgumentError collectPl(0, lmax = -1)
end

@testset "dnPl" begin
    @test dnPl(0.5, 3, 10) == 0
    @test dnPl(0.5, 3, 10, zeros(100)) == 0
    @testset "double factorial" begin
        @test dnPl(0.5, 2, 2) == 3
        @test dnPl(0.5, 3, 3) == 15
        @test dnPl(0.5, 3, 3, zeros(1)) == 15
    end
    @test_throws ArgumentError dnPl(0.5, 10, 3, zeros(1))
    @test_throws ArgumentError dnPl(0.5, -1, 3)
    @test_throws ArgumentError dnPl(0.5, 1, -3)
    @test_throws DomainError dnPl(2.5, 1, 3)
end

@testset "collectPl!" begin
    v = zeros(10);
    x = 0.6

    collectPl!(v, x)
    for l in 0:length(v) - 1
        @test v[firstindex(v) + l] == Pl(x, l)
    end

    v .= 0
    lmax = 6
    collectPl!(v, x, lmax = lmax)
    for l in 0:lmax
        @test v[firstindex(v) + l] == Pl(x, l)
    end

    @test_throws DomainError collectPl!(v, -2, lmax = lmax)
    @test_throws ArgumentError collectPl!(v, 0, lmax = -1)
end

@testset "collectdnPl!" begin
    v = zeros(10);
    @test_throws DomainError collectdnPl!(v, -2, lmax = 3, n = 2)
    @test_throws ArgumentError collectdnPl!(v, 0, lmax = -1, n = 3)
    @test_throws ArgumentError collectdnPl!(v, 0.5, lmax = 1, n = -3)
    @test_throws ArgumentError collectdnPl!(zeros(1), 0, lmax = 3, n = 3)
end

@testset "collectdnPl" begin
    @test_throws DomainError collectdnPl(-2, lmax = 3, n = 2)
    @test_throws ArgumentError collectdnPl(0, lmax = -1, n = 3)
    @test_throws ArgumentError collectdnPl(0.5, lmax = 1, n = -3)
end

@testset "dPl/dx" begin
    x = 2rand() - 1
    xh = tohyper(x)
    lmax = 5
    P = eps1.(collectPl(xh, lmax = lmax))
    P2 = collectdnPl(x, lmax = lmax, n = 1)
    @test P[0] == P2[0] == 0
    @test P[1] == P2[1] == 1
    @test P[2] ≈ P2[2] ≈ 3x
    @test P[3] ≈ P2[3] ≈ (-3 + 15x^2)/2
    @test P[4] ≈ P2[4] ≈ 1/8*(-60x + 140x^3)
    @test P[5] ≈ P2[5] ≈ 1/8*(15 - 210x^2 + 315x^4)

    @testset "x = ±1" begin
        lmax = 10
        x = 1
        P = collectPl(tohyper(x), lmax = lmax)
        dP1h = eps1.(P)
        dP1 = collectdnPl(x, lmax = lmax, n = 1)

        for l in axes(dP1,1)
            @test dP1h[l] == dP1[l] == l*(l+1)/2
        end

        x = -1
        P = collectPl(tohyper(x), lmax = lmax)
        dPm1h = eps1.(P)
        dPm1 = collectdnPl(x, lmax = lmax, n = 1)

        for l in axes(dPm1,1)
            @test dPm1h[l] == dPm1[l] == (-1)^(l+1) * l*(l+1)/2
        end
    end
end

@testset "d^2Pl/dx^2" begin
    x = 2rand() - 1
    xh = tohyper(x)
    lmax = 5
    P = eps1eps2.(collectPl(xh, lmax = lmax))
    P2 = collectdnPl(x, lmax = lmax, n = 2)

    @test P[0] == P2[0] == 0
    @test P[1] == P2[1] == 0

    @test P[2] ≈ P2[2] ≈ 3
    @test P[3] ≈ P2[3] ≈ 15x
    @test P[4] ≈ P2[4] ≈ 1/8*(-60 + 420x^2)
    @test P[5] ≈ P2[5] ≈ 1/8*(-420x + 1260x^3)
end

@testset "Askey-Gasper" begin
    x = 0.5
    for lmax=1:100:1_000
        @test sum(collectPl(x, lmax = lmax)) >= 0
        @test sum(collectdnPl(x, lmax = lmax, n = 0)) >= 0
    end
end


text:
using MacroTools: postwalk, @capture

@testset "Tfunc" begin
    dict = OrderedDict(
        :age_t0 => :(30),
        :duration_t0 => :(6),
        :premium_t0 => :(100),
        :sum_assured_t0 => :(10000),
        :age => :(floor(Int, age_t0 + t/12)),
        :duration => :(duration_t0 + t),
        :qx => :(age / 1000),
        :lapse_rate_m => :(1 / qx / 1000),
        :no_pols_som => :(if t==1 1 else no_pols_eom[t-1] end),
        :no_deaths => :(no_pols_som*qx/12*(1-0.5*lapse_rate_m)),
        :no_surrenders => :(no_pols_som*lapse_rate_m*(1-0.5*qx)),
        :no_pols_eom => :(if t==1 1 else no_pols_eom[t-1] end - no_deaths - no_surrenders),
        :premium_income => :(premium_t0*no_pols_eom/12),
        :claims_outgo => :(sum_assured_t0*no_deaths),
        :profit => :(premium_income-claims_outgo),
        :yield_curve => :(0.1),
        :profit_discounted => :(profit*(1+yield_curve)^(-t)),
        :liability => :(-profit_discounted+ if t<600 liability[t+1] else 0 end)
    )
    @test transform_expr_Tfunc(OrderedDict(i=>CodePoint(j, nothing) for (i,j) ∈ dict)) isa Expr

    T = collect(1:600)
    f₁ = @Tfunc begin
        age_t0 = 30
        duration_t0 = 6
        premium_t0 = 100
        sum_assured_t0 = 10000
        age = floor(Int, age_t0 + t/12)
        duration = duration_t0 + t
        qx = age / 1000
        lapse_rate_m = 1 / qx / 1000
        no_pols_som = if t==1 1 else no_pols_eom[t-1] end
        no_deaths = no_pols_som*qx/12*(1-0.5*lapse_rate_m)
        no_surrenders = no_pols_som*lapse_rate_m*(1-0.5*qx)
        no_pols_eom = if t==1 1 else no_pols_eom[t-1] end - no_deaths - no_surrenders
        premium_income = premium_t0*no_pols_eom/12
        claims_outgo = sum_assured_t0*no_deaths
        profit = premium_income-claims_outgo
        yield_curve = 0.1
        profit_discounted  = profit*(1+yield_curve)^(-t)
        liability = -profit_discounted+ if t<600 liability[t+1] else 0 end
        testme = no_deaths*t*profit*0 + age_t0
    end

    @test [f₁(;age_t0).testme[end] for age_t0 ∈ 10:5:1000] == collect(10:5:1000)
    @test f₁().age == @.(floor(Int, 30 + T/12))

end

transform_expr_Tfunc(formulas::OrderedDict{Symbol, CodePoint}; fname::Union{Symbol, Nothing}=nothing) where T = begin
    variables = collect(keys(formulas))
    fguts = transform_expr_all(formulas)

    # "send" input as kwargs and "recieve" output as named tuples
    fkwargs = Expr(:parameters, [Expr(:kw, var, nothing) for var ∈ variables]...)
    freturn = Expr(:tuple, [Expr(:(=), var, var) for var ∈ variables]...)

    # TODO: build seperate pipeline not to have to mutate
    # Mutate the function to suppress variable allocation when kwargs are passed, I think this is quite expensive
    varinits = Dict(var=>gensymx(string(var)*"_init") for var in variables)

    finits = Expr(:block,
        [Expr(:(=), varinits[var], Expr(:call, :(===), var, :nothing)) for var ∈ variables]...
    )
    for var in variables
        fguts = postwalk(x->@capture(x, $var[_] = __) ? Expr(:if, varinits[var], Expr(:block, formulas[var].line, x)) : x, fguts)
        fguts = postwalk(x->@capture(x, $var = __) ? Expr(:if, varinits[var], Expr(:block, formulas[var].line, x)) : x, fguts)
    end

    # Different form for named and unnamed functions
    funcbody = Expr(:block, finits, fguts, freturn)
    if fname === nothing
        Expr(:->, Expr(:tuple, fkwargs), funcbody)
    else
        Expr(:(=), Expr(:call, fname, fkwargs), funcbody)
    end

end

macro Tfunc(expr)
    formulas = expr_to_assignment_dict(expr; line=__source__)
    esc(transform_expr_Tfunc(formulas))
end

macro Tfunc(fname::Symbol, expr)
    formulas = expr_to_assignment_dict(expr; line=__source__)
    esc(transform_expr_Tfunc(formulas; fname))
end


text:
#provide a pair of points in latitude and longitude, in 'degrees' to get
#the distance back in Km
function getDistanceFromLatLongInKm(lat1,long1,lat2,long2)
  RadiusEarth = 6371; # Earth's radius in km
  dLat = deg2rad(lat2-lat1); # deg2rad, built in fn in Julia
  dLong = deg2rad(long2-long1);
  h = sin(dLat/2)*sin(dLat/2) + cos(deg2rad(lat1))*cos(deg2rad(lat2))*sin(dLong/2)*sin(dLong/2)
  d = RadiusEarth * 2 * atan2(sqrt(h),sqrt(1-h)); # Distance in km
  return d;
end

max_stars_repo_name: mantzaris/distLatLong

max_issues_repo_name: mantzaris/distLatLong

max_forks_repo_name: mantzaris/distLatLong


text:
using JSON;

@enum MessageType begin
    RUN_COMMAND
    GET_WORKSPACES
    SUBSCRIBE
    GET_OUTPUTS
    GET_TREE
    GET_MARKS
    GET_BAR_CONFIG
    GET_VERSION
    GET_BINDING_MODES
    GET_CONFIG
    SEND_TICK
    SYNC
    GET_BINDING_STATE
end

struct Message
    type::MessageType
    payload::String
    function Message(type, payload)
        return vcat(b"i3-ipc",
                    reinterpret(UInt8, UInt32[length(payload)]),
                    reinterpret(UInt8, UInt32[UInt32(type)]),
                    Vector{UInt8}(payload))
    end
end


text:
# Warning: setupLLCgrid no longer exists
# Example: read monthly average temperature/salinity fields on poseidon
# Here compare reading with read_mdsio versus MeshArrays built-in method.
# At bottom, also comparing statistics computed via two types of files.
using Revise
using ECCOtour, MeshArrays, MITgcmTools, Statistics

# other useful packages in Reemergence project
# using Statistics, PyPlot, Distributions, FFTW, LinearAlgebra, StatsBase
path_grid="../inputs/GRID_LLC90/"
γ = setupLLCgrid(path_grid)

path = "/poseidon/ecco/ECCOv4r4/MITgcm/exps/iter129_bulkformula/run/diags/"
fileroot = "state_3d_set1" # contains θ, S, dρ/dz, see accompanying meta file

# get a list of all files in directory.
filelist = searchdir(path,fileroot) # first filter for state_3d_set1

# just save the list of ".data" files, eliminate ".meta"
datafilelist  = filter(x -> occursin("data",x),filelist) # second filter for "data"

# should read this from the meta file.
nz = 50  # number of depths
nc = 3   # number of variables

# read the first one into a 4D array: lon, lat, depth, variable type
@time state1 = read_mdsio(path,datafilelist[1]);
@time state2 = γ.read(path*datafilelist[1],MeshArray(γ,Float32,nz*nc));

#compare salinity at depth level 20.
level = 48
propertyno = 2;
ntiles = size(state2,1)
S1 = state1[:,:,level,propertyno]
S2 = state2[:,(propertyno-1)*nz+level]

# compare maxima
xcount = [sum(count(!iszero,S2[i])) for i = 1:ntiles]
@time maxtest = maximum([maximum(filter(!iszero,S2[i])) for i = 1:ntiles if xcount[i] > 0])
@assert maxtest == maximum(filter(!iszero,S1)) # filter out zeroes

# compare minimum
@time mintest = minimum([minimum(filter(!iszero,S2[i])) for i = 1:ntiles if xcount[i] > 0])
@assert mintest == minimum(filter(!iszero,S1)) # filter out zeroes

# compute mean of S1, S2
xsum = sum([sum(S2[i]) for i = 1:ntiles if xcount[i] > 0]) # works b.c. 0 on land
xbar = xsum/sum(xcount)
@assert xbar == mean(filter(!iszero,S1))

x²sum = sum([sum((filter(!iszero,S2[i]).-xbar).^2) for i=1:ntiles if xcount[i]>0])
σ = sqrt(x²sum/(sum(xcount)-1))
@assert σ == std(filter(!iszero,S1))


text:
# Convenience function that calls the more complicated simulatedata
# function with fixed arguments.
"""
    t, fluxes, errors, truegalaxyvector = mockdata()

Returns 50 regular (i.e. co-occuring/simultaneous) observations in three filters
and the true galaxy vector.
`fluxes` and `errors` have both dimensions 3 × 50.
`truegalaxyvector` is a 3 dimensional vector.
"""
function mockdata()

    PyPlot.ioff()

    Temp = 5000.0
    Sc   =   10.0

    _, _, _, _, regflux, regσ, _, obstimeco = @suppress simulatedata(N=50, SN=50.0, Temperature=Temp, offset=5.0, scale=Sc, logs=-6.0);
    close("all")
    PyPlot.ion()

    return obstimeco, Matrix(regflux[:,1:3]'), Matrix(regσ[:,1:3]'), Sc * vec(unitgalaxyvector(Temp)[1:3])

end


text:
"""
The LaTeXStrings module exists mainly to make LaTeX equations easier to type as
literal strings, and so that the resulting strings display as formatted equations
in supporting environments like IJulia.

See in particular the `LaTeXString` type and the `L"..."` constructor macro.
"""
module LaTeXStrings
export LaTeXString, latexstring, @L_str

# IJulia supports LaTeX output for any object with a text/latex
# writemime method, but these are annoying to type as string literals
# in Julia because of all the escaping required, e.g. "\$\\alpha +
# \\beta\$".  To simplify this, we add a new string type with a macro
# constructor, so that one can simply do L"$\alpha + \beta$".

@doc raw"""
A `LaTeXString` is a string type whose contents represent a fragment of LaTeX code,
typically containing an equation (`$...$`).   In certain environments (e.g. IJulia)
this will display with LaTeX-like formatting.   For the most part, you can use
a `LaTeXString` object in any context that expects an `AbstractString` object.

The `L"..."` macro is convenient for constructing `LaTeXString` objects, because
it eliminates the need to escape backslashes and dollar signs, and implicitly inserts
dollar signs around the string if none are present.  For example, `L"$\alpha$"`, `L"\alpha"`,
and `LaTeXString("\$\\alpha\$")` are all equivalent.
"""
struct LaTeXString <: AbstractString
    s::String
end

# the only point of using LaTeXString to represent equations, since
# IJulia doesn't support LaTeX output other than equations, so add $'s
# around the string if they aren't there (ignoring \$ and %$)
_maybe_wrap_equation(s) = occursin(r"[^\\%]\$|^\$", s) ? s : string('\$', s, '\$')

"""
    latexstring(args...)

Similar to `string(args...)`, but generates a `LaTeXString` instead of a `String`.
"""
latexstring(args...) = latexstring(string(args...))
function latexstring(s::String)
    return LaTeXString(_maybe_wrap_equation(s))
end
latexstring(s::AbstractString) = latexstring(String(s))

if isdefined(Meta, :parseatom)
    const parseatom = Meta.parseatom
else
    parseatom(s, i; filename=nothing) = Meta.parse(s, i; greedy=false)
end

@doc raw"""
    L"..."

Creates a `LaTeXString` and is equivalent to `latexstring(raw"...")`, except that
`%$` can be used for interpolation.

```jldoctest
julia> L"x = \sqrt{2}"
L"$x = \sqrt{2}$"

julia> L"x = %$(sqrt(2))"
L"$x = 1.4142135623730951$"
```
"""
macro L_str(s::String)
    i = firstindex(s)
    buf = IOBuffer(maxsize=ncodeunits(s))
    ex = Expr(:call, GlobalRef(LaTeXStrings, :latexstring))
    while i <= ncodeunits(s)
        c = @inbounds s[i]
        i = nextind(s, i)
        if c === '%' && i <= ncodeunits(s)
            c = @inbounds s[i]
            if c === '$'
                position(buf) > 0 && push!(ex.args, String(take!(buf)))
                atom, i = parseatom(s, nextind(s, i), filename=string(__source__.file))
                Meta.isexpr(atom, :incomplete) && error(atom.args[1])
                atom !== nothing && push!(ex.args, atom)
                continue
            else
                print(buf, '%')
            end
        else
            print(buf, c)
        end
    end
    position(buf) > 0 && push!(ex.args, String(take!(buf)))
    return esc(ex)
end

Base.write(io::IO, s::LaTeXString) = write(io, s.s)
Base.show(io::IO, ::MIME"application/x-latex", s::LaTeXString) = print(io, s.s)
Base.show(io::IO, ::MIME"text/latex", s::LaTeXString) = print(io, s.s)
function Base.show(io::IO, s::LaTeXString)
    @static if isdefined(Base, :escape_raw_string)  # Julia ≥ 1.4
        print(io,"L\"")
        Base.escape_raw_string(io, s.s)
        print(io,'"')
    else
        print(io, 'L')
        Base.print_quoted_literal(io, s.s)   # Julia < 1.6
    end
end

Base.firstindex(s::LaTeXString) = firstindex(s.s)
Base.lastindex(s::LaTeXString) = lastindex(s.s)
Base.iterate(s::LaTeXString, i::Int) = iterate(s.s, i)
Base.iterate(s::LaTeXString) = iterate(s.s)
Base.nextind(s::LaTeXString, i::Int) = nextind(s.s, i)
Base.prevind(s::LaTeXString, i::Int) = prevind(s.s, i)
Base.eachindex(s::LaTeXString) = eachindex(s.s)
Base.length(s::LaTeXString) = length(s.s)
Base.getindex(s::LaTeXString, i::Integer) = getindex(s.s, i)
Base.getindex(s::LaTeXString, i::Int) = getindex(s.s, i) # for method ambig in Julia 0.6
Base.getindex(s::LaTeXString, i::UnitRange{Int}) = getindex(s.s, i)
Base.getindex(s::LaTeXString, i::UnitRange{<:Integer}) = getindex(s.s, i)
Base.getindex(s::LaTeXString, i::AbstractVector{<:Integer}) = getindex(s.s, i)
Base.codeunit(s::LaTeXString, i::Integer) = codeunit(s.s, i)
Base.codeunit(s::LaTeXString) = codeunit(s.s)
Base.ncodeunits(s::LaTeXString) = ncodeunits(s.s)
Base.codeunits(s::LaTeXString) = codeunits(s.s)
Base.sizeof(s::LaTeXString) = sizeof(s.s)
Base.isvalid(s::LaTeXString, i::Integer) = isvalid(s.s, i)
Base.pointer(s::LaTeXString) = pointer(s.s)
Base.IOBuffer(s::LaTeXString) = IOBuffer(s.s)
Base.unsafe_convert(T::Union{Type{Ptr{UInt8}},Type{Ptr{Int8}},Cstring}, s::LaTeXString) = Base.unsafe_convert(T, s.s)

end # module


text:
# Raw BPF calls

const SYS_bpf = 321

function bpf(cmd::API.BPFCmd, attr::Ref{T}) where {T}
    vec = zeros(UInt8, 120)
    GC.@preserve vec attr begin
        p_attr = Base.unsafe_convert(Ptr{Cvoid}, Base.cconvert(Ptr{Cvoid}, attr))
        ccall(:memcpy, Ptr{Cvoid}, (Ptr{UInt8},   Ptr{Cvoid}, Csize_t),
                                    pointer(vec), p_attr,     sizeof(T))
        ptr = Ptr{Cvoid}(pointer(vec))
        res = ccall(:syscall, Cint, (Clong,    Cint,      Ptr{Cvoid}, Cuint),
                                     SYS_bpf,  Cint(cmd), ptr,        Cuint(120)) # FIXME: sizeof(T)
        ccall(:memcpy, Ptr{Cvoid}, (Ptr{Cvoid}, Ptr{UInt8},   Csize_t),
                                    p_attr,     pointer(vec), sizeof(T))
        return res
    end
end


text:
using DataStructures
using ProgressMeter

mutable struct NGram{S<:AbstractString, T}
    n::Int
    text::S
    inds::Dict{T, Set{Int}}
    counts::Dict{T, Float64}
    function NGram{S, T}(text) where {S<:AbstractString, T}
        new(length(text), text, Dict{T, Set{Int}}(), Dict{T, Int}())
    end
end

entropy(ngram::NGram) = sum(values(ngram.counts)) * (ngram.n - 1)

struct Utterance{S<:AbstractString}
    text::S
    count::Int
    ngrams::Set{NGram{S, Utterance{S}}}
    function Utterance{S}(text::S, count::Int) where {S<:AbstractString}
        new(text, count, Set{NGram{S, Utterance{S}}}())
    end
end

NGram(text::S) where {S} = NGram{S, Utterance{S}}(text)
Utterance(text::S, count::Int) where {S} = Utterance{S}(text, count)

function overlaps(utterance, ngram1, ngram2)
    inds1, inds2 = ngram1.inds[utterance], ngram2.inds[utterance]
    nb1, nb2 = sizeof(ngram1.text), sizeof(ngram2.text)
    ret = 0
    for i1 in inds1
        for i2 in inds2
            # if ngram1.text == " other " || ngram2.text == " other "
            #     println("$(ngram1.text), $(ngram2.text), $i1, $i2, $(nb1), $(nb2)")
            # end
            if i2 <= i1 <= i1 + nb1 <= i2 + nb2
                ret += 1
            elseif i1 <= i2 <= i2 + nb2 <= i1 + nb1
            elseif (i1 <= i2 < i1 + nb1
                 || i2 <= i1 < i2 + nb2)
                ret += 1
            end
            # in " other "
            # " other ": 1, 7
            # "r ": 6, 2
        end
    end
    ret = ret / (length(inds2))# * length(inds1))
    #display((ret, utterance.text, ngram1.text, ngram2.text))
    return ret
end

function addinstance!(ngram::NGram, utterance::Utterance, i::Int)
    ngram.counts[utterance] = utterance.count
    inds = get!(Set{Int}, ngram.inds, utterance)
    push!(inds, i)
    push!(utterance.ngrams, ngram)
end

function buildngrams!(utterance::Utterance{S}, ngrams) where {S}
    utext = utterance.text
    i = 1
    while i < endof(utext)
        j = i
        while j < endof(utext)
            j = nextind(utext, j)
            text = utterance.text[i:j]
            ngram = get!(()->NGram(text), ngrams, text)
            addinstance!(ngram, utterance, i)
        end
        i = nextind(utext, i)
    end
end

function buildngrams(counts::Accumulator{S}) where {S}
    ngrams = Dict{S, NGram{S, Utterance{S}}}()
    @showprogress 1 "building ngrams: " for (text, count) in counts
        utterance = Utterance(text, count)
        buildngrams!(utterance, ngrams)
    end
    return values(ngrams)
end

function fix_double_counts!(oldng::NGram, newng::NGram, utterance::Utterance, ngrams::PriorityQueue)
    #println("fixing dc for $((oldng.text, newng.text)) in \"$(utterance.text)\"")
    seen = Set([oldng, newng])
    for ngram in utterance.ngrams
        if !(ngram in seen)
            olfrac1 = overlaps(utterance, ngram, oldng)
            olfrac2 = overlaps(utterance, ngram, newng)
            olfrac = olfrac1 * olfrac2
            if ngram.text == " other "
                println("fixing dc for $((oldng.text, newng.text)) in \"$(utterance.text)\"")
                print("$olfrac1, $olfrac2, $(sum(values(ngram.counts))), $(utterance.count) -- from $(entropy(ngram)) ")
            end
            ngram.counts[utterance] += utterance.count * olfrac
            if ngram.text == " other "
                println("to $(entropy(ngram))")
            end
            if olfrac > 0 && ngram in keys(ngrams)
                ngrams[ngram] = entropy(ngram)
            end
            push!(seen, ngram)
        end
    end
end

function decrement_overlaps!(ngram::NGram, ngrams::PriorityQueue)
    for utterance in keys(ngram.inds)
        seen = Set([ngram])
        for ngram2 in utterance.ngrams
            if !(ngram2 in seen)
                #display((ngram.text,ngram2.text=>ngram2.count, utterance.text=>utterance.count, overlaps(utterance, ngram2, ngram)))
                olfrac = overlaps(utterance, ngram2, ngram)
                # if ngram2.text == " other "
                #     print("in \"$(utterance.text)\": $olfrac, $(sum(values(ngram2.counts))), $(utterance.count) -- from $(entropy(ngram2)) ")
                # end
                ngram2.counts[utterance] -= utterance.count * olfrac
                ngram2.counts[utterance] = max(ngram2.counts[utterance], 0)
                # if ngram2.text == " other "
                #     println("to $(entropy(ngram2))")
                # end
                if olfrac > 0 && ngram2 in keys(ngrams)
                    #otherwise maybe update?
                    ngrams[ngram2] = (entropy(ngram2), ngram2.text)
                elseif olfrac > 0
                    #fix_double_counts!(ngram, ngram2, utterance, ngrams)
                end
                # if ngram2.text in vocab
                #     pop!(vocab, ngram2.text)
                # end
                # ngrams[ngram2] = entropy(ngram2)
                push!(seen, ngram2)
            end
        end
    end
end

function buildvocab(counts::Accumulator{S}, maxsize::Int) where {S}
    #display(counts)
    vocab = OrderedDict{S, Float64}(counter(string.(convert(Vector{Char}, join(keys(counts), "")))))
    #vocab = OrderedDict{S, Float64}()
    ngrams = Dict(ng => (entropy(ng), ng.text) for ng in buildngrams(counts))
    ngrams = PriorityQueue(ngrams, Base.Order.Reverse)
    p = Progress(maxsize - length(vocab), 1, "building subword vocab: ")
    while length(vocab) < maxsize
        #display(structure(ngrams))
        best, (bestentropy, besttext) = dequeue_pair!(ngrams)
        if bestentropy <= 0
            break
        end
        #println("decrementing \"$(best.text)\" with entropy $bestentropy")
        decrement_overlaps!(best, ngrams)
        vocab[besttext] = bestentropy
        next!(p)
    end
    return vocab
end

function buildvocab(counts::Associative, maxsize::Number)
    buildvocab(counter(convert(Dict{String, Int}, counts)), convert(Int, maxsize))
end

function segment(utext::S, vocab::Associative{S, Float64}) where {S <: AbstractString}
    if length(utext) == 1
        return [utext]
    end
    vocab = keys(vocab)
    i, segments = 1, OrderedDict(1 => Vector{S}())
    while true
        j = i
        while j <= endof(utext)
            nextj = j < endof(utext) ? nextind(utext, j): j + 1
            subword = utext[i:j]
            #display(subword => subword in vocab)
            if subword in vocab
                curlen = if nextj in keys(segments)
                    length(segments[nextj])
                else
                    length(utext) + 1
                end
                if length(segments[i]) + 1 < curlen
                    segments[nextj] = segments[i][:]
                    push!(segments[nextj], subword)
                end
            end
            j = nextj
        end
        inds = collect(keys(segments))
        sort!(inds)
        firstind = findfirst(inds .== i)
        if firstind < length(inds)
            i = inds[firstind + 1]
        else
            break
        end
    end
    #display(segments)
    return segments[endof(utext) + 1]
end

function segment(utext::AbstractString, vocab::Associative)
    segment(utext, OrderedDict{String, Float64}(vocab))
end

function segment(utexts::AbstractVector, vocab::Associative)
    [tok for utext in utexts for tok in segment(utext, vocab)]
end


text:
struct PosteriorLTISDE{Tprior<:LTISDE, Tdata} <: AbstractGP
    prior::Tprior
    data::Tdata
end

# Avoids method ambiguity.
posterior(fx::FiniteLTISDE, y::AbstractVector) = _posterior(fx, y)
posterior(fx::FiniteLTISDE, y::AbstractVector{<:Real}) = _posterior(fx, y)

_posterior(fx, y) = PosteriorLTISDE(fx.f, (y=y, x=fx.x, Σy=fx.Σy))

const FinitePosteriorLTISDE = FiniteGP{<:PosteriorLTISDE}

function AbstractGPs.cov(fx::FinitePosteriorLTISDE)
    @error "Intentionally not implemented. Please don't try to explicitly compute this cov. matrix."
end

function AbstractGPs.marginals(fx::FinitePosteriorLTISDE)
    if fx.x != fx.f.data.x
        x, Σys, ys, _, pr_indices = build_inference_data(fx.f, fx.x)

        model = build_lgssm(fx.f.prior, x, Σys)
        Σys_new = noise_var_to_time_form(fx.x, fx.Σy)
        σ²s_pr_full = build_prediction_obs_vars(pr_indices, x, Σys_new)
        model_post = replace_observation_noise_cov(posterior(model, ys), σ²s_pr_full)
        return destructure(x, map(marginals, marginals(model_post))[pr_indices])
    else
        f = Zygote.literal_getfield(fx, Val(:f))
        prior = Zygote.literal_getfield(f, Val(:prior))
        x = Zygote.literal_getfield(fx, Val(:x))
        data = Zygote.literal_getfield(f, Val(:data))
        Σy = Zygote.literal_getfield(data, Val(:Σy))
        Σy_diag = Zygote.literal_getfield(Σy, Val(:diag))
        y = Zygote.literal_getfield(data, Val(:y))

        Σy_new = Zygote.literal_getfield(fx, Val(:Σy))

        model = build_lgssm(AbstractGPs.FiniteGP(prior, x, Σy))
        Σys_new = noise_var_to_time_form(x, Σy_new)
        ys = observations_to_time_form(x, y)
        model_post = replace_observation_noise_cov(posterior(model, ys), Σys_new)
        return destructure(x, map(marginals, marginals(model_post)))
    end
end

function AbstractGPs.mean_and_var(fx::FinitePosteriorLTISDE)
    ms = marginals(fx)
    return map(mean, ms), map(var, ms)
end

AbstractGPs.mean(fx::FinitePosteriorLTISDE) = mean_and_var(fx)[1]

AbstractGPs.var(fx::FinitePosteriorLTISDE) = mean_and_var(fx)[2]

function AbstractGPs.rand(rng::AbstractRNG, fx::FinitePosteriorLTISDE)
    x, Σys, ys, _, pr_indices = build_inference_data(fx.f, fx.x)

    # Get time-form observation vars and observations.
    Σys_pr = noise_var_to_time_form(fx.x, fx.Σy)

    model = build_lgssm(fx.f.prior, x, Σys)
    σ²s_pr_full = build_prediction_obs_vars(pr_indices, x, Σys_pr)
    model_post = replace_observation_noise_cov(posterior(model, ys), σ²s_pr_full)
    return destructure(x, rand(rng, model_post)[pr_indices])
end

AbstractGPs.rand(fx::FinitePosteriorLTISDE) = rand(Random.GLOBAL_RNG, fx)

function AbstractGPs.logpdf(fx::FinitePosteriorLTISDE, y_pr::AbstractVector{<:Real})

    x, Σys, ys, tr_indices, pr_indices = build_inference_data(
        fx.f, fx.x, fx.Σy, fill(missing, length(y_pr)),
    )

    # Get time-form observation vars and observations.
    Σys_pr = noise_var_to_time_form(fx.x, fx.Σy)
    ys_pr = observations_to_time_form(fx.x, y_pr)

    Σys_pr_full = build_prediction_obs_vars(pr_indices, x, Σys_pr)
    ys_pr_full = build_prediction_obs(tr_indices, pr_indices, x, ys_pr)

    model = build_lgssm(fx.f.prior, x, Σys)
    model_post = replace_observation_noise_cov(posterior(model, ys), Σys_pr_full)
    return logpdf(model_post, ys_pr_full)
end

# Join the dataset used to construct `f` and the one specified by `x_pred`, `σ²s_pred`, and
# `y_pred`. This is used in all of the inference procedures. Also provide a collection of
# indices that can be used to obtain the requested prediction locations from the joint data.

# This is the most naive way of going about this.
# The present implementation assumes that there are no overlapping data, which will lead to
# numerical issues if violated.
function build_inference_data(
    f::PosteriorLTISDE,
    x_pred::AbstractVector,
    Σy_pred::AbstractMatrix{<:Real},
    y_pred::AbstractVector{<:Union{Missing, Real}},
)
    d = f.data
    return merge_datasets((d.x, x_pred), (d.Σy, Σy_pred), (d.y, y_pred))
end

function merge_datasets(
    (x1, x2)::Tuple{AbstractVector, AbstractVector},
    (Σy1, Σy2)::Tuple{AbstractMatrix{<:Real}, AbstractMatrix{<:Real}},
    (y1, y2)::Tuple{
        AbstractVector{<:Union{Missing, Real}}, AbstractVector{<:Union{Missing, Real}},
    },
)
    # Merge and sort the inputs temporally.
    x_raw = merge_inputs(x1, x2)
    sort_indices, x = sort_in_time(x_raw)

    # Merge and sort the observation covariances temporally.
    Σys_raw = vcat(noise_var_to_time_form(x1, Σy1), noise_var_to_time_form(x2, Σy2))
    Σys = Σys_raw[sort_indices]

    # Merge and sort the observations temporally.
    ys_1 = observations_to_time_form(x1, y1)
    ys_2 = observations_to_time_form(x2, y2)
    ys_raw = vcat(ys_1, ys_2)
    ys = ys_raw[sort_indices]

    # The last length(x_pred) indices belong to the predictions.
    x1_indices = sortperm(sort_indices)[1:length(ys_1)]
    x2_indices = sortperm(sort_indices)[end-length(ys_2) + 1:end]

    return x, Σys, ys, x1_indices, x2_indices
end

# If no observations or variances are provided, make the observations arbitrary and the
# variances very large to simulate missing data.
function build_inference_data(f::PosteriorLTISDE, x_pred::AbstractVector)
    Σy_pred = Diagonal(fill(_large_var_const(), length(x_pred)))
    y_pred = fill(missing, length(x_pred))
    return build_inference_data(f, x_pred, Σy_pred, y_pred)
end

# Functions that make predictions at new locations require missings to be placed at the
# locations of the training data.
function build_prediction_obs_vars(
    pr_indices::AbstractVector{<:Integer},
    x_full::AbstractVector,
    Σys_pr::AbstractVector,
)
    σ²s_pr_full = get_zeros(x_full)
    σ²s_pr_full[pr_indices] .= Σys_pr
    return σ²s_pr_full
end

get_zeros(x::AbstractVector{T}) where {T<:Real} = zeros(T, length(x))

function build_prediction_obs(
    tr_indices::AbstractVector{<:Integer},
    pr_indices::AbstractVector{<:Integer},
    x_full::AbstractVector,
    y_pr::AbstractVector{T},
) where {T}
    y_pr_full = Vector{Union{Missing, T}}(undef, length(get_times(x_full)))
    y_pr_full[tr_indices] .= missing
    y_pr_full[pr_indices] .= y_pr
    return y_pr_full
end


text:
module EDF

using Dates

#####
##### Types
#####

"""
    EDF.PatientID

Type representing the local patient identification field of an EDF header, assuming the
file is EDF+ compliant.
`EDF.File`s which are parsed from files which are not EDF+ compliant do not contain this
type; the corresponding field is instead a `String`.

# Fields

* `code` (`String` or `Missing`): The code by which a patient is referred, if known
* `sex` (`Char` or `Missing`): Patient sex, `'M'`, `'F'`, or `missing` if unknown
* `birthdate` (`Date` or `Missing`): Patient date of birth, if known
* `name` (`String` or `Missing`): Patient name, if known
"""
struct PatientID
    code::Union{String,Missing}
    sex::Union{Char,Missing}
    birthdate::Union{Date,Missing}
    name::Union{String,Missing}
end

"""
    EDF.RecordingID

Type representing the local recording identification field of an EDF header, assuming
the file is EDF+ compliant.
`EDF.File`s which are parsed from files which are not EDF+ compliant do not contain
this type; the corresponding field is instead a `String`.

# Fields

* `startdate` (`Date` or `Missing`): Start date of the recording
* `admincode` (`String` or `Missing`): Administration code for the recording
* `technician` (`String` or `Missing`): Identifier for the technician or investigator
* `equipment` (`String` or `Missing`): Identifier for the equipment used
"""
struct RecordingID
    startdate::Union{Date,Missing}
    admincode::Union{String,Missing}
    technician::Union{String,Missing}
    equipment::Union{String,Missing}
end

"""
    EDF.AnnotationsList

Type representing a time-stamp annotations list (TAL).

# Fields

* `offset` (`Float64`): Offset from the recording start time (specified in the header)
  at which the event in this TAL starts
* `duration` (`Float64` or `Nothing`): Duration of the event, if specified
* `event` (`Vector{String}`): List of events for this TAL
"""
struct AnnotationsList
    offset::Float64
    duration::Union{Float64,Nothing}
    event::Vector{String}
end

"""
    EDF.RecordAnnotation

Type containing all annotations applied to a particular data record.

# Fields

* `offset` (`Float64`): Offset from the recording start time (specified in the header)
  at which the current data record starts
* `event` (`Vector{String}` or `Nothing`): The event that marks the start of the data
  record, if applicable
* `annotations` (`Vector{EDF.AnnotationsList}`): The time-stamped annotations lists (TALs)
  in the current data record
* `n_bytes` (`Int`): The number of raw bytes per data record in the "EDF Annotation" signal
"""
mutable struct RecordAnnotation
    offset::Float64
    event::Union{Vector{String},Nothing}
    annotations::Vector{AnnotationsList}
    n_bytes::Int

    RecordAnnotation() = new()
    RecordAnnotation(offset, event, annotations, n_bytes) = new(offset, event, annotations, n_bytes)
end

"""
    EDF.Header

Type representing the header record for an EDF file.

## Fields

* `version` (`String`): Version of the data format
* `patient` (`String` or `EDF.PatientID`): Local patient identification
* `recording` (`String` or `EDF.RecordingID`): Local recording identification
* `start` (`DateTime`): Date and time the recording started
* `n_records` (`Int`): Number of data records
* `duration` (`Float64`): Duration of a data record in seconds
* `n_signals` (`Int`): Number of signals in a data record
* `nb_header` (`Int`): Total number of raw bytes in the header record
"""
struct Header
    version::String
    patient::Union{String,PatientID}
    recording::Union{String,RecordingID}
    continuous::Bool
    start::DateTime
    n_records::Int
    duration::Float64
    n_signals::Int
    nb_header::Int
end

# TODO: Make the vector of samples mmappable
# Also TODO: Refactor to make signals immutable
"""
    EDF.Signal

Type representing a single signal extracted from an EDF file.

# Fields

* `label` (`String`): The name of the signal, e.g. `F3-M2`
* `transducer` (`String`): Transducer type
* `physical_units` (`String`): Units of measure for the signal, e.g. `uV`
* `physical_min` (`Float32`): The physical minimum value of the signal
* `physical_max` (`Float32`): The physical maximum value of the signal
* `digital_min` (`Float32`): The minimum value of the signal that could occur in a data record
* `digital_max` (`Float32`): The maximum value of the signal that could occur in a data record
* `prefilter` (`String`): Description of any prefiltering done to the signal
* `n_samples` (`Int16`): The number of samples in a data record (NOT overall)
* `samples` (`Vector{Int16}`): The encoded sample values of the signal

!!! note
    Samples are stored in a `Signal` object in the same encoding as they appear in raw
    EDF files. See [`decode`](@ref) for decoding signals to their physical values.
"""
mutable struct Signal
    label::String
    transducer::String
    physical_units::String
    physical_min::Float32
    physical_max::Float32
    digital_min::Float32
    digital_max::Float32
    prefilter::String
    n_samples::Int16
    samples::Vector{Int16}

    Signal() = new()
end

"""
    EDF.File

Type representing a parsed EDF file.
All data defined in the file is accessible from this type by inspecting its fields
and the fields of the types of those fields.

# Fields

* `header` (`EDF.Header`): File-level metadata extracted from the file header
* `signals` (`Vector{EDF.Signal}`): All signals extracted from the data records
* `annotations` (`Vector{EDF.RecordAnnotation}` or `Nothing`): A vector of length
  `header.n_records` where each element contains annotations for the corresponding
  data record, if annotations are present in the file
"""
struct File
    header::Header
    signals::Vector{Signal}
    annotations::Union{Vector{RecordAnnotation},Nothing}
end

function Base.show(io::IO, edf::File)
    print(io, "EDF.File with ", length(edf.signals), " signals")
end

#####
##### Utilities
#####

"""
    EDF.decode(signal::Signal)

Decode the sample values in the given signal and return a `Vector` of the physical values.
"""
function decode(signal::Signal)
    digital_range = signal.digital_max - signal.digital_min
    physical_range = signal.physical_max - signal.physical_min
    return @. ((signal.samples - signal.digital_min) / digital_range) * physical_range + signal.physical_min
end

#####
##### The rest
#####

include("read.jl")
include("write.jl")

end # module


text:
#=
AdvancedPrimTest:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-10
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using Test

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/AdvancedPrim.jl")
matrix =  GBMatrix([[0.0, 7.0, 5.0, 5.0] [7.0, 0.0, 0.0, 0.0] [5.0, 0.0, 0.0, 0.0] [5.0, 0.0, 0.0, 0.0]])
m = GBVector([Inf, 0.0, 0.0, 0.0]) #source matrix
weight, mst = a_prim(matrix, 4, m)
@test weight == 17.0
@test mst[5] == 5
@test mst[6] == 5
@test mst[7] == 7