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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ 3fea381c-8b4f-11eb-0db0-9701a3447ddb
using Distributions, Plots, StatsPlots, Survival, Statistics, PlutoUI

# ╔═╡ f8988e20-8ce1-11eb-096d-45c3fb657e83
md" # RWRI#15 Technical Summary"

# ╔═╡ e026e85e-8ce2-11eb-24ef-4dd84cb657be
html"<button onclick='present()'>present</button>"

# ╔═╡ d4a856ea-8b55-11eb-1ab3-7347511e1d5b
import Pkg

# ╔═╡ d89bd6fa-8b55-11eb-000e-cd190c432d25
Pkg.add("PlutoUI")

# ╔═╡ 300d7882-8c17-11eb-27f2-8df170245c92
gr()

# ╔═╡ 38b63284-8b50-11eb-0100-3d640819af62
K = 1_000

# ╔═╡ 985a22de-8ce2-11eb-268d-e13e945f4aa6
md"## Central Limit"

# ╔═╡ 36030816-8b51-11eb-3bbe-99f9c372a692
plot(
	histogram(rand(Uniform(), K), label="single uniform"),
	histogram(( rand(Uniform(), K) + rand(Uniform(), K) ) / 2, label="2 avg uniform"),
	histogram(( rand(Uniform(), K) + rand(Uniform(), K) + rand(Uniform(), K) ) / 3, label="3 avg uniform"),
	histogram(( rand(Uniform(), K) + rand(Uniform(), K) + rand(Uniform(), K) + rand(Uniform(), K)  ) / 4, label="4 avg uniform"),
	layout=(2,2)
)

# ╔═╡ 7a718cd2-8b50-11eb-101b-8178142e592f
plot(
	plot(rand(Normal(), K), title="Behavioral Differences", ylims=(-100, 100), label="gaussian"),
	plot(rand(TDist(1.1), K), label="power law", ylims=(-100, 100)),
	layout=(2,1)
)

# ╔═╡ adbc0b80-8b50-11eb-023c-f3e60c2888cc
begin
	plot(
		histogram(rand(Normal(), K), title="Normal"),
		histogram(rand(TDist(1.1), K), title="Power Law"),
		layout = (2, 1), legend = false
	)
end

# ╔═╡ ee16db52-8b56-11eb-0d84-47c12457033d
begin
	pl = TDist(1.1); x_pl = rand(pl, K)
	normal = Normal(); x_normal = rand(normal, K)
	plot(
		plot(x_normal, ccdf(normal, x_normal),
			seriestype=:scatter, xlims=(0.1,100), xaxis=:log, yaxis=:log),
		plot(x_pl, ccdf(pl, x_pl),
			seriestype=:scatter, xlims=(0.1,100), xaxis=:log, yaxis=:log),
		layout=(1,2)
	)
end

# ╔═╡ a46a5422-8ce2-11eb-2f13-db6dc08f691b
md"## Law of Large Numbers"

# ╔═╡ a958d146-8bf4-11eb-3fb5-59fdcaebd563
plot([cor(rand(Normal(), n), rand(Normal(), n)) for n=1:500],
	 ylims=(-1,1), seriestype=:line,
	 ylabel="correlation", xlabel="sample size"
)

# ╔═╡ 66638b3e-8be7-11eb-1f72-919dd4c2d00d
function corr(n)
	return cor(rand(Normal(), n), rand(Normal(), n))
end

# ╔═╡ 7db4ec02-8bf2-11eb-302c-83a265e29552
@bind samplesize Slider(10:1000;show_value=true)

# ╔═╡ b3714786-8be7-11eb-3550-9968b2c626df
begin
	n = 1000;
	histogram(map(1:n) do i corr(samplesize) end, bins=25, xlims=(-0.75,0.75),
			  xlabel="correlation", ylabel="frequency", label="")
end

# ╔═╡ 26472e82-8bf5-11eb-1ae5-63a5546fbd5c
plot(
	plot([mean(rand(Normal(0,1), n)) for n=1:10^3], seriestype=:line, title="Law of Large Numbers, NOT! Mean",label="Normal", ylims=(-5,5), ylabel="mean"),
	plot([mean(rand(TDist(1.1), n)) for n=1:10^3], seriestype=:line, label="PowerLaw", ylims=(-50,50), xlabel="sample size", ylabel="mean"),
	layout=(2,1)
)

# ╔═╡ d95ac7ce-8bf6-11eb-019d-81056d6ba663
plot(
	plot([var(rand(Normal(0,1), n)) for n=1:10^3], seriestype=:line, title="Law of Large Numbers, NOT! Variance",label="Normal", yaxis=:log,
		 ylabel="log(variance)"),
	plot([var(rand(TDist(1.1), n)) for n=1:10^3], seriestype=:line, yaxis=:log, label="PowerLaw", xlabel="sample size", ylabel="log(variance)"),
	layout=(2,1)
)

# ╔═╡ Cell order:
# ╠═f8988e20-8ce1-11eb-096d-45c3fb657e83
# ╠═e026e85e-8ce2-11eb-24ef-4dd84cb657be
# ╠═d4a856ea-8b55-11eb-1ab3-7347511e1d5b
# ╠═d89bd6fa-8b55-11eb-000e-cd190c432d25
# ╠═3fea381c-8b4f-11eb-0db0-9701a3447ddb
# ╠═300d7882-8c17-11eb-27f2-8df170245c92
# ╠═38b63284-8b50-11eb-0100-3d640819af62
# ╟─985a22de-8ce2-11eb-268d-e13e945f4aa6
# ╠═36030816-8b51-11eb-3bbe-99f9c372a692
# ╠═7a718cd2-8b50-11eb-101b-8178142e592f
# ╠═adbc0b80-8b50-11eb-023c-f3e60c2888cc
# ╠═ee16db52-8b56-11eb-0d84-47c12457033d
# ╟─a46a5422-8ce2-11eb-2f13-db6dc08f691b
# ╠═a958d146-8bf4-11eb-3fb5-59fdcaebd563
# ╠═66638b3e-8be7-11eb-1f72-919dd4c2d00d
# ╠═7db4ec02-8bf2-11eb-302c-83a265e29552
# ╠═b3714786-8be7-11eb-3550-9968b2c626df
# ╠═26472e82-8bf5-11eb-1ae5-63a5546fbd5c
# ╠═d95ac7ce-8bf6-11eb-019d-81056d6ba663


text:
function initialize_without_derivatives(u0, f, p, t0, order, var=1e-3)
    q = order
    d = length(u0)

    m0 = zeros(d*(q+1))
    m0[1:d] = u0
    if !isinplace(f)
        m0[d+1:2d] = f(u0, p, t0)
    else
        f(m0[d+1:2d], u0, p, t0)
    end
    P0 = [zeros(d, d) zeros(d, d*q);
          zeros(d*q, d) diagm(0 => var .* ones(d*q))]
    return m0, P0
end

function __initialize_with_derivatives_forwarddiff(u0, f, p, t0, order::Int)
    error("Legacy code! ForwardDiff is not included in the dependencies anymore. Use the TaylorSeries.jl approach instead, implemented in `initialize_with_derivatives`.")
    @warn "Better to use the TaylorSeries.jl approach"
    f = isinplace(f) ? iip_to_oop(f) : f

    d = length(u0)
    q = order

    uElType = eltype(u0)
    m0 = zeros(uElType, d*(q+1))
    P0 = zeros(uElType, d*(q+1), d*(q+1))

    m0[1:d] .= u0
    m0[d+1:2d] .= f(u0, p, t0)

    f_derivatives = Function[f]
    for o in 2:q
        _curr_f_deriv = f_derivatives[end]
        dfdu(u, p, t) = ForwardDiff.jacobian(u -> _curr_f_deriv(u, p, t), u)
        dfdt(u, p, t) = ForwardDiff.derivative(t -> _curr_f_deriv(u, p, t), t)
        df(u, p, t) = dfdu(u, p, t) * f(u, p, t) + dfdt(u, p, t)
        push!(f_derivatives, df)
        m0[o*d+1:(o+1)*d] = df(u0, p, t0)
    end

    return m0, P0
end

function initialize_with_derivatives(u0, f, p, t0, order::Int)
    f = isinplace(f) ? iip_to_oop(f) : f

    d = length(u0)
    q = order

    set_variables("u", numvars=d, order=order+1)

    uElType = eltype(u0)
    m0 = zeros(uElType, d*(q+1))
    P0 = zeros(uElType, d*(q+1), d*(q+1))

    m0[1:d] .= u0
    m0[d+1:2d] .= f(u0, p, t0)

    # Make sure that the vector field f does not depend on t
    f_t_taylor = taylor_expand(t -> f(u0, p, t), t0)
    @assert !(eltype(f_t_taylor) <: TaylorN)

    fp = taylor_expand(u -> f(u, p, t0), u0)
    f_derivatives = [fp]
    for o in 2:q
        _curr_f_deriv = f_derivatives[end]
        dfdu = stack([derivative.(_curr_f_deriv, i) for i in 1:d])'
        # dfdt(u, p, t) = ForwardDiff.derivative(t -> _curr_f_deriv(u, p, t), t)
        # df(u, p, t) = dfdu(u, p, t) * f(u, p, t) + dfdt(u, p, t)
        df = dfdu * fp
        push!(f_derivatives, df)
        m0[o*d+1:(o+1)*d] = evaluate(df)
    end

    return m0, P0
end

function initialize_with_derivatives(u0, f::DynamicalODEFunction, p, t0, order::Int)
    f = isinplace(f) ? iip_to_oop(f) : f
    function _f(u, p, t)
        _u = ArrayPartition(u[1:1], u[2:2])
        _du = f(_u, p, t)
        du = vcat(_du...)
        return du
    end

    _u0 = vcat(u0...)
    d = length(u0)
    q = order

    set_variables("u", numvars=d, order=order+1)

    uElType = eltype(u0)
    m0 = zeros(uElType, d*(q+1))
    P0 = zeros(uElType, d*(q+1), d*(q+1))

    m0[1:d] .= u0
    _du0 = _f(_u0, p, t0)
    m0[d+1:2d] .= _f(_u0, p, t0)

    # Make sure that the vector field f does not depend on t
    f_t_taylor = taylor_expand(t -> _f(_u0, p, t), t0)
    @assert !(eltype(f_t_taylor) <: TaylorN)

    fp = taylor_expand(u -> _f(u, p, t0), _u0)
    f_derivatives = [fp]
    for o in 2:q
        _curr_f_deriv = f_derivatives[end]
        dfdu = stack([derivative.(_curr_f_deriv, i) for i in 1:d])'
        # dfdt(u, p, t) = ForwardDiff.derivative(t -> _curr_f_deriv(u, p, t), t)
        # df(u, p, t) = dfdu(u, p, t) * f(u, p, t) + dfdt(u, p, t)
        df = dfdu * fp
        push!(f_derivatives, df)
        m0[o*d+1:(o+1)*d] = evaluate(df)
    end

    return m0, P0
end

function iip_to_oop(f!)
    function f(u, p, t)
        du = copy(u)
        f!(du, u, p, t)
        return du
    end
    return f
end

ext: jl


text:
using Documenter
using ModelicaScriptingTools

makedocs(
    sitename = "HH-Modelica",
    format = Documenter.HTML(),
    modules = Module[]
)

deploydocs(
    repo = "github.com/CSchoel/hh-modelica.git",
    devbranch = "main",
    versions = ["v^", "v#.#", "stable" => "v^"]
)

max_stars_repo_name: CSchoel/hh-modelica

max_issues_repo_name: CSchoel/hh-modelica

max_forks_repo_name: CSchoel/hh-modelica


text:
using ParameterizedFunctions, OrdinaryDiffEq, DiffEqParamEstim
using BlackBoxOptim, NLopt, Plots, QuadDIRECT
gr(fmt=:png)

Xiang2015Bounds = Tuple{Float64, Float64}[(9, 11), (20, 30), (2, 3)] # for local optimizations
xlow_bounds = [9.0,20.0,2.0]
xhigh_bounds = [11.0,30.0,3.0]
LooserBounds = Tuple{Float64, Float64}[(0, 22), (0, 60), (0, 6)] # for global optimization
GloIniPar = [0.0, 0.5, 0.1] # for global optimizations
LocIniPar = [9.0, 20.0, 2.0] # for local optimization

g1 = @ode_def LorenzExample begin
  dx = σ*(y-x)
  dy = x*(ρ-z) - y
  dz = x*y - β*z
end σ ρ β
p = [10.0,28.0,2.66] # Parameters used to construct the dataset
r0 = [1.0; 0.0; 0.0]                #[-11.8,-5.1,37.5] PODES Initial values of the system in space # [0.1, 0.0, 0.0]
tspan = (0.0, 30.0)                 # PODES sample of 3000 observations over the (0,30) timespan
prob = ODEProblem(g1, r0, tspan,p)
tspan2 = (0.0, 3.0)                 # Xiang test sample of 300 observations with a timestep of 0.01
prob_short = ODEProblem(g1, r0, tspan2,p)

dt = 30.0/3000
tf = 30.0
tinterval = 0:dt:tf
t  = collect(tinterval)

h = 0.01
M = 300
tstart = 0.0
tstop = tstart + M * h
tinterval_short = 0:h:tstop
t_short = collect(tinterval_short)

# Generate Data
data_sol_short = solve(prob_short,Vern9(),saveat=t_short,reltol=1e-9,abstol=1e-9)
data_short = convert(Array, data_sol_short) # This operation produces column major dataset obs as columns, equations as rows
data_sol = solve(prob,Vern9(),saveat=t,reltol=1e-9,abstol=1e-9)
data = convert(Array, data_sol)

