text:
using MosaicViews
using Test
using Colors, ColorVectorSpace

@testset "MosaicView" begin
    @test_throws ArgumentError MosaicView(rand(2))
    @test_throws ArgumentError MosaicView(rand(2,2))
    @test_throws ArgumentError MosaicView(rand(2,2,2,2,2))

    @testset "3D input" begin
        A = zeros(Int,2,2,2)
        A[:,:,1] = [1 2; 3 4]
        A[:,:,2] = [5 6; 7 8]
        mv = @inferred MosaicView(A)
        @test parent(mv) === A
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (4, 2)
        @test @inferred(getindex(mv,1,1)) === 1
        @test @inferred(getindex(mv,2,1)) === 3
        @test_throws BoundsError mv[0,1]
        @test_throws BoundsError mv[1,0]
        @test_throws BoundsError mv[1,3]
        @test_throws BoundsError mv[5,1]
        @test all(mv .== vcat(A[:,:,1],A[:,:,2]))
        # singleton dimension doesn't change anything
        @test mv == MosaicView(reshape(A,2,2,2,1))
    end

    @testset "4D input" begin
        A = zeros(Int,2,2,1,2)
        A[:,:,1,1] = [1 2; 3 4]
        A[:,:,1,2] = [5 6; 7 8]
        mv = @inferred MosaicView(A)
        @test parent(mv) === A
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (2, 4)
        @test @inferred(getindex(mv,1,1)) === 1
        @test @inferred(getindex(mv,2,1)) === 3
        @test_throws BoundsError mv[0,1]
        @test_throws BoundsError mv[1,0]
        @test_throws BoundsError mv[3,1]
        @test_throws BoundsError mv[1,5]
        @test all(mv .== hcat(A[:,:,1,1],A[:,:,1,2]))
        A = zeros(Int,2,2,2,3)
        A[:,:,1,1] = [1 2; 3 4]
        A[:,:,1,2] = [5 6; 7 8]
        A[:,:,1,3] = [9 10; 11 12]
        A[:,:,2,1] = [13 14; 15 16]
        A[:,:,2,2] = [17 18; 19 20]
        A[:,:,2,3] = [21 22; 23 24]
        mv = @inferred MosaicView(A)
        @test parent(mv) === A
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (4, 6)
        @test @inferred(getindex(mv,1,1)) === 1
        @test @inferred(getindex(mv,2,1)) === 3
        @test all(mv .== vcat(hcat(A[:,:,1,1],A[:,:,1,2],A[:,:,1,3]), hcat(A[:,:,2,1],A[:,:,2,2],A[:,:,2,3])))
    end
end

@testset "mosaicview" begin
    @test_throws ArgumentError mosaicview(rand(2))
    @test_throws ArgumentError mosaicview(rand(2,2))

    @testset "3D input" begin
        A = [(k+1)*l-1 for i in 1:2, j in 1:3, k in 1:2, l in 1:2]
        B = reshape(A, 2, 3, :)
        @test_throws ArgumentError mosaicview(B, nrow=0)
        @test_throws ArgumentError mosaicview(B, ncol=0)
        @test_throws ArgumentError mosaicview(B, nrow=1, ncol=1)
        mv = mosaicview(B)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(B)
        @test size(mv) == (8, 3)
        @test @inferred(getindex(mv,3,1)) === 2
        @test mv == [
            1  1  1
            1  1  1
            2  2  2
            2  2  2
            3  3  3
            3  3  3
            5  5  5
            5  5  5
        ]
        mv = mosaicview(B, nrow=2)
        @test mv == MosaicView(A)
        @test typeof(mv) != typeof(MosaicView(A))
        @test parent(parent(mv)).data == B
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(B)
        @test size(mv) == (4, 6)
    end

    @testset "4D input" begin
        A = [(k+1)*l-1 for i in 1:2, j in 1:3, k in 1:2, l in 1:2]
        @test_throws ArgumentError mosaicview(A, nrow=0)
        @test_throws ArgumentError mosaicview(A, ncol=0)
        @test_throws ArgumentError mosaicview(A, nrow=1, ncol=1)
        mv = mosaicview(A)
        @test mv == MosaicView(A)
        @test typeof(mv) != typeof(MosaicView(A))
        @test parent(parent(mv)).data == reshape(A, 2, 3, :)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (4, 6)
        mv = mosaicview(A, npad=1)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (5, 7)
        @test mv == mosaicview(A, nrow=2, npad=1)
        @test mv == mosaicview(A, ncol=2, npad=1)
        @test @inferred(getindex(mv,3,1)) === 0
        @test @inferred(getindex(mv,2,5)) === 3
        @test mv == [
            1  1  1  0  3  3  3
            1  1  1  0  3  3  3
            0  0  0  0  0  0  0
            2  2  2  0  5  5  5
            2  2  2  0  5  5  5
        ]
        mv = mosaicview(A, ncol=3, npad=1)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (5, 11)
        @test mv == [
            1  1  1  0  3  3  3  0  0  0  0
            1  1  1  0  3  3  3  0  0  0  0
            0  0  0  0  0  0  0  0  0  0  0
            2  2  2  0  5  5  5  0  0  0  0
            2  2  2  0  5  5  5  0  0  0  0
        ]
        mv = mosaicview(A, rowmajor=true, ncol=3, npad=1)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (5, 11)
        @test @inferred(getindex(mv,3,1)) === 0
        @test @inferred(getindex(mv,2,5)) === 2
        @test mv == [
            1  1  1  0  2  2  2  0  3  3  3
            1  1  1  0  2  2  2  0  3  3  3
            0  0  0  0  0  0  0  0  0  0  0
            5  5  5  0  0  0  0  0  0  0  0
            5  5  5  0  0  0  0  0  0  0  0
        ]
        mv = mosaicview(A, rowmajor=true, ncol=3, npad=2)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (6, 13)
        @test mv == [
            1  1  1  0  0  2  2  2  0  0  3  3  3
            1  1  1  0  0  2  2  2  0  0  3  3  3
            0  0  0  0  0  0  0  0  0  0  0  0  0
            0  0  0  0  0  0  0  0  0  0  0  0  0
            5  5  5  0  0  0  0  0  0  0  0  0  0
            5  5  5  0  0  0  0  0  0  0  0  0  0
        ]
        mv = mosaicview(A, -1.0, rowmajor=true, ncol=3, npad=1)
        @test typeof(mv) <: MosaicView
        @test eltype(mv) == eltype(A)
        @test size(mv) == (5, 11)
        @test mv == [
             1   1   1  -1   2   2   2  -1   3   3   3
             1   1   1  -1   2   2   2  -1   3   3   3
            -1  -1  -1  -1  -1  -1  -1  -1  -1  -1  -1
             5   5   5  -1  -1  -1  -1  -1  -1  -1  -1
             5   5   5  -1  -1  -1  -1  -1  -1  -1  -1
        ]
    end

    @testset "Colorant Array" begin
        A = rand(RGB{Float32}, 2, 3, 2, 2)
        mv = mosaicview(A)
        @test eltype(mv) == eltype(A)
        @test mv == @inferred(MosaicView(A))
        mv = mosaicview(A, rowmajor=true, ncol=3)
        @test eltype(mv) == eltype(A)
        @test @inferred(getindex(mv, 3, 4)) == RGB(0,0,0)
        mv = mosaicview(A, colorant"white", rowmajor=true, ncol=3)
        @test eltype(mv) == eltype(A)
        @test @inferred(getindex(mv, 3, 4)) == RGB(1,1,1)
    end
end


text:
@testset "growcut3d" begin
    img = generate_3d_img()
    display_3D(img)
    @test true
    clicks=zeros(Int,axes(img))
    clicks[50,50,10] = 1
    clicks[50,50,90] = 1
    clicks[12,12,12] = 2
    clicks[38,38,38] = 2
    clicks[5,5,5] = 3
    clicks[90,90,90] = 3
    clicks[90,5,90] = 3
    clicks[5,90,5] = 3
    l = segment_image(GrowCut(),img,clicks,27,27)
    @test l[50,50,50] == 1
    @test l[55,45,98] == 1
    @test l[25,25,25] == 2
    l = segment_image(GrowCut(),img,clicks,26,26)
    @test l[11,11,30] == 2
    @test l[5,5,90] == 3
    @test l[90,90,5] == 3
end

ext: jl


text:
struct TwoLoopVars{T,TR,M}
    d::T # search direction
    S::TR # change in x
    Y::TR # change in successive gradients
    α::M
    ρ::M
end
function TwoLoopVars(x, memory)
    d = copy(x)

    S = [copy(x) for i = 1:memory]
    Y = [copy(x) for i = 1:memory]
    α = similar(x, memory)
    ρ = similar(x, memory)
    TwoLoopVars(d, S, Y, α, ρ)
end
lbfgs_vars(method, x) = TwoLoopVars(x, method.memory)
function solve(
    problem::OptimizationProblem,
    x0,
    approach::LineSearch{<:LBFGS{<:Inverse,<:TwoLoop},<:Any,<:QNScaling},
    options::OptimizationOptions,
    cache = preallocate_qn_caches(mstyle(problem), x0),
)
    solve(problem, (x0, nothing), approach, options, cache)
end
function solve(
    problem::OptimizationProblem,
    s0::Tuple,
    approach::LineSearch{<:LBFGS{<:Inverse,<:TwoLoop},<:Any,<:QNScaling},
    options::OptimizationOptions,
    cache = preallocate_qn_caches(mstyle(problem), first(s0)),
)
    _solve(mstyle(problem), problem, s0, approach, options, cache)
end

function _solve(
    mstyle,
    problem::OptimizationProblem,
    s0::Tuple,
    approach::LineSearch{<:LBFGS{<:Inverse,<:TwoLoop},<:Any,<:QNScaling},
    options::OptimizationOptions,
    cache,
)
    t0 = time()

    #==============
         Setup
    ==============#
    x0, B0 = s0
    T = eltype(x0)
    # Remove B from all of this
    objvars = prepare_variables(problem, approach, x0, copy(x0), B0)

    qnvars = lbfgs_vars(modelscheme(approach), objvars.x)
    ∇f0 = norm(objvars.∇fz, Inf)
    f0 = objvars.fz
    P = initial_preconditioner(approach, x0)

    #========================
         First iteration
    ========================#
    objvars, qnvars, P =
        iterate(mstyle, 1, qnvars, objvars, P, approach, problem, problem, options)
    iter = 1
    # Check for gradient convergence
    is_converged = converged(approach, objvars, ∇f0, options)
    while iter <= options.maxiter && !any(is_converged)
        iter += 1
        # take a step and update approximation
        objvars, qnvars, P = iterate(
            mstyle,
            iter,
            qnvars,
            objvars,
            P,
            approach,
            problem,
            problem,
            options,
            false,
        )
        # Check for gradient convergence
        is_converged = converged(approach, objvars, ∇f0, options)
    end
    x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
    return ConvergenceInfo(
        approach,
        (
            P = P,
            B = B,
            ρs = norm(x .- z),
            ρx = norm(x),
            solution = z,
            fx = fx,
            minimum = fz,
            ∇fz = ∇fz,
            f0 = f0,
            ∇f0 = ∇f0,
            iter = iter,
            time = time() - t0,
        ),
        options,
    )
end
function iterate(
    mstyle::InPlace,
    iter::Integer,
    qnvars::TwoLoopVars,
    objvars,
    P,
    approach::LineSearch{<:LBFGS{<:Inverse,<:TwoLoop},<:Any,<:QNScaling},
    prob::OptimizationProblem,
    problem::OptimizationProblem,
    options::OptimizationOptions,
    is_first = nothing,
)
    # split up the approach into the hessian approximation scheme and line search
    scheme, linesearch, K = modelscheme(approach), algorithm(approach), approach.scaling # make grabber and get a better name than K
    current_memory = min(iter - 1, scheme.memory)
    x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
    Tf = typeof(fz)

    # Move nexts into currs
    fx = fz
    copyto!(x, z)
    copyto!(∇fx, ∇fz)

    # Update preconditioner
    P = update_preconditioner(scheme, x, P) # returns nothing????

    # Update current gradient and calculate the search direction
    d = find_direction!(scheme, copy(∇fz), qnvars, current_memory, K, P) # solve Bd = -∇fx
    φ = _lineobjective(mstyle, problem, ∇fz, z, x, d, fx, real(dot(∇fx, d))) # real is needed to convert complex dots to actually being real

    # Perform line search along d
    α, f_α, ls_success = find_steplength(mstyle, linesearch, φ, Tf(1))

    s = current_memory == length(qnvars.S) ? qnvars.S[1] : qnvars.S[1+current_memory]
    @. s = α * d
    z = retract(problem, z, x, s)

    # Update approximation
    fz, ∇fz, qnvars =
        update_obj!(problem, qnvars, α, x, ∇fx, z, ∇fz, current_memory, scheme)
    return (
        x = x,
        fx = fx,
        ∇fx = ∇fx,
        z = z,
        fz = fz,
        ∇fz = ∇fz,
        B = nothing,
        Pg = nothing,
    ),
    qnvars,
    P
end

