text:
using LinearAlgebra: norm
using OrdinaryDiffEq
using Trixi

# This elixir calculates the sound generated by a pair of co-rotating vortices in a quiescent
# medium.
# Note that this elixir can take several hours to run. For optimal results consider increasing the
# initial_refinement_level to 8.
#
# References:
# - Duck Joo Lee, Sam Ok Koo (1995)
#   Numerical study of sound generation due to a spinning vortex
#   https://doi.org/10.2514/3.12327
# - Michael Schlottke-Lakemper (2017)
#   A direct-hybrid method for aeroacoustic analysis
#   https://doi.org/10.18154/RWTH-2017-04082
module VortexPairSetup

using LinearAlgebra: norm
using Trixi

# Parameters that describe the co-rotating vortex pair
struct VortexPair{RealT<:Real}
  r0::RealT # Distance between origin and each vortex center
  rc::RealT # Vortex core radius
  c0::RealT # Speed of sound
  circulation::RealT # Circulation of the vortices
  rho0::RealT # Density
end

# Analytical flow solution, used for the initial condition of the flow simulation
function velocity(x, t, vortex_pair::VortexPair)
  @unpack r0, rc, circulation = vortex_pair

  omega = circulation / (4 * pi * r0^2)
  si, co = sincos(omega * t)
  b = SVector(r0 * co, r0 * si) # vortex centers are b and -b
  z_plus = x - b
  z_minus = x + b

  # Transform to polar coordinates
  r_plus = norm(z_plus)
  r_minus = norm(z_minus)
  theta_plus = atan(z_plus[2], z_plus[1])
  theta_minus = atan(z_minus[2], z_minus[1])

  si_plus, co_plus = sincos(theta_plus)
  si_minus, co_minus = sincos(theta_minus)

  v1 = -circulation/(2 * pi) * ( r_plus /(rc^2 + r_plus^2)  * si_plus +
                                r_minus/(rc^2 + r_minus^2) * si_minus)
  v2 = circulation/(2 * pi) * ( r_plus /(rc^2 + r_plus^2)  * co_plus +
                               r_minus/(rc^2 + r_minus^2) * co_minus )

  return SVector(v1, v2)
end

# Initial condition of the flow simulation. Uses constant density rho0 and analytical velocities.
# The pressure is calculated using the given speed of sound c0 and Bernoulli's principle
struct InitialCondition{RealT<:Real}
  vortex_pair::VortexPair{RealT}
end

function (initial_condition::InitialCondition)(x, t, equations::CompressibleEulerEquations2D)
  @unpack vortex_pair = initial_condition
  @unpack rho0, c0 = vortex_pair
  gamma = equations.gamma

  v = velocity(x, t, vortex_pair)
  p0 = rho0 * c0^2 / gamma
  p = p0 - 0.5 * (gamma-1)/gamma * sum(v.^2) # Bernoulli's principle

  prim = SVector(rho0, v[1], v[2], p)
  return prim2cons(prim, equations)
end

# For both the flow and acoustics solvers, a sponge layer is used to dampen the density
# and pressure towards the freestream values (for the flow solver) and the perturbed pressure
# to zero (for the acoustics solver).
struct SpongeLayer{RealT<:Real, uEltype<:Real, N, SourceTerms}
  sponge_layer_min::NTuple{4, RealT} # (-x,+x,-y,+y) min coordinates of sponge layer per direction
  sponge_layer_max::NTuple{4, RealT} # (-x,+x,-y,+y) max coordinates of sponge layer per direction
  reference_values::NTuple{N, uEltype} # reference values for variables affected by sponge layer
  source_terms::SourceTerms # source terms to be used outside the sponge zone
end

function SpongeLayer(; sponge_layer_min, sponge_layer_max, reference_values, source_terms=nothing)
  return SpongeLayer(sponge_layer_min, sponge_layer_max, reference_values, source_terms)
end

function (sponge_layer::SpongeLayer)(u, x, t, equations)
  @unpack sponge_layer_min, sponge_layer_max, reference_values, source_terms = sponge_layer

  if lies_in_sponge_layer(x, sponge_layer_min, sponge_layer_max)
    return source_term_sponge_layer(u, x, t, equations, sponge_layer_min, sponge_layer_max,
                                    reference_values)
  elseif source_terms !== nothing
    return source_terms(u, x, t, equations)
  else
    return SVector(ntuple(v -> zero(eltype(u)), Val(nvariables(equations))))
  end
end

function lies_in_sponge_layer(x, sponge_layer_min, sponge_layer_max)
  return (sponge_layer_min[1] <= x[1] <= sponge_layer_max[1]) || # -x direction
         (sponge_layer_min[2] <= x[1] <= sponge_layer_max[2]) || # +x direction
         (sponge_layer_min[3] <= x[2] <= sponge_layer_max[3]) || # -y direction
         (sponge_layer_min[4] <= x[2] <= sponge_layer_max[4])    # +y direction
end

function source_term_sponge_layer(u, x, t, equations::AcousticPerturbationEquations2D,
                                  sponge_layer_min::NTuple{4}, sponge_layer_max::NTuple{4},
                                  reference_values)
  # Perturbed pressure source is -alpha^2 * (u - reference_value) where alpha in [0,1] is a damping
  # factor depending on the position inside the sponge layer

  # Calculate the damping factors for each direction (=0 if x is not in the corresponding sponge
  # zone) and take the maximum. This ensures proper damping if x lies in a corner where the sponge
  # zones for two directions overlap
  alphas = ntuple(i -> calc_damping_factor(x, i, sponge_layer_min, sponge_layer_max),
                                           Val(2*ndims(equations)))
  alpha_square = maximum(alphas)^2

  return SVector(0, 0, -alpha_square*(u[3] - reference_values[1]/u[6]^2), 0, 0, 0, 0)
end

function source_term_sponge_layer(u, x, t, equations::CompressibleEulerEquations2D,
                                  sponge_layer_min::NTuple{4}, sponge_layer_max::NTuple{4},
                                  reference_values)
  # Calculate the damping factors for each direction (=0 if x is not in the corresponding sponge
  # zone) and take the maximum. This ensures proper damping if x lies in a corner where the sponge
  # zones for two directions overlap
  alphas = ntuple(i -> calc_damping_factor(x, i, sponge_layer_min, sponge_layer_max),
                                           Val(2*ndims(equations)))
  alpha_square = maximum(alphas)^2

  u_prim = cons2prim(u, equations)
  s = SVector(-alpha_square*(u_prim[1] - reference_values[1]), 0, 0,
              -alpha_square*(u_prim[4] - reference_values[2]))

  return prim2cons(s, equations)
end

function calc_damping_factor(x, direction, sponge_layer_min, sponge_layer_max)
  # Damping factor alpha grows linearly from 0 to 1 depending on how deep x lies in the sponge layer
  # If x does not lie in the sponge layer, this returns 0

  # Get the coordinate that determines how deep we are in the sponge zone
  if direction in (1, 2)
    pos = x[1]
  else
    pos = x[2]
  end

  # Determine where the sponge layer begins/ends to allow calculating the damping factor
  if iseven(direction)
    sponge_begin = sponge_layer_min[direction]
    sponge_end = sponge_layer_max[direction]
  else
    sponge_begin = sponge_layer_max[direction]
    sponge_end = sponge_layer_min[direction]
  end

  alpha = (pos - sponge_begin) / (sponge_end - sponge_begin)

  # alpha lies in [0, 1] if and only if x lies in the sponge zone
  if 0 <= alpha <= 1
    return alpha
  else
    return zero(alpha)
  end
end

# Boundary condition for the flow problem: The sponge layer dampens density and pressure towards the
# freestream values. The freestream values (converted into conservative variables) are therefore
# used as a Dirichlet boundary
struct BoundaryCondition{uEltype}
  rho::uEltype
  rho_e::uEltype
end

function (bc::BoundaryCondition)(u_inner, orientation, direction, x, t,
                                 surface_flux_function,
                                 equations::CompressibleEulerEquations2D)
  u_boundary = SVector(bc.rho, zero(bc.rho), zero(bc.rho), bc.rho_e)

  # Calculate boundary flux
  if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
  end

  return flux
end

end # module

import .VortexPairSetup

###############################################################################
# shared parameters, mesh and solver for both semidiscretizations

# Parameters of the vortex pair
Mach = 1/9
c0 = 1.0
r0 = 1.0
circulation = 4 * pi * r0 * c0 * Mach
rho = 1.0

rc = 2/9 * r0 * 1.0

T_r = 8 * pi^2 * r0^2 / circulation # Rotational period of the vortex pair
T_a = T_r / 2 # Acoustic period of the vortex pair

vortex_pair = VortexPairSetup.VortexPair(r0, rc, c0, circulation, rho)

# Shared mesh for both semidiscretizations
coordinates_min = (-135*r0, -135*r0) # minimum coordinates (min(x), min(y))
coordinates_max = ( 135*r0,  135*r0) # maximum coordinates (max(x), max(y))
refinement_patches = (
  (type="sphere", center=(0.0, 0.0), radius=85.0*r0),
  (type="sphere", center=(0.0, 0.0), radius=20.0*r0),
  (type="sphere", center=(0.0, 0.0), radius=10.0*r0),
  (type="sphere", center=(0.0, 0.0), radius=5.0*r0)
)
initial_refinement_level=7
n_cells_max=500_000
mesh = TreeMesh(coordinates_min, coordinates_max,
                initial_refinement_level=initial_refinement_level,
                refinement_patches=refinement_patches,
                n_cells_max=n_cells_max, # set maximum capacity of tree data structure
                periodicity=false)

# Create DG solver with polynomial degree = 3 and (local) Lax-Friedrichs/Rusanov flux as surface flux
solver = DGSEM(polydeg=3, surface_flux=flux_lax_friedrichs)

###############################################################################
# semidiscretization Euler equations

gamma = 1.4
equations_euler = CompressibleEulerEquations2D(gamma)

initial_condition_euler = VortexPairSetup.InitialCondition(vortex_pair)

sponge_layer_euler = VortexPairSetup.SpongeLayer(sponge_layer_min=(-135*r0, 115*r0, -135*r0, 115*r0),
                                                 sponge_layer_max=(-115*r0, 135*r0, -115*r0, 135*r0),
                                                 reference_values=(rho, rho * c0^2 / gamma)) # (rho0, p0)

boundary_condition_euler = VortexPairSetup.BoundaryCondition(rho, (rho * c0^2 / gamma) / (gamma-1))

semi_euler = SemidiscretizationHyperbolic(mesh, equations_euler, initial_condition_euler, solver,
                                          boundary_conditions=boundary_condition_euler,
                                          source_terms=sponge_layer_euler)

###############################################################################
# semidiscretization acoustic perturbation equations
equations_acoustics = AcousticPerturbationEquations2D(v_mean_global=(13.0, 26.0), c_mean_global=39.0,
                                                      rho_mean_global=52.0) # global mean values will be overwritten

sponge_layer_acoustics = VortexPairSetup.SpongeLayer(sponge_layer_min=(-135*r0, 100*r0, -135*r0, 100*r0),
                                                     sponge_layer_max=(-100*r0, 135*r0, -100*r0, 135*r0),
                                                     reference_values=(0.0,))

"""
    boundary_condition_zero(u_inner, orientation, direction, x, t, surface_flux_function,
                            equations::AcousticPerturbationEquations2D)

Boundary condition that uses a boundary state where the state variables are zero and the mean
variables are the same as in `u_inner`.
"""
function boundary_condition_zero(u_inner, orientation, direction, x, t, surface_flux_function,
                                 equations::AcousticPerturbationEquations2D)
  value = zero(eltype(u_inner))
  u_boundary = SVector(value, value, value, cons2mean(u_inner, equations)...)

  # Calculate boundary flux
  if iseven(direction) # u_inner is "left" of boundary, u_boundary is "right" of boundary
    flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
  else # u_boundary is "left" of boundary, u_inner is "right" of boundary
    flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
  end

  return flux
end

semi_acoustics = SemidiscretizationHyperbolic(mesh, equations_acoustics, initial_condition_constant,
                                              solver, boundary_conditions=boundary_condition_zero,
                                              source_terms=sponge_layer_acoustics)

###############################################################################
# ODE solvers, callbacks etc. for averaging the flow field

# Create ODE problem
tspan1 = (0.0, 400.0)
ode_averaging = semidiscretize(semi_euler, tspan1)

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

analysis_interval = 5000
alive_callback = AliveCallback(analysis_interval=analysis_interval)

tspan_averaging = (50.0, 400.0)
averaging_callback = AveragingCallback(semi_euler, tspan_averaging)

stepsize_callback = StepsizeCallback(cfl=0.8)

callbacks_averaging = CallbackSet(summary_callback, alive_callback, averaging_callback,
                                  stepsize_callback)

