text:
export  Turtle,
        Forward,
        Turn,
        Pendown,
        Penup,
        Circle,
        Orientation,
        Rectangle,
        Push,
        Pop,
        Message,
        HueShift,
        Randomize_saturation,
        Pencolor,
        Penwidth,
        Pen_opacity_random,
        Reposition

"""
    Turtle()
    Turtle(O)
    Turtle(0, 0)
    Turtle(O, pendown=true, orientation=0, pencolor=(1.0, 0.25, 0.25))

Create a Turtle. You can command a turtle to move and draw "turtle graphics".

The commands (unusually for Julia) start with a capital letter, and angles are specified in degrees.

Basic commands are `Forward()`, `Turn()`, `Pendown()`, `Penup()`, `Pencolor()`,
`Penwidth()`, `Circle()`, `Orientation()`, `Rectangle()`, and `Reposition()`.

Others include `Push()`, `Pop()`, `Message()`, `HueShift()`, `Randomize_saturation()`,
`Reposition()`, and `Pen_opacity_random()`.
"""
mutable struct Turtle
    xpos::Float64
    ypos::Float64
    pendown::Bool
    orientation::Float64 # stored in radians but set in degrees
    pencolor::Tuple{Float64, Float64, Float64} # it's an RGB turtle
    function Turtle(xpos, ypos, pendown::Bool, orientation, pencolor=(0.0, 0.0, 0.0))
        new(xpos, ypos, pendown, orientation, pencolor)
    end
end

Turtle(x, y) = Turtle(x, y, true, 0, (0, 0, 0))
Turtle(pos::Point=O) = Turtle(pos.x, pos.y, true, 0, (0, 0, 0))
Turtle(pos::Point, pendown::Bool) = Turtle(pos.x, pos.y, pendown, 0, (0, 0, 0))
Turtle(pos::Point, pendown::Bool, orientation) = Turtle(pos.x, pos.y, pendown, orientation, (0, 0, 0))
Turtle(pos::Point, pendown::Bool, orientation, col::NTuple{3, Number}) = Turtle(pos.x, pos.y, pendown, orientation, col)
Turtle(pos::Point, pendown::Bool, orientation, r, g, b) = Turtle(pos.x, pos.y, pendown, orientation, (r, g, b))

Base.broadcastable(t::Turtle) = Ref(t)

# a stack to hold pushed/popped turtle positions
const queue = Array{Array{Float64, 1}, 1}()

"""
    Forward(t::Turtle, d=1)

Move the turtle forward by `d` units. The stored position is updated.
"""
function Forward(t::Turtle, d=1)
    oldx, oldy = t.xpos, t.ypos
    t.xpos += (d * cos(t.orientation))
    t.ypos += (d * sin(t.orientation))
    if t.pendown
        gsave()
        sethue(t.pencolor...)
        line(Point(oldx, oldy), Point(t.xpos, t.ypos), :stroke)
        grestore()
    end
end

"""
    Turn(t::Turtle, r=5.0)

Increase the turtle's rotation by `r` degrees. See also `Orientation`.
"""
function Turn(t::Turtle, r=5.0)
    t.orientation = mod2pi(t.orientation + deg2rad(r))
end

"""
    Orientation(t::Turtle, r=0.0)

Set the turtle's orientation to `r` degrees. See also `Turn`.
"""
function Orientation(t::Turtle, r=0.0)
    t.orientation = mod2pi(deg2rad(r))
end

"""
    Pendown(t::Turtle)

Put that pen down and start drawing.
"""
function Pendown(t::Turtle)
    t.pendown = true
end

"""
    Penup(t::Turtle)

Pick that pen up and stop drawing.
"""
function Penup(t::Turtle)
    t.pendown = false
end

"""
    Circle(t::Turtle, radius=1.0)

Draw a filled circle centered at the current position with the given radius.
"""
function Circle(t::Turtle, radius=1.0)
    gsave()
    sethue(t.pencolor...)
    circle(t.xpos, t.ypos, radius, :fill)
    grestore()
end

"""
    Rectangle(t::Turtle, width=10.0, height=10.0)

Draw a filled rectangle centered at the current position with the given radius.
"""
function Rectangle(t::Turtle, width=10.0, height=10.0)
    gsave()
    sethue(t.pencolor...)
    rect(t.xpos-width/2, t.ypos-height/2, width, height, :fill)
    grestore()
end

"""
    Push(t::Turtle)

Save the turtle's position and orientation on a stack.
"""
function Push(t::Turtle)
# save xpos, ypos, turn in a queue
    global queue
    push!(queue, [t.xpos, t.ypos, t.orientation])
end

"""
    Pop(t::Turtle)

Lift the turtle's position and orientation off a stack.
"""
function Pop(t::Turtle)
# restore xpos, ypos, turn from queue
    global queue
    if isempty(queue)
        @info("stack was already empty")
        t.xpos, t.ypos, t.orientation = 0.0, 0.0, 0.0
    else
        t.xpos, t.ypos, t.orientation = pop!(queue)
    end
end

"""
    Message(t::Turtle, txt)

Write some text at the current position.
"""
function Message(t::Turtle, txt)
    gsave()
    sethue(t.pencolor...)
    text(txt, t.xpos, t.ypos)
    grestore()
end

"""
    HueShift(t::Turtle, inc=1.0)

Shift the Hue of the turtle's pen forward by `inc`. Hue values range between
0 and 360. (Don't start with black, otherwise the saturation and brightness values
will be black.)
"""
function HueShift(t::Turtle, inc=1.0)
    r, g, b = t.pencolor
    hsv = convert(Colors.HSV, Colors.RGB(r, g, b))
    newhsv = (mod(hsv.h + inc, 360), hsv.s, hsv.v)
    newrgb = convert(Colors.RGB, Colors.HSV(newhsv[1], newhsv[2], newhsv[3]))
    sethue(newrgb)
    t.pencolor = (newrgb.r, newrgb.g, newrgb.b)
end

"""
    Randomize_saturation(t::Turtle)

Randomize the saturation of the turtle's pen color.
"""
function Randomize_saturation(t::Turtle)
    r, g, b = t.pencolor
    hsv = convert(Colors.HSV, Colors.RGB(r, g, b))
    newhsv = (hsv.h, rand(0.5:0.05:1.0), hsv.v)
    newrgb= convert(Colors.RGB, Colors.HSV(newhsv[1], newhsv[2], newhsv[3]))
    sethue(newrgb)
    t.pencolor = (newrgb.r, newrgb.g, newrgb.b)
end

"""
    Pencolor(t::Turtle, r, g, b)

Set the Red, Green, and Blue colors of the turtle.
"""
function Pencolor(t::Turtle, r, g, b)
    sethue(r, g, b)
    t.pencolor = (r, g, b)
end

function Pencolor(t::Turtle, col::ColorTypes.Colorant)
    temp = convert(RGBA, col)
    Pencolor(t::Turtle, temp.r, temp.g, temp.b)
end

function Pencolor(t::Turtle, col::AbstractString)
    temp = parse(RGBA, col)
    Pencolor(t::Turtle, temp.r, temp.g, temp.b)
end

Pencolor(t::Turtle, col::NTuple{3, Number}) = Pencolor(t, col...)

"""
    Reposition(t::Turtle, pos::Point)
    Reposition(t::Turtle, x, y)

Reposition: pick the turtle up and place it at another position.
"""
function Reposition(t::Turtle, x, y)
    t.xpos = x
    t.ypos = y
end

Reposition(t::Turtle, pos::Point) = Reposition(t, pos.x, pos.y)

"""
    Penwidth(t::Turtle, w)

Set the width of the line drawn.
"""
Penwidth(t::Turtle, w) = setline(w)

"""
    Pen_opacity_random(t::Turtle)

Change the opacity of the pen to some value at random.
"""
Pen_opacity_random(t::Turtle) = setopacity(rand())

# end


text:
# data
" Name of substance, temperature range, coefficent"
solidThermoData = Dict(
    "Co" => [[298., 700., 1394., 1768.],
    [10.9943 54.375 -55.5132 25.817 0.164533 -4.7033 30.3258 0;
    -204.576 515.582 -421.55 129.558 17.9267 94.6155 -272.856 0;
    -12418.9 15326.4 -7087.15 1167.31 3319.77 10138.7 -10472.8 0]],
    "CoO" => [[298., 1600., 3000.],
    [43.65 22.38373 -16.89386 6.556161 0.532263 -249.8643 102.6892 -237.7353;
    35.9629 14.1274 1.164304 -0.116488 5.161759 -237.8181 106.4037 -237.7353]],
    "Co3O4" => [[298., 1500. ],
    [177.0953 -100.4361 185.7792 -49.04861 -3.483951 -971.5415 331.2828 -910.02]],
    "Cu" => [[298., 1358.],
    [17.72891 28.0987 -31.25289 13.97243 0.068611 -6.056591 47.89592 0] ],
    "CuO" => [[298, 2000.],
    [48.56494 7.498607 -0.05598 0.013851 -0.760082 -173.4272 94.85128 -156.0632]],
    "Fe" => [ [298., 700., 1042., 1100., 1809.],
    [18.42868 24.64301 -8.91372 9.664706 -0.012643 -6.573022 42.51488 0;
    -57767.65 137919.7 -122773.2 38682.42 3993.08 24078.67 -87364.01 0;
    -325.8859 28.92876 0 0 411.9629 745.8231 241.8766 0;
    -776.7387 919.4005 -383.7184 57.08148 242.1369 697.6234 -558.3674 0] ],
    "FeO" => [ [298., 1650. ], #, 5000. ],
    [45.7512 18.78553 -5.952201 0.852779 -0.081265 -286.7429 110.312 -272.0441 #;
    #68.19920 -4.501232E-10 1.195227E-10 -1.064302E-11 -3.092680E-10 -281.4326 137.8377 -249.5321
    ] ],
    "Fe3O4" => [ [298., 900., 3000.],
    [104.2096 178.5108 10.6151 1.132534 -0.994202 -1163.336 212.0585 -1120.894;
    200.832 1.586435E-7 -6.661682E-8 9.452452E-9 3.186020E-8 -1174.135 388.079 -1120.894] ],
    "Fe2O3" => [ [298., 950., 1050., 2500.],
    [93.43834 108.3577 -50.86447 25.58683 -1.61133 -863.2094 161.0719 -825.5032;
    150.624 0 0 0 0 -875.6066 252.8814 -825.5032;
    110.9362 32.04714 -9.192333 0.901506 5.433677 -843.1471 228.3548 -825.5032] ],
    "MgAl2O4" => [ [298., 2408.],
    [146.5551 35.91295 -0.183978 0.031409 -3.659941 -2356.684 234.7416 -2299.112 ] ],
    "Ni" => [ [298., 600., 700., 1728.],
    [13.6916 82.49509 -174.9548 161.6011 -0.092417 -6.833644 27.669 0;
    1248.045 -1257.51 0 0 -165.1266 -788.8263 1212.891 0;
    16.49839 18.74913 -6.639841 1.717278 1.872051 -0.467675 51.67575 0] ],
    "NiO" => [ [298., 525., 565., 2228.],
    [54.0207 -61.579 20.700 299.55 -0.09328 -254.135 118.9857 -239.701;
    54.686 6.3974 0.0 0.0 0.0 -256.996 99.924 -239.701;
    44.20243 11.7264 -1.03832 -0.099521 0.64577 -250.8932 90.9702 -239.701] ]
    )

" Thermodynamic data coefficent for gases"
gasThermoData = Dict(
    "CH4" => [[298., 1300, 6000],
    [-0.703029 108.4773 -42.52157 5.862788 0.678565 -76.84376 158.7163 -74.87310;
    85.81217 11.26467 -2.114146 0.138190 -26.42221 -153.5327 224.4143 -74.87310]],
    "CO" => [[298., 1300, 6000],
    [25.56759 6.096130 4.054656 -2.671301 0.131021 -118.0089 227.3665 -110.5271;
    35.15070 1.300095 -0.205921 0.013550 -3.282780 -127.8375 231.7120 -110.5271]],
    "CO2" => [[298., 1200, 6000],
    [24.99735 55.18696 -33.69137 7.948387 -0.136638 -403.6075 228.2431 -393.5224;
    58.16639 2.720074 -0.492289 0.038844 -6.447293 -425.9186 263.6125 -393.5224]],
    "H2" => [[298., 1000., 2500.],
    [33.066178 -11.363417 11.432816 -2.772874 -0.158558 -9.980797 172.707974 0;
    18.563083 12.257357 -2.859786 0.268238 1.977990 -1.147438 156.288133 0.]],
    "H2O" => [[500., 1700],
    [30.09200 6.832514 6.793435 -2.534480 0.082139 -250.8810 223.3967 -241.8264]],
    "N2" => [ [100., 500., 2000., 6000.],
    [28.98641 1.853978 -9.647459 16.63537 0.000117 -8.671914 226.4168 0;
    19.50583 19.88705 -8.598535 1.369784 0.527601 -4.935202 212.39 0;
    35.51872 1.128728 -0.196103 0.014662 -4.55376 -18.97091 224.981 0] ],
    "O2" => [[100., 700., 2000., 6000.],
    [31.32234 -20.23531 57.86644 -36.50624 -0.007374 -8.903471 246.7945 0;
    30.03235 8.772972 -3.988133 0.788313 -0.741599 -11.32468 236.1663 0;
    20.91111 10.72071 -2.020498 0.146449 9.245722 5.337651 237.6185 0]]
)

" Molar mass of substance, [g/mol]"
molarMassData = Dict(
    "Co" => 58.93,
    "CoO" => 74.93,
    "Co3O4" => 240.80,
    "Cu" => 63.55,
    "CuO" => 79.55,
    "Fe" => 55.85,
    "FeO" => 71.85,
    "Fe3O4" => 231.56,
    "Fe2O3" => 159.69,
    "MgAl2O4" => 142.27,
    "N2" => 28,
    "Ni" => 58.69,
    "NiO" => 74.69,
    "O2" => 32,
)

ext: jl


text:
include("utils.jl")

TEST_GROUP = get(ENV, "MANIFOLDS_TEST_GROUP", "all")

@info "Manifolds.jl Test settings:\n\n" *
      "Testing Float32:  $(TEST_FLOAT32)\n" *
      "Testing Double64: $(TEST_DOUBLE64)\n" *
      "Testing Static:   $(TEST_STATIC_SIZED)\n\n" *
      "Test group:       $(TEST_GROUP)\n\n" *
      "Check test/utils.jl if you wish to change these settings."

