text:
using Plots
using LinearAlgebra
using Test
using FLoops

pyplot()

using Images

function Base.isapprox(n1::Nothing, n2::Nothing)
    return true
end

# For now, all cameras are aligned on the z axis
struct Camera
    # Set of all pixels, counts as scene resolution
    pixels

    # physical size of aperture
    size::Vector{Float64}

    # camera's distance from screen
    focal_length::Float64

    # camera's position
    p::Vector{Float64}
end

struct Ray
    # Velocity vector
    l::Vector{Float64}

    # Position vector
    p::Vector{Float64}

    # Color
    c::RGB

end

struct Surface

    # Reflectivity
    r::Float64

    # Transmission
    t::Float64

    # Color
    c::RGBA

    # index of refraction
    ior::Float64

    function Surface(in_r, in_t, in_c, in_ior)
        if !isapprox(in_r+in_t+in_c.alpha, 1)
            error("invalid surface definition, RTC < 1")
        end
        new(in_r,in_t,in_c, in_ior)
    end

    Surface(in_r, in_t, in_c::Float64, in_ior) =
         new(in_r, in_t, RGBA(0,0,0,0), in_ior)
end

abstract type Object end

struct Sphere <: Object
    # Lens position
    p::Vector{Float64}

    # Lens radius
    r::Float64

    s::Surface
end

function Lens(p, r, ior)
    return Sphere(p, r, Surface(0,1,RGBA(0,0,0,0),ior))
end

function ReflectingSphere(p, r)
    return Sphere(p,r,Surface(1,0,RGBA(0,0,0,0),0))
end

function ColoredSphere(p, r, c::RGB)
    return Sphere(p, r, Surface(0,0,RGBA(c), 0))
end

struct SkyBox <: Object
    # Skybox position
    p::Vector{Float64}

    # Skybox radius
    r::Float64
end

function sphere_normal_at(ray, sphere)
    n = normalize(ray.p .- sphere.p)

    return n
end

# note: light moves as a particular speed with respect to the medium it is
#       moving through, so...
#       n_2*v = n_1*l + (n_1*cos(theta_1) - n_2*cos(theta_2))*n
#       Other approximations: ior = n_1/n_2, c = -n*l
function refract(ray, lens::Sphere, ior)
    n = sphere_normal_at(ray, lens)

    if dot(n, ray.l) > 0
        n .*= -1
    end
    c = dot(-n, ray.l);
    d = 1.0 - ior^2 * (1.0 - c^2);

    if (d < 0.0)
        return reflect(ray, n)
    end

    ray_vel = ior * ray.l + (ior * c - sqrt(d)) * n;
    return Ray(ray_vel, ray.p, ray.c)
end

# note: for reflection, l_x -> l_x, but l_y -> -l_y
#       In this case, the y component of l = cos(theta)*n
#       so new vector: v = l + 2cos(theta)*n
function reflect(ray, n)
    ray_vel = ray.l .- 2*dot(ray.l, n).*n
    ray_pos = ray.p .+ 0.001*ray.l
    return Ray(ray_vel, ray_pos, ray.c)
end

function intersection(ray::Ray, sphere::S;
                      threshold = 0.01) where
                      {S <: Union{Sphere, SkyBox}}
    relative_dist = ray.p-sphere.p
    a = dot(ray.l, ray.l)
    b = 2.0 * dot(relative_dist, ray.l)
    c = dot(relative_dist, relative_dist) - sphere.r*sphere.r
    discriminant = b*b - 4*a*c

    if discriminant < 0
        return nothing
    elseif discriminant > 0
        roots = [(-b + sqrt(discriminant)) / (2*a),
                 (-b - sqrt(discriminant)) / (2*a)]
        min = minimum(roots)
        max = maximum(roots)

        if min > threshold
            return (min)*ray.l
        elseif max > threshold
            return (max)*ray.l
        else
            return nothing
        end
    else
        # Returns nothing if tangential
        return nothing
        #return (-b/(2*a))*ray.l
    end
end

function propagate(rays::Array{Ray}, objects::Vector{O},
                    num_intersections) where {O <: Object}
    @floop ThreadedEx() for j = 1:length(rays)
        rays[j] = propagate(rays[j], objects, num_intersections)
    end

    return rays
end

function propagate(ray::Ray, objects::Vector{O},
                   num_intersections) where {O <: Object}

    for i = 1:num_intersections
        if ray.l != zeros(length(ray.l))
            intersect_final = [Inf, Inf]
            intersected_object = nothing
            for object in objects
                intersect = intersection(ray, object)
                if intersect != nothing &&
                   sum(intersect[:].^2) < sum(intersect_final[:].^2)
                    intersect_final = intersect
                    intersected_object = object
                end
            end

            if intersect_final != nothing
                ray = Ray(ray.l, ray.p .+ intersect_final, ray.c)
                if typeof(intersected_object) == Sphere
                    reflected_ray = ray
                    refracted_ray = ray
                    colored_ray = ray
                    if !isapprox(intersected_object.s.t, 0)
                        ior = 1/intersected_object.s.ior
                        if dot(ray.l,
                               sphere_normal_at(ray,
                                                intersected_object)) > 0
                            ior = intersected_object.s.ior
                        end

                        refracted_ray = refract(ray, intersected_object, ior)
                        refracted_ray = propagate(refracted_ray, objects,
                                                  num_intersections-i)
                    end

                    if !isapprox(intersected_object.s.r, 0)
                        n = sphere_normal_at(ray, intersected_object)
                        reflected_ray = reflect(ray, n)
                        reflected_ray = propagate(reflected_ray, objects,
                                                  num_intersections-i)
                    end

                    if !isapprox(intersected_object.s.c.alpha, 0)
                        ray_color = RGB(intersected_object.s.c)
                        ray_vel = zeros(length(ray.l))
                        colored_ray = Ray(ray_vel, ray.p, ray_color)
                    end

                    ray_color = intersected_object.s.t*refracted_ray.c +
                                intersected_object.s.r*reflected_ray.c +
                                intersected_object.s.c.alpha*colored_ray.c

                    ray = Ray(zeros(length(ray.l)), ray.p, ray_color)

                elseif typeof(intersected_object) == Mirror
                    ray = reflect(ray, intersected_object.n)
                elseif typeof(intersected_object) == SkyBox
                    ray_color = pixel_color(ray.p, 1000)
                    ray_vel = zeros(length(ray.l))
                    ray = Ray(ray_vel, ray.p, ray_color)
                end
            else
                println("hit nothing")
            end
        end
    end

    return ray
end

function pixel_color(position, extents)
    c = RGB(0)
    if position[1] < extents && position[1] > -extents
        c += RGB((position[1]+extents)/(2.0*extents), 0, 0)
    else
        println(position)
    end

    if position[2] < extents && position[2] > -extents
        c += RGB(0,0,(position[2]+extents)/(2.0*extents))
    else
        println(position)
    end

    if position[3] < extents && position[3] > -extents
        c += RGB(0,(position[3]+extents)/(2.0*extents), 0)
    else
        println(position)
    end

    return c
end

function convert_to_img(rays::Array{Ray}, filename)
    color_array = Array{RGB}(undef, size(rays)[2], size(rays)[1])
    for i = 1:length(color_array)
         color_array[i] = rays[i].c
    end

    save(filename, color_array)
end

function init_rays(cam::Camera)

    res = size(cam.pixels)
    dim = cam.size

    pixel_width = dim ./ res

    # create a set of rays that go through every pixel in our grid.
    rays = Array{Ray}(undef, res[1], res[2])
    for i = 1:res[1]
        for j = 1:res[2]
            pixel_loc = [cam.p[1] + 0.5*dim[1] - i*dim[1]/res[1] +
                         0.5*pixel_width[1],
                         cam.p[2] + 0.5*dim[2] - j*dim[2]/res[2] +
                         0.5*pixel_width[2],
                         cam.p[3]+cam.focal_length]
            l = normalize(pixel_loc - cam.p)
            rays[res[2]*(i-1) + j] = Ray(l, pixel_loc, RGB(0))
        end
    end

    return rays

end

function ray_trace(objects::Vector{O}, cam::Camera; filename="check.png",
                   num_intersections = 10) where {O <: Object}

    rays = init_rays(cam)

    rays = propagate(rays, objects, num_intersections)

    convert_to_img(rays, filename)

    return rays

end

function main()
    sky = [SkyBox([0.0, 0.0, 0.0], 1000)]
    spheres = [ColoredSphere([-50,0,-25], 20, RGB(0, 0, 1))]
    #spheres = [Lens([0,0,-25], 20, 1.5), ReflectingSphere([0,50,-100],20),
    #           ColoredSphere([-50,0,-25], 20, RGB(0, 0, 1))]
#=
    spheres = [Lens([0,0,-25], 20, 1.5), ReflectingSphere([0,50,-100],20),
               ColoredSphere([-50,0,-25], 20, RGB(0, 0, 1)),
               Sphere([30, 25, -60], 20,
                      Surface(0.0, 0.75, RGBA(1,0,0,0.25), 1.5)),
               Sphere([50, 0, -25], 20,
                      Surface(0.5, 0.0, RGBA(0,1,0,0.5), 1.5)),
               Sphere([-30, 25, -60], 20,
                      Surface(0.5, 0.5, RGBA(1,1,1,0), 1.5))]
=#
    #bg_spheres = sphere_rand(10, [-400 400; -400 400; -350 -400], 50)

    objects = vcat(sky, spheres)

    blank_img = Array{RGB}(undef, 1920, 1080)
    #blank_img = Array{RGB}(undef, 20, 10)
    blank_img[:] .= RGB(0)

    cam = Camera(blank_img, [160,90], -100, [0,20,100])

    return ray_trace(objects, cam; num_intersections=20)

end


text:
using NESSie: potprefactor

@testset "potprefactor" begin
    for T in testtypes
        @test typeof(potprefactor(T)) == T
    end
end

@testset "defaultopt" begin
    for T in testtypes
        @test typeof(defaultopt(T)) == Option{T}
    end
end

ext: jl


text:
@testset "Utilities" begin
    # create visualizer
    vis = Dojo.Visualizer();

    # default background
    Dojo.set_background!(vis)

    # add floor
    Dojo.set_floor!(vis)

    # add quadratic surface
    Dojo.set_surface!(vis, x -> x' * x,
        color=Dojo.RGBA(1.0, 0.0, 0.0, 1.0))

    # lighting
    Dojo.set_light!(vis)

    # camera
    Dojo.set_camera!(vis,
        zoom=0.25)

    # test that methods don't fail
    @test true
end

@testset "MeshCat mechanism" begin
    # create visualizer
    vis = Dojo.Visualizer();
    # get environment and simulate
    env = DojoEnvironments.get_environment("halfcheetah";
        timestep=0.1)
    reset(env)
    storage = Dojo.simulate!(env.mechanism, 0.25,
        record=true,
        verbose=true);
    # visualize simulation
    Dojo.visualize(env.mechanism, storage,
        vis=vis)

    # test that methods don't fail
    @test true
end

@testset "URDF mesh" begin
    # create visualizer
    vis = Dojo.Visualizer();
    # get environment and simulate
    env = DojoEnvironments.get_environment("quadruped"; timestep=0.1)
    reset(env)
    storage = Dojo.simulate!(env.mechanism, 0.25,
        record=true,
        verbose=false)
    # visualize simulation
    Dojo.visualize(env.mechanism, storage,
        vis=vis)
    # visualize w/ contact
    Dojo.visualize(env.mechanism, storage,
        vis=vis,
        show_contact=true)

    # test that methods don't fail
    @test true
end


text:
# This file was generated, do not modify it. # hide
box_model = UnivariateBoxCoxTransformer()
box = machine(box_model, ys)
z = transform(box, ys)

ridge_model = RidgeRegressor(lambda=0.1)
ridge = machine(ridge_model, W, z)
ẑ = predict(ridge, W)

ŷ = inverse_transform(box, ẑ)


text:
include("../utils.jl")
include("group_utils.jl")

