text:
"""
    HydrostaticBoussinesqModel <: BalanceLaw

A `BalanceLaw` for ocean modeling.

write out the equations here

ρₒ = reference density of sea water
cʰ = maximum horizontal wave speed
cᶻ = maximum vertical wave speed
αᵀ = thermal expansitivity coefficient
νʰ = horizontal viscosity
νᶻ = vertical viscosity
κʰ = horizontal diffusivity
κᶻ = vertical diffusivity
fₒ = first coriolis parameter (constant term)
β  = second coriolis parameter (linear term)

# Usage

    HydrostaticBoussinesqModel(problem)

"""
struct HydrostaticBoussinesqModel{C, PS, P, MA, TA, F, FT, I} <: BalanceLaw
    param_set::PS
    problem::P
    coupling::C
    momentum_advection::MA
    tracer_advection::TA
    forcing::F
    state_filter::I
    ρₒ::FT
    cʰ::FT
    cᶻ::FT
    αᵀ::FT
    νʰ::FT
    νᶻ::FT
    κʰ::FT
    κᶻ::FT
    κᶜ::FT
    fₒ::FT
    β::FT
    function HydrostaticBoussinesqModel{FT}(
        param_set::PS,
        problem::P;
        coupling::C = Uncoupled(),
        momentum_advection::MA = nothing,
        tracer_advection::TA = NonLinearAdvectionTerm(),
        forcing::F = Forcing(),
        state_filter::I = nothing,
        ρₒ = FT(1000),  # kg / m^3
        cʰ = FT(0),     # m/s
        cᶻ = FT(0),     # m/s
        αᵀ = FT(2e-4),  # (m/s)^2 / K
        νʰ = FT(5e3),   # m^2 / s
        νᶻ = FT(5e-3),  # m^2 / s
        κʰ = FT(1e3),   # m^2 / s # horizontal diffusivity
        κᶻ = FT(1e-4),  # m^2 / s # background vertical diffusivity
        κᶜ = FT(1e-1),  # m^2 / s # diffusivity for convective adjustment
        fₒ = FT(1e-4),  # Hz
        β = FT(1e-11), # Hz / m
    ) where {FT <: AbstractFloat, PS, P, C, MA, TA, F, I}
        return new{C, PS, P, MA, TA, F, FT, I}(
            param_set,
            problem,
            coupling,
            momentum_advection,
            tracer_advection,
            forcing,
            state_filter,
            ρₒ,
            cʰ,
            cᶻ,
            αᵀ,
            νʰ,
            νᶻ,
            κʰ,
            κᶻ,
            κᶜ,
            fₒ,
            β,
        )
    end
end

HBModel = HydrostaticBoussinesqModel

boundary_conditions(ocean::HBModel) = ocean.problem.boundary_conditions

@inline noforcing(args...) = 0

function Forcing(; u = noforcing, v = noforcing, η = noforcing, θ = noforcing)
    return (u = u, v = v, η = η, θ = θ)
end

"""
    vars_state(::HBModel, ::Prognostic)

prognostic variables evolved forward in time

u = (u,v) = (zonal velocity, meridional velocity)
η = sea surface height
θ = temperature
"""
function vars_state(m::HBModel, ::Prognostic, T)
    @vars begin
        u::SVector{2, T}
        η::T # real a 2-D variable TODO: should be 2D
        θ::T
    end
end

"""
    init_state_prognostic!(::HBModel)

sets the initial value for state variables
dispatches to ocean_init_state! which is defined in a problem file such as SimpleBoxProblem.jl
"""
function init_state_prognostic!(m::HBModel, Q::Vars, A::Vars, local_geometry, t)
    return ocean_init_state!(m, m.problem, Q, A, local_geometry, t)
end

"""
    vars_state(::HBModel, ::Auxiliary)
helper variables for computation

second half is because there is no dedicated integral kernels
these variables are used to compute vertical integrals
w = vertical velocity
wz0 = w at z = 0
pkin = bulk hydrostatic pressure contribution

first half of these are fields that are used for computation
y = north-south coordinate

"""
function vars_state(m::HBModel, ::Auxiliary, T)
    @vars begin
        y::T     # y-coordinate of the box
        w::T     # ∫(-∇⋅u)
        pkin::T  # ∫(-αᵀθ)
        wz0::T   # w at z=0
        uᵈ::SVector{2, T}    # velocity deviation from vertical mean
        ΔGᵘ::SVector{2, T}   # vertically averaged tendency
    end
end

function ocean_init_aux! end

"""
    init_state_auxiliary!(::HBModel)

sets the initial value for auxiliary variables (those that aren't related to vertical integrals)
dispatches to ocean_init_aux! which is defined in a problem file such as SimpleBoxProblem.jl
"""
function init_state_auxiliary!(
    m::HBModel,
    state_auxiliary::MPIStateArray,
    grid,
    direction,
)
    init_state_auxiliary!(
        m,
        (m, A, tmp, geom) -> ocean_init_aux!(m, m.problem, A, geom),
        state_auxiliary,
        grid,
        direction,
    )
end

"""
    vars_state(::HBModel, ::Gradient)

variables that you want to take a gradient of
these are just copies in our model
"""
function vars_state(m::HBModel, ::Gradient, T)
    @vars begin
        ∇u::SVector{2, T}
        ∇uᵈ::SVector{2, T}
        ∇θ::T
    end
end

"""
    compute_gradient_argument!(::HBModel)

copy u and θ to var_gradient
this computation is done pointwise at each nodal point

# arguments:
- `m`: model in this case HBModel
- `G`: array of gradient variables
- `Q`: array of state variables
- `A`: array of aux variables
- `t`: time, not used
"""
@inline function compute_gradient_argument!(m::HBModel, G::Vars, Q::Vars, A, t)
    G.∇θ = Q.θ

    velocity_gradient_argument!(m, m.coupling, G, Q, A, t)

    return nothing
end

@inline function velocity_gradient_argument!(
    m::HBModel,
    ::Uncoupled,
    G,
    Q,
    A,
    t,
)
    G.∇u = Q.u

    return nothing
end

"""
    vars_state(::HBModel, ::GradientFlux, FT)

the output of the gradient computations
multiplies ∇u by viscosity tensor and ∇θ by the diffusivity tensor
"""
function vars_state(m::HBModel, ::GradientFlux, T)
    @vars begin
        ∇ʰu::T
        ν∇u::SMatrix{3, 2, T, 6}
        κ∇θ::SVector{3, T}
    end
end

"""
    compute_gradient_flux!(::HBModel)

copy ∇u and ∇θ to var_diffusive
this computation is done pointwise at each nodal point

# arguments:
- `m`: model in this case HBModel
- `D`: array of diffusive variables
- `G`: array of gradient variables
- `Q`: array of state variables
- `A`: array of aux variables
- `t`: time, not used
"""
@inline function compute_gradient_flux!(
    m::HBModel,
    D::Vars,
    G::Grad,
    Q::Vars,
    A::Vars,
    t,
)
    # store ∇ʰu for continuity equation (convert gradient to divergence)
    D.∇ʰu = G.∇u[1, 1] + G.∇u[2, 2]

    velocity_gradient_flux!(m, m.coupling, D, G, Q, A, t)

    κ = diffusivity_tensor(m, G.∇θ[3])
    D.κ∇θ = -κ * G.∇θ

    return nothing
end

@inline function velocity_gradient_flux!(m::HBModel, ::Uncoupled, D, G, Q, A, t)
    ν = viscosity_tensor(m)
    D.ν∇u = -ν * G.∇u

    return nothing
end

"""
    viscosity_tensor(::HBModel)

uniform viscosity with different values for horizontal and vertical directions

# Arguments
- `m`: model object to dispatch on and get viscosity parameters
"""
@inline viscosity_tensor(m::HBModel) = Diagonal(@SVector [m.νʰ, m.νʰ, m.νᶻ])

"""
    diffusivity_tensor(::HBModel)

uniform diffusivity in the horizontal direction
applies convective adjustment in the vertical, bump by 1000 if ∂θ∂z < 0

# Arguments
- `m`: model object to dispatch on and get diffusivity parameters
- `∂θ∂z`: value of the derivative of temperature in the z-direction
"""
@inline function diffusivity_tensor(m::HBModel, ∂θ∂z)
    ∂θ∂z < 0 ? κ = m.κᶜ : κ = m.κᶻ

    return Diagonal(@SVector [m.κʰ, m.κʰ, κ])
end

"""
    vars_integral(::HBModel)

location to store integrands for bottom up integrals
∇hu = the horizontal divegence of u, e.g. dw/dz
"""
function vars_state(m::HBModel, ::UpwardIntegrals, T)
    @vars begin
        ∇ʰu::T
        αᵀθ::T
    end
end

"""
    integral_load_auxiliary_state!(::HBModel)

copy w to var_integral
this computation is done pointwise at each nodal point

arguments:
m -> model in this case HBModel
I -> array of integrand variables
Q -> array of state variables
A -> array of aux variables
"""
@inline function integral_load_auxiliary_state!(
    m::HBModel,
    I::Vars,
    Q::Vars,
    A::Vars,
)
    I.∇ʰu = A.w # borrow the w value from A...
    I.αᵀθ = -m.αᵀ * Q.θ # integral will be reversed below

    return nothing
end

"""
    integral_set_auxiliary_state!(::HBModel)

copy integral results back out to aux
this computation is done pointwise at each nodal point

arguments:
m -> model in this case HBModel
A -> array of aux variables
I -> array of integrand variables
"""
@inline function integral_set_auxiliary_state!(m::HBModel, A::Vars, I::Vars)
    A.w = I.∇ʰu
    A.pkin = I.αᵀθ

    return nothing
end

"""
    vars_reverse_integral(::HBModel)

location to store integrands for top down integrals
αᵀθ = density perturbation
"""
function vars_state(m::HBModel, ::DownwardIntegrals, T)
    @vars begin
        αᵀθ::T
    end
end

"""
    reverse_integral_load_auxiliary_state!(::HBModel)

copy αᵀθ to var_reverse_integral
this computation is done pointwise at each nodal point

arguments:
m -> model in this case HBModel
I -> array of integrand variables
A -> array of aux variables
"""
@inline function reverse_integral_load_auxiliary_state!(
    m::HBModel,
    I::Vars,
    Q::Vars,
    A::Vars,
)
    I.αᵀθ = A.pkin

    return nothing
end

"""
    reverse_integral_set_auxiliary_state!(::HBModel)

copy reverse integral results back out to aux
this computation is done pointwise at each nodal point

arguments:
m -> model in this case HBModel
A -> array of aux variables
I -> array of integrand variables
"""
@inline function reverse_integral_set_auxiliary_state!(
    m::HBModel,
    A::Vars,
    I::Vars,
)
    A.pkin = I.αᵀθ

    return nothing
end

"""
    flux_first_order!(::HBModel)

calculates the hyperbolic flux contribution to state variables
this computation is done pointwise at each nodal point

# arguments:
m -> model in this case HBModel
F -> array of fluxes for each state variable
Q -> array of state variables
A -> array of aux variables
t -> time, not used

# computations
∂ᵗu = ∇⋅(g*η + g∫αᵀθdz + v⋅u)
∂ᵗθ = ∇⋅(vθ) where v = (u,v,w)
"""
@inline function flux_first_order!(
    m::HBModel,
    F::Grad,
    Q::Vars,
    A::Vars,
    t::Real,
    direction,
)
    # ∇ʰ • (g η)
    hydrostatic_pressure!(m, m.coupling, F, Q, A, t)

    # ∇ʰ • (- ∫(αᵀ θ))
    kinematic_pressure!(m, F, Q, A, t)

    # ∇ʰ • (v ⊗ u)
    momentum_advection!(m, m.momentum_advection, F, Q, A, t)

    # ∇ • (u θ)
    tracer_advection!(m, m.tracer_advection, F, Q, A, t)

    return nothing
end

@inline function hydrostatic_pressure!(m::HBModel, ::Uncoupled, F, Q, A, t)
    η = Q.η
    Iʰ = @SMatrix [
        1 -0
        -0 1
        -0 -0
    ]

    F.u += grav(parameter_set(m)) * η * Iʰ

    return nothing
end

@inline function kinematic_pressure!(m::HBModel, F, Q, A, t)
    pkin = A.pkin
    Iʰ = @SMatrix [
        1 -0
        -0 1
        -0 -0
    ]
    F.u += grav(parameter_set(m)) * pkin * Iʰ

    return nothing
end

momentum_advection!(::HBModel, ::Nothing, _...) = nothing
@inline function momentum_advection!(
    ::HBModel,
    ::NonLinearAdvectionTerm,
    F,
    Q,
    A,
    t,
)
    u = Q.u
    @inbounds v = @SVector [Q.u[1], Q.u[2], A.w]

    F.u += v * u'

    return nothing
end

tracer_advection!(::HBModel, ::Nothing, _...) = nothing
@inline function tracer_advection!(
    ::HBModel,
    ::NonLinearAdvectionTerm,
    F,
    Q,
    A,
    t,
)
    θ = Q.θ
    @inbounds v = @SVector [Q.u[1], Q.u[2], A.w]

    F.θ += v * θ

    return nothing
end

"""
    flux_second_order!(::HBModel)

calculates the parabolic flux contribution to state variables
this computation is done pointwise at each nodal point

# arguments:
- `m`: model in this case HBModel
- `F`: array of fluxes for each state variable
- `Q`: array of state variables
- `D`: array of diff variables
- `A`: array of aux variables
- `t`: time, not used

# computations
∂ᵗu = -∇⋅(ν∇u)
∂ᵗθ = -∇⋅(κ∇θ)
"""
@inline function flux_second_order!(
    m::HBModel,
    F::Grad,
    Q::Vars,
    D::Vars,
    HD::Vars,
    A::Vars,
    t::Real,
)
    F.u += D.ν∇u
    F.θ += D.κ∇θ

    return nothing
end

"""
    source!(::HBModel)

Calculates the source term contribution to state variables.
This computation is done pointwise at each nodal point.

Arguments:
    m -> model in this case HBModel
    F -> array of fluxes for each state variable
    Q -> array of state variables
    A -> array of aux variables
    t -> time, not used

Computations:
    ∂ᵗu = -f × u
    ∂ᵗη = w|(z=0)
"""
@inline function source!(
    m::HBModel,
    S::Vars,
    Q::Vars,
    D::Vars,
    A::Vars,
    t::Real,
    direction,
)
    # explicit forcing for SSH
    wz0 = A.wz0
    S.η += wz0

    coriolis_force!(m, m.coupling, S, Q, A, t)

    # Arguments for forcing functions
    # args = y, t, u, v, w, η, θ
    args = tuple(A.y, t, Q.u[1], Q.u[2], A.w, Q.η, Q.θ)

    Su = m.forcing.u(args...)
    Sv = m.forcing.v(args...)

    S.u += @SVector [Su, Sv]
    S.η += m.forcing.η(args...)
    S.θ += m.forcing.θ(args...)

    return nothing
