text:
# Robust Optimization
# 1. Function for MILP with Box Uncertainty and Budget of Uncertainty
# 2. Fucntion for Two-Stage Stochastic LP with Box Uncertainty
# Version: 3.0
# Author: Edward J. Xu, edxu96@outlook.com
# Date: April 8th, 2019
module RobustOptim
using JuMP
using GLPKMathProgInterface
# using PrettyTables

# Function for MILP with Box Uncertainty and Budget of Uncertainty
function milpBoxBudget(; num_x, num_y, vec_min_y, vec_max_y, vec_c, vec_f, vec_b, mat_a, mat_b,
    mat_a_bar, mat_a_hat, mat_b_bar, vec_b_bar, vec_gammaCap)
    m = length(mat_a_bar[:,1])
    n = length(mat_a_bar[1,:])
    for i = 1: m
        if vec_gammaCap[i] > n
            println("Error: vec_gammaCap[$i] is greater than m in nominal and hat data!")
        end
    end
    println("---------------------------------------------------------\n",
            "    1/2. Robust MILP with Box Uncertainty and Budget\n",
            "---------------------------------------------------------\n")
    println("Input: vec_gammaCap = $vec_gammaCap.")
    #
    model = Model(solver = GLPKSolverMIP())
    @variable(model, vec_y[1: num_y], Int)
    @variable(model, vec_x[1: num_x] >= 0)
    @objective(model, Min, (transpose(vec_c) * vec_x + transpose(vec_f) * vec_y)[1])
    @constraint(model, vec_y[1: num_y] .<= vec_max_y)
    @constraint(model, vec_y[1: num_y] .>= vec_min_y)
    @constraint(model, mat_a * vec_x + mat_b * vec_y .>= vec_b)
    # Transformation of Box Uncertainty.
    @variable(model, vec_lambda[1: m] >= 0)
    @variable(model, mat_mu[1: m, 1: num_x] >= 0)
    @variable(model, vec_z[1: num_x] >= 0)
    @constraint(model, - vec_z .<= vec_x)
    @constraint(model, vec_x .<= vec_z)
    for i = 1: m
        # sum() of variables without coefficients can be used directly
        @constraint(model, transpose(mat_a_bar[i, :]) * vec_x + vec_gammaCap[i] * vec_lambda[i] +
            sum(mat_mu[i, :]) + transpose(mat_b_bar[i, :]) * vec_y <= vec_b_bar[i])
        @constraint(model, vec_lambda[i] .+ hcat(mat_mu[i, :]) .>=  hcat(mat_a_hat[i, :]) .* vec_z[:])
    end
    solve(model)
    obj_result = getobjectivevalue(model)
    vec_result_y = getvalue(vec_y)
    println("Result: vec_y = $vec_result_y.")
    vec_result_x = getvalue(vec_x)
    println("Result: vec_x = $vec_result_x.")
    vec_result_z = getvalue(vec_z)
    mat_result_mu = getvalue(mat_mu)
    vec_result_lambda = getvalue(vec_lambda)
    println("---------------------------------------------------------\n",
            "   2/2. Nominal Ending\n",
            "---------------------------------------------------------\n")
    return (obj_result, vec_result_y, vec_result_x, vec_result_z, mat_result_mu, vec_result_lambda)
end

# Function for LP with Box Uncertainty and Budget of Uncertainty
function lpBoxBudget(; num_x, vec_c, vec_b, mat_a, mat_a_bar, mat_a_hat, vec_b_bar, vec_gammaCap)
    m = length(mat_a_bar[:,1])
    n = length(mat_a_bar[1,:])
    for i = 1: m
        if vec_gammaCap[i] > n
            println("Error: vec_gammaCap[$i] is greater than m in nominal and hat data!")
        end
    end
    println("---------------------------------------------------------\n",
            "    1/2. Robust MILP with Box Uncertainty and Budget\n",
            "---------------------------------------------------------\n")
    println("Input: vec_gammaCap = $vec_gammaCap.")
    #
    model = Model(solver = GLPKSolverMIP())
    @variable(model, vec_x[1: num_x] >= 0)
    @objective(model, Min, (transpose(vec_c) * vec_x)[1])
    # Transformation of Box Uncertainty.
    @variable(model, vec_lambda[1: m] >= 0)
    @variable(model, mat_mu[1: m, 1: num_x] >= 0)
    @variable(model, vec_z[1: num_x] >= 0)
    @constraint(model, - vec_z .<= vec_x)
    @constraint(model, vec_x .<= vec_z)
    for i = 1: m
        # sum() of variables without coefficients can be used directly
        @constraint(model, transpose(mat_a_bar[i, :]) * vec_x + vec_gammaCap[i] * vec_lambda[i] +
            sum(mat_mu[i, :]) <= vec_b_bar[i])
        @constraint(model, vec_lambda[i] .+ hcat(mat_mu[i, :]) .>=  hcat(mat_a_hat[i, :]) .* vec_z[:])
    end
    solve(model)
    obj_result = getobjectivevalue(model)
    vec_result_x = getvalue(vec_x)
    # println("Result: vec_x = $vec_result_x.")
    vec_result_z = getvalue(vec_z)
    mat_result_mu = getvalue(mat_mu)
    vec_result_lambda = getvalue(vec_lambda)
    println("---------------------------------------------------------\n",
            "   2/2. Nominal Ending\n",
            "---------------------------------------------------------\n")
    return (obj_result, vec_result_x, vec_result_z, mat_result_mu, vec_result_lambda)
end

