text:
"""
a quantum link model of qed
"""
function nonsym_qed_qlm_ham(;m=0.8,a=1,e=1.0,l=1)
    (Xl,Yl,Zl,Il) = nonsym_spintensors(l);
    Pl = Xl + 1im*Yl; Ml = Xl - 1im*Yl;
    (Xf,Yf,Zf,If) = nonsym_spintensors(1//2);
    Pf = Xf + 1im*Yf; Mf = Xf - 1im*Yf;

    data = Array{Union{Missing,typeof(Pf)},3}(missing,2,4,4)

    data[1,1,1] = Il; data[1,end,end] = Il;
    data[2,1,1] = If; data[2,end,end] = If;

    data[1,1,end] = a*e*e*0.5*Zl*Zl
    data[2,1,end] = m*Zf

    pref = -1im/(a*2*sqrt(l*(l+1)))

    data[2,1,2] = pref*Pf
    data[1,2,2] = Pl
    data[2,2,end] = Pf

    data[2,1,3] = -pref*Mf
    data[1,3,3] = Ml
    data[2,3,end] = Mf

    return MPOHamiltonian(data)
end

"""
qed_qlm has a local gauge transform G_j, here we construct G_j^2
useful to find the groundstate of qed_qlm in the G = 0 sector
"""
function qed_qlm_G2(;l=1)
    (Xl,Yl,Zl,Il) = nonsym_spintensors(l);
    Pl = Xl + 1im*Yl; Ml = Xl - 1im*Yl;

    (Xf,Yf,Zf,If) = nonsym_spintensors(1//2);
    Pf = Xf + 1im*Yf; Mf = Xf - 1im*Yf;

    @tensor todecomp[-1 -2 -3;-4 -5 -6]:= Zl[-1;-4]*If[-2;-5]*Il[-3;-6]
    @tensor todecomp[-1 -2 -3;-4 -5 -6]+= Il[-1;-4]*If[-2;-5]*Zl[-3;-6]
    @tensor todecomp[-1 -2 -3;-4 -5 -6]-= Il[-1;-4]*(Zf+If*0.5)[-2;-5]*Il[-3;-6]
    todecomp = todecomp*todecomp

    string = MPSKit.decompose_localmpo(add_util_leg(todecomp))

    nIl = add_util_leg(Il);
    nIf = add_util_leg(If);

    data = Array{Union{Missing,typeof(string[1])},3}(missing,2,3,3)

    data[1,1,1] = nIl; data[1,end,end] = nIl;
    data[2,1,1] = nIf; data[2,end,end] = nIf;

    data[1,1,2] = string[1];
    data[2,2,2] = string[2];
    data[1,2,end] = string[3];

    return MPOHamiltonian(data)
end

"""
the lattice qed hamiltonian, using a seperate U(1) charge for even and odd matter sites, allowing one to construct manifestly guage invariant mps's
"""
function u1_qed_ham(;m=1,g=2,max_U=3)
    # 2 U_1 charges, even and odd
    link_space = Rep[U₁×U₁]((i,-i)=>1 for i in -max_U:max_U);
    odd_matter_space = Rep[U₁×U₁]((0,0)=>1,(0,1)=>1);
    even_matter_space = Rep[U₁×U₁]((-1,0)=>1,(0,0)=>1);

    #the E operator is gauge invariant, so we can make it
    E = TensorMap(ComplexF64,link_space,link_space)
    for (f1,f2) in fusiontrees(E) #there are easier ways than iterating over the fusiontrees, but I'm not sure how
        curent_charge = f1.coupled.sectors[1].charge;
        E[f1,f2] = fill(curent_charge,1,1)
    end

    #even_matter_N
    even_matter_N = TensorMap(ComplexF64,even_matter_space,even_matter_space)
    for (f1,f2) in fusiontrees(even_matter_N)
        curent_charge = f1.coupled.sectors[1].charge;
        even_matter_N[f1,f2] = fill(curent_charge == 0 ? 0 : 1,1,1);
    end

    #odd_matter_N
    odd_matter_N = TensorMap(ComplexF64,odd_matter_space,odd_matter_space)
    for (f1,f2) in fusiontrees(odd_matter_N)
        curent_charge = f1.coupled.sectors[2].charge;
        odd_matter_N[f1,f2] = fill(curent_charge == 1 ? 0 : 1,1,1);
    end

    #the U operator is not gauge invariant, so it has to also act on the virtual level!
    virtual_link_space = Rep[U₁×U₁]((0,1)=>1,(0,-1)=>1,(1,0)=>1,(-1,0)=>1);
    ou = oneunit(virtual_link_space);
    U = TensorMap(ones,ComplexF64,virtual_link_space*link_space,link_space*virtual_link_space);
    #I think - by construction - that U will act like a creation/anihilation operator

    #odd_matter_flip
    odd_matter_flip_1 = TensorMap(ones,ComplexF64,ou*odd_matter_space,odd_matter_space*virtual_link_space)
    odd_matter_flip_2 = TensorMap(ones,ComplexF64,virtual_link_space*odd_matter_space,odd_matter_space*ou)

    even_matter_flip_1 = TensorMap(ones,ComplexF64,ou*even_matter_space,even_matter_space*virtual_link_space)
    even_matter_flip_2 = TensorMap(ones,ComplexF64,virtual_link_space*even_matter_space,even_matter_space*ou)

    #think this fixes jordan wigner?
    blocks(odd_matter_flip_1)[U₁(0)⊠U₁(1)] .*= -1;
    blocks(even_matter_flip_1)[U₁(0)⊠U₁(0)] .*= -1;

    #matter - even - matter - odd
    data = Array{Any,3}(missing,4,3,3);
    data[:,1,1] .= 1.0
    data[:,3,3] .= 1.0

    #even matter site:
    data[1,1,2] = even_matter_flip_1;
    data[1,2,3] = even_matter_flip_2;
    data[1,1,3] = m * even_matter_N;

    #link
    data[2,1,3] = g^2/2 * E^2
    data[2,2,2] = U * 1im/2;

    #odd matter site
    data[3,1,2] = odd_matter_flip_1;
    data[3,2,3] = odd_matter_flip_2;
    data[3,1,3] = -m * odd_matter_N;

    #link
    data[4,1,3] = g^2/2 * E^2;
    data[4,2,2] = U* 1im/2;

    MPOHamiltonian(data)
end

max_stars_repo_name: maartenvd/MPSKitModels

max_issues_repo_name: maartenvd/MPSKitModels

max_forks_repo_name: maartenvd/MPSKitModels


text:
"""
    Client(endpoint; headers=Dict(), introspect=true)
    Client(endpoint, ws_endpoint; headers=Dict(), introspect=true)

GraphQL client. If just `endpoint` is provided, `ws_endpoint` is
assumed to be the same as endpoint with "http" replaced by "ws".

By default, introspection will be performed on the client. This can be turned
off by setting the `introspect` keyword argument to `false`.

# Fields - Public Interface
- `endpoint::String`: endpoint for queries and mutations.
- `ws_endpoint::String`: endpoint for subscriptions.
- `headers::Dict`: contains client specific headers.
- `introspection_complete::Bool`: set to `true` once introspection has been performed.
- `queries::Vector{String}`: list of available queries.
- `mutations::Vector{String}`: list of available mutations.
- `subscriptions::Vector{String}`: list of available subscriptions.

# Fields - Internal Use
- `type_to_fields_map::Dict{String, Dict{String, Dict{String, Any}}}`: maps GQL types
    to their fields.
- `query_to_type_map::Dict{String, String}`: maps GQL queries, mutations and
    subscriptions to the type(s) of their outputs(s).
- `query_to_args_map::Dict{String, Dict{String, String}}`: maps GQL queries, mutations
    and subscriptions to their arguments and types.
- `input_object_fields_to_type_map::Dict{String, Dict{String, String}}`: maps GQL input
    objects to their fields/types.
- `schema::Dict{String, Any}`: full schema of server.
- `introspected_types::Dict{String, DataType}`: dictionary containing all introspected types.

# Examples
```julia
julia> client = Client("https://countries.trevorblades.com")
GraphQLClient Client
       endpoint: https://countries.trevorblades.com
    ws_endpoint: wss://countries.trevorblades.com

julia> client = Client("https://countries.trevorblades.com", "wss://countries.trevorblades.com")
GraphQLClient Client
       endpoint: https://countries.trevorblades.com
    ws_endpoint: wss://countries.trevorblades.com
```
"""
mutable struct Client
    endpoint::String
    ws_endpoint::String
    headers::Dict
    introspection_complete::Bool
    type_to_fields_map::Dict{String, Dict{String, Dict{String, Any}}}
    query_to_type_map::Dict{String, String}
    query_to_args_map::Dict{String, Dict{String, String}}
    input_object_fields_to_type_map::Dict{String, Dict{String, String}}
    queries::Vector{String}
    mutations::Vector{String}
    subscriptions::Vector{String}
    schema::Dict{String, Any}
    introspected_types::Dict{String, DataType}
end
Client(endpoint; kwargs...) = Client(endpoint, replace(endpoint, r"^(http)" => "ws"); kwargs...)
function Client(endpoint, ws_endpoint; headers=Dict(), introspect=true)
    client = Client(
        endpoint,
        ws_endpoint,
        headers,
        false,
        Dict{String, Any}(),
        Dict{String, Any}(),
        Dict{String, Any}(),
        Dict{String, Any}(),
        String[],
        String[],
        String[],
        Dict{String, Any}(),
        Dict{String, DataType}(),
    )
    if introspect
        return full_introspection!(client)
    else
        return client
    end
end

function Base.show(io::IO, ::MIME"text/plain", client::Client)
    println(io, "GraphQLClient Client")
    println(io, "       endpoint: $(client.endpoint)")
    print(io, "    ws_endpoint: $(client.ws_endpoint)")
end

#############
# Accessors #
#############

"""
    get_queries(client::Client)

Returns all queries available from GraphQL server. If introspection
has not been performed, will run `full_introspection!(client)`.
"""
function get_queries(client::Client)
    !client.introspection_complete && full_introspection!(client)
    return client.queries
end

"""
    get_mutations(client::Client)

Returns all mutations available from GraphQL server. If introspection
has not been performed, will run `full_introspection!(client)`.
"""
function get_mutations(client::Client)
    !client.introspection_complete && full_introspection!(client)
    return client.mutations
end

"""
    get_subscriptions(client::Client)

Returns all subscriptions available from GraphQL server. If introspection
has not been performed, will run `full_introspection!(client)`.
"""
function get_subscriptions(client::Client)
    !client.introspection_complete && full_introspection!(client)
    return client.subscriptions
end

"""
    CLIENT

`Ref` which contains the global `Client`.
"""
const CLIENT = Ref{Client}()

"""
    global_graphql_client()

Retrieve the global `Client`.
"""
function global_graphql_client()
    if !isassigned(CLIENT)
        throw(GraphQLClientException("Global client not set. Use global_graphql_client(client) to set."))
    end
    return CLIENT[]
end

"""
    global_graphql_client(client::Client)

Set the global `Client.`
"""
function global_graphql_client(client::Client)
    return CLIENT[] = client
end

ext: jl


text:
# The formulas for `initialize_dual_potential!` are in Section 6.1.2 of SFVTP19.
function initialize_dual_potential!(::KL{ρ}, a::DiscreteMeasure) where {ρ}
    c = -ρ * log(length(a.log_density))
    a.dual_potential .= c
end

function initialize_dual_potential!(::TV{ρ}, a::DiscreteMeasure) where {ρ}
    c = -ρ * sign(log(length(a.log_density)))
    a.dual_potential .= c
end

# Specialized implementation for the KL-divergence
# Prop. 12 of SFVTP19.
function sinkhorn_divergence!(
    D::KL{ρ},
    C,
    a::DiscreteMeasure,
    b::DiscreteMeasure,
    ϵ = 1e-1;
    kwargs...,
) where {ρ}
    f = a.dual_potential
    g = b.dual_potential
    scaled_exp = z -> exp(-z / ρ)
    unbalanced_sinkhorn!(D, C, a, a, ϵ; kwargs...)
    term_aa = fdot(scaled_exp, a.density, f)
    unbalanced_sinkhorn!(D, C, b, b, ϵ; kwargs...)
    term_bb = fdot(scaled_exp, b.density, g)
    unbalanced_sinkhorn!(D, C, a, b, ϵ; kwargs...)
    term_ab = fdot(scaled_exp, a.density, f) + fdot(scaled_exp, b.density, g)
    return -(ρ + ϵ / 2) * (term_ab - term_aa - term_bb)
end

# Needed to avoid ambiguity errors
sinkhorn_divergence!(
    D::KL{ρ},
    C::AbstractMatrix,
    a::DiscreteMeasure,
    b::DiscreteMeasure,
    ϵ = 1e-1;
    kwargs...,
) where {ρ} = throw(ArgumentError("Must pass a cost function `C`, not a cost matrix."))


text:
using Flux, OptimizationFlux, Zygote, StochasticDiffEq
using LinearAlgebra, Statistics
println("NNPDENS_tests")
using Test, NeuralPDE

using Random
Random.seed!(100)

println("one-dimensional heat equation")
x0 = Float32[11.] # initial points
tspan = (0.0f0,5.0f0)
dt = 0.5 # time step
d = 1 # number of dimensions
m = 10 # number of trajectories (batch size)

g(X) = sum(X.^2)   # terminal condition
f(X,u,σᵀ∇u,p,t) = Float32(0.0)
μ_f(X,p,t) = zero(X) #Vector d x 1
σ_f(X,p,t) = Diagonal(ones(Float32,d)) #Matrix d x d
prob = TerminalPDEProblem(g, f, μ_f, σ_f, x0, tspan)

hls = 10 + d #hidden layer size
opt = Flux.ADAM(0.005)  #optimizer
#sub-neural network approximating solutions at the desired point
u0 = Flux.Chain(Dense(d,hls,relu),
                Dense(hls,hls,relu),
                Dense(hls,1))
# sub-neural network approximating the spatial gradients at time point
σᵀ∇u = Flux.Chain(Dense(d+1,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,d))
pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)

ans = solve(prob, pdealg, verbose=true, maxiters=200, trajectories=m,
                            alg=EM(), dt=dt, pabstol = 1f-6)

u_analytical(x,t) = sum(x.^2) .+ d*t
analytical_ans = u_analytical(x0, tspan[end])

error_l2 = sqrt((ans-analytical_ans)^2/ans^2)

println("one-dimensional heat equation")
# println("numerical = ", ans)
# println("analytical = " ,analytical_ans)
println("error_l2 = ", error_l2, "\n")
@test error_l2 < 0.1

println("high-dimensional heat equation")
d = 50 # number of dimensions
x0 = fill(8.0f0,d)
tspan = (0.0f0,2.0f0)
dt = 0.5
m = 50 # number of trajectories (batch size)

g(X) = sum(X.^2)
f(X,u,σᵀ∇u,p,t) = Float32(0.0)
μ_f(X,p,t) = zero(X) #Vector d x 1
σ_f(X,p,t) = Diagonal(ones(Float32,d)) #Matrix d x d
prob = TerminalPDEProblem(g, f, μ_f, σ_f, x0, tspan)

hls = 10 + d #hidden layer size
#sub-neural network approximating solutions at the desired point
u0 = Flux.Chain(Dense(d,hls,relu),
                Dense(hls,hls,relu),
                Dense(hls,1))
# sub-neural network approximating the spatial gradients at time point
σᵀ∇u = Flux.Chain(Dense(d+1,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,d))
pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)

