text:
##
cd(@__DIR__())
lines = readlines("input.txt")
tm = parse(Int, lines[1])
@show tm
ids = map(split(lines[2], ",")) do x
    x == "x" && return missing
    parse(Int128, x)
end |> skipmissing |> collect

##
ts = (tm .÷ ids .+ 1) .* ids
@show idn = ids[argmin(ts)]
@show wt = minimum(ts) - tm
@show idn * wt

## PARTE 2
lines = readlines("input.txt")
idsts = map(split(lines[2], ",")) do x
    x == "x" && return missing
    parse(Int128, x)
end |> enumerate
idsts = Iterators.filter(idsts) do (i,d)
    !ismissing(d)
end
ts = map(first, idsts) .- 1
ids = map(last, idsts)

##
n = 1068781
function find_tm(ids, ts; start = 1)
    n = start
    seqts = zeros(Int128, length(ids))
    ds = similar(ids)
    ds[1] = 0
    ds[2:end] .= 1
    @. seqts = ((n ÷ ids) + ds) * ids - ts
    while !all(==(seqts[1]), seqts[2:end])
        n += 1
        @. seqts = ((n ÷ ids) + ds) * ids - ts
    end
    n
end
function extended_gcd(a::T, b::T) where T
    old_r::T, r::T = a, b
    old_s::T, s::T = 1, 0
    old_t::T, t::T = 0, 1
    while r != 0
        q = old_r ÷ r
        old_r, r = r, old_r - q * r
        old_s, s = s, old_s - q * s
        old_t, t = t, old_t - q * t
    end
    (old_s, old_t), old_r, (t, s)
end
function cong_x(y::Union{Int128,BigInt}, p::Union{Int128,BigInt}, N::Integer = 1_000_000_000)
    c = zero(typeof(y))
    k = y ÷ p |> abs
    if y < 0
        y += k * p
    end
    while y < 0 && c <= N
        y += p
        c += 1
    end
    c = zero(typeof(y))
    k = y ÷ p |> abs
    if y > p
        y -= k * p
    end
    while y > p && y > 0 && c <= N
        y -= p
        c += 1
    end
    y
end
function find_tm_2(a::T, b, t1, t2, N::Integer = 1_000_000_000) where T
    (s, t), _, _ = extended_gcd(a, b)
    @show (s, t)
    x::T = t1 * t * b + t2 * s * a
    cong_x(x, a * b, N)
end
function find_tm_opt(ids, ts, N::Integer = 1_000_000_000)
    foldl(zip(ids[2:end], ts[2:end]); init=(ids[1],zero(Int128))) do (a, x), (b, t)
        @show (a, x), (b, t)
        a * b, find_tm_2(a, b, x, -t, N)
    end |> last
end

##
k = 9
@time x = find_tm_opt(ids[1:k], ts[1:k], 100_000_000_000)
@show x
@show (x .+ ts[1:k]) .% ids[1:k];


text:
export parseToAUTO

#Custom type to hold f90-code components
type f90_subroutine
    name::AbstractString
    args::AbstractString
    header::AbstractString
end

#Function that returns a sequence of all numerical arguments of the model
function declareNumericalArguments(M::Model)
    declare = "\tDOUBLE PRECISION "
    for var in keys(M.odes)
        declare *= var * ", "
    end
    for alg in keys(M.alg)
        declare *= alg * ", "
    end
    for par in keys(M.pars)
        declare *= par * ", "
    end
    declare *= "\n"
    return(declare)
end

#Parse algebraic equations
function parseAlgs(M::Model)
    parsed = ""
    for a in keys(M.alg)
        eq = M.alg[a]
        parsed *= "\t$a = $eq\n"
    end
    return(parsed)
end

#Parse Variables (declaration)
function parseVars(M::Model)
    parsed = ""
    for var in enumerate(keys(M.odes))
        var_name = var[2]
        var_i = var[1]
        parsed *= "\t$var_name = U($var_i)\n"
    end
    return(parsed)
end

#Parse singular point solutions for variables
function parseVarsPt(M::Model, values)
    parsed = ""
    for var in enumerate(keys(M.odes))
        var_i = var[1]
        var_val = values[var_i]
        parsed *= "\tU($var_i)=$var_val\n"
    end
    return(parsed)
end

#Parse parameters (declaration)
function parsePars(M::Model)
    parsed = ""
    for par in enumerate(keys(M.pars))
        par_name = par[2]
        par_i = par[1]
        parsed *= "\t$par_name = P($par_i)\n"
    end
    return(parsed)
end

#Parse parameter set corresponding to singular point solutions
function parseParsPt(M::Model)
    parsed = ""
    for par in enumerate(keys(M.pars))
        par_i = par[1]
        par_val = M.pars[par[2]]
        parsed *= "\tP($par_i)=$par_val\n"
    end
    return(parsed)
end

#Assemble f90 subroutines (generic)
function assembleSubroutine(routine::f90_subroutine, body::AbstractString)
    name = routine.name
    args = routine.args
    header = routine.header
    footer = "END SUBROUTINE $name"
    if header == ""
        R = "SUBROUTINE $name\n$footer\n"
    else
        R = "SUBROUTINE $name($args)\n$header\n$body\n$footer"
    end
    return(R)
end

#Assemble subroutines specifically
function parseSubroutine(routine_name, header_table, M::Model, values)
    if routine_name == "FUNC"
        body = "\n" * declareNumericalArguments(M) * "\n" * parseAlgs(M) * "\n" *  parseVars(M) * "\n" * parsePars(M)
    elseif routine_name == "STPNT"
        body = "\n" * parseParsPt(M) * "\n" *parseVarsPt(M, values)
    else
        body = ""
    end
    return(assembleSubroutine(header_table[routine_name], body))
end

#Top-level wrapper function for parsing model to f90-file
function parseToAUTO(M::Model, values, header_table = f90_header_table)
    parsed = "!--------------------------------------\n!\tGenerated with XPPjl\n!--------------------------------------\n\n"
    for routine_name in keys(header_table)
        parsed *= parseSubroutine(routine_name, header_table, M, values) * "\n\n"
    end
    model_name = M.name
    f = open(model_name * ".f90", "w")
    print(parsed)
    write(f, parsed)
    close(f)

    #Make parameter file
    pars = ""
    for p in keys(M.auto_specs)
        pars *= p * "= " * string(M.auto_specs[p]) * ", "
    end
    println(pars)
    f = open("c."*model_name, "w")
    write(f, pars)
    close(f)
end

f90_header_table= Dict(
    "FUNC" => f90_subroutine(
        "FUNC",
        "NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP",
        "\tIMPLICIT NONE
        INTEGER, INTENT(IN) :: NDIM, ICP(*), IJAC
        DOUBLE PRECISION, INTENT(IN) :: U(NDIM), PAR(*)
        DOUBLE PRECISION, INTENT(OUT) :: F(NDIM)
        DOUBLE PRECISION, INTENT(INOUT) :: DFDU(NDIM,NDIM), DFDP(NDIM,*)",
        ),
    "STPNT" => f90_subroutine(
        "STPNT",
        "NDIM,U,PAR,T",
        "\tIMPLICIT NONE
        INTEGER, INTENT(IN) :: NDIM
        DOUBLE PRECISION, INTENT(INOUT) :: U(NDIM),PAR(*)
        DOUBLE PRECISION, INTENT(IN) :: T"
        ),
    "BCND" => f90_subroutine("BCND", "", ""),
    "ICND" => f90_subroutine("ICND", "", ""),
    "FOPT" => f90_subroutine("FOPT", "", ""),
    "PVLS" => f90_subroutine("PVLS", "", ""),
    )

ext: jl

max_stars_repo_name: novakgroupoxford/XPPjl

max_issues_repo_name: novakgroupoxford/XPPjl

max_forks_repo_name: novakgroupoxford/XPPjl


text:
##
include("main_sim_parallel.jl")