###############################################################################
# run simulation for averaging the flow field

# OrdinaryDiffEq's `solve` method evolves the solution in time and executes the passed callbacks
sol_averaging = solve(ode_averaging, CarpenterKennedy2N54(williamson_condition=false),
                      dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
                      save_everystep=false, callback=callbacks_averaging);

# Print the timer summary
summary_callback()

###############################################################################
# set up coupled semidiscretization

source_region(x) = sum(abs2, x) < 6.0^2 # calculate sources within radius 6 around origin
# gradually reduce acoustic source term amplitudes to zero, starting at radius 5
weights(x) = sum(abs2, x) < 5.0^2 ? 1.0 : cospi(0.5 * (norm(x) - 5.0))

semi = SemidiscretizationEulerAcoustics(semi_acoustics, semi_euler,
                                        source_region=source_region, weights=weights)

###############################################################################
# ODE solvers, callbacks etc. for the coupled simulation

# Create ODE problem
tspan = (0.0, 7.0*T_a)
ode = semidiscretize(semi, tspan)
# We need an additional ODE for the pure flow problem
ode_euler = semidiscretize(semi.semi_euler, tspan)

# Set up coupling callback
cfl_acoustics = 0.8
cfl_euler = 0.8
euler_acoustics_coupling = EulerAcousticsCouplingCallback(
  ode_euler, "out/averaging.h5", CarpenterKennedy2N54(williamson_condition=false),
  cfl_acoustics, cfl_euler, callback=SaveRestartCallback(interval=2300, output_directory="out/euler/"))

# At the beginning of the main loop, the SummaryCallback prints a summary of the simulation setup
# and resets the timers
summary_callback = SummaryCallback()

analysis_interval = 5000
analysis_callback = AnalysisCallback(semi, interval=analysis_interval)
alive_callback = AliveCallback(analysis_interval=analysis_interval)

output_directory = "out/"
save_solution = SaveSolutionCallback(interval=2300, output_directory=output_directory)
save_restart = SaveRestartCallback(interval=2300, output_directory=output_directory)

callbacks = CallbackSet(summary_callback, alive_callback, analysis_callback, save_solution,
                        save_restart, euler_acoustics_coupling)

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);

# Print the timer summary
summary_callback()


text:
for OP in (:checked_abs, :checked_neg)
    @eval begin
       @inline function $OP(x::T) where T<:SafeInteger
            ix = baseint(x)
            result = $OP(ix)
            return safeint(result)
        end
    end
end

for OP in (:checked_add, :checked_sub,
           :checked_mul, :checked_div, :checked_fld, :checked_cld,
           :checked_rem, :checked_mod)
    @eval begin
       @inline function $OP(x::T, y::T) where T<:SafeSigned
            ix = baseint(x)
            iy = baseint(y)
            result = $OP(ix, iy)
            return safeint(result)
        end

       @inline function $OP(x::T, y::T) where T<:SafeUnsigned
            ix = baseint(x)
            iy = baseint(y)
            result = $OP(ix, iy)
            return safeint(result)
        end

        @inline function $OP(x::T1, y::T2) where {T1<:SafeSigned, T2<:SafeSigned}
            xx, yy = promote(x, y)
            ix = baseint(xx)
            iy = baseint(yy)
            result = $OP(ix, iy)
            return safeint(result)
        end

        @inline function $OP(x::T1, y::T2) where {T1<:SafeUnsigned, T2<:SafeUnsigned}
            xx, yy = promote(x, y)
            ix = baseint(xx)
            iy = baseint(yy)
            result = $OP(ix, iy)
            return safeint(result)
        end

        @inline function $OP(x::T1, y::T2) where {T1<:SafeInteger, T2<:Base.BitInteger}
            xx, yy = promote(x, y)
            return $OP(xx, yy)
        end

        @inline function $OP(x::T1, y::T2) where {T1<:Base.BitInteger, T2<:SafeInteger}
            xx, yy = promote(x, y)
            return $OP(xx, yy)
        end
    end
end

for OP in (:add_with_overflow, :sub_with_overflow, :mul_with_overflow)
    @eval begin
       @inline function $OP(x::T, y::T) where T<:SafeSigned
            ix = baseint(x)
            iy = baseint(y)
            value, bool = $OP(ix, iy)
            return safeint(value), bool
        end

        @inline function $OP(x::T, y::T) where T<:SafeUnsigned
            ix = baseint(x)
            iy = baseint(y)
            value, bool = $OP(ix, iy)
            return safeint(value), bool
        end

        @inline function $OP(x::T1, y::T2) where {T1<:SafeSigned, T2<:SafeSigned}
            xx, yy = promote(x, y)
            ix = baseint(xx)
            iy = baseint(yy)
            value, bool = $OP(ix, iy)
            return safeint(value), bool
        end

        @inline function $OP(x::T1, y::T2) where {T1<:SafeUnsigned, T2<:SafeUnsigned}
            xx, yy = promote(x, y)
            ix = baseint(xx)
            iy = baseint(yy)
            value, bool = $OP(ix, iy)
            return safeint(value), bool
        end

    end
end

ext: jl


text:
function Base.append!(
    krg::Kriging{N,M},
    x::AbstractMatrix,
    y::AbstractMatrix,
) where {N,M}
    xs = krg.xscaler * x
    ys = krg.yscaler * y
    for i in 1:M
        append!(krg.gps[i], xs, ys[i, :])
    end
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ROCmLLVM_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ROCmLLVM")
JLLWrappers.@generate_main_file("ROCmLLVM", UUID("866f66b9-32a5-5756-8490-10a13b5dc930"))
end  # module ROCmLLVM_jll


text:
using JLD2

mutable struct TestType1
    x::Int
end
mutable struct TestType2
    x::Int
end
struct TestType3
    x::TestType2
end

mutable struct TestType4
    x::Int
end
mutable struct TestType5
    x::TestType4
end
mutable struct TestType6 end
primitive type TestType7 8 end
struct TestType8
    a::TestType4
    b::TestType5
    c::TestType6
    d::TestType7
end
struct TestTypeContainer{T}
    a::T
    b::Int
end
primitive type TestType9 16 end
struct TestType10
    a::Int
    b::UInt8
end
struct TestType11
    a::Int
    b::UInt8
end
struct TestType12
    x::String
end
mutable struct TestType13
    a
    TestType13() = new()
end
struct TestType14{T}
    x::T
    y::Int
end
struct TestType15{T,S}
    x::T
    y::S
end
struct TestType16{T}
    x::T
    y::Int
end
struct TestTypeContainer2{T}
    a::T
    b::Int
end
struct TestTypeContainer3{T,S}
    a::T
    b::S
end
struct TestTypeContainer4{T}
    a::T
    b::Int
end
struct TestTypeContainer5{T,S}
    a::T
    b::S
end
struct TestType17
    x::Int
end
struct TestType18
    x::Int
end
struct TestType19{T}
    x::T
end
primitive type TestType20 16 end
struct TestType21 end

fn = joinpath(tempdir(),"test.jld")
file = jldopen(fn, "w")
write(file, "x1", TestType1(57))
write(file, "x2", TestType3(TestType2(1)))
write(file, "x3", TestType4(1))
write(file, "x4", TestType5(TestType4(2)))
write(file, "x5", [TestType5(TestType4(i)) for i = 1:5])
write(file, "x6", TestType6())
write(file, "x7", reinterpret(TestType7, 0x77))
write(file, "x8", TestType8(TestType4(2), TestType5(TestType4(3)),
                            TestType6(), reinterpret(TestType7, 0x12)))
write(file, "x9", (TestType4(1),
                   (TestType5(TestType4(2)),
                    [TestType5(TestType4(i)) for i = 1:5]),
                   TestType6()))
write(file, "x10", TestTypeContainer(TestType4(3), 4))
write(file, "x11", reinterpret(TestType9, 0x1234))
write(file, "x12", TestType10(1234, 0x56))
write(file, "x13", TestType11(78910, 0x11))
write(file, "x14", TestType12("abcdefg"))
write(file, "x15", TestType13())
write(file, "x16", TestType14(1.2345, 67))
write(file, "x17", TestType15(8.91011, 12))
write(file, "x18", TestType16(12.131415, 1617))
write(file, "x19", TestTypeContainer2(TestType4(3), 4))
write(file, "x20", TestTypeContainer3(TestType4(3), 4))
write(file, "x21", TestTypeContainer4(TestType4(3), 4))
write(file, "x22", TestTypeContainer5(TestType4(3), 4))
write(file, "x23", TestType17(5678910))
write(file, "x24", Dict(TestType4(1) => TestType5(TestType4(2))))
write(file, "x25", Dict(TestType17(1) => reinterpret(TestType9, 0x4567)))
write(file, "x26", TestType18(1337))
write(file, "x27", TestType19(31337))
write(file, "x28", reinterpret(TestType20, 0x1337))
write(file, "x29", TestType21())

close(file)

workspace()

# workspace doesn't work anymore unless we call @eval Core.Main afterwards. Unfortunately, we can't
# just put all this code in a block, because then it won't be evaluated at top-level. So
# we need a bunch of @eval Core.Mains here.
@eval Core.Main using .LastMain.JLD2, Compat, Test
@eval Core.Main mutable struct TestType1
    x::Float64
end
@eval Core.Main mutable struct TestType2
    x::Int
end
@eval Core.Main struct TestType3
    x::TestType1
end
@eval Core.Main struct TestTypeContainer{T}
    a::T
    b::Int
end
@eval Core.Main primitive type TestType9 8 end
@eval Core.Main struct TestType10
    a::Int
    b::UInt8
    c::UInt8
end
@eval Core.Main struct TestType11
    b::UInt8
end
@eval Core.Main struct TestType12
    x::Int
end
@eval Core.Main mutable struct TestType13
    a
end
@eval Core.Main struct TestType14{T,S}
    x::T
    y::S
end
@eval Core.Main struct TestType15{T}
    x::T
    y::Float64
end
@eval Core.Main struct TestType16{T<:Integer}
    x::T
    y::Int
end
@eval Core.Main struct TestTypeContainer2{T}
    a::T
    b::String
end
@eval Core.Main struct TestTypeContainer3{T,S}
    a::T
    b::S
end
@eval Core.Main struct TestTypeContainer4{T,S}
    a::T
    b::S
end
@eval Core.Main struct TestTypeContainer5{T}
    a::T
    b::Int
end
@eval Core.Main abstract type TestType18 end
@eval Core.Main struct TestType19 end
@eval Core.Main struct TestType20
    x::UInt16
end
@eval Core.Main struct TestType21
    x::Int
end
@eval Core.Main begin
const TestType17 = 5

file = jldopen(LastMain.fn, "r")
x = read(file, "x1")
@test isa(x, TestType1)
@test x.x === 57.0
@test_throws MethodError read(file, "x2")
println("The following missing type warnings are a sign of normal operation.")
@test read(file, "x3").x == 1
@test read(file, "x4").x.x == 2

x = read(file, "x5")
for i = 1:5
    @test x[i].x.x == i
end
@test isempty(fieldnames(typeof(read(file, "x6"))))
@test reinterpret(UInt8, read(file, "x7")) == 0x77

x = read(file, "x8")
@test x.a.x == 2
@test x.b.x.x == 3
@test isempty(fieldnames(typeof(x.c)))
@test reinterpret(UInt8, x.d) == 0x12

x = read(file, "x9")
@test isa(x, Tuple)
@test length(x) == 3
@test x[1].x == 1
@test isa(x[2], Tuple)
@test length(x[2]) == 2
@test x[2][1].x.x == 2
for i = 1:5
    @test x[2][2][i].x.x == i
end
@test isempty(fieldnames(typeof(x[3])))

x = read(file, "x10")
@test isa(x, TestTypeContainer{Any})
@test x.a.x === 3
@test x.b === 4

x = read(file, "x11")
@test !isa(x, TestType9)
@test reinterpret(UInt16, x) === 0x1234

x = read(file, "x12")
@test !isa(x, TestType10)
@test x.a === 1234
@test x.b === 0x56

x = read(file, "x13")
@test isa(x, TestType11)
@test x.b === 0x11

x = read(file, "x14")
@test !isa(x, TestType12)
@test x.x == "abcdefg"

@test_throws JLD2.UndefinedFieldException x = read(file, "x15")

x = read(file, "x16")
@test !isa(x, TestType14)
@test x.x === 1.2345
@test x.y === 67

x = read(file, "x17")
@test !isa(x, TestType15)
@test x.x === 8.91011
@test x.y === 12

x = read(file, "x18")
@test !isa(x, TestType16)
@test x.x === 12.131415
@test x.y === 1617

x = read(file, "x19")
@test !isa(x, TestTypeContainer2)
@test x.a.x === 3
@test x.b === 4