function iterate(
    mstyle::OutOfPlace,
    cache,
    objvars,
    P,
    approach::LineSearch{<:LBFGS{<:Inverse,<:TwoLoop},<:Any,<:QNScaling},
    prob::OptimizationProblem,
    problem::OptimizationProblem,
    options::OptimizationOptions,
    is_first = nothing,
)
    # split up the approach into the hessian approximation scheme and line search
    scheme, linesearch = modelscheme(approach), algorithm(approach)
    # Move nexts into currs
    x, fx, ∇fx, z, fz, ∇fz, B, Pg = objvars
    Tf = typeof(fz)
    fx = fz
    x = copy(z)
    ∇fx = copy(∇fz)
    # Update preconditioner
    P = update_preconditioner(scheme, x, P)
    # Update current gradient and calculate the search direction
    d = find_direction(B, P, ∇fx, scheme) # solve Bd = -∇fx
    φ = _lineobjective(mstyle, problem, problem, ∇fz, z, x, d, fx, dot(∇fx, d))

    # # Perform line search along d
    α, f_α, ls_success = find_steplength(linesearch, φ, Tf(1))

    # # Calculate final step vector and update the state
    s = @. α * d
    z = retract(problem, z, x, s)

    # Update approximation
    fz, ∇fz, B = update_obj(problem, s, ∇fx, z, ∇fz, B, scheme, is_first)

    return x, fx, ∇fx, z, fz, ∇fz, B, P
end

max_stars_repo_name: JuliaNLSolvers/NLSolvers

max_issues_repo_name: JuliaNLSolvers/NLSolvers

max_forks_repo_name: JuliaNLSolvers/NLSolvers


text:
export tmapfoldl;

"""
    tmapfoldl(f::Function, op::Function, itr::AbstractArray; init)::AbstractArray

Multi-threaded version of [mapfoldl(f, op, itr; init)](https://docs.julialang.org/en/v1.2/base/collections/#Base.mapfoldl-Tuple{Any,Any,Any}). Only the **mapping** part is parallelised.
"""
function tmapfoldl(f::Function, op::Function, itr::T; init) where T<:AbstractArray
	return foldl(op, tmap(f, itr); init = init);
end

"""
    tmapfoldl(f::Function, op::Function, itr::AbstractArray)::AbstractArray

Multi-threaded version of [mapfoldl(f, op, itr)](https://docs.julialang.org/en/v1.2/base/collections/#Base.mapfoldl-Tuple{Any,Any,Any}). Only the **mapping** part is parallelised.
"""
function tmapfoldl(f::Function, op::Function, itr::T) where T<:AbstractArray
	return foldl(op, tmap(f, itr));