@testset "Manifolds.jl" begin
    if TEST_GROUP ∈ ["all", "test_manifolds"]
        include_test("differentiation.jl")

        include_test("ambiguities.jl")

        @testset "utils test" begin
            Random.seed!(42)
            @testset "usinc_from_cos" begin
                @test Manifolds.usinc_from_cos(-1) == 0
                @test Manifolds.usinc_from_cos(-1.0) == 0.0
            end
            @testset "log_safe!" begin
                n = 8
                Q = qr(randn(n, n)).Q
                A1 = Matrix(Hermitian(Q * Diagonal(rand(n)) * Q'))
                @test exp(Manifolds.log_safe!(similar(A1), A1)) ≈ A1 atol = 1e-6
                A1_fail = Matrix(Hermitian(Q * Diagonal([-1; rand(n - 1)]) * Q'))
                @test_throws DomainError Manifolds.log_safe!(similar(A1_fail), A1_fail)

                T = triu!(randn(n, n))
                T[diagind(T)] .= rand.()
                @test exp(Manifolds.log_safe!(similar(T), T)) ≈ T atol = 1e-6
                T_fail = copy(T)
                T_fail[1] = -1
                @test_throws DomainError Manifolds.log_safe!(similar(T_fail), T_fail)

                A2 = Q * T * Q'
                @test exp(Manifolds.log_safe!(similar(A2), A2)) ≈ A2 atol = 1e-6
                A2_fail = Q * T_fail * Q'
                @test_throws DomainError Manifolds.log_safe!(similar(A2_fail), A2_fail)

                A3 = exp(SizedMatrix{n,n}(randn(n, n)))
                @test A3 isa SizedMatrix
                @test exp(Manifolds.log_safe!(similar(A3), A3)) ≈ A3 atol = 1e-6
                @test exp(Manifolds.log_safe(A3)) ≈ A3 atol = 1e-6

                A3_fail = Float64[1 2; 3 1]
                @test_throws DomainError Manifolds.log_safe!(similar(A3_fail), A3_fail)

                A4 = randn(ComplexF64, n, n)
                @test exp(Manifolds.log_safe!(similar(A4), A4)) ≈ A4 atol = 1e-6
            end
            @testset "isnormal" begin
                @test !Manifolds.isnormal([1.0 2.0; 3.0 4.0])
                @test !Manifolds.isnormal(complex.(reshape(1:4, 2, 2), reshape(5:8, 2, 2)))

                # diagonal
                @test Manifolds.isnormal(diagm(randn(5)))
                @test Manifolds.isnormal(diagm(randn(ComplexF64, 5)))
                @test Manifolds.isnormal(Diagonal(randn(5)))
                @test Manifolds.isnormal(Diagonal(randn(ComplexF64, 5)))

                # symmetric/hermitian
                @test Manifolds.isnormal(Symmetric(randn(3, 3)))
                @test Manifolds.isnormal(Hermitian(randn(3, 3)))
                @test Manifolds.isnormal(Hermitian(randn(ComplexF64, 3, 3)))
                x = Matrix(Symmetric(randn(3, 3)))
                x[3, 1] += eps()
                @test !Manifolds.isnormal(x)
                @test Manifolds.isnormal(x; atol=sqrt(eps()))

                # skew-symmetric/skew-hermitian
                skew(x) = x - x'
                @test Manifolds.isnormal(skew(randn(3, 3)))
                @test Manifolds.isnormal(skew(randn(ComplexF64, 3, 3)))

                # orthogonal/unitary
                @test Manifolds.isnormal(Matrix(qr(randn(3, 3)).Q); atol=sqrt(eps()))
                @test Manifolds.isnormal(
                    Matrix(qr(randn(ComplexF64, 3, 3)).Q);
                    atol=sqrt(eps()),
                )
            end
            @testset "realify/unrealify!" begin
                # round trip real
                x = randn(3, 3)
                @test Manifolds.realify(x, ℝ) === x
                @test Manifolds.unrealify!(similar(x), x, ℝ) == x

                # round trip complex
                x2 = randn(ComplexF64, 3, 3)
                x2r = Manifolds.realify(x2, ℂ)
                @test eltype(x2r) <: Real
                @test size(x2r) == (6, 6)
                x2c = Manifolds.unrealify!(similar(x2), x2r, ℂ)
                @test x2c ≈ x2

                # matrix multiplication is preserved
                x3 = randn(ComplexF64, 3, 3)
                x3r = Manifolds.realify(x3, ℂ)
                @test x2 * x3 ≈ Manifolds.unrealify!(similar(x2), x2r * x3r, ℂ)
            end
            @testset "allocation" begin
                @test allocate([1 2; 3 4], Float64, Size(3, 3)) isa Matrix{Float64}
                @test allocate(SA[1 2; 3 4], Float64, Size(3, 3)) isa MMatrix{3,3,Float64}
                @test allocate(SA[1 2; 3 4], Size(3, 3)) isa MMatrix{3,3,Int}
            end
            @testset "eigen_safe" begin
                @test Manifolds.eigen_safe(SA[1.0 0.0; 0.0 1.0]) isa
                      Eigen{Float64,Float64,<:SizedMatrix{2,2},<:SizedVector{2}}
            end
        end

        @test Manifolds.is_metric_function(flat)
        @test Manifolds.is_metric_function(sharp)

        include_test("groups/group_utils.jl")
        include_test("notation.jl")
        # starting with tests of simple manifolds
        include_test("manifolds/centered_matrices.jl")
        include_test("manifolds/circle.jl")
        include_test("manifolds/cholesky_space.jl")
        include_test("manifolds/elliptope.jl")
        include_test("manifolds/euclidean.jl")
        include_test("manifolds/fixed_rank.jl")
        include_test("manifolds/generalized_grassmann.jl")
        include_test("manifolds/generalized_stiefel.jl")
        include_test("manifolds/grassmann.jl")
        include_test("manifolds/hyperbolic.jl")
        include_test("manifolds/lorentz.jl")
        include_test("manifolds/multinomial_doubly_stochastic.jl")
        include_test("manifolds/multinomial_symmetric.jl")
        include_test("manifolds/positive_numbers.jl")
        include_test("manifolds/probability_simplex.jl")
        include_test("manifolds/projective_space.jl")
        include_test("manifolds/rotations.jl")
        include_test("manifolds/skewhermitian.jl")
        include_test("manifolds/spectrahedron.jl")
        include_test("manifolds/sphere.jl")
        include_test("manifolds/sphere_symmetric_matrices.jl")
        include_test("manifolds/stiefel.jl")
        include_test("manifolds/symmetric.jl")
        include_test("manifolds/symmetric_positive_definite.jl")
        include_test("manifolds/symmetric_positive_semidefinite_fixed_rank.jl")
        include_test("manifolds/symplectic.jl")
        include_test("manifolds/symplecticstiefel.jl")
        include_test("manifolds/tucker.jl")

        include_test("manifolds/essential_manifold.jl")
        include_test("manifolds/multinomial_matrices.jl")
        include_test("manifolds/oblique.jl")
        include_test("manifolds/torus.jl")

        #meta manifolds
        include_test("manifolds/product_manifold.jl")
        include_test("manifolds/power_manifold.jl")
        include_test("manifolds/vector_bundle.jl")
        include_test("manifolds/graph.jl")

        include_test("metric.jl")
        include_test("statistics.jl")
        include_test("approx_inverse_retraction.jl")
    end

    if TEST_GROUP ∈ ["test_lie_groups", "all"]
        # Lie groups and actions
        include_test("groups/groups_general.jl")
        include_test("groups/validation_group.jl")
        include_test("groups/circle_group.jl")
        include_test("groups/translation_group.jl")
        include_test("groups/general_linear.jl")
        include_test("groups/special_linear.jl")
        include_test("groups/special_orthogonal.jl")
        include_test("groups/heisenberg.jl")
        include_test("groups/product_group.jl")
        include_test("groups/semidirect_product_group.jl")
        include_test("groups/special_euclidean.jl")
        include_test("groups/group_operation_action.jl")
        include_test("groups/rotation_action.jl")
        include_test("groups/translation_action.jl")
        include_test("groups/connections.jl")
        include_test("groups/metric.jl")
    end
    if TEST_GROUP ∈ ["all", "test_manifolds"]
        include_test("recipes.jl")
    end
end


text:
struct DimensionalityReduction <: AbstractProcessingStep
    table::Observable{SimpleTable}
    card::ProcessingCard
end

function columns_out(step::DimensionalityReduction)
    card = step.card
    isempty(columns_in(card)) && return Symbol[]
    name = extract_positional_argument(card.outputs; strict=false)
    kwargs = extract_named_arguments(card.method; strict=false)
    dims = 0
    if all(!isnothing, (name, kwargs))
        for (k, v) in kwargs
            k == "dims" && (dims = parse(Int, v))
        end
    end
    return Symbol.(name, '_', 1:dims)
end

# Add custom type to represent multi-dimensional scaling
struct MDS end

# Unify interface on how to pass dims
_fit(::Type{PCA}, X; dims, kwargs...) = fit(PCA, X; maxoutdim=dims, pratio=1, kwargs...)
_fit(::Type{PPCA}, X; dims, kwargs...) = fit(PPCA, X; maxoutdim=dims, kwargs...)
_fit(::Type{FactorAnalysis}, X; dims, kwargs...) = fit(FactorAnalysis, X; maxoutdim=dims, kwargs...)
_fit(::Type{ICA}, X; dims, kwargs...) = fit(ICA, X, dims; kwargs...)

function test_dims(output; dims)
    _dims = size(output, 1)
    (_dims < dims) && error("Could only find $_dims dimensions. $dims were requested.")
    (_dims > dims) && error("Found too many dimensions.") # Should never happen
    return output
end

function project(an, data; dims, kwargs...)
    anres = _fit(an, data; dims, kwargs...)
    output = transform(anres, data)
    return test_dims(output; dims)
end

function project(::Type{MDS}, data; dims, distance=Euclidean(), kwargs...)
    output = classical_mds(pairwise(distance, data, dims=2), dims; kwargs...)
    return test_dims(output; dims)
end

const dimensionalityreductions = (
    pca=PCA,
    ppca=PPCA,
    factoranalysis=FactorAnalysis,
    ica=ICA,
    mds=MDS,
)

function DimensionalityReduction(table::Observable{SimpleTable})

    analysis_names = collect(map(string, keys(dimensionalityreductions)))

    analysis_options = vecmap(analysis_names) do name
        return name * " dims" => [string(i) for i in 1:100]
    end
    pushfirst!(analysis_options, "+" => String[])

    default_names = ":projection"

    wdgs = (
        inputs=RichTextField("Inputs", data_options(table, keywords=[""]), ""),
        outputs=RichTextField("Outputs", ["" => ["projection"]], default_names),
        method=RichTextField("Method", analysis_options, ""),
    )

    card = ProcessingCard(:Project; wdgs...)
    return DimensionalityReduction(table, card)
end

function (dimres::DimensionalityReduction)(data)
    card = dimres.card
    args, kwargs = extract_all_arguments(card.inputs)
    name = extract_positional_argument(card.outputs)

    cols = Tables.getcolumn.(Ref(data), Symbol.(args))
    X = reduce(vcat, transpose.(cols))
    options = Pair{Symbol, Any}[Symbol(k) => Tables.getcolumn(data, Symbol(v)) for (k, v) in kwargs]
    for (k, v) in extract_named_arguments(card.method)
        if k == "dims"
            push!(options, :dims => parse(Int, v))
        else
            push!(options, Symbol(k) => v)
        end
    end

    f = extract_function(card.method)
    an = dimensionalityreductions[Symbol(f)]
    projected_data = project(an, X; options...)
    rows = eachrow(projected_data)
    return SimpleTable(
        Symbol.(name, '_', eachindex(rows)),
        collect(AbstractVector, rows)
    )
end


text:
## Might be more elegant and concise to have Fermi terms and Pauli terms be
## OpTerm{FermiOp} and OpTerm{Pauli}, instead of FermiTerm, etc.
## We implement this idea here.

"""
    struct OpTerm{W<:AbstractOp, T<:AbstractVector{W}, CoeffT} <: AbstractTerm{W, CoeffT}

The operator term
"""
struct OpTerm{W<:AbstractOp, T<:AbstractVector{W}, CoeffT} <: AbstractTerm{W, CoeffT}
    ops::T
    coeff::CoeffT
end

"""
    op_string(t::OpTerm)

Return the string of operators in `t`, omitting the coefficient.

# Examples
```jldocstring
julia> op_string(PauliTerm("XY", 2.0))
2-element Vector{Pauli}:
 Pauli: X
 Pauli: Y
```
"""
op_string(t::OpTerm) = t.ops

"""
    term_type(::Type{T}) where T <: AbstractOp

Return the concrete type `OpTerm{T}` associated with `T`.
This is used in generic code to construct an `OpTerm` when `T` is known.
"""
term_type(::Type{T}) where {T <: AbstractOp} = OpTerm{T}

"""
    strip_typeof(::OpTerm{W, T, CoeffT}) where {W, T, CoeffT}

Return the type of the argument with concrete `T` and `CoeffT` replaced
by variables. `strip_typeof` is useful when a constructor is available only
for the simpler type returned by `strip_typeof`.
# Examples
```jldoctest
julia> typeof(PauliTerm("XX"))
PauliTerm{Vector{Pauli}, Complex{Int64}} (alias for OpTerm{Pauli, Array{Pauli, 1}, Complex{Int64}})

julia> strip_typeof(PauliTerm("XX"))
PauliTerm (alias for OpTerm{Pauli, T} where T<:AbstractArray{Pauli, 1})

julia> typeof(OpTerm{PauliI}("XX"))
OpTerm{PauliI, Vector{PauliI}, Complex{Int64}}

julia> QuantumOps.strip_typeof(OpTerm{PauliI}("XX"))
OpTerm{PauliI, T} where T<:AbstractVector{PauliI}
```

The following attempt to construct an operator will throw an exception because the requested
method does not exist.
```julia
julia> typeof(OpTerm{PauliI}("XX"))("YY")
ERROR: MethodError: no method matching OpTerm{PauliI, Vector{PauliI}, Complex{Int64}}(::String)
```

But, this works as intended.
```jldoctest
julia> QuantumOps.strip_typeof(OpTerm{PauliI}("XX"))("YY")
2-factor OpTerm{PauliI, Vector{PauliI}, Complex{Int64}}:
YY * (1 + 0im)
```
"""
strip_typeof(::OpTerm{W, T, CoeffT}) where {W, T, CoeffT} = OpTerm{W}

"""
    DenseOpTerm

An alias for `OpTerm` with dense storage.
"""
const DenseOpTerm = OpTerm{W, T} where {W<:AbstractOp, T<:Union{Vector{W}, StaticArrays.SVector}}

## TODO: Commented out and tests pass. Is this necessary?
# OpTerm(term::AbstractVector, coeff=_DEFAULT_COEFF) = OpTerm(term, coeff)

"""
    OpTerm{T}(term::V, coeff=_DEFAULT_COEFF) where {T, V<:AbstractVector{T}}

Construct an `OpTerm{T}` from a vector of `T`.

# Examples
```jldoctest
julia> OpTerm{Pauli}([X, Y, Z])
3-factor PauliTerm{Vector{Pauli}, Complex{Int64}}:
XYZ * (1 + 0im)
```
"""
OpTerm{T}(term::V, coeff=_DEFAULT_COEFF) where {T, V<:AbstractVector{T}} = OpTerm(term, coeff)

OpTerm(term::V, coeff=_DEFAULT_COEFF) where {T, V<:AbstractVector{T}} = OpTerm(term, coeff)

"""
    OpTerm{T}(s::AbstractString, coeff=_DEFAULT_COEFF)

Construct `OpTerm{T}` from a string representation.

# Examples
```jldoctest
julia> OpTerm{Pauli}("XYZI")
4-factor PauliTerm{Vector{Pauli}, Complex{Int64}}:
XYZI * (1 + 0im)

julia> OpTerm{FermiOp}("++--")
4-factor FermiTerm{Vector{FermiOp}, Complex{Int64}}:
++-- * (1 + 0im)
```

Recall that `FermiTerm` is an alias for `OpTerm{FermiOp}`.

```jldoctest
julia> FermiTerm("++--", 2.0)
4-factor FermiTerm{Vector{FermiOp}, Complex{Int64}}:
++-- * 2.0
```
"""
OpTerm{T}(s::AbstractString, coeff) where T = OpTerm(Vector{T}(s), coeff)

## TODO: Splitting is expensive because the entire string is scanned if there is no space.
## instead we may split on the space before the first operator character.
## But, FermiOp uses characters that may also be in numbers... Hmm...
## I'm not sure I want to support this. Parsing the string is tricky and I don't know
## what the use case is.
## It would not be hard t do this efficiently for Pauli, which has op chars very different
## from the numbers.
function OpTerm{T}(s::AbstractString) where T
    return OpTerm{T}(s, _DEFAULT_COEFF)
    # strs = split(s, ' '; limit=2)
    # if length(strs) == 1
    #     return OpTerm{T}(s, _DEFAULT_COEFF)
    # elseif length(strs) == 2
    # end
end

# This works. But OpSum{Pauli} will fail because it will not
# call this method.
# function OpTerm{T}(s::AbstractString) where T <: AbstractPauli
#     j = 0
#     for c in s
#         if c in "IXYZ"
#             break
#         end
#         j += 1
#     end
#     if j == 0
#         return OpTerm{T}(s, _DEFAULT_COEFF)
#     end
#     strs = split(s, ' '; limit=2)
#     if length(strs) == 2
#         return OpTerm{T}(strs[2], parse(ComplexF64, strs[1]))
#     else
#         throw(ArgumentError("Unable to parse string for PauliTerm"))
#     end
# end

## TODO: document and give examples for all signatures below
"""
    OpTerm{OpT}(inds::AbstractVector{<:Integer}, coeff=_DEFAULT_COEFF) where OpT <: AbstractOp
    OpTerm{OpT}(ops::AbstractOp...; coeff=_DEFAULT_COEFF)
    OpTerm(ops::AbstractOp...; coeff=_DEFAULT_COEFF)
    OpTerm{OpT}(s::Symbol, coeff=_DEFAULT_COEFF) where OpT

Construct an `OpTerm{OpT}` from an array of integer indices (codes) `inds` for `OpT`,
from arguments `ops`, etc.

# Examples

```jldoctest
julia> PauliTerm([0,1,2,3])
4-factor PauliTerm{Vector{Pauli}, Complex{Int64}}:
IXYZ * (1 + 0im)

julia> FermiTerm([0,1,2,3,4,5])
6-factor FermiTerm{Vector{FermiOp}, Complex{Int64}}:
INE+-0 * (1 + 0im)

julia> OpTerm{PauliI}([3, 2, 1, 0])
4-factor OpTerm{PauliI, Vector{PauliI}, Complex{Int64}}:
ZYXI * (1 + 0im)

julia> PauliTerm(0:3)
4-factor PauliTerm{Vector{Pauli}, Complex{Int64}}:
IXYZ * (1 + 0im)
```
"""
function OpTerm{OpT}(inds::AbstractVector{<:Integer}, coeff=_DEFAULT_COEFF) where OpT <: AbstractOp
#    return OpTerm(OpT.(inds), coeff)
    return OpTerm([OpT(i) for i in inds], coeff)
end

## Documented elsewhere (above, when this comment was written)
function OpTerm{OpT}(ops::AbstractOp...; coeff=_DEFAULT_COEFF) where OpT
    if isempty(ops)
        return OpTerm(OpT[], coeff)
    end
    return OpTerm([ops...], coeff)
end

## Documented elsewhere (above, when this comment was written)
function OpTerm(ops::AbstractOp...; coeff=_DEFAULT_COEFF)
    return OpTerm([ops...], coeff)
end

function OpTerm{OpT}(s::Symbol, coeff=_DEFAULT_COEFF) where OpT
    return OpTerm{OpT}(Vector{OpT}(String(s)), coeff)
end

function _op_term_macro_helper(::Type{T}, str) where T
    v = split(str, "*")
    if length(v) == 1
        return :($T($str))
    end
    if length(v) == 2
        num = parse(ComplexF64, v[1])
        return :($T($(strip(v[2])), $num))
    end
    throw(ErrorException("bad"))
end

"""
    rand_op_term(::Type{OpT}, n::Integer; coeff=_DEFAULT_COEFF) where {OpT <: AbstractOp}

Return a random `OpTerm` of `n` tensor factors.
"""
rand_op_term(::Type{OpT}, n::Integer; coeff=_DEFAULT_COEFF) where {OpT <: AbstractOp} =
    term_type(OpT)(rand(OpT, n), coeff)

####
#### Algebra
####

function multiply_keeping_phase!(target::AbstractArray{<:AbstractOp},
                                s1::AbstractArray{<:AbstractOp}, s2::AbstractArray{<:AbstractOp})
    length(s1) != length(s2) && throw(DimensionMismatch())
    _phase_data = AbstractOps.phase_data(eltype(s1))
    @inbounds for i in eachindex(s1)
        AbstractOps.accumulate_phase!(_phase_data, s1[i], s2[i])
        _prod = s1[i] * s2[i]
        if iszero(_prod)  # If any factor is zero, the entire product is zero
            fill!(target, Zero)
            return (target, 1)
        end
        target[i] = _prod
    end
    cum_phase = AbstractOps.compute_phase(_phase_data, s1, s2)
    return(target, cum_phase)
end

function mul!(target::AbstractArray{<:AbstractOp}, ps1::DenseOpTerm, ps2::DenseOpTerm)
    s_new, phase = multiply_keeping_phase!(target, op_string(ps1), op_string(ps2))
    if iszero(s_new)
        return strip_typeof(ps1)(s_new, zero(ps1.coeff))
    end
    return strip_typeof(ps1)(s_new, ps1.coeff * ps2.coeff * phase)
end

function Base.:*(ps1::DenseOpTerm, ps2::DenseOpTerm)
    return mul!(similar(op_string(ps1)), ps1, ps2)
end

function Base.zero(term::OpTerm{T}) where T
    new_string = similar(op_string(term))
    fill!(new_string, zero(T))
    return strip_typeof(term)(new_string, zero(term.coeff))
end


text:
# Automatic differentiation on sn.jl:
include("sn.jl")

using ForwardDiff
using DiffResults

function sn_jac(l_max::Int64,r::T,b::T) where {T <: Real}
  # Computes the derivative of s_n(r,b) with respect to r, b.
  # Create a vector for use with ForwardDiff
  x=[r,b]
  # Compute the length of the vector s_n:
  n_max = l_max^2+2*l_max
  # Allocate an array for s_n:
  sn = zeros(typeof(r),n_max+1)

  # Now, define a wrapper of s_n! for use with ForwardDiff:
  function diff_sn(x::Array{T,1}) where {T <: Real}
  # x should be a two-element vector with values [r,b]
  r_tmp,b_tmp = x
  sn = zeros(typeof(r_tmp),n_max+1)
  s_n!(l_max,r_tmp,b_tmp,sn)
  return sn
  end

  # Set up a type to store s_n and it's Jacobian with respect to x:
  out = DiffResults.JacobianResult(sn,x)
  # Compute the Jacobian (and value):
  out = ForwardDiff.jacobian!(out,diff_sn,x)
  # Place the value in the s_n vector:
  sn = DiffResults.value(out)
  # And, place the Jacobian in an array:
  sn_jacobian = DiffResults.jacobian(out)
  # Add in s2:
  s2_grad=zeros(typeof(r),2)
  s_2,Eofk,Em1mKdm=s2!(r,b,s2_grad)
  sn_jacobian[3,:]=s2_grad
return sn,sn_jacobian
end

# Here are some example lines of how to use it:
#r = 100.0; b=100.5
#l_max = 10
#sn,sn_jacobian= sn_jac(l_max,r,b)

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
# ===============================
# Written by AAE
# TU Delft, Spring 2014
# simulkade.com
# ===============================

# ====================== 1D CARTESIAN MESH =======================
"""
m = createMesh1D(Nx::Int, Width::Real)
It creates a uniform mesh on a 1D Cartesian domain.

Inputs:

   + Nx: number of cells in the domain (integer)
	 + Width: width or length of the domain (Real)

Outputs:

   + m: a mesh structure

Usage:
```
Nx=10
Lx=1.0
m=createMesh1D(Nx, Lx)
```
"""
function createMesh1D(Nx::Int, Width::Real)
# builds a uniform 1D mesh:
# Nx is the number of cells in x (horizontal) direction
# Width is the domain length in x direction

# mesh dimension

# cell size is dx
dx = Width/Nx

# numbering system of cells, like the single index numbering of Matlab
# +2 is added to account for the ghost cells that are added at
# the boundaries
MeshStructure(1,
		[Nx],
		CellSize(dx*ones(Nx+2), [0.0], [0.0]),
		CellLocation([1:Nx;]*dx.-dx/2,[0.0],[0.0]),
		FaceLocation([0:Nx;]*dx,[0.0],[0.0]),
		[1],
		[1])
end

"""
m = createMesh1D(facelocationX::Array{Real,1})
It creates a non/uniform mesh on a 1D Cartesian domain.

Inputs:

   + facelocationX: location of the the cell boundaries on
	 the x-axis.

Outputs:

   + m: a mesh structure

Usage:
```
	x= [0.0, 1.0, 1.4, 2.5, 4.1, 6.0, 10.0]
	m=createMesh1D(x)
```
"""
function createMesh1D(facelocationX::Array{T,1}) where T<:Real
# builds a uniform 1D mesh:
# facelocationX is the location of each cell face
# mesh dimension
# cell size is dx
Nx = length(facelocationX)-1
# numbering system of cells, like the single index numbering of Matlab
# +2 is added to account for the ghost cells that are added at
# the boundaries
MeshStructure(1,
		[Nx],
		CellSize([facelocationX[2]-facelocationX[1];
		facelocationX[2:end]-facelocationX[1:end-1];
		facelocationX[end]-facelocationX[end-1]], [0.0], [0.0]),
		CellLocation(0.5*(facelocationX[2:end]+facelocationX[1:end-1]),[0.0],[0.0]),
		FaceLocation(facelocationX,[0.0],[0.0]),
		[1],
		[1])
end

# ================= 1D CYLINDRICAL MESH ==========================
"""
m = createMeshCylindrical1D(Nr::Int, Radius::Real)
It creates a uniform mesh on a 1D Radial domain.

Inputs:

   + Nr: number of cells in the domain (integer)
	 + Radius: Radius of the domain (Real)

Outputs:

   + m: a mesh structure

Usage:

Nr=10
R=1.0
m=createMesh1D(Nr, R)
"""
function createMeshCylindrical1D(Nr::Int, Radius::Real)
# builds a uniform 1D cylindrical mesh:
# Nx is the number of cells in r (radial) direction
# Radius is the domain length in r direction
# mesh dimension
# cell size is dr
dr = Radius/Nr
# numbering system of cells, like the single index numbering of Matlab
# +2 is added to account for the ghost cells that are added at
# the boundaries
MeshStructure(1.5,
		[Nr],
		CellSize(dr*ones(Nr+2), [0.0], [0.0]),
		CellLocation([1:Nr;]*dr.-dr/2.0,[0.0],[0.0]),
		FaceLocation([0:Nr;]*dr,[0.0],[0.0]),
		[1],
		[1])
end

"""
m = createMeshCylindrical1D(facelocationR::Array{Real,1})
It creates a non/uniform mesh on a 1D radial domain.

Inputs:

   + facelocationX: location of the the cell boundaries on
	 the r-axis.

Outputs:

   + m: a mesh structure

Usage:

	r= [0.1, 1.0, 1.4, 2.5, 4.1, 6.0, 10.0]
	m=createMeshCylindrical1D(r)
"""
function createMeshCylindrical1D(facelocationR::Array{T,1}) where T<:Real
# builds a uniform 1D cylindrical mesh:
# Nx is the number of cells in r (radial) direction
# Radius is the domain length in r direction
# mesh dimension
# cell size is dr
Nr = length(facelocationR)-1
# numbering system of cells, like the single index numbering of Matlab
# +2 is added to account for the ghost cells that are added at
# the boundaries
MeshStructure(1.5,
		[Nr],
		CellSize([facelocationR[2]-facelocationR[1];
		facelocationR[2:end]-facelocationR[1:end-1];
		facelocationR[end]-facelocationR[end-1]], [0.0], [0.0]),
		CellLocation(0.5*(facelocationR[2:end]+facelocationR[1:end-1]),[0.0],[0.0]),
		FaceLocation(facelocationR,[0.0],[0.0]),
		[1],
		[1])
end

# ========================= 2D CARTESIAN MESH ==========================
"""
m = createMesh2D(Nx::Int, Ny::Int, Width::Real, Height::Real)
It creates a uniform mesh on a 2D Cartesian domain.

Inputs:

   + Nx: number of cells in the x-direction (integer)
   + Ny: number of cells in the y-direction (integer)
	 + Width: length of the domain in the x-direction (Real)
   + Height: length of the domain in the y-direction (Real)

Outputs:

   + m: a mesh structure

Usage:

Nx=10
Ny=15
Lx=1.0
Ly=2.5
m=createMesh2D(Nx, Ny, Lx, Ly)
"""
function createMesh2D(Nx::Int, Ny::Int, Width::Real, Height::Real)
# builds a uniform 2D mesh:
# Nx is the number of cells in x (horizontal) direction
# Width is the domain length in x direction

# numbering system of cells, like the single index numbering of Matlab
# +2 is added to account for the ghost cells that are added at
# the boundaries
# cell size is dx
dx = Width/Nx
dy = Height/Ny
G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)
MeshStructure(2,
	[Nx, Ny],
	CellSize(dx*ones(Nx+2), dy*ones(Ny+2), [0.0]),
	CellLocation([1:Nx;]*dx.-dx/2.0,[1:Ny;]*dy.-dy/2.0,[0.0]),
	FaceLocation([0:Nx;]*dx,[0:Ny;]*dy,[0.0]),
	G[[1,end],[1,end]][:],
	[1])
end

function createMesh2D(facelocationX::Array{T,1}, facelocationY::Array{T,1}) where T<:Real
Nx = length(facelocationX)-1
Ny = length(facelocationY)-1
G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)
MeshStructure(2,
	[Nx, Ny],
	CellSize([facelocationX[2]-facelocationX[1]; facelocationX[2:end]-facelocationX[1:end-1]; facelocationX[end]-facelocationX[end-1]],
	[facelocationY[2]-facelocationY[1]; facelocationY[2:end]-facelocationY[1:end-1]; facelocationY[end]-facelocationY[end-1]],
	[0.0]),
	CellLocation(0.5*(facelocationX[2:end]+facelocationX[1:end-1]), 0.5*(facelocationY[2:end]+facelocationY[1:end-1]), [0.0]),
	FaceLocation(facelocationX, facelocationY, [0.0]),
	G[[1,end],[1,end]][:],
	[1])
end

# ============================== 2D RADIAL MESH =======================================
function createMeshRadial2D(Nr::Int, Ntheta::Int, Radius::Real, Angle::Real)
# builds a uniform radial 2D mesh:
# Nr is the number of cells in r (radial) direction
# Radius is the domain length in r direction
# Ntheta is the number of cells in angular direction (not really a direction)
# Angle is the angular domain size
# numbering system of cells, like the single index numbering of Matlab
# +2 is added to account for the ghost cells that are added at
# the boundaries
# cell size is dx, dy
if Angle>2*pi
	Angle = 2*pi
	println("The domain size adjusted to match a maximum of 2*pi.")
end
dr = Radius/Nr
dtheta = Angle/Ntheta
G=reshape([1:(Nr+2)*(Ntheta+2);], Nr+2, Ntheta+2)
MeshStructure(2.8,
	[Nr, Ntheta],
	CellSize(dr*ones(Nr+2), dtheta*ones(Ntheta+2), [0.0]),
	CellLocation([1:Nr;]*dr.-dr/2.0,[1:Ntheta;]*dtheta.-dtheta/2.0,[0.0]),
	FaceLocation([0:Nr;]*dr,[0:Ntheta;]*dtheta,[0.0]),
	G[[1,end],[1,end]][:],
	[1])
end

function createMeshRadial2D(facelocationR::Array{T,1}, facelocationTheta::Array{T,1}) where T<:Real
if facelocationTheta[end]>2.0*pi
	facelocationTheta = facelocationTheta/facelocationTheta[end]*(2.0*pi)
	println("The domain size adjusted to match a maximum of 2*pi.")
end
Nx = length(facelocationR)-1
Ny = length(facelocationTheta)-1
G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)
MeshStructure(2.8,
	[Nx, Ny],
	CellSize([facelocationR[2]-facelocationR[1]; facelocationR[2:end]-facelocationR[1:end-1]; facelocationR[end]-facelocationR[end-1]],
	[facelocationTheta[2]-facelocationTheta[1]; facelocationTheta[2:end]-facelocationTheta[1:end-1]; facelocationTheta[end]-facelocationTheta[end-1]],
	[0.0]),
	CellLocation(0.5*(facelocationR[2:end]+facelocationR[1:end-1]), 0.5*(facelocationTheta[2:end]+facelocationTheta[1:end-1]), [0.0]),
	FaceLocation(facelocationR, facelocationTheta, [0.0]),
	G[[1,end],[1,end]][:],
	[1])
end

# ===================== 2D CYLINDRICAL MESH ==========================
function createMeshCylindrical2D(Nr::Int, Nz::Int, Radius::Real, Height::Real)
# builds a uniform cylindrical 2D mesh:
# Nr is the number of cells in r (radial) direction
# Radius is the domain length in r direction

# numbering system of cells, like the single index numbering of Matlab
# +2 is added to account for the ghost cells that are added at
# the boundaries
# cell size is dr, dz
dr = Radius/Nr
dz = Height/Nz
G=reshape([1:(Nr+2)*(Nz+2);], Nr+2, Nz+2)
MeshStructure(2.5,
	[Nr, Nz],
	CellSize(dr*ones(Nr+2), dz*ones(Nz+2), [0.0]),
	CellLocation([1:Nr;]*dr.-dr/2.0,[1:Nz;]*dz.-dz/2.0,[0.0]),
	FaceLocation([0:Nr;]*dr,[0:Nz;]*dz,[0.0]),
	G[[1,end],[1,end]][:],
	[1])
end

function createMeshCylindrical2D(facelocationR::Array{T,1}, facelocationY::Array{T,1}) where T<:Real
Nx = length(facelocationR)-1
Ny = length(facelocationY)-1
G=reshape([1:(Nx+2)*(Ny+2);], Nx+2, Ny+2)
MeshStructure(2.5,
	[Nx, Ny],
	CellSize([facelocationR[2]-facelocationR[1]; facelocationR[2:end]-facelocationR[1:end-1]; facelocationR[end]-facelocationR[end-1]],
	[facelocationY[2]-facelocationY[1]; facelocationY[2:end]-facelocationY[1:end-1]; facelocationY[end]-facelocationY[end-1]],
	[0.0]),
	CellLocation(0.5*(facelocationR[2:end]+facelocationR[1:end-1]), 0.5*(facelocationY[2:end]+facelocationY[1:end-1]), [0.0]),
	FaceLocation(facelocationR, facelocationY, [0.0]),
	G[[1,end],[1,end]][:],
	[1])
end

# ================================= 3D CARTESIAN MESH ========================================
function createMesh3D(Nx::Int, Ny::Int, Nz::Int, Width::Real, Height::Real, Depth::Real)
# builds a uniform 3D mesh:
# Nx is the number of cells in x direction
# Ny is the number of cells in y direction
# Nz is the number of cells in z direction
# Width is the domain length in x direction
# Height is the domain length in y direction
# Depth is the domain length in z direction
# numbering system of cells, like the single index numbering of Matlab
# +2 is added to account for the ghost cells that are added at
# the boundaries
dx = Width/Nx
dy = Height/Ny
dz = Depth/Nz
G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)
MeshStructure(3,
	[Nx, Ny, Nz],
	CellSize(dx*ones(Nx+2), dy*ones(Ny+2), dz*ones(Nz+2)),
	CellLocation([1:Nx;]*dx.-dx/2.0, [1:Ny;]*dy.-dy/2.0, [1:Nz;]*dz.-dz/2.0),
	FaceLocation([0:Nx;]*dx, [0:Ny;]*dy, [0:Nz;]*dz),
	G[[1,end],[1,end],[1,end]][:],
	[G[[1, end], [1, end], 2:Nz+1][:];
	G[[1, end], 2:Ny+1, [1, end]][:];
	G[2:Nx+1, [1, end], [1, end]][:]])