plot(data_sol_short,vars=(1,2,3)) # the short solution
plot(data_sol,vars=(1,2,3)) # the longer solution
interpolation_sol = solve(prob,Vern7(),saveat=t,reltol=1e-12,abstol=1e-12)
plot(interpolation_sol,vars=(1,2,3))

xyzt = plot(data_sol_short, plotdensity=10000,lw=1.5)
xy = plot(data_sol_short, plotdensity=10000, vars=(1,2))
xz = plot(data_sol_short, plotdensity=10000, vars=(1,3))
yz = plot(data_sol_short, plotdensity=10000, vars=(2,3))
xyz = plot(data_sol_short, plotdensity=10000, vars=(1,2,3))
plot(plot(xyzt,xyz),plot(xy, xz, yz, layout=(1,3),w=1), layout=(2,1), size=(800,600))

xyzt = plot(data_sol, plotdensity=10000,lw=1.5)
xy = plot(data_sol, plotdensity=10000, vars=(1,2))
xz = plot(data_sol, plotdensity=10000, vars=(1,3))
yz = plot(data_sol, plotdensity=10000, vars=(2,3))
xyz = plot(data_sol, plotdensity=10000, vars=(1,2,3))
plot(plot(xyzt,xyz),plot(xy, xz, yz, layout=(1,3),w=1), layout=(2,1), size=(800,600))

obj_short = build_loss_objective(prob_short,Tsit5(),L2Loss(t_short,data_short),tstops=t_short)
res1 = bboptimize(obj_short;SearchRange = LooserBounds, MaxSteps = 7e3)
# Tolernace is still too high to get close enough

obj_short = build_loss_objective(prob_short,Tsit5(),L2Loss(t_short,data_short),tstops=t_short,reltol=1e-9)
res1 = bboptimize(obj_short;SearchRange = LooserBounds, MaxSteps = 7e3)
# With the tolerance lower, it achieves the correct solution in 3.5 seconds.

obj_short = build_loss_objective(prob_short,Vern9(),L2Loss(t_short,data_short),tstops=t_short,reltol=1e-9,abstol=1e-9)
res1 = bboptimize(obj_short;SearchRange = LooserBounds, MaxSteps = 7e3)
# With the more accurate solver Vern9 in the solution of the ODE, the convergence is less efficient!

# Fastest BlackBoxOptim: 3.5 seconds

obj_short = build_loss_objective(prob_short,Vern9(),L2Loss(t_short,data_short),tstops=t_short,reltol=1e-9,abstol=1e-9)