@testset "Special Euclidean group" begin
    @testset "SpecialEuclidean($n)" for n in (2, 3, 4)
        G = SpecialEuclidean(n)
        @test isa(G, SpecialEuclidean{n})
        @test repr(G) == "SpecialEuclidean($n)"
        M = base_manifold(G)
        @test M === TranslationGroup(n) × SpecialOrthogonal(n)
        @test submanifold(G, 1) === TranslationGroup(n)
        @test submanifold(G, 2) === SpecialOrthogonal(n)
        Rn = Rotations(n)
        x = Matrix(I, n, n)

        if n == 2
            t = Vector{Float64}.([1:2, 2:3, 3:4])
            ω = [[1.0, 2.0], [2.0, 1.0], [1.0, 3.0]]
            tuple_pts = [(ti, exp(Rn, x, hat(Rn, x, ωi))) for (ti, ωi) in zip(t, ω)]
            tuple_v = ([-1.0, 2.0], hat(Rn, x, [1.0, -0.5]))
        elseif n == 3
            t = Vector{Float64}.([1:3, 2:4, 4:6])
            ω = [[1.0, 2.0, 3.0], [3.0, 2.0, 1.0], [1.0, 3.0, 2.0]]
            tuple_pts = [(ti, exp(Rn, x, hat(Rn, x, ωi))) for (ti, ωi) in zip(t, ω)]
            tuple_v = ([-1.0, 2.0, 1.0], hat(Rn, x, [1.0, 0.5, -0.5]))
        else # n == 4
            t = Vector{Float64}.([1:4, 2:5, 3:6])
            ω = [
                [1.0, 2.0, 3.0, 1.0, 2.0, 3.0],
                [3.0, 2.0, 1.0, 1.0, 2.0, 3.0],
                [1.0, 3.0, 2.0, 1.0, 2.0, 3.0],
            ]
            tuple_pts = [(ti, exp(Rn, x, hat(Rn, x, ωi))) for (ti, ωi) in zip(t, ω)]
            tuple_v = ([-1.0, 2.0, 1.0, 3.0], hat(Rn, x, [1.0, 0.5, -0.5, 0.0, 2.0, 1.0]))
        end

        @testset "product point" begin
            reshapers = (Manifolds.ArrayReshaper(), Manifolds.StaticReshaper())
            for reshaper in reshapers
                shape_se = Manifolds.ShapeSpecification(reshaper, M.manifolds...)
                pts = [Manifolds.prod_point(shape_se, tp...) for tp in tuple_pts]
                v_pts = [Manifolds.prod_point(shape_se, tuple_v...)]

                g1, g2 = pts[1:2]
                t1, R1 = g1.parts
                t2, R2 = g2.parts
                g1g2 = Manifolds.prod_point(shape_se, R1 * t2 + t1, R1 * R2)
                @test isapprox(G, compose(G, g1, g2), g1g2)
                @test affine_matrix(G, g1g2) ≈ affine_matrix(G, g1) * affine_matrix(G, g2)
                tmp = copy(g1)
                Manifolds._padpoint!(G, tmp)
                @test tmp == g1
                tmp = copy(v_pts[1])
                Manifolds._padvector!(G, tmp)
                @test tmp == v_pts[1]

                w = translate_diff(G, pts[1], make_identity(G, pts[1]), v_pts[1])
                w2 = allocate(w)
                w2.parts[1] .= w.parts[1]
                w2.parts[2] .= pts[1].parts[2] * w.parts[2]
                @test screw_matrix(G, w2) ≈
                      affine_matrix(G, pts[1]) * screw_matrix(G, v_pts[1])

                test_group(
                    G,
                    pts,
                    v_pts,
                    v_pts;
                    test_diff = true,
                    diff_convs = [(), (LeftAction(),)],
                )
            end
        end

        @testset "product repr" begin
            pts = [ProductRepr(tp...) for tp in tuple_pts]
            v_pts = [ProductRepr(tuple_v...)]

            g1, g2 = pts[1:2]
            t1, R1 = g1.parts
            t2, R2 = g2.parts
            g1g2 = ProductRepr(R1 * t2 + t1, R1 * R2)
            @test isapprox(G, compose(G, g1, g2), g1g2)
            g1g2mat = affine_matrix(G, g1g2)
            @test g1g2mat ≈ affine_matrix(G, g1) * affine_matrix(G, g2)
            @test affine_matrix(G, g1g2mat) === g1g2mat
            @test affine_matrix(G, make_identity(G, pts[1])) isa SDiagonal{n,Float64}
            @test affine_matrix(G, make_identity(G, pts[1])) == SDiagonal{n,Float64}(I)

            w = translate_diff(G, pts[1], make_identity(G, pts[1]), v_pts[1])
            w2 = allocate(w)
            w2.parts[1] .= w.parts[1]
            w2.parts[2] .= pts[1].parts[2] * w.parts[2]
            w2mat = screw_matrix(G, w2)
            @test w2mat ≈ affine_matrix(G, pts[1]) * screw_matrix(G, v_pts[1])
            @test screw_matrix(G, w2mat) === w2mat

            test_group(
                G,
                pts,
                v_pts,
                v_pts;
                test_diff = true,
                diff_convs = [(), (LeftAction(),)],
            )
        end

        @testset "affine matrix" begin
            pts = [affine_matrix(G, ProductRepr(tp...)) for tp in tuple_pts]
            v_pts = [screw_matrix(G, ProductRepr(tuple_v...))]
            test_group(
                G,
                pts,
                v_pts,
                v_pts;
                test_diff = true,
                diff_convs = [(), (LeftAction(),)],
            )
        end

        @testset "hat/vee" begin
            shape_se =
                Manifolds.ShapeSpecification(Manifolds.ArrayReshaper(), M.manifolds...)
            x = Manifolds.prod_point(shape_se, tuple_pts[1]...)
            V = Manifolds.prod_point(shape_se, tuple_v...)
            vexp = [V.parts[1]; vee(Rn, x.parts[2], V.parts[2])]
            v = vee(G, x, V)
            @test v ≈ vexp
            @test hat(G, x, v) ≈ V

            v = vee(G, affine_matrix(G, x), screw_matrix(G, V))
            @test v ≈ vexp
            @test hat(G, affine_matrix(G, x), v) ≈ screw_matrix(G, V)
        end
    end

    G = SpecialEuclidean(11)
    @test affine_matrix(G, make_identity(G, ones(12, 12))) isa
          Diagonal{Float64,Vector{Float64}}
    @test affine_matrix(G, make_identity(G, ones(12, 12))) == Diagonal(ones(11))
end


text:
function asterplot(; x::AbstractVector = [], y::AbstractVector = Int[], weight::AbstractVector = Float64[], holesize = 0)

    v = piechart(x = x, y = y, holesize = holesize)

    #Shim weights into y2 slot
    for d in enumerate(v.data[1].values)
        if weight != Float64[]
            d[2]["y2"] = weight[d[1]]
        else
            d[2]["y2"] = 1
        end
    end

    v.marks[1].from.transform = [VegaTransform(Dict{Any, Any}("type" => "pie", "field" => "y2")),
                             VegaTransform(Dict{Any, Any}("type" => "formula", "expr" => "($(v.height) - $(v.marks[1].properties.enter.innerRadius.value)) * (datum.y/100.0) + $(v.marks[1].properties.enter.innerRadius.value)", "field" => "calcOuterRadius"))
                            ]

    v.marks[1].properties.enter.outerRadius = VegaValueRef(field = "calcOuterRadius")

    #make a copy of the entire properties, then change to make gray background
    push!(v.marks, deepcopy(v.marks[1]))

    #Change stroke to gray, remove fill, make outerRadius total height
    v.marks[2].properties.enter.stroke = VegaValueRef(value = "gray")
    v.marks[2].properties.enter.fill = nothing
    v.marks[2].properties.enter.outerRadius = VegaValueRef(value = v.height)

    return v

end

max_stars_repo_name: jpthor/vega


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

using Base.Test
using FEMBasis

@testset "interpolate scalar field" begin
    # in unit square: T(X,t) = t*(1 + X[1] + 3*X[2] - 2*X[1]*X[2])
    B = Quad4()
    X = ((0.0, 0.0), (1.0, 0.0), (1.0, 1.0), (0.0, 1.0))
    T = (1.0, 2.0, 3.0, 4.0)
    T_known(X) = 1 + X[1] + 3*X[2] - 2*X[1]*X[2]
    T_interpolated = interpolate(B, T, (0.0, 0.0))
    @test isapprox(T_interpolated, T_known((0.5, 0.5)))
end

@testset "calculate gradient dB/dX" begin
    B = Quad4()
    X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
    dBdX = grad(B, X, (0.0, 0.0))
    @test isapprox(dBdX, 1/2*[-1 1 1 -1; -1 -1 1 1])
end

@testset "interpolate gradient of scalar field" begin
    # in unit square: grad(T)(X) = [1-2X[2], 3-2*X[1]]
    B = Quad4()
    X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
    T = (1.0, 2.0, 3.0, 4.0)
    dTdX = grad(B, T, X, (0.0, 0.0))
    dTdX_expected(X) = [1-2*X[2] 3-2*X[1]]
    @test isapprox(dTdX, dTdX_expected([0.5, 0.5]))
end

@testset "interpolate vector field" begin
    # in unit square, u(X,t) = [1/4*t*X[1]*X[2], 0, 0]
    B = Quad4()
    u = Vector[[0.0, 0.0], [0.0, 0.0], [1/4, 0.0], [0.0, 0.0]]
    u_known(X) = [1/4*X[1]*X[2], 0]
    u_interpolated = interpolate(B, u, (0.0, 0.0))
    @test isapprox(u_interpolated, u_known((0.5, 0.5)))
end

@testset "interpolate gradient of vector field" begin
    # in unit square, u(X) = t*[X[1]*(X[2]+1), X[1]*(4*X[2]-1)]
    # => u_i,j = t*[X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
    B = Quad4()
    X = ([0.0, 0.0], [1.0, 0.0], [1.0, 1.0], [0.0, 1.0])
    u = ([0.0, 0.0], [1.0, -1.0], [2.0, 3.0], [0.0, 0.0])
    dudX = grad(B, u, X, (0.0, 0.0))
    dudX_expected(X) = [X[2]+1 X[1]; 4*X[2]-1 4*X[1]]
    @test isapprox(dudX, dudX_expected([0.5, 0.5]))
end

@testset "test BasisInfo" begin
    B = BasisInfo(Seg2)
    @test length(B) == 2
    @test size(B) == (1, 2)
    X = ((0.0,), (1.1))
    xi = (0.0,)
    eval_basis!(B, X, xi)
    @test isapprox(B.invJ, inv(B.J))
    @test isapprox(B.detJ, det(B.J))
    B = BasisInfo(Quad4)
    X = ((0.0,0.0), (1.1,0.0), (1.0,1.0), (0.0,1.0))
    xi = (0.0,0.0)
    eval_basis!(B, X, xi)
    @test isapprox(B.invJ, inv(B.J))
    @test isapprox(B.detJ, det(B.J))
    B = BasisInfo(Hex8)
    X = ((0.0,0.0,0.0), (1.1,0.0,0.0), (1.0,1.0,0.0), (0.0,1.0,0.0),
         (0.0,0.0,1.0), (1.0,0.0,1.0), (1.0,1.0,1.0), (0.0,1.0,1.0))
    xi = (0.0,0.0,0.0)
    eval_basis!(B, X, xi)
    @test isapprox(B.invJ, inv(B.J))
    @test isapprox(B.detJ, det(B.J))
end

@testset "evaluate gradient using BasisInfo" begin
    B = BasisInfo(Quad4)
    X = ((0.0,0.0), (1.0,0.0), (1.0,1.0), (0.0,1.0))
    xi = (0.0, 0.0)
    eval_basis!(B, X, xi)
    u = ((0.0, 0.0), (1.0, -1.0), (2.0, 3.0), (0.0, 0.0))
    gradu = zeros(2, 2)
    grad!(B, gradu, u)
    gradu_expected = [1.5 0.5; 1.0  2.0]
    @test isapprox(gradu, gradu_expected)
end

@testset "test curves" begin
    bi = BasisInfo(Seg2)
    X1 = ((0.0,0.0,0.0), (1.0,0.0,0.0))
    X2 = ((0.0,0.0,0.0), (0.0,1.0,0.0))
    X3 = ((0.0,0.0,0.0), (0.0,0.0,1.0))
    xi = (0.0, 0.0)
    eval_basis!(bi, X1, xi)
    @test isapprox(bi.detJ, 0.5)
    eval_basis!(bi, X2, xi)
    @test isapprox(bi.detJ, 0.5)
    eval_basis!(bi, X3, xi)
    @test isapprox(bi.detJ, 0.5)
    @test isapprox(jacobian(Seg2(), X1, xi), [0.5 0.0 0.0])
    @test isapprox(jacobian(Seg2(), X2, xi), [0.0 0.5 0.0])
    @test isapprox(jacobian(Seg2(), X3, xi), [0.0 0.0 0.5])
end

@testset "test manifolds" begin
    bi = BasisInfo(Quad4)
    X1 = ((0.0,0.0,0.0), (1.0,0.0,0.0), (1.0,1.0,0.0), (0.0,1.0,0.0))
    X2 = ((0.0,0.0,0.0), (0.0,1.0,0.0), (0.0,1.0,1.0), (0.0,0.0,1.0))
    X3 = ((0.0,0.0,0.0), (0.0,0.0,1.0), (1.0,0.0,1.0), (1.0,0.0,0.0))
    xi = (0.0, 0.0)
    eval_basis!(bi, X1, xi)
    @test isapprox(bi.detJ, 0.25)
    eval_basis!(bi, X2, xi)
    @test isapprox(bi.detJ, 0.25)
    eval_basis!(bi, X3, xi)
    @test isapprox(bi.detJ, 0.25)
    J1 = jacobian(Quad4(), X1, xi)
    J2 = jacobian(Quad4(), X2, xi)
    J3 = jacobian(Quad4(), X3, xi)
    @test isapprox(J1, [0.5 0.0 0.0; 0.0 0.5 0.0])
    @test isapprox(J2, [0.0 0.5 0.0; 0.0 0.0 0.5])
    @test isapprox(J3, [0.0 0.0 0.5; 0.5 0.0 0.0])
end


text:
@testset "char freq from text" begin
    test_text = "十七道身影，十七道白色的身影，宛如星丸跳跃一般从山腰处朝山顶方向而来，这十七道身影的主人，年纪最小的也超过了五旬，一个个神色凝重，他们身穿的白袍代表的是内门，而胸前那金色的唐字则是唐门长老的象征。"
    cf = charfreq(test_text)

    @test "十" in keys(cf) && cf["十"] == 3
    @test length(cf) == 65
end

@testset "char freq from iterable" begin
    tokens = ["a", "b", "a", "c", "b", "a", "c", "b", "a"]
    cf = charfreq(tokens)

    @test "a" in keys(cf) && cf["a"] == 4
    @test "b" in keys(cf) && cf["b"] == 3
    @test "c" in keys(cf) && cf["c"] == 2

    @test length(cf) == 3
    @test size(cf) == length(tokens)
end


text:
# Copyright 2017-2018 Gandalf Software, Inc. (Scott Paul Jones)
# Licensed under MIT License, see LICENSE.md