end

function createMesh3D(facelocationX::Array{T,1}, facelocationY::Array{T,1}, facelocationZ::Array{T,1}) where T<:Real
Nx = length(facelocationX)-1
Ny = length(facelocationY)-1
Nz = length(facelocationZ)-1
G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)
MeshStructure(3,
	[Nx, Ny, Nz],
	CellSize([facelocationX[2]-facelocationX[1]; facelocationX[2:end]-facelocationX[1:end-1]; facelocationX[end]-facelocationX[end-1]],
	  [facelocationY[2]-facelocationY[1]; facelocationY[2:end]-facelocationY[1:end-1]; facelocationY[end]-facelocationY[end-1]],
	  [facelocationZ[2]-facelocationZ[1]; facelocationZ[2:end]-facelocationZ[1:end-1]; facelocationZ[end]-facelocationZ[end-1]]),
	CellLocation(0.5*(facelocationX[2:end]+facelocationX[1:end-1]),
	  0.5*(facelocationY[2:end]+facelocationY[1:end-1]),
	  0.5*(facelocationZ[2:end]+facelocationZ[1:end-1])),
	FaceLocation(facelocationX, facelocationY, facelocationZ),
	G[[1,end],[1,end],[1,end]][:],
	[G[[1, end], [1, end], 2:Nz+1][:];
	G[[1, end], 2:Ny+1, [1, end]][:];
	G[2:Nx+1, [1, end], [1, end]][:]])