x = read(file, "x20")
@test isa(x, TestTypeContainer3{Any,Int})
@test x.a.x === 3
@test x.b === 4

x = read(file, "x21")
@test !isa(x, TestTypeContainer4)
@test x.a.x === 3
@test x.b === 4

x = read(file, "x22")
@test !isa(x, TestTypeContainer5)
@test x.a.x === 3
@test x.b === 4

x = read(file, "x23")
@test x.x === 5678910

x = read(file, "x24")
@test first(x).first.x === 1
@test first(x).second.x.x === 2

x = read(file, "x25")
@test first(x).first.x === 1
@test reinterpret(UInt16, first(x).second) === 0x4567

x = read(file, "x26")
@test x.x === 1337

x = read(file, "x27")
@test x.x === 31337

x = read(file, "x28")
@test !isa(x, TestType20)
@test reinterpret(UInt16, x) == 0x1337

x = read(file, "x29")
@test sizeof(x) == 0

close(file)
end


text:
## path to TRPL folder --> set it before running the code!
path_to_trpl = "D:\\Projects\\TRPL" # some path to the TRPL folder (just an example here)

## activate package environment
cd(@__DIR__)
using Pkg
Pkg.activate("Project.toml")

## loading Packages
using OrdinaryDiffEq, DiffEqSensitivity, ForwardDiff, LinearAlgebra,
      Random, Printf, Distributions, BlackBoxOptim, NLopt, DataFrames, CSV

## input parameters
# Parameters here are setup for a quick test!
# For a real optimization set n_iter >= 2500000, max_time_loc_opt = 3600 to be save!
job_dir = string(path_to_trpl,"\\Optimizations\\Test\\Job2") # job directory (just for tests)
name_syn = "syn_data_1" # name of synthetic data set
n_iter = 2500 # number of optimization iteration and runs
n_runs =  1 # number of optimization runs
ftol_loc_opt = 1e-14 # tolerance for local optimization
max_time_loc_opt = 15 # maximum time in s for local optimization
#          τ_r0,  τ_nr0, τ_2t, τ_1t, τ_1e, τ_1d,    α,    β,    r,   C
lb = vcat([   1,      1,    1,    1,    1,    1,    1,    1,    1,   1],
           repeat([1e-6],10)) # lower parameter bounds
#          τ_r0,  τ_nr0, τ_2t, τ_1t, τ_1e, τ_1d,    α,    β,    r,   C
ub = vcat([5000,   5000, 1000, 1000,  1e9, 1000,  1e5,  1e5, 1000, 1e6],
           repeat([ 1.0],10)) # upper parameter bounds

## data directory
cd(joinpath(path_to_trpl,"Data\\Syn\\CSV"))

## load measurement data
data_syn = Matrix(CSV.read(name_syn*".csv",DataFrame,skipto=2))

## job directory
cd(job_dir)

## function for data processing (data structre is explained in README.md)
function data_processing(data)
    t = data[4:end,1]
    num_traj = size(data,2)-1
    num_pars = 10
    dims = num_pars + num_traj
    C_scaling = data[1,2]
    ND = 10 .^ data[2,2:end]
    background = data[3,2:end]
    IPL = data[4:end,2:num_traj+1]
    tspan = (t[1], t[end])
    N = length(t)
    idx_IPL = zeros(Bool,size(IPL))
    for i = 1:num_traj
	idx_IPL[:,i] = IPL[:,i] .> 0
    end
    return t, tspan, num_traj, num_pars, dims, IPL, idx_IPL, ND, background, N, C_scaling
end

## process synthetic data
t, tspan, num_traj, num_pars, dims, IPL, idx_IPL, ND, background, N, C_scaling = data_processing(data_syn)

## lower and upper bound dimension
lb = lb[1:dims]
ub = ub[1:dims]

## definition rate equations
function rate_equations(du, u, p, t)
    du[1] = -((1 - u[2]) / p[3]) * u[1] - ((1 - u[3]) / p[4]) * u[1] + (1 / p[5]) * u[3] * 1 / p[7] -
            ((1 + u[1]) / p[1]) * u[1] - ((1 + u[1]) / (p[2] * (1 + p[9] * u[1]))) * u[1]
    du[2] = ((1 - u[2]) / (p[3])) * p[8] * u[1]
    du[3] = ((1 - u[3]) / p[4]) * p[7] * u[1] - u[3] / p[5] - u[3] / p[6]
end

## ODEProblem definition with temporary u0 and parameters
u0 = [0.1, 0, 0]
prob = ODEProblem(rate_equations, u0, tspan, lb)

## definition loss function
function loss(par)
    l = 0.0
    function prob_func(prob, i, repeat)
        remake(prob, u0 = [par[num_pars+i], 0, 0])
    end
    ensemble_simu = EnsembleProblem(prob, prob_func = prob_func)
    sol = Array(solve(ensemble_simu, Rosenbrock23(), p = par[1:num_pars-1], EnsembleThreads(),
		              trajectories = num_traj, saveat = t,save_idxs = [1],
		              reltol=1e-6,abstol=1e-8))[1,:,:]
    for i = 1:num_traj
        model = ((C_scaling * par[num_pars] / ND[i]) .* sol[:,i] .* (sol[:,i] .+ 1)) .+ background[i]
        idx_used = idx_IPL[:,i] .& (model .> 0) # we ignore zeros for the log
        l += sum((model[idx_used]) .- (IPL[idx_used,i] .* log.(model[idx_used])))
    end
    return l
end

## definition nlopt loss function
function loss_nlopt(par::Vector,grad::Vector)
    if length(grad) > 0
	    grad[:] = ForwardDiff.gradient(loss,par)
    end
    return loss(par)
end

## multi run optimization
function multirun_opt(n_iter,n_runs)
    results = Array{Any}(undef,(n_runs+2,dims+3))
    results[1:n_runs,1] = ["Run $i" for i = 1:n_runs]
    results[end-1:end,[1,2,dims+3]] = ["lb" "x" "x";"UP" "x" "x"]
    results[end-1:end,3:end-1] = hcat(lb,ub)'
    for i = 1:n_runs
        seed_num = abs(rand(Int))
        println("$(seed_num) RNG Seed")
        Random.seed!(seed_num)
        println("Run $(i)")
        flush(stdout)
        optim_setup = bbsetup(loss,Method = :de_rand_1_bin_radiuslimited,
                              SearchRange = collect(zip(lb, ub)),
                              NumDimensions = dims, MaxFuncEvals = n_iter,
                              TraceMode = :compact, TraceInterval = 30)
        optim_result = bboptimize(optim_setup)
        p_optim = best_candidate(optim_result)
        optl = Opt(:LN_NEWUOA_BOUND,dims)
        optl.min_objective = loss_nlopt
        optl.lower_bounds = lb
        optl.upper_bounds = ub
        optl.ftol_abs = ftol_loc_opt
        optl.maxtime = max_time_loc_opt
        (optfl,optpl,retl) = optimize(optl,p_optim)
        println(optfl)
        println("")
        println(optpl)
        println("")
        results[i,2:end] = vcat(optfl,optpl,Int(seed_num))
    end
    return results
end
## multi run optimization
res_opt = multirun_opt(n_iter,n_runs)

## output preparation optimization results
df_opt = DataFrame(res_opt,:auto)
var_names_opt = vcat(["Run/Bounds","F_min", "τ_r0", "τ_nr0", "τ_2t",
		              "τ_1t", "τ_1e", "τ_1d", "α", "β", "r","C"],
		              [repeat("η_0$i", 1) for i = 1:num_traj],"Seed")
rename!(df_opt, var_names_opt)

## output results
CSV.write(name_syn*"_opt_res.csv", df_opt)

## message
println("Script has finished")

max_stars_repo_name: MaxGrossmann/TRPL

max_issues_repo_name: MaxGrossmann/TRPL

max_forks_repo_name: MaxGrossmann/TRPL


text:
@testset "ParametersEstimation.jl" begin
    when = [1,2,3]
    times = [1,3]
    marks = [2, 3]

    ts = TimeSeries(times)
    mts = MarkedTimeSeries(times, marks)

    μ = 2.0
    ϕ = 1
    γ = 1
    δ = 1.0

    @testset "volfunc" begin
        @test_throws ErrorException SEMPP.volfunc(when, ts, -γ, δ)
        @test_logs (:warn,"no marks but δ non zero") SEMPP.volfunc(when, ts, γ, δ)
        @test size(SEMPP.volfunc(when, mts, γ, δ)) == size(when)
        @test SEMPP.volfunc(when, mts, γ, δ) == [0, 3*exp(-1), 3*exp(-2)]
    end
end

include("ParametersEstimation_tests/DiscreteProcessEstimation_tests.jl")
include("ParametersEstimation_tests/ContinuousProcessEstimation_tests.jl")


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    groupby(data, var)

Partition spatial `data` into groups of constant value
for spatial variable `var`.

### Notes

Missing values are grouped into a separate group.
"""
function groupby(data::Data, var::Symbol)
  vars = Tables.schema(values(data)).names
  @assert var ∈ vars "invalid variable name"

  # data for variable
  vdata = data[var]

  # partition function with missings
  function f(i, j)
    vi, vj = vdata[i], vdata[j]
    mi, mj = ismissing(vi), ismissing(vj)
    (mi && mj) || ((!mi && !mj) && (vi == vj))
  end

  # partition function without missings
  g(i, j) = vdata[i] == vdata[j]

  # select the appropriate predicate function
  pred = Missing <: eltype(vdata) ? f : g

  # perform partition
  p = partition(data, PredicatePartition(pred))

  # retrieve value from each subset
  vals = [vdata[s[1]] for s in subsets(p)]

  # save metadata
  metadata = Dict(:values => vals)

  Partition(data, subsets(p), metadata)
end


text:
module Drums
using MusicManipulations
using MIDI: NOTEON, NOTEOFF
# all pitches of digitally connected instruments that support extra velocities
const DIGITAL = [0x26,0x28,0x33,0x35,0x3b]

"""
    getnotes_td50(track::MIDITrack, tpq = 960)

Get notes from midi track. Take care of Roland TD-50's ability to have velocities up to 159 for snare and ride.
"""
function getnotes_td50(track::MIDI.MIDITrack, tpq = 960)
    notes = MoreVelNote[]
    tracktime = UInt(0)
    extravel = 0
    for (i, event) in enumerate(track.events)
        tracktime += event.dT
        # Read through events until a noteon with velocity higher tha 0 is found
        if isa(event, MIDIEvent) && event.status & 0xF0 == NOTEON && event.data[2] > 0
            duration = UInt(0)
            #Test if the next event is an extra velocity event and modify velocity if needed.
            if event.data[1] in DIGITAL && event.data[2]==0x7f && track.events[i+1].status==0xb0 && track.events[i+1].data[1]==0x58
                extravel = floor(UInt8,track.events[i+1].data[2]/2)
                if extravel > 32
                    extravel = 32
                end
            end
            #Test if the previous event is an extra velocity event and modify velocity if needed.
            if i>2 #first event is alwas METAEvent
                if event.data[1] in DIGITAL && event.data[2]==0x7f && track.events[i-1].status==0xb0 && track.events[i-1].data[1]==0x58
                    extravel = floor(UInt8,track.events[i-1].data[2]/2)
                    if extravel > 32
                        extravel = 32
                    end
                end
            end
            for event2 in track.events[i+1:length(track.events)]
                duration += event2.dT
                # If we have a MIDI event & it's a noteoff (or a note on with 0 velocity), and it's for the same note as the first event we found, make a note
                # Many MIDI files will encode note offs as note ons with velocity zero
                if isa(event2, MIDI.MIDIEvent) && (event2.status & 0xF0 == MIDI.NOTEOFF || (event2.status & 0xF0 == MIDI.NOTEON && event2.data[2] == 0)) && event.data[1] == event2.data[1]
                    push!(notes, MoreVelNote(event.data[1], event.data[2]+extravel, tracktime, duration, event.status & 0x0F))
                    break
                end
            end
            extravel = 0
        end
    end
    sort!(notes, lt=((x, y)->x.position<y.position))
    return Notes(notes, tpq)
end

"""
    rm_hihatfake(notes::MIDI.Notes, BACK = 100, FORW = 100, CUTOFF = 0x16)