"""
# ----------------------------------------------------------------------------------------------
## simulation with fixed parameters (single run)
# χ_M_0 = deg2rad(60)
# v_M_0 = V_M_0 * [sin(χ_M_0); cos(χ_M_0)]

# χ_f_d = deg2rad(180)
# s_BPNG.v̂_f_d = [cos(γ_f_d)*sin(χ_f_d); cos(γ_f_d)*cos(χ_f_d); sin(γ_f_d)]
# df, f = sim_plot(p_M_0, v_M_0, p_T_0, v_T_0, s_BPNG, A_M_max);
# display(f)

# ----------------------------------------------------------------------------------------------
## simulation with one parameter variation
# Variable 1
χ_M_0_list  = deg2rad.(30:5:90)
v_M_0_list  = χ_M_0_list |> Map(χ_M_0 -> V_M_0 * [sin(χ_M_0); cos(χ_M_0)]) |> collect

# Variable 2
χ_f_d_list  = deg2rad.(90:5:270)
v̂_f_d_list  = χ_f_d_list    |> Map(χ_f_d -> [cos(γ_f_d)*sin(χ_f_d); cos(γ_f_d)*cos(χ_f_d); sin(γ_f_d)])             |> collect
s_BPNG_list = v̂_f_d_list    |> Map(v̂_f_d -> BPNG(N, dim, σ_M_lim, v̂_f_d, Bias_IACG_StationaryTarget_2D, s_Bias))    |> collect

# size(df_list)

# df_list     = v_M_0_list |> Map(v_M_0 -> main( p_M_0, v_M_0, p_T_0, v_T_0, s_BPNG; A_M_max = A_M_max) ) |> tcollect
df_list     = s_BPNG_list   |> Map(s_BPNG -> main( p_M_0, v_M_0, p_T_0, v_T_0, s_BPNG; A_M_max = A_M_max))          |> tcollect

t_list      = df_list   |> Map(df -> df.time)     |> collect
# sol_list    = df_list   |> Map(df -> df.sol)      |> collect
x_list      = df_list   |> Map(df -> df.sol |> Map(datum -> datum.pursuer.p[1])     |> collect)    |> collect
y_list      = df_list   |> Map(df -> df.sol |> Map(datum -> datum.pursuer.p[2])     |> collect)    |> collect
χ_list      = df_list   |> Map(df -> df.sol |> Map(datum -> atan(datum.pursuer.v[1], datum.pursuer.v[2]) )     |> collect)    |> collect
A_list      = df_list   |> Map(df -> df.sol |> Map(datum -> norm(datum.pursuer.u))  |> collect)    |> collect
r_list      = df_list   |> Map(df -> df.sol |> Map(datum -> datum.r)                |> collect)    |> collect
σ_M_list    = df_list   |> Map(df -> df.sol |> Map(datum -> datum.σ_M)              |> collect)    |> collect
miss_list   = r_list    |> Map(r  -> r[end])    |> collect

f_xy = plot(size = (800,600))
for i in 1:length(t_list)
    plot!(f_xy, x_list[i], y_list[i], label = :false, xlabel = "x [m]", ylabel = "y [m]")
end
f_xy = fig_print([], [], "xy"; fig_handle = f_xy, ar_val = :equal) 	# To just save the result
display(f_xy)

f_r = plot()
for i in 1:length(t_list)
    plot!(f_r, t_list[i], r_list[i], label = :false, xlabel = "t [s]", ylabel = "r [m]")
end
f_r = fig_print([], [], "r"; fig_handle = f_r) 	# To just save the result
display(f_r)

# f_miss = fig_print(rad2deg.(χ_f_d_list), miss_list, "miss", :false, "χ_f_d [deg]", "r_f [m]")
# display(f_miss)

f_σ = plot()
for i in 1:length(t_list)
    plot!(f_σ, t_list[i], rad2deg.(σ_M_list[i]), label = :false, xlabel = "t [s]", ylabel = "σ_M [deg]")
end
f_σ = fig_print([], [], "sigma_M"; fig_handle = f_σ) 	# To just save the result
display(f_σ)

f_χ = plot()
for i in 1:length(t_list)
    plot!(f_χ, t_list[i], rad2deg.(χ_list[i]), label = :false, xlabel = "t [s]", ylabel = "χ_M [deg]")
    plot!(f_χ, t_list[i][[1,end],1], rad2deg.(χ_f_d_list[i]*ones(2)), label = :false, linestyle = :dot, linewidth = 0.5)
end
f_χ = fig_print([], [], "chi_M"; fig_handle = f_χ) 	# To just save the result
display(f_χ)

f_A = plot()
for i in 1:length(t_list)
    plot!(f_A, t_list[i], A_list[i], label = :false, xlabel = "t [s]", ylabel = "A_M [m/s^2]")
end
f_A = fig_print([], [], "A_M"; fig_handle = f_A) 	# To just save the result
display(f_A)

# ----------------------------------------------------------------------------------------------
## simulation with two parameters variation
df_list_list  = v_M_0_list |> Map( v_M_0 -> s_BPNG_list |> Map(s_BPNG -> main( p_M_0, v_M_0, p_T_0, v_T_0, s_BPNG; A_M_max = A_M_max)) |> tcollect  ) |> tcollect

t_list      = df_list_list |> Map( df_ic -> df_ic   |> Map(df_param -> df_param.time)      |> collect) |> collect
# sol_list    = df_list_list |> Map( df_ic -> df_ic   |> Map(df_param -> df_param.sol)       |> collect) |> collect

x_list      = df_list_list |> Map( df_ic -> df_ic  |> Map(df_param -> df_param.sol |> Map(datum -> datum.pursuer.p[1])                              |> collect  )  |> collect )    |> collect
y_list      = df_list_list |> Map( df_ic -> df_ic  |> Map(df_param -> df_param.sol |> Map(datum -> datum.pursuer.p[2])                              |> collect  )  |> collect )    |> collect
χ_list      = df_list_list |> Map( df_ic -> df_ic  |> Map(df_param -> df_param.sol |> Map(datum -> atan(datum.pursuer.v[1], datum.pursuer.v[2]) )   |> collect  )  |> collect )    |> collect
A_list      = df_list_list |> Map( df_ic -> df_ic  |> Map(df_param -> df_param.sol |> Map(datum -> norm(datum.pursuer.u))                           |> collect  )  |> collect )    |> collect
r_list      = df_list_list |> Map( df_ic -> df_ic  |> Map(df_param -> df_param.sol |> Map(datum -> datum.r)                                         |> collect  )  |> collect )    |> collect
σ_M_list    = df_list_list |> Map( df_ic -> df_ic  |> Map(df_param -> df_param.sol |> Map(datum -> datum.σ_M)                                       |> collect  )  |> collect )    |> collect
miss_list   = r_list       |> Map( r_ic_list -> r_ic_list   |> Map(r  -> r[end])    |> collect ) |> collect

miss_list   = hcat(miss_list...)'

#
f_miss_surf = contour(rad2deg.(χ_M_0_list), rad2deg.(χ_f_d_list), miss_list', xlabel = "χ_M_0 [deg]", ylabel = "χ_f_d [deg]", fill = true)
# plot(rad2deg.(χ_M_0_list), rad2deg.(χ_f_d_list), miss_list, st = [:surface, :contourf], layout = 2)
f_miss_surf = fig_print([], [], "miss"; fig_handle = f_miss_surf) 	# To just save the result

i = 2 # in 1:length(v_M_0_list)

for i in 1:length(v_M_0_list)
    f_xy = plot(size = (800,600))

    for j in 1:length(s_BPNG_list)
        plot!(f_xy, x_list[i][j], y_list[i][j], label = :false, xlabel = "x [m]", ylabel = "y [m]")
    end
    f_xy = fig_print([], [], "xy_ic_"*string(i); fig_handle = f_xy, ar_val = :equal) 	# To just save the result
    display(f_xy)
end

# f_r = plot()
# for j in 1:length(s_BPNG_list)
#     plot!(f_r, t_list[i][j], r_list[i][j], label = :false, xlabel = "t [s]", ylabel = "r [m]")
# end
# f_r = fig_print([], [], "r_ic_"*string(i); fig_handle = f_r) 	# To just save the result
# display(f_r)

# f_σ = plot()
# for j in 1:length(s_BPNG_list)
#     plot!(f_σ, t_list[i][j], rad2deg.(σ_M_list[i][j]), label = :false, xlabel = "t [s]", ylabel = "σ_M [deg]")
# end
# f_σ = fig_print([], [], "sigma_M_ic_"*string(i); fig_handle = f_σ) 	# To just save the result
# display(f_σ)

# f_χ = plot()
# for j in 1:length(s_BPNG_list)
#     plot!(f_χ, t_list[i][j], rad2deg.(χ_list[i][j]), label = :false, xlabel = "t [s]", ylabel = "χ_M [deg]")
#     plot!(f_χ, t_list[i][j][[1,end],1], rad2deg.(χ_f_d_list[j]*ones(2)), label = :false, linestyle = :dot, linewidth = 0.5)
# end
# f_χ = fig_print([], [], "chi_M_ic_"*string(i); fig_handle = f_χ) 	# To just save the result
# display(f_χ)

# f_A = plot()
# for j in 1:length(s_BPNG_list)
#     plot!(f_A, t_list[i][j], A_list[i][j], label = :false, xlabel = "t [s]", ylabel = "A_M [m/s^2]")
# end
# f_A = fig_print([], [], "A_M_ic_"*string(i); fig_handle = f_A) 	# To just save the result
# display(f_A)

# ---------------------------------------------------------------------------------------------------
## 3D guidance
include("main_missiles.jl")
include("sim_missiles.jl")

##
# Parallel execution (works only in REPL mode)
# Threads.nthreads()

dim = 3

γ_f_d_list  = deg2rad.(-90:5:0)
χ_f_d_list  = deg2rad.(90:5:270)
v̂_f_d_list  = γ_f_d_list    |> Map( γ_f_d -> χ_f_d_list |> Map(χ_f_d -> [cos(γ_f_d)*sin(χ_f_d); cos(γ_f_d)*cos(χ_f_d); sin(γ_f_d)])             |> collect ) |> collect
s_BPNG_list = v̂_f_d_list    |> Map( v̂_f_d_γ -> v̂_f_d_γ  |> Map(v̂_f_d -> BPNG(N, dim, σ_M_lim, v̂_f_d, Bias_IACG_StationaryTarget, s_Bias))       |> collect ) |> collect

df_list_list  = s_BPNG_list |> Map( s_BPNG_γ -> s_BPNG_γ |> Map(s_BPNG -> main( p_M_0, v_M_0, p_T_0, v_T_0, s_BPNG; A_M_max = A_M_max)) |> tcollect  ) |> tcollect

x_list      = df_list_list |> Map( df_ic -> df_ic  |> Map(df_param -> df_param.sol |> Map(datum -> datum.pursuer.p[1])  |> collect  )  |> collect )    |> collect
y_list      = df_list_list |> Map( df_ic -> df_ic  |> Map(df_param -> df_param.sol |> Map(datum -> datum.pursuer.p[2])  |> collect  )  |> collect )    |> collect
z_list      = df_list_list |> Map( df_ic -> df_ic  |> Map(df_param -> df_param.sol |> Map(datum -> datum.pursuer.p[3])  |> collect  )  |> collect )    |> collect

f_3D = plot(aspect_ratio = :equal)
for i in 1:length(γ_f_d_list)
    for j in 1:length(χ_f_d_list)
        plot!(f_3D, x_list[i][j], y_list[i][j], z_list[i][j], label=:false)
    end
end
plot!(f_3D, xlabel = "x [m]", ylabel = "y [m]", zlabel = "z [m]", camera = (45,30), size = (800,800),
xlims = (0,5.5E3), ylims = (-1E3, 4.5E3), zlims = (-1E3,4E3))
f_3D = fig_print([], [], "Trajectory_3D"; fig_handle = f_3D, ar_val = :equal) 	# To just save the result

"""

## ---------------------------------------------------------------------------------------------------
# massive numerical analysis example
# ---------------------------------------------------------------------------------------------------
using LaTeXStrings
include("main_sim_parallel.jl")

# ---------------------------------------------------------------------------------------------------
## List of parameters
# Variable 1
γ_xref_M_0_list = deg2rad.(0:5:60)
χ_M_0_list  = pi / 2 .- γ_xref_M_0_list # deg2rad.(30:5:90)
v_M_0_list  = χ_M_0_list |> Map(χ_M_0 -> V_M_0 * [sin(χ_M_0); cos(χ_M_0)]) |> collect

# Variable 2
γ_xref_f_d_list = deg2rad.(-180:5:0)
χ_f_d_list  = pi / 2 .- γ_xref_f_d_list # deg2rad.(90:5:270)
v̂_f_d_list  = χ_f_d_list    |> Map(χ_f_d -> [cos(γ_f_d) * sin(χ_f_d); cos(γ_f_d) * cos(χ_f_d); sin(γ_f_d)])             |> collect
# s_BPNG_list = v̂_f_d_list    |> Map(v̂_f_d -> BPNG(N, dim, σ_M_lim, v̂_f_d, Bias_IACG_StationaryTarget_2D, s_Bias))    |> collect

# Variable 3-4
α_list      = 0.6:0.2:1
n_list      = 0:1:3
# s_BPNG_list   = α_list      |> Map( α -> n_list  |> Map( n -> begin
#                     s_Bias  = ComponentArray(α = α, δ = 0.01, n = n, r_ref = 10E3, k = 4, m = 10);
#                     BPNG(N, dim, σ_M_lim, v̂_f_d, Bias_IACG_StationaryTarget, s_Bias)
#                 end) |> collect ) |> collect

main(p_M_0, v_M_0, p_T_0, v_T_0, s_BPNG; A_M_max=A_M_max)
# --------> the above line is to execute main() once for compilation before passing to multiple threads

# ---------------------------------------------------------------------------------------------------
## simulation for paper: gain design parameter variation
mutable struct Results
    α
    n
    r_f_result
    r_min_result
end
s_sim = Array{Results}(undef, length(α_list), length(n_list))