ans = solve(prob, pdealg, verbose=true, maxiters=150, trajectories=m,
                            alg=EM(), dt=dt, pabstol = 1f-6)

u_analytical(x,t) = sum(x.^2) .+ d*t
analytical_ans = u_analytical(x0, tspan[end])
error_l2 = sqrt((ans - analytical_ans)^2/ans^2)

println("high-dimensional heat equation")
# println("numerical = ", ans)
# println("analytical = " ,analytical_ans)
println("error_l2 = ", error_l2, "\n")
@test error_l2 < 1.0

println("Black-Scholes-Barenblatt equation")
d = 30 # number of dimensions
x0 = repeat([1.0f0, 0.5f0], div(d,2))
tspan = (0.0f0,1.0f0)
dt = 0.2
m = 30 # number of trajectories (batch size)

r = 0.05f0
sigma = 0.4f0
f(X,u,σᵀ∇u,p,t) = r * (u - sum(X.*σᵀ∇u))
g(X) = sum(X.^2)
μ_f(X,p,t) = zero(X) #Vector d x 1
σ_f(X,p,t) = Diagonal(sigma*X) #Matrix d x d
prob = TerminalPDEProblem(g, f, μ_f, σ_f, x0, tspan)

hls  = 10 + d #hide layer size
opt = Flux.ADAM(0.001)
u0 = Flux.Chain(Dense(d,hls,relu),
                Dense(hls,hls,relu),
                Dense(hls,1))
σᵀ∇u = Flux.Chain(Dense(d+1,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,d))
pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)

ans = solve(prob, pdealg, verbose=true, maxiters=150, trajectories=m,
                            alg=EM(), dt=dt, pabstol = 1f-6)

u_analytical(x, t) = exp((r + sigma^2).*(tspan[end] .- tspan[1])).*sum(x.^2)
analytical_ans = u_analytical(x0, tspan[1])
error_l2 = sqrt((ans .- analytical_ans)^2/ans^2)

println("Black Scholes Barenblatt equation")
# println("numerical ans= ", ans)
# println("analytical ans = " , analytical_ans)
println("error_l2 = ", error_l2, "\n")
@test error_l2 < 1.0

# Allen-Cahn Equation
d = 10 # number of dimensions
x0 = fill(0.0f0,d)
tspan = (0.3f0,0.6f0)
dt = 0.015 # time step
m = 20 # number of trajectories (batch size)

g(X) = 1.0 / (2.0 + 0.4*sum(X.^2))
f(X,u,σᵀ∇u,p,t) = u .- u.^3
μ_f(X,p,t) = zero(X) #Vector d x 1
σ_f(X,p,t) = Diagonal(ones(Float32,d)) #Matrix d x d
prob = TerminalPDEProblem(g, f, μ_f, σ_f, x0, tspan)

hls = 20 + d #hidden layer size
opt = Flux.ADAM(5^-3)  #optimizer
#sub-neural network approximating solutions at the desired point
u0 = Flux.Chain(Dense(d,hls,relu),
                Dense(hls,hls,relu),
                Dense(hls,1))

# sub-neural network approximating the spatial gradients at time point
σᵀ∇u = Flux.Chain(Dense(d+1,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,d))
pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)

ans = solve(prob, pdealg, verbose=true, maxiters=150, trajectories=m,
                            alg=EM(), dt=dt, pabstol = 1f-6)

prob_ans = 0.30879
error_l2 = sqrt((ans - prob_ans)^2/ans^2)

println("Allen-Cahn equation")
# println("numerical = ", ans)
# println("prob_ans = " , prob_ans)
println("error_l2 = ", error_l2, "\n")
@test error_l2 < 1.0

# Hamilton Jacobi Bellman Equation
d = 30 # number of dimensions
x0 = fill(0.0f0,d)
tspan = (0.0f0, 1.0f0)
dt = 0.2
m = 30 # number of trajectories (batch size)
λ = 1.0f0
#
g(X) = log(0.5f0 + 0.5f0*sum(X.^2))
f(X,u,σᵀ∇u,p,t) = -λ*sum(σᵀ∇u.^2)
μ_f(X,p,t) = zero(X)  #Vector d x 1 λ
σ_f(X,p,t) = Diagonal(sqrt(2.0f0)*ones(Float32,d)) #Matrix d x d
prob = TerminalPDEProblem(g, f, μ_f, σ_f, x0, tspan)

hls = 256 #hidden layer size
opt = Flux.ADAM(0.1)  #optimizer
#sub-neural network approximating solutions at the desired point
u0 = Flux.Chain(Dense(d,hls,relu),
                Dense(hls,hls,relu),
                Dense(hls,hls,relu),
                Dense(hls,1))
# sub-neural network approximating the spatial gradients at time point
σᵀ∇u = Flux.Chain(Dense(d+1,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,d))
pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)

@time ans = solve(prob, pdealg, verbose=true, maxiters=150, trajectories=m,
                            alg=EM(), dt=dt, pabstol = 1f-4)

T = tspan[2]
MC = 10^5
W() = randn(d,1)
u_analytical(x, t) = -(1/λ)*log(mean(exp(-λ*g(x .+ sqrt(2.0)*abs.(T-t).*W())) for _ = 1:MC))
analytical_ans = u_analytical(x0, tspan[1])

error_l2 = sqrt((ans - analytical_ans)^2/ans^2)

println("Hamilton Jacobi Bellman Equation")
# println("numerical = ", ans)
# println("analytical = " , analytical_ans)
println("error_l2 = ", error_l2, "\n")
@test error_l2 < 1.0

# Nonlinear Black-Scholes Equation with Default Risk
d = 20 # number of dimensions
x0 = fill(100.0f0,d)
tspan = (0.0f0,1.0f0)
dt = 0.125 # time step
m = 20 # number of trajectories (batch size)

g(X) = minimum(X)
δ = 2.0f0/3
R = 0.02f0
f(X,u,σᵀ∇u,p,t) = -(1 - δ)*Q(u)*u - R*u

vh = 50.0f0
vl = 70.0f0
γh = 0.2f0
γl = 0.02f0
function Q(u)
    Q = 0
    if u < vh
        Q = γh
    elseif  u >= vl
        Q = γl
    else  #if  u >= vh && u < vl
        Q = ((γh - γl) / (vh - vl)) * (u - vh) + γh
    end
end

µc = 0.02f0
σc = 0.2f0

μ_f(X,p,t) = µc*X #Vector d x 1
σ_f(X,p,t) = σc*Diagonal(X) #Matrix d x d
prob = TerminalPDEProblem(g, f, μ_f, σ_f, x0, tspan)

hls = 256 #hidden layer size
opt = Flux.ADAM(0.008)  #optimizer
#sub-neural network approximating solutions at the desired point
u0 = Flux.Chain(Dense(d,hls,relu),
                Dense(hls,hls,relu),
                Dense(hls,1))

# sub-neural network approximating the spatial gradients at time point
σᵀ∇u = Flux.Chain(Dense(d+1,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,hls,relu),
                  Dense(hls,d))
pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)

@time ans = solve(prob, pdealg, verbose=true, maxiters=100, trajectories=m,
                            alg=EM(), dt=dt, pabstol = 1f-6)

prob_ans = 57.3
error_l2 = sqrt((ans - prob_ans)^2/ans^2)

println("Nonlinear Black-Scholes Equation with Default Risk")
# println("numerical = ", ans)
# println("prob_ans = " , prob_ans)
println("error_l2 = ", error_l2, "\n")
@test error_l2 < 1.0


text:
#=
Inputs:
	bs = 50 : no. of strains in each bin
	ml = false : whether to merge last bin with previous one
=#

include("linreg.jl")
include("bin.jl")

lin(x;m=0,c=0) = m*x + c # linear function

function generate_scaling_plot(repsizes, mspacerlens; bs=50, ml=false, alldata=true,
	x="N",
	ytext="mean spacer length",
	cas=true,
	pvalue=0,
	fs=12,
	saveplot=false)

	# plot all data or not
	if alldata == true
		figsize = (6.4,7); ylim1 = 25; ylim2 = 56.5; ytickvals=25:5:55
	else
	    figsize = (6.4,4.8); ylim1 = 31; ylim2 = 37.5; ytickvals = 31:37
    end

    # x-axis is N or L
    if x == "N"
    	xrange = collect(0:0.1:3); xlim1 = 0; xlim2 = 3; textx = 2.7
    	xtick_vals = log10.(vcat(collect(1:10),collect(20:10:100),collect(200:100:1000)))
		xtick_labels = string.(vcat(["1"],fill("",8),["10"],fill("",8),["100"],fill("",8),["1000"]))
		xtext="spacer repertoire size"
    elseif x == "L"
    	xrange = collect(log10(7e5):0.01:log10(1.7e7)); xlim1 = log10(7e5); xlim2 = log10(1.7e7); textx=7.09
    	xtick_vals = log10.(vcat(collect(7e5:1e5:1e6),collect(2e6:1e6:1e7)))
    	xtick_labels = vcat(fill("",3),["\$10^6\$"],fill("",8),["\$10^7\$"])
    	xtext = "genome length"
    end

	### linear regression ###

	m, c, dm, dc = linreg(x=log10.(repsizes), y=mspacerlens)

	m1 = round(m, digits=2)
	dm1 = round(dm, digits=2)

	propconst = log(10)/m
	dpropconst = log(10)*dm/m^2

	propconst1 = round(propconst, digits=1)
	dpropconst1 = round(dpropconst, digits=1);

	if cas == false
		slopetext = "fit, slope = $m1 \$\\pm\$ $dm1"
	else
		slopetext = "fit, \$\\ln\\,$x = ($propconst1 \\pm $dpropconst1)\\ l_s\$ + const."
	end

	### bin strains by repsize, and find mean and standard error in bins ###
	# Sort repsizes and partition them into bins of size bs.
	# If `mergelast=true`, the last bin (which may contain < bs strains) is merged with the penultimate one.
	# Strains with equal repsize are randomly shuffled and may fall in either bin if they lie on the boundary between 2 bins.

	inds_inbins = bin_equalheight(repsizes, binsize=bs, mergelast=ml)

	mean_logN_inbins = [mean(log10.(repsizes[inds])) for inds in inds_inbins]
	stderr_logN_inbins = [std(log10.(repsizes[inds]))/sqrt(length(inds)) for inds in inds_inbins]

	mean_l_inbins = [mean(mspacerlens[inds]) for inds in inds_inbins]
	stderr_l_inbins = [std(mspacerlens[inds])/sqrt(length(inds)) for inds in inds_inbins];

	### correlation coefficient ###

	rho = round(cor(log10.(repsizes),mspacerlens), digits=2)

	### plot binned data and linear fit ###

	fig, ax = subplots(figsize=figsize, dpi=200)

	# plot mean and standard error of binned data
	errorbar(mean_logN_inbins, mean_l_inbins, linestyle="none", marker="o", markersize=1,
	    xerr=stderr_logN_inbins, yerr=stderr_l_inbins, elinewidth=1, capsize=0)

	# plot linear fit
	plot(xrange, lin.(xrange, m=m, c=c))

	# plot unbinned data
	if alldata == true
		scatter(log10.(repsizes), mspacerlens, color="green", s=0.4, alpha=0.5)

		if pvalue == 0
			pvaluetext = "<10^{-6}"
		else
			pvaluetext = string("=",round(pvalue, digits=2))
		end
		text(textx,49.8,"Pearson's \$r = $rho,\\ p$pvaluetext\$", fontsize=fs)

		legendtext = [slopetext, "individual datapoints", "data (binned)"]
    	legendloc = (0.9,0.83)
    	filename = "fit_alldata.svg"
    else
    	legendtext = [slopetext, "data (binned)"]
        legendloc = "upper left"
        filename = "fit.svg"
	end

	xlabel(xtext, fontsize=fs)
	xlim(xlim1,xlim2)
	xticks(xtick_vals, xtick_labels, fontsize=fs)

	ylabel(ytext, fontsize=fs)
	ylim(ylim1,ylim2)
	yticks(ytickvals, fontsize=fs)

	legend(legendtext, fontsize=fs, loc=legendloc)

	ax.spines["top"].set_visible(false)
	ax.spines["right"].set_visible(false)

	if saveplot==true
		savefig(filename,format="svg")
	end

	return rho, m, c, dm, dc, propconst1, dpropconst1, inds_inbins, mean_logN_inbins, mean_l_inbins, stderr_logN_inbins, stderr_l_inbins
end


text:
const bulma_grid_section_template = mt"""
{{{description}}}

```@raw html
<div class="columns is-multiline bulma-grid-card-section">
```

{{{cards}}}

```@raw html
</div>
```
"""

const bulma_grid_card_template = mt"""
```@raw html
<div class="column is-half">
    <div class="card bulma-grid-card">
        <div class="card-content">
            <h3 class="is-size-5 bulma-grid-card-text">
                {{{title}}}
            </h3>
            <p class="is-size-6 bulma-grid-card-description">
                {{{description}}}
            </p>
        </div>
        <div class="card-image bulma-grid-card-cover">
```
[![card image]({{{coverpath}}})](@ref {{id}})
```@raw html
        </div>
    </div>
</div>
```
"""

function cardtheme(::Val{:bulmagrid})
    templates = Dict(
        "card" => bulma_grid_card_template,
        "section" => bulma_grid_section_template
    )
    return templates, abspath(@__DIR__, "style.css")
end


text:
"""
Approximate radius of the Earth (km) used in geoemetry functions.
"""
const RADIUS_OF_EARTH_KM = 6371.0

"""
Factor used to convert speed from mph to kph units.
"""
const KMH_PER_MPH = 1.60934

"""
URLs to OpenStreetMap APIs.
"""
const OSM_URLS = Dict(
    :nominatim_search => "https://nominatim.openstreetmap.org/search",
    :overpass_map => "http://overpass-api.de/api/interpreter",
    :overpass_status => "https://overpass-api.de/api/status"
)

