text:
# MIT license
# Copyright (c) Microsoft Corporation. All rights reserved.
# See LICENSE in the project root for full license information.

using Test
using OpticSim
using OpticSim.Emitters
using OpticSim.Geometry

using Random
using StaticArrays

@testset "Emitters" begin
    @testset "AngularPower" begin
        @testset "Lambertian" begin
            @test typeof(AngularPower.Lambertian()) === AngularPower.Lambertian{Float64}
            @test_throws MethodError AngularPower.Lambertian(String)

            @test apply(AngularPower.Lambertian(), Transform(), 1, Ray(zeros(3), ones(3))) === 1
        end

        @testset "Cosine" begin
            @test AngularPower.Cosine().cosine_exp === one(Float64)
            @test_throws MethodError AngularPower.Cosine(String)

            @test apply(AngularPower.Cosine(), Transform(), 1, Ray(zeros(3), ones(3))) === 0.5773502691896258
        end

        @testset "Gaussian" begin
            @test AngularPower.Gaussian(0, 1).gaussianu === 0
            @test AngularPower.Gaussian(0, 1).gaussianv === 1

            @test apply(AngularPower.Gaussian(0, 1), Transform(), 1, Ray(zeros(3), ones(3))) === 0.7165313105737892
        end
    end

    @testset "Directions" begin
        @testset "Constant" begin
            @test Directions.Constant(0, 0, 0).direction === Vec3(Int)
            @test Directions.Constant(Vec3()).direction === Vec3()
            @test Directions.Constant().direction === unitZ3()

            @test Base.length(Directions.Constant()) === 1
            @test Emitters.generate(Directions.Constant(), 0) === unitZ3()

            @test Base.iterate(Directions.Constant()) === (unitZ3(), 2)
            @test Base.iterate(Directions.Constant(), 2) === nothing
            @test Base.getindex(Directions.Constant(), 1) === unitZ3()
            @test Base.getindex(Directions.Constant(), 2) === unitZ3()
            @test Base.firstindex(Directions.Constant()) === 0
            @test Base.lastindex(Directions.Constant()) === 0
            @test Base.copy(Directions.Constant()) === Directions.Constant()
        end

        @testset "RectGrid" begin
            @test Directions.RectGrid(unitX3(), 0., 0., 0, 0).uvec === unitY3()
            @test Directions.RectGrid(unitX3(), 0., 0., 0, 0).vvec === unitZ3()
            @test Directions.RectGrid(0., 0., 0, 0).uvec === unitX3()
            @test Directions.RectGrid(0., 0., 0, 0).vvec === unitY3()

            @test Base.length(Directions.RectGrid(0., 0., 2, 3)) === 6
            @test collect(Directions.RectGrid(ones(Vec3), π/4, π/4, 2, 3)) == [
                [0.6014252485829169, 0.7571812496798511, 1.2608580392339483],
                [1.1448844430815608, 0.4854516524305293, 0.9891284419846265],
                [0.643629619770125, 1.0798265148205617, 1.0798265148205617],
                [1.225225479837374, 0.7890285847869373, 0.7890285847869373],
                [0.6014252485829169, 1.2608580392339483, 0.7571812496798511],
                [1.1448844430815608, 0.9891284419846265, 0.4854516524305293],
            ]
        end

        @testset "UniformCone" begin
            @test Directions.UniformCone(unitX3(), 0., 0).uvec === unitY3()
            @test Directions.UniformCone(unitX3(), 0., 0).vvec === unitZ3()
            @test Directions.UniformCone(0., 0).uvec === unitX3()
            @test Directions.UniformCone(0., 0).vvec === unitY3()

            @test Base.length(Directions.UniformCone(0., 1)) === 1

            @test collect(Directions.UniformCone(π/4, 2, rng=Random.MersenneTwister(0))) == [
                [0.30348115383395624, -0.6083405920618145, 0.7333627433388552],
                [0.16266571675478964, 0.2733479444418462, 0.9480615833700192],
            ]
        end

        @testset "HexapolarCone" begin
            @test Directions.HexapolarCone(unitX3(), 0., 0).uvec === unitY3()
            @test Directions.HexapolarCone(unitX3(), 0., 0).vvec === unitZ3()
            @test Directions.HexapolarCone(0., 0).uvec === unitX3()
            @test Directions.HexapolarCone(0., 0).vvec === unitY3()

            @test Base.length(Directions.HexapolarCone(0., 1)) === 7
            @test Base.length(Directions.HexapolarCone(0., 2)) === 19
            @test collect(Directions.HexapolarCone(π/4, 1)) == [
                [0.0, 0.0, 1.0],
                [0.7071067811865475, 0.0, 0.7071067811865476],
                [0.3535533905932738, 0.6123724356957945, 0.7071067811865476],
                [-0.3535533905932736, 0.6123724356957946, 0.7071067811865477],
                [-0.7071067811865475, 8.659560562354932e-17, 0.7071067811865476],
                [-0.35355339059327406, -0.6123724356957944, 0.7071067811865476],
                [0.3535533905932738, -0.6123724356957945, 0.7071067811865476],
            ]
        end
    end

    @testset "Origins" begin
        @testset "Point" begin
            @test Origins.Point(Vec3()).origin === Vec3()
            @test Origins.Point(0, 0, 0).origin === Vec3(Int)
            @test Origins.Point().origin === Vec3()

            @test Base.length(Origins.Point()) === 1
            @test Emitters.visual_size(Origins.Point()) === 1
            @test Emitters.visual_size(Origins.Point()) === 1
            @test Emitters.generate(Origins.Point(), 1) === Vec3()

            @test Base.iterate(Origins.Point(), 1) === (Vec3(), 2)
            @test Base.iterate(Origins.Point(), 2) === nothing
            @test Base.getindex(Origins.Point(), 1) === Vec3()
            @test Base.getindex(Origins.Point(), 2) === Vec3()
            @test Base.firstindex(Origins.Point()) === 0
            @test Base.lastindex(Origins.Point()) === 0
            @test Base.copy(Origins.Point()) === Origins.Point()
        end

        @testset "RectUniform" begin
            @test Origins.RectUniform(1, 2, 3).width === 1
            @test Origins.RectUniform(1, 2, 3).height === 2
            @test Origins.RectUniform(1, 2, 3).samples_count === 3

            @test Base.length(Origins.RectUniform(1, 2, 3)) === 3
            @test Emitters.visual_size(Origins.RectUniform(1, 2, 3)) === 2

            @test collect(Origins.RectUniform(1, 2, 3, rng=Random.MersenneTwister(0))) == [
                [0.3236475079774124, 0.8207130758528729, 0.0],
                [-0.3354342018663148, -0.6453423070674709, 0.0],
                [-0.221119890668799, -0.5930468839161547, 0.0],
            ]
        end

        @testset "RectGrid" begin
            @test Origins.RectGrid(1., 2., 3, 4).ustep === .5
            @test Origins.RectGrid(1., 2., 3, 4).vstep === 2/3

            @test Base.length(Origins.RectGrid(1., 2., 3, 4)) === 12
            @test Emitters.visual_size(Origins.RectGrid(1., 2., 3, 4)) === 2.

            @test collect(Origins.RectGrid(1., 2., 3, 4)) == [
                [-0.5, -1.0, 0.0],
                [0.0, -1.0, 0.0],
                [0.5, -1.0, 0.0],
                [-0.5, -0.33333333333333337, 0.0],
                [0.0, -0.33333333333333337, 0.0],
                [0.5, -0.33333333333333337, 0.0],
                [-0.5, 0.33333333333333326, 0.0],
                [0.0, 0.33333333333333326, 0.0],
                [0.5, 0.33333333333333326, 0.0],
                [-0.5, 1.0, 0.0],
                [0.0, 1.0, 0.0],
                [0.5, 1.0, 0.0],
            ]
        end

        @testset "RectJitterGrid" begin
            @test Origins.RectJitterGrid(1., 2., 3, 4, 1).width === 1.0
            @test Origins.RectJitterGrid(1., 2., 3, 4, 1).height === 2.0
            @test Origins.RectJitterGrid(1., 2., 2, 2, 3).ustep === 0.5
            @test Origins.RectJitterGrid(1., 2., 2, 2, 3).vstep === 1.0

            @test Base.length(Origins.RectJitterGrid(1., 2., 5, 6, 7)) === 210
            @test Emitters.visual_size(Origins.RectJitterGrid(1., 2., 5, 6, 7)) === 2.0

            @test collect(Origins.RectJitterGrid(1., 2., 2, 2, 2, rng=Random.MersenneTwister(0))) == [
                [-0.08817624601129381, -0.08964346207356355, 0.0],
                [-0.4177171009331574, -0.8226711535337354, 0.0],
                [0.1394400546656005, -0.7965234419580773, 0.0],
                [0.021150832966014832, -0.9317307444943552, 0.0],
                [-0.3190858046118913, 0.9732164043865108, 0.0],
                [-0.2070942241283379, 0.5392892841426182, 0.0],
                [0.13001792513452393, 0.910046541351011, 0.0],
                [0.08351809722107484, 0.6554484126999125, 0.0],
            ]
        end

        @testset "Hexapolar" begin
            @test Origins.Hexapolar(1, 0, 0).nrings === 1
            @test_throws MethodError Origins.Hexapolar(1., 0, 0)

            @test Base.length(Origins.Hexapolar(1, 0, 0)) === 7
            @test Base.length(Origins.Hexapolar(2, 0, 0)) === 19
            @test Emitters.visual_size(Origins.Hexapolar(0, 1, 2)) === 4

            @test collect(Origins.Hexapolar(1, π/4, π/4)) == [
                [0.0, 0.0, 0.0],
                [0.7853981633974483, 0.0, 0.0],
                [0.39269908169872425, 0.6801747615878316, 0.0],
                [-0.392699081698724, 0.6801747615878317, 0.0],
                [-0.7853981633974483, 9.618353468608949e-17, 0.0],
                [-0.3926990816987245, -0.6801747615878315, 0.0],
                [0.39269908169872425, -0.6801747615878316, 0.0],
            ]
        end
    end

    @testset "Spectrum" begin
        @testset "Uniform" begin
            @test Spectrum.UNIFORMSHORT === 0.450
            @test Spectrum.UNIFORMLONG === 0.680
            @test Spectrum.Uniform(0, 1).low_end === 0
            @test Spectrum.Uniform(0, 1).high_end === 1
            @test Spectrum.Uniform().low_end === 0.450
            @test Spectrum.Uniform().high_end === 0.680

            @test Emitters.generate(Spectrum.Uniform(rng=Random.MersenneTwister(0))) === (1.0, 0.6394389268348049)
        end

        @testset "DeltaFunction" begin
            @test Emitters.generate(Spectrum.DeltaFunction(2)) === (1, 2)
            @test Emitters.generate(Spectrum.DeltaFunction(2.)) === (1., 2.)
            @test Emitters.generate(Spectrum.DeltaFunction(π)) === (true, π) # hopefully this is ok!
        end

        @testset "Measured" begin
            # TODO
        end
    end

    @testset "Sources" begin
        expected_rays = [
            OpticalRay([0., 0., 0.], [0., 0., 1.], 1., 0.6394389268348049),
            OpticalRay([0., 0., 0.], [0., 0., 1.], 1., 0.6593820037230804),
            OpticalRay([0., 0., 0.], [0., 0., 1.], 1., 0.48785013357074763),
        ]

        @testset "Source" begin
            @test Sources.Source().transform === Transform()
            @test Sources.Source().spectrum === Spectrum.Uniform()
            @test Sources.Source().origins === Origins.Point()
            @test Sources.Source().directions === Directions.Constant()
            @test Sources.Source().power_distribution === AngularPower.Lambertian()
            @test Sources.Source().sourcenum === 0

            @test Sources.Source() === Sources.Source(
                Transform(),
                Spectrum.Uniform(),
                Origins.Point(),
                Directions.Constant(),
                AngularPower.Lambertian())

            @test Base.length(Sources.Source()) === 1
            @test Base.length(Sources.Source(; origins=Origins.Hexapolar(1, 0., 0.))) === 7
            @test Base.length(Sources.Source(; directions=Directions.HexapolarCone(0., 1))) === 7
            @test Base.length(Sources.Source(;
                origins=Origins.Hexapolar(1, 0., 0.),
                directions=Directions.HexapolarCone(0., 1)
            )) === 49

            @test Base.iterate(Sources.Source(spectrum=Spectrum.Uniform(rng=Random.MersenneTwister(0)))) === (expected_rays[1], Sources.SourceGenerationState(2, 0, Vec3()))

            @test Base.getindex(Sources.Source(spectrum=Spectrum.Uniform(rng=Random.MersenneTwister(0))), 0) === expected_rays[1]
            @test Emitters.generate(Sources.Source(spectrum=Spectrum.Uniform(rng=Random.MersenneTwister(0))), 0) === expected_rays[1]

            @test Emitters.generate(Sources.Source(spectrum=Spectrum.Uniform(rng=Random.MersenneTwister(0)))) === (expected_rays[1], Sources.SourceGenerationState(0, -2, Vec3()))

            @test Base.firstindex(Sources.Source()) === 0
            @test Base.lastindex(Sources.Source()) === 0
            @test Base.copy(Sources.Source()) === Sources.Source()
        end

        @testset "CompositeSource" begin
            s() = Sources.Source(spectrum=Spectrum.Uniform(rng=Random.MersenneTwister(0)))
            tr = Transform()
            cs1 = Sources.CompositeSource(tr, [s()])
            cs2 = Sources.CompositeSource(tr, [s(), s()])
            cs3 = Sources.CompositeSource(tr, [s(), cs2])

            @test cs1.transform === tr

            @test cs1.uniform_length === 1
            @test cs2.uniform_length === 1
            @test cs3.uniform_length === -1

            @test cs1.total_length === 1
            @test cs2.total_length === 2
            @test cs3.total_length === 3

            @test cs1.start_indexes == []
            @test cs2.start_indexes == []
            @test cs3.start_indexes == [0, 1, 3]

            @test Base.length(cs1) === 1
            @test Base.length(cs2) === 2
            @test Base.length(cs3) === 3

            @test Base.iterate(cs1) === (expected_rays[1], Sources.SourceGenerationState(2, 0, Vec3()))

            @test collect(cs2) == vcat(expected_rays[1:1], expected_rays[1:1])

            @test collect(cs3) == vcat(expected_rays[1:1], expected_rays[2:2], expected_rays[2:2])
        end
    end
end # testset Emitters


text:
include("obj_io.jl")
include("surface_mesh.jl")

ObjIO.tests()
SurfaceMesh.tests()