for i_α in 1:length(α_list)
    for i_n in 1:length(n_list)
        α       = α_list[i_α]
        n       = n_list[i_n]
        s_Bias  = ComponentArray(α=α, δ=0.01, n=n, r_ref=5E3, k=4, m=10)
        s_BPNG_list = v̂_f_d_list    |> Map(v̂_f_d -> BPNG(N, dim, σ_M_lim, v̂_f_d, Bias_IACG_StationaryTarget_2D, s_Bias))    |> collect

        # execution
        df_list_list  = v_M_0_list |> Map(v_M_0 -> s_BPNG_list |> Map(s_BPNG -> main(p_M_0, v_M_0, p_T_0, v_T_0, s_BPNG; A_M_max=A_M_max)) |> tcollect) |> tcollect

        # data collection
        t_list      = df_list_list |> Map(df_ic -> df_ic   |> Map(df_param -> df_param.time)      |> collect) |> collect
        # sol_list    = df_list_list |> Map( df_ic -> df_ic   |> Map(df_param -> df_param.sol)       |> collect) |> collect
        x_list      = df_list_list |> Map(df_ic -> df_ic  |> Map(df_param -> df_param.sol  |> Map(datum -> datum.pursuer.p[1])                              |> collect)  |> collect)    |> collect
        y_list      = df_list_list |> Map(df_ic -> df_ic  |> Map(df_param -> df_param.sol  |> Map(datum -> datum.pursuer.p[2])                              |> collect)  |> collect)    |> collect
        χ_list      = df_list_list |> Map(df_ic -> df_ic  |> Map(df_param -> df_param.sol  |> Map(datum -> atan(datum.pursuer.v[1], datum.pursuer.v[2]))    |> collect)  |> collect)    |> collect
        γ_xref_list = df_list_list |> Map(df_ic -> df_ic  |> Map(df_param -> df_param.sol  |> Map(datum -> atan(datum.pursuer.v[2], datum.pursuer.v[1]))    |> collect)  |> collect)    |> collect
        λ_xref_list = df_list_list |> Map(df_ic -> df_ic  |> Map(df_param -> df_param.sol  |> Map(datum -> atan(datum.evador.p[2] - datum.pursuer.p[2], datum.evador.p[1] - datum.pursuer.p[1]))    |> collect)  |> collect)    |> collect
        A_list      = df_list_list |> Map(df_ic -> df_ic  |> Map(df_param -> df_param.sol  |> Map(datum -> norm(datum.pursuer.u))                           |> collect)  |> collect)    |> collect
        r_list      = df_list_list |> Map(df_ic -> df_ic  |> Map(df_param -> df_param.sol  |> Map(datum -> datum.r)                                         |> collect)  |> collect)    |> collect
        σ_M_list    = df_list_list |> Map(df_ic -> df_ic  |> Map(df_param -> df_param.sol  |> Map(datum -> datum.σ_M)                                       |> collect)  |> collect)    |> collect

        r_f_result   = r_list       |> Map(r_ic_list -> r_ic_list   |> Map(r  -> r[end])    |> collect) |> collect
        r_min_result = r_list       |> Map(r_ic_list -> r_ic_list   |> Map(r  -> minimum(r)) |> collect) |> collect

        σ_M_signed_list  = γ_xref_list .- λ_xref_list
        γ_xref_pred_list = γ_xref_list .- N / (N - 1) .* σ_M_signed_list

        r_f_result      = hcat(r_f_result...)'
        r_min_result    = hcat(r_min_result...)'
        s_sim[i_α, i_n] = Results(α, n, r_f_result, r_min_result)

        γ_xref_M_0_list = 90 .- rad2deg.(χ_M_0_list)
        γ_xref_f_d_list = 90 .- rad2deg.(χ_f_d_list)
        K_r_max         = (N - 1 + 0.01) * 2^n
        γ_xref_f_d_ub = min(rad2deg.(((N - 1) * tan(σ_M_lim) / K_r_max)^(1 / α)) .- γ_xref_M_0_list / (N - 1), zeros(length(γ_xref_M_0_list)))
        γ_xref_f_d_lb =      rad2deg.(-((N - 1) * tan(σ_M_lim) / K_r_max)^(1 / α)) .- γ_xref_M_0_list / (N - 1)

        # r_f
        f_r_f = contour(γ_xref_M_0_list, γ_xref_f_d_list, r_f_result',
                        xlabel=L"\gamma_0 \textrm{~[deg]}", ylabel=L"\gamma_f_d \textrm{~[deg]}", fill=true, seriescolor=:coolwarm, ylims=extrema(γ_xref_f_d_list))
        plot!(f_r_f, γ_xref_M_0_list, [γ_xref_f_d_lb γ_xref_f_d_ub], linecolor=:red, linestyle=:dash, linewidth=2,
                    label=["LB" "UB"], ylims=extrema(γ_xref_f_d_list)) # title="r_f: α = $α, n = $n"
        f_r_f = fig_print([], [], "r_f_alpha_$(i_α)_n_$(i_n)", [], [], [], f_r_f) 	# To just save the result
        display(f_r_f)

        # r_min
        f_r_min = contour(γ_xref_M_0_list, γ_xref_f_d_list, r_min_result',
        xlabel=L"\gamma_0 \textrm{~[deg]}", ylabel=L"\gamma_f_d \textrm{~[deg]}", fill=true, seriescolor=:coolwarm, ylims=extrema(γ_xref_f_d_list))
        plot!(f_r_min, γ_xref_M_0_list, [γ_xref_f_d_lb γ_xref_f_d_ub], linecolor=:red, linestyle=:dash, linewidth=2,
                    label=["LB" "UB"], ylims=extrema(γ_xref_f_d_list)) # title="r_f: α = $α, n = $n"
        f_r_min = fig_print([], [], "r_min_alpha_$(i_α)_n_$(i_n)", [], [], [], f_r_min) 	# To just save the result
        display(f_r_min)

        # OX
        f_OX = plot()
        for i in 1:length(γ_xref_M_0_list)
            for j in 1:length(γ_xref_f_d_list)
                if r_min_result[i,j] < 4 && r_f_result[i,j] < 4
                    plot!(f_OX, γ_xref_M_0_list[i] * ones(2), γ_xref_f_d_list[j] * ones(2), markershape=:circle, markercolor=:blue, label=:false)
                else
                    plot!(f_OX, γ_xref_M_0_list[i] * ones(2), γ_xref_f_d_list[j] * ones(2), markershape=:xcross, markercolor=:red, label=:false)
                end
            end
        end
        plot!(f_OX, xlabel=L"\gamma_0 \textrm{~[deg]}", ylabel=L"\gamma_f_d \textrm{~[deg]}")
        f_OX = fig_print([], [], "OX_alpha_$(i_α)_n_$(i_n)", [], [], [], f_OX) 	# To just save the result
        display(f_OX)

        # e_γ_f
        for i in 1:length(γ_xref_M_0_list)
            f_e_γ_f = plot!()

            for j in 1:length(γ_xref_f_d_list)
                γ_xref_f_d = γ_xref_f_d_list[j]
                e_γ_f_list = 180 / pi .* γ_xref_pred_list[i][j] .- γ_xref_f_d
                plot!(f_e_γ_f, t_list[i][j], abs.(e_γ_f_list), label=:false, xlabel="t [s]", ylabel="|e_γ_f| [deg]", ylims=(0, 60))
            end
            f_e_γ_f = fig_print([], [], "e_gamma_f_alpha_$(i_α)_n_$(i_n)_ic_$i", [], [], [], f_e_γ_f) 	# To just save the result
            # display(f_e_γ_f)
        end

        # xy
        for i in 1:length(v_M_0_list)
            f_xy = plot(size=(800, 600))

            for j in 1:length(s_BPNG_list)
                plot!(f_xy, x_list[i][j], y_list[i][j], label=:false, xlabel="x [m]", ylabel="y [m]", ylims=(-500, 4.5E3))
            end
            f_xy = fig_print([], [], "xy_alpha_$(i_α)_n_$(i_n)_ic_$i", [], [], [], f_xy, ar_val=:equal) 	# To just save the result
            # display(f_xy)
        end

        # f_r = plot()
        # for j in 1:length(s_BPNG_list)
        #     plot!(f_r, t_list[i][j], r_list[i][j], label = :false, xlabel = "t [s]", ylabel = "r [m]")
        # end
        # f_r = fig_print([], [], "r", [], [], [], f_r) 	# To just save the result
        # display(f_r)

        # f_σ = plot()
        # for j in 1:length(s_BPNG_list)
        #     plot!(f_σ, t_list[i][j], rad2deg.(σ_M_list[i][j]), label = :false, xlabel = "t [s]", ylabel = "σ_M [deg]")
        # end
        # f_σ = fig_print([], [], "sigma_M", [], [], [], f_σ) 	# To just save the result
        # display(f_σ)

        # f_χ = plot()
        # for j in 1:length(s_BPNG_list)
        #     plot!(f_χ, t_list[i][j], rad2deg.(χ_list[i][j]), label = :false, xlabel = "t [s]", ylabel = "χ_M [deg]")
        #     plot!(f_χ, t_list[i][j][[1,end],1], rad2deg.(χ_f_d_list[j]*ones(2)), label = :false, linestyle = :dot, linewidth = 0.5)
        # end
        # f_χ = fig_print([], [], "chi_M", [], [], [], f_χ) 	# To just save the result
        # display(f_χ)

        # f_A = plot()
        # for j in 1:length(s_BPNG_list)
        #     plot!(f_A, t_list[i][j], A_list[i][j], label = :false, xlabel = "t [s]", ylabel = "A_M [m/s^2]")
        # end
        # f_A = fig_print([], [], "A_M", [], [], [], f_A) 	# To just save the result
        # display(f_A)

    end
end


text:
function logistic_map(r)
    x = [0.5]

    # Your code goes here!

    return x
end

max_stars_repo_name: basimr/lovelace

max_issues_repo_name: basimr/project-lovelace