# These use the "Holy Trait Trick", which was created by @timholy (Tim Holy),
# which was made possible by the type system in Julia from @jeff.bezanson (Jeff Bezanson).

# One is proud to stand on the shoulders of such people!

## Traits for string and character types ##

"""
    ValidatedStyle(A)
    ValidatedStyle(typeof(A))

`ValidatedStyle` specifies the whether a string or character type is always valid or not
When you define a new `AbstractString` or `AbstractChar` type, you can choose to implement it
as always validated, or validation state unknown.

    Strs.ValidatedStyle(::Type{<:MyString}) = AwaysValid()

The default is `UnknownValidity`
"""
abstract type ValidatedStyle end
struct AlwaysValid     <: ValidatedStyle end
struct UnknownValidity <: ValidatedStyle end

ValidatedStyle(::Type{<:AbstractChar}) = UnknownValidity()
@static isdefined(Base, :AbstractChar) || (ValidatedStyle(::Type{Char}) = UnknownValidity())
ValidatedStyle(::Type{<:Chr}) = AlwaysValid()

ValidatedStyle(::Type{<:AbstractString}) = UnknownValidity()
ValidatedStyle(::Type{<:Str}) = AlwaysValid()

ValidatedStyle(A::T) where {T<:Union{AbsChar,AbstractString}} = ValidatedStyle(T)

# single or multiple codeunits per codepoint

"""
    CodePointStyle(A)
    CodePointStyle(typeof(A))

`CodePointStyle` specifies the whether a string type uses one or multiple codeunits to
encode a single codepoint.
When you define a new `AbstractString` type, you can choose to implement it with either
single or multi-codeunit indexing.

    Strs.CodePointStyle(::Type{<:MyString}) = CodeUnitMulti()

The default is `CodeUnitMulti()`
"""
abstract type CodePointStyle end

struct CodeUnitSingle <: CodePointStyle end
struct CodeUnitMulti  <: CodePointStyle end

CodePointStyle(::Type{<:CSE}) = CodeUnitSingle()
CodePointStyle(::Type{<:Union{UTF8_CSEs,UTF16CSE}}) = CodeUnitMulti()

CodePointStyle(::Type{T}) where {T<:AbstractString} = CodePointStyle(cse(T))

CodePointStyle(::Type{<:SubString{T}}) where {T<:AbstractString} = CodePointStyle(T)
CodePointStyle(v::AbstractString) = CodePointStyle(typeof(v))

is_multi(str::AbstractString) = CodePointStyle(str) === CodeUnitMulti()

# Type of character set

"""
    CharSetStyle(A)
    CharSetStyle(typeof(A))

`CharSetStyle` specifies the information about the character set used by the string or
characters.
When you define a new `AbstractString` or `AbstractChar` type,
you can choose to implement it with

    Strs.CharSetStyle(::Type{<:MyString}) = CharSetISOCompat()

The default is `CharSetUnicode()`
"""
abstract type CharSetStyle end

"""Codepoints are not in Unicode compatible order"""
struct CharSetOther         <: CharSetStyle end
"""Characters 0-0x7f same as ASCII"""
struct CharSetASCIICompat   <: CharSetStyle end
"""Characters 0-0x9f follows ISO 8859"""
struct CharSetISOCompat     <: CharSetStyle end
"""Characters 0-0xd7ff, 0xe000-0xffff follow Unicode BMP"""
struct CharSetBMPCompat     <: CharSetStyle end
"""Full Unicode character set, no additions"""
struct CharSetUnicode       <: CharSetStyle end
"""Unicode character set, plus encodings of invalid characters"""
struct CharSetUnicodePlus   <: CharSetStyle end
"""8-bit Binary string, not text"""
struct CharSetBinary        <: CharSetStyle end
"""Raw bytes, words, or character string, unknown encoding/character set"""
struct CharSetUnknown       <: CharSetStyle end

CharSetStyle(A::AbstractString) = CharSetStyle(typeof(A))

CharSetStyle(::Type{<:AbstractString})    = CharSetUnicodePlus()
CharSetStyle(::Type{<:Str{<:Union{Text1CSE, Text2CSE, Text4CSE}}}) = CharSetUnknown()
CharSetStyle(::Type{<:Str{BinaryCSE}})    = CharSetBinary()
CharSetStyle(::Type{<:Str{ASCIICSE}})     = CharSetASCIICompat()
CharSetStyle(::Type{<:Str{<:Latin_CSEs}}) = CharSetISOCompat()
CharSetStyle(::Type{<:Str{<:UCS2_CSEs}})  = CharSetBMPCompat()

CharSetStyle(A::AbstractChar)   = CharSetStyle(typeof(A))

CharSetStyle(::Type{<:AbstractChar}) = CharSetUnicode()
CharSetStyle(::Type{Char})      = CharSetUnicodePlus() # Encodes invalid characters also
CharSetStyle(::Type{UInt8})     = CharSetBinary()
CharSetStyle(::Type{ASCIIChr})  = CharSetASCIICompat()
CharSetStyle(::Type{LatinChr})  = CharSetISOCompat()
CharSetStyle(::Type{_LatinChr}) = CharSetISOCompat()
CharSetStyle(::Type{UCS2Chr})   = CharSetBMPCompat()
CharSetStyle(::Type{UInt16})    = CharSetUnknown()
CharSetStyle(::Type{UInt32})    = CharSetUnknown()

# must check range if CS1 is smaller than CS2, even if CS2 is valid for it's range
_isvalid(::ValidatedStyle, ::Type{ASCIICharSet}, ::Type{T}, val) where {T<:CharSet} = is_ascii(val)

(_isvalid(::ValidatedStyle, ::Type{LatinCharSet}, ::Type{T}, val)
 where {T<:Union{Text2CharSet, Text4CharSet, UCS2CharSet, UTF32CharSet}}) =
     is_latin(val)

(_isvalid(::ValidatedStyle, ::Type{UCS2CharSet}, ::Type{T}, val)
 where {T<:Union{Text2CharSet, Text4CharSet, UTF32CharSet}}) =
     is_bmp(val)

_isvalid(::ValidatedStyle, ::Type{UTF32CharSet}, ::Type{<:CharSet}, val) =
    is_unicode(val)

# no checking needed for cases where it is a superset of T
(_isvalid(::AlwaysValid, ::Type{LatinCharSet}, ::Type{T}, val)
  where {T<:Union{Text1CharSet,ASCIICharSet,LatinSubSet}}) = true

(_isvalid(::AlwaysValid, ::Type{UCS2CharSet}, ::Type{T}, val)
  where {T<:Union{Text1CharSet,ASCIICharSet,LatinCharSet,LatinSubSet,UCS2SubSet}}) = true

(_isvalid(::AlwaysValid, ::Type{UTF32CharSet}, ::Type{T}, val)
  where {T<:Union{Text1CharSet,ASCIICharSet,LatinCharSet,UCS2CharSet,LatinSubSet,UCS2SubSet,UTF32SubSet}}) =
    true

# no subsets allowed for these
_isvalid(::AlwaysValid, ::Type{LatinSubSet}, ::Type{ASCIICharSet}, val) = false
(_isvalid(::AlwaysValid, ::Type{UCS2SubSet}, ::Type{T}, val)
  where {T<:Union{Text1CharSet,ASCIICharSet,LatinCharSet,LatinSubSet}}) = false
(_isvalid(::AlwaysValid, ::Type{UTF32SubSet}, ::Type{T}, val)
 where {T<:Union{Text1CharSet,Text2CharSet,ASCIICharSet,LatinCharSet,
                 UCS2CharSet,LatinSubSet,UCS2SubSet}}) = false

(_isvalid(::AlwaysValid, ::Type{S}, ::Type{<:CodeUnitTypes}, chr)
 where {S<:Union{Text1CharSet,Text2CharSet,Text4CharSet}}) =
     chr <= typemax(S)

# Different character sets
is_valid(::Type{Str{<:CSE{CS1}}}, str::T) where {T<:Str{<:CSE{CS2}},CS1} where {CS2} =
    _isvalid(ValidatedStyle(T), CS1, CS2, str)

# Same character set
is_valid(::Type{Str{<:CSE{CS}}}, str::T) where {T<:Str{<:CSE{CS}}} where {CS} =
    _isvalid(ValidatedStyle(T), str)

_isvalid(::AlwaysValid, v) = true

_isvalid(::UnknownValidity, str::T) where {T<:Str} = _isvalid(T, pointer(str), ncodeunits(str))

# By default, check that it is valid Unicode codepoint
_isvalid(::UnknownValidity, v) = _isvalid(UnknownValidity(), UTF32CharSet, charset(v), v)

is_valid(::Type{T}, str::T) where {T<:Str}       = _isvalid(ValidatedStyle(T), str)
is_valid(::Type{T}, chr::T) where {T<:Chr} = _isvalid(ValidatedStyle(T), chr)

is_valid(str::T) where {T<:Str}       = _isvalid(ValidatedStyle(T), str)
is_valid(chr::T) where {T<:Chr} = _isvalid(ValidatedStyle(T), chr)

# Different character set
function is_valid(::Type{S}, chr::T) where {S<:Chr, T<:Chr}
    CS = charset(S)
    CT = charset(T)
    CS == CT ? _isvalid(ValidatedStyle(T), chr) : _isvalid(ValidatedStyle(T), CS, CT, chr)
end

# Not totally sure how to get rid of some of these, they really should be handled
# by the compiler, using the ValidatedStyle trait along with the character sets

_isvalid_chr(::Type{ASCIICharSet}, v)  = is_ascii(v)
_isvalid_chr(::Type{LatinCharSet}, v)  = is_latin(v)
_isvalid_chr(::Type{UCS2CharSet},  v)  = is_bmp(v)
_isvalid_chr(::Type{UTF32CharSet}, v)  = is_unicode(v)
_isvalid_chr(::Type{LatinSubSet}, v)   = is_latin(v)
_isvalid_chr(::Type{UCS2SubSet},  v)   = is_bmp(v)
_isvalid_chr(::Type{UTF32SubSet}, v)   = is_unicode(v)
_isvalid_chr(::Type{Text1CharSet}, v)  = v <= typemax(UInt8)
_isvalid_chr(::Type{Text2CharSet}, v)  = v <= typemax(UInt16)
_isvalid_chr(::Type{Text4CharSet}, v)  = v <= typemax(UInt32)
_isvalid_chr(::Type{BinaryCharSet}, v) = v <= typemax(UInt8)

# Not totally sure about this, base Char is rather funky in v0.7
_isvalid_chr(::Type{UniPlusCharSet}, v) = v <= typemax(UInt32)

is_valid(::Type{T}, v::Unsigned) where {T<:Chr} =
    _isvalid_chr(charset(T), v)
is_valid(::Type{T}, v::Signed) where {T<:Chr} =
    0 <= v <= typemax(UInt32) && _isvalid_chr(charset(T), v%UInt32)

is_valid(::Type{Char}, ch::Union{Text1Chr, ASCIIChr, LatinChars, UCS2Chr, UTF32Chr}) = true
is_valid(::Type{Char}, ch::Text2Chr) = is_bmp(ch)
is_valid(::Type{Char}, ch::Text4Chr) = is_unicode(ch)
is_valid(::Type{T},    ch::Char) where {T<:Chr} = Base.isvalid(ch) && is_valid(T, ch%UInt32)

# For now, there is only support for immutable `Str`s, when mutable `Str`s are added.

"""
    MutableStyle(A)
    MutableStyle(typeof(A))

`MutableStyle` specifies the whether a string type is mutable or not

    Strs.MutableStyle(::Type{<:MyString}) = MutableStr()

The default is `ImmutableStr`
"""
abstract type MutableStyle end
struct ImmutableStr <: MutableStyle end
struct MutableStr   <: MutableStyle end

MutableStyle(A::AbstractString) = MutableStyle(typeof(A))

MutableStyle(A::Str) = ImmutableStr()

_ismutable(::ImmutableStr, str::Type{T}) where {T<:Str} = false
_ismutable(::MutableStr, str::Type{T}) where {T<:Str} = true
is_mutable(::Type{T}) where {T<:Str} = _ismutable(MutableStyle(T))

isimmutable(str::T) where {T<:Str} = !is_mutable(T)

# Comparison traits

"""
    CompareStyle(Union{A, typeof(A)}, Union{B, typeof(B)})

`CompareStyle` specifies how to compare two strings with character set encodings A and B

    Strs.CompareStyle(::Type{<:MyString}, ::Type{String}) = ByteCompare()

The default is `CodePointCompare`
"""
abstract type CompareStyle end

# Todo: merge these with EqualsStyle

struct NoCompare        <: CompareStyle end # For equality checks, can't be equal
struct ByteCompare      <: CompareStyle end # Compare bytewise
struct ASCIICompare     <: CompareStyle end # Compare bytewise for ASCII subset, else codepoint
struct WordCompare      <: CompareStyle end # Compare first not equal word with <
struct UTF16Compare     <: CompareStyle end # Compare first not equal word adjusting > 0xd7ff
struct WidenCompare     <: CompareStyle end # Narrower can be simply widened for comparisons
struct CodePointCompare <: CompareStyle end # Compare CodePoints

CompareStyle(::Type{<:CSE}, ::Type{<:CSE}) = CodePointCompare()

CompareStyle(::Type{C}, ::Type{C}) where {C<:CSE} =
    ByteCompare()
CompareStyle(::Type{C}, ::Type{C}) where {C<:Union{Word_CSEs,Quad_CSEs}} =
    WordCompare()