"""
Exclusion filters used when querying different OpenStreetMap networks.
"""
const WAY_EXCLUSION_FILTERS = Dict(
     :drive_mainroads => Dict(
        "area" => ["yes"],
        "highway" => ["cycleway", "footway", "path", "pedestrian", "steps", "track", "corridor", "elevator", "escalator", "proposed", "construction", "bridleway", "abandoned", "platform", "raceway", "service", "residential"],
        "motor_vehicle" => ["no"],
        "motorcar" => ["no"],
        "access" => ["private"],
        "service" => ["parking", "parking_aisle", "driveway", "private", "emergency_access"]
    ),
    :drive => Dict(
        "area" => ["yes"],
        "highway" => ["cycleway", "footway", "path", "pedestrian", "steps", "track", "corridor", "elevator", "escalator", "proposed", "construction", "bridleway", "abandoned", "platform", "raceway", "service"],
        "motor_vehicle" => ["no"],
        "motorcar" => ["no"],
        "access" => ["private"],
        "service" => ["parking", "parking_aisle", "driveway", "private", "emergency_access"]
    ),
    :drive_service => Dict(
        "area" => ["yes"],
        "highway" => ["cycleway", "footway", "path", "pedestrian", "steps", "track", "corridor", "elevator", "escalator", "proposed", "construction", "bridleway", "abandoned", "platform", "raceway"],
        "motor_vehicle" => ["no"],
        "motorcar" => ["no"],
        "access" => ["private"],
        "service" => ["parking", "parking_aisle", "private", "emergency_access"]
    ),
    :walk => Dict(
        "area" => ["yes"],
        "highway" => ["cycleway", "motor", "proposed", "construction", "abandoned", "platform", "raceway"],
        "foot" => ["no"],
        "access" => ["private"],
        "service" => ["private"]
    ),
    :bike => Dict(
        "area" => ["yes"],
        "highway" => ["footway", "steps", "corridor", "elevator", "escalator", "motor", "proposed", "construction", "abandoned", "platform", "raceway"],
        "bicycle" => ["no"],
        "access" => ["private"],
        "service" => ["private"]
    ),
    :all => Dict(
        "area" => ["yes"],
        "highway" => ["proposed", "construction", "abandoned", "platform", "raceway"],
        "access" => ["private"],
        "service" => ["private"]
    ),
    :all_private => Dict(
        "area" => ["yes"],
        "highway" => ["proposed", "construction", "abandoned", "platform", "raceway"]
    ),
    :none => Dict{AbstractString,Vector{AbstractString}}(),
    :rail => Dict("highway" => ["proposed", "platform"])
)

"""
Concantenates OpenStreetMap exclusion filters into a query string.
"""
function concatenate_exclusions(exclusions::AbstractDict{S,Vector{S}})::S where {S <: AbstractString}
    filters = ""

    for (k, v) in exclusions
        filters *= """["$k"!~"$(join(v, '|'))"]"""
    end

    return filters
end

"""
OpenStreetMap query strings used for different transport networks, to test queries see `https://overpass-api.de/query_form.html`.
"""
const WAY_FILTERS_QUERY = Dict(
    :drive_mainroads => """["highway"]$(concatenate_exclusions(WAY_EXCLUSION_FILTERS[:drive_mainroads]))""",
    :drive => """["highway"]$(concatenate_exclusions(WAY_EXCLUSION_FILTERS[:drive]))""",
    :drive_service => """["highway"]$(concatenate_exclusions(WAY_EXCLUSION_FILTERS[:drive_service]))""",
    :walk => """["highway"]$(concatenate_exclusions(WAY_EXCLUSION_FILTERS[:walk]))""",
    :bike => """["highway"]$(concatenate_exclusions(WAY_EXCLUSION_FILTERS[:bike]))""",
    :all => """["highway"]$(concatenate_exclusions(WAY_EXCLUSION_FILTERS[:all]))""",
    :all_private => """["highway"]$(concatenate_exclusions(WAY_EXCLUSION_FILTERS[:all_private]))""",
    :none => """["highway"]$(concatenate_exclusions(WAY_EXCLUSION_FILTERS[:none]))""",
    :rail => """["railway"]$(concatenate_exclusions(WAY_EXCLUSION_FILTERS[:rail]))"""
)

"""
OpenStreetMap query strings used for getting relation data in addition to nodes and ways that are contained within.
"""
const RELATION_FILTERS_QUERY = Dict(
    :drive_mainroads => """["type"="restriction"]["restriction"][!"conditional"][!"hgv"]""",
    :drive => """["type"="restriction"]["restriction"][!"conditional"][!"hgv"]""",
    :drive_service => """["type"="restriction"]["restriction"][!"conditional"][!"hgv"]""",
    :walk => nothing,
    :bike => nothing,
    :all => """["type"="restriction"]["restriction"][!"conditional"][!"hgv"]""",
    :all_private => """["type"="restriction"]["restriction"][!"conditional"][!"hgv"]""",
    :none => nothing,
    :rail => nothing
)

"""
OpenStreetMap metadata tags.
"""
const OSM_METADATA = Set(["version", "timestamp", "changeset", "uid", "user", "generator", "note", "meta", "count"])

"""
OpenStreetMap download formats.
"""
const OSM_DOWNLOAD_FORMAT = Dict(
    :osm => "xml",
    :xml => "xml",
    :json => "json"
)

"""
Default maxspeed based on highway type.
"""
const DEFAULT_MAXSPEEDS = Dict(
    "motorway" => 100,
    "trunk" => 100,
    "primary" => 100,
    "secondary" => 100,
    "tertiary" => 50,
    "unclassified" => 50,
    "residential" => 50,
    "other" => 50
)

"""
Default number of lanes based on highway type.
"""
const DEFAULT_LANES = Dict(
    "motorway" => 3,
    "trunk" => 3,
    "primary" => 2,
    "secondary" => 2,
    "tertiary" => 1,
    "unclassified" => 1,
    "residential" => 1,
    "other" => 1
)

"""
Default oneway attribute based on highway type.
"""
const DEFAULT_ONEWAY = Dict(
    "motorway" => false,
    "trunk" => false,
    "primary" => false,
    "secondary" => false,
    "tertiary" => false,
    "unclassified" => false,
    "residential" => false,
    "roundabout" => true,
    "other" => false
)

"""
Values that a determine a road is oneway.
"""
const ONEWAY_TRUE = Set(["true", "yes", "1", "-1", 1, -1])

"""
Values that a determine a road is not oneway.
"""
const ONEWAY_FALSE = Set(["false", "no", "0", 0])

"""
Default factor applied to maxspeed when the `lane_efficiency` weight is used to contruct OSMGraph object.
"""
const LANE_EFFICIENCY = Dict(
    1 => 0.7,
    2 => 0.8,
    3 => 0.9,
    4 => 1.0
)

"""
Default data types used to construct OSMGraph object.
"""
const DEFAULT_OSM_ID_TYPE = Int64
const DEFAULT_OSM_INDEX_TYPE = Int32
const DEFAULT_OSM_EDGE_WEIGHT_TYPE = Float64
const DEFAULT_OSM_MAXSPEED_TYPE = Int16
const DEFAULT_OSM_LANES_TYPE = Int8

"""
Default height of buildings in metres.
"""
const DEFAULT_BUILDING_HEIGHT_PER_LEVEL = 4

"""
Default maximum levels of buildings.
"""
const DEFAULT_MAX_BUILDING_LEVELS = 3

"""
Delimiters used to clean maxspeed and lanes data.
"""
const COMMON_OSM_STRING_DELIMITERS = r"[+^:;,|-]"


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/FEMBeam.jl/blob/master/LICENSE

using FEMBeam
using FEMBase.Test

# ABAQUS code to test:
# *NODE
#  1, 0.0, 0.0, 0.0
#  2, 0.0, 1.0, 0.0
# *ELEMENT, TYPE=B33, ELSET=EALL
#  1, 1,2
# *BEAM GENERAL SECTION, ELSET=EALL, DENSITY=100.0, SECTION=GENERAL
#  2.5E-2, 1.0E-4, 0.0E-4, 1.5E-4, 3.0E-5
#  1.0, 0.0, 0.0
#  100.0E2, 80.0E2
# *STEP
# *MATRIX GENERATE, STIFFNESS, MASS, LOAD
# *MATRIX OUTPUT, STIFFNESS, MASS, LOAD
# *DLOAD
#  1, P1, 0.0E3
#  1, P2, 0.0E3
# ** 1, PX, 54.0E3
#  1, PY, 72.0E3
# ** 1, PZ, 108.0E3
# *END STEP
# *STEP
# *STATIC
# *BOUNDARY
#  1, 1,6
# *DLOAD
#  1, P1, 0.0E3
#  1, P2, 0.0E3
# ** 1, PX, 54.0E3
#  1, PY, 72.0E3
# ** 1, PZ, 108.0E3
# *NODE PRINT
#  U, COORD
# *END STEP

L = 1.0
X = Dict(
    1 => [0.0, 0.0, 0.0],
    2 => [L, 0.0, 0.0])
beam = Element(Seg2, [1, 2])
E = 210.0e6
A = 20.0e-2
G = 84.0e6
Ix = J = 5.0e-5
I2 = Iy = 10.0e-5
I1 = Iz = 20.0e-5
update!(beam, "geometry", X)
update!(beam, "youngs modulus", E)
update!(beam, "shear modulus", G)
update!(beam, "cross-section area", A)
update!(beam, "torsional moment of inertia 1", I1)
update!(beam, "torsional moment of inertia 2", I2)
update!(beam, "polar moment of inertia", J)
update!(beam, "orientation", [1.0 0.0 0.0; 0.0 1.0 0.0; 0.0 0.0 1.0])
problem = Problem(Beam, "beam", 6)
add_elements!(problem, [beam])
assemble!(problem, 0.0)
K = Matrix(sparse(problem.assembly.K))
K_expected = [
    E*A/L 0 0 0 0 0 -E*A/L 0 0 0 0 0
    0 12*E*Iz/L^3 0 0 0 6*E*Iz/L^2 0 -12*E*Iz/L^3 0 0 0 6*E*Iz/L^2
    0 0 12*E*Iy/L^3 0 -6*E*Iy/L^2 0 0 0 -12*E*Iy/L^3 0 -6*E*Iy/L^2 0
    0 0 0 G*Ix/L 0 0 0 0 0 -G*Ix/L 0 0
    0 0 -6*E*Iy/L^2 0 4*E*Iy/L 0 0 0 6*E*Iy/L^2 0 2*E*Iy/L 0
    0 6*E*Iz/L^2 0 0 0 4*E*Iz/L 0 -6*E*Iz/L^2 0 0 0 2*E*Iz/L
    -E*A/L 0 0 0 0 0 E*A/L 0 0 0 0 0
    0 -12*E*Iz/L^3 0 0 0 -6*E*Iz/L^2 0 12*E*Iz/L^3 0 0 0 -6*E*Iz/L^2
    0 0 -12*E*Iy/L^3 0 6*E*Iy/L^2 0 0 0 12*E*Iy/L^3 0 6*E*Iy/L^2 0
    0 0 0 -G*Ix/L 0 0 0 0 0 G*J/L 0 0
    0 0 -6*E*Iy/L^2 0 2*E*Iy/L 0 0 0 6*E*Iy/L^2 0 4*E*Iy/L 0
    0 6*E*Iz/L^2 0 0 0 2*E*Iz/L 0 -6*E*Iz/L^2 0 0 0 4*E*Iz/L
    ]

if !isapprox(K, K_expected)
    K_diff = K-K_expected
    K_diff[abs.(K_diff) .< 1.0e-9] = 0
    println("K")
    display(K)
    println("K_expected")
    display(K_expected)
    println("difference")
    display(K_diff)
end

function solution(X2, n, q1, q2, qx, qy, qz; T=nothing)
    X = Dict(1 => [0.0, 0.0, 0.0], 2 => X2)
    beam = Element(Seg2, [1, 2])
    update!(beam, "geometry", X)
    update!(beam, "normal", n)
    update!(beam, "youngs modulus", 100.0e2)
    update!(beam, "shear modulus", 80.0e2)
    update!(beam, "cross-section area", 3.0e-2)
    update!(beam, "torsional moment of inertia 1", 1.0e-4)
    update!(beam, "torsional moment of inertia 2", 1.5e-4)
    update!(beam, "polar moment of inertia", 3.0e-5)
    update!(beam, "distributed load 1", q1)
    update!(beam, "distributed load 2", q2)
    update!(beam, "distributed load x", qx)
    update!(beam, "distributed load y", qy)
    update!(beam, "distributed load z", qz)
    if T != nothing
        update!(beam, "orientation", T)
    end
    time = 0.0
    X1, X2 = beam("geometry", time)
    L = norm(X2-X1)
    n1 = beam("normal", time)
    t = (X2-X1)/L
    n2 = cross(t,n1)
    #println("t = $t, n1  = $n1, n2 = $n2")
    problem = Problem(Beam, "example 1", 6)
    add_elements!(problem, [beam])
    assemble!(problem, time)
    K = Matrix(sparse(problem.assembly.K))
    f = Vector(sparse(problem.assembly.f)[:])
    u = zeros(12)
    fd = 7:12
    u[fd] = K[fd,fd]\f[fd]
    return K, u, f
end

include("read_abaqus.jl")

KA1, uA1, fA1 = solution([1.0, 0.0, 0.0], [0.0, 1.0, 0.0], 0.0, 0.0, 0.0, 72.0e3, 108.0e3)
KA1_expected = read_mtx_from_str("""
-1,1, -1,1,  3.333333333333318e+01
-1,2, -1,1, -8.333333333333497e+00
1,1, -1,1,  2.500000000000001e+01
2,1, -1,1, -2.500000000000001e+01
-1,2, -1,2,  3.333333333333317e+01
1,1, -1,2,  2.500000000000002e+01
2,1, -1,2, -2.500000000000002e+01
1,1, 1,1,  3.000000000000020e+02
2,1, 1,1, -3.000000000000020e+02
1,2, 1,2,  1.799999999999999e+01
1,6, 1,2,  9.000000000000000e+00
2,2, 1,2, -1.799999999999999e+01
2,6, 1,2,  8.999999999999996e+00
1,3, 1,3,  1.200000000000000e+01
1,5, 1,3, -5.999999999999999e+00
2,3, 1,3, -1.200000000000000e+01
2,5, 1,3, -5.999999999999998e+00
1,4, 1,4,  2.400000000000012e-01
2,4, 1,4, -2.400000000000012e-01
1,5, 1,5,  4.000000000000005e+00
2,3, 1,5,  5.999999999999999e+00
2,5, 1,5,  1.999999999999994e+00
1,6, 1,6,  6.000000000000007e+00
2,2, 1,6, -8.999999999999998e+00
2,6, 1,6,  2.999999999999991e+00
2,1, 2,1,  3.000000000000020e+02
2,2, 2,2,  1.799999999999999e+01
2,6, 2,2, -8.999999999999996e+00
2,3, 2,3,  1.200000000000000e+01
2,5, 2,3,  5.999999999999998e+00
2,4, 2,4,  2.400000000000012e-01
2,5, 2,5,  4.000000000000004e+00
2,6, 2,6,  6.000000000000007e+00""")
@test isapprox(KA1, KA1_expected)