max_forks_repo_name: basimr/project-lovelace


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -3.892086739821988 <= q <= 1.133112254487742)

                     add_NL_constraint(m, :(swish(-0.9215225812660202 + -0.8659264200816277*swish(0.9754315885766114 + 0.9182647481159649*$(x[1]) + 0.2963572993024828*$(x[2])) + -0.7655479829179535*swish(0.7516798050016731 + -0.4187259277447679*$(x[1]) + -0.991762240929563*$(x[2]))) + swish(0.9152212659843197 + -0.09622429015475253*swish(0.9754315885766114 + 0.9182647481159649*$(x[1]) + 0.2963572993024828*$(x[2])) + 0.18727818139193309*swish(0.7516798050016731 + -0.4187259277447679*$(x[1]) + -0.991762240929563*$(x[2]))) - $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:
# TODO Add VPRKpSecondary tableau.

"Parameters for right-hand side function of variational partitioned Runge-Kutta methods."
mutable struct ParametersVPRKpSecondary{DT, TT, ET <: VODE{DT,TT}, D, S} <: AbstractParametersVPRK{DT,TT,ET,D,S}
    equ::ET
    tab::TableauVPRK{TT}
    Δt::TT

    ω::Matrix{TT}
    R::Vector{TT}

    t::TT
    q::Vector{DT}
    p::Vector{DT}
end

function ParametersVPRKpSecondary(equ::ET, tab::TableauVPRK{TT}, Δt::TT) where {DT, TT, ET <: VODE{DT,TT}}
    D = equ.d
    S = tab.s

    # compute reduction matrix
    ω = zeros(TT, S-1, S)

    for i in 1:(S-1)
        for j in 1:S
            ω[i,j] = tab.q.b[j] * tab.q.c[j]^(i-1)
        end
    end

    R = convert(Vector{TT}, [1, tab.R∞])
    q = zeros(DT,D)
    p = zeros(DT,D)

    ParametersVPRKpSecondary{DT,TT,ET,D,S}(equ, tab, Δt, ω, R, 0, q, p)
end

function compute_stages_vprk!(x, q̅, p̅, Q, V, Λ, P, F, R, Φ, params)
    # copy x to V
    compute_stages_v_vprk!(x, V, params)

    # copy x to Λ
    compute_stages_λ_vprk!(x, Λ, params)

    # compute Q
    compute_stages_q_vprk!(q̅, Q, V, Λ, params)

    # compute p̅, R and Ψ
    compute_projection_vprk!(q̅, p̅, Q, V, Λ, R, Φ, params)

    # compute P and F
    compute_stages_p_vprk!(Q, V, P, F, params)
end

function compute_stages_q_vprk!(q̅::Vector{ST}, Q::Matrix{ST}, V::Matrix{ST}, Λ::Matrix{ST},
                params::ParametersVPRKpSecondary{DT,TT,ET,D,S}) where {ST,DT,TT,ET,D,S}

    @assert D == size(Q,1) == size(V,1) == size(Λ,1) == length(q̅)
    @assert S == size(Q,2) == size(V,2) == size(Λ,2)

    local y1::ST
    local y2::ST
    local y3::ST
    local y4::ST

    # compute Q
    for i in 1:S
        for k in 1:D
            y1 = y2 = y3 = y4 = 0
            for j in 1:S
                y1 += params.tab.q.a[i,j] * V[k,j]
                y2 += params.tab.q.â[i,j] * V[k,j]
                y3 += params.tab.q.a[i,j] * Λ[k,j]
                y4 += params.tab.q.â[i,j] * Λ[k,j]
            end
            Q[k,i] = params.q[k] + params.Δt * (y1 + y2 + y3 + y4)
        end
    end

    # compute q̅
    for k in 1:D
        y1 = y2 = y3 = y4 = 0
        for j in 1:S
            y1 += params.tab.q.b[j] * V[k,j]
            y2 += params.tab.q.b̂[j] * V[k,j]
            y3 += params.tab.q.b[j] * Λ[k,j]
            y4 += params.tab.q.b̂[j] * Λ[k,j]
        end
        q̅[k] = params.q[k] + params.Δt * (y1 + y2 + y3 + y4)
    end
end

@generated function compute_projection_vprk!(q̅::Vector{ST}, p̅::Vector{ST},
                Q::Matrix{ST}, V::Matrix{ST}, Λ::Matrix{ST}, R::Matrix{ST}, Ψ::Matrix{ST},
                params::ParametersVPRKpSecondary{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    # create temporary vectors
    tQ = zeros(ST,D)
    tV = zeros(ST,D)
    tΛ = zeros(ST,D)
    tR = zeros(ST,D)
    tΨ = zeros(ST,D)
    tH = zeros(ST,D)
    tΩ = zeros(ST,D,D)

    quote
        local t₀::TT = params.t
        local t₁::TT = params.t + params.Δt
        local tᵢ::TT
        local v = zeros(q̅)

        # compute p̅=ϑ(q̅)
        params.equ.ϑ(t₁, q̅, v, p̅)

        for i in 1:S
            simd_copy_xy_first!($tQ, Q, i)
            simd_copy_xy_first!($tV, V, i)
            simd_copy_xy_first!($tΛ, Λ, i)

            tᵢ = t₀ + params.Δt * params.tab.p.c[i]

            params.equ.g(tᵢ, $tQ, $tΛ, $tR)
            params.equ.ω(tᵢ, $tQ, $tΩ)
            params.equ.dH(tᵢ, $tQ, $tH)

            # TODO Check if ω() returns Ω or Ω^T -> this decides the sign on dH below
            #      Seems to be Ω^T -> then dH should be subtracted

            simd_mult!($tΨ, $tΩ, $tV)

            $tΨ .-= $tH

            simd_copy_yx_first!($tR, R, i)
            simd_copy_yx_first!($tΨ, Ψ, i)
        end
    end
end

function compute_rhs_vprk!(b::Vector{ST}, P::Matrix{ST}, F::Matrix{ST}, R::Matrix{ST},
                params::ParametersVPRKpSecondary{DT,TT,ET,D,S}) where {ST,DT,TT,ET,D,S}

    local z1::ST
    local z2::ST
    local z3::ST
    local z4::ST

    # compute b = - [(P-AF-AR)]
    for i in 1:S
        for k in 1:D
            z1 = z2 = z3 = z4 = 0
            for j in 1:S
                z1 += params.tab.p.a[i,j] * F[k,j]
                z2 += params.tab.p.â[i,j] * F[k,j]
                z3 += params.tab.p.a[i,j] * R[k,j]
                z4 += params.tab.p.â[i,j] * R[k,j]
            end
            b[D*(i-1)+k] = (P[k,i] - params.p[k]) - params.Δt * (z1 + z2 + z3 + z4)
        end
    end
end

function compute_rhs_vprk_projection!(b::Vector{ST}, p̅::Vector{ST},
                F::Matrix{ST}, R::Matrix{ST}, Ψ::Matrix{ST}, offset::Int,
                params::ParametersVPRKpSecondary{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    local ψ::ST
    local z1::ST
    local z2::ST
    local z3::ST
    local z4::ST

    # for i in 1:S
    #     for k in 1:D
    #         b[offset+D*(i-1)+k] = Ψ[k,i]
    #     end
    # end

    for i in 1:(S-1)
        for k in 1:D
            ψ = 0
            for j in 1:S
                ψ += params.ω[i,j] * Ψ[k,j]
            end
            b[offset+D*(i-1)+k] = ψ
        end
    end

    for k in 1:D
        z1 = z2 = z3 = z4 = 0
        for j in 1:S
            z1 += params.tab.p.b[j] * F[k,j]
            z2 += params.tab.p.b̂[j] * F[k,j]
            z3 += params.tab.p.b[j] * R[k,j]
            z4 += params.tab.p.b̂[j] * R[k,j]
        end
        b[offset+D*(S-1)+k] = (p̅[k] - params.p[k]) - params.Δt * (z1 + z2 + z3 + z4)
    end
end

"Compute stages of variational partitioned Runge-Kutta methods."
@generated function function_stages!(x::Vector{ST}, b::Vector{ST},
                params::ParametersVPRKpSecondary{DT,TT,ET,D,S}
            ) where {ST,DT,TT,ET,D,S}

    scache = NonlinearFunctionCacheVPRK{ST}(D,S)
    pcache = NonlinearFunctionCacheVPRKprojection{ST}(D,S)

    function_stages = quote
        compute_stages_vprk!(x, $pcache.q̅, $pcache.p̅,
                                $scache.Q, $scache.V, $pcache.Λ,
                                $scache.P, $scache.F, $pcache.R,
                                $pcache.Φ, params)

        # compute b = [P-AF-AR]
        compute_rhs_vprk!(b, $scache.P, $scache.F, $pcache.R, params)

        # compute b = Φ
        compute_rhs_vprk_projection!(b, $pcache.p̅, $scache.F, $pcache.R, $pcache.Φ, D*S, params)

        compute_rhs_vprk_correction!(b, $scache.V, params)
    end

    return function_stages
end

@doc raw"""
Variational partitioned Runge-Kutta integrator with projection on secondary constraint.

The VPRK integrator solves the following system of equations for the internal stages,
```math
\begin{align*}
Q_{n,i} &= q_{n} + h \sum \limits_{j=1}^{s} a_{ij} \, \big( V_{n,j} + \Lambda_{n,j} \big) , \\
P_{n,i} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{a}_{ij} \, \big( F_{n,j} + \nabla \vartheta (Q_{n,j})^{T} \cdot \Lambda_{n,j} \big) - d_i \lambda , \\
0 &= \sum \limits_{i=1}^{s} d_i V_i , \\
0 &= \sum \limits_{j=1}^{s} \omega_{ij} \Psi_{n,j} ,
\end{align*}
```
with definitions
```math
\begin{align*}
P_{n,i} &= \dfrac{\partial L}{\partial v} (Q_{n,i}, V_{n,i}) , \\
F_{n,i} &= \dfrac{\partial L}{\partial q} (Q_{n,i}, V_{n,i}) , \\
\Psi_{n,i} &= \psi(Q_{n,i}, V_{n,i}, P_{n,i}, F_{n,i}) ,
\end{align*}
```
and update rule
```math
\begin{align*}
q_{n+1} &= q_{n} + h \sum \limits_{i=1}^{s} b_{i} \, \big( V_{n,i} + \Lambda_{n,i} \big) , \\
p_{n+1} &= p_{n} + h \sum \limits_{i=1}^{s} \bar{b}_{i} \, \big( F_{n,i} + \nabla \vartheta (Q_{n,j})^{T} \cdot \Lambda_{n,j} \big) , \\
0 &= \phi (q_{n+1}, p_{n+1}) ,
\end{align*}
```
satisfying the symplecticity conditions
```math
\begin{align*}
b_{i} \bar{a}_{ij} + \bar{b}_{j} a_{ji} &= b_{i} \bar{b}_{j} , &
\bar{b}_i &= b_i ,
\end{align*}
```
the primary constraint,
```math
\begin{align*}
\phi(q,p) = p - \vartheta (q) = 0 ,
\end{align*}
```
at the final solution ``(q_{n+1}, p_{n+1})``,
and super positions of the secondary constraints,
```math
\begin{align*}
\psi(q,\dot{q},p,\dot{p})
= \dot{p} - \dot{q} \cdot \nabla \vartheta (q)
= \big( \nabla \vartheta (q) - \nabla \vartheta^{T} (q) \big) \cdot \dot{q} - \nabla H (q)
= 0,
\end{align*}
```
which, evaluated at the internal stages, read
```math
\begin{align*}
\Psi_{n,j} = \big( \nabla \vartheta (Q_{n,j}) - \nabla \vartheta^{T} (Q_{n,j}) \big) \cdot V_{n,j} - \nabla H (Q_{n,j}) .
\end{align*}
```
Here, ``\omega`` is a ``(s-1) \times s`` matrix, chosen such that the resulting
method has optimal order.
The vector ``d`` is zero for Gauss-Legendre methods and needs to be chosen
appropriately for Gauss-Lobatto methods (for details see documentation of
VPRK methods).
"""
struct IntegratorVPRKpSecondary{DT, TT, PT <: ParametersVPRKpSecondary{DT,TT},
                                        ST <: NonlinearSolver{DT},
                                        IT <: InitialGuessPODE{DT,TT}} <: AbstractIntegratorVPRK{DT,TT}
    params::PT
    solver::ST
    iguess::IT

    scache::NonlinearFunctionCacheVPRK{DT}
    pcache::NonlinearFunctionCacheVPRKprojection{DT}

    q::Vector{Vector{TwicePrecision{DT}}}
    p::Vector{Vector{TwicePrecision{DT}}}
end

function IntegratorVPRKpSecondary(equation::ET, tableau::TableauVPRK{TT}, Δt::TT) where {DT, TT, ET <: VODE{DT,TT}}
    D = equation.d
    M = equation.n
    S = tableau.s

    # create params
    params = ParametersVPRKpSecondary(equation, tableau, Δt)

    # create solver
    solver = create_nonlinear_solver(DT, 2*D*S, params)

    # create initial guess
    iguess = InitialGuessPODE(get_config(:ig_interpolation), equation, Δt)

    # create cache for internal stage vectors and update vectors
    scache = NonlinearFunctionCacheVPRK{DT}(D,S)
    pcache = NonlinearFunctionCacheVPRKprojection{DT}(D,S)

    # create solution vectors
    q = create_solution_vector(DT, D, M)
    p = create_solution_vector(DT, D, M)

    # create integrator
    IntegratorVPRKpSecondary{DT, TT, typeof(params), typeof(solver), typeof(iguess)}(
                params, solver, iguess, scache, pcache, q, p)
end

function initial_guess!(int::IntegratorVPRKpSecondary{DT,TT}, m::Int) where {DT,TT}
    for i in 1:int.params.tab.s
        evaluate!(int.iguess, m, int.scache.y, int.scache.z, int.scache.v, int.params.tab.q.c[i], int.params.tab.p.c[i])
        for k in 1:int.params.equ.d
            int.solver.x[int.params.equ.d*(0*int.params.tab.s+i-1)+k] = int.scache.v[k]
            int.solver.x[int.params.equ.d*(1*int.params.tab.s+i-1)+k] = 0
        end
    end
end

"Integrate ODE with variational partitioned Runge-Kutta integrator."
function integrate_step!(int::IntegratorVPRKpSecondary{DT,TT}, sol::SolutionPDAE{DT,TT}, m::Int, n::Int) where {DT,TT}
    # check if m and n are compatible with solution dimensions
    check_solution_dimension_asserts(sol, m, n)

    # set time and solution for nonlinear solver
    int.params.t = sol.t[0] + (n-1)*int.params.Δt
    int.params.q .= int.q[m]
    int.params.p .= int.p[m]

    # compute initial guess
    initial_guess!(int, m)

    # call nonlinear solver
    solve!(int.solver)

    # print solver status
    print_solver_status(int.solver.status, int.solver.params, n)

    # check if solution contains NaNs or error bounds are violated
    check_solver_status(int.solver.status, int.solver.params, n)

    # compute final update
    compute_stages_vprk!(int.solver.x, int.pcache.q̅, int.pcache.p̅,
                         int.scache.Q, int.scache.V, int.pcache.Λ,
                         int.scache.P, int.scache.F, int.pcache.R,
                         int.pcache.Φ, int.params)

    # compute unprojected solution
    update_solution!(int.q[m], int.scache.V, int.params.tab.q.b, int.params.tab.q.b̂, int.params.Δt)
    update_solution!(int.p[m], int.scache.F, int.params.tab.p.b, int.params.tab.p.b̂, int.params.Δt)

    # add projection to solution
    update_solution!(int.q[m], int.pcache.Λ, int.params.tab.q.b, int.params.tab.q.b̂, int.params.Δt)
    update_solution!(int.p[m], int.pcache.R, int.params.tab.p.b, int.params.tab.p.b̂, int.params.Δt)

    # copy solution to initial guess for next time step
    update!(int.iguess, m, sol.t[0] + n*int.params.Δt, int.q[m], int.p[m])

    # take care of periodic solutions
    cut_periodic_solution!(int.q[m], int.params.equ.periodicity)

    # copy to solution
    copy_solution!(sol, int.q[m], int.p[m], int.pcache.λ, n, m)
end


text:
""" DSLs for defining categories, diagrams, and related structures.

Here "diagram" means diagram in the ordinary sense of category theory, not
wiring diagram.
"""
module DiagrammaticPrograms
export @graph, @category, @diagram

using Base.Iterators: repeated
using MLStyle: @match

using ...Present, ...Graphs, ...CategoricalAlgebra
using ...Graphs.BasicGraphs: TheoryGraph

# Data types
############

@present TheoryNamedGraph <: TheoryGraph begin
  Name::AttrType
  vname::Attr(V, Name)
  ename::Attr(E, Name)
end

""" Default graph type for [`@category`](@ref) macro and related macros.
"""
@acset_type NamedGraph(TheoryNamedGraph, index=[:src,:tgt],
                       unique_index=[:vname,:ename]) <: AbstractGraph

@present TheoryMaybeNamedGraph <: TheoryGraph begin
  VName::AttrType
  EName::AttrType
  vname::Attr(V, VName)
  ename::Attr(E, EName)
end

@acset_type _MaybeNamedGraph(TheoryMaybeNamedGraph, index=[:src,:tgt],
                             unique_index=[:vname]) <: AbstractGraph

""" Default graph type for [`@graph`](@ref) macro.

Vertex names are uniquely indexed and edge names are optional and unindexed.
"""
const MaybeNamedGraph{Name} = _MaybeNamedGraph{Name,Union{Nothing,Name}}

vertex_named(g, name) = incident(g, name, :vname)
edge_named(g, name)= incident(g, name, :ename)

# Graphs
########

""" Construct a graph in a simple, declarative style.

The syntax is reminiscent of Graphviz. Each line a declares a vertex or set of
vertices, or an edge. For example, the following defines a directed triangle:

```julia
@graph begin
  v0, v1, v2
  fst: v0 → v1
  snd: v1 → v2
  comp: v0 → v2
end
```

Vertices in the graph must be uniquely named, whereas edges names are optional
and need not be unique.
"""
macro graph(graph_type, body)
  :(parse_graph_dsl($(esc(graph_type)), $(Meta.quot(body))))
end
macro graph(body)
  :(parse_graph_dsl(MaybeNamedGraph{Symbol}, $(Meta.quot(body))))
end

function parse_graph_dsl(G::Type, body::Expr)
  g = G()
  body = Base.remove_linenums!(body)
  for expr in body.args
    @match reparse_arrow(expr) begin
      # v
      v::Symbol => add_vertex!(g, vname=v)
      # u, v, ...
      Expr(:tuple, vs...) => add_vertices!(g, length(vs), vname=vs)
      # u → v
      Expr(:call, :(→), u::Symbol, v::Symbol) => parse_edge!(g, u, v)
      # e : u → v
      Expr(:call, :(:), e::Symbol, Expr(:call, :(→), u::Symbol, v::Symbol)) =>
        parse_edge!(g, u, v, ename=e)
      # (e, f, ...) : u → v
      Expr(:call, (:), Expr(:tuple, es...), Expr(:call, :(→), u::Symbol, v::Symbol)) =>
        for e in es; parse_edge!(g, u, v, ename=e) end
      _ => error("@graph macro cannot parse line: $expr")
    end
  end
  return g
end

function parse_edge!(g, s::Symbol, t::Symbol;
                     ename::Union{Symbol,Nothing}=nothing)
  e = add_edge!(g, vertex_named(g, s), vertex_named(g, t))
  if has_subpart(g, :ename)
    g[e,:ename] = ename
  end
  return e
end

function reparse_arrow(expr)
  # `f : x → y` is parsed by Julia as `(f : x) → y`, not `f : (x → y)`.
  @match expr begin
    Expr(:call, :(→), Expr(:call, :(:), f, x), y) =>
      Expr(:call, :(:), f, Expr(:call, :(→), x, y))
    Expr(head, args...) => Expr(head, (reparse_arrow(arg) for arg in args)...)
    _ => expr
  end
end

# Categories
############

""" Present a category by generators and relations.

The result is a finitely presented category (`FinCat`) represented by a graph,
possibly with path equations. For example, the simplex category truncated to one
dimension is:

```julia
@category begin
  V, E
  (δ₀, δ₁): V → E
  σ₀: E → V

  σ₀ ∘ δ₀ == id(V)
  σ₀ ∘ δ₁ == id(V)
end
```

The objects and morphisms must be uniquely named.
"""
macro category(body)
  parse_category_dsl(NamedGraph{Symbol}, body)
end

function parse_category_dsl(G::Type, body::Expr)
  g, eqs = G(), Pair[]
  body = Base.remove_linenums!(body)
  for expr in body.args
    @match reparse_arrow(expr) begin
      # v
      v::Symbol => add_vertex!(g, vname=v)
      # u, v, ...
      Expr(:tuple, vs...) => add_vertices!(g, length(vs), vname=vs)
      # e : u → v
      Expr(:call, :(:), e::Symbol, Expr(:call, :(→), u::Symbol, v::Symbol)) =>
        parse_edge!(g, u, v, ename=e)
      # (e, f, ...) : u → v
      Expr(:call, (:), Expr(:tuple, es...), Expr(:call, :(→), u::Symbol, v::Symbol)) =>
        for e in es; parse_edge!(g, u, v, ename=e) end
      # f == g
      Expr(:call, :(==), lhs, rhs) => push!(eqs, parse_equation(g, lhs, rhs))
      _ => error("@category macro cannot parse line: $expr")
    end
  end
  isempty(eqs) ? FinCat(g) : FinCat(g, eqs)
end

""" Given generating graph, parse equation between morphisms.
"""
parse_equation(g, lhs, rhs) = parse_hom(g, lhs) => parse_hom(g, rhs)

""" Given generating graph, parse morphism expression.
"""
function parse_hom(g, expr)
  function parse(expr)
    @match expr begin
      Expr(:call, :compose, args...) => mapreduce(parse, vcat, args)
      Expr(:call, :(⋅), f, g) => vcat(parse(f), parse(g))
      Expr(:call, :(∘), f, g) => vcat(parse(g), parse(f))
      Expr(:call, :id, x) => empty(Path, g, vertex_named(g, x))
      f::Symbol => Path(g, edge_named(g, f))
      _ => error("Invalid morphism expression $expr")
    end
  end
  parse(expr)
end

# Diagrams
##########

""" Present a diagram in a given category.

Recall that a *diagram* in a category ``C`` is a functor ``F: J → C`` from a
small category ``J`` into ``C``. Given the category ``C``, this macro presents a
diagram in ``C``, i.e., constructs a finitely presented indexing category ``J``
together with a functor ``F: J → C``. This method of simultaneous definition is
often more convenient than defining ``J`` and ``F`` separately.

For example, the following diagram specifies the paths of length two in a graph:

```julia
@diagram FinCat(TheoryGraph) begin
  v::V
  (e1, e2)::E
  t::tgt : e1 → v
  s::src : e2 → v
end
```
"""
macro diagram(category, body)
  :(parse_diagram_dsl($(esc(category)), $(Meta.quot(body))))
end

function parse_diagram_dsl(C::Cat, body::Expr)
  g, eqs = NamedGraph{Symbol}(), Pair[]
  F_ob, F_hom = [], []
  body = Base.remove_linenums!(body)
  for expr in body.args
    @match reparse_arrow(expr) begin
      # x => X
      Expr(:call, :(=>), x::Symbol, X::Symbol) ||
      # x::X
      Expr(:(::), x::Symbol, X::Symbol) => begin
        add_vertex!(g, vname=x)
        push!(F_ob, X)
      end
      # (x, y, ...) => X
      Expr(:call, :(=>), Expr(:tuple, xs...), X::Symbol) ||
      # (x, y, ...)::X
      Expr(:(::), Expr(:tuple, xs...), X::Symbol) => begin
        add_vertices!(g, length(xs), vname=xs)
        append!(F_ob, repeated(X, length(xs)))
      end
      # (e: x → y) => f
      Expr(:call, :(=>), Expr(:call, :(:), e::Symbol,
                              Expr(:call, :(→), x::Symbol, y::Symbol)), f::Symbol) ||
      # e::f : x → y
      Expr(:call, :(:), Expr(:(::), e::Symbol, f::Symbol),
           Expr(:call, :(→), x::Symbol, y::Symbol)) => begin
        parse_edge!(g, x, y, ename=e)
        push!(F_hom, f)
      end
      # f == g
      Expr(:call, :(==), lhs, rhs) => push!(eqs, parse_equation(g, lhs, rhs))
      _ => error("@diagram macro cannot parse line: $expr")
    end
  end
  J = isempty(eqs) ? FinCat(g) : FinCat(g, eqs)
  F = FinDomFunctor(F_ob, F_hom, J, C)
  is_functorial(F) || error("Constructed @diagram is not functorial: $F")
  return F
end

end


text:
@testset "Dual Type" begin
    h = LPHRepresentation(spzeros(Int, 2, 2), [1, 2], [3, 4], [4, 5], [6, 7])
    p = @inferred polyhedron(h)
    @test p isa SimplePolyhedron{2, Rational{BigInt}, LPHRepresentation{2, Rational{BigInt}, SparseMatrixCSC{Rational{BigInt},Int}}, Polyhedra.Hull{2, Rational{BigInt}, Vector{Rational{BigInt}}}}
    h = hrep(zeros(2, 2), zeros(2))
    p = @inferred polyhedron(h)
    @test p isa SimplePolyhedron{2, Float64, MixedMatHRep{2, Float64, Matrix{Float64}}, MixedMatVRep{2, Float64, Matrix{Float64}}}
    h = hrep(spzeros(2, 2), zeros(2))
    p = @inferred polyhedron(h)
    @test p isa SimplePolyhedron{2, Float64, MixedMatHRep{2, Float64, SparseMatrixCSC{Float64, Int}}, MixedMatVRep{2, Float64, Matrix{Float64}}}
    v = vrep(zeros(2, 3))
    p = @inferred polyhedron(v)
    @test p isa SimplePolyhedron{3, Float64, MixedMatHRep{3, Float64, Matrix{Float64}}, MixedMatVRep{3, Float64, Matrix{Float64}}}
    v = vrep(spzeros(2, 3))
    p = @inferred polyhedron(v)
    @test p isa SimplePolyhedron{3, Float64, MixedMatHRep{3, Float64, Matrix{Float64}}, MixedMatVRep{3, Float64, SparseMatrixCSC{Float64, Int}}}
end

struct BadPoly{N, T} <: Polyhedra.Polyhedron{N, T}
end

@testset "Unimplemented methods" begin
    p = BadPoly{2, Int}()
    @test_throws ErrorException detectvlinearity!(p)
    @test_throws ErrorException detecthlinearity!(p)
    @test_throws ErrorException eliminate(p, [1], FourierMotzkin())
    @test_throws ErrorException eliminate(p, [1], BlockElimination())
    @test_throws ErrorException hrep(p)
    @test_throws ErrorException vrep(p)
    @test_throws ErrorException loadpolyhedron!(p, "a", "b")
    @test_throws ErrorException loadpolyhedron!(p, "a", "ine")
    @test_throws ErrorException loadpolyhedron!(p, "a", "ext")
end

@testset "SimplePolyhedron constructor with nothing" begin
    vr = convexhull([-1, 0], [0, -1]) + conichull([1, 1], [-1, 1])
    hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1)
    p = SimplePolyhedron{2, Int, typeof(hr), typeof(vr)}(hr, nothing, lp_solver)
    @test hrep(p) === hr
    @test !vrepiscomputed(p)
    @test p.solver === lp_solver
    p = SimplePolyhedron{2, Int, typeof(hr), typeof(vr)}(nothing, vr, lp_solver)
    @test !hrepiscomputed(p)
    @test vrep(p) === vr
    @test p.solver === lp_solver
end


text:
using Plots
using BenchmarkTools
#2(a)
function vfsolvex(vnew, kgrid, tolerance, imax, σ=1.5)
    β=0.9

    v=vnew .+ 2*tolerance
    cartesianindex = Array{CartesianIndex{2}, length(v)}
    i=1

    c = kgrid .- kgrid'
    c[c .< 0] .= 0
    if σ==1
      u = log.(c)
    else
      u = (c.^(1-σ).-1)./(1-σ)
    end
    u[kgrid .- kgrid' .< 0] .= -Inf

    while maximum(abs.(v-vnew)) > tolerance && i <= imax
      v = vnew;

      (vnew, cartesianindex) = findmax(u .+ β*v', dims = 2);
      i += 1;
    end
    kprimeindex = getindex.(cartesianindex, 2)
    return (v=vnew, kprime=kgrid[kprimeindex], kprimeindex=kprimeindex)
  end

  kupper = 5
  klower = 0.01
  n = 500
  kgrid = collect(range(klower, stop=kupper, length = n))
  (v, kprime, kprimeindex) = vfsolvex(zeros(n), kgrid, 0.001, 1000)
  #(b)
  (v, kprime, kprimeindex) = vfsolvex(zeros(500), kgrid, 0.001, 1000);
   plot(kgrid, v, label = "v, σ=1.5")

  #(c)
  plot()
  for x in [1.25, 1.1, 1.001, 1.00001, 1]
    kgrid = collect(range(klower, stop=kupper, length=500))
    (v, kprime, kprimeindex) = vfsolvex(zeros(500), kgrid, 0.001, n, x);
    display(plot!(kgrid, v, label = "v, σ=$x"))
  end
#as σ->1, the value function is getting closer to log utility value function
#the vfsolvex autometically compute log utility solution when σ=1.
#We also solve the log utility problem analytically in PS 5, and the numerical and
#analytical solution coinside

#(d)
function policy(x)
    e=getindex(findall(kgrid .== x),1)
    return kprime[e]
end

function findindex(x)
    a=findmin(abs.(x .- kgrid), dims=1)
    return a[2]
end

function findkpath(T, k_0)
    w = findindex(k_0)
    kpath = zeros(T)
    k = kgrid[w]
        for i in 1:T
            k = policy(k)
            kpath[i] = k
        end
    return kpath
end

function plotkpath(T, k_0)
    w = findindex(k_0)
    kpath = zeros(T)
    k = kgrid[w]
        for i in 1:T
            k = policy(k)
            kpath[i] = k
        end
    time = collect(range(1, T, length = T))
    plot(time, kpath, label = "kpath")
end

findkpath(100, 5)
plotkpath(100, 5)
# this is a cake eating problem with no production, so capital always decline

max_stars_repo_name: paradoxuu/Dynamic-Programming-Problems

max_issues_repo_name: paradoxuu/Dynamic-Programming-Problems

max_forks_repo_name: paradoxuu/Dynamic-Programming-Problems


text:
## A template for defining a new derivative evaluation method

# Define the method struct inheriting from either `GenerativeDerivativeMethod` or `NonGenerativeDerivativeMethod`
# Here we'll use `GenerativeDerivativeMethod` assuming our method will add extra supports
# This should contain any information that is needed to implement the method
struct MyDerivMethod <: GenerativeDerivativeMethod
    my_attr::Bool # REPLACE WITH ACTUAL ATTRIBUTE
    # ADD ANY MORE INFORMATION THAT IS NEEDED
end

# Extend `support_label` (only needed for generative methods)
function InfiniteOpt.support_label(method::MyDerivMethod)::DataType
    return InternalLabel # REPLACE WITH DESIRED SUPPORT LABEL FOR ANY EXTRA SUPPORTS THAT ARE ADDED
end

# Extend `generate_derivative_supports` (only needed for generative methods)
function InfiniteOpt.generate_derivative_supports(
    pref::IndependentParameterRef,
    method::MyDerivMethod
    )::Vector{Float64}
    # collect the supports already added to `pref`
    curr_supps = supports(pref, label = All)
    # generate the extra supports and return them (REPLACE BELOW WITH ACTUAL)
    min_s = minimum(curr_supps)
    max_s = maximum(curr_supps)
    return rand(2) * (max_s - min_s) .+ min_s
end

# Extend `evaluate_derivative`
# Generative methods must include calls of `add_derivative_supports`
# It will likely also be convenient to use `make_reduced_expr`
function InfiniteOpt.evaluate_derivative(
    dref::GeneralVariableRef,
    method::MyDerivMethod,
    write_model::JuMP.AbstractModel
    )::Vector{JuMP.AbstractJuMPScalar}
    # get the basic derivative information
    vref = derivative_argument(dref)
    pref = operator_parameter(dref)
    # make sure internal supports are added to the model
    InfiniteOpt.add_derivative_supports(pref)
    # generate the derivative expressions h_i corresponding to the equations of
    # the form h_i = 0 (REPLACE BELOW WITH ACTUAL IMPLEMENTATION)
    supps = supports(pref, label = All)
    exprs = Vector{JuMP.AbstractJuMPScalar}(undef, length(supps) - 1)
    for i in eachindex(exprs)
        d = InfiniteOpt.make_reduced_expr(dref, pref, supps[i], write_model)
        v1 = InfiniteOpt.make_reduced_expr(vref, pref, supps[i], write_model)
        v2 = InfiniteOpt.make_reduced_expr(vref, pref, supps[i + 1], write_model)
        change = supps[i + 1] - supps[i]
        exprs[i] = JuMP.@expression(write_model, -change * d + v2 - v1)
    end
    return exprs
end


text:
abstract type LayerContent end
abstract type GeometricObject <: LayerContent end

export BBox, bbox

struct Angle
    rad::Float64
end

struct Transform
    scale::Float64
    rotation::Angle
    translation::SVector{2, Float64}
end

struct Point <: GeometricObject
    xy::SVector{2, Float64}

    Point(x::Real, y::Real) = new(SVector(convert(Float64, x), convert(Float64, y)))
    function Point(xy)
        if Base.length(xy) != 2
            error("need 2 elements")
        end
        new(SVector{2, Float64}(xy...))
    end
end

const P = Point

struct Line <: GeometricObject
    from::Point
    to::Point
end

struct Arc <: GeometricObject
    center::Point
    radius::Float64
    start_angle::Angle
    end_angle::Angle
end

struct Bezier <: GeometricObject
    from::Point
    c1::Point
    c2::Point
    to::Point
end

struct Rect <: GeometricObject
    center::Point
    width::Float64
    height::Float64
    angle::Angle
end

struct Polygon <: GeometricObject
    points::Vector{Point}
end

struct Circle <: GeometricObject
    center::Point
    radius::Float64
end

struct TextExtent
    xbearing::Float64
    ybearing::Float64
    width::Float64
    height::Float64
    xadvance::Float64
    yadvance::Float64
end

struct Txt <: GeometricObject
    pos::Point
    text::String
    size::Float64
    halign::Symbol
    valign::Symbol
    angle::Angle
    font::String
    extent::Union{Nothing, TextExtent}
end

abstract type PathCommand end

export Move, Lineto, Curveto, Close

struct Move <: PathCommand
    p::Point
end

struct Lineto <: PathCommand
    p::Point
end

struct CurveTo <: PathCommand
    c1::Point
    c2::Point
    p::Point
end

struct Close <: PathCommand
end

struct Path <: GeometricObject
    commands::Vector{PathCommand}
end

geoms = (:Point, :Line, :Arc, :Circle, :Rect, :Bezier, :Path, :Polygon, :Txt)

abstract type Attribute end

mutable struct Attributes
    attrs::Dict{Type{<:Attribute}, Union{Function, Attribute}}
end

struct Clip
    shape
    # shape::Union{Nothing, Shape, Tuple{Shape, Shape}}
end

Base.Broadcast.broadcastable(c::Clip) = Ref(c)

struct Opacity
    opacity::Float64
end

Base.Broadcast.broadcastable(o::Opacity) = Ref(o)

struct Operator
    operator::Symbol
end

Base.Broadcast.broadcastable(o::Operator) = Ref(o)

mutable struct Layer <: LayerContent
    transform::Union{Tuple{Function, Vector}, Transform}
    content::Vector{LayerContent}
    parent::Union{Layer, Nothing}
    attrs::Attributes
    clip::Clip
    opacity::Opacity
    operator::Operator
end

mutable struct Shape{T <: GeometricObject} <: LayerContent
    # make this nicer and more julian
    content::Union{Tuple{Function, Vector}, T}
    parent::Union{Layer, Nothing}
    solved::Union{T, Nothing}
    attrs::Attributes
    clip::Clip
    opacity::Opacity
    operator::Operator
end

mutable struct Shapes{T <: GeometricObject} <: LayerContent
    # make this nicer and more julian
    content::Union{Tuple{Function, Vector}, Array{T}}
    parent::Union{Layer, Nothing}
    solved::Union{Array{T}, Nothing}
    attrs::Attributes
    clip::Clip
    opacity::Opacity
    operator::Operator
end

struct Gradient
    from::Point
    to::Point
    stops::Vector{Float64}
    colors::Vector{<:Colors.Colorant}
end

struct RadialGradient
    from::Circle
    to::Circle
    stops::Vector{Float64}
    colors::Vector{<:Colors.Colorant}
end

Clip(s1::Shape, s2::Shape) = Clip((s1, s2))

const FillContent = Union{Nothing, Colors.Colorant, Gradient, RadialGradient, Array{<:Colors.Colorant}, Pair{Symbol, <:Function}}
struct Fill{T<:FillContent} <: Attribute
    val::T
end

const TextfillContent = FillContent
struct Textfill{T<:TextfillContent} <: Attribute
    val::T
end

const StrokeContent = Union{Nothing, Colors.Colorant, Gradient, RadialGradient, Array{<:Colors.Colorant}, Pair{Symbol, <:Function}}
struct Stroke{T<:StrokeContent} <: Attribute
    val::T
end

const LinestyleContent = Symbol
struct Linestyle{T<:LinestyleContent} <: Attribute
    val::T
end

const LinewidthContent = Real
struct Linewidth{T<:LinewidthContent} <: Attribute
    val::T
end

const MarkersizeContent = Real
struct Markersize{T<:MarkersizeContent} <: Attribute
    val::T
end

const MarkerContent = Symbol
struct Marker{T<:MarkerContent} <: Attribute
    val::T
end

const VisibleContent = Bool
struct Visible{T<:VisibleContent} <: Attribute
    val::T
end

struct BBox
    from::Point
    to::Point
end


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

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

Cancels a contact with a specified contact ID.

# Arguments
- `contact_id`: UUID of a contact.

"""
cancel_contact(contactId; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("DELETE", "/contact/$(contactId)"; aws_config=aws_config)
cancel_contact(contactId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("DELETE", "/contact/$(contactId)", params; aws_config=aws_config)

"""
    create_config(config_data, name)
    create_config(config_data, name, params::Dict{String,<:Any})

Creates a Config with the specified configData parameters. Only one type of configData can
be specified.

# Arguments
- `config_data`: Parameters of a Config.
- `name`: Name of a Config.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"tags"`: Tags assigned to a Config.
"""
create_config(configData, name; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/config", Dict{String, Any}("configData"=>configData, "name"=>name); aws_config=aws_config)
create_config(configData, name, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/config", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("configData"=>configData, "name"=>name), params)); aws_config=aws_config)

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

Creates a DataflowEndpoint group containing the specified list of DataflowEndpoint objects.
The name field in each endpoint is used in your mission profile DataflowEndpointConfig to
specify which endpoints to use during a contact. When a contact uses multiple
DataflowEndpointConfig objects, each Config must match a DataflowEndpoint in the same group.

# Arguments
- `endpoint_details`: Endpoint details of each endpoint in the dataflow endpoint group.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"tags"`: Tags of a dataflow endpoint group.
"""
create_dataflow_endpoint_group(endpointDetails; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/dataflowEndpointGroup", Dict{String, Any}("endpointDetails"=>endpointDetails); aws_config=aws_config)
create_dataflow_endpoint_group(endpointDetails, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/dataflowEndpointGroup", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("endpointDetails"=>endpointDetails), params)); aws_config=aws_config)

"""
    create_mission_profile(dataflow_edges, minimum_viable_contact_duration_seconds, name, tracking_config_arn)
    create_mission_profile(dataflow_edges, minimum_viable_contact_duration_seconds, name, tracking_config_arn, params::Dict{String,<:Any})

Creates a mission profile.  dataflowEdges is a list of lists of strings. Each lower level
list of strings has two elements: a from ARN and a to ARN.

# Arguments
- `dataflow_edges`: A list of lists of ARNs. Each list of ARNs is an edge, with a from
  Config and a to Config.
- `minimum_viable_contact_duration_seconds`: Smallest amount of time in seconds that
  you’d like to see for an available contact. AWS Ground Station will not present you with
  contacts shorter than this duration.
- `name`: Name of a mission profile.
- `tracking_config_arn`: ARN of a tracking Config.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"contactPostPassDurationSeconds"`: Amount of time after a contact ends that you’d like
  to receive a CloudWatch event indicating the pass has finished.
- `"contactPrePassDurationSeconds"`: Amount of time prior to contact start you’d like to
  receive a CloudWatch event indicating an upcoming pass.
- `"tags"`: Tags assigned to a mission profile.
"""
create_mission_profile(dataflowEdges, minimumViableContactDurationSeconds, name, trackingConfigArn; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/missionprofile", Dict{String, Any}("dataflowEdges"=>dataflowEdges, "minimumViableContactDurationSeconds"=>minimumViableContactDurationSeconds, "name"=>name, "trackingConfigArn"=>trackingConfigArn); aws_config=aws_config)
create_mission_profile(dataflowEdges, minimumViableContactDurationSeconds, name, trackingConfigArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/missionprofile", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("dataflowEdges"=>dataflowEdges, "minimumViableContactDurationSeconds"=>minimumViableContactDurationSeconds, "name"=>name, "trackingConfigArn"=>trackingConfigArn), params)); aws_config=aws_config)

"""
    delete_config(config_id, config_type)
    delete_config(config_id, config_type, params::Dict{String,<:Any})

Deletes a Config.

# Arguments
- `config_id`: UUID of a Config.
- `config_type`: Type of a Config.

"""
delete_config(configId, configType; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("DELETE", "/config/$(configType)/$(configId)"; aws_config=aws_config)
delete_config(configId, configType, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("DELETE", "/config/$(configType)/$(configId)", params; aws_config=aws_config)

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

Deletes a dataflow endpoint group.

# Arguments
- `dataflow_endpoint_group_id`: UUID of a dataflow endpoint group.

"""
delete_dataflow_endpoint_group(dataflowEndpointGroupId; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("DELETE", "/dataflowEndpointGroup/$(dataflowEndpointGroupId)"; aws_config=aws_config)
delete_dataflow_endpoint_group(dataflowEndpointGroupId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("DELETE", "/dataflowEndpointGroup/$(dataflowEndpointGroupId)", params; aws_config=aws_config)

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

Deletes a mission profile.

# Arguments
- `mission_profile_id`: UUID of a mission profile.

"""
delete_mission_profile(missionProfileId; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("DELETE", "/missionprofile/$(missionProfileId)"; aws_config=aws_config)
delete_mission_profile(missionProfileId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("DELETE", "/missionprofile/$(missionProfileId)", params; aws_config=aws_config)

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

Describes an existing contact.

# Arguments
- `contact_id`: UUID of a contact.

"""
describe_contact(contactId; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/contact/$(contactId)"; aws_config=aws_config)
describe_contact(contactId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/contact/$(contactId)", params; aws_config=aws_config)

"""
    get_config(config_id, config_type)
    get_config(config_id, config_type, params::Dict{String,<:Any})

Returns Config information. Only one Config response can be returned.

# Arguments
- `config_id`: UUID of a Config.
- `config_type`: Type of a Config.

"""
get_config(configId, configType; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/config/$(configType)/$(configId)"; aws_config=aws_config)
get_config(configId, configType, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/config/$(configType)/$(configId)", params; aws_config=aws_config)

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

Returns the dataflow endpoint group.

# Arguments
- `dataflow_endpoint_group_id`: UUID of a dataflow endpoint group.

"""
get_dataflow_endpoint_group(dataflowEndpointGroupId; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/dataflowEndpointGroup/$(dataflowEndpointGroupId)"; aws_config=aws_config)
get_dataflow_endpoint_group(dataflowEndpointGroupId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/dataflowEndpointGroup/$(dataflowEndpointGroupId)", params; aws_config=aws_config)

"""
    get_minute_usage(month, year)
    get_minute_usage(month, year, params::Dict{String,<:Any})

Returns the number of minutes used by account.

# Arguments
- `month`: The month being requested, with a value of 1-12.
- `year`: The year being requested, in the format of YYYY.

"""
get_minute_usage(month, year; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/minute-usage", Dict{String, Any}("month"=>month, "year"=>year); aws_config=aws_config)
get_minute_usage(month, year, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/minute-usage", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("month"=>month, "year"=>year), params)); aws_config=aws_config)

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

Returns a mission profile.

# Arguments
- `mission_profile_id`: UUID of a mission profile.

"""
get_mission_profile(missionProfileId; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/missionprofile/$(missionProfileId)"; aws_config=aws_config)
get_mission_profile(missionProfileId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/missionprofile/$(missionProfileId)", params; aws_config=aws_config)

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

Returns a satellite.

# Arguments
- `satellite_id`: UUID of a satellite.

"""
get_satellite(satelliteId; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/satellite/$(satelliteId)"; aws_config=aws_config)
get_satellite(satelliteId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/satellite/$(satelliteId)", params; aws_config=aws_config)

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

Returns a list of Config objects.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Maximum number of Configs returned.
- `"nextToken"`: Next token returned in the request of a previous ListConfigs call. Used to
  get the next page of results.
"""
list_configs(; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/config"; aws_config=aws_config)
list_configs(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/config", params; aws_config=aws_config)

"""
    list_contacts(end_time, start_time, status_list)
    list_contacts(end_time, start_time, status_list, params::Dict{String,<:Any})

Returns a list of contacts. If statusList contains AVAILABLE, the request must include
groundStation, missionprofileArn, and satelliteArn.

# Arguments
- `end_time`: End time of a contact.
- `start_time`: Start time of a contact.
- `status_list`: Status of a contact reservation.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"groundStation"`: Name of a ground station.
- `"maxResults"`: Maximum number of contacts returned.
- `"missionProfileArn"`: ARN of a mission profile.
- `"nextToken"`: Next token returned in the request of a previous ListContacts call. Used
  to get the next page of results.
- `"satelliteArn"`: ARN of a satellite.
"""
list_contacts(endTime, startTime, statusList; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/contacts", Dict{String, Any}("endTime"=>endTime, "startTime"=>startTime, "statusList"=>statusList); aws_config=aws_config)
list_contacts(endTime, startTime, statusList, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/contacts", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("endTime"=>endTime, "startTime"=>startTime, "statusList"=>statusList), params)); aws_config=aws_config)

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

Returns a list of DataflowEndpoint groups.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Maximum number of dataflow endpoint groups returned.
- `"nextToken"`: Next token returned in the request of a previous
  ListDataflowEndpointGroups call. Used to get the next page of results.
"""
list_dataflow_endpoint_groups(; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/dataflowEndpointGroup"; aws_config=aws_config)
list_dataflow_endpoint_groups(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/dataflowEndpointGroup", params; aws_config=aws_config)

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

Returns a list of ground stations.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Maximum number of ground stations returned.
- `"nextToken"`: Next token that can be supplied in the next call to get the next page of
  ground stations.
- `"satelliteId"`: Satellite ID to retrieve on-boarded ground stations.
"""
list_ground_stations(; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/groundstation"; aws_config=aws_config)
list_ground_stations(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/groundstation", params; aws_config=aws_config)

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

Returns a list of mission profiles.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Maximum number of mission profiles returned.
- `"nextToken"`: Next token returned in the request of a previous ListMissionProfiles call.
  Used to get the next page of results.
"""
list_mission_profiles(; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/missionprofile"; aws_config=aws_config)
list_mission_profiles(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/missionprofile", params; aws_config=aws_config)

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

Returns a list of satellites.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"maxResults"`: Maximum number of satellites returned.
- `"nextToken"`: Next token that can be supplied in the next call to get the next page of
  satellites.
"""
list_satellites(; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/satellite"; aws_config=aws_config)
list_satellites(params::AbstractDict{String, Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/satellite", params; aws_config=aws_config)

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

Returns a list of tags for a specified resource.

# Arguments
- `resource_arn`: ARN of a resource.

"""
list_tags_for_resource(resourceArn; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/tags/$(resourceArn)"; aws_config=aws_config)
list_tags_for_resource(resourceArn, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("GET", "/tags/$(resourceArn)", params; aws_config=aws_config)

"""
    reserve_contact(end_time, ground_station, mission_profile_arn, satellite_arn, start_time)
    reserve_contact(end_time, ground_station, mission_profile_arn, satellite_arn, start_time, params::Dict{String,<:Any})

Reserves a contact using specified parameters.

# Arguments
- `end_time`: End time of a contact.
- `ground_station`: Name of a ground station.
- `mission_profile_arn`: ARN of a mission profile.
- `satellite_arn`: ARN of a satellite
- `start_time`: Start time of a contact.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"tags"`: Tags assigned to a contact.
"""
reserve_contact(endTime, groundStation, missionProfileArn, satelliteArn, startTime; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/contact", Dict{String, Any}("endTime"=>endTime, "groundStation"=>groundStation, "missionProfileArn"=>missionProfileArn, "satelliteArn"=>satelliteArn, "startTime"=>startTime); aws_config=aws_config)
reserve_contact(endTime, groundStation, missionProfileArn, satelliteArn, startTime, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/contact", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("endTime"=>endTime, "groundStation"=>groundStation, "missionProfileArn"=>missionProfileArn, "satelliteArn"=>satelliteArn, "startTime"=>startTime), params)); aws_config=aws_config)

"""
    tag_resource(resource_arn, tags)
    tag_resource(resource_arn, tags, params::Dict{String,<:Any})

Assigns a tag to a resource.

# Arguments
- `resource_arn`: ARN of a resource tag.
- `tags`: Tags assigned to a resource.

"""
tag_resource(resourceArn, tags; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/tags/$(resourceArn)", Dict{String, Any}("tags"=>tags); aws_config=aws_config)
tag_resource(resourceArn, tags, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("POST", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tags"=>tags), params)); aws_config=aws_config)

"""
    untag_resource(resource_arn, tag_keys)
    untag_resource(resource_arn, tag_keys, params::Dict{String,<:Any})

Deassigns a resource tag.

# Arguments
- `resource_arn`: ARN of a resource.
- `tag_keys`: Keys of a resource tag.

"""
untag_resource(resourceArn, tagKeys; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("DELETE", "/tags/$(resourceArn)", Dict{String, Any}("tagKeys"=>tagKeys); aws_config=aws_config)
untag_resource(resourceArn, tagKeys, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("DELETE", "/tags/$(resourceArn)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("tagKeys"=>tagKeys), params)); aws_config=aws_config)

"""
    update_config(config_data, config_id, config_type, name)
    update_config(config_data, config_id, config_type, name, params::Dict{String,<:Any})

Updates the Config used when scheduling contacts. Updating a Config will not update the
execution parameters for existing future contacts scheduled with this Config.

# Arguments
- `config_data`: Parameters of a Config.
- `config_id`: UUID of a Config.
- `config_type`: Type of a Config.
- `name`: Name of a Config.

"""
update_config(configData, configId, configType, name; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("PUT", "/config/$(configType)/$(configId)", Dict{String, Any}("configData"=>configData, "name"=>name); aws_config=aws_config)
update_config(configData, configId, configType, name, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("PUT", "/config/$(configType)/$(configId)", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("configData"=>configData, "name"=>name), params)); aws_config=aws_config)

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

Updates a mission profile. Updating a mission profile will not update the execution
parameters for existing future contacts.

# Arguments
- `mission_profile_id`: UUID of a mission profile.

# Optional Parameters
Optional parameters can be passed as a `params::Dict{String,<:Any}`. Valid keys are:
- `"contactPostPassDurationSeconds"`: Amount of time after a contact ends that you’d like
  to receive a CloudWatch event indicating the pass has finished.
- `"contactPrePassDurationSeconds"`: Amount of time after a contact ends that you’d like
  to receive a CloudWatch event indicating the pass has finished.
- `"dataflowEdges"`: A list of lists of ARNs. Each list of ARNs is an edge, with a from
  Config and a to Config.
- `"minimumViableContactDurationSeconds"`: Smallest amount of time in seconds that you’d
  like to see for an available contact. AWS Ground Station will not present you with contacts
  shorter than this duration.
- `"name"`: Name of a mission profile.
- `"trackingConfigArn"`: ARN of a tracking Config.
"""
update_mission_profile(missionProfileId; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("PUT", "/missionprofile/$(missionProfileId)"; aws_config=aws_config)
update_mission_profile(missionProfileId, params::AbstractDict{String, <:Any}; aws_config::AbstractAWSConfig=global_aws_config()) = groundstation("PUT", "/missionprofile/$(missionProfileId)", params; aws_config=aws_config)


text:
# Simple tests for TaylorSeries implementation
using TaylorSeries
using Base.Test

# Tests for 1-d Tylor expansions
xT(a) = Taylor([a,one(a)],15)
xT0 = Taylor(xT(0),15)
xTI = im*xT0
z = zero(xT0)
u = 1.0*one(xT0)
tol1 = eps(1.0)

@test eltype(convert(Taylor{Complex128},u)) == Complex128
@test eltype(convert(Taylor{Complex128},1)) == Complex128
@test promote(0,Taylor(1.0,0)) == (z,u)
@test eltype(promote(0,Taylor(u))[1]) == Float64
@test eltype(promote(1.0+im, z)[1]) == Complex{Float64}
@test eltype(TaylorSeries.fixshape(z,u)[1]) == Float64

@test length(Taylor(0)) == 0
@test length(TaylorSeries.fixshape(z,convert(Taylor{Int64},[0]))[1]) == 15
@test TaylorSeries.firstnonzero(xT0) == 1
@test TaylorSeries.firstnonzero(z) == z.order+1

@test u == 1
@test 0.0 == z
@test xT0.coeffs[2] == 1
@test z+1 == u
@test xT0+xT0 == 2xT0
@test xT0-xT0 == z

xsquare = Taylor([0,0,1],15)
@test xT0^0 == xT0^0.0 == one(xT0)
@test xT0*xT0 == xsquare
@test (-xT0)^2 == xsquare
@test xsquare/xT0 == xT0
@test xT0/(xT0*3) == (1/3)*u
@test xT0/3im == -xTI/3
@test 1/(1-xT0) == Taylor(ones(xT0.order+1))
@test Taylor([0,1,1])/xT0 == xT0+1
@test (xT0+im)^2 == xsquare+2im*xT0-1
@test (xT0+im)^3 == Taylor([-1im,-3,3im,1],15)
@test (xT0+im)^4 == Taylor([1,-4im,-6,4im,1],15)
@test imag(xsquare+2im*xT0-1) == 2xT0
@test (Rational(1,2)*xsquare).coeffs[3] == 1//2
@test xT0^2/xsquare == u
@test ((1+xT0)^(1/3)).coeffs[3]+1/9 <= tol1
@test 1-xsquare == (1+xT0)-xT0*(1+xT0)
@test (1-xsquare)^2 == (1+xT0)^2 * (1-xT0)^2
@test (sqrt(1+xT0)).coeffs[3] == -1/8
@test ((1-xsquare)^(1//2))^2 == 1-xsquare
@test ((1-xT0)^(1//4)).coeffs[15] == -4188908511//549755813888
@test abs(((1+xT0)^3.2).coeffs[14] + 5.4021062656e-5) < tol1

@test isapprox( rem(4.1 + xT0,4).coeffs[1], (0.1 + xT0).coeffs[1] )
@test isapprox( mod(4.1 + xT0,4).coeffs[1], (0.1 + xT0).coeffs[1] )
@test isapprox( mod2pi(2pi + 0.1 + xT0).coeffs[1], (0.1 + xT0).coeffs[1] )

@test log(exp(xsquare)) == xsquare
@test exp(log(1-xsquare)) == 1-xsquare
@test log((1-xT0)^2) == 2*log(1-xT0)
@test real(exp(xTI)) == cos(xT0)
@test imag(exp(xTI)) == sin(xT0)
@test exp(xTI') == cos(xT0)-im*sin(xT0)
@test (exp(xT0))^(2im) == cos(2xT0)+im*sin(2xT0)
@test (exp(xT0))^Taylor(-5.2im) == cos(5.2xT0)-im*sin(5.2xT0)
@test abs((tan(xT0)).coeffs[8]- 17/315) < tol1
@test abs((tan(xT0)).coeffs[14]- 21844/6081075) < tol1
@test evalTaylor(exp(Taylor([0,1],17)),1.0) == e
@test evalTaylor(exp(Taylor([0,1],1))) == 1.0
@test evalTaylor(exp(xT0),xT0^2) == exp(xT0^2)

@test deriv( exp(xT(1.0)), 5 ) == exp(1.0)
@test deriv( exp(xT(1.0pi)), 3 ) == exp(1.0pi)
@test isapprox( deriv(exp(xT(1.0pi)), 10) , exp(1.0pi) )
@test integTaylor(diffTaylor(exp(xT0)),1) == exp(xT0)
@test integTaylor(cos(xT0)) == sin(xT0)

@test_throws ErrorException 1/xT0
@test_throws ErrorException z/z
@test_throws ErrorException xT0^1.5
@test_throws ErrorException sqrt(xT0)
@test_throws ErrorException log(xT0)
@test_throws ErrorException cos(xT0)/sin(xT0)
@test_throws ErrorException deriv( exp(xT(1.0pi)), 30 )

# Tests for 1-d Tylor expansions
set_numVars(2)
set_maxOrder(17)
xTN = TaylorN([0,1,0],17)
yTN = TaylorN([0,0,1])
zN = zero(xTN)
uN = convert(TaylorN{Float64},one(yTN))

@test TaylorN(zN,5) == 0
@test TaylorN(uN) == convert(TaylorN{Complex},1)
@test get_numVars(zN) == zN.numVars
@test length(TaylorN(0,0)) == binomial( get_numVars() + get_maxOrder(), get_maxOrder(zN) )
@test eltype(convert(TaylorN{Complex128},[1])) == Complex128

@test TaylorN(1)+xTN+yTN == TaylorN([1,1,1])
@test xTN-yTN == TaylorN([0,1,-1])
@test xTN*yTN == TaylorN([0,0,0,0,1,0])
@test (1/(1-xTN)).coeffs[4] == 1.0
@test (yTN/(1-xTN)).coeffs[5] == 1.0
@test mod(1+xTN,1) == +xTN
@test (rem(1+xTN,1)).coeffs[1] == 0
@test diffTaylor(mod2pi(2pi+yTN),2) == diffTaylor(yTN,2)
@test diffTaylor(yTN) == zN

@test diffTaylor(2xTN*yTN^2,1) == 2yTN^2
@test xTN*xTN^3 == xTN^4
@test (1+xTN)^(3//2) == ((1+xTN)^0.5)^3
@test real( exp(1im * xTN)) == cos(xTN)
@test imag((exp(yTN))^(-1im)) == -sin(yTN)
@test evalTaylor(exp( xTN+yTN )) == 1
@test isapprox(evalTaylor(exp( xTN+yTN ), [1,1]), e^2)

println("    \033[32;1mSUCCESS\033[0m")


text:
"""
    Pendulum <: Environment

    classic system with one rotational degree of freedom
    https://underactuated.mit.edu/pend.html
"""
struct Pendulum end

function pendulum(;
    representation=:minimal,
    max_speed=8.0,
    max_torque=8.0,
    timestep=0.05,
    gravity=-10.0,
    mass=1.0,
    len=1.0,
    damper=0.0,
    seed=1,
    vis=Visualizer(),
    name=:robot,
    opts_step=SolverOptions(),
    opts_grad=SolverOptions(),
    T=Float64)

    mechanism = get_mechanism(:pendulum,
        timestep=timestep,
        gravity=gravity,
        mass=mass,
        len=len,
        damper=damper)

    initialize!(mechanism, :pendulum)

    if representation == :minimal
        nx = minimal_dimension(mechanism)
        no = 3
    elseif representation == :maximal
        nx = maximal_dimension(mechanism)
        no = 13
    end
    nu = input_dimension(mechanism)

    high = [1.0, 1.0, max_speed]
    aspace = BoxSpace(input_dimension(mechanism),
        low=[-timestep * max_torque],
        high=[timestep * max_torque])
    ospace = BoxSpace(no,
        low=-high,
        high=high)
    rng = [MersenneTwister(seed),]

    x = Inf * ones(nx)
    fx = zeros(nx, nx)
    fu = zeros(nx, nu)

    u_prev = Inf * ones(nu)
    control_mask = ones(1,1)
    control_scaling = Diagonal(ones(nu))
    build_robot(mechanism,
        vis=vis,
        name=name)

    info = Dict(:maximal_speed => max_speed, :maximal_torque => max_torque)

    TYPES = [T, typeof(mechanism), typeof(aspace), typeof(ospace), typeof(info)]
    env = Environment{Pendulum, TYPES...}(mechanism, representation, aspace, ospace,
        x, fx, fu,
        u_prev, control_mask' * control_scaling, nx, nu, no, info,
        rng, vis, opts_step, opts_grad)
    return env
end

function Base.reset(env::Environment{Pendulum};
    x=nothing)

    initialize!(env.mechanism, :pendulum)

    if x != nothing
        env.state .= x
    else
        if env.representation == :minimal
            high = [π, 1.0]
            low = -high
            env.state .= rand(env.rng[1], env.num_states) .* (high .- low) .+ low
        elseif env.representation == :maximal
            env.state .= pendulum_nominal_max()
        end
        env.input_previous .= Inf
    end
    return get_observation(env)
end

function get_observation(env::Environment{Pendulum})
    if env.representation == :minimal
        θ, ω = env.state
        return [cos(θ), sin(θ), ω]
    else env.representation == :maximal
        return env.state
    end
end

function Base.step(env::Environment{Pendulum}, x, u;
    gradients=false,
    attitude_decompress=false)
    mechanism = env.mechanism
    timestep= mechanism.timestep
    max_torque = env.info[:maximal_torque]

    x0 = x
    u0 = clamp.(u, -max_torque, max_torque)
    env.input_previous .= u0

    z0 = env.representation == :minimal ? minimal_to_maximal(mechanism, x0) : x0
    z1 = step!(mechanism, z0, timestep * u0; opts = env.opts_step)
    env.state .= env.representation == :minimal ? maximal_to_minimal(mechanism, z1) : z1

    # Compute cost function
    costs = cost(env, x0, u0)

    # Gradients
    if gradients
        if env.representation == :minimal
            fx, fu = get_minimal_gradients!(env.mechanism, z0, timestep * u0,
                opts=env.opts_grad)
        elseif env.representation == :maximal
            fx, fu = get_maximal_gradients!(env.mechanism, z0, timestep * u0,
                opts=env.opts_grad)
            if attitude_decompress
                A0 = attitude_jacobian(z0, length(env.mechanism.bodies))
                A1 = attitude_jacobian(z1, length(env.mechanism.bodies))
                fx = A1 * fx * A0'
                fu = A1 * fu
            end
        end
        env.dynamics_jacobian_state .= fx
        env.dynamics_jacobian_input .= timestep * fu
    end

    info = Dict()
    return get_observation(env), -costs, false, info
end

function angle_normalize(x)
    return ((x + 101π) % (2 * π)) - π
end

function pendulum_nominal_max()
    x1 = [0.0; 0.0; -0.5]
    v1 = [0.0; 0.0; 0.0]
    q1 = [1.0; 0.0; 0.0; 0.0]
    ω1 = [0.0; 0.0; 0.0]
    z1 = [x1; v1; q1; ω1]
end

function pendulum_goal_max()
    xT = [0.0; 0.0; 0.5]
    vT = [0.0; 0.0; 0.0]
    qT = [0.0; 1.0; 0.0; 0.0]
    ωT = [0.0; 0.0; 0.0]
    zT = [xT; vT; qT; ωT]
end

function cost(env::Environment{Pendulum}, x, u)
    if env.representation == :minimal
        θ, ω = x
        c = angle_normalize(θ - π)^2 + 1e-1 * ω^2 + 1e-3 * (u[1])^2 # angle_normalize enforces angle ∈ [-π, π]
        c = angle_normalize(θ - π)^2 + 1e-3 * ω^2 + 1e-3 * (env.mechanism.timestep * u[1])^2 # angle_normalize enforces angle ∈ [-π, π]
    else
        c = Inf
    end
    return c * 0.1
end


text:
module CommunalHelperDreamCrumbleWallOnRumble

using ..Ahorn, Maple
using Ahorn.CommunalHelper
using Ahorn.CommunalHelperEntityPresets: CustomDreamBlockData

const entityData = appendkwargs(CustomDreamBlockData, :(
    permanent::Bool=false,
))
@mapdefdata Entity "CommunalHelper/DreamCrumbleWallOnRumble" DreamCrumbleWallOnRumble entityData

const placements = Ahorn.PlacementDict(
    "Dream Crumble Wall On Rumble (Communal Helper)" => Ahorn.EntityPlacement(
        DreamCrumbleWallOnRumble,
        "rectangle",
    ),
)

Ahorn.minimumSize(entity::DreamCrumbleWallOnRumble) = 8, 8
Ahorn.resizable(entity::DreamCrumbleWallOnRumble) = true, true

function Ahorn.selection(entity::DreamCrumbleWallOnRumble)
    x, y = Ahorn.position(entity)

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    return Ahorn.Rectangle(x, y, width, height)
end

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::DreamCrumbleWallOnRumble, room::Maple.Room)
    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    renderDreamBlock(ctx, 0, 0, width, height, entity.data)
end

end