CompareStyle(::Type{UTF16CSE},    ::Type{UTF16CSE})    = UTF16Compare()
CompareStyle(::Type{UTF16CSE},    ::Type{<:UCS2_CSEs}) = UTF16Compare()
CompareStyle(::Type{<:UCS2_CSEs}, ::Type{UTF16CSE})    = UTF16Compare()

CompareStyle(::Type{ASCIICSE}, ::Type{<:Union{Binary_CSEs,Latin_CSEs,UTF8_CSEs}}) =
    ByteCompare()
CompareStyle(::Type{ASCIICSE}, ::Type{<:Union{Word_CSEs,Quad_CSEs}}) =
    WidenCompare()

CompareStyle(::Type{<:Latin_CSEs}, ::Type{<:Latin_CSEs}) =
    ByteCompare()
CompareStyle(::Type{<:Latin_CSEs}, ::Type{UTF8_CSEs}) =
    ASCIICompare()
CompareStyle(::Type{<:Latin_CSEs}, ::Type{<:Union{Word_CSEs,Quad_CSEs}}) =
    WidenCompare()

CompareStyle(::Type{<:UCS2_CSEs}, ::Type{<:Union{ASCIICSE,Binary_CSEs,Latin_CSEs,Quad_CSEs}}) =
    WidenCompare()
CompareStyle(::Type{<:UCS2_CSEs},  ::Type{<:UCS2_CSEs}) =
    WordCompare()

CompareStyle(::Type{<:UTF32_CSEs},
             ::Type{<:Union{ASCIICSE,Binary_CSEs,Latin_CSEs,Text2CSE,UCS2_CSEs}}) =
    WidenCompare()
CompareStyle(::Type{<:UTF32_CSEs},  ::Type{<:UTF32_CSEs}) =
    WordCompare()

CompareStyle(::Type{S}, ::Type{T}) where {S<:AbstractString, T<:AbstractString} =
    CompareStyle(cse(S), cse(T))
CompareStyle(::Type{T}, ::Type{T}) where {T<:AbstractString} = ByteCompare()

CompareStyle(A::AbstractString, B::AbstractString) = CompareStyle(typeof(A), typeof(B))

Base.promote_rule(::Type{NoCompare},    ::Type{<:CompareStyle})   = NoCompare
Base.promote_rule(::Type{ByteCompare},  ::Type{CodePointCompare}) = ByteCompare
Base.promote_rule(::Type{WordCompare},  ::Type{CodePointCompare}) = WordCompare
Base.promote_rule(::Type{UTF16Compare}, ::Type{CodePointCompare}) = UTF16Compare
Base.promote_rule(::Type{ASCIICompare}, ::Type{CodePointCompare}) = ASCIICompare
Base.promote_rule(::Type{WidenCompare}, ::Type{CodePointCompare}) = WidenCompare

"""
    CanContain(Union{A, typeof(A)}, Union{B, typeof(B)})

`CanContainStyle` specifies whether the first string can contain a substring of the second type,
    and if so, what is the most efficient method of comparison

    Strs.CanContain(::Type{<:MyString}, ::Type{String}) = ByteCompare()

Returns an instance of type `CompareStyle`, default `CodePointCompare`
"""
CanContain(::Type{<:CSE}, ::Type{<:CSE}) = CodePointCompare()

CanContain(::Type{C}, ::Type{C}) where {C<:CSE} = ByteCompare()

CanContain(::Type{<:Binary_CSEs}, ::Type{<:Union{_UCS2CSE,_UTF32CSE}}) =
    NoCompare()
CanContain(::Type{<:Binary_CSEs}, ::Type{<:Byte_CSEs}) =
    ByteCompare()
CanContain(::Type{<:Binary_CSEs}, ::Type{<:Union{Word_CSEs, Quad_CSEs}}) =
    WidenCompare()

CanContain(::Type{ASCIICSE}, ::Type{<:SubSet_CSEs}) =
    NoCompare()
CanContain(::Type{ASCIICSE}, ::Type{<:Union{Binary_CSEs, LatinCSE, UTF8_CSEs}}) =
    ByteCompare()
CanContain(::Type{ASCIICSE}, ::Type{<:Union{Word_CSEs, Quad_CSEs}}) =
    WidenCompare()

CanContain(::Type{<:Latin_CSEs}, ::Type{<:Union{_UCS2CSE,_UTF32CSE}}) =
    NoCompare()
CanContain(::Type{<:Latin_CSEs}, ::Type{<:Union{Binary_CSEs,ASCIICSE,Latin_CSEs}}) =
    ByteCompare()
CanContain(::Type{<:Latin_CSEs}, ::Type{<:UTF8_CSEs}) =
    ASCIICompare()
CanContain(::Type{<:Latin_CSEs}, ::Type{<:Union{Word_CSEs, Quad_CSEs}}) =
    WidenCompare()

CanContain(::Type{<:UTF8_CSEs}, ::Type{<:Union{ASCIICSE,Binary_CSEs}}) =
    ByteCompare()
CanContain(::Type{<:UTF8_CSEs}, ::Type{<:Latin_CSEs}) =
    ASCIICompare()

CanContain(::Type{<:Union{Text2CSE,UCS2CSE}}, ::Type{_UTF32CSE}) =
    NoCompare()
CanContain(::Type{<:Union{Text2CSE,UCS2_CSEs}}, ::Type{<:Word_CSEs}) =
    ByteCompare()
CanContain(::Type{<:Union{Text2CSE,UCS2_CSEs}},
           ::Type{<:Union{ASCIICSE, Binary_CSEs, Latin_CSEs, Quad_CSEs}}) =
    WidenCompare()

CanContain(::Type{UTF16CSE}, ::Type{<:Union{ASCIICSE,Binary_CSEs,Latin_CSEs}}) =
    WidenCompare()
CanContain(::Type{UTF16CSE}, ::Type{<:Union{Text2CSE, UCS2_CSEs}}) =
    ByteCompare()

CanContain(::Type{<:Quad_CSEs}, ::Type{<:Quad_CSEs}) =
    ByteCompare()
CanContain(::Type{<:Quad_CSEs},
           ::Type{<:Union{Binary_CSEs,ASCIICSE,Latin_CSEs,Text2CSE,UCS2_CSEs}}) =
    WidenCompare()

CanContain(::Type{S}, ::Type{T}) where {S<:AbstractString, T<:AbstractString} =
    CanContain(cse(S), cse(T))
CanContain(::Type{T}, ::Type{T}) where {T<:AbstractString} = ByteCompare()

CanContain(A::AbstractString, B::AbstractString) = CanContain(typeof(A), typeof(B))

"""
    EqualsStyle(Union{A, typeof(A)}, Union{B, typeof(B)})

`EqualsStyle` specifies whether the first string can equal a substring of the second type,
    and if so, what is the most efficient method of comparison

This is determined by the Strs.CanContain trait
Returns an instance of type `CompareStyle`, default `CodePointCompare`
"""
EqualsStyle(::Type{S}, ::Type{T}) where {S<:CSE,T<:CSE} =
    Base.promote_typeof(CanContain(S, T), CanContain(T, S))()
EqualsStyle(::Type{T}, ::Type{T}) where {T<:CSE} = ByteCompare()

EqualsStyle(::Type{S}, ::Type{T}) where {S<:AbstractString, T<:AbstractString} =
    EqualsStyle(cse(S), cse(T))
EqualsStyle(::Type{T}, ::Type{T}) where {T<:AbstractString} = ByteCompare()

EqualsStyle(A::AbstractString, B::AbstractString) =  EqualsStyle(typeof(A), typeof(B))


text:
"""
*loglike*

Computes the log likelihood for the TAX model.
* `δ`: transfer of attention parameter
* `γ`: probability weighting parameter
* `β`: utility curvature parameter
* `θ`: decision consistency parameter
* `pa`: outcome probabilities for gamble a
* `va`: outcome values for gamble a
* `pb`: outcome probabilities for gamble b
* `vb`: outcome values for gamble b

Function Signature
```julia
loglike(δ, β, γ, θ, pa, va, pb, vb, data)
```
"""
function loglike(δ, β, γ, θ, pa, va, pb, vb, data)
    eua,eub = expected_utilities(δ, β, γ, θ, pa, va, pb, vb)
    p = choice_prob(eua, eub, θ)
    p = max(p, eps())
    p = min(1-eps(), p)
    return logpdf(Bernoulli(p), data)
end

function choice_prob(eua, eub, θ)
    return 1/(1 + exp(-θ*(eua - eub)))
end

function simulate(δ, β, γ, θ, pa, va, pb, vb)
    eua,eub = expected_utilities(δ, β, γ, θ, pa, va, pb, vb)
    p = choice_prob(eua, eub, θ)
    return rand() ≤ p ? true : false
end

function expected_utilities(δ, β, γ, θ, pa, va, pb, vb)
    model = TAX(;γ, δ, β)
    gambleA = Gamble(;p=pa, v=va)
    gambleB = Gamble(;p=pb, v=vb)
    eua = mean(model, gambleA)
    eub = mean(model, gambleB)
    return eua,eub
end

function random_design(outcome_dist, n_vals, n_choices)
    design = Vector{Vector{Float64}}()
    for _ in 1:n_choices
        p = rand(Dirichlet(n_vals, 1))
        push!(design, p)
        v = rand(outcome_dist, n_vals)
        push!(design, v)
    end
    return design
end

function abs_zscore(x)
    model = ExpectedUtility(1)
    gambleA = Gamble(;p=x[1], v=x[2])
    gambleB = Gamble(;p=x[3], v=x[4])
    μa = mean(model, gambleA)
    μb = mean(model, gambleB)
    σa = var(model, gambleA)
    σb = var(model, gambleB)
    return abs((μa - μb)/sqrt(σa + σb))
end

# using Revise, Turing, Random, MCMCChains

# # Set the true probability of heads in a coin.
# p_true = 0.5

# # Iterate from having seen 0 observations to 100 observations.
# Ns = 0:100

# # Draw data from a Bernoulli distribution, i.e. draw heads or tails.
# Random.seed!(12)
# data = rand(Bernoulli(p_true), last(Ns))

# # Declare our Turing model.
# @model function coinflip(y)
#     # Our prior belief about the probability of heads in a coin.
#     p ~ Beta(1, 1)

#     # The number of observations.
#     N = length(y)
#     for n in 1:N
#         # Heads or tails of a coin are drawn from a Bernoulli distribution.
#         y[n] ~ Bernoulli(p)
#     end
# end

# # Settings of the Hamiltonian Monte Carlo (HMC) sampler.
# iterations = 1000
# ϵ = 0.05
# τ = 10

# # Start sampling.
# chain = sample(coinflip(data), HMC(ϵ, τ), iterations)

# lpfun = function f(chain::Chains, dist, chains) # function to compute the logpdf values
#     niter, nparams, nchains = size(chain)
#     lp = zeros(niter + nchains) # resulting logpdf values
#     for i = 1:nparams
#         lp += map(p ->logpdf(dist(p...), data), Array(chain[:,i,:]))
#     end
#     return lp
# end

# lpfun = function f(chain::Chains) # function to compute the logpdf values
#     niter, _, nchains = size(chain)
#     parms = (Symbol.(chain.name_map.parameters)...,)
#     lp = zeros(niter*nchains, 1) # resulting logpdf values
#     for parm in parms
#         lp .+= map(p -> logpdf(Bernoulli(p...), data), Array(chain[parm]))
#     end
#     return lp
# end

# DIC, pD = dic(chain, lpfun)
#