end

# ======================================= 3D CYLINDRICAL MESH ==============================================
function createMeshCylindrical3D(Nr::Int, Ntheta::Int, Nz::Int, Radius::Real, Angle::Real, Height::Real)
# builds a uniform 3D mesh:
# Nr is the number of cells in x direction
# Ntheta is the number of cells in y direction
# Nz is the number of cells in z direction
# Radius is the domain length in x direction
# Angle is the domain length in y direction
# Depth is the domain length in z direction

# numbering system of cells, like the single index numbering of Matlab
# +2 is added to account for the ghost cells that are added at
# the boundaries
if Angle>2*pi
	Angle = 2*pi
	println("The domain size adjusted to match a maximum of 2*pi.")
end
dr = Radius/Nr
dtheta = Angle/Ntheta
dz = Height/Nz
G=reshape([1:(Nr+2)*(Ntheta+2)*(Nz+2);], Nr+2, Ntheta+2, Nz+2)
MeshStructure(3.2,
	[Nr, Ntheta, Nz],
	CellSize(dr*ones(Nr+2), dtheta*ones(Ntheta+2), dz*ones(Nz+2)),
	CellLocation([1:Nr;]*dr.-dr/2.0, [1:Ntheta;]*dtheta.-dtheta/2.0, [1:Nz;]*dz.-dz/2.0),
	FaceLocation([0:Nr;]*dr, [0:Ntheta;]*dtheta, [0:Nz;]*dz),
	G[[1,end],[1,end],[1,end]][:],
	[G[[1, end], [1, end], 2:Nz+1][:];
	G[[1, end], 2:Ntheta+1, [1, end]][:];
	G[2:Nr+1, [1, end], [1, end]][:]])
end

function createMeshCylindrical3D(facelocationR::Array{T,1}, facelocationTheta::Array{T,1}, facelocationZ::Array{T,1}) where T<:Real
if facelocationTheta[end]>2*pi
	facelocationTheta = facelocationTheta/facelocationTheta[end]*2.0*pi
	println("The domain size adjusted to match a maximum of 2*pi.")
end
Nx = length(facelocationR)-1
Ny = length(facelocationTheta)-1
Nz = length(facelocationZ)-1
G=reshape([1:(Nx+2)*(Ny+2)*(Nz+2);], Nx+2, Ny+2, Nz+2)
MeshStructure(3.2,
	[Nx, Ny, Nz],
	CellSize([facelocationR[2]-facelocationR[1]; facelocationR[2:end]-facelocationR[1:end-1]; facelocationR[end]-facelocationR[end-1]],
	  [facelocationTheta[2]-facelocationTheta[1]; facelocationTheta[2:end]-facelocationTheta[1:end-1]; facelocationTheta[end]-facelocationTheta[end-1]],
	  [facelocationZ[2]-facelocationZ[1]; facelocationZ[2:end]-facelocationZ[1:end-1]; facelocationZ[end]-facelocationZ[end-1]]),
	CellLocation(0.5*(facelocationR[2:end]+facelocationR[1:end-1]),
	  0.5*(facelocationTheta[2:end]+facelocationTheta[1:end-1]),
	  0.5*(facelocationZ[2:end]+facelocationZ[1:end-1])),
	FaceLocation(facelocationR, facelocationTheta, facelocationZ),
	G[[1,end],[1,end],[1,end]][:],
	[G[[1, end], [1, end], 2:Nz+1][:];
	G[[1, end], 2:Ny+1, [1, end]][:];
	G[2:Nx+1, [1, end], [1, end]][:]])
end


text:
function EventState(arg0::EventDetector)
    return EventState((EventDetector,), arg0)
end

function do_event(obj::EventState, arg0::SpacecraftState)
    return jcall(obj, "doEvent", EventState_EventOccurrence, (SpacecraftState,), arg0)
end

function evaluate_step(obj::EventState, arg0::OrekitStepInterpolator)
    return jcall(obj, "evaluateStep", jboolean, (OrekitStepInterpolator,), arg0)
end

function get_event_date(obj::EventState)
    return jcall(obj, "getEventDate", AbsoluteDate, ())
end

function get_event_detector(obj::EventState)
    return jcall(obj, "getEventDetector", EventDetector, ())
end

function init(obj::EventState, arg0::SpacecraftState, arg1::AbsoluteDate)
    return jcall(obj, "init", void, (SpacecraftState, AbsoluteDate), arg0, arg1)
end

function reinitialize_begin(obj::EventState, arg0::OrekitStepInterpolator)
    return jcall(obj, "reinitializeBegin", void, (OrekitStepInterpolator,), arg0)
end

function try_advance(obj::EventState, arg0::SpacecraftState, arg1::OrekitStepInterpolator)
    return jcall(obj, "tryAdvance", jboolean, (SpacecraftState, OrekitStepInterpolator), arg0, arg1)
end


text:
__precompile__()
module NewsvendorModel

using Distributions
using QuadGK
using Printf

import Base: show

# export my_round
export AbstractNewsvendorProblem, underage_cost, overage_cost, profit_shift,
distr, q_opt, q_scarf, lost_sales, sales, leftover, profit, mismatch_cost,
solve, penalty, critical_fractile, q_max, q_min
export NVModel, rounded, nvmodel

"""
An abstract newsvendor problem is essentially described by having a
  - demand distribution,
  - cost of overage,
  - cost of underage, and
  - profit_shift term.
Its standard concrete type might come as a NewsvendorModel, defined by
unit cost, uni selling price etc.
"""
abstract type AbstractNewsvendorProblem end

include("./newsvendorproblem.jl")
include("./leftover.jl")
include("./nvmodel.jl")

end # module


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
name = "libass"
version = v"0.14.0"