# Fucntion for Two-Stage Stochastic MILP with Box Uncertainty
function milpAdjustBox(; num_x, num_y, num_z, vec_min_y, vec_max_y, vec_c, vec_f, vec_g, vec_b,
    mat_a, mat_b, mat_d, mat_a_bar, mat_a_hat, mat_b_bar, mat_d_bar, vec_b_bar)
    println("---------------------------------------------------------\n",
            "   1/2. Adjustable Robust LP with Box Uncertainty\n",
            "---------------------------------------------------------\n")
    model = Model(solver = GLPKSolverMIP())
    # 1. Standard LP
    @variable(model, vec_y[1: num_y], Int)
    @constraint(model, vec_y >= vec_min_y)
    @constraint(model, vec_max_y >=  vec_y)
    @variable(model, vec_x[1: num_x] >= 0)
    @variable(model, gamma)
    @variable(model, vec_alpha[1: num_z])
    @variable(model, vec_theta1[1: num_z])
    @objective(model, Min, (transpose(vec_f) * vec_y + transpose(vec_c) * vec_x)[1] + gamma)
    @constraint(model, mat_b * vec_y + mat_a * vec_x + mat_d * (vec_alpha + vec_theta1) .>= vec_b)
    # 2. Get rid of the uncertainty in objective function
    @variable(model, vec_beta[1: num_z])
    @constraint(model, gamma >= (transpose(vec_g) * (vec_alpha + vec_theta1))[1])
    @constraint(model, vec_theta1 .<= vec_beta)
    @constraint(model, - vec_beta .<= vec_theta1)
    @constraint(model, vec_alpha .>= - vec_theta1)  # z(zeta) must be greater than 0
    # 3. Transformation of Box Uncertainty.
    @variable(model, vec_theta2[1: num_z])
    vec_theta2_rep = (ones(num_y, num_z) .* vec_theta2')[:]  # Repeat every theta2 num_y times and then concatenate.
    @constraint(model, mat_b_bar * vec_y + mat_a_bar * vec_x + mat_a_hat * vec_theta2_rep +
        mat_d_bar * (vec_alpha + vec_theta1) .<= vec_b_bar)
    @constraint(model, vec_theta2_rep .<= vec_x)
    @constraint(model, - vec_x .<= vec_theta2_rep)
    # @constraint(model, [m = 1: num_y], sum([mat_a_bar[m, (p-1) * num_y + m] * vec_x[(p-1) * num_y + m] + vec_theta2[p]]
    #     for p = 1: 10)[1] + mat_b_bar[m, m] * vec_y[m] <= vec_b_bar[m])
    # @constraint(model, [m = 1: num_y, p = 1: num_z], - mat_a_hat[m, (p-1) * num_y + m] * vec_x[(p-1) * num_y + m]
    #     <= vec_theta2[p])
    # @constraint(model, [m = 1: num_y, p = 1: num_z], vec_theta2[p] <=
    #     mat_a_hat[m, (p-1) * num_y + m] * vec_x[(p-1) * num_y + m])
    # Solve the model
    solve(model)
    obj_result = getobjectivevalue(model)
    vec_result_y = getvalue(vec_y)
    println("Result: vec_y = $vec_result_y.")
    vec_result_x = getvalue(vec_x)
    println("Result: vec_x = $vec_result_x.")
    vec_result_theta1 = getvalue(vec_theta1)
    vec_result_theta2 = getvalue(vec_theta2)
    vec_result_gamma = getvalue(gamma)
    println("---------------------------------------------------------\n",
            "   2/2. Nominal Ending\n",
            "---------------------------------------------------------\n")
    return (obj_result, vec_result_y, vec_result_x, vec_result_gamma, vec_result_theta1, vec_result_theta2)
end

# Fucntion for Two-Stage Stochastic LP with Box Uncertainty
function lpAdjustBox(; num_x, m, num_z, vec_c, vec_g, vec_b,
    mat_a, mat_d, mat_a_bar, mat_a_hat, mat_d_bar, vec_b_bar)
    println("---------------------------------------------------------\n",
            "   1/2. Adjustable Robust LP with Box Uncertainty\n",
            "---------------------------------------------------------\n")
    model = Model(solver = GLPKSolverMIP())
    # 1. Standard LP
    @variable(model, vec_x[1: num_x] >= 0)
    @variable(model, gamma)
    @variable(model, vec_alpha[1: num_z])
    @variable(model, vec_theta1[1: num_z])
    @objective(model, Min, (transpose(vec_c) * vec_x)[1] + gamma)
    @constraint(model, mat_a * vec_x + mat_d * (vec_alpha + vec_theta1) .>= vec_b)
    # 2. Get rid of the uncertainty in objective function
    @variable(model, vec_beta[1: num_z])
    @constraint(model, gamma >= (transpose(vec_g) * (vec_alpha + vec_theta1))[1])
    @constraint(model, vec_theta1 .<= vec_beta)
    @constraint(model, - vec_beta .<= vec_theta1)
    @constraint(model, vec_alpha .>= - vec_theta1)  # z(zeta) must be greater than 0
    # 3. Transformation of Box Uncertainty.
    @variable(model, vec_theta2[1: num_z])
    vec_theta2_rep = (ones(m, num_z) .* vec_theta2')[:]  # Repeat every theta2 num_y times and then concatenate.
    @constraint(model, mat_a_bar * vec_x + mat_a_hat * vec_theta2_rep +
        mat_d_bar * (vec_alpha + vec_theta1) .<= vec_b_bar)
    @constraint(model, vec_theta2_rep .<= vec_x)
    @constraint(model, - vec_x .<= vec_theta2_rep)
    # Solve the model
    solve(model)
    obj_result = getobjectivevalue(model)
    # println("Result: obj = $(obj).")
    vec_result_x = getvalue(vec_x)
    # println("Result: vec_x = $vec_result_x.")
    vec_result_theta1 = getvalue(vec_theta1)
    vec_result_theta2 = getvalue(vec_theta2)
    vec_result_gamma = getvalue(gamma)
    println("---------------------------------------------------------\n",
            "   2/2. Nominal Ending\n",
            "---------------------------------------------------------\n")
    return (obj_result, vec_result_x, vec_result_gamma, vec_result_theta1, vec_result_theta2)
end

# Fucntion for LP with Box Uncertainty
function lpBox(; num_x, m, num_z, vec_c, vec_b, mat_a, mat_a_bar, mat_a_hat, vec_b_bar)
    println("---------------------------------------------------------\n",
            "   1/2. Adjustable Robust LP with Box Uncertainty\n",
            "---------------------------------------------------------\n")
    model = Model(solver = GLPKSolverMIP())
    # 1. Standard LP
    @variable(model, vec_x[1: num_x] >= 0)
    @objective(model, Min, (transpose(vec_c) * vec_x)[1])
    @constraint(model, mat_a * vec_x .>= vec_b)
    # 2. Transformation of Box Uncertainty.
    @variable(model, vec_theta2[1: m])
    vec_theta2_rep = (ones(m, num_z) .* vec_theta2')[:]  # Repeat every theta2 num_y times and then concatenate.
    @constraint(model, mat_a_bar * vec_x + mat_a_hat * vec_theta2_rep .<= vec_b_bar)
    @constraint(model, vec_theta2_rep .<= vec_x)
    @constraint(model, - vec_x .<= vec_theta2_rep)
    # Solve the model
    solve(model)
    obj_result = getobjectivevalue(model)
    # println("Result: obj = $(obj).")
    vec_result_x = getvalue(vec_x)
    # println("Result: vec_x = $vec_result_x.")
    vec_result_theta2 = getvalue(vec_theta2)
    println("---------------------------------------------------------\n",
            "   2/2. Nominal Ending\n",
            "---------------------------------------------------------\n")
    return (obj_result, vec_result_x, vec_result_theta2)
end

end


text:
function getInteractionHamIPicture(
    Ham_S::Operator,
    Ham_int::Operator,
    t::AbstractFloat,
)::Operator
    U_op = evolutionOperator(Ham_S, t)
    return U_op' * Ham_int * U_op
end

function getInteractionHamIPictureA(
    Ham_S::Operator,
    Ham_int::Operator,
    t::AbstractFloat,
)::Array
    U_op = evolutionOperator(Ham_S, t)
    return (U_op' * Ham_int * U_op).data
end

function getInteractionHamIPictureA(
    Ham_int::Array,
    H_lambda::Array,
    H_S::Array,
    H_Sinv::Array,
    t::AbstractFloat,
)
    U_op = evolutionOperatorA(H_lambda, H_S, H_Sinv, t)
    U_opd = adjoint(U_op)
    return U_opd * Ham_int * U_op
end

ext: jl


text:
function generate_contact_expressions(model::Model, env::Environment;
		T = Float64, jacobians = false, nv = model.nq)

	# model dimensions
	nq = model.nq
	nu = model.nu
	nw = model.nw
	nc = model.nc

	# friction dimensions
	nb = nc * friction_dim(env)
	ncf = nc * dim(env)

	# simulator dimensions
	nz = num_var(model, env)
	nθ = num_data(model)

	# Declare variables
	@variables q0[1:nq]
	@variables q1[1:nq]
	@variables γ1[1:nc]
	@variables b1[1:nb]
	@variables λ1[1:ncf]
	@variables u1[1:nu]
	@variables w1[1:nw]
	@variables q2[1:nq]
	@variables vt[1:nb]
	@variables ψ1[1:nc]
	@variables η1[1:nb]
	@variables h[1:1]
	@variables k[1:ncf]

	@variables z[1:nz]
	@variables θ[1:nθ]
	@variables κ[1:1]

	# Expressions
	expr = Dict{Symbol, Expr}()

	# Contact Jacobian
	J = J_func(model, env, q2)
	J = reshape(J, (ncf, nv))
	J = Symbolics.simplify.(J)

	# Signed distance
	ϕ = ϕ_func(model, env, q2)
	ϕ = Symbolics.simplify.(ϕ)[1:nc]

	# Contact forces
	cf = contact_forces(model, env, γ1, b1, q2, k)
	cf = Symbolics.simplify.(cf)

	# Velocity stack
	vs = velocity_stack(model, env, q1, q2, k, h)
	vs = Symbolics.simplify.(vs)

	# Dynamics
	d = dynamics(model, env, h, q0, q1, u1, w1, λ1, q2)
	d = Symbolics.simplify.(d)

	# Functions
	expr[:J]    = build_function(J, q2, checkbounds=true)[1]
	expr[:ϕ]    = build_function(ϕ, q2, checkbounds=true)[1]
	expr[:cf]   = build_function(cf, γ1, b1, q2, k, checkbounds=true)[1]
	expr[:vs]   = build_function(vs, q1, q2, k, h, checkbounds=true)[1]
	expr[:d]    = build_function(d, h, q0, q1, u1, w1, λ1, q2, checkbounds=true)[1]

	if jacobians
		dJ = Symbolics.jacobian(J'*λ1, q2, simplify = true) # dq2 = ∂q2  + dJ(λ1, q2)
		dλ1 = J # this is an analytical result, might be wrong but should be okay

		dcf  = Symbolics.jacobian(cf, [q2; γ1; b1], simplify=true) # NOTE: input order change

		vsq2 = Symbolics.jacobian(-vs, q2, simplify = true)
		vsq1h = Symbolics.jacobian(-vs, [q1; h], simplify = true)

		# Maximum dissipation (eq.)
		ψ_stack = transpose(E_func(model, env)) * ψ1
		md = η1 - vt - ψ_stack
		md = Symbolics.simplify.(md)
		mdvs = Symbolics.jacobian(-md, vt, simplify = true)
		mdψη = Symbolics.jacobian(md, [ψ1; η1], simplify = true)

		# # Residual constraints
		rc = res_con(model, env, z, θ, κ)
		rc = Symbolics.simplify.(rc)
		rcz = Symbolics.jacobian(rc, z, simplify = true)
		rcθ = Symbolics.jacobian(rc, θ, simplify = true)

		# Functions
		expr[:dJ]   = build_function(dJ, λ1, q2, checkbounds=true)[1]
		expr[:dλ1]  = build_function(dλ1, q2, checkbounds=true)[1]

		expr[:dcf]  = build_function(dcf, γ1, b1, q2, k, checkbounds=true)[1]

		expr[:vsq2]  = build_function(vsq2, q1, q2, k, h, checkbounds=true)[1]
		expr[:vsq1h] = build_function(vsq1h, q1, q2, k, h, checkbounds=true)[1]

		expr[:mdvs] = build_function(mdvs, vt, ψ1, η1, checkbounds=true)[1]
		expr[:mdψη] = build_function(mdψη, vt, ψ1, η1, checkbounds=true)[1]

		expr[:rc]  = build_function(rc, z, θ, κ, checkbounds=true)[2]
		expr[:rcz] = build_function(rcz, z, θ, checkbounds=true)[2]
		expr[:rcθ] = build_function(rcθ, z, θ, checkbounds=true)[2]
	end
	return expr
end

"""
	generate_residual_expressions(model::Model)
Generate fast residual methods using Symbolics symbolic computing tools.
"""
function generate_residual_expressions(model::Model, env::Environment;
		mapping = (a,b,c) -> Diagonal(ones(num_var(model, env))), jacobians = :full,
		T = Float64, nv = model.nq)

	nq = model.nq
	nu = model.nu
	nc = model.nc
	nb = nc * friction_dim(env)
	nz = num_var(model, env)
	nθ = num_data(model)

	# Declare variables
	@variables z[1:nz]
	@variables θ[1:nθ]
	@variables κ[1:1]

	# Residual
	r = residual(model, env, z, θ, κ)
	r = Symbolics.simplify.(r)

	# Build function
	expr = Dict{Symbol, Expr}()
	expr[:r]  = build_function(r, z, θ, κ, checkbounds=true)[2]

	if jacobians == :full
		# contact expressions
		expr_contact = generate_contact_expressions(model, env,
			T = T, jacobians = false, nv = nv)

		m = mapping(model, env, z)
		m = simplify.(m)
		rz = Symbolics.jacobian(r, z, simplify = true)
		rz = rz * m
		rz = simplify.(rz)

		rθ = Symbolics.jacobian(r, θ, simplify = true) # TODO: sparse version

		rz_sp = similar(rz, T)
		rθ_sp = similar(rθ, T)

		expr[:rz] = build_function(rz, z, θ, checkbounds=true)[2]
		expr[:rθ] = build_function(rθ, z, θ, checkbounds=true)[2]
	else
		expr_contact = generate_contact_expressions(model, env,
			T = T, jacobians = true, nv = nv)

		rz_sp = zeros(nz, nz)
		rθ_sp = zeros(nz, nθ)
	end

	expr = merge(expr, expr_contact)

	return expr, rz_sp, rθ_sp
end

function instantiate_residual!(s::Simulation, path_res, path_jac;
	jacobians = :full)

	# load expressions
	expr = load_expressions(path_res)
	instantiate_contact_methods!(s.con, expr, jacobians = jacobians)

	# residual
	s.res.r!  = eval(expr[:r])

	# Jacobians
	if jacobians == :full
		s.res.rz! = eval(expr[:rz])
		s.res.rθ! = eval(expr[:rθ])
	else
		_rz(rz, z, θ) = rz_approx!(s, rz, z, θ)
		_rθ(rθ, z, θ) = rθ_approx!(s, rθ, z, θ)
		s.res.rz! = _rz
		s.res.rθ! = _rθ
	end

	# pre-allocated memory
	jac = load(path_jac)
	s.rz = jac["rz_sp"]
	s.rθ = jac["rθ_sp"]

	return nothing
end

"""
	instantiate_residual!(model,
		path::AbstractString="model/residual.jld2")
Evaluates the residual expressions to generate functions, stores them into the model.
"""
function instantiate_residual!(fct::ResidualMethods, expr::Dict{Symbol,Expr};
	jacobians = :full)

	fct.r!  = eval(expr[:r])

	if jacobians == :full
		fct.rz! = eval(expr[:rz])
		fct.rθ! = eval(expr[:rθ])
	else
		# fct.rz! = rz_approx!
		# fct.rθ! = rθ_approx!
	end

	return nothing
end

function instantiate_contact_methods!(fct::ContactMethods, expr::Dict{Symbol,Expr};
	jacobians = :full)

	fct.J = eval(expr[:J])
	fct.ϕ = eval(expr[:ϕ])
	fct.cf = eval(expr[:cf])
	fct.vs = eval(expr[:vs])
	fct.d = eval(expr[:d])

	if jacobians == :approx
		fct.dJ = eval(expr[:dJ])
		fct.dλ1 = eval(expr[:dλ1])
		fct.dcf = eval(expr[:dcf])
		fct.vsq2 = eval(expr[:vsq2])
		fct.vsq1h = eval(expr[:vsq1h])
		fct.mdvs = eval(expr[:mdvs])
		fct.mdψη = eval(expr[:mdψη])
		fct.rc = eval(expr[:rc])
		fct.rcz = eval(expr[:rcz])
		fct.rcθ = eval(expr[:rcθ])
	end

	return nothing
end


text:
using Documenter, BioSequences

makedocs(
    format = Documenter.HTML(),
    sitename = "BioSequences.jl",
    pages = [
        "Home"                           => "index.md",
        "Biological Symbols"             => "symbols.md",
        "BioSequences Types"             => "types.md",
        "Constructing sequences"         => "construction.md",
        "Indexing & modifying sequences" => "transforms.md",
        "Predicates"                     => "predicates.md",
        "Random sequences"               => "random.md",
        "Pattern matching and searching" => "sequence_search.md",
        "Sequence Demultipllexing"       => "demultiplexer.md",
        "Sequence composition"           => "composition.md",
        "Iteration"                      => "iteration.md",
        "Counting"                       => "counting.md",
        "I/O"                            => "io.md"
    ],
    authors = "Ben J. Ward, D.C.Jones, Kenta Sato, The BioJulia Organisation and other contributors."
)

deploydocs(
    repo = "github.com/BioJulia/BioSequences.jl.git",
    deps = nothing,
    make = nothing
)


text:
using AmbientForcing
using OrdinaryDiffEq
using Test
using LinearAlgebra

using BenchmarkTools

# Simple spherical dynamics

function g(du, u, p, t)
    du .= 0
    du[1] = 1.0
    du[end] = 1.0 - sum(u .^ 2)
    return nothing
end

N = 100
M = Diagonal([ones(N - 1); zeros(1)])
u0 = [zeros(N - 1); 1.0]
Frand = [1.0; zeros(N - 1)]

odeg = ODEFunction(g, mass_matrix=M)
prob = ODEProblem(odeg, u0, (0.0, 1.0))

solve(prob, Rodas4())
afoprob = ambient_forcing_problem(odeg, u0, 10.0, Frand)
sol = solve(afoprob, Tsit5(), save_everystep=false)
@btime solve(afoprob, Tsit5(), save_everystep=false)

afoprob_fd = ambient_forcing_problem(odeg, u0, 10.0, Frand, method=:ForwardDiff)
sol_fd = solve(afoprob_fd, Tsit5(), save_everystep=false)
@btime solve(afoprob_fd, Tsit5(), save_everystep=false)

@test Array(sol) == Array(sol_fd)


text:
function track(model::PopulationModel, Z::Array{T,2}, M::BitArray{2}) where T
    track(Z, weights(model), M)
end

# all subtypes of PopulationModel must have at least two methods
# PopulationModel(opts...)  # a constructor
# fit(p::Type{PopulationModel}, Z::Array{T,N), opts) # fit the model to an input

# For "spatial" models we might additionally have
# track(p::PopulationModel, Z)      # track the expression of populations over time
# bootstrap(p::PopulationModel, Z, fitopts)  # recalulate the model using a resampled Z
#                                   (applicable where Z is a 2-order tensor)

# ^ these methods should also work with the basic dimensionality reduction techniques in MultivariateStats

# Where Z is a 3rd or higher order tensor track & bootstrap will have to take a dimension


text:
function compute_update0!(H2, dHdtau, H, min_dxy2, _dx, _dy, dx, dy)
  @inbounds @fastmath for iy=1:size(dHdtau,2)
    for ix=1:size(dHdtau,1)
      dHdtau[ix, iy] = (-((-((H[(ix + 1) + 1, iy + 1] - H[ix + 1, iy + 1])) * _dx - -((H[ix + 1, iy + 1] - H[ix, iy + 1])) * _dx)) * _dx - (-((H[ix + 1, (iy + 1) + 1] - H[ix + 1, iy + 1])) * _dy - -((H[ix + 1, iy + 1] - H[ix + 1, iy])) * _dy) * _dy) + 0.6 * dHdtau[ix, iy]
      H2[ix+1,iy+1] = H[ix+1,iy+1] + min_dxy2/4.1*dHdtau[ix,iy]  # sets the BC as H[1]=H[end]=0
    end
  end
  return
end
function compute_update1!(H2, dHdtau, H, min_dxy2, _dx, _dy, dx, dy)
  @turbo for iy=1:size(dHdtau,2)
    for ix=1:size(dHdtau,1)
      dHdtau[ix, iy] = (-((-((H[(ix + 1) + 1, iy + 1] - H[ix + 1, iy + 1])) * _dx - -((H[ix + 1, iy + 1] - H[ix, iy + 1])) * _dx)) * _dx - (-((H[ix + 1, (iy + 1) + 1] - H[ix + 1, iy + 1])) * _dy - -((H[ix + 1, iy + 1] - H[ix + 1, iy])) * _dy) * _dy) + 0.6 * dHdtau[ix, iy]
      H2[ix+1,iy+1] = H[ix+1,iy+1] + min_dxy2/4.1*dHdtau[ix,iy]  # sets the BC as H[1]=H[end]=0
    end
  end
  return
end

@testset "conv" begin
  lx, ly = 10.0, 10.0   # domain size
  nx, ny = 12, 12     # number of grid points
  nt     = 1e3          # max number of iterations
  dx, dy = lx/nx, ly/ny
  xc, yc = LinRange(dx/2, lx-dx/2, nx), LinRange(dy/2, ly-dy/2, ny)
  dHdtau = zeros(nx-2, ny-2)
  dHdtau2 = zeros(nx-2, ny-2)
  H      = exp.(.-(xc.-lx/2).^2 .-(yc.-ly/2)'.^2)
  H1     = copy(H)
  H2     = copy(H)
  min_dxy2 = min(dx,dy)^2
  _dx, _dy = 1.0/dx, 1.0/dy
  compute_update0!(H1, dHdtau, H, min_dxy2, _dx, _dy, dx, dy)
  compute_update1!(H2, dHdtau2, H, min_dxy2, _dx, _dy, dx, dy)
  @test H1 ≈ H2
  @test dHdtau ≈ dHdtau2
end


text:
using Parameters
import Base: length,show,print,println
using StaticArrays
abstract type AbstractReaction end
export AbstractReaction

@with_kw struct ElementaryReaction{T<:AbstractRate,Q<:Integer,V1<:AbstractArray,V2<:AbstractArray,V3<:AbstractArray,V4<:AbstractArray,V5<:AbstractArray} <: AbstractReaction
    index::Q
    reactants::V1
    reactantinds::V2
    products::V3
    productinds::V4
    kinetics::T
    electronchange::Int64 = 0
    radicalchange::Int64 = -100
    reversible::Bool = true
    pairs::V5 = @SArray [@SArray [""]]
end
export ElementaryReaction

getrxnstr(rxn::T) where {T<:AbstractReaction} = join([join(getfield.(rxn.reactants,:name),"+"),join(getfield.(rxn.products,:name),"+")],"<=>")
export getrxnstr
show(io::IO,rxn::T) where {T<:AbstractReaction} = print(io,getrxnstr(rxn))
print(rxn::T) where {T<:AbstractReaction} = print(getrxnstr(rxn))
println(rxn::T) where {T<:AbstractReaction} = print(string(getrxnstr(rxn),"\n"))

"""
matches reactants to the products that most resemble them
returns a StaticArray of length 2 static arrays containing the names of each pair
"""
function getpairs(rxn::T) where {T<:AbstractReaction}
    if length(rxn.reactants) == 1
        if length(rxn.products) == 1
            return [ [rxn.reactants[1].name,rxn.products[1].name ] ]
        elseif length(rxn.products) == 2
            return [ [rxn.reactants[1].name,rxn.products[1].name ],  [rxn.reactants[1].name,rxn.products[2].name ]]
        elseif length(rxn.products) == 3
            return [ [rxn.reactants[1].name,rxn.products[1].name ], [rxn.reactants[1].name,rxn.products[2].name ], [rxn.reactants[1].name,rxn.products[3].name ]]
        end
    elseif length(rxn.reactants) == 2
        if length(rxn.products) == 1
            return [ [rxn.reactants[1].name,rxn.products[1].name ], [rxn.reactants[2].name,rxn.products[1].name ]]
        else
            return choosepairs(rxn)
        end
    elseif length(rxn.reactants) == 3
        if length(rxn.products) == 1
            return [ [rxn.reactants[1].name,rxn.products[1].name ], [rxn.reactants[2].name,rxn.products[1].name ], [rxn.reactants[3].name,rxn.products[1].name ]]
        else
            return choosepairs(rxn)
        end
    else
        return choosepairs(rxn)
    end
end

"""
returns a negative number representing compositional similarity between two species
this only uses the number of bonds and number of each atom type
"""
function getsimilarity(spc1::T,spc2::T2) where {T<:AbstractSpecies,T2<:AbstractSpecies}
    s = -abs(spc1.bondnum - spc2.bondnum)
    keys1 = keys(spc1.atomnums)
    keys2 = keys(spc2.atomnums)
    for x in keys1
        if x in keys2
            s -= abs(spc1.atomnums[x]-spc2.atomnums[x])
        else
            s -= spc1.atomnums[x]
        end
    end
    for y in keys2
        if y in keys1
            continue
        else
            s -= spc2.atomnums[y]
        end
    end
    return s
end

function choosepairs(rxn::T) where {T<:AbstractReaction}
    pairs = []
    for (i,r) in enumerate(rxn.reactants)
        ind = argmax([getsimilarity(r,p) for p in rxn.products])
        push!(pairs,[r.name,rxn.products[ind].name])
    end
    return pairs
end

length(r::T) where {T<:AbstractReaction}= 1
export length


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# Uniform random generation

## from types: rand(::Type, [dims...])

### GLOBAL_RNG fallback for all types

@inline rand(T::Type) = rand(GLOBAL_RNG, T)

### random floats

# CloseOpen(T) is the fallback for an AbstractFloat T
@inline rand(r::AbstractRNG=GLOBAL_RNG, ::Type{T}=Float64) where {T<:AbstractFloat} =
    rand(r, CloseOpen(T))

# generic random generation function which can be used by RNG implementors
# it is not defined as a fallback rand method as this could create ambiguities

rand_generic(r::AbstractRNG, ::CloseOpen{Float16}) =
    Float16(reinterpret(Float32,
                        (rand_ui10_raw(r) % UInt32 << 13) & 0x007fe000 | 0x3f800000) - 1)

rand_generic(r::AbstractRNG, ::CloseOpen{Float32}) =
    reinterpret(Float32, rand_ui23_raw(r) % UInt32 & 0x007fffff | 0x3f800000) - 1

rand_generic(r::AbstractRNG, ::Close1Open2_64) =
    reinterpret(Float64, 0x3ff0000000000000 | rand(r, UInt64) & 0x000fffffffffffff)

rand_generic(r::AbstractRNG, ::CloseOpen_64) = rand(r, Close1Open2()) - 1.0

#### BigFloat

const bits_in_Limb = sizeof(Limb) << 3
const Limb_high_bit = one(Limb) << (bits_in_Limb-1)

struct BigFloatRandGenerator
    prec::Int
    nlimbs::Int
    limbs::Vector{Limb}
    shift::UInt

    function BigFloatRandGenerator(prec::Int=precision(BigFloat))
        nlimbs = (prec-1) ÷ bits_in_Limb + 1
        limbs = Vector{Limb}(nlimbs)
        shift = nlimbs * bits_in_Limb - prec
        new(prec, nlimbs, limbs, shift)
    end
end

function _rand(rng::AbstractRNG, gen::BigFloatRandGenerator)
    z = BigFloat()
    limbs = gen.limbs
    rand!(rng, limbs)
    @inbounds begin
        limbs[1] <<= gen.shift
        randbool = iszero(limbs[end] & Limb_high_bit)
        limbs[end] |= Limb_high_bit
    end
    z.sign = 1
    unsafe_copy!(z.d, pointer(limbs), gen.nlimbs)
    (z, randbool)
end

function rand(rng::AbstractRNG, gen::BigFloatRandGenerator, ::Close1Open2{BigFloat})
    z = _rand(rng, gen)[1]
    z.exp = 1
    z
end

function rand(rng::AbstractRNG, gen::BigFloatRandGenerator, ::CloseOpen{BigFloat})
    z, randbool = _rand(rng, gen)
    z.exp = 0
    randbool &&
        ccall((:mpfr_sub_d, :libmpfr), Int32,
              (Ptr{BigFloat}, Ptr{BigFloat}, Cdouble, Int32),
              &z, &z, 0.5, Base.MPFR.ROUNDING_MODE[])
    z
end

# alternative, with 1 bit less of precision
# TODO: make an API for requesting full or not-full precision
function rand(rng::AbstractRNG, gen::BigFloatRandGenerator, ::CloseOpen{BigFloat}, ::Void)
    z = rand(rng, Close1Open2(BigFloat), gen)
    ccall((:mpfr_sub_ui, :libmpfr), Int32, (Ptr{BigFloat}, Ptr{BigFloat}, Culong, Int32),
          &z, &z, 1, Base.MPFR.ROUNDING_MODE[])
    z
end

rand_generic(rng::AbstractRNG, I::FloatInterval{BigFloat}) =
    rand(rng, BigFloatRandGenerator(), I)

### random integers

rand_ui10_raw(r::AbstractRNG) = rand(r, UInt16)
rand_ui23_raw(r::AbstractRNG) = rand(r, UInt32)

@inline rand_ui52_raw(r::AbstractRNG) = reinterpret(UInt64, rand(r, Close1Open2()))
@inline rand_ui52(r::AbstractRNG) = rand_ui52_raw(r) & 0x000fffffffffffff

### random complex numbers

rand(r::AbstractRNG, ::Type{Complex{T}}) where {T<:Real} = complex(rand(r, T), rand(r, T))

### random characters

# returns a random valid Unicode scalar value (i.e. 0 - 0xd7ff, 0xe000 - # 0x10ffff)
function rand(r::AbstractRNG, ::Type{Char})
    c = rand(r, 0x00000000:0x0010f7ff)
    (c < 0xd800) ? Char(c) : Char(c+0x800)
end

### arrays of random numbers

function rand!(r::AbstractRNG, A::AbstractArray{T}, ::Type{X}=T) where {T,X}
    for i in eachindex(A)
        @inbounds A[i] = rand(r, X)
    end
    A
end

rand!(A::AbstractArray, ::Type{X}) where {X} = rand!(GLOBAL_RNG, A, X)
# NOTE: if the second parameter above is defaulted to eltype(A) and the
# method below is removed, then some specialized methods (e.g. for
# rand!(::Array{Float64})) will fail to be called
rand!(A::AbstractArray) = rand!(GLOBAL_RNG, A)

rand(r::AbstractRNG, dims::Dims)       = rand(r, Float64, dims)
rand(                dims::Dims)       = rand(GLOBAL_RNG, dims)
rand(r::AbstractRNG, dims::Integer...) = rand(r, Dims(dims))
rand(                dims::Integer...) = rand(Dims(dims))

rand(r::AbstractRNG, T::Type, dims::Dims)                   = rand!(r, Array{T}(dims))
rand(                T::Type, dims::Dims)                   = rand(GLOBAL_RNG, T, dims)
rand(r::AbstractRNG, T::Type, d::Integer, dims::Integer...) = rand(r, T, Dims((d, dims...)))
rand(                T::Type, d::Integer, dims::Integer...) = rand(T, Dims((d, dims...)))
# note: the above methods would trigger an ambiguity warning if d was not separated out:
# rand(r, ()) would match both this method and rand(r, dims::Dims)
# moreover, a call like rand(r, NotImplementedType()) would be an infinite loop

#### arrays of floats

rand!(r::AbstractRNG, A::AbstractArray, ::Type{T}) where {T<:AbstractFloat} =
    rand!(r, A, CloseOpen{T}())

function rand!(r::AbstractRNG, A::AbstractArray, I::FloatInterval)
    for i in eachindex(A)
        @inbounds A[i] = rand(r, I)
    end
    A
end

function rand!(rng::AbstractRNG, A::AbstractArray, I::FloatInterval{BigFloat})
    gen = BigFloatRandGenerator()
    for i in eachindex(A)
        @inbounds A[i] = rand(rng, gen, I)
    end
    A
end

rand!(A::AbstractArray, I::FloatInterval) = rand!(GLOBAL_RNG, A, I)

## Generate random integer within a range

abstract type RangeGenerator end

### RangeGenerator for BitInteger

# remainder function according to Knuth, where rem_knuth(a, 0) = a
rem_knuth(a::UInt, b::UInt) = a % (b + (b == 0)) + a * (b == 0)
rem_knuth(a::T, b::T) where {T<:Unsigned} = b != 0 ? a % b : a

# maximum multiple of k <= 2^bits(T) decremented by one,
# that is 0xFFFF...FFFF if k = typemax(T) - typemin(T) with intentional underflow
# see http://stackoverflow.com/questions/29182036/integer-arithmetic-add-1-to-uint-max-and-divide-by-n-without-overflow
maxmultiple(k::T) where {T<:Unsigned} =
    (div(typemax(T) - k + oneunit(k), k + (k == 0))*k + k - oneunit(k))::T

# maximum multiple of k within 1:2^32 or 1:2^64 decremented by one, depending on size
maxmultiplemix(k::UInt64) = k >> 32 != 0 ?
    maxmultiple(k) :
    (div(0x0000000100000000, k + (k == 0))*k - oneunit(k))::UInt64

struct RangeGeneratorInt{T<:Integer,U<:Unsigned} <: RangeGenerator
    a::T   # first element of the range
    k::U   # range length or zero for full range
    u::U   # rejection threshold
end

# generators with 32, 128 bits entropy
RangeGeneratorInt(a::T, k::U) where {T,U<:Union{UInt32,UInt128}} =
    RangeGeneratorInt{T,U}(a, k, maxmultiple(k))

# mixed 32/64 bits entropy generator
RangeGeneratorInt(a::T, k::UInt64) where {T} =
    RangeGeneratorInt{T,UInt64}(a, k, maxmultiplemix(k))

function RangeGenerator(r::UnitRange{T}) where T<:Unsigned
    isempty(r) && throw(ArgumentError("range must be non-empty"))
    RangeGeneratorInt(first(r), last(r) - first(r) + oneunit(T))
end

for (T, U) in [(UInt8, UInt32), (UInt16, UInt32),
               (Int8, UInt32), (Int16, UInt32), (Int32, UInt32),
               (Int64, UInt64), (Int128, UInt128), (Bool, UInt32)]

    @eval RangeGenerator(r::UnitRange{$T}) = begin
        isempty(r) && throw(ArgumentError("range must be non-empty"))
        # overflow ok:
        RangeGeneratorInt(first(r), convert($U, unsigned(last(r) - first(r)) + one($U)))
    end
end

### RangeGenerator for BigInt

struct RangeGeneratorBigInt <: RangeGenerator
    a::BigInt         # first
    m::BigInt         # range length - 1
    nlimbs::Int       # number of limbs in generated BigInt's (z ∈ [0, m])
    nlimbsmax::Int    # max number of limbs for z+a
    mask::Limb        # applied to the highest limb
end

function RangeGenerator(r::UnitRange{BigInt})
    m = last(r) - first(r)
    m < 0 && throw(ArgumentError("range must be non-empty"))
    nd = ndigits(m, 2)
    nlimbs, highbits = divrem(nd, 8*sizeof(Limb))
    highbits > 0 && (nlimbs += 1)
    mask = highbits == 0 ? ~zero(Limb) : one(Limb)<<highbits - one(Limb)
    nlimbsmax = max(nlimbs, abs(last(r).size), abs(first(r).size))
    return RangeGeneratorBigInt(first(r), m, nlimbs, nlimbsmax, mask)
end

### rand(::RangeGenerator)

# this function uses 32 bit entropy for small ranges of length <= typemax(UInt32) + 1
# RangeGeneratorInt is responsible for providing the right value of k
function rand(rng::AbstractRNG, g::RangeGeneratorInt{T,UInt64}) where T<:Union{UInt64,Int64}
    local x::UInt64
    if (g.k - 1) >> 32 == 0
        x = rand(rng, UInt32)
        while x > g.u
            x = rand(rng, UInt32)
        end
    else
        x = rand(rng, UInt64)
        while x > g.u
            x = rand(rng, UInt64)
        end
    end
    return reinterpret(T, reinterpret(UInt64, g.a) + rem_knuth(x, g.k))
end

function rand(rng::AbstractRNG, g::RangeGeneratorInt{T,U}) where {T<:Integer,U<:Unsigned}
    x = rand(rng, U)
    while x > g.u
        x = rand(rng, U)
    end
    (unsigned(g.a) + rem_knuth(x, g.k)) % T
end

function rand(rng::AbstractRNG, g::RangeGeneratorBigInt)
    x = MPZ.realloc2(g.nlimbsmax*8*sizeof(Limb))
    limbs = unsafe_wrap(Array, x.d, g.nlimbs)
    while true
        rand!(rng, limbs)
        @inbounds limbs[end] &= g.mask
        MPZ.mpn_cmp(x, g.m, g.nlimbs) <= 0 && break
    end
    # adjust x.size (normally done by mpz_limbs_finish, in GMP version >= 6)
    x.size = g.nlimbs
    while x.size > 0
        @inbounds limbs[x.size] != 0 && break
        x.size -= 1
    end
    MPZ.add!(x, g.a)
end

#### arrays

function rand!(rng::AbstractRNG, A::AbstractArray, g::RangeGenerator)
    for i in eachindex(A)
        @inbounds A[i] = rand(rng, g)
    end
    return A
end

### random values from UnitRange

rand(rng::AbstractRNG, r::UnitRange{<:Union{Signed,Unsigned,BigInt,Bool}}) =
    rand(rng, RangeGenerator(r))

rand!(rng::AbstractRNG, A::AbstractArray,
      r::UnitRange{<:Union{Signed,Unsigned,BigInt,Bool,Char}}) =
          rand!(rng, A, RangeGenerator(r))

## random values from AbstractArray

rand(rng::AbstractRNG, r::AbstractArray) = @inbounds return r[rand(rng, 1:length(r))]
rand(                  r::AbstractArray) = rand(GLOBAL_RNG, r)

### arrays

function rand!(rng::AbstractRNG, A::AbstractArray, r::AbstractArray)
    g = RangeGenerator(1:(length(r)))
    for i in eachindex(A)
        @inbounds A[i] = r[rand(rng, g)]
    end
    return A
end

rand!(A::AbstractArray, r::AbstractArray) = rand!(GLOBAL_RNG, A, r)

rand(rng::AbstractRNG, r::AbstractArray{T}, dims::Dims) where {T} =
    rand!(rng, Array{T}(dims), r)
rand(                  r::AbstractArray, dims::Dims)       = rand(GLOBAL_RNG, r, dims)
rand(rng::AbstractRNG, r::AbstractArray, dims::Integer...) = rand(rng, r, Dims(dims))
rand(                  r::AbstractArray, dims::Integer...) = rand(GLOBAL_RNG, r, Dims(dims))

## random values from Dict, Set, IntSet

function rand(r::AbstractRNG, t::Dict)
    isempty(t) && throw(ArgumentError("collection must be non-empty"))
    rg = RangeGenerator(1:length(t.slots))
    while true
        i = rand(r, rg)
        Base.isslotfilled(t, i) && @inbounds return (t.keys[i] => t.vals[i])
    end
end

rand(r::AbstractRNG, s::Set) = rand(r, s.dict).first

function rand(r::AbstractRNG, s::IntSet)
    isempty(s) && throw(ArgumentError("collection must be non-empty"))
    # s can be empty while s.bits is not, so we cannot rely on the
    # length check in RangeGenerator below
    rg = RangeGenerator(1:length(s.bits))
    while true
        n = rand(r, rg)
        @inbounds b = s.bits[n]
        b && return n
    end
end

function nth(iter, n::Integer)::eltype(iter)
    for (i, x) in enumerate(iter)
        i == n && return x
    end
end
nth(iter::AbstractArray, n::Integer) = iter[n]

rand(r::AbstractRNG, s::Union{Associative,AbstractSet}) = nth(s, rand(r, 1:length(s)))

rand(s::Union{Associative,AbstractSet}) = rand(GLOBAL_RNG, s)

### arrays

function rand!(r::AbstractRNG, A::AbstractArray, s::Union{Dict,Set,IntSet})
    for i in eachindex(A)
        @inbounds A[i] = rand(r, s)
    end
    A
end

# avoid linear complexity for repeated calls with generic containers
rand!(r::AbstractRNG, A::AbstractArray, s::Union{Associative,AbstractSet}) =
    rand!(r, A, collect(s))

rand!(A::AbstractArray, s::Union{Associative,AbstractSet}) = rand!(GLOBAL_RNG, A, s)

rand(r::AbstractRNG, s::Associative{K,V}, dims::Dims) where {K,V} =
    rand!(r, Array{Pair{K,V}}(dims), s)

rand(r::AbstractRNG, s::AbstractSet{T}, dims::Dims) where {T} = rand!(r, Array{T}(dims), s)
rand(r::AbstractRNG, s::Union{Associative,AbstractSet}, dims::Integer...) =
    rand(r, s, Dims(dims))
rand(s::Union{Associative,AbstractSet}, dims::Integer...) = rand(GLOBAL_RNG, s, Dims(dims))
rand(s::Union{Associative,AbstractSet}, dims::Dims) = rand(GLOBAL_RNG, s, dims)

## random characters from a string

isvalid_unsafe(s::String, i) = !Base.is_valid_continuation(unsafe_load(pointer(s), i))
isvalid_unsafe(s::AbstractString, i) = isvalid(s, i)
_endof(s::String) = sizeof(s)
_endof(s::AbstractString) = endof(s)

function rand(rng::AbstractRNG, s::AbstractString)::Char
    g = RangeGenerator(1:_endof(s))
    while true
        pos = rand(rng, g)
        isvalid_unsafe(s, pos) && return s[pos]
    end
end

rand(s::AbstractString) = rand(GLOBAL_RNG, s)

### arrays

# we use collect(str), which is most of the time more efficient than specialized methods
# (except maybe for very small arrays)
rand!(rng::AbstractRNG, A::AbstractArray, str::AbstractString) = rand!(rng, A, collect(str))
rand!(A::AbstractArray, str::AbstractString) = rand!(GLOBAL_RNG, A, str)
rand(rng::AbstractRNG, str::AbstractString, dims::Dims) =
    rand!(rng, Array{eltype(str)}(dims), str)

rand(rng::AbstractRNG, str::AbstractString, d::Integer, dims::Integer...) =
    rand(rng, str, Dims((d, dims...)))

rand(str::AbstractString, dims::Dims) = rand(GLOBAL_RNG, str, dims)
rand(str::AbstractString, d::Integer, dims::Integer...) = rand(GLOBAL_RNG, str, d, dims...)

max_stars_repo_name: baggepinnen/julia

max_issues_repo_name: baggepinnen/julia

max_forks_repo_name: afniedermayer/julia


text:
#########################################################################
#
#   Base rules
#
#########################################################################

# derivation neutral functions
@deriv_rule colon(x,y)   x     0.
@deriv_rule colon(x,y)   y     0.

@deriv_rule length(x)    x     0.

@deriv_rule size(x)      x     0.
@deriv_rule size(x,y)    x     0.
@deriv_rule size(x,y)    y     0.

@deriv_rule fill(x,y)    x     0.
@deriv_rule fill(x,y)    y     0.

@deriv_rule similar(x,y) x     0.
@deriv_rule similar(x,y) y     0.

@deriv_rule zeros(x)     x     0.

@deriv_rule ones(x)      x     0.

@deriv_rule cell(x)      x     0.

@deriv_rule sign(x)      x     0.

@deriv_rule reverse(x)   x     0.

#  tuple  TODO : specialized macro for this kind of function
@deriv_rule tuple(x)        x     ds[1]
@deriv_rule tuple(x,y)      x     ds[1]
@deriv_rule tuple(x,y)      y     ds[2]
@deriv_rule tuple(x,y,z)    x     ds[1]
@deriv_rule tuple(x,y,z)    y     ds[2]
@deriv_rule tuple(x,y,z)    z     ds[3]
@deriv_rule tuple(x,y,z,t)  x     ds[1]
@deriv_rule tuple(x,y,z,t)  y     ds[2]
@deriv_rule tuple(x,y,z,t)  x     ds[3]
@deriv_rule tuple(x,y,z,t)  t     ds[4]

#  vcat
@deriv_rule vcat(x)        x     ds[1]
@deriv_rule vcat(x,y)      x     ds[1]
@deriv_rule vcat(x,y)      y     ds[2]
@deriv_rule vcat(x,y,z)    x     ds[1]
@deriv_rule vcat(x,y,z)    y     ds[2]
@deriv_rule vcat(x,y,z)    z     ds[3]
@deriv_rule vcat(x,y,z,t)  x     ds[1]
@deriv_rule vcat(x,y,z,t)  y     ds[2]
@deriv_rule vcat(x,y,z,t)  x     ds[3]
@deriv_rule vcat(x,y,z,t)  t     ds[4]

# reshape
@deriv_rule reshape(x::AbstractArray, a, b)        x    reshape(ds, size(x))
@deriv_rule reshape(x::AbstractArray, a, b)        a    0.
@deriv_rule reshape(x::AbstractArray, a, b)        b    0.
@deriv_rule reshape(x::AbstractArray, d::Tuple)    x    reshape(ds, size(x))
@deriv_rule reshape(x::AbstractArray, d::Tuple)    d    0.

# square root
@deriv_rule sqrt(x::Real)              x     0.5 * x ^ (-0.5) * ds
@deriv_rule sqrt(x::AbstractVector)    x     0.5 .* x .^ (-0.5) .* ds

# addition
@deriv_rule +(x::Real         , y::Real )            x     ds
@deriv_rule +(x::AbstractArray, y::AbstractArray)    x     ds
@deriv_rule +(x::Real         , y::Real )            y     ds
@deriv_rule +(x::AbstractArray, y::AbstractArray)    y     ds

@deriv_rule +(x::Real         , y::AbstractArray)    x     sum(ds)
@deriv_rule +(x::AbstractArray, y       )            x     ds
@deriv_rule +(x::AbstractArray, y::Real )            y     sum(ds)
@deriv_rule +(x               , y::AbstractArray)    y     ds

# dot addition
@deriv_rule .+(x::Real         , y::Real )            x     ds
@deriv_rule .+(x::Real         , y::AbstractArray)    x     sum(ds)
@deriv_rule .+(x::AbstractArray, y       )            x     ds
@deriv_rule .+(x::Real         , y::Real )            y     ds
@deriv_rule .+(x::AbstractArray, y::Real )            y     sum(ds)
@deriv_rule .+(x               , y::AbstractArray)    y     ds

# unary substraction
@deriv_rule -(x::Real )                              x     -ds
@deriv_rule -(x::AbstractArray)                      x     -ds

# binary substraction
@deriv_rule -(x::Real         , y::Real )            x     ds
@deriv_rule -(x::AbstractArray, y::AbstractArray)    x     ds
@deriv_rule -(x::Real         , y::Real )            y     -ds
@deriv_rule -(x::AbstractArray, y::AbstractArray)    y     -ds

# dot binary substraction
@deriv_rule .-(x::Real         , y::Real )            x     ds
@deriv_rule .-(x::Real         , y::AbstractArray)    x     sum(ds)
@deriv_rule .-(x::AbstractArray, y       )            x     ds
@deriv_rule .-(x::Real         , y::Real )            y     -ds
@deriv_rule .-(x::AbstractArray, y::Real )            y     -sum(ds)
@deriv_rule .-(x               , y::AbstractArray)    y     -ds

# sum()
@deriv_rule sum(x::Real )                            x     ds
@deriv_rule sum(x::AbstractArray)                    x     ones(size(x)).*ds

# dot()
@deriv_rule dot(x::Real         , y::Real )          x     y * ds
@deriv_rule dot(x::Real         , y::Real )          y     x * ds

@deriv_rule dot(x::AbstractArray, y::AbstractArray)  x     y.*ds
@deriv_rule dot(x::AbstractArray, y::AbstractArray)  y     x.*ds

# log() and exp()
@deriv_rule log(x::Real )                            x     ds / x
@deriv_rule log(x::AbstractArray)                    x     ds ./ x

@deriv_rule exp(x::Real )                            x     exp(x) * ds
@deriv_rule exp(x::AbstractArray)                    x     exp(x) .* ds

@deriv_rule log1p(x::Real)                           x     ds  / (1 + x)
@deriv_rule log1p(x::AbstractArray)                  x     ds ./ (1 + x)

@deriv_rule expm1(x::Real)                           x     (1. + expm1(x))  * ds
@deriv_rule expm1(x::AbstractArray)                  x     (1. + expm1(x)) .* ds
# note : derivative uses expm1() and not exp() to reuse the
#   already calculated expm1()

# trig functions
@deriv_rule sin(x::Real )                            x     cos(x) * ds
@deriv_rule sin(x::AbstractArray)                    x     cos(x) .* ds

@deriv_rule cos(x::Real )                            x     -sin(x) * ds
@deriv_rule cos(x::AbstractArray)                    x     -sin(x) .* ds

@deriv_rule tan(x::Real )                            x     (1. + tan(x)  * tan(x))  * ds
@deriv_rule tan(x::AbstractArray)                    x     (1. + tan(x) .* tan(x)) .* ds

@deriv_rule sinh(x::Real )                           x     cosh(x) * ds
@deriv_rule sinh(x::AbstractArray)                   x     cosh(x) .* ds

@deriv_rule cosh(x::Real )                           x     sinh(x) * ds
@deriv_rule cosh(x::AbstractArray)                   x     sinh(x) .* ds

@deriv_rule tanh(x::Real )                           x     (1. - tanh(x)  * tanh(x))  * ds
@deriv_rule tanh(x::AbstractArray)                   x     (1. - tanh(x) .* tanh(x)) .* ds

@deriv_rule asin(x::Real )                           x     ds  / sqrt(1 - x *x)
@deriv_rule asin(x::AbstractArray)                   x     ds ./ sqrt(1 - x.*x)

@deriv_rule acos(x::Real )                           x     -ds  / sqrt(1 - x *x)
@deriv_rule acos(x::AbstractArray)                   x     -ds ./ sqrt(1 - x.*x)

@deriv_rule atan(x::Real )                           x     ds  / (1 + x *x)
@deriv_rule atan(x::AbstractArray)                   x     ds ./ (1 + x.*x)

# round, floor, ceil, trunc, mod2pi
@deriv_rule round(x::Real)                           x     0.
@deriv_rule round(x::AbstractArray)                  x     0.

@deriv_rule floor(x::Real)                           x     0.
@deriv_rule floor(x::AbstractArray)                  x     0.

@deriv_rule ceil(x::Real)                            x     0.
@deriv_rule ceil(x::AbstractArray)                   x     0.

@deriv_rule trunc(x::Real)                           x     0.
@deriv_rule trunc(x::AbstractArray)                  x     0.

@deriv_rule mod2pi(x::Real)                          x     ds

# abs, max(), min()
@deriv_rule abs(x::Real )                            x     sign(x) * ds
@deriv_rule abs(x::AbstractArray)                    x     sign(x) .* ds

@deriv_rule max(x::Real         , y::Real )          x     (x > y) * ds
@deriv_rule max(x::Real         , y::AbstractArray)  x     sum((x .> y) .* ds)
@deriv_rule max(x::AbstractArray, y::Real )          x     (x .> y) .* ds
@deriv_rule max(x::AbstractArray, y::AbstractArray)  x     (x .> y) .* ds

@deriv_rule max(x::Real         , y::Real )          y     (x < y) * ds
@deriv_rule max(x::Real         , y::AbstractArray)  y     (x .< y) .* ds
@deriv_rule max(x::AbstractArray, y::Real )          y     sum((x .< y) .* ds)
@deriv_rule max(x::AbstractArray, y::AbstractArray)  y     (x .< y) .* ds

@deriv_rule min(x::Real         , y::Real )          x     (x < y) * ds
@deriv_rule min(x::Real         , y::AbstractArray)  x     sum((x .< y) .* ds)
@deriv_rule min(x::AbstractArray, y::Real )          x     (x .< y) .* ds
@deriv_rule min(x::AbstractArray, y::AbstractArray)  x     (x .< y) .* ds

@deriv_rule min(x::Real         , y::Real )          y     (x > y) * ds
@deriv_rule min(x::Real         , y::AbstractArray)  y     (x .> y) .* ds
@deriv_rule min(x::AbstractArray, y::Real )          y     sum((x .> y) .* ds)
@deriv_rule min(x::AbstractArray, y::AbstractArray)  y     (x .> y) .* ds

# maximum, minimum
@deriv_rule maximum(x::Real         )     x     ds
@deriv_rule maximum(x::AbstractArray)     x     (x .== maximum(x)) .* ds

@deriv_rule minimum(x::Real         )     x     ds
@deriv_rule minimum(x::AbstractArray)     x     (x .== minimum(x)) .* ds

# multiplication
@deriv_rule *(x::Real         , y::Real )            x     y * ds
@deriv_rule *(x::Real         , y::AbstractArray)    x     sum(y .* ds)
@deriv_rule *(x::AbstractArray, y::Real )            x     y .* ds
@deriv_rule *(x::AbstractArray, y::AbstractArray)    x     ds * y'

@deriv_rule *(x::Real         , y::Real )            y     x * ds
@deriv_rule *(x::Real         , y::AbstractArray)    y     x .* ds
@deriv_rule *(x::AbstractArray, y::Real )            y     sum(x .* ds)
@deriv_rule *(x::AbstractArray, y::AbstractArray)    y     x' * ds

# dot multiplication
@deriv_rule .*(x::Real         , y::Real )           x     y .* ds
@deriv_rule .*(x::Real         , y::AbstractArray)   x     sum(y .* ds)
@deriv_rule .*(x::AbstractArray, y::Real )           x     y .* ds
@deriv_rule .*(x::AbstractArray, y::AbstractArray)   x     y .* ds

@deriv_rule .*(x::Real         , y::Real )           y     x * ds
@deriv_rule .*(x::Real         , y::AbstractArray)   y     x .* ds
@deriv_rule .*(x::AbstractArray, y::Real )           y     sum(x .* ds)
@deriv_rule .*(x::AbstractArray, y::AbstractArray)   y     x .* ds

# power  (both args reals)
@deriv_rule ^(x::Real, y::Real)                      x     y * x ^ (y-1) * ds
@deriv_rule ^(x::Real, y::Real)                      y     log(x) * x ^ y * ds

# dot power
@deriv_rule .^(x::Real         , y::Real )           x     y * x ^ (y-1) * ds
@deriv_rule .^(x::Real         , y::AbstractArray)   x     sum(y .* x .^ (y-1) .* ds)
@deriv_rule .^(x::AbstractArray, y::Real )           x     y * x .^ (y-1) .* ds
@deriv_rule .^(x::AbstractArray, y::AbstractArray)   x     y .* x .^ (y-1) .* ds

@deriv_rule .^(x::Real         , y::Real )           y     log(x) * x ^ y * ds
@deriv_rule .^(x::AbstractArray, y::Real )           y     sum( log(x) .* x .^ y .* ds)
@deriv_rule .^(x::Real         , y::AbstractArray)   y     log(x) .* x .^ y .* ds
@deriv_rule .^(x::AbstractArray, y::AbstractArray)   y     log(x) .* x .^ y .* ds

# division
@deriv_rule /(x::Real          , y::Real )           x     ds / y
@deriv_rule /(x::AbstractArray , y::Real )           x     ds ./ y

@deriv_rule /(x::Real          , y::Real )           y     -x * ds / (y * y)
@deriv_rule /(x::AbstractArray , y::Real )           y     sum(-x .* ds) / (y * y)

# dot division
@deriv_rule ./(x::Real         , y::Real )           x     ds / y
@deriv_rule ./(x::Real         , y::AbstractArray)   x     sum(ds ./ y)
@deriv_rule ./(x::AbstractArray, y::Real )           x     ds ./ y
@deriv_rule ./(x::AbstractArray, y::AbstractArray)   x     ds ./ y

@deriv_rule ./(x::Real         , y::Real )           y     -x * ds / (y * y)
@deriv_rule ./(x::Real         , y::AbstractArray)   y     -x * ds ./ (y .* y)
@deriv_rule ./(x::AbstractArray, y::Real )           y     -sum(x .* ds) / (y * y)
@deriv_rule ./(x::AbstractArray, y::AbstractArray)   y     -x .* ds ./ (y .* y)

# transpose
@deriv_rule transpose(x::Real )                      x     ds
@deriv_rule transpose(x::AbstractArray)              x     transpose(ds)

# erf, erfc, gamma, beta, lbeta, lgamma
@deriv_rule erf(x::Real)                       x     2/sqrt(π) * exp(-x  * x)  * ds
@deriv_rule erf(x::AbstractArray)              x     2/sqrt(π) * exp(-x .* x) .* ds

@deriv_rule erfc(x::Real)                      x     -2/sqrt(π) * exp(-x  * x)  * ds
@deriv_rule erfc(x::AbstractArray)             x     -2/sqrt(π) * exp(-x .* x) .* ds

@deriv_rule gamma(x::Real)                     x     polygamma(0,x)  * gamma(x)  * ds
@deriv_rule gamma(x::AbstractArray)            x     polygamma(0,x) .* gamma(x) .* ds

@deriv_rule lgamma(x::Real)                    x     polygamma(0,x)  * ds
@deriv_rule lgamma(x::AbstractArray)           x     polygamma(0,x) .* ds

@deriv_rule beta(x::Real         , y::Real)            x   beta(x,y)  * (digamma(x)-digamma(x+y))  * ds
@deriv_rule beta(x::AbstractArray, y::AbstractArray)   x   beta(x,y) .* (digamma(x)-digamma(x+y)) .* ds
@deriv_rule beta(x::Real         , y::Real)            y   beta(x,y)  * (digamma(y)-digamma(x+y))  * ds
@deriv_rule beta(x::AbstractArray, y::AbstractArray)   y   beta(x,y) .* (digamma(y)-digamma(x+y)) .* ds

@deriv_rule lbeta(x::Real         , y::Real)            x   (polygamma(0,x)-polygamma(0,x+y))  * ds
@deriv_rule lbeta(x::AbstractArray, y::AbstractArray)   x   (polygamma(0,x)-polygamma(0,x+y)) .* ds
@deriv_rule lbeta(x::Real         , y::Real)            y   (polygamma(0,y)-polygamma(0,x+y))  * ds
@deriv_rule lbeta(x::AbstractArray, y::AbstractArray)   y   (polygamma(0,y)-polygamma(0,x+y)) .* ds


text:
struct MidpointPropagator{E<:Exponentiator} <: MagnusPropagator
    A::Function
    a::Number
    Exp::E
end
MidpointPropagator(A::Function, Exp::E) where E<:Exponentiator =
    MidpointPropagator(A, 1, Exp)
MidpointPropagator(B::LinearMap,
                   f::Function,
                   C::LinearMap,
                   a::Number,
                   Exp::E) where E<:Exponentiator =
                       MidpointPropagator(t -> B + f(t)*C, a, Exp)

(p::MidpointPropagator{E})(t::Real, τ::Real,
                           v::AbstractVector, w::AbstractVector) where E<:Exponentiator =
                               p.Exp(p.A(t+τ/2), p.a*τ, v, w)

export MidpointPropagator


text:
export apply!, pointcrop!, linecrop!, linepool!, heatmappool!, aggregate!

function apply!(f::Function, buf::AbstractMatrix, c::Canvas)
  for j ∈ 1:size(buf,2)
    for i ∈ 1:size(buf,1)
      x, y = ij2xy(i, j, c)
      buf[i,j] = f(x, y)
    end
  end
end

function apply!(f::Function, buf::AbstractVector{Point2f0}, c::Canvas)
  xyrect = c.xyrect[]
  n = width(c.ijhome)
  p = 0
  for x ∈ range(left(xyrect), right(xyrect); length=n+1)
    y = f(x)
    if y isa AbstractVector
      for k ∈ 1:length(y)
        ij = xy2ij(Point2f0(x, y[k]), c)
        p = set!(buf, p, ij)
      end
    else
      ij = xy2ij(Point2f0(x, y), c)
      p = set!(buf, p, ij)
    end
  end
  truncate!(buf, p)
  length(buf) < 1 && push!(buf, Point2f0(Inf32, Inf32))  # needed because Makie does not like 0-point lines
end

function pointcrop!(buf::AbstractVector{Point2f0}, c::Canvas)
  data = c.datasrc.data
  data === missing && return
  xyrect = c.xyrect[]
  p = 0
  for k in 1:size(data,1)
    if left(xyrect) <= data[k,1] <= right(xyrect) && bottom(xyrect) <= data[k,2] <= top(xyrect)
      ij = xy2ij(Point2f0(data[k,1], data[k,2]), c)
      p = set!(buf, p, ij)
    end
  end
  truncate!(buf, p)
  length(buf) < 1 && push!(buf, Point2f0(Inf32, Inf32))  # needed because Makie does not like 0-point lines
end

function linecrop!(buf::AbstractVector{Point2f0}, c::Canvas)
  data = c.datasrc.data
  data === missing && return
  p = 0
  n = size(data,1)
  if n > 0
    xyrect = c.xyrect[]
    lastpt = data[1,:]
    for k in 1:n
      nextpt = data[k < n ? k + 1 : n, :]
      lastpt[1] < left(xyrect) && data[k,1] < left(xyrect) && nextpt[1] < left(xyrect) && continue
      lastpt[1] > right(xyrect) && data[k,1] > right(xyrect) && nextpt[1] > right(xyrect) && continue
      lastpt[2] < bottom(xyrect) && data[k,2] < bottom(xyrect) && nextpt[2] < bottom(xyrect) && continue
      lastpt[2] > top(xyrect) && data[k,2] > top(xyrect) && nextpt[2] > top(xyrect) && continue
      lastpt = data[k,:]
      ij = xy2ij(Point2f0(data[k,1], data[k,2]), c)
      p = set!(buf, p, ij)
    end
  end
  truncate!(buf, p)
  length(buf) < 1 && push!(buf, Point2f0(Inf32, Inf32))  # needed because Makie does not like 0-point lines
end

# TODO: check if average time of pooling is correct
function linepool!(buf::AbstractVector{Point2f0}, c::Canvas; pooling=orderedextrema)
  data = c.datasrc.data
  data === missing && return
  n = size(data,1)
  n < 2 && return linecrop!(buf, c)
  x0 = data[1,1]
  dx = data[2,1] - x0
  x2ndx(x) = round(Int, (x - x0) / dx)
  xyrect = c.xyrect[]
  n1 = x2ndx(left(xyrect))
  n1 < 1 && (n1 = 1)
  n2 = x2ndx(right(xyrect))
  n2 > n && (n2 = n)
  p = 0
  if n1 <= n && n2 >= 1
    m = width(c.ijrect)
    b = fld(n2 - n1, m)
    if b <= 1
      for k ∈ n1:n2
        ij = xy2ij(Point2f0(data[k,1], data[k,2]), c)
        p = set!(buf, p, ij)
      end
    else
      for k ∈ n1:b:n2
        k2 = k + b - 1
        k2 > n && (k2 = n)
        blk = @view data[k:k2,2]
        x = (data[k,1] + data[k2,1]) / 2
        for y ∈ pooling(blk)
          ij = xy2ij(Point2f0(x, y), c)
          p = set!(buf, p, ij)
        end
      end
    end
  end
  truncate!(buf, p)
  length(buf) < 1 && push!(buf, Point2f0(Inf32, Inf32))  # needed because Makie does not like 0-point lines
end

function heatmappool!(buf::AbstractMatrix, c::Canvas; pooling=mean)
  data = c.datasrc.data
  data === missing && return
  pwidth, qwidth = size(data)
  xylims = c.datasrc.xyrect
  xpp = width(xylims) / pwidth
  ypq = height(xylims) / qwidth
  iwidth, jwidth = size(buf)
  xyrect = c.xyrect[]
  xpi = width(xyrect) / iwidth
  ypj = height(xyrect) / jwidth
  ppi = xpi / xpp
  qpj = ypj / ypq
  for j ∈ 1:jwidth
    for i ∈ 1:iwidth
      p = ((i - 1) * xpi + left(xyrect) - left(xylims)) / xpp
      q = ((j - 1) * ypj + bottom(xyrect) - bottom(xylims)) / ypq
      p1 = round(Int, p - ppi/2) + 1
      p2 = round(Int, p + ppi/2) + 1
      q1 = round(Int, q - qpj/2) + 1
      q2 = round(Int, q + qpj/2) + 1
      if p1 > pwidth || p2 < 1 || q1 > qwidth || q2 < 1
        buf[i,j] = 0
      else
        p1 < 1 && (p1 = 1)
        p2 > pwidth && (p2 = pwidth)
        q1 < 1 && (q1 = 1)
        q2 > qwidth && (q2 = qwidth)
        blk = @view data[p1:p2,q1:q2]
        buf[i,j] = pooling(blk)
      end
    end
  end
end

# TODO: add support for spreading?
function aggregate!(buf::AbstractMatrix, c::Canvas)
  m, n = size(buf)
  buf .= zero(buf[1,1])
  data = c.datasrc.data
  for k ∈ 1:size(data,1)
    i, j = xy2ij(Int, data[k,1], data[k,2], c)
    0 <= i < m && 0 <= j < n && (buf[i+1, j+1] += 1)
  end
end

### helpers

function set!(buf, p, pt)
  if p < length(buf)
    p += 1
    buf[p] = pt
  else
    push!(buf, pt)
    p = length(buf)
  end
  p
end

function truncate!(buf, p)
  p < length(buf) && splice!(buf, p+1:length(buf))
end

orderedextrema(x) = argmin(x) < argmax(x) ? (minimum(x), maximum(x)) : (maximum(x), minimum(x))


text:
####
#### Reference state
####

"""
    compute_ref_state!(
        state,
        grid::Grid,
        param_set::PS;
        Pg::FT,
        Tg::FT,
        qtg::FT,
    ) where {PS, FT}

TODO: add better docs once the API converges

The reference profiles, given
 - `grid` the grid
 - `param_set` the parameter set
"""
function compute_ref_state!(state, grid::Grid, param_set::PS; Pg::FT, Tg::FT, qtg::FT) where {PS, FT}

    p0_c = center_ref_state(state).p0
    ρ0_c = center_ref_state(state).ρ0
    α0_c = center_ref_state(state).α0
    p0_f = face_ref_state(state).p0
    ρ0_f = face_ref_state(state).ρ0
    α0_f = face_ref_state(state).α0

    q_pt_g = TD.PhasePartition(qtg)
    ts_g = TD.PhaseEquil_pTq(param_set, Pg, Tg, qtg)
    θ_liq_ice_g = TD.liquid_ice_pottemp(ts_g)

    # We are integrating the log pressure so need to take the log of the
    # surface pressure
    logp = log(Pg)

    # Form a right hand side for integrating the hydrostatic equation to
    # determine the reference pressure
    function rhs(logp, u, z)
        p_ = exp(logp)
        ts = thermo_state_pθq(param_set, p_, θ_liq_ice_g, qtg)
        R_m = TD.gas_constant_air(ts)
        T = TD.air_temperature(ts)
        return -FT(CPP.grav(param_set)) / (T * R_m)
    end

    # Perform the integration
    z_span = (grid.zmin, grid.zmax)
    @show z_span
    prob = ODE.ODEProblem(rhs, logp, z_span)
    sol = ODE.solve(prob, ODE.Tsit5(), reltol = 1e-12, abstol = 1e-12)

    parent(p0_f) .= sol.(vec(grid.zf))
    parent(p0_c) .= sol.(vec(grid.zc))

    p0_f .= exp.(p0_f)
    p0_c .= exp.(p0_c)

    # Compute reference state thermodynamic profiles
    @inbounds for k in real_center_indices(grid)
        ts = thermo_state_pθq(param_set, p0_c[k], θ_liq_ice_g, qtg)
        α0_c[k] = TD.specific_volume(ts)
    end

    @inbounds for k in real_face_indices(grid)
        ts = thermo_state_pθq(param_set, p0_f[k], θ_liq_ice_g, qtg)
        α0_f[k] = TD.specific_volume(ts)
    end

    ρ0_f .= 1 ./ α0_f
    ρ0_c .= 1 ./ α0_c
end


text:
using FGES
using Documenter

DocMeta.setdocmeta!(FGES, :DocTestSetup, :(using FGES); recursive=true)

makedocs(;
    modules=[FGES],
    authors="RobertGregg <robwgregg@gmail.com> and contributors",
    repo="https://github.com/RobertGregg/FGES.jl/blob/{commit}{path}#{line}",
    sitename="FGES.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)


text:
struct Simulation

    u::Vector{Vector{Float64}}
    t::Vector{Float64}

end


text:
export FixedPolynomialBasis, ScaledMonomialBasis, MonomialBasis

"""
    abstract type AbstractPolynomialBasis end

Polynomial basis of a subspace of the polynomials [Section~3.1.5, BPT12].

[BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R.
*Semidefinite Optimization and Convex Algebraic Geometry*.
Society for Industrial and Applied Mathematics, **2012**.
"""
abstract type AbstractPolynomialBasis end

"""
    struct FixedPolynomialBasis{PT<:MultivariatePolynomials.AbstractPolynomialLike, PV<:AbstractVector{PT}} <: AbstractPolynomialBasis
        polynomials::PV
    end

Polynomial basis with the polynomials of the vector `polynomials`.
For instance, `FixedPolynomialBasis([1, x, 2x^2-1, 4x^3-3x])` is the Chebyshev
polynomial basis for cubic polynomials in the variable `x`.
"""
struct FixedPolynomialBasis{PT<:MultivariatePolynomials.AbstractPolynomialLike, PV<:AbstractVector{PT}} <: AbstractPolynomialBasis
    polynomials::PV
end

function MultivariatePolynomials.polynomialtype(mb::FixedPolynomialBasis{PT}, T::Type) where PT
    C = MultivariatePolynomials.coefficienttype(PT)
    U = typeof(zero(C) * zero(T) + zero(C) * zero(T))
    MultivariatePolynomials.polynomialtype(PT, U)
end
function MultivariatePolynomials.polynomial(f::Function, fpb::FixedPolynomialBasis)
    sum(ip -> f(ip[1]) * ip[2], enumerate(fpb.polynomials))
end

"""
    struct MonomialBasis{MT<:MultivariatePolynomials.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPolynomialBasis
        monomials::MV
    end

Monomial basis with the monomials of the vector `monomials`.
For instance, `MonomialBasis([1, x, y, x^2, x*y, y^2])` is the monomial basis
for the subspace of quadratic polynomials in the variables `x`, `y`.
"""
struct MonomialBasis{MT<:MultivariatePolynomials.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPolynomialBasis
    monomials::MV
end
MonomialBasis(monomials) = MonomialBasis(monovec(monomials))

MultivariatePolynomials.polynomialtype(mb::MonomialBasis{MT}, T::Type) where MT = MultivariatePolynomials.polynomialtype(MT, T)
MultivariatePolynomials.polynomial(f::Function, mb::MonomialBasis) = polynomial(f, mb.monomials)
function MultivariatePolynomials.coefficients(p, ::Type{<:MonomialBasis})
    return coefficients(p)
end

"""
    struct ScaledMonomialBasis{MT<:MultivariatePolynomials.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPolynomialBasis
        monomials::MV
    end

*Scaled monomial basis* (see [Section 3.1.5, BPT12]) with the monomials of the vector `monomials`.
Given a monomial ``x^\\alpha = x_1^{\\alpha_1} \\cdots x_n^{\\alpha_n}`` of degree ``d = \\sum_{i=1}^n \\alpha_i``,
the corresponding polynomial of the basis is
```math
{d \\choose \\alpha}^{\\frac{1}{2}} x^{\\alpha} \\quad \\text{ where } \\quad
{d \\choose \\alpha} = \\frac{d!}{\\alpha_1! \\alpha_2! \\cdots \\alpha_n!}.
```

For instance, create a polynomial with the basis ``[xy^2, xy]`` creates the polynomial
``\\sqrt{3} a xy^2 + \\sqrt{2} b xy`` where `a` and `b` are new JuMP decision variables.
Constraining the polynomial ``axy^2 + bxy`` to be zero with the scaled monomial basis constrains
`a/√3` and `b/√2` to be zero.

[BPT12] Blekherman, G.; Parrilo, P. A. & Thomas, R. R.
*Semidefinite Optimization and Convex Algebraic Geometry*.
Society for Industrial and Applied Mathematics, **2012**.
"""
struct ScaledMonomialBasis{MT<:MultivariatePolynomials.AbstractMonomial, MV<:AbstractVector{MT}} <: AbstractPolynomialBasis
    monomials::MV
end
ScaledMonomialBasis(monomials) = ScaledMonomialBasis(monovec(monomials))

MultivariatePolynomials.polynomialtype(mb::ScaledMonomialBasis{MT}, T::Type) where MT = MultivariatePolynomials.polynomialtype(MT, promote_type(T, Float64))
scaling(m::MultivariatePolynomials.AbstractMonomial) = √(factorial(degree(m)) / prod(factorial, exponents(m)))
MultivariatePolynomials.polynomial(f::Function, mb::ScaledMonomialBasis) = polynomial(i -> scaling(mb.monomials[i]) * f(i), mb.monomials)
unscale_coef(t::MultivariatePolynomials.AbstractTerm) = coefficient(t) / scaling(monomial(t))
function MultivariatePolynomials.coefficients(p, ::Type{<:ScaledMonomialBasis})
    return unscale_coef.(terms(p))
end


text:
using Test
include("../src/classification.jl")
using Main.Classification

@testset "One hot encoding" begin
    @test onehot(0, [0, 1]) == Int8[1 0]
    @test onehot(7, 5:8) == Int8[0 0 1 0]
    @test onehot(:a, [:c, :g, :a, :5]) == Int8[0 0 1 0]
    @test onehot([2 3 4 2 2], 1:4) == Int8[0 0 0 0 0; 1 0 0 1 1; 0 1 0 0 0; 0 0 1 0 0]
    @test onehot([2 3 4 2 1]) == Int8[0 0 0 0 1; 1 0 0 1 0; 0 1 0 0 0; 0 0 1 0 0]
end

@testset "One cold decoding" begin
    @test onecold(onehot([2 3 4 2 1]), 1:4) == [2 3 4 2 1]
end


text:
using Kalman, GaussianDistributions, LinearAlgebra
using GaussianDistributions: ⊕ # independent sum of Gaussian r.v.
using Statistics
using GLMakie
using StaticArrays

function onestep(p, v)
    v += @SVector randn(2)
    w = 0.1 * @SVector randn(2)
    p += v + w
    p, v
end

x0 = zero(SVector{4, Float64})  #start point
P = SMatrix{4,4,Float64}(I)
# H: observation matrix
H = SMatrix{2,4,Float64}(I) # as only location is observed, not speed
# F: state-trasition
F = SMatrix{4,4,Float64}([
    1 0 1 0
    0 1 0 1
    0 0 1 0
    0 0 0 1
   ])
# Q: the covariance of the process noise
Q = SMatrix{4,4}(Diagonal([0.01, 0.01, 0.1, 0.1]))
# R: the covariance of the observation noise
R = SMatrix{2,2,Float64}(I)

p = @SVector zeros(2) # initial poistion
v = @SVector ones(2) # initial velocity
n = 100

pp = Gaussian(x0, P)
ps = [pp] # vector of filtered Gaussians
obs = Point2f0[]
for i in 1:n
    global pp, p, v
    # predict
    pp = F*pp ⊕ Gaussian(zero(x0), Q) #same as Gaussian(Φ*p.μ, Φ*p.Σ*Φ' + Q)
    # observe
    p, v = onestep(p, v)
    xobs = p + @SVector randn(2)
    push!(obs, xobs)
    # correct
    pp, yres, _ = Kalman.correct(Kalman.JosephForm(), pp, (Gaussian(xobs, R), H))
    push!(ps, pp) # save filtered density
end

path = [Point2f0(mean(p)[1:2]) for p in ps]
vel = [Point2f0(mean(p)[3:4]) for p in ps]

fig = Figure()
ax = Axis(fig[1, 1], aspect = DataAspect())
scatter!(obs)
lines!(obs)
sl = Slider(fig[2, 1], range = 1:n, startvalue = 1)
p = @lift([Point2f0(path[$(sl.value)])])
d = @lift([Point2f0(vel[$(sl.value)])])
arrows!(p, d, arrowcolor = :red, linecolor = :red)#, linewidth = 3, )
w = 0.01only(diff([extrema(vcat(obs...))...]))
on(p) do pp
    pp = p[]
    ps = [pp; pp .+ d[]]
    xlim = extrema(first, ps) .+ (-w, w)
    ylim = extrema(last, ps) .+ (-w, w)
    limits!(ax, xlim..., ylim...)
end
fig


text:
# Train an normalizing flow for seismic imaging
# Authors: Ali Siahkoohi, alisk@gatech.edu
# Date: February 2021
"""
Example:
julia train_hint_imaging.jl
"""

using DrWatson
@quickactivate :ReliabilityAwareImaging

using InvertibleNetworks
using HDF5
using JLD2
using Random
using Statistics
using ArgParse
using ProgressMeter
using Flux: gpu, cpu
using Flux.Data: DataLoader
using Flux.Optimise: Optimiser, ExpDecay, update!, ADAM

# Random seed
Random.seed!(19)

s = ArgParseSettings()
@add_arg_table! s begin
    "--max_epoch"
        help = "Maximum # of epochs"
        arg_type = Int
        default = 100
    "--lr"
        help = "starting learning rate"
        arg_type = Float32
        default = 1f-4
    "--lr_step"
        help = "lr scheduler applied at evey lr_step epoch"
        arg_type = Int
        default = 2
    "--batchsize"
        help = "batch size"
        arg_type = Int
        default = 16
    "--n_hidden"
        help = "# hidden channels"
        arg_type = Int
        default = 128
    "--depth"
        help = "depth of the network"
        arg_type = Int
        default = 12
    "--sim_name"
        help = "simulation name"
        arg_type = String
        default = "seismic-imaging-posterior-learning"
end
parsed_args = parse_args(s)

max_epoch = parsed_args["max_epoch"]
lr = parsed_args["lr"]
lr_step = parsed_args["lr_step"]
batchsize = parsed_args["batchsize"]
n_hidden = parsed_args["n_hidden"]
depth = parsed_args["depth"]
sim_name = parsed_args["sim_name"]

# Define raw data directory
mkpath(datadir("training-data"))
data_path = datadir("training-data", "training-pairs.h5")

# Download the dataset into the data directory if it does not exist
if isfile(data_path) == false
    run(`wget https://www.dropbox.com/s/53u8ckb9aje8xv4/'
        'training-pairs.h5 -q -O $data_path`)
end

# Load seismic images and create training and testing data
file = h5open(data_path, "r")
X_train = file["dm"][:, :, :, :]
Y_train = file["rtm"][:, :, :, :]

Y_train[1:10, :, :, :] .= 0f0

nx, ny, nc, nsamples = size(X_train)

AN_x = ActNorm(nsamples)
AN_y = ActNorm(nsamples)
X_train = AN_x.forward(X_train)
Y_train = AN_y.forward(Y_train)

AN_params_x = get_params(AN_x)
AN_params_y = get_params(AN_y)

# Split in training/testing
ntrain = Int(floor((nsamples*.9)))
train_idx = randperm(nsamples)[1:ntrain]

# Dimensions after wavelet squeeze to increase no. of channels
nx = Int(nx/2)
ny = Int(ny/2)
n_in = Int(nc*4)

# Create network
CH = NetworkConditionalHINT(
    n_in, n_hidden, depth, k1=3, k2=3, p1=1, p2=1
) |> gpu

X = wavelet_squeeze(X_train[:, :, :, 1:64]) |> gpu
Y = wavelet_squeeze(Y_train[:, :, :, 1:64]) |> gpu
# CH.forward(X, Y)

# Training
# Batch extractor
train_loader = DataLoader(train_idx, batchsize=batchsize, shuffle=true)
num_batches = length(train_loader)

# Optimizer
opt = Optimiser(ExpDecay(lr, 9f-1, num_batches*lr_step, 1f-6), ADAM(lr))

# Training log keeper
fval = zeros(Float32, num_batches*max_epoch)
p = Progress(num_batches*max_epoch)

for epoch=1:max_epoch
    for (itr, idx) in enumerate(train_loader)
        Base.flush(Base.stdout)

        # Apply wavelet squeeze (change dimensions to -> n1/2 x n2/2 x nc*4 x ntrain)
        X = wavelet_squeeze(X_train[:, :, :, idx])
        Y = wavelet_squeeze(Y_train[:, :, :, idx])

        X = X |> gpu
        Y = Y |> gpu

        fval[(epoch-1)*num_batches + itr] = loss_supervised(CH, X, Y)[1]

        ProgressMeter.next!(
            p;
            showvalues = [
                (:Epoch, epoch),
                (:Itreration, itr),
                (:NLL, fval[(epoch-1)*num_batches + itr])
            ]
        )

        # Update params
        for p in get_params(CH)
            update!(opt, p.data, p.grad)
        end
        clear_grad!(CH)
    end

    # Saving parameters and logs
    Params = get_params(CH) |> cpu
    save_dict = @strdict epoch max_epoch lr lr_step batchsize n_hidden depth sim_name Params fval train_idx AN_params_x AN_params_y
    @tagsave(
        datadir(sim_name, savename(save_dict, "jld2"; digits=6)),
        save_dict;
        safe=true
    )

end


text:
# these are the default solver packages in order of decreasing precedence
# for each supported problem class.

const LPsolvers = [(:Clp,:ClpSolver),
                   (:GLPKMathProgInterface,:GLPKSolverLP),
                   (:Gurobi,:GurobiSolver),
                   (:CPLEX,:CplexSolver),
                   (:Mosek,:MosekSolver),
                   (:Xpress,:XpressSolver)]

const MIPsolvers = [(:Cbc,:CbcSolver),
                    (:GLPKMathProgInterface,:GLPKSolverMIP),
                    (:Gurobi,:GurobiSolver),
                    (:CPLEX,:CplexSolver),
                    (:Mosek,:MosekSolver),
                    (:Xpress,:XpressSolver)]

const QPsolvers = [(:Gurobi,:GurobiSolver),
                   (:CPLEX,:CplexSolver),
                   (:Mosek,:MosekSolver),
                   (:Ipopt,:IpoptSolver),
                   (:Xpress,:XpressSolver)]

const SDPsolvers = [(:Mosek,:MosekSolver),
                    (:SCS,:SCSSolver)]

const NLPsolvers = [(:Ipopt,:IpoptSolver),
                    (:KNITRO,:KnitroSolver),
                    (:Mosek,:MosekSolver)]

const Conicsolvers = [(:ECOS,:ECOSSolver),
                      (:SCS,:SCSSolver),
                      (:Mosek,:MosekSolver)]

using Base.Meta

# Don't load packages for default solvers until needed.
# This reduces the startup time for MathProgBase.
for solvertype in ["LP", "MIP", "QP", "SDP", "NLP", "Conic"]
    typename = Symbol("Default"*solvertype*"Solver")
    @eval begin
        type $typename <: SolverInterface.AbstractMathProgSolver
        end
    end
    defaultname = Symbol("default"*solvertype*"solver")
    @eval $defaultname = ($typename)()

    solvers = Symbol(solvertype*"solvers")

    for t in (:LinearQuadraticModel,:ConicModel,:NonlinearModel)
        @eval function SolverInterface.$t(s::$typename)
            for (pkgname, solvername) in $solvers
                try
                    eval(Expr(:import,pkgname))
                catch
                    if isdir(Pkg.dir((string(pkgname))))
                        warn("Package ",string(pkgname),
                             " is installed but couldn't be loaded. ",
                             "You may need to run `Pkg.build(\"$pkgname\")`")
                    end
                    continue
                end
                ex = Expr(:(=), $(quot(defaultname)),
                          Expr(:call, Expr(:., pkgname, quot(solvername))))
                eval(ex)
                ex = Expr(:call, $(quot(t)), $(quot(defaultname)))
                return eval(ex)
            end
            suggestions = join(["\"$(pkgname)\", " for (pkgname,solvername) in $solvers], ' ')
            error("No ",$solvertype, " solver detected. Try installing one of the following packages: ", suggestions, " and restarting Julia")
        end
    end
end


text:
#=
test_timeevolution:
- Julia version: 1.5.2
- Author: qcancam
- Date: 2021-01-03
=#

include("../src/MatrixProductStates.jl")

using .MatrixProductStates, SparseArrays, Arpack, PyCall, PyPlot

#np=pyimport("numpy")
#@pyimport matplotlib.pyplot as plt

g = 1.0; L = 3
T=100.0
N=1000.0
τ = T/N

function H_LZ(g, L)
    H = Array{Complex{Float64}, 2}[]
    id = [1  0;
          0  1]
    σˣ = [0  1;
          1  0]
    σᶻ = [1  0;
          0 -1]
    for i in 1:L-2
        push!(H,g*σᶻ⊗id+σˣ⊗id)
    end
    #=L-1=#
    push!(H,g*σᶻ⊗id+σˣ⊗id+g*id⊗σᶻ+id⊗σˣ)

    return H
end

function H_Z_MPO(L)
    id = [1  0;
          0  1]
    σˣ = [0  1;
          1  0]
    σᶻ = [1  0;
          0 -1]
    W_tnsr = zeros(Complex{Float64}, 2, 2, 2, 2)
    W_tnsr[1, 1, :, :] = id
    W_tnsr[2, 1, :, :] = σᶻ
    W_tnsr[2, 2, :, :] = id

    return MPO(W_tnsr, L)
end

#H = H_LZ(g, L)
ψ = randn(MPS{L, Complex{Float64}}, 100, 2)

glist=[-T/2+T/N*i for i in 1:N]
H1=[H_LZ(g, L)[1] for g in glist]
H2=[H_LZ(g, L)[2] for g in glist]
HLm1=[H_LZ(g, L)[L-1] for g in glist]

ϕ,time_list,expect_list= real_timedep_time_evolution(ψ, H1,H2,HLm1,H_Z_MPO(L), T,N,50)

#ϕ= real_time_evolution(ψ, H[1],H[2],H[L-1],H_TFIM_MPO(g,L), T,N,50)

fig = plt.figure()
ax = fig.add_subplot(1, 1, 1)
ax.plot(time_list, expect_list, "-", c="b", ms=10)
ax.set_xlabel("t")
ax.set_ylabel("⟨-∑σᶻ⟩")
#ax.set_title("2D Line Plot Sample")
#ax.grid(true)
plt.show()


text:
"""
    ecdf(y::Vector{<:Real})::Tuple{Vector{<:Real}, Vector{<:Real}}

Compute the empirical cumulative distribution function using the Gumbel formula.

The empirical quantiles are computed using the Gumbel plotting positions as
as recommended by [Makkonen (2006)](https://journals.ametsoc.org/jamc/article/45/2/334/12668/Plotting-Positions-in-Extreme-Value-Analysis).

# Example
```julia-repl
julia> (x, F̂) = Extremes.ecdf(y)
```

# Reference
Makkonen, L. (2006). Plotting positions in extreme value analysis. Journal of
Applied Meteorology and Climatology, 45(2), 334-340.
"""
function ecdf(y::Vector{<:Real})::Tuple{Vector{<:Real}, Vector{<:Real}}
    ys = sort(y)
    n = length(ys)
    p = collect(1:n)/(n+1)

    return ys, p
end

"""
    checkstationarity(model::EVA)

Check if the extreme value model `model` is stationary.
"""
function checkstationarity(m::EVA)

    if getcovariatenumber(m) > 0
        @info "The graph is optimized for stationary models and the model provided is not."
    end

end

"""
    checknonstationarity(model::EVA)

Check if the extreme value model `model` is nonstationary.
"""
function checknonstationarity(m::EVA)

    if getcovariatenumber(m) == 0
        @info "The graph is optimized for non-stationary models and the model provided is not."
    end

end

include(joinpath("validationplots", "plots.jl"))
include(joinpath("validationplots", "plots_data.jl"))
include(joinpath("validationplots", "plots_std_data.jl"))


text:
using Knet: logsoftmax, softmax

function kl_div(input_logits, target_logits)
    logsoft_input = logsoftmax(input_logits; dims=1)
    soft_target = softmax(target_logits; dims=1)
    logsoft_target = logsoftmax(target_logits; dims=1)
    return mean(soft_target .* (logsoft_target .- logsoft_input))
end

function mse_loss(y_hat, y)
    return mean((y - y_hat)^2)
end

function dice_loss(score, target)
    smooth = 1e-5
    intersect = sum(score .* target)
    y_sum = sum(target .* target)
    z_sum = sum(score .* score)
    loss = (2 .* intersect .+ smooth) / (z_sum + y_sum .+ smooth)
    loss = 1 .- loss
    return loss
end

max_stars_repo_name: skoc/julia-project

max_issues_repo_name: skoc/julia-project

max_forks_repo_name: skoc/julia-project


text:
"""
    drawTagBox!(image, tag)

Draw a box around the tag.
`imageCol = RGB.(image)
foreach(tag->drawTagBox!(imageCol, tag), tags)
`
"""
function drawTagBox!(image::Array{T,2}, tag::AprilTag; width::Int = 1, drawReticle::Bool = true) where T <: RGB
    cpoints = map(tag -> Point(round.(Int,tag)...),tag.p)
    # Reticle
    if drawReticle
        midx = (cpoints[3].x + cpoints[1].x) / 2
        midy = (cpoints[3].y + cpoints[1].y) / 2
        x1 = Point(round(midx - (cpoints[2].x - cpoints[1].x) / 4), round(midy - (cpoints[2].y - cpoints[1].y) / 4))
        x2 = Point(round(midx + (cpoints[2].x - cpoints[1].x) / 4), round(midy + (cpoints[2].y - cpoints[1].y) / 4))
        y1 = Point(round(midx - (cpoints[4].x - cpoints[1].x) / 4), round(midy - (cpoints[4].y - cpoints[1].y) / 4))
        y2 = Point(round(midx + (cpoints[4].x - cpoints[1].x) / 4), round(midy + (cpoints[4].y - cpoints[1].y) / 4))
        draw!(image, LineSegment( x1, x2), RGB{N0f8}(1.0, 0.0, 0.0), boundedBresenham)
        draw!(image, LineSegment( y1, y2), RGB{N0f8}(1.0, 0.0, 0.0), boundedBresenham)
    end
    # Make a box around the tag
    if width <= 1
        draw!(image, LineSegment( cpoints[1], cpoints[2]), RGB{N0f8}(0.0, 0.0, 1.0), boundedBresenham)
        draw!(image, LineSegment( cpoints[2], cpoints[3]), RGB{N0f8}(0.0, 0.0, 1.0), boundedBresenham)
        draw!(image, LineSegment( cpoints[3], cpoints[4]), RGB{N0f8}(1.0, 0.0, 0.0), boundedBresenham)
        draw!(image, LineSegment( cpoints[4], cpoints[1]), RGB{N0f8}(0.0, 1.0, 0.0), boundedBresenham)
    else
        drawThickLine!(image, cpoints[1], cpoints[2], RGB{N0f8}(0.0, 0.0, 1.0), width)
        drawThickLine!(image, cpoints[2], cpoints[3], RGB{N0f8}(0.0, 0.0, 1.0), width)
        drawThickLine!(image, cpoints[3], cpoints[4], RGB{N0f8}(1.0, 0.0, 0.0), width)
        drawThickLine!(image, cpoints[4], cpoints[1], RGB{N0f8}(0.0, 1.0, 0.0), width)
    end

    return nothing
end

function boundedBresenham(img::AbstractArray{T, 2}, y0::Int, x0::Int, y1::Int, x1::Int, color::T) where T<:Colorant
    dx = abs(x1 - x0)
    dy = abs(y1 - y0)

    sx = x0 < x1 ? 1 : -1
    sy = y0 < y1 ? 1 : -1;

    err = (dx > dy ? dx : -dy) / 2

    (row,col) = size(img)
    while true
        if 1 <= x0 <= col && 1 <= y0 <= row
            img[y0, x0] = color
        end
        (x0 != x1 || y0 != y1) || break
        e2 = err
        if e2 > -dx
            err -= dy
            x0 += sx
        end
        if e2 < dy
            err += dx
            y0 += sy
        end
    end

    img
end

function drawThickLine!(image, startpoint, endpoint, colour, thickness)
    (row, col) = size(image)
    x1 = startpoint.x
    y1 = startpoint.y
    x2 = endpoint.x
    y2 = endpoint.y

    draw!(image, LineSegment(x1,y1,x2,y2), colour, boundedBresenham)

    if x1 != x2 && abs((y2-y1)/(x2-x1)) < 1
        for tn = 1:thickness
            i = tn ÷ 2
            iseven(tn) && draw!(image, LineSegment(x1,y1-i,x2,y2-i), colour, boundedBresenham)
            isodd(tn)  && draw!(image, LineSegment(x1,y1+i,x2,y2+i), colour, boundedBresenham)
        end
    else
        for tn = 1:thickness
            i = tn ÷ 2
            iseven(tn) && draw!(image, LineSegment(x1-i,y1,x2-i,y2), colour, boundedBresenham)
            isodd(tn)  && draw!(image, LineSegment(x1+i,y1,x2+i,y2), colour, boundedBresenham)
        end
    end
end

"""
    drawTagAxes!(image, tag, CameraMatrix)

Draw the tag x, y, and z axes to show the orientation.
`imageCol = RGB.(image)
foreach(tag->drawTagAxes!(imageCol, tag, K), tags)
`
"""
function drawTagAxes!(image::Array{T,2}, tag::AprilTag, K::Array{Float64,2}) where T <: RGB

    Kp = [K [0;0]; 0.0 0.0 1.0 0.0]
    pose = AprilTags.homography_to_pose(tag.H, K[1,1], K[2,2], K[1,3], K[2,3])

    # calculate and project
    p0 = Kp*pose[:,4]
    p0 /= p0[3]

    p1 = Kp*(pose[:,4] + pose[:,1])
    p1 /= p1[3]

    p2 = Kp*(pose[:,4] + pose[:,2])
    p2 /= p2[3]

    p3 = Kp*(pose[:,4] + pose[:,3])
    p3 /= p3[3]

    ip0 = Point(round.(Int,p0[1:2])...)
    ip1 = Point(round.(Int,p1[1:2])...)
    ip2 = Point(round.(Int,p2[1:2])...)
    ip3 = Point(round.(Int,p3[1:2])...)
    #TODO check that the line stays on the image to avoid bounds errors
    draw!(image, LineSegment(ip0, ip1), RGB{N0f8}(1.0, 0.0, 0.0), boundedBresenham)
    draw!(image, LineSegment(ip0, ip2), RGB{N0f8}(0.0, 1.0, 0.0), boundedBresenham)
    draw!(image, LineSegment(ip0, ip3), RGB{N0f8}(0.0, 0.0, 1.0), boundedBresenham)
    image
end


text:
function update_lambda!(state::State, data::Data, prior::Prior,
                        flags::Vector{Flag})
  update_lambdaC!(state, data, prior, flags)
  update_lambdaT!(state, data, prior, flags)
end

function sample_lambda(y::AbstractVector{<:Real}, mu::AbstractVector{<:Real},
                       psi::AbstractVector{<:Real},
                       omega::AbstractVector{<:Real},
                       nu::AbstractVector{<:Real},
                       eta::AbstractVector{<:Real},
                       zeta::AbstractVector{<:Real},
                       v::AbstractVector{<:Real}, flags::Vector{Flag})
  logeta = log.(eta)

  if UpdateLambdaWithSkewT() in flags
    st = [Util.AltSkewT(loc=mu[k], df=nu[k], altscale=sqrt(omega[k]),
                        altskew=psi[k]) for k in eachindex(eta)]
    lambda = [let
                logmix = logpdf.(st, y[n]) + logeta
                MCMC.wsample_logprob(logmix)
              end for n in eachindex(y)]
  else  # update with latent var representation.
    lambda = [let
                loc = mu + psi * zeta[n]
                scale = sqrt.(omega) / sqrt(v[n])
                logmix = normlogpdf.(loc, scale, y[n]) +
                         gammalogpdf.(nu/2, 2 ./ nu, v[n]) + logeta
                MCMC.wsample_logprob(logmix)
              end for n in eachindex(y)]
  end
  return lambda
end

function update_lambdaC!(state::State, data::Data, prior::Prior,
                         flags::Vector{Flag})
  state.lambdaC .= sample_lambda(data.yC_finite,
                                 state.mu, state.psi, state.omega, state.nu,
                                 state.etaC, state.zetaC, state.vC, flags)
end

function update_lambdaT!(state::State, data::Data, prior::Prior,
                         flags::Vector{Flag})
  etaT_star = state.beta ? state.etaT : state.etaC
  state.lambdaT .= sample_lambda(data.yT_finite,
                                 state.mu, state.psi, state.omega, state.nu,
                                 etaT_star, state.zetaT, state.vT, flags)
end

max_stars_repo_name: luiarthur/CytofDensityEstimation

max_issues_repo_name: luiarthur/CytofDensityEstimation

max_forks_repo_name: luiarthur/CytofDensityEstimation


text:
abstract type AbstractTree{T} end

### Interface (TODO finish later)

# Return height of the tree.
# height(t::AbstractTree) = nothing

# Print a tree
# Base.show(io::IO, t::AbstractTree) = nothing


text:
function make_vandermonte_real(c::Vector{T}) where {T}
    n = length(c)
    V = zeros(T, n, n)
    for k = 1:n
        for kk = 1:n
            V[k, kk] = c[kk] ^ (k-1)
        end
    end
    return V
end

function find_real_eigenvalue(A::Matrix{BigFloat})
    n, _ = size(A)
    isodd(n) || error("A needs to be odd sized")
    ev = eigvals(A)
    min_val, min_ind = findmin(abs.(imag.(ev)))
    (min_val != 0.0) && error("something is wrong")
    return real.(ev[min_ind])
end

function calc_b_hat(A, c)
    gamma_0 = find_real_eigenvalue(A)
    n, _ = size(A)
    rhs = BigFloat.(collect(1 .// (1:n)))
    rhs[1] = 1 - gamma_0
    V = make_vandermonte_real(c)
    return vcat([gamma_0], V \ rhs)
end

function radau_butcher_table_plus(s::Int64)
    # ABC returns coefficients for S-stage Gauss-Legendre methods. (From RKGL, DiffMan package v.2 by K. Engo, A. Marthinsen & H. Munthe-Kaas.)
    P_Radau = output_p_Radau(s)
    c = roots(BigFloat.(coeffs(P_Radau)))
    c = sort(real.(c))
    V = make_vandermonte_real(c)'
    J = diagm(0 => 1 .// (1:s) )
    A = diagm(0 => c) * V * J / V
    b = (ones(BigFloat, 1, s) * J / V)[:]
    bi = get_dense_radau(J, V, c)
    if isodd(s)
        b̂ = calc_b_hat(A, c)
    else
        b̂ = zeros(BigFloat, s + 1) .+ NaN
    end
    return A, b, c, bi, b̂
end

function make_legendre_poly_by_order(n::Int64)
    P_x2_minus_1 = Poly([-1, 0, 1])
    return (1 // (2^n * factorial(n))) * polyder(P_x2_minus_1^n, n)
end

function output_p_Radau(n::Int64)
    Pn = make_legendre_poly_by_order(n)
    Pn_m = make_legendre_poly_by_order(n-1)
    P_2x_minus_1 = Poly([-1, 2])
    P_Radau = polyval(Pn, P_2x_minus_1) - polyval(Pn_m, P_2x_minus_1)
    return P_Radau
end

function get_dense_radau(J, V, c)
    s = length(c)
    interp_poly = (J / V)' .* ((1 ./ c) * collect(1:s)')
    return interp_poly
end


text:
using JSON

export Scope,
       withscope,
       Observable,
       @private,
       setobservable!,
       on, onjs,
       evaljs,
       onmount,
       onimport,
       ondependencies,
       adddeps!,
       import!,
       addconnection!

import Compat.Sockets: send
import Observables: Observable, AbstractObservable, listeners

# bool marks if it is synced
const ObsDict = Dict{String, Tuple{AbstractObservable, Union{Nothing,Bool}}}

mutable struct Scope
    dom::Any
    observs::ObsDict
    private_obs::Set{String}
    systemjs_options::Any
    # This should be an array of scopes, but there currently exist some circular
    # relationships between imports and scopes.
    imports::Vector{Asset}
    # A collection of handler functions associated with various observables in
    # this scope. Of the form
    # "observable-name" => ["array", "of", "JS", "strings"]
    # where each JS-string is a function that is invoked when the observable
    # changes.
    jshandlers::Any
    pool::ConnectionPool

    mount_callbacks::Vector{JSString}

    function Scope(
            dom::Any,
            observs::ObsDict, # bool marks if it is synced
            private_obs::Set{String},
            systemjs_options::Any,
            imports::Vector{Asset},
            jshandlers::Any,
            pool::ConnectionPool,
            mount_callbacks::Vector{JSString}
        )
        scope = new(
            dom, observs, private_obs,
            systemjs_options, imports, jshandlers, pool,
            mount_callbacks
        )
        register_scope!(scope)
        return scope
    end
end

function Base.getproperty(scope::Scope, property::Symbol)
    if property === :id
        Base.depwarn(
            "Accessing `scope.id` is deprecated, use `scopeid(scope)` instead.",
            :webio_scope_id,
        )
        return scopeid(scope)
    end
    return getfield(scope, property)
end

"""
    scopeid(scope)

Get a unique identifier for the specified scope.
This is the identifier that is used to reference scopes when communicating
between Julia and the Browser.

This defers to `Base.objectid` under the hood, though is cast to a string for
mostly historical reasons.
"""
scopeid(scope::Scope) = string(objectid(scope))

const global_scope_registry = Dict{String, Scope}()

function register_scope!(scope::Scope)
    global_scope_registry[scopeid(scope)] = scope
end

function deregister_scope!(scope::Scope)
    delete!(global_scope_registry, scopeid(scope))
end

function lookup_scope(uuid::String)
    get(global_scope_registry, uuid) do
        error("Could not find scope $(uuid)")
    end
end

"""
    Scope(<keyword arguments>)

An object which can send and receive messages.

# Arguments
- `dom`: The DOM node that will be rendered when the scope is displayed in the
    browser.
- `imports`: An collection of `Asset`s to load (see [`Asset`](@ref) for
    more documentation) when the scope is mounted. If the entry is not an
    `Asset` then it should be an argument to construct an `Asset`.

# Examples
```julia
myscope = Scope(
    dom=node(:p, "I'm a little scope!"),
    imports=[Asset("foo.js"), "bar.css", "spam" => "spam.js"],
)
```
"""
function Scope(;
        dom = dom"span"(),
        outbox::Union{Channel, Nothing} = nothing,
        observs::Dict = ObsDict(),
        private_obs::Set{String} = Set{String}(),
        systemjs_options = nothing,
        imports = Asset[],
        jshandlers::Dict = Dict(),
        mount_callbacks::Vector{JSString} = JSString[],
        # Deprecated
        id=nothing,
    )
    if id !== nothing
        Base.depwarn(
            "`Scope(; id, kwargs...)` is deprecated, "
                * "use `Scope(kwargs...)` instead.",
            :webio_scope_id_keyword_arg,
        )
    end
    imports = Asset[Asset(i) for i in imports]
    pool = outbox !== nothing ? ConnectionPool(outbox) : ConnectionPool()
    return Scope(
        dom, observs, private_obs, systemjs_options,
        imports, jshandlers, pool, mount_callbacks
    )
end

function Scope(id; kwargs...)
    Base.depwarn(
        "Scope(id; kwargs...) is deprecated, use Scope(kwargs...) instead.",
        :webio_scope_id_positional_arg,
    )
    return Scope(; kwargs...)
end

(w::Scope)(arg) = (w.dom = arg; w)
Base.wait(scope::Scope) = ensure_connection(scope.pool)

function Observables.on(f, w::Scope, key)
    key = string(key)
    listener, _ = Base.@get! w.observs key (Observable{Any}(w, key, nothing), nothing)
    on(f, listener)
end

private(scope::Scope, props...) = foreach(p->push!(scope.private_obs, string(p)), props)

macro private(ex)
    if ex.head == :(=)
        ref, val = ex.args
        scope, key = ref.args
        quote
            x=$(esc(ex))
            private($(esc(scope)), $(esc(key)))
            x
        end
    else
        error("Wrong usage of @private. Use as in `@private scope[\"key\"] = observable`")
    end
end

# we need to maintain mapping from an ID to observable
# in order to allow interpolation of observables.
const observ_id_dict = WeakKeyDict()

function setobservable!(ctx, key, obs; sync=nothing)
    key = string(key)
    if haskey(ctx.observs, key)
        @warn("An observable named $key already exists in scope $(scopeid(ctx)).
             Overwriting.")
    end

    ctx.observs[key] = (obs, sync)

    # the following metadata is stored for use in interpolation
    # of observables into DOM trees and `@js` expressions
    observ_id_dict[obs]  = (WeakRef(ctx), key)
    obs
end

# Ask JS to send stuff
function setup_comm(f, ob::AbstractObservable)
    if haskey(observ_id_dict, ob)
        scope, key = observ_id_dict[ob]
        # if !(key in scope.value.private_obs)
        #     evaljs(scope.value, js"""
        #            console.log(this)
        #            this.observables[$key].sync = true
        #     """)
        # end
    end
end

# TODO: hook `off` up

function Base.getindex(w::Scope, key)
    key = string(key)
    if haskey(w.observs, key)
        w.observs[key][1]
    else
        Observable{Any}(w, key, nothing)
    end
end

function Base.setindex!(w::Scope, obs, key)
    setobservable!(w, key, obs)
end

function Observable{T}(ctx::Scope, key, value; sync=nothing) where {T}
    setobservable!(ctx, key, Observable{T}(value), sync=sync)
end

function Observable(ctx::Scope, key, val::T; sync=nothing) where {T}
    Observable{T}(ctx, key, val; sync=sync)
end

function JSON.lower(scope::Scope)
    obs_dict = Dict()
    for (k, ob_) in scope.observs
        ob, sync = ob_
        if k in scope.private_obs
            continue
        end
        if sync === nothing
            # by default, we sync if there are any listeners
            # other than the JS back edge
            sync = any(f-> !isa(f, SyncCallback), listeners(ob))
        end
        obs_dict[k] = Dict(
            "sync" => sync,
            "value" => ob[],
            "id" => obsid(ob)
        )
    end
    Dict(
        "id" => scopeid(scope),
        "systemjs_options" => scope.systemjs_options,
        "imports" => lowerassets(scope.imports),
        "handlers" => scope.jshandlers,
        "mount_callbacks" => scope.mount_callbacks,
        "observables" => obs_dict
    )
end

function render(s::Scope)
    node(s, s.dom)
end

"""
    send_command(scope, command[, "key" => "value", ...])

Send a command message for a scope. A command is essentially a fire-and-forget
style message; no response or acknowledgement is expected.
"""
function send_command(scope::Scope, command, data::Pair...)
    message = Dict(
        "type" => "command",
        "command" => command,
        "scope" => scopeid(scope),
        data...
    )
    send(scope.pool, message)
    nothing
end

"""
Send a command to update the frontend's value for an observable.
"""
send_update_observable(scope::Scope, name::AbstractString, value) = send_command(
    scope,
    "update_observable",
    "name" => name,
    "value" => value,
)

function send_request(scope::Scope, request, data::Pair...)
    send_request(scope.pool, request, "scope" => scopeid(scope), data...)
end

macro evaljs(ctx, expr)
    @warn("@evaljs is deprecated, use evaljs function instead")
    :(send_request($(esc(ctx)), "eval", "expression" => $(esc(expr))))
end

function evaljs(ctx, expr)
    send_request(ctx, "eval", "expression" => expr)
end

function onmount(scope::Scope, f::JSString)
    push!(scope.mount_callbacks, f)
    return scope
end

function onimport(scope::Scope, f::JSString)
    onmount(scope, js"""
        function () {
            var handler = ($(f));
            ($(Async(scope.imports))).then((imports) => handler.apply(this, imports));
        }
        """)
end

onmount(scope::Scope, f) = onmount(scope, JSString(f))
onimport(scope::Scope, f) = onimport(scope, JSString(f))

Base.@deprecate ondependencies(ctx, jsf) onimport(ctx, jsf)

"""
A callable which updates the frontend
"""
struct SyncCallback
    ctx
    f
end

(s::SyncCallback)(xs...) = s.f(xs...)
"""
Set observable without synchronizing with the counterpart on the browser
"""
function set_nosync(ob, val)
    Observables.setexcludinghandlers(ob, val, x -> !(x isa SyncCallback))
end

const lifecycle_commands = ["scope_created"]

function dispatch(ctx, key, data)
    if haskey(ctx.observs, string(key))
        # this message has come from the browser
        # so don't update the browser back!
        set_nosync(ctx.observs[key][1], data)
    else
        if !(key in lifecycle_commands)
            @warn("$key does not have a handler for scope id $(scopeid(ctx))")
        end
    end
end

function onjs(ctx, key, f)
    push!(Base.@get!(ctx.jshandlers, key, []), f)
end

function offjs(ctx, key, f)
    if f in get(ctx.jshandlers, key, [])
        keys = ctx.jshandlers[key]
        deleteat!(keys, findall(in(f), keys))
    end
    nothing
end

function ensure_sync(ctx, key)
    ob = ctx.observs[key][1]
    # have at most one synchronizing handler per observable
    if !any(x->isa(x, SyncCallback) && x.ctx==ctx, listeners(ob))
        f = SyncCallback(ctx, (msg) -> send_update_observable(ctx, key, msg))
        on(SyncCallback(ctx, f), ob)
    end
end

function onjs(ob::AbstractObservable, f)
    if haskey(observ_id_dict, ob)
        ctx, key = observ_id_dict[ob]
        ctx = ctx.value
        # make sure updates are set up to propagate to JS
        ensure_sync(ctx, key)
        onjs(ctx, key, f)
    else
        error("This observable is not associated with any scope.")
    end
end

function Base.show(io::IO, ::WEBIO_NODE_MIME, x::Scope)
    write(io, JSON.json(render(x)))
end
function Base.show(io::IO, m::MIME"text/html", x::Scope)
    show(io, m, render(x))
end

function _show(io::IO, el::Scope, indent_level=0)
    showindent(io, indent_level)
    _show(io, el.dom, indent_level)
end

Base.@deprecate_binding Context Scope
Base.@deprecate handle! on
Base.@deprecate handlejs! onjs

# Adding assets to scopes
import!(scope::Scope, x) = push!(scope.imports, Asset(x))
Base.@deprecate adddeps!(scope, x) import!(scope, x)


text:
module IntoTheJungleCodeModSetColorGradeIfFlagOnSpawnController

using ..Ahorn, Maple

@mapdef Entity "IntoTheJungleCodeMod/SetColorGradeIfFlagOnSpawnController" SetColorGradeIfFlagOnSpawnController(x::Integer, y::Integer, flag::String="", colorGrade::String="")

const placements = Ahorn.PlacementDict(
    "Set Color Grade If Flag On Spawn Controller\n(Into The Jungle code mod)" => Ahorn.EntityPlacement(
        SetColorGradeIfFlagOnSpawnController
    )
)

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

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::SetColorGradeIfFlagOnSpawnController, room::Maple.Room) = Ahorn.drawImage(ctx, Ahorn.Assets.northernLights, -12, -12)

end

max_stars_repo_name: Jackal-Celeste/JungleHelper

max_issues_repo_name: Jackal-Celeste/JungleHelper

max_forks_repo_name: Jackal-Celeste/JungleHelper


text:
abstract type AbstractMethod end

export BSLSpline

struct BSLSpline <: AbstractMethod

    p :: Int

end