# some checker bump-map. It can be any noise function for instance
bump_map_w = 64
bump_map_h = 64
bump_map = zeros(Float64, bump_map_w, bump_map_h)
for i in 1:bump_map_h
  for j in 1:bump_map_w
    bump_map[i, j] = xor(((i - 1) % 16 < 8), ((j - 1) % 16 < 8)) ? 1.0 : 0.0
  end
end

bump_from_uv = function(uv)
  # bump-map coordinates
  j = 1 + uv[1] * bump_map_w
  i = 1 + uv[2] * bump_map_h
  j0 = floor(Int, j)
  j1 = min(floor(Int, j) + 1, bump_map_w)
  i0 = floor(Int, i)
  i1 = min(floor(Int, i) + 1, bump_map_h)
#  bump = bump_map[i0, j0]
  bump = bump_map[i1, j1] * (j - j0) * (i - i0)   # syntax thing, you can't just split a string into several
  bump += bump_map[i1, j0] * (j1 - j) * (i - i0)
  bump += bump_map[i0, j0] * (j1 - j) * (i1 - i)
  bump += bump_map[i0, j1] * (j - j0) * (i1 - i)
  return bump
end

if max([bump_from_uv((0.1, v)) for v in 0.1:0.001:0.15]...) != 1.
  println("bilinear interpolation is broken")
end

if min([bump_from_uv((u, 0.1)) for u in 0.1:0.001:0.15]...) != 0.
  println("bilinear interpolation is broken")
end

distance_function = function(x)
  d = sqrt(x[1]^2 + x[2]^2 + x[3]^2) - 10
  uv = (x[1] / 10., x[2] / 10.)
#  if uv[1] > 0. && uv[1] < 1. && uv[2] > 0. && uv[2] < 1.
#    d -= bump_from_uv(uv)
#  end
  return d
end
gradient_function = SurfaceMesh.distance_to_gradient_operator(distance_function)
(vertexes, faces) = SurfaceMesh.build_3d_mesh((-10, -10, -10), (10, 10, 10), 1., distance_function, gradient_function)

open("test.obj", "w") do file
  write(file, ObjIO.str_from_vertexes(vertexes))
  write(file, ObjIO.str_from_faces(faces))
end

ext: jl


text:
module TestMmap
using Test
using PhaseSpaceIO
using PhaseSpaceIO.Testing
using PhaseSpaceIO: Format, FormatEGS, FormatIAEA

@testset "mmap $P" for P in [EGSParticle{Nothing},
                          EGSParticle{Float32},
                          IAEAParticle{0,0},
                          IAEAParticle{1,0},
                          IAEAParticle{0,1},
                          IAEAParticle{rand(0:3), rand(0:3)},
                         ]

    mktempdir() do dir
        for len in [1,2,3,4,10,100,rand(1000:2000)]
            ps = [arbitrary(P) for _ in 1:len]
            if Format(first(ps)) isa FormatEGS
                path = joinpath(dir, "phsp.egsphsp")
            else
                path = joinpath(dir, "phsp.IAEAheader")
            end
            phsp_write(path, ps)

            ps1 = phsp_iterator(collect, path)
            ps2 = PhspVector(path)
            @test all(ps1 .≈ ps)
            @test all(ps2 .≈ ps)

            p1 = first(ps)
            p2 = first(ps2)
            @test p1.w == p2.w
            @test length(ps2) == length(ps)
            @test first(ps) == first(ps2)
            @test last(ps) == last(ps2)
            @test ps1 == ps2
            # for (p1, p2) in zip(ps1, ps2)
            #     @test p1.u == p2.u
            #     @test p1.v == p2.v
            #     @test p1.w == p2.w
            #     @test p1 == p2
            # end
        end
    end
end

end#module


text:
include("actorinterfacestests.jl")
include("stack.jl")
include("basics.jl")


text:
function resample(n_samples, lbs, ubs, prev_X, x_sol; mean_shift=0.75, std_shrink = 0.3)
    """ Gaussian resamples free variables depending on previous solution.
    Arguments:
    - n_samples: number of samples
    - lbs: lower bounds
    - bs: upper bounds
    - prev_X: previous data
    - x_sol: solution array with same # of columns as lbs/ubs
    - sigma: by how much to move the mean to new optimum, and to shrink the std.
    Returns:
    - new_samples
    """
    X = deepcopy(prev_X);
    prev_samples, n_features = size(prev_X);
    prev_mean = [sum(X[:,i])/prev_samples for i=1:n_features];
    prev_std = [sum((X[:,i]-ones(prev_samples)* prev_mean[i]).^2)/prev_samples for i=1:n_features];
    new_std = prev_std*(1-std_shrink);
    new_mean = mean_shift*x_sol + (1-mean_shift)*prev_mean;
    dist = MvNormal(new_mean, new_std);
    while size(X,1) < prev_samples + n_samples
        new_row = rand(dist)';
        if all(lbs .<= new_row) && all(ubs .>= new_row);
            X = vcat(X, Matrix(new_row));
        end
    end
    return X
end

# Sampling based on GaussianProcesses
#         plan = randomLHC(Int(round(get_param(bbf, :n_samples)*ratio)), n_dims);
#         random_samples = scaleLHC(plan,[(bbf.lbs[i], bbf.ubs[i]) for i in vks]);
#         μ, σ = predict_f(bbf.gp, random_samples');
#         cdf_0 = [Distributions.cdf(Distributions.Normal(μ[i], σ[i]),0) for i=1:size(random_samples,1)];
#          #TODO: add criterion for information as well (something like sortperm(σ))
#         # Sample places with high probability of being near boundary (0),
#         # but also balance feasibility ratio.
#         p = bbf.feas_ratio
# #         balance_fn = x -> -1*tan(2*atan(-0.5)*(x-0.5)) + 0.5
#         balance_fn = x -> -1/0.5^2*(x-0.5)^3 + 0.5;
#         indices = sortperm(abs.(cdf_0 .- balance_fn(p)));
#         samples = DataFrame(random_samples[indices[1:get_param(bbf, :n_samples)],:], vks);
#         eval!(bbf, samples);

# function predict(bbf::BlackBoxFunction, X::AbstractArray)
#     μ, σ = predict_f(bbf.gp, transpose(X))
#     return μ, σ
# end
#
# function optimize_gp!(bbf::BlackBoxFunction)
#     """ Optimizes a GaussianProcess over a BlackBoxFunction,
#     with adaptively changing kernel. """
# #         bbf.gp = ElasticGPE(length(bbf.idxs), # data
# #         mean = MeanConst(sum(bbf.Y)/length(bbf.Y)), logNoise = -10)
#     lbs = [bbf.lbs[key] for key in bbf.vks];
#     ubs = [bbf.ubs[key] for key in bbf.vks];
#     bbf.gp = GPE(transpose(Array(bbf.X)), bbf.Y, # data
#     MeanConst(sum(bbf.Y)/length(bbf.Y)),
#     SEArd(log.((ubs-lbs)./(2*sqrt(length(bbf.Y)))), -5.))
#     optimize!(bbf.gp); #TODO: optimize GP
#                        # Instead of regenerating at every run, figure out
#                        # how to update.
# end


text:
using OrdinaryDiffEq, Test, DiffEqOperators, Random, LinearAlgebra, SparseArrays
let N = 20
Random.seed!(0); u0 = normalize(randn(N))
dd = -2 * ones(N); du = ones(N-1)
A = diagm(-1 => du, 0 => dd, 1 => du)
_f = (u,p,t) -> A*u - u.^3
_f_ip = (du,u,p,t) -> (mul!(du, A, u); du .-= u.^3)
_jac = (u,p,t) -> A - 3 * diagm(0 => u.^2)
_jac_ip = (J,u,p,t) -> begin
  copyto!(J, A)
  @inbounds for i = 1:N
    J[i, i] -= 3 * u[i]^2
  end
end
# f = ODEFunction(_f; jac=_jac)
# f_ip = ODEFunction(_f_ip; jac=_jac_ip, jac_prototype=zeros(N,N))
jac_prototype = DiffEqArrayOperator(zeros(N,N); update_func=_jac_ip)
f = ODEFunction(_f; jac_prototype=jac_prototype)
f_ip = ODEFunction(_f_ip; jac_prototype=jac_prototype)
prob = ODEProblem(f, u0, (0.0, 1.0))
prob_ip = ODEProblem(f_ip, u0, (0.0, 1.0))

@testset "Classical ExpRK - Low Order" begin
  dt = 0.01; tol=1e-3
  Algs = [LawsonEuler,NorsettEuler,ETDRK2]
  for Alg in Algs
    sol = solve(prob, Alg(krylov=true, m=20); dt=dt, reltol=tol)
    sol_ref = solve(prob, Tsit5(); reltol=tol)
    @test isapprox(sol(1.0), sol_ref(1.0); rtol=tol)

    sol = solve(prob_ip, Alg(krylov=true, m=20); dt=dt, reltol=tol)
    sol_ref = solve(prob_ip, Tsit5(); reltol=tol)
    @test isapprox(sol(1.0), sol_ref(1.0); rtol=tol)

    println(Alg) # prevent Travis hanging
  end
end

@testset "Classical ExpRK - High Order" begin
  dt = 0.05; tol=1e-5
  Algs = [ETDRK3,ETDRK4,HochOst4]
  for Alg in Algs
    sol = solve(prob, Alg(krylov=true, m=20); dt=dt, reltol=tol)
    sol_ref = solve(prob, Tsit5(); reltol=tol)
    @test isapprox(sol(1.0), sol_ref(1.0); rtol=tol)

    sol = solve(prob_ip, Alg(krylov=true, m=20); dt=dt, reltol=tol)
    sol_ref = solve(prob_ip, Tsit5(); reltol=tol)
    @test isapprox(sol(1.0), sol_ref(1.0); rtol=tol)

    println(Alg) # prevent Travis hanging
  end
end

@testset "EPIRK" begin
  dt = 0.05; tol=1e-5
  Algs = [Exp4, EPIRK4s3A, EPIRK4s3B, EXPRB53s3, EPIRK5P1, EPIRK5P2]
  for Alg in Algs
    sol = solve(prob, Alg(); dt=dt, reltol=tol)
    sol_ref = solve(prob, Tsit5(); reltol=tol)
    @test isapprox(sol(1.0), sol_ref(1.0); rtol=tol)

    sol = solve(prob_ip, Alg(); dt=dt, reltol=tol)
    sol_ref = solve(prob_ip, Tsit5(); reltol=tol)
    @test isapprox(sol(1.0), sol_ref(1.0); rtol=tol)
    println(Alg) # prevent Travis hanging
  end

  sol = solve(prob, EPIRK5s3(); dt=dt, reltol=tol)
  sol_ref = solve(prob, Tsit5(); reltol=tol)
  @test_broken isapprox(sol(1.0), sol_ref(1.0); rtol=tol)

  sol = solve(prob_ip, EPIRK5s3(); dt=dt, reltol=tol)
  sol_ref = solve(prob_ip, Tsit5(); reltol=tol)
  @test_broken isapprox(sol(1.0), sol_ref(1.0); rtol=tol)
  println(EPIRK5s3) # prevent Travis hanging
end

@testset "Adaptive exponential Rosenbrock" begin
  # Regression tests adapted from ode_dense_tests.jl
  interp_points = 0.0:1/16:1.0
  interp_results = [zeros(N) for _ = 1:length(interp_points)]
  function regression_test(prob, alg, tol)
    sol1 = solve(prob, alg, dt=1/4, dense=true, adaptive=true)
    sol1(interp_results, interp_points)
    sol2 = solve(prob, alg, dt=1/16, dense=true, adaptive=false)
    for i in eachindex(sol2)
      err = maximum(abs.(sol2[i] - interp_results[i]))
      @test err < tol
    end
  end

  println("Exprb32, out-of-place")
  regression_test(prob, Exprb32(m=N), 3e-4)
  println("Exprb32, inplace")
  regression_test(prob_ip, Exprb32(m=N), 3e-4)
  println("Exprb43, out-of-place")
  regression_test(prob, Exprb43(m=N), 3e-4)
  println("Exprb43, inplace")
  regression_test(prob_ip, Exprb43(m=N), 3e-4)
end
end

@testset "ExpRK with custom jacobian" begin
  N = 10
  # Sparse Jacobian
  Random.seed!(0); u0 = normalize(randn(N))
  dd = -2 * ones(N); du = ones(N-1)
  A = spdiagm(-1 => du, 0 => dd, 1 => du)
  f = (u,p,t) -> A*u
  exp_fun = ODEFunction(f;
                    jac=(u,p,t) -> A,
                    analytic=(u,p,t) -> exp(t*Matrix(A)) * u)
  prob = ODEProblem(exp_fun, u0, (0.0,1.0))
  sol = solve(prob, LawsonEuler(krylov=true, m=N); dt=0.1)
  @test sol(1.0) ≈ exp_fun.analytic(u0,nothing,1.0)
  # Matrix-free Jacobian
  # Need to implement the missing interface for DerivativeOperator first
  @test_broken begin
    L = DerivativeOperator{Float64}(2,2,1.0,N,:Dirichlet0,:Dirichlet0)
    exp_fun2 = ODEFunction(L;
                      jac_prototype=L,
                      analytic=(u,p,t) -> exp(t*full(L)) * u)
    prob = ODEProblem(exp_fun2, u0, (0.0,1.0))
    sol = solve(prob, LawsonEuler(krylov=true, m=N); dt=0.1)
    @test sol(1.0) ≈ exp_fun2.analytic(u0,nothing,1.0)
  end
end

@testset "ExpRK with default jacobian" begin
  N = 10
  Random.seed!(0); u0 = normalize(randn(N))
  dd = -2 * ones(N); du = ones(N-1)
  A = diagm(-1 => du, 0 => dd, 1 => du)
  f = (du,u,p,t) -> mul!(du,A,u)
  jac = (J,u,p,t) -> (J .= A; nothing)
  exp_fun = ODEFunction(f; jac=jac, analytic=(u,p,t) -> exp(t*A)*u)
  prob = ODEProblem(exp_fun, u0, (0.0,1.0))
  sol = solve(prob, LawsonEuler(krylov=true, m=N); dt=0.1)
  @test sol(1.0) ≈ exp_fun.analytic(u0,nothing,1.0)
end


text:
"""
  RestInterfaces

Library for building REST APIs.
"""
module RestInterfaces

include("./HttpErrors.jl")
include("./HttpMethods.jl")
include("./Resources.jl")
include("./Middleware.jl")
include("./Utils.jl")

end # module Rest


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

using Documenter, Korg