# Collection of sources required to build libass
sources = [
    ArchiveSource("https://github.com/libass/libass/releases/download/$(version)/libass-$(version).tar.xz",
                  "881f2382af48aead75b7a0e02e65d88c5ebd369fe46bc77d9270a94aa8fd38a2"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/libass-*
apk add nasm
./configure --prefix=${prefix} --build=${MACHTYPE} --host=${target} --disable-require-system-font-provider --disable-static
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("FreeType2_jll"),
    Dependency("FriBidi_jll"),
    # Future versions of bzip2 should allow a more relaxed compat because the
    # soname of the macOS library shouldn't change at every patch release.
    Dependency("Bzip2_jll", v"1.0.6"; compat="=1.0.6"),
    Dependency("Zlib_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
using PMDPs

function linear_params(objective)
    [
        (params = Dict(pairs((nᵣ=3, c=1, T=10, expected_res=6., res_budget_μ=5.,
            objective=objective))),
         vi = true, name="linear_problem"),
        (params = Dict(pairs((nᵣ=5, c=2, T=40, expected_res=30., res_budget_μ=5.,
            objective=objective))),
         vi = true, name="linear_problem"), # last VI size to run
        (params = Dict(pairs((nᵣ=10, c=5, T=100, expected_res=100., res_budget_μ=5.,
            objective=objective))),
         vi = false, name="linear_problem"),
        (params = Dict(pairs((nᵣ=10, c=40, T=1000, expected_res=Float64(800), res_budget_μ=5.,
            objective=objective))),
         vi = false, name="linear_problem"),
        # Dict(pairs((nᵣ=50, c=40, T=1000, expected_res=Float64(4000), res_budget_μ=5.,
        #     objective=objective)))
    ]
end

function graph_params(objective)
    [
        (params = Dict(pairs((NV=5, NE=8, seed=1, NP=20, c=5, T=100, expected_res=Float64(80),
            res_budget_μ=5., objective=objective))),
         vi = false, name="graph_problem"),
        (params = Dict(pairs((NV=8, NE=20, seed=1, NP=50, c=10, T=1000, expected_res=Float64(400),
            res_budget_μ=5., objective=objective))),
         vi = false, name="graph_problem"),
        (params = Dict(pairs((NV=15, NE=30, seed=1, NP=100, c=10, T=1000, expected_res=Float64(600),
            res_budget_μ=5., objective=objective))),
         vi = false, name="graph_problem"),
        # Dict(pairs((NV=30, NE=45, seed=1, NP=100, c=10, T=1000, expected_res=Float64(900),
        #     res_budget_μ=5., objective=objective)))
    ]
end

function get_linear_problems(objective::Symbol)
    pps = [(pp=PMDPs.linear_pp(;p.params...), p...)
        for p in linear_params(objective)]
end
function get_graph_problems(objective::Symbol)
    pps = [(pp=PMDPs.graph_pp(;p.params...), p...)
        for p in graph_params(objective)]
end

"""
n_lp ... how many lp problems to use (e.g. 3 - use the first three problems)
n_gp ... how many gp problems to use (e.g. 3 - use the first three problems).
By default, first 100 problems will be used.
"""
function get_benchmarks(;objectives=(:revenue, :utilization), n_lp=100, n_gp=100)
    problems=[]
    for obj in objectives
        lpps = get_linear_problems(obj)
        gpps = get_graph_problems(obj)
        push!(problems, lpps[1:minimum([n_lp, length(lpps)])]...)
        push!(problems, gpps[1:minimum([n_gp, length(gpps)])]...)
    end
    problems
end

function get_fast_benchmarks()
    get_benchmarks(; n_lp=3, n_gp=0)
end

function get_tiny_benchmarks()
    get_benchmarks(; n_lp=2, n_gp=1)
end

# linear_params(:revenue)
# graph_params(:revenue)

# get_linear_problems(:revenue)
# get_graph_problems(:revenue)

# get_benchmarks()
# fb = get_fast_benchmarks()

# fb[1].params

# get_linear_problems(:revenue)
# get_graph_problems(:revenue)


text:
module TestMetadata

using Test
using MLJModels
import MLJBase

metadata_file = joinpath(@__DIR__, "..", "src",
                         "registry", "Metadata.toml")
pca = MLJModels.Handle("PCA", "MultivariateStats")
cnst = MLJModels.Handle("ConstantRegressor", "MLJModels")
i = MLJModels.info_given_handle(metadata_file)[pca]

@testset "Handle constructors" begin
    @test MLJModels.Handle("PCA") ==
        MLJModels.Handle("PCA", "MultivariateStats")
    # TODO: add tests here when duplicate model names enter registry
end

@testset "building INFO_GIVEN_HANDLE" begin
    @test isempty(MLJModels.localmodeltypes(MLJBase))
    @test issubset(Set([DeterministicConstantClassifier,
                        DeterministicConstantRegressor,
                        ConstantClassifier,
                        ConstantRegressor,
                        FeatureSelector,
                        OneHotEncoder,
                        Standardizer,
                        UnivariateBoxCoxTransformer,
                        UnivariateStandardizer]),
                   MLJModels.localmodeltypes(MLJModels))
    @test MLJModels.info_given_handle(metadata_file)[pca][:name] == "PCA"
    d1 = MLJModels.info_given_handle(metadata_file)[cnst]
    d2 = MLJBase.info_dict(ConstantRegressor)
    for (k, v) in d1
        if v isa Vector
            @test Set(v) == Set(d2[k])
        else
            @test v == d2[k]
        end
    end
end

h = Vector{Any}(undef, 7)
h[1] = MLJModels.Handle("1", "a")
h[3] = MLJModels.Handle("1", "b")
h[2] = MLJModels.Handle("2", "b")
h[4] = MLJModels.Handle("3", "b")
h[5] = MLJModels.Handle("4", "c")
h[6] = MLJModels.Handle("3", "d")
h[7] = MLJModels.Handle("5", "e")
info_given_handle = Dict([h[j]=>i for j in 1:7]...)

@testset "building AMBIGUOUS_NAMES" begin
    @test Set(MLJModels.ambiguous_names(info_given_handle)) == Set(["1", "3"])
end

@testset "building PKGS_GIVEN_NAME" begin
    d = MLJModels.pkgs_given_name(info_given_handle)
    @test Set(d["1"]) == Set(["a", "b"])
    @test d["2"]==["b",]
    @test Set(d["3"]) == Set(["b", "d"])
    @test d["4"] == ["c",]
    @test d["5"] == ["e",]
end

@testset "building NAMES" begin
    model_names = MLJModels.model_names(info_given_handle)
    @test Set(model_names) == Set(["1", "2", "3", "4", "5"])
end

end
true


text:
include("D_gain.jl")
include("D_inv_gain.jl")
include("A.jl")
include("A_gain.jl")
include("A_inv_gain.jl")


text:
# https://gist.github.com/silgon/0ba43e00e0749cdf4f8d244e67cd9d6a
# https://riptutorial.com/julia-lang/example/19443/parsing-json
# https://github.com/JuliaIO/JSON.jl

import Pkg;
Pkg.add("JSON")

import JSON
###################
### Write data ####
###################
# dictionary to write
dict1 = Dict(
                "param1" => 1,
                "param2" => 2,
                "dict" => Dict(
                                "d1" => 1.,
                                "d2" => 1.,
                                "d3" => 1.
                            )
            )

# pass data as a json string (how it shall be displayed in a file)
stringdata = JSON.json(dict1)

# write the file with the stringdata variable information
open("write_read.json", "w") do f
    write(f, stringdata)
end

###################
### Read data #####
###################
# create variable to write the information
dict2 = Dict()
open("write_read.json", "r") do f
    global dict2
    dicttxt = read(f, String)  # file information to string
    dict2 = JSON.parse(dicttxt)  # parse and transform data
end

# print both dictionaries
println(dict1)
println(dict2)


text:
module MoneyMarkets

tBillR(t, P0, F) = (1.0 - P0/F) * 360.0/t
tBillD(r, t, F) = r*t*F/360.0

holdingPerYield(P0, P1, D1) = (P1+D1)/P0 - 1.0
effAnnYield(t, P0, P1, D1) = ((P1+D1)/P0)^(365/t) - 1.0
moneyMktYield(t, P0, P1, D1) = ((P1+D1)/P0 - 1.0)*360/t

function twrr(bv, b_inf)
  n = length(bv)-1; r = 1.0
  for i in 1:n r *= bv[i+1]/(bv[i]+b_inf[i]) end
  r^(1/n) - 1.0
end

function twrr(n, bv, b_inf)
  r = 1.0
  for i in 1:(length(bv)-1) r *= bv[i+1]/(bv[i]+b_inf[i]) end
  r^(1/n) - 1.0
end

end


text:
import Pkg

@static if VERSION >= v"1.4.0"
    function _get_package_version(name)
        for (uuid, dep) in Pkg.dependencies()
            if dep.name == name
                return dep.version
            end
        end
        return nothing
    end
else
    function _get_package_version(name)
        return get(Pkg.installed(), name, nothing)
    end
end

const MLJ_VERSION = _get_package_version("MLJ")


text:
# ---
# title: Flip
# cover: flip.gif
# description: flip the input image horizontally or vertically
# ---

# [`FlipX`](@ref)/[`FlipY`](@ref) can be used to flip the input image horizontally/vertically.

using Augmentor
using ImageShow, ImageCore

img_in = testpattern(RGB, ratio=0.5)

mosaicview(
    img_in,
    augment(img_in, FlipX()),
    augment(img_in, FlipY());
    fillvalue=colorant"white", nrow=1, npad=10
)

# To perform a random flip, you can also pass the probablity to the constructor. For example, `FlipX(0.5)`
# flips the image with half chance.

# ## References

#md # ```@docs
#md # FlipX
#md # FlipY
#md # ```

## save covers #src
using ImageMagick #src
using FileIO #src
include(joinpath("..", "assets", "utilities.jl")) #src
cover = make_gif(testpattern(RGB, ratio=0.5), Either(FlipX(), FlipY()), 4) #src
ImageMagick.save("flip.gif", cover; fps=1) #src


text:
using Plots
using Wavelets
WT.paul10
supertype(typeof(WT.paul10))
wavelet(WT.morl)
wavelet(WT.morl,4)
plotlyjs()
J=11; n = 2^J
x = testfunction(n,"Bumps")
supertype(typeof(WT.dog0))
supertype(typeof(WT.paul5))
supertype(typeof(WT.morl))
EmphasizeFrequencyInfo = WT.Morlet(20)
y = cwt(x, wavelet(WT.morl))
heatmap(abs.(y)); plot!(x+90,label="")
y = cwt(x, wavelet(WT.dog1))
heatmap(abs.(y)); plot!(x+90,label="")
y = cwt(x, wavelet(WT.paul4))
heatmap(abs.(y)); plot!(x+90,label="")


text:
using BinDeps
using Compat
@BinDeps.setup

version = "1.2.0"
url = "https://github.com/liberfa/erfa/releases/download/v$version/erfa-$version.tar.gz"

# This function returns true if a library satisfies our
# requirements. BinDeps passes to this function the library path and
# the dlopen'ed handle. We test is certain symbols we need exist.
validate(l, h) = (Libdl.dlsym_e(h, "eraAb") != C_NULL &&
                  Libdl.dlsym_e(h, "eraG2icrs") != C_NULL &&
                  Libdl.dlsym_e(h, "eraIcrs2g") != C_NULL)

erfa = library_dependency("liberfa"; validate = validate)

provides(Sources, URI(url), erfa)
provides(BuildProcess, Autotools(libtarget="src/liberfa.la"), erfa)

# Windows binaries cross-compiled in erfa with
#
# mkdir -p /usr/lib
# i686-w64-mingw32-gcc -Isrc -o usr/lib/liberfa.dll -O3 -shared \
#     -static-libgcc src/*.c
# 7za a erfa-win32.7z usr
# x86_64-w64-mingw32-gcc -Isrc -o usr/lib/liberfa.dll -O3 -shared \
#     -static-libgcc src/*.c
# 7za a erfa-win64.7z usr

provides(Binaries, URI("https://bintray.com/artifact/download/kbarbary/generic/erfa-win$(Sys.WORD_SIZE).7z"), erfa, os = :Windows)

@BinDeps.install @compat Dict(:liberfa => :liberfa)


text: include("mlogloss.jl")


text:
struct Angle

   i :: Int64
   j :: Int64
   k :: Int64

   ktheta :: Float64
   theta0 :: Float64

   kub :: Float64
   s0 :: Float64

end


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

using JuliaFEM, Test

meshfile = "/geometry/2d_block/BLOCK_1elem.med"
mesh = aster_read_mesh(dirname(@__DIR__)*meshfile)

# field problem
block = Problem(Elasticity, "BLOCK", 2)
block.properties.formulation = :plane_strain
block.properties.finite_strain = false
block.properties.geometric_stiffness = false

block_elements = create_elements(mesh, "BLOCK")
update!(block_elements, "youngs modulus", 288.0)
update!(block_elements, "poissons ratio", 1/3)
traction_elements = create_elements(mesh, "TOP")
update!(traction_elements, "displacement traction force 2", 288.0*9/8)
add_elements!(block, block_elements, traction_elements)

# boundary conditions
bc = Problem(Dirichlet, "symmetry bc 23", 2, "displacement")
bc_sym_23_elements = create_elements(mesh, "LEFT")
bc_sym_13_elements = create_elements(mesh, "BOTTOM")
update!(bc_sym_23_elements, "displacement 1", 0.0)
update!(bc_sym_13_elements, "displacement 2", 0.0)
add_elements!(bc, bc_sym_23_elements, bc_sym_13_elements)

analysis = Analysis(Linear, block, bc)
run!(analysis)

u = block("displacement", 0.0)
u3_expected = [-0.5, 1.0]
@debug("displacement", u, u3_expected)
@test isapprox(u[3], u3_expected)


text:
function exit_reward(e::Episode; reward::Float64=1.0)
    linds = LinearIndices(e.grid.cells)
    if linds[e.robot.y, e.robot.x] in e.grid.exits
        return reward
    end
    0.0
end

function food_reward(e::Episode; reward::Float64=1.0)
    if object(e.grid.cells[e.robot.y, e.robot.x]) == "food"
        e.grid.cells[e.robot.y, e.robot.x].obj = type_to_float("empty")
        return reward
    end
    0.0
end

function coverage_reward(e::Episode; reward::Float64=1.0)
    if ~("visited" in keys(e.meta))
        e.meta["visited"] = falses(size(e.grid.cells))
    end
    if ~e.meta["visited"][e.robot.y, e.robot.x]
        e.meta["visited"][e.robot.y, e.robot.x] = true
        return reward
    end
    0.0
end


text:
export allele_avg, richness, summary, summary_stats

"""
    allele_avg(data::PopData; rounding::Bool = true)
Returns a NamedTuple of the average number of alleles ('mean') and
standard deviation (`stdev`) of a `PopData`. Use `rounding = false` to
not round results. Default (`true`) rounds to 4 digits.
"""
function allele_avg(data::PopData; rounding::Bool = true, population = false)
    tmp = richness(data)
    if rounding
        (mean = round(mean(tmp.richness), digits = 4), stdev = round(variation(tmp.richness), digits = 4))
    else
        (mean = mean(tmp.richness), stdev = variation(tmp.richness))
    end
end

"""
    richness(data::PopData; by::String = "locus")
Calculates various allelic richness and returns a table of per-locus
allelic richness. Use `by = "population"` to calculate richness by
locus by population.
"""
function richness(data::PopData; by::String = "locus")
    if by == "locus"
        DataFrames.combine(
            groupby(data.loci, :locus),
            :genotype => (geno -> length(unique_alleles(geno))) => :richness
        )
    elseif by == "population"
        DataFrames.combine(
            groupby(data.loci, [:locus, :population]),
            :genotype => (geno -> length(unique_alleles(geno))) => :richness
        )
    else
        throw(ArgumentError("Please use by = \"locus\" (default) or \"population\""))
    end
end

#TODO add citations to docstring
"""
    summary(data::PopData; by::String = "global")
Provides summary statistics for a `PopData` object. Use `by = "locus"` for
summary information by locus. Global values are given as unweighted means of
the per-locus parameters.

### Het_obs
observed heterozygosity given as:\n
1 - ∑ₖ ∑ᵢ Pₖᵢᵢ/np \n
where Pkii represents the proportion of homozygote `i` in sample `k` and `np`
is the number of samples in that population

### HT
overall gene diversity given as: \n
1 - ∑ᵢ(p̄ᵢ² + (HS / (ñ × np)) - Het_obs / (2 × ñ × np)) \n
where p̄ᵢ = ∑ₖpₖᵢ / np

### HS
within population gene diversity given as: \n
1 - ∑ᵢ(pᵢ² + HS / (ñ × np) - Het_obs / (2 × ñ × np)) \n
where ñ = np / ∑ₖ(1/nₖ) \n
where p̄ᵢ² = ∑ₖ(pᵢₖ² / np)

### DST
amount of gene diversity among samples given as: \n
HT - HS

### DST′
amount of gene diversity among samples adjusted for sample size given as: \n
(np / (np-1)) × Dst

### HT′
overall gene diversity adjusted for sample size given as: \n
HS + DST′

### FST
proportion of the total genetic variance in subpopulations relative to the total genetic variance  given as: \n
DST / HT

### FST′
proportion of the total genetic variance in subpopulations relative to the total genetic variance, adjusted for heterozygosity given as: \n
DST′ / HT′

### FIS
proportion of the genetic variance in a locus relative to the genetic variance within subpopulations given as: \n
1 - (Het_obs / HS)

### DEST
population differentiation (Jost 2008) given as: \n
(np/(np-1)) × (Ht'-Hs) / (1-Hs)
"""
function Base.summary(data::PopData; by::String = "global")
    # observed/expected het per locus per pop
    het_df = DataFrames.combine(
        groupby(data.loci, [:locus, :population]),
        :genotype => nonmissing => :n,
        :genotype => hetero_o => :het_obs,
        :genotype => (i -> hetero_e(i)) => :het_exp,
        :genotype => allele_freq => :alleles
    )
    # collapse down to retrieve averages and counts
    n_df = DataFrames.combine(
        groupby(het_df, :locus),
        :n => (n -> sum(.!iszero.(n)))=> :count,
        :n => (n -> sum(.!iszero.(n)) ./ sum(reciprocal.(n))) => :mn,
        [:het_obs, :het_exp, :n] => ((o,e,n) -> mean(skipmissing(gene_diversity_nei87.(e,o,sum(.!iszero.(n)) ./ sum(reciprocal.(n)))))) => :HS,
        :het_obs => (o -> mean(skipmissing(o)))=> :Het_obs,
        :alleles => (alleles ->  sum(values(avg_allele_freq(alleles, 2))))=> :avg_freq
        )

    Ht = 1.0 .- n_df.avg_freq .+ (n_df.HS ./ n_df.mn ./ n_df.count) - (n_df.Het_obs ./ 2.0 ./ n_df.mn ./ n_df.count)
    DST = Ht .- n_df.HS
    DST′ = n_df.count ./ (n_df.count .- 1) .* DST
    HT′ = n_df.HS .+ DST′

    if lowercase(by) == "locus"
        FIS =  1.0 .- (n_df.Het_obs ./ n_df.HS)
        FST = DST ./ Ht
        DEST = DST′ ./ (1 .- n_df.HS)
        FST′ = DST′ ./ HT′
        insertcols!(
            select!(
                n_df,
                :locus,
                :Het_obs => (i -> round.(i, digits = 4)) => :Het_obs,
                :HS => (i -> round.(i, digits = 4)) => :HS,
            ),
            :HT => round.(Ht, digits = 4),
            :DST => round.(DST, digits = 4),
            :HT′ => round.(HT′, digits = 4),
            :DST′ => round.(DST′, digits = 4),
            :FST => round.(FST, digits = 4),
            :FST′ => round.(FST′, digits = 4),
            :FIS => round.(FIS, digits = 4),
            :DEST => round.(DEST, digits = 4)
        )
    elseif lowercase(by) == "global"
        Het_obs = mean(n_df.Het_obs)
        HS = mean(n_df.HS)
        Ht = mean(Ht)
        DST = mean(DST)
        DST′ = mean(DST′)
        HT′ = mean(HT′)

        DataFrame(
        :Het_obs => round.(Het_obs, digits = 4),
        :HS => round.(HS, digits = 4),
        :HT => round.(Ht , digits = 4),
        :DST => round.(DST, digits = 4),
        :HT′ => round.(HT′, digits = 4),
        :DST′ => round.(DST′, digits = 4),
        :FST => round(DST / Ht, digits = 4),
        :FST′ => round(DST′ / HT′, digits = 4),
        :FIS => round(1 - (Het_obs / HS), digits = 4),
        :DEST => round(DST′ / (1 - HS), digits = 4)
        )
    else
        throw(ArgumentError("Use by = \"global\" or \"locus\""))
    end
end

function _summary(data::PopData; by::String = "global")
    # observed/expected het per locus per pop
    het_df = DataFrames.combine(
        groupby(data.loci, [:locus, :population]),
        :genotype => nonmissing => :n,
        :genotype => hetero_o => :het_obs,
        :genotype => (i -> hetero_e(i)) => :het_exp,
        :genotype => allele_freq => :alleles
    )
    # collapse down to retrieve averages and counts
    n_df = DataFrames.combine(
        groupby(het_df, :locus),
        :n => (n -> sum(.!iszero.(n)))=> :count,
        :n => (n -> sum(.!iszero.(n)) ./ sum(reciprocal.(n))) => :mn,
        [:het_obs, :het_exp, :n] => ((o,e,n) -> mean(skipmissing(gene_diversity_nei87.(e,o,sum(.!iszero.(n)) ./ sum(reciprocal.(n)))))) => :HS,
        :het_obs => (o -> mean(skipmissing(o)))=> :Het_obs,
        :alleles => (alleles ->  sum(values(avg_allele_freq(alleles)).^2))=> :avg_freq
    )

    Ht = 1.0 .- n_df.avg_freq .+ (n_df.HS ./ n_df.mn ./ n_df.count) - (n_df.Het_obs ./ 2.0 ./ n_df.mn ./ n_df.count)
    DST = Ht .- n_df.HS
    DST′ = n_df.count ./ (n_df.count .- 1) .* DST
    HT′ = n_df.HS .+ DST′

    if lowercase(by) != "global"
        FIS =  1.0 .- (n_df.Het_obs ./ n_df.HS)
        FST = DST ./ Ht
        DEST = DST′ ./ (1 .- n_df.HS)
        FST′ = DST′ ./ HT′
        DataFrame(
            :locus => n_df.locus,
            :Het_obs => n_df.Het_obs,
            :HS => n_df.HS,
            :HT => Ht,
            :DST => DST,
            :HT′ => HT′,
            :DST′ =>DST′,
            :FST => FST,
            :FST′ => FST′,
            :FIS => FIS,
            :DEST => DEST
        )
    else
        Het_obs = mean(n_df.Het_obs)
        HS = mean(n_df.HS)
        Ht = mean(Ht)
        DST = mean(DST)
        DST′ = mean(DST′)
        HT′ = mean(HT′)

        DataFrame(
        :Het_obs => Het_obs,
        :HS => HS,
        :HT => Ht,
        :DST => DST,
        :HT′ => HT′,
        :DST′ => DST′,
        :FST => DST / Ht,
        :FST′ => DST′ / HT′,
        :FIS => 1 - (Het_obs / HS),
        :DEST => DST′ / (1 - HS)
        )
    end
end

const summary_stats = summary


text:
struct GenericStorage{T,scalar} <: AbstractStorage{T, scalar}
  data::Array{T,4}
  norb::Int
  N::Int
  M::Int

  function GenericStorage{T, scalar}(data::AbstractArray{T,4}) where {T, scalar}
    norb = size(data,1)
    N = size(data,3)
    M = size(data,4)
    @assert size(data) == (norb, norb, N, M)
    @assert !(scalar && norb > 1)
    return new(data, norb, N, M)
  end
end

function GenericStorage(data::AbstractArray{T,4}, scalar=false) where T
  GenericStorage{T,scalar}(data)
end

function GenericStorage(::Type{T}, N::Integer, M::Integer, norb=1, scalar=false) where T <: Number
  data = zeros(T, norb, norb, N, M)
  return GenericStorage(data, scalar)
end

GenericStorage(N::Integer, M::Integer, norb=1, scalar=false) = GenericStorage(ComplexF64, N, M, norb, scalar)

@inline function Base.getindex(X::GenericStorage{T,scalar}, k, l, i, j) where {T, scalar}
  return X.data[k,l,i,j]
end

@inline function Base.setindex!(X::GenericStorage{T,scalar}, v, k, l, i, j) where {T, scalar}
  return X.data[k,l,i,j] = v
end


text:
module EntropyM

export entropy

# include("Utils.jl")

using Main.UtilsM: removeExtras

function entropy(text::String, alphabet::String)

    refined = removeExtras(text, alphabet)
    return sum([
        (-i) * log2(i)
        for i in [
            count(string(letter), refined) / length(refined)
            for letter in alphabet if occursin(letter, text)
        ]
    ])

end # function

end # module

max_stars_repo_name: Squalm/Cipher-Tools

max_issues_repo_name: Squalm/Cipher-Tools

max_forks_repo_name: Squalm/Cipher-Tools


text: # generate training data for ML model used in Deka, Misra


text:
domain_type(::Type{<:AbstractPICDataStructure{T,N,G}}) where {T,N,G} = parameterless_type(G)

function domain_discretization(::Type{<:AbstractPICDataStructure{T,N,G}}) where {T,N,G}
    domain_discretization(G)
end

domain_discretization(::Type{<:AbstractAxisGrid{N}}) where N = LatticeGrid{N}()
domain_discretization(::Type{<:ParticlePositions}) = ParticleGrid()


text:
using PyFenicsAD
using FiniteDifferences
using Zygote
using PyCall
using Test

fenics = pyimport("fenics")
fa = pyimport("fenics_adjoint")
ufl = pyimport("ufl")

mesh = fa.UnitSquareMesh(3, 2)
V = fenics.FunctionSpace(mesh, "P", 1)

function assemble_fenics(u, kappa0, kappa1)

    f = fa.Expression(
        "10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)",
        degree = 2,
    )

    inner, grad, dx = ufl.inner, ufl.grad, ufl.dx
    J_form = 0.5 * inner(kappa0 * grad(u), grad(u)) * dx - kappa1 * f * u * dx
    # J = fa.assemble(J_form)
    # Avoid PyCall's automatic conversion to Julia types
    J = pycall(fa.assemble, PyObject, J_form)
    return J
end

templates = (fa.Function(V), fa.Constant(0.0), fa.Constant(0.0))
inputs = (ones(V.dim()), ones(1) * 0.5, ones(1) * 0.6)

@register_fem_function(zygote_assemble_fenics, templates, assemble_fenics)

hh(inputs...) = zygote_assemble_fenics(inputs...)
hh0(x) = hh(x, inputs[2], inputs[3])
hh1(y) = hh(inputs[1], y, inputs[3])
hh2(z) = hh(inputs[1], inputs[2], z)

@testset "chainrules_zygote_assemble_primal" begin
    out = zygote_assemble_fenics(inputs...)

    u1 = fa.interpolate(fa.Constant(1.0), V)
    J = assemble_fenics(u1, fa.Constant(0.5), fa.Constant(0.6))
    @test out[1] == PyCall.convert(Float64, J)
end

@testset "chainrules_zygote_assemble_pullback" begin
    out, vjp_fun = Zygote.pullback(hh, inputs...)

    g = ones(size(out))
    vjp_out = vjp_fun(g)

    fdm = FiniteDifferences.central_fdm(2, 1)

    fdm_jac1 = FiniteDifferences.jacobian(fdm, hh0, inputs[1])
    fdm_jac2 = FiniteDifferences.jacobian(fdm, hh1, inputs[2])
    fdm_jac3 = FiniteDifferences.jacobian(fdm, hh2, inputs[3])

    @test isapprox(vjp_out[1], fdm_jac1[1]', atol = 1e-5)
    @test isapprox(vjp_out[2], fdm_jac2[1]', atol = 1e-5)
    @test isapprox(vjp_out[3], fdm_jac3[1]', atol = 1e-5)
end

mesh = fa.UnitSquareMesh(6, 5)
V = fenics.FunctionSpace(mesh, "P", 1)

function solve_fenics(kappa0, kappa1)

    f = fa.Expression(
        "10*exp(-(pow(x[0] - 0.5, 2) + pow(x[1] - 0.5, 2)) / 0.02)",
        degree = 2,
    )

    u = fa.Function(V)
    bcs = [fa.DirichletBC(V, fa.Constant(0.0), "on_boundary")]

    inner, grad, dx = ufl.inner, ufl.grad, ufl.dx
    JJ = 0.5 * inner(kappa0 * grad(u), grad(u)) * dx - kappa1 * f * u * dx
    v = fenics.TestFunction(V)
    F = fenics.derivative(JJ, u, v)
    fa.solve(F == 0, u, bcs = bcs)
    return u
end

templates = (fa.Constant(0.0), fa.Constant(0.0))
inputs = (ones(1) * 0.5, ones(1) * 0.6)

@register_fem_function(zygote_solve_fenics, templates, solve_fenics)

@testset "chainrules_zygote_solve_primal" begin
    out = zygote_solve_fenics(inputs...)
    u = solve_fenics(fa.Constant(0.5), fa.Constant(0.6))
    @test isapprox(out, to_numpy(u))
end

@testset "chainrules_zygote_solve_pullback" begin
    out, vjp_fun = Zygote.pullback(zygote_solve_fenics, inputs...)
    # g = np.ones_like(numpy_output)
    g = ones(size(out))
    vjp_out = vjp_fun(g)
    @test isapprox(vjp_out[1], [-2.91792642])
    @test isapprox(vjp_out[2], [2.43160535])
end


text:
# This file was generated, do not modify it. # hide
tuned_iterated_pipe = TunedModel(model=iterated_pipe,
                                 range=[r1, r2],
                                 tuning=tuning,
                                 measures=[brier_loss, auc, accuracy],
                                 resampling=StratifiedCV(nfolds=6, rng=123),
                                 acceleration=CPUThreads(),
                                 n=40)

max_stars_repo_name: alan-turing-institute/MLJTutorials

max_issues_repo_name: alan-turing-institute/MLJTutorials

max_forks_repo_name: alan-turing-institute/MLJTutorials


text:
"""
    struct SignSymmetry <: Sparsity end

Sign symmetry as developed in [Section III.C, L09].

[L09] Lofberg, Johan.
*Pre-and post-processing sum-of-squares programs in practice*.
IEEE transactions on automatic control 54, no. 5 (2009): 1007-1011.
"""
struct SignSymmetry <: Sparsity end

function binary_exponent(exponents, ::Type{T}) where T
    cur = zero(T)
    for exponent in exponents
        cur <<= 1
        if isodd(exponent)
            cur |= one(T)
        end
    end
    return cur
end
bin_dot(x, y) = isodd(count_ones(x & y))
function buckets_sign_symmetry(monos, r, ::Type{T}, ::Type{U}) where {T, U}
    # We store vector of indices instead of vector of monos in the bucket
    # as `monos` is already sorted and the buckets content will be sorted
    # so it will remain sorted and `DynamicPolynomials` can keep the
    # `MonomialVector` struct in `monos[I]`.
    buckets = Dict{U, Vector{Int}}()
    for (idx, mono) in enumerate(monos)
        exp = binary_exponent(MP.exponents(mono), T)
        i = 1 + binary_exponent([bin_dot(ri, exp) for ri in r], U)
        if !haskey(buckets, i)
            buckets[i] = Int[]
        end
        push!(buckets[i], idx)
    end
    return [monos[I] for I in values(buckets)]
end
function sign_symmetry(monos::AbstractVector{<:MP.AbstractMonomial}, n, ::Type{T}, gram_monos) where T
    r = xor_complement((binary_exponent(MP.exponents(mono), T) for mono in monos), n, T)
    return buckets_sign_symmetry(gram_monos,
                                 r, T, appropriate_type(length(r)))
end
function sparsity(monos::AbstractVector{<:MP.AbstractMonomial}, sp::SignSymmetry,
                  gram_monos::AbstractVector = SumOfSquares.Certificate.monomials_half_newton_polytope(monos, tuple()))
    n = MP.nvariables(monos)
    return sign_symmetry(monos, n, appropriate_type(n), gram_monos)
end


text:
function create_empty_model()
    sc = SupplyChain()

    return sc
end

@testset "Happy Path" begin
    @test haversine(0, 0, 0, 0) == 0

    @test haversine(51.510357, -0.116773, 38.889931, -77.009003) ≈ 5897658.289

    @test add_lane!(SupplyChain(),
                             Lane(Customer("c1", Seattle),
                                  Customer("c2", Seattle);
                                  unit_cost=1.0,
                                  minimum_quantity=0,
                                  time=0)) isa Lane

    @test !isnothing(create_empty_model())
end


text:
#!/usr/bin/env julia
using Documenter
using HypertextLiteral

# Setup for doctests embedded in docstrings.
DocMeta.setdocmeta!(HypertextLiteral, :DocTestSetup, :(using HypertextLiteral))

# Highlight indented code blocks as Julia code.
using Documenter.Expanders: ExpanderPipeline, Selectors, Markdown, iscode
abstract type DefaultLanguage <: ExpanderPipeline end
Selectors.order(::Type{DefaultLanguage}) = 99.0
Selectors.matcher(::Type{DefaultLanguage}, node, page, doc) =
    iscode(node, "")
Selectors.runner(::Type{DefaultLanguage}, node, page, doc) =
    page.mapping[node] = Markdown.Code("julia", node.code)

makedocs(
    sitename = "HypertextLiteral.jl",
    format = Documenter.HTML(prettyurls=(get(ENV, "CI", nothing) == "true")),
    pages = [
        "Usage" => "index.md",
        "Extend" => "extend.md",
        "Design" => "design.md",
        "Notation" => "notation.md",
        "Primitives" => "primitives.md",
        "Reference" => "reference.md",
    ],
    modules = [HypertextLiteral])

deploydocs(
    repo = "github.com/clarkevans/HypertextLiteral.jl.git",
)


text:
# tests performance compared with discrete value iteration
using MCTS
using POMDPModels
using DiscreteValueIteration
using POMDPs

mdp = LegacyGridWorld()
N = 50
rng = MersenneTwister(1)
init_states = [GridWorldState(rand(rng, 1:mdp.size_x), rand(rng, 1:mdp.size_y)) for i in 1:N]

solvers = Dict(:vi => ValueIterationSolver(),
               :mcts1k => MCTSSolver(n_iterations=1000,
                                   depth=20,
                                   exploration_constant=10.0),
               :mcts10k => MCTSSolver(n_iterations=10000,
                                   depth=20,
                                   exploration_constant=10.0),
               :dpw => DPWSolver(n_iterations=1000,
                                 depth=20,
                                 exploration_constant=10.0,
                                 k_action=100.0, # these are all just set to large values so that behavior should be equivalent to MCTS
                                 alpha_action=1.0,
                                 k_state=100.0,
                                 alpha_action=1.0,
                                 ))

policies = Dict([(k, solve(s,mdp)) for (k,s) in solvers])
rewards = SharedArray(Float64, length(solvers), N)
index = Dict([(k, i) for (i, k) in enumerate(keys(solvers))])

for (k,p) in policies
    println("simulating $k")
    # @time @sync @parallel for i in 1:N
    @time for i in 1:N
        sim = MCTS.RolloutSimulator(rng=MersenneTwister(i))
        rewards[index[k],i] = simulate(sim, mdp, deepcopy(p), init_states[i])
    end
end

r = sdata(rewards)
for (k,i) in index
    println("$k: $(mean(r[i,:]))")
end


text:
function pull_request_build(data::GitHubAutoMergeData, ::NewVersion; check_license)::Nothing
    # Rules:
    # 1. Registry consistency tests pass
    # 2. Only changes a subset of the following files:
    #     - `E/Example/Compat.toml`
    #     - `E/Example/Deps.toml`
    #     - `E/Example/Versions.toml`
    # 3. Sequential version number
    #     - if the last version was 1.2.3 then the next can be 1.2.4, 1.3.0 or 2.0.0
    #     - does not apply to JLL packages
    # 4. Compat for all dependencies
    #     - there should be a [compat] entry for Julia
    #     - all [deps] should also have [compat] entries
    #     - all [compat] entries should have upper bounds
    #     - dependencies that are standard libraries do not need [compat] entries
    #     - dependencies that are JLL packages do not need [compat] entries
    # 5. If it is a patch release, then it does not narrow the Julia compat range
    # 6. (only applies to JLL packages) The only dependencies of the package are:
    #     - Pkg
    #     - Libdl
    #     - other JLL packages
    # 7. Version can be installed
    #     - given the proposed changes to the registry, can we resolve and install the new version of the package?
    #     - i.e. can we run `Pkg.add("Foo")`
    # 8. Package repository contains only OSI-approved license(s) in the license file at toplevel in the version being registered
    # 9. Version can be loaded
    #     - once it's been installed (and built?), can we load the code?
    #     - i.e. can we run `import Foo`

    this_is_jll_package = is_jll_name(data.pkg)
    @info("This is a new package pull request",
          pkg = data.pkg,
          version = data.version,
          this_is_jll_package)

    update_status(data;
                  state = "pending",
                  context = "automerge/decision",
                  description = "New version. Pending.")

    this_pr_can_use_special_jll_exceptions =
        this_is_jll_package && data.authorization == :jll

    # Each element is a tuple of a guideline and whether it's
    # applicable. Instead of a guideline there can be the symbol
    # `:update_status` in which case the PR status will be updated
    # with the results so far before continuing to the following
    # guidelines.
    guidelines =
        [
            (guideline_registry_consistency_tests_pass, true),     # 1
            (guideline_pr_only_changes_allowed_files, true),       # 2
            (guideline_sequential_version_number,                  # 3
                !this_pr_can_use_special_jll_exceptions),
            (guideline_compat_for_all_deps, true),                 # 4
            (guideline_patch_release_does_not_narrow_julia_compat, # 5
                !this_pr_can_use_special_jll_exceptions),
            (guideline_allowed_jll_nonrecursive_dependencies,      # 6
                this_is_jll_package),
            (:update_status, true),
            (guideline_version_can_be_pkg_added, true),            # 7
            # `guideline_version_has_osi_license` must be run after
            # `guideline_version_can_be_pkg_added` so that it can use the downloaded code!
            (guideline_version_has_osi_license, check_license),    # 8
            (guideline_version_can_be_imported, true),             # 9
        ]

    checked_guidelines = Guideline[]

    for (guideline, applicable) in guidelines
        applicable || continue
        if guideline == :update_status
            if !all(passed, checked_guidelines)
                update_status(data;
                              state = "failure",
                              context = "automerge/decision",
                              description = "New version. Failed.")
            end
        else
            check!(guideline, data)
            @info(guideline.info,
                  meets_this_guideline = passed(guideline),
                  message = message(guideline))
            push!(checked_guidelines, guideline)
        end
    end

    if all(passed, checked_guidelines) # success
        description = "New version. Approved. name=\"$(data.pkg)\". sha=\"$(data.current_pr_head_commit_sha)\""
        update_status(data;
                      state = "success",
                      context = "automerge/decision",
                      description = description)
        this_pr_comment_pass = comment_text_pass(data.registration_type,
                                                 data.suggest_onepointzero,
                                                 data.version)
        my_retry(() -> update_automerge_comment!(data, this_pr_comment_pass))
    else # failure
        update_status(data;
                      state = "failure",
                      context = "automerge/decision",
                      description = "New version. Failed.")
        failing_messages = message.(filter(!passed, checked_guidelines))
        this_pr_comment_fail = comment_text_fail(data.registration_type,
                                                 failing_messages,
                                                 data.suggest_onepointzero,
                                                 data.version)
        my_retry(() -> update_automerge_comment!(data, this_pr_comment_fail))
        throw(AutoMergeGuidelinesNotMet("The automerge guidelines were not met."))
    end
    return nothing
end


text:
using GNSSTools
using FFTW
using AllanDeviations
using Statistics
using LsqFit
using NNLS
using PyPlot
pygui(true)

function psd_at_f(h_parms, f, v_0)
    h₋₂, h₋₁, h₀, h₁, h₂ = h_parms
    return v_0^2 * (h₋₂*f^(-4) + h₋₁*f^(-3) + h₀*f^(-2) + h₁*f^(-1) + h₂*f^0)
end

function fit_h_parms_to_signal_PSD(psd_dBc, f_c=10e6,
                                   freqs=[1,10,100,1e3,1e4,1e5])
    N = length(freqs)
    @assert N == length(psd_dBc)
    b = zeros(N-1)
    A = zeros(N-1)
    for i in 1:N-1
        j = N - i
        A[i] = (psd_dBc[j] + psd_dBc[j+1])/2 + 10*log10(freqs[j+1] - freqs[j])
        if i == 1
            b[i] = 2*10^(A[i]/10) / (freqs[j+1] - freqs[j])  # b₀
        elseif i == 2
            b[i] = 2*10^(A[i]/10) / log(freqs[j+1] - freqs[j])  # b₋₁
        else
            b[i] = 2*10^(A[i]/10) / ((2-i)*(freqs[j+1]^(2-i) - freqs[j]^(2-i)))
        end
    end
    # order is h2, h1, h0, h-1, h-2
    h = b ./ f_c.^2
    # reverse order to h-2, h-1, h₀, h₁, h₂
    return reverse(h)
end

function calc_Sy(x, f_s)
    N = length(x)
    X = fft(x)[1:floor(Int, N/2)+1]
    X = (1/(f_s*N)) .* abs2.(X)
    X[2:end-1] .= 2 .* X[2:end-1]
    X_log10 = 10 .* log10.(X)
    frequencies = Array(range(0, f_s/2, length=length(X)))
    return (X, X_log10, frequencies)
end

t_length = 1.  # second
f_s = 5e6  # Hz
v_0 = L1_freq  # Hz
h_parms = h_parms_tcxo[4]
h₋₂, h₋₁, h₀, h₁, h₂ = h_parms

phase_noise = generate_phase_noise(t_length, f_s, v_0, h_parms)

truth = zeros(floor(Int, length(phase_noise)/2)+1)
freqs = Array(range(0, f_s/2, length=length(truth)))
for i in 2:length(truth)
    f = freqs[i]
    truth[i] = v_0^2 * (h₋₂*f^(-4) + h₋₁*f^(-3) + h₀*f^(-2) + h₁*f^(-1) + h₂*f^0)
end

# Take the difference
N = length(phase_noise)
X = fft(phase_noise)[1:floor(Int, N/2)+1]
X = (1/(f_s*N)) .* abs2.(X)
X[2:end-1] .= 2 .* X[2:end-1]

freqs2 = freqs[2:end]
difference = 10 .* log10.(abs.(truth .- X))  # dB
ax1 = fig.add_subplot(1, 1, 1)
ax1.plot(freqs, difference, "k")

fig = figure()
ax1 = fig.add_subplot(3, 1, 1)
ax1.plot(freqs2, 10 .* log10.(X[2:end]), "b", label="Generated Phase Noise PSD",
         linewidth=5)
ax1.plot(freqs2, 10 .* log10.(truth[2:end]), "r", label="Truth PSD",
         linestyle="--")
ax1.set_xscale("log")
xlabel("Frequency (Hz)")
ylabel("Power")
title("Generated and Truth PSDs")
legend()
ax2 = fig.add_subplot(3, 1, 2)
ax2.plot(freqs2, difference[2:end], "k")
ax2.set_xscale("log")
xlabel("Frequency (Hz)")
ylabel("Power")
title("PSD of the Difference Between the Generated and Truth")
ax3 = fig.add_subplot(3, 1, 3)
ax3.plot(Array(range(0, t_length, length=N)), phase_noise, "k")
xlabel("Time (s)")
ylabel("Amplitude (rad)")
title("Generated Phase Noise in Time Domain")
fig.subplots_adjust(hspace=0.9)

prn = 1
channel = "I"
signal_type = define_l1ca_code_type()
signal = definesignal(signal_type, f_s, t_length; f_d=800, CN0=38.5,
                      receiver_h_parms=h_parms, code_start_idx=100)
generatesignal!(signal)
acqresults, trackresults = process(signal, signal_type, prn, channel;
                                   show_plot=false)
# plotresults(trackresults)
phase_residuals = 2π.*signal.f_d.*trackresults.t .- trackresults.phi
fig = figure()
ax1 = fig.add_subplot(1, 1, 1)
# ax1.plot(trackresults.t, trackresults.phi, "k")
# ax1.plot(trackresults.t, 2π.*signal.f_d.*trackresults.t, "r")
ax1.plot(trackresults.t, phase_residuals, "k")
ax1.plot(signal.t, signal.phase_noise, "b")

################################################################################
f_carrier = 10e6
f_s = 25e6
t_length = 1.
t = Array(0:1/f_s:t_length-1/f_s)
carrier = cis.(2π.*f_carrier.*t)
N = length(t)
N_over_2 = floor(Int, length(carrier)/2)
frequencies = Array(range(0, f_s/2, length=N_over_2+1))
fig = figure()
ax1 = fig.add_subplot(1, 1, 1)
cm_red = get_cmap("Reds")
for i in 1:length(h_parms_tcxo)
    h_parms = h_parms_tcxo[i]
    ϕ_noise = real.(generate_phase_noise(t_length, f_s, f_carrier, h_parms))
    oscillator_noise = cis.(ϕ_noise)
    x = carrier .* oscillator_noise
    X = fft(x)[1:floor(Int, N/2)+1]
    X = (1/(f_s*N)) .* abs2.(X)
	X[2:end-1] .= 2 .* X[2:end-1]
    X_log10 = 10 .* log10.(X)
    max_val = maximum(X_log10)
    color = cm_red(1. * i / length(h_parms_tcxo))
    ax1.plot(frequencies .- f_carrier, X_log10 .- max_val, color=color)
end

cm_blue = get_cmap("Blues")
for i in 1:length(h_parms_ocxo)
    h_parms = h_parms_ocxo[i]
    ϕ_noise = real.(generate_phase_noise(t_length, f_s, f_carrier, h_parms))
    oscillator_noise = cis.(ϕ_noise)
    x = carrier .* oscillator_noise
    X = fft(x)[1:floor(Int, N/2)+1]
    X = (1/(f_s*N)) .* abs2.(X)
	X[2:end-1] .= 2 .* X[2:end-1]
    X_log10 = 10 .* log10.(X)
    max_val = maximum(X_log10)
    color = cm_blue(1. * i / length(h_parms_ocxo))
    ax1.plot(frequencies .- f_carrier, X_log10 .- max_val, color=color)
end

ϕ_noise = real.(generate_phase_noise(t_length, f_s, f_carrier, h_parms_cesium))
oscillator_noise = cis.(ϕ_noise)
x = carrier .* oscillator_noise
X = fft(x)[1:floor(Int, N/2)+1]
X = (1/(f_s*N)) .* abs2.(X)
X[2:end-1] .= 2 .* X[2:end-1]
X_log10 = 10 .* log10.(X)
max_val = maximum(X_log10)
ax1.plot(frequencies .- f_carrier, X_log10 .- max_val, "k--")

xscale("log")
xlim([1, 1e6])
ylim([-200, 0])
xlabel("Freguency (Hz)")
ylabel("dBc/Hz")

# Allan Deviation plot of first four TCXOs
fig = figure()
ax1 = fig.add_subplot(1, 1, 1)
for i in 1:4
    ϕ_noise = real.(generate_phase_noise(t_length, f_s, f_carrier, h_parms_tcxo[i]))
    oscillator_noise = cis.(ϕ_noise)
    x = carrier .* oscillator_noise
    result = allandev(oscillator_noise, f_s)
    ax1.plot(result.tau, result.deviation)
end
τ = Array(range(1e-7, 1, length=1000))
h₋₂, h₋₁, h₀, h₁, h₂ = h_parms_tcxo[4]
σₜ = zeros(length(τ))
for i in 1:length(τ)
    σₜ[i] = sqrt(h₋₂*τ[i]*(2π)^2/6 + h₋₁*2*log(2) + h₀/(2*τ[i]))
end
xscale("log")
yscale("log")

################################################################################
file = "/home/sjbilardi/Documents/SeNSeLab/GNSSTools.jl/test/data/known_oscillator_data/20210327_013800_000005_Aero-North-TCXO_1575.42_25.sc16"
data = loaddata(:sc16, file, 25e6, L1_freq, L1_freq, 1., skip_to=1)
replica = definesignal(define_l1ca_code_type(), 25e6, 1e-3)
for i in 1:32
    fd_est, n0_est, SNR_est = courseacquisition(data, replica, i)
    println("PRN $i\t$SNR_est")
end

# acqresults, trackresults = process(data, define_l1ca_code_type(), 1);

fig = figure()
for i in 1:32
    ax = fig.add_subplot(8, 4, i)
    acqresults, trackresults = process(data, define_l1ca_code_type(), i;
                                       show_plot=false, replica_t_length=1e-3,
                                       fine_acq_method=:carrier, state_num=3,
                                       dynamickf=false)
    ax.plot(trackresults.t, real.(trackresults.ZP))
    ax.plot(trackresults.t, imag.(trackresults.ZP))
    tick_params(left=false, bottom=false, labelleft=false, labelbottom=false)
    title("PRN $i")
end
subplots_adjust(hspace=0.4, wspace=0.4,  top=0.91, left=0.08, right=0.93,
                bottom=0.06)

# Track PRN 25
acqresults, trackresults = process(data, define_l1ca_code_type(), 25;
                                   show_plot=false, replica_t_length=1e-3,
                                   fine_acq_method=:carrier)

# Detrend phase and plot PSD
phi = trackresults.phi .- 2π.*median(trackresults.fds).*trackresults.t
# plot_spectrum(phi, 1/replica.t_length)

f_s = 1e3
N = length(phi)
X = fft(phi)[1:floor(Int, N/2)+1]
X = (1/(f_s*N)) .* abs2.(X)
X[2:end-1] .= 2 .* X[2:end-1]
X_phi = deepcopy(X)
frequencies = Array(range(0, f_s/2, length=floor(Int, N/2)+1))
frequencies_phi = deepcopy(frequencies)
v_0 = L1_freq
@. model(f, p) = v_0^2 * (p[1]*f^(-4) + p[2]*f^(-3) + p[3]*f^(-2) + p[4]*f^(-1) + p[5]*f^0)
@. model_truncated(f, p) = v_0^2 * (p[1]*f^(-4) + p[2]*f^(-3) + p[3]*f^(-2))
p₀ = h_parms_tcxo[4]
p₀2 = h_parms_tcxo[1]
# fit = curve_fit(model, frequencies[2:end], X[2:end], p₀)
fit = curve_fit(model_truncated, frequencies[2:end], X_phi[2:end], p₀2)

f_s = 5e6
t_length = 1.
# f_carrier = 10e6
# t_length = 4

# ϕ_noise = real.(generate_phase_noise(t_length, f_s, f_carrier, fit.param))
# t = Array(range(0, t_length, length=length(ϕ_noise)))
# carrier = cis.(2π.*f_carrier.*t)
# oscillator_noise = cis.(ϕ_noise)
# N = length(ϕ_noise)
# x = carrier .* oscillator_noise
# X = fft(x)[1:floor(Int, N/2)+1]
# X = (1/(f_s*N)) .* abs2.(X)
# X[2:end-1] .= 2 .* X[2:end-1]
# X_log10 = 10 .* log10.(X)
# max_val = maximum(X_log10)
frequencies = Array(range(0, f_s/2, length=floor(Int, f_s*t_length)))

truth = zeros(length(frequencies))
truth2 = zeros(length(frequencies))
for i in 2:length(frequencies)
    truth[i] = psd_at_f(fit.param, frequencies[i], L1_freq)
    truth2[i] = psd_at_f(fit.param, frequencies[i], 10e6)
end
truth = 10 .* log10.(truth[2:end])
truth2 = 10 .* log10.(truth2[2:end])

fig = figure()
ax1 = fig.add_subplot(1,1,1)
ax1.plot(frequencies_phi, 10 .* log10.(X_phi), label="Phase measurements from PLL")
ax1.plot(frequencies[2:end], truth, label="Three h-parameter fit")
ax1.plot(frequencies[2:end], truth2, label="Est. of phase noise at 10MHz")
xlim([1, 1e6])
xscale("log")
xlabel("Freguency (Hz)")
ylabel("dB/Hz")
legend()
# ax2 = fig.add_subplot(2,1,2)
# ax2.plot(frequencies .- f_carrier, X_log10 .- max_val, "k-")
# xscale("log")
# xlim([1, 1e6])
# ylim([-200, 0])
# xlabel("Freguency (Hz)")
# ylabel("dBc/Hz")

################################################################################\
# f_s = 36e6
f_carrier1 = 16.3676e6
f_carrier2 = 10e6
f_s = 1*f_carrier1
t_length = 1
N = floor(Int, t_length*f_s)
N_over_2 = floor(Int, N/2)
t = Array(range(0, t_length, length=N))
frequencies = Array(range(0, f_s/2, length=N_over_2+1))
# carrier1 = cis.(2π.*f_carrier1.*t)
# carrier2 = cis.(2π.*f_carrier2.*t)
phase_noise = zeros(Complex{Float64}, N)
# signal = zeros(Complex{Float64}, N)
dBc_freq = [1., 1e1, 1e2, 1e3, 1e4]
X = zeros(Complex{Float64}, N_over_2+1)
x = zeros(Complex{Float64}, N)
rakon_it5300B = Dict("h_parms"=>[1.9e-6, 2e-6, 9.7e-8, 9.1e-11, 9.9e-15],
                     "dBc"=>[-57, -88, -112, -130, -140],
                     "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
                     "name"=>"TCXO, Rakon IT5300B",
                     "f_0"=>f_carrier1)
agilent_e4424b = Dict("h_parms"=>[2.8e-9, 1.8e-11, 1.2e-10, 1e-11, 1.8e-15],
                      "dBc"=>[-88, -118, -128, -135, -147],
                      "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
                      "name"=>"OCXO, Agilent E4424B",
                      "f_0"=>10e6)
isotemp_91_1 = Dict("h_parms"=>[2.2e-9, 1.6e-10, 1.6e-10, 7.1e-13, 7.2e-16],
                    "dBc"=>[-89, -120, -140, -151, -154],
                    "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
                    "name"=>"OCXO, ISOTEMP 91-1",
                    "f_0"=>f_carrier2)
oscilloquartz_8607 = Dict("h_parms"=>[3.9e-14, 8.7e-13, 2.2e-14, 3.2e-13, 6.4e-15],
                          "dBc"=>[-122, -137, -143, -145, -145],
                          "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
                          "name"=>"OCXO, Oscilloquartz 8607",
                          "f_0"=>f_carrier2)
srs_prs10 = Dict("h_parms"=>[5.3e-11, 4.7e-11, 1.6e-15, 9.6e-14, 9.9e-16],
                 "dBc"=>[-103, -135, -150, -152, -153],
                 "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
                 "name"=>"Rubidum, SRS PRS10",
                 "f_0"=>f_carrier2)
cesium = Dict("h_parms"=>[1.5e-12, 6.2e-11, 7.6e-23, 1.4e-14, 6.2e-19],
              "dBc"=>[-105, -135, -160, -170, -190],
              "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
              "name"=>"Cesium, Typical Values",
              "f_0"=>f_carrier2)
mti_250l = Dict("h_parms"=>[1.7e-9, 1e-10, 1.6e-10, 3.5e-13, 1.8e-16],
              "dBc"=>[-90, -120, -140, -155, -160],
              "dBc_freq"=>[1., 1e1, 1e2, 1e3, 1e4],
              "name"=>"OCXO, MTI 250L",#, Low g-sensitivity",
              "f_0"=>f_carrier2)

oscillators = [rakon_it5300B, agilent_e4424b, isotemp_91_1, oscilloquartz_8607,
               srs_prs10, cesium, mti_250l]
# figsize = (6.5, 4)
colorz = ["k", "b", "g", "r", "y", "m", "c", "gray"]
fig = figure()
ax1 = fig.add_subplot(1,1,1)
for i in 1:length(oscillators)
    color = colorz[i]
    oscillator = oscillators[i]
    if oscillator["name"] == "TCXO, Rakon IT5300B"
        # carrier = carrier1
        f_carrier = f_carrier1
    else
        # carrier = carrier2
        f_carrier = f_carrier2
    end
    h_parms = oscillator["h_parms"] ./ oscillator["f_0"]^2
    generate_phase_noise!(phase_noise, t_length, f_carrier, h_parms)
    ϕ = real.(phase_noise)
    # signal = carrier .* cis.(ϕ)
    # signal .= carrier .* cis.(real.(phase_noise))
    # x = signal
    x .= real.(phase_noise) .+ 0im
    # X, X_log, freqs = calc_Sy(x, f_s)
    fft!(x);
    X .= abs2.(x[1:N_over_2+1]) / (f_s*N);
    X[2:end-1] .= 2 .* X[2:end-1];
    X .= 10 .* log10.(X);
    # X_max = maximum(real.(X));
    # X = X .- X_max;
    # ax1.plot(frequencies[2:end] .- f_carrier, X[2:end], color,
    #          label=string(oscillator["name"], " from h Parms"))
    ax1.plot(frequencies[2:end], real.(X[2:end]), color,
             label=string(oscillator["name"], " from h Parms"))
    ax1.scatter(dBc_freq, oscillator["dBc"].+(10*log10(2)),
                c=color, label=string(oscillator["name"], " Truth"))
end
xscale("log")
xlim([1, 1e6])
xlabel("Frequency (Hz)")
ylabel("dBc/Hz")
legend()

@. function model(x, p)
    f = x[1]
    y = x[2]
    return  y - p[1]*f^(-4) + p[2]*f^(-3) + p[3]*f^(-2) + p[4]*f^(-1) + p[5]*f^0
end

@. model(x, p) = v_0^2 * (p[1]*f^(-4) + p[2]*f^(-3) + p[3]*f^(-2) + p[4]*f^(-1) + p[5]*f^0)