KA2, uA2, fA2 = solution([1.0, 0.0, 0.0], [0.0, -1.0, 0.0], 0.0, 0.0, 0.0, 72.0e3, 108.0e3)
KA2_expected = read_mtx_from_str("""
-1,1, -1,1,  3.333333333333318e+01
-1,2, -1,1, -8.333333333333497e+00
1,1, -1,1,  2.500000000000001e+01
2,1, -1,1, -2.500000000000001e+01
-1,2, -1,2,  3.333333333333317e+01
1,1, -1,2,  2.500000000000002e+01
2,1, -1,2, -2.500000000000002e+01
1,1, 1,1,  3.000000000000020e+02
2,1, 1,1, -3.000000000000020e+02
1,2, 1,2,  1.799999999999999e+01
1,6, 1,2,  9.000000000000000e+00
2,2, 1,2, -1.799999999999999e+01
2,6, 1,2,  8.999999999999996e+00
1,3, 1,3,  1.200000000000000e+01
1,5, 1,3, -5.999999999999999e+00
2,3, 1,3, -1.200000000000000e+01
2,5, 1,3, -5.999999999999998e+00
1,4, 1,4,  2.400000000000012e-01
2,4, 1,4, -2.400000000000012e-01
1,5, 1,5,  4.000000000000005e+00
2,3, 1,5,  5.999999999999999e+00
2,5, 1,5,  1.999999999999994e+00
1,6, 1,6,  6.000000000000007e+00
2,2, 1,6, -8.999999999999998e+00
2,6, 1,6,  2.999999999999991e+00
2,1, 2,1,  3.000000000000020e+02
2,2, 2,2,  1.799999999999999e+01
2,6, 2,2, -8.999999999999996e+00
2,3, 2,3,  1.200000000000000e+01
2,5, 2,3,  5.999999999999998e+00
2,4, 2,4,  2.400000000000012e-01
2,5, 2,5,  4.000000000000004e+00
2,6, 2,6,  6.000000000000007e+00
""")
@test isapprox(KA2, KA2_expected)

KA3, uA3, fA3 = solution([1.0, 0.0, 0.0], [0.0, 0.0, 1.0], 0.0, 0.0, 0.0, 72.0e3, 108.0e3)
KA3_expected = read_mtx_from_str("""
-1,1, -1,1,  3.333333333333318e+01
-1,2, -1,1, -8.333333333333497e+00
1,1, -1,1,  2.500000000000001e+01
2,1, -1,1, -2.500000000000001e+01
-1,2, -1,2,  3.333333333333317e+01
1,1, -1,2,  2.500000000000002e+01
2,1, -1,2, -2.500000000000002e+01
1,1, 1,1,  3.000000000000020e+02
2,1, 1,1, -3.000000000000020e+02
1,2, 1,2,  1.200000000000000e+01
1,6, 1,2,  5.999999999999999e+00
2,2, 1,2, -1.200000000000000e+01
2,6, 1,2,  5.999999999999998e+00
1,3, 1,3,  1.799999999999999e+01
1,5, 1,3, -9.000000000000000e+00
2,3, 1,3, -1.799999999999999e+01
2,5, 1,3, -8.999999999999996e+00
1,4, 1,4,  2.400000000000012e-01
2,4, 1,4, -2.400000000000012e-01
1,5, 1,5,  6.000000000000007e+00
2,3, 1,5,  8.999999999999998e+00
2,5, 1,5,  2.999999999999991e+00
1,6, 1,6,  4.000000000000005e+00
2,2, 1,6, -5.999999999999999e+00
2,6, 1,6,  1.999999999999994e+00
2,1, 2,1,  3.000000000000020e+02
2,2, 2,2,  1.200000000000000e+01
2,6, 2,2, -5.999999999999998e+00
2,3, 2,3,  1.799999999999999e+01
2,5, 2,3,  8.999999999999996e+00
2,4, 2,4,  2.400000000000012e-01
2,5, 2,5,  6.000000000000007e+00
2,6, 2,6,  4.000000000000004e+00
""")
@test isapprox(KA3, KA3_expected)

KA4, uA4, fA4 = solution([1.0, 0.0, 0.0], [0.0, 0.0, -1.0], 0.0, 0.0, 0.0, 72.0e3, 108.0e3)
KA4_expected = read_mtx_from_str("""
-1,1, -1,1,  3.333333333333318e+01
-1,2, -1,1, -8.333333333333497e+00
1,1, -1,1,  2.500000000000001e+01
2,1, -1,1, -2.500000000000001e+01
-1,2, -1,2,  3.333333333333317e+01
1,1, -1,2,  2.500000000000002e+01
2,1, -1,2, -2.500000000000002e+01
1,1, 1,1,  3.000000000000020e+02
2,1, 1,1, -3.000000000000020e+02
1,2, 1,2,  1.200000000000000e+01
1,6, 1,2,  5.999999999999999e+00
2,2, 1,2, -1.200000000000000e+01
2,6, 1,2,  5.999999999999998e+00
1,3, 1,3,  1.799999999999999e+01
1,5, 1,3, -9.000000000000000e+00
2,3, 1,3, -1.799999999999999e+01
2,5, 1,3, -8.999999999999996e+00
1,4, 1,4,  2.400000000000012e-01
2,4, 1,4, -2.400000000000012e-01
1,5, 1,5,  6.000000000000007e+00
2,3, 1,5,  8.999999999999998e+00
2,5, 1,5,  2.999999999999991e+00
1,6, 1,6,  4.000000000000005e+00
2,2, 1,6, -5.999999999999999e+00
2,6, 1,6,  1.999999999999994e+00
2,1, 2,1,  3.000000000000020e+02
2,2, 2,2,  1.200000000000000e+01
2,6, 2,2, -5.999999999999998e+00
2,3, 2,3,  1.799999999999999e+01
2,5, 2,3,  8.999999999999996e+00
2,4, 2,4,  2.400000000000012e-01
2,5, 2,5,  6.000000000000007e+00
2,6, 2,6,  4.000000000000004e+00
""")
@test isapprox(KA4, KA4_expected)

KB1, uB1, fB1 = solution([0.0, 1.0, 0.0], [1.0, 0.0, 0.0], 0.0, 0.0, 0.0, 72.0e3, 108.0e3)
KB1_expected = read_mtx_from_str("""
-1,1, -1,1,  3.333333333333318e+01
-1,2, -1,1, -8.333333333333497e+00
1,2, -1,1,  2.500000000000001e+01
2,2, -1,1, -2.500000000000001e+01
-1,2, -1,2,  3.333333333333317e+01
1,2, -1,2,  2.500000000000002e+01
2,2, -1,2, -2.500000000000002e+01
1,1, 1,1,  1.799999999999999e+01
1,6, 1,1, -9.000000000000000e+00
2,1, 1,1, -1.799999999999999e+01
2,6, 1,1, -8.999999999999996e+00
1,2, 1,2,  3.000000000000020e+02
2,2, 1,2, -3.000000000000020e+02
1,3, 1,3,  1.200000000000000e+01
1,4, 1,3,  5.999999999999999e+00
2,3, 1,3, -1.200000000000000e+01
2,4, 1,3,  5.999999999999998e+00
1,4, 1,4,  4.000000000000005e+00
2,3, 1,4, -5.999999999999999e+00
2,4, 1,4,  1.999999999999994e+00
1,5, 1,5,  2.400000000000012e-01
2,5, 1,5, -2.400000000000012e-01
1,6, 1,6,  6.000000000000007e+00
2,1, 1,6,  8.999999999999998e+00
2,6, 1,6,  2.999999999999991e+00
2,1, 2,1,  1.799999999999999e+01
2,6, 2,1,  8.999999999999996e+00
2,2, 2,2,  3.000000000000020e+02
2,3, 2,3,  1.200000000000000e+01
2,4, 2,3, -5.999999999999998e+00
2,4, 2,4,  4.000000000000004e+00
2,5, 2,5,  2.400000000000012e-01
2,6, 2,6,  6.000000000000007e+00
""")
@test isapprox(KB1, KB1_expected)

KB2, uB2, fB2 = solution([0.0, 1.0, 0.0], [-1.0, 0.0, 0.0], 0.0, 0.0, 0.0, 72.0e3, 108.0e3)
KB2_expected = read_mtx_from_str("""
-1,1, -1,1,  3.333333333333318e+01
-1,2, -1,1, -8.333333333333497e+00
1,2, -1,1,  2.500000000000001e+01
2,2, -1,1, -2.500000000000001e+01
-1,2, -1,2,  3.333333333333317e+01
1,2, -1,2,  2.500000000000002e+01
2,2, -1,2, -2.500000000000002e+01
1,1, 1,1,  1.799999999999999e+01
1,6, 1,1, -9.000000000000000e+00
2,1, 1,1, -1.799999999999999e+01
2,6, 1,1, -8.999999999999996e+00
1,2, 1,2,  3.000000000000020e+02
2,2, 1,2, -3.000000000000020e+02
1,3, 1,3,  1.200000000000000e+01
1,4, 1,3,  5.999999999999999e+00
2,3, 1,3, -1.200000000000000e+01
2,4, 1,3,  5.999999999999998e+00
1,4, 1,4,  4.000000000000005e+00
2,3, 1,4, -5.999999999999999e+00
2,4, 1,4,  1.999999999999994e+00
1,5, 1,5,  2.400000000000012e-01
2,5, 1,5, -2.400000000000012e-01
1,6, 1,6,  6.000000000000007e+00
2,1, 1,6,  8.999999999999998e+00
2,6, 1,6,  2.999999999999991e+00
2,1, 2,1,  1.799999999999999e+01
2,6, 2,1,  8.999999999999996e+00
2,2, 2,2,  3.000000000000020e+02
2,3, 2,3,  1.200000000000000e+01
2,4, 2,3, -5.999999999999998e+00
2,4, 2,4,  4.000000000000004e+00
2,5, 2,5,  2.400000000000012e-01
2,6, 2,6,  6.000000000000007e+00
""")
@test isapprox(KB2, KB2_expected)

KB3, uB3, fB3 = solution([0.0, 1.0, 0.0], [0.0, 0.0, 1.0], 0.0, 0.0, 0.0, 72.0e3, 108.0e3)
KB3_expected = read_mtx_from_str("""
-1,1, -1,1,  3.333333333333318e+01
-1,2, -1,1, -8.333333333333497e+00
1,2, -1,1,  2.500000000000001e+01
2,2, -1,1, -2.500000000000001e+01
-1,2, -1,2,  3.333333333333317e+01
1,2, -1,2,  2.500000000000002e+01
2,2, -1,2, -2.500000000000002e+01
1,1, 1,1,  1.200000000000000e+01
1,6, 1,1, -5.999999999999999e+00
2,1, 1,1, -1.200000000000000e+01
2,6, 1,1, -5.999999999999998e+00
1,2, 1,2,  3.000000000000020e+02
2,2, 1,2, -3.000000000000020e+02
1,3, 1,3,  1.799999999999999e+01
1,4, 1,3,  9.000000000000000e+00
2,3, 1,3, -1.799999999999999e+01
2,4, 1,3,  8.999999999999996e+00
1,4, 1,4,  6.000000000000007e+00
2,3, 1,4, -8.999999999999998e+00
2,4, 1,4,  2.999999999999991e+00
1,5, 1,5,  2.400000000000012e-01
2,5, 1,5, -2.400000000000012e-01
1,6, 1,6,  4.000000000000005e+00
2,1, 1,6,  5.999999999999999e+00
2,6, 1,6,  1.999999999999994e+00
2,1, 2,1,  1.200000000000000e+01
2,6, 2,1,  5.999999999999998e+00
2,2, 2,2,  3.000000000000020e+02
2,3, 2,3,  1.799999999999999e+01
2,4, 2,3, -8.999999999999996e+00
2,4, 2,4,  6.000000000000007e+00
2,5, 2,5,  2.400000000000012e-01
2,6, 2,6,  4.000000000000004e+00
""")
@test isapprox(KB3, KB3_expected)

KB4, uB4, fB4 = solution([0.0, 1.0, 0.0], [0.0, 0.0, -1.0], 0.0, 0.0, 0.0, 72.0e3, 108.0e3)
KB4_expected = read_mtx_from_str("""
-1,1, -1,1,  3.333333333333318e+01
-1,2, -1,1, -8.333333333333497e+00
1,2, -1,1,  2.500000000000001e+01
2,2, -1,1, -2.500000000000001e+01
-1,2, -1,2,  3.333333333333317e+01
1,2, -1,2,  2.500000000000002e+01
2,2, -1,2, -2.500000000000002e+01
1,1, 1,1,  1.200000000000000e+01
1,6, 1,1, -5.999999999999999e+00
2,1, 1,1, -1.200000000000000e+01
2,6, 1,1, -5.999999999999998e+00
1,2, 1,2,  3.000000000000020e+02
2,2, 1,2, -3.000000000000020e+02
1,3, 1,3,  1.799999999999999e+01
1,4, 1,3,  9.000000000000000e+00
2,3, 1,3, -1.799999999999999e+01
2,4, 1,3,  8.999999999999996e+00
1,4, 1,4,  6.000000000000007e+00
2,3, 1,4, -8.999999999999998e+00
2,4, 1,4,  2.999999999999991e+00
1,5, 1,5,  2.400000000000012e-01
2,5, 1,5, -2.400000000000012e-01
1,6, 1,6,  4.000000000000005e+00
2,1, 1,6,  5.999999999999999e+00
2,6, 1,6,  1.999999999999994e+00
2,1, 2,1,  1.200000000000000e+01
2,6, 2,1,  5.999999999999998e+00
2,2, 2,2,  3.000000000000020e+02
2,3, 2,3,  1.799999999999999e+01
2,4, 2,3, -8.999999999999996e+00
2,4, 2,4,  6.000000000000007e+00
2,5, 2,5,  2.400000000000012e-01
2,6, 2,6,  4.000000000000004e+00
""")
@test isapprox(KB4, KB4_expected)