Remove fake tip notes generated by Hihat at certain actions by spotting them and
just removing every Hihat Head event with velocity less than `CUTOFF`and position
maximum `BACK` ticks before or `FORW` ticks after foot close.
"""
function rm_hihatfake!(notes::MIDI.Notes, BACK = 100, FORW = 100, CUTOFF = 0x16)

    #first map special closed notes
    for note in notes
        if note.pitch == 0x16
            note.pitch = 0x1a
        elseif note.pitch == 0x2a
            note.pitch = 0x2e
        end
    end
    #then look for fake notes
   i = 1
   deleted = 0
   len = length(notes)
   while i <= len
      #find foot close
      if notes[i].pitch == 0x2c || notes[i].pitch == 0x1a
         #go back and remove all fake tip strokes
         j = i-1
         #search all notes in specified BACK region
         while j>0 && notes[i].position-notes[j].position < BACK
            #if they are quiet enough
            if notes[j].pitch == 0x2e && notes[j].velocity <= CUTOFF
               #remove them
               deleteat!(notes.notes,j)
               deleted += 1
               i-=1
               len-=1
            else
               j-=1
            end
         end
         #go forward and remove all fake tip strokes
         j=i+1
         #search all notes in specified FORW region
         while j<=len && notes[j].position-notes[i].position < FORW
            #if they are quiet enough
            if notes[j].pitch == 0x2e && notes[j].velocity <= CUTOFF
               #remove them
               deleteat!(notes.notes,j)
               deleted += 1
               len-=1
            else
               j+=1
            end
         end
      end
      i+=1
   end
   println("deleted $(deleted) fake notes")
end

# Map the pitches of midinotes to Instruments of Roland TD-50
const MAP_TD50 = Dict{UInt8,String}(
    0x16=>"Hihat Rim (closed)",
    0x1a=>"Hihat Rim",
    0x24=>"Kick",
    0x25=>"Snare RimClick",
    0x26=>"Snare",
    0x27=>"Tom 4 Rimshot",
    0x28=>"Snare Rimshot",
    0x29=>"Tom 4",
    0x2a=>"Hihat Head (closed)",
    0x2b=>"Tom 3",
    0x2c=>"Hihat Foot Close",
    0x2d=>"Tom 2",
    0x2e=>"Hihat Head",
    0x2f=>"Tom 2 Rimshot",
    0x30=>"Tom 1",
    0x31=>"Cymbal 1",
    0x32=>"Tom 1 Rimshot",
    0x33=>"Ride Head",
    0x34=>"Cymbal 2",
    0x35=>"Ride Bell",
    0x37=>"Cymbal 1",
    0x39=>"Cymbal 2",
    0x3a=>"Tom 3 Rimshot",
    0x3b=>"Ride Rim")

# All posible pitches in an Array
const ALLPITCHES_TD50 = collect(keys(MAP_TD50))

# Map the pitches to numbers for plotting in a graph
const REORDER_TD50 = Dict{UInt8,UInt8}(
    0x16=>8,
    0x1a=>5,
    0x24=>0,
    0x25=>3,
    0x26=>1,
    0x27=>19,
    0x28=>2,
    0x29=>18,
    0x2a=>7,
    0x2b=>16,
    0x2c=>6,
    0x2d=>14,
    0x2e=>4,
    0x2f=>15,
    0x30=>12,
    0x31=>20,
    0x32=>13,
    0x33=>9,
    0x34=>21,
    0x35=>11,
    0x37=>20,
    0x39=>21,
    0x3a=>17,
    0x3b=>10)

###############################################################################
#velocity quantization
###############################################################################

"""
    td50_velquant_interval(notes::MIDI.Notes, numintervals::Int)

Divide the velocity range in `numintervals` intervals and quantize the
velocities of each `Note` to the mean value of all notes of the corresponding
instrument in this interval.
"""
function td50_velquant_interval(notes::MIDI.Notes{N}, numintervals::Int) where {N}
    #get notes separated by pitches
    sep = separatepitches(notes)
    newnotes = MIDI.Notes{N}(Vector{N}(), notes.tpq)

    for pitch in keys(sep)
        #short acces to needed notes
        pitchnotes = sep[pitch].notes

        #take care of different maximum velocities
        maxvel = 0
        pitch in DIGITAL ? maxvel = 160 : maxvel = 128

        #do a histogram and weight it with the velocities
        hist = zeros(maxvel)
        for note in pitchnotes
            hist[note.velocity] += 1
        end
        whist = copy(hist)
        for i = 1:length(hist)
            whist[i] *= i
        end

        #create the partitioning and compute corresponding means
        intlength = ceil(Int, maxvel/numintervals)
        meanvals = zeros(Int, numintervals)
        for i = 0:numintervals-1
            start = i*intlength+1
            ende = i*intlength+intlength
            if ende > maxvel
                ende = maxvel
            end
            piece = whist[start:ende]
            hits = sum(hist[start:ende])
            if hits != 0
                piece ./= hits
            end
            meanvals[i+1] = round(Int,mean(piece)*intlength)
        end

        #quantize notes
        for note in pitchnotes
            quant = ceil(Int, note.velocity/intlength)
            note.velocity = meanvals[quant]
        end

        # append to field of quantized notes
        append!(newnotes.notes, pitchnotes)
    end

    #restore temporal order
    sort!(newnotes.notes, lt=((x, y)->x.position<y.position))
    return newnotes
end

end


text:
function (eos::PressureEquation{<:Murnaghan1st})(v)
    @unpack v0, b0, b′0 = getparam(eos)
    return b0 / b′0 * ((v0 / v)^b′0 - 1)
end
function (eos::PressureEquation{<:Murnaghan2nd})(v)
    @unpack v0, b0, b′0, b″0 = getparam(eos)
    h = b′0^2 - 2b0 * b″0
    r = (v0 / v)^h
    return 2b0 / b′0 / (h / b′0 * ((r + 1) / (r - 1)) - 1)
end
function (eos::PressureEquation{<:BirchMurnaghan2nd})(v)
    @unpack v0, b0 = getparam(eos)
    f = ToEulerianStrain(v0)(v)
    return 3b0 * f * (2f + 1)^_2½
end
function (eos::PressureEquation{<:BirchMurnaghan3rd})(v)
    @unpack v0, b0, b′0 = getparam(eos)
    f = ToEulerianStrain(v0)(v)
    return 3f / 2 * b0 * (2f + 1)^_2½ * (2 + 3f * (b′0 - 4))
end
function (eos::PressureEquation{<:BirchMurnaghan4th})(v)
    @unpack v0, b0, b′0, b″0 = getparam(eos)
    f = ToEulerianStrain(v0)(v)
    h = b″0 * b0 + b′0^2
    return b0 / 2 * (2f + 1)^_2½ * ((9h - 63b′0 + 143) * f^2 + 9f * (b′0 - 4) + 6)
end
function (eos::PressureEquation{<:PoirierTarantola2nd})(v)
    @unpack v0, b0 = getparam(eos)
    f = ToNaturalStrain(v0)(v)
    return -3b0 * f * exp(-3f)
end
function (eos::PressureEquation{<:PoirierTarantola3rd})(v)
    @unpack v0, b0, b′0 = getparam(eos)
    f = ToNaturalStrain(v0)(v)
    return -3b0 / 2 * f * exp(-3f) * (3f * (b′0 - 2) + 1)
end
function (eos::PressureEquation{<:PoirierTarantola4th})(v)
    @unpack v0, b0, b′0, b″0 = getparam(eos)
    f = ToNaturalStrain(v0)(v)
    h = b″0 * b0 + b′0^2
    return -3b0 / 2 * f * exp(-3f) * (3f^2 * (h + 3b′0 + 3) + 3f * (b′0 - 2) + 2)
end
function (eos::PressureEquation{<:Vinet})(v)
    @unpack v0, b0, b′0 = getparam(eos)
    x, y = (v / v0)^_⅓, 3 / 2 * (b′0 - 1)
    return 3b0 / x^2 * (1 - x) * exp(y * (1 - x))
end
function (eos::PressureEquation{<:AntonSchmidt})(v)
    @unpack v0, b0, b′0 = getparam(eos)
    x, n = v / v0, -b′0 / 2
    return -b0 * x^n * log(x)
end
function (eos::PressureEquation{<:Holzapfel{Z}})(v) where {Z}
    @unpack v0, b0, b′0 = getparam(eos)
    η = (v / v0)^_⅓
    p0 = FERMI_GAS_CONSTANT * (Z / v0)^(5 / 3)
    c0 = -log(3b0 / p0 |> NoUnits)
    c2 = 3 / 2 * (b′0 - 3) - c0
    return 3b0 * (1 - η) / η^5 * exp(c0 * (1 - η)) * (1 + c2 * η * (1 - η))
end


text:
using HttpServer
using AccessControl
using EscapeString    # For escaping the user's name before displaying to the user
using LoggedDicts     # Data store for sessions as well as access control data (users, login credentials, permissions)

# Server-side sessions stored in a LoggedDict
sessions = LoggedDict("sessions", "sessions.log", true)

# Configure access control
function member_redirect(username::AbstractString)              # Function that defines a redirect location
    role = AccessControl.get_role(username)
    "/members/$(role)"
end
acdata          = LoggedDict("acdata", "acdata.log")            # Initialize access control data
session_config  = Dict(:datastore => sessions)
login_config    = Dict(:success_redirect => member_redirect)    # Use member_redirect(username) after login
logout_config   = Dict(:redirect => "/home")
pwdreset_config = Dict(:success_redirect => member_redirect)
AccessControl.update_config!(acdata, session = session_config, login = login_config, logout = logout_config, pwdreset = pwdreset_config)

# Add users to acdata.
# INSECURE!!! Because usernames and passwords are visible in plain text. See Admin Access for the secure way to do this.
AccessControl.create_user!("Alice", "pwd_alice", Set(["gold"]))
AccessControl.create_user!("Bob",   "pwd_bob",   Set(["silver"]))

# Handlers
function home!(req, res)
    res.data = "This is the home page. Anyone can visit here.

                $(login_form())"
end

function members_gold!(req, res)
    !has_role(req, "gold") && (notfound!(res); return)        # Check whether user is logged in and has the "gold" role
    session_id = read_sessionid(req)
    username   = session_get(session_id, "username")
    username   = escapestring(username, :html_text)
    res.data   = "Welcome $(username)! This page displays information for GOLD members only.

                  $(logout_pwdreset_links())"
end

function members_silver!(req, res)
    !has_role(req, "silver") && (notfound!(res); return)        # Check whether user is logged in and has the "silver" role
    session_id = read_sessionid(req)
    username   = session_get(session_id, "username")
    username   = escapestring(username, :html_text)
    res.data   = "Welcome $(username)! This page displays information for SILVER members only.

                  $(logout_pwdreset_links())"
end

# App
function app(req::Request)
    res  = Response()
    rsrc = req.resource
    if rsrc == "/home"
        home!(req, res)
    elseif rsrc == "/members/gold"
        members_gold!(req, res)
    elseif rsrc == "/members/silver"
        members_silver!(req, res)
    elseif haskey(acpaths, rsrc)        # Note this new clause in the if statement for handling access control
        acpaths[rsrc](req, res)
    else
        notfound!(res)
    end
    res
end

# Run the app under HTTPS rather than HTTP
rel(filename::AbstractString, p::AbstractString) = joinpath(dirname(filename), p)
if !isfile("keys/server.crt")
    @unix_only begin
	run(`mkdir -p $(rel(@__FILE__, "keys"))`)
	run(`openssl req -x509 -nodes -days 365 -newkey rsa:2048 -keyout $(rel(@__FILE__, "keys/server.key")) -out $(rel(@__FILE__, "keys/server.crt"))`)
    end
end

server = Server((req, res) -> app(req))
cert   = MbedTLS.crt_parse_file(rel(@__FILE__, "keys/server.crt"))
key    = MbedTLS.parse_keyfile(rel(@__FILE__, "keys/server.key"))
@async run(server, port = 8000, ssl = (cert, key))


text:
# Intersection methods
# TODO: Change intersection methods -- This breaks with infinite slope
function intersect(l1::LineSegment, l2::LineSegment)
    # Check for all cases in which they don't intersect
    # Will check for: 1) Whether l1.p1 and l1.p2 are on opposite sides of
    # l2.p1->l2.p2 2) Whether l2.p1 and l2.p2 are on opposite sides of
    # l1.p1->l2.p2 3) Whether they have same slope
    # See: http://jeffe.cs.illinois.edu/teaching/373/notes/x06-sweepline.pdf
    p1, p2 = l1.p1, l1.p2
    p3, p4 = l2.p1, l2.p2
    if ccw(p1, p3, p4) == ccw(p2, p3, p4)
        return (false, Point(Inf, Inf))
    elseif ccw(p3, p1, p2) == ccw(p4, p1, p2)
        return (false, Point(Inf, Inf))
    elseif l1.slope == l2.slope
        return (false, Point(Inf, Inf))
    end

    # If we make it this far, then line segments should intersect
    # Pull out points
    x1 = l1.p1[1]
    y1 = l1.p1[2]
    x2 = l1.p2[1]
    y2 = l1.p2[2]
    x3 = l2.p1[1]
    y3 = l2.p1[2]
    x4 = l2.p2[1]
    y4 = l2.p2[2]

    # Get intersection
    denom = (x1-x2)*(y3-y4) - (y1-y2)*(x3-x4)
    x_intersect = (x1*y2 - y1*x2)*(x3-x4) - (x1-x2)*(x3*y4 - y3*x4)/denom
    y_intersect = (x1*y2 - y1*x2)*(y3-y4) - (y1-y2)*(x3*y4 - y3*x4)/denom
    intersection = Point(x_intersect, y_intersect)

    # Check if intersection is on both intervals
    if in(x_intersect, Interval(x1, x2)) && in(x_intersect, Interval(x3, x4))
        out = (true, intersection)
    # If not, return intersection but negative status
    else
        out = (false, intersection)
    end
in
    return out
end

function intersect(ch1::ConvexHull, ch2::ConvexHull)
    error("Intersection between convex sets is not currently implemented")
end

function intersect(ch::ConvexHull, q::Quadrant)

    # Get information about the set of extreme points
    ext_pts = ch.extremepoints
    xmin, _, xmax, _ = extremax(ext_pts)
    ymin, _, ymax, _ = extremay(ext_pts)
    npts = length(ch)

    # Get information about quadrant
    o = q.origin
    lx = LineSegment(Point(clamp(q.lx.p1[1], xmin, xmax), clamp(q.lx.p1[2], ymin, ymax)),
                     Point(clamp(q.lx.p2[1], xmin, xmax), clamp(q.lx.p2[2], ymin, ymax)))
    ly = LineSegment(Point(clamp(q.ly.p1[1], xmin, xmax), clamp(lx.p1[2], ymin, ymax)),
                     Point(clamp(q.ly.p2[1], xmin, xmax), clamp(lx.p2[2], ymin, ymax)))

    # Start a convex hull to store outpoints in
    prev_pt = ext_pts[end]
    outhull = in(prev_pt, q) ? [prev_pt] : []

    # Will check between each point individually
    # Points need to be ordered either clockwise
    # or counter-clockwise
    for i=1:npts
        # Current extreme point and corresponding line
        c_ep = ext_pts[i]
        c_ls = LineSegment(prev_pt, c_ep)

        # Compute intersections with lines
        stat_x, c_inter_x = intersect(lx, c_ls)
        stat_y, c_inter_y = intersect(ly, c_ls)

        # If there is intersection, push to hull
        stat_x ? push!(outhull, c_inter_x) : nothing
        stat_y ? push!(outhull, c_inter_y) : nothing

        # Add current point if it is in quadrant
        in(c_ep, q) ? push!(outhull, c_ep) : nothing
        prev_pt = c_ep

    end

    return outhull
end

# function intersect(ch::ConvexHull, q::Quadrant)
#     # TODO: Need to consider special cases where it is not binding on one side
#     # but is on the other. It should be as simple as returning all extrema to
#     # the right side of the binding extreme. Right thing to do here is break
#     # this function into cases and write the cases out, that will make it
#     # look much less messy

#     # Get information about the set of extreme points
#     ext_pts = ch.extremepoints
#     _, xmin, _, xmax = extremax(ext_pts)
#     _, ymin, _, ymax = extremay(ext_pts)
#     npts = length(ep)

#     # Get information about quadrant
#     o = q.origin

#     # Counters to determine whether points belong in the
#     # set or not. vcount and hcount are for number of intersections
#     # the line makes with the convex hull each way. If either are
#     # even then origin is not in intersection
#     checknext = true
#     storecurr = false
#     vcount = 0
#     hcount = 0
#     # Start a convex hull to store outpoints in
#     outhull = []

#     # Will check between each point individually
#     # Points need to be ordered either clockwise
#     # or counter-clockwise
#     for i=2:npts+1
#         # Pull out points
#         if i<=npts
#             x1 = ep[i-1]
#             x2 = ep[i]
#         else
#             x1 = ep[npts]
#             x2 = ep[1]
#         end

#         # First add point if we are supposed to store
#         if storecurr
#             push!(outhull, x1)
#         end

#         # Check whether intersects
#         w1intersect = evalliney(x1, x2, o.x)
#         w2intersect = evallinex(x1, x2, o.y)
#         cond1 = ((x2.x - o.x)*(x1.x - o.x) < 0) && (w1intersect >= o.y)
#         cond2 = ((x2.y - o.y)*(x1.y - o.y) < 0) && (w2intersect >= o.x)

#         # Check first condition
#         if cond1
#             storecurr = !storecurr
#             push!(outhull, Point(w1, w1intersect))
#             vcount += 1
#         # Account for if point is right on top of value
#         elseif isapprox(x2.x, o.x) && isapprox(x2.y, w1intersect)
#             if storecurr
#                 push!(outhull, x2)
#             end
#             storecurr = !storecurr
#             vcount += 1
#         end

#         # Check second condition
#         if cond2
#             storecurr = !storecurr
#             # Make sure we don't double count if both conditions hold
#             push!(outhull, Point(w2intersect, w2))
#             hcount += 1
#         elseif isapprox(x2.y, o.y) && isapprox(x2.x, w2intersect)
#             if storecurr
#                 push!(outhull, x2)
#             end
#             storecurr = !storecurr
#             hcount += 1
#         end
#     end

#     if (hcount*vcount)%2!=0 || (o.x<minx(ep)[1].x && o.y<miny(ep)[1].y)
#         push!(outhull, o)
#     end

#     return outhull
# end


text:
using LightGraphs, MetaGraphs
using LightGraphs: AbstractGraph, inneighbors
using MetaGraphs: AbstractMetaGraph

export adjlist, getvprops, geteprops, getvalue, getpairindices, cleangraph!
function adjlist(g::AbstractGraph)
    n = nv(g)
    adj = Vector{Int}[inneighbors(g, i) for i = 1:n]
    return adj
end

function getvprops(g::AbstractMetaGraph, p::Array{<:Pair,1})
    propsinds = getpairindices(p,[])
    nf = length(propsinds)
    n = nv(g)
    features = Array{Float64,2}(undef, nf, n)
    for (i,v) in enumerate(vertices(g))
        feature = [getvalue(keys, props(g,v)) for keys in propsinds]
        features[:,i] .= feature
    end
    return features
end

function geteprops(g::AbstractMetaGraph, p::Array{<:Pair,1}, rev::Bool=false)
    propsinds = getpairindices(p,[])
    nf = length(propsinds)
    nofv = nv(g)
    nofe = ne(g)
    features = fill(Float64(0.0), nf, nofv, nofv)
    for (i,e) in enumerate(edges(g))
        r,s = dst(e),src(e)
        feature = [getvalue(keys, props(g,Edge(s,r))) for keys in propsinds]
        rev && (s = dst(e); r = src(e))
        features[:,s,r] .= feature
    end
    return features
end

function getvalue(indices, dict)
    v = dict[indices[1]]
    if length(indices) > 1
        for i in indices[2:end]
            v = v[i]
        end
    end
    v
end

function getpairindices(p::Array{<:Pair,1}, op)
    out = []
    for i in p
        oc = vcat(op, i[1])
        if isa(i[2], Array{<:Pair,1})
            push!(out, getpairindices(i[2], oc)...)
        elseif isa(i[2], Nothing)
            push!(out, oc)
        end
    end
    return out
end

function cleangraph!(g)
    rnodes = []
    redges = []
    for e in edges(g)
        edgeinfo = props(g, e)
        maininfo = props(g, src(e))
        neiginfo = props(g, dst(e))
        (edgeinfo[:dist] == -1.0) && push!(redges, e)
        ((maininfo[:nodeev][:ucRsrpRslt] == -1) | (maininfo[:cellsinfo] == -1 )) && push!(rnodes, src(e))
        ((neiginfo[:nodeev][:ucRsrpRslt] == -1) | (neiginfo[:cellsinfo] == -1)) && push!(rnodes, dst(e))
    end
    for n in vertices(g)
        maininfo = props(g,n)
        ((maininfo[:nodeev][:ucRsrpRslt] == -1) | (maininfo[:cellsinfo] == -1 )) && push!(rnodes, n)
    end
    for e in union(redges)
        rem_edge!(g,e)
    end
    for n in union(rnodes)
        rem_vertex!(g, n)
    end
end


text:
const AbstractMatrixBlock{T, M} = SubArray{T,2,M,Tuple{UnitRange{Int},UnitRange{Int}},false}

"""
$(TYPEDEF)