end


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/MLJTutorials/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("MLJTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

using MLJ

@load LinearRegressor pkg=MLJLinearModels

using RDatasets, DataFrames
boston = dataset("MASS", "Boston")
first(boston, 3)

describe(boston, :mean, :std, :eltype)

data = coerce(boston, autotype(boston, :discrete_to_continuous));

y = data.MedV
X = select(data, Not(:MedV));

mdl = LinearRegressor()

X_uni = select(X, :LStat) # only a single feature
mach_uni = machine(mdl, X_uni, y)
fit!(mach_uni)

fp = fitted_params(mach_uni)
@show round.(fp.coefs, sigdigits=3)
@show round(fp.intercept, sigdigits=3)

using PyPlot

figure(figsize=(8,6))
plot(X.LStat, y, ls="none", marker="o")
Xnew = (LStat = collect(range(extrema(X.LStat)..., length=100)),)
plot(Xnew.LStat, predict(mach_uni, Xnew))

mach = machine(mdl, X, y)
fit!(mach)

fp = fitted_params(mach)
@show round.(fp.coefs[1:3], sigdigits=3)
@show round(fp.intercept, sigdigits=3)

println(rpad(" Feature", 11), "| ", "Coefficient")
println("-"^24)
for (i, name) in enumerate(names(X))
    println(rpad("$name", 11), "| ", round(fp.coefs[i], sigdigits=3))
end
println(rpad("Intercept", 11), "| ", round(fp.intercept, sigdigits=3))

ŷ = predict(mach, X)
round(rms(ŷ, y), sigdigits=4)

figure(figsize=(8,6))
res = ŷ .- y
stem(res)

figure(figsize=(8,6))
hist(res, density=true)

X2 = hcat(X, X.LStat .* X.Age);

rename!(X2, :x1 => :interaction);

mach = machine(mdl, X2, y)
fit!(mach)
ŷ = predict(mach, X2)
round(rms(ŷ, y), sigdigits=4)

X3 = hcat(X.LStat, X.LStat.^2)
mach = machine(mdl, X3, y)
fit!(mach)
ŷ = predict(mach, X3)
round(rms(ŷ, y), sigdigits=4)

Xnew = (LStat = Xnew.LStat, LStat2 = Xnew.LStat.^2)

figure(figsize=(8,6))
plot(X.LStat, y, ls="none", marker="o")
plot(Xnew.LStat, predict(mach, Xnew))

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials


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

using Markdown
using InteractiveUtils

# ╔═╡ d1f4afa8-2501-4447-96a6-7fb93f463e2c
using IterativeSolvers, DataFrames

# ╔═╡ 57c64cca-55ee-4c2b-ad27-2942a262df7b
md" ### Ex-02-08.jl"

# ╔═╡ c7b3e4fb-8f17-4e0e-82a5-18d30f8b198a
begin
	A = Float64[
	  16 4 8;
	  4 5 -4;
	  8 -4 22;
	]
	a = Float64[4; 2; 5]
	r = vcat(collect(1:6), collect(10:10:100))
end;

# ╔═╡ 21a0af0a-7873-45f5-bd59-31c43cbb0f0c
md" ##### Jacobi method:"

# ╔═╡ 137e76a3-ea70-49e0-a1f7-a2eadaac9963
x_j = [(i=i, j=jacobi(A, a; maxiter=i)) for i in r]

# ╔═╡ 31d9aa57-e96b-409b-b1ab-9e5966af6bf5
x_j[1]

# ╔═╡ 70d35d30-1798-4982-9b9e-3812a575d179
A * x_j[11].j

# ╔═╡ 954c3dec-86ad-4a29-9528-a790a9ebf2c0
md" ##### Gauss-Seidel method:"

# ╔═╡ 53491c26-f5d1-48d0-a795-7c20a31b00df
x_gs = [(i=i, gs=gauss_seidel(A, a; maxiter=i)) for i in r]

# ╔═╡ b2a3e4de-cd59-466a-9872-f72e78bb91ba
A * x_gs[11].gs

# ╔═╡ df81e70f-ec63-4be6-8517-d962b315dc0f
md" ##### Succesive Over-Relaxation method:"

# ╔═╡ cebd114d-8a7c-44c0-bbca-e39b124a263e
x_sor = [(i=i, sor=sor(A, a, 1.0; maxiter=i)) for i in r]

# ╔═╡ 9b86fca1-6867-4c70-b7a2-6e648bbfc26e
A * x_sor[11].sor

# ╔═╡ c1dd5a8b-5256-410f-9abc-58afd36cb1af
md" ##### Sequential Succesive Over-Relaxation method:"

# ╔═╡ 09fd7541-fe59-435b-8f6d-73b72b0c01f0
x_ssor = [(i=i, ssor=ssor(A, a, 1.0; maxiter=i)) for i in r]

# ╔═╡ 96d4ef12-7a35-47c5-a255-5056e582e618
A * x_ssor[11].ssor

# ╔═╡ Cell order:
# ╟─57c64cca-55ee-4c2b-ad27-2942a262df7b
# ╠═d1f4afa8-2501-4447-96a6-7fb93f463e2c
# ╠═c7b3e4fb-8f17-4e0e-82a5-18d30f8b198a
# ╟─21a0af0a-7873-45f5-bd59-31c43cbb0f0c
# ╠═137e76a3-ea70-49e0-a1f7-a2eadaac9963
# ╠═31d9aa57-e96b-409b-b1ab-9e5966af6bf5
# ╠═70d35d30-1798-4982-9b9e-3812a575d179
# ╟─954c3dec-86ad-4a29-9528-a790a9ebf2c0
# ╠═53491c26-f5d1-48d0-a795-7c20a31b00df
# ╠═b2a3e4de-cd59-466a-9872-f72e78bb91ba
# ╟─df81e70f-ec63-4be6-8517-d962b315dc0f
# ╠═cebd114d-8a7c-44c0-bbca-e39b124a263e
# ╠═9b86fca1-6867-4c70-b7a2-6e648bbfc26e
# ╟─c1dd5a8b-5256-410f-9abc-58afd36cb1af
# ╠═09fd7541-fe59-435b-8f6d-73b72b0c01f0
# ╠═96d4ef12-7a35-47c5-a255-5056e582e618


text:
using Test
using CUDA
using Random
using Rectangle

include("../environments/dummy_app.jl")

function kernel_test_is_point_in_rect(
    point_x,
    point_y,
    rect_x,
    rect_y,
    width,
    height,
    result,
)

    threadID = threadIdx().x

    result[threadID] = is_point_in_rect(point_x, point_y, rect_x, rect_y, width, height)

    return
end

@testset "Point in Rect" begin

    app_width = 600
    app_height = 400

    result_gpu = CUDA.fill(true, 1)

    for i = 1:10000

        point_x = rand(1:app_width)
        point_y = rand(1:app_width)

        rect_x = rand(1:app_width)
        rect_y = rand(1:app_width)
        rect_width = rand(100:400)
        rect_height = rand(50:200)

        CUDA.@cuda threads = 1 blocks = 1 kernel_test_is_point_in_rect(
            point_x,
            point_y,
            rect_x,
            rect_y,
            rect_width,
            rect_height,
            result_gpu,
        )
        CUDA.synchronize()

        result = Array(result_gpu)

        rect = Rect(rect_x, rect_y, rect_x + rect_width, rect_y + rect_height)
        @test inside((point_x, point_y), rect) == result[1]

    end

end


text:
using Test, LinearMaps, LinearAlgebra, SparseArrays, BenchmarkTools

@testset "basic functionality" begin
    A = 2 * rand(ComplexF64, (20, 10)) .- 1
    v = rand(ComplexF64, 10)
    w = rand(ComplexF64, 20)
    V = rand(ComplexF64, 10, 3)
    W = rand(ComplexF64, 20, 3)
    α = rand()
    β = rand()
    M = @inferred LinearMap(A)
    N = @inferred LinearMap(M)

    @testset "LinearMaps.jl" begin
        @test eltype(M) == eltype(A)
        @test size(M) == size(A)
        @test size(N) == size(A)
        @test !isreal(M)
        @test ndims(M) == 2
        @test_throws ErrorException size(M, 3)
        @test length(M) == length(A)
    end

    Av = A * v
    AV = A * V

    @testset "mul! and *" begin
        @test M * v == Av
        @test N * v == Av
        @test @inferred mul!(copy(w), M, v) == mul!(copy(w), A, v)
        b = @benchmarkable mul!($w, $M, $v)
        @test run(b, samples=3).allocs == 0
        @test @inferred mul!(copy(w), N, v) == Av

        # mat-vec-mul
        @test @inferred mul!(copy(w), M, v, 0, 0) == zero(w)
        @test @inferred mul!(copy(w), M, v, 0, 1) == w
        @test @inferred mul!(copy(w), M, v, 0, β) == β * w
        @test @inferred mul!(copy(w), M, v, 1, 1) ≈ Av + w
        @test @inferred mul!(copy(w), M, v, 1, β) ≈ Av + β * w
        @test @inferred mul!(copy(w), M, v, α, 1) ≈ α * Av + w
        @test @inferred mul!(copy(w), M, v, α, β) ≈ α * Av + β * w
        @test @inferred mul!(copy(w), M, v, α)    ≈ α * Av

        # test mat-mat-mul!
        @test @inferred mul!(copy(W), M, V, α, β) ≈ α * AV + β * W
        @test @inferred mul!(copy(W), M, V, α) ≈ α * AV
        @test @inferred mul!(copy(W), M, V) ≈ AV
        @test typeof(M * V) <: LinearMap
    end
end


text:
import Aswing
using Test

mydir,_ = splitdir(@__FILE__)

# tolerance used for tests
atol = 1e-5

# Read and input .asw file
inputasw = joinpath(mydir,"hawk.asw")
asw = Aswing.readasw(inputasw)
Aswing.setgeom(asw)

# Construct .pnt file from scratch
constraints = Aswing.Constraints(
    linear_acceleration = [Aswing.IPFORX,Aswing.IPFORY,Aswing.IPFORZ],
    angular_acceleration = [Aswing.IPMOMX,Aswing.IPMOMY,Aswing.IPMOMZ],
    velocity = [Aswing.IPVIAS,Aswing.IPBETA,Aswing.IPVACZ],
    rotation_rate = [Aswing.IPROTX,Aswing.IPROTY,Aswing.IPROTZ],
    position = [Aswing.IPPOSX,Aswing.IPPOSY,Aswing.IPPOSZ],
    phi = Aswing.IPBANK,
    theta = Aswing.IPVACX,
    psi = Aswing.IPHEAD,
    flap_defl_ctrl_var = [Aswing.IPRACX, Aswing.IPRACY, Aswing.IPRACZ, Aswing.IPFLP4],
    eng_pwr_ctrl_var = [Aswing.IPENG1])

parameters = Aswing.Parameters(
    linear_acceleration = zeros(3),
    angular_acceleration = zeros(3),
    velocity = 60.0,
    beta = 0.0,
    alpha = 1.968,
    rotation_rate = zeros(3),
    position = zeros(3),
    phi = 0.0,
    theta = 0.29,
    psi = 0.0,
    flap_defl_ctrl_var = [0.0, -1.8050, 0.0, 0.0],
    eng_pwr_ctrl_var = [0.0],
    sum_force = zeros(3),
    sum_mom = zeros(3),
    lift = 332.158,
    climb_angle = -1.678,
    radial_acceleration = 0.0,
    usr1 = 0.0,
    usr2 = 0.0,
    least_squared = 0.0)

pnt = [Aswing.OperatingPoint(
    constraints = constraints,
    parameters = parameters,
    mach_from_airspeed = false,
    machpg = 0.0)]

Aswing.setpnt(pnt)

# Save .asw and .pnt files
outputasw = joinpath(mydir,"out.asw")
Aswing.aswsav(outputasw)

outputpnt = joinpath(mydir,"out.pnt")
Aswing.pntsav(outputpnt)

# test reading/writing .asw files
newasw = Aswing.readasw(outputasw)
@testset "asw" begin
    @test asw.namein == newasw.namein
    @test asw.si_units == newasw.si_units
    @test asw.constants == newasw.constants
    @testset "reference values" begin
        for field in fieldnames(Aswing.ReferenceValues)
            @test isapprox(getfield(asw.refvals, field), getfield(newasw.refvals, field))
        end
    end
    @testset "weights" begin
        for i = 1:length(asw.weights)
            for field in fieldnames(Aswing.Weight)
                @test isapprox(getfield(asw.weights[i], field), getfield(newasw.weights[i], field))
            end
        end
    end
    @testset "sensors" begin
        for i = 1:length(asw.sensors)
            for field in fieldnames(Aswing.Sensor)
                @test isapprox(getfield(asw.sensors[i], field), getfield(newasw.sensors[i], field))
            end
        end
    end
    @testset "engines" begin
        for i = 1:length(asw.engines)
            for field in fieldnames(Aswing.Engine)
                @test isapprox(getfield(asw.engines[i], field), getfield(newasw.engines[i], field))
            end
        end
    end
    @test asw.struts == newasw.struts
    @test asw.joints == newasw.joints
    @test asw.jangles == newasw.jangles
    @test asw.grounds == newasw.grounds
    @testset "beams" begin
        for i = 1:length(asw.beams)
            asw.beams[i].number == newasw.beams[i].number
            asw.beams[i].physical_index == newasw.beams[i].physical_index
            asw.beams[i].name == newasw.beams[i].name
            for field in fieldnames(Aswing.SpanwiseVariables)
                @test isapprox(getfield(asw.beams[i].spanwise_variables, field),
                    getfield(newasw.beams[i].spanwise_variables, field))
            end
        end
    end
end

# test reading/writing operating points
newpnt = Aswing.readpnt(outputpnt)
@testset "pnt" begin
    for i = 1:length(pnt)
        @testset "constraints" begin
            for field in fieldnames(Aswing.Constraints)
                pntvar = getfield(pnt[i].constraints, field)
                newpntvar = getfield(newpnt[i].constraints, field)
                if field == :flap_defl_ctrl_var
                    @test pntvar[1] == newpntvar[1]
                    @test pntvar[2] == newpntvar[2]
                    @test pntvar[3] == newpntvar[3]
                    @test pntvar[4] == newpntvar[4]
                elseif field == :eng_pwr_ctrl_var
                    @test pntvar[1] == newpntvar[1]
                else
                    @test pntvar == newpntvar
                end
            end
        end
        @testset "parameters" begin
            for field in fieldnames(Aswing.Parameters)
                pntvar = getfield(pnt[i].parameters, field)
                newpntvar = getfield(newpnt[i].parameters, field)
                if field == :flap_defl_ctrl_var
                    @test pntvar[1] == newpntvar[1]
                    @test pntvar[2] == newpntvar[2]
                    @test pntvar[3] == newpntvar[3]
                    @test pntvar[4] == newpntvar[4]
                elseif field == :eng_pwr_ctrl_var
                    @test pntvar[1] == newpntvar[1]
                else
                    @test isapprox(pntvar, newpntvar, atol=1e-2)
                end
            end
        end
        @test pnt[i].mach_from_airspeed == newpnt[i].mach_from_airspeed
        @test pnt[i].machpg == newpnt[i].machpg
        @test pnt[i].ground_image == newpnt[i].ground_image
        @test pnt[i].wflap == newpnt[i].wflap
        @test pnt[i].wpeng == newpnt[i].wpeng
    end
end

# test steady solution
steady_solution = Aswing.solvesteady(1)
@test isapprox(steady_solution.pnt_param.lift_force/steady_solution.pnt_param.drag_force, 33.6092341750539, atol=atol)

# get stability derivatives
stabderivs = Aswing.stabilityderivatives()
@testset "stabderivs" begin
    @test isapprox(stabderivs.CXu, -0.03652047310144528, atol=atol)
    @test isapprox(stabderivs.CYu, 0.0, atol=atol)
    @test isapprox(stabderivs.CLu, 1.216904723395558, atol=atol)
    @test isapprox(stabderivs.CXa, 0.514992210591996, atol=atol)
    @test isapprox(stabderivs.CYa, 0.0, atol=atol)
    @test isapprox(stabderivs.CLa, 6.2730761556092975, atol=atol)
    @test isapprox(stabderivs.CXb, 0.0, atol=atol)
    @test isapprox(stabderivs.CYb, -0.42126215334348716, atol=atol)
    @test isapprox(stabderivs.CLb, 0.0, atol=atol)
    @test isapprox(stabderivs.Clu, 0.0, atol=atol)
    @test isapprox(stabderivs.Cmu, 0.04767457583546996, atol=atol)
    @test isapprox(stabderivs.Cnu, 0.0, atol=atol)
    @test isapprox(stabderivs.Cla, 0.0, atol=atol)
    @test isapprox(stabderivs.Cma, -1.166473545760255, atol=atol)
    @test isapprox(stabderivs.Cna, 0.0, atol=atol)
    @test isapprox(stabderivs.Clb, -0.17672244009048493, atol=atol)
    @test isapprox(stabderivs.Cmb, 0.0, atol=atol)
    @test isapprox(stabderivs.Cnb, 0.033382702016544044, atol=atol)
    @test isapprox(stabderivs.CXp, 0.0, atol=atol)
    @test isapprox(stabderivs.CYp, -0.6745224799824928, atol=atol)
    @test isapprox(stabderivs.CLp, 0.0, atol=atol)
    @test isapprox(stabderivs.CXq, -0.45680201241386725, atol=atol)
    @test isapprox(stabderivs.CYq, 0.0, atol=atol)
    @test isapprox(stabderivs.CLq, 7.284291585430452, atol=atol)
    @test isapprox(stabderivs.CXr, 0.0, atol=atol)
    @test isapprox(stabderivs.CYr, 0.25951346743318265, atol=atol)
    @test isapprox(stabderivs.CLr, 0.0, atol=atol)
    @test isapprox(stabderivs.Clp, -0.7078626211816347, atol=atol)
    @test isapprox(stabderivs.Cmp, 0.0, atol=atol)
    @test isapprox(stabderivs.Cnp, -0.06251261778613514, atol=atol)
    @test isapprox(stabderivs.Clq, 0.0, atol=atol)
    @test isapprox(stabderivs.Cmq, -16.26577215749712, atol=atol)
    @test isapprox(stabderivs.Cnq, 0.0, atol=atol)
    @test isapprox(stabderivs.Clr, 0.1623276548641932, atol=atol)
    @test isapprox(stabderivs.Cmr, 0.0, atol=atol)
    @test isapprox(stabderivs.Cnr, -0.03741744012260867, atol=atol)
    @test isapprox(stabderivs.CXdF, [0.0, 5.16115826001186e-6, 0.0, -0.0010600668088926565], atol=atol)
    @test isapprox(stabderivs.CYdF, [0.0065501671250972815, 0.0, 0.0019375002936612249, 0.0], atol=atol)
    @test isapprox(stabderivs.CLdF, [0.0, 0.0037351659479835007, 0.0, 0.0490397612326285], atol=atol)
    @test isapprox(stabderivs.CXdP, [0.0018661778027707714], atol=atol)
    @test isapprox(stabderivs.CYdP, [-5.292539460008798e-9], atol=atol)
    @test isapprox(stabderivs.CLdP, [-0.00010449565952022606], atol=atol)
end

# test setting new paneling
newintervals = [Aswing.getintervals(i) for i  in 1:4]*2
Aswing.nodeset(newintervals)
@test all(newintervals .== [Aswing.getintervals(i) for i in 1:4])

# make sure solution remains similar
steady_solution = Aswing.solvesteady(1)
@test isapprox(steady_solution.pnt_param.lift_force/steady_solution.pnt_param.drag_force, 33.6573231422232, atol=atol)

# get static margin (note that this is using the increased paneling)
sm_qs, sm_rigid, fail = Aswing.getstaticmargin(1)
@test isapprox(sm_qs/100, 0.26230193526999124, atol=atol)
@test isapprox(sm_rigid/100, 0.31478364312678236, atol=atol)

# recover from static margin test (restores default paneling)
Aswing.setgeom(asw)
Aswing.setpnt(pnt)
steady_solution = Aswing.solvesteady(1)

# set up and test unsteady solution
Aswing.retainparam()
Aswing.setcons(Aswing.Constraints("free"))
unsteady_solution = Aswing.solveunsteady(0.1, 10)
@test isapprox(unsteady_solution[11].pnt_param.lift_force/unsteady_solution[11].pnt_param.drag_force, 33.6092341750539, atol=atol)

# test eigenvalues
asweigs = Aswing.geteigs(16, -0.1+0.0im, 1)
@test isapprox(asweigs,
[
-0.02665193416160748 + 0.0im               ,
-0.01689216581284332 + 0.4242720980364713im,
 -0.5286716816061375 + 1.521907733573002im ,
  -7.053040642626703 + 0.0im               ,
 -3.8200059691829984 + 6.2769874595045065im,
  -8.726782350618162 + 8.588264563476661im ,
 -11.749382348992544 + 7.338095974651897im,
 ], atol=atol)


text:
#!/bin/bash
# -*- mode: julia -*-
#=
thisdir="$(dirname "${BASH_SOURCE[0]}")"
exec "$thisdir/julia.sh" --startup-file=no "$@" ${BASH_SOURCE[0]}
=#

pkgid = Base.PkgId(Base.UUID("438e738f-606a-5dbb-bf0a-cddfbfd45ab0"), "PyCall")
sysimage_path = unsafe_string(Base.JLOptions().image_file)
if haskey(Base.loaded_modules, pkgid)
    @info "PyCall is compiled in: `$(sysimage_path)`"
else
    error("PyCall is not compiled in: ", sysimage_path)
end


text:
using JuliaPackageTest
using Test

@testset "JuliaPackageTest.jl" begin
    # Write your tests here.
    x = 2
    y = 3
    @test JuliaPackageTest.sum_values(x,y) == 5
end


text:
export BrickletHallEffectV2MagneticFluxDensityCallbackConfiguration
struct BrickletHallEffectV2MagneticFluxDensityCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
    option::Char
    min::Integer
    max::Integer
end

export BrickletHallEffectV2CounterConfig
struct BrickletHallEffectV2CounterConfig
    high_threshold::Integer
    low_threshold::Integer
    debounce::Integer
end

export BrickletHallEffectV2CounterCallbackConfiguration
struct BrickletHallEffectV2CounterCallbackConfiguration
    period::Integer
    value_has_to_change::Bool
end

export BrickletHallEffectV2SPITFPErrorCount
struct BrickletHallEffectV2SPITFPErrorCount
    error_count_ack_checksum::Integer
    error_count_message_checksum::Integer
    error_count_frame::Integer
    error_count_overflow::Integer
end

export BrickletHallEffectV2Identity
struct BrickletHallEffectV2Identity
    uid::String
    connected_uid::String
    position::Char
    hardware_version::Vector{Integer}
    firmware_version::Vector{Integer}
    device_identifier::Integer
end

export BrickletHallEffectV2
"""
Measures magnetic flux density between -7mT and +7mT
"""
mutable struct BrickletHallEffectV2 <: TinkerforgeDevice
    ipcon::IPConnection
    deviceInternal::PyObject

    """
    Creates an object with the unique device ID *uid* and adds it to
    the IP Connection *ipcon*.
    """
    function BrickletHallEffectV2(uid::String, ipcon::IPConnection)
        package = pyimport("tinkerforge.bricklet_hall_effect_v2")
        deviceInternal = package.BrickletHallEffectV2(uid, ipcon.ipconInternal)

        return new(ipcon, deviceInternal)
    end
end

export get_magnetic_flux_density
"""
    $(SIGNATURES)

Returns the `magnetic flux density (magnetic induction) <https://en.wikipedia.org/wiki/Magnetic_flux>`__.

If you want to get the value periodically, it is recommended to use the
:cb:`Magnetic Flux Density` callback. You can set the callback configuration
with :func:`Set Magnetic Flux Density Callback Configuration`.
"""
function get_magnetic_flux_density(device::BrickletHallEffectV2)
    return device.deviceInternal.get_magnetic_flux_density()
end

export set_magnetic_flux_density_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Magnetic Flux Density` callback is triggered
periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after the value has changed. If the value didn't change
within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.

It is furthermore possible to constrain the callback with thresholds.

The `option`-parameter together with min/max sets a threshold for the :cb:`Magnetic Flux Density` callback.

The following options are possible:

.. csv-table::
 :header: "Option", "Description"
 :widths: 10, 100

 "'x'",    "Threshold is turned off"
 "'o'",    "Threshold is triggered when the value is *outside* the min and max values"
 "'i'",    "Threshold is triggered when the value is *inside* or equal to the min and max values"
 "'<'",    "Threshold is triggered when the value is smaller than the min value (max is ignored)"
 "'>'",    "Threshold is triggered when the value is greater than the min value (max is ignored)"

If the option is set to 'x' (threshold turned off) the callback is triggered with the fixed period.
"""
function set_magnetic_flux_density_callback_configuration(device::BrickletHallEffectV2, period, value_has_to_change, option, min, max)
    device.deviceInternal.set_magnetic_flux_density_callback_configuration(period, value_has_to_change, option, min, max)
end

export get_magnetic_flux_density_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by :func:`Set Magnetic Flux Density Callback Configuration`.
"""
function get_magnetic_flux_density_callback_configuration(device::BrickletHallEffectV2)
    return device.deviceInternal.get_magnetic_flux_density_callback_configuration()
end

export get_counter
"""
    $(SIGNATURES)

Returns the current value of the counter.

You can configure the low/high thresholds and the debounce time
with :func:`Set Counter Config`.

If you set reset counter to *true*, the count is set back to 0
directly after it is read.

If you want to get the count periodically, it is recommended to use the
:cb:`Counter` callback. You can set the callback configuration
with :func:`Set Counter Callback Configuration`.
"""
function get_counter(device::BrickletHallEffectV2, reset_counter)
    return device.deviceInternal.get_counter(reset_counter)
end

export set_counter_config
"""
    $(SIGNATURES)

Sets a high and a low threshold as well as a debounce time.

If the measured magnetic flux density goes above the high threshold or
below the low threshold, the count of the counter is increased by 1.

The debounce time is the minimum time between two count increments.
"""
function set_counter_config(device::BrickletHallEffectV2, high_threshold, low_threshold, debounce)
    device.deviceInternal.set_counter_config(high_threshold, low_threshold, debounce)
end

export get_counter_config
"""
    $(SIGNATURES)

Returns the counter config as set by :func:`Set Counter Config`.
"""
function get_counter_config(device::BrickletHallEffectV2)
    return device.deviceInternal.get_counter_config()
end

export set_counter_callback_configuration
"""
    $(SIGNATURES)

The period is the period with which the :cb:`Counter`
callback is triggered periodically. A value of 0 turns the callback off.

If the `value has to change`-parameter is set to true, the callback is only
triggered after at least one of the values has changed. If the values didn't
change within the period, the callback is triggered immediately on change.

If it is set to false, the callback is continuously triggered with the period,
independent of the value.
"""
function set_counter_callback_configuration(device::BrickletHallEffectV2, period, value_has_to_change)
    device.deviceInternal.set_counter_callback_configuration(period, value_has_to_change)
end

export get_counter_callback_configuration
"""
    $(SIGNATURES)

Returns the callback configuration as set by
:func:`Set Counter Callback Configuration`.
"""
function get_counter_callback_configuration(device::BrickletHallEffectV2)
    return device.deviceInternal.get_counter_callback_configuration()
end

export get_spitfp_error_count
"""
    $(SIGNATURES)

Returns the error count for the communication between Brick and Bricklet.

The errors are divided into

* ACK checksum errors,
* message checksum errors,
* framing errors and
* overflow errors.

The errors counts are for errors that occur on the Bricklet side. All
Bricks have a similar function that returns the errors on the Brick side.
"""
function get_spitfp_error_count(device::BrickletHallEffectV2)
    return device.deviceInternal.get_spitfp_error_count()
end

export set_bootloader_mode
"""
    $(SIGNATURES)

Sets the bootloader mode and returns the status after the requested
mode change was instigated.

You can change from bootloader mode to firmware mode and vice versa. A change
from bootloader mode to firmware mode will only take place if the entry function,
device identifier and CRC are present and correct.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_bootloader_mode(device::BrickletHallEffectV2, mode)
    return device.deviceInternal.set_bootloader_mode(mode)
end

export get_bootloader_mode
"""
    $(SIGNATURES)

Returns the current bootloader mode, see :func:`Set Bootloader Mode`.
"""
function get_bootloader_mode(device::BrickletHallEffectV2)
    return device.deviceInternal.get_bootloader_mode()
end

export set_write_firmware_pointer
"""
    $(SIGNATURES)

Sets the firmware pointer for :func:`Write Firmware`. The pointer has
to be increased by chunks of size 64. The data is written to flash
every 4 chunks (which equals to one page of size 256).

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function set_write_firmware_pointer(device::BrickletHallEffectV2, pointer)
    device.deviceInternal.set_write_firmware_pointer(pointer)
end

export write_firmware
"""
    $(SIGNATURES)

Writes 64 Bytes of firmware at the position as written by
:func:`Set Write Firmware Pointer` before. The firmware is written
to flash every 4 chunks.

You can only write firmware in bootloader mode.

This function is used by Brick Viewer during flashing. It should not be
necessary to call it in a normal user program.
"""
function write_firmware(device::BrickletHallEffectV2, data)
    return device.deviceInternal.write_firmware(data)
end

export set_status_led_config
"""
    $(SIGNATURES)

Sets the status LED configuration. By default the LED shows
communication traffic between Brick and Bricklet, it flickers once
for every 10 received data packets.

You can also turn the LED permanently on/off or show a heartbeat.

If the Bricklet is in bootloader mode, the LED is will show heartbeat by default.
"""
function set_status_led_config(device::BrickletHallEffectV2, config)
    device.deviceInternal.set_status_led_config(config)
end

export get_status_led_config
"""
    $(SIGNATURES)

Returns the configuration as set by :func:`Set Status LED Config`
"""
function get_status_led_config(device::BrickletHallEffectV2)
    return device.deviceInternal.get_status_led_config()
end

export get_chip_temperature
"""
    $(SIGNATURES)

Returns the temperature as measured inside the microcontroller. The
value returned is not the ambient temperature!

The temperature is only proportional to the real temperature and it has bad
accuracy. Practically it is only useful as an indicator for
temperature changes.
"""
function get_chip_temperature(device::BrickletHallEffectV2)
    return device.deviceInternal.get_chip_temperature()
end

export reset
"""
    $(SIGNATURES)

Calling this function will reset the Bricklet. All configurations
will be lost.

After a reset you have to create new device objects,
calling functions on the existing ones will result in
undefined behavior!
"""
function reset(device::BrickletHallEffectV2)
    device.deviceInternal.reset()
end

export write_uid
"""
    $(SIGNATURES)

Writes a new UID into flash. If you want to set a new UID
you have to decode the Base58 encoded UID string into an
integer first.

We recommend that you use Brick Viewer to change the UID.
"""
function write_uid(device::BrickletHallEffectV2, uid)
    device.deviceInternal.write_uid(uid)
end

export read_uid
"""
    $(SIGNATURES)

Returns the current UID as an integer. Encode as
Base58 to get the usual string version.
"""
function read_uid(device::BrickletHallEffectV2)
    return device.deviceInternal.read_uid()
end

export get_identity
"""
    $(SIGNATURES)

Returns the UID, the UID where the Bricklet is connected to,
the position, the hardware and firmware version as well as the
device identifier.

The position can be 'a', 'b', 'c', 'd', 'e', 'f', 'g' or 'h' (Bricklet Port).
A Bricklet connected to an :ref:`Isolator Bricklet <isolator_bricklet>` is always at
position 'z'.

The device identifier numbers can be found :ref:`here <device_identifier>`.
|device_identifier_constant|
"""
function get_identity(device::BrickletHallEffectV2)
    return device.deviceInternal.get_identity()
end

export register_callback
"""
Registers the given *function* with the given *callback_id*.
"""
function register_callback(device::BrickletHallEffectV2, callback_id, function_)
    if isnothing(function_)
        device.registered_callbacks.pop(callback_id, None)
    else
        device.registered_callbacks[callback_id] = function_
    end
end


text:
"""
decurse recursive patterns
"""
struct MemoTreeChildren{P}
    visited::Dict
    child::P
    descend::Bool
end

function MemoTreeChildren(children::Union{Vector, Tuple}, visited::Dict=Dict())
    children_ = []
    for x in children
        push!(children_,
              MemoTreeChildren(visited, x,
                               !haskey(visited, x)))
        visited[x] = true
    end
    children_
end

children(x::CombinedParser, visited::Dict) =
    children(x)

children(x::WrappedAssertion, visited::Dict) =
    children(x.parser, visited)

children(x::MemoTreeChildren) =
    if x.descend
        MemoTreeChildren(children(x.child), x.visited)
    else
        tuple()
    end

function printnode(io::IO, x::MemoTreeChildren)
    printnode(io, x.child)
    x.descend || isempty(children(x.child)) || printstyled(io, " # branches hidden", color=:light_black)
end

Base.show(io::IO, x::MemoTreeChildren) =
    show(io,x.child)

hasregex(x::CombinedParser) = false
hasregex(x::ConstantParser) = true
hasregex(x::WrappedAssertion) = hasregex(x.parser)
hasregex(x::WrappedParser) = false
hasregex(x::Union{AtStart,AtEnd}) = true
function Base.show(io::IO, x::CombinedParser)
    if hasregex(x)
        print(io, "re\"",regex_string(x),"\"") ##!!
    else
        mc = MemoTreeChildren(Dict{Any,Bool}(x=>true),x, true)
        print_tree(io, mc, indicate_truncation=false)
        get(io,:compact,false) || println(io,"::",result_type(x))
    end
end

printnode(io::IO,x::Bytes{N}) where N =
    print(io, "$(N) TypedBytes::$(result_type(x))")

printnode(io::IO, x::CombinedParser) =
    _printnode(io, x)

function _printnode(io::IO, x::CombinedParser)
    printstyled(io, regex_prefix(x), bold=true, color=:cyan)
    if isempty(children(x))
        printstyled(io, regex_inner(x), bold=true, color=:cyan)
        ##printnode(io, x.parser)
    else
        printstyled(io, children_char,bold=true)
    end
    printstyled(io, regex_suffix(x), bold=true, color=:cyan)
    printstyled(io, " ")
    print_constructor(io,x) # , color=:yellow)
end


text:
function ACQR_em_s_run(y,nx::Int;max_iter::Int=100,tol::Float64=1e-6,txo::Bool=true)
	#=
	estimate A, C, Q, R, m1, P1 using the EM algorithm
	and the subspace algorithm for the starting point
	for model

	x_{t+1} = A*x_{t} + w_{t}
	y_{t}   = C*x_{t} + v_{t}

	cov(w_{t},v_{t}) = [Q 0;0 R]
	x1 -> N(m1,P1)

	javier.cara@upm.es, 2107-05

	=#

	# data as a matrix and by rows
	y,ny,nt = byrow(y)

	# starting values from the ssi algorithm
	i = nx+1
	Ai,Ci,Qi,Ri,Si = ACQRS_sub(y,nx,i)
	# x11 = C^{-1}(y_1 - v_1)
	m1i = Ci\y[:,1:1]
	P1i = zeros(nx,nx)

	# em
	A,C,Q,R,m1,P1,loglikv,aic = ACQR_em_s(y,Ai,Ci,Qi,Ri,m1i,max_iter,tol,txo)

	# output
	m = Dict("A"=>A,"C"=>C,"Q"=>Q,"R"=>R,"m1"=>m1,"P1"=>P1,
	"loglikv"=>loglikv[end],"aic"=>aic)
	return m

end


text:
# See Julia Dates
# https://docs.julialang.org/en/v1/stdlib/Dates/#stdlib-dates-api-1

# TODO: remove logging using test instead
using Test
using Logging
logger = ConsoleLogger(stdout, Logging.Info)
global_logger(logger)

using Dates
# https://docs.julialang.org/en/v1/stdlib/Dates/#Dates.DateFormat
# y	        1996, 96	    Returns year of 1996, 0096
# Y	        1996, 96	    Returns year of 1996, 0096. Equivalent to y
# m	        1, 01	        Matches 1 or 2-digit months
# u	        Jan	            Matches abbreviated months according to the locale keyword
# U	        January	        Matches full month names according to the locale keyword
# d	        1, 01	        Matches 1 or 2-digit days
# H	        00	            Matches hours (24-hour clock)
# I	        00	            For outputting hours with 12-hour clock
# M	        00	            Matches minutes
# S	        00	            Matches seconds
# s	        .500	        Matches milliseconds
# e	        Mon, Tues	    Matches abbreviated days of the week
# E	        Monday	        Matches full name days of the week
# p	        AM	            Matches AM/PM (case-insensitive)
# yyyymmdd	19960101	    Matches fixed-width year, month, and day

# https://docs.julialang.org/en/v1/stdlib/Dates/#Constructors
function learn_constructors()
    @test DateTime(2013,1,1) == DateTime(2013)   # partial
    @test DateTime(2013,7,1) == DateTime(2013,7) # partial

    @test_throws ArgumentError DateTime(2013,2,29) # 29 out of range (1:28)

    @test now()              isa DateTime
    @test Date(2020, 12, 23) isa Date
    @test Date(now())        isa Date

    # Construct DateTime from Date
    @test DateTime(2020,1,1) == DateTime(Date(2020,1,1))
end
# https://docs.julialang.org/en/v1/stdlib/Dates/#Durations/Comparisons
function learn_durations_comparisons()
    date_1 = Date(2000,2,1)
    date_2 = Date(2012,2,29)

    @test date_1 <  date_2
    @test date_1 != date_2

    @test_throws MethodError date_1 + date_2
    @test_throws MethodError date_1 * date_2
    @test_throws MethodError date_1 / date_2

    # Duration in Day()
    @test Day( 4411) == date_2 - date_1
    @test Day(-4411) == date_1 - date_2

    dt_1 = DateTime(date_1)
    dt_2 = DateTime(date_2)

    # Duration in Millisecond()
    @test Millisecond(381110400000) == dt_2 - dt_1
end

# https://docs.julialang.org/en/v1/stdlib/Dates/#Accessor-Functions
function learn_accessor()
    date = Date(2014,1,31)

    @test 2014 == year(date)
    @test 1    == month(date)
    @test 5    == week(date)
    @test 31   == day(date)

    @test Year(2014) == Year(date)
    @test 2014       != Year(date)
    @test Day(31)    == Day(date)
    @test 31         != Day(date)

    @test (2014, 1)     == yearmonth(date)
    @test (1, 31)       == monthday(date)
    @test (2014, 1, 31) == yearmonthday(date)
end

# https://docs.julialang.org/en/v1/stdlib/Dates/#Query-Functions
function learn_query()
    date = Date(2014,1,29)

    @test 3           == dayofweek(date) # Wednesday
    @test "Wednesday" == dayname(date)
    @test "Wed"       == dayabbr(date)
    @test "January"   == monthname(date)
    @test "Jan"       == monthabbr(date)
    @test 31          == daysinmonth(date)
    @test false       == isleapyear(date)
    @test 29          == dayofyear(date)
    @test 1           == quarterofyear(date)
    @test 29          == dayofquarter(date)
    # There are a lot more.
end

function learn_time_type_period_arithmetic()
    @test false == isleapyear(2014)

    @test Date(2014,2,28) == Date(2014,1,29) + Month(1)
    @test Date(2014,2,28) == Date(2014,1,30) + Month(1)
    @test Date(2014,2,28) == Date(2014,1,31) + Month(1)

    # Tricky stuff
    @test Date(2014,2,28)  == (Date(2014,1,29) + Day(1)) + Month(1)
    @test Date(2014,3,1)  ==  Date(2014,1,29) + Day(1)  + Month(1)
    @test Date(2014,3,1)  ==  Date(2014,1,29) + Month(1) + Day(1)

    # How to test this?
    ref_date = Date(2014,1,29)
    for date in Date(2014,1,29):Day(1):Date(2014,2,3)
        @test ref_date == date
        ref_date += Day(1)
    end

    # Trick
    ref_date = Date(2014,1,31)
    for date in Date(2014,1,31):Month(1):Date(2014,7,31)
        #@test ref_date == date # does not work
        ref_date += Month(1)
    end
end

# https://docs.julialang.org/en/v1/stdlib/Dates/#Adjuster-Functions
function learn_adjuster()
    @test Date(2014,7,14) ==         firstdayofweek(Date(2014,7,16))
    @test "Monday"        == dayname(firstdayofweek(Date(2014,7,16)))

    @test Date(2014,7,31) == lastdayofmonth(Date(2014,7,16))
    @test Date(2014,9,30) == lastdayofquarter(Date(2014,7,16))

    istuesday = x -> Tuesday == dayofweek(x)

    # 2014-07-13 is a Sunday, the following Tuesday is 2014-07-15
    @test Date(2014,7,15) == tonext(istuesday, Date(2014,7,13))
    # Convenience method to do above
    @test Date(2014,7,15) == tonext(Date(2014,7,13), Tuesday)

    # Find the following Thanksgiving
    @test Date(2014,11,27) == tonext(Date(2014,7,13)) do date
        # Return true when it is the 4th Thursday of November,
        # that is Thanksgiving.
        Thursday == dayofweek(date) &&
        4        == dayofweekofmonth(date) && # 4th Thursday
        November == month(date)
    end

    # filter
    # Pittsburgh street cleaning: Every 2nd Tuesday from April to November
    Ans = filter(Date(2014):Day(1):Date(2015)) do date
        Tuesday == dayofweek(date) &&
        April <= month(date) <= November &&
        2 == dayofweekofmonth(date)
    end
    @test Ans == [Date("2014-04-08"), Date("2014-05-13"), Date("2014-06-10"), Date("2014-07-08"), Date("2014-08-12"), Date("2014-09-09"), Date("2014-10-14"), Date("2014-11-11")] == [Date("2014-04-08"), Date("2014-05-13"), Date("2014-06-10"), Date("2014-07-08"), Date("2014-08-12"), Date("2014-09-09"), Date("2014-10-14"), Date("2014-11-11")]
end

# https://docs.julialang.org/en/v1/stdlib/Dates/#Period-Types
function learn_period_types()
    @test 1 == Year(1).value
    @test Year(3) == Year(1) + Year(2)

    @test 5       == div(Year(10), Year(2))
    @test Year(5) == div(Year(10), 2) # Year(10) ÷ 2
    @test Year(5) == Year(10) / 2
    @test 5       == Year(10) / Year(2)

    @test Year(3) == div(Year(10), 3)
    @test_throws InexactError Year(10) / 3

    @test 10 == Millisecond(10).value
    @test 10 == Dates.value(Millisecond(10))
end

# https://docs.julialang.org/en/v1/stdlib/Dates/#Rounding
function learn_rounding()
    @test Date(1985,8,1) == floor(Date(1985,8,16), Month)
    @test DateTime(2013,2,13,0,45,0) == ceil(DateTime(2013,2,13,0,31,20), Minute(15))
    @test DateTime(2016,8,7) == round(DateTime(2016,8,6,20,15), Day)
    # ISO 8601 standard, 0000-01-01T00:00:00 was chosen as base (or "rounding epoch")
end

function learn_parse()
    @test Date(2021,6,18) == Date("06/18/2021", DateFormat("m/d/y"))
    @test Date(2021,6,18) == Date("06/18/2021", dateformat"m/d/y")
    @test Date(2021,6,18) == Date("6/18/2021",  dateformat"m/d/y")
    @test Date(0021,6,18) == Date("6/18/21",    dateformat"m/d/y")
    @test Date(0021,6,18) == Date("6/18/21",    dateformat"m/d/y")

    @test Date(2021,6,18) == Date("2021-Jun-18", dateformat"y-u-d")

    @test DateTime(2021,6,18,15,38, 17) == DateTime("2021-Jun-18T15:38:17", DateFormat("y-u-dTHH:MM:SS"))
end

function learn_print()
    @test "Friday"      == Dates.format(Date(2021,6,18), "E")
    @test "Saturday"    == Dates.format(Date(2021,6,19), "E")
    @test "Friday-2021" == Dates.format(Date(2021,6,18), "E-yyyy")
    @test "Friday-21"   == Dates.format(Date(2021,6,18), "E-yy")

    @test "2021-Jun-18" == Dates.format(Date(2021,6,18), "yyyy-u-d")
    @test "2021-Jun-18" == Dates.format(Date(2021,6,18), "yyyy-u-dd")
    @test "2021-Jun-1"  == Dates.format(Date(2021,6,1), "yyyy-u-d")
    @test "2021-Jun-01" == Dates.format(Date(2021,6,1), "yyyy-u-dd")

    # Built-in formats
    # ISODateFormat     dateformat"yyyy-mm-dd"
    # ISODateTimeFormat dateformat"yyyy-mm-ddTHH:MM:SS.s"

    # How to compare?
    #@test isequal(ISODateTimeFormat, dateformat"yyyy-mm-ddTHH:MM:SS.s")
end

function learn_arithemtic()
    # +
    @test Date(2021,6,21) == Date(2021,6,18) + Day(3)
    @test_throws MethodError Date(2021,6,18) + Millisecond(3*86400*1000)
    #@test Date(2021,6,21) == Date(2021,6,18) + Millisecond(3*86400*1000)
    @test Date(2021,6,21) == DateTime(2021,6,18) + Millisecond(3*86400*1000)

    # -
    @test Day(11) == Date(2021,10,31) - Date(2021, 10,20)
    @test Day(11) == convert(Day, DateTime(2021,10,31) - DateTime(2021, 10,20)) # milliseconds
    @test Millisecond(11*86400*1000) == DateTime(2021,10,31) - DateTime(2021, 10,20)

    # >
    @test Date(2021,6,21) > Date(2021,6,18)

    # Access value as Int64
    @test 11 == Day(11).value
    @test Day(11).value isa Int64

    # Comparing between Date and DateTime
    @test Date(2021,10,31) == DateTime(2021,10,31)
    @test Date(2021,10,31) > DateTime(2021,10,30, 23)
    @test Date(2021,10,31) < DateTime(2021,10,31, 1)
end

"""
Use to generate a list of dates for use in production schedule.
"""
function print_dates(io::IO)
    for date in Date(2014,1,29):Day(1):Date(2014,2,3)
        str = Dates.format(date, "dd-u-Y")
        println(io, "DATE \"$str\"")
    end
end

# generate a list of date

@test DateTime("2019-2-1T00:00:00") == DateTime(2019, 2, 1)

@test DateTime(2020, 12, 29) > DateTime(2020, 12, 25)
io = stdout
io = devnull
print_dates(io)

@testset "Base" begin
    learn_constructors()
    learn_durations_comparisons()
    learn_accessor()
    learn_query()
    learn_time_type_period_arithmetic()
    learn_adjuster()
    learn_period_types()
    learn_rounding()
    learn_parse()
    learn_print()
    learn_arithemtic()

end
nothing # suppress last line printout


text:
module WarpedRegressors

using AbstractGPs
using Bijectors
using Random

using AbstractGPs: AbstractGP

include("warped_regressor.jl")
include("extras.jl")

export warp

end # module


text:
module ImplicitBnd

using Compat
using Compat.Test
using StaticArrays
using IntervalArithmetic
using EAGOSmoothMcCormickGrad

function XaboutY(x,y,tol)
    return abs(x-y) <= tol
end

EAGOSmoothMcCormickGrad.set_diff_relax(0)

# sets up initial options
opts1 = mc_opts(Float64)          # sets options for relaxtion
opts1.kmax = 3             # sets number of iterations
opts1.style = "KrawczykCW"   # sets style of contractor

#=
generates the expansion point parameters for the function using the opts
options using inverse preconditioner
=#
f(x,p) = x[1]*p[1]+p[1]
g(x,p) = [x[1]*p[1]+p[1];
          x[1]*p[1]+2*p[1]]
function h1(x,p)
    t1 = x[1]^2
    t2 = x[1]*p[1]
    t3 = 4.0
    t4 = t1 + t2
    t5 = t2 + t3
    return [t5]
end
hj1(x,p) = [2*x[1]+p[1]]
P = [Interval(6.0,9.0)]
X = [Interval(-0.78,-0.4)]
p = [7.5]
pmid = mid.(P)

#=
relaxes the equality h(x,p)
=#
np = 1
szero = @SVector zeros(np)
sone = @SVector ones(np)
p_mc = [SMCg{np,Float64}(p[i],p[i],sone,sone,@interval(P[i].lo,P[i].hi),false,[∅],[1.0]) for i=1:np]

#=
generates the expansion point parameters for the function using the opts
options using in place LDU full pivot conditioner with sparse calc
=#
f1(x,p) = x[1]*p[1]+p[1]
g1(x,p) = [x[1]*p[1]+p[1];
          x[1]*p[1]+2*p[1]]
function h1!(hout,x,p)
    t1 = x[1]^2
    t2 = x[1]*p[1]
    t3 = 4.0
    t4 = t1 + t2
    t5 = t2 + t3
    hout[1] = t5
end
function hj1!(hout,x,p)
    hout[1] = 2*x[1]+p[1]
end
P1 = [Interval(6.0,9.0)]
X1 = [Interval(-0.78,-0.4)]
p1 = [7.5]
pmid1 = mid.(P1)

np1 = 1
szero1 = @SVector zeros(np1)
sone1 = @SVector ones(np1)
p_mc1 = [SMCg{np,Float64}(p1[i],p1[i],sone1,sone1,@interval(P1[i].lo,P1[i].hi),false,[∅],[1.0]) for i=1:np1]

param = GenExpansionParams(h1,hj1,X,P,pmid,opts1)
@test XaboutY(param[2][1].cc,-0.46623417721518967,1E-8)
@test XaboutY(param[2][1].cv,-0.5792721518987343,1E-8)
hbnds = MC_impRelax(h1,hj1,p_mc,pmid,X,P,opts1,param)
@test XaboutY(hbnds[1].cc,-0.4907369110006179,1E-8)
@test XaboutY(hbnds[1].cv,-0.5758134311908059,1E-8)
fbnds = impRelax_f(f,h1,hj1,X,P,p,pmid,opts1,param)
@test XaboutY(fbnds.cc,3.955578533996293,1E-8)
@test XaboutY(fbnds.cv,2.9947278349146904,1E-8)
fgbnds = impRelax_fg(f,g,h1,hj1,X,P,p,pmid,opts1,param)
@test XaboutY(fgbnds[1].cc,3.955578533996293,1E-8)
@test XaboutY(fgbnds[1].cv,2.9947278349146904,1E-8)
@test XaboutY(fgbnds[2][1].cc,3.955578533996293,1E-8)
@test XaboutY(fgbnds[2][1].cv,2.9947278349146904,1E-8)
@test XaboutY(fgbnds[2][2].cc,11.455578533996293,1E-8)
@test XaboutY(fgbnds[2][2].cv,10.494727834914691,1E-8)

param1 = InGenExpansionParams(h1!,hj1!,X1,P1,pmid1,opts1)
@test XaboutY(param1[2][1].cc,-0.46623417721518967,1E-8)
@test XaboutY(param1[2][1].cv,-0.5792721518987343,1E-8)
# relaxation of f(x,p) & g(x,p) at (x(p),p) SPARSE SOLVER
hbnds1 = MC_NimpRelax(h1!,hj1!,p_mc1,pmid1,X1,P1,opts1,param1)
@test XaboutY(hbnds1[1].cc,-0.4907369110006179,1E-8)
@test XaboutY(hbnds1[1].cv,-0.5758134311908059,1E-8)
fbnds1 = NimpRelax_f(f1,h1!,hj1!,X1,P1,p,pmid1,opts1,param1)
@test XaboutY(fbnds1.cc,3.955578533996293,1E-8)
@test XaboutY(fbnds1.cv,2.9947278349146904,1E-8)
fgbnds1 = NimpRelax_fg(f1,g1,h1!,hj1!,X1,P1,p,pmid1,opts1,param1)
@test XaboutY(fgbnds1[1].cc,3.955578533996293,1E-8)
@test XaboutY(fgbnds1[1].cv,2.9947278349146904,1E-8)
@test XaboutY(fgbnds1[2][1].cc,3.955578533996293,1E-8)
@test XaboutY(fgbnds1[2][1].cv,2.9947278349146904,1E-8)
@test XaboutY(fgbnds1[2][2].cc,11.455578533996293,1E-8)
@test XaboutY(fgbnds1[2][2].cv,10.494727834914691,1E-8)

# relaxation of f(x,p) & g(x,p) at (x(p),p) DIRECT SOLVER
param2 = IndGenExpansionParams(h1!,hj1!,X1,P1,pmid1,opts1)
@test XaboutY(param2[2][1].cc,-0.46623417721518967,1E-8)
@test XaboutY(param2[2][1].cv,-0.5792721518987343,1E-8)
hbnds2 = MC_NdimpRelax(h1!,hj1!,p_mc1,pmid1,X1,P1,opts1,param2)
@test XaboutY(hbnds2[1].cc,-0.4907369110006179,1E-8)
@test XaboutY(hbnds2[1].cv,-0.5758134311908059,1E-8)
fbnds2 = NdimpRelax_f(f1,h1!,hj1!,X1,P1,p,pmid1,opts1,param2)
@test XaboutY(fbnds2.cc,3.955578533996293,1E-8)
@test XaboutY(fbnds2.cv,2.9947278349146904,1E-8)
fgbnds2 = NdimpRelax_fg(f1,g1,h1!,hj1!,X1,P1,p,pmid1,opts1,param2)
@test XaboutY(fgbnds2[1].cc,3.955578533996293,1E-8)
@test XaboutY(fgbnds2[1].cv,2.9947278349146904,1E-8)
@test XaboutY(fgbnds2[2][1].cc,3.955578533996293,1E-8)
@test XaboutY(fgbnds2[2][1].cv,2.9947278349146904,1E-8)
@test XaboutY(fgbnds2[2][2].cc,11.455578533996293,1E-8)
@test XaboutY(fgbnds2[2][2].cv,10.494727834914691,1E-8)

opts2 = mc_opts(Float64)   # sets options for relaxtion
opts2.kmax = 3             # sets number of iterations
opts2.style = "NewtonGS"   # sets style of contractor

param3 = GenExpansionParams(h1,hj1,X,P,pmid,opts2)
@test XaboutY(param3[2][1].cc,-0.4593323920935859,1E-8)
@test XaboutY(param3[2][1].cv,-0.5667248956273347,1E-8)
hbnds3 = MC_impRelax(h1,hj1,p_mc,pmid,X,P,opts2,param)
@test XaboutY(hbnds3[1].cc,-0.5014435013618904,1E-8)
@test XaboutY(hbnds3[1].cv,-0.5667248956273347,1E-8)
fbnds3 = impRelax_f(f,h1,hj1,X,P,p,pmid,opts2,param)
@test XaboutY(fbnds3.cc,3.891338991828657,1E-8)
@test XaboutY(fbnds3.cv,3.0592452208305847,1E-8)
fgbnds3 = impRelax_fg(f,g,h1,hj1,X,P,p,pmid,opts2,param)
@test XaboutY(fgbnds3[1].cc,3.891338991828657,1E-8)
@test XaboutY(fgbnds3[1].cv,3.0592452208305847,1E-8)
@test XaboutY(fgbnds3[2][1].cc,3.891338991828657,1E-8)
@test XaboutY(fgbnds3[2][1].cv,3.0592452208305847,1E-8)
@test XaboutY(fgbnds3[2][2].cc,11.391338991828658,1E-8)
@test XaboutY(fgbnds3[2][2].cv,10.559245220830585,1E-8)

param4 = InGenExpansionParams(h1!,hj1!,X1,P1,pmid1,opts2)
@test XaboutY(param4[2][1].cc,-0.4593323920935859,1E-8)
@test XaboutY(param4[2][1].cv,-0.5667248956273347,1E-8)
hbnds4 = MC_NimpRelax(h1!,hj1!,p_mc1,pmid1,X1,P1,opts2,param1)
@test XaboutY(hbnds4[1].cc,-0.5014435013618905,1E-8)
@test XaboutY(hbnds4[1].cv,-0.5667248956273346,1E-8)
fbnds4 = NimpRelax_f(f1,h1!,hj1!,X1,P1,p,pmid1,opts2,param1)
@test XaboutY(fbnds4.cc,3.891338991828657,1E-8)
@test XaboutY(fbnds4.cv,3.0592452208305865,1E-8)
fgbnds4 = NimpRelax_fg(f1,g1,h1!,hj1!,X1,P1,p,pmid1,opts2,param1)
@test XaboutY(fgbnds4[1].cc,3.891338991828657,1E-8)
@test XaboutY(fgbnds4[1].cv,3.0592452208305865,1E-8)
@test XaboutY(fgbnds4[2][1].cc,3.891338991828657,1E-8)
@test XaboutY(fgbnds4[2][1].cv,3.0592452208305865,1E-8)
@test XaboutY(fgbnds4[2][2].cc,11.391338991828658,1E-8)
@test XaboutY(fgbnds4[2][2].cv,10.559245220830586,1E-8)

# relaxation of f(x,p) & g(x,p) at (x(p),p) DIRECT SOLVER
param5 = IndGenExpansionParams(h1!,hj1!,X1,P1,pmid1,opts2)
@test XaboutY(param5[2][1].cc,-0.4593323920935859,1E-8)
@test XaboutY(param5[2][1].cv,-0.5667248956273347,1E-8)
hbnds5 = MC_NdimpRelax(h1!,hj1!,p_mc1,pmid1,X1,P1,opts2,param2)
@test XaboutY(hbnds5[1].cc,-0.5014435013618904,1E-8)
@test XaboutY(hbnds5[1].cv,-0.5667248956273347,1E-8)
fbnds5 = NdimpRelax_f(f1,h1!,hj1!,X1,P1,p,pmid1,opts2,param2)
@test XaboutY(fbnds5.cc,3.891338991828657,1E-8)
@test XaboutY(fbnds5.cv,3.0592452208305847,1E-8)
fgbnds5 = NdimpRelax_fg(f1,g1,h1!,hj1!,X1,P1,p,pmid1,opts2,param2)
@test XaboutY(fgbnds5[1].cc,3.891338991828657,1E-8)
@test XaboutY(fgbnds5[1].cv,3.0592452208305847,1E-8)
@test XaboutY(fgbnds5[2][1].cc,3.891338991828657,1E-8)
@test XaboutY(fgbnds5[2][1].cv,3.0592452208305847,1E-8)
@test XaboutY(fgbnds5[2][2].cc,11.391338991828658,1E-8)
@test XaboutY(fgbnds5[2][2].cv,10.559245220830585,1E-8)

end


text:
module SPRING
using Plots
    toto(x) = cos(x)

    struct MyStruct
        myfield::Float64
    end

    function go()
        x=[Float64(i) for i ∈ 1:40]
        y = toto.(x)
        plot(x,y)

        mf=MyStruct(3)
        @show mf

    end
end # module


text:
Pkg.add("Images")
Pkg.add("DataFrames")
using Images
using DataFrames

images = imread(nameFile)
temp = float32sc(img)
if ndims(temp) == 3
  temp = mean(temp.data, 1)
end
x[i, :] = reshape(temp, 1 imageSize)

function read_data(typeData, labelsInfo, imageSize, path)
  x = zeros(size(labelsInfo, 1), imageSize)
  for(index, idImage) in enumerate(labelsInfo["ID"])
    nameFile = "$(path)/$(typeData)Resized/$(idImage).Bmp
    img = imread(nameFile)

    temp = float32sc(img)

    if ndims(temp) == 3
      temp = mean(temp.data, 1)
    end

    x[index, :] = reshape(temp, 1, imageSize)
  end
  return x
end

max_stars_repo_name: danielbrownjr/super-disco

max_issues_repo_name: danielbrownjr/super-disco

max_forks_repo_name: danielbrownjr/super-disco


text:
# include("../src/CALiPPSO.jl")
# using .CALiPPSO
using CALiPPSO
using Documenter

makedocs(sitename="CALiPPSO",
    authors = "Rafael Diaz Hernandez Rojas",
    pages = [
        "Home" => "index.md",
        "Introduction: theory and few details" => "theory.md",
        "User Guide" => Any["Installation" => "installation.md", "Basic usage" =>"basic_usage.md", "How the main function works" => "mainfunction.md", "Changing the default behaviour" =>"changing_default.md"],
        "Examples" => "tests.md",
        "`Structs` of CALiPPSO" => "types.md",
        "Problem solving and possible issues" => "issues.md",
        "To Do's"=> "todos.md",
        "API Reference" => "api.md"],
        format = Documenter.HTML(prettyurls = false)
    )

deploydocs(
    repo = "github.com/rdhr/CALiPPSO.jl.git",
)

max_stars_repo_name: rdhr/CALiPPSO

max_issues_repo_name: rdhr/CALiPPSO

max_forks_repo_name: rdhr/CALiPPSO


text:
# RELOG: Reverse Logistics Optimization
# Copyright (C) 2020-2021, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using JuMP

function resolve(model_old, filename::AbstractString; kwargs...)::OrderedDict
    @info "Reading $filename..."
    instance = RELOG.parsefile(filename)
    return resolve(model_old, instance; kwargs...)
end

function resolve(model_old, instance::Instance; optimizer = nothing)::OrderedDict
    milp_optimizer = lp_optimizer = optimizer
    if optimizer === nothing
        milp_optimizer = _get_default_milp_optimizer()
        lp_optimizer = _get_default_lp_optimizer()
    end

    @info "Building new graph..."
    graph = build_graph(instance)
    _print_graph_stats(instance, graph)

    @info "Building new optimization model..."
    model_new = RELOG.build_model(instance, graph, milp_optimizer)

    @info "Fixing decision variables..."
    _fix_plants!(model_old, model_new)
    JuMP.set_optimizer(model_new, lp_optimizer)

    @info "Optimizing MILP..."
    JuMP.optimize!(model_new)

    if !has_values(model_new)
        @warn("No solution available")
        return OrderedDict()
    end

    @info "Extracting solution..."
    solution = get_solution(model_new, marginal_costs = true)

    return solution
end

function _fix_plants!(model_old, model_new)::Nothing
    T = model_new[:instance].time

    # Fix open_plant variables
    for ((node_old, t), var_old) in model_old[:open_plant]
        value_old = JuMP.value(var_old)
        node_new = model_new[:graph].name_to_process_node_map[(
            node_old.location.plant_name,
            node_old.location.location_name,
        )]
        var_new = model_new[:open_plant][node_new, t]
        JuMP.unset_binary(var_new)
        JuMP.fix(var_new, value_old)
    end

    # Fix is_open variables
    for ((node_old, t), var_old) in model_old[:is_open]
        value_old = JuMP.value(var_old)
        node_new = model_new[:graph].name_to_process_node_map[(
            node_old.location.plant_name,
            node_old.location.location_name,
        )]
        var_new = model_new[:is_open][node_new, t]
        JuMP.unset_binary(var_new)
        JuMP.fix(var_new, value_old)
    end

    # Fix plant capacities
    for ((node_old, t), var_old) in model_old[:capacity]
        value_old = JuMP.value(var_old)
        node_new = model_new[:graph].name_to_process_node_map[(
            node_old.location.plant_name,
            node_old.location.location_name,
        )]
        var_new = model_new[:capacity][node_new, t]
        JuMP.delete_lower_bound(var_new)
        JuMP.delete_upper_bound(var_new)
        JuMP.fix(var_new, value_old)
    end

    # Fix plant expansion
    for ((node_old, t), var_old) in model_old[:expansion]
        value_old = JuMP.value(var_old)
        node_new = model_new[:graph].name_to_process_node_map[(
            node_old.location.plant_name,
            node_old.location.location_name,
        )]
        var_new = model_new[:expansion][node_new, t]
        JuMP.delete_lower_bound(var_new)
        JuMP.delete_upper_bound(var_new)
        JuMP.fix(var_new, value_old)
    end
end

max_stars_repo_name: ANL-CEEESA/RELOG

max_issues_repo_name: ANL-CEEESA/RELOG

max_forks_repo_name: ANL-CEEESA/RELOG


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

## semantic version numbers (https://semver.org/)

const VInt = UInt32
"""
    VersionNumber

Version number type which follow the specifications of
[semantic versioning](https://semver.org/), composed of major, minor
and patch numeric values, followed by pre-release and build
alpha-numeric annotations. See also [`@v_str`](@ref).

# Examples
```jldoctest
julia> VersionNumber("1.2.3")
v"1.2.3"

julia> VersionNumber("2.0.1-rc1")
v"2.0.1-rc1"
```
"""
struct VersionNumber
    major::VInt
    minor::VInt
    patch::VInt
    prerelease::Tuple{Vararg{Union{UInt64,String}}}
    build::Tuple{Vararg{Union{UInt64,String}}}

    function VersionNumber(major::VInt, minor::VInt, patch::VInt,
            pre::Tuple{Vararg{Union{UInt64,String}}},
            bld::Tuple{Vararg{Union{UInt64,String}}})
        major >= 0 || throw(ArgumentError("invalid negative major version: $major"))
        minor >= 0 || throw(ArgumentError("invalid negative minor version: $minor"))
        patch >= 0 || throw(ArgumentError("invalid negative patch version: $patch"))
        for ident in pre
            if ident isa Integer
                ident >= 0 || throw(ArgumentError("invalid negative pre-release identifier: $ident"))
            else
                if !occursin(r"^(?:|[0-9a-z-]*[a-z-][0-9a-z-]*)$"i, ident) ||
                    isempty(ident) && !(length(pre)==1 && isempty(bld))
                    throw(ArgumentError("invalid pre-release identifier: $(repr(ident))"))
                end
            end
        end
        for ident in bld
            if ident isa Integer
                ident >= 0 || throw(ArgumentError("invalid negative build identifier: $ident"))
            else
                if !occursin(r"^(?:|[0-9a-z-]*[a-z-][0-9a-z-]*)$"i, ident) ||
                    isempty(ident) && length(bld)!=1
                    throw(ArgumentError("invalid build identifier: $(repr(ident))"))
                end
            end
        end
        new(major, minor, patch, pre, bld)
    end
end
VersionNumber(major::Integer, minor::Integer = 0, patch::Integer = 0,
        pre::Tuple{Vararg{Union{Integer,AbstractString}}} = (),
        bld::Tuple{Vararg{Union{Integer,AbstractString}}} = ()) =
    VersionNumber(VInt(major), VInt(minor), VInt(patch),
        map(x->x isa Integer ? UInt64(x) : String(x), pre),
        map(x->x isa Integer ? UInt64(x) : String(x), bld))

VersionNumber(v::Tuple) = VersionNumber(v...)

function print(io::IO, v::VersionNumber)
    v == typemax(VersionNumber) && return print(io, "∞")
    print(io, v.major)
    print(io, '.')
    print(io, v.minor)
    print(io, '.')
    print(io, v.patch)
    if !isempty(v.prerelease)
        print(io, '-')
        join(io, v.prerelease,'.')
    end
    if !isempty(v.build)
        print(io, '+')
        join(io, v.build,'.')
    end
end
show(io::IO, v::VersionNumber) = print(io, "v\"", v, "\"")

Broadcast.broadcastable(v::VersionNumber) = Ref(v)

const VERSION_REGEX = r"^
    v?                                      # prefix        (optional)
    (\d+)                                   # major         (required)
    (?:\.(\d+))?                            # minor         (optional)
    (?:\.(\d+))?                            # patch         (optional)
    (?:(-)|                                 # pre-release   (optional)
    ([a-z][0-9a-z-]*(?:\.[0-9a-z-]+)*|-(?:[0-9a-z-]+\.)*[0-9a-z-]+)?
    (?:(\+)|
    (?:\+((?:[0-9a-z-]+\.)*[0-9a-z-]+))?    # build         (optional)
    ))
$"ix

function split_idents(s::AbstractString)
    idents = split(s, '.')
    ntuple(length(idents)) do i
        ident = idents[i]
        occursin(r"^\d+$", ident) ? parse(UInt64, ident) : String(ident)
    end
end

function VersionNumber(v::AbstractString)
    v == "∞" && return typemax(VersionNumber)
    m = match(VERSION_REGEX, v)
    m === nothing && throw(ArgumentError("invalid version string: $v"))
    major, minor, patch, minus, prerl, plus, build = m.captures
    major = parse(VInt, major)
    minor = minor !== nothing ? parse(VInt, minor) : VInt(0)
    patch = patch !== nothing ? parse(VInt, patch) : VInt(0)
    if prerl !== nothing && !isempty(prerl) && prerl[1] == '-'
        prerl = prerl[2:end] # strip leading '-'
    end
    prerl = prerl !== nothing ? split_idents(prerl) : minus !== nothing ? ("",) : ()
    build = build !== nothing ? split_idents(build) : plus  !== nothing ? ("",) : ()
    return VersionNumber(major, minor, patch, prerl, build)
end

"""
    @v_str

String macro used to parse a string to a [`VersionNumber`](@ref).

# Examples
```jldoctest
julia> v"1.2.3"
v"1.2.3"

julia> v"2.0.1-rc1"
v"2.0.1-rc1"
```
"""
macro v_str(v); VersionNumber(v); end

typemin(::Type{VersionNumber}) = v"0-"

function typemax(::Type{VersionNumber})
    ∞ = typemax(VInt)
    VersionNumber(∞, ∞, ∞, (), ("",))
end

ident_cmp(a::Integer, b::Integer) = cmp(a, b)
ident_cmp(a::Integer, b::String ) = isempty(b) ? +1 : -1
ident_cmp(a::String,  b::Integer) = isempty(a) ? -1 : +1
ident_cmp(a::String,  b::String ) = cmp(a, b)

function ident_cmp(
    A::Tuple{Vararg{Union{Integer,String}}},
    B::Tuple{Vararg{Union{Integer,String}}},
)
    for (a, b) in zip(A, B)
       c = ident_cmp(a,b)
       (c != 0) && return c
    end
    length(A) < length(B) ? -1 :
    length(B) < length(A) ? +1 : 0
end

function ==(a::VersionNumber, b::VersionNumber)
    (a.major != b.major) && return false
    (a.minor != b.minor) && return false
    (a.patch != b.patch) && return false
    (ident_cmp(a.prerelease, b.prerelease) != 0) && return false
    (ident_cmp(a.build, b.build) != 0) && return false
    return true
end

issupbuild(v::VersionNumber) = length(v.build)==1 && isempty(v.build[1])

function isless(a::VersionNumber, b::VersionNumber)
    (a.major < b.major) && return true
    (a.major > b.major) && return false
    (a.minor < b.minor) && return true
    (a.minor > b.minor) && return false
    (a.patch < b.patch) && return true
    (a.patch > b.patch) && return false
    (!isempty(a.prerelease) && isempty(b.prerelease)) && return true
    (isempty(a.prerelease) && !isempty(b.prerelease)) && return false
    c = ident_cmp(a.prerelease,b.prerelease)
    (c < 0) && return true
    (c > 0) && return false
    (!issupbuild(a) && issupbuild(b)) && return true
    (issupbuild(a) && !issupbuild(b)) && return false
    c = ident_cmp(a.build,b.build)
    (c < 0) && return true
    return false
end

function hash(v::VersionNumber, h::UInt)
    h += 0x8ff4ffdb75f9fede % UInt
    h = hash(v.major, h)
    h = hash(v.minor, h)
    h = hash(v.patch, h)
    h = hash(v.prerelease, ~h)
    h = hash(v.build, ~h)
end

lowerbound(v::VersionNumber) = VersionNumber(v.major, v.minor, v.patch, ("",), ())
upperbound(v::VersionNumber) = VersionNumber(v.major, v.minor, v.patch, (), ("",))

thispatch(v::VersionNumber) = VersionNumber(v.major, v.minor, v.patch)
thisminor(v::VersionNumber) = VersionNumber(v.major, v.minor, 0)
thismajor(v::VersionNumber) = VersionNumber(v.major, 0, 0)

nextpatch(v::VersionNumber) = v < thispatch(v) ? thispatch(v) : VersionNumber(v.major, v.minor, v.patch+1)
nextminor(v::VersionNumber) = v < thisminor(v) ? thisminor(v) : VersionNumber(v.major, v.minor+1, 0)
nextmajor(v::VersionNumber) = v < thismajor(v) ? thismajor(v) : VersionNumber(v.major+1, 0, 0)

function check_new_version(existing::Vector{VersionNumber}, ver::VersionNumber)
    if isempty(existing)
        for v in [v"0.0.1", v"0.1", v"1"]
            lowerbound(v) <= ver <= v && return
        end
        error("version $ver is invalid initial version (try 0.0.1, 0.1, 1.0)")
    end
    issorted(existing) || (existing = sort(existing))
    idx = searchsortedlast(existing, ver)
    idx > 0 || error("version $ver less than least existing version $(existing[1])")
    prv = existing[idx]
    ver == prv && error("version $ver already exists")
    nxt = thismajor(ver) != thismajor(prv) ? nextmajor(prv) :
          thisminor(ver) != thisminor(prv) ? nextminor(prv) : nextpatch(prv)
    ver <= nxt || error("version $ver skips over $nxt")
    thispatch(ver) <= ver && return # regular or build release
    idx < length(existing) && thispatch(existing[idx+1]) <= nxt &&
        error("version $ver is pre-release of existing version $(existing[idx+1])")
    return # acceptable new version
end

## julia version info

"""
    VERSION

A `VersionNumber` object describing which version of Julia is in use. For details see
[Version Number Literals](@ref man-version-number-literals).
"""
const VERSION = try
    ver = VersionNumber(VERSION_STRING)
    if !isempty(ver.prerelease)
        if GIT_VERSION_INFO.build_number >= 0
            ver = VersionNumber(ver.major, ver.minor, ver.patch, (ver.prerelease..., GIT_VERSION_INFO.build_number), ver.build)
        else
            println("WARNING: no build number found for pre-release version")
            ver = VersionNumber(ver.major, ver.minor, ver.patch, (ver.prerelease..., "unknown"), ver.build)
        end
    elseif GIT_VERSION_INFO.build_number > 0
        println("WARNING: ignoring non-zero build number for VERSION")
    end
    ver
catch e
    println("while creating Base.VERSION, ignoring error $e")
    VersionNumber(0)
end

const libllvm_version = VersionNumber("6.0.1")

function banner(io::IO = stdout)
    if GIT_VERSION_INFO.tagged_commit
        commit_string = TAGGED_RELEASE_BANNER
    elseif isempty(GIT_VERSION_INFO.commit)
        commit_string = ""
    else
        days = Int(floor((ccall(:jl_clock_now, Float64, ()) - GIT_VERSION_INFO.fork_master_timestamp) / (60 * 60 * 24)))
        days = max(0, days)
        unit = days == 1 ? "day" : "days"
        distance = GIT_VERSION_INFO.fork_master_distance
        commit = GIT_VERSION_INFO.commit_short

        if distance == 0
            commit_string = "Commit $(commit) ($(days) $(unit) old master)"
        else
            branch = GIT_VERSION_INFO.branch
            commit_string = "$(branch)/$(commit) (fork: $(distance) commits, $(days) $(unit))"
        end
    end

    commit_date = isempty(Base.GIT_VERSION_INFO.date_string) ? "" : " ($(split(Base.GIT_VERSION_INFO.date_string)[1]))"

    if get(io, :color, false)
        c = text_colors
        tx = c[:normal] # text
        jl = c[:normal] # julia
        d1 = c[:bold] * c[:blue]    # first dot
        d2 = c[:bold] * c[:red]     # second dot
        d3 = c[:bold] * c[:green]   # third dot
        d4 = c[:bold] * c[:magenta] # fourth dot

        print(io,"""               $(d3)_$(tx)
           $(d1)_$(tx)       $(jl)_$(tx) $(d2)_$(d3)(_)$(d4)_$(tx)     |  Documentation: https://docs.julialang.org
          $(d1)(_)$(jl)     | $(d2)(_)$(tx) $(d4)(_)$(tx)    |
           $(jl)_ _   _| |_  __ _$(tx)   |  Type \"?\" for help, \"]?\" for Pkg help.
          $(jl)| | | | | | |/ _` |$(tx)  |
          $(jl)| | |_| | | | (_| |$(tx)  |  Version $(VERSION)$(commit_date)
         $(jl)_/ |\\__'_|_|_|\\__'_|$(tx)  |  $(commit_string)
        $(jl)|__/$(tx)                   |

        """)
    else
        print(io,"""
                       _
           _       _ _(_)_     |  Documentation: https://docs.julialang.org
          (_)     | (_) (_)    |
           _ _   _| |_  __ _   |  Type \"?\" for help, \"]?\" for Pkg help.
          | | | | | | |/ _` |  |
          | | |_| | | | (_| |  |  Version $(VERSION)$(commit_date)
         _/ |\\__'_|_|_|\\__'_|  |  $(commit_string)
        |__/                   |

        """)
    end
end

max_stars_repo_name: mhaldrich/julia

max_issues_repo_name: mhaldrich/julia

max_forks_repo_name: mhaldrich/julia


text:
using ProjAssistant
@quickactivate "Chemostat_Rath2017"

@time begin

    using ProgressMeter
    using Base.Threads

    import Chemostat
    import Chemostat.MetNets
    import Chemostat.MetLP
    import Chemostat.MetEP
    const Ch = Chemostat

    import Chemostat_Rath2017
    const ChR = Chemostat_Rath2017
    const Rd = ChR.RathData
    const H1 = ChR.Human1
    const HG = H1.HumanGEM
    const AG = H1.HartGEMs

    using Plots

    using Statistics
end

## ---------------------------------------------------------------------
# DESCRIPTION
# This converge MaxEnt assuming the culture is glucose limited

## ---------------------------------------------------------------------
let
    mets_map = HG.load_mets_map()
    for rath_met in ["GLC", "GLN", "GAL"]
        model_met = mets_map[rath_met]
        @show rath_met model_met
    end
end


text:
@testset "EDM" begin
    n = 17
    r = 3
    ijs = [(3,1), (3,2), (4,1), (4,2), (4,3), (5,1), (5,3), (5,4), (8,7),
        (9,5), (9,6), (9,7), (9,8), (11,10), (13,10), (13,11), (13,12),
        (14,10), (14,11), (14,12), (14,13),
        (15, 1), (15,3), (15,4), (15,5), (15,7), (15,8), (15,9),
        (16,5), (16,6), (16,9), (16,12), (16,13), (16,14), (16,15),
        (17,10), (17,11), (17,12), (17,13), (17,14), (17,15), (17,16)]
    append!(ijs, [(i,i) for i in 1:n])
    II, JJ = Chordal.unzip(ijs)
    sp = sparse(II, JJ, ones(length(II)))
    sp = sp + tril(sp)'

    V = randn(n, r)
    VV = V*V'
    Afull = diag(VV)*ones(n)' + ones(n)*diag(VV)' - 2VV
    A = sp .* Afull

    W = edm_completion(A)
    @test Chordal.is_edm(W)
end


text:
using AlgebraicRelations.DB
using AlgebraicRelations.Queries
using SQLite
using DataFrames
using Catlab

@present WorkplaceSchema <: TheorySQL begin
  # Data tables
  employees::Ob
  employees_1_person::Attr(employees, Int)
  employees_2_id::Attr(employees, Int)

  names::Ob
  names_1_person::Hom(names, employees)
  names_2_full_name::Attr(names, String)

  salary::Ob
  salary_1_person::Hom(salary, employees)
  salary_2_salary::Attr(salary, Real)

  # Relation tables
  manager::Ob
  manager_1_person::Hom(manager, employees)
  manager_2_manager::Hom(manager, employees)

  relation::Ob
  relation_1_person1::Hom(relation, employees)
  relation_2_person2::Hom(relation, employees)
  relation_3_relationship::Attr(relation, Real)
end;

@db_schema Workplace(WorkplaceSchema)
schema = Workplace()

db = SQLite.DB()
splt_stmts = split(generate_schema_sql(schema), "\n")

@testset "Generate DB Schema" begin
  for stmt in splt_stmts
    @test DBInterface.execute(db, stmt) isa SQLite.Query
  end
end

# Fill out the table
insert_stmts = ["INSERT INTO employees    VALUES (1, 1);",
                "INSERT INTO employees    VALUES (2, 2);",
                "INSERT INTO employees    VALUES (3, 3);",
                "INSERT INTO employees    VALUES (4, 4);",
                "INSERT INTO names        VALUES (1, 'Alice Smith');",
                "INSERT INTO names        VALUES (2, 'Bob Jones');",
                "INSERT INTO names        VALUES (3, 'Eve Johnson');",
                "INSERT INTO names        VALUES (4, 'John Doe');",
                "INSERT INTO manager      VALUES (1, 1);",
                "INSERT INTO manager      VALUES (2, 1);",
                "INSERT INTO manager      VALUES (3, 4);",
                "INSERT INTO manager      VALUES (4, 1);",
                "INSERT INTO salary       VALUES (1, 150000);",
                "INSERT INTO salary       VALUES (2, 50000);",
                "INSERT INTO salary       VALUES (3, 80000);",
                "INSERT INTO salary       VALUES (4, 90000);"]

for stmt in insert_stmts
  DBInterface.execute(db, stmt)
end

@testset "Generate SQL Queries" begin
  q = @query schema (p, n) where (p, n, m, m1) begin
    manager(p,m)
    manager(m, m1)
    manager(m1, m)
    names(p, n)
  end

  A = DBInterface.execute(db, to_sql(q)) |> DataFrame
  @test A isa DataFrame
  @test ["Alice Smith",
         "Bob Jones",
         "John Doe"] == A[!, "n"]
  @test draw_query(q) isa Catlab.Graphics.Graphviz.Graph
end


text:
mutable struct FastaIterator{A}
    reader::FASTA.Reader # UnionAll type, but whatever
    record::FASTA.Record
    seq::LongSequence{A}

    function FastaIterator{A}(reader, record, seq) where {A <: BioSequences.Alphabet, R}
        new(reader, record, seq)
    end
end

function FastaIterator{A}(reader::FASTA.Reader) where {A <: BioSequences.Alphabet}
    return FastaIterator{A}(reader, FASTA.Record(), LongSequence{A}())
end

function FastaIterator{A}(io::IO) where {A <: BioSequences.Alphabet}
    return FastaIterator{A}(FASTA.Reader(io))
end

function FastaIterator{A}(path::AbstractString) where {A <: BioSequences.Alphabet}
    return FastaIterator{A}(FASTA.Reader(open(path)))
end

Base.eltype(::Type{FastaIterator{A}}) where A = LongSequence{A}
Base.IteratorEltype(::Type{FastaIterator{A}}) where A = Base.HasEltype()
Base.IteratorSize(::Type{<:FastaIterator}) = Base.SizeUnknown()

function Base.iterate(it::FastaIterator, ::Nothing=nothing)
    if eof(it.reader)
        close(it.reader)
        return nothing
    end
    read!(it.reader, it.record)
    copy!(it.seq, it.record)
    return it.seq, nothing
end

struct CanonicalKmerIterator{M}
    it::M
end

function CanonicalKmerIterator(it::M) where {M <: BioSequences.AbstractMerIterator}
    return CanonicalKmerIterator{M}(it)
end

Base.IteratorEltype(::Type{CanonicalKmerIterator{M}}) where M = Base.HasEltype()
Base.IteratorSize(::Type{CanonicalKmerIterator{M}}) where M = Base.IteratorSize(M)
Base.length(x::CanonicalKmerIterator{M}) where M = length(x.it)

@inline function Base.iterate(it::CanonicalKmerIterator, s...)
    itval = iterate(it.it, s...)
    itval === nothing && return nothing
    itval = itval::Tuple{Any, Any} # hack from Base.generator to help inference
    return canonical(itval[1]), itval[2]
end

mutable struct KmerSketcher{M,F}
    hasher::MinHasher{F}
    bases::Int

    function KmerSketcher{M,F}(h::MinHasher{F}) where {M <: Mer, F}
        new{M,F}(h, 0)
    end
end

KmerSketcher{M}(h::MinHasher{F}) where {M,F} = KmerSketcher{M,F}(h)
function KmerSketcher{M,F}(s::Integer) where {M <: Mer, F}
    KmerSketcher{M,F}(MinHasher(s))
end

KmerSketcher{M}(s::Integer) where {M <: Mer} = KmerSketcher{M,hash}(s)

function Base.empty!(sk::KmerSketcher)
    empty!(sk.hasher)
    sk.bases = 0
    return sk
end

function update!(sketcher::KmerSketcher{M}, seq::BioSequence) where M
    it = CanonicalKmerIterator(each(M, seq))
    MinHash.update!(sketcher.hasher, it)
    sketcher.bases += length(seq)
    return sketcher
end

function update!(sketcher::KmerSketcher, io::IO)
    it = FastaIterator{DNAAlphabet{4}}(io)
    for seq in it
        update!(sketcher, seq)
    end
    return sketcher
end

# Parameterize to allow different hash functions. Rarely used, so there are few convenience
# functions for this feature.
struct KmerHashes{K,F}
    bases::Int
    sketch::MinHash.MinHashSketch
end

# This also exists in BioSequences, but is an internal method
ksize(::Type{<:AbstractMer{A,K}}) where {A,K} = K

function KmerHashes(sk::KmerSketcher{M,F}) where {M,F}
    return KmerHashes{ksize(M),F}(sk.bases, MinHash.MinHashSketch(sk.hasher))
end

function kmer_minhash(seq::Union{IO, BioSequence}, s::Integer, ::Val{K}) where K
    sketcher = KmerSketcher{DNAMer{K}}(s)
    update!(sketcher, seq)
    return KmerHashes(sketcher)
end

function kmer_minhash_each(it, s::Integer, ::Val{K}) where K
    sketcher = KmerSketcher{DNAMer{K}}(s)
    result = KmerHashes{K,hash}[]
    for i in it
        update!(sketcher, i)
        push!(result, KmerHashes(sketcher))
        empty!(sketcher)
    end
    return result
end

function kmer_minhash_each(io::IO, s::Integer, v::Val)
    it = FastaIterator{DNAAlphabet{4}}(io)
    return kmer_minhash_each(it, s, v)
end

function kmer_minhash_paths(paths::Vector{String}, s::Integer, v::Val{K}) where K
    processes = Task[]
    for path in paths
        push!(processes, Threads.@spawn open(x -> kmer_minhash(x, s, v), path))
    end
    result = KmerHashes{K,hash}[]
    for process in processes
        push!(result, fetch(process))
    end
    return result
end

function intersectionlength(x::T, y::T) where {T<:KmerHashes}
    return MinHash.intersectionlength(x.sketch, y.sketch)
end

function intersectionlength(v::Vector{KmerHashes{K,F}}) where {K,F}
    v2 = [i.sketch for i in v]
    return MinHash.intersectionlength(v2)
end

function mash_length(h1::KmerHashes{K,F}, h2::KmerHashes{K,F}) where {K,F}
    w = intersectionlength(h1, h2)
    j = w / (length(h1.sketch) + length(h2.sketch) - w)
    return iszero(w) ? 1.0 : -log(2j / (j+1)) / k
end


text:
module Models
export ModelParams

using FactorizationMachines: FMMatrix, FMFloat

abstract ModelParams

immutable GaussianModelParams <: ModelParams
    k0::Bool
    k1::Bool
    num_factors::Int

    mean::Float64
    stddev::Float64

    GaussianModelParams(; k0 = true, k1 = true, num_factors = 8, mean = .0, stddev = .01) = new(k0, k1, num_factors, mean, stddev)
end
const gauss = GaussianModelParams

function init(params::GaussianModelParams, X::FMMatrix, y::Vector{FMFloat})
    # initialization
    num_attributes, num_samples = size(X)
    # sanity check
    assert(length(y) == num_samples)

    # create initial model
    w0 = .0
    w = zeros(num_attributes)
    V = randn(params.num_factors, num_attributes) .* params.stddev + params.mean

    # new model
    model = FMModel(params.k0, params.k1, w0, w, V, params.num_factors)
end

type FMModel
    k0::Bool
    k1::Bool

    w0::FMFloat
    w::Vector{FMFloat}
    V::Matrix{FMFloat}

    num_factors::Int
end

function predict_instance!(model::FMModel,
                           idx::StridedVector{Int64}, x::StridedVector{FMFloat},
                           f_sum::Vector{FMFloat}, sum_sqr::Vector{FMFloat})
    fill!(f_sum, .0)
    fill!(sum_sqr, .0)
    result = zero(FMFloat)
    if model.k0
        result += model.w0
    end
    if model.k1
        for i in 1:length(idx)
            result += model.w[idx[i]] * x[i]
        end
    end
    for f in 1:model.num_factors
        for i in 1:length(idx)
            d = model.V[f,idx[i]] * x[i]
            f_sum[f] += d
            sum_sqr[f] += d * d
        end
        result += 0.5 * (f_sum[f] * f_sum[f] - sum_sqr[f])
    end
    result
end

end


text:
function clone{T<:AbstractMolecularContainer}(x::T, new_props::Dict=Dict())
    return T(;merge(features(x), new_props)...)
end


text:
using ContMechTensors
using Base.Test

include("test_misc.jl")
include("test_ops.jl")
include("test_ad.jl")

# Build the docs
include("../docs/make.jl")