KC1, uC1, fC1 = solution([0.0, 0.0, 1.0], [1.0, 0.0, 0.0], 0.0, 0.0, 0.0, 72.0e3, 108.0e3)
KC1_expected = read_mtx_from_str("""
-1,1, -1,1,  3.333333333333318e+01
-1,2, -1,1, -8.333333333333497e+00
1,3, -1,1,  2.500000000000001e+01
2,3, -1,1, -2.500000000000001e+01
-1,2, -1,2,  3.333333333333317e+01
1,3, -1,2,  2.500000000000002e+01
2,3, -1,2, -2.500000000000002e+01
1,1, 1,1,  1.799999999999999e+01
1,5, 1,1,  9.000000000000000e+00
2,1, 1,1, -1.799999999999999e+01
2,5, 1,1,  8.999999999999996e+00
1,2, 1,2,  1.200000000000000e+01
1,4, 1,2, -5.999999999999999e+00
2,2, 1,2, -1.200000000000000e+01
2,4, 1,2, -5.999999999999998e+00
1,3, 1,3,  3.000000000000020e+02
2,3, 1,3, -3.000000000000020e+02
1,4, 1,4,  4.000000000000005e+00
2,2, 1,4,  5.999999999999999e+00
2,4, 1,4,  1.999999999999994e+00
1,5, 1,5,  6.000000000000007e+00
2,1, 1,5, -8.999999999999998e+00
2,5, 1,5,  2.999999999999991e+00
1,6, 1,6,  2.400000000000012e-01
2,6, 1,6, -2.400000000000012e-01
2,1, 2,1,  1.799999999999999e+01
2,5, 2,1, -8.999999999999996e+00
2,2, 2,2,  1.200000000000000e+01
2,4, 2,2,  5.999999999999998e+00
2,3, 2,3,  3.000000000000020e+02
2,4, 2,4,  4.000000000000004e+00
2,5, 2,5,  6.000000000000007e+00
2,6, 2,6,  2.400000000000012e-01
""")
@test isapprox(KC1, KC1_expected)

KC2, uC2, fC2 = solution([0.0, 0.0, 1.0], [-1.0, 0.0, 0.0], 0.0, 0.0, 0.0, 72.0e3, 108.0e3)
KC2_expected = read_mtx_from_str("""
-1,1, -1,1,  3.333333333333318e+01
-1,2, -1,1, -8.333333333333497e+00
1,3, -1,1,  2.500000000000001e+01
2,3, -1,1, -2.500000000000001e+01
-1,2, -1,2,  3.333333333333317e+01
1,3, -1,2,  2.500000000000002e+01
2,3, -1,2, -2.500000000000002e+01
1,1, 1,1,  1.799999999999999e+01
1,5, 1,1,  9.000000000000000e+00
2,1, 1,1, -1.799999999999999e+01
2,5, 1,1,  8.999999999999996e+00
1,2, 1,2,  1.200000000000000e+01
1,4, 1,2, -5.999999999999999e+00
2,2, 1,2, -1.200000000000000e+01
2,4, 1,2, -5.999999999999998e+00
1,3, 1,3,  3.000000000000020e+02
2,3, 1,3, -3.000000000000020e+02
1,4, 1,4,  4.000000000000005e+00
2,2, 1,4,  5.999999999999999e+00
2,4, 1,4,  1.999999999999994e+00
1,5, 1,5,  6.000000000000007e+00
2,1, 1,5, -8.999999999999998e+00
2,5, 1,5,  2.999999999999991e+00
1,6, 1,6,  2.400000000000012e-01
2,6, 1,6, -2.400000000000012e-01
2,1, 2,1,  1.799999999999999e+01
2,5, 2,1, -8.999999999999996e+00
2,2, 2,2,  1.200000000000000e+01
2,4, 2,2,  5.999999999999998e+00
2,3, 2,3,  3.000000000000020e+02
2,4, 2,4,  4.000000000000004e+00
2,5, 2,5,  6.000000000000007e+00
2,6, 2,6,  2.400000000000012e-01
""")
@test isapprox(KC2, KC2_expected)

KC3, uC3, fC3 = solution([0.0, 0.0, 1.0], [0.0, 1.0, 0.0], 0.0, 0.0, 0.0, 72.0e3, 108.0e3)
KC3_expected = read_mtx_from_str("""
-1,1, -1,1,  3.333333333333318e+01
-1,2, -1,1, -8.333333333333497e+00
1,3, -1,1,  2.500000000000001e+01
2,3, -1,1, -2.500000000000001e+01
-1,2, -1,2,  3.333333333333317e+01
1,3, -1,2,  2.500000000000002e+01
2,3, -1,2, -2.500000000000002e+01
1,1, 1,1,  1.200000000000000e+01
1,5, 1,1,  5.999999999999999e+00
2,1, 1,1, -1.200000000000000e+01
2,5, 1,1,  5.999999999999998e+00
1,2, 1,2,  1.799999999999999e+01
1,4, 1,2, -9.000000000000000e+00
2,2, 1,2, -1.799999999999999e+01
2,4, 1,2, -8.999999999999996e+00
1,3, 1,3,  3.000000000000020e+02
2,3, 1,3, -3.000000000000020e+02
1,4, 1,4,  6.000000000000007e+00
2,2, 1,4,  8.999999999999998e+00
2,4, 1,4,  2.999999999999991e+00
1,5, 1,5,  4.000000000000005e+00
2,1, 1,5, -5.999999999999999e+00
2,5, 1,5,  1.999999999999994e+00
1,6, 1,6,  2.400000000000012e-01
2,6, 1,6, -2.400000000000012e-01
2,1, 2,1,  1.200000000000000e+01
2,5, 2,1, -5.999999999999998e+00
2,2, 2,2,  1.799999999999999e+01
2,4, 2,2,  8.999999999999996e+00
2,3, 2,3,  3.000000000000020e+02
2,4, 2,4,  6.000000000000007e+00
2,5, 2,5,  4.000000000000004e+00
2,6, 2,6,  2.400000000000012e-01
""")
@test isapprox(KC3, KC3_expected)

KC4, uC4, fC4 = solution([0.0, 0.0, 1.0], [0.0, -1.0, 0.0], 0.0, 0.0, 0.0, 72.0e3, 108.0e3)
KC4_expected = read_mtx_from_str("""
-1,1, -1,1,  3.333333333333318e+01
-1,2, -1,1, -8.333333333333497e+00
1,3, -1,1,  2.500000000000001e+01
2,3, -1,1, -2.500000000000001e+01
-1,2, -1,2,  3.333333333333317e+01
1,3, -1,2,  2.500000000000002e+01
2,3, -1,2, -2.500000000000002e+01
1,1, 1,1,  1.200000000000000e+01
1,5, 1,1,  5.999999999999999e+00
2,1, 1,1, -1.200000000000000e+01
2,5, 1,1,  5.999999999999998e+00
1,2, 1,2,  1.799999999999999e+01
1,4, 1,2, -9.000000000000000e+00
2,2, 1,2, -1.799999999999999e+01
2,4, 1,2, -8.999999999999996e+00
1,3, 1,3,  3.000000000000020e+02
2,3, 1,3, -3.000000000000020e+02
1,4, 1,4,  6.000000000000007e+00
2,2, 1,4,  8.999999999999998e+00
2,4, 1,4,  2.999999999999991e+00
1,5, 1,5,  4.000000000000005e+00
2,1, 1,5, -5.999999999999999e+00
2,5, 1,5,  1.999999999999994e+00
1,6, 1,6,  2.400000000000012e-01
2,6, 1,6, -2.400000000000012e-01
2,1, 2,1,  1.200000000000000e+01
2,5, 2,1, -5.999999999999998e+00
2,2, 2,2,  1.799999999999999e+01
2,4, 2,2,  8.999999999999996e+00
2,3, 2,3,  3.000000000000020e+02
2,4, 2,4,  6.000000000000007e+00
2,5, 2,5,  4.000000000000004e+00
2,6, 2,6,  2.400000000000012e-01
""")
@test isapprox(KC4, KC4_expected)


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "qwt"
version = v"6.0.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/ig-or/qwt.git", "30812fc4a2630ec9956f4daa6c4f6d4b331c74a6")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/qwt
qmake
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "windows"),
    Platform("x86_64", "linux"; libc="glibc"),
]

# The products that we will ensure are always built
products = [
    LibraryProduct(["libqwt", "qwt"], :qwt),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    # TODO: in the future replace `Qt_jll` with `Qt5Base_jll` + `Qt5Svg_jll`
    Dependency(PackageSpec(name="Qt_jll", uuid="ede63266-ebff-546c-83e0-1c6fb6d0efc8"))
    Dependency(PackageSpec(name="GLEW_jll", uuid="bde7f898-03f7-559e-8810-194d950ce600"))
    Dependency(PackageSpec(name="GLFW_jll", uuid="0656b61e-2033-5cc2-a64a-77c0f6c09b89"))
    Dependency(PackageSpec(name="GLU_jll", uuid="bd17208b-e95e-5925-bf81-e2f59b3e5c61"))
]

include("../../fancy_toys.jl")

platforms_linux = filter(p -> Sys.islinux(p), platforms)
platforms_win = filter(p -> Sys.iswindows(p), platforms)

if any(should_build_platform.(triplet.(platforms_linux)))
    build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"7")
end
if any(should_build_platform.(triplet.(platforms_win)))
    build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version = v"8")
end

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
using StatsBase, Plots ; gr()

function proportionFished(gF,sF,n,N,withReplacement = false)
    function fishing()
        fishInPond = [ones(Int64,gF); zeros(Int64,sF)]
        fishCaught = Int64[]

        for fish in 1:n
            fished = rand(fishInPond)
            push!(fishCaught,fished)
            if withReplacement == false
                deleteat!(fishInPond, findfirst(x->x==fished, fishInPond))
            end
        end
        sum(fishCaught)
    end

    simulations = [fishing() for _ in 1:N]
    proportions = counts(simulations,0:n)/N

    if withReplacement
        plot!(0:n, proportions,
             line=:stem, marker=:circle, c=:blue, ms=6, msw=0,
             label="With replacement",
             xlabel="n",
             ylims=(0, 0.6), ylabel="Probability")
    else
        plot!(0:n, proportions,
            line=:stem, marker=:xcross, c=:red, ms=6, msw=0,
            label="Without replacement")
    end
end

N = 10^6
goldFish, silverFish, n = 3, 4, 3
plot()
proportionFished(goldFish, silverFish, n, N)
proportionFished(goldFish, silverFish, n, N, true)

max_stars_repo_name: newptcai/StatsWithJuliaBook

max_issues_repo_name: newptcai/StatsWithJuliaBook

max_forks_repo_name: newptcai/StatsWithJuliaBook


text:
using BotCoreLCMTypes
using LCMCore
using LCMCore: fingerprint

using Test

#=
The following Python code was used to retrieve all of the
fingerprints for the various bot_core lcmtypes automatically:

import sys
import inspect
import struct
import bot_core

for name, obj in inspect.getmembers(bot_core):
    if inspect.isclass(obj):
        print("{:s} => {:d},".format(name, struct.unpack('q', obj._get_packed_fingerprint())[0]))

import sys
import inspect
import struct
import drake

for name, obj in inspect.getmembers(drake):
    if inspect.isclass(obj):
        print("{:s} => {:d},".format(name, struct.unpack('q', obj._get_packed_fingerprint())[0]))

=#

"""
Equality test for most types, but we compare LCMType instances
by looking at their fields, since otherwise == just falls back
to === which compares object identity and is thus not very
informative.
"""
closeenough(x, y) = x == y
function closeenough(x::LCMType, y::LCMType)
    (typeof(x) == typeof(y)) || return false
    for name in fieldnames(typeof(x))
        closeenough(getfield(x, name), getfield(y, name)) || return false
    end
    true
end

@testset "BotCoreLCMTypes" begin
    expected_fingerprints_network_order = Dict(
        atlas_command_t => 1312111223708868662,
        force_torque_t => -575550166132800226,
        gps_data_t => 1591358810433180375,
        gps_satellite_info_list_t => 6630954057076318457,
        gps_satellite_info_t => -3841388522418257571,
        image_metadata_t => -8143784636512515174,
        image_sync_t => 3436534526547922253,
        image_t => -1083725853459778796,
        images_t => 4243505650783468064,
        ins_t => 4613985382350956424,
        joint_angles_t => -8005155243827168194,
        joint_state_t => -6587704215521773762,
        kvh_raw_imu_batch_t => 8667107086997737736,
        kvh_raw_imu_t => 2860203241911737633,
        planar_lidar_t => -8260152468248931869,
        pointcloud2_t => 5339961563394395403,
        pointcloud_t => 7588611463407503629,
        pointfield_t => 7874274974017539762,
        pose_t => 6778111308394862126,
        position_3d_t => 8736894522790027246,
        quaternion_t => 2236618827848571702,
        raw_t => -8160235860338714832,
        rigid_transform_t => -5853054794096467990,
        robot_state_t => 8515145021206371399,
        robot_urdf_t => 6244330983731234563,
        sensor_status_t => -5959806439264436958,
        six_axis_force_torque_array_t => -6289052062373553992,
        six_axis_force_torque_t => -1425131401932437513,
        system_status_t => -5287778015234963678,
        twist_t => 5454792568601052517,
        utime_t => 3436534526547922253,
        vector_3d_t => 2207305183986155182,
        viewer_command_t => 3337564407834538480,
        viewer_draw_t => 4918545815192686401,
        viewer_geometry_data_t => 566142053834305117,
        viewer_link_data_t => 7145691636786734417,
        viewer_load_robot_t => 4120136222412670857,
    )

    for (lcmtype, fingerprint_network_order) in expected_fingerprints_network_order
        @testset "$lcmtype" begin
            @testset "constructor" begin
                msg = lcmtype()
                for name in fieldnames(typeof(msg))
                    @test closeenough(getfield(msg, name), LCMCore.defaultval(typeof(getfield(msg, name))))
                end
            end
            @testset "fingerprint" begin
                @test hton(fingerprint(lcmtype)) == fingerprint_network_order
            end
        end
    end
end


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control
#
# This file is part of module
#   Modia3D.Composition (Modia3D/Composition/_module.jl)
#

mutable struct TreeJointFixed <: Modia3D.AbstractJoint
   _internal::ModiaMath.ComponentInternal  # Data common to assembly component classes
   obj1::Object3D                          # obj2 is fixed to obj1
   obj2::Object3D
end

function Base.show(io::IO, joint::TreeJointFixed)
   print(io,"Fixed(", ModiaMath.fullName(joint.obj1),
                   ", ", ModiaMath.fullName(joint.obj2))
end

"""
    joint = Modia3D.Fixed(obj1, obj2; r=zeros(3), R=nothing, q=nothing)

Return a Fixed `joint` that fixes `obj2::`[`Object3D`](@ref) relatively to
`obj1::`[`Object3D`](@ref) via the relative position vector from `obj1` to `obj2`
(resolved in `obj1`) and the relative transformation matrix `R` from `obj1` to `obj2`
or alternatively the relative quaternion `q`.

# Arguments
- `obj1::Object3D`: Parent object.

- `obj2::Object3D`: Object fixed relative to `obj1`.

- `r::AbstractVector`: Relative position vector from `obj1` to `obj2`
   resolved in `obj1`.

- `R::Union{ModiaMath.RotationMatrix,Nothing}`: Rotation matrix defining the rotation
   from `obj1` to `obj`. If both `R = nothing` and `q = nothing`,
   a null rotation is defined.

- `q::Union{ModiaMath.Quaternion,Nothing}`: Quaternion defining the rotation
   from `obj1` to `obj2`. If both `R = nothing` and `q = nothing`,
   a null rotation is defined.
"""
function Fixed(obj1::Object3D, obj2::Object3D;
               r::AbstractVector = ModiaMath.ZeroVector3D,
               R::Union{ModiaMath.RotationMatrix,NOTHING} = nothing,
               q::Union{ModiaMath.Quaternion,NOTHING} = nothing)::TreeJointFixed

    connect(obj1, obj2, r=r, R=R, q=q)
    if obj2.parent ≡ obj1
		obj2.joint = TreeJointFixed(ModiaMath.ComponentInternal(), obj1, obj2)
        return obj2.joint
    elseif obj1.parent ≡ obj2
		obj1.joint = TreeJointFixed(ModiaMath.ComponentInternal(), obj2, obj1)
        return obj1.joint
    else
        error("Fixed: unexpected error. obj1 = ", ModiaMath.fullName(obj1),
              ", obj2 = ", ModiaMath.fullName(obj2))
    end
