text:
using Plots
using Random
Random.seed!(1)

include_ddp()

include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)

function visualize!(vis, model, q;
      r = 0.025, Δt = 0.1)

	default_background!(vis)

	torso = Cylinder(Point3f0(0.0, 0.0, 0.0), Point3f0(0.0, 0.0, model.l_torso),
		convert(Float32, 0.035))
	setobject!(vis["torso"], torso,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	thigh_1 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l_thigh1),
		convert(Float32, 0.0175))
	setobject!(vis["thigh1"], thigh_1,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	calf_1 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l_calf1),
		convert(Float32, 0.0125))
	setobject!(vis["leg1"], calf_1,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	thigh_2 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l_thigh2),
		convert(Float32, 0.0175))
	setobject!(vis["thigh2"], thigh_2,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	calf_2 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l_calf2),
		convert(Float32, 0.0125))
	setobject!(vis["leg2"], calf_2,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	thigh_3 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l_thigh3),
		convert(Float32, 0.0175))
	setobject!(vis["thigh3"], thigh_3,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	calf_3 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l_calf3),
		convert(Float32, 0.0125))
	setobject!(vis["leg3"], calf_3,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	thigh_4 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l_thigh4),
		convert(Float32, 0.0175))
	setobject!(vis["thigh4"], thigh_4,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	calf_4 = Cylinder(Point3f0(0.0,0.0,0.0), Point3f0(0.0, 0.0, model.l_calf4),
		convert(Float32, 0.0125))
	setobject!(vis["leg4"], calf_4,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

	hip1 = setobject!(vis["hip1"], Sphere(Point3f0(0),
        convert(Float32, 0.035)),
        MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))

	hip2 = setobject!(vis["hip2"], Sphere(Point3f0(0),
        convert(Float32, 0.035)),
        MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))

	knee1 = setobject!(vis["knee1"], Sphere(Point3f0(0),
        convert(Float32, 0.025)),
        MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))

	knee2 = setobject!(vis["knee2"], Sphere(Point3f0(0),
		convert(Float32, 0.025)),
		MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))

	knee3 = setobject!(vis["knee3"], Sphere(Point3f0(0),
		convert(Float32, 0.025)),
		MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))

	knee4 = setobject!(vis["knee4"], Sphere(Point3f0(0),
        convert(Float32, 0.025)),
        MeshPhongMaterial(color = RGBA(0, 0, 0, 1.0)))

	feet1 = setobject!(vis["feet1"], Sphere(Point3f0(0),
        convert(Float32, r)),
        MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

	feet2 = setobject!(vis["feet2"], Sphere(Point3f0(0),
		convert(Float32, r)),
		MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

	feet3 = setobject!(vis["feet3"], Sphere(Point3f0(0),
		convert(Float32, r)),
		MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

	feet4 = setobject!(vis["feet4"], Sphere(Point3f0(0),
        convert(Float32, r)),
        MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0, 1.0)))

	T = length(q)
	p_shift = [0.0, 0.0, r]
	for t = 1:T
		MeshCat.atframe(anim, t) do
			p = [q[t][1]; 0.0; q[t][2]] + p_shift

			k_torso = kinematics_1(model, q[t], body = :torso, mode = :ee)
			p_torso = [k_torso[1], 0.0, k_torso[2]] + p_shift

			k_thigh_1 = kinematics_1(model, q[t], body = :thigh_1, mode = :ee)
			p_thigh_1 = [k_thigh_1[1], 0.0, k_thigh_1[2]] + p_shift

			k_calf_1 = kinematics_2(model, q[t], body = :calf_1, mode = :ee)
			p_calf_1 = [k_calf_1[1], 0.0, k_calf_1[2]] + p_shift

			k_thigh_2 = kinematics_1(model, q[t], body = :thigh_2, mode = :ee)
			p_thigh_2 = [k_thigh_2[1], 0.0, k_thigh_2[2]] + p_shift

			k_calf_2 = kinematics_2(model, q[t], body = :calf_2, mode = :ee)
			p_calf_2 = [k_calf_2[1], 0.0, k_calf_2[2]] + p_shift

			k_thigh_3 = kinematics_2(model, q[t], body = :thigh_3, mode = :ee)
			p_thigh_3 = [k_thigh_3[1], 0.0, k_thigh_3[2]] + p_shift

			k_calf_3 = kinematics_3(model, q[t], body = :calf_3, mode = :ee)
			p_calf_3 = [k_calf_3[1], 0.0, k_calf_3[2]] + p_shift

			k_thigh_4 = kinematics_2(model, q[t], body = :thigh_4, mode = :ee)
			p_thigh_4 = [k_thigh_4[1], 0.0, k_thigh_4[2]] + p_shift

			k_calf_4 = kinematics_3(model, q[t], body = :calf_4, mode = :ee)
			p_calf_4 = [k_calf_4[1], 0.0, k_calf_4[2]] + p_shift

			settransform!(vis["thigh1"], cable_transform(p, p_thigh_1))
			settransform!(vis["leg1"], cable_transform(p_thigh_1, p_calf_1))
			settransform!(vis["thigh2"], cable_transform(p, p_thigh_2))
			settransform!(vis["leg2"], cable_transform(p_thigh_2, p_calf_2))
			settransform!(vis["thigh3"], cable_transform(p_torso, p_thigh_3))
			settransform!(vis["leg3"], cable_transform(p_thigh_3, p_calf_3))
			settransform!(vis["thigh4"], cable_transform(p_torso, p_thigh_4))
			settransform!(vis["leg4"], cable_transform(p_thigh_4, p_calf_4))
			settransform!(vis["torso"], cable_transform(p, p_torso))
			settransform!(vis["hip1"], Translation(p))
			settransform!(vis["hip2"], Translation(p_torso))
			settransform!(vis["knee1"], Translation(p_thigh_1))
			settransform!(vis["knee2"], Translation(p_thigh_2))
			settransform!(vis["knee3"], Translation(p_thigh_3))
			settransform!(vis["knee4"], Translation(p_thigh_4))
			settransform!(vis["feet1"], Translation(p_calf_1))
			settransform!(vis["feet2"], Translation(p_calf_2))
			settransform!(vis["feet3"], Translation(p_calf_3))
			settransform!(vis["feet4"], Translation(p_calf_4))
		end
	end

	settransform!(vis["/Cameras/default"],
	    compose(Translation(0.0, 0.0, -1.0), LinearMap(RotZ(-pi / 2.0))))

	MeshCat.setanimation!(vis, anim)
end

contact_control_path = "/home/taylor/Research/ContactControl.jl/src"

using Parameters
# Utilities
include(joinpath(contact_control_path, "utils.jl"))

# Solver
include(joinpath(contact_control_path, "solver/cones.jl"))
include(joinpath(contact_control_path, "solver/interior_point.jl"))
include(joinpath(contact_control_path, "solver/lu.jl"))

# Environment
include(joinpath(contact_control_path, "simulator/environment.jl"))

# Dynamics
include(joinpath(contact_control_path, "dynamics/model.jl"))

# Simulator
include(joinpath(contact_control_path, "simulation/contact_methods.jl"))
include(joinpath(contact_control_path, "simulation/simulation.jl"))
include(joinpath(contact_control_path, "simulator/trajectory.jl"))

include(joinpath(contact_control_path, "dynamics/code_gen_dynamics.jl"))
include(joinpath(contact_control_path, "dynamics/fast_methods_dynamics.jl"))

# Models
include(joinpath(contact_control_path, "dynamics/quaternions.jl"))
include(joinpath(contact_control_path, "dynamics/mrp.jl"))
include(joinpath(contact_control_path, "dynamics/euler.jl"))

# include("dynamics/particle_2D/model.jl")
# include("dynamics/particle/model.jl")
# include(joinpath(contact_control_path, "dynamics/hopper_2D/model.jl"))
# include("dynamics/hopper_3D/model.jl")
# include("dynamics/hopper_3D_quaternion/model.jl")
include(joinpath(contact_control_path, "dynamics/quadruped/model.jl"))
# include("dynamics/quadruped_simple/model.jl")
# include("dynamics/biped/model.jl")
# include(joinpath(contact_control_path, "dynamics/flamingo/model.jl"))
# include("dynamics/pushbot/model.jl")
# include("dynamics/planarpush/model.jl")
# include("dynamics/planarpush_2D/model.jl")
# include("dynamics/rigidbody/model.jl")
# include("dynamics/box/model.jl")

# Simulation
include(joinpath(contact_control_path, "simulation/environments/flat.jl"))
# include("simulation/environments/piecewise.jl")
# include("simulation/environments/quadratic.jl")
# include("simulation/environments/slope.jl")
# include("simulation/environments/sinusoidal.jl")
# include("simulation/environments/stairs.jl")

include(joinpath(contact_control_path, "simulation/residual_approx.jl"))
include(joinpath(contact_control_path, "simulation/code_gen_simulation.jl"))

# Visuals
using MeshCatMechanisms
include(joinpath(contact_control_path, "dynamics/visuals.jl"))
include(joinpath(contact_control_path, "dynamics/visual_utils.jl"))

s = get_simulation("quadruped", "flat_2D_lc", "flat")

# @load joinpath(pwd(), "examples/contact_implicit", "flamingo_stand_100hz.jld2") q u γ b ψ η μ h
@load joinpath(pwd(), "examples/contact_implicit", "quadruped_gait_100Hz.jld2") qm um γm bm ψm ηm μm hm
visualize!(vis, s.model, qm, Δt = h)

nq = s.model.dim.q
m = s.model.dim.u

function initial_configuration(model::Quadruped, θ1, θ2, θ3)
    q1 = zeros(nq)
    q1[3] = pi / 2.0
    q1[4] = -θ1
    q1[5] = θ2

    q1[8] = -θ1
    q1[9] = θ2

    q1[2] = model.l_thigh1 * cos(q1[4]) + model.l_calf1 * cos(q1[5])

    q1[10] = -θ3
    q1[11] = acos((q1[2] - model.l_thigh2 * cos(q1[10])) / model.l_calf2)

    q1[6] = -θ3
    q1[7] = acos((q1[2] - model.l_thigh2 * cos(q1[6])) / model.l_calf2)

    return q1
end

θ1 = pi / 3.5
θ2 = pi / 3.5
θ3 = pi / 3.5

q1 = initial_configuration(s.model, θ1, θ2, θ3)
q0 = copy(q1)
visualize!(vis, s.model, [q1])

strd_diff = kinematics_2(s.model, qm[1], body = :calf_2, mode = :ee)[1] - kinematics_2(s.model, q1, body = :calf_2, mode = :ee)[1]

q0[1] += strd_diff
q1[1] += strd_diff

T = 5
h = 0.01

qT1 = qm[1]
qT = qm[2]

x1 = [q0; q1]
# xT = [q0; q1]#
xT = [qT1; qT]
visualize!(vis, s.model, [q0, q1], Δt = h)

visualize!(vis, s.model, [q1, qT], Δt = 1.0)

visualize!(vis, s.model, [qT1, qT], Δt = h)

struct Dynamics{T}
	s::Simulation
	ip_dyn::InteriorPoint
	ip_jac::InteriorPoint
	h::T
end

function gen_dynamics(s::Simulation, h;
		dyn_opts =  InteriorPointOptions{Float64}(
						r_tol = 1.0e-8,
						κ_tol = 1.0e-4,
						κ_init = 0.1,
						diff_sol = true),
		jac_opts =  InteriorPointOptions{Float64}(
						r_tol = 1.0e-8,
						κ_tol = 1.0e-4,
						κ_init = 0.1,
						diff_sol = true))

	z = zeros(num_var(s.model, s.env))
	θ = zeros(num_data(s.model))

	ip_dyn = interior_point(z, θ,
		idx_ineq = inequality_indices(s.model, s.env),
		r! = s.res.r!,
		rz! = s.res.rz!,
		rθ! = s.res.rθ!,
		rz = s.rz,
		rθ = s.rθ,
		opts = dyn_opts)

	ip_dyn.opts.diff_sol = false

	ip_jac = interior_point(z, θ,
		idx_ineq = inequality_indices(s.model, s.env),
		r! = s.res.r!,
		rz! = s.res.rz!,
		rθ! = s.res.rθ!,
		rz = s.rz,
		rθ = s.rθ,
		opts = jac_opts)

	ip_jac.opts.diff_sol = true

	Dynamics(s, ip_dyn, ip_jac, h)
end

d = gen_dynamics(s, h,
	dyn_opts = InteriorPointOptions{Float64}(κ_tol = 1.0e-5, κ_init = 0.1),
	jac_opts = InteriorPointOptions{Float64}(κ_tol = 1.0e-5, κ_init = 0.1))

function f!(d::Dynamics, q0, q1, u1, mode = :dynamics)
	s = d.s
	ip = (mode == :dynamics ? d.ip_dyn : d.ip_jac)
	h = d.h

	z_initialize!(ip.z, s.model, s.env, copy(q1))
	θ_initialize!(ip.θ, s.model, copy(q0), copy(q1), copy(u1), zeros(s.model.dim.w), s.model.μ_world, h)

	status = interior_point_solve!(ip)

	!status && (@warn "dynamics failure")
end

function f(d::Dynamics, q0, q1, u1)
	f!(d, q0, q1, u1, :dynamics)
	return copy(d.ip_dyn.z[1:d.s.model.dim.q])
end

f(d, q0, q1, zeros(m))

function fq0(d::Dynamics, q0, q1, u1)
	f!(d, q0, q1, u1, :jacobian)
	return copy(d.ip_jac.δz[1:d.s.model.dim.q, 1:d.s.model.dim.q])
end

fq0(d, q0, q1, zeros(m))

function fq1(d::Dynamics, q0, q1, u1)
	f!(d, q0, q1, u1, :jacobian)
	return copy(d.ip_jac.δz[1:d.s.model.dim.q, d.s.model.dim.q .+ (1:d.s.model.dim.q)])
end

fq1(d, q0, q1, zeros(m))

function fx1(d::Dynamics, q0, q1, u1)
	f!(d, q0, q1, u1, :jacobian)
	return copy(d.ip_jac.δz[1:d.s.model.dim.q, 1:(2 * d.s.model.dim.q)])
end

fx1(d, q0, q1, zeros(m))

function fu1(d::Dynamics, q0, q1, u1)
	f!(d, q0, q1, u1, :jacobian)
	return copy(d.ip_jac.δz[1:d.s.model.dim.q, 2 * d.s.model.dim.q .+ (1:d.s.model.dim.u)])
end

fu1(d, q0, q1, zeros(m))

struct QuadrupedCI{I, T} <: Model{I, T}
    n::Int
    m::Int
    d::Int
	dynamics::Dynamics
end

model = QuadrupedCI{Midpoint, FixedTime}(2 * s.model.dim.q, s.model.dim.u, 0, d)
u1 = zeros(model.m)

function fd(model::QuadrupedCI{Midpoint, FixedTime}, x, u, w, h, t)
	nq = model.dynamics.s.model.dim.q
	q0 = x[1:nq]
	q1 = x[nq .+ (1:nq)]

	q2 = f(model.dynamics, q0, q1, u)

	return [q1; q2]
end

fd(model, x1, u1, zeros(0), h, 1)

function fdx(model::QuadrupedCI{Midpoint, FixedTime}, x, u, w, h, t)
	nq = model.dynamics.s.model.dim.q
	q0 = x[1:nq]
	q1 = x[nq .+ (1:nq)]
	dq2dx1 = fx1(model.dynamics, q0, q1, u)

	return [zeros(nq, nq) I; dq2dx1]
end

fdx(model, x1, u1, zeros(0), h, 1)

function fdu(model::QuadrupedCI{Midpoint, FixedTime}, x, u, w, h, t)
	nq = model.dynamics.s.model.dim.q
	q0 = x[1:nq]
	q1 = x[nq .+ (1:nq)]
	dq2du1 = fu1(model.dynamics, q0, q1, u)
	return [zeros(nq, model.m); dq2du1]
end

fdu(model, x1, u1, zeros(0), h, 1)

n = model.n
m = model.m

T = 20
# ū = [u[1] + 0.0 * randn(model.m) for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Rollout
# @save joinpath(pwd(), "examples/contact_implicit", "flamingo_stand_100hz_v2.jld2") u_stand
@load joinpath(pwd(), "examples/contact_implicit", "quadruped_stand_100Hz.jld2") u_stand
u_stand = u_stand[1:m]
ū = [u_stand + 0.00 * randn(model.m) for t = 1:T-1]

x̄ = rollout(model, x1, ū, w, h, T)

q̄ = state_to_configuration(x̄)
visualize!(vis, s.model, q̄, Δt = h)

# Objective
Q = [t < T ? 1.0e-1 * Diagonal(ones(model.n)) : 10.0 * Diagonal(ones(model.n)) for t = 1:T]
q = [-2.0 * Q[t] * xT for t = 1:T]
R = [Diagonal(1.0 * ones(model.m)) for t = 1:T-1]
r = [-2.0 * R[t] * u_stand for t = 1:T-1]

obj = StageCosts([QuadraticCost(Q[t], q[t],
	t < T ? R[t] : nothing, t < T ? r[t] : nothing) for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
    if t < T
		Q = obj.cost[t].Q
		q = obj.cost[t].q
	    R = obj.cost[t].R
		r = obj.cost[t].r
        return x' * Q * x + q' * x + u' * R * u + r' * u
    elseif t == T
		Q = obj.cost[T].Q
		q = obj.cost[T].q
        return x' * Q * x + q' * x
    else
        return 0.0
    end
end

# Constraints
ul = -1.0 * ones(model.m)
uu = 1.0 * ones(model.m)
p = [t < T ? 2 * m : n for t = 1:T]
info_t = Dict(:ul => ul, :uu => uu, :inequality => (1:2 * m))
info_T = Dict(:xT => xT)
con_set = [StageConstraint(p[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	p = cons.con[t].p

	if t < T
		ul = cons.con[t].info[:ul]
		uu = cons.con[t].info[:uu]
		c .= [ul - u; u - uu]
	else
		c .= x - cons.con[T].info[:xT]
	end
end

prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T,
	analytical_dynamics_derivatives = true)

# Solve
@time constrained_ddp_solve!(prob,
	max_iter = 100, max_al_iter = 2,
	ρ_init = 1.0, ρ_scale = 10.0,
	con_tol = 1.0e-2)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

q̄ = state_to_configuration(x̄)

vis = Visualizer()
render(vis)
visualize!(vis, s.model, q̄, Δt = h)

plot(hcat(ū...)', linetype = :steppost)

norm(x̄[end] - xT, Inf)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using Random

"""

Arguments:
    - `f`: Function to be optimized
    - `g`: Gradient function for `f`
    - `c`: Constraints function. Evaluates all the constraint equations and return a Vector of values
    - `x0`: (Vector) Initial position to start from
    - `n`: (Int) Number of evaluations allowed. Remember `g` costs twice of `f`
    - `prob`: (String) Name of the problem. So you can use a different strategy for each problem
"""
function optimize(f, g, c, x0, n, prob)
    p = squared_sum_penalty_function(c)
    x_best, xlist = penalty_method(f, p, x0, n)
    return x_best
end

function sum_penalty_function(c)
    f = function(x)
        ceval = c(x)
        sum(ceval[ceval .> 0])
    end
    return f
end

function squared_sum_penalty_function(c)
    f = function(x)
        ceval = c(x)
        sum(ceval[ceval .> 0].^2)
    end
    return f
end

function penalty_method(f, p, x, n; ρ=1, γ=2)
    numEvalsLeft = n
    hj_max_eval = 10;
    xlist = Array{Float64}[]
    push!(xlist, x)
    while numEvalsLeft >= hj_max_eval
        x, evalsUsed, feasFound = hooke_jeeves(f, ρ, p, x, hj_max_eval)
        push!(xlist, x)
        numEvalsLeft = numEvalsLeft - evalsUsed
        ρ *= γ
    end
    return x,xlist
end

function last_feasible_x(currXFeas, x, feasFound, last_x_feas)
    if currXFeas
        return x
    else
        if !feasFound
            return x
        else
            return last_x_feas
        end
    end
end

basis(i, n) = [k == i ? 1.0 : 0.0 for k in 1 : n]

function hooke_jeeves(f, ρ, p, x, maxEval; α=1, γ=0.5)
    xf = f(x)
    xp = p(x)
    y, m = xf + ρ*xp, length(x)
    evalsUsed = 1
    last_x_feas = x
    if xp==0
        feasFound = true
        currXFeas = true
    else
        feasFound = false
        currXFeas = false
    end
    while evalsUsed < maxEval
        improved = false
        for i in 1 : m
            x′ = x + α*basis(i, m)
            if evalsUsed >= maxEval
                return last_feasible_x(currXFeas, x, feasFound, last_x_feas), evalsUsed, feasFound
            end
            xf′ = f(x′)
            xp′ = p(x′)
            y′ = xf′ + ρ*xp′
            evalsUsed = evalsUsed + 1
            if y′ < y
                x, y, improved = x′, y′, true
                if xp′==0
                    last_x_feas = x
                    feasFound = true
                    currXFeas = true
                else
                    currXFeas = false
                end
            else
                x′ = x - α*basis(i, m)
                if evalsUsed >= maxEval
                    return last_feasible_x(currXFeas, x, feasFound, last_x_feas), evalsUsed, feasFound
                end
                xf′ = f(x′)
                xp′ = p(x′)
                y′ = xf′ + ρ*xp′
                evalsUsed = evalsUsed + 1
                if y′ < y
                    x, y, improved = x′, y′, true
                    if xp′==0
                        last_x_feas = x
                        feasFound = true
                        currXFeas = true
                    else
                        currXFeas = false
                    end
                end
            end
        end
        if !improved
            α *= γ
        end
    end
    return last_x_feas, evalsUsed, feasFound
end

ext: jl


text:
# Testing
# =======
#
# Utilities to assist testing of BioJulia packages.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/BioCore.jl/blob/master/LICENSE.md

module Testing

"""
    random_array(n::Integer, elements, probs)

Create a random array of length `n`, composed of possible `elements`, which are
selected according to their `prob`.
"""
function random_array(n::Integer, elements, probs)
    @assert sum(probs) == 1.0 "the probabilities in probs, don't sum to 1.0!"
    cumprobs = cumsum(probs)
    x = Vector{eltype(elements)}(undef, n)
    for i in 1:n
        x[i] = elements[searchsorted(cumprobs, rand()).start]
    end
    return x
end

"""
    random_seq(n::Integer, nts, probs)

Create a random sequence of length `n`, composed of `nts`, according to their
`prob`.
"""
function random_seq(n::Integer, nts, probs)
    return String(random_array(n, nts, probs))
end

"""
    random_dna(n, probs = [0.24, 0.24, 0.24, 0.24, 0.04])

Create a random DNA sequence of length `n`, by sampling the nucleotides
A, C, G, T, and N, according to their probability in `probs`.
"""
function random_dna(n, probs = [0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'T', 'N'], probs)
end

"""
    random_rna(n, probs = [0.24, 0.24, 0.24, 0.24, 0.04])

Create a random RNA sequence of length `n`, by sampling the nucleotides
A, C, G, U, and N, according to their probability in `probs`.
"""
function random_rna(n, probs = [0.24, 0.24, 0.24, 0.24, 0.04])
    return random_seq(n, ['A', 'C', 'G', 'U', 'N'], probs)
end

"""
    random_aa(n, probs = [0.24, 0.24, 0.24, 0.24, 0.04])

Create a random amino acid sequence of length `n`, by sampling the possible
amino acid characters.
"""
function random_aa(len)
    return random_seq(len,
        ['A', 'R', 'N', 'D', 'C', 'Q', 'E', 'G', 'H', 'I',
         'L', 'K', 'M', 'F', 'P', 'S', 'T', 'W', 'Y', 'V', 'X' ],
        push!(fill(0.049, 20), 0.02))
end

"""
    intempdir(fn::Function, parent = tempdir())

Execute some function `fn` in a temporary directory.
After the function is executed, the directory is cleaned by doing the
equivalent of `rm -r`.
"""
function intempdir(fn::Function, parent = tempdir())
    dirname = mktempdir(parent)
    try
        cd(fn, dirname)
    finally
        rm(dirname, recursive = true)
    end
end

"""
    random_interval(minstart, maxstop)

Create a random interval between a minimum starting point, and a maximum stop
point.
"""
function random_interval(minstart, maxstop)
    start = rand(minstart:maxstop)
    return start:rand(start:maxstop)
end

end # Module Testing

max_forks_repo_name: BioJulia/BioGenerics


text:
using Distributions, Plots; pyplot()

N = 10^5
nMin, nStep, nMax = 10, 10, 200
nn = Int(nMax/nStep)
sampleSizes = nMin:nStep:nMax
trueB = 5
trueDist = Uniform(-2, trueB)

MLEest(data) = maximum(data)
MMest(data)  = mean(data) + sqrt(3)*std(data)

res = Dict{Symbol,Array{Float64}}(
    ((sym) -> sym => Array{Float64}(undef,nn)).(
        [:MSeMLE,:MSeMM, :VarMLE,:VarMM,:BiasMLE,:BiasMM]))

for (i, n) in enumerate(sampleSizes)
    mleEst, mmEst = Array{Float64}(undef, N), Array{Float64}(undef, N)
    for j in 1:N
        sample    = rand(trueDist,n)
        mleEst[j] = MLEest(sample)
        mmEst[j]  = MMest(sample)
    end
    meanMLE, meanMM = mean(mleEst), mean(mmEst)
    varMLE, varMM = var(mleEst), var(mmEst)

    res[:MSeMLE][i] = varMLE + (meanMLE - trueB)^2
    res[:MSeMM][i] = varMM + (meanMM - trueB)^2
    res[:VarMLE][i] = varMLE
    res[:VarMM][i] = varMM
    res[:BiasMLE][i] = meanMLE - trueB
    res[:BiasMM][i] = meanMM - trueB
end

p1 = scatter(sampleSizes, [res[:MSeMLE] res[:MSeMM]], c=[:blue :red],
    label=["Mean sq.err (MLE)" "Mean sq.err (MM)"])
p2 = scatter(sampleSizes, [res[:VarMLE] res[:VarMM]], c=[:blue :red],
    label=["Variance (MLE)" "Variance (MM)"])
p3 = scatter(sampleSizes, [res[:BiasMLE] res[:BiasMM]], c=[:blue :red],
    label=["Bias (MLE)" "Bias (MM)"])

plot(p1, p2, p3, ms=10, shape=:xcross, xlabel="n",
    layout=(1,3), size=(1200, 400))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using BinDeps

depsdir = @__DIR__

# Build TreeSitter
BinDeps.run(@build_steps begin
    ChangeDirectory(joinpath(depsdir, "tree-sitter"))
    `script/build-lib`
end)

# # Build Your Grammar File
# BinDeps.run(@build_steps begin
#     ChangeDirectory(depsdir)
#     `./generate.sh`
#     `make clean all`
# end)


text:
# PerceptronLoss

using LaTeXStrings, LossFunctions, Measures, Latexify

Plots.plot(PerceptronLoss(), c=:olive, linewidth=2, annotation=(1, 1.6, L"L\left( a \right) = max\{0, - a \}", 18))

max_stars_repo_name: Ellon-M/visualizations

max_issues_repo_name: Ellon-M/visualizations

max_forks_repo_name: Ellon-M/visualizations


text:
module TestEffectiveNumbers
using Test
using LinearAlgebra
using Statistics

using Diversity
using Diversity: powermean

numbers = [1.0, 2, 4, 8, 16];
numspecies = 100;
fragments = rand(numspecies);
weights = rand(numspecies);
weights /= sum(weights);
Z1 = ones(typeof(weights[1]), (length(weights), length(weights)));
numcommunities = 8;
manyweights = rand(numspecies, numcommunities);
manyweights *= Diagonal(reshape(mapslices(v -> 1. / sum(v), manyweights, dims=1),
                                numcommunities));

# Simple power means - we no longer export these, but we should check
# them anyway as everything relies on them

@testset "powermean" begin
    # Check that an exception is thrown when 'values' and weights are different lengths
    @test_throws DimensionMismatch powermean(numbers, 0, weights)

    # Some simple values
    @test powermean([1.0], 0.0, [1.0]) ≈ 1.0
    @test powermean(numbers, 0.0) ≈ 4.0
    @test powermean(numbers, [-Inf]) ≈ [1]
    @test powermean(numbers, [1.0, -1.0]) ≈ [31.0/5.0, 80.0/31.0]
    @test powermean(numbers, Inf, [1.0, 1.0, 1.0, 1.0, 0.0]) ≈ 8
    @test isnan(powermean(numbers, 0.0, 0.0 * numbers))

    # Power mean with some random numbers
    @test powermean(fragments, 0) ≈ prod(fragments .^ (1. / numspecies))
    @test powermean(fragments, 1) ≈ mean(fragments)
    @test powermean(fragments, Inf) ≈ maximum(fragments)
    @test powermean(fragments, 0, weights) ≈ prod(fragments .^ weights)
    @test powermean(fragments, 1, weights) ≈ sum(fragments .* weights)
    @test powermean(manyweights, -1, manyweights) .^ -1 ≈ numspecies * ones(size(manyweights, 2))
end

@testset "qD" begin
    # Basic qD diversity calculation
    @test qD(weights, 0) ≈ mapreduce((x) -> x ≈ 0 ? 0 : 1, +, weights)
    @test qD(Metacommunity(weights), 0) == qD(weights, 0)
    @test qD(Metacommunity(weights, UniqueTypes(numspecies)), 0) ==
        qD(weights, 0)
    @test_throws ErrorException qD(Metacommunity(weights,
                                   GeneralTypes(rand(numspecies, numspecies))),
                                   0)
    @test qD(weights, 1) ≈ prod(weights .^ -weights)
    @test qD(weights, 2) ≈ 1.0 / sum(weights .^ 2)
    @test qD(weights, Inf) ≈ 1.0 / maximum(weights)

    @test qD(weights, [1, 2]) ≈ [qD(weights, 1), qD(weights, 2)]

    @test typeof(qD(manyweights[:,1], 0)) <: AbstractFloat
    @test typeof(qD(manyweights[:,1], [0])) <: Vector

    for i in 1:size(manyweights, 2)
        @test qD(manyweights[:,i], [0]) ≈ numspecies * ones((1, size(manyweights[:,i], 2)))
    end
end

@testset "qDZ" begin
    # General Leinster-Cobbold diversity calculation
    @test qDZ(weights, [1, 2]) ≈ qD(weights, [1, 2])
    @test qDZ(weights, [0, 1, 2, 3, Inf], Z1) ≈ [1, 1, 1, 1, 1]

    for i in 1:size(manyweights, 2)
        @test qDZ(manyweights[:,i], [0, 1, 2, Inf],
                  ones((size(manyweights[:,i], 1),
                        size(manyweights[:,i], 1)))) ≈ ones((4, size(manyweights[:,i], 2)))
    end
end

end


text:
module DPHTTest
using DataProcessingHierarchyTools
const DPHT = DataProcessingHierarchyTools
import DataProcessingHierarchyTools:level, filename, shash, _hash
using Test
using LinearAlgebra

import Base:hcat, zero

struct TestData <: DPHData
    data::String
    setid::Vector{Int64}
end

zero(::Type{TestData}) = TestData("", Int64[])

DPHT.level(::Type{TestData}) = "session"
DPHT.filename(::Type{TestData}) = "data.mat"

function Base.hcat(X1::TestData, X2::TestData)
    ndata = "$(X1.data)$(X2.data)"
    setid = [X1.setid;X2.setid .+ X1.setid[end]]
    TestData(ndata, setid)
end

function TestData()
    open(filename(TestData), "r") do f
        data = readline(f)
        ss = readline(f)
        setid = [parse(Int64,s) for s in split(ss, ',')]
        TestData(data, setid)
    end
end

function TestData(x)
    X = TestData("test", [1])
    save_data(X)
    X
end

function save_data(X::TestData)
    open(filename(TestData),"w") do f
        write(f, "$(X.data)\n")
        for ss in X.setid[1:end-1]
            write(f, "$(ss),")
        end
        write(f, "$(X.setid[end])\n")
    end
end

dirs =["20140903/session01", "20140903/session02"]

@testset "Level functions" begin
    thislevel = DPHT.level("newWorkingMemory/Pancake/20130923")
    @test thislevel == "day"
    _name = DPHT.get_level_name("day","newWorkingMemory/Pancake/20130923")
    @test _name == "20130923"
    _name = DPHT.get_level_name("subject","newWorkingMemory/Pancake/20130923")
    @test _name == "Pancake"
    _pth = DPHT.process_level("session", "Pancake/20130923/session01/array01/channel001")
    @test _pth == "./../.."
    rpath = DPHT.get_relative_path("subjects", "newWorkingMemory/Pancake/20130923")
    @test rpath == "Pancake/20130923"
    rpath = DPHT.get_relative_path("day", "newWorkingMemory/Pancake/20130923")
    @test rpath == "."
    @test_throws ErrorException DPHT.get_relative_path("channel", "newWorkingMemory/Pancake/20130923")

    thislevel = "newWorkingMemory/Pancake/20130923/session01/array01/channel001/cell01"
    @test DPHT.get_level_path("array", thislevel) == "newWorkingMemory/Pancake/20130923/session01/array01"
    @test DPHT.get_level_path("session", thislevel) == "newWorkingMemory/Pancake/20130923/session01"

    dd = tempdir()
    cd(dd) do
        for d2 in dirs
            mkpath(d2)
        end
        func() = touch("test.txt")
        DPHT.process_dirs(func, dirs)
        for d2 in dirs
            @test isfile("$(d2)/test.txt")
        end
        level_dirs = DPHT.get_level_dirs("session", dirs)
        @test level_dirs == dirs
        #cleanup
        cd(dirs[1]) do
            thislevel = DPHT.level()
            @test thislevel == "session"
            #test getting all level directories belonging to a particular level
            session_dirs = DPHT.get_level_dirs("session","..")
            @test session_dirs[1] == "../session01"
        end
        dirsn = DPHT.get_level_dirs("session",dirs[1])
        @test dirsn[1] == "20140903/session01/."
        dirsn = DPHT.get_level_dirs("session")
        @test dirsn[1] == dirs[1]
        @test dirsn[2] == dirs[2]
        for d2 in dirs
            rm(d2;recursive=true)
        end
    end
    ss = "/tmp/temp2/Pancake/20130923/session01/array01/channel001/cell01"
    pp = DPHT.get_fullname(ss)
    @test pp == joinpath("Pancake","20130923","session01","array01","channel001","cell01")
    ss = "/tmp/temp2/Pancake"
    pp = DPHT.get_fullname(ss)
    @test pp == "Pancake"
end

@testset "Types" begin
    dd = tempdir()
    cd(dd) do
        for d2 in dirs
            mkpath(d2)
        end
        Y = DPHT.process_dirs(TestData, dirs, 1)
        for d2 in dirs
            cd(d2) do
                @test isfile(filename(TestData))
                ll = DPHT.process_level(TestData)
                @test ll == "."
            end
        end
        @test Y.data == "testtest"
        @test Y.setid == [1,2]
    end
end

@testset "Graceful fail" begin
    dd = tempdir()
    mkpath(dirs[1])
    cd(dirs[1]) do
        @test_throws ArgumentError DPHT.process_level("rubbish")
    end
    rm(dirs[1];recursive=true)
    @test_throws ErrorException DPHT.level(DPHT.DPHData)
    @test_throws ErrorException DPHT.filename(DPHT.DPHData)
end

@testset "Shortnames" begin
    sn = DPHT.get_shortname("cell01")
    @test sn == "c01"
    sn = DPHT.get_shortname("newWorkingMemory/James/20140904/session01/array01/channel030/cell01")
    @test sn == "J20140904s01a01g030c01"
end

struct MyArgs <: DPHT.DPHDataArgs
    f1::Float64
    f2::Int64
    f3::AbstractVector{Float64}
end

struct MyData <: DPHT.DPHData
    dd::Float64
    args::MyArgs
end

DPHT.filename(::Type{MyData}) = "mydata.mat"
DPHT.datatype(::Type{MyArgs}) = MyData
DPHT.level(::Type{MyData}) = "session"

function MyData(args::MyArgs;force_redo=false, do_save=true)
    fname = DPHT.filename(args)
    redo = force_redo || !isfile(fname)
    if !redo
        X = DPHT.load(args)
    else
        x = args.f1*args.f2*sum(args.f3)
        X = MyData(x, args)
        if do_save
            DPHT.save(X)
        end
    end
    X
end

@testset "ArgsHash" begin
    args = MyArgs(1.0, 3, -1.0:0.5:10.0)
    h = shash(args)
    @test h == 0x34c727e2f750c253
    @test DPHT.filename(args) == "mydata_34c727e2f750c253.mat"
	h2 = _hash(args)
	if Base.VERSION < v"1.5.0"
		@test h2 == h
	elseif Base.VERSION < v"1.6.0"
		@test h2 == 0x34c727e2f750c253
	end
end

@testset "ArgsCheck" begin
    args = MyArgs(1.0, 3, -1.0:0.5:10.0)
    @test DPHT.check_args(args, 1.0, 3, -1.0:0.5:10.0)
    @test DPHT.check_args(args, 1.0, 3, collect(-1.0:0.5:10.0))
    @test DPHT.check_args(args, 2.0, 3, -1.0:0.5:10.0) == false
    @test DPHT.check_args(args, 1.0, 2, -1.0:0.5:10.0) == false
    @test DPHT.check_args(args, 1.0, 3, -1.0:1.0:10.0) == false
    args2 = MyArgs(1.0, 3, -1.0:0.5:10.0)
    @test DPHT.check_args(args, args2)
    args3 = MyArgs(1.0, 3, -1.0:0.1:10.0)
    @test DPHT.check_args(args, args3) == false
    args4 = MyArgs(1.0, 4, -1.0:0.1:10.0)
    @test DPHT.check_args(args, args4) == false
end

@testset "Save and load" begin
    dd = tempdir()
    args = MyArgs(1.0, 3, -1.0:0.5:10.0)
    cd(dd) do
        DPHT.visit_dirs(MyData, dirs, args)
        cd(dirs[1]) do
            X = MyData(args;do_save=true, force_redo=false)
            mkdir("array01")
            cd("array01") do
                cd(DPHT.process_level(DPHT.level(MyData))) do
                    X2 = MyData(args;do_save=true, force_redo=false)
                    @test X.dd ≈ X2.dd
                end
            end
        end
        cd(dirs[2]) do
            rm(DPHT.filename(args))
            mkdir("array01")
            cd("array01") do
                cd(DPHT.process_level(DPHT.level(MyData))) do
                    @test_throws ErrorException q = DPHT.load(TestData)
                end
            end
        end

        for d2 in dirs
            rm(d2;recursive=true)
        end
        args2 = MyArgs(1.0, 3, [1.0])
        X = MyData(args2)
        fname = DPHT.filename(args2)
        @test isfile(fname)
        X2 = MyData(args2)
        @test X.args.f3 == X2.args.f3
        rm(fname)
    end
end

struct SArgs <: DPHT.DPHDataArgs
    a::Int64
end

struct S <: DPHT.DPHData
    μ::Vector{Float64}
    Σ::Matrix{Float64}
    α0::Float64
    args::SArgs
end

DPHT.filename(::Type{S}) = "mydata.mat"
DPHT.datatype(::Type{SArgs}) = S

@testset "Variable sanitasion" begin
    ss = S([0.0, 0.0], Matrix{Float64}(I,2,2), 0.05, SArgs(1.0))
    Q = convert(Dict{String,Any}, ss)
    @test Q["mu"] ≈ ss.μ
    @test Q["Sigma"] ≈ ss.Σ
    @test Q["alpha0"] ≈ ss.α0

    ss2 = convert(S, Q)
    @test ss2.μ ≈ ss.μ
    @test ss2.Σ ≈ ss.Σ
    @test ss2.α0 ≈ ss.α0

    Q = Dict{String,Any}("mu" => 1.0, "Sigma" => 0.5,
                         "alpha0" => 0.05,
                         "args" => Dict{String,Any}("a" => 3))
    ss3 = convert(S, Q)
    @test ss3.μ ≈ [1.0]
    @test ss3.Σ ≈ fill(0.5, 1,1)
    @test ss3.α0 ≈ 0.05

    ss2 = DPHT.desanitise("mu0")
    @test ss2 == "μ0"
    ss2 = DPHT.desanitise("mu_f")
    @test ss2 == "μ_f"
    ss2 = DPHT.desanitise("alpha_0f")
    @test ss2 == "α_0f"
end

@testset "Array of objects" begin
    args = [SArgs(1.0), SArgs(2.0), SArgs(3.0)]
    ss = S[]
    for a in args
        push!(ss, S([0.0, 0.0], [[0.1 0.1];[0.3 0.1]], 0.05, a))
    end
    dd = tempdir()
    cd(dd) do
        DPHT.save(ss)
        ss2 = DPHT.load(args)
        @test ss2[1].args == args[1]
        @test ss2[1].μ == ss[1].μ
        @test ss2[1].Σ == ss[1].Σ
        @test ss2[2].args == args[2]
        @test ss2[2].μ == ss[2].μ
        @test ss2[2].Σ == ss[2].Σ
        @test ss2[3].args == args[3]
        @test ss2[3].μ == ss[3].μ
        @test ss2[3].Σ == ss[3].Σ
        fname = DPHT.filename(args)
        h = shash(args)
        @test h == 0xb6f003559185097d
        rm(fname)
		for _ss in ss
			DPHT.save(_ss, overwrite=true)
		end
		args2 = DPHT.findargs(S)
		@test length(args2) == length(args)
		@test findfirst(x->x.a==args[1].a,args2) != nothing
		@test findfirst(x->x.a==args[2].a,args2) != nothing
		@test findfirst(x->x.a==args[3].a,args2) != nothing
		fargs = DPHT.findargs(S;a=2)
		@test length(fargs) == 1
    end
end

@testset "git-annex" begin
    tdir = tempdir()
    s_args = SArgs(2)
    ss = S([0.1], [0.0 1.0; 1.0 0.0], 0.05, s_args)
    cd(tdir) do
        mkpath("testdata")
        cd("testdata") do
            DPHT.save(ss, overwrite=true)
            if DPHT.git_annex() != nothing
                run(`git init`)
                run(`$(DPHT.git_annex()) init .`)
                run(`$(DPHT.git_annex()) add .`)
                run(`$(DPHT.git_annex()) lock`)
                ss2 = S([0.1], [0.0 1.0; 1.0 0.0], 0.03, s_args)
                mkpath("../testdata2")
                cd("../testdata2") do
                    run(`git init --bare `)
                end
                run(`git remote add origin ../testdata2`)
                run(`$(DPHT.git_annex()) sync origin`)
                run(`$(DPHT.git_annex()) copy -t origin .`)
                @test isfile(DPHT.filename(ss2.args))
                @test_throws ErrorException DPHT.save(ss2)
                run(`$(DPHT.git_annex()) drop .`)
                DPHT.reset!(s_args)
                DPHT.save(ss2, overwrite=true)
                run(`$(DPHT.git_annex()) sync origin`)
                run(`$(DPHT.git_annex()) copy -t origin .`)
                run(`$(DPHT.git_annex()) drop .`)
            end
            @test DPHT.computed(s_args)
            ss2 = DPHT.load(s_args)
            @test ss2.args.a == ss.args.a
            #cleanup
            if DPHT.git_annex() != nothing
                fname = DPHT.filename(s_args)
                run(`$(DPHT.git_annex()) drop --force --from origin $fname`)
                run(`$(DPHT.git_annex()) drop --force $fname`)
            end
        end
        #cleanup
        chmod("testdata", 0o777, recursive=true)
        rm("testdata", recursive=true)
        if isdir("testdata2")
            chmod("testdata2", 0o777, recursive=true)
            rm("testdata2", recursive=true)
        end
    end
end
end#module


text:
"""
    LookupArrayTrait

Abstract supertype of all traits of a [`LookupArray`](@ref).

These modify the behaviour of the lookup index.

The term "Trait" is used loosely - these may be fields of an object
of traits hard-coded to specific types.
"""
abstract type LookupArrayTrait end

"""
    Order <: LookupArrayTrait

Traits for the order of a [`LookupArray`](@ref). These determine how
`searchsorted` finds values in the index, and how objects are plotted.
"""
abstract type Order <: LookupArrayTrait end

"""
    Ordered <: Order

Supertype for the order of an ordered [`LookupArray`](@ref),
including [`ForwardOrdered`](@ref) and [`ReverseOrdered`](@ref).
"""
abstract type Ordered <: Order end

"""
    AutoOrder <: Order

    AutoOrder()

Specifies that the `Order` of a `LookupArray` will be found automatically
where possible.
"""
struct AutoOrder <: Order end

"""
    ForwardOrdered <: Ordered

    ForwardOrdered()

Indicates that the `LookupArray` index is in the normal forward order.
"""
struct ForwardOrdered <: Ordered end

"""
    ReverseOrdered <: Ordered

    ReverseOrdered()

Indicates that the `LookupArray` index is in the reverse order.
"""
struct ReverseOrdered <: Ordered end

"""
    Unordered <: Order

    Unordered()

Indicates that `LookupArray` is unordered.

This means the index cannot be searched with `searchsortedfirst`
or similar optimised methods - instead it will use `findfirst`.
"""
struct Unordered <: Order end

isrev(x) = isrev(typeof(x))
isrev(::Type{<:ForwardOrdered}) = false
isrev(::Type{<:ReverseOrdered}) = true

"""
   Locus <: LookupArrayTrait

Abstract supertype of types that indicate the position of index values
where they represent [`Intervals`](@ref).

These allow for values array cells to align with the [`Start`](@ref),
[`Center`](@ref), or [`End`](@ref) of values in the lookup index.

This means they can be plotted with correct axis markers, and allows automatic
converrsions to between formats with different standards (such as NetCDF and GeoTiff).
"""
abstract type Locus <: LookupArrayTrait end

"""
    Center <: Locus

    Center()

Indicates a lookup value is for the center of its corresponding array cell.
"""
struct Center <: Locus end

"""
    Start <: Locus

    Start()

Indicates a lookup value is for the start of its corresponding array cell,
in the direction of the lookup index order.
"""
struct Start <: Locus end

"""
    End <: Locus

    End()

Indicates a lookup value is for the end of its corresponding array cell,
in the direction of the lookup index order.
"""
struct End <: Locus end

"""
    AutoLocus <: Locus

    AutoLocus()

Indicates a interval where the index position is not yet known.
This will be filled with a default value on object construction.
"""
struct AutoLocus <: Locus end

"""
    Sampling <: LookupArrayTrait

Indicates the sampling method used by the index: [`Points`](@ref)
or [`Intervals`](@ref).
"""
abstract type Sampling <: LookupArrayTrait end

struct NoSampling <: Sampling end
struct AutoSampling <: Sampling end

locus(sampling::AutoSampling) = AutoLocus()

"""
    Points <: Sampling

    Points()

[`Sampling`](@ref) lookup where single samples at exact points.

These are always plotted at the center of array cells.
"""
struct Points <: Sampling end

locus(sampling::Points) = Center()

"""
    Intervals <: Sampling

    Intervals(locus::Locus)

[`Sampling`](@ref) specifying that sampled values are the mean (or similar)
value over an _interval_, rather than at one specific point.

Intervals require a [`Locus`](@ref) of [`Start`](@ref), [`Center`](@ref) or
[`End`](@ref) to define the location in the interval that the index values refer to.
"""
struct Intervals{L} <: Sampling
    locus::L
end
Intervals() = Intervals(AutoLocus())

locus(sampling::Intervals) = sampling.locus

"""
    rebuild(::Intervals, locus::Locus) => Intervals

Rebuild `Intervals` with a new Locus.
"""
rebuild(::Intervals, locus) = Intervals(locus)

"""
    Span <: LookupArrayTrait

Defines the type of span used in a [`Sampling`](@ref) index.
These are [`Regular`](@ref) or [`Irregular`](@ref).
"""
abstract type Span <: LookupArrayTrait end

struct NoSpan <: Span end

"""
    AutoSpan <: Span

    AutoSpan()

The span will be guessed and replaced in `format` or `set`.
"""
struct AutoSpan <: Span end

struct AutoStep end
struct AutoBounds end

"""
    Regular <: Span

    Regular(step=AutoStep())

`Points` or `Intervals` that have a fixed, regular step.
"""
struct Regular{S} <: Span
    step::S
end
Regular() = Regular(AutoStep())

val(span::Regular) = span.step

Base.step(span::Regular) = span.step

"""
    Irregular <: Span

    Irregular(bounds::Tuple)
    Irregular(lowerbound, upperbound)

`Points` or `Intervals` that have an `Irrigular` step size. To enable bounds tracking
and accuract selectors, the starting bounds are provided as a 2 tuple, or 2 arguments.
`(nothing, nothing)` is acceptable input, the bounds will be guessed from the index,
but may be innaccurate.
"""
struct Irregular{B<:Union{<:Tuple{<:Any,<:Any},AutoBounds}} <: Span
    bounds::B
end
Irregular() = Irregular(AutoBounds())
Irregular(lowerbound, upperbound) = Irregular((lowerbound, upperbound))

bounds(span::Irregular) = span.bounds
val(span::Irregular) = span.bounds

"""
    Explicit(bounds::AbstractMatix)

Intervals where the span is explicitly listed for every interval.

This uses a matrix where with length 2 columns for each index value,
holding the lower and upper bounds for that specific index.
"""
struct Explicit{B} <: Span
    val::B
end
Explicit() = Explicit(AutoBounds())

val(span::Explicit) = span.val

"""
    AutoIndex

Detect a `LookupArray` index from the context. This is used in `NoLookup` to simply
use the array axis as the index when the array is constructed, and in `set` to
change the `LookupArray` type without changing the index values.
"""
struct AutoIndex <: AbstractVector{Int} end

Base.size(::AutoIndex) = (0,)


text:
using DashBootstrapComponents

pagination = html_div(dbc_pagination(max_value = 10))

max_issues_repo_name: sthagen/dash-bootstrap-components

max_forks_repo_name: sthagen/dash-bootstrap-components


text:
using Test

@testset "Hierarchical Sampling" begin

maxT = 100
nSamples = 10
its = 5000

testEvents = Array{Array{Float64}}(undef, nSamples)
for i in 1:nSamples
    testEvents[i] = HawkesProcesses.simulate(0.5, 0.5, Distributions.Exponential(1/0.5), maxT)
end

testSamples = HawkesProcesses.hierarchical_fit(testEvents, maxT, its)

@test isa(testSamples, Tuple)
@test length(testSamples) == 3
@test length(testSamples[1]) == its
@test length(testSamples[2]) == its
@test length(testSamples[3]) == its

end


text:
using Weave

const DEST_DIR = joinpath(@__DIR__, "..", "src", "experiments")

for (root, dirs, files) in walkdir(joinpath(@__DIR__, "experiments"))
    for f in files
        if splitext(f)[2] == ".jl"
            tangle(joinpath(root, f); informat="script", out_path=DEST_DIR)
        end
    end
end


text:
"""
    mcppv(ts::AbstractVector; shuffle::Int = 500)

return p-value of a change point in mean in `ts`.\\
`shuffle` is number of randomization used to generate p-value,
uses global random number generator.
"""
function mcppv(ts::AbstractVector; shuffle::Int = 500)
    length(ts) <= 1 && return 1.0                   # large pvalue is no change

    dev = ts .- mean(ts)                            # centered
    cusum = Vector{Float64}(undef, length(ts))
    cumsum!(cusum, dev)                             # observed cumsum
    cmax_obs = maximum(abs, cusum)                  # max abs cumsum observed

    n = 0
    for _ = 1:shuffle
        shuffle!(dev)
        cumsum!(cusum, dev)
        cmax = maximum(abs, cusum)
        if cmax >= cmax_obs
            n += 1
        end
    end

    (n + 1) / (shuffle + 1)                         # p-value
end

"""
    mcptime(ts::AbstractVector; minlen::Int = 1)

return index of start of new segment via minimum ssq in `ts`.\\
`minlen` is minimum length of a segment you want.\\
0 is returned if `minlen` is not least twice length of `ts`.
"""
function mcptime(ts::AbstractVector; minlen::Int = 3)
    minlen < 1 && return 0
    len = length(ts)
    len < 2 * minlen && return 0

    chgpoint = 0
    minssq = Inf
    @inbounds for brk = minlen+1:len+1-minlen
        @views ssqtot = ssq(ts[1:brk-1]) + ssq(ts[brk:len])
        if ssqtot < minssq
            minssq = ssqtot
            chgpoint = brk
        end
    end

    chgpoint
end

"""
ssq(v) returns sum of squares
"""
@inline ssq(v::AbstractVector) = var(v) * (length(v) - 1)

"""
    mcpoint(ts::AbstractVector; pcut=0.05, shuffle=500, minlen=3)

return index of change point in mean in `ts` if it exists, 0 otherwise.\\
`pcut` is p-value threshold of existence of a change point, default is 0.05.\\
`shuffle` is number of randomizations to use,\\
`minlen` is minimum allowed length for a segment (prevent short segments).
"""
function mcpoint(ts::AbstractVector; pcut = 0.05, shuffle::Int = 500, minlen::Int = 3)
    minlen < 1 && return 0
    len = length(ts)
    len < 2 * minlen && return 0

    chgpoint = 0
    if mcppv(ts; shuffle = shuffle) <= pcut
        brk = mcptime(ts; minlen = minlen)
        if minlen < brk <= (len - minlen + 1)
            chgpoint = brk
        end
    end

    chgpoint
end

"""
    mcplast(ts::AbstractVector; pcut=0.05, shuffle=500, minlen=3)

return index of rightmost change point in `ts` if any, 0 otherwise.\\
`pcut` is p-value threshold of existence of a change point, default is 0.05.\\
`shuffle` is number of randomizations to use.\\
`minlen` is minimum allowed length for a segment (avoids short segments).
"""
function mcplast(ts::AbstractVector; pcut = 0.05, shuffle::Int = 500, minlen::Int = 3)
    minlen < 1 && return 0
    len = length(ts)
    len < 2 * minlen && return 0

    head = 1
    brk = mcpoint(ts[head:len]; pcut = pcut, minlen = minlen, shuffle = shuffle)
    while brk > 0
        head += brk - 1         # keep checking right hand segment
        brk = mcpoint(ts[head:len]; pcut = pcut, minlen = minlen, shuffle = shuffle)
    end

    head == 1 ? 0 : head
end

##---   all change points if they exist

"""
    mcpall(ts::AbstractVector; pcut = 0.05, shuffle = 500, minlen = 3)

return vector of index of all change points in `ts` via recursive partitioning.\\
`pcut` is p-value threshold of existence of a change point, default is 0.05.\\
`shuffle` is number of randomizations to use,\\
`minlen` is minimum allowed length for a segment (avoid short segments).
"""
function mcpall(ts::AbstractVector; pcut = 0.05, shuffle::Int = 500, minlen::Int = 3)
    chgpts = Vector{Int}(undef, 0)
    minlen < 1 && return chgpts
    len = length(ts)
    len < 2 * minlen && return chgpts

    _mcpall!(chgpts, 1, ts; pcut = pcut, shuffle = shuffle, minlen = minlen)
    sort!(chgpts)
end

function _mcpall!(
    chgpts::AbstractVector,
    head::Int,
    ts::AbstractVector;
    pcut = 0.05,
    shuffle::Int = 500,
    minlen::Int = 3,
)
    len = length(ts)
    chg = mcpoint(ts; pcut = pcut, shuffle = shuffle, minlen = minlen)
    chg == 0 && return                      # no change point, stop
    append!(chgpts, head + chg - 1)         # remember change point
    _mcpall!(chgpts, head, ts[1:chg-1]; pcut = pcut, shuffle = shuffle, minlen = minlen)
    _mcpall!(
        chgpts,
        head + chg - 1,
        ts[chg:len];
        pcut = pcut,
        shuffle = shuffle,
        minlen = minlen,
    )

    nothing
end

##---   mcplot

"""
    mcplot(ts::AbstractVector; chgpts=Int[], palette=:seaborn_bright)

returns plot of `ts` with `chgpts` as change points (default is none),\\
`palette` is the default color scheme to use.

`mcplot(ts)` is a plain plot of the time series,
`mcplot(ts, chgpts=mcpall(ts))` is plot of ts with all change points.
"""
function mcplot(ts::AbstractVector; chgpts = Int[], palette = :seaborn_bright)
    len = length(ts)
    plt = plot(
        size = (600, 400),
        bar_width = 0.9,
        xlims = (0.5, len + 0.5),
        titlefont = (8),
        xtickfontsize = 6,
        palette = palette,
    )

    cp = ndims(chgpts) == 0 ? [chgpts] : vec(chgpts)    # always a vector

    vline!(-0.5 .+ cp, lc = :blue, ls = :dash, lw = 1, label = false)

    n = 0
    for (t1, t2) in zip(reverse([1; cp]), reverse([cp .- 1; length(ts)]))
        n += 1
        m = mean(ts[t1:t2])
        plot!(t1:t2, ts[t1:t2], st = :bar, lc = nothing, label = false, fc = n)
        plot!([t1 - 0.5, t2 + 0.5], [m, m], lc = :blue, ls = :dot, lw = 2, label = false)
    end

    plt
end


text:
using OrdinaryDiffEq, Test, ParameterizedFunctions

d_alembert = @ode_def DAlembert begin
    dx = a - b*x + c*t
end a b c
function d_alembert_analytic(u0,p,t::Number)
    a,b,c = p
    ebt = exp(b*t)
    @. exp(-b*t)*(-a*b + c + ebt*(a*b + c*(b*t - 1)) + b^2 * u0)/(b^2)
end

p = (1., 2., 3.)
u0 = [1.0]
tspan = (0.0,10.0)
prob = ODEProblem(ODEFunction(d_alembert.f,
                              jac = d_alembert.jac,
                              analytic=d_alembert_analytic),
                              u0,tspan,p)

sol = solve(prob,Tsit5(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 1e-7
sol = solve(prob,Rosenbrock23(),abstol=1e-8,reltol=1e-8)
@test sol.errors[:l2] < 1e-7
sol = solve(prob,Rodas4(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 1e-7
sol = solve(prob,Veldd4(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 1e-7
sol = solve(prob,Rodas5(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 1e-7
sol = solve(prob,TRBDF2(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 1e-7
sol = solve(prob,Trapezoid(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 1e-7
sol = solve(prob,KenCarp3(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 1e-7
sol = solve(prob,KenCarp4(),abstol=1e-10,reltol=1e-10)
@test sol.errors[:l2] < 1e-7


text:
# This file contains some real first-rate type piracy. Sorry in advance.

import LinearAlgebra: \

UpperTriangular(D::Diagonal) = D
LowerTriangular(D::Diagonal) = D


text:
const USE_GPU = false
using ParallelStencil
using ParallelStencil.FiniteDifferences2D
@static if USE_GPU
    @init_parallel_stencil(CUDA, Float64, 2)
else
    @init_parallel_stencil(Threads, Float64, 2)
end
using Plots, Printf, ImplicitGlobalGrid
import MPI

macro H_dx(ix, iy)
    esc(:(H[$ix+1, $iy] - H[$ix, $iy]))
end
macro H_dy(ix, iy)
    esc(:(H[$ix, $iy+1] - H[$ix, $iy]))
end
macro Hv_dy(ix, iy)
    esc(:(H[$ix, $iy+1] * v[$ix, $iy+1] - H[$ix, $iy] * v[$ix, $iy]))
end
macro Hv_py(ix, iy)
    esc(:(H[$ix, $iy+1] * v[$ix, $iy+1] + H[$ix, $iy] * v[$ix, $iy]))
end
macro Huv_px(ix, iy)
    esc(:(H[$ix+1, $iy] * u[$ix+1, $iy] * v[$ix+1, $iy] + H[$ix, $iy] * u[$ix, $iy] * v[$ix, $iy]))
end
macro Huv_dx(ix, iy)
    esc(:(H[$ix, $iy+1] * v[$ix, $iy+1] * u[$ix, $iy+1] - H[$ix, $iy] * v[$ix, $iy] * u[$ix, $iy]))
end

"""
	compute_H!(H_new, H, G, u, v, dt_dx, dt_dy, size_H1_2, λx_2, λy_2)

ParallelStencil kernel that computes a time step of the shallow water 2D
process.
"""
@parallel_indices (ix, iy) function compute_H!(
    H_new,
    H,
    u,
    v,
    dt_dx,
    dt_dy,
    size_H1_2,
    λx_2,
    λy_2,
)
    if (ix <= size_H1_2 && iy <= size(H, 2))
        temp = H[ix+1, iy] -
               dt_dx * (
            1 / 2 * (H[ix+2, iy] * u[ix+2, iy] - H[ix, iy] * u[ix, iy]) -
            λx_2 * (@H_dx(ix + 1, iy) - @H_dx(ix, iy))
        )
        if (iy == size(H, 2))
            H_new[ix+1, iy] = (
                temp + dt_dy * (1 / 2 * @Hv_py(ix + 1, iy - 1) - λy_2 * @H_dy(ix + 1, iy - 1))
            )
        elseif (iy == 1)
            H_new[ix+1, iy] = (
                temp - dt_dy * (1 / 2 * @Hv_py(ix + 1, iy) - λy_2 * @H_dy(ix + 1, iy))
            )
        else
            H_new[ix+1, iy] = (
                temp - dt_dy * (1 / 2 * (@Hv_dy(ix + 1, iy) + @Hv_dy(ix + 1, iy - 1)) - λy_2 * (@H_dy(ix + 1, iy) - @H_dy(ix + 1, iy - 1)))
            )
        end
    end
    return
end

"""
	compute_u!(u_new, u, H_new, H, dt_dx, dt_dy, size_H1_2, g)

ParallelStencil kernel that computes the Lax-Friedrichs for the momentum in x-direction of the shallow water
equations.
"""
@parallel_indices (ix, iy) function compute_u!(
    u_new,
    u,
    H_new,
    H,
    dt_dx,
    dt_dy,
    size_H1_2,
    g,
)
    if (ix <= size_H1_2 && iy <= size(H, 2))
        u_new[ix+1, iy] =
            (
                1 / 2 * (H[ix, iy] * u[ix, iy] + H[ix+2, iy] * u[ix+2, iy]) -
                dt_dx * 1 / 2 * (
                    H[ix+2, iy] * u[ix+2, iy]^2 + 0.5 * g * H[ix+2, iy]^2 -
                    H[ix, iy] * u[ix, iy]^2 - 0.5 * g * H[ix, iy]^2
                )
            ) / H_new[ix+1, iy]
    end
    return
end

"""
	compute_v!(v_new, v, H_new, H, u, dt_dx, dt_dy, size_H2_2, g)

ParallelStencil kernel that computes the Lax-Friedrichs for the momentum in y-direction of the shallow water
equations.
"""
@parallel_indices (ix, iy) function compute_v!(
    v_new,
    v,
    H_new,
    H,
    u,
    dt_dx,
    dt_dy,
    size_H2_2,
    g,
)
    if (ix <= size(H, 1) && iy <= size_H2_2)
        temp = 1 / 2 * (H[ix, iy] * v[ix, iy] + H[ix, iy+2] * v[ix, iy+2]) - dt_dy * 1 / 2 * (
            H[ix, iy+2] * v[ix, iy+2]^2 - H[ix, iy] * v[ix, iy]^2 +
            0.5 * g * H[ix, iy+2]^2 - 0.5 * g * H[ix, iy]^2
        )
        if (ix == size(H, 1))
            v_new[ix, iy+1] =
                (
                    temp + dt_dx * 1 / 2 * @Huv_px(ix - 1, iy + 1)
                ) / H_new[ix, iy+1]
        elseif (ix == 1)
            v_new[ix, iy+1] =
                (
                    temp - dt_dx * 1 / 2 * @Huv_px(ix, iy + 1)
                ) / H_new[ix, iy+1]
        else
            v_new[ix, iy+1] =
                (
                    temp - dt_dx * 1 / 2 * (@Huv_dx(ix + 1, iy + 1) + @Huv_dx(ix - 1, iy + 1))
                ) / H_new[ix, iy+1]

        end
    end
    return
end

"""
	bc_x!(A)

Copy boundary conditions in x dimension of array `A`.
"""
@parallel_indices iy function bc_x!(A::Data.Array)
    A[1, iy] = A[2, iy]
    A[end, iy] = A[end-1, iy]
    return
end

"""
	shallow_water_2D_xpu_mpi(; nx, ny, dam1D_x=true, dam2D=false, do_visu=false)

2D shallow water equations solver for an instantaneous dam break.
The Lax-Friedrichs Method was applied to the continuity equation.
Geometry (length of 40 meters, width of 20 meters) and initial conditions
(half of domain has initial water level of 5 meters, other half is dry)
matches BASEMENT version 2.8 Test Case H_1 "Dam break in a closed channel."
Solution of momentum equations requires division by H. Therefore a minimum
H must be set throughout the domain, to avoid numerical instability.

# Arguments
    - `nx`: number of discretised cells for x dimension.
    - `ny`: number of discretised cells for y dimension.
    - `dam1D_x`: if true, 1D dam break in x-direction, else 1D dam break in y-direction. dam2D needs to be false to activate 1D dam break.
    - `dam2D`: if true, 2D dam break
    - `do_visu`: if true, each physical time step will be ploted.

# Return values
    - `H`: The solution arrray (Water surface height in m).
    - `xc`: The x-coord vector.
"""
@views function shallow_water_2D_xpu_mpi(;
    # Numerics
    nx = 512,
    ny = 256,
    # Modeling
    dam1D_x = true,
    dam2D = false,
    # Visualisation
    do_visu = false,
)
    # Physics
    Lx, Ly = 40.0, 20.0
    g = 9.81
    u_max = 5 #from review of results with very small timesteps
    H_init = 5.0
    # Numerics
    nout = 100
    H_min = 1.0
    # Derived numerics
    me, dims = init_global_grid(nx, ny, 1)  # Initialization of MPI and more
    @static if USE_GPU
        select_device()
    end  # select one GPU per MPI local rank (if >1 GPU per node)
    dx, dy = Lx / nx_g(), Ly / ny_g()
    nt = 10000
    xc, yc = LinRange(dx / 2, Lx - dx / 2, nx), LinRange(dy / 2, Ly - dy / 2, ny)
    # Array initialisation
    H = @zeros(nx, ny) .+ H_min
    if (dam1D_x && !dam2D)
        H[1:round(Int64(nx / 2)), :] = @zeros(round(Int64(nx / 2)), ny) .+ H_init #1D dam break in x-direction
    elseif (!dam1D_x && !dam2D)
        H[:, 1:round(Int64(ny / 2))] = zeros(Float64, nx, round(Int64(ny / 2))) .+ H_init #1D dam break in y direction
    else
        H[1:round(Int64(nx/2)), 1:round(Int64(ny/2))] = zeros(Float64, round(Int64(nx/2)), round(Int64(ny/2))) .+ H_init #2D dam break
    end
    H_new = copy(H)
    u = @zeros(nx, ny)
    u_new = copy(u)
    v = @zeros(nx, ny)
    v_new = copy(v)
    size_H1_2, size_H2_2 = size(H, 1) - 2, size(H, 2) - 2

    # Preparation of visualisation
    if do_visu
        if (me == 0)
            ENV["GKSwstype"] = "nul"
            if isdir("plots/part-2/temp") == false
                mkdir("plots/part-2/temp")
            end
            loadpath = "./plots/part-2/temp"
            anim = Animation(loadpath, String[])
            println("Animation directory: $(anim.dir)")
        end
        nx_v, ny_v = (nx - 2) * dims[1], (ny - 2) * dims[2]
        if (nx_v * ny_v * sizeof(Data.Number) > 0.8 * Sys.free_memory())
            error("Not enough memory for visualization.")
        end
        H_v = zeros(nx_v, ny_v) # global array for visu
        H_inn = zeros(nx - 2, ny - 2) # no halo local array for visu
        Xi_g, Yi_g = LinRange(dx + dx / 2, Lx - dx - dx / 2, nx_v), LinRange(dy + dy / 2, Ly - dy - dy / 2, ny_v) # inner points only
    end
    t_tic = 0.0
    niter = 0
    time = 0
    # Time loop
    for it = 0:nt
        if (it == 11)
            t_tic = Base.time()
            niter = 0
        end
        dt = dx / (max(maximum(abs.(u)), maximum(abs.(v)), u_max)) / 4.1
        time = time + dt
        λx = dx / 2 / dt #x-direction signal speed for Lax-Friedrichs scheme
        λy = dy / 2 / dt #y-direction signal speed for Lax-Friedrichs scheme
        dt_dx, dt_dy = dt / dx, dt / dy
        λx_2, λy_2 = λx / 2, λy / 2

        #Lax-Friedrichs for Continuity
        @parallel compute_H!(H_new, H, u, v, dt_dx, dt_dy, size_H1_2, λx_2, λy_2)
        #Lax-Friedrichs for Momentum in x-direction
        @parallel compute_u!(u_new, u, H_new, H, dt_dx, dt_dy, size_H1_2, g)
        #Lax-Friedrichs for Momentum in y-direction
        @parallel compute_v!(v_new, v, H_new, H, u, dt_dx, dt_dy, size_H2_2, g)

        #No-flux boundary conditions
        @parallel 1:size(H, 2) bc_x!(H_new)
        u_new[1, :] .= 0
        u_new[end, :] .= 0
        v_new[:, 1] .= 0
        v_new[:, end] .= 0

        # pointer swap
        H .= H_new
        u .= u_new
        v .= v_new

        update_halo!(H)
        update_halo!(u)
        update_halo!(v)

        niter += 1
        # Visualize
        if do_visu && (it % nout == 0)
            H_inn .= H[2:end-1, 2:end-1]
            gather!(H_inn, H_v)
            z_max = maximum(H_v)
            z_min = minimum(H_v)
            if (me == 0)
                plot(Xi_g, Yi_g, Array(H_v)', st = :surface, xlims = (Xi_g[1], Xi_g[end]), ylims = (Yi_g[1], Yi_g[end]), zlims = (z_min, z_max), c = :blues, colorbar = false, xlabel = "Lx", ylabel = "Ly", zlabel = "water surface elevation (m)", title = "time = $(round(time, sigdigits=3)) s")
                frame(anim)
            end
        end
    end
    t_toc = Base.time() - t_tic
    A_eff = (3 * 2) / 1e9 * nx * ny * sizeof(Float64)  # Effective main memory access per iteration [GB]
    t_it = t_toc / niter                      # Execution time per iteration [s]
    T_eff = A_eff / t_it                       # Effective memory throughput [GB/s]
    if (me == 0)
        @printf(
            "Time = %1.3f sec, T_eff = %1.3f GB/s, Center %s = %1.6f m (niter = %d)\n",
            t_toc,
            round(T_eff, sigdigits = 3),
            (round(Int64(nx / 2)), round(Int64(ny / 2))),
            H[round(Int64(nx / 2)), round(Int64(ny / 2))],
            niter
        )
    end
    if (do_visu && me == 0)
        if (dam1D_x && !dam2D)
            gif(anim, "plots/part-2/shallow_water_2D_xpu_mpi_1D-dam-x.gif", fps = 5)
        elseif (!dam1D_x && !dam2D)
            gif(anim, "plots/part-2/shallow_water_2D_xpu_mpi_1D-dam-y.gif", fps = 5)
        else
            gif(anim, "plots/part-2/shallow_water_2D_xpu_mpi_2D-dam.gif", fps = 15)
        end
    end
    # Retrieve the global domain
    finalize_global_grid()

    return Array(H), xc
end

# shallow_water_2D_xpu_mpi(; nx=256, ny=256, dam2D = false, dam1D_x = false, do_visu = true)


text:
using Plots
using GDAL
using Shapefile
using GBIF
using StatsBase
using CSV
using DataFrames
using Statistics

cd("/Users/kiristern/Documents/GitHub/SDM/")
include("../BioClim/src/lib/SDMLayer.jl")
include("../BioClim/src/lib/gdal.jl")
include("../BioClim/src/lib/worldclim.jl")
include("../BioClim/src/lib/bioclim.jl")
include("../BioClim/src/lib/shapefiles.jl")
include("explo_fnc.jl")

# Import CSV file
sp_table = CSV.read("data/preds/originals/0010800-190621201848488.csv")
## View sp_table in consol: values(sp_table)
# select only scientificName, decimalLatitude, decimalLongitude, and year
sp_table = sp_table[:, [:species, :infraspecificEpithet, :taxonRank,
                        :decimalLatitude, :decimalLongitude, :year]]
## if subspecies, specify in new_sp name
sp_table.new_sp = copy(sp_table.species)
for i in 1:length(sp_table)
    if sp_table.taxonRank[i] == "SUBSPECIES"
        sp_table.new_sp[i] = string(sp_table.species[i], "_", sp_table.infraspecificEpithet[i])
    end
end

# merge two columns
# sp_table.new_sp = sp_table[:species] .* "_" .* sp_table[:infraspecificEpithet]

# Get number of observatons per species/subspecies
# newdf = by(df, :new_sp, n = :new_sp => length)
# show(sort(newdf, order(:new_sp)), allrows=true)
#
# # Show result for subspecies observations
# show(sp_table[sp_table.taxonRank .== "SUBSPECIES", [:species, :infraspecificEpithet, :new_sp]], allrows=true)

# # select only occurences with included year date
# sp_table = sp_table[sp_table.year .!= "NA", :]
# # convert year "string" to "integer"
# sp_table.year = parse.(Int64, sp_table.year)

# Prepare csv data using fnc
df = prepare_csvdata(sp_table)

## From Tim's BioClim/main.jl
# Get the worldclim data by their layer number
cd("/Users/kiristern/Documents/GitHub/BioClim/")

@info "Extract and crop bioclim variables"
@time wc_vars = [clip(worldclim(i), df) for i in 1:19];
# Make the prediction for each layer
@info "Predictions for each layer"
@time predictions = [bioclim(wc_vars[i], df) for i in 1:length(wc_vars)];
# Make the final prediction by taking the minimum
@info "Minimum-consensus aggregation"
@time prediction = reduce(minimum, predictions);
# Get the threshold for NaN given a percentile
@info "Threshold estimation"
@time threshold = first(quantile(prediction[df], [0.05]))
@info "5% threshold:\t$(round(threshold; digits=3))"
# Filter the predictions based on the threshold
@info "Final prediction filtering"
@time for i in eachindex(prediction.grid)
    prediction.grid[i] < threshold && (prediction.grid[i] = NaN)
end

worldmap = clip(worldshape(50), prediction)

sdm_plot = plot([0.0], lab = " ", msw = 0.0, ms = 0.0, size = (900, 450), frame=:box,
                title = first(unique(df.species)))
xaxis!(sdm_plot, (prediction.left, prediction.right), "Longitude")
yaxis!(sdm_plot, (prediction.bottom, prediction.top), "Latitude")

for p in worldmap
    sh = Shape([pp.x for pp in p.points], [pp.y for pp in p.points])
    plot!(sdm_plot, sh, c = :lightgrey, lab = " ")
end

heatmap!(
    sdm_plot,
    longitudes(prediction), latitudes(prediction), prediction.grid,
    aspectratio = 92.60/60.75, c = :BuPu,
    clim = (0.0, maximum(filter(!isnan, prediction.grid)))
    )

for p in worldmap
    xy = map(x -> (x.x, x.y), p.points)
    plot!(sdm_plot, xy, c=:grey, lab = " ", leg=false)
end

sdm_plot

savefig("../SDM/figures/sdm.png")
savefig("../SDM/figures/sdm_map-$(first(unique(sp_table.species))).pdf")

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
#=##############################################################################
# DESCRIPTION
    Testing of an isolated, planar, 45-deg swept-back wing (Bertin's planar wing
    in Example 7.2, pp. 343 of Bertin's Aerodynamics for Engineers). Validated
    with experimental data from Weber and Brebner (1958), Low-speed tests on
    45-deg swept-back wings, part I, Tables 3 and 4.

# AUTHORSHIP
  * Author    : Eduardo J. Alvarez
  * Email     : Edo.AlvarezR@gmail.com
  * Created   : Oct 2019
  * License   : MIT
=###############################################################################

"""
    Test FLOWVLM solver with an isolated, planar, swept wing.
"""
function bertin_VLM(;   # TEST OPTIONS
                        tol=0.025,
                        wake_coupled=true,
                        nsteps=200,
                        vlm_fsgm=-1,
                        surf_fsgm=0.05,
                        # OUTPUT OPTIONS
                        save_path=nothing,
                        run_name="bertins",
                        prompt=true,
                        verbose=true, verbose2=true, v_lvl=1,
                        disp_plot=true, figsize_factor=5/6
                        )

    if verbose; println("\t"^(v_lvl)*"Running Bertin's wing test..."); end;

    # ------------- GENERATE BERTIN'S WING -------------------------------------
    if verbose; println("\t"^(v_lvl+1)*"Generating geometry..."); end;
    # Experimental conditions
    magVinf = 163*0.3048            # (m/s) freestream
    rhoinf = 9.093/10^1             # (kg/m^3) air density
    alpha = 4.2                     # (deg) angle of attack
    qinf = 0.5*rhoinf*magVinf^2     # (Pa) static pressure

    # Geometry
    twist = 0.0                     # (deg) root twist
    lambda = 45.0                   # (deg) sweep
    gamma = 0.0                     # (deg) Dihedral
    b = 98*0.0254                   # (m) span
    ar = 5.0                        # Aspect ratio
    tr = 1.0                        # Taper ratio

    # Discretization
    n = 4*2^4                       # Number of horseshoes
    r = 12.0                        # Geometric expansion
    central = false                 # Central expansion

    # Freestream function
    Vinf(X, t) = magVinf*[cos(alpha*pi/180), 0.0, sin(alpha*pi/180)]

    # Generate wing
    wing = vlm.simpleWing(b, ar, tr, twist, lambda, gamma;
                                                    n=n, r=r, central=central)

    # ------------- SIMULATION SETUP -------------------------------------------
    if verbose; println("\t"^(v_lvl+1)*"Simulation setup..."); end;

    wake_len = 2*b              # (m) length to develop the wake
    lambda_vpm = 2.0            # target core overlap of vpm wake

    # Simulation options
    telapsed = wake_len/magVinf # (s) total time to perform maneuver
    # nsteps = 2000             # Number of time steps
    Vcruise = 0.0               # (m/s) aircraft velocity during cruise (dummy)
    RPMh_w = 0.0                # Rotor RPM during hover (dummy)

    # Solver options
    p_per_step = 1              # Number of particle sheds per time steps (dummy)
    overwrite_sigma = lambda_vpm * magVinf * (telapsed/nsteps)/p_per_step # Smoothing core size
    # vlm_sigma = -1            # VLM regularization core size (deactivated with -1)
    vlm_sigma = vlm_fsgm*b
    surf_sigma = surf_fsgm*b    # Smoothing radius of lifting surface on VPM
    # wake_coupled = true       # Coupled VPM wake with VLM solution
    shed_unsteady = true        # Whether to shed unsteady-loading wake
    # shed_unsteady = false
    vlm_init = true             # Initialize with the VLM semi-infinite wake solution

    # Maneuver definition (dummy)
    Vaircraft(t) = zeros(3)     # Translational velocity of system
    angle_wing(t) = zeros(3)    # Angle of the system

    angle = ()                  # Angle of each tilting system
    RPM = ()                    # RPM of each rotor system
    Vvehicle = Vaircraft        # Velocity of the vehicle
    anglevehicle = angle_wing   # Angle of the vehicle

    maneuver = uns.KinematicManeuver(angle, RPM, Vvehicle, anglevehicle)

    # System definitions
    system = vlm.WingSystem()   # System of all FLOWVLM objects
    vlm.addwing(system, "BertinsWing", wing)

    vlm_system = system         # System solved through VLM solver
    wake_system = system        # System that will shed a VPM wake

    # Vehicle definition
    vehicle = uns.VLMVehicle(   system;
                                vlm_system=vlm_system,
                                wake_system=wake_system
                             )

    if verbose
        println("\t"^(v_lvl+1)*"Core overlap:\t\t$(lambda_vpm)")
        println("\t"^(v_lvl+1)*"Core size:\t\t$(round(overwrite_sigma/b, digits=3))*b")
        println("\t"^(v_lvl+1)*"Time step translation:\t$(round(magVinf * (telapsed/nsteps)/b, digits=3))*b")
    end

    # ------------- SIMULATION MONITOR -----------------------------------------
    y2b = 2*wing._ym/b

    # Weber's lift distribution data (Table 3)
    web_2yb = [0.0, 0.041, 0.082, 0.163, 0.245, 0.367, 0.510, 0.653, 0.898, 0.949]
    web_Cl = [0.235, 0.241, 0.248, 0.253, 0.251, 0.251, 0.251, 0.246, 0.192, 0.171]
    web_CL = 0.238
    web_ClCL = web_Cl/web_CL

    # Weber's drag distribution data (Table 3)
    web_Cd = [0.059, 0.025, 0.016, 0.009, 0.007, 0.006, 0.006, 0.004, -0.002, -0.007]
    web_CD = 0.005
    web_CdCD = web_Cd/web_CD

    function monitor(sim, PFIELD, T, DT; figname="monitor_$(save_path)", nsteps_plot=1)

        aux = PFIELD.nt/nsteps
        clr = (1-aux, 0, aux)

        if PFIELD.nt==0 && disp_plot
            figure(figname, figsize=[7*2, 5*2]*figsize_factor)
            subplot(221)
            xlim([0,1])
            xlabel(L"$\frac{2y}{b}$")
            ylabel(L"$\frac{Cl}{CL}$")
            title("Spanwise lift distribution")

            subplot(222)
            xlim([0,1])
            xlabel(L"$\frac{2y}{b}$")
            ylabel(L"$\frac{Cd}{CD}$")
            title("Spanwise drag distribution")

            subplot(223)
            xlabel("Simulation time (s)")
            ylabel(L"Lift Coefficient $C_L$")

            subplot(224)
            xlabel("Simulation time (s)")
            ylabel(L"Drag Coefficient $C_D$")

            figure(figname*"_2", figsize=[7*2, 5*1]*figsize_factor)
            subplot(121)
            xlabel(L"$\frac{2y}{b}$")
            ylabel(L"Circulation $\Gamma$")
            subplot(122)
            xlabel(L"$\frac{2y}{b}$")
            ylabel(L"Effective velocity $V_\infty$")
        end

        if PFIELD.nt%nsteps_plot==0 && disp_plot
            figure(figname)

            vlm.calculate_field(wing, "Ftot"; S=b^2/ar, qinf=qinf, rhoinf=rhoinf)
            vlm.calculate_field(wing, "CFtot"; S=b^2/ar, qinf=qinf, rhoinf=rhoinf)
            vlm.calculate_field(wing, "Cftot/CFtot"; S=b^2/ar, qinf=qinf, rhoinf=rhoinf)
            ClCL1 = wing.sol["Cl/CL"]
            CdCD1 = wing.sol["Cd/CD"]

            info = vlm.fields_summary(wing)
            CL = info["CL"]
            CD = info["CD"]

            subplot(221)
            plot(web_2yb, web_ClCL, "ok", label="Weber's experimental data")
            plot(y2b, ClCL1, "-", label="FLOWVLM", alpha=0.5, color=clr)

            subplot(222)
            plot(web_2yb, web_CdCD, "ok", label="Weber's experimental data")
            plot(y2b, CdCD1, "-", label="FLOWVLM", alpha=0.5, color=clr)

            subplot(223)
            plot([0, T], web_CL*ones(2), ":k", label="Weber's experimental data")
            plot([T], [CL], "o", label="FLOWVLM", alpha=0.5, color=clr)

            subplot(224)
            plot([0, T], web_CD*ones(2), ":k", label="Weber's experimental data")
            plot([T], [CD], "o", label="FLOWVLM", alpha=0.5, color=clr)

            figure(figname*"_2")
            subplot(121)
            plot(y2b, wing.sol["Gamma"], "-", label="FLOWVLM", alpha=0.5, color=clr)
            if wake_coupled && PFIELD.nt!=0
                subplot(122)
                plot(y2b, norm.(wing.sol["Vkin"]), "-", label="FLOWVLM", alpha=0.5, color=[clr[1], 1, clr[3]])
                plot(y2b, norm.(wing.sol["Vvpm"]), "-", label="FLOWVLM", alpha=0.5, color=clr)
                plot(y2b, [norm(Vinf(vlm.getControlPoint(wing, i), T)) for i in 1:vlm.get_m(wing)]/magVinf,
                                                            "-k", label="FLOWVLM", alpha=0.5)
            end
        end

        return false
    end

    # ------------- RUN SIMULATION ---------------------------------------------
    # Simulation setup
    Vref = Vcruise                  # Reference velocity
    RPMref = RPMh_w                 # Reference RPM
    ttot = telapsed                 # Total time to perform maneuver
    Vinit = Vref*Vaircraft(0)       # Initial vehicle velocity
                                    # Maximum number of particles
    max_particles = ceil(Int, (nsteps+2)*(2*vlm.get_m(vehicle.vlm_system)+1)*p_per_step)

    simulation = uns.Simulation(vehicle, maneuver, Vref, RPMref, ttot;
                                                                    Vinit=Vinit)

    if verbose; println("\t"^(v_lvl+1)*"Running simulation..."); end;
    pfield = uns.run_simulation(simulation, nsteps;
                                      # SIMULATION OPTIONS
                                      Vinf=Vinf,
                                      # SOLVERS OPTIONS
                                      p_per_step=p_per_step,
                                      overwrite_sigma=overwrite_sigma,
                                      vlm_sigma=vlm_sigma,
                                      surf_sigma=surf_sigma,
                                      vlm_init=vlm_init,
                                      max_particles=max_particles,
                                      wake_coupled=wake_coupled,
                                      shed_unsteady=shed_unsteady,
                                      extra_runtime_function=monitor,
                                      # OUTPUT OPTIONS
                                      save_path=save_path,
                                      run_name=run_name,
                                      prompt=prompt,
                                      verbose=verbose2, v_lvl=v_lvl+1,
                                      save_horseshoes=!wake_coupled
                                      )

    # ------------- POST-PROCESSING --------------------------------------------
    if verbose; println("\t"^(v_lvl+1)*"Postprocessing..."); end;
    # Simulation fift and drag
    vlm.calculate_field(wing, "CFtot"; S=b^2/ar, qinf=qinf, rhoinf=rhoinf)
    info = vlm.fields_summary(wing)
    CLsim = info["CL"]
    CDsim = info["CD"]

    # Weber's experimental lift and drag (Table 4)
    CLexp = 0.238
    CDexp = 0.005

    # Error
    CLerr = abs(CLexp-CLsim)/CLexp
    CDerr = abs(CDexp-CDsim)/CDexp

    res = CLerr<tol && (CDerr<tol || true)

    if verbose
        t = "\t"^(v_lvl+1)
        @printf "%0s%10s\t%-11s\t%-11s\t%-11s\n"    t "PARAMETER"   "Experimental"  "Simulation"    "Error %"
        @printf "%0s%10s\t%11.5e\t%11.5e\t%11.5e\n" t "CL"          CLexp           CLsim           100*CLerr
        @printf "%0s%10s\t%11.5e\t%11.5e\t%11.5e\n" t "CD"          CDexp           CDsim           100*CDerr
        println("\t"^(v_lvl+1)*"TEST RESULT:\t$res")
    end

    return res
end

"""
    Test FLOWVLM solver on kinematics of an isolated, planar, swept wing.
"""
function bertin_kinematic(;   # TEST OPTIONS
                        tol=0.025,
                        wake_coupled=true,
                        nsteps=150,
                        vlm_fsgm=-1,
                        surf_fsgm=0.046,
                        p_per_step = 1,
                        vlm_rlx = -1,
                        VehicleType=uns.VLMVehicle,
                        # OUTPUT OPTIONS
                        save_path=nothing,
                        run_name="bertins",
                        prompt=true,
                        verbose=true, verbose2=true, v_lvl=1,
                        disp_plot=true, figsize_factor=5/6
                        )

    if verbose; println("\t"^(v_lvl)*"Running Bertin's wing test..."); end;

    # ------------- GENERATE BERTIN'S WING -------------------------------------
    if verbose; println("\t"^(v_lvl+1)*"Generating geometry..."); end;
    # Experimental conditions
    magVinf = 163*0.3048            # (m/s) freestream
    rhoinf = 9.093/10^1             # (kg/m^3) air density
    alpha = 4.2                     # (deg) angle of attack
    qinf = 0.5*rhoinf*magVinf^2     # (Pa) static pressure

    # Geometry
    twist = 0.0                     # (deg) root twist
    lambda = 45.0                   # (deg) sweep
    gamma = 0.0                     # (deg) Dihedral
    b = 98*0.0254                   # (m) span
    ar = 5.0                        # Aspect ratio
    tr = 1.0                        # Taper ratio

    # Discretization
    n = 4*2^4                       # Number of horseshoes
    r = 12.0                        # Geometric expansion
    central = false                 # Central expansion

    # Freestream function
    # Vinf(X, t) = magVinf*[cos(alpha*pi/180), 0.0, sin(alpha*pi/180)]
    # Vinf(X, t) = 1e-12*ones(3)      # (Don't make this zero or things will break)
    # Here I had to give it an initial freestream or the unsteady shedding would
    # in the first step
    Vinf(X, t) = t==0 ? magVinf*[1,0,0] : 1e-12*ones(3)

    # Generate wing
    wing = vlm.simpleWing(b, ar, tr, twist, lambda, gamma;
                                                    n=n, r=r, central=central)

    # Pitch wing to corresponding angle of attack
    O = zeros(3)                    # Coordinate system origin
    Oaxis = gt.rotation_matrix2(0.0, -alpha, 0.0) # Coordinate system axes
    vlm.setcoordsystem(wing, O, Oaxis)

    # ------------- SIMULATION SETUP -------------------------------------------
    if verbose; println("\t"^(v_lvl+1)*"Simulation setup..."); end;

    wake_len = 2*b              # (m) length to develop the wake
    lambda_vpm = 2.0            # target core overlap of vpm wake

    # Simulation options
    telapsed = wake_len/magVinf # (s) total time to perform maneuver
    # nsteps = 2000             # Number of time steps
    Vcruise = magVinf           # (m/s) aircraft velocity during cruise
    RPMh_w = 0.0                # Rotor RPM during hover (dummy)

    # Solver options
    # p_per_step = 1              # Number of particle sheds per time steps
    overwrite_sigma = lambda_vpm * magVinf * (telapsed/nsteps)/p_per_step # Smoothing core size
    # vlm_sigma = -1            # VLM regularization core size (deactivated with -1)
    vlm_sigma = vlm_fsgm*b
    surf_sigma = surf_fsgm*b    # Smoothing radius of lifting surface on VPM
    # wake_coupled = true       # Coupled VPM wake with VLM solution
    shed_unsteady = true        # Whether to shed unsteady-loading wake
    # shed_unsteady = false
    # vlm_rlx = -1                # VLM relaxation (deactivated with -1)
    vlm_init = true             # Initialize with the VLM semi-infinite wake solution

    # Maneuver definition
    Vaircraft(t) = [-1,0,0]     # Translational velocity of system
    angle_wing(t) = zeros(3)    # Angle of the system

    angle = ()                  # Angle of each tilting system
    RPM = ()                    # RPM of each rotor system
    Vvehicle = Vaircraft        # Velocity of the vehicle
    anglevehicle = angle_wing   # Angle of the vehicle

    maneuver = uns.KinematicManeuver(angle, RPM, Vvehicle, anglevehicle)

    # System definitions
    system = vlm.WingSystem()   # System of all FLOWVLM objects
    vlm.addwing(system, "BertinsWing", wing)

    vlm_system = system         # System solved through VLM solver
    wake_system = system        # System that will shed a VPM wake

    # Vehicle definition
    vehicle = VehicleType(      system;
                                vlm_system=vlm_system,
                                wake_system=wake_system
                             )

    if verbose
        println("\t"^(v_lvl+1)*"Core overlap:\t\t$(lambda_vpm)")
        println("\t"^(v_lvl+1)*"Core size:\t\t$(round(overwrite_sigma/b, digits=3))*b")
        println("\t"^(v_lvl+1)*"Time step translation:\t$(round(magVinf * (telapsed/nsteps)/b, digits=3))*b")
    end

    # ------------- SIMULATION MONITOR -----------------------------------------
    y2b = 2*wing._ym/b

    # Weber's lift distribution data (Table 3)
    web_2yb = [0.0, 0.041, 0.082, 0.163, 0.245, 0.367, 0.510, 0.653, 0.898, 0.949]
    web_Cl = [0.235, 0.241, 0.248, 0.253, 0.251, 0.251, 0.251, 0.246, 0.192, 0.171]
    web_CL = 0.238
    web_ClCL = web_Cl/web_CL

    # Weber's drag distribution data (Table 3)
    web_Cd = [0.059, 0.025, 0.016, 0.009, 0.007, 0.006, 0.006, 0.004, -0.002, -0.007]
    web_CD = 0.005
    web_CdCD = web_Cd/web_CD

    prev_wing = nothing

    function monitor(sim, PFIELD, T, DT; figname="monitor_$(save_path)", nsteps_plot=1)

        aux = PFIELD.nt/nsteps
        clr = (1-aux, 0, aux)

        if PFIELD.nt==0 && disp_plot
            figure(figname, figsize=[7*2, 5*2]*figsize_factor)
            subplot(221)
            xlim([0,1])
            xlabel(L"$\frac{2y}{b}$")
            ylabel(L"$\frac{Cl}{CL}$")
            title("Spanwise lift distribution")

            subplot(222)
            xlim([0,1])
            xlabel(L"$\frac{2y}{b}$")
            ylabel(L"$\frac{Cd}{CD}$")
            title("Spanwise drag distribution")

            subplot(223)
            xlabel("Simulation time (s)")
            ylabel(L"Lift Coefficient $C_L$")

            subplot(224)
            xlabel("Simulation time (s)")
            ylabel(L"Drag Coefficient $C_D$")

            figure(figname*"_2", figsize=[7*2, 5*1]*figsize_factor)
            subplot(121)
            xlabel(L"$\frac{2y}{b}$")
            ylabel(L"Circulation $\Gamma$")
            subplot(122)
            xlabel(L"$\frac{2y}{b}$")
            ylabel(L"Effective velocity $V_\infty$")
        end

        if PFIELD.nt!=0 && PFIELD.nt%nsteps_plot==0 && disp_plot
            figure(figname)

            # Force at each VLM element
            Ftot = uns.calc_aerodynamicforce(wing, prev_wing, PFIELD, Vinf, DT,
                                                            rhoinf; t=PFIELD.t)
            L, D, S = uns.decompose(Ftot, [0,0,1], [-1,0,0])
            vlm._addsolution(wing, "L", L)
            vlm._addsolution(wing, "D", D)
            vlm._addsolution(wing, "S", S)

            # Force per unit span at each VLM element
            ftot = uns.calc_aerodynamicforce(wing, prev_wing, PFIELD, Vinf, DT,
                                        rhoinf; t=PFIELD.t, per_unit_span=true)
            l, d, s = uns.decompose(ftot, [0,0,1], [-1,0,0])

            # Lift of the wing
            Lwing = norm(sum(L))
            CLwing = Lwing/(qinf*b^2/ar)
            ClCL = norm.(l) / (Lwing/b)

            # Drag of the wing
            Dwing = norm(sum(D))
            CDwing = Dwing/(qinf*b^2/ar)
            CdCD = [sign(dot(this_d, [1,0,0])) for this_d in d].*norm.(d) / (Dwing/b) # Preserves the sign of drag

            vlm._addsolution(wing, "Cl/CL", ClCL)
            vlm._addsolution(wing, "Cd/CD", CdCD)

            subplot(221)
            plot(web_2yb, web_ClCL, "ok", label="Weber's experimental data")
            plot(y2b, ClCL, "-", label="FLOWVLM", alpha=0.5, color=clr)

            subplot(222)
            plot(web_2yb, web_CdCD, "ok", label="Weber's experimental data")
            plot(y2b, CdCD, "-", label="FLOWVLM", alpha=0.5, color=clr)

            subplot(223)
            plot([0, T], web_CL*ones(2), ":k", label="Weber's experimental data")
            plot([T], [CLwing], "o", label="FLOWVLM", alpha=0.5, color=clr)

            subplot(224)
            plot([0, T], web_CD*ones(2), ":k", label="Weber's experimental data")
            plot([T], [CDwing], "o", label="FLOWVLM", alpha=0.5, color=clr)

            figure(figname*"_2")
            subplot(121)
            plot(y2b, wing.sol["Gamma"], "-", label="FLOWVLM", alpha=0.5, color=clr)
            if wake_coupled && PFIELD.nt!=0
                subplot(122)
                plot(y2b, norm.(wing.sol["Vkin"])/magVinf, "-", label="FLOWVLM", alpha=0.5, color=[clr[1], 1, clr[3]])
                if VehicleType==uns.VLMVehicle
                    plot(y2b, norm.(wing.sol["Vvpm"]), "-", label="FLOWVLM", alpha=0.5, color=clr)
                end
                plot(y2b, [norm(Vinf(vlm.getControlPoint(wing, i), T)) for i in 1:vlm.get_m(wing)],
                                                            "-k", label="FLOWVLM", alpha=0.5)
            end
        end

        prev_wing = deepcopy(wing)

        return false
    end

    # ------------- RUN SIMULATION ---------------------------------------------
    # Simulation setup
    Vref = Vcruise                  # Reference velocity
    RPMref = RPMh_w                 # Reference RPM
    ttot = telapsed                 # Total time to perform maneuver
    Vinit = Vref*Vaircraft(0)       # Initial vehicle velocity
                                    # Maximum number of particles
    max_particles = ceil(Int, (nsteps+2)*(2*vlm.get_m(vehicle.vlm_system)+1)*p_per_step)

    simulation = uns.Simulation(vehicle, maneuver, Vref, RPMref, ttot;
                                                                    Vinit=Vinit)

    if verbose; println("\t"^(v_lvl+1)*"Running simulation..."); end;
    pfield = uns.run_simulation(simulation, nsteps;
                                      # SIMULATION OPTIONS
                                      Vinf=Vinf,
                                      # SOLVERS OPTIONS
                                      p_per_step=p_per_step,
                                      overwrite_sigma=overwrite_sigma,
                                      vlm_sigma=vlm_sigma,
                                      surf_sigma=surf_sigma,
                                      vlm_rlx=vlm_rlx,
                                      vlm_init=vlm_init,
                                      max_particles=max_particles,
                                      wake_coupled=wake_coupled,
                                      shed_unsteady=shed_unsteady,
                                      extra_runtime_function=monitor,
                                      # OUTPUT OPTIONS
                                      save_path=save_path,
                                      run_name=run_name,
                                      prompt=prompt,
                                      verbose=verbose2, v_lvl=v_lvl+1,
                                      save_horseshoes=!wake_coupled
                                      )

    # ------------- POST-PROCESSING --------------------------------------------
    if verbose; println("\t"^(v_lvl+1)*"Postprocessing..."); end;
    # Simulation fift and drag
    vlm.calculate_field(wing, "CFtot"; S=b^2/ar, qinf=qinf, rhoinf=rhoinf)
    info = vlm.fields_summary(wing)
    CLsim = info["CL"]
    CDsim = info["CD"]

    # Weber's experimental lift and drag (Table 4)
    CLexp = 0.238
    CDexp = 0.005

    # Error
    CLerr = abs(CLexp-CLsim)/CLexp
    CDerr = abs(CDexp-CDsim)/CDexp

    res = CLerr<tol && (CDerr<tol || true)

    if verbose
        t = "\t"^(v_lvl+1)
        @printf "%0s%10s\t%-11s\t%-11s\t%-11s\n"    t "PARAMETER"   "Experimental"  "Simulation"    "Error %"
        @printf "%0s%10s\t%11.5e\t%11.5e\t%11.5e\n" t "CL"          CLexp           CLsim           100*CLerr
        @printf "%0s%10s\t%11.5e\t%11.5e\t%11.5e\n" t "CD"          CDexp           CDsim           100*CDerr
        println("\t"^(v_lvl+1)*"TEST RESULT:\t$res")
    end

    return res
end

max_stars_repo_name: dendisuhubdy/FLOWUnsteady

max_issues_repo_name: dendisuhubdy/FLOWUnsteady

max_forks_repo_name: dendisuhubdy/FLOWUnsteady


text:
using Base.Test
import WaveletScattering: PointwiseLayer

modulus = PointwiseLayer(bottoms = [:bottom], tops = [:top])

@test !Mocha.can_do_bp(modulus)
@test !Mocha.has_neuron(modulus)
@test !Mocha.has_param(modulus)


text:
"""
Julia program to find the Delannoy number for a given rectangular grid.
A Delannoy number describes the number of paths from the southwest corner
(0, 0) of a rectangular grid to the northeast corner (m, n) using only single steps north, northeast, or east.
Delannoy number is named after the French mathematican Henri Delannoy.
"""

function delannoy_number(r, c)

    # Initialize the dp array with '0' as value.
    dp = zeros(Int, c + 1, r + 1)

    for i in 1:c
        dp[1,i] = 1
    end
    for i in 2:c+1
        dp[i,1] = 1
    end
    # From Each point calculate the number of paths, to the north-east point, that can be reached by
    # traversing through the immediate right or immediate top or immediate top-right point.
    for i in 2:c+1
        for j in 2:r+1
            dp[i,j] = dp[i-1,j] + dp[i-1,j-1] + dp[i,j-1]
        end
    end
    return dp[c+1,r+1]
end

print("Enter the co-ordinates of the north-east corner.\nEnter the row co-ordinate: ")
x = readline()
x = parse(Int, x)
print("Enter the column co-ordinate: ")
y = readline()
y = parse(Int, y)
res = delannoy_number(x, y)
print("The Delannoy number of the given grid is $res.")

"""
Time Complexity - O(x * y), where `x`, `y` is the given co-ordinates.
Space Complexity - O(1)

SAMPLE INPUT AND OUTPUT

SAMPLE I

Enter the co-ordinates of the north-east corner.
Enter the row co-ordinate: 5
Enter the column co-ordinate: 6
The Delannoy number of the given grid is 3653.
"""

max_issues_repo_name: adarshnjena/NeoAlgo

max_forks_repo_name: adarshnjena/NeoAlgo


text:
# # Maximum Likelihood Estimation: The Normal Linear Model
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`maxlikenlm.ipynb`](@__NBVIEWER_ROOT_URL__examples/generated/maxlikenlm.ipynb)
#-
#
# The following tutorial will introduce maximum likelihood estimation
# in Julia for the normal linear model.
#
# The normal linear model (sometimes referred to as the OLS model) is
# the workhorse of regression modeling and is utilized across a number
# of diverse fields.  In this tutorial, we will utilize simulated data
# to demonstrate how Julia can be used to recover the parameters of
# interest.
#
# The first order of business is to use the `Optim` package
# and also include the `NLSolversBase` routine:
#

using Optim, NLSolversBase, Random
using LinearAlgebra: diag
Random.seed!(0);                            # Fix random seed generator for reproducibility

#md # !!! tip
#md #     Add Optim with the following command at the Julia command prompt:
#md #     `Pkg.add("Optim")`
#
# The first item that needs to be addressed is the data generating process or DGP.
# The following code will produce data from a normal linear model:

n = 500                             # Number of observations
nvar = 2                            # Number of variables
β = ones(nvar) * 3.0                # True coefficients
x = [ones(n) randn(n, nvar - 1)]    # X matrix of explanatory variables plus constant
ε = randn(n) * 0.5                  # Error variance
y = x * β + ε;                      # Generate Data

# In the above example, we have 500 observations, 2 explanatory
# variables plus an intercept, an error variance equal to 0.5,
# coefficients equal to 3.0, and all of these are subject to change by
# the user. Since we know the true value of these parameters, we
# should obtain these values when we maximize the likelihood function.
#
# The next step in our tutorial is to define a Julia function for the
# likelihood function. The following function defines the likelihood
# function for the normal linear model:

function Log_Likelihood(X, Y, β, log_σ)
    σ = exp(log_σ)
    llike = -n/2*log(2π) - n/2* log(σ^2) - (sum((Y - X * β).^2) / (2σ^2))
    llike = -llike
end

# The log likelihood function accepts 4 inputs: the matrix of
# explanatory variables (X), the dependent variable (Y), the β's, and
# the error varicance. Note that we exponentiate the error variance in
# the second line of the code because the error variance cannot be
# negative and we want to avoid this situation when maximizing the
# likelihood.
#
# The next step in our tutorial is to optimize our function. We first
# use the `TwiceDifferentiable` command in order to obtain the Hessian
# matrix later on, which will be used to help form t-statistics:

func = TwiceDifferentiable(vars -> Log_Likelihood(x, y, vars[1:nvar], vars[nvar + 1]),
                           ones(nvar+1); autodiff=:forward);

# The above statment accepts 4 inputs: the x matrix, the dependent
# variable y, and a vector of β's and the error variance.  The
# `vars[1:nvar]` is how we pass the vector of β's and the `vars[nvar +
# 1]` is how we pass the error variance. You can think of this as a
# vector of parameters with the first 2 being β's and the last one is
# the error variance.
#
# The `ones(nvar+1)` are the starting values for the parameters and
# the `autodiff=:forward` command performs forward mode automatic
# differentiation.
#
# The actual optimization of the likelihood function is accomplished
# with the following command:

opt = optimize(func, ones(nvar+1))

## Test the results                #src
using Test                    #src
@test Optim.converged(opt)         #src
@test Optim.g_residual(opt) < 1e-8 #src

# The first input to the command is the function we wish to optimize
# and the second input are the starting values.
#
# After a brief period of time, you should see output of the
# optimization routine, with the parameter  estimates being very close
# to our simulated values.
#
# The optimization routine stores several quantities and we can obtain
# the maximim likelihood estimates with the following command:

parameters = Optim.minimizer(opt)
@test parameters ≈ [3.00279, 2.96455, -0.648693] atol=1e-5 #src

# !!! Note
#     Fieldnames for all of the quantities can be obtained with the following command:
#     fieldnames(opt)
#
# Since we paramaterized our likelihood to use the exponentiated
# value, we need to exponentiate it to get back to our original log
# scale:

parameters[nvar+1] = exp(parameters[nvar+1])

# In order to obtain the correct Hessian matrix, we have to "push" the
# actual parameter values that maximizes the likelihood function since
# the `TwiceDifferentiable` command uses the next to last values to
# calculate the Hessian:

numerical_hessian = hessian!(func,parameters)

# We can now invert our Hessian matrix to obtain the variance-covariance matrix:

var_cov_matrix = inv(numerical_hessian)

# In this example, we are only interested in the statistical
# significance of the coefficient estimates so we obtain those with
# the following command:

β = parameters[1:nvar]
@test β ≈ [3.00279, 2.96455] atol=1e-5 #src

# We now need to obtain those elements of the variance-covariance
# matrix needed to obtain our t-statistics, and we can do this with
# the following commands:

temp = diag(var_cov_matrix)
temp1 = temp[1:nvar]

# The t-statistics are formed by dividing element-by-element the
# coefficients by their standard errors, or the square root of the
# diagonal elements of the variance-covariance matrix:

t_stats = β./sqrt.(temp1)
@test t_stats ≈ [39.7191, 39.9506] atol=1e-4 #src

# From here, one may examine other statistics of interest using the
# output from the optimization routine.

#md # ## [Plain Program](@id maxlikenlm-plain-program)
#md #
#md # Below follows a version of the program without any comments.
#md # The file is also available here: [maxlikenlm.jl](maxlikenlm.jl)
#md #
#md # ```julia
#md # @__CODE__
#md # ```