`SegmentedBlockDiagonalMatrix` is an `AbstractMatrix` backed by a parent `AbstractMatrix`, which
additionally stores a sequence of views into the diagonal blocks of the parent matrix. This
type is useful for storing and updating block-diagonal matrices whose block contents
may change but whose overall structure is fixed, such as configuration derivative <-> velocity
jacobians.
"""
struct SegmentedBlockDiagonalMatrix{T, M<:AbstractMatrix{T}} <: AbstractMatrix{T}
    parent::M
    blocks::Vector{AbstractMatrixBlock{T, M}}

    function SegmentedBlockDiagonalMatrix{T, M}(parent::M, block_indices) where {T, M<:AbstractMatrix{T}}
        check_contiguous_block_ranges(parent, block_indices)
        blocks = collect(view(parent, indices...) for indices in block_indices)
        new{T, M}(parent, blocks)
    end
end

SegmentedBlockDiagonalMatrix(parent::M, block_indices) where {T, M<:AbstractMatrix{T}} = SegmentedBlockDiagonalMatrix{T, M}(parent, block_indices)

function SegmentedBlockDiagonalMatrix{T}(initializer, rows::Integer, cols::Integer, block_indices) where T
    parent = Matrix{T}(initializer, rows, cols)
    SegmentedBlockDiagonalMatrix{T}(parent, block_indices)
end

Base.parent(m::SegmentedBlockDiagonalMatrix) = m.parent
Base.size(m::SegmentedBlockDiagonalMatrix) = size(m.parent)
@propagate_inbounds Base.getindex(v::SegmentedBlockDiagonalMatrix, i::Int) = v.parent[i]
@propagate_inbounds Base.setindex!(v::SegmentedBlockDiagonalMatrix, value, i::Int) = v.parent[i] = value
Base.IndexStyle(::Type{<:SegmentedBlockDiagonalMatrix}) = IndexLinear()
blocks(m::SegmentedBlockDiagonalMatrix) = m.blocks

function check_contiguous_block_ranges(parent::AbstractMatrix, block_indices)
    if !_is_contiguous_and_diagonal(parent, block_indices)
        throw(ArgumentError("The `block_indices` should be a vector of index ranges corresponding to non-overlapping contiguous diagonal blocks"))
    end
end

function _is_contiguous_and_diagonal(parent::AbstractMatrix, block_indices)
    expected_starts = first.(axes(parent))
    for inds in block_indices
        if first.(inds) !== expected_starts
            return false
        end
        expected_starts = last.(inds) .+ 1
    end
    if expected_starts !== last.(axes(parent)) .+ 1
        return false
    end
    return true
end

function LinearAlgebra.mul!(C::Matrix, A::Matrix, B::SegmentedBlockDiagonalMatrix)
    # TODO: coordination with BLAS threads
    @boundscheck size(C) == (size(A, 1), size(B, 2)) || throw(DimensionMismatch("Output size mismatch."))
    A′ = Base.unalias(C, A)
    Threads.@threads for block in blocks(B) # allocates 32 bytes (see https://github.com/JuliaLang/julia/issues/29748)
        @inbounds begin
            Acols, Ccols = parentindices(block)
            Aview = uview(A′, :, Acols)
            Cview = uview(C, :, Ccols)
            if block == I
                copyto!(Cview, Aview)
            else
                mul!(Cview, Aview, block)
            end
        end
    end
    return C
end

function LinearAlgebra.mul!(C::Matrix, A::SegmentedBlockDiagonalMatrix, B::Matrix)
    # TODO: coordination with BLAS threads
    @boundscheck size(C) == (size(A, 1), size(B, 2)) || throw(DimensionMismatch("Output size mismatch."))
    B′ = Base.unalias(C, B)
    Threads.@threads for block in blocks(A) # allocates 32 bytes (see https://github.com/JuliaLang/julia/issues/29748)
        @inbounds begin
            Crows, Brows = parentindices(block)
            Bview = uview(B, Brows, :)
            Cview = uview(C, Crows, :)
            if block == I
                copyto!(Cview, Bview)
            else
                mul!(Cview, block, Bview)
            end
        end
    end
    return C
end


text:
export cursor, focused, frameCount, frameRate
export height, noCursor, noSmooth, smooth, width

# start with arrow cursor by default, as expected
currCursor = GLFW.CreateStandardCursor(GLFW.ARROW_CURSOR)

function cursor(window)
	GLFW.SetInputMode(window, GLFW.CURSOR, GLFW.CURSOR_NORMAL)
	GLFW.SetCursor(window, currCursor)
end

function cursor(window, cursorType)
	GLFW.SetInputMode(window, GLFW.CURSOR, GLFW.CURSOR_NORMAL)
	if cursorType == "HAND"
		currCursor = GLFW.CreateStandardCursor(GLFW.HAND_CURSOR)
	elseif cursorType == "ARROW"
		currCursor = GLFW.CreateStandardCursor(GLFW.ARROW_CURSOR)
	elseif cursorType == "CROSS"
		currCursor = GLFW.CreateStandardCursor(GLFW.CROSSHAIR_CURSOR)
	elseif cursorType == "MOVE"
		currCursor = GLFW.CreateStandardCursor(GLFW.IBEAM_CURSOR)
	elseif cursorType == "TEXT"
		currCursor = GLFW.CreateStandardCursor(GLFW.IBEAM_CURSOR)
	elseif cursorType == "WAIT"
		# currCursor = GLFW.CreateStandardCursor(GLFW.HAND_CURSOR)
	end
	GLFW.SetCursor(window, currCursor)
end

function focused(window)
    if GLFW.GetWindowAttrib(window, GLFW.FOCUSED)
        return true
    else
        return false
    end
end

function frameCount()
	return state.frameCount
end

function frameRate()
	return state.frameRate
end

function frameRate(fRate)
	state.frameRate = fRate
	return nothing
end

function height()
	return state.height
end

function noCursor(window)
	GLFW.SetInputMode(window, GLFW.CURSOR, GLFW.CURSOR_HIDDEN)
end

function noSmooth()
	glDisable(GL_MULTISAMPLE)
end

function smooth()
	glEnable(GL_MULTISAMPLE)
end

function width()
	return state.width
end


text:
@kwdef struct Page{OP, T}
    size::Int
    page::Int
    pagelen::Int
    next::Union{Nothing, String}
    values::Vector{T}
    _extras::NamedTuple
end
@json_struct Page

struct BBPaginator{FRG, P}
    api::FRG
    page::P
    BBPaginator(api::FRG, page::Page) where {FRG <: Forge} = new{FRG, P}(api, page)
end

paginate(api::BitbucketAPI, f::Function, args...) = paginate(api, f(api, args...))

paginate(api::BitbucketAPI, (page, _resp)::Tuple{Page, HTTP.Response}) = BBPaginator(api, page)

Base.iterate(p::BBPaginator{BitbucketAPI}) =
    isempty(p.page.values) ? nothing : (p.page.values[1], (p.page, 2))

function Base.iterate(pgn::BBPaginator{BitbucketAPI, P}, (p, offset)::Tuple{P, Integer}) where {OP, P <: Page{OP}}
    api = pgn.api
    offset <= length(p.values) && return p.page.values[offset], (p, offset + 1)
    p.next === nothing && return nothing
    rate_limit(api, OP)
    resp = try
        HTTP.request(
            :GET, p.next, request_headers(api, OP);
            status_exception=false,
        )
    catch e
        rethrow(HTTPError(e, stacktrace(catch_backtrace())))
    end
    has_rate_limits(api, OP) && rate_limit_update!(api, OP, resp)
    resp.status != 200 && throw(HTTPError(resp, HTTP.StatusError(resp.status, resp), stacktrace()))
    try
        p2 = postprocess(postprocessor(api, OP), resp, into(api, OP))
        p2.values[1], (p2, 2)
    catch e
        rethrow(PostProcessorError(resp, e, stacktrace(catch_backtrace())))
    end
end


text:
######################################
### BOILERPLATE FROM LECTURE NOTES ###
######################################

abstract type Event end
abstract type State end

# Captures an event and the time it takes place
struct TimedEvent
    event::Event
    time::Float64
    location_ID::Int
    arrival_time::Float64
end

# Comparison of two timed events - this will allow us to use them in a heap/priority-queue
isless(te1::TimedEvent, te2::TimedEvent) = te1.time < te2.time

"""
Object which represents state of overall system, e.g. total jobs orbiting,
total jobs within system.
"""

mutable struct NetworkState <: State
    number_in_system::Int # Total number of jobs within the network
    queues::Array{Int64} #Int64 used because scenario 4 has an infinite queue
    number_in_system_decreased::Bool
    orbiting_jobs::Int
    NetworkState() = new(0)
    NetworkState(x::Int, y::Array{Int8}, z::Int) = new(x,y, false, z)
end

struct ArrivalEvent <: Event end
struct EndOfServiceEvent <: Event end
struct OverflowEvent <: Event end

# Generic events that we can always use
"""
    EndSimEvent()