text:
"""
    computeDelta!(p::QuadraticEigenvalueProblem{N,T}, c::AIMCache{N,Polynomial{T}}) where {N <: Unsigned, T <: Number}

Compute and return the AIM "quantization condition".

# Input
- `m::AIMSteppingMethod`: The stepping method to use.
- `p::QuadraticEigenvalueProblem`: A quadratic frequency problem.
- `c::AIMCache`: An AIM cache object created from p.

# Output
An object of type `Polynomial{T}` whose roots are the problem's eigenvalues.
"""
function computeDelta!(m::AIMSteppingMethod, p::QuadraticEigenvalueProblem{N,T}, c::AIMCache{N,Polynomial{T}}) where {N <: Unsigned, T <: Number}

    if (get_niter(p) + one(N)) != c.size
        error("The provided cache cannot hold data for $(get_niter(p)) AIM iterations.
        Make sure that the passed AIMCache corresponds to the passed AIM problem")
    end

    # Initialize the current data arrays with the initial data for the first step
    copy!(c.ccda, c.icda)
    copy!(c.cdda, c.idda)

    # Perform the aim steps
    for i in one(N):get_niter(p)
        AIMStep!(m, p, c)
    end

    # Compute and return the AIM "quantization condition"
    return c.cdda[1]*c.pcda[1] - c.pdda[1]*c.ccda[1]
end

"""
    computeDelta!(m::AIMSteppingMethod, p::NumericAIMProblem{N,T}, c::AIMCache{N,T}, ω::T) where {N <: Unsigned, T <: Number}

Compute and return the AIM "quantization condition".

# Input
- `m::AIMSteppingMethod`: The stepping method to use.
- `p::QuadraticEigenvalueProblem`: A quadratic frequency problem.
- `c::AIMCache`: An AIM cache object created from p.
- `ω::T`: Point to evaluate the quantization condition.

# Output
An object of type `T` which represents the AIM quantization condition at point ω.
"""
function computeDelta!(m::AIMSteppingMethod, p::NumericAIMProblem{N,T}, c::AIMCache{N,T}, ω::T) where {N <: Unsigned, T <: Number}

    if (get_niter(p) + one(N)) != c.size
        error("The provided cache cannot hold data for $(get_niter(p)) AIM iterations.
        Make sure that the passed AIMCache corresponds to the passed AIM problem")
    end

    # Fill iitial arrays with data corresponding to the point ω
    recomputeInitials!(p,c,ω)

    # Initialize the current data arryas with the initial data for the first step
    copy!(c.ccda, c.icda)
    copy!(c.cdda, c.idda)

    # Perform the aim steps
    for i in one(N):get_niter(p)
        AIMStep!(m, p, c)
    end

    # Compute and return the AIM "quantization condition" at point ω
    return c.cdda[1]*c.pcda[1] - c.pdda[1]*c.ccda[1]
end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

include("model_AdditionalUnattendContent.jl")
include("model_ApiEntityReference.jl")
include("model_ApiError.jl")
include("model_ApiErrorBase.jl")
include("model_AvailabilitySetListResult.jl")
include("model_AvailabilitySetProperties.jl")
include("model_BootDiagnostics.jl")
include("model_BootDiagnosticsInstanceView.jl")
include("model_Caching.jl")
include("model_ComputeLongRunningOperationProperties.jl")
include("model_CreateOption.jl")
include("model_DataDisk.jl")
include("model_DataDiskImage.jl")
include("model_DiagnosticsProfile.jl")
include("model_DiskEncryptionSettings.jl")
include("model_DiskInstanceView.jl")
include("model_HardwareProfile.jl")
include("model_ImageDataDisk.jl")
include("model_ImageListResult.jl")
include("model_ImageOSDisk.jl")
include("model_ImageProperties.jl")
include("model_ImageStorageProfile.jl")
include("model_InnerError.jl")
include("model_InstanceViewStatus.jl")
include("model_KeyVaultKeyReference.jl")
include("model_KeyVaultSecretReference.jl")
include("model_LinuxConfiguration.jl")
include("model_ListUsagesResult.jl")
include("model_MaintenanceRedeployStatus.jl")
include("model_NetworkInterfaceReferenceProperties.jl")
include("model_NetworkProfile.jl")
include("model_OSDisk.jl")
include("model_OSDiskImage.jl")
include("model_OSProfile.jl")
include("model_OperationStatusResponse.jl")
include("model_Plan.jl")
include("model_PurchasePlan.jl")
include("model_Resource.jl")
include("model_ResourceSku.jl")
include("model_ResourceSkuCapabilities.jl")
include("model_ResourceSkuCapacity.jl")
include("model_ResourceSkuCosts.jl")
include("model_ResourceSkuRestrictions.jl")
include("model_ResourceSkusResult.jl")
include("model_RollingUpgradePolicy.jl")
include("model_RollingUpgradeProgressInfo.jl")
include("model_RollingUpgradeRunningStatus.jl")
include("model_RollingUpgradeStatusInfoProperties.jl")
include("model_Sku.jl")
include("model_SshConfiguration.jl")
include("model_SshPublicKey.jl")
include("model_StorageAccountType.jl")
include("model_StorageProfile.jl")
include("model_SubResource.jl")
include("model_SubResourceReadOnly.jl")
include("model_UpdateResource.jl")
include("model_UpgradePolicy.jl")
include("model_Usage.jl")
include("model_UsageName.jl")
include("model_VaultCertificate.jl")
include("model_VaultSecretGroup.jl")
include("model_VirtualHardDisk.jl")
include("model_VirtualMachineAgentInstanceView.jl")
include("model_VirtualMachineCaptureParameters.jl")
include("model_VirtualMachineCaptureResultProperties.jl")
include("model_VirtualMachineExtensionHandlerInstanceView.jl")
include("model_VirtualMachineExtensionImageProperties.jl")
include("model_VirtualMachineExtensionInstanceView.jl")
include("model_VirtualMachineExtensionProperties.jl")
include("model_VirtualMachineHealthStatus.jl")
include("model_VirtualMachineIdentity.jl")
include("model_VirtualMachineImageProperties.jl")
include("model_VirtualMachineInstanceView.jl")
include("model_VirtualMachineListResult.jl")
include("model_VirtualMachineProperties.jl")
include("model_VirtualMachineScaleSetDataDisk.jl")
include("model_VirtualMachineScaleSetExtensionListResult.jl")
include("model_VirtualMachineScaleSetExtensionProfile.jl")
include("model_VirtualMachineScaleSetExtensionProperties.jl")
include("model_VirtualMachineScaleSetIPConfigurationProperties.jl")
include("model_VirtualMachineScaleSetIdentity.jl")
include("model_VirtualMachineScaleSetInstanceView.jl")
include("model_VirtualMachineScaleSetInstanceViewStatusesSummary.jl")
include("model_VirtualMachineScaleSetListResult.jl")
include("model_VirtualMachineScaleSetListSkusResult.jl")
include("model_VirtualMachineScaleSetListWithLinkResult.jl")
include("model_VirtualMachineScaleSetManagedDiskParameters.jl")
include("model_VirtualMachineScaleSetNetworkConfigurationDnsSettings.jl")
include("model_VirtualMachineScaleSetNetworkConfigurationProperties.jl")
include("model_VirtualMachineScaleSetNetworkProfile.jl")
include("model_VirtualMachineScaleSetOSDisk.jl")
include("model_VirtualMachineScaleSetOSProfile.jl")
include("model_VirtualMachineScaleSetProperties.jl")
include("model_VirtualMachineScaleSetPublicIPAddressConfiguration.jl")
include("model_VirtualMachineScaleSetPublicIPAddressConfigurationDnsSettings.jl")
include("model_VirtualMachineScaleSetPublicIPAddressConfigurationProperties.jl")
include("model_VirtualMachineScaleSetSku.jl")
include("model_VirtualMachineScaleSetSkuCapacity.jl")
include("model_VirtualMachineScaleSetStorageProfile.jl")
include("model_VirtualMachineScaleSetUpdateIPConfigurationProperties.jl")
include("model_VirtualMachineScaleSetUpdateNetworkConfigurationProperties.jl")
include("model_VirtualMachineScaleSetUpdateNetworkProfile.jl")
include("model_VirtualMachineScaleSetUpdateOSDisk.jl")
include("model_VirtualMachineScaleSetUpdateOSProfile.jl")
include("model_VirtualMachineScaleSetUpdateProperties.jl")
include("model_VirtualMachineScaleSetUpdatePublicIPAddressConfiguration.jl")
include("model_VirtualMachineScaleSetUpdatePublicIPAddressConfigurationProperties.jl")
include("model_VirtualMachineScaleSetUpdateStorageProfile.jl")
include("model_VirtualMachineScaleSetUpdateVMProfile.jl")
include("model_VirtualMachineScaleSetVMExtensionsSummary.jl")
include("model_VirtualMachineScaleSetVMInstanceIDs.jl")
include("model_VirtualMachineScaleSetVMInstanceRequiredIDs.jl")
include("model_VirtualMachineScaleSetVMInstanceView.jl")
include("model_VirtualMachineScaleSetVMListResult.jl")
include("model_VirtualMachineScaleSetVMProfile.jl")
include("model_VirtualMachineScaleSetVMProperties.jl")
include("model_VirtualMachineSize.jl")
include("model_VirtualMachineSizeListResult.jl")
include("model_VirtualMachineStatusCodeCount.jl")
include("model_WinRMConfiguration.jl")
include("model_WinRMListener.jl")
include("model_WindowsConfiguration.jl")
include("model_AvailabilitySet.jl")
include("model_Image.jl")
include("model_ImageReference.jl")
include("model_ManagedDiskParameters.jl")
include("model_NetworkInterfaceReference.jl")
include("model_RollingUpgradeStatusInfo.jl")
include("model_VirtualMachine.jl")
include("model_VirtualMachineCaptureResult.jl")
include("model_VirtualMachineExtension.jl")
include("model_VirtualMachineExtensionImage.jl")
include("model_VirtualMachineImageResource.jl")
include("model_VirtualMachineScaleSet.jl")
include("model_VirtualMachineScaleSetExtension.jl")
include("model_VirtualMachineScaleSetIPConfiguration.jl")
include("model_VirtualMachineScaleSetNetworkConfiguration.jl")
include("model_VirtualMachineScaleSetUpdate.jl")
include("model_VirtualMachineScaleSetUpdateIPConfiguration.jl")
include("model_VirtualMachineScaleSetUpdateNetworkConfiguration.jl")
include("model_VirtualMachineScaleSetVM.jl")
include("model_VirtualMachineImage.jl")


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

using Markdown
using InteractiveUtils

# ╔═╡ c930a0fa-4ef8-4963-b4b2-0d576ecdde7e
using Pkg, DrWatson

# ╔═╡ 919b41e5-7e43-45e0-8676-373b77ac6208
begin
	@quickactivate "NumericalMethodsforEngineers"
	using Roots, Plots
	using Symbolics, SymbolicUtils, Latexify, LaTeXStrings, ModelingToolkit
	using NumericalMethodsforEngineers, Test
end

# ╔═╡ c34bdc36-5013-4096-9cb2-581ad9e02239
md" ## f-03-01.jl"

# ╔═╡ 0b1514c5-eded-4f60-ac31-6f27a4deff1d
begin
	@variables f x
	f = exp(x) - x^4
	L(f)
end

# ╔═╡ 7c4288f0-38d9-4c39-b6b1-789aba5465b6
f_j = build_function(f, [x], expression=Val{false});

# ╔═╡ c7c87623-9c41-4a68-8673-47a71c67f083
begin
	y = -10.0:0.1:10.0
	plot(y, f_j.(y), leg=false)
	hline!([0.0])
	vline!([fzeros(f_j, [-1, 2])])
end

# ╔═╡ 2708dfb1-09dd-458c-8b5f-2e2ec3d800c4
begin
	D = Differential(x)
	fp_s = expand_derivatives(D(f))
	fp_j = build_function(fp_s, [x], expression=Val{false})
	L(fp_s)
end

# ╔═╡ 1dd78a44-b910-4bb2-8af0-9a1ab3ee4a51
Roots.newton(f_j, fp_j, 8)        # 8.6131694564414

# ╔═╡ 8054c3fd-3ecb-4c2a-b606-7fd518d750a3
begin
	fpp_s = expand_derivatives(D(D(f)))
	fpp_j = build_function(fpp_s, [x], expression=Val{false})
	L(fpp_s)
end

# ╔═╡ 7e2f4bd6-4f9a-4a72-8e70-e5fc56f806ee
halley(f_j, fp_j, fpp_j, 8)

# ╔═╡ 3e022154-6933-4ba7-8cf9-2d48d3cbb27f
secant_method(f_j, [8, 8.5])

# ╔═╡ 26cd9127-ae9e-4f09-97e4-d10f514e2164
md" ## End of f-03-01.jl"

# ╔═╡ Cell order:
# ╟─c34bdc36-5013-4096-9cb2-581ad9e02239
# ╠═c930a0fa-4ef8-4963-b4b2-0d576ecdde7e
# ╠═919b41e5-7e43-45e0-8676-373b77ac6208
# ╠═0b1514c5-eded-4f60-ac31-6f27a4deff1d
# ╠═7c4288f0-38d9-4c39-b6b1-789aba5465b6
# ╠═c7c87623-9c41-4a68-8673-47a71c67f083
# ╠═2708dfb1-09dd-458c-8b5f-2e2ec3d800c4
# ╠═1dd78a44-b910-4bb2-8af0-9a1ab3ee4a51
# ╠═8054c3fd-3ecb-4c2a-b606-7fd518d750a3
# ╠═7e2f4bd6-4f9a-4a72-8e70-e5fc56f806ee
# ╠═3e022154-6933-4ba7-8cf9-2d48d3cbb27f
# ╟─26cd9127-ae9e-4f09-97e4-d10f514e2164


text:
@testset "High-level wrapper" begin
    @testset "Generated structs" begin
        test(StructDefinition{true}, struct_by_name, :VkPhysicalDeviceProperties, :(
            @auto_hash_equals struct PhysicalDeviceProperties <: HighLevelStruct
                api_version::VersionNumber
                driver_version::VersionNumber
                vendor_id::UInt32
                device_id::UInt32
                device_type::PhysicalDeviceType
                device_name::String
                pipeline_cache_uuid::NTuple{Int(VK_UUID_SIZE), UInt8}
                limits::PhysicalDeviceLimits
                sparse_properties::PhysicalDeviceSparseProperties
            end
        ))

        test(StructDefinition{true}, struct_by_name, :VkExternalBufferProperties, :(
            @auto_hash_equals struct ExternalBufferProperties <: HighLevelStruct
                next::Any
                external_memory_properties::ExternalMemoryProperties
            end
        ))

        test(StructDefinition{true}, struct_by_name, :VkPipelineExecutableInternalRepresentationKHR, :(
            @auto_hash_equals struct PipelineExecutableInternalRepresentationKHR <: HighLevelStruct
                next::Any
                name::String
                description::String
                is_text::Bool
                data_size::UInt
                data::OptionalPtr{Ptr{Cvoid}}
            end
        ))

        test(StructDefinition{true}, struct_by_name, :VkApplicationInfo, :(
            @auto_hash_equals struct ApplicationInfo <: HighLevelStruct
                next::Any
                application_name::String
                application_version::VersionNumber
                engine_name::String
                engine_version::VersionNumber
                api_version::VersionNumber
            end
        ))

        test(StructDefinition{true}, struct_by_name, :VkInstanceCreateInfo, :(
            @auto_hash_equals struct InstanceCreateInfo <: HighLevelStruct
                next::Any
                flags::UInt32
                application_info::OptionalPtr{ApplicationInfo}
                enabled_layer_names::Vector{String}
                enabled_extension_names::Vector{String}
            end
        ))

        test(StructDefinition{true}, struct_by_name, :VkXcbSurfaceCreateInfoKHR, :(
            @auto_hash_equals struct XcbSurfaceCreateInfoKHR <: HighLevelStruct
                next::Any
                flags::UInt32
                connection::Ptr{vk.xcb_connection_t}
                window::vk.xcb_window_t
            end
        ))
    end

    @testset "Conversion to low-level structs" begin
        test_ex(
            Constructor(
                StructDefinition{false}(struct_by_name(:VkApplicationInfo)),
                StructDefinition{true}(struct_by_name(:VkApplicationInfo)),
            ),
            :(
                _ApplicationInfo(x::ApplicationInfo) = _ApplicationInfo(x.application_version, x.engine_version, x.api_version; x.next, x.application_name, x.engine_name)
            ),
        )
        test_ex(
            Constructor(
                StructDefinition{false}(struct_by_name(:VkInstanceCreateInfo)),
                StructDefinition{true}(struct_by_name(:VkInstanceCreateInfo)),
            ),
            :(
                _InstanceCreateInfo(x::InstanceCreateInfo) = _InstanceCreateInfo(x.enabled_layer_names, x.enabled_extension_names; x.next, x.flags, application_info = convert_nonnull(_ApplicationInfo, x.application_info))
            ),
        )
        test_ex(
            Constructor(
                StructDefinition{false}(struct_by_name(:VkRenderingAttachmentInfo)),
                StructDefinition{true}(struct_by_name(:VkRenderingAttachmentInfo)),
            ),
            :(
                _RenderingAttachmentInfo(x::RenderingAttachmentInfo) = _RenderingAttachmentInfo(x.image_layout, x.resolve_image_layout, x.load_op, x.store_op, convert_nonnull(_ClearValue, x.clear_value); x.next, x.image_view, x.resolve_mode, x.resolve_image_view)
            ),
        )
    end

    @testset "Additional constructor" begin
        test(Constructor, StructDefinition{true}(StructDefinition{false}(struct_by_name(:VkApplicationInfo))), :(
            ApplicationInfo(application_version::VersionNumber, engine_version::VersionNumber, api_version::VersionNumber; next = C_NULL, application_name = "", engine_name = "") = ApplicationInfo(next, application_name, application_version, engine_name, engine_version, api_version)
        ))
    end

    @testset "`Core -> High-level constructors`" begin
        function test_constructor_core_to_hl(name, with_next_types, body)
            spec = struct_by_name(name)
            def = StructDefinition{true}(struct_by_name(name))
            cons = Constructor(def, spec)
            expected = :($(VulkanGen.struct_name(name, true))(x::$name) = $body)
            with_next_types && push!(expected.args[1].args, :(next_types::Type...))
            test_ex(cons, expected)
        end

        test_constructor_core_to_hl(:VkLayerProperties, false, :(
            LayerProperties(from_vk(String, x.layerName), from_vk(VersionNumber, x.specVersion), from_vk(VersionNumber, x.implementationVersion), from_vk(String, x.description))
        ))

        test_constructor_core_to_hl(:VkQueueFamilyProperties, false, :(
            QueueFamilyProperties(x.queueFlags, x.queueCount, x.timestampValidBits, Extent3D(x.minImageTransferGranularity))
        ))

        test_constructor_core_to_hl(:VkPhysicalDeviceMemoryProperties, false, :(
            PhysicalDeviceMemoryProperties(x.memoryTypeCount, MemoryType.(x.memoryTypes), x.memoryHeapCount, MemoryHeap.(x.memoryHeaps))
        ))

        test_constructor_core_to_hl(:VkDisplayPlaneCapabilities2KHR, true, :(
            DisplayPlaneCapabilities2KHR(load_next_chain(x.pNext, next_types...), DisplayPlaneCapabilitiesKHR(x.capabilities))
        ))

        test_constructor_core_to_hl(:VkDrmFormatModifierPropertiesListEXT, true, :(
            DrmFormatModifierPropertiesListEXT(load_next_chain(x.pNext, next_types...), unsafe_wrap(Vector{DrmFormatModifierPropertiesEXT}, x.pDrmFormatModifierProperties, x.drmFormatModifierCount; own = true))
        ))

        test_constructor_core_to_hl(:VkPhysicalDeviceGroupProperties, true, :(
            PhysicalDeviceGroupProperties(load_next_chain(x.pNext, next_types...), x.physicalDeviceCount, PhysicalDevice.(x.physicalDevices), from_vk(Bool, x.subsetAllocation))
        ))

        test_constructor_core_to_hl(:VkAccelerationStructureVersionInfoKHR, true, :(
            AccelerationStructureVersionInfoKHR(load_next_chain(x.pNext, next_types...), unsafe_wrap(Vector{UInt8}, x.pVersionData, 2VK_UUID_SIZE; own = true))
        ))
    end

    @testset "`Low-level -> High-level constructors`" begin
        function test_constructor_ll_to_hl(name, ex)
            s = struct_by_name(name)
            def_hl = StructDefinition{true}(s)
            def_ll = StructDefinition{false}(s)
            c = Constructor(def_hl, def_ll)
            test_ex(c, ex)
        end
        test_constructor_ll_to_hl(:VkPhysicalDeviceFeatures2, :(
            function PhysicalDeviceFeatures2(x::_PhysicalDeviceFeatures2, next_types::Type...)
                (; deps) = x
                GC.@preserve deps PhysicalDeviceFeatures2(x.vks, next_types...)
            end
        ))
        test_constructor_ll_to_hl(:VkPhysicalDeviceFeatures, :(
            PhysicalDeviceFeatures(x::_PhysicalDeviceFeatures) = PhysicalDeviceFeatures(x.vks)
        ))
    end
end

max_stars_repo_name: JuliaGPU/VulkanAbstraction

max_issues_repo_name: JuliaGPU/VulkanAbstraction

max_forks_repo_name: JuliaGPU/VulkanAbstraction


text:
# === Design Notes ===
#
# `DataValueArray{T, N}` is a struct-of-arrays representation of
# `Array{DataValue{T}, N}`. This makes it easy to define complicated operations
# (e.g. matrix multiplication) by reusing the existing definition for
# `Array{T}`.
#
# One complication when defining functions that operate on the internal fields
# of a `DataValueArray` is that developers must take care to ensure that they
# do not index into an undefined entry in the `values` field. This is not a
# problem for `isbits` types, which are never `#undef`, but will trigger
# an exception for any other type.
#
# TODO: Ensure that size(values) == size(isnull) using inner constructor.
"""
`DataValueArray{T, N}` is an efficient alternative to `Array{DataValue{T}, N}`.
It allows users to easily define operations on arrays with null values by
reusing operations that only work on arrays without any null values.
"""
immutable DataValueArray{T, N} <: AbstractArray{DataValue{T}, N}
    values::Array{T, N}
    isnull::Array{Bool, N}

    function DataValueArray{T, N}(d::AbstractArray{T, N}, m::AbstractArray{Bool, N}) where {T, N}
        if size(d) != size(m)
            msg = "values and missingness arrays must be the same size"
            throw(ArgumentError(msg))
        end
        new(d, m)
    end
end

DataValueVector{T} = DataValueArray{T, 1}
DataValueMatrix{T} = DataValueArray{T, 2}


text:
#===============================================================================

                        Definitions of Workspaces
        Defines *structs*:
        - StandardGlobalSubworkspace - convenience container that is most likely
                             going to be used even with custom global Workspaces
        - StandardLocalSubworkspace - convenience container that is most likely
                             going to be used even with custom local Workspaces
        - GenericGlobalWorkspace - simple global Workspace for simple problems
        - GenericLocalWorkspace - simple local Workspace for simple problems
        As well as functions:
        - init_global_workspace
        - create_workspaces

===============================================================================#

#                              GLOBAL WORKSPACE
#-------------------------------------------------------------------------------

#=        |-------------------------------------------------------|
          |    each custom GlobalWorkspace MUST overwrite these:  |
          |-------------------------------------------------------|           =#
"""
    init_global_workspace(
        ::MCMCBackend,
        num_mcmc_steps,
        updates::Vector{<:MCMCUpdate},
        data,
        θinit::Vector{T};
        kwargs...
    ) where T

Initialize the `<custom>GlobalWorkspace`. `<custom>MCMCBackend` points to which
`GlobalWorkspace` constructors to use, `updates` is a list of MCMC updates,
`θinit` is the initial guess for the parameter and `kwargs` are the named
arguments passed to `run!`.
"""
function init_global_workspace(
        ::MCMCBackend,
        num_mcmc_steps,
        updates::Vector{<:MCMCUpdate},
        data,
        θinit::Vector{T};
        kwargs...
    ) where T
    error("init_global_workspace not implemented")
end

"""
    loglikelihood(ws::GlobalWorkspace, ::Proposal)

Evaluate the loglikelihood for the proposal Law and observations stored in a
global workspace.
"""
function Distributions.loglikelihood(ws::GlobalWorkspace, ::Proposal)
    error(
        string(
            "loglikelihood(ws::GlobalWorkspace, ::Proposal) for a custom ",
            "GlobalWorkspace is not implemented"
        )
    )
end

"""
    loglikelihood(ws::GlobalWorkspace, ::Previous)

Evaluate the loglikelihood for the accepted Law and observations stored in a
global workspace.
"""
function Distributions.loglikelihood(ws::GlobalWorkspace, ::Previous)
    error(
        string(
            "loglikelihood(ws::GlobalWorkspace, ::Previous) for a custom ",
            "GlobalWorkspace is not implemented"
        )
    )
end

#=
        for custom GlobalWorkspaces these methods will work so long as
        they contain sub_ws::StandardGlobalSubworkspace, otherwise they
        need to be overwritten.
                                                                              =#

"""
    num_mcmc_steps(ws::GlobalWorkspace)

Return the total set number of MCMC iterations.
"""
num_mcmc_steps(ws::GlobalWorkspace) = num_mcmc_steps(ws.sub_ws)

"""
    state(ws::GlobalWorkspace)

Return currently accepted state of the chain.
"""
state(ws::GlobalWorkspace) = state(ws.sub_ws)

"""
    state(ws::GlobalWorkspace, step)

Return the state of the chain accepted at the the `step` iteration of the Markov
chain.
"""
state(ws::GlobalWorkspace, step) = state(ws.sub_ws, step)

"""
    state°(ws::GlobalWorkspace, step)

Return the state of the chain proposed at the the `step` iteration of the Markov
chain.
"""
state°(ws::GlobalWorkspace, step) = state°(ws.sub_ws, step)

"""
    num_updt(ws::GlobalWorkspace)

Return the total set number of MCMC updates that may be performed at each MCMC
iteration.
"""
num_updt(ws::GlobalWorkspace) = num_updt(ws.sub_ws)

"""
    estim_mean(ws::GlobalWorkspace)

Return the empirical mean of the parameter.
"""
estim_mean(ws::GlobalWorkspace) = estim_mean(ws.sub_ws)

"""
    estim_cov(ws::GlobalWorkspace)

Return the empirical covariance of the parameter.
"""
estim_cov(ws::GlobalWorkspace) = estim_cov(ws.sub_ws)

#=
                            standard sub-Workspace
                                                                              =#

"""
    struct StandardGlobalSubworkspace{T,TD} <: GlobalWorkspace{T}
        state::Vector{T}
        state_history::Vector{Vector{Vector{T}}}
        state_proposal_history::Vector{Vector{Vector{T}}}
        data::TD
        stats::GenericChainStats{T}
    end

Standard containers expected to be present in every global workspace. `state` is
the currently accepted parameter θ, `state_history` is a chain of `state`s that
have been accepted and `state_proposal_history` is a chain of `state`s that have
been proposed. `data` are the data passed to an MCMC sampler (usually just a
pointer) and `stats` gathers some basic online information about the chain.
"""
struct StandardGlobalSubworkspace{T,TD} <: GlobalWorkspace{T}
    state::Vector{T}
    state_history::Vector{Vector{Vector{T}}}
    state_proposal_history::Vector{Vector{Vector{T}}}
    data::TD
    stats::GenericChainStats{T}

    function StandardGlobalSubworkspace(
            num_mcmc_steps,
            num_updates,
            data::TD,
            θinit::Vector{T},
        ) where {T,TD}
        M, NU = num_mcmc_steps, num_updates # name alias
        ∅s = [Vector{T}(undef, length(θinit)) for _=1:NU]
        new{T,TD}(
            deepcopy(θinit),
            [deepcopy(∅s) for _=1:M],
            [deepcopy(∅s) for _=1:M],
            data,
            GenericChainStats(θinit, NU, M),
        )
    end
end

num_mcmc_steps(ws::StandardGlobalSubworkspace) = length(ws.state_history)
state(ws::StandardGlobalSubworkspace) = ws.state
function state(ws::StandardGlobalSubworkspace, step)
    ws.state_history[step.mcmciter][step.pidx]
end
function state°(ws::StandardGlobalSubworkspace, step)
    ws.state_proposal_history[step.mcmciter][step.pidx]
end
num_updt(ws::StandardGlobalSubworkspace) = length(first(ws.state_history))
estim_mean(ws::StandardGlobalSubworkspace) = ws.stats.mean
estim_cov(ws::StandardGlobalSubworkspace) = ws.stats.cov

#=
                        simple, custom GlobalWorkspace
                                                                              =#

"""
    struct GenericGlobalWorkspace{T,TD,TL} <: GlobalWorkspace{T}
        sub_ws::StandardGlobalSubworkspace{T,TD}
        P::TL
        P°::TL
    end

Generic global workspace with `sub_ws` containing current `state` and keeping
track of its history and some basic statistics. `P` and `P°` are the target laws
with accepted and proposal `state` set as parameters.
"""
struct GenericGlobalWorkspace{T,TD,TL} <: GlobalWorkspace{T}
    sub_ws::StandardGlobalSubworkspace{T,TD}
    P::TL
    P°::TL
end

function init_global_workspace(
        ::GenericMCMCBackend,
        num_mcmc_steps,
        updates::Vector{<:MCMCUpdate},
        data,
        θinit::Vector{T};
        kwargs...
    ) where T
    sub_ws = StandardGlobalSubworkspace(
        num_mcmc_steps,
        length(updates),
        data,
        θinit,
    )
    GenericGlobalWorkspace{T,typeof(data),typeof(data.P)}(
        sub_ws,
        deepcopy(data.P),
        deepcopy(data.P),
    )
end

function Distributions.loglikelihood(ws::GenericGlobalWorkspace, ::Proposal)
    loglikelihood(ws.P°, ws.sub_ws.data.obs)
end

#=
                generic methods for all GlobalWorkspaces
                                                                              =#

function Base.summary(ws::GlobalWorkspace; init=false)
    summary(Base.stdout, ws, Val(init))
end

Base.summary(io::IO, ws::GlobalWorkspace; init=false) = summary(io, ws, Val(init))

function Base.summary(io::IO, ws::GlobalWorkspace, ::Val{true})
    println(io, "Number of MCMC iterations: ", num_mcmc_steps(ws))
    println(io, "Initial θ guess: ", state(ws))
    println(io, "Number of updates at each MCMC iteration: ", num_updt(ws))
    println(io, "Chosen `GlobalWorkspace:` ", remove_curly(typeof(ws)))
end

function Base.summary(io::IO, ws::GlobalWorkspace, ::Val{false})
    println(io, "Number of MCMC iterations: ", num_mcmc_steps(ws))
    println(io, "Estimated E[θ]: ", estim_mean(ws))
    println(io, "Estimated Cov[θ]: ")
    show(io, "text/plain", estim_cov(ws))
end

#                              LOCAL WORKSPACE
#-------------------------------------------------------------------------------

#=
        for custom LocalWorkspaces these methods MUST be overwritten:
                                                                              =#
"""
    create_workspace(
        ::MCMCBackend,
        mcmcupdate,
        global_ws::GlobalWorkspace,
        num_mcmc_steps
    ) where {T}

Create a local workspace for a given `mcmcupdate`.
"""
function create_workspace(
        ::MCMCBackend,
        mcmcupdate,
        global_ws::GlobalWorkspace,
        num_mcmc_steps
    ) where {T}
    error("create_workspace (localworkspace) is not defined.")
end

"""
    accepted(ws::LocalWorkspace, i::Int)

Return boolean for whether the ith update has been accepted.
"""
function accepted(ws::LocalWorkspace, i::Int)
    error("accepted(ws::LocalWorkspace, i) is not implemented")
end

"""
    set_accepted!(ws::LocalWorkspace, i::Int, v)

Set boolean for whether the ith update has been accepted.
"""
set_accepted!(ws::LocalWorkspace, i::Int, v) = (ws.acceptance_history[i] = v)

#=
        for custom LocalWorkspaces these methods will work so long as
        they contain sub_ws::StandardGlobalSubworkspace and
        sub_ws°::StandardGlobalSubworkspace, otherwise they need to be
        overwritten.
                                                                              =#
"""
    ll(ws::LocalWorkspace)

Return log-likelihood of the currently accepted parameter.
"""
ll(ws::LocalWorkspace) = ll(ws.sub_ws)

"""
    ll°(ws::LocalWorkspace)

Return log-likelihood of the currently proposed parameter.
"""
ll°(ws::LocalWorkspace) = ll(ws.sub_ws°)

"""
    ll(ws::LocalWorkspace, i::Int)

Return log-likelihood of the ith accepted parameter.
"""
ll(ws::LocalWorkspace, i::Int) = ll(ws.sub_ws, i)

"""
    ll°(ws::LocalWorkspace, i::Int)

Return log-likelihood of the ith proposed parameter.
"""
ll°(ws::LocalWorkspace, i::Int) = ll(ws.sub_ws°)

"""
    state(ws::LocalWorkspace)

Return the last accepted state.
"""
state(ws::LocalWorkspace) = state(ws.sub_ws)

"""
    state°(ws::LocalWorkspace)

Return the last proposed state.
"""
state°(ws::LocalWorkspace) = state(ws.sub_ws°)

#=
        for custom LocalWorkspaces these are well defined, but may
        not make sense:
                                                                              =#
"""
    create_workspaces(v::MCMCBackend, mcmc::MCMC)

Create local workspaces, one for each update.
"""
function create_workspaces(backend::MCMCBackend, mcmc::MCMC)
    [
        create_workspace(
            backend,
            mcmc.updates[i],
            mcmc.workspace,
            mcmc.schedule.num_mcmc_steps
        ) for i in 1:mcmc.schedule.num_updates
    ]
end

"""
    llr(ws::LocalWorkspace, i::Int)

Compute log-likelihood ratio at ith mcmc iteration.
"""
llr(ws::LocalWorkspace, i::Int) = sum(ll°(ws, i) .- ll(ws, i))

"""
    name_of_update(ws::LocalWorkspace)

Return the name of the update.
"""
name_of_update(ws::LocalWorkspace) = "unknown name"

#=
                        standard local sub-Workspace
                                                                              =#

"""
    struct StandardLocalSubworkspace{T} <: LocalWorkspace{T}
        state::Vector{T}
        ll::Vector{Float64}
        ll_history::Vector{Vector{Float64}}
        ∇ll::Vector{Vector{Float64}}
        momentum::Vector{Float64}
    end

Standard containers likely to be present in every local workspace. `state` is
the currently accepted *subset* of parameter θ that the corresponding update
operates on, `ll` is the corresponding log-likelihood and `ll_history` is the
chain of log-likelihoods. A single `ll` is of the type `Vector{Float64}` to
reflect the fact that a problem might admit a natural factorisation into
independent components that may be operated on independently, in parallel with
each entry in `ll` corresponding to a separate component. `∇ll` is the gradient
of log-likelihood (needed by gradient-based algorithms) and `momentum` is the
variable needed for the Hamiltionan dynamics. `∇ll` and `momentum` may be simply
left untouched if the problem does not need them.
"""
struct StandardLocalSubworkspace{T} <: LocalWorkspace{T}
    state::Vector{T}
    ll::Vector{Float64}
    ll_history::Vector{Vector{Float64}}
    ∇ll::Vector{Vector{Float64}}
    momentum::Vector{Float64}

    function StandardLocalSubworkspace(
            state::Vector{T}, M; num_indep_blocks=1
        ) where T
        new{T}(
            deepcopy(state),
            Float64[-Inf for _=1:num_indep_blocks],
            [zeros(Float64, num_indep_blocks) for _=1:M],
            [zeros(Float64, length(state)) for _=1:num_indep_blocks],
            zero(state),
        )
    end
end

ll(ws::StandardLocalSubworkspace) = ws.ll
ll(ws::StandardLocalSubworkspace, i::Int) = ws.ll_history[i]
state(ws::StandardLocalSubworkspace) = ws.state

#=
                        simple, custom LocalWorkspace
                                                                              =#

"""
    struct GenericLocalWorkspace{T} <: LocalWorkspace{T}
        sub_ws::StandardLocalSubworkspace{T}
        sub_ws°::StandardLocalSubworkspace{T}
        acceptance_history::Vector{Bool}
    end

Generic local workspace with `sub_ws` containing a subset of `state` that the
corresponding updates operates on, with currently accepted value of `state` as
well as its log-likelihood. `sub_ws°` corresponds to a proposal `state`.
`acceptance_history` keeps track of accept/reject decisions.
"""
struct GenericLocalWorkspace{T} <: LocalWorkspace{T}
    sub_ws::StandardLocalSubworkspace{T}
    sub_ws°::StandardLocalSubworkspace{T}
    acceptance_history::Vector{Bool}
    updt_name::String
end

function create_workspace(
        ::GenericMCMCBackend,
        mcmcupdate,
        global_ws::GenericGlobalWorkspace{T},
        num_mcmc_steps
    ) where {T}
    _state = state(global_ws, mcmcupdate)[:]
    GenericLocalWorkspace{T}(
        StandardLocalSubworkspace(_state, num_mcmc_steps),
        StandardLocalSubworkspace(_state, num_mcmc_steps),
        Vector{Bool}(undef, num_mcmc_steps),
        string(remove_curly(typeof(mcmcupdate))),
    )
end

accepted(ws::GenericLocalWorkspace, i::Int) = ws.acceptance_history[i]
name_of_update(ws::GenericLocalWorkspace) = ws.updt_name


text:
using ControlSystems, RobustAndOptimalControl, MonteCarloMeasurements
using RobustAndOptimalControl: bisect_a

# Example from the diskmargin paper

L = tf(25, [1,10,10,10])
dm = diskmargin(L, 0)
@test dm.ω0 ≈ 1.94   atol=0.02
@test dm.γmin ≈ 0.63    atol=0.02
@test dm.γmax ≈ 1.59    atol=0.02
@test dm.α ≈ 0.46   atol=0.02
show(dm)
plot(dm)
nyquistplot(L)
plot!(dm, nyquist=true)
plot!(Disk(dm), nyquist=true)

@test dm.alpha == dm.α
@test length(dm.gainmargin) == 2
@test dm.phasemargin == dm.ϕm
@test :gainmargin ∈ propertynames(dm)

## Frequency-dependent margin
w = exp10.(LinRange(-2, 2, 500))
dms = diskmargin(L, 0, w)
plot(dms)

##
s = tf("s")
L = 6.25*(s + 3)*(s + 5) / (s*(s + 1)^2 *(s^2 + 0.18s + 100))

## αmax > 2
dm = diskmargin(L, 0, 200)
@test dm.γmax < dm.γmin
@test dm.γmin ≈ 0 atol = 1e-6
@test dm.ϕm ≈ 90 atol = 1e-4

w = exp10.(LinRange(-1, 2, 500))
dms = diskmargin(L, 0, w)
plot(dms)

# using IntervalArithmetic
# δ(a=1) = Complex(-a..a, -a..a)
# Δ(n, a) = diagm([δ(a) for _ in 1:n])
# M = [0 1; -0.1 -0.1]
# D = Δ(2, 1)
# 0 ∈ det(I-M*D)

## Loop at a time
a = 10
P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0)
K = ss(I(2))

Li = K*P
Lo = P*K

@test tf(minreal(broken_feedback(Li, 1))) ≈ tf(1, [1, 0])
@test tf(minreal(broken_feedback(Lo, 1))) ≈ tf(1, [1, 0])

dm = loop_diskmargin(Li)[1]
@test dm.α ≈ 2
@test dm.γmax > 1e10
@test dm.ϕm ≈ 90

dmm = loop_diskmargin(P, K)
dmm.input[1].α == dm.α

##
L3 = let
    tempA = [1.0 0.0 9.84 0.0 0.0; 0.0 1.0 0.01 2.14634e-6 0.0; 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 1.0 -1.73983959887; 0.0 0.0 0.0 0.0 0.56597684805]
    tempB = [0.0 4.901416e-5 0.00019883877999999998; 0.0 0.0 0.0; 0.0 0.0 4.0414389999999996e-5; 0.0 -0.02004649371 0.0; 0.0 -0.00490141631 0.0]
    tempC = [0.0 -0.83516488404 0.0 0.0 0.0; 186.74725411661 0.0 0.0 0.0 0.0; -7.44299057498 0.0 7035.08410814126 0.0 0.0]
    tempD = [0.0 0.0 0.0; 34875.36444283988 0.0 0.0; 48304.01940122544 0.0 0.0]
    ss(tempA, tempB, tempC, tempD, 0.01)
end

w = 2π .* exp10.(LinRange(-2, 2, 300))

dm = loop_diskmargin(L3, 0, 4.05)
@test dm[1].α ≈ 0.794418036911981 rtol=1e-3

dm = diskmargin(L3, ss(I(3), L3.Ts), 0, w)
plot(dm.simultaneous)
plot!(dm.simultaneous_input)
plot!(dm.simultaneous_output)

plot(dm.input)
plot!(dm.output)

## MIMO

a = 10
P = ss([0 a; -a 0], I(2), [1 a; -a 1], 0)
K = ss(I(2))

w = 2π .* exp10.(LinRange(-2, 2, 300))
# @time bisect_a(P, K, w)
##

# break at input (pass outputs through)
a = bisect_a(P, K, w; Z = [], tol=1e-4)
@test minimum(a) < 0.0998 # from DM paper

# break at output (pass inputs through)
a = bisect_a(P, K, w; W = [], tol=1e-4)
@test minimum(a) < 0.0998

# break at both input and output
a = bisect_a(P, K, w; tol=1e-4)
au = bisect_a(P, K, w; tol=1e-4, N=640, upper=true)
@test minimum(a) < 0.0499

ar = bisect_a(P, K, w; tol=1e-4, δ=δr)
aur = bisect_a(P, K, w; tol=1e-4, N=640, upper=true, δ=δr)

dm = diskmargin(P, K, 0, w)
adm = [dm.α for dm in dm.simultaneous]

if isinteractive()
    plot(w,  a,   xscale=:log10, xlabel="Frequency", ylims=(0,3), lab="Lower Complex")
    plot!(w, au,  xscale=:log10, xlabel="Frequency", ylims=(0,3), lab="Upper Complex")
    plot!(w, ar,  xscale=:log10, xlabel="Frequency", ylims=(0,3), lab="Lower Real")
    plot!(w, aur, xscale=:log10, xlabel="Frequency", ylims=(0,3), lab="Upper Real")
    plot!(w, adm, xscale=:log10, xlabel="Frequency", ylims=(0,3), lab="μ")
    display(current())
end

##
w = 2π .* exp10.(LinRange(-2, 2, 300))

L3 = let
    tempA = [1.0 0.0 9.84 0.0 0.0; 0.0 1.0 0.01 2.14634e-6 0.0; 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 1.0 -1.73983959887; 0.0 0.0 0.0 0.0 0.56597684805]
    tempB = [0.0 4.901416e-5 0.00019883877999999998; 0.0 0.0 0.0; 0.0 0.0 4.0414389999999996e-5; 0.0 -0.02004649371 0.0; 0.0 -0.00490141631 0.0]
    tempC = [0.0 -0.83516488404 0.0 0.0 0.0; 186.74725411661 0.0 0.0 0.0 0.0; -7.44299057498 0.0 7035.08410814126 0.0 0.0]
    tempD = [0.0 0.0 0.0; 34875.36444283988 0.0 0.0; 48304.01940122544 0.0 0.0]
    ss(tempA, tempB, tempC, tempD, 0.01)
end

a = bisect_a(L3, ss(I(3), L3.Ts), w; tol=2e-3)
au = bisect_a(L3, ss(I(3), L3.Ts), w; tol=2e-3, upper=true, N=256)
plot(w, a, xscale=:log10, xlabel="Frequency", ylims=(0,Inf))
plot!(w, au, xscale=:log10, xlabel="Frequency", ylims=(0,Inf))

dm = diskmargin(L3, 0, w)
plot(dm)

## Test diskmargin with particles in the system

L3 = let
    tempA = [1.0 0.0 9.84 0.0 0.0; 0.0 1.0 0.01 2.14634e-6 0.0; 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 1.0 -1.73983959887; 0.0 0.0 0.0 0.0 0.56597684805]
    tempB = [0.0 4.901416e-5 0.00019883877999999998; 0.0 0.0 0.0; 0.0 0.0 4.0414389999999996e-5; 0.0 -0.02004649371 0.0; 0.0 -0.00490141631 0.0]*(1 + 0.1*Particles(32))
    tempC = [0.0 -0.83516488404 0.0 0.0 0.0; 186.74725411661 0.0 0.0 0.0 0.0; -7.44299057498 0.0 7035.08410814126 0.0 0.0]
    tempD = [0.0 0.0 0.0; 34875.36444283988 0.0 0.0; 48304.01940122544 0.0 0.0]
    ss(tempA, tempB, tempC, tempD, 0.01)
end

unsafe_comparisons(true)
dm = loop_diskmargin(L3, 0, 4.05)
@test dm[1].α isa Particles
@test pmean(dm[1].α) ≈ 0.794418036911981 rtol=0.01

L3 = let
    tempA = [1.0 0.0 9.84 0.0 0.0; 0.0 1.0 0.01 2.14634e-6 0.0; 0.0 0.0 1.0 0.0 0.0; 0.0 0.0 0.0 1.0 -1.73983959887; 0.0 0.0 0.0 0.0 0.56597684805]
    tempB = [0.0 4.901416e-5 0.00019883877999999998; 0.0 0.0 0.0; 0.0 0.0 4.0414389999999996e-5; 0.0 -0.02004649371 0.0; 0.0 -0.00490141631 0.0]
    tempC = [0.0 -0.83516488404 0.0 0.0 0.0; 186.74725411661 0.0 0.0 0.0 0.0; -7.44299057498 0.0 7035.08410814126 0.0 0.0]
    tempD = [0.0 0.0 0.0; 34875.36444283988 0.0 0.0; 48304.01940122544 0.0 0.0]
    ss(tempA, tempB, tempC, tempD, 0.01)
end

P = L3
K = ss(1.0I(3), L3.timeevol)

w = [0.01, 0.1, 1, 10, 100]
mu = 1e5*[2.961061403927230
2.961017564832289
0.029884614923336
0.029880921297303
0.000453111209111
0.000452684368322
0.000033811562301
0.000033712186847
0.000003489834242
0.000003488286208]

mu2 = structured_singular_value(permutedims(freqresp(L3, w), (2,3,1)))
@test mu2 ≈ mu[2:2:end]

##
w = exp10.(LinRange(-2, 2, 300))
K = ss(I(3), L3.timeevol)
L = feedback(L3, K)
mu = structured_singular_value(permutedims(freqresp(L, w), (2,3,1)))
plot(mu)
@test mu[1] ≈ 1 rtol=1e-3
@test maximum(mu) ≈ 2.503148529400597 rtol=1e-3
@test mu[end] ≈ 0.397718223553195 rtol=0.03

##
w = 2π .* exp10.(LinRange(-2, 2, 300))
# break at input (pass outputs through)
# M,D = RobustAndOptimalControl.get_M(P, K, w; W = [])
M,D = RobustAndOptimalControl.get_M(L3, w)

# Test stability of the computation
# @time mu = [structured_singular_value(M) for i = 1:10]
# plot(w, mu, xscale=:log10)

## compare structured_singular_value and bisect_a
# mu = structured_singular_value(M)
# a = bisect_a(P, K, w; W = [], tol=1e-4)
# au = bisect_a(P, K, w; W = [], tol=1e-4, upper=true, N=2560)

# plot(w, mu, xscale=:log10, lab="mu")
# plot!(w, inv.(a), xscale=:log10, lab="lower")
# plot!(w, inv.(au), xscale=:log10, lab="upper")

## ==================================================

##
w = 2π .* exp10.(LinRange(-2, 2, 300))
a = [-0.2 10;-10 -0.2]; b = I(2); c = [1 8;-10 1];
P = ss(a,b,c,0);
K = ss([1 -2;0 1]);
dm = diskmargin(K*P) # disk margins at plant inputs
dm = diskmargin(P*K); # disk margins at plant outputs
MMIO = diskmargin(P,K,0,w)
plot(MMIO.simultaneous, lab="simultaneous")
plot!(MMIO.simultaneous_output, lab="output")
plot!(MMIO.simultaneous_input, lab="input")

## Compare mu and diskmargin
# w = 2π .* exp10.(LinRange(-3, 3, 300))
# M,D = RobustAndOptimalControl.get_M(P, K, w; σ=0)
# mu = structured_singular_value(M)
# a = bisect_a(P, K, w; tol=1e-3, σ=0)
# au = bisect_a(P, K, w; tol=1e-3, upper=true, N=1024, σ=0)
# MMIO = diskmargin(P,K,0,w)

# dma = [dm.α for dm in MMIO.simultaneous]

# plot(w, dma)
# plot!(w, inv.(mu), xscale=:log10, lab="mu", sp=1)
# plot!(w, a, xscale=:log10, lab="lower", sp=1)
# plot!(w, au, xscale=:log10, lab="upper", sp=1)


text:
using ExprTools
using Test

@testset "ExprTools.jl" begin
    include("function.jl")
end


text:
# albums.jl
## https://developer.spotify.com/documentation/web-api/reference/albums/

## https://developer.spotify.com/documentation/web-api/reference/albums/get-album/
@doc """
# Get an Album(s)
**Summary**: Get Spotify catalog information for an album(s)

`album_id` _Required_: The Spotify ID for the album. Up to 20 albums can be passed by
comma separating the ID's\n
`market` _Optional_: An ISO 3166-1 alpha-2 country code. Default "US"\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/albums/get-album/)
""" ->
function album_get(album_id, market="US")
    return spotify_request("albums/$album_id?market=$market")
end

## https://developer.spotify.com/documentation/web-api/reference/albums/get-albums-tracks/
@doc """
# Get an Album's Tracks
**Summary**: Get Spotify catalog information about an album's tracks. Optional parameters can be used to limit the number of tracks returned

`album_id` _Required_: The Spotify ID for the album. Up to 20 albums can be passed by
`limit` _Optional_: The maximum number of tracks to return. Default 20\n
`offset` _Optional_: The index of the first track to return. Default 0\n
`market` _Optional_: An ISO 3166-1 alpha-2 country code. Default "US"\n

[Reference](https://developer.spotify.com/documentation/web-api/reference/albums/get-albums-tracks/)
""" ->
function album_get_tracks(album_id, limit=20, offset=0, market="US")
    return spotify_request("albums/$album_id/tracks?limit=$limit&offset=$offset&market=$market")
end


text:
# To build the documentation:
#    - julia --project="." make.jl
#    - empty!(ARGS); include("make.jl")
# To build the documentation without running the tests:
#    - julia --project="." make.jl preview
#    - push!(ARGS,"preview"); include("make.jl")

using Documenter, Literate, BetaML, Test

if "preview" in ARGS
    println("*** Attention: code in the tutorial will not be run/tested")
else
    println("*** Building documentation and testing tutorials...")
end

push!(LOAD_PATH,"../src/")

const _TUTORIAL_DIR = joinpath(@__DIR__, "src", "tutorials")
# Important: If some tutorial is removed but the md file is left, this may still  be used by Documenter
_TUTORIAL_SUBDIR = [
    #"Getting started",
    "Regression - bike sharing",
    "Classification - cars",
    "Clustering - Iris"
]

function link_example(content)
    edit_url = match(r"EditURL = \"(.+?)\"", content)[1]
    footer = match(r"^(---\n\n\*This page was generated using)"m, content)[1]
    content = replace(
        content, footer => "[View this file on Github]($(edit_url)).\n\n" * footer
    )
    return content
end

function _file_list(full_dir, relative_dir, extension)
    return map(
        file -> joinpath(relative_dir, file),
        filter(file -> endswith(file, extension), sort(readdir(full_dir))),
    )
end

"""
    _include_sandbox(filename)
Include the `filename` in a temporary module that acts as a sandbox. (Ensuring
no constants or functions leak into other files.)
"""
function _include_sandbox(filename)
    mod = @eval module $(gensym()) end
    return Base.include(mod, filename)
end

function literate_directory(dir)
    rm.(_file_list(dir, dir, ".md"))
    for filename in _file_list(dir, dir, ".jl")
        # `include` the file to test it before `#src` lines are removed. It is
        # in a testset to isolate local variables between files.
        if ! ("preview" in ARGS)
            @testset "$(filename)" begin
               _include_sandbox(filename)
             end
             Literate.markdown(
                 filename,
                 dir;
                 documenter = true,
                 postprocess = link_example
             )
        else
            Literate.markdown(
                filename,
                dir;
                documenter = true,
                postprocess = link_example,
                codefence =  "```text" => "```"
            )
        end

    end
    return nothing
end

println("Starting literating tutorials (.jl --> .md)...")
literate_directory.(joinpath.(_TUTORIAL_DIR, _TUTORIAL_SUBDIR))

println("Starting making the doculentation...")
makedocs(sitename="BetaML.jl Documentation",
         authors = "Antonello Lobianco",
         pages = [
            "Index" => "index.md",

            "Tutorial" => vcat("tutorials/Betaml_tutorial_getting_started.md",map(
                    subdir -> subdir => map(
                        file -> joinpath("tutorials", subdir, file),
                        filter(
                            file -> endswith(file, ".md"),
                            sort(readdir(joinpath(_TUTORIAL_DIR, subdir))),
                        ),
                    ),
                    _TUTORIAL_SUBDIR,
                ),
            ),
            "API (Reference manual)" => [
              "Perceptron" => "Perceptron.md",
              "Trees" => "Trees.md",
              "Nn"   => "Nn.md",
              "Clustering" => "Clustering.md",
              "Utils" => "Utils.md",
            ],
         ],
         format = Documenter.HTML(prettyurls = false),
         #strict = true,
         #doctest = false
)
println("Starting deploying the documentation...")
deploydocs(
    repo = "github.com/sylvaticus/BetaML.jl.git",
)


text:
using DelimitedFiles

function ReadData(fname::AbstractString)
    dataStr = readdlm(fname, String)

    data = zeros(Int, length(dataStr), length(dataStr[1]))
    for i = 1:length(dataStr)
        data[i,:] = map(x -> parse(Int, x), split(dataStr[i],""))
    end

    return data
end

function StepData!(data::Matrix)
    # need to do slicing to modify data in place
    data[:,:] += ones(Int, size(data)...)

    flashes = findall(x -> x == 10, data)
    while ! isempty(flashes)
        newFlashes = Array{CartesianIndex{2},1}(undef, 0)

        for flash in flashes
            x, y = flash[1], flash[2]
            xMin, xMax = x == 1 ? x : x-1, x == size(data)[1] ? x : x+1
            yMin, yMax = y == 1 ? y : y-1, y == size(data)[2] ? y : y+1
            data[xMin:xMax, yMin:yMax] += ones(Int, xMax-xMin+1, yMax-yMin+1)

            for i=xMin:xMax, j=yMin:yMax
                if data[i,j] == 10
                    push!(newFlashes, CartesianIndex(i,j))
                end
            end
        end

        flashes = newFlashes
    end

    map!(x -> x>=10 ? 0 : x, data, data)
end

function Part1(data::Matrix)
    nSteps = 100
    flashes = 0

    for i = 1:nSteps
        StepData!(data)
        flashes += count(iszero, data)
    end

    return flashes
end

function Part2(data::Matrix)
    steps = 0
    while ! all(iszero, data)
        StepData!(data)
        steps += 1
    end

    return steps
end

#datafile = "../data/test11.txt"
datafile = "../data/day11.txt"

data = ReadData(datafile)
println(Part1(data))

# need to re-read data since Part1 modifies it
data = ReadData(datafile)
println(Part2(data))


text:
module AtomicGraphNets

export AGNConv, AGNPool, AGNConvDEQ
include("layers.jl")
using .Layers: AGNConv, AGNPool, AGNConvDEQ

include("models.jl")
export build_CGCNN, build_SGCNN, build_CGCNN_DEQ

end