makedocs(sitename="Korg",
         modules=[Korg],
         pages=[
                "Quickstart" => "index.md"
                "Function Reference" => "API.md"
                "Developer Documentation" => "devdocs.md"
                "Install" => "install.md"
               ],
        authors="Adam Wheeler and Matthew Abruzzo",
        format=Documenter.HTML(assets=["assets/favicon.ico"])
       )

deploydocs(repo = "github.com/ajwheeler/Korg.jl.git",
           devbranch="main")


text:
using Mongo
using PropertyGraph
using PropertyGraphMongo

using LibBSON
using Base.Test
using UUID
using Compat

mongoDBDir = "/tmp/Mongo.jl-test.db"
mkpath(mongoDBDir)
mongod = spawn(`mongod --dbpath $mongoDBDir`)
sleep(1) # wait for listen on port

try
	testfiles = [
				"testsavechanges",
				"testquery",
				]

	println("Running tests:")

	for t in testfiles
		testfilepath = joinpath(Pkg.dir("PropertyGraphMongo"),"test","$(t).jl")
		println("running $(testfilepath) ...")
		include(testfilepath)
	end

finally
	kill(mongod)
	rm(mongoDBDir, recursive=true)
end


text:
export DMD_PG_Options, solveDMD_withPG

# This file is available under the terms of the MIT License

mutable struct DMD_PG_Options{T <: AbstractFloat}
    itm::Integer
    tol::T
    ptf::Integer

    prox::Function
end

function DMD_PG_Options(;itm=100,tol=1e-5,ptf=10,
                        prox=prox_null::Function)

    return DMD_PG_Options(itm, tol, ptf, prox)
end

function solveDMD_withPG(params, opts)
    # iteration parameters
    itm = opts.itm;
    tol = opts.tol;
    ptf = opts.ptf;
    # prox function
    prox = opts.prox;
    # problem dimensions
    m = params.m;
    n = params.n;
    k = params.k;

    T = typeof(real(params.X[1]))

    # pre-allocate the variables
    gar  = zeros(T, k << 1);
    garp1 = zeros(T, k << 1);
    ar  = copy(params.ar);
    prox(ar)
    arm1 = copy(ar);

    function a_fg!(gar, arnew)
        copyto!(params.ar, arnew);
        return aBFuncGrad(gar, params)
    end

    function a_f(arnew)
        copyto!(params.ar, arnew);
        return aBFunc(params);
    end

    # initialize variables
    obj = a_fg!(gar, ar);
    err = T(1.0);
    noi = 0;

    if ( isnan(obj) || isinf(obj) )
        println("BAD INITIAL GUESS ", ar)
        println("ABORT")
        return
    end

    obj_his = zeros(itm)
    err_his = zeros(itm)
    lns_his = zeros(itm)

    alphastart = 10.0

    while err >= tol
        eta, lns_noi = descent_PG!(arm1, ar, gar, obj, a_f, prox,
                                   alphastart=alphastart);
        prox(ar);

        obj = a_fg!(gar, ar);
        err = sqrt(sum(abs2, ar - arm1)) / eta;
        # err = sqrt(sum(abs2,ar - arm1));
        # err = vecnorm(ar - arm1)
        copyto!(arm1, ar);
        noi += 1;
        # @show gar;
        # @show ar;
        (noi % ptf == 0) && @printf("PG: iter %4d, obj %1.6e, err %1.6e\n",
                                                    noi, obj, err);

        obj_his[noi] = obj
        err_his[noi] = err
        lns_his[noi] = lns_noi
        noi >= itm && break;
        alphastart = min(2 * eta, 10.0)
    end

    return obj_his[1:noi], err_his[1:noi], lns_his[1:noi]

end

# --------------------------
# Descent Line Search
# --------------------------
function descent_PG!(xm1, x, g, obj, f, prox;
                     alphamin=1e1 * eps(real(eltype((x)))),
                     tol=1e2 * eps(real(eltype(x))),
                     alphastart=10.0)
    T = real(eltype(x))
    n = length(x);
    alpha = T(alphastart);
    noi = 0;
    for i = 1:n
        x[i] = xm1[i] - alpha * g[i];
    end
    prox(x);
    objp1 = f(x);
    while objp1 > obj || isnan(objp1)
        alpha *= T(0.5);
        for i = 1:n
            x[i] = xm1[i] - alpha * g[i];
        end
        prox(x);
        objp1 = f(x);
        noi += 1;
        # noi ≥ 5 && break;
        alpha < alphamin && break;
    end
    # @show noi, a
    if (objp1 > obj || isnan(objp1))
        alpha = T(0)
        for i = 1:n
            x[i] = xm1[i]
        end
    end
    return alpha, noi;
end


text:
module Jrill

using HTTP
using BufferedStreams
using DataFrames
import JSON

#TODO jdbc
# TODO organize repo into a module structure
# with different functions
# TODO HTTPS
# TODO <IP address or host name>
# TODO actual query
# TODO results to Dataframe

struct DrillConnection
   host::String
   port::Int
   ssl::Bool
end

"""
	make_url(conn::DrillConnection)

Create base url for connection to Drill REST API

# Arguments
- `conn::DrillConnection`: Drill Connection info

# Returns
- base url string for HTTP requests
"""
function make_url(conn::DrillConnection)
    proto = conn.ssl ? "https" : "http"
    return "$proto://$(conn.host):$(conn.port)"
end

"""
   drill_query(conn::DrillConnection, query::String; to_df::Bool=false, limit::Int=nothing)

Submit a query and return results.

# Arguments
- `conn::DrillConnection`: Drill Connection info
- `query::String`: query to submit
- `to_df::Bool=false`: to return as Dataframe otherwise return JSON.parse of the result
- `limit::Int=nothing`: sets autoLimit in the HTTP request to limit number of rows (Drill 1.16+)

# Returns
- Dataframe or the result of JSON.parse depending on to_df
"""
function drill_query(conn::DrillConnection, query::String; to_df::Bool=false, limit::Int=-1)
    url = make_url(conn)
    io = Base.BufferStream()
    headers = ["Content-Type" => "application/json"]
    params = Dict(
                  "queryType"=>"SQL",
                  "query" => query
                 )
    if limit != -1
       params["autoLimit"] = limit
    end

    r = HTTP.request("POST", "$url/query.json",
             headers,
             JSON.json(params),
             response_stream=io)
    # TODO create Dataframe
    if to_df
       df = DataFrame()
       retval = JSON.parse(io)
       for (index,row) in enumerate(retval["rows"])
           if index == 1
	       # Create dataframe defn
	       for key in keys(row)
                  df[Symbol(key)] = row[key]
	       end
	   else
	       # push dict as new row
	       push!(df, collect(values(row)))
	   end
       end
       return df
    else
      return JSON.parse(io)
    end
end

"""
	drill_cluster(conn::DrillConnection)

Get Drillbit information, such as port numbers.

# Arguments
- `conn::DrillConnection`: Drill Connection info

# Returns
- the result of JSON.parse

"""
function drill_cluster(conn::DrillConnection)
    url = make_url(conn)
    io = Base.BufferStream()
    r = HTTP.request("GET", "$url/cluster.json", response_stream=io)
    return JSON.parse(io)
end

"""
	Get the status of Drill.

TODO
"""
#function drill_status(conn::DrillConnection)
#end

"""
	drill_metrics(conn::DrillConnection)

Get the current memory metrics.

# Arguments
- `conn::DrillConnection`: Drill Connection info

# Returns
- the result of JSON.parse

"""
function drill_metrics(conn::DrillConnection)
    url = make_url(conn)
    io = Base.BufferStream()
    r = HTTP.request("GET", "$url/status/metrics", response_stream=io)
    return JSON.parse(io)
end

"""
	drill_threads(conn::DrillConnection)

Get the status of threads.

# Arguments
- `conn::DrillConnection`: Drill Connection info

# Returns
- String of the request body
"""
function drill_threads(conn::DrillConnection)
    url = make_url(conn)
    r = HTTP.request("GET", "$url/status/threads")
    #TODO make this parasable/browsable in someway?
    return String(r.body)
end

"""
	drill_options(conn::DrillConnection)

List the name, default, and data type of the system and session options.

# Examples
```julia-repl
julia> dc = drill_connection("localhost", 8047)
julia> drill_options(drill_conn)
```
# Arguments
- `conn::DrillConnection`: Drill Connection info

# Returns
- the result of JSON.parse
"""
function drill_options(conn::DrillConnection)
    url = make_url(conn)
    io = Base.BufferStream()
    r = HTTP.request("GET", "$url/options.json", response_stream=io)
    return JSON.parse(io)
end

function drill_connection(host::String, port::Integer; ssl::Bool=false)
    return DrillConnection(host,port, ssl)
end

## Test block #
#println("hello world")
#drill_conn = drill_connection("192.168.0.20", 8047)
#query = "SELECT * FROM sys.version"
#query = "show databases"
#query = "SELECT * FROM sys.profiles_json"
#ret = drill_query(drill_conn, query, to_df=true)
#println(ret)

end # module


text:
@testset "121.best-time-to-buy-and-sell-stock.jl" begin
    @test max_profit_121([7, 6, 4, 3, 1]) == 0
    @test max_profit_121([7, 1, 5, 3, 6, 4]) == 5
    @test max_profit_121([3, 2, 6, 5, 0, 3]) == 4
end


text:
function df_to_dict_with_id(df::DataFrame)
    Dict(col[1] => Dict(zip(Symbol.(df[!, 1]), col[2])) for col in eachcol(df, true))
end

function df_to_dict(df::DataFrame)
    Dict(Symbol.(col[1]) => col[2] for col in eachcol(df, true))
end

function dictzip(df::DataFrame, which::Tuple{Symbol, Symbol})
    return Dict(zip(df[which[1]], df[which[2]]))
end

function insert_default!(dict::Dict, default, keys)
	for k in keys
		haskey(dict, k) || (dict[k] = default)
	end
end

max_issues_repo_name: bobbyxiong/redispatch-ptg

max_forks_repo_name: bobbyxiong/redispatch-ptg


text:
"""
fluxonium.jl - do some calculations for fluxonium
"""

using Dierckx
using FFTW
using HDF5
using LaTeXStrings
using LinearAlgebra
import Plots
using StaticArrays
using Zygote

WDIR = get(ENV, "RBQOC_PATH", "../../")
include(joinpath(WDIR, "src", "spin", "spin.jl"))

# Construct paths.
const EXPERIMENT_META = "spin"
const EXPERIMENT_NAME = "figures"
const SAVE_PATH = joinpath(WDIR, "out", EXPERIMENT_META, EXPERIMENT_NAME)
const DFQ_PLOT_FILE_PATH = joinpath(SAVE_PATH, "dfq.png")
const DFQ_DATA_FILE_PATH = joinpath(SAVE_PATH, "dfq.h5")
const PINK_PLOT_FILE_PATH = joinpath(SAVE_PATH, "pink.png")
const PINKFW_PLOT_FILE_PATH = joinpath(SAVE_PATH, "pinkfw.png")

## SYSTEM DEFINITION ##
# FLUXONIUM_STATE_COUNT is the state count used in the T1 calculations
const FLUXONIUM_STATE_COUNT = 110
const FLUXONIUM_LEVELS = Array(range(0., stop=FLUXONIUM_STATE_COUNT - 1,
                               length=FLUXONIUM_STATE_COUNT))