Return an event that ends the simulation.
"""
struct EndSimEvent <: Event end

"""
    LogStateEvent()

Return an event that prints a log of the current simulation state.
"""
struct LogStateEvent <: Event end

"""
The main simulation function gets an initial state and an initial event
that gets things going. Optional arguments are the maximal time for the
simulation, times for logging events, and a call-back function.
"""

function simulate(init_state::State, init_timed_event::TimedEvent
                    ;
                    max_time::Float64 = 10.0,
                    log_times::Vector{Float64} = Float64[],
                    callback = (time, state) -> nothing)

    # The event queue
    priority_queue = BinaryMinHeap{TimedEvent}()

    # Put the standard events in the queue
    push!(priority_queue, init_timed_event)
    # FIXME: TESTING NAN FOR THE FINAL VALUE OF TIMED EVENT
    push!(priority_queue, TimedEvent(EndSimEvent(), max_time, 0, NaN))
    for log_time in log_times
        push!(priority_queue, TimedEvent(LogStateEvent(), log_time, 0, NaN))
    end

    # initilize the state
    state = deepcopy(init_state)
    time = 0.0

    # Callback at simulation start
    callback(time, state, TimedEvent[])

    # The main discrete event simulation loop - SIMPLE!
    while true
        # Get the next event
        timed_event = pop!(priority_queue)

        # Advance the time
        time = timed_event.time

        # Act on the event
        new_timed_events = process_event(time, state, timed_event.event)

        # If the event was an end of simulation then stop
        if timed_event.event isa EndSimEvent
            break
        end

        # The event may spawn 0 or more events which we put in the priority queue
        for nte in new_timed_events
            push!(priority_queue,nte)
        end

        # Callback for each simulation event
        callback(time, state, timed_event)
    end
end;

"""
Overload Simulate with the Project 2 specific Simulate function, differs by taking in the premade NetworkParameters object as an arg.
"""

function simulate(init_state::State, init_timed_event::TimedEvent, scenario::NetworkParameters
    ;
    max_time::Float64 = 10.0,
    log_times::Vector{Float64} = Float64[],
    callback = (time, state) -> nothing)

    # The event queue
    priority_queue = BinaryMinHeap{TimedEvent}()

    # Put the standard events in the queue
    push!(priority_queue, init_timed_event)
    # FIXME: TESTING NAN FOR THE FINAL VALUE OF TIMED EVENT
    push!(priority_queue, TimedEvent(EndSimEvent(), max_time, 0, NaN))
    for log_time in log_times
        push!(priority_queue, TimedEvent(LogStateEvent(), log_time, 0, NaN))
    end

    # Initialize the network state
    state = deepcopy(init_state)
    time = 0.0

    # # initialize the queue states of the stations
    # state.queues = zeros(Int8, scenario.L) This is probably a bad approach to take when taking a "generic" system so Im commenting it out for the time being

    # Callback at simulation start
    callback(time, state, init_timed_event)

    # The main discrete event simulation loop
    while true
        # Get the next event
        timed_event = pop!(priority_queue)

        # Advance the time
        time = timed_event.time

        # Act on the event
        new_timed_events = process_event(time, state, timed_event.location_ID, timed_event.event, scenario, timed_event.arrival_time)

        # If the event was an end of simulation then stop
        if timed_event.event isa EndSimEvent
            break
        end

        # The event may spawn 0 or more events which we put in the priority queue
        for nte in new_timed_events
            push!(priority_queue,nte)
        end

        # Callback for each simulation event
        callback(time, state, timed_event)

        # Reset state changed variable
        state.number_in_system_decreased = false
    end
    return state
end


text:
using SciMLBenchmarks, Test

@testset "weave_file" begin
    benchmarks_dir = joinpath(dirname(@__DIR__), "benchmarks")
    SciMLBenchmarks.weave_file(joinpath(benchmarks_dir, "Testing"), "test.jmd")

    @test isfile(joinpath(dirname(@__DIR__), "script", "Testing", "test.jl"))
    @test isfile(joinpath(dirname(@__DIR__), "html", "Testing", "test.html"))
    @test isfile(joinpath(dirname(@__DIR__), "markdown", "Testing", "test.md"))
end


text:
# This file is auto-generated by AWSMetadata.jl
using AWS
using AWS.AWSServices: translate
using AWS.Compat
using AWS.UUIDs

"""
    CreateParallelData()

Creates a parallel data resource in Amazon Translate by importing an input file from Amazon S3. Parallel data files contain examples of source phrases and their translations from your translation memory. By adding parallel data, you can influence the style, tone, and word choice in your translation output.

# Required Parameters
- `ClientToken`: A unique identifier for the request. This token is automatically generated when you use Amazon Translate through an AWS SDK.
- `Name`: A custom name for the parallel data resource in Amazon Translate. You must assign a name that is unique in the account and region.
- `ParallelDataConfig`: Specifies the format and S3 location of the parallel data input file.

# Optional Parameters
- `Description`: A custom description for the parallel data resource in Amazon Translate.
- `EncryptionKey`:
"""
create_parallel_data(ClientToken, Name, ParallelDataConfig; aws_config::AWSConfig=global_aws_config()) = translate("CreateParallelData", Dict{String, Any}("ClientToken"=>ClientToken, "Name"=>Name, "ParallelDataConfig"=>ParallelDataConfig); aws_config=aws_config)
create_parallel_data(ClientToken, Name, ParallelDataConfig, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("CreateParallelData", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClientToken"=>ClientToken, "Name"=>Name, "ParallelDataConfig"=>ParallelDataConfig), args)); aws_config=aws_config)

"""
    DeleteParallelData()

Deletes a parallel data resource in Amazon Translate.

# Required Parameters
- `Name`: The name of the parallel data resource that is being deleted.

"""
delete_parallel_data(Name; aws_config::AWSConfig=global_aws_config()) = translate("DeleteParallelData", Dict{String, Any}("Name"=>Name); aws_config=aws_config)
delete_parallel_data(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("DeleteParallelData", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Name"=>Name), args)); aws_config=aws_config)

"""
    DeleteTerminology()

A synchronous action that deletes a custom terminology.

# Required Parameters
- `Name`: The name of the custom terminology being deleted.

"""
delete_terminology(Name; aws_config::AWSConfig=global_aws_config()) = translate("DeleteTerminology", Dict{String, Any}("Name"=>Name); aws_config=aws_config)
delete_terminology(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("DeleteTerminology", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Name"=>Name), args)); aws_config=aws_config)

"""
    DescribeTextTranslationJob()

Gets the properties associated with an asycnhronous batch translation job including name, ID, status, source and target languages, input/output S3 buckets, and so on.

# Required Parameters
- `JobId`: The identifier that Amazon Translate generated for the job. The StartTextTranslationJob operation returns this identifier in its response.

"""
describe_text_translation_job(JobId; aws_config::AWSConfig=global_aws_config()) = translate("DescribeTextTranslationJob", Dict{String, Any}("JobId"=>JobId); aws_config=aws_config)
describe_text_translation_job(JobId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("DescribeTextTranslationJob", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("JobId"=>JobId), args)); aws_config=aws_config)

"""
    GetParallelData()

Provides information about a parallel data resource.

# Required Parameters
- `Name`: The name of the parallel data resource that is being retrieved.

"""
get_parallel_data(Name; aws_config::AWSConfig=global_aws_config()) = translate("GetParallelData", Dict{String, Any}("Name"=>Name); aws_config=aws_config)
get_parallel_data(Name, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("GetParallelData", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Name"=>Name), args)); aws_config=aws_config)

"""
    GetTerminology()

Retrieves a custom terminology.

# Required Parameters
- `Name`: The name of the custom terminology being retrieved.
- `TerminologyDataFormat`: The data format of the custom terminology being retrieved, either CSV or TMX.

"""
get_terminology(Name, TerminologyDataFormat; aws_config::AWSConfig=global_aws_config()) = translate("GetTerminology", Dict{String, Any}("Name"=>Name, "TerminologyDataFormat"=>TerminologyDataFormat); aws_config=aws_config)
get_terminology(Name, TerminologyDataFormat, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("GetTerminology", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("Name"=>Name, "TerminologyDataFormat"=>TerminologyDataFormat), args)); aws_config=aws_config)

"""
    ImportTerminology()