end

@inline function coriolis_force!(m::HBModel, ::Uncoupled, S, Q, A, t)
    # f × u
    f = coriolis_parameter(m, m.problem, A.y)
    u, v = Q.u # Horizontal components of velocity
    S.u -= @SVector [-f * v, f * u]

    return nothing
end

"""
    wavespeed(::HBModel)

calculates the wavespeed for rusanov flux
"""
@inline wavespeed(m::HBModel, n⁻, _...) = abs(SVector(m.cʰ, m.cʰ, m.cᶻ)' * n⁻)

"""
    update_penalty(::HBModel)
    set Δη = 0 when computing numerical fluxes
"""
# We want not have jump penalties on η (since not a flux variable)
function update_penalty!(
    ::RusanovNumericalFlux,
    ::HBModel,
    n⁻,
    λ,
    ΔQ::Vars,
    Q⁻,
    A⁻,
    Q⁺,
    A⁺,
    t,
)
    ΔQ.η = -0

    return nothing
end

filter_state!(Q, filter::Nothing, grid) = nothing
filter_state!(Q, filter, grid) = apply!(Q, UnitRange(1, size(Q, 2)), grid)

"""
    update_auxiliary_state!(::HBModel)

Applies the vertical filter to the zonal and meridional velocities to preserve numerical incompressibility
Applies an exponential filter to θ to anti-alias the non-linear advective term

Doesn't actually touch the aux variables any more, but we need a better filter interface than this anyways
"""
function update_auxiliary_state!(
    dg::DGModel,
    m::HBModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    FT = eltype(Q)
    MD = dg.modeldata

    # `update_aux!` gets called twice, once for the real elements and once for
    # the ghost elements.  Only apply the filters to the real elems.
    if elems == dg.grid.topology.realelems
        # required to ensure that after integration velocity field is divergence free
        vert_filter = MD.vert_filter
        apply!(Q, (:u,), dg.grid, vert_filter, direction = VerticalDirection())

        exp_filter = MD.exp_filter
        apply!(Q, (:θ,), dg.grid, exp_filter, direction = VerticalDirection())

        filter_state!(Q, m.state_filter, dg.grid)
    end

    compute_flow_deviation!(dg, m, m.coupling, Q, t)

    return true
end

@inline compute_flow_deviation!(dg, ::HBModel, ::Uncoupled, _...) = nothing

"""
    update_auxiliary_state_gradient!(::HBModel)

    ∇hu to w for integration
    performs integration for w and pkin (should be moved to its own integral kernels)
    copies down w and wz0 because we don't have 2D structures
"""
function update_auxiliary_state_gradient!(
    dg::DGModel,
    m::HBModel,
    Q::MPIStateArray,
    t::Real,
    elems::UnitRange,
)
    FT = eltype(Q)
    A = dg.state_auxiliary
    D = dg.state_gradient_flux

    # load -∇ʰu as ∂ᶻw
    index_w = varsindex(vars_state(m, Auxiliary(), FT), :w)
    index_∇ʰu = varsindex(vars_state(m, GradientFlux(), FT), :∇ʰu)
    @views @. A.data[:, index_w, elems] = -D.data[:, index_∇ʰu, elems]

    # compute integrals for w and pkin
    indefinite_stack_integral!(dg, m, Q, A, t, elems) # bottom -> top
    reverse_indefinite_stack_integral!(dg, m, Q, A, t, elems) # top -> bottom

    # We are unable to use vars (ie A.w) for this because this operation will
    # return a SubArray, and adapt (used for broadcasting along reshaped arrays)
    # has a limited recursion depth for the types allowed.
    number_aux = number_states(m, Auxiliary())
    index_wz0 = varsindex(vars_state(m, Auxiliary(), FT), :wz0)
    info = basic_grid_info(dg)
    Nqh, Nqk = info.Nqh, info.Nqk
    nelemv, nelemh = info.nvertelem, info.nhorzelem
    nrealelemh = info.nhorzrealelem

    # project w(z=0) down the stack
    data = reshape(A.data, Nqh, Nqk, number_aux, nelemv, nelemh)
    flat_wz0 = @view data[:, end:end, index_w, end:end, 1:nrealelemh]
    boxy_wz0 = @view data[:, :, index_wz0, :, 1:nrealelemh]
    boxy_wz0 .= flat_wz0

    return true
end

"""
    boundary_state!(nf, bc, ::HBModel, args...)

applies boundary conditions for the hyperbolic fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""
@inline boundary_state!(nf, bc, ocean::HBModel, args...) =
    _ocean_boundary_state!(nf, bc, ocean, args...)

#=
"""
    boundary_state!(nf, ::HBModel, args...)

applies boundary conditions for the hyperbolic fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""
@inline function boundary_state!(nf, ocean::HBModel, args...)
    boundary_conditions = ocean.problem.boundary_conditions
    return ocean_boundary_state!(nf, boundary_conditions, ocean, args...)
end
=#

"""
    ocean_boundary_state!(nf, bc::OceanBC, ::HBModel, args...)

splits boundary condition application into velocity and temperature conditions
"""
@inline function ocean_boundary_state!(nf, bc::OceanBC, ocean::HBModel, args...)
    ocean_velocity_boundary_state!(nf, bc.velocity, ocean, args...)
    ocean_temperature_boundary_state!(nf, bc.temperature, ocean, args...)

    return nothing
end

"""
    ocean_boundary_state!(nf, boundaries::Tuple, ::HBModel,
                          Q⁺, A⁺, n, Q⁻, A⁻, bctype)

applies boundary conditions for the first-order and gradient fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""
@generated function ocean_boundary_state!(
    nf::Union{NumericalFluxFirstOrder, NumericalFluxGradient},
    boundaries::Tuple,
    ocean,
    Q⁺,
    A⁺,
    n,
    Q⁻,
    A⁻,
    bctype,
    t,
    args...,
)
    N = fieldcount(boundaries)
    return quote
        Base.Cartesian.@nif(
            $(N + 1),
            i -> bctype == i, # conditionexpr
            i -> ocean_boundary_state!(
                nf,
                boundaries[i],
                ocean,
                Q⁺,
                A⁺,
                n,
                Q⁻,
                A⁻,
                t,
            ), # expr
            i -> error("Invalid boundary tag")
        ) # elseexpr
        return nothing
    end
end

"""
    ocean_boundary_state!(nf, boundaries::Tuple, ::HBModel,
                          Q⁺, A⁺, D⁺, n, Q⁻, A⁻, D⁻, bctype)

applies boundary conditions for the second-order fluxes
dispatches to a function in OceanBoundaryConditions.jl based on bytype defined by a problem such as SimpleBoxProblem.jl
"""
@generated function ocean_boundary_state!(
    nf::NumericalFluxSecondOrder,
    boundaries::Tuple,
    ocean,
    Q⁺,
    D⁺,
    HD⁺,
    A⁺,
    n,
    Q⁻,
    D⁻,
    HD⁻,
    A⁻,
    bctype,
    t,
    args...,
)
    N = fieldcount(boundaries)
    return quote
        Base.Cartesian.@nif(
            $(N + 1),
            i -> bctype == i, # conditionexpr
            i -> ocean_boundary_state!(
                nf,
                boundaries[i],
                ocean,
                Q⁺,
                D⁺,
                A⁺,
                n,
                Q⁻,
                D⁻,
                A⁻,
                t,
            ), # expr
            i -> error("Invalid boundary tag")
        ) # elseexpr
        return nothing
    end
end


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

using Markdown
using InteractiveUtils

# ╔═╡ dc645c7a-98c1-4922-afc6-502cb29a3ffa
using Pkg; Pkg.activate("..")

# ╔═╡ 13c11517-1774-4ba7-add1-b84051860488
using STMOZOO.SudokuSolver, PlutoUI

# ╔═╡ 171fee18-20f6-11eb-37e5-2d04caea8c35
md"""
# Sudoku: Simulated annealing

By Sarah Laperre

## Introduction

Sudoku is probably the best known puzle game. Because simple rules and a variety of puzzles evrywhere, most people know hot solve them. In this project, I will take a look at one algorithm used for solving sudokus: simulated annealing. But first a short introduction to sudoku.

### A short history
Contrary to popular belief, sudoku does not originate from Japan, but from France. In the 19th centuryfirst sudoku-like puzzles appeared in newspapers. It was called "carré magique diabolique", and was a magic square with certain numbers removed. A magic square is a number-filled square where ueach of the rows, columns and diagonals adds up to the same number. It was later renamed to "Number Place" and simplified, resulting in the current rules. The puzzle became popular when it was renamed to sudoku, meaning single number, by a japanese puzzle company. In 2004, a program was made to genrate sudokus with unique solutions. Even now, sudokus are still linked to newspapers.

### The rules
Many sudoku variants exist, with different characters, grid sizes, additional constraints or blocks with different shapes. Here, only the original 9x9 sudoku will be considered.

A sudoku is a 9x9 grid, divided in 9 3x3 regions, called blocks. within the grid, some positions will already be filled. These filled positions will be referred to as fixed positions. The goal is to fill all empty positions, in such a way that every row, column and block contains the numbers 1 till 9 once.

It has been proven that there are at least 17 fixed positions needed for a sudoku to have a uniqe solution. There are $$6.67 * 10^{21}$$ different ways to fill a 9x9 sudoku grid.
"""

# ╔═╡ 38c3a299-b597-40aa-abbf-71afde17f26d
md"""
## Exploring sudokus

Sudokus level 1 till 3 will be used as explanation. Level 4 till 6 can be used to try some different sudokus, and blank can be used as a template for different sudokus.

To use as input, sudokus are represented by an array, with blanks filled by a zero.
"""

# ╔═╡ 7a25aac1-c7b2-4a94-81da-9325c30a0336
lvl1 =  [4 9 2 8 1 7 3 6 5
	8 1 3 5 6 4 7 9 0
	6 7 5 9 2 3 4 8 1
	9 5 4 0 8 1 2 3 7
	2 0 1 3 0 9 6 0 4
	7 3 6 2 4 0 8 1 9
	5 2 9 7 3 6 1 4 8
	0 4 8 1 5 2 9 7 6
	1 6 7 4 9 8 5 2 0]

# ╔═╡ ec046206-68a9-406d-a2ac-11ffdb2dee82
lvl2 =  [7 4 1 5 9 0 3 8 2
	3 0 6 2 4 7 1 9 5
	2 9 5 1 8 3 6 0 7
	9 7 2 0 0 8 0 1 4
	0 0 4 0 5 0 8 0 0
	1 5 0 4 0 0 7 3 6
	5 0 3 9 7 2 4 6 8
	4 6 9 8 3 5 2 0 1
	8 2 7 0 1 4 9 5 3]

# ╔═╡ cd1397a7-30ed-4542-993a-67548895691d
lvl3 =  [0 5 1 0 0 0 8 4 0
	3 8 4 0 5 0 2 7 6
	6 7 2 0 0 0 5 1 9
	1 3 9 7 0 8 4 2 5
	2 0 0 3 1 5 0 0 7
	5 6 7 4 0 2 1 3 8
	4 9 3 0 0 0 6 5 2
	8 2 6 0 3 0 7 9 1
	0 1 5 0 0 0 3 8 0]

# ╔═╡ 51464e74-c439-49f4-802c-9dba8e5fcadc
lvl4 =  [0 7 4 0 0 0 1 6 0
	8 0 6 0 0 0 2 0 5
    0 1 5 0 7 0 8 3 0
	0 0 3 8 0 6 7 0 0
	0 2 0 0 0 0 0 5 0
	0 0 9 2 0 3 4 0 0
	0 3 2 0 6 0 5 8 0
	9 0 8 0 0 0 6 0 1
	0 6 7 0 0 0 9 2 0]

# ╔═╡ 9cc0129e-3122-41b6-8d3f-ed0dd5420da4
lvl5 =  [7 6 3 9 4 1 0 5 8
	9 0 5 0 0 0 0 0 7
	0 0 0 0 0 5 0 6 9
	2 0 4 0 3 0 0 0 5
	5 0 0 7 0 8 0 0 3
	6 0 0 0 5 0 8 0 1
	3 7 0 5 0 0 0 0 0
	8 0 0 0 0 0 1 0 4
	4 9 0 3 8 2 5 7 6]

# ╔═╡ f991f1de-a92d-45e0-8294-44715e616e6c
lvl6 =  [2 0 0 0 5 0 0 0 6
	0 9 0 7 6 0 0 2 0
	0 0 6 0 0 0 3 0 0
	0 0 0 4 3 7 0 5 0
	7 6 0 9 0 5 0 3 8
	0 5 0 6 8 2 0 0 0
	0 0 1 0 0 0 8 0 0
	0 2 0 0 4 9 0 7 0
	3 0 0 0 7 0 0 0 9]

# ╔═╡ 45587b5d-86c9-436f-a573-03a8fdb68ebb
blank = [0 0 0 0 0 0 0 0 0
	0 0 0 0 0 0 0 0 0
	0 0 0 0 0 0 0 0 0
	0 0 0 0 0 0 0 0 0
	0 0 0 0 0 0 0 0 0
	0 0 0 0 0 0 0 0 0
	0 0 0 0 0 0 0 0 0
	0 0 0 0 0 0 0 0 0
	0 0 0 0 0 0 0 0 0]

# ╔═╡ 24d2d6e2-5283-4700-9459-6b84b3f80c4f
md"""
The fuction show_sudoku(sudoku) can be used to make this look more like an actual sudoku.
"""

# ╔═╡ eb18ffab-3c45-4bd6-af6a-1c7b422f14bc
with_terminal() do
	show_sudoku(lvl1)
end

# ╔═╡ 9556bf81-9486-4b00-8f82-18febb29e22e
with_terminal() do
	show_sudoku(lvl2)
end

# ╔═╡ f307ed1d-28a8-4896-9931-f004e11cffd5
with_terminal() do
	show_sudoku(lvl3)
end

# ╔═╡ fb4aeb8c-20f7-11eb-0444-259de7b76883
md"""

## Simulated annealing

Simulated annealing is used to find global optima in large, but discrete search spaces. Its uses include solving the travelling salesman problemSince the goal with sudoku is to find a single solution among many possibilities, simulated annealing can be used as a solver. However, it is better at finding an aproximate global optimum than a pricize local optimum. Simulated annealing falls into somewhat of a grey zone between local and global seach. As such it seemed interesting to find out how well SA performs when solving a sudoku.

Annealing is a method in matalwork, where the speed by which the metal cools is controlled. This results in fewer defects, because the material forms a larger crystal structure.

Similarly, in simulated annealing, you start with a certain temperature. This temperature represents the chance that a step with a worse score (more errors) is accepted. A high temperature means a higher acceptance for mistakes. As the program runs, the temperature gets lower (cools down) and fewer errors are accepted. This means that the search is more global with higher temperatures, but becomes more local as the temperature decreases.

## SA for sudokus

SA works by comparing 2 different neighbouring states. For sudokus, these are randomly filled out grids, with 2 positions swapped in one of the states. To compare the states, the number of errors in rows and columns are counted. If there are less errors in the new state, it is accepted automatically. When the new state has more errors, it is accepted or rejected depending on the temperature.

### Getting started
As such, the first step was to fill the sudoku with random numbers.

The filled sudoku looks as follows:
"""

# ╔═╡ f40d46d6-c022-4ee5-aa06-99439e69f5f9
filledlvl2 = fill_full(lvl2)

# ╔═╡ fdd4e550-20f8-11eb-227b-25f36708484d
with_terminal() do
	show_sudoku(filledlvl2)
end

# ╔═╡ ca79bf63-af4d-46a3-9ef7-ce04c404fcfd
md"""
Notice that each block contains all numbers from 1 till 9. Because it would be cumbersome to count errors within the blocks, the sudoku is randomly filled block by block and swaps occur within the blocks. This way, blocks contain all numbers from the start, and this remains throughout the program.

Since it is not clear in the filled sudoku which positions are fixed, a separate array is formed to keep track of these positions.
"""

# ╔═╡ acb21fa7-64e8-4a6b-ad00-d64b33252cb0
fixedlvl2 = fixed_values(lvl2)

# ╔═╡ 025fd6e8-20f9-11eb-3e7d-3519f3c4b58f
with_terminal() do
	show_sudoku(fixedlvl2)
end

# ╔═╡ 096eff98-20f9-11eb-1e61-99d5714895ba
md"""
Since there are blocks with less than 2 non-fixed positions, we need to identify these positions from the start. In these cases, the block will already be solved by now. From the blocks with more non-fixed positions, 1 will be chosen at random in each iteration. From this block 2 non-fixed positions will be chosen, again at random, and swapped to create a new neighbouring state.
"""

# ╔═╡ 165509ca-20f9-11eb-107c-550cbba0f0e9
fixed_blocks(fixedlvl2)

# ╔═╡ 1fffc82a-20f9-11eb-198c-c160d7dac87d
md"""
For the lvl2 sudoku, there are 3 blocks where positions can be swapped. Blocks are numbered as follows:

1 4 7

2 5 8

3 4 9

Block 1 has no positions that can be swapped and will be solved in het filled sudoku.
"""

# ╔═╡ 49832a8e-20f9-11eb-0841-19a40a12db18
block(1, 1, lvl2)

# ╔═╡ 9c369269-b5af-4fcf-955c-7609f55f43a6
block(1, 1, filledlvl2)

# ╔═╡ 55e0e274-20f9-11eb-36c0-753f228f7e9b
md"""
Comparing fixed positions in blocks 2 and 5 reveals a disadvantage of picking a random block: even though block 2 has fewer positions that can be switched, the block is still picked just as often. This results in a longer runtime.
"""

# ╔═╡ 2a520041-4a9e-458d-8d06-f67b51143aab
block(4, 1, fixedlvl2) # block 2

# ╔═╡ 6731c4f5-c6d3-4efc-a982-80cbe1ffd357
block(4, 4, fixedlvl2) # block 5

# ╔═╡ b1551758-20f9-11eb-3e8f-ff9a7127d7f8
md"""
### Calculating errors
The following function is used to caculate the number of errors in a given state.
"""

# ╔═╡ 45189a82-20fa-11eb-0423-05ce1b84639d
function nr_errors(sudoku)
	errors = 0
  	for i in 1:9
    	errors += (9 - length(freq_table(sudoku[i,:])))
	    errors += (9 - length(freq_table(sudoku[:,i])))
   	end
    return errors
end

# ╔═╡ 41d8f1dc-20fa-11eb-3586-a989427c1fd6
SudokuSolver.nr_errors(filledlvl2)

# ╔═╡ 4ed4215e-20fa-11eb-11ee-f7741591163c
md"""
### Determining starting temperature

A good starting temperature is calculated from the standard deviation of 50 starting states. Since 50 starting states need to be generated anyway, I opted to keep track of these starting states, and return the best one to start from. For easy sudokus like lvl1, this can already solve the sudoku.
"""

# ╔═╡ 6c5473b4-20fa-11eb-327b-51ac560530eb
T0lvl1, bestlvl1 = start_temp(lvl1)

# ╔═╡ d3b42073-d0de-43f3-8627-f9b350fa5167
SudokuSolver.nr_errors(bestlvl1)

# ╔═╡ 252903c4-8aa3-40aa-b9cd-8fbd00c13953
md" For lvl2"

# ╔═╡ 91a5cbbf-dd94-4786-b65b-6755c7f76804
T0lvl2, bestlvl2 = start_temp(lvl2)

# ╔═╡ 3639da42-6a95-4df9-b777-9cf68758a6b0
SudokuSolver.nr_errors(bestlvl2)

# ╔═╡ 1d3d599f-78c3-4d9f-ac2b-e9d81b086836
md" For lvl3"

# ╔═╡ 8b4ea346-9bed-4889-bc7e-e0795a212812
T0lvl3, bestlvl3 = start_temp(lvl3)

# ╔═╡ 134bf884-c7eb-4fb3-b993-71ae40133cd2
SudokuSolver.nr_errors(bestlvl3)

# ╔═╡ 7518c2c0-20fa-11eb-32c0-a9db2a91cbc5
md"""
### Deciding if a worse state is kept
Following function is the core of how SA works. If a state has more errors than the last, this dicides whether or not it is accepted. It is only dependent on the difference in errors (diff) and the current temperature.
"""

# ╔═╡ 34027942-20fb-11eb-261e-3b991ce4c9f8
function decide(diff, T)
   	chance = exp(-diff/T)
   	if rand() < chance
    	return true
   	end
    return false
end

# ╔═╡ 3bbeb85c-20fc-11eb-04d0-fb12d8ace50a
md"""
### Actual solving

The final function looks as follows.
"""

# ╔═╡ 8623ac1a-20fa-11eb-2d45-49cce0fdac86
function sudoku_solver(sudoku)
    if !valid(sudoku)
        error("Input is not a valid sudoku")
    end

	solved = 0
	cooling = 0.99 # Determines the rate by which the temperature decreases.
	stuck = 0
    restart = 0
	show_sudoku(sudoku)
	fixed = fixed_values(sudoku)
	T0, filled = start_temp(sudoku)
    T = T0
    blocks, possible = possible_swap(fixed)
	errors = nr_errors(filled)
	iterations = 81 - sum(fixed) # T0 and iterations are dependent on how many non-fixed values there are.
	state = filled
    best = filled # Keep track of best encountered result.
    bestscore = 90

	if errors == 0
		solved = 1
	end

	while solved == 0
		lastscore = errors
		for i in 1:iterations # Temperature does not decrease for a certain number of iterations.
			newstate = swap(state, fixed, blocks, possible) # Create a new neighbouring state.
			newerrors = nr_errors(newstate)
			if newerrors == 0
				solved = 1
				break
			end

			diff = newerrors - errors
			if diff <= 0 # Better score is always accepted.
				state = newstate
				errors = newerrors
                if errors < bestscore # Remember best result.
                    best = state
                    bestscore = errors
                end
			else
				kept = decide(diff, T) # SA at work.
				if kept
					state = newstate
					errors = newerrors
				end
			end
		end

		T *= cooling # Temperature decreases.
		if errors == 0
			solved = 1
			break
		end

		if errors >= lastscore
			stuck += 1
		else
			stuck = 0
		end

		if stuck > 1000 # When running too long:
            restart +=1
            if restart > 20 # Return to best solution (restart).
                state = best
                errors = bestscore
                restart = 0
				stuck = 0
            else
                T += T0/3 # Or increase temperature.
				stuck = 0
            end
		end
	end

	show_sudoku(state)
	return state
end

# ╔═╡ 48de2108-1677-4277-a33b-a02650636a76
md"""
Over all, it takes a long time to run, depending on how many positions are fixed. SA seems to be quik to find good results, but finding the actual solution is more difficult. The restart seems to be necessary to keep the program on track. Due to the more global searching, sudokus with few errors are quickly found, but for the same reason, this good solution is also quickly discarded. A solution with as few as 2 errors is found quikly, but the actual sulution is harder to find. The cooling rate needs to kept low in order not to miss the solution.
"""

# ╔═╡ 1c4bdb0f-52e8-4eca-a10d-85fee6fa3e1a
with_terminal() do
	SudokuSolver.sudoku_solver(lvl1)
end

# ╔═╡ 702e4676-c461-4588-bccb-b63e30d410e5
with_terminal() do
	SudokuSolver.sudoku_solver(lvl2)
end

# ╔═╡ 67fd3250-f165-4201-9426-7fe89825404d
md"""
## Conclusion

SA is a quik way to solve simple sudokus, but is not a good fit for sudokus with fewer fixed values. Since SA does not keep track of which solutions have been tried already, it ends up trying more (identical) combinations than other algorithms use for sudoku solving. Better options would be a branch-and-bound method, or brute-force search with backtracking. Both keep track of arlready tried solutions and as such, dont try the same one repeatedly. Given that SA does reach a low number of errors relatively quikly, perhaps a hybrid approach would work?
"""

# ╔═╡ bac46fbd-5c17-4418-b59d-136bac701b63
# with_terminal() do
# 	SudokuSolver.sudoku_solver(lvl3)
# end

# ╔═╡ Cell order:
# ╟─171fee18-20f6-11eb-37e5-2d04caea8c35
# ╠═dc645c7a-98c1-4922-afc6-502cb29a3ffa
# ╠═13c11517-1774-4ba7-add1-b84051860488
# ╟─38c3a299-b597-40aa-abbf-71afde17f26d
# ╟─7a25aac1-c7b2-4a94-81da-9325c30a0336
# ╟─ec046206-68a9-406d-a2ac-11ffdb2dee82
# ╟─cd1397a7-30ed-4542-993a-67548895691d
# ╟─51464e74-c439-49f4-802c-9dba8e5fcadc
# ╟─9cc0129e-3122-41b6-8d3f-ed0dd5420da4
# ╟─f991f1de-a92d-45e0-8294-44715e616e6c
# ╟─45587b5d-86c9-436f-a573-03a8fdb68ebb
# ╟─24d2d6e2-5283-4700-9459-6b84b3f80c4f
# ╠═eb18ffab-3c45-4bd6-af6a-1c7b422f14bc
# ╠═9556bf81-9486-4b00-8f82-18febb29e22e
# ╠═f307ed1d-28a8-4896-9931-f004e11cffd5
# ╟─fb4aeb8c-20f7-11eb-0444-259de7b76883
# ╠═f40d46d6-c022-4ee5-aa06-99439e69f5f9
# ╠═fdd4e550-20f8-11eb-227b-25f36708484d
# ╟─ca79bf63-af4d-46a3-9ef7-ce04c404fcfd
# ╠═acb21fa7-64e8-4a6b-ad00-d64b33252cb0
# ╠═025fd6e8-20f9-11eb-3e7d-3519f3c4b58f
# ╟─096eff98-20f9-11eb-1e61-99d5714895ba
# ╠═165509ca-20f9-11eb-107c-550cbba0f0e9
# ╟─1fffc82a-20f9-11eb-198c-c160d7dac87d
# ╠═49832a8e-20f9-11eb-0841-19a40a12db18
# ╠═9c369269-b5af-4fcf-955c-7609f55f43a6
# ╟─55e0e274-20f9-11eb-36c0-753f228f7e9b
# ╠═2a520041-4a9e-458d-8d06-f67b51143aab
# ╠═6731c4f5-c6d3-4efc-a982-80cbe1ffd357
# ╟─b1551758-20f9-11eb-3e8f-ff9a7127d7f8
# ╠═45189a82-20fa-11eb-0423-05ce1b84639d
# ╠═41d8f1dc-20fa-11eb-3586-a989427c1fd6
# ╟─4ed4215e-20fa-11eb-11ee-f7741591163c
# ╠═6c5473b4-20fa-11eb-327b-51ac560530eb
# ╠═d3b42073-d0de-43f3-8627-f9b350fa5167
# ╟─252903c4-8aa3-40aa-b9cd-8fbd00c13953
# ╠═91a5cbbf-dd94-4786-b65b-6755c7f76804
# ╠═3639da42-6a95-4df9-b777-9cf68758a6b0
# ╟─1d3d599f-78c3-4d9f-ac2b-e9d81b086836
# ╠═8b4ea346-9bed-4889-bc7e-e0795a212812
# ╠═134bf884-c7eb-4fb3-b993-71ae40133cd2
# ╟─7518c2c0-20fa-11eb-32c0-a9db2a91cbc5
# ╠═34027942-20fb-11eb-261e-3b991ce4c9f8
# ╟─3bbeb85c-20fc-11eb-04d0-fb12d8ace50a
# ╠═8623ac1a-20fa-11eb-2d45-49cce0fdac86
# ╟─48de2108-1677-4277-a33b-a02650636a76
# ╠═1c4bdb0f-52e8-4eca-a10d-85fee6fa3e1a
# ╠═702e4676-c461-4588-bccb-b63e30d410e5
# ╟─67fd3250-f165-4201-9426-7fe89825404d
# ╠═bac46fbd-5c17-4418-b59d-136bac701b63

ext: jl


text:
"""
    JobQueue

An object representing and AWS batch job queue.

See [`AWSBatch.create_job_queue`](@ref).
"""
struct JobQueue
    arn::String

    function JobQueue(queue::AbstractString; aws_config::AbstractAWSConfig=global_aws_config())
        arn = job_queue_arn(queue; aws_config=aws_config)
        arn === nothing && error("No job queue ARN found for: $queue")
        new(arn)
    end
end

Base.:(==)(a::JobQueue, b::JobQueue) = a.arn == b.arn
function describe(queue::JobQueue; aws_config::AbstractAWSConfig=global_aws_config())
    return describe_job_queue(queue; aws_config=aws_config)
end
function describe_job_queue(queue::JobQueue; aws_config::AbstractAWSConfig=global_aws_config())
    return describe_job_queue(queue.arn; aws_config=aws_config)
end
function max_vcpus(queue::JobQueue; aws_config::AbstractAWSConfig=global_aws_config())
    sum(max_vcpus(ce; aws_config=aws_config) for ce in compute_environments(queue; aws_config=aws_config))
end

function _create_compute_environment_order(envs)
    map(enumerate(envs)) do (i, env)
        Dict{String,Any}("computeEnvironment"=>env, "order"=>i)
    end
end

"""
    create_job_queue(name, envs, priority=1; aws_config=global_aws_config())

Create a job queue with name `name` and priority `priority` returning the associated `JobQueue` object.
`envs` must be an iterator of compute environments given by ARN.

See the AWS docs [here](https://docs.aws.amazon.com/batch/latest/APIReference/API_CreateJobQueue.html).
"""
function create_job_queue(name::AbstractString, envs, priority::Integer=1;
                          enabled::Bool=true,
                          tags::AbstractDict=Dict{String,Any}(),
                          aws_config::AbstractAWSConfig=global_aws_config())
    env = _create_compute_environment_order(envs)
    args = Dict{String,Any}()
    enabled || (args["state"] = "DISABLED")
    isempty(tags) || (args["tags"] = tags)
    return @mock Batch.create_job_queue(env, name, priority, args; aws_config=aws_config)
end

"""
    list_job_queues(;aws_config=global_aws_config())

Get a list of `JobQueue` objects as returned by `Batch.describe_job_queues()`.
"""
function list_job_queues(;aws_config::AbstractAWSConfig=global_aws_config())
    [JobQueue(q["jobQueueArn"]) for q ∈ Batch.describe_job_queues(;aws_config=aws_config)["jobQueues"]]
end

"""
    compute_environments(queue::JobQueue; aws_config=global_aws_config())

Get a list of `ComputeEnvironment` objects associated with the `JobQueue`.
"""
function compute_environments(queue::JobQueue; aws_config::AbstractAWSConfig=global_aws_config())
    ce_order = describe(queue; aws_config=aws_config)["computeEnvironmentOrder"]

    compute_envs = Vector{ComputeEnvironment}(undef, length(ce_order))
    for ce in ce_order
        i, arn = ce["order"], ce["computeEnvironment"]
        compute_envs[i] = ComputeEnvironment(arn)
    end

    return compute_envs
end

function job_queue_arn(queue::AbstractString; aws_config::AbstractAWSConfig=global_aws_config())
    startswith(queue, "arn:") && return queue
    json = describe_job_queue(queue; aws_config=aws_config)
    isempty(json) ? nothing : json["jobQueueArn"]
end

function describe_job_queue(queue::AbstractString;
                            aws_config::AbstractAWSConfig=global_aws_config())::OrderedDict
    json = @mock Batch.describe_job_queues(Dict("jobQueues" => [queue]); aws_config=aws_config)
    queues = json["jobQueues"]
    len = length(queues)::Int
    @assert len <= 1
    return len == 1 ? first(queues) : OrderedDict()
end


text:
# K-mer
# =====
#
# Compact k-mer sequence type.
#
# A Kmer is a sequence ≤ 32nt, without any 'N's, packed in a single 64 bit
# value.  While BioSequence is an efficient general-purpose sequence
# representation, Kmer is useful for applications like assembly, k-mer counting,
# k-mer based quantification in RNA-Seq, etc that rely on manipulating many
# short sequences as efficiently (space and time) as possible.
#
# This file is a part of BioJulia.
# License is MIT: https://github.com/BioJulia/Bio.jl/blob/master/LICENSE.md

# Representation
# --------------
#
# Four kinds of nucleotides are encoded as follows:
#
#   nucleotide | binary
#   ---------- | ------
#       A      |   00
#       C      |   01
#       G      |   10
#     T / U    |   11
#
# NucleicAcids are filled from MSBs to LSBs and right-aligned so that all k-mers
# are lexicographically ordered. For example, the memory layout of "TACG" is:
#   64-bit: 0b 00 00 … 00 11 00 01 10
#    4-mer:                T  A  C  G

bitstype 64 Kmer{T<:NucleicAcid, K} <: Sequence

typealias DNAKmer{K} Kmer{DNA, K}
typealias RNAKmer{K} Kmer{RNA, K}
typealias DNACodon DNAKmer{3}
typealias RNACodon RNAKmer{3}

function Kmer{T<:NucleicAcid}(nts::T...)
    return make_kmer(nts)
end

function DNACodon(x::DNA, y::DNA, z::DNA)
    return make_kmer((x, y, z))
end

function RNACodon(x::RNA, y::RNA, z::RNA)
    return make_kmer((x, y, z))
end

# Conversion
# ----------

function Base.convert{K}(::Type{DNAKmer{K}}, x::UInt64)
    mask = ~UInt64(0) >> (64 - 2K)
    return reinterpret(DNAKmer{K}, x & mask)
end
function Base.convert{K}(::Type{RNAKmer{K}}, x::UInt64)
    mask = ~UInt64(0) >> (64 - 2K)
    return reinterpret(RNAKmer{K}, x & mask)
end
Base.convert{K}(::Type{UInt64}, x::DNAKmer{K}) = reinterpret(UInt64, x)
Base.convert{K}(::Type{UInt64}, x::RNAKmer{K}) = reinterpret(UInt64, x)
Base.convert{K}(::Type{DNAKmer{K}}, x::RNAKmer{K}) = reinterpret(DNAKmer{K}, x)
Base.convert{K}(::Type{RNAKmer{K}}, x::DNAKmer{K}) = reinterpret(RNAKmer{K}, x)

function Base.convert{T,K}(::Type{Kmer{T,K}}, seq::AbstractString)
    return make_kmer(Kmer{T,K}, seq)
end

function Base.convert{T,K,A<:DNAAlphabet}(::Type{Kmer{T,K}}, seq::BioSequence{A})
    return make_kmer(Kmer{DNA,K}, seq)
end

function Base.convert{T,K,A<:RNAAlphabet}(::Type{Kmer{T,K}}, seq::BioSequence{A})
    return make_kmer(Kmer{RNA,K}, seq)
end

Base.convert{T}(::Type{Kmer{T}}, seq::AbstractString) = convert(Kmer{T,length(seq)}, seq)
Base.convert{A<:DNAAlphabet}(::Type{Kmer}, seq::BioSequence{A}) = convert(Kmer{DNA,length(seq)}, seq)
Base.convert{A<:RNAAlphabet}(::Type{Kmer}, seq::BioSequence{A}) = convert(Kmer{RNA,length(seq)}, seq)
Base.convert{A<:DNAAlphabet}(::Type{DNAKmer}, seq::BioSequence{A}) = convert(DNAKmer{length(seq)}, seq)
Base.convert{A<:RNAAlphabet}(::Type{RNAKmer}, seq::BioSequence{A}) = convert(RNAKmer{length(seq)}, seq)

# create a kmer from a sequence whose elements are convertible to a nucleotide
function make_kmer{T,K}(::Type{Kmer{T,K}}, seq)
    seqlen = length(seq)
    if seqlen > 32
        throw(ArgumentError("cannot create a k-mer loger than 32nt"))
    elseif seqlen != K
        throw(ArgumentError("cannot create a $(K)-mer from a sequence of length $(seqlen)"))
    end

    x = UInt64(0)
    for c in seq
        nt = convert(T, c)
        if isambiguous(nt)
            error("A K-mer cannot contain an ambiguous nucleotide in its sequence")
        end
        x = (x << 2) | UInt64(trailing_zeros(nt))
    end

    return Kmer{T,K}(x)
end

make_kmer{K,T}(seq::NTuple{K,T}) = make_kmer(Kmer{T,K}, seq)

Base.convert{K}(::Type{BioSequence}, x::DNAKmer{K}) = DNASequence(x)
Base.convert{K}(::Type{BioSequence}, x::RNAKmer{K}) = RNASequence(x)
Base.convert{A<:DNAAlphabet,K}(::Type{BioSequence{A}}, x::DNAKmer{K}) = BioSequence{A}([nt for nt in x])
Base.convert{A<:RNAAlphabet,K}(::Type{BioSequence{A}}, x::RNAKmer{K}) = BioSequence{A}([nt for nt in x])
Base.convert{S<:AbstractString}(::Type{S}, seq::Kmer) = convert(S, [Char(x) for x in seq])

# Basic Functions
# ---------------

alphabet{k}(::Type{DNAKmer{k}}) = (DNA_A, DNA_C, DNA_G, DNA_T)
alphabet{k}(::Type{RNAKmer{k}}) = (RNA_A, RNA_C, RNA_G, RNA_U)

Base.hash(x::Kmer, h::UInt) = hash(UInt64(x), h)

Base.length{T,K}(x::Kmer{T, K}) = K
Base.eltype{T,k}(::Type{Kmer{T,k}}) = T

@inline function inbounds_getindex{T,K}(x::Kmer{T,K}, i::Integer)
    return reinterpret(T, 0x01 << ((UInt64(x) >> 2(K - i)) & 0b11))
end

Base.summary{k}(x::DNAKmer{k}) = string("DNA ", k, "-mer")
Base.summary{k}(x::RNAKmer{k}) = string("RNA ", k, "-mer")

Base.:-{T,K}(x::Kmer{T,K}, y::Integer) = Kmer{T,K}(UInt64(x) - y % UInt64)
Base.:+{T,K}(x::Kmer{T,K}, y::Integer) = Kmer{T,K}(UInt64(x) + y % UInt64)
Base.:+{T,K}(x::Integer, y::Kmer{T,K}) = y + x
Base.:(==){T,k}(x::Kmer{T,k}, y::Kmer{T,k}) = UInt64(x) == UInt64(y)
Base.isless{T,K}(x::Kmer{T,K}, y::Kmer{T,K}) = isless(UInt64(x), UInt64(y))

# Other functions
# ---------------

"""
    complement(kmer::Kmer)

Return the complement of `kmer`.
"""
complement{T,k}(x::Kmer{T,k}) = Kmer{T,k}(~UInt64(x))

"""
    reverse(kmer::Kmer)

Return the reverse of `kmer`.
"""
Base.reverse{T,k}(x::Kmer{T,k}) = Kmer{T,k}(nucrev2(UInt64(x)) >> (64 - 2k))

"""
    reverse_complement(kmer::Kmer)

Return the reverse complement of `kmer`
"""
reverse_complement(x::Kmer) = complement(reverse(x))

"""
    mismatches(a::Kmer, b::Kmer)

Return the number of mismatches between `a` and `b`.
"""
mismatches{T,k}(a::Kmer{T,k}, b::Kmer{T,k}) = bitpar_mismatches2(UInt64(a), UInt64(b))

"""
    canonical(kmer::Kmer)

Return the canonical k-mer of `x`.

A canonical k-mer is the numerical lesser of a k-mer and its reverse complement.
This is useful in hashing/counting k-mers in data that is not strand specific,
and thus observing k-mer is equivalent to observing its reverse complement.
"""
function canonical(x::Kmer)
    y = reverse_complement(x)
    return x < y ? x : y
end

function Base.rand{T,k}(::Type{Kmer{T,k}})
    return convert(Kmer{T,k}, rand(UInt64))
end

function Base.rand{T,k}(::Type{Kmer{T,k}}, size::Integer)
    return [rand(Kmer{T,k}) for _ in 1:size]
end

# K-mer neighbor
# --------------

# neighbors on a de Bruijn graph
immutable KmerNeighborIterator{T, K}
    x::Kmer{T, K}
end

"""
    neighbors(kmer::Kmer)

Return an iterator through k-mers neighboring `kmer` on a de Bruijn graph.
"""
neighbors(x::Kmer) = KmerNeighborIterator(x)

Base.length(::KmerNeighborIterator) = 4
Base.eltype{T,k}(::Type{KmerNeighborIterator{T,k}}) = Kmer{T,k}
Base.start(it::KmerNeighborIterator) = UInt64(0)
Base.done(it::KmerNeighborIterator, i) = i == 4
function Base.next{T, K}(it::KmerNeighborIterator{T, K}, i)
    return Kmer{T,K}((UInt64(it.x) << 2) | i), i + 1
end

# Counters
# --------

function gc_content{T,k}(kmer::Kmer{T,k})
    if k == 0
        return 0.0
    else
        return (count_g(kmer) + count_c(kmer)) / k
    end
end

function count_a{T,k}(kmer::Kmer{T,k})
    return count_a(reinterpret(UInt64, kmer)) - (32 - k)
end

function count_c{T,k}(kmer::Kmer{T,k})
    return count_c(reinterpret(UInt64, kmer))
end

function count_g{T,k}(kmer::Kmer{T,k})
    return count_g(reinterpret(UInt64, kmer))
end

function count_t{T,k}(kmer::Kmer{T,k})
    return count_t(reinterpret(UInt64, kmer))
end

# Count A, C, T/U, G respectively in a kmer stored in a UInt64
function count_a(x::UInt64)
    xinv = ~x
    return count_ones(((xinv >>> 1) & xinv) & 0x5555555555555555)
end
count_c(x::UInt64) = count_ones((((~x) >>> 1) & x) & 0x5555555555555555)
count_g(x::UInt64) = count_ones(((x >>> 1) & (~x)) & 0x5555555555555555)
count_t(x::UInt64) = count_ones((x    & (x >>> 1)) & 0x5555555555555555)

# Shuffle
# -------

function Base.shuffle{T,k}(kmer::Kmer{T,k})
    # Fisher-Yates shuffle
    for i in 1:k-1
        j = rand(i:k)
        kmer = swap(kmer, i, j)
    end
    return kmer
end

# Swap two nucleotides at `i` and `j`.
function swap{T,k}(kmer::Kmer{T,k}, i, j)
    i = 2k - 2i
    j = 2k - 2j
    b = convert(UInt64, kmer)
    x = ((b >> i) $ (b >> j)) & UInt64(0x03)
    return Kmer{T,k}(b $ ((x << i) | (x << j)))
end


text:
using Test
using Dorothy

@DorothyAll()

@testset "Models" begin

	@testset "Header" begin

		h1 = MolecularModelHeader()
		@test isempty(keys(h1))
		@test_throws Exception h1.smurf = true
		h1.time = 10.0
		@test h1.time == 10.0
		@test get(h1, :time, nothing) == 10.0
		@test length(keys(h1)) == 1
		h1.time = 10
		@test h1.time == 10.0
		@test length(keys(h1)) == 1
		@test_throws Exception h1.time = "ten dot zero"
		@test length(keys(h1)) == 1
		delete!(h1, :time)
		@test isempty(keys(h1))
		h1.title = "something"
		@test_throws Exception h1.title = nothing
		@test length(keys(h1)) == 1
		h1.step = 5
		@test_throws Exception h1.step = 2.5
		@test length(keys(h1)) == 2
		h1.lambda = 0.6
		@test length(keys(h1)) == 3
		h1.cell = OrthorhombicCell(25.0)
		@test dims(h1.cell) ≈ [25.0, 25.0, 25.0]
		@test_throws Exception h1.cell = zeros(3,4)
		h1.virial = zeros(3,3)
		h1.pressure = ones(3,3)
		show(IOBuffer(), h1)
		show(IOBuffer(), MIME"text/plain"(), h1)
		@test h1.lambda == h1[:lambda]
		@test !haskey(h1, :smurf)
		@test get(h1, :lambda, nothing) == 0.6
		@test get(h1, :smurf, nothing) == nothing
		@test get(() -> nothing, h1, :smurf) == nothing
		delete!(h1, :time)
		@test get!(h1, :time, 5.0) == 5.0
		@test h1.time == 5.0
		delete!(h1, :lambda)
		@test get!(() -> 1.0, h1, :lambda) == 1.0
		@test h1.lambda == 1.0

	end

	@testset "MolecularModel" begin

		m1 = MolecularModel()
		m2 = MolecularModel(10)
		@test isempty(m1)
		@test length(m2) == 10
		@test isempty(keys(m1))
		@test isempty(keys(m2))
		show(IOBuffer(), m2)
		show(IOBuffer(), MIME"text/plain"(), m2)
		@test firstindex(m2) == 1
		@test lastindex(m2) == 10
		@test collect(eachindex(m2)) == collect(1:10)
		@test_throws Exception m2.charges = fill("smurf", 10)
		m2.charges = fill(10.0, 10)
		@test m2.charges[1] == 10.0
		@test m2[1].charge == 10.0
		@test length(keys(m2)) == 1
		@test_throws Exception m2.occupancies = fill("tapir", 10)
		m2.occupancies = fill(true, 10)
		@test_throws Exception m2.bfactors = fill(nothing, 10)
		m2.bfactors = fill(10, 10)
		@test length(keys(m2)) == 3
		delete!(m2, :bfactors)
		@test length(keys(m2)) == 2
		delete!(m2, :charges)
		delete!(m2, :occupancies)
		@test isempty(keys(m2))
		m2.R = fill([0,0,0], 10)
		@test_throws Exception m2.R = fill(nothing, 10)
		m2.V = fill([0,0,0], 10)
		@test_throws Exception m2.V = fill(nothing, 10)
		m2.F = fill([0,0,0], 10)
		@test_throws Exception m2.F = fill(nothing, 10)
		@test_throws Exception m2.ids = fill(0, (10,3))
		@test_throws Exception m2.ids = fill(nothing, 10)
		m2.ids = 1:10
		@test_throws Exception m2.resids = fill(0, (10,3))
		@test_throws Exception m2.resids = fill(nothing, 10)
		m2.resids = fill(1, 10)
		@test_throws Exception m2.names = fill(nothing, 10)
		m2.names = fill("cat", 10)
		@test_throws Exception m2.resnames = fill(nothing, 10)
		m2.resnames = fill("kitten", 10)
		@test_throws Exception m2.chainids = fill(nothing, 10)
		m2.chainids = fill("A", 10)
		@test_throws Exception m2.invalid = fill(nothing, 10)
		@test haskey(m2, :resnames) == true
		@test haskey(m2, :tapir) == false
		@test getkey(m2, :resnames, nothing) == :resnames
		@test getkey(m2, :tapir, nothing) == nothing
		@test get(m2, :resnames, nothing) == fill("kitten", 10)
		@test get(m2, :tapir, nothing) == nothing
		@test_throws Exception m2.tapir = fill("kitten", 10)
		@test get!(m2, :resnames, nothing) == fill("kitten", 10)
		@test get!(m2, :occupancies, fill(1.0, 10)) == fill(1.0, 10)
		get!(m2, :R, undef)
	end

	@testset "Particle" begin
		m1 = MolecularModel(10)
		for (i, p) in enumerate(m1)
			@test p.i == i
		end
		m1.R = fill([0,0,0], 10)
		m1.ids = 1:10
		@test m1[1].R == [0.0, 0.0, 0.0]
		@test m1[2].id == 2
		show(IOBuffer(), m1[3])
		p1 = m1[4]
		@test length(keys(p1)) == 2
		@test_throws Exception p1.name = "CA"
		@test haskey(p1, :id) == true
		@test haskey(p1, :tapir) == false
		@test getkey(p1, :id, nothing) == :id
		@test getkey(p1, :tapir, nothing) == nothing
		@test get(p1, :id, nothing) == 4
		@test get(p1, :tapir, nothing) == nothing
		@test p1.id == 4
		p1.id = 3
	end

	@testset "Views" begin
		m1 = MolecularModel(10)
		v1 = view(m1, 1:5)
		@test parent(v1) == m1
		@test parentindices(v1) == 1:5
		m1.R = fill([0,0,0], 10)
		m1.ids = 1:10
		@test v1.ids == 1:5
		@test v1.R == fill([0,0,0], 5)
		fill!(v1.R, [1,1,1])
		@test v1.R[1:5] == fill([1,1,1], 5)
		@test length(keys(v1)) == 2
		@test_throws Exception v1.resids = 1:5
		@test haskey(v1, :ids) == true
		@test haskey(v1, :tapir) == false
		@test getkey(v1, :ids, nothing) == :ids
		@test getkey(v1, :tapir, nothing) == nothing
		@test get(v1, :ids, nothing) == 1:5
		@test get(v1, :tapir, nothing) == nothing
	end

	@testset "Collection" begin

		@testset "Copy" begin
			m1 = MolecularModel(20)
			m1.ids = 1:20
			m1.header.title = "capybara"
			m2 = MolecularModel(m1)
			@test m2.header.title == "capybara"
			m3 = MolecularModel(view(m1, 6:10))
			@test m3.ids == 6:10
			@test m3.header.title == "capybara"
		end

		@testset "Delete" begin
			m1 = MolecularModel(10)
			m1.ids = 1:10
			@test_throws Exception deleteat!(m, 11)
			m1.topology = []
			pair!(m1.topology, (3, 4))
			pair!(m1.topology, (5, 6))
			deleteat!(m1, 3)
			@test pairs(m1.topology) == [(4,5)]
			@test m1.ids == [1:2..., 4:10...]
			deleteat!(m1, 3:5)
			@test isempty(pairs(m1.topology))
			@test m1.ids == [1:2..., 7:10...]
			pair!(m1.topology, (1, 2))
			empty!(m1)
			@test isempty(m1)
			@test isempty(pairs(m1.topology))
		end

		@testset "Splice" begin
			m1 = MolecularModel(10)
			m1.ids = 1:10
			m1.names = fill("A", 10)
			m2 = MolecularModel(10)
			m2.ids = 11:20
			splice!(m1, 2:3)
			@test length(m1) == 8
			@test m1.ids == [1, 4:10...]
			@test_throws Exception insert!(m1, 2, m2)
			m2.names = fill("B", 10)
			insert!(m1, 2, m2)
			@test length(m1) == 18
			@test m1.ids == [1, 11:20..., 4:10...]
			m1.topology = []
			pair!(m1.topology, (1, 2))
			m2.topology = []
			prepend!(m1, m2[1])
			@test length(m1) == 19
			@test ! ((1, 2) in m1.topology)
			@test (2, 3) in m1.topology
			pair!(m2.topology, (1, 2))
			pair!(m2.topology, (2, 3))
			append!(m1, view(m2, 1:5))
			@test length(m1) == 24
			@test (20, 21) in m1.topology
			@test (21, 22) in m1.topology
		end

		@testset "Resize" begin
			m1 = MolecularModel()
			resize!(m1, 10)
			m1.topology = []
			pair!(m1.topology, (6, 7))
			m1.ids = 1:10
			resize!(m1, 20)
			@test length(m1) == 20
			@test length(pairs(m1.topology)) == 1
			@test m1.ids[5] == 5
			resize!(m1, 5)
			@test length(m1) == 5
			@test isempty(pairs(m1.topology))
		end

	end

end # @testset


text:
using Test
using MPI
using CLIMA.Mesh.BrickMesh
using Random

function main()
  MPI.Init()
  comm = MPI.COMM_WORLD
  crank = MPI.Comm_rank(comm)
  csize = MPI.Comm_size(comm)

  Nelemtotal = 113
  Random.seed!(1234);
  globalcode = randperm(Nelemtotal)

  @assert csize > 1

  bs = [(i==1) ? (1:0) :
        BrickMesh.linearpartition(Nelemtotal, i-1, csize-1) for i=1:csize]
  as = [BrickMesh.linearpartition(Nelemtotal, i, csize) for i=1:csize]

  codeb = globalcode[bs[crank+1]]

  (so, ss, rs) = BrickMesh.getpartition(comm, codeb)

  codeb = codeb[so]
  codec = []
  for r = 0:csize-1
    sendrange = ss[r+1]:ss[r+2]-1
    rcounts = MPI.Allgather(Cint(length(sendrange)), comm)
    c = MPI.Gatherv(view(codeb, sendrange), rcounts, r, comm)
    if r == crank
      codec = c
    end
  end

  codea = (1:Nelemtotal)[as[crank+1]]

  @test sort(codec) == codea
end

main()

max_stars_repo_name: akshaysridhar/CLIMA

max_issues_repo_name: akshaysridhar/CLIMA

max_forks_repo_name: akshaysridhar/CLIMA


text:
# Action
# ======

struct Action
    name::Symbol
    order::Int
end

# An action list type.
# * Action names are unique.
# * Actions are sorted by its order.
struct ActionList
    actions::Vector{Action}

    function ActionList(actions::Vector{Action}=Action[])
        list = new(Action[])
        union!(list, actions)
        return list
    end
end

function Base.:(==)(l1::ActionList, l2::ActionList)
    if length(l1) != length(l2)
        return false
    end
    for i in 1:lastindex(l1.actions)
        if l1.actions[i].name != l2.actions[i].name
            return false
        end
    end
    return true
end

function Base.hash(list::ActionList, h::UInt)
    for a in list
        h = xor(h, hash(a.name))
    end
    return h
end

function Base.push!(list::ActionList, action::Action)
    i = findfirst(a -> a.name == action.name, list.actions)
    if i !== nothing
        if action.order < list.actions[i].order
            list.actions[i] = action  # replace
        end
    else
        push!(list.actions, action)
    end
    sort!(list.actions, by=a->a.order)
    return list
end

function Base.union!(list::ActionList, actions::Union{ActionList,Vector{Action}})
    for a in actions
        push!(list, a)
    end
    return list
end

function Base.length(list::ActionList)
    return length(list.actions)
end

function Base.iterate(list::ActionList, i=1)
    if i > length(list.actions)
        return nothing
    end
    return list.actions[i], i + 1
end

function action_names(list::ActionList)
    return [a.name for a in list.actions]
end


text:
#!/usr/bin/julia

function err()
   println("Usage: please input a non-negative integer")
end

function fac(n)
    if n < 0
        err()
        exit()
    end

    out = 1
    for i = 1:n
        out = out * i
    end
    return out
end

try
    println(fac(parse(Int, ARGS[1])))
catch e
    err()
end


text:
function drawcircle!(mat :: Matrix, xc, yc, r)

    x = 0;
    y = r;
    p = 1 - r;

     while x < y

         mat[round(xc+x), round(yc+y)] = 1
         mat[round(xc-x), round(yc+y)] = 1
         mat[round(xc+x), round(yc-y)] = 1
         mat[round(xc-x), round(yc-y)] = 1
         mat[round(xc+y), round(yc+x)] = 1
         mat[round(xc-y), round(yc+x)] = 1
         mat[round(xc+y), round(yc-x)] = 1
         mat[round(xc-y), round(yc-x)] = 1

         x += 1

         if p < 0
            p = p + 2 * x + 1
         else
            y = y - 1
            p = p + 2 * (x - y) + 1
         end

         mat[round(xc+x), round(yc+y)] = 1
         mat[round(xc-x), round(yc+y)] = 1
         mat[round(xc+x), round(yc-y)] = 1
         mat[round(xc-x), round(yc-y)] = 1
         mat[round(xc+y), round(yc+x)] = 1
         mat[round(xc-y), round(yc+x)] = 1
         mat[round(xc+y), round(yc-x)] = 1
         mat[round(xc-y), round(yc-x)] = 1
     end
end


text:
module ThreeBody

    ######
    using TaylorSeries, TaylorIntegration

    export M_1, M_2, μ, mu,
        L_1,L_2,L_3,L_x4,L_y4,L_x5,L_y5,
        L1,L2,L3,Lx4,Lx5,
        C_L1,C_L2,C_L3,C_L4,C_L5,
        mass_moon, mass_sun, mass_earth,
        _V_pseudo, _CR3BP_ode!, _C_Jacobi, _vy,
        collision_approach, JT_accuracy

    ######

    ### Earth, Moon, Sun masses ###
    const mass_moon  = 7.347673e22 #kg
    const mass_earth = 5.9722e24 #kg
    const mass_sun = 1.9855e30 #kg
    # mass_satelite ≈ 1e3 kg
    ### Earth, Moon, Sun distances ###
    const radius_earth_moon = 3.85e5 #km
    const radius_earth_sun  = 1.496e8 #km

    ### Mass parameters (with distance normalized to 1)
    M_1 = mass_earth
    M_2 = mass_moon
    mu(M_1,M_2) = M_2/(M_1 + M_2)
    μ = mu(M_1,M_2)

    ### Lagrangian Points ###
    # L1, L2, L3 can be obtained numerically more precisely
    """
    \$ 1 - ( μ /3 )^{1/3} \$
    """
    L1(μ)  =  1 - (μ/3)^(1/3)
    """
    \$ 1 + ( μ /3 )^{1/3} \$
    """
    L2(μ)  =  1 + (μ/3)^(1/3)
    """
    \$ - 1 - 5μ /12  \$
    """
    L3(μ)  =  -(1 + (5μ/12))
    """
    \$  (1 - 2μ)/2  \$
    """
    Lx4(μ) = 1/2 * (1 - 2μ)
    """
    \$  (1 - 2μ)/2  \$
    """
    Lx5(μ) = 1/2 * (1 - 2μ)

    L_1  = L1(μ)
    L_2  = L2(μ)
    L_3  = L3(μ)
    L_x4 = Lx4(μ)
    L_x5 = Lx5(μ)
    L_y4 = √3/2
    L_y5 = -√3/2

    ####  Potentials ###
    V_G(r, μ) = -(1-μ)/norm(r-[-μ,0.0]) - μ/norm(r-[1-μ,0.0])
    V_Ω(r) = -norm(r)^2/2

    """
    Returns the no coriolis potential for the CR3BP at value `r`.
    """
    _V_pseudo(r, μ) = V_G(r,μ) + V_Ω(r)

    # V_pseudo(r) = _V_pseudo(r,μ)

    ### Jacobi Constants ###
    C_L1 = _V_pseudo([L_1,0.0],μ)
    C_L2 = _V_pseudo([L_2,0.0],μ)
    C_L3 = _V_pseudo([L_3,0.0],μ)
    C_L4 = _V_pseudo([L_x4,L_y4],μ)
    C_L5 = _V_pseudo([L_x5,L_y5],μ);

    ### Eqs of motion ###
    """
    Equations of motion for CR3BP ready for integrating via `TaylorIntegration`
    """
    function _CR3BP_ode!(t, x, dx, μ)

        β = 1.0-μ
        r_13_cubed = ( (x[1]+μ)^2 + x[2]^2 )^(3/2)
        r_23_cubed = ( (x[1]-β)^2 + x[2]^2 )^(3/2)

        dx[1] = x[3]
        dx[2] = x[4]
        dx[3] = -( β*(x[1]+μ)/r_13_cubed + μ*(x[1]-β)/r_23_cubed ) + x[1] + 2x[4]
        dx[4] = -(β/r_13_cubed + μ/r_23_cubed - 1)*x[2] - 2x[3]

        nothing
    end

    # CR3BP_ode!(t, x, dx) = _CR3BP_ode!(t, x, dx, μ)

    """
    Distance between two 2-dimn solutions.
    """
    function JT_accuracy(ϕ,ϕ_JT)
        dims = size(ϕ)[2]
        diffs = zero(ϕ[:,1])
        for dim in 1:dims
            diffs += (ϕ[:,dim] - ϕ_JT[:,dim]).^2
        end
        return sqrt.(diffs)
    end

    """
    Returns the indexes for which two trajectories differ in less than `d_min`.
    """
    function collision_approach(ϕ1, ϕ2; d_min=0)
        d_ϕ1ϕ2 = JT_accuracy(ϕ1[:,1:2], ϕ2[:,1:2])
        if d_min == 0
            d_min = minimum(d_ϕ1ϕ2)
        end
        #approached_distance = Vector{eltype(d_A1_A2)}()
        #approached_time = Vector{eltype(tv)}()
        approached_index = Vector{Int64}()
        for i in eachindex(d_ϕ1ϕ2)
            if d_ϕ1ϕ2[i] <= d_min
                #push!(approached_distance, d_ϕ1ϕ2[i])
                #push!(approached_time, tv[i])
                push!(approached_index, i)
            end
        end

        return approached_index
    end

    ### Jacobi Constant ###
    """
    Computes Jacobi Energy \$ V(r) + v^2/2 \$.
    """
    _C_Jacobi(r,v,μ) = _V_pseudo(r,μ) + (v⋅v)/2.0

    ### implicit y-velocity ###
    """
    Obtains y-velocity for initial conditions to be at energy level `E`
    """
    _vy(r, vx, μ; E=C_L1) = √( 2*( E - _V_pseudo(r,μ) ) - vx^2 )

    println("M_1: $M_1 , M_2: $M_2, μ: $μ")
    println("L_1: $(round(L_1,3)), L_2: $(round(L_2,3)), L_3: $(round(L_1,3)), L_4x: $(round(L_x4,3)), L_4y: $(round(L_y4,3))")
    println("C_L1: $(round(C_L1,3))")
end

max_stars_repo_name: blas-ko/tesis

max_issues_repo_name: blas-ko/tesis

max_forks_repo_name: blas-ko/tesis


text:
# Helper funs
#
# (c) A. Dominik, 2020

# relu(), sigm(), used from Knet!

leaky_sigm(x; l=0.01) = Knet.sigm(x) + l*x
leaky_tanh(x; l=0.01) = tanh(x) + l*x
leaky_relu(x; l=0.01) = Knet.relu(x) + l*x


text:
module Coulomb
    using CxxWrap
    @wrapmodule(joinpath(@__DIR__, "../deps/libcoulomb.so"))
    function __init__()
        @initcxx
        initialize()
    end
end

import .Coulomb

Gaussian_size(::Gaussian_s) = 1
Gaussian_size(::Gaussian_p) = 3
Gaussian_size(s::STO3G) = Gaussian_size(s.g[1])

function make_shells(basis::Vector{STO3G})
    l = Vector{Cint}()
    alpha = Vector()
    d = Vector()
    R = Vector()
    i = 1
    Nshell = 0
    while i <= length(basis)
        sto = basis[i]
        g = sto.g[1]
        step = Gaussian_size(sto)
        push!(l, step == 1 ? 0 : 1)
        push!(alpha, [sto.g[1].α, sto.g[2].α, sto.g[3].α])
        push!(d, [sto.d...])
        push!(R, [g.R[1], g.R[2], g.R[3]])
        i += step
        Nshell += 1
    end
    return Nshell, l, vcat(alpha...), vcat(d...), vcat(R...)
end

function Ve!(ve::Array, basis::Vector{STO3G})
    Nshell, l, alpha, d, R = make_shells(basis)
    Coulomb.Ve(
        pointer(ve), Nshell, pointer(l), pointer(alpha), pointer(d), pointer(R)
    )
end

function Veff!(veff::Array, ρ::Array, Nshell::Int, l::Vector{Cint}, alpha::Vector{Float64}, d::Vector{Float64}, R::Vector{Float64})
    Coulomb.Veff(
        pointer(veff), pointer(ρ), Nshell, pointer(l), pointer(alpha), pointer(d), pointer(R)
    )
end


text:
# A pseudo-VGG19 model
# This model skips the max-pool layers and uses 2x2 strides to control image size
# Classification layers are left out.
# 16-layers (8 down, 8 up)
function pVGG_layers(image_size::Tuple{Int, Int})
    nlayers = 16
    kk = repeat([(3, 3)], nlayers);  # kernel size, all equal
    ss = repeat([(1, 1)], nlayers);  # stride size
    ss[1] = ss[3] = ss[5] = ss[9] = ss[13] = (2, 2)
    d = (1, 1);  # dilation size

    pc = padchain(nlayers, image_size, kk, ss)

    ae = Chain(
               # Encoder section
               # Block 1
               Conv(kk[1], 3=>64, relu, pad=pc[1], stride=ss[1], dilation=d),
               Conv(kk[2], 64=>64, relu, pad=pc[2], stride=ss[2], dilation=d),
               # Block 2
               Conv(kk[3], 64=>128, relu, pad=pc[3], stride=ss[3], dilation=d),
               Conv(kk[4], 128=>128, relu, pad=pc[4], stride=ss[4], dilation=d),
               # Block 3
               Conv(kk[5], 128=>256, relu, pad=pc[5], stride=ss[5], dilation=d),
               Conv(kk[6], 256=>256, relu, pad=pc[6], stride=ss[6], dilation=d),
               Conv(kk[7], 256=>256, relu, pad=pc[7], stride=ss[7], dilation=d),
               Conv(kk[8], 256=>256, relu, pad=pc[8], stride=ss[8], dilation=d),
               # Block 4
               Conv(kk[9], 256=>512, relu, pad=pc[9], stride=ss[9], dilation=d),
               Conv(kk[10], 512=>512, relu, pad=pc[10], stride=ss[10], dilation=d),
               Conv(kk[11], 512=>512, relu, pad=pc[11], stride=ss[11], dilation=d),
               Conv(kk[12], 512=>512, relu, pad=pc[12], stride=ss[12], dilation=d),
               # Block 5
               Conv(kk[13], 512=>512, relu, pad=pc[13], stride=ss[13], dilation=d),
               Conv(kk[14], 512=>512, relu, pad=pc[14], stride=ss[14], dilation=d),
               Conv(kk[15], 512=>512, relu, pad=pc[15], stride=ss[15], dilation=d),
               Conv(kk[16], 512=>512, relu, pad=pc[16], stride=ss[16], dilation=d),
               # Decoder section
               # Block 5
               ConvTranspose(kk[16], 512=>512, relu, pad=pc[16], stride=ss[16], dilation=d),
               ConvTranspose(kk[15], 512=>512, relu, pad=pc[15], stride=ss[15], dilation=d),
               ConvTranspose(kk[14], 512=>512, relu, pad=pc[14], stride=ss[14], dilation=d),
               ConvTranspose(kk[13], 512=>512, relu, pad=pc[13], stride=ss[13], dilation=d),
               # Block 4
               ConvTranspose(kk[12], 512=>512, relu, pad=pc[12], stride=ss[12], dilation=d),
               ConvTranspose(kk[11], 512=>512, relu, pad=pc[11], stride=ss[11], dilation=d),
               ConvTranspose(kk[10], 512=>512, relu, pad=pc[10], stride=ss[10], dilation=d),
               ConvTranspose(kk[9], 512=>256, relu, pad=pc[9], stride=ss[9], dilation=d),
               # Block 3
               ConvTranspose(kk[8], 256=>256, relu, pad=pc[8], stride=ss[8], dilation=d),
               ConvTranspose(kk[7], 256=>256, relu, pad=pc[7], stride=ss[7], dilation=d),
               ConvTranspose(kk[6], 256=>256, relu, pad=pc[6], stride=ss[6], dilation=d),
               ConvTranspose(kk[5], 256=>128, relu, pad=pc[5], stride=ss[5], dilation=d),
               # Block 2
               ConvTranspose(kk[4], 128=>128, relu, pad=pc[4], stride=ss[4], dilation=d),
               ConvTranspose(kk[3], 128=>64, relu, pad=pc[3], stride=ss[3], dilation=d),
               # Block 1
               ConvTranspose(kk[2], 64=>64, relu, pad=pc[2], stride=ss[2], dilation=d),
               ConvTranspose(kk[1], 64=>3, relu, pad=pc[1], stride=ss[1], dilation=d)
               );

    return ae
end

struct pVGG <: AutoEncoderModel
    layers::Chain
end

pVGG(image_size::Tuple{Int, Int}) = pVGG(pVGG_layers(image_size))

Base.show(io::IO, ::pVGG) = print(io, "pVGG()")

Flux.@functor pVGG

(m::pVGG)(x) = m.layers(x)


text:
include("interps.jl")

f(x) = 1 / (25 * x * x + 1)

function myerror(w)
  N = 1000
  res = abs(f(-1.) - w(-1.))
  for i=1:N
    temp = abs(f(-1. + (2*i/N)) - w(-1. + (2*i/N)))
    if temp > res
      res = temp
    end
  end
  return res
end

function sredniokwadratowa(w)
  N = 1000
  res = abs(f(-1.) - w(-1.))
  for i=1:N
    temp = abs(f(-1. + (2*i/N)) - w(-1. + (2*i/N)))
    if temp > res
      res = temp
    end
  end
  return res
end

# Rownoodlegle 10
xss = linspace(-1,1,10)
w = interpNewton(xss, map(f,xss))
myerror(w)

#zera czebyszewa
xss = czebyszew_zeroes(10)
w = interpNewton(xss, map(f,xss))
myerror(w)

#punkty ekstyrmalne czebyszewa
xss = czebyszew_extreme(9)
w = interpNewton(xss, map(f,xss))
myerror(w)

#lagrande
w = x -> lagrandePolymonial(9,x)
w(1)
myerror(w)

#czebyszew optymalny
w = x -> czebyszewPolymonial(9,x)
w(0.15)
myerror(w)

``

max_stars_repo_name: Magikis/Uniwersity

max_issues_repo_name: Magikis/Uniwersity

max_forks_repo_name: Magikis/Uniwersity


text:
module LaTeXFormatTests

using Test

using Documenter

# Documenter package docs
info("Building Documenter's docs with LaTeX.")
const Documenter_root = normpath(joinpath(dirname(@__FILE__), "..", "..", "docs"))
doc = makedocs(
    debug = true,
    root = Documenter_root,
    modules = [Documenter],
    clean = false,
    format = :latex,
    sitename = "Documenter.jl",
    authors = "Michael Hatherly, Morten Piibeleht, and contributors.",
    pages = [
        "Home" => "index.md",
        "Manual" => Any[
            "Guide" => "man/guide.md",
            "man/examples.md",
            "man/syntax.md",
            "man/doctests.md",
            "man/hosting.md",
            "man/latex.md",
            "man/internals.md",
            "man/contributing.md",
        ],
        "Library" => Any[
            "Public" => "lib/public.md",
            "Internals" => Any[
                "Internals" => "lib/internals.md",
                "lib/internals/anchors.md",
                "lib/internals/builder.md",
                "lib/internals/cross-references.md",
                "lib/internals/docchecks.md",
                "lib/internals/docsystem.md",
                "lib/internals/documents.md",
                "lib/internals/dom.md",
                "lib/internals/expanders.md",
                "lib/internals/formats.md",
                "lib/internals/generator.md",
                "lib/internals/mdflatten.md",
                "lib/internals/selectors.md",
                "lib/internals/utilities.md",
                "lib/internals/walkers.md",
                "lib/internals/writers.md",
            ]
        ]
    ]
)

@testset "LaTeX" begin
    @test isa(doc, Documenter.Documents.Document)
end

end


text:
"""
Generates a random value from a truncated normal distribution that represents
a new item. It's used to simulate incoming application to the registry.
"""
function getRegistryCandidate()
    rand(TruncatedNormal(50, 20, 0, 100))
end

"""
Calculates the mean quality of the items in the given registry
"""
function registryMean(registry)
    length(registry) == 0 ? 0 : mean(registry)
end

max_stars_repo_name: mlodykaz/tcr-simulation

max_issues_repo_name: mlodykaz/tcr-simulation

max_forks_repo_name: mlodykaz/tcr-simulation


text:
# ======================================
# Solve function for continuous systems
# ======================================

"""
    solve(ivp::IVP{<:AbstractContinuousSystem}, tspan, alg; kwargs...)

Solves the initial-value problem defined by `ivp` over the time span `tspan`,
using the algorihm `alg`. If no algorithm is given, a default algorithm is chosen.

### Input

- `ivp`   -- initial-value problem
- `tspan` -- time span for this initial-value problem
- `alg`   -- reachability algorithm

Additional options are passed as arguments or keyword arguments; see the notes
below for details. See the online documentation for examples.

### Output

The solution of a reachability problem, as an instance of a `ReachSolution`.

### Notes

- Use the `alg`, `algorithm` or `opC` keyword arguments to specify the algorithm
  to solve the initial-value problem. Algorithm-specific options should be passed
  to the algorithm constructor as well.

- Use the `tspan` keyword argument to specify the time span; it can be:
    - a tuple,
    - an interval, or
    - a vector with two components.

- Use the `T` keyword argument to specify the time horizon; the initial time is
  then assumed to be zero.

- Use the `static` keyword argument to force conversion to static arrays in the
  algorithm (should be supported by the algorithm).

- Use the `NSTEPS` keyword argument to specify the number of discrete steps
  solved in the set-based recurrence.

- Use the `threading` option to use multi-threading parallelism. This option applies
  for initial-value problems whose initial condition is a vector of sets.
"""
function solve(ivp::IVP{<:AbstractContinuousSystem}, args...; kwargs...)

    # preliminary checks
    _check_dim(ivp)

    # get time span (or the emptyset if NSTEPS was specified)
    tspan = _get_tspan(args...; kwargs...)

    # get the continuous post or find a default one
    cpost = _get_cpost(ivp, args...; kwargs...)
    if cpost == nothing
        cpost = _default_cpost(ivp, tspan; kwargs...)
    end

    # run the continuous-post operator
    F = post(cpost, ivp, tspan; kwargs...)

    got_ensemble = get(kwargs, :ensemble, false)
    if got_ensemble
        @requires DifferentialEquations
        # compute trajectories using ensemble simulation
        ensemble_sol = _solve_ensemble(ivp, args...; kwargs...)
        dict = Dict{Symbol,Any}(:ensemble=>ensemble_sol)
        sol = ReachSolution(F, cpost, dict)
    else
        # wrap the flowpipe and algorithm in a solution structure
        sol = ReachSolution(F, cpost)
    end

    return sol
end

# solve for distributed initial conditions; uses multi-threaded implementation by default
function solve(ivp::IVP{AT, VT}, args...; kwargs...) where {AT<:AbstractContinuousSystem,
                                                            ST<:Union{LazySet, IntervalBox}, VT<:AbstractVector{ST}}
    _check_dim(ivp)
    tspan = _get_tspan(args...; kwargs...)
    cpost = _get_cpost(ivp, args...; kwargs...)
    if cpost == nothing
        cpost = _default_cpost(ivp, tspan; kwargs...)
    end

    X0 = initial_state(ivp)
    S = system(ivp)

    threading = haskey(kwargs, :threading) ? kwargs[:threading] : (Threads.nthreads() > 1)

    F = _solve_distributed(cpost, S, X0, tspan, Val(threading); kwargs...)
    return ReachSolution(MixedFlowpipe(F), cpost)
end

function _solve_distributed(cpost, S, X0, tspan, threading::Val{false}; kwargs...)
    return [post(cpost, IVP(S, X0i), tspan; kwargs...) for X0i in X0]
end

function _solve_distributed(cpost, S, X0, tspan, threading::Val{true}; kwargs...)
    nsets = length(X0)
    FT = Flowpipe{numtype(cpost), rsetrep(cpost)} # TODO add third parameter
    sol_tot = Vector{FT}(undef, nsets)

    Threads.@threads for i in 1:length(X0)
        sol_i = ReachabilityAnalysis.post(cpost, IVP(S, X0[i]), tspan; kwargs...)
        sol_tot[i] = sol_i
    end
    return sol_tot
end

# ====================================
# Continuous post operator interface
# ====================================

"""
    AbstractPost

Abstract supertype of all post operator types.
"""
abstract type AbstractPost end

"""
    AbstractContinuousPost

Abstract supertype of all continuous post operators.
"""
abstract type AbstractContinuousPost <: AbstractPost end

setrep(::InitialValueProblem{HS, ST}) where {HS, ST} = ST

# ==================
# Argument handling
# ==================

# extend dimension for common IntervalArithmetic types
LazySets.dim(::IA.Interval) = 1
LazySets.dim(::IA.IntervalBox{D, N}) where {D, N} = D

# lazy sets (or sets that behave like such -- TODO update once UnionSet <: LazySet)
_dim(X::Union{<:LazySet, <:IA.Interval, <:IA.IntervalBox, <:UnionSet, <:UnionSetArray}) = dim(X)

# singleton elements
_dim(X::Number) = 1
_dim(X::AbstractVector{N}) where {N<:Number} = length(X)

# vector of sets
function  _dim(X::AbstractVector{UT}) where {UT<:Union{<:LazySet,
                <:IA.Interval, <:IA.IntervalBox, <:UnionSet, <:UnionSetArray}}
    n = _dim(first(X))
    all(X -> _dim(X) == n, X)  || throw(ArgumentError("dimension mismatch between " *
            "the initial sets in this array; expected only sets of dimension $n"))
    return n
end

_dim(X) = throw(ArgumentError("the type of the initial condition, $(typeof(X)), cannot be handled"))

function _check_dim(S, X0; throw_error::Bool=true)
    n = statedim(S)
    d = _dim(X0)
    n == d && return true

    if throw_error
        throw(ArgumentError("the state-space dimension should match " *
                            "the dimension of the initial state, but they are " *
                            "$n and $d respectively"))
    end
    return false
end

_dim(X0::Tuple{VT, VT}) where {N, VT<:AbstractVector{N}} = length(X0[1]) + length(X0[2])
_dim(X0::Tuple{<:LazySet{N}, <:LazySet{N}}) where {N} = dim(X0[1]) + dim(X0[2])

function _check_dim(S::Union{SecondOrderLinearContinuousSystem,
                             SecondOrderAffineContinuousSystem,
                             SecondOrderConstrainedLinearControlContinuousSystem,
                             SecondOrderConstrainedAffineControlContinuousSystem},
                    X0; throw_error::Bool=true)
    n = statedim(S)
    d = _dim(X0)
    d == 2*n && return true

    if throw_error
        throw(ArgumentError("the dimension of the initial state should be " *
                            "twice the state-space dimension, but they are " *
                            "$d and $n respectively"))
    end
    return false
end

function _check_dim(ivp::InitialValueProblem{<:AbstractContinuousSystem}; throw_error::Bool=true)
    S = system(ivp)
    X0 = initial_state(ivp)
    _check_dim(S, X0, throw_error=throw_error)
end

# promotion from tuple-like arguments
@inline _promote_tspan((t1, t2)::Tuple{T, T}) where {T} = TimeInterval(t1, t2)
@inline _promote_tspan((t1, t2)::Tuple{T, S}) where {T, S} = TimeInterval(promote(t1, t2))

# no-op, corresponds to (inf(tspan), sup(tspan))
@inline _promote_tspan(tspan::IA.Interval) = tspan

# no-op, takes interval wrapped data; corresponds to (min(tspan), max(tspan))
@inline _promote_tspan(tspan::Interval) = tspan.dat

# number T defaults to a time interval of the form [0, T]
@inline _promote_tspan(tspan::Number) = TimeInterval(zero(tspan), tspan)

# catch-all array like tspan
@inline function _promote_tspan(tspan::AbstractArray)
    @assert length(tspan) == 2 throw(ArgumentError("the length of `tspan` must be two, but " *
                                                   "it is of length $(length(tspan))"))
    return TimeInterval(tspan[1], tspan[2])
end

function _get_tspan(args...; kwargs...)
    got_tspan = haskey(kwargs, :tspan)
    got_T = haskey(kwargs, :T)
    got_NSTEPS = haskey(kwargs, :NSTEPS)

    # throw error for repeated arguments
    if got_tspan && got_T
        throw(ArgumentError("cannot parse the time horizon `T` and the " *
            "time span `tspan` simultaneously; use one or the other"))
    elseif got_tspan && got_NSTEPS
        throw(ArgumentError("cannot parse the time span `tspan` and the " *
            "number of steps `NSTEPS` simultaneously; use one or the other"))
    elseif got_T && got_NSTEPS
        throw(ArgumentError("cannot parse the time horizon `T` and the " *
            "number of steps `NSTEPS` simultaneously; use one or the other"))
    end

    # parse time span
    if got_tspan
        tspan =_promote_tspan(kwargs[:tspan])

    elseif got_T
        tspan = _promote_tspan(kwargs[:T])

    elseif got_NSTEPS
        tspan = emptyinterval() # defined a posteriori

    elseif !isempty(args) && applicable(_promote_tspan, args[1])
        # applies to tuples, vectors, LazySets.Interval, IA.Interval, and numbers
        tspan = _promote_tspan(args[1])
#=
        (args[1] isa Tuple{Float64, Float64} || # time span given as tuple
                              args[1] isa Vector{Float64}         || # time span given as vector
                              args[1] isa Interval ||  # time span given as LazySets.Interval
                              args[1] isa IntervalArithmetic.Interval || # time span given as IA.Interval
                              args[1] isa Real) # got time horizon as first argument
=#

    elseif  haskey(kwargs, :max_jumps)
        # got maximum number of jumps, applicable to hybrid systems
        tspan = emptyinterval()

    else

        # couldn't find time span => error
        throw(ArgumentError("the time span has not been specified, but is required " *
                            "for `solve`; you should specify either the time horizon " *
                            "`T=...`, the time span `tspan=...`, or the number of steps, `NSTEPS=...`"))
    end
    return tspan
end

# return the time horizon given a time span
# the check_positive flag is used for algorithms that do not support negative
# times
function _get_T(tspan::TimeInterval; check_zero::Bool=true, check_positive::Bool=true)
    t0 = inf(tspan)
    if check_zero
        @assert iszero(t0) "this algorithm can only handle zero initial time"
    end
    T = sup(tspan)
    if check_positive
        @assert T > 0 "the time horizon should be positive"
    end
    return T
end

tstart(Δt::TimeInterval) = inf(Δt)
tend(Δt::TimeInterval) = sup(Δt)
tspan(Δt::TimeInterval) = Δt

function _get_cpost(ivp, args...; kwargs...)
    got_alg = haskey(kwargs, :alg)
    got_algorithm = haskey(kwargs, :algorithm)
    got_opC = haskey(kwargs, :opC)
    no_args = isempty(args) || args[1] === nothing

    # continous post was specified
    if got_alg
        cpost = kwargs[:alg]
    elseif got_algorithm
        cpost = kwargs[:algorithm]
    elseif got_opC
        cpost = kwargs[:opC]

    # check if either args[1] or args[2] are the post
    elseif !no_args && typeof(args[1]) <: AbstractContinuousPost
        cpost = args[1]
    elseif !no_args && length(args) == 2 && typeof(args[2]) <: AbstractContinuousPost
        cpost = args[2]
    elseif !no_args && length(args) > 2
        throw(ArgumentError("the number of arguments, $(length(args)), is not valid"))

    # no algorithm found => compute default
    else
        cpost = nothing
    end
    return cpost
end

# ===================
# Algorithm defaults
# ===================

# number of discrete steps used if only the time horizonis given
const DEFAULT_NSTEPS = 100

function _default_cpost(ivp::AbstractContinuousSystem, tspan; kwargs...)
    if isaffine(ivp)
        if haskey(kwargs, :δ)
            δ = kwargs[:δ]
        elseif haskey(kwargs, :N)
            N = kwargs[:N]
            δ = diam(tspan) / N
        elseif haskey(kwargs, :NSTEPS)
            NSTEPS = kwargs[:NSTEPS]
            δ = diam(tspan) / NSTEPS
        elseif haskey(kwargs, :num_steps)
            num_steps = kwargs[:num_steps]
            δ = diam(tspan) / num_steps
        else
            δ = diam(tspan) / DEFAULT_NSTEPS
        end
        if statedim(ivp) == 1 && !is_second_order(ivp)
            opC = INT(δ=δ)
        else
            static = haskey(kwargs, :static) ? kwargs[:static] : false
            opC = GLGM06(δ=δ, static=static)
        end
    else
        # check additional kwargs options if they exist, allowing some aliases
        orderQ = get(kwargs, :orderQ, DEFAULT_ORDER_Q_TMJETS)
        orderT = get(kwargs, :orderT, DEFAULT_ORDER_T_TMJETS)

        if haskey(kwargs, :max_steps)
            maxsteps = kwargs[:max_steps]
        elseif haskey(kwargs, :maxsteps)
            maxsteps = kwargs[:maxsteps]
        else
            maxsteps = DEFAULT_MAX_STEPS_TMJETS
        end

        if haskey(kwargs, :abs_tol)
            abstol = kwargs[:abs_tol]
        elseif haskey(kwargs, :abstol)
            abstol = kwargs[:abstol]
        else
            abstol = DEFAULT_ABS_TOL_TMJETS
        end

        opC = TMJets(orderQ=orderQ, orderT=orderT, maxsteps=maxsteps, abstol=abstol)
    end
    return opC
end

"""
    _default_cpost(ivp::IVP{<:AbstractContinuousSystem}, tspan; kwargs...)

### Input

- `ivp`   -- initial-value problem
- `tspan` -- time-span

### Output

A continuous post-operator that can handle the given initial-value problem.

### Notes

If the system is affine, then:

- If it is one-dimensional, algorithm `INT` is used, otherwise,
- Algorithm `GLGM06` is used.

If the system is not affine, then the algorithm `TMJets` is used.
"""
function _default_cpost(ivp::IVP{<:AbstractContinuousSystem}, tspan; kwargs...)
    _default_cpost(system(ivp), tspan; kwargs...)
end

# return a default algorithm, assuming that the first mode is representative
function _default_cpost(ivp::IVP{<:AbstractHybridSystem}, tspan; kwargs...)
    first_mode = mode(system(ivp), 1)
    _default_cpost(first_mode, tspan; kwargs...)
end


text:
using Test

using Cbc
using JuMP

using MIPVerify: set_log_level!
using MIPVerify: get_max_index, get_norm, get_default_tightening_options

@isdefined(TestHelpers) || include("TestHelpers.jl")

@testset "MIPVerify" begin
    TestHelpers.reset_timer!()
    set_log_level!("warn")

    include("integration.jl")
    include("net_components.jl")
    include("utils.jl")
    include("models.jl")
    include("batch_processing_helpers.jl")

    @testset "get_default_tightening_options" begin
        @test get_default_tightening_options(() -> Cbc.Optimizer()) ==
              Dict("logLevel" => 0, "seconds" => 20)
        @test get_default_tightening_options(Cbc.Optimizer) ==
              Dict("logLevel" => 0, "seconds" => 20)
        @test get_default_tightening_options(() -> nothing) == Dict()
    end

    @testset "get_max_index" begin
        @test_throws MethodError get_max_index([])
        @test get_max_index([3]) == 1
        @test get_max_index([3, 1, 4]) == 3
        @test get_max_index([3, 1, 4, 1, 5, 9, 2]) == 6
    end

    @testset "get_norm" begin
        @testset "real-valued arrays" begin
            xs = [1, -2, 3]
            @test get_norm(1, xs) == 6
            @test get_norm(2, xs) == sqrt(14)
            @test get_norm(Inf, xs) == 3
            @test_throws DomainError get_norm(3, xs)
        end
        @testset "variable-valued arrays" begin
            @testset "l1" begin
                m = TestHelpers.get_new_model()
                x1 = @variable(m, lower_bound = 1, upper_bound = 5)
                x2 = @variable(m, lower_bound = -8, upper_bound = -2)
                x3 = @variable(m, lower_bound = 3, upper_bound = 10)
                xs = [x1, x2, x3]
                n_1 = get_norm(1, xs)
                n_2 = get_norm(2, xs)
                n_inf = get_norm(Inf, xs)

                @objective(m, Min, n_1)
                optimize!(m)
                @test JuMP.objective_value(m) ≈ 6

                if Base.find_package("Gurobi") !== nothing
                    # Skip these tests if Gurobi is not installed.
                    # Cbc does not solve problems with quadratic objectives
                    @objective(m, Min, n_2)
                    optimize!(m)
                    @test JuMP.objective_value(m) ≈ 14
                end

                @objective(m, Min, n_inf)
                optimize!(m)
                @test JuMP.objective_value(m) ≈ 3

                @test_throws DomainError get_norm(3, xs)
            end
        end
    end

    println()
    TestHelpers.print_timer()
    println()
    println()
end


text:
module Evolutionary
    using Random, LinearAlgebra, Statistics
    using Base: @kwdef
    using UnPack: @unpack
    using NLSolversBase: AbstractObjective, NonDifferentiable, ConstraintBounds,
                         value, value!, nconstraints_x, nconstraints, AbstractConstraints

    import NLSolversBase: NonDifferentiable, f_calls, value, value!
    import Base: show, copy, minimum, summary, identity, getproperty, rand, getindex, length, copyto!, setindex!

    export AbstractStrategy, strategy, mutationwrapper,
           IsotropicStrategy, AnisotropicStrategy, NoStrategy,
           isfeasible, BoxConstraints, apply!, penalty, penalty!,
           PenaltyConstraints, WorstFitnessConstraints, MixedTypePenaltyConstraints,
           # ES mutations
           gaussian, cauchy,
           # GA mutations
           flip, bitinversion, domainrange, inversion, insertion, swap2, scramble, shifting, PM, MIPM,
           # ES recombinations
           average, marriage,
           # GA recombinations
           singlepoint, twopoint, uniform,
           discrete, waverage, intermediate, line, HX, LX, MILX,
           PMX, OX1, CX, OX2, POS,
           # GA selections
           ranklinear, uniformranking, roulette, rouletteinv, sus, susinv, tournament, truncation,
           # DE selections
           random, permutation, randomoffset, best,
           # DE recombinations
           uniformbin, exponential,
           # GP exports
           Terminal, crosstree, mutatetree,
           # Optimization methods
           ES, CMAES, GA, DE, TreeGP,
           # re-export
           NonDifferentiable, value, value!

    # optimize API
    include("api/types.jl")
    include("api/results.jl")
    include("api/utilities.jl")
    include("api/constraints.jl")
    include("api/optimize.jl")

    # ES & GA recombination functions
    include("recombinations.jl")

    # ES & GA mutation functions
    include("mutations.jl")

    # GA selection functions
    include("selections.jl")

    # Evolution Strategy
    include("es.jl")
    include("cmaes.jl")

    # Genetic Algorithms
    include("ga.jl")

    # Differential Evolution
    include("de.jl")

    # Genetic Programming
    include("gp.jl")

    # deprecations
    @deprecate isotropic(recombinant, strategy) gaussian(recombinant, strategy)
    @deprecate anisotropic(recombinant, strategy) gaussian(recombinant, strategy)
    @deprecate isotropicSigma(strategy) gaussian(strategy)
    @deprecate anisotropicSigma(strategy) gaussian(strategy)
    @deprecate averageSigma(strategy) average(strategy)
    @deprecate strategy(kwargs...) IsotropicStrategy(N) false
end


text:
using Millboard
using Base.Test

board = Dict()
@test """
+==+""" == table(board) |> string

board = Dict(:a=>:b, "a"=>"b")
@test """
+=====+=======+
| key | value |
+=====+=======+
|  :a |    :b |
+-----+-------+
|   a |     b |
+-----+-------+""" == table(board) |> string

board = Dict("1x3"=>[1 2 3], "2x3"=>[1 2 3; 4 5 6], "3x1"=> [1; 2; 3])
@test """
+=====+=======+
| key | value |
+=====+=======+
| 1x3 | 1 2 3 |
+-----+-------+
|     | 1 2 3 |
| 2x3 | 4 5 6 |
+-----+-------+
|     |     1 |
|     |     2 |
| 3x1 |     3 |
+-----+-------+""" == table(board) |> string


text:
using NeidArchive
using Test
using Dates

@testset "NeidArchive.jl" begin
    user_nexsci = "pyneidprop"
    passwd_nexsci = "pielemonquietyellow"
    cookie_nexsci = "./neidadmincookie.txt"
    query_result_file = "./criteria.csv"

    @testset "login" begin
        NeidArchive.login(userid=user_nexsci, password=passwd_nexsci, cookiepath=cookie_nexsci)
    end

    @testset "Query" begin
        param = Dict{String,String}()
        param["datalevel"] = "l0"
        param["piname"] = "Logsdon"
        param["datetime"] = NeidArchive.datetime_one_day(Date(2021,1,31))
        NeidArchive.query(param, cookiepath=cookie_nexsci, outpath=query_result_file)
        num_lines = countlines(query_result_file) - 1
        #println("# Query resulted in file with ", num_lines, " entries.")
        #
        #df = CSV.read(query_result_file, DataFrame)
        #num_files = size(df, 1)
        @test num_lines == 238
    end

end


text:
function pilot(length, freq; samplerate=1)
    L = round(Int, uconvert(NoUnits, length*samplerate))
    dθ = uconvert(NoUnits, 2π*freq/samplerate)

    cos.((0:L-1) * dθ)
end

function pilotfilt(sig, ω)
    # DSP.jl normalizes to half-cycle/sample, so that 1 is the nyquist,
    # but the rest of this code uses 1 cycle/sample (so 0.5 is nyquist)
    bpf = digitalfilter(Bandpass(2ω*0.97, 2ω/0.97),
                        FIRWindow(;transitionwidth=2ω*0.1))
    filt(bpf, sig)[length(bpf):end]
end

"""
Find a pilot tone in the given signal.
"""
function findpilot(sig, freq, len; samplerate=1)
    # TODO: this isn't robust to a wideband energy increase. It looks for an
    # energy pulse in the target band, but that could also get triggered if
    # everything just gets louder
    ω = freqnorm(freq, samplerate)
    L = durnorm(len, samplerate)
    filtered = pilotfilt(sig, ω)

    env = abs2.(hilbert(filtered))
    _, onset = findmax(xcorr(env, ones(L); padmode=:none))
    onset -= L

    prerollpower = sum(x->x^2, filtered[1:onset-1]) / (onset-1)
    pilotpower = sum(x->x^2, filtered[(0:L-1).+onset]) / L

    snr = pow2db(pilotpower/prerollpower)
    if snr < 6
        @warn @sprintf "Narrowband Pilot SNR only %.2f dB" snr
    end
    onset
end

"""
Compute the average phase difference of `sig` (which should correspond to the
frequency). The frequency is given in cycles/sample
"""
function avgphasediff(sig)
    # we remove the ends to get rid of the ringing from the hilbert transform
    an = hilbert(sig)[2000:end-2000]
    diff = @view(an[2:end]) ./ @view an[1:end-1]
    mean(angle, diff)/2π
end

"""
Takes the given signal and pilot frequency and returns the resample factor
necessary to synchronize the signal.
"""
function pilotsync(sig, freq; samplerate=1)
    ω = freqnorm(freq, samplerate)
    # bpf = gaussian(N, 0.20) .* cos.((0:N-1) * ω)
    # filtered = filt(bpf, sig)[length(bpf):end]
    filtered = pilotfilt(sig, ω)
    ω̂₁ = avgphasediff(filtered)
    @info "frequency estimated from hilbert transform: $(ω̂₁)"

    # estimate amplitude, frequency, and phase
    model(n, p) = @. p[1] * cos(2π*p[2]*n + p[3])

    fit = curve_fit(model, (0:length(filtered)-1), filtered,
                    [1.0, ω̂₁, 0.0], autodiff=:forwarddiff)

    # @info "frequency estimated from lsq fit: $(fit.param[2]/2π*samplerate)"
    @info "fit amplitude: $(fit.param[1])"
    @info "fit frequency: $(fit.param[2] * samplerate)"
    @info "fit phase (rad): $(fit.param[3])"

    fit.param[2] / ω
end

"""
    phaselock(x; alpha=0.05)

Simple PLL (phase-locked loop). Assumes `x` to be a real or complex sinusoidal
signal. Returns a tuple of (freq, phase, err, x2), where each is a vector with
the same length as `x`. `x2` will be real or complex to match `x`.

`alpha` controls how agressive the PLL is in correcting its frequency and phase
estimates. Lower values are more stable but converge more slowly.

## Example

```
using Plots

in_p = 0:0.1:500
x = cos.(in_p) + randn(length(in_p))*0.3

f, p, err, x2 = phaselock(x)
```

Algorithm ported from C code at http://liquidsdr.org/blog/pll-simple-howto/
Released under the MIT License, Copyright (c) 2007 - 2016 Joseph Gaeddert
"""
function phaselock(x::Vector{<:Complex}; alpha=0.05)
    beta = 0.5*alpha^2
    freq = 0.0;
    phase = 0.0;
    freqs = zeros(length(x))
    phases = zeros(length(x))
    errs = zeros(length(x))
    out = zeros(Complex{Float64}, length(x))

    for i in eachindex(x)
        out[i] = exp(im*phase)
        err = angle(x[i]*conj(out[i]))

        freq += beta * err
        phase += freq + alpha * err

        freqs[i] = freq/2π
        phases[i] = phase/2π
        errs[i] = err
    end

    (freqs, phases, errs, out)
end

function phaselock(x::Vector{<:Real}; kwargs...)
    freqs, phases, errs, out = phaselock(hilbert(x); kwargs...)

    (freqs, phases, errs, real.(out))
end


text:
abstract type SPIRType end

struct VoidType <: SPIRType end

abstract type ScalarType <: SPIRType end

struct IntegerType <: ScalarType
    width::Int
    signed::Bool
end

function IntegerType(inst::Instruction)
    (width, signed) = inst.arguments
    IntegerType(width, signed)
end

struct FloatType <: ScalarType
    width::Int
end

FloatType(inst::Instruction) = FloatType(inst.arguments[])

struct BooleanType <: ScalarType end

struct VectorType <: SPIRType
    eltype::ScalarType
    n::Int
end

function VectorType(inst::Instruction, eltype::ScalarType)
    VectorType(eltype, last(inst.arguments))
end

struct MatrixType <: SPIRType
    eltype::VectorType
    n::Int
end

function MatrixType(inst::Instruction, eltype::VectorType)
    MatrixType(eltype, last(inst.arguments))
end

struct ImageType <: SPIRType
    sampled_type::SPIRType
    dim::Dim
    format::ImageFormat
    arrayed::Int
    multisampled::Bool
    depth::Optional{Bool}
    access_qualifier::Optional{AccessQualifier}
    sampled::Optional{Bool}
    function ImageType(sampled_type::SPIRType, dim::Dim, format::ImageFormat, arrayed::Int, multisampled::Bool, depth::Optional{Bool}, sampled::Optional{Bool}, access_qualifier::Optional{AccessQualifier})
        if dim == DimSubpassData
            !sampled || error("Subpass image cannot be sampled")
            format == ImageFormatUnknown || error("Subpass image format must be unknown")
        end
        new(sampled_type, dim, format, arrayed, multisampled, depth, access_qualifier)
    end
    function ImageType(sampled_type, dim, format, arrayed, multisampled, depth, sampled, access_qualifier)
        ImageType(convert(SPIRType, sampled_type), convert(Dim, dim), convert(ImageFormat, format), convert(Int, arrayed), convert(Bool, multisampled), convert(Optional{Bool}, depth), convert(Optional{Bool}, sampled), convert(Optional{AccessQualifier}, access_qualifier))
    end
end

function ImageType(inst::Instruction, sampled_type::SPIRType)
    @assert inst.opcode == OpTypeImage
    (sampled_type_id, dim, depth, arrayed, multisampled, sampled, format) = inst.arguments[1:7]
    access_qualifier = length(inst.arguments) > 7 ? last(inst.arguments) : nothing
    sampled = sampled == 0 ? nothing : Bool(sampled - 1)
    if depth == 2
        depth = nothing
    end
    ImageType(sampled_type, dim, format, arrayed, multisampled, depth, sampled, access_qualifier)
end

struct SamplerType <: SPIRType end

struct SampledImageType <: SPIRType
    image_type::ImageType
end

struct ArrayType <: SPIRType
    eltype::SPIRType
    "Constant expression giving the size of the array, if determined at compile-time. SPIR-V runtime arrays will have this field set to `nothing`."
    size::Optional{Instruction}
end

struct OpaqueType <: SPIRType
    name::Symbol
end

struct StructType <: SPIRType
    members::Vector{SPIRType}
end

struct PointerType <: SPIRType
    storage_class::StorageClass
    type::SPIRType
end

PointerType(inst::Instruction, type::SPIRType) = PointerType(first(inst.arguments), type)

struct Constant
    value::Any
    is_spec_const::Bool
end

struct FunctionType <: SPIRType
    rettype::ID
    argtypes::Vector{ID}
end

function parse_type(inst::Instruction, types::SSADict{SPIRType}, results::SSADict{Any})::SPIRType
    @match op = inst.opcode begin
        &OpTypeVoid => VoidType()
        &OpTypeInt => IntegerType(inst)
        &OpTypeFloat => FloatType(inst)
        &OpTypeBool => BooleanType()
        &OpTypeVector => VectorType(inst, types[first(inst.arguments)])
        &OpTypeMatrix => MatrixType(inst, types[first(inst.arguments)])
        &OpTypeImage => ImageType(inst, types[first(inst.arguments)])
        &OpTypeSampler => SamplerType()
        &OpTypeSampledImage => SampledImageType(types[inst.arguments[]])
        &OpTypeArray => ArrayType(types[first(inst.arguments)], results[last(inst.arguments)])
        &OpTypeRuntimeArray => ArrayType(types[inst.arguments[]], nothing)
        &OpTypeStruct => StructType([types[id] for id in inst.arguments[]])
        &OpTypeOpaque => OpaqueType(Symbol(inst.arguments[]))
        &OpTypePointer => PointerType(inst, types[last(inst.arguments)])
        _ => error("$op does not represent a SPIR-V type or the corresponding type is not implemented.")
    end
end