end

computeKinematics!(joint::TreeJointFixed, obj::Object3D, analysis::ModiaMath.AnalysisType, time::Float64) = nothing
computeForceTorqueAndResidue!(joint::TreeJointFixed, obj::Object3D, analysis::ModiaMath.AnalysisType, time::Float64) = nothing


text:
using BenchmarkTools

using Interpolations
using Interpol

Atest = rand((10^2, 10^2, 10^2));

@benchmark @inbounds($(Atest)[LinearInterpolation, x,y,z]) setup = (x=$(1+99*rand()); y=$(1+99*rand()); z=$(1+99*rand()))
@benchmark @inbounds($(Atest)[x,y,z]) setup = (x=$(round(Int, 1+99*rand())); y=$(round(Int, 1+99*rand())); z=$(round(Int, 1+99*rand())))
@benchmark @inbounds(inplaceadd($(Atest), LinearInterpolation, 1., x,y,z)) setup = (x=$(1+99*rand()); y=$(1+99*rand()); z=$(1+99*rand()))
@benchmark $(interpolate(Atest, BSpline(Linear()), OnGrid()))[x,y,z] setup = (x=$(1+99*rand()); y=$(1+99*rand()); z=$(1+99*rand()))

@benchmark for (x,y,z) in xyz @inbounds(@inbounds($(Atest)[LinearInterpolation, x,y,z])) end setup = (xyz=$(zip((1+99*rand(1000)), (1+99*rand(1000)), (1+99*rand(1000)))))
@benchmark for (x,y,z) in xyz @inbounds($(interpolate(Atest, BSpline(Linear()), OnGrid()))[x,y,z]) end setup = (xyz=$(zip((1+99*rand(1000)), (1+99*rand(1000)), (1+99*rand(1000)))))

@benchmark for (x,y) in xy; ($Atest[ConstantBoundary, x,y]) end setup = (xy=$(zip(round(Int, -1000+2000*rand(1000)), round(Int, -1000+2000*rand(1000)))))

max_stars_repo_name: JKrehl/Interpol

max_issues_repo_name: JKrehl/Interpol

max_forks_repo_name: JKrehl/Interpol


text:
import CLBLAS
import OpenCL.cl

const clblas = CLBLAS

const COUNT = int(ARGS[1])

function seq_mat_mul_sdot{T}(Mdim::Int, Ndim::Int, Pdim::Int,
                             A::Array{T}, B::Array{T}, C::Array{T})
    for i in 1:Ndim
        for j in 1:Mdim
            tmp = zero(Float32)
            for k in 1:Pdim
                @inbounds tmp += A[(i-1)*Ndim+k] * B[(k-1)*Pdim+j]
            end
            @inbounds C[(i-1)*Ndim+j] = tmp
        end
    end
end

A = rand(Float32, 10, 10)
B = rand(Float32, 10, 10)
C = rand(Float32, 10, 10)
# Force compilation
seq_mat_mul_sdot(10, 10, 10, A, B, C)

clblas.setup(true)
device, ctx, queue = clblas.get_next_compute_context()

matSizes = Array(Tuple, COUNT)
for i in 1:COUNT
   matSizes[i] = (rand(20:100) * 32, rand(20:100) * 32, rand(20:100) * 32)
end

const alpha = 5.0
const beta  = 2.0

cpu_times        =  zeros(Float32, COUNT)
cpu_seq_times    =  zeros(Float32, COUNT)
gpu_times        =  zeros(Float32, COUNT)
gpu_times_block  =  zeros(Float32, COUNT)
gpu_times_clblas =  zeros(Float32, COUNT)

kernel_source = open(readall, "GEMM-NVIDIA.cl")
prg  = cl.Program(ctx, source=kernel_source) |> cl.build!
mmul = cl.Kernel(prg, "mmul")

kernel_source = open(readall, "block_form.cl")
prg2  = cl.Program(ctx, source=kernel_source) |> cl.build!
block_mmul = cl.Kernel(prg2, "mmul")

blocksize = 16
localmem1 = cl.LocalMem(Float32, blocksize^2)
localmem2 = cl.LocalMem(Float32, blocksize^2)

for id in 1:COUNT
   m, k, n = matSizes[id]
   A = rand(Float32, m, k)
   B = rand(Float32, k, n)
   C = rand(Float32, m, n)

   info("-----CPU------")
   startT = time()
     C_RES = alpha*A*B + beta*C
   runtime = time() - startT
   cpu_times[id] = runtime

   info("------ CPU SEQ ---------")
   startT = time()
   #seq_mat_mul_sdot(n, m, k, A, B, C)
   runtime = time() - startT
   cpu_seq_times[id] = runtime

   info("-----GPU------")
   A_buf = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=A)
   B_buf = cl.Buffer(Float32, ctx, (:r, :copy), hostbuf=B)
   C_buf = cl.Buffer(Float32, ctx, (:rw, :copy), hostbuf=C)

   evt = mmul[queue, (m, n), (blocksize, blocksize)](convert(Csize_t, m), convert(Csize_t, k),
                                                     alpha, beta, C_buf, A_buf, B_buf, localmem1, localmem2)
   runtime = evt[:profile_duration] / 1e9
   gpu_times[id] = runtime

   info("-----GPU BLOCK------")
   evt = cl.call(queue, block_mmul, (m, n), (16,16),
                 int32(m), int32(k),
                 C_buf, A_buf, B_buf, localmem1, localmem2)
   runtime = evt[:profile_duration] / 1e9
   gpu_times_block[id] = runtime

   info("-----GPU CLBLAS------")
   m = convert(Csize_t, m)
   k = convert(Csize_t, k)
   n = convert(Csize_t, n)

   lda = m
   ldb = k
   ldc = m

   event = cl.UserEvent(ctx)
   ptr_event = [pointer(event)]

   startT = time()
   clblas.clblasSgemm(1, 0, 0, m, n, k,
                           alpha, pointer(A_buf), 0, lda,
                           pointer(B_buf), 0, ldb, beta,
                           pointer(C_buf), 0, ldc,
                           1, [pointer(queue)],
                           0, C_NULL, pointer(ptr_event))

   cl.finish(queue)
   endT = time()
#   startT = cl.CL_ulong[0]
#   cl.api.clGetEventProfilingInfo( event.id, cl.CL_PROFILING_COMMAND_START, sizeof(cl.CL_long), startT, C_NULL)
#   println("Start:", startT[1])

#   endT = cl.CL_ulong[0]
#   cl.api.clGetEventProfilingInfo(event.id, cl.CL_PROFILING_COMMAND_END, sizeof(cl.CL_long), endT, C_NULL)
#   println("End:", endT[1])

#   runtime = (endT[1] - startT[1] ) / 1e9
   gpu_times_clblas[id] = endT - startT
end

using PyCall
@pyimport matplotlib.pyplot as plt
#plt.plot([1:COUNT], cpu_seq_times, color="blue", linewidth=2.0, linestyle="--", label="CPU")
plt.plot([1:COUNT], cpu_times, color="red", linewidth=2.0, linestyle="--", label="CPU")
plt.plot([1:COUNT], gpu_times, color="orange", linewidth=2.0, linestyle="--", label="GPU")
plt.plot([1:COUNT], gpu_times_block, color="green", linewidth=2.0, linestyle="--", label="GPU Block")
plt.plot([1:COUNT], gpu_times_clblas, color="black", linewidth=2.0, linestyle="--", label="GPU BLAS")

plt.ylabel("time in seconds")
plt.title("CPU vs GPU times")
plt.show()


text:
# This file was generated, do not modify it. # hide
train, test = partition(eachindex(targets), 0.70, shuffle=true,
                        rng=StableRNG(52))


text:
mutable struct LogNormalSwitches{T1} <: AbstractDistributionSwitches{T1}
    pvec :: ParamVector
end

mutable struct LogNormal{T1} <: AbstractDistributionLH{T1}
    objId :: ObjectId
    switches :: LogNormalSwitches{T1}
    lb :: T1
    mean :: T1
    std :: T1
end

Lazy.@forward LogNormal.switches (
    ModelParams.get_pvector
);

function draw(u :: LogNormal{T1}, nDims, rng :: AbstractRNG) where T1
    drawV = u.lb .+ exp.(u.mean .+ u.std .* randn(rng, T1, nDims...));
    # is_lognormal(u)  &&  exp_draws!(drawV);
    return drawV
end

function draw(u :: LogNormal{T1}, rng :: AbstractRNG)  where T1
    drawV = u.lb + exp(u.mean + u.std * rand(rng, T1));
    # is_lognormal(u)  &&  (drawV = exp(drawV));
    return drawV
end

function quantiles(u :: LogNormal{T1}, pctV)  where T1
    normalDistr = Distributions.Normal(u.mean, u.std);
    qV = u.lb .+ exp.(Distributions.quantile.(normalDistr, pctV));
    # is_lognormal(u)  &&  exp_draws!(qV);
    return qV
end

Base.eltype(::LogNormal{T1}) where T1 = T1;

isbounded(u :: LogNormal{T1}) where T1 = false;

Base.show(io :: IO, um :: LogNormal{T1}) where T1 =
    print(io, "LogNormal{$T1}($(um.mean), $(um.std))");

# Validate draws
function check_draws(u :: LogNormal{T1}, drawM) where T1
    isValid = true;
    isValid = isValid  &&  all(x -> x >= u.lb, drawM);
    return isValid
end

## ----------  Construction

ModelParams.get_pvector(switches :: LogNormalSwitches{T1}) where T1 = switches.pvec;

init_distribution(objId :: ObjectId, switches :: LogNormalSwitches{T1}) where T1 =
    init_log_normal(switches);

function init_log_normal(switches :: LogNormalSwitches{T1}) where T1
    pvec = get_pvector(switches);
    return LogNormal{T1}(switches.pvec.objId, switches,
        param_default_value(pvec, :lb),
        param_default_value(pvec, :mean),
        param_default_value(pvec, :std))
end

function init_log_normal_switches(objId :: ObjectId;
    lbVal = 1.0, meanVal = 0.0, stdVal = 1.0,
    calLb = false, calMean = true, calStd = true
    )
    T1 = typeof(meanVal);
    pLb = init_lb(; lbVal, calLb);
    pMean = init_mean(; meanVal, calMean);
    pStd = init_std(; stdVal, calStd);
    pvec = ParamVector(objId, [pLb, pMean, pStd]);
    return LogNormalSwitches{T1}(pvec)
end

function init_lb(; descr = "Lower bound", lsym = "lb", lbVal = 0.0,
    lbLb = -10.0, lbUb = 10.0, calLb = false)
    return Param(:lb, descr, lsym, lbVal, lbVal, lbLb, lbUb, calLb);
end

make_test_log_normal_switches() = init_log_normal_switches(ObjectId(:none);
    meanVal = 1.2, stdVal = 2.3);
make_test_log_normal() = init_log_normal(make_test_log_normal_switches());

# ---------------------


text:
using Dates
using Test

using MarketData

using TimeSeries

@testset "modify" begin

@testset "rename" begin
    re_ohlc = rename(ohlc, [:a, :b, :c, :d])
    re_cl   = rename(cl, [:vector])
    re_cls  = rename(cl, :symbol)

    @testset "change colnames with multi-member vector" begin
        @test colnames(re_ohlc) == [:a, :b, :c, :d]
        @test_throws ArgumentError rename(ohlc, [:a])
    end

    @testset "change colnames with single-member vector" begin
        @test colnames(re_cl) == [:vector]
        @test_throws ArgumentError rename(cl, [:a, :b])
    end

    @testset "change colnames with pair" begin
        re_ohlc_2 = rename(ohlc, :Open => :a)
        @test colnames(re_ohlc_2) == [:a, :High, :Low, :Close]
        @test_throws ArgumentError rename(ohlc, :Unknown => :A)
    end

    @testset "change colnames with several pairs" begin
        re_ohlc_2 = rename(ohlc, :Open => :a, :Close => :d)
        @test colnames(re_ohlc_2) == [:a, :High, :Low, :d]
        @test_throws MethodError rename(ohlc)
    end

    @testset "change colnames with dict" begin
        re_ohlc_2 = rename(ohlc, Dict(:Open => :a, :Close => :d)...)
        @test colnames(re_ohlc_2) == [:a, :High, :Low, :d]
    end

    @testset "change colnames with function" begin
        @testset "lambda function" begin
            f = colname -> Symbol(uppercase(string(colname)))
            re_ohlc_2 = rename(f, ohlc)
            @test colnames(re_ohlc_2) == [:OPEN, :HIGH, :LOW, :CLOSE]
        end

        @testset "function composition" begin
            f = Symbol ∘ uppercase ∘ string
            re_ohlc_2 = rename(f, ohlc)
            @test colnames(re_ohlc_2) == [:OPEN, :HIGH, :LOW, :CLOSE]
        end

        @testset "do block" begin
            re_ohlc_2 = rename(ohlc) do x
                x |> string |> uppercase |> Symbol
            end
            @test colnames(re_ohlc_2) == [:OPEN, :HIGH, :LOW, :CLOSE]
        end

        @testset "automatic string/symbol" begin
            re_ohlc_2 = rename(uppercase, ohlc, String)
            @test colnames(re_ohlc_2) == [:OPEN, :HIGH, :LOW, :CLOSE]
        end

    end
end  # @testset "rename"

@testset "rename!" begin
    let
        ta    = first(ohlc)
        cols  = [:Open, :High, :Low, :Close]
        cols′ = [:A, :B, :C, :D]
        rename!(ta, [:A, :B, :C, :D])
        @test colnames(ohlc) == cols
        @test colnames(ta)   == cols′
    end

    let
        ta = first(cl)
        rename!(ta, [:A])
        @test colnames(cl) == [:Close]
        @test colnames(ta) == [:A]
    end

    let
        ta = first(cl)
        rename!(ta, :A)
        @test colnames(cl) == [:Close]
        @test colnames(ta) == [:A]
    end

    let
        ta    = first(ohlc)
        cols  = [:Open,  :High, :Low,  :Close]
        cols′ = [:Open′, :High, :Low′, :Close]
        rename!(ta, :Open => :Open′, :Low => :Low′)
        @test colnames(ohlc) == cols
        @test colnames(ta)   == cols′
        @test_throws MethodError rename!(ta)
        @test_throws ArgumentError rename!(ohlc, :Unknown => :A)
    end

    let
        ta    = first(ohlc)
        cols  = [:Open, :High, :Low, :Close]
        cols′ = [:open, :high, :low, :close]
        rename!(Symbol ∘ lowercase ∘ string, ta)
        @test colnames(ohlc)  == cols
        @test colnames(ta)    == cols′
    end

    let
        ta    = first(ohlc)
        cols  = [:Open, :High, :Low, :Close]
        cols′ = [:open, :high, :low, :close]
        rename!(lowercase, ta, String)
        @test colnames(ohlc)  == cols
        @test colnames(ta)    == cols′
    end