Creates or updates a custom terminology, depending on whether or not one already exists for the given terminology name. Importing a terminology with the same name as an existing one will merge the terminologies based on the chosen merge strategy. Currently, the only supported merge strategy is OVERWRITE, and so the imported terminology will overwrite an existing terminology of the same name. If you import a terminology that overwrites an existing one, the new terminology take up to 10 minutes to fully propagate and be available for use in a translation due to cache policies with the DataPlane service that performs the translations.

# Required Parameters
- `MergeStrategy`: The merge strategy of the custom terminology being imported. Currently, only the OVERWRITE merge strategy is supported. In this case, the imported terminology will overwrite an existing terminology of the same name.
- `Name`: The name of the custom terminology being imported.
- `TerminologyData`: The terminology data for the custom terminology being imported.

# Optional Parameters
- `Description`: The description of the custom terminology being imported.
- `EncryptionKey`: The encryption key for the custom terminology being imported.
"""
import_terminology(MergeStrategy, Name, TerminologyData; aws_config::AWSConfig=global_aws_config()) = translate("ImportTerminology", Dict{String, Any}("MergeStrategy"=>MergeStrategy, "Name"=>Name, "TerminologyData"=>TerminologyData); aws_config=aws_config)
import_terminology(MergeStrategy, Name, TerminologyData, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("ImportTerminology", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("MergeStrategy"=>MergeStrategy, "Name"=>Name, "TerminologyData"=>TerminologyData), args)); aws_config=aws_config)

"""
    ListParallelData()

Provides a list of your parallel data resources in Amazon Translate.

# Optional Parameters
- `MaxResults`: The maximum number of parallel data resources returned for each request.
- `NextToken`: A string that specifies the next page of results to return in a paginated response.
"""
list_parallel_data(; aws_config::AWSConfig=global_aws_config()) = translate("ListParallelData"; aws_config=aws_config)
list_parallel_data(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("ListParallelData", args; aws_config=aws_config)

"""
    ListTerminologies()

Provides a list of custom terminologies associated with your account.

# Optional Parameters
- `MaxResults`: The maximum number of custom terminologies returned per list request.
- `NextToken`: If the result of the request to ListTerminologies was truncated, include the NextToken to fetch the next group of custom terminologies.
"""
list_terminologies(; aws_config::AWSConfig=global_aws_config()) = translate("ListTerminologies"; aws_config=aws_config)
list_terminologies(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("ListTerminologies", args; aws_config=aws_config)

"""
    ListTextTranslationJobs()

Gets a list of the batch translation jobs that you have submitted.

# Optional Parameters
- `Filter`: The parameters that specify which batch translation jobs to retrieve. Filters include job name, job status, and submission time. You can only set one filter at a time.
- `MaxResults`: The maximum number of results to return in each page. The default value is 100.
- `NextToken`: The token to request the next page of results.
"""
list_text_translation_jobs(; aws_config::AWSConfig=global_aws_config()) = translate("ListTextTranslationJobs"; aws_config=aws_config)
list_text_translation_jobs(args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("ListTextTranslationJobs", args; aws_config=aws_config)

"""
    StartTextTranslationJob()

Starts an asynchronous batch translation job. Batch translation jobs can be used to translate large volumes of text across multiple documents at once. For more information, see async. Batch translation jobs can be described with the DescribeTextTranslationJob operation, listed with the ListTextTranslationJobs operation, and stopped with the StopTextTranslationJob operation.  Amazon Translate does not support batch translation of multiple source languages at once.

# Required Parameters
- `ClientToken`: A unique identifier for the request. This token is auto-generated when using the Amazon Translate SDK.
- `DataAccessRoleArn`: The Amazon Resource Name (ARN) of an AWS Identity Access and Management (IAM) role that grants Amazon Translate read access to your input data. For more nformation, see identity-and-access-management.
- `InputDataConfig`: Specifies the format and S3 location of the input documents for the translation job.
- `OutputDataConfig`: Specifies the S3 folder to which your job output will be saved.
- `SourceLanguageCode`: The language code of the input language. For a list of language codes, see what-is-languages. Amazon Translate does not automatically detect a source language during batch translation jobs.
- `TargetLanguageCodes`: The language code of the output language.

# Optional Parameters
- `JobName`: The name of the batch translation job to be performed.
- `ParallelDataNames`: The names of the parallel data resources to use in the batch translation job. For a list of available parallel data resources, use the ListParallelData operation.
- `TerminologyNames`: The name of the terminology to use in the batch translation job. For a list of available terminologies, use the ListTerminologies operation.
"""
start_text_translation_job(ClientToken, DataAccessRoleArn, InputDataConfig, OutputDataConfig, SourceLanguageCode, TargetLanguageCodes; aws_config::AWSConfig=global_aws_config()) = translate("StartTextTranslationJob", Dict{String, Any}("ClientToken"=>ClientToken, "DataAccessRoleArn"=>DataAccessRoleArn, "InputDataConfig"=>InputDataConfig, "OutputDataConfig"=>OutputDataConfig, "SourceLanguageCode"=>SourceLanguageCode, "TargetLanguageCodes"=>TargetLanguageCodes); aws_config=aws_config)
start_text_translation_job(ClientToken, DataAccessRoleArn, InputDataConfig, OutputDataConfig, SourceLanguageCode, TargetLanguageCodes, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("StartTextTranslationJob", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClientToken"=>ClientToken, "DataAccessRoleArn"=>DataAccessRoleArn, "InputDataConfig"=>InputDataConfig, "OutputDataConfig"=>OutputDataConfig, "SourceLanguageCode"=>SourceLanguageCode, "TargetLanguageCodes"=>TargetLanguageCodes), args)); aws_config=aws_config)

"""
    StopTextTranslationJob()

Stops an asynchronous batch translation job that is in progress. If the job's state is IN_PROGRESS, the job will be marked for termination and put into the STOP_REQUESTED state. If the job completes before it can be stopped, it is put into the COMPLETED state. Otherwise, the job is put into the STOPPED state. Asynchronous batch translation jobs are started with the StartTextTranslationJob operation. You can use the DescribeTextTranslationJob or ListTextTranslationJobs operations to get a batch translation job's JobId.

# Required Parameters
- `JobId`: The job ID of the job to be stopped.

"""
stop_text_translation_job(JobId; aws_config::AWSConfig=global_aws_config()) = translate("StopTextTranslationJob", Dict{String, Any}("JobId"=>JobId); aws_config=aws_config)
stop_text_translation_job(JobId, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("StopTextTranslationJob", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("JobId"=>JobId), args)); aws_config=aws_config)

"""
    TranslateText()

Translates input text from the source language to the target language. For a list of available languages and language codes, see what-is-languages.

# Required Parameters
- `SourceLanguageCode`: The language code for the language of the source text. The language must be a language supported by Amazon Translate. For a list of language codes, see what-is-languages. To have Amazon Translate determine the source language of your text, you can specify auto in the SourceLanguageCode field. If you specify auto, Amazon Translate will call Amazon Comprehend to determine the source language.
- `TargetLanguageCode`: The language code requested for the language of the target text. The language must be a language supported by Amazon Translate.
- `Text`: The text to translate. The text string can be a maximum of 5,000 bytes long. Depending on your character set, this may be fewer than 5,000 characters.

# Optional Parameters
- `TerminologyNames`: The name of the terminology list file to be used in the TranslateText request. You can use 1 terminology list at most in a TranslateText request. Terminology lists can contain a maximum of 256 terms.
"""
translate_text(SourceLanguageCode, TargetLanguageCode, Text; aws_config::AWSConfig=global_aws_config()) = translate("TranslateText", Dict{String, Any}("SourceLanguageCode"=>SourceLanguageCode, "TargetLanguageCode"=>TargetLanguageCode, "Text"=>Text); aws_config=aws_config)
translate_text(SourceLanguageCode, TargetLanguageCode, Text, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("TranslateText", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("SourceLanguageCode"=>SourceLanguageCode, "TargetLanguageCode"=>TargetLanguageCode, "Text"=>Text), args)); aws_config=aws_config)

"""
    UpdateParallelData()

Updates a previously created parallel data resource by importing a new input file from Amazon S3.

# Required Parameters
- `ClientToken`: A unique identifier for the request. This token is automatically generated when you use Amazon Translate through an AWS SDK.
- `Name`: The name of the parallel data resource being updated.
- `ParallelDataConfig`: Specifies the format and S3 location of the parallel data input file.

# Optional Parameters
- `Description`: A custom description for the parallel data resource in Amazon Translate.
"""
update_parallel_data(ClientToken, Name, ParallelDataConfig; aws_config::AWSConfig=global_aws_config()) = translate("UpdateParallelData", Dict{String, Any}("ClientToken"=>ClientToken, "Name"=>Name, "ParallelDataConfig"=>ParallelDataConfig); aws_config=aws_config)
update_parallel_data(ClientToken, Name, ParallelDataConfig, args::AbstractDict{String, <:Any}; aws_config::AWSConfig=global_aws_config()) = translate("UpdateParallelData", Dict{String, Any}(mergewith(_merge, Dict{String, Any}("ClientToken"=>ClientToken, "Name"=>Name, "ParallelDataConfig"=>ParallelDataConfig), args)); aws_config=aws_config)


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Functions related to the conversion from Euler angles to DCM.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export angle_to_dcm

"""
    angle_to_dcm(θ₁::Number[, θ₂::Number[, θ₃::Number]], rot_seq::Symbol = :ZYX)
    angle_to_dcm(Θ::EulerAngles)

Create a direction cosine matrix that perform a set of rotations (`θ₁`, `θ₂`,
`θ₃`) about the coordinate axes specified in `rot_seq`.

The input values of the origin Euler angles can also be passed inside the
structure `Θ` (see [`EulerAngles`](@ref)).

The rotation sequence is defined by a `Symbol` specifing the rotation axes. The
possible values depends on the number of rotations as follows:

- **1 rotation** (`θ₁`): `:X`, `:Y`, or `:Z`.
- **2 rotations** (`θ₁`, `θ₂`): `:XY`, `:XZ`, `:YX`, `:YZ`, `:ZX`, or `:ZY`.
- **3 rotations** (`θ₁`, `θ₂`, `θ₃`): `:XYX`, `XYZ`, `:XZX`, `:XZY`, `:YXY`,
    `:YXZ`, `:YZX`, `:YZY`, `:ZXY`, `:ZXZ`, `:ZYX`, or `:ZYZ`

# Remarks

This function assigns `dcm = A3 * A2 * A1` in which `Ai` is the DCM related with
the *i*-th rotation, `i Є [1,2,3]`. If the *i*-th rotation is not specified,
then `Ai = I`.

# Example