opt = Opt(:GN_ORIG_DIRECT_L, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Accurate 3.2 seconds

opt = Opt(:GN_CRS2_LM, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Accurate 3.0 seconds

opt = Opt(:GN_ISRES, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Accurate to single precision 8.2 seconds

opt = Opt(:GN_ESCH, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Approximatively accurate, good starting values for local optimization

opt = Opt(:LN_BOBYQA, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # 0.1 seconds

opt = Opt(:LN_NELDERMEAD, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.29 sec

opt = Opt(:LD_SLSQP, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.21 sec

opt = Opt(:LN_COBYLA, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 1.84 sec

opt = Opt(:LN_NEWUOA_BOUND, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.18 sec ROUNDOFF LIMITED

opt = Opt(:LN_PRAXIS, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.18 sec

opt = Opt(:LN_SBPLX, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.65 sec

opt = Opt(:LD_MMA, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.7 sec

opt = Opt(:LD_LBFGS, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.12 sec

opt = Opt(:LD_TNEWTON_PRECOND_RESTART, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj_short.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Accurate 0.15 sec

# BB with Vern9 converges very slowly. The final values are within the NarrowBounds.
obj = build_loss_objective(prob,Vern9(),L2Loss(t,data),tstops=t,reltol=1e-9,abstol=1e-9)

res1 = bboptimize(obj;SearchRange = LooserBounds, MaxSteps = 4e3) # Default adaptive_de_rand_1_bin_radiuslimited 33 sec [10.2183, 24.6711, 2.28969]
#res1 = bboptimize(obj;SearchRange = LooserBounds, Method = :adaptive_de_rand_1_bin, MaxSteps = 4e3) # Method 32 sec [13.2222, 25.8589, 2.56176]
#res1 = bboptimize(obj;SearchRange = LooserBounds, Method = :dxnes, MaxSteps = 2e3) # Method dxnes 119 sec  [16.8648, 24.393, 2.29119]
#res1 = bboptimize(obj;SearchRange = LooserBounds, Method = :xnes, MaxSteps = 2e3) # Method xnes 304 sec  [19.1647, 24.9479, 2.39467]
#res1 = bboptimize(obj;SearchRange = LooserBounds, Method = :de_rand_1_bin_radiuslimited, MaxSteps = 2e3) # Method 44 sec  [13.805, 24.6054, 2.37274]
#res1 = bboptimize(obj;SearchRange = LooserBounds, Method = :generating_set_search, MaxSteps = 2e3) # Method 195 sec [19.1847, 24.9492, 2.39412]

# using Evolutionary
# N = 3
# @time result, fitness, cnt = cmaes(obj, N; μ = 3, λ = 12, iterations = 1000) # cmaes( rastrigin, N; μ = 15, λ = P, tol = 1e-8)

opt = Opt(:GN_ORIG_DIRECT_L, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Fail to converge

opt = Opt(:GN_CRS2_LM, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 20000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Hit and miss. converge approximately accurate values for local opt.91 seconds

opt = Opt(:GN_ISRES, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 50000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Approximately accurate within local bounds

opt = Opt(:GN_ESCH, 3)
lower_bounds!(opt,[0.0,0.0,0.0])
upper_bounds!(opt,[22.0,60.0,6.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 20000)
@time (minf,minx,ret) = NLopt.optimize(opt,GloIniPar) # Approximately accurate

opt = Opt(:LN_BOBYQA, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Claims SUCCESS but does not iterate to the true values.

opt = Opt(:LN_NELDERMEAD, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-9)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Inaccurate final values

opt = Opt(:LD_SLSQP, 3)
lower_bounds!(opt,[9.0,20.0,2.0])
upper_bounds!(opt,[11.0,30.0,3.0])
min_objective!(opt, obj.cost_function2)
xtol_rel!(opt,1e-12)
maxeval!(opt, 10000)
@time (minf,minx,ret) = NLopt.optimize(opt,LocIniPar) # Inaccurate final values

obj_short = build_loss_objective(prob_short,Tsit5(),L2Loss(t_short,data_short),tstops=t_short)
lower = [0.0,0.0,0.0]
upper = [50.0,50.0,50.0]
splits = ([1.0,5.0,15.0],[0,10,20],[0,10,20])
@time root, x0 = analyze(obj_short,splits,lower,upper)

minimum(root)

obj = build_loss_objective(prob,Vern9(),L2Loss(t,data),tstops=t,reltol=1e-9,abstol=1e-9)
lower = [0.0,0.0,0.0]
upper = [50.0,50.0,50.0]
splits = ([0,5.0,15.0],[0,15,30],[0,2,5])
@time root, x0 = analyze(obj,splits,lower,upper)

minimum(root)


text:
using Libdl
"""
    PetscLibType{PetscScalar, PetscInt}(petsc_library)

A container for specific PETSc libraries.

All other containers for PETSc objects should be typed on this to ensure that
dispatch is correct.
"""
struct PetscLibType{PetscScalar, PetscInt, LibType}
    petsc_library::LibType
end
function PetscLibType{ST, IT}(petsc_library) where {ST, IT}
    LT = typeof(petsc_library)
    return PetscLibType{ST, IT, LT}(petsc_library)
end

"""
    scalartype(petsclib::PetscLibType)

return the scalar type for the associated `petsclib`
"""
scalartype(::PetscLibType{ST}) where {ST} = ST
scalartype(::Type{PetscLib}) where {PetscLib <: PetscLibType{ST}} where {ST} =
    ST

"""
    realtype(petsclib::PetscLibType)

return the real type for the associated `petsclib`
"""
realtype(::PetscLibType{ST}) where {ST} = real(ST)
realtype(::Type{PetscLib}) where {PetscLib <: PetscLibType{ST}} where {ST} =
    real(ST)

"""
    inttype(petsclib::PetscLibType)

return the int type for the associated `petsclib`
"""
inttype(::PetscLibType{ST, IT}) where {ST, IT} = IT
inttype(
    ::Type{PetscLib},
) where {PetscLib <: PetscLibType{ST, IT}} where {ST, IT} = IT

function initialize(libhdl::Ptr{Cvoid})
    PetscInitializeNoArguments_ptr = dlsym(libhdl, :PetscInitializeNoArguments)
    @chk ccall(PetscInitializeNoArguments_ptr, PetscErrorCode, ())
end

const petsclibs = map(libs) do lib
    libhdl = dlopen(lib...)

    # initialize petsc
    PetscInitializeNoArguments_ptr =
        dlsym(libhdl, :PetscInitializeNoArguments)
    @chk ccall(PetscInitializeNoArguments_ptr, PetscErrorCode, ())

    PETSC_REAL = DataTypeFromString(libhdl, "Real")
    PETSC_SCALAR = DataTypeFromString(libhdl, "Scalar")
    PETSC_INT_SIZE = PetscDataTypeGetSize(libhdl, PETSC_INT)

    PetscReal =
        PETSC_REAL == PETSC_DOUBLE ? Cdouble :
        PETSC_REAL == PETSC_FLOAT ? Cfloat :
        error("PETSC_REAL = $PETSC_REAL not supported.")

    PetscScalar =
        PETSC_SCALAR == PETSC_REAL ? PetscReal :
        PETSC_SCALAR == PETSC_COMPLEX ? Complex{PetscReal} :
        error("PETSC_SCALAR = $PETSC_SCALAR not supported.")

    PetscInt =
        PETSC_INT_SIZE == 4 ? Int32 :
        PETSC_INT_SIZE == 8 ? Int64 :
        error("PETSC_INT_SIZE = $PETSC_INT_SIZE not supported.")

    PetscFinalize_ptr = dlsym(libhdl, :PetscFinalize)
    @chk ccall(PetscFinalize_ptr, PetscErrorCode, ())

    # TODO: PetscBLASInt, PetscMPIInt ?
    return PetscLibType{PetscScalar, PetscInt}(lib[1])
end

# New macro is really to track the update
macro for_petsc(expr)
    quote
        for petsclib in petsclibs
            # String for the library
            petsc_library = petsclib.petsc_library

            # types we dispatch on
            PetscLib = typeof(petsclib)
            UnionPetscLib = Union{PetscLib, Type{PetscLib}}

            PetscScalar = scalartype(petsclib)
            PetscReal = realtype(petsclib)
            PetscInt = inttype(petsclib)

            @eval esc($expr)
        end
    end
end

# TODO: Remove this macro once we have switch over to @for_petsc macro
macro for_libpetsc(expr)
    quote
        for petsclib in petsclibs
            # String for the library
            libpetsc = petsclib.petsc_library

            # types we dispatch on
            PetscLib = typeof(petsclib)
            UnionPetscLib = Union{PetscLib, Type{PetscLib}}

            PetscScalar = scalartype(petsclib)
            PetscReal = realtype(petsclib)
            PetscInt = inttype(petsclib)

            @eval esc($expr)
        end
    end
end

# TODO: These should be removed!
const scalar_types = map(x -> scalartype(x), petsclibs)
@assert length(scalar_types) == length(unique(scalar_types))
@for_libpetsc inttype(::Type{$PetscScalar}) = $PetscInt


text:
import Base.zero, Base.rand, Base.convert, Base.abs, Base.typemax, Base.typemin
"""
Type to specify what bootstrap to run, and to hold the resulting bootstrapped values
"""
abstract type AbstractBootstrapper{T} end

mutable struct GenericBootstrapper{T} <: AbstractBootstrapper{T}
    nbootstrap::Int64
    func::Function
    values::Array{T,1}
end

mutable struct Bootstrapper{T<:Real} <: AbstractBootstrapper{T}
    nbootstrap::Int64
    func::Function
    values::Array{T,1}
    μ::T
    σ::T
end

function Bootstrapper(::Type{T},nbootstrap::Int64,func::Function) where T <: Real
    #test the function
    a = func(rand(T,2))
    if length(a) > 1
        throw(ArgumentError("func should return a scalar value"))
    end
    Bootstrapper(nbootstrap,func,zeros(T,nbootstrap),zero(T),zero(T))
end

function summarize!(B::AbstractBootstrapper{T}) where T<:Real
    B.μ = mean(B.values)
    B.σ = std(B.values)
end

"""
Compute bootstrap statistics `nbootstrap` times using `func`

    function run_bootstrap{T}(nbootstrap::Int64,func::Function, X::Array{T,1},args...)
"""
function run_bootstrap(nbootstrap::Int64,func::Function, X::Array{T,1},args...;RNG=MersenneTwister(rand(UInt32))) where T <: Any
    B = Bootstrapper(T, nbootstrap, func)
    run_bootstrap!(B,X,args...;RNG=RNG)
    B
end

function run_bootstrap!(B::AbstractBootstrapper{T}, X::Array{T,1},args...;RNG=MersenneTwister(rand(UInt32))) where T <: Any
    n = length(X)
    idx = zeros(Int64,n)
    for i in 1:B.nbootstrap
        rand!(RNG, idx,1:n)
        B.values[i] = B.func(X[idx],args...)
    end
    summarize!(B)
end

struct ZScore
	v::Float64
end

zero(::Type{ZScore}) = ZScore(zero(Float64))
rand(::Type{ZScore}) = randn()
convert(::Type{ZScore}, X::Float64) = ZScore(X)
convert(::Type{Float64}, X::ZScore) = X.v
abs(z::ZScore) = abs(z.v)

struct Percentile
    v::Float64
    Percentile(v) = 0.0 < v < 100.0 ? new(v) : throw(ArgumentError("Percentiles must be between 0.0 and 100.0"))
end

Base.typemax(::Type{Percentile}) = 100.0
Base.typemin(::Type{Percentile}) = 0.0


text:
module Iris_Tests
using Test
using DataDeps
using MLDatasets

data_dir = withenv("DATADEPS_ALWAY_ACCEPT"=>"true") do
    datadep"Iris"
end

@testset "Iris" begin
    X  = Iris.features()
    Y  = Iris.labels()
    @test X isa Matrix{Float64}
    @test Y isa Vector{String}
    @test size(X) == (4, 150)
    @test size(Y) == (150,)
end

end


text:
function match(cur_reactor::AbstractReactor, cur_constraint::Symbol, cur_value::Union{Void, Number}=nothing, is_singular::Bool=true)
  cur_func = function (cur_T_bar::Number)
    tmp_reactor = deepcopy(cur_reactor)

    tmp_reactor.T_bar = real(cur_T_bar)

    tmp_reactor.T_bar <= 0 && return NaN

    tmp_reactor.sigma_v = calc_sigma_v(tmp_reactor)

    isnan(tmp_reactor.sigma_v) && return NaN

    cur_error = calc_I_P(tmp_reactor)

    cur_error /= symbols(:I_P)

    cur_error -= 1.0

    secondary_constraint = getfield(Fussy, Symbol("$(tmp_reactor.constraint)_equation"))(tmp_reactor)

    if cur_value == nothing
      tertiary_constraint = getfield(Fussy, Symbol("$(cur_constraint)_equation"))(tmp_reactor)
    else
      tertiary_constraint = getfield(Fussy, Symbol("$(cur_constraint)_equation"))(tmp_reactor, cur_value)
    end

    tmp_dict = equation_set_dict(
      EquationSet(tmp_reactor, secondary_constraint, tertiary_constraint)
    )

    any(is_nan, values(tmp_dict)) && return NaN

    @assert all(is_finite, values(tmp_dict))

    cur_error = subs(cur_error, tmp_dict...)

    cur_error = real(float(cur_error))

    cur_error
  end

  matched_reactors = _match(cur_reactor, cur_func, cur_constraint, cur_value)
  matched_reactors = filter!(work_reactor -> work_reactor.is_good, matched_reactors)

  isempty(matched_reactors) && return ( is_singular ? nothing : [] )

  if length(matched_reactors) > 1
    work_reactors = filter(work_reactor -> work_reactor.is_valid, matched_reactors)

    if !isempty(work_reactors)
      is_singular || return work_reactors

      if length(work_reactors) > 1
        custom_log("Bad match! Too many reactors!")
        custom_log(cur_constraint, " - ", cur_value)
        custom_log(work_reactors)
      end

      min_index = indmin(
        map(work_reactor -> (cur_constraint == :cost ? work_reactor.W_M : work_reactor.cost), work_reactors)
      )

      return work_reactors[min_index]
    end
  end

  is_singular && return first(matched_reactors)

  matched_reactors
end

function _match(cur_reactor::AbstractReactor, cur_function::Function, cur_constraint::Symbol, cur_value::Union{Void, Number}=nothing)
  cur_atol = 10 * eps()

  cur_T_bar_list = find_roots(
    cur_function, min_T_bar, max_T_bar,
    no_pts = Int(ceil(no_pts_T_bar/2))+1, reltol = 1e-2, abstol = 1e-4
  )

  isempty(cur_T_bar_list) && return []

  matched_reactors = []

  for (cur_index, cur_T_bar) in enumerate(cur_T_bar_list)
    isapprox(0.0, cur_T_bar, atol=cur_atol) && continue

    tmp_reactor = deepcopy(cur_reactor)

    tmp_reactor.T_bar = cur_T_bar

    tmp_reactor.sigma_v = calc_sigma_v(tmp_reactor)

    secondary_constraint = getfield(Fussy, Symbol("$(tmp_reactor.constraint)_equation"))(tmp_reactor)

    if cur_value == nothing
      tertiary_constraint = getfield(Fussy, Symbol("$(cur_constraint)_equation"))(tmp_reactor)
    else
      tertiary_constraint = getfield(Fussy, Symbol("$(cur_constraint)_equation"))(tmp_reactor, cur_value)
    end

    cur_I_P = EquationSet(tmp_reactor, secondary_constraint, tertiary_constraint).I_P

    isnan(cur_I_P) && continue

    cur_I_P = float(cur_I_P)

    isa(cur_I_P, Complex) && continue

    tmp_reactor.I_P = float(cur_I_P)

    push!(matched_reactors, update!(tmp_reactor))
  end

  matched_reactors
end


text:
function compare_to(::Type{Precision}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(Precision, "compareTo", jint, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function compare_to(::Type{Precision}, arg0::jdouble, arg1::jdouble, arg2::jint)
    return jcall(Precision, "compareTo", jint, (jdouble, jdouble, jint), arg0, arg1, arg2)
end

function equals(::Type{Precision}, arg0::jdouble, arg1::jdouble)
    return jcall(Precision, "equals", jboolean, (jdouble, jdouble), arg0, arg1)
end

function equals(::Type{Precision}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(Precision, "equals", jboolean, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function equals(::Type{Precision}, arg0::jdouble, arg1::jdouble, arg2::jint)
    return jcall(Precision, "equals", jboolean, (jdouble, jdouble, jint), arg0, arg1, arg2)
end

function equals(::Type{Precision}, arg0::jfloat, arg1::jfloat)
    return jcall(Precision, "equals", jboolean, (jfloat, jfloat), arg0, arg1)
end

function equals(::Type{Precision}, arg0::jfloat, arg1::jfloat, arg2::jfloat)
    return jcall(Precision, "equals", jboolean, (jfloat, jfloat, jfloat), arg0, arg1, arg2)
end

function equals(::Type{Precision}, arg0::jfloat, arg1::jfloat, arg2::jint)
    return jcall(Precision, "equals", jboolean, (jfloat, jfloat, jint), arg0, arg1, arg2)
end

function equals_including_na_n(::Type{Precision}, arg0::jdouble, arg1::jdouble)
    return jcall(Precision, "equalsIncludingNaN", jboolean, (jdouble, jdouble), arg0, arg1)
end

function equals_including_na_n(::Type{Precision}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(Precision, "equalsIncludingNaN", jboolean, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function equals_including_na_n(::Type{Precision}, arg0::jdouble, arg1::jdouble, arg2::jint)
    return jcall(Precision, "equalsIncludingNaN", jboolean, (jdouble, jdouble, jint), arg0, arg1, arg2)
end

function equals_including_na_n(::Type{Precision}, arg0::jfloat, arg1::jfloat)
    return jcall(Precision, "equalsIncludingNaN", jboolean, (jfloat, jfloat), arg0, arg1)
end

function equals_including_na_n(::Type{Precision}, arg0::jfloat, arg1::jfloat, arg2::jfloat)
    return jcall(Precision, "equalsIncludingNaN", jboolean, (jfloat, jfloat, jfloat), arg0, arg1, arg2)
end

function equals_including_na_n(::Type{Precision}, arg0::jfloat, arg1::jfloat, arg2::jint)
    return jcall(Precision, "equalsIncludingNaN", jboolean, (jfloat, jfloat, jint), arg0, arg1, arg2)
end

function equals_with_relative_tolerance(::Type{Precision}, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(Precision, "equalsWithRelativeTolerance", jboolean, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function is_mathematical_integer(::Type{Precision}, arg0::jdouble)
    return jcall(Precision, "isMathematicalInteger", jboolean, (jdouble,), arg0)
end

function is_mathematical_integer(::Type{Precision}, arg0::jfloat)
    return jcall(Precision, "isMathematicalInteger", jboolean, (jfloat,), arg0)
end

function representable_delta(::Type{Precision}, arg0::jdouble, arg1::jdouble)
    return jcall(Precision, "representableDelta", jdouble, (jdouble, jdouble), arg0, arg1)
end

function round(::Type{Precision}, arg0::jdouble, arg1::jint)
    return jcall(Precision, "round", jdouble, (jdouble, jint), arg0, arg1)
end

function round(::Type{Precision}, arg0::jdouble, arg1::jint, arg2::jint)
    return jcall(Precision, "round", jdouble, (jdouble, jint, jint), arg0, arg1, arg2)
end

function round(::Type{Precision}, arg0::jfloat, arg1::jint)
    return jcall(Precision, "round", jfloat, (jfloat, jint), arg0, arg1)
end

function round(::Type{Precision}, arg0::jfloat, arg1::jint, arg2::jint)
    return jcall(Precision, "round", jfloat, (jfloat, jint, jint), arg0, arg1, arg2)
end


text:
type State <: AbsState; id; end
type Action <: AbsAction; id; end

==(lhs::State, rhs::State) = lhs.id == rhs.id
isequal(lhs::State, rhs::State) = lhs.id == rhs.id
hash(s::State) = hash(s.id)

==(lhs::Action, rhs::Action) = lhs.id == rhs.id
isequal(lhs::Action, rhs::Action) = lhs.id == rhs.id
hash(a::Action) = hash(a.id)

mutable struct MDP <: AbsEnvironment
    ns#number of states
    na#number of actions
    graph::Dict{State, Dict{Action,Array{Tuple{State, Float64, Float64}}}}#state + action ->  [state, reward, prob]
    MDP(ns, na, graph) = new(ns, na, graph)
    function MDP(ns=3, na=3)
        states = [State(i) for i in 1:ns]
        actions = [Action(i) for i in 1:na]
        graph = Dict{State, Dict{Action,Array{Tuple{State, Float64, Float64}}}}()
        for s in states
            numA = rand(1:na)
            indices = randperm(na)
            graph[s] = Dict{Action, Array{Tuple{State, Float64, Float64}}}()

            for k=1:numA
                a = actions[indices[k]]
                numS = rand(1:ns)
                indicesS = randperm(ns)
                probs = rand(numS,)
                probs = probs ./ sum(probs)

                for j=1:numS
                    sprime = states[indicesS[j]]
                    if j == 1
                        graph[s][a] = [(sprime, rand(-100:100), probs[j])]
                    else
                        push!(graph[s][a], (sprime, rand(-100:100), probs[j]))
                    end
                end
            end
        end
        new(ns, na, graph)
    end
end

getSuccessors(s::State, a::Action, env::MDP) = env.graph[s][a]
getAllStates(env::MDP) = keys(env.graph)
getActions(s::State, env::MDP) = keys(env.graph[s])

function printmdp(mdp::MDP)
    println("Number of states: $(mdp.ns)")
    println("Number of actions: $(mdp.na)\n")

    for s in keys(mdp.graph)
        for a in keys(mdp.graph[s])
            for t in mdp.graph[s][a]
                println("$(s) + $(a) => $(t[1]), reward: $(t[2]), probability: $(t[3])")
            end
            println("")
        end
    end
end

max_stars_repo_name: ozanarkancan/ReinforcementLearning

max_issues_repo_name: ozanarkancan/ReinforcementLearning

max_forks_repo_name: ozanarkancan/ReinforcementLearning


text:
function correlate(w)

del = w[:,1]
E = w[:,2]./del .- 1

F = 4.8274*E.^4 - 5.9816*E.^3 + 4.0274*E.^2 + 0.23247.*E .+ 0.15174

B = 0.5*(-225.86.*E.^3 - 3016.6.*E.^2 - 208.68*E .- 17.915
        + 131.9*E.^3 - 167.32*E.^2 + 76.642.*E .- 11.068)
B[E.<-0.0616] = -225.86.*E[E.<-0.0616].^3 - 3016.6.*E[E.<-0.0616].^2 - 208.68*E[E.<-0.0616] .- 17.915
B[E.>-0.0395] =  131.9*E[E.>-0.0395].^3 - 167.32*E[E.>-0.0395].^2 + 76.642.*E[E.>-0.0395] .- 11.068

S = 0.5*(451.55*E.^3 + 2010*E.^2 + 138.96*E .+ 11.296
                        - 96.739*E.^3 + 117.74*E.^2 - 46.432*E .+ 6.8074)
S[E.<-0.0582] = 451.55*E[E.<-0.0582].^3 + 2010*E[E.< -0.0582].^2 + 138.96*E[E.< -0.0582] .+ 11.296
S[E.>-0.042]  = -96.739*E[E.>-0.042].^3 + 117.74*E[E.> -0.042].^2 - 46.432*E[E.> -0.042] .+ 6.8074

dfde = 4*4.8274*E.^3 - 3*5.9816*E.^2 + 2*4.0274*E .+ 0.23247

return del, E, F ,B, S, dfde

end

max_stars_repo_name: desanga/IBMFVM

max_issues_repo_name: desanga/IBMFVM

max_forks_repo_name: desanga/IBMFVM


text:
@safetestset "isfqdn" begin
    using Validator

    @test Validator.isfqdn("domain.com") === true
    @test Validator.isfqdn("dom.plato") === true
    @test Validator.isfqdn("a.domain.co") === true
    @test Validator.isfqdn("foo--bar.com") === true
    @test Validator.isfqdn("xn--froschgrn-x9a.com") === true
    @test Validator.isfqdn("rebecca.blackfriday") === true
    @test Validator.isfqdn("1337.com") === true

    @test Validator.isfqdn("abc") === false
    @test Validator.isfqdn("256.0.0.0") === false
    @test Validator.isfqdn("_.com") === false
    @test Validator.isfqdn("*.some.com") === false
    @test Validator.isfqdn("s!ome.com") === false
    @test Validator.isfqdn("domain.com/") === false
    @test Validator.isfqdn("/more.com") === false
    @test Validator.isfqdn("domain.com�") === false
    @test Validator.isfqdn("domain.co\u00A0m") === false
    @test Validator.isfqdn("domain.co\u1680m") === false
    @test Validator.isfqdn("domain.co\u2006m") === false
    @test Validator.isfqdn("domain.co\u2028m") === false
    @test Validator.isfqdn("domain.co\u2029m") === false
    @test Validator.isfqdn("domain.co\u202Fm") === false
    @test Validator.isfqdn("domain.co\u205Fm") === false
    @test Validator.isfqdn("domain.co\u3000m") === false
    @test Validator.isfqdn("domain.com\uDC00") === false
    @test Validator.isfqdn("domain.co\uEFFFm") === false
    @test Validator.isfqdn("domain.co\uFDDAm") === false
    @test Validator.isfqdn("domain.co\uFFF4m") === false
    @test Validator.isfqdn("domain.com©") === false
    @test Validator.isfqdn("example.0") === false
    @test Validator.isfqdn("192.168.0.9999") === false
    @test Validator.isfqdn("192.168.0") === false

    # allowTrailingDot = true
    @test Validator.isfqdn("example.com.", true, false, true) === true

    @test Validator.isfqdn("example.0", true, false, true) === false
    @test Validator.isfqdn("192.168.0", true, false, true) === false
    @test Validator.isfqdn("192.168.0.9999", true, false, true) === false

    # requireTld = false
    @test Validator.isfqdn("example.0", false) === false
    @test Validator.isfqdn("192.168.0", false) === false
    @test Validator.isfqdn("192.168.0.9999", false) === false

    # requireTld = false, allowNumericTld = true
    @test Validator.isfqdn("example.0", false, false, false, true) === true
    @test Validator.isfqdn("192.168.0", false, false, false, true) === true
    @test Validator.isfqdn("192.168.0.9999", false, false, false, true) === true

    # allowWildcard = true
    @test Validator.isfqdn("*.example.com", true, false, false, false, true) === true
    @test Validator.isfqdn("*.shop.example.com", true, false, false, false, true) === true

end


text:
"""
    FlipSignBridge{T, S1, S2}

Bridge constrained variables in `S1` into constrained variables in `S2` by
multiplying the variables by `-1` and taking the point reflection of the set
across the origin. The flipped `MOI.VectorOfVariables`-in-`S` constraint is
stored in the `flipped_constraint` field by convention.
"""
abstract type FlipSignBridge{
    T, S1<:MOI.AbstractSet, S2<:MOI.AbstractSet} <: AbstractBridge end

function supports_constrained_variable(
    ::Type{<:FlipSignBridge{T, S1}}, ::Type{S1}) where {T, S1<:MOI.AbstractVectorSet}
    return true
end
function MOIB.added_constrained_variable_types(
    ::Type{<:FlipSignBridge{T, S1, S2}}) where {T, S1, S2}
    return [(S2,)]
end
function MOIB.added_constraint_types(::Type{<:FlipSignBridge})
    return Tuple{DataType, DataType}[]
end

# Attributes, Bridge acting as a model
function MOI.get(bridge::FlipSignBridge, ::MOI.NumberOfVariables)
    return length(bridge.flipped_variables)
end
function MOI.get(bridge::FlipSignBridge, ::MOI.ListOfVariableIndices)
    return bridge.flipped_variables
end
function MOI.get(::FlipSignBridge{T, S1, S2},
                 ::MOI.NumberOfConstraints{MOI.VectorOfVariables, S2}) where {T, S1, S2<:MOI.AbstractVectorSet}
    return 1
end
function MOI.get(bridge::FlipSignBridge{T, S1, S2},
                 ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables, S2}) where {T, S1, S2<:MOI.AbstractVectorSet}
    return [bridge.flipped_constraint]
end

# References
function MOI.delete(model::MOI.ModelLike, bridge::FlipSignBridge)
    MOI.delete(model, bridge.flipped_variables)
end

function MOI.delete(model::MOI.ModelLike, bridge::FlipSignBridge, i::IndexInVector)
    MOI.delete(model, bridge.flipped_variables[i.value])
    deleteat!(bridge.flipped_variables, i.value)
end

# Attributes, Bridge acting as a constraint
function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintSet, bridge::FlipSignBridge{T, S1}) where {T, S1<:MOI.AbstractVectorSet}
    return S1(length(bridge.flipped_variables))
end

function MOI.get(model::MOI.ModelLike,
                 attr::Union{MOI.ConstraintPrimal, MOI.ConstraintDual},
                 bridge::FlipSignBridge)
    return -MOI.get(model, attr, bridge.flipped_constraint)
end

function MOI.get(model::MOI.ModelLike,
                 attr::Union{MOI.VariablePrimal, MOI.VariablePrimalStart},
                 bridge::FlipSignBridge, i::IndexInVector)
    return -MOI.get(model, attr, bridge.flipped_variables[i.value])
end

function MOIB.bridged_function(bridge::FlipSignBridge{T}, i::IndexInVector) where T
    func = MOI.SingleVariable(bridge.flipped_variables[i.value])
    return MOIU.operate(-, T, func)
end
function unbridged_map(bridge::FlipSignBridge{T}, vi::MOI.VariableIndex,
                       i::IndexInVector) where T
    func = MOIU.operate(-, T, MOI.SingleVariable(vi))
    return (bridge.flipped_variables[i.value] => func,)
end

function MOI.supports(model::MOI.ModelLike, attr::MOI.VariablePrimalStart,
                      ::Type{<:FlipSignBridge})
    return MOI.supports(model, attr, MOI.VariableIndex)
end
function MOI.set(model::MOI.ModelLike, attr::MOI.VariablePrimalStart,
                 bridge::FlipSignBridge, value, i::IndexInVector)
    MOI.set(model, attr, bridge.flipped_variables[i.value], -value)
end

"""
    NonposToNonnegBridge{T, F<:MOI.AbstractVectorFunction, G<:MOI.AbstractVectorFunction} <:
        FlipSignBridge{T, MOI.Nonpositives, MOI.Nonnegatives, F, G}

Transforms constrained variables in `Nonpositives` into constrained variables in
`Nonnegatives`.
"""
struct NonposToNonnegBridge{T} <: FlipSignBridge{T, MOI.Nonpositives, MOI.Nonnegatives}
    flipped_variables::Vector{MOI.VariableIndex}
    flipped_constraint::MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.Nonnegatives}
end
function bridge_constrained_variable(::Type{NonposToNonnegBridge{T}},
                                     model::MOI.ModelLike,
                                     set::MOI.Nonpositives) where T
    flipped_variables, flipped_constraint = MOI.add_constrained_variables(
        model, MOI.Nonnegatives(set.dimension))
    return NonposToNonnegBridge{T}(flipped_variables, flipped_constraint)
end


text:
qmax_default(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = 9//8
qmin_default(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = 1//5

delta_default(alg) = 1//1
delta_default(alg::SRIW1) = 1//6

isadaptive(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = false
isadaptive(alg::Union{StochasticDiffEqAdaptiveAlgorithm,StochasticDiffEqRODEAdaptiveAlgorithm}) = true
isadaptive(alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = all(isadaptive.(alg.algs))

alg_order(alg::EM) = 1//2
alg_order(alg::LambaEM) = 1//2
alg_order(alg::ImplicitEM) = 1//2
alg_order(alg::ImplicitEulerHeun) = 1//2
alg_order(alg::ImplicitRKMil) = 1//1
alg_order(alg::WangLi3SMil_A) = 1//1
alg_order(alg::WangLi3SMil_B) = 1//1
alg_order(alg::WangLi3SMil_C) = 1//1
alg_order(alg::WangLi3SMil_D) = 1//1
alg_order(alg::WangLi3SMil_E) = 1//1
alg_order(alg::WangLi3SMil_F) = 1//1
alg_order(alg::ISSEM) = 1//2
alg_order(alg::ISSEulerHeun) = 1//2
alg_order(alg::SplitEM) = 1//2
alg_order(alg::PCEuler) = 1//2
alg_order(alg::IIF1M) = 1//2
alg_order(alg::IIF2M) = 1//2
alg_order(alg::IIF1Mil) = 1//1
alg_order(alg::EulerHeun) = 1//2
alg_order(alg::LambaEulerHeun) = 1//2
alg_order(alg::RandomEM) = 1//2
alg_order(alg::RKMil) = 1//1
alg_order(alg::RKMilCommute) = 1//1
alg_order(alg::RKMil_General) = 1//1

# Generalised version of SROCK1, both Ito ans Stratonovich, will have strong order of 1//2
# and weak order of 1 for Multidimensional Weiner process
# Stratonovich version strong order 1 for 1 dimensional Weiner Process or if noise is commutative
# Ito version can have strong order version for 1 dimensional Weiner Process,
# diagonal noise or commutative noise
alg_order(alg::SROCK1) = 1//2
alg_order(alg::SROCK2) = 1//1
alg_order(alg::KomBurSROCK2) = 1//1
alg_order(alg::SROCKC2) = 1//1
alg_order(alg::SROCKEM) = alg.strong_order_1 ? 1//1 : 1//2
alg_order(alg::SKSROCK) = 1//2
alg_order(alg::TangXiaoSROCK2) = 1//1

alg_order(alg::SRI) = alg.tableau.order
alg_order(alg::SRIW1) = 3//2
alg_order(alg::SRIW2) = 3//2
alg_order(alg::SOSRI) = 3//2
alg_order(alg::SOSRI2) = 3//2
alg_order(alg::SRA) = alg.tableau.order
alg_order(alg::SRA1) = 2//1
alg_order(alg::SRA2) = 2//1
alg_order(alg::SRA3) = 2//1
alg_order(alg::SOSRA) = 2//1
alg_order(alg::SOSRA2) = 2//1
alg_order(alg::SKenCarp) = 2//1
alg_order(alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = maximum(alg_order.(alg.algs))
get_current_alg_order(alg::StochasticDiffEqAlgorithm,cache) = alg_order(alg)
get_current_alg_order(alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm},cache) = alg_order(alg.algs[cache.current])

beta2_default(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = 2//(5alg_order(alg))
beta1_default(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm},beta2) = 7//(10alg_order(alg))

isdtchangeable(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = true

alg_interpretation(alg::StochasticDiffEqAlgorithm) = :Ito
alg_interpretation(alg::EulerHeun) = :Stratonovich
alg_interpretation(alg::LambaEulerHeun) = :Stratonovich
alg_interpretation(alg::KomBurSROCK2) = :Stratonovich
alg_interpretation(alg::RKMil{interpretation}) where {interpretation} = interpretation
alg_interpretation(alg::SROCK1{interpretation}) where {interpretation} = interpretation
alg_interpretation(alg::RKMilCommute{interpretation}) where {interpretation} = interpretation
alg_interpretation(alg::RKMil_General) = alg.interpretation
alg_interpretation(alg::ImplicitRKMil{CS,AD,F,S,N,T2,Controller,interpretation}) where {CS,AD,F,S,N,T2,Controller,interpretation} = interpretation

alg_compatible(prob,alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = true

alg_compatible(prob,alg::StochasticDiffEqAlgorithm) = false
alg_compatible(prob,alg::SRI) = is_diagonal_noise(prob)
alg_compatible(prob,alg::SRIW1) = is_diagonal_noise(prob)
alg_compatible(prob,alg::SRIW2) = is_diagonal_noise(prob)
alg_compatible(prob,alg::SOSRI) = is_diagonal_noise(prob)
alg_compatible(prob,alg::SOSRI2) = is_diagonal_noise(prob)
alg_compatible(prob,alg::SRA) = true
alg_compatible(prob,alg::SRA1) = true
alg_compatible(prob,alg::SRA2) = true
alg_compatible(prob,alg::SRA3) = true
alg_compatible(prob,alg::SOSRA) = true
alg_compatible(prob,alg::SOSRA2) = true
alg_compatible(prob,alg::SKenCarp) = true
alg_compatible(prob,alg::EM) = true
alg_compatible(prob,alg::LambaEM) = true
alg_compatible(prob,alg::WangLi3SMil_A) = true
alg_compatible(prob,alg::WangLi3SMil_B) = true
alg_compatible(prob,alg::WangLi3SMil_C) = true
alg_compatible(prob,alg::WangLi3SMil_D) = true
alg_compatible(prob,alg::WangLi3SMil_E) = true
alg_compatible(prob,alg::WangLi3SMil_F) = true
alg_compatible(prob,alg::SROCK1) = true
alg_compatible(prob,alg::SROCK2) = true
alg_compatible(prob,alg::KomBurSROCK2) = true
alg_compatible(prob,alg::SROCKC2) = true
alg_compatible(prob,alg::SROCKEM) = true
alg_compatible(prob,alg::SKSROCK) = true
alg_compatible(prob,alg::TangXiaoSROCK2) = true
alg_compatible(prob,alg::EulerHeun) = true
alg_compatible(prob,alg::LambaEulerHeun) = true
alg_compatible(prob,alg::SplitEM) = true
alg_compatible(prob,alg::PCEuler) = true
alg_compatible(prob,alg::ImplicitEM) = true
alg_compatible(prob,alg::ImplicitEulerHeun) = true
alg_compatible(prob,alg::ISSEM) = true
alg_compatible(prob,alg::ISSEulerHeun) = true
alg_compatible(prob,alg::RKMil) = is_diagonal_noise(prob)
alg_compatible(prob,alg::ImplicitRKMil) = is_diagonal_noise(prob)
alg_compatible(prob,alg::RKMilCommute) = true # No good check for commutative noise
alg_compatible(prob,alg::RKMil_General) = true
alg_compatible(prob,alg::IIF1M) = true
alg_compatible(prob,alg::IIF2M) = true
alg_compatible(prob,alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = max((alg_compatible(prob,a) for a in alg.algs)...)

alg_needs_extra_process(alg::Union{StochasticDiffEqAlgorithm,StochasticDiffEqRODEAlgorithm}) = false
alg_needs_extra_process(alg::Union{StochasticDiffEqCompositeAlgorithm,StochasticDiffEqRODECompositeAlgorithm}) = max((alg_needs_extra_process(a) for a in alg.algs)...)
alg_needs_extra_process(alg::SRI) = true
alg_needs_extra_process(alg::SRIW1) = true
alg_needs_extra_process(alg::SRIW2) = true
alg_needs_extra_process(alg::SOSRI) = true
alg_needs_extra_process(alg::SOSRI2) = true
alg_needs_extra_process(alg::SRA) = true
alg_needs_extra_process(alg::SRA1) = true
alg_needs_extra_process(alg::SRA2) = true
alg_needs_extra_process(alg::SRA3) = true
alg_needs_extra_process(alg::SOSRA) = true
alg_needs_extra_process(alg::SOSRA2) = true
alg_needs_extra_process(alg::SKenCarp) = true

alg_autodiff(alg::StochasticDiffEqNewtonAlgorithm{CS,AD,Controller}) where {CS,AD,Controller} = AD
alg_autodiff(alg::StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}) where {CS,AD,Controller} = AD

get_current_alg_autodiff(alg, cache) = alg_autodiff(alg)
get_current_alg_autodiff(alg::StochasticDiffEqCompositeAlgorithm, cache) = alg_autodiff(alg.algs[cache.current])

get_chunksize(alg::StochasticDiffEqNewtonAlgorithm{CS,AD,Controller}) where {CS,AD,Controller} = CS
get_chunksize(alg::StochasticDiffEqNewtonAdaptiveAlgorithm{CS,AD,Controller}) where {CS,AD,Controller} = CS

alg_mass_matrix_compatible(alg::StochasticDiffEqAlgorithm) = false

alg_can_repeat_jac(alg::StochasticDiffEqAlgorithm) = true

function alg_mass_matrix_compatible(alg::Union{StochasticDiffEqNewtonAlgorithm,StochasticDiffEqNewtonAdaptiveAlgorithm})
    if alg.symplectic
        return true
    elseif alg.theta == 1
        return true
    else
        error("Algorithm must be set as symplectic or theta=1 for mass matrices")
    end
end

is_split_step(::StochasticDiffEqAlgorithm) = false
is_split_step(::EM{split}) = split
is_split_step(::LambaEM{split}) = split

alg_stability_size(alg::SOSRI2) = 10.6
alg_stability_size(alg::SOSRA2) = 5.3

is_composite(alg) = false
is_composite(alg::StochasticDiffEqCompositeAlgorithm) = true
is_composite(alg::StochasticDiffEqRODECompositeAlgorithm) = true
function unwrap_alg(integrator, is_nlsolve)
  alg = integrator.alg
  if !is_composite(alg)
    return alg
  elseif typeof(alg.choice_function) <: AutoSwitch
    num = is_nlsolve ? 2 : 1
    return alg.algs[num]
  else
    return alg.algs[integrator.cache.current]
  end
end

issplit(::StochasticDiffEqAlgorithm) = false
issplit(::SplitSDEAlgorithms) = true


text:
#=
    Recommendation API
    ------------------

Recommendations are performed by calling:

    recommend(env, request)

where:
    `env::SearchEnv` is the search environment object
    `request::InternalRequest` is the request

The recommender function signatures should be of the form:

    some_recommender(request; environment=nothing)

where:
    `some_recommender` is the name of the recommender
    `request` is the recommendation request, containing all needed parameters
    `environment` a named tuple with at least two fields
        • `dbdata` - an `IndexedTable` or `NDSparse` object containint the data
        • `id_key` - the name of the data primary key
        • `ranker` - the ranking function

The arguments above should be enough to implement any ranker.
=#
function recommend(env::SearchEnv, request)
    recommender = safe_symbol_eval(request.recommender, DEFAULT_RECOMMENDER_NAME)
    if recommender !== search_recommender
        # Generic recommender, does not need full environment i.e. strip indexes
        env = build_data_env(env)
    end
    ## Call recommender
    recommender(request; environment=env)
    ##
end


text:
for N in [Float64, Rational{Int}, Float32]
    B1 = BallInf(zeros(N, 2), N(1))
    B2 = Ball1(ones(N, 2), N(1))
    B3 = Hyperrectangle(low=N[-1, -1], high=N[2, 2])
    S = Singleton(N[10, 10])
    UXY = UnionSet(B1, B2)

    # type alias
    @test B1 ∪ B2 == UXY

    # array type (union of a finite number of convex sets)
    Uarr = UnionSetArray([B1, B2])

    # constructor without argumenst
    UnionSetArray()

    # getindex & length
    @test Uarr[1] == B1 && Uarr[2] == B2
    @test length(Uarr) == 2

    # swap
    U2 = swap(UXY)
    @test UXY.X == U2.Y && UXY.Y == U2.X

    # neutral and absorbing elements
    E = EmptySet{N}(2)
    @test neutral(UnionSet) == neutral(UnionSetArray) == EmptySet
    @test UnionSet(B1, E) == UnionSet(E, B1) == B1
    U = Universe{N}(2)
    @test absorbing(UnionSet) == absorbing(UnionSetArray) == Universe
    @test UnionSet(B1, U) == UnionSet(U, B1) == U

    for U in [UXY, Uarr]
        # dimension
        @test dim(U) == dim(B1)

        # support vector (default algorithm)
        d = N[1, 0]
        @test σ(d, U) == [N(2), N(1)]
        @test σ(d, U, algorithm="support_vector") == [N(2), N(1)]

        # support vector (support function algorithm)
        @test σ(d, U, algorithm="support_function") == [N(2), N(1)]

        # support function
        @test ρ(d, U) == N(2)

        # an_element and membership
        @test an_element(U) ∈ U

        # emptiness
        @test !isempty(U)

        # boundedness
        @test isbounded(U) && isboundedtype(typeof(U))

        # inclusion
        subset, point = ⊆(U, B3, true)
        @test U ⊆ B3 && subset && point == N[]
        subset, point = ⊆(U, B2, true)
        @test U ⊈ B2 && !subset && point ∈ U && point ∉ B2

        # isdisjoint
        disjoint1, point1 = isdisjoint(U, B1, true)
        disjoint2, point2 = isdisjoint(B1, U, true)
        @test !isdisjoint(U, B1) && !isdisjoint(B1, U) && !disjoint1 &&
              !disjoint2 && point1 ∈ U && point1 ∈ B1 && point2 ∈ U &&
              point2 ∈ B1
        disjoint1, point1 = isdisjoint(U, S, true)
        disjoint2, point2 = isdisjoint(S, U, true)
        @test isdisjoint(U, S) && isdisjoint(S, U) && disjoint1 &&
              disjoint2 && point1 == point2 == N[]

        # vertices list
        @test ispermutation(vertices_list(U),
            [N[1, 1], N[1, -1], N[-1, 1], N[-1, -1],
            N[1, 2], N[1, 0], N[2, 1], N[0, 1]])
        @test ispermutation(vertices_list(U, apply_convex_hull=true),
            [N[1, -1], N[-1, 1], N[-1, -1], N[1, 2], N[2, 1]])
    end

    # linear map
    # (only a sufficient test because there is currently no equivalence check:
    #  test that the map of the union is equivalent to the union of each map)
    A = N[1 0; 0 1]
    Y = linear_map(A, UXY)
    @test isequivalent(Y.X, linear_map(A, UXY.X))
    @test isequivalent(Y.Y, linear_map(A, UXY.Y))
    Y = linear_map(A, Uarr)
    for k in 1:length(array(Uarr))
        @test isequivalent(array(Y)[k], linear_map(A, array(Uarr)[k]))
    end

    # an_element and membership
    X = HPolyhedron([HalfSpace(N[1, 0], N(0)), HalfSpace(N[-1, 0], N(-1))])
    UXY = UnionSet(X, B1)
    @test an_element(UXY) ∈ UXY
    Uarr = UnionSetArray([EmptySet{N}(2), B1])
    @test an_element(Uarr) ∈ Uarr

    # boundedness of unbounded unions
    @test !isbounded(UXY) && !isboundedtype(typeof(UXY))
    Uarr = UnionSetArray([X, B1])
    @test !isbounded(Uarr) && !isboundedtype(typeof(Uarr))

    # is_polyhedral
    @test is_polyhedral(UXY)
    if N isa AbstractFloat
        U2 = UnionSet(B1, Ball2(zeros(N, 2), N(1)))
        @test !is_polyhedral(U2)
    end

    # volume
    X = Interval(N(1), N(2))
    Y = Interval(N(3), N(4))
    @test volume(UnionSet(X, Y)) == N(2)
end

for N in [Float64]
    B1 = BallInf(zeros(N, 2), N(1))
    B2 = Ball1(ones(N, 2), N(1))
    B3 = Hyperrectangle(low=N[-1, -1], high=N[2, 2])
    S = Singleton(N[10, 10])
    UXY = UnionSet(B1, B2)
    Uarr = UnionSetArray([B1, B2])

    for U in [UXY, Uarr]
        # intersection
        @test !isempty(intersection(U, B3)) && !isempty(intersection(B3, U))
        @test isempty(intersection(U, S)) && isempty(intersection(S, U))
    end

    # emptiness
    emptyP = HPolyhedron([HalfSpace(N[1, 0], N(0)), HalfSpace(N[-1, 0], N(-1))])
    @test !isempty(emptyP ∪ B1) && !isempty(B1 ∪ emptyP) &&
          isempty(emptyP ∪ emptyP)
    @test !isempty(UnionSetArray([emptyP, B1])) &&
          !isempty(UnionSetArray([B1, B2, emptyP])) &&
          isempty(UnionSetArray([emptyP, emptyP]))

    # boundedness
    unboundedP = HPolyhedron([HalfSpace(N[1, 0], N(0))])
    @test !isbounded(unboundedP ∪ B1) && !isbounded(B1 ∪ unboundedP)
    @test !isbounded(UnionSetArray([unboundedP, B1])) &&
          !isbounded(UnionSetArray([B1, B2, unboundedP]))

    # isdisjoint
    disjoint1, point1 = isdisjoint(UXY, UXY, true)
    disjoint2, point2 = isdisjoint(UXY, Uarr, true)
    disjoint3, point3 = isdisjoint(Uarr, UXY, true)
    disjoint4, point4 = isdisjoint(Uarr, Uarr, true)
    @test !isdisjoint(UXY, UXY) && !isdisjoint(UXY, Uarr) &&
          !isdisjoint(Uarr, UXY) && !isdisjoint(Uarr, Uarr) && !disjoint1 &&
          !disjoint2 && !disjoint3 && !disjoint4 && point1 ∈ UXY &&
          point2 ∈ UXY && point2 ∈ Uarr && point3 ∈ UXY && point3 ∈ Uarr &&
          point4 ∈ Uarr
end


text:
struct LazyArrayStyle{N} <: AbstractArrayStyle{N} end
LazyArrayStyle(::Val{N}) where N = LazyArrayStyle{N}()
LazyArrayStyle{M}(::Val{N}) where {N,M} = LazyArrayStyle{N}()
"""
    BroadcastLayout{F}()

is returned by `MemoryLayout(A)` if a matrix `A` is a `BroadcastArray`.
`F` is the typeof function that broadcast operation is applied.
"""
struct BroadcastLayout{F} <: AbstractLazyLayout end

@inline tuple_type_memorylayouts(::Type{I}) where I<:Tuple = tuple(MemoryLayout(Base.tuple_type_head(I)), tuple_type_memorylayouts(Base.tuple_type_tail(I))...)
@inline tuple_type_memorylayouts(::Type{Tuple{A}}) where {A} = (MemoryLayout(A),)
@inline tuple_type_memorylayouts(::Type{Tuple{A,B}}) where {A,B} = (MemoryLayout(A),MemoryLayout(B))
@inline tuple_type_memorylayouts(::Type{Tuple{A,B,C}}) where {A,B,C} = (MemoryLayout(A),MemoryLayout(B),MemoryLayout(C))
@inline tuple_type_memorylayouts(::Type{Tuple{A,B,C,D}}) where {A,B,C,D} = (MemoryLayout(A),MemoryLayout(B),MemoryLayout(C),MemoryLayout(D))
@inline tuple_type_memorylayouts(::Type{Tuple{A,B,C,D,E}}) where {A,B,C,D,E} = (MemoryLayout(A),MemoryLayout(B),MemoryLayout(C),MemoryLayout(D),MemoryLayout(E))

@inline broadcastlayout(::Type{F}, _...) where F = BroadcastLayout{F}()

function _copyto!(_, ::BroadcastLayout, dest::AbstractArray{<:Any,N}, bc::AbstractArray{<:Any,N}) where N
    materialize!(dest, _broadcastarray2broadcasted(bc))
    dest
end

struct BroadcastArray{T, N, F, Args} <: LazyArray{T, N}
    f::F
    args::Args
end

const BroadcastVector{T,F,Args} = BroadcastArray{T,1,F,Args}
const BroadcastMatrix{T,F,Args} = BroadcastArray{T,2,F,Args}

LazyArray(bc::Broadcasted) = BroadcastArray(bc)

BroadcastArray{T,N,F,Args}(bc::Broadcasted) where {T,N,F,Args} = BroadcastArray{T,N,F,Args}(bc.f,bc.args)
BroadcastArray{T,N}(bc::Broadcasted{Style,Axes,F,Args}) where {T,N,Style,Axes,F,Args} = BroadcastArray{T,N,F,Args}(bc.f,bc.args)
BroadcastArray{T}(bc::Broadcasted{<:Union{Nothing,BroadcastStyle},<:Tuple{Vararg{Any,N}},<:Any,<:Tuple}) where {T,N} =
    BroadcastArray{T,N}(bc)

BroadcastVector(bc::Broadcasted) = BroadcastVector{combine_eltypes(bc.f, bc.args)}(bc)
BroadcastMatrix(bc::Broadcasted) = BroadcastMatrix{combine_eltypes(bc.f, bc.args)}(bc)

@inline MemoryLayout(::Type{BroadcastArray{T,N,F,Args}}) where {T,N,F,Args} =
    broadcastlayout(F, tuple_type_memorylayouts(Args)...)

_broadcast2broadcastarray() = ()
_broadcast2broadcastarray(a, b...) = tuple(a, _broadcast2broadcastarray(b...)...)
_broadcast2broadcastarray(a::Broadcasted{DefaultArrayStyle{0}}, b...) = tuple(materialize(a), _broadcast2broadcastarray(b...)...)
_broadcast2broadcastarray(a::Broadcasted, b...) = tuple(BroadcastArray(a), _broadcast2broadcastarray(b...)...)

_BroadcastArray(bc::Broadcasted) = BroadcastArray{combine_eltypes(bc.f, bc.args)}(bc)
BroadcastArray(bc::Broadcasted{S}) where S =
    _BroadcastArray(instantiate(Broadcasted{S}(bc.f, _broadcast2broadcastarray(bc.args...))))

BroadcastArray(f, A, As...) = BroadcastArray(broadcasted(f, A, As...))
BroadcastArray{T}(f, A, As...) where T = BroadcastArray{T}(instantiate(broadcasted(f, A, As...)))
BroadcastMatrix(f, A...) = BroadcastMatrix(broadcasted(f, A...))
BroadcastVector(f, A...) = BroadcastVector(broadcasted(f, A...))

BroadcastArray{T,N}(f, A...) where {T,N} = BroadcastArray{T,N,typeof(f),typeof(A)}(f, A)

BroadcastArray(b::BroadcastArray) = b
BroadcastVector(A::BroadcastVector) = A
BroadcastMatrix(A::BroadcastMatrix) = A

@inline __broadcastarray2broadcasted() = ()
@inline __broadcastarray2broadcasted(a, b...) = tuple(_broadcastarray2broadcasted(a), __broadcastarray2broadcasted(b...)...)
@inline _broadcastarray2broadcasted(lay::BroadcastLayout, a) = broadcasted(call(lay, a), __broadcastarray2broadcasted(arguments(lay, a)...)...)
@inline _broadcastarray2broadcasted(lay::BroadcastLayout, a::BroadcastArray) = broadcasted(call(lay, a), __broadcastarray2broadcasted(arguments(lay, a)...)...)
@inline _broadcastarray2broadcasted(_, a) = a
@inline _broadcastarray2broadcasted(lay, a::BroadcastArray) = error("Overload LazyArrays._broadcastarray2broadcasted(::$(lay), _)")
@inline _broadcastarray2broadcasted(::DualLayout{ML}, a) where ML = _broadcastarray2broadcasted(ML(), a)
@inline _broadcastarray2broadcasted(a) = _broadcastarray2broadcasted(MemoryLayout(a), a)
@inline _broadcasted(A) = instantiate(_broadcastarray2broadcasted(A))
broadcasted(A::BroadcastArray) = _broadcasted(A)
broadcasted(A::SubArray{<:Any,N,<:BroadcastArray}) where N = _broadcasted(A)
Broadcasted(A::BroadcastArray) = broadcasted(A)::Broadcasted
Broadcasted(A::SubArray{<:Any,N,<:BroadcastArray}) where N = broadcasted(A)::Broadcasted

@inline BroadcastArray(A::AbstractArray) = BroadcastArray(call(A), arguments(A)...)

axes(A::BroadcastArray) = axes(broadcasted(A))
size(A::BroadcastArray) = map(length, axes(A))

@propagate_inbounds getindex(A::BroadcastArray{T,N}, kj::Vararg{Int,N}) where {T,N} = convert(T,broadcasted(A)[kj...])::T

converteltype(::Type{T}, A::AbstractArray) where T = convert(AbstractArray{T}, A)
converteltype(::Type{T}, A) where T = convert(T, A)
sub_materialize(::BroadcastLayout, A) = converteltype(eltype(A), materialize(_broadcasted(A)))

copy(bc::Broadcasted{<:LazyArrayStyle}) = BroadcastArray(bc)

# BroadcastArray are immutable
copy(bc::BroadcastArray) = bc
map(::typeof(copy), bc::BroadcastArray) = bc
copy(bc::AdjOrTrans{<:Any,<:BroadcastArray}) = bc

# Replacement for #18.
# Could extend this to other similar reductions in Base... or apply at lower level?
# for (fname, op) in [(:sum, :add_sum), (:prod, :mul_prod),
#                     (:maximum, :max), (:minimum, :min),
#                     (:all, :&),       (:any, :|)]
function Base._sum(f, A::BroadcastArray, ::Colon)
    bc = broadcasted(A)
    T = Broadcast.combine_eltypes(f ∘ bc.f, bc.args)
    out = zero(T)
    @simd for I in eachindex(bc)
        @inbounds out += f(bc[I])
    end
    out
end
function Base._prod(f, A::BroadcastArray, ::Colon)
    bc = broadcasted(A)
    T = Broadcast.combine_eltypes(f ∘ bc.f, bc.args)
    out = one(T)
    @simd for I in eachindex(bc)
        @inbounds out *= f(bc[I])
    end
    out
end

BroadcastStyle(::Type{<:LazyArray{<:Any,N}}) where N = LazyArrayStyle{N}()
BroadcastStyle(::Type{<:Adjoint{<:Any,<:LazyVector{<:Any}}}) where N = LazyArrayStyle{2}()
BroadcastStyle(::Type{<:Transpose{<:Any,<:LazyVector{<:Any}}}) where N = LazyArrayStyle{2}()
BroadcastStyle(::Type{<:Adjoint{<:Any,<:LazyMatrix{<:Any}}}) where N = LazyArrayStyle{2}()
BroadcastStyle(::Type{<:Transpose{<:Any,<:LazyMatrix{<:Any}}}) where N = LazyArrayStyle{2}()
BroadcastStyle(::Type{<:SubArray{<:Any,1,<:LazyMatrix,<:Tuple{Slice,Any}}}) = LazyArrayStyle{1}()
BroadcastStyle(L::LazyArrayStyle{N}, ::StaticArrayStyle{N}) where N = L
BroadcastStyle(L::LazyArrayStyle{N}, ::StructuredMatrixStyle)  where N = L

## scalar-range broadcast operations ##
# Ranges already support smart broadcasting
for op in (+, -, big)
    @eval begin
        broadcasted(::LazyArrayStyle{1}, ::typeof($op), r::AbstractRange) =
            broadcast(DefaultArrayStyle{1}(), $op, r)
    end
end

for op in (-, +, *, /)
    @eval broadcasted(::LazyArrayStyle{1}, ::typeof($op), r::AbstractRange, x::Real) = broadcast(DefaultArrayStyle{1}(), $op, r, x)
end

for op in (-, +, *, \)
    @eval broadcasted(::LazyArrayStyle{1}, ::typeof($op), x::Real, r::AbstractRange) = broadcast(DefaultArrayStyle{1}(), $op, x, r)
end

broadcasted(::LazyArrayStyle{N}, op, r::AbstractFill{T,N}) where {T,N} = broadcast(DefaultArrayStyle{N}(), op, r)
broadcasted(::LazyArrayStyle{N}, op, r::AbstractFill{T,N}, x::Number) where {T,N} = broadcast(DefaultArrayStyle{N}(), op, r, x)
broadcasted(::LazyArrayStyle{N}, op, x::Number, r::AbstractFill{T,N}) where {T,N} = broadcast(DefaultArrayStyle{N}(), op, x, r)
broadcasted(::LazyArrayStyle{N}, op, r::AbstractFill{T,N}, x::Ref) where {T,N} = broadcast(DefaultArrayStyle{N}(), op, r, x)
broadcasted(::LazyArrayStyle{N}, op, x::Ref, r::AbstractFill{T,N}) where {T,N} = broadcast(DefaultArrayStyle{N}(), op, x, r)
broadcasted(::LazyArrayStyle{N}, op, r1::AbstractFill{T,N}, r2::AbstractFill{V,N}) where {T,V,N} = broadcast(DefaultArrayStyle{N}(), op, r1, r2)
broadcasted(::LazyArrayStyle{1}, ::typeof(*), a::AbstractFill, b::AbstractRange) = broadcast(DefaultArrayStyle{1}(), *, a, b)
broadcasted(::LazyArrayStyle{1}, ::typeof(*), a::AbstractRange, b::AbstractFill) = broadcast(DefaultArrayStyle{1}(), *, a, b)
broadcasted(::LazyArrayStyle{1}, ::typeof(*), a::Zeros{<:Any,1}, b::AbstractRange) = broadcast(DefaultArrayStyle{1}(), *, a, b)
broadcasted(::LazyArrayStyle{1}, ::typeof(*), a::AbstractRange, b::Zeros{<:Any,1}) = broadcast(DefaultArrayStyle{1}(), *, a, b)

for op in (:*, :/, :\)
    @eval broadcasted(::LazyArrayStyle{N}, ::typeof($op), a::Zeros{T,N}, b::Zeros{V,N}) where {T,V,N} = broadcast(DefaultArrayStyle{N}(), $op, a, b)
end

for op in (:*, :/)
    @eval begin
        broadcasted(::LazyArrayStyle{N}, ::typeof($op), a::Zeros{T,N}, b::AbstractArray{V,N}) where {T,V,N} = broadcast(DefaultArrayStyle{N}(), $op, a, b)
        broadcasted(::LazyArrayStyle{N}, ::typeof($op), a::Zeros{T,N}, b::Broadcasted) where {T,N} = broadcast(DefaultArrayStyle{N}(), $op, a, b)
    end
end
for op in (:*, :\)
    @eval begin
        broadcasted(::LazyArrayStyle{N}, ::typeof($op), a::AbstractArray{T,N}, b::Zeros{V,N}) where {T,V,N} = broadcast(DefaultArrayStyle{N}(), $op, a, b)
        broadcasted(::LazyArrayStyle{N}, ::typeof($op), a::Broadcasted, b::Zeros{V,N}) where {V,N} = broadcast(DefaultArrayStyle{N}(), $op, a, b)
    end
end

###
# support
###

_broadcast_colsupport(ax, ::Tuple{}, A, j) = ax[1]
_broadcast_colsupport(ax, ::Tuple{<:Any}, A, j) = colsupport(A,j)
_broadcast_colsupport(ax, Aax::Tuple{OneTo{Int},<:Any}, A, j) = length(Aax[1]) == 1 ? ax[1] : colsupport(A,j)
_broadcast_colsupport(ax, ::Tuple{<:Any,<:Any}, A, j) = colsupport(A,j)
_broadcast_rowsupport(ax, ::Tuple{}, A, j) = ax[2]
_broadcast_rowsupport(ax, ::Tuple{<:Any}, A, j) = ax[2]
_broadcast_rowsupport(ax, Aax::Tuple{<:Any,OneTo{Int}}, A, j) = length(Aax[2]) == 1 ? ax[2] : rowsupport(A,j)
_broadcast_rowsupport(ax, ::Tuple{<:Any,<:Any}, A, j) = rowsupport(A,j)

colsupport(lay::BroadcastLayout{typeof(*)}, A, j) = intersect(_broadcast_colsupport.(Ref(axes(A)), axes.(arguments(lay,A)), arguments(lay,A), Ref(j))...)
rowsupport(lay::BroadcastLayout{typeof(*)}, A, j) = intersect(_broadcast_rowsupport.(Ref(axes(A)), axes.(arguments(lay,A)), arguments(lay,A), Ref(j))...)

for op in (:+, :-)
    @eval begin
        rowsupport(lay::BroadcastLayout{typeof($op)}, A, j) = convexunion(_broadcast_rowsupport.(Ref(axes(A)), axes.(arguments(lay,A)), arguments(lay,A), Ref(j))...)
        colsupport(lay::BroadcastLayout{typeof($op)}, A, j) = convexunion(_broadcast_colsupport.(Ref(axes(A)), axes.(arguments(lay,A)), arguments(lay,A), Ref(j))...)
    end
end

###
# SubArray
###

# TODO: special case adjtrans to skip the `isone` check and return numbers instead of 1-vectors.
#

sublayout(b::BroadcastLayout, _) = b

@inline _broadcastviewinds(::Tuple{}, inds) = ()
@inline _broadcastviewinds(ax::Tuple{OneTo{Int},Vararg{Any}}, inds::Tuple{Number,Vararg{Any}}) =
    tuple(isone(length(ax[1])) ? 1 : convert(Int,inds[1]), _broadcastviewinds(tail(ax), tail(inds))...)
@inline _broadcastviewinds(ax::Tuple{OneTo{Int},Vararg{Any}}, inds::Tuple{AbstractVector{<:Integer},Vararg{Any}}) =
    tuple(isone(length(ax[1])) ? convert(typeof(inds[1]),Base.OneTo(1)) : inds[1], _broadcastviewinds(tail(ax), tail(inds))...)
@inline function _broadcastviewinds(ax::Tuple{OneTo{Int},Vararg{Any}}, inds::Tuple{Any,Vararg{Any}})
    @assert isone(length(ax[1]))
    tuple(Base.OneTo(1), _broadcastviewinds(tail(ax), tail(inds))...)
end

@inline _broadcastviewinds(ax, inds) = # don't support special broacasting
    tuple(inds[1], _broadcastviewinds(tail(ax), tail(inds))...)

_viewifmutable(a, inds::Number...) = a[inds...]
@inline _viewifmutable(a, inds...) = view(a, inds...)
@inline _viewifmutable(a::AbstractFill, inds...) = a[inds...]
@inline _viewifmutable(a::AbstractRange, inds...) = a[inds...]
@inline _viewifmutable(a::AbstractRange, inds::Number...) = a[inds...]
# _viewifmutable(a::BroadcastArray, inds...) = a[inds...]
_viewifmutable(a::AdjOrTrans{<:Any,<:AbstractVector}, k::Integer, j::Integer) = a[k,j]
function _viewifmutable(a::AdjOrTrans{<:Any,<:AbstractVector}, k::Integer, j)
    @assert k == 1
    _viewifmutable(parent(a), j)
end
@inline _broadcastview(a, inds) = _viewifmutable(a, _broadcastviewinds(axes(a), inds)...)
@inline _broadcastview(a::Number, inds) = a
@inline _broadcastview(a::Base.RefValue, inds) = a

@inline __broadcastview(inds) = ()
@inline __broadcastview(inds, a, b...) = (_broadcastview(a, inds), __broadcastview(inds, b...)...)

@inline function _broadcast_sub_arguments(lay, P, V)
    args = arguments(lay, P)
    __broadcastview(parentindices(V), args...)
end

@inline _broadcast_sub_arguments(lay::DualLayout{ML}, P, V::AbstractVector) where ML =
    arguments(ML(), view(P', parentindices(V)[2]))

@inline _broadcast_sub_arguments(A, V) = _broadcast_sub_arguments(MemoryLayout(A), A, V)
@inline _broadcast_sub_arguments(V) =  _broadcast_sub_arguments(parent(V), V)
@inline arguments(b::BroadcastLayout, V::SubArray) = _broadcast_sub_arguments(V)
@inline call(b::BroadcastLayout, a::SubArray) = call(b, parent(a))

###
# Transpose
###

call(b::BroadcastLayout, a::AdjOrTrans) = call(b, parent(a))

transposelayout(b::BroadcastLayout) = b

_adjoint(a) = adjoint(a)
_adjoint(a::Ref) = a
_transpose(a) = transpose(a)
_transpose(a::Ref) = a

arguments(b::BroadcastLayout, A::Adjoint) = map(_adjoint, arguments(b, parent(A)))
arguments(b::BroadcastLayout, A::Transpose) = map(_transpose, arguments(b, parent(A)))

###
# Show
###

_broadcastarray_summary(io, A) = _broadcastarray_summary(io, A.f, arguments(A)...)
function _broadcastarray_summary(io, f, args...)
    print(io, "$f.(")
    summary(io, first(args))
    for a in tail(args)
        print(io, ", ")
        summary(io, a)
    end
    print(io, ")")
end

for op in (:+, :-, :*, :\, :/)
    @eval begin
        function _broadcastarray_summary(io::IO, ::typeof($op), args...)
            if length(args) == 1
                print(io, "($($op)).(")
                summary(io, first(args))
                print(io, ")")
            else
                print(io, "(")
                summary(io, first(args))
                print(io, ")")
                for a in tail(args)
                    print(io, " .$($op) (")
                    summary(io, a)
                    print(io, ")")
                end
            end
        end
    end
end

function _broadcastarray_summary(io::IO, ::typeof(Base.literal_pow), ::Base.RefValue{typeof(^)}, x, ::Base.RefValue{Val{K}}) where {N,K}
    print(io, "(")
    summary(io, x)
    print(io, ") .^ $K")
end

function _broadcastarray_summary(io::IO, ::typeof(^), x, y)
    print(io, "(")
    summary(io, x)
    print(io, ") .^ ")
    summary(io, y)
end

Base.array_summary(io::IO, C::BroadcastArray, inds::Tuple{Vararg{OneTo}}) = _broadcastarray_summary(io, C)
function Base.array_summary(io::IO, C::BroadcastArray, inds)
    _broadcastarray_summary(io, C)
    print(io, " with indices ", Base.inds2string(inds))
end

function Base.array_summary(io::IO, C::Adjoint{<:Any,<:LazyArray}, inds::Tuple{Vararg{OneTo}})
    print(io, "(")
    summary(io, parent(C))
    print(io, ")'")
end
function Base.array_summary(io::IO, C::Adjoint{<:Any,<:LazyArray}, inds)
    print(io, "(")
    summary(io, parent(C))
    print(io, ")' with indices ", Base.inds2string(inds))
end

function Base.array_summary(io::IO, C::Transpose{<:Any,<:LazyArray}, inds::Tuple{Vararg{OneTo}})
    print(io, "transpose(")
    summary(io, parent(C))
    print(io, ")")
end
function Base.array_summary(io::IO, C::Transpose{<:Any,<:LazyArray}, inds)
    print(io, "transpose(")
    summary(io, parent(C))
    print(io, ") with indices ", Base.inds2string(inds))
end

###
# Mul
###

_broadcast_mul_mul(A, B) = simplify(Mul(broadcast(*, A...), B))
_broadcast_mul_mul((a,B)::Tuple{AbstractVector,AbstractMatrix}, C) = a .* (B*C)
_broadcast_mul_mul((A,b)::Tuple{AbstractMatrix,AbstractVector}, C) = b .* (A*C)
@inline copy(M::Mul{BroadcastLayout{typeof(*)}}) = _broadcast_mul_mul(arguments(BroadcastLayout{typeof(*)}(), M.A), M.B)
@inline copy(M::Mul{BroadcastLayout{typeof(*)},<:AbstractLazyLayout}) = _broadcast_mul_mul(arguments(BroadcastLayout{typeof(*)}(), M.A), M.B)
@inline copy(M::Mul{BroadcastLayout{typeof(*)},ApplyLayout{typeof(*)}}) = _broadcast_mul_mul(arguments(BroadcastLayout{typeof(*)}(), M.A), M.B)

getindex(A::BroadcastMatrix{<:Any,typeof(*),<:Tuple{AbstractVector,AbstractMatrix}}, ::Colon, j::Integer) = A.args[1] .* A.args[2][:,j]
getindex(A::BroadcastMatrix{<:Any,typeof(*),<:Tuple{AbstractMatrix,AbstractVector}}, ::Colon, j::Integer) = A.args[1][:,j] .* A.args[2]


text:
# Matern 5/2 ARD covariance Function

@doc """
# Description
Constructor for the ARD Matern 5/2 kernel (covariance)

k(x,x') = σ²(1+√3*d/L + 5d²/3L²)exp(-√5*d/L), where d = |x-x'| and L = diag(ℓ₁,ℓ₂,...)
# Arguments:
* `ll::Vector{Float64}`: Log of the length scale ℓ
* `lσ::Float64`: Log of the signal standard deviation σ
""" ->
type Mat52Ard <: StationaryARD
    ℓ2::Vector{Float64}   # Log of Length scale
    σ2::Float64           # Log of signal std
    Mat52Ard(ll::Vector{Float64}, lσ::Float64) = new(exp(2.0*ll), exp(2.0*lσ))
end

function set_params!(mat::Mat52Ard, hyp::Vector{Float64})
    length(hyp) == num_params(mat) || throw(ArgumentError("Mat52 kernel only has $(num_params(mat)) parameters"))
    d = length(mat.ℓ2)
    mat.ℓ2 = exp(2.0*hyp[1:d])
    mat.σ2 = exp(2.0*hyp[d+1])
end

get_params(mat::Mat52Ard) = [log(mat.ℓ2)/2.0; log(mat.σ2)/2.0]
get_param_names(mat::Mat52Ard) = [get_param_names(mat.ℓ2, :ll); :lσ]
num_params(mat::Mat52Ard) = length(mat.ℓ2) + 1

metric(mat::Mat52Ard) = WeightedEuclidean(1.0./(mat.ℓ2))
cov(mat::Mat52Ard, r::Float64) = mat.σ2*(1+sqrt(5)*r+5/3*r^2)*exp(-sqrt(5)*r)

function grad_kern(mat::Mat52Ard, x::Vector{Float64}, y::Vector{Float64})
    #r = distance(mat,x,y)
    wdiff = (x-y).^2./mat.ℓ2
    r = sqrt(sum(wdiff))
    exp_r = exp(-sqrt(5)*r)

    g1 = mat.σ2*(5/3)*(1+sqrt(5)*r)*exp_r.*wdiff # dK_d(log ℓ)
    g2 = 2.0*mat.σ2*(1+sqrt(5)*r+(5/3)*r^2)*exp_r  # dK_d(log σ)

    return [g1; g2]
end

function grad_stack!(stack::AbstractArray, mat::Mat52Ard, X::Matrix{Float64}, data::StationaryARDData)
    d = size(X,1)
    R = distance(mat,X)
    exp_R = exp(-sqrt(5)*R)

    stack[:,:,d+1] = 2.0*cov(mat, X)
    part = (5/3) * mat.σ2 .* exp_R .* (1.0 + sqrt(5)*R)
    broadcast!(*, view(stack, :, :, 1:d), (5/3)*mat.σ2, 1 + sqrt(5)*R, exp_R, data.dist_stack, reshape(1.0./mat.ℓ2, (1,1,d)))

    return stack
end


text:
using Pkg
Pkg.add("Coverage")
using Coverage
cov_res = process_folder()
haskey(ENV, "COVERALLS_URL") && Coveralls.submit(cov_res)


text:
export ARC_SecA_ls
"""
A one dimensionnal ARC algorithm using the exact second derivative in
the second order Taylor expansion enriched with a cubic term.
For more documentation referer to ARC_generic_ls.
"""
function ARC_SecA_ls(h :: LineModel, stop_ls :: LS_Stopping; kwargs...)
    ls_meta = LS_Function_Meta(dir = :Sec)
    (state, stop_ls.meta.optimal) = ARC_generic_ls(h, stop_ls, f_meta = ls_meta; kwargs...)

    return (state, stop_ls.meta.optimal)
end


text:
using Plots, LaTeXStrings
using DelimitedFiles
using Revise
using JLD2, FileIO

cd("julia/")
push!(LOAD_PATH, ".")
using VirusLoadCurve

tdata = [4.0/24.0, 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.]
RAWDATA = readdlm("../../data/Smith2018/Virus_Best10.txt", ' ', Float64, '\n')

k = 6
ttdata = repeat(tdata, inner=1)
tend = ttdata[end]
vvdata = reshape(RAWDATA[k, :], 1, :)[:]

data = VirusLoadCurve.VirusLoadData(ttdata, vvdata)

par0 = [0.70, 2.88, 6.00, 7.60, 0.20, 5.0]
resultVLF = VirusLoadCurve.fitVLF(data; ϵ=0.1)
summary(resultVLF)
save("test.jld2", Dict("result" => resultVLF, "data" => data))

resultVLF = load("test.jld2", "result")

plot(data)
plot!(resultVLF, empirical=true, ylims=(-2, 9))
ylims!((-2, 9))

pt = VirusLoadCurve.Boxplots(resultVLF)


text:
module LineExtraction

using Color
using Images
using ImageView

include("extraction.jl")

export extract

end # module


text:
module PlotlyLight

using Random
using JSON3
using EasyConfig
using Cobweb
using Artifacts

export Plot, Config, collectrows

const cdn_url = "https://cdn.plot.ly/plotly-2.11.0.min.js"
const plotlyjs = joinpath(artifact"plotly.min.js", basename(cdn_url))
const templates_dir = artifact"plotly_templates"
const templates = map(x -> replace(x, ".json" => ""), readdir(templates_dir))

#-----------------------------------------------------------------------------# defaults
module Defaults
using EasyConfig: Config
using JSON3
export src, class, style, parent_class, parent_style, config, layout

src             = Ref(:cdn)
class           = Ref("")
style           = Ref("height: 100%;")
parent_class    = Ref("")
parent_style    = Ref("height: 100vh;")
config          = Ref(Config(displaylogo=false, responsive=true))
layout          = Ref(Config())

function reset!()
    src[]           = :cdn
    class[]         = ""
    style[]         = "height: 100%;"
    parent_class[]  = ""
    parent_style[]  = "height: 100vh;"
    config[]        = Config(displaylogo=false, responsive=true)
    layout[]        = Config()
end
end

#-----------------------------------------------------------------------------# src!
src_opts = [:cdn, :local, :standalone, :none]
"""
    src!(x::Symbol) # `x` must be one of: $src_opts

- `:cdn` → Use PlotlyJS CDN.
- `:local` → Use local artifact.
- `:standalone` → Write JS into the HTML file directly (can be shared and viewed offline).
- `:none` → For when inserting into a page with Plotly.js already included.
"""
src!(x::Symbol) = (x in src_opts || error("src must be one of: $src_opts"); Defaults.src[] = x)

#-----------------------------------------------------------------------------# template!
"""
    template(t)

Load the template `t`, which must be one of:

```
$(join(templates, "\n"))
```
"""
function template(t)
    string(t) in templates || error("$t not found.  Options are one of: $(join(templates, ", ")).")
    open(io -> JSON3.read(io, Config), joinpath(templates_dir, string(t) * ".json"))
end

"""
    template!(t)

Replace the `template` key of the default layout with `t`.  See also : `PlotlyLight.template`.
"""
function template!(t)
    Defaults.layout[].template = template(t)
end

#-----------------------------------------------------------------------------# Plot
"""
    Plot(data, layout, config; id, js)

- A Plotly.js plot with components `data`, `layout`, and `config`.
  - `data = Config()`: A `Config` (single trace) or `Vector{Config}` (multiple traces).
  - `layout = Config()`.
  - `config = Config(displaylogo=false, responsive=true)`.
- Each of the three components are converted to JSON via `JSON3.write`.
- See the Plotly Javascript docs here: https://plotly.com/javascript/.
- Keyword Args:
  - `id`: The `id` of the `<div>` the plot will be created in.  Default: `randstring(10)`.
  - `js`:  `Cobweb.Javascript` to add after the creation of the plot.  Default:
    - `Cobweb.Javascript("console.log('plot created!')")`

### Example

    p = Plot(Config(x=1:10, y=randn(10)))
    p.layout.title.text = "My Title!"
    p
"""
Base.@kwdef mutable struct Plot
    data::Vector{Config}    = Config[]
    layout::Config          = Defaults.layout[]
    config::Config          = Defaults.config[]
    id::String              = randstring(10)            # id of graphDiv
    js::Cobweb.Javascript   = Cobweb.Javascript("console.log('plot created!')")
end
function Plot(traces, layout=Defaults.layout[], config=Defaults.config[]; kw...)
    data = traces isa Config ? [traces] : traces
    Plot(; kw..., data,
        layout = merge(layout, Defaults.layout[]),
        config = merge(config, Defaults.config[])
    )
end

#-----------------------------------------------------------------------------# Display
Base.display(::Cobweb.CobwebDisplay, o::Plot) = display(Cobweb.CobwebDisplay(), Cobweb.Page(o))

function Base.show(io::IO, M::MIME"text/html", o::Plot)
    class, style = Defaults.class, Defaults.style
    parent_class, parent_style = Defaults.parent_class, Defaults.parent_style
    parent_style = if get(io, :is_pluto, false)
        s = replace(parent_style[], r"height.*;" => "")
        "height: 400px;" * s
    else
        parent_style[]
    end
    src = Defaults.src[]
    src in [:cdn, :standalone, :none, :local] || error("`src` must be :cdn, :standalone, :none, or :local")
    println(io, "
")
    println(io, "

")
    println(io, "
")

    if src === :cdn
        println(io, "")
    elseif src === :standalone
        print(io, "")
    elseif src === :local
        println(io, "")
    else
        # :none
    end

    println(io, "\n")
end

#-----------------------------------------------------------------------------# vecvec
collectrows(x::AbstractMatrix) = collect.(eachrow(x))

end # module


text:
using JLD2, FileIO

data = zeros(Int8, 2, 2)
data[2,2] =1

save("file.jld2","data", data)
#
uu = load("file.jld2", "data")
uu


text:
"""
ULMFiT - Sentiment Analyzer
"""

struct SentimentAnalyzer
    vocab::Vector
    rnn_layers::Flux.Chain
    linear_layers::Flux.Chain
end

function SentimentAnalyzer(weights=load_weights(), clsfr_hidden_sz::Integer=50)
    inRNN = size(weights[2], 2)
    hidRNN = size(weights[7], 2)
    outRNN = size(weights[12], 1)/4
    clsfr_hidden_sz = size(weights[end-1], 2)/3
    clsfr_out_sz = size(weights[end-1], 1)
    vocab = intern.(string.(readdlm("vocab.csv",',', header=false)[:, 1]))
    sa = SentimentAnalyzer(
        vocab,
        Chain(
            DroppedEmbeddings(size(weights[1])...),
            AWD_LSTM(inRNN, hidRNN),
            AWD_LSTM(hidRNN, hidRNN),
            AWD_LSTM(hidRNN, outRNN)
        ),
        Chain(
            PooledDense(outRNN, clsfr_hidden_sz, relu)),
            BatchNorm(clsfr_hidden_sz, relu),
            Dense(clsfr_hidden_sz, clsfr_out_sz, sigmoid),
        )
    )
    test_mode!(sa)
    return sa
end

function get_sentiment(sa::SentimentAnlayzer, text::String)
    tokens = tokenize(text)
    idxs = map(x -> indices([x], sa.vocab, "_unk_"), tokens)
    H = sa.rnn_layers(idxs)
    sentiment_score = sa.linear_layers(H)
    return sentiment_score
end


text:
@testset "Testing memcpy" begin
    # Test error checking first.
    x = Vector{Vector{Int}}()
    y = Vector{Vector{Int}}()

    @test_throws ArgumentError CachedArrays.memcpy!(x, y)

    # Source/Dest length mismatch
    x = rand(Float32, 1024)
    y = rand(Float32, 1025)
    @test_throws ArgumentError CachedArrays.memcpy!(x, y)

    # Test alignment
    popfirst!(y)
    #@test_throws ArgumentError CachedArrays.memcpy!(x, y)
    #@test_throws ArgumentError CachedArrays.memcpy!(y, x)

    # Now, see if the memcpy works correctly.
    x = rand(Float32, 1024 * 10)
    y = rand(Float32, 1024 * 10)
    @test x != y
    CachedArrays.memcpy!(x, y)
    @test x == y

    # Make this kind of a wierd number of bytes.
    x = rand(Float32, 1024 * 10 + 5)
    y = rand(Float32, 1024 * 10 + 5)
    @test x != y
    CachedArrays.memcpy!(x, y)
    @test x == y
end


text:
using Base.Test

println("testing MPI...")
run(`julia mpi_test.jl`)
run(`mpirun -np 2 julia mpi_test.jl`)
println("done testing MPI.")

println("testing data source/sink...")
include("source_test.jl")
include("sink_test.jl")
run(`mpirun -np 2 julia sink_test.jl`)
println("done testing data source/sink.")

println("testing filter...")
include("filter_test1.jl")
println("done testing filter.")

println("testing join...")
include("join_test1.jl")
println("done testing join.")

println("testing aggregate...")
include("aggregate_test1.jl")
include("aggregate_rename_test.jl")
println("done testing aggregate...")

println("testing push filter up...")
include("push_filter_test1.jl")
println("done testing push filter up...")

include("q26_test.jl")
include("q05_test.jl")
# run manually since boost causes travis error
#include("length_unique_test.jl")
include("tables_cat.jl")
include("tables_transpose_hcat.jl")
include("stencil_test1.jl")
run(`mpirun -np 2 julia stencil_test2.jl`)
run(`mpirun -np 2 julia cumsum_test.jl`)
include("index_test.jl")

# commented out since MKL BLACS and SCALAPACK libraries are not linked automatically currently
#run(`mpirun -np 4 julia --depwarn=no 2d_test.jl`)
# mpiicpc -O0 -g -std=c++11  -DMKL_ILP64 -qopenmp -I${MKLROOT}/include -I/usr/local/hdf5/include  -g -fpic  -o main1 main1.cc -L/usr/local/hdf5/lib -lhdf5  -lm  -Wl,--start-group -L${MKLROOT}/lib/intel64 -lmkl_scalapack_ilp64 -lmkl_intel_ilp64 -lmkl_core -lmkl_intel_thread -lmkl_blacs_intelmpi_ilp64 -lpthread -lm -ldl -Wl,--end-group
# mpirun -np 4 ./main1


text:
@test length(Dojo.module_dir()) > 0
@test Dojo.scn(1.0) == " 1.0e+0"
@test Dojo.scn(0.0) == " 0.0e+0"
@test Dojo.scn(Inf) == " Inf"
@test Dojo.scn(123.1,
    digits=0) == " 1e+2"

# show functions
test_origin = Origin{Float64}()
display(test_origin)
@test true

test_body = Box(1.0,2.0,3.0,4.0)
display(test_body)
display(test_body.shape)
display(test_body.state)
@test true

test_joint = JointConstraint(Fixed(test_origin,test_body))
display(test_joint)
display(test_joint.translational)
display(test_joint.rotational)
@test true

test_mechanism = Mechanism(test_origin,[test_body],[test_joint])
display(test_mechanism)
display(test_mechanism.matrix_entries[1])
display(test_mechanism.system)
@test true

test_storage = Storage(3,4)
display(test_storage)
@test true

test_options = SolverOptions()
display(test_options)
@test true


text:
using Genie.Renderer.Json
import JSON3
import HTTP

yq = randn(3)

z = json(yq)

h = String(z.body)
out_j3 = JSON3.read(h)
out = convert(Vector{Float64}, out_j3)

max_stars_repo_name: RoyCCWang/WebAppRKHS

max_issues_repo_name: RoyCCWang/WebAppRKHS

max_forks_repo_name: RoyCCWang/WebAppRKHS


text:
# -*- coding: utf-8 -*-
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.10.3
#   kernelspec:
#     display_name: Julia 1.7.0
#     language: julia
#     name: julia-1.7
# ---

# %%
using BenchmarkTools

# %%
using Random

function mcpi(N, rng = MersenneTwister())
    c = 0
    for _ in 1:N
        c += rand(rng)^2 + rand(rng)^2 ≤ 1
    end
    4c/N
end

println("Julia v", VERSION)
print("mcpi(10^8):")
@btime mcpi(10^8)

# %%
using Base.Threads
using Random

function mcpi_atomic(N, rng = MersenneTwister())
    a = Atomic{Int}(0)
    @threads for _ in 1:N
        atomic_add!(a, Int(rand(rng)^2 + rand(rng)^2 ≤ 1))
    end
    4a[]/N
end

println("Julia v", VERSION)
@show Threads.nthreads()
print("mcpi_atomic(10^8):")
@btime mcpi_atomic(10^8)

# %% [markdown]
# https://gist.github.com/TsuMakoto/4138d3c2fd05a400d02eff0b91e34658

# %%
using Base.Threads
using Random

mutable struct MCAtomic{T}; @atomic n::T; end

function mcpi_tm(N, rng = MersenneTwister())
    mc = MCAtomic(0)
    @threads for _ in 1:N
        @atomic mc.n += Int(rand(rng)^2 + rand(rng)^2 ≤ 1)
    end
    4mc.n/N
end

println("Julia v", VERSION)
@show Threads.nthreads()
print("mcpi_tm(10^8):")
@btime mcpi_tm(10^8)

# %%
using Base.Threads
using Distributed: splitrange
using Random

function mcpi_splittange_atomic(N)
    ranges = splitrange(1, N, nthreads())
    a = Atomic{Int}(0)
    @threads for ran in ranges
        rng = MersenneTwister()
        c = 0
        for _ in ran
            c += rand(rng)^2 + rand(rng)^2 ≤ 1
        end
        atomic_add!(a, c)
    end
    4a[]/N
end

println("Julia v", VERSION)
@show Threads.nthreads()
print("mcpi_splittange_atomic(10^8):")
@btime mcpi_splittange_atomic(10^8)

# %%
# The following code is a modified version of
#
#    function _threadsfor(iter, lbody, schedule)
#    macro threads(args...)
#
# in https://github.com/JuliaLang/julia/blob/9f3265399227fbfc4f0160ec3592a9262bd3eb5f/base/threadingconstructs.jl
#
# Its license is MIT: https://julialang.org/license

using Base.Threads
using Base.Threads: threading_run

function _my_threadsfor(iter, lbody, prebody, postbody, schedule)
    lidx = iter.args[1]         # index
    range = iter.args[2]
    quote
        local threadsfor_fun
        let range = $(esc(range))
        function threadsfor_fun(onethread=false)
            r = range # Load into local variable
            lenr = length(r)
            # divide loop iterations among threads
            if onethread
                tid = 1
                len, rem = lenr, 0
            else
                tid = threadid()
                len, rem = divrem(lenr, nthreads())
            end
            # not enough iterations for all the threads?
            if len == 0
                if tid > rem
                    return
                end
                len, rem = 1, 0
            end
            # compute this thread's iterations
            f = firstindex(r) + ((tid-1) * len)
            l = f + len - 1
            # distribute remaining iterations evenly
            if rem > 0
                if tid <= rem
                    f = f + (tid-1)
                    l = l + tid
                else
                    f = f + rem
                    l = l + rem
                end
            end
            # run prebody
            $(esc(prebody))
            # run this thread's iterations
            for i = f:l
                local $(esc(lidx)) = @inbounds r[i]
                $(esc(lbody))
            end
            # run postbody
            $(esc(postbody))
        end
        end
        if threadid() != 1 || ccall(:jl_in_threaded_region, Cint, ()) != 0
            $(if schedule === :static
              :(error("`@my_threads :static` can only be used from thread 1 and not nested"))
              else
              # only use threads when called from thread 1, outside @threads
              :(Base.invokelatest(threadsfor_fun, true))
              end)
        else
            threading_run(threadsfor_fun)
        end
        nothing
    end
end

"""
    @my_threads
A macro to parallelize a `for` loop to run with multiple threads.
It splits the iteration space among multiple tasks with `prebody` and `postbody`.
It runs those tasks on threads according to a scheduling policy.
Usage:
```julia
@my_threads [schedule] begin
    prebody
end for ...
    ...
end begin
    postbody
end
```
"""
macro my_threads(args...)
    na = length(args)
    if na == 4
        sched, prebody, ex, bostbody = args
        if sched isa QuoteNode
            sched = sched.value
        elseif sched isa Symbol
            # for now only allow quoted symbols
            sched = nothing
        end
        if sched !== :static
            throw(ArgumentError("unsupported schedule argument in @threads"))
        end
    elseif na == 3
        sched = :default
        prebody, ex, postbody = args
    else
        throw(ArgumentError("wrong number of arguments in @my_threads"))
    end
    if !(isa(ex, Expr) && ex.head === :for)
        throw(ArgumentError("@my_threads requires a `for` loop expression"))
    end
    if !(ex.args[1] isa Expr && ex.args[1].head === :(=))
        throw(ArgumentError("nested outer loops are not currently supported by @my_threads"))
    end
    return _my_threadsfor(ex.args[1], ex.args[2], prebody, postbody, sched)
end

@doc @my_threads

# %%
using Random

function mcpi_my_threads(N)
    a = Atomic{Int}(0)
    @my_threads begin
        rng = MersenneTwister()
        c = 0
    end for _ in 1:N
        c += rand(rng)^2 + rand(rng)^2 ≤ 1
    end begin
        atomic_add!(a, c)
    end
    4a[]/N
end

println("Julia v", VERSION)
@show Threads.nthreads()
print("mcpi_my_threads(10^8):")
@btime mcpi_my_threads(10^8)

# %%

max_stars_repo_name: genkuroki/public

max_issues_repo_name: genkuroki/public

max_forks_repo_name: genkuroki/public