const SQRT_FLUXONIUM_LEVELS_TRUNC = map(sqrt, FLUXONIUM_LEVELS[2:FLUXONIUM_STATE_COUNT])
const ANNIHILATE = diagm(1 => SQRT_FLUXONIUM_LEVELS_TRUNC)
const CREATE = diagm(-1 => SQRT_FLUXONIUM_LEVELS_TRUNC)
const PHI_OP = PHI_OSC * 2^(-1//2) * (CREATE + ANNIHILATE)
const H_EXP_RAW = exp(1im * PHI_OP)
const H_LC = diagm(E_PLASMA * FLUXONIUM_LEVELS)

## DOMEGA FIT ##
const FBFQ_MIN = 0.3987965
const FBFQ_MAX = 1 - FBFQ_MIN
const FBFQ_SAMPLE_COUNT = Integer(1e3)
const FBFQ_SAMPLES = Array(range(FBFQ_MIN, stop=FBFQ_MAX, length=FBFQ_SAMPLE_COUNT))
const FBFQ_DFQ_POLYDEG = 10

function fbfq_hamiltonian(fbfq)
    reduced_flux = 2 * pi * fbfq
    h_exp = H_EXP_RAW * exp(1im * reduced_flux)
    h_cos = 0.5 * (h_exp + h_exp')
    h_fluxonium = real(H_LC - EJ * h_cos)
    return h_fluxonium
end

function fbfq_gphie(fbfq)
    hamiltonian = fbfq_hamiltonian(fbfq)
    evals, evecs = eigen(Hermitian(hamiltonian))
    gphie = evecs[:, 1]' * PHI_OP * evecs[:, 2]
    return gphie
end

# derivative of the 2-pi reduced amplitude w.r.t the external flux by flux quantum
const DA_DFBFQ = 4 * pi * EL * fbfq_gphie(0.5)

function fbfq_fq(fbfq)
    hamiltonian = fbfq_hamiltonian(fbfq)
    eigvals_ = eigvals(Hermitian(hamiltonian))
    fq = (eigvals_[2] - eigvals_[1])
    return fq
end

function fbfq_dfq_helin(fbfq)
    dfbfq = fbfq - 0.5
    delta = 0.014 * 2 * pi / 2
    A = 0.15 * 2 * pi
    B = 0.06
    dfq = (2 * A^2 * dfbfq / B / sqrt(A^2 * dfbfq^2 + B^2 * delta^2))
    return dfq / (2 * pi)
end

function fit_dfq(;plot=false, save=false, pull=false)
    if pull
        (fqs, dfqs) = h5open(DFQ_DATA_FILE_PATH, "r") do data_file
            fqs = read(data_file, "fqs")
            dfqs = read(data_file, "dfqs")
            return(fqs, dfqs)
        end
    else
        fqs = zeros(size(FBFQ_SAMPLES)[1])
        dfqs = zeros(size(FBFQ_SAMPLES)[1])
        dfqhs = zeros(size(FBFQ_SAMPLES)[1])
        for (i, fbfq) in enumerate(FBFQ_SAMPLES)
            fqs[i] = fbfq_fq(fbfq)
            (dfq,) = Zygote.gradient(fbfq_fq, fbfq)
            dfqs[i] = real(dfq)
            dfqhs[i] = fbfq_dfq_helin(fbfq)
        end
    end

    # fbfq_dfq_poly = Polynomials.fit(FBFQ_SAMPLES, dfqs, FBFQ_DFQ_POLYDEG)
    # dfqs_poly = map(fbfq_dfq_poly, FBFQ_SAMPLES)

    fbfq_dfq_dierckx = Dierckx.Spline1D(FBFQ_SAMPLES, dfqs)
    dfqs_dierckx = map(fbfq_dfq_dierckx, FBFQ_SAMPLES)

    if save
        h5open(DFQ_DATA_FILE_PATH, "cw") do data_file
            write(data_file, "fbfqs", FBFQ_SAMPLES)
            write(data_file, "fqs", fqs)
            write(data_file, "dfqs", dfqs)
        end
    end

    if plot
        fig = Plots.plot(dpi=DPI, legend=:bottomright)
        Plots.scatter!(fig, FBFQ_SAMPLES, fqs, label="fq", alpha=ALPHA, markersize=MS)
        Plots.scatter!(fig, FBFQ_SAMPLES, dfqs, label="dfq", alpha=ALPHA, markersize=MS)
        # Plots.plot!(fig, FBFQ_SAMPLES, dfqs_poly, label="dfq poly")
        Plots.plot!(fig, FBFQ_SAMPLES, dfqs_dierckx, label="dfq dierckx")
        Plots.plot!(fig, FBFQ_SAMPLES, dfqhs, label="dfqh")
        Plots.xlabel!(L"$\Phi / \Phi_{0}$")
        Plots.ylabel!("Amp (GHz)")
        Plots.savefig(fig, DFQ_PLOT_FILE_PATH)
    end
end

function plot_pink_noise_from_white(count; ndist=FBFQ_NDIST, dt_inv=1e2, plot=false, seed=0,
                                    namp=NAMP_PREFACTOR)
    Random.seed!(seed)
    freqs = fftfreq(count, dt_inv)
    times = (0:count-1) ./ dt_inv
    time = count / dt_inv

    white_noise = rand(ndist, count)
    pink_noise_ = Array{Complex{Float64}, 1}(white_noise)
    # transform white noise to frequency domain
    fft!(pink_noise_)
    white_fft = Array{Float64, 1}(map(abs, pink_noise_)) / count
    # square root of the spectral density is the
    # fourier transform of the noise
    # normalize by count
    for i in 2:length(pink_noise_)
        pink_noise_[i] = pink_noise_[i] / (sqrt(abs(freqs[i])) * count)
    end
    # normalize to dt_inv, this is the fft value at f=0
    pink_noise_[1] = dt_inv / count
    pink_fft = Array{Float64, 1}(map(abs, pink_noise_))
    # transform to time domain
    ifft!(pink_noise_)
    # take modulus, normalize by count
    for i = 1:length(pink_noise_)
        pink_noise_[i] = abs(pink_noise_[i]) * count
    end
    pink_noise_ = Array{Float64, 1}(pink_noise_)

    # get delta_fq for plotting
    fnoise = pink_noise_ * namp

    if plot
        # inds = 1:length(pink_fft)
        inds = 1:Int(1e4):length(pink_fft)
        subfig1 = Plots.plot()
        Plots.plot!(subfig1, freqs[inds], pink_fft[inds], label="pink")
        Plots.plot!(subfig1, freqs[inds], white_fft[inds], label="white")
        Plots.xlabel!(L"f \; \textrm{(GHz)}")
        Plots.ylabel!(L"|\hat{x}(f)| \; \textrm{(a.u.)}")
        subfig2 = Plots.plot()
        Plots.plot!(subfig2, times[inds], pink_noise_[inds], label="pink")
        Plots.plot!(subfig2, times[inds], white_noise[inds], label="white")
        Plots.xlabel!(L"t \; \textrm{(ns)}")
        Plots.ylabel!(L"x(t) \; \textrm{(a.u.)}")
        subfig3 = Plots.plot()
        Plots.plot!(subfig3, times[inds], fnoise[inds], label=:none)
        Plots.xlabel!(L"t \; \textrm{(ns)}")
        Plots.ylabel!(L"\Delta f_{q} \; \textrm{(GHz)}")
        layout = @layout [a; b; c]
        fig = Plots.plot(subfig1, subfig2, subfig3, dpi=DPI, layout=layout)
        Plots.savefig(fig, PINKFW_PLOT_FILE_PATH)
    end

    return fnoise
end

function plot_pink_noise_from_density(count; dt_inv=1e2, namp=NAMP_PREFACTOR, plot=false, seed=0,
                                      ndist=FBFQ_NDIST)
    Random.seed!(seed)
    times = Array(0:1:count-1) / dt_inv
    time = count / dt_inv
    freqs = fftfreq(count, dt_inv)

    white_noise = rand(ndist, count)
    pink_noise_ = Array{Complex{Float64}, 1}(freqs)

    pink_noise_[1] = dt_inv
    for i = 2:length(pink_noise_)
        pink_noise_[i] = 1 / sqrt(abs(pink_noise_[i]))
    end
    pink_fft = Array{Float64, 1}(pink_noise_)

    ifft!(pink_noise_)
    for i = 1:length(pink_noise_)
        pink_noise_[i] = abs(pink_noise_[i]) / count
    end
    pink_noise_ = Array{Float64, 1}(pink_noise_)

    fnoise = pink_noise * namp

    if plot
        subfig1 = Plots.plot()
        Plots.plot!(subfig1, freqs, pink_fft, label="pink")
        Plots.xlabel!(L"f \; \textrm{(GHz)}")
        Plots.ylabel!(L"|\hat{x}(f)| \; \textrm{(a.u.)}")
        subfig2 = Plots.plot()
        Plots.plot!(subfig2, times, pink_noise, label="pink")
        Plots.xlabel!(L"t \; \textrm{(ns)}")
        Plots.ylabel!(L"x(t) \; \textrm{(a.u.)}")
        subfig3 = Plots.plot()
        Plots.plot!(subfig3, times, fnoise, label=:none)
        Plots.xlabel!(L"t \; \textrm{(ns)}")
        Plots.ylabel!(L"\Delta f_{q} \; \textrm{(GHz)}")
        layout = @layout [a; b; c]
        fig = Plots.plot(subfig1, subfig2, subfig3, dpi=DPI, layout=layout)
        Plots.savefig(fig, PINK_PLOT_FILE_PATH)
    end

    return fnoise
end

function plot_pink_noise_timmer(count; dt_inv=DT_NOISE_INV, namp=NAMP_PREFACTOR,
                                plot=false, seed=0, ndist=STD_NORMAL)
    countby2 = Int(floor(count / 2))
    sqrt_count = sqrt(count)
    Random.seed!(seed)
    times = Array(0:1:count-1) / dt_inv
    time = count / dt_inv
    freqs = fftfreq(count, dt_inv)

    noise_fft = Array{Complex{Float64}}(zeros(count))
    for i = 2:countby2
        noise = (rand(ndist) + 1im * rand(ndist)) / freqs[i]^(1//2)
        noise_fft[i] = noise
        noise_fft[end - i + 2] = conj(noise)
    end
    if iseven(count)
        noise_fft[countby2 + 1] = rand(ndist) / abs(freqs[countby2 + 1])^(1//2)
    else
        noise = (rand(ndist) + 1im * rand(ndist)) / freqs[countby2 + 1]^(1//2)
        noise_fft[countby2 + 1] = noise
        noise_fft[end - (countby2 + 1) + 2] = conj(noise)
    end
    noise_fft_ = copy(noise_fft)
    fft!(noise_fft)
    pink_noise = Array{Float64}(zeros(count))
    for i = 1:count
        pink_noise[i] = real(noise_fft[i]) / sqrt_count
    end

    if plot
        subfig1 = Plots.plot()
        Plots.scatter!(subfig1, freqs, map(x -> abs(x)^2, noise_fft_), label="pink")
        Plots.xlabel!(L"f \; \textrm{(GHz)}")
        Plots.ylabel!(L"|\hat{x}(f)|^2 \; \textrm{(a.u.)}")
        subfig2 = Plots.plot()
        Plots.scatter!(subfig2, times, pink_noise, label="pink")
        Plots.xlabel!(L"t \; \textrm{(ns)}")
        Plots.ylabel!(L"x(t) \; \textrm{(a.u.)}")
        layout = @layout [a; b]
        fig = Plots.plot(subfig1, subfig2, dpi=DPI, layout=layout)
        Plots.savefig(fig, PINK_PLOT_FILE_PATH)
    end

    return pink_noise
end

function plot_pink_noise_timmer(count; dt_inv=DT_NOISE_INV, namp=NAMP_PREFACTOR,
                                plot=false, seed=0, ndist=STD_NORMAL)
    countby2 = Int(floor(count / 2))
    sqrt_count = sqrt(count)
    Random.seed!(seed)
    times = Array(0:1:count-1) / dt_inv
    time = count / dt_inv
    freqs = fftfreq(count, dt_inv)

    noise_fft = Array{Complex{Float64}}(zeros(count))
    for i = 2:countby2
        noise = (rand(ndist) + 1im * rand(ndist)) / freqs[i]^(1//2)
        noise_fft[i] = noise
        noise_fft[end - i + 2] = conj(noise)
    end
    if iseven(count)
        noise_fft[countby2 + 1] = rand(ndist) / abs(freqs[countby2 + 1])^(1//2)
    else
        noise = (rand(ndist) + 1im * rand(ndist)) / freqs[countby2 + 1]^(1//2)
        noise_fft[countby2 + 1] = noise
        noise_fft[end - (countby2 + 1) + 2] = conj(noise)
    end
    noise_fft_ = copy(noise_fft)
    fft!(noise_fft)
    pink_noise = Array{Float64}(zeros(count))
    for i = 1:count
        pink_noise[i] = real(noise_fft[i]) / sqrt_count
    end

    if plot
        subfig1 = Plots.plot()
        Plots.scatter!(subfig1, freqs, map(x -> abs(x)^2, noise_fft_), label="pink")
        Plots.xlabel!(L"f \; \textrm{(GHz)}")
        Plots.ylabel!(L"|\hat{x}(f)|^2 \; \textrm{(a.u.)}")
        subfig2 = Plots.plot()
        Plots.scatter!(subfig2, times, pink_noise, label="pink")
        Plots.xlabel!(L"t \; \textrm{(ns)}")
        Plots.ylabel!(L"x(t) \; \textrm{(a.u.)}")
        layout = @layout [a; b]
        fig = Plots.plot(subfig1, subfig2, dpi=DPI, layout=layout)
        Plots.savefig(fig, PINK_PLOT_FILE_PATH)
    end

    return pink_noise
end

function pink_noise_filter(count; dt_inv=DT_NOISE_INV, namp=NAMP_PREFACTOR,
                           ndist=STD_NORMAL, seed=0, plot=false)
    b1 = 0.049922035
    b2 = -0.095993537
    b3 = 0.050612699
    b4 = -0.004408786
    a1 = 1
    a2 = -2.494956002
    a3 = 2.017265875
    a4 = -0.522189400
    Random.seed!(seed)
    xs = rand(ndist, count)
    ys = zeros(count)
    ys[1] = b1 * xs[1]
    ys[2] = b1 * xs[2] + b2 * xs[1] - a2 * ys[1]
    ys[3] = b1 * xs[3] + b2 * xs[2] + b3 * xs[1] - a2 * ys[2] - a3 * ys[1]
    for i = 4:count
        ys[i] = (b1 * xs[i] + b2 * xs[i - 1] + b3 * xs[i - 2] + b4 * xs[i - 3]
                 - a2 * ys[i - 1] - a3 * ys[i - 2] - a4 * ys[i - 3])
    end
    ys = ys * dt_inv
    times = Array(0:1:count - 1) / dt_inv
    freqs = fftfreq(count, dt_inv)
    fmin = minimum(freqs)
    fmax = maximum(freqs)
    pl_freqs = range(fmin, stop=fmax, length=Int(3e2))
    pink_lines = [freq == 0 ? 0 : count / abs(freq) for freq in pl_freqs]
    ys_fft = fft(ys)

    if plot
        subfig1 = Plots.plot()
        Plots.scatter!(subfig1, freqs, map(x -> abs(x)^2, ys_fft))
        Plots.plot!(subfig1, pl_freqs, pink_lines)
        Plots.xlabel!(L"f \; \textrm{(GHz)}")
        Plots.ylabel!(L"|\hat{x}(f)|^2 \; \textrm{(a.u.)}")
        subfig2 = Plots.plot()
        Plots.scatter!(subfig2, times, ys)
        Plots.xlabel!(L"t \; \textrm{(ns)}")
        Plots.ylabel!(L"x(t) \; \textrm{(a.u.)}")
        layout = @layout [a; b]
        fig = Plots.plot(subfig1, subfig2, dpi=DPI, layout=layout)
        Plots.savefig(fig, PINK_PLOT_FILE_PATH)
    end

    return ys
end

max_stars_repo_name: SchusterLab/rbqoc

max_issues_repo_name: SchusterLab/rbqoc

max_forks_repo_name: SchusterLab/rbqoc


text:
@system Nitrogen begin
    SPAD: SPAD_greenness => 60 ~ preserve(parameter)
    _a: SPAD_N_coeff_a => 0.0004 ~ preserve(u"g/m^2", parameter)
    _b: SPAD_N_coeff_b => 0.012 ~ preserve(u"g/m^2", parameter)
    _c: SPAD_N_coeff_c => 0 ~ preserve(u"g/m^2", parameter)
    N(SPAD, _a, _b, _c): leaf_nitrogen_content => begin
        a*SPAD^2 + b*SPAD + c
    end ~ preserve(u"g/m^2", parameter)

    Np(N, SLA) => N * SLA ~ track(u"percent")
    SLA: specific_leaf_area => 200 ~ preserve(u"cm^2/g")
end


text:
function abseigs1d(abso::Vector{Float64},eigs::Vector{Float64},J::Matrix{Float64},
    Q::Float64,A::Float64,beta::Float64,rho::Float64,ts::Float64,xs::Float64,
    Qs::Float64,As::Float64)
    J1d(J,Q,A,beta,rho,ts,xs,Qs,As);
    eigs[1] = 0.5*(J[1,1] + J[2,2]) + sqrt.(0.25*(J[1,1] + J[2,2])^2 -
        J[1,1]*J[2,2] + J[1,2]*J[2,1]);
    eigs[2] = 0.5*(J[1,1] + J[2,2]) - sqrt.(0.25*(J[1,1] + J[2,2])^2 -
        J[1,1]*J[2,2] + J[1,2]*J[2,1]);
    abso[1] = abs(eigs[1]);
    abso[2] = abs(eigs[2]);
end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
# Laplacian

"""
    laplacian!(v,w)

Evaluate the discrete Laplacian of `w` and return it as `v`. The data `w` can be
of type dual/primary nodes or edge components or edges; `v` must be of the same type.

# Example

```jldoctest
julia> w = Nodes(Dual,(8,6));

julia> v = deepcopy(w);

julia> w[4,3] = 1.0;

julia> laplacian!(v,w)
Nodes{Dual,8,6,Float64} data
Printing in grid orientation (lower left is (1,1))
6×8 Array{Float64,2}:
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   1.0  0.0  0.0  0.0  0.0
 0.0  0.0  1.0  -4.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0   1.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0   0.0  0.0  0.0  0.0  0.0
```
"""
function laplacian!(out::Nodes{Dual,NX, NY}, w::Nodes{Dual,NX, NY}) where {NX, NY}
    view(out,2:NX-1,2:NY-1) .= view(w,2:NX-1,1:NY-2) .+ view(w,1:NX-2,2:NY-1) .+
                               view(w,3:NX,2:NY-1) .+ view(w,2:NX-1,3:NY) .-
                               4.0.*view(w,2:NX-1,2:NY-1)
    #@inbounds for y in 2:NY-1, x in 2:NX-1
    #    out[x,y] = w[x,y-1] + w[x-1,y] - 4w[x,y] + w[x+1,y] + w[x,y+1]
    #end
    out
end

function laplacian!(out::Nodes{Primal,NX, NY}, w::Nodes{Primal,NX, NY}) where {NX, NY}
    view(out,2:NX-2,2:NY-2) .= view(w,2:NX-2,1:NY-3) .+ view(w,1:NX-3,2:NY-2) .+
                               view(w,3:NX-1,2:NY-2) .+ view(w,2:NX-2,3:NY-1) .-
                               4.0.*view(w,2:NX-2,2:NY-2)
    #@inbounds for y in 2:NY-2, x in 2:NX-2
    #    out[x,y] = w[x,y-1] + w[x-1,y] - 4w[x,y] + w[x+1,y] + w[x,y+1]
    #end
    out
end

function laplacian!(out::XEdges{Dual,NX, NY}, w::XEdges{Dual,NX, NY}) where {NX, NY}
    view(out,2:NX-2,2:NY-1) .= view(w,2:NX-2,1:NY-2) .+ view(w,1:NX-3,2:NY-1) .+
                               view(w,3:NX-1,2:NY-1) .+ view(w,2:NX-2,3:NY) .-
                               4.0.*view(w,2:NX-2,2:NY-1)
    out
end

function laplacian!(out::YEdges{Dual,NX, NY}, w::YEdges{Dual,NX, NY}) where {NX, NY}
    view(out,2:NX-1,2:NY-2) .= view(w,2:NX-1,1:NY-3) .+ view(w,1:NX-2,2:NY-2) .+
                               view(w,3:NX,2:NY-2)   .+ view(w,2:NX-1,3:NY-1) .-
                               4.0.*view(w,2:NX-1,2:NY-2)
    out
end

function laplacian!(out::XEdges{Primal,NX, NY}, w::XEdges{Primal,NX, NY}) where {NX, NY}
   view(out,2:NX-1,2:NY-2) .= view(w,2:NX-1,1:NY-3) .+ view(w,1:NX-2,2:NY-2) .+
                              view(w,3:NX,2:NY-2) .+ view(w,2:NX-1,3:NY-1) .-
                              4.0.*view(w,2:NX-1,2:NY-2)
    out
end

function laplacian!(out::YEdges{Primal,NX, NY}, w::YEdges{Primal,NX, NY}) where {NX, NY}
   view(out,2:NX-2,2:NY-1) .= view(w,2:NX-2,1:NY-2) .+ view(w,1:NX-3,2:NY-1) .+
                                view(w,3:NX-1,2:NY-1) .+ view(w,2:NX-2,3:NY) .-
                               4.0.*view(w,2:NX-2,2:NY-1)
    out
end

function laplacian!(out::Edges{C,NX, NY}, w::Edges{C,NX, NY}) where {C,NX, NY}
  laplacian!(out.u,w.u)
  laplacian!(out.v,w.v)
  out
end

function laplacian!(out::EdgeGradient{C,D,NX, NY}, w::EdgeGradient{C,D,NX, NY}) where {C,D,NX, NY}
  laplacian!(out.dudx,w.dudx)
  laplacian!(out.dvdx,w.dvdx)
  laplacian!(out.dudy,w.dudy)
  laplacian!(out.dvdy,w.dvdy)
  out
end

#=
function laplacian!(out::Edges{Dual,NX, NY}, w::Edges{Dual,NX, NY}) where {NX, NY}
  view(out.u,2:NX-2,2:NY-1) .= view(w.u,2:NX-2,1:NY-2) .+ view(w.u,1:NX-3,2:NY-1) .+
                          view(w.u,3:NX-1,2:NY-1) .+ view(w.u,2:NX-2,3:NY) .-
                          4.0.*view(w.u,2:NX-2,2:NY-1)
  view(out.v,2:NX-1,2:NY-2) .= view(w.v,2:NX-1,1:NY-3) .+ view(w.v,1:NX-2,2:NY-2) .+
                               view(w.v,3:NX,2:NY-2) .+ view(w.v,2:NX-1,3:NY-1) .-
                               4.0.*view(w.v,2:NX-1,2:NY-2)
  #@inbounds for y in 2:NY-1, x in 2:NX-2
  #    out.u[x,y] = w.u[x,y-1] + w.u[x-1,y] - 4w.u[x,y] + w.u[x+1,y] + w.u[x,y+1]
  #end
  #@inbounds for y in 2:NY-2, x in 2:NX-1
  #    out.v[x,y] = w.v[x,y-1] + w.v[x-1,y] - 4w.v[x,y] + w.v[x+1,y] + w.v[x,y+1]
  #end
  out
end
=#

#=
function laplacian!(out::Edges{Primal,NX, NY}, w::Edges{Primal,NX, NY}) where {NX, NY}
  view(out.u,2:NX-1,2:NY-2) .= view(w.u,2:NX-1,1:NY-3) .+ view(w.u,1:NX-2,2:NY-2) .+
                               view(w.u,3:NX,2:NY-2) .+ view(w.u,2:NX-1,3:NY-1) .-
                              4.0.*view(w.u,2:NX-1,2:NY-2)
  view(out.v,2:NX-2,2:NY-1) .= view(w.v,2:NX-2,1:NY-2) .+ view(w.v,1:NX-3,2:NY-1) .+
                                view(w.v,3:NX-1,2:NY-1) .+ view(w.v,2:NX-2,3:NY) .-
                               4.0.*view(w.v,2:NX-2,2:NY-1)
  #@inbounds for y in 2:NY-2, x in 2:NX-1
  #    out.u[x,y] = w.u[x,y-1] + w.u[x-1,y] - 4w.u[x,y] + w.u[x+1,y] + w.u[x,y+1]
  #end
  #@inbounds for y in 2:NY-1, x in 2:NX-2
  #    out.v[x,y] = w.v[x,y-1] + w.v[x-1,y] - 4w.v[x,y] + w.v[x+1,y] + w.v[x,y+1]
  #end
  out
end
=#

"""
    laplacian(w)

Evaluate the discrete Laplacian of `w`. The data `w` can be of type dual/primary
nodes or edges. The returned result is of the same type as `w`.

# Example

```jldoctest
julia> q = Edges(Primal,(8,6));

julia> q.u[2,2] = 1.0;

julia> laplacian(q)
Edges{Primal,8,6,Float64} data
u (in grid orientation)
5×8 Array{Float64,2}:
 0.0   0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0   1.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  -4.0  1.0  0.0  0.0  0.0  0.0  0.0
 0.0   0.0  0.0  0.0  0.0  0.0  0.0  0.0
v (in grid orientation)
6×7 Array{Float64,2}:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
function laplacian(w::Nodes{C,NX,NY}) where {C<:CellType,NX,NY}
  laplacian!(Nodes(C,w), w)
end

function laplacian(w::XEdges{C,NX,NY}) where {C<:CellType,NX,NY}
  laplacian!(XEdges(C,w), w)
end

function laplacian(w::YEdges{C,NX,NY}) where {C<:CellType,NX,NY}
  laplacian!(YEdges(C,w), w)
end

function laplacian(w::Edges{C,NX,NY}) where {C<:CellType,NX,NY}
  laplacian!(Edges(C,w), w)
end

function laplacian(w::EdgeGradient{C,D,NX,NY}) where {C<:CellType,D<:CellType,NX,NY}
  laplacian!(EdgeGradient(C,w), w)
end

"""
    laplacian_symm!(v,w)

Evaluate the symmetric 5-point discrete Laplacian of `w` and return it as `v`. The data `w` can be
of type dual nodes only for now. This symmetric Laplacian also evaluates the
partial Laplacians (using only available stencil data) on the ghost nodes.
"""
function laplacian_symm!(out::Nodes{Dual,NX, NY}, w::Nodes{Dual,NX, NY}) where {NX, NY}
    @inbounds for y in 2:NY-1, x in 2:NX-1
        out[x,y] = w[x,y-1] + w[x-1,y] - 4w[x,y] + w[x+1,y] + w[x,y+1]
    end
    @inbounds for y in 2:NY-1
        out[1,y]  = w[1,y-1]            - 4w[1,y] + w[2,y] + w[1,y+1]
        out[NX,y] = w[NX,y-1] + w[NX-1,y]- 4w[NX,y]        + w[NX,y+1]
    end
    @inbounds for x in 2:NX-1
        out[x,1]  = w[x-1,1] + w[x+1,1] - 4w[x,1] + w[x,2]
        out[x,NY] = w[x-1,NY]+ w[x+1,NY]- 4w[x,NY] + w[x,NY-1]
    end
    out[1,1] = -4w[1,1] + w[1,2] + w[2,1]
    out[NX,1] = -4w[NX,1] + w[NX-1,1] + w[NX,2]
    out[1,NY] = -4w[1,NY] + w[1,NY-1] + w[2,NY]
    out[NX,NY] = -4w[NX,NY] + w[NX,NY-1] + w[NX-1,NY]
    out
end

"""
    plan_laplacian(dims::Tuple,[with_inverse=false],[fftw_flags=FFTW.ESTIMATE],
                          [factor=1.0],[dx=1.0],[dtype=Float64])

Constructor to set up an operator for evaluating the discrete Laplacian on
dual or primal nodal data of dimension `dims`. If the optional keyword
`with_inverse` is set to `true`, then it also sets up the inverse Laplacian
(the lattice Green's function, LGF). These can then be applied, respectively, with
`*` and `\\` operations on data of the appropriate size. The optional
parameter `factor` is a scalar used to multiply the result of the operator and
divide the inverse. The optional parameter
`dx` is used in adjusting the uniform value of the LGF to match the behavior
of the continuous analog at large distances; this is set to 1.0 by default. The
type of data on which to act is floating point by default, but can also be ComplexF64.
This is specified with the optional parameter `dtype`

Instead of the first argument, one can also supply `w::Nodes` to specify the
size of the domain.

# Example

```jldoctest
julia> w = Nodes(Dual,(5,5));

julia> w[3,3] = 1.0;

julia> L = plan_laplacian(5,5;with_inverse=true)
Discrete Laplacian (and inverse) on a (nx = 5, ny = 5) grid acting on Float64 data with spacing 1.0

julia> s = L\\w
Nodes{Dual,5,5,Float64} data
Printing in grid orientation (lower left is (1,1))
5×5 Array{Float64,2}:
 0.16707    0.129276     0.106037     0.129276    0.16707
 0.129276   0.0609665   -0.00734343   0.0609665   0.129276
 0.106037  -0.00734343  -0.257343    -0.00734343  0.106037
 0.129276   0.0609665   -0.00734343   0.0609665   0.129276
 0.16707    0.129276     0.106037     0.129276    0.16707

julia> L*s ≈ w
true
```
"""
function plan_laplacian end

"""
    plan_laplacian!(dims::Tuple,[with_inverse=false],[fftw_flags=FFTW.ESTIMATE],
                          [factor=1.0][,nthreads=length(Sys.cpu_info())])

Same as [`plan_laplacian`](@ref), but operates in-place on data. The number of threads `threads` defaults to the number of
logical CPU cores on the system.
"""
function plan_laplacian! end

struct Laplacian{NX, NY, T, R,inplace}
    factor::T
    dx::Float64
    conv::Union{CircularConvolution{NX, NY, T},Nothing}
end

Base.eltype(::Laplacian{NX,NY,T}) where {NX,NY,T} = T

for (lf,inplace) in ((:plan_laplacian,false),
                     (:plan_laplacian!,true))
    @eval function $lf(dims::Tuple{Int,Int};
                   with_inverse = false, fftw_flags = FFTW.ESTIMATE, factor::Real = 1.0, dx = 1.0, dtype = Float64, nthreads = length(Sys.cpu_info()))
        NX, NY = dims
        if !with_inverse
            return Laplacian{NX, NY, dtype, false, $inplace}(convert(dtype,factor),convert(Float64,dx),nothing)
        end

        G = view(LGF_TABLE, 1:NX, 1:NY)
        Laplacian{NX, NY, dtype, true, $inplace}(convert(dtype,factor),convert(Float64,dx),CircularConvolution(G, fftw_flags,dtype=dtype,nthreads=nthreads))
    end

    @eval function $lf(nx::Int, ny::Int;
        with_inverse = false, fftw_flags = FFTW.ESTIMATE, factor = 1.0, dx = 1.0, dtype = Float64, nthreads = length(Sys.cpu_info()))
        $lf((nx, ny), with_inverse = with_inverse, fftw_flags = fftw_flags, factor = factor, dx = dx, dtype = dtype, nthreads = nthreads)
    end

    @eval function $lf(w::ScalarGridData;
        with_inverse = false, fftw_flags = FFTW.ESTIMATE, factor = 1.0, dx = 1.0, dtype = Float64, nthreads = length(Sys.cpu_info())) where {T<:CellType,NX,NY}
        $lf(size(w), with_inverse = with_inverse, fftw_flags = fftw_flags, factor = factor, dx = dx, dtype = dtype, nthreads = nthreads)
    end
end

function Base.show(io::IO, L::Laplacian{NX, NY, T, R, inplace}) where {NX, NY, T, R, inplace}
    nodedims = "(nx = $NX, ny = $NY)"
    inverse = R ? " (and inverse)" : ""
    isinplace = inplace ? " in-place" : ""
    print(io, "Discrete$isinplace Laplacian$inverse on a $nodedims grid acting on $T data with
               factor $(L.factor) and spacing $(L.dx)")
end

mul!(out::T, L::Laplacian, s::T) where T<:GridData = (laplacian!(out, s); out .*= L.factor)

*(L::Laplacian{MX,MY,T,R,false}, s::GridData) where {MX,MY,T,R} =
      L.factor*laplacian(s)

function (*)(L::Laplacian{MX,MY,T,R,true}, s::GridData) where {MX,MY,T,R}
    mul!(s,L,deepcopy(s))
end

#=
function ldiv!(out::Nodes{C,NX, NY,T},
                   L::Laplacian{MX, MY, T, true, inplace},
                   s::Nodes{C, NX, NY,T}) where {C <: CellType, NX, NY, MX, MY, T, inplace}

    mul!(out.data, L.conv, s.data)
    inv_factor = 1.0/L.factor

    # Adjust the behavior at large distance to match continuous kernel
    out.data .-= (sum(s.data)/2π)*(GAMMA+log(8)/2-log(L.dx))
    out.data .*= inv_factor
    out
end
=#
for (datatype) in (:Nodes, :XEdges, :YEdges)
  @eval function ldiv!(out::$datatype{C,NX, NY,T},
                   L::Laplacian{MX, MY, T, true, inplace},
                   s::$datatype{C, NX, NY,T}) where {C <: CellType, NX, NY, MX, MY, T, inplace}

    mul!(out.data, L.conv, s.data)
    inv_factor = 1.0/L.factor

    # Adjust the behavior at large distance to match continuous kernel
    out.data .-= (sum(s.data)/2π)*(GAMMA+log(8)/2-log(L.dx))
    out.data .*= inv_factor
    out
  end

  #@eval \(L::Laplacian{MX,MY,T,R,false},s::$datatype{C,NX,NY}) where {MX,MY,T,R,C <: CellType,NX,NY} =
  #  ldiv!($datatype(C,s), L, s)

end

function ldiv!(out::Edges{C,NX,NY},L::Laplacian,s::Edges{C,NX,NY}) where {C,NX,NY}
  ldiv!(out.u,L,s.u)
  ldiv!(out.v,L,s.v)
  out
end

function ldiv!(out::EdgeGradient{C,D,NX,NY},L::Laplacian,s::EdgeGradient{C,D,NX,NY}) where {C,D,NX,NY}
  ldiv!(out.dudx,L,s.dudx)
  ldiv!(out.dvdx,L,s.dvdx)
  ldiv!(out.dudy,L,s.dudy)
  ldiv!(out.dvdy,L,s.dvdy)
  out
end

\(L::Laplacian{MX,MY,T,R,false},s::G) where {MX,MY,T,R,G<:GridData} =
  ldiv!(G(), L, s)

\(L::Laplacian{MX,MY,T,R,true},s::GridData{NX,NY}) where {MX,MY,T,R,NX,NY} =
    ldiv!(s, L, deepcopy(s))

#=
\(L::Laplacian{MX,MY,T,R,false},s::Nodes{C,NX,NY}) where {MX,MY,T,R,C <: CellType,NX,NY} =
  ldiv!(Nodes(C,s), L, s)

\(L::Laplacian{MX,MY,T,R,true},s::Nodes{C,NX,NY}) where {MX,MY,T,R,C <: CellType,NX,NY} =
  ldiv!(s, L, deepcopy(s))
=#


text:
module CODATA

using PhysicalConstant, Unitful, Measurements

import PhysicalConstants: @constant, @derived_constant

@constant(m_e, "Electron mass", 9.109_383_56e-31,
          BigFloat(910_938_356)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 1.1e-38,
          BigFloat(11)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA ")
@constant(m_h, "Hydrogen mass", 1.673_723_600e-27,
          BigFloat(1673_723_600)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 2.1e-35,
          BigFloat(21)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA ")
@constant(m_n, "Neutron mass", 1.674_927_471e-27,
          BigFloat(1674_927_471)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 2.1e-35,
          BigFloat(21)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA ")
@constant(m_p, "Proton mass", 1.672_621_898e-27,
          BigFloat(1672_621_898)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 2.1e-35,
          BigFloat(21)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA ")
@constant(m_m, "Muon mass", 1.836_152_67e-28,
          BigFloat(1836_152_67)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 2.1e-35,
          BigFloat(21)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA ")
@constant(m_t, "Tau mass", 3.167_735_502e-27,
          BigFloat(3167_735_502)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 2.1e-35,
          BigFloat(21)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA 2018")
@constant(m_u, "Atomic mass constant", 1.660_539_040e-27,
          BigFloat(1660_539_040)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          u"kg", 2.0e-35,
          BigFloat(20)/BigFloat(1000_000_000_000_000_000_000_000_000_000_000_000),
          "CODATA ")

end


text:
@testset "h-for" begin
    F.def_LOCAL_VARS!()
    s = """
        @def list = [1,2,3]
        """ |> fd2html_td
    hs = raw"""
        ABC
        {{for x in list}}
          {{fill x}}
        {{end}}
        """
    tokens = F.find_tokens(hs, F.HTML_TOKENS, F.HTML_1C_TOKENS)
    hblocks, tokens = F.find_all_ocblocks(tokens, F.HTML_OCB)
    qblocks = F.qualify_html_hblocks(hblocks)
    @test qblocks[1] isa F.HFor
    @test qblocks[1].vname == "x"
    @test qblocks[1].iname == "list"
    @test qblocks[2] isa F.HFun
    @test qblocks[3] isa F.HEnd

    content, head, i = F.process_html_for(hs, qblocks, 1)
    @test isapproxstr(content, "1 2 3")
end

@testset "h-for2" begin
    F.def_LOCAL_VARS!()
    s = """
        @def list = ["path/to/badge1.png", "path/to/badge2.png"]
        """ |> fd2html_td
    h = raw"""
        ABC
        {{for x in list}}
            {{fill x}}
        {{end}}
        """ |> F.convert_html
    @test isapproxstr(h, """
        ABC
        path/to/badge1.png
        path/to/badge2.png
        """)
end

@testset "h-for3" begin
    F.def_LOCAL_VARS!()
    s = """
        @def iter = (("a", 1), ("b", 2), ("c", 3))
        """ |> fd2html_td
    h = raw"""
        ABC
        {{for (n, v) in iter}}
            name:{{fill n}}
            value:{{fill v}}
        {{end}}
        """ |> F.convert_html
    @test isapproxstr(h, """
        ABC
        name:a
        value:1
        name:b
        value:2
        name:c
        value:3
        """)

    s = """
        @def iter2 = ("a"=>10, "b"=>7, "c"=>3)
        """ |> fd2html_td
    h = raw"""
        ABC
        {{for (n, v) in iter2}}
            name:{{fill n}}
            value:{{fill v}}
        {{end}}
        """ |> F.convert_html
    @test isapproxstr(h, """
        ABC
        name:a
        value:10
        name:b
        value:7
        name:c
        value:3
        """)
end

# read from file, see FranklinFAQ-001
@testset "for-file" begin
    gotd()
    write("members.csv", """
    name,github
    Eric Mill,konklone
    Parker Moore,parkr
    Liu Fengyun,liufengyun
    """)
    s = """
        @def members = eachrow(readdlm("members.csv", ',', skipstart=1))
        ~~~

        {{for (name, alias) in members}}

            {{name}}

        {{end}}

        ~~~
        """ |> fd2html
    @test isapproxstr(s, """

                    Eric Mill

                    Parker Moore

                    Liu Fengyun

""")
end


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(csc(c+d*x)^3*(A+A*csc(c+d*x))*(a+a*csc(c+d*x)), x) == :(-2*A*a*d^-1*cot(c+d*x)+-7//8*A*a*d^-1*arctanh(cos(c+d*x))+-2//3*A*a*d^-1*cot(c+d*x)^3+-7//8*A*a*d^-1*cot(c+d*x)*csc(c+d*x)+-1//4*A*a*d^-1*csc(c+d*x)^3*cot(c+d*x))
@test integrate(csc(c+d*x)^2*(A+A*csc(c+d*x))*(a+a*csc(c+d*x)), x) == :(-1*A*a*d^-1*arctanh(cos(c+d*x))+-5//3*A*a*d^-1*cot(c+d*x)+-1*A*a*d^-1*cot(c+d*x)*csc(c+d*x)+-1//3*A*a*d^-1*csc(c+d*x)^2*cot(c+d*x))
@test integrate((A+A*csc(c+d*x))*(a+a*csc(c+d*x))*csc(c+d*x), x) == :(-2*A*a*d^-1*cot(c+d*x)+-3//2*A*a*d^-1*arctanh(cos(c+d*x))+-1//2*A*a*d^-1*cot(c+d*x)*csc(c+d*x))
@test integrate((A+A*csc(c+d*x))*(a+a*csc(c+d*x))*sin(c+d*x), x) == :(2*A*a*x+-1*A*a*d^-1*arctanh(cos(c+d*x))+-1*A*a*d^-1*cos(c+d*x))
@test integrate(sin(c+d*x)^2*(A+A*csc(c+d*x))*(a+a*csc(c+d*x)), x) == :(3//2*A*a*x+-2*A*a*d^-1*cos(c+d*x)+-1//2*A*a*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate(sin(c+d*x)^3*(A+A*csc(c+d*x))*(a+a*csc(c+d*x)), x) == :(A*a*x+-2*A*a*d^-1*cos(c+d*x)+1//3*A*a*d^-1*cos(c+d*x)^3+-1*A*a*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate(csc(c+d*x)^3*(A+A*csc(c+d*x))*(a+-1*a*csc(c+d*x)), x) == :(-1//8*A*a*d^-1*arctanh(cos(c+d*x))+-1//8*A*a*d^-1*cot(c+d*x)*csc(c+d*x)+1//4*A*a*d^-1*csc(c+d*x)^3*cot(c+d*x))
@test integrate(csc(c+d*x)^2*(A+A*csc(c+d*x))*(a+-1*a*csc(c+d*x)), x) == :(1//3*A*a*d^-1*cot(c+d*x)^3)
@test integrate((A+A*csc(c+d*x))*(a+-1*a*csc(c+d*x))*csc(c+d*x), x) == :(-1//2*A*a*d^-1*arctanh(cos(c+d*x))+(1/2)*A*a*d^-1*cot(c+d*x)*csc(c+d*x))
@test integrate((A+A*csc(c+d*x))*(a+-1*a*csc(c+d*x))*sin(c+d*x), x) == :(A*a*d^-1*arctanh(cos(c+d*x))+-1*A*a*d^-1*cos(c+d*x))
@test integrate(sin(c+d*x)^2*(A+A*csc(c+d*x))*(a+-1*a*csc(c+d*x)), x) == :(-1//2*A*a*x+-1//2*A*a*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate(sin(c+d*x)^3*(A+A*csc(c+d*x))*(a+-1*a*csc(c+d*x)), x) == :(1//3*A*a*d^-1*cos(c+d*x)^3)
@test integrate(csc(c+d*x)^3*(A+-1*A*csc(c+d*x))*(a+a*csc(c+d*x)), x) == :(-1//8*A*a*d^-1*arctanh(cos(c+d*x))+-1//8*A*a*d^-1*cot(c+d*x)*csc(c+d*x)+1//4*A*a*d^-1*csc(c+d*x)^3*cot(c+d*x))
@test integrate(csc(c+d*x)^2*(A+-1*A*csc(c+d*x))*(a+a*csc(c+d*x)), x) == :(1//3*A*a*d^-1*cot(c+d*x)^3)
@test integrate((A+-1*A*csc(c+d*x))*(a+a*csc(c+d*x))*csc(c+d*x), x) == :(-1//2*A*a*d^-1*arctanh(cos(c+d*x))+(1/2)*A*a*d^-1*cot(c+d*x)*csc(c+d*x))
@test integrate(csc(c+d*x)^-1*(A+-1*A*csc(c+d*x))*(a+a*csc(c+d*x)), x) == :(A*a*d^-1*arctanh(cos(c+d*x))+-1*A*a*d^-1*cos(c+d*x))
@test integrate(csc(c+d*x)^-2*(A+-1*A*csc(c+d*x))*(a+a*csc(c+d*x)), x) == :(-1//2*A*a*x+-1//2*A*a*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate(csc(c+d*x)^-3*(A+-1*A*csc(c+d*x))*(a+a*csc(c+d*x)), x) == :(1//3*A*a*d^-1*cos(c+d*x)^3)
@test integrate(csc(c+d*x)^3*(A+-1*A*csc(c+d*x))*(a+-1*a*csc(c+d*x)), x) == :(2*A*a*d^-1*cot(c+d*x)+-7//8*A*a*d^-1*arctanh(cos(c+d*x))+2//3*A*a*d^-1*cot(c+d*x)^3+-7//8*A*a*d^-1*cot(c+d*x)*csc(c+d*x)+-1//4*A*a*d^-1*csc(c+d*x)^3*cot(c+d*x))
@test integrate(csc(c+d*x)^2*(A+-1*A*csc(c+d*x))*(a+-1*a*csc(c+d*x)), x) == :(A*a*d^-1*arctanh(cos(c+d*x))+-5//3*A*a*d^-1*cot(c+d*x)+A*a*d^-1*cot(c+d*x)*csc(c+d*x)+-1//3*A*a*d^-1*csc(c+d*x)^2*cot(c+d*x))
@test integrate((A+-1*A*csc(c+d*x))*(a+-1*a*csc(c+d*x))*csc(c+d*x), x) == :(2*A*a*d^-1*cot(c+d*x)+-3//2*A*a*d^-1*arctanh(cos(c+d*x))+-1//2*A*a*d^-1*cot(c+d*x)*csc(c+d*x))
@test integrate(csc(c+d*x)^-1*(A+-1*A*csc(c+d*x))*(a+-1*a*csc(c+d*x)), x) == :(-2*A*a*x+-1*A*a*d^-1*arctanh(cos(c+d*x))+-1*A*a*d^-1*cos(c+d*x))
@test integrate(csc(c+d*x)^-2*(A+-1*A*csc(c+d*x))*(a+-1*a*csc(c+d*x)), x) == :(3//2*A*a*x+2*A*a*d^-1*cos(c+d*x)+-1//2*A*a*d^-1*cos(c+d*x)*sin(c+d*x))
@test integrate(csc(c+d*x)^-3*(A+-1*A*csc(c+d*x))*(a+-1*a*csc(c+d*x)), x) == :(-1*A*a*x+-2*A*a*d^-1*cos(c+d*x)+1//3*A*a*d^-1*cos(c+d*x)^3+A*a*d^-1*cos(c+d*x)*sin(c+d*x))

max_stars_repo_path: test/4 Trig functions/4.6 Cosecant/4.6.3.1 (a+b csc)^m (d csc)^n (A+B csc).jl


text:
module GeoRegions

## Modules Used
using Dates
using DelimitedFiles
using GeometryBasics
using Logging
using PolygonOps
using PrettyTables

import Base: show, read

## Exporting the following functions:
export
        GeoRegion, RectRegion, PolyRegion, RegionGrid, RectGrid, PolyGrid,
        resetGeoRegions, templateGeoRegions, listGeoRegions,
        isGeoRegion, addGeoRegions, removeGeoRegion, coordGeoRegion, isinGeoRegion,
        Point2

## Abstract types
"""
    GeoRegion

Abstract supertype for geographical regions.
"""
abstract type GeoRegion end

"""
    RectRegion{ST<:AbstractString, FT<:Real} <: GeoRegion

Structure containing information on a rectilinear (i.e. rectangular shape on a lon-lat grid) region `regID`, with parent ID `parID`, name `name`, and bounds `N`, `S`, `E`, `W`.

If `W` < 0º, then `is180` is true, and if `E` > 180, then `is360` is true.  Both `is180` and `is360` can simultaneously true, but they cannot both be false.
"""
struct RectRegion{ST<:AbstractString, FT<:Real} <: GeoRegion
    regID :: ST
    parID :: ST
    name  :: ST
    N     :: FT
    S     :: FT
    E     :: FT
    W     :: FT
    is180 :: Bool
    is360 :: Bool
end

"""
    PolyRegion{ST<:AbstractString, FT<:Real} <: GeoRegion

Structure containing information on a polygonal (and non-rectangular) GeoRegion with ID `regID`, with parent ID `parID`, name `name`, and rectangular bounds `N`, `S`, `E`, `W`.  The actual shape within the rectangular bound is given by `shape`.

If `W` < 0º, then `is180` is true, and if `E` > 180, then `is360` is true.  Both `is180` and `is360` can simultaneously true, but they cannot both be false.
"""
struct PolyRegion{ST<:AbstractString, FT<:Real} <: GeoRegion
    regID :: ST
    parID :: ST
    name  :: ST
    N     :: FT
    S     :: FT
    E     :: FT
    W     :: FT
    shape :: Vector{Point2{FT}}
    is180 :: Bool
    is360 :: Bool
end

"""
    RegionGrid

Abstract supertype for geographical regions.
"""
abstract type RegionGrid end

"""
    RectGrid{FT<:Real}

Structure containing information needed to extract data from a given lon-lat grid for a given RectRegion.
"""
struct RectGrid{FT<:Real} <: RegionGrid
    igrid :: Vector{Int}
    ilon  :: Vector{Int}
    ilat  :: Vector{Int}
    glon  :: Vector{FT}
    glat  :: Vector{FT}
end

"""
    PolyGrid{FT<:Real}

Structure containing information needed to extract data from a given lon-lat grid for a given PolyRegion with numeric elements of type `FT`. In addition to the information stored by the RectGrid type, it also has a `mask` that allows for non-rectilinear shape extraction.
"""
struct PolyGrid{FT<:Real} <: RegionGrid
    igrid :: Vector{Int}
    ilon  :: Vector{Int}
    ilat  :: Vector{Int}
    glon  :: Vector{FT}
    glat  :: Vector{FT}
    mask  :: Array{FT,2}
end

## Including other files in the module
include("Read.jl")
include("Create.jl")
include("Query.jl")
include("IsIn.jl")
include("IsInGeoRegion.jl")
include("Extract.jl")
include("Show.jl")

end # module


text:
using NeuralVerification
using Test

macro no_error(ex)
    quote
        try $(esc(ex))
            true
        catch
            false
        end
    end
end

at = @__DIR__

# Small MNIST Network

mnist_small = read_nnet("$at/../examples/networks/mnist_small.nnet")

# entry 23 in MNIST datset
input_center = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,121,254,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13,230,253,248,99,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,118,253,253,225,42,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61,253,253,253,74,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,32,206,253,253,186,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,211,253,253,239,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,254,253,253,133,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,142,255,253,186,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,149,229,254,207,21,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,54,229,253,254,105,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,152,254,254,213,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,112,251,253,253,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,29,212,253,250,149,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36,214,253,253,137,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,75,253,253,253,59,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,93,253,253,189,17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,224,253,253,84,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43,235,253,126,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,99,248,253,119,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,225,235,49,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
output_center = [-1311.1257826380004,4633.767704436501,-654.0718535670002,-1325.349417307,1175.2361184373997,-1897.8607293569007,-470.3405972940001,830.8337987382,-377.7467076115001,572.3674015264198]

in_epsilon = 1 #0-255
out_epsilon = 10 #logit domain

input_low = input_center .- in_epsilon
input_high = input_center .+ in_epsilon

output_low = output_center .- out_epsilon
output_high = output_center .+ out_epsilon

inputSet = Hyperrectangle(low=input_low, high=input_high)
outputSet = Hyperrectangle(low=output_low, high=output_high)

problem_hyperrect_small = Problem(mnist_small, inputSet, outputSet)
solver_reluVal = ReluVal(max_iter = 2)

@test @no_error solve(solver_reluVal, problem_hyperrect_small)

solver_reluplex=Reluplex()
@test @no_error solve(solver_reluplex, problem_hyperrect_small)

# Deep MNIST network
mnist_large = read_nnet("$at/../examples/networks/mnist_large.nnet")

#mnist_large = read_nnet("/Users/christopherlazarus/Desktop/jsand/NNet/readNNet/python/mnist_large.nnet")

# entry 23 in MNIST datset
input_center = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,121,254,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13,230,253,248,99,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,118,253,253,225,42,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61,253,253,253,74,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,32,206,253,253,186,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,211,253,253,239,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,254,253,253,133,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,142,255,253,186,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,149,229,254,207,21,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,54,229,253,254,105,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,152,254,254,213,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,112,251,253,253,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,29,212,253,250,149,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36,214,253,253,137,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,75,253,253,253,59,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,93,253,253,189,17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,224,253,253,84,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43,235,253,126,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,99,248,253,119,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,225,235,49,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
output_center = [131.8781755,134.8987015,141.6166255,158.34307,129.8803525,104.8286425,98.64196,133.6358395,131.1716215,105.10621]

in_epsilon = 1 #0-255
out_epsilon = 10 #logit domain

input_low = input_center .- in_epsilon
input_high = input_center .+ in_epsilon

output_low = output_center .- out_epsilon
output_high = output_center .+ out_epsilon

inputSet = Hyperrectangle(low=input_low, high=input_high)
outputSet = Hyperrectangle(low=output_low, high=output_high)

problem_hyperrect_deep = Problem(mnist_large, inputSet, outputSet)
solver_reluVal = ReluVal(max_iter = 2)

@test @no_error solve(solver_reluVal, problem_hyperrect_deep)

solver_reluplex=Reluplex()
@test @no_error solve(solver_reluplex, problem_hyperrect_deep)

# Wide MNIST network

mnist_wide = read_nnet("$at/../examples/networks/mnist-1-100.nnet")
input_center = [0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,121,254,136,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,13,230,253,248,99,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,4,118,253,253,225,42,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,61,253,253,253,74,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,32,206,253,253,186,9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,211,253,253,239,69,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,254,253,253,133,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,142,255,253,186,8,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,149,229,254,207,21,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,54,229,253,254,105,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,152,254,254,213,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,112,251,253,253,26,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,29,212,253,250,149,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,36,214,253,253,137,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,75,253,253,253,59,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,93,253,253,189,17,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,224,253,253,84,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,43,235,253,126,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,99,248,253,119,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,225,235,49,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0]
output_center = [-1219.6781446362043,2713.019467032775,-3826.9579099907123,-3306.084501236369,2641.9250628243003,-1829.0475478735848,-2477.0765438475523,84.03007306734946,-2072.196780418711,145.0876643197917]

in_epsilon = 1 #0-255
out_epsilon = 10 #logit domain

input_low = input_center .- in_epsilon
input_high = input_center .+ in_epsilon

output_low = output_center .- out_epsilon
output_high = output_center .+ out_epsilon

inputSet = Hyperrectangle(low=input_low, high=input_high)
outputSet = Hyperrectangle(low=output_low, high=output_high)

problem_hyperrect_wide = Problem(mnist_wide, inputSet, outputSet)
solver_reluVal = ReluVal(max_iter = 2)

@test @no_error solve(solver_reluVal, problem_hyperrect_wide)

solver_reluplex=Reluplex()
@test @no_error solve(solver_reluplex, problem_hyperrect_wide)


text:
@testset "Pezz" begin
    # @test pezz(buzz(fizz(7))) == "Pezz"
    @test 7 |> fizz |> buzz |> pezz == "Pezz"
    @test 21 |> fizz |> buzz |> pezz == "FizzPezz"
    @test 35 |> fizz |> buzz |> pezz == "BuzzPezz"
    @test 105 |> fizz |> buzz |> pezz == "FizzBuzzPezz"
    @test 105 |> fizz |> pezz |> buzz == "FizzPezzBuzz"
    @test 105 |> pezz |> buzz |> fizz == "PezzBuzzFizz"
    # 以下も念のため
    @test 1 |> fizz |> buzz |> pezz == 1
    @test 3 |> fizz |> buzz |> pezz == "Fizz"
    @test 5 |> fizz |> buzz |> pezz == "Buzz"
    @test 15 |> fizz |> buzz |> pezz == "FizzBuzz"
    @test 15 |> buzz |> fizz |> pezz == "BuzzFizz"
    @test 104 |> fizz |> buzz |> pezz == 104
end


text:
#=
  Package: FishEBM
  File: FileIO.jl
  Devin Rose
  Functions for summarizing the results of a simulation
  Created: June 2016
=#

"""
  INPUT: ageDataFrame = DataFrame of yearly adult age-specific population.
  OUTPUT: ageSUMMARY.csv: file containing yearly age-specific population levels.

  Last update: August 2016
"""
function ageData(adf::DataFrame, path::String)
  separateDirChar = getDirChar()

  file = string(path,"$(separateDirChar)ageSpecificSUMMARY.csv")
  writetable(file, adf);

  file = string(path,"$(separateDirChar)yearlyAgeSpecificSUMMARY.csv")
  yearlyData = DataFrame(Year = 0, Age2 = 0, Age3 = 0, Age4 = 0, Age5 = 0, Age6 = 0, Age7 = 0, Age8Plus = 0, Total = 0)
  writetable(file, getYearlyPop(adf, yearlyData))

end

"""
  INPUT: popDataFrame = DataFrame of weekly stage population.
  OUTPUT: simSUMMARY.csv: file containing weekly population levels.

  Last update: August 2016
"""
function aliveData(popDataFrame::DataFrame, path::String)
  separateDirChar = getDirChar()

  file = string(path,"$(separateDirChar)stageSUMMARY.csv")
  writetable(file, popDataFrame)

  file = string(path,"$(separateDirChar)yearlyStageSUMMARY.csv")
  yearlyData = DataFrame(Year = 0, Stage1 = 0, Stage2 = 0, Stage3 = 0, Stage4 = 0, Total = 0)
  writetable(file, getYearlyPop(popDataFrame, yearlyData))
end

"""
  Description: Creates a directory based on the user defined path and has
    dependency on exists. It is reconmended that one passes isdir(path) into
    exists.

  Returns: N/A

  Last update: June 2016
"""
function createDir(path::String, exists::Bool)
  if exists == true
    print("A directory already exists in the path $path.\n")
  else
    mkdir(path)
  end
end

"""
  INPUT:  path = user specified path.
  OUTPUT: datePath = dir. for date folder
          & dirExist = True/False.
"""
function dateDir(path::String)
  currentDate = string(Dates.today())
  datePath = string(path,string(getDirChar(), currentDate))
  dirExist = isdir(datePath)

  return (datePath,dirExist)
end

"""
  Description: Finds the character used to separate the path based on current OS

  Returns: String

  Last update: January 2017
"""
function getDirChar()
  @assert(is_windows() || is_unix(), "There is currently no functionality for
  the operating system: "$(Sys.KERNEL)", now aborting.")

  if is_windows() # For a windows file system
    return "\\"
  elseif is_unix() # For apple and linux distributions
    return "/"
  end
end

"""
  Description: Generic function to sum weekly data of dataframes into yearly data

  Last update: August 2016
"""
function getYearlyData(weeklyData::DataFrame, yearlyDataFrame::DataFrame)
  numYears = ceil((size(weeklyData)[1] - 1)/52)
  weekMin = 1
  weekMax = 52

  for year = 1:numYears
    yearData = weeklyData[(weeklyData[:Week] .>= weekMin)&(weeklyData[:Week] .<= weekMax), :]

    #Use size of weekly dataframe - 2 to account for the week and total columns that aren't needed here
    dataVector = fill(0, (size(weeklyData)[2] - 2))
    for i = 1:length(dataVector)
      dataVector[i] = sum(yearData[i+1])
    end
    push!(yearlyDataFrame, vcat(year, dataVector..., sum(dataVector)))
    weekMin = weekMin + 52
    weekMax = weekMax + 52
  end

  return yearlyDataFrame
end

"""
  Description: Generic function to return yearly snapshot of population

  Last update: August 2016
"""
function getYearlyPop(weeklyData::DataFrame, yearlyDataFrame::DataFrame)
  numYears = ceil((size(weeklyData)[1] - 1)/52)
  week = 1

  for year = 1:numYears
    yearData = weeklyData[(weeklyData[:Week] .== week), :]

    #Use size of weekly dataframe - 2 to account for the week and total columns that aren't needed here
    dataVector = fill(0, (size(weeklyData)[2] - 2))
    for i = 1:length(dataVector)
      dataVector[i] = sum(yearData[i+1])
    end
    push!(yearlyDataFrame, vcat(year, dataVector..., sum(dataVector)))
    week = week + 52
  end

  return yearlyDataFrame
end

"""
  INPUT: hdf = DataFrame of weekly age-specific harvest levels and total harvest.
  OUTPUT: harvestSUMMARY.csv: file containing weekly harvest levels.

  Last update: August 2016
"""
function harvestData(hdf::DataFrame, zoneData::DataFrame, path::String)
  #Output weekly harvest by age
  file = string(path,"$(getDirChar())harvestSUMMARY.csv")
  writetable(file, hdf)

  #Output weekly harvest by zone
  file = string(path,"$(getDirChar())harvestZoneSUMMARY.csv")
  writetable(file, zoneData)

  #Output yearly harvest by age
  file = string(path,"$(getDirChar())yearlyHarvestSUMMARY.csv")
  yearHarvest = DataFrame(Year = 0, Age2 = 0, Age3 = 0, Age4 = 0, Age5 = 0, Age6 = 0, Age7 = 0, Age8Plus = 0, Total = 0)
  writetable(file, getYearlyData(hdf, yearHarvest))

  #Output yearly harvest by zone
  file = string(path,"$(getDirChar())yearlyHarvestZoneSUMMARY.csv")
  yearZoneHarvest = DataFrame(Year = 0, z1 = 0, z2 = 0, z3 = 0, z4 = 0, z5 = 0, z6 = 0, z7 = 0, z8 = 0, z9 = 0, z10 = 0,
                              z11 = 0, z12 = 0, z13 = 0, z14 = 0, z15 = 0, z16 = 0, z17 = 0, z18 = 0, Total = 0)
  writetable(file, getYearlyData(zoneData, yearZoneHarvest))

end

"""
  INPUT: kdf = DataFrame of weekly killed data by natural and extra mortality.
  OUTPUT: killedSUMMARY.csv: file containing weekly mortality levels.

  Last update: August 2016
"""
function killedData(kdf::DataFrame, path::String)
  file = string(path,"$(getDirChar())killedSUMMARY.csv")

  #Sum natural, extra, and Compensatory mortalities before outputting to .csv
  for i = 1:size(kdf)[1]
    totalKill = 0
    for j = 2:4
      totalKill += kdf[i,j]
    end #for j
    kdf[i,5] = totalKill
  end #for i
  #Output weekly mortalities
  writetable(file, kdf)

  #Output yearly mortalities
  file = string(path,"$(getDirChar())yearlyKilledSUMMARY.csv")
  yearMortality = DataFrame(Year = 0, Natural = 0, Extra = 0, Compensatory = 0, Total = 0)
  writetable(file, getYearlyData(kdf, yearMortality))
end

"""
  INPUT:  path = user specified path.

  OUTPUT: resultsPath = dir. for results folder
    & dirExist = True/False.
"""
function resultsDir(path::String)
    resultsPath = string(path, "$(getDirChar())results")
    dirExist = isdir(resultsPath)

    return (resultsPath,dirExist)
end

"""
  INPUT:  path = user specified path.
  OUTPUT: newRunPath = dir. for new run folder,
          currentRunPath = dir. for current run folder,
          & dirExist = True/False.
"""
function runDir(path::String)
  dirCounter = 0
  stopCriteria = 0

  while (stopCriteria == 0)
    newRunPath = string(path,string("$(getDirChar())run_", dirCounter))
    dirExist = isdir(newRunPath)

    if dirExist == false
      stopCriteria = 1
    elseif dirExist == true
      dirCounter = dirCounter + 1
    end
  end

  if dirCounter == 0
    currentRunPath = string(path,string("$(getDirChar())run_", dirCounter))
  else
    currentRunPath = string(path,string("$(getDirChar())run_", dirCounter-1))
  end

  newRunPath = string(path,string("$(getDirChar())run_", dirCounter))
  dirExist = isdir(newRunPath)

  return (newRunPath, currentRunPath, dirExist)
end

"""
    Description: Finds the project directory based on the by splitting the
      string of the current directory pathway at the first and second intstance
      of '\\' (in reverse order).

    Returns: String

    Last update: June 2016
"""
function setProjPath()
  projPath = Base.source_path() #find current project path
  separateDirChar = getDirChar()
  pathEndSearch = rsearch(projPath, separateDirChar)[1]
  endOfDirPath = projPath[pathEndSearch:length(projPath)]
  projPath = split(projPath, endOfDirPath)[1]
  pathEndSearch = rsearch(projPath, separateDirChar)[1]
  endOfDirPath = projPath[pathEndSearch:length(projPath)]
  projPath = string(split(projPath, endOfDirPath)[1], "$(getDirChar())simulations")

  return ascii(projPath) #returns the path to FishEBM.jl
end

"""
  Description: Creates the newest directory for the current simulation.

  Returns: N/A

  Last update: June 2016
"""
function simDir()

  path = setProjPath()
  results = resultsDir(path)
  date = dateDir(results[1])
  run = runDir(date[1])

  createDir(results[1],results[2])
  createDir(date[1],date[2])
  createDir(run[1],run[3])
end

"""
  INPUT: runDir, userInput, k, effort, initStock, adultAssumpt,
    agentAssumpt.
  OUTPUT: simREADME.txt

    Description :Function creats a formatted textfile containing information
      particular to a simulation's run. The file saves in the same directory as
      the current date's run number. For descriptions of the INPUT's please
      refer to types.jl. Note, runDir and userInput are user given PATH and
      description respectively.

  Last update: September 2016
"""

function simReadme(adultAssumpt::AdultAssumptions, agentAssumpt::AgentAssumptions, effort::Vector{Float64}, initStock::Vector{Int64}, carryingCap::Vector{Float64}, path::String, userInput::String)

  file_name = string(path,"$(getDirChar())simREADME.txt")
  output_file = open(file_name, "w")

  # Version number
  write(output_file,"-------------------------\n")
  write(output_file, "FishEBM.jl Version Number: "*FishEBM.VersionNumber*"\n-------------------------\n")

  # Simulation description
  write(output_file,"-------------------------\n")
  write(output_file, "Simulation Description: \n-------------------------\n")
  write(output_file, userInput)
  write(output_file, "\n")

  write(output_file,"\n-------------------------\n")
  write(output_file, "General Assumptions: \n-------------------------\n")

  carryingCap_string = string("Carrying capacity: ",carryingCap,"\n")
  write(output_file, carryingCap_string)

  effort_string = "Effort vector: "
  write(output_file, effort_string)
  show(output_file, effort)
  write(output_file, "\n")

  initStock_string = "Initial stock values are: "
  write(output_file,initStock_string)
  show(output_file, initStock)
  write(output_file, "\n")

  write(output_file,"\n-------------------------\n")
  adultAssumpt_string = "Adult Assumptions: \n-------------------------\n"
  write(output_file,adultAssumpt_string)

  write(output_file,"Natural Mortality vector: ")
  show(output_file,adultAssumpt.naturalmortality)
  write(output_file, "\n")

  halfmature_string = string("Age at half maturity: ",adultAssumpt.halfmature,"\n")
  write(output_file,halfmature_string)

  write(output_file,"Broodsize (Age specific fecundity) vector: ")
  show(output_file,adultAssumpt.broodsize)
  write(output_file, "\n")

  fc_string = string("Fecundity compensation value: ",adultAssumpt.fecunditycompensation,"\n")
  write(output_file,fc_string)

  matC_string = string("Maturity compensation value: ",adultAssumpt.maturitycompensation,"\n")
  write(output_file,matC_string)

  morC_string = string("Mortality compensation value: ",adultAssumpt.mortalitycompensation,"\n")
  write(output_file,morC_string)

  catch_string = "Catchability: "
  write(output_file, catch_string)
  show(output_file,adultAssumpt.catchability)
  write(output_file, "\n")

  write(output_file,"\n-------------------------\n")
  agentAssumpt_string = "Agent Assumptions: \n-------------------------\n"
  write(output_file,agentAssumpt_string)

  natMor_string = "Natural mortality array: \n"
  write(output_file, natMor_string)
  show(output_file, agentAssumpt.naturalmortality)
  write(output_file, "\n")

  extraMor_string = "Extra mortality vector: "
  write(output_file, extraMor_string)
  show(output_file, agentAssumpt.extramortality)
  write(output_file, "\n")

  growth_string = "Growth vector: "
  write(output_file, growth_string)
  show(output_file, agentAssumpt.growth)
  write(output_file, "\n")

  move_string = "Movement matricies: \n"
  write(output_file, move_string)
  show(output_file, agentAssumpt.movement)
  write(output_file, "\n")

  auto_string = "Autonomy vector: "
  write(output_file, auto_string)
  show(output_file, agentAssumpt.autonomy)

  close(output_file)
end

"""
  INPUT: sdf = DataFrame of weekly age-specific spawn levels and total spawn size.
  OUTPUT: spawnSUMMARY.csv: file containing weekly spawn levels.

  Last update: August 2016
"""
function spawnData(sdf::DataFrame, path::String)
  #Output weekly spawn summary
  file = string(path,"$(getDirChar())spawnSUMMARY.csv")
  writetable(file, sdf)

  #Output yearly spawn summary
  file = string(path,"$(getDirChar())yearlySpawnSUMMARY.csv")
  yearSpawn = DataFrame(Year = 0, Age2 = 0, Age3 = 0, Age4 = 0, Age5 = 0, Age6 = 0, Age7 = 0, Age8Plus = 0, Total = 0)
  writetable(file, getYearlyData(sdf, yearSpawn))
end


text:
using Catalyst, DelaySSAToolkit
# 0 -->A1, => A2, ..., => A5 => 0 with delay 1 for each cascade delay

rn = @reaction_network begin
    5, 0 --> A1
end

delay_trigger_affect! = []
chain_len = 10
delay_trigger_affect! = function (integrator, rng)
    append!(integrator.de_chan[1], 1.)
end

u0 = zeros(chain_len)
de_chan0 = []
for _ in 1:chain_len
    push!(de_chan0, [])
end
tspan = (0.,50.)
delay_complete_affect! = []
for i in 1:chain_len-1
    push!(delay_complete_affect!, function (integrator, rng)
        integrator.u[i] -= 1 # A_prev minus 1
        integrator.u[i+1] += 1 # A plus 1
        append!(integrator.de_chan[i+1], 1.) # add to the delay channel
    end
    )
end
push!(delay_complete_affect!, function (integrator, rng)
    integrator.u[chain_len] -= 1 # A_prev minus 1
end)

delay_trigger = Dict(Pair(1, delay_trigger_affect!))
delay_complete = Dict(i=>delay_complete_affect![i] for i in 1:chain_len)
delay_interrupt = Dict()

delayjumpset = DelayJumpSet(delay_trigger, delay_complete, delay_interrupt)
jumpsys = convert(JumpSystem, rn, combinatoric_ratelaws = false)
dprob = DiscreteProblem(jumpsys, u0, tspan)

algo_list =[DelayDirect(),DelayRejection(),DelayMNRM(),DelayDirectCR()]

meanlist=[]
medianlist=[]
minlist=[]
maxlist=[]
stdlist=[]
for algo in algo_list
    djprob = DelayJumpProblem(jumpsys,dprob,algo,delayjumpset, de_chan0, save_positions=(false,false))
    a=@benchmark solve(djprob, SSAStepper())
    push!(meanlist,copy(mean(a).time/1e9))
    push!(medianlist,copy(median(a).time/1e9))
    push!(minlist,copy(minimum(a).time/1e9))
    push!(maxlist,copy(maximum(a).time/1e9))
    push!(stdlist,copy(std(a).time/1e9))
    print(djprob.aggregator," OK","\n")
end

df=DataFrame(mean=meanlist,median=medianlist,min=minlist,max=maxlist,std=stdlist)
CSV.write("test/zxytest/results/cascade_save_positions=f.csv",df)

df=CSV.read("test/zxytest/results/cascade_save_positions=f.csv",DataFrame)
medianlist=df.median

medianvalue=[string(round(mt,digits=6),"s") for mt in medianlist]

algo_name = ["DelayDirect","DelayRejection","DelayMNRM","DelayDirectCR"]
using Plots
p1=bar(algo_name,medianlist,legend=:false,title="cascade",ylabel="mediantime")
scatter!(algo_name, 0.0005 .+ medianlist , markeralpha=0, series_annotations=medianvalue)
savefig(p1,"test/zxytest/results/cascade_median.png")


text:
"""
      Quadrotor

      Dynamics model from "Trajectory generation and control for precise
      aggressive maneuvers with quadrotors". Orientation represented with
      Modified Rodrigues Parameters.
"""

struct Quadrotor{I, T} <: Model{I, T}
      n::Int
      m::Int
      d::Int

      mass # mass
      J       # inertia matrix
      Jinv    # inertia matrix inverse
      g       # gravity
      L    # length
      kf   # coefficient
      km   # coefficient
end

function f(model::Quadrotor, z, u, w)
      # states
      x = view(z,1:3)
      r = view(z,4:6)
      v = view(z,7:9)
      ω = view(z,10:12)

      # controls
      w1 = u[1]
      w2 = u[2]
      w3 = u[3]
      w4 = u[4]

      # forces
      F1 = model.kf * w1
      F2 = model.kf * w2
      F3 = model.kf * w3
      F4 = model.kf * w4

      F = @SVector [0.0,
                    0.0,
                    F1 + F2 + F3 + F4] # total rotor force in body frame

      # moments
      M1 = model.km * w1
      M2 = model.km * w2
      M3 = model.km * w3
      M4 = model.km * w4

      τ = @SVector [model.L * (F2 - F4),
                    model.L * (F3 - F1),
                    (M1 - M2 + M3 - M4)] # total rotor torque in body frame

      SVector{12}([v;
                   0.25 * ((1.0 - r' * r) * ω - 2.0 * cross(ω, r) + 2.0*(ω' * r) * r);
                   model.g + (1.0 / model.mass) * MRP(r[1], r[2], r[3]) * F;
                   model.Jinv * (τ - cross(ω, model.J * ω))])
end

n, m, d = 12, 4, 12

mass = 0.5
J = Diagonal(@SVector[0.0023, 0.0023, 0.004])
Jinv = Diagonal(@SVector[1.0 / 0.0023, 1.0 / 0.0023, 1.0 / 0.004])
L = 0.175
kf = 1.0
km = 0.0245

model = Quadrotor{Midpoint, FixedTime}(n, m, d,
                  mass,
                  J,
                  Jinv,
                  @SVector[0.0, 0.0, -9.81],
                  L,
                  kf,
                  km)

function visualize!(vis, p::Quadrotor, q; Δt = 0.1)
    default_background!(vis)

    obj_path = joinpath(pwd(),
      "models/quadrotor/drone.obj")
    mtl_path = joinpath(pwd(),
      "models/quadrotor/drone.mtl")

    ctm = ModifiedMeshFileObject(obj_path, mtl_path, scale = 1.0)
    setobject!(vis["drone"], ctm)
    settransform!(vis["drone"], LinearMap(RotZ(pi) * RotX(pi / 2.0)))

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

    for t = 1:length(q)
        MeshCat.atframe(anim, t) do
            settransform!(vis["drone"],
                  compose(Translation(q[t][1:3]),
                        LinearMap(MRP(q[t][4:6]...) * RotX(pi / 2.0))))
        end
    end
    # settransform!(vis["/Cameras/default"], compose(Translation(0.0, 0.0, 0.0),
    # LinearMap(RotZ(pi/2))))
    MeshCat.setanimation!(vis, anim)
end

max_issues_repo_name: thowell/DirectMotionPlanning