end  # @testset "rename!"

end  # @testset "modify"


text:
type Vehicle
    speed::Int8
    position::Int64
    tipo::Int8
    change::Int8
    num::Int64
    len::Int64
end

type Highway1D
	count::Int64
	N::Int64
	highway::Vector{Vehicle}
end

function Highway1D(N::Int64)
    highway = Vehicle[]
    for i in 1:N
        cell = Vehicle(round(Int8, -1), i, round(Int8, -1), zero(Int8), 0, 0)
        push!(highway, cell)
    end
    return Highway1D(0, N, highway)
end

type Highway2D
  lanes::Int8
  highway::Vector{Highway1D}
end

function Highway2D(lanes::Int64, N::Int64)
    C = Highway1D[]
    if lanes > 1
        for i = 1:lanes
            push!(C, Highway1D(N))
        end
        return Highway2D(round(Int8,lanes), C)
    end
    if lanes == 1
        return Highway1D(N)
    end
    if lanes < 1
      println("There must be at least one lane!")
    end
end


text:
"""
    reliability(probabilities, frequencies; deviation=true, kwargs...)
    reliability(probabilities, frequencies, low_high; deviation=true, kwargs...)
    reliability(probabilities, frequencies, low, high; deviation=true, kwargs...)

Plot a reliability diagram of the observed `frequencies` versus the predicted
`probabilities`, optionally with consistency bars ranging from `low` to `high` for
each probability.

If `deviation` is `true` (default), the difference `frequencies - probabilities` is
plotted versus `probabilities`. This can be helpful for inspecting how closely the
observed frequencies match the predicted probabilities.

## Attributes

$(Makie.ATTRIBUTES)

# Examples

```jldoctest
julia> probabilities = sort(rand(10));

julia> frequencies = rand(10);

julia> reliability(probabilities, frequencies);

julia> # plot frequencies instead of deviations
       reliability(probabilities, frequencies; deviation=false);

julia> # use different colors for markers
       reliability(probabilities, frequencies; markercolor=:red);

julia> # random consistency bars
       low = max.(0, probabilities .- rand.());

julia> high = min.(1, probabilities .+ rand.());

julia> # plot pre-computed consistency bars
       reliability(probabilities, frequencies, low, high);

julia> # alternatively consistency bars can be specified as vector of tuples
       reliability(probabilities, frequencies, map(tuple, low, high));
```

# References

Bröcker, J. and Smith, L.A. (2007). *Increasing the reliability of reliability diagrams*.
Weather and forecasting, 22(3), pp. 651-661.
"""
@recipe(Reliability, x_y_low_high) do scene
    return Makie.Attributes(;
        color=Makie.theme(scene, :linecolor),
        colormap=Makie.theme(scene, :colormap),
        colorrange=Makie.automatic,
        linestyle=Makie.theme(scene, :linestyle),
        linewidth=Makie.theme(scene, :linewidth),
        visible=Makie.theme(scene, :visible),
        marker=Makie.theme(scene, :marker),
        markercolor=Makie.theme(scene, :markercolor),
        markercolormap=Makie.theme(scene, :colormap),
        markercolorrange=Makie.automatic,
        markersize=Makie.theme(scene, :markersize),
        markervisible=Makie.theme(scene, :visible),
        strokecolor=Makie.theme(scene, :markerstrokecolor),
        strokewidth=Makie.theme(scene, :markerstrokewidth),
        deviation=Makie.Node(true),
        barswhiskerwidth=0,
        barscolor=Makie.theme(scene, :linecolor),
        barscolormap=Makie.theme(scene, :colormap),
        barscolorrange=Makie.automatic,
        barslinewidth=Makie.theme(scene, :linewidth),
        barsvisible=Makie.theme(scene, :visible),
        cycle=[:color],
        inspectable=Makie.theme(scene, :inspectable),
    )
end

"""
    reliability(
        probabilities::AbstractVector{<:Real},
        outcomes::AbstractVector{Bool};
        binning=EqualMass(),
        consistencybars=ConsistencyBars(),
        kwargs...,
    )

Plot a reliability diagram of the observed frequencies versus the mean probabilities
in different clusters together with a set of consistency bars.

Optionally, the binning algorithm that is used to form the clusters and the consistency
bars can be configured with the `binning` and `consistencybars` keyword arguments. If
`consistencybars === nothing`, then no consistency bars are computed.

# Examples

```jldoctest
julia> probabilities = rand(10);

julia> outcomes = rand(Bool, 10);

julia> # default binning and consistency bars
       reliability(probabilities, outcomes);

julia> # custom options: without consistency bars
       reliability(probabilities, outcomes; consistencybars=nothing);
```

See also: [`EqualMass`](@ref), [`EqualSize`](@ref), [`ConsistencyBars`](@ref)
"""
reliability(::AbstractVector{<:Real}, ::AbstractVector{Bool})

function Makie.plot!(plot::Reliability)
    # extract
    x_y_low_high = plot.x_y_low_high
    deviation = plot.deviation

    # compute pairs of mean predictions and (deviations of) mean outcomes
    x_y = Makie.lift(x_y_low_high, deviation) do x_y_low_high, deviation
        return map(x_y_low_high) do (x, y, _, _)
            return deviation ? Makie.Point2f0(x, y - x) : Makie.Point2f0(x, y)
        end
    end

    # compute pairs of mean predictions and (deviations of) consistency bounds
    x_low_high = Makie.lift(x_y_low_high, deviation) do x_y_low_high, deviation
        return map(x_y_low_high) do (x, _, low, high)
            return if deviation
                Makie.Point3f0(x, low - x, high - x)
            else
                Makie.Point3f0(x, low, high)
            end
        end
    end

    # plot lines
    Makie.lines!(
        plot,
        x_y;
        color=plot.color,
        colormap=plot.colormap,
        colorrange=plot.colorrange,
        linestyle=plot.linestyle,
        linewidth=plot.linewidth,
        visible=plot.visible,
    )

    # plot consistency bars
    Makie.rangebars!(
        plot,
        x_low_high;
        color=plot.barscolor,
        colormap=plot.barscolormap,
        linewidth=plot.barslinewidth,
        visible=plot.barsvisible,
        whiskerwidth=plot.barswhiskerwidth,
    )

    # plot points
    Makie.scatter!(
        plot,
        x_y;
        color=plot.markercolor,
        colormap=plot.colormap,
        colorrange=plot.colorrange,
        marker=plot.marker,
        markersize=plot.markersize,
        strokecolor=plot.strokecolor,
        strokewidth=plot.strokewidth,
        visible=plot.markervisible,
    )

    return plot
end


text:
println("# Test Distributions")

# Multinomial-Dirichlet
println(" * Dirichlet - Multinomial")

d = DirichletMultinomial(5, 1.0)

println(" * - add data")

add_data!(d, [2, 3, 3])

@test d.n == 3
@test d.counts[2] == 1
@test d.counts[3] == 2

println(" * - add sparse data")

add_data!(d, sparse([1,0,1,0,0]))

@test d.n == 4
@test d.counts[1] == 1
@test d.counts[3] == 3

println(" * - remove data")

remove_data!(d, [3])

@test d.n == 3
@test d.counts[3] == 2

println(" * - remove sparse data")

remove_data!(d, sparse([1,0,0,0,0]))

@test d.n == 2
@test d.counts[1] == 0

println(" * - log pdf")
@test logpred(d, [3])[1] == logpred(d, sparse([0,0,1,0,0]))[1]

# student-t Distribution
println(" * Generalized student-t Distribution")

df = 1
mean = 0
sigma = 1

d = Distributions.TDist(df)
@test_approx_eq_eps BNP.tlogpdf(0, df, mean, sigma) Distributions.logpdf(d, 0) 1e-4

df = 1
mean = 1
sigma = 1

d = Distributions.TDist(df)
@test_approx_eq_eps BNP.tlogpdf(0, df, mean, sigma) Distributions.logpdf(d, 0 - mean) 1e-4

df = 1
mean = 0
sigma = 10

d = Distributions.TDist(df)

# Normal-Gamma
println(" * Gamma-Normal")

d = GammaNormal()

(" * - add data")

x = 1.0

add_data!(d, x)
@test d.n == 1

add_data!(d, vec([1.0, 2.0, 1.5, 1.0]))
@test d.n == 5

println(" * - remove data")

remove_data!(d, 2.0)
@test d.n == 4

println(" * - log pdf")

x2 = [2.0]
@test logpred(d, x)[1] > logpred(d, x2)[1]

# Normal-Gamma
println(" * GaussianDiagonal")

d = GaussianDiagonal{NormalNormal}(NormalNormal[NormalNormal() for i in 1:2])

println(" * - add data")
add_data!(d, rand(2))
@test d.n == 1

println("# FINISHED")


text:
# function visualize!(vis, model, q;
#        Δt = 0.1,
# 	   color = RGBA(1,0,0,1.0))
#
# 	default_background!(vis)
#
# 	l2 = Cylinder(Point3f0(-model.l * 4.0, 0.0, 0.0),
# 		Point3f0(model.l * 4.0, 0.0, 0.0),
# 		convert(Float32, 0.025))
#
# 	setobject!(vis["slider"], l2, MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))
#
#     l1 = Cylinder(Point3f0(0.0, 0.0, 0.0),
# 		Point3f0(0.0, 0.0, model.l),
# 		convert(Float32, 0.025))
#
#     setobject!(vis["arm"], l1,
# 		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))
#
#     setobject!(vis["base"], HyperSphere(Point3f0(0.0),
#         convert(Float32, 0.1)),
#         MeshPhongMaterial(color = color))
#
#     setobject!(vis["ee"], HyperSphere(Point3f0(0.0),
#         convert(Float32, 0.05)),
#         MeshPhongMaterial(color = color))
#
# 	anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))
#
# 	for t = 1:length(q)
# 	    MeshCat.atframe(anim,t) do
# 			x = q[t]
# 			px = x[1] + model.l * sin(x[2])
# 			pz = -model.l * cos(x[2])
#             settransform!(vis["arm"], cable_transform([x[1]; 0;0], [px; 0.0; pz]))
#             settransform!(vis["base"], Translation([x[1]; 0.0; 0.0]))
#             settransform!(vis["ee"], Translation([px; 0.0; pz]))
# 	    end
# 	end
#
# 	settransform!(vis["/Cameras/default"],
# 		compose(Translation(0.0, 0.0, -1.0), LinearMap(RotZ(- pi / 2))))
#
# 	MeshCat.setanimation!(vis,anim)
# end

function _create_cartpole!(vis, model;
        i = 0,
        tl = 1.0,
	    color = RGBA(0, 0, 0, tl))

	l2 = Cylinder(Point3f0(-model.l * 10.0, 0.0, 0.0),
		Point3f0(model.l * 10.0, 0.0, 0.0),
		convert(Float32, 0.0125))

	setobject!(vis["slider_$i"], l2, MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, tl)))

    l1 = Cylinder(Point3f0(0.0, 0.0, 0.0),
		Point3f0(0.0, 0.0, model.l),
		convert(Float32, 0.025))

    setobject!(vis["arm_$i"], l1,
		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, tl)))

    setobject!(vis["base_$i"], HyperSphere(Point3f0(0.0),
        convert(Float32, 0.1)),
        MeshPhongMaterial(color = color))

    setobject!(vis["ee_$i"], HyperSphere(Point3f0(0.0),
        convert(Float32, 0.05)),
        MeshPhongMaterial(color = color))
end

function _set_cartpole!(vis, model, x;
    i = 0)

	px = x[1] + model.l * sin(x[2])
	pz = -model.l * cos(x[2])
    settransform!(vis["arm_$i"], cable_transform([x[1]; 0;0], [px; 0.0; pz]))
    settransform!(vis["base_$i"], Translation([x[1]; 0.0; 0.0]))
    settransform!(vis["ee_$i"], Translation([px; 0.0; pz]))
end

function visualize!(vis, model, q;
        i = 0,
        tl = 1.0,
        Δt = 0.1,
        color = RGBA(0,0,0,1.0))

	default_background!(vis)
    _create_cartpole!(vis, model, i = i, color = color, tl = tl)

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

	for t = 1:length(q)
	    MeshCat.atframe(anim,t) do
			_set_cartpole!(vis, model, q[t], i = i)
	    end
	end

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

	MeshCat.setanimation!(vis,anim)
end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
using LinearAlgebra
using BenchmarkTools
ny = 1000
nx = 100
Sₐ = (Diagonal(abs.(rand(nx))));
Sₑ = (Diagonal((3rand(ny))));
#Sₐ = (abs.(rand(nx,nx)));
#Sₑ = ((10*rand(ny,ny)));
K = rand(ny,nx);
y = randn(ny);

prob = OptimalEstimation.OptimalEstimationProblem(K, Sₑ, y, Sₐ)

#OptimalEstimation.solve(prob) ≈ OptimalEstimation.solve_stable(prob)
@time OptimalEstimation.solve(prob, errorAnalysis=true) ;
@time OptimalEstimation.solve_stable(prob,errorAnalysis=true);

max_stars_repo_name: RadiativeTransfer/OptimalEstimation

max_issues_repo_name: RadiativeTransfer/OptimalEstimation

max_forks_repo_name: RadiativeTransfer/OptimalEstimation


text:
# Mostly a blunt copy from CometLogger.jl
# https://github.com/rejuvyesh/CometLogger.jl/blob/master/src/CometLogger.jl
CoreLogging.catch_exceptions(lg::WandbLogger) = false

CoreLogging.min_enabled_level(lg::WandbLogger) = lg.min_level

# For now, log everything that is above the lg.min_level
CoreLogging.shouldlog(lg::WandbLogger, level, _module, group, id) = true

function preprocess(name, val::T, data) where {T}
    if isstructtype(T) && !(val isa PyObject)
        fn = logable_propertynames(val)
        for f in fn
            prop = getfield(val, f)
            preprocess(name * "/$f", prop, data)
        end
    else
        push!(data, name => val)
    end
    return data
end