```jldoctest
julia> angle_to_dcm(pi / 2, :X)
3×3 StaticArrays.SArray{Tuple{3,3},Float64,2,9}:
 1.0   0.0          0.0
 0.0   6.12323e-17  1.0
 0.0  -1.0          6.12323e-17

julia> angle_to_dcm(pi / 5, pi / 7, :YZ)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  0.728899  0.433884  -0.529576
 -0.351019  0.900969   0.25503
  0.587785  0.0        0.809017

julia> angle_to_dcm(pi / 5, pi / 7, 0, :YZY)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  0.728899  0.433884  -0.529576
 -0.351019  0.900969   0.25503
  0.587785  0.0        0.809017

julia> dcm = angle_to_dcm(pi / 2, pi / 3, pi / 4, :ZYX)
3×3 StaticArrays.SMatrix{3, 3, Float64, 9} with indices SOneTo(3)×SOneTo(3):
  3.06162e-17  0.5       -0.866025
 -0.707107     0.612372   0.353553
  0.707107     0.612372   0.353553

```
"""
function angle_to_dcm(θ::Number, rot_seq::Symbol)
    sa, ca = sincos(θ)

    if rot_seq == :X
        return DCM(
            1,   0,   0,
            0, +ca, +sa,
            0, -sa, +ca
        )'
    elseif rot_seq == :Y
        return DCM(
            +ca, 0, -sa,
              0, 1,   0,
            +sa, 0, +ca
        )'
    elseif rot_seq == :Z
        return DCM(
            +ca, +sa, 0,
            -sa, +ca, 0,
              0,   0, 1
        )'
    else
        throw(ArgumentError("rot_seq must be :X, :Y, or :Z"))
    end
end

function angle_to_dcm(
    θ₁::T1,
    θ₂::T2,
    rot_seq::Symbol
) where {T1<:Number, T2<:Number}
    T = promote_type(T1, T2)

    # Compute the sines and cosines.
    s₁, c₁ = sincos(T(θ₁))
    s₂, c₂ = sincos(T(θ₂))

    if rot_seq == :XY
        return DCM(
              c₂,  s₁ * s₂, -c₁ * s₂,
            T(0),       c₁,       s₁,
              s₂, -s₁ * c₂,  c₁ * c₂
        )'
    elseif rot_seq == :XZ
        return DCM(
              c₂, c₁ * s₂, s₁ * s₂,
             -s₂, c₁ * c₂, s₁ * c₂,
            T(0),     -s₁,      c₁
        )'
    elseif rot_seq == :YX
        return DCM(
                 c₁, T(0),      -s₁,
            s₁ * s₂,   c₂,  c₁ * s₂,
            s₁ * c₂,  -s₂,  c₁ * c₂
        )'
    elseif rot_seq == :YZ
        return DCM(
             c₁ * c₂,   s₂, -s₁ * c₂,
            -c₁ * s₂,   c₂,  s₁ * s₂,
                  s₁, T(0),       c₁
        )'
    elseif rot_seq == :ZX
        return DCM(
                  c₁,       s₁, T(0),
            -c₂ * s₁,  c₂ * c₁,   s₂,
             s₂ * s₁, -s₂ * c₁,   c₂
        )'
    elseif rot_seq == :ZY
        return DCM(
            c₂ * c₁, c₂ * s₁, -s₂ ,
                -s₁,      c₁, T(0),
            s₂ * c₁, s₂ * s₁,  c₂
        )'
    else
        throw(ArgumentError("The rotation sequence :$rot_seq is not valid."))
    end
end

function angle_to_dcm(
    θ₁::T1,
    θ₂::T2,
    θ₃::T3,
    rot_seq::Symbol = :ZYX
) where {T1<:Number, T2<:Number, T3<:Number}
    T = promote_type(T1, T2, T3)

    # Compute the sines and cosines.
    s₁, c₁ = sincos(T(θ₁))
    s₂, c₂ = sincos(T(θ₂))
    s₃, c₃ = sincos(T(θ₃))

    if rot_seq == :ZYX
        return DCM(
                 c₂ * c₁,                c₂ * s₁,             -s₂ ,
            s₃ * s₂ * c₁ - c₃ * s₁, s₃ * s₂ * s₁ + c₃ * c₁, s₃ * c₂,
            c₃ * s₂ * c₁ + s₃ * s₁, c₃ * s₂ * s₁ - s₃ * c₁, c₃ * c₂
        )'
    elseif rot_seq == :XYX
        return DCM(
              c₂,               s₁ * s₂,               -c₁ * s₂,
            s₂ * s₃, -s₁ * c₂ * s₃ + c₁ * c₃, c₁ * c₂ * s₃ + s₁ * c₃,
            s₂ * c₃, -s₁ * c₃ * c₂ - c₁ * s₃, c₁ * c₃ * c₂ - s₁ * s₃
        )'
    elseif rot_seq == :XYZ
        return DCM(
             c₂ * c₃,  s₁ * s₂ * c₃ + c₁ * s₃, -c₁ * s₂ * c₃ + s₁ * s₃,
            -c₂ * s₃, -s₁ * s₂ * s₃ + c₁ * c₃,  c₁ * s₂ * s₃ + s₁ * c₃,
                s₂,             -s₁ * c₂,                 c₁ * c₂
        )'
    elseif rot_seq == :XZX
        return DCM(
               c₂,               c₁ * s₂,                 s₁ * s₂,
            -s₂ * c₃,  c₁ * c₃ * c₂ - s₁ * s₃,  s₁ * c₃ * c₂ + c₁ * s₃,
             s₂ * s₃, -c₁ * c₂ * s₃ - s₁ * c₃, -s₁ * c₂ * s₃ + c₁ * c₃
        )'

    elseif rot_seq == :XZY
        return DCM(
            c₃ * c₂, c₁ * c₃ * s₂ + s₁ * s₃, s₁ * c₃ * s₂ - c₁ * s₃,
              -s₂,             c₁ * c₂,                s₁ * c₂,
            s₃ * c₂, c₁ * s₂ * s₃ - s₁ * c₃, s₁ * s₂ * s₃ + c₁ * c₃
        )'
    elseif rot_seq == :YXY
        return DCM(
            -s₁ * c₂ * s₃ + c₁ * c₃,  s₂ * s₃, -c₁ * c₂ * s₃ - s₁ * c₃,
                       s₁ * s₂,         c₂,              c₁ * s₂,
             s₁ * c₃ * c₂ + c₁ * s₃, -s₂ * c₃,  c₁ * c₃ * c₂ - s₁ * s₃
        )'
    elseif rot_seq == :YXZ
        return DCM(
             c₁ * c₃ + s₂ * s₁ * s₃, c₂ * s₃, -s₁ * c₃ + s₂ * c₁ * s₃,
            -c₁ * s₃ + s₂ * s₁ * c₃, c₂ * c₃,  s₁ * s₃ + s₂ * c₁ * c₃,
                  s₁ * c₂,             -s₂,         c₂ * c₁
        )'
    elseif rot_seq == :YZX
        return DCM(
                       c₁ * c₂,         s₂,              -s₁ * c₂,
            -c₃ * c₁ * s₂ + s₃ * s₁,  c₂ * c₃,  c₃ * s₁ * s₂ + s₃ * c₁,
             s₃ * c₁ * s₂ + c₃ * s₁, -s₃ * c₂, -s₃ * s₁ * s₂ + c₃ * c₁
        )'
    elseif rot_seq == :YZY
        return DCM(
            c₁ * c₃ * c₂ - s₁ * s₃, s₂ * c₃, -s₁ * c₃ * c₂ - c₁ * s₃,
                -c₁ * s₂,             c₂,               s₁ * s₂,
            c₁ * c₂ * s₃ + s₁ * c₃, s₂ * s₃, -s₁ * c₂ * s₃ + c₁ * c₃
        )'
    elseif rot_seq == :ZXY
        return DCM(
            c₃ * c₁ - s₂ * s₃ * s₁, c₃ * s₁ + s₂ * s₃ * c₁, -s₃ * c₂,
                -c₂ * s₁,                c₂ * c₁,              s₂,
            s₃ * c₁ + s₂ * c₃ * s₁, s₃ * s₁ - s₂ * c₃ * c₁,  c₂ * c₃
        )'
    elseif rot_seq == :ZXZ
        return DCM(
            -s₁ * c₂ * s₃ + c₁ * c₃, c₁ * c₂ * s₃ + s₁ * c₃, s₂ * s₃,
            -s₁ * c₃ * c₂ - c₁ * s₃, c₁ * c₃ * c₂ - s₁ * s₃, s₂ * c₃,
                       s₁ * s₂,               -c₁ * s₂,         c₂
        )'
    elseif rot_seq == :ZYZ
        return DCM(
             c₁ * c₃ * c₂ - s₁ * s₃,  s₁ * c₃ * c₂ + c₁ * s₃, -s₂ * c₃,
            -c₁ * c₂ * s₃ - s₁ * c₃, -s₁ * c₂ * s₃ + c₁ * c₃,  s₂ * s₃,
                       c₁ * s₂,                 s₁ * s₂,          c₂
        )'
    else
        throw(ArgumentError("The rotation sequence :$rot_seq is not valid."))
    end
end

angle_to_dcm(Θ::EulerAngles) = angle_to_dcm(Θ.a1, Θ.a2, Θ.a3, Θ.rot_seq)


text:
using CurveProximityQueries
using ConvexBodyProximityQueries
using IntervalArithmetic
using Plots
using StaticArrays

curve = Bernstein([
    [-0.6, 0.0],
    [-0.6, 0.0],
    [0.1, -0.4],
    [0.4, 1.6],
    [-0.2, -1.6],
    [-0.1, 0.4],
    [0.8, 0.0],
    [0.8, 0.0],
])
#ell = Ellipse(curve, 0.2..0.8)

curve2 = Bernstein([[0.6, 0.6], [0.4, 1.0], [0.2, 0.2], [0.0, 0.6]])
#ell2 = Ellipse(curve2, 0.2..0.9)

pts = (curve(0.5), curve2(0.55))

points = curve.(collect(0.0:0.01:1.0))
points2 = curve2.(collect(0.0:0.01:1.0))

plot(first.(points), last.(points), label = "1")
plot!(first.(points2), last.(points2), label = "2")

closest = closest_points(curve,curve2)
scatter!([i for i in first.(closest)], [i for i in last.(closest)])
# plot!(curve, 0.2..0.8,
#       linewidth=2.0,
#       linecolor=:black,
#      )
# plot!(curve, 0.8..1.0,
#       linewidth=2.0,
#       linestyle=:dash,
#       linecolor=Gray(0.5),
#      )
# scatter!(curve(0.2), markersize=4, markercolor=:black)
# scatter!(curve(0.8), markersize=4, markercolor=:black)
#
# plot!(curve2, 0.0..0.2,
#       linewidth=2.0,
#       linestyle=:dash,
#       linecolor=Gray(0.5),
#      )
# plot!(curve2, 0.2..0.9,
#       linewidth=2.0,
#       linecolor=:black,
#      )
# plot!(curve2, 0.9..1.0,
#       linewidth=2.0,
#       linestyle=:dash,
#       linecolor=Gray(0.5),
#      )
# scatter!(curve2(0.2), markersize=4, markercolor=:black)
# scatter!(curve2(0.9), markersize=4, markercolor=:black)
#
# plot!(pts,
#       linestyle=:dash,
#       markersize=8.0,
#       linewidth=2.0,
#      )
# savefig("out/fig5b.pdf")

max_stars_repo_name: QPImaging/xfromProjections

max_issues_repo_name: QPImaging/xfromProjections

max_forks_repo_name: QPImaging/xfromProjections


text:
using QuasiGeostrophy, FFTW, Test
using LinearAlgebra
include(pwd() * "/test/test_utils.jl")
const tol = 1e1

## Define 1D Test
Ω = S¹(0, 2π)
Nx = 2^4;
fourier_grid = FourierGrid(Nx, Ω)
fieldnames = ("ϕ1", "ϕ2", "ϕ3", "ϕ4", "ϕ5")
ϕ1, ϕ2, ϕ3, ϕ4, ϕ5 = create_fields(names = fieldnames, grid = fourier_grid)
# initialize fields with nontrivial data
x = fourier_grid.grid[1]
ϕ1(sin.(x)); ϕ2(cos.(x));
ϕ3( @. sin(x) * cos(x) );
ϕ4( @. sin(x) + cos(x) );
ϕ5( @. cos(x)^2 - sin(x)^2 );

@testset "1D Field Algebra Test" begin
    println("checking ϕ1 * ϕ2 - ϕ3")
    @test norm(ϕ1 * ϕ2 - ϕ3)/norm(ϕ3) < eps(tol)
    println("checking ϕ1 + ϕ2 - ϕ4")
    @test norm(ϕ1 + ϕ2 - ϕ4)/norm(ϕ4) < eps(tol)
end
println(" ")
∂x = create_operators(fourier_grid)

@testset "1D Field Calculus Test" begin
    println("checking ∂x(ϕ1)  - ϕ2")
    @test norm(∂x(ϕ1)  - ϕ2)/norm(ϕ2) < eps(tol)
    println("checking ∂x(ϕ3)  - ϕ5")
    @test norm(∂x(ϕ3)  - ϕ5)/norm(ϕ5) < eps(tol)
end

## Define 2D Test
Ω = S¹(0,2π) × S¹(0,2π)
Nx = 2^4; Ny = 2^4;
fourier_grid = FourierGrid((Nx, Ny), Ω)

fieldnames = ("ϕ1", "ϕ2", "ϕ3", "ϕ4", "ϕ5", "ϕ6")
ϕ1, ϕ2, ϕ3, ϕ4, ϕ5, ϕ6 = create_fields(names = fieldnames, grid = fourier_grid)

x, y = fourier_grid.grid
ϕ1(sin.(x) .+ 0 .* y)
ϕ2(sin.(y) .+ 0 .* x)
ϕ3(sin.(x) .* sin.(y))
ϕ4(sin.(x) .+ sin.(y))
ϕ5(cos.(x) .* sin.(y))
ϕ6(sin.(x) .* cos.(y))

@testset "2D Field Algebra Test" begin
    @test norm(ϕ1 * ϕ2 - ϕ3)/norm(ϕ3) < eps(tol)
    @test norm(ϕ1 + ϕ2 - ϕ4)/norm(ϕ4) < eps(tol)
end

println(" ")
∂x, ∂y = create_operators(fourier_grid)

@testset "2D Field Calculus Test" begin
    println("checking ∂y(ϕ1) ")
    @test norm(∂y(ϕ1)) < eps(tol)
    println("checking ∂x(ϕ3) - ϕ5 ")
    @test norm( ∂x(ϕ3) - ϕ5 )/norm(ϕ5) < eps(tol)
    println("checking ∂y(ϕ3) - ϕ6 ")
    @test norm( ∂y(ϕ3) - ϕ6 )/norm(ϕ6) < eps(tol)
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
using PkgUtility
using PlantHydraulics
using Test

include("test_PlantHydraulics.jl")