"""
    logable_propertynames(val::Any)
Returns a tuple with the name of the fields of the structure `val` that
should be logged to Wandb. This function should be overridden when
you want Wandb to ignore some fields in a structure when logging
it. The default behaviour is to return the  same result as `propertynames`.
See also: [`Base.propertynames`](@ref)
"""
logable_propertynames(val::Any) = propertynames(val)

## Default unpacking of key-value dictionaries
function preprocess(name, dict::AbstractDict, data)
    for (key, val) in dict
        # convert any key into a string, via interpolating it
        preprocess("$name/$key", val, data)
    end
    return data
end

# Split complex numbers into real/complex pairs
preprocess(name, val::Complex, data) = push!(data, name * "/re" => real(val), name * "/im" => imag(val))

process(lg::WandbLogger, name::AbstractString, obj, step::Int) = log(lg, Dict(name => obj); step=step)

function CoreLogging.handle_message(lg::WandbLogger, level, message, _module, group, id, file, line; kwargs...)
    i_step = lg.step_increment # :log_step_increment default value

    if !isempty(kwargs)
        data = Vector{Pair{String,Any}}()
        # ∀ (k-v) pairs, decompose values into objects that can be serialized
        for (key, val) in pairs(kwargs)
            # special key describing step increment
            if key == :log_step_increment
                i_step = val
                continue
            end
            preprocess(message * "/$key", val, data)
        end
        iter = increment_step!(lg, i_step)
        for (name, val) in data
            process(lg, name, val, iter)
        end
    end
end


text:
if !isdefined(Main, :CustomStencil)
    include("../module/CustomStencil.jl")
    using Main.CustomStencil
end
include("misc.jl")
coefs = [-5 0.1875 0.0625]

dense_stencil = @def_stencil_expression(
            @sum(i,-2,2,
                @sum(j,-2,2,
                    @sum(k,-2,2, c[max(abs(i),abs(j),abs(k))]*D[x+i,y+j,z+k]))))

st_inst = NewStencilInstance(dense_stencil, 2, coefs)

radius = 2
nx = 8
ny = 8
nz = 8
dx = 1<<(nx)
dy = 1<<(ny)
dz = 1<<(nz)

t_steps = 5
## Actual time steps
data = CreateData(dx,dy,dz,radius)
dout = ApplyStencil(st_inst, data, t_steps)
#heatmap(dout1[:,:,Int(round(dz/2))])
o = Array(dout)

isdefined

isdefined(Main, :CustomStencil)


text:
function environment(
        f::Function;
        globalconfig::Vector=[],
        threadconfig::Vector=[]
    )
    # Save the current settings
    #
    # CPLGetConfigOption() will return the value of the config option, be it
    #     either defined through environment variable, CPLSetConfigOption() or
    #     CPLSetThreadLocalConfigOption() (from the same thread).
    # CPLGetThreadLocalConfigOption() will return the value of the config
    #     option, but only if it has been set with
    #     CPLSetThreadLocalConfigOption()
    #
    # (ref https://github.com/mapbox/rasterio/pull/997#issuecomment-287117289)
    globalsettings = Dict(k => getconfigoption(k) for (k,v) in globalconfig)
    localsettings = Dict(k => getthreadconfigoption(k) for (k,v) in threadconfig)
    for (k,v) in threadconfig; setthreadconfigoption(k, v) end
    for (k,v) in globalconfig; setconfigoption(k, v) end

    try
        f()
    finally
        # Restore previous settings
        for (k,v) in globalsettings
            if v == ""
                clearconfigoption(k)
            else
                setconfigoption(k, v)
            end
        end
        for (k,v) in localsettings
            if v == ""
                clearthreadconfigoption(k)
            else
                setthreadconfigoption(k, v)
            end
        end
    end
end

function executesql(f::Function, dataset::Dataset, args...)
    result = unsafe_executesql(dataset, args...)
    try
        f(result)
    finally
        releaseresultset(dataset, result)
    end
end

function createfeature(f::Function, featuredefn::FeatureDefn)
    feature = unsafe_createfeature(featuredefn)
    reference(featuredefn)
    # the additional reference & dereference here is to deal with the case
    # where you start by (1) creating a featuredefn (0 references)
    # before (2) using it to create a feature here.
    # if we do not artificially increase the reference, then destroy(feature)
    # will release the featuredefn, when we're going to handle it ourselves
    # later. Therefore we dereference (rather than release) the featuredefn.
    try
        f(feature)
    finally
        destroy(feature)
        dereference(featuredefn)
    end
end

"""
    addfielddefn!(layer::AbstractFeatureLayer, name, etype::OGRFieldType; <keyword arguments>)

Create a new field on a layer.

This function should not be called while there are feature objects in existence
that were obtained or created with the previous layer definition.

Not all drivers support this function. You can query a layer to check if it
supports it with the OLCCreateField capability. Some drivers may only support
this method while there are still no features in the layer. When it is
supported, the existing features of the backing file/database should be updated
accordingly.

Drivers may or may not support not-null constraints. If they support creating
fields with not-null constraints, this is generally before creating any feature
to the layer.

### Parameters
* `layer`:  the layer to write the field definition.
* `name`:   name of the field definition to write to disk.
* `etype`:  type of the field definition to write to disk.

### Keyword arguments
* `nwidth`:     the preferred formatting width. 0 (default) indicates undefined.
* `nprecision`: number of decimals for formatting. 0 (default) for undefined.
* `justify`:    the formatting justification ([OJUndefined], OJLeft or OJRight)
* `approx`:     If `true` (default `false`), the field may be created in a
                slightly different form depending on the limitations of the
                format driver.
"""
function addfielddefn!(
        layer::AbstractFeatureLayer,
        name::AbstractString,
        etype::OGRFieldType;
        nwidth::Integer             = 0,
        nprecision::Integer         = 0,
        justify::OGRJustification   = GDAL.OJUndefined,
        approx::Bool                = false
    )
    fielddefn = unsafe_createfielddefn(name, etype)
    setparams!(fielddefn, name, etype, nwidth = nwidth, nprecision = nprecision,
        justify = justify)
    addfielddefn!(layer, fielddefn)
    destroy(fielddefn)
    layer
end

function addfielddefn(
        f::Function,
        layer::AbstractFeatureLayer,
        name::AbstractString,
        etype::OGRFieldType;
        nwidth::Integer             = 0,
        nprecision::Integer         = 0,
        justify::OGRJustification   = GDAL.OJUndefined,
        approx::Bool                = false
    )
    fielddefn = unsafe_createfielddefn(name, etype)
    setparams!(fielddefn, name, etype, nwidth = nwidth, nprecision = nprecision,
        justify = justify)
    try
        f(fielddefn)
        addfielddefn!(layer, fielddefn)
    finally
        destroy(fielddefn)
    end
end

"""
    writegeomdefn!(layer::AbstractFeatureLayer, name, etype::OGRwkbGeometryType, approx=false)

Write a new geometry field on a layer.

This function should not be called while there are feature objects in existence
that were obtained or created with the previous layer definition.

Not all drivers support this function. You can query a layer to check if it
supports it with the OLCCreateField capability. Some drivers may only support
this method while there are still no features in the layer. When it is
supported, the existing features of the backing file/database should be updated
accordingly.

Drivers may or may not support not-null constraints. If they support creating
fields with not-null constraints, this is generally before creating any feature
to the layer.

### Parameters
* `layer`:  the layer to write the field definition.
* `name`:   name of the field definition to write to disk.
* `etype`:  type of the geometry field defintion to write to disk.

### Keyword arguments
* `approx`: If `true` (default `false`), the geometry field may be created in a
            slightly different form depending on the limitations of the driver.
"""
function writegeomdefn!(
        layer::AbstractFeatureLayer,
        name::AbstractString,
        etype::OGRwkbGeometryType;
        approx::Bool = false
    )
    geomdefn = unsafe_creategeomdefn(name, etype)
    addgeomdefn!(layer, geomdefn)
    destroy(geomdefn)
    return layer
end

function writegeomdefn(
        f::Function,
        layer::AbstractFeatureLayer,
        name::AbstractString,
        etype::OGRwkbGeometryType;
        approx::Bool = false
    )
    geomdefn = unsafe_creategeomdefn(name, etype)
    try
        f(geomdefn)
        addgeomdefn!(layer, geomdefn)
    finally
        destroy(geomdefn)
    end
end

for gdalfunc in (
        :boundary, :buffer, :centroid, :clone, :convexhull, :create,
        :createcolortable, :createcoordtrans, :copy, :createfeaturedefn,
        :createfielddefn, :creategeom, :creategeomcollection,
        :creategeomfieldcollection, :creategeomdefn, :createlayer,
        :createlinearring, :createlinestring, :createmultilinestring,
        :createmultipoint, :createmultipolygon, :createmultipolygon_noholes,
        :createpoint, :createpolygon, :createRAT, :createstylemanager,
        :createstyletable, :createstyletool, :curvegeom, :delaunaytriangulation,
        :difference, :forceto, :fromGML, :fromJSON, :fromWKB, :fromWKT,
        :gdalbuildvrt, :gdaldem, :gdalgrid, :gdalnearblack, :gdalrasterize,
        :gdaltranslate, :gdalvectortranslate, :gdalwarp, :getband,
        :getcolortable, :getfeature, :getgeom, :getlayer, :getmaskband,
        :getoverview, :getpart, :getspatialref, :importCRS, :intersection, :importEPSG,
        :importEPSGA, :importESRI, :importPROJ4, :importWKT, :importXML,
        :importURL, :lineargeom, :newspatialref, :nextfeature, :pointalongline,
        :pointonsurface, :polygonfromedges, :polygonize, :read, :sampleoverview,
        :simplify, :simplifypreservetopology, :symdifference, :union, :update,
        :readraster,
    )
    eval(quote
        function $(gdalfunc)(f::Function, args...; kwargs...)
            obj = $(Symbol("unsafe_$gdalfunc"))(args...; kwargs...)
            try
                f(obj)
            finally
                destroy(obj)
            end
        end
    end)
end


text:
# default choices for g and dg
g(u::T, j::Int64, zl = 0.0, zr = 1.0) where {T<:Real} =  - ( u + zl )^ j * (zr -u)
dg(u::T, j::Int64, zl = 0.0, zr = 1.0) where {T<:Real} = ( u + zl )^(j-1) * ( (j+1) * u + zl - j * zr )


text:
"""
    remove_duplicates(x::Vector, y::Vector)

Remove duplicate tuples `(x[i],y[i])` from the vectors `x` and `y`.
"""
function remove_duplicates(x::Vector, y::Vector)
    points = hcat(x, y)
    unique_points = unique(points, dims = 1)

    return unique_points[:, 1], unique_points[:, 2]
end

"""
    voronoiarea(x::Vector, y::Vector, rw) -> Vector

Compute the area of each Voronoi cell of the generators `(x[i], y[i])` in the vectors `x` and `y`.

`rw` is the boundary window.
"""
function voronoiarea(x::Vector, y::Vector, rw::Vector=[0.0; 1.0; 0.0; 1.0], epsilon = 1e-9)
    da = DeldirArguments(x, y, rw, epsilon)
    deldirwrapper!(da)

    npd = Int64(da.npd[1])
    unordered_area = @view da.dirsum[2*npd + 1:3*npd]

    area = unordered_area[da.indices]
    return area
end

"""
    edges(D) -> Vector, Vector

Collect the edges of a dataframe in vectors that are ready to be plotted.
"""
function edges(D::DataFrame)
    x1 = D[!, "x1"]
    y1 = D[!, "y1"]
    x2 = D[!, "x2"]
    y2 = D[!, "y2"]

    N = DataFrames.nrow(D)
    x = Vector{Float64}(undef, 3*N)
    y = similar(x)

    nx = 0
    for n in 1:N
        nx += 1
    	x[nx] = x1[n]
    	y[nx] = y1[n]

        nx += 1
    	x[nx] = x2[n]
    	y[nx] = y2[n]

        nx += 1
    	x[nx] = NaN
    	y[nx] = NaN
    end

    return x, y
end


text:
module SignalAnalysis

using Requires
using DocStringExtensions

using SignalBase
using SignalBase.Units
export nframes, nchannels, sampletype, framerate, duration
export 𝓈

const 𝓈 = Units.s

include("signals.jl")
include("basic.jl")
include("dsp.jl")
include("tfa.jl")
include("generate.jl")
include("array.jl")
include("rand.jl")

function __init__()
  @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" include("plot.jl")
  @require InteractiveViz="d14badfc-0adb-4d57-980e-37858d990fa5" include("iplot.jl")
end

end # module


text:
"""
    PyRef([x])

A reference to a Python object converted from `x`, or a null reference if `x` is not given.

This is baically just a mutable reference to a pointer to a Python object.
It owns the reference (if non-NULL) and automatically decrefs it when finalized.

Building block for more complex wrapper types such as `PyObject` and `PyList`.
"""
mutable struct PyRef
    ptr::CPyPtr
    PyRef(::Val{:new}, ptr::Ptr, borrowed::Bool) = begin
        x = new(CPyPtr(ptr))
        borrowed && C.Py_IncRef(ptr)
        finalizer(x) do x
            if CONFIG.isinitialized
                ptr = x.ptr
                if !isnull(ptr)
                    with_gil(false) do
                        C.Py_DecRef(ptr)
                    end
                    x.ptr = C_NULL
                end
            end
        end
        x
    end
end
export PyRef

pynewref(x::Ptr, check::Bool = false) =
    (check && isnull(x)) ? pythrow() : PyRef(Val(:new), x, false)
pyborrowedref(x::Ptr, check::Bool = false) =
    (check && isnull(x)) ? pythrow() : PyRef(Val(:new), x, true)

ispyreftype(::Type{PyRef}) = true
pyptr(x::PyRef) = x.ptr
isnull(x::PyRef) = isnull(pyptr(x))
Base.unsafe_convert(::Type{CPyPtr}, x::PyRef) = pyptr(x)
C.PyObject_TryConvert__initial(o, ::Type{PyRef}) = C.putresult(pyborrowedref(o))

PyRef(x) = begin
    ptr = C.PyObject_From(x)
    isnull(ptr) && pythrow()
    pynewref(ptr)
end
PyRef() = pynewref(C_NULL)

Base.convert(::Type{PyRef}, x::PyRef) = x
Base.convert(::Type{PyRef}, x) = PyRef(x)

# Cache some common standard modules
for name in [
    "os",
    "io",
    "sys",
    "pprint",
    "collections",
    "collections.abc",
    "numbers",
    "fractions",
    "datetime",
    "numpy",
    "pandas",
]
    f = Symbol("py", replace(name, "." => ""), "module")
    rf = Symbol("_", f)
    @eval $rf = PyRef()
    @eval $f(::Type{T}) where {T} = begin
        r = $rf
        if isnull(r.ptr)
            m = pyimport(PyRef, $name)
            C.Py_IncRef(m.ptr)
            r.ptr = m.ptr
        end
        (r isa T) ? r : pyconvert(T, r)
    end
    @eval $f() = $f(PyObject)
end