text:
"""
    CameraState(;
        position::VecE2      = VecE2(0.,0.),
        zoom::Real           = 1.,
        rotation::Real       = 0.,
        canvas_width::Int64  = DEFAULT_CANVAS_WIDTH
        canvas_height::Int64 = DEFAULT_CANVAS_HEIGHT
    )

Representation of the internal camera state.

 - `position::VecE2`: the (x,y) position of the camera in [metre].
     The x-direction is measured to the east, y direction to the north.
 - `zoom::Real`: the zoom level of the camera, expressed in [pixels / metre]
 - `rotation::Real`: the rotation angle of the camera in [radian]
 - `canvas_width::Int64`: the canvas width in [pixel]
 - `canvas_height::Int64`: the canvas width in [pixel]
"""
@with_kw mutable struct CameraState
    position  :: VecE2 = VecE2(0.,0.)
    zoom      :: Real = 1.
    rotation  :: Real = 0.
    canvas_width  :: Int64 = DEFAULT_CANVAS_WIDTH
    canvas_height :: Int64 = DEFAULT_CANVAS_HEIGHT
end
position(cs::CameraState) = cs.position
zoom(cs::CameraState) = cs.zoom
rotation(cs::CameraState) = cs.rotation
canvas_width(cs::CameraState) = cs.canvas_width
canvas_height(cs::CameraState) = cs.canvas_height

camera_move!(cs::CameraState, dx::Real, dy::Real) = cs.position = cs.position + VecE2(dx, dy)
camera_move!(cs::CameraState, Δ::VecE2) = cs.position = cs.position + Δ
camera_move_pix!(cs::CameraState, dx::Real, dy::Real) = cs.position = cs.position + VecE2(dx/cs.zoom, dy/cs.zoom)
camera_move_pix!(cs::CameraState, Δ::VecE2) = cs.position = cs.position + VecE2(Δ.x/cs.zoom, Δ.y/cs.zoom)
camera_rotate!(cs::CameraState, θ::Real) = cs.rotation += θ # [radians]
camera_zoom!(cs::CameraState, factor::Real) = cs.zoom *= factor

function set_camera!(
    cs::CameraState;
    x::Real=cs.position.x,
    y::Real=cs.position.y,
    zoom::Real=cs.zoom,
    rotation::Real=cs.rotation
)
    cs.position = VecE2(x,y)
    cs.zoom = zoom
    cs.rotation = rotation
    return cs
end
reset_camera!(cs::CameraState) = set_camera!(cs, x=0., y=0., zoom=1., rotation=0.)

"""
Camera abstract type
"""
abstract type Camera end
position(c::Camera) = position(c.state)
zoom(c::Camera) = zoom(c.state)
rotation(c::Camera) = rotation(c.state)
canvas_width(c::Camera) = canvas_width(c.state)
canvas_height(c::Camera) = canvas_height(c.state)

"""
    StaticCamera <: Camera

Fix the position and the zoom as specified in the constructor.

# Constructor
`StaticCamera(position::VecE2=(0,0), zoom::Real=4)`
"""
struct StaticCamera <: Camera
    state::CameraState
end
StaticCamera(;kwargs...) = StaticCamera(CameraState(;kwargs...))
update_camera!(camera::StaticCamera, ::Scene) = camera.state

"""
    TargetFollowCamera <: Camera
Camera which follows the vehicle with ID `target_id`.
By default, the target vehicle is tracked in x and y direction.
Tracking in either direction can be disabled by setting the
`x` or `y` keys to a desired value.

# Constructor

`TargetFollowCamera(target_id; x=NaN, y=NaN, kwargs...)`

"""
mutable struct TargetFollowCamera{I} <: Camera
    state::CameraState
    target_id::I
    x::Float64
    y::Float64
end
function TargetFollowCamera(target_id; x=NaN, y=NaN, kwargs...)
    TargetFollowCamera(CameraState(;kwargs...), target_id, x, y)
end

function update_camera!(camera::TargetFollowCamera{I}, scene::Scene{E}) where {I,E<:Entity}
    target = get_by_id(scene, camera.target_id)
    x, y = posg(target.state)[1:2]
    x = isnan(camera.x) ? x : camera.x
    y = isnan(camera.y) ? y : camera.y
    set_camera!(camera.state, x=x, y=y)
    return camera.state
end

"""
    ZoomingCamera <: Camera
Camera which gradually changes the zoom level of the scene to `zoom_target` with step size `dz`.
"""
mutable struct ZoomingCamera <: Camera
    state::CameraState
    zoom_target::Float64
    dz::Float64
end
function ZoomingCamera(;zoom_target=20., dz=.5, kwargs...)
    ZoomingCamera(CameraState(;kwargs...), zoom_target, dz)
end

function update_camera!(camera::ZoomingCamera, scene::Scene{E}) where {E<:Entity}
    zt, zc = camera.zoom_target, zoom(camera)
    if zt < zc  # zooming in
        set_camera!(camera.state, zoom=max(zt, zc-camera.dz))
    elseif zt > zc  # zooming out
        set_camera!(camera.state, zoom=min(zt, zc+camera.dz))
    end
    return camera.state
end

"""
    SceneFollowCamera

Camera centered over all vehicles.

By default, the scene is tracked in x and y direction and the zoom level
is adapted to fit all vehicles in the scene. Tracking in either direction
can be disabled by setting the `x` or `y` keys to a desired value. The
zoom level can be fixed by passing a value to `zoom`.
The value of `padding` specifies the width of the additional border around
the zoomed-in area.
"""
struct SceneFollowCamera <: Camera
    state::CameraState
    x::Float64
    y::Float64
    zoom::Float64
    padding::Float64
    min_width::Float64
    min_height::Float64
end
SceneFollowCamera(; x=NaN, y=NaN, zoom=NaN, padding=4., kwargs...) = SceneFollowCamera(CameraState(;kwargs...), x, y, zoom, padding, 10, 10)
function update_camera!(camera::SceneFollowCamera, scene::Scene{E}) where {E<:Entity}
    if isnan(camera.zoom)
        pos = [posg(veh.state) for veh in scene]
        X = [p.x for p in pos]
        Y = [p.y for p in pos]
        p = camera.padding
        x_min, x_max = minimum(X)-p, maximum(X)+p
        y_min, y_max = minimum(Y)-p, maximum(Y)+p
        width = max(x_max-x_min, camera.min_width)
        height = max(y_max-y_min, camera.min_height)
        x_zoom = canvas_width(camera) / width
        y_zoom = canvas_height(camera) / height
        x_zoom, y_zoom
        x = isnan(camera.x) ? (x_min+x_max)/2 : camera.x
        y = isnan(camera.y) ? (y_min+y_max)/2 : camera.y
        zoom = min(x_zoom, y_zoom)
    else
        # TODO: is following the center of mass really the best thing to do?
        C = sum([posg(veh.state)[1:2] for veh in scene])/length(scene)  # center of mass
        x = isnan(camera.x) ? C[1] : camera.x
        y = isnan(camera.y) ? C[2] : camera.y
        zoom = camera.zoom
    end

    set_camera!(camera.state, x=x, y=y, zoom=zoom)
    return camera.state
end

"""
    ComposedCamera <: Camera
Composition of several cameras. The `update_camera` actions of the individual cameras are applied in the order in which they are saved in the `cameras` array.
States of individual cameras are ignored, the state of the composed camera is the one that will be used for rendering.

Example Usage

    cam = ComposedCamera(cameras=[SceneFollowCamera(), ZoomingCamera()])
"""
mutable struct ComposedCamera <: Camera
    state::CameraState
    cameras::Array{Camera}
end
ComposedCamera(cameras; kwargs...) = ComposedCamera(CameraState(;kwargs...), cameras)

function AutomotiveVisualization.update_camera!(camera::ComposedCamera, scene::Scene{E}) where {E<:Entity}
    for cam in camera.cameras
        cam.state = camera.state
        camera.state = update_camera!(cam, scene)
    end
    return camera.state
end

update_camera!(::Nothing, ::Scene) = nothing


text:
function veroneselift(x,liftscale,d::Integer)
    X = collect((prod(y) for y in with_replacement_combinations(x, d)))
    X = X.*liftscale

    return X
end

function veroneseliftscale(n::Integer,d::Integer)
    power = collect(multiexponents(n, d))
    df = factorial(d)
    scaling(m) = sqrt(df / prod(factorial, m))
    liftscale = collect(scaling(poly) for poly in power)
    return liftscale
end

ext: jl

max_stars_repo_name: zhemingwang/DataDrivenSwitchControl

max_issues_repo_name: zhemingwang/DataDrivenSwitchControl

max_forks_repo_name: zhemingwang/DataDrivenSwitchControl


text:
module SnowFg

using ..Ahorn, Maple

placements = Maple.SnowFg

function Ahorn.canFgBg(effect::Maple.SnowFg)
    return true, false
end

end


text:
module TestRanges

using Test
using MLJBase
using MLJTuning
using Random

# `in` for MLJType is overloaded to be `===` based. For purposed of
# testing here, we need `==` based:
function _in(x, itr)::Union{Bool,Missing}
    for y in itr
        ismissing(y) && return missing
        y == x && return true
    end
    return false
end
_issubset(itr1, itr2) = all(_in(x, itr2) for x in itr1)

mutable struct DummyModel <: Deterministic
    lambda::Float64
    metric::Float64
    kernel::Char
end

dummy_model = DummyModel(4, 9.5, 'k')

mutable struct SuperModel <: Deterministic
    K::Int64
    model1::DummyModel
    model2::DummyModel
end

dummy_model = DummyModel(1.2, 9.5, 'k')
super_model = SuperModel(4, dummy_model, deepcopy(dummy_model))

r1 = range(super_model, :(model1.kernel), values=['c', 'd'])
r2 = range(super_model, :K, lower=1, upper=10, scale=:log10)

@testset "models from cartesian range and resolutions" begin

    # with method:
    m1 = MLJTuning.grid(super_model, [r1, r2], [nothing, 7])
    m1r = MLJTuning.grid(MersenneTwister(123), super_model, [r1, r2],
                         [nothing, 7])

    # generate all models by hand:
    models1 = [SuperModel(1, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(1, DummyModel(1.2, 9.5, 'd'), dummy_model),
               SuperModel(2, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(2, DummyModel(1.2, 9.5, 'd'), dummy_model),
               SuperModel(3, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(3, DummyModel(1.2, 9.5, 'd'), dummy_model),
               SuperModel(5, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(5, DummyModel(1.2, 9.5, 'd'), dummy_model),
               SuperModel(7, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(7, DummyModel(1.2, 9.5, 'd'), dummy_model),
               SuperModel(10, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(10, DummyModel(1.2, 9.5, 'd'), dummy_model)]

    @test _issubset(models1, m1) && _issubset(m1, models1)
    @test m1r != models1
    @test _issubset(models1, m1r) && _issubset(m1, models1)

    # with method:
    m2 = MLJTuning.grid(super_model, [r1, r2], [1, 7])

    # generate all models by hand:
    models2 = [SuperModel(1, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(2, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(3, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(5, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(7, DummyModel(1.2, 9.5, 'c'), dummy_model),
               SuperModel(10, DummyModel(1.2, 9.5, 'c'), dummy_model)]

    @test _issubset(models2, m2) && _issubset(m2, models2)

end

@testset "processing user specification of range" begin
    r1 = range(Int, :h1, lower=1, upper=10)
    r2 = range(Int, :h2, lower=20, upper=30)
    s = range(Char, :j1, values = ['x', 'y'])
    @test_throws ArgumentError MLJTuning.process_user_range("junk", 42, 1)
    @test(@test_logs((:warn, r"Ignoring"),
                     MLJTuning.process_user_range((s, 3), 42, 1)) ==
          ((s, ), (2, )))
    @test MLJTuning.process_user_range(r1, 42, 1) == ((r1, ), (42, ))
    @test MLJTuning.process_user_range((r1, 3), 42, 1) == ((r1, ), (3, ))
    @test MLJTuning.process_user_range(s, 42, 1) == ((s, ), (2,))
    @test MLJTuning.process_user_range([(r1, 3), r2, s], 42, 1) ==
        ((r1, r2, s), (3, 42, 2))
end

end
true


text:
mutable struct GLFunc
    p::Ptr{Cvoid}
end

include(joinpath("..", "deps", "deps.jl"))

gl_represent(x::GLenum) = GLENUM(x).name
gl_represent(x) = repr(x)

function debug_opengl_expr(func_name, args)
    if enable_opengl_debugging && func_name != :glGetError
        quote
            err = glGetError()
            if err != GL_NO_ERROR
                arguments = gl_represent.(tuple($(args...)))
                warn("OpenGL call to $($func_name), with arguments: $(arguments)
                Failed with error: $(GLENUM(err).name).")
            end
        end
    else
        :()
    end
end

# based on getCFun macro
macro glfunc(opengl_func)
    arguments = map(opengl_func.args[1].args[2:end]) do arg
        isa(arg, Symbol) && return Expr(:(::), arg)
        arg
    end
    # Get info out of arguments of `opengl_func`
    arg_names       = map(arg->arg.args[1], arguments)
    return_type     = opengl_func.args[2]
    input_types     = map(arg->arg.args[2], arguments)
    func_name       = opengl_func.args[1].args[1]
    func_name_sym   = Expr(:quote, func_name)
    func_name_str   = string(func_name)
    ptr_expr        = :(getprocaddress_e($func_name_str))

    if iswindows() # windows has some function pointers statically available and some not, this is how we deal with it:
        ptr = Libdl.dlsym_e(gl_lib, func_name)
        if (ptr != C_NULL)
            ptr_expr = :(($func_name_sym, "opengl32"))
            ret = quote
                function $func_name($(arg_names...))
                    result = ccall($ptr_expr, $return_type, ($(input_types...),), $(arg_names...))
                    $(debug_opengl_expr(func_name, arg_names))
                    result
                end
                $(Expr(:export, func_name))
            end
            return esc(ret)
        end
    end
    ptr_sym = gensym("$(func_name)_func_pointer")
    ret = quote
        $ptr_sym = GLFunc(C_NULL)
        function $func_name($(arg_names...))
            if $ptr_sym.p::Ptr{Cvoid} == C_NULL
                $ptr_sym.p::Ptr{Cvoid} = $ptr_expr
            end
            result = ccall($ptr_sym.p::Ptr{Cvoid}, $return_type, ($(input_types...),), $(arg_names...))
            $(debug_opengl_expr(func_name, arg_names))
            result
        end
        $(Expr(:export, func_name))
        end
    return esc(ret)
end

if iswindows()
    const gl_lib = Libdl.dlopen("opengl32")
end
include("glFunctions.jl")
if iswindows()
    Libdl.dlclose(gl_lib)
end


text:
function install_env_path(
    m::Module,
    options::Configs.Comonicon = Configs.read_options(m);
    shell::String = basename(ENV["SHELL"]),
    home_dir = homedir(),
    env = ENV,
    yes::Bool = false,
)

    rcfile = detect_rcfile(shell, home_dir)
    install_path = expanduser(options.install.path)
    msg = "cannot detect $(options.install.path)/bin in PATH, do you want to add it in PATH?"
    if !contains_path(rcfile, install_path, env) && Tools.prompt(msg, yes)
        write_path(rcfile, install_path)
    end

    msg = "cannot detect $(options.install.path)/completions in FPATH, do you want to add it in FPATH?"
    if !contains_fpath(rcfile, install_path, env) && Tools.prompt(msg, yes)
        write_fpath(rcfile, install_path)
    end

    @info "open a new terminal, or source $rcfile to enable the new PATH."
    return
end

function contains_path(rcfile::String, path::String, env = ENV)
    bin = joinpath(path, "bin")
    return contains_shell_var(rcfile, "PATH", bin, env)
end

function contains_fpath(rcfile::String, path::String, env = ENV)
    completions = joinpath(path, "completions")
    return contains_shell_var(rcfile, "FPATH", completions, env)
end

function write_path(rcfile::String, path::String)
    bin = joinpath(path, "bin")
    return export_shell_var(rcfile, "PATH", bin)
end

function write_fpath(rcfile::String, path::String)
    completions = joinpath(path, "completions")
    return export_shell_var(rcfile, "FPATH", completions)
end

function contains_shell_var(rcfile, var, value, env)
    if !haskey(env, var)
        for line in readlines(rcfile)
            if strip(line) == "export $var=\"$(relhome(value)):\$$var\"" ||
               strip(line) == "export $var=\"$value:\$$var\""
                return true
            end
        end
        return false
    end

    for each in split(env[var], ":")
        (each == relhome(value) || each == value) && return true
    end
    return false
end

function export_shell_var(rcfile, var, path)
    isempty(rcfile) && return
    script = """
    # generated by Comonicon
    export $var="$(relhome(path)):\$$var"
    """
    open(rcfile, "a") do io
        write(io, "\n" * script * "\n")
    end
end

relhome(path) = joinpath("\$HOME", relpath(path, homedir()))

function detect_rcfile(shell::String = basename(ENV["SHELL"]), home = homedir())
    if shell == "zsh"
        rcfile = joinpath((haskey(ENV, "ZDOTDIR") ? ENV["ZDOTDIR"] : home), ".zshrc")
    elseif shell == "bash"
        rcfile = joinpath(home, ".bashrc")
        if !isfile(rcfile)
            rcfile = joinpath(home, ".bash_profile")
        end
        return rcfile
    else
        error(
            "auto installation for $shell is not supported, " *
            "please open an issue under $COMONICON_URL",
        )
    end
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdpageheader

"""
    'feffery'_antdpageheader(;kwargs...)
    'feffery'_antdpageheader(children::Any;kwargs...)
    'feffery'_antdpageheader(children_maker::Function;kwargs...)

An AntdPageHeader component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): The content of the tab - will only be displayed if this tab is selected
- `id` (String; optional)
- `backClicks` (Real; optional)
- `className` (String; optional)
- `ghost` (Bool; optional)
- `historyBackDisabled` (Bool; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `showBackIcon` (Bool; optional)
- `style` (Dict; optional)
- `subTitle` (String; optional)
- `title` (String; optional)
"""
function 'feffery'_antdpageheader(; kwargs...)
        available_props = Symbol[:children, :id, :backClicks, :className, :ghost, :historyBackDisabled, :loading_state, :showBackIcon, :style, :subTitle, :title]
        wild_props = Symbol[]
        return Component("'feffery'_antdpageheader", "AntdPageHeader", "feffery_antd_components", available_props, wild_props; kwargs...)
end

'feffery'_antdpageheader(children::Any; kwargs...) = 'feffery'_antdpageheader(;kwargs..., children = children)
'feffery'_antdpageheader(children_maker::Function; kwargs...) = 'feffery'_antdpageheader(children_maker(); kwargs...)

max_stars_repo_name: ReilYoung/feffery-antd-components

max_issues_repo_name: ReilYoung/feffery-antd-components

max_forks_repo_name: ReilYoung/feffery-antd-components


text:
function replaceind_indval(IV::Tuple, iĩ::Pair)
  i, ĩ = iĩ
  return ntuple(n -> first(IV[n]) == i ? ĩ => last(IV[n]) : IV[n], length(IV))
end

# atol controls the tolerance cutoff for determining which eigenvectors are in the null
# space of the isometric MPS tensors. Setting to 1e-2 since we only want to keep
# the eigenvectors corresponding to eigenvalues of approximately 1.
function subspace_expansion(
  ψ::InfiniteCanonicalMPS, H, b::Tuple{Int,Int}; maxdim, cutoff, atol=1e-2, kwargs...
)
  n1, n2 = b
  lⁿ¹ = commoninds(ψ.AL[n1], ψ.C[n1])
  rⁿ¹ = commoninds(ψ.AR[n2], ψ.C[n1])

  dˡ = dim(lⁿ¹)
  dʳ = dim(rⁿ¹)
  @assert dˡ == dʳ
  if dˡ ≥ maxdim
    println(
      "Current bond dimension at bond $b is $dˡ while desired maximum dimension is $maxdim, skipping bond dimension increase",
    )
    return (ψ.AL[n1], ψ.AL[n2]), ψ.C[n1], (ψ.AR[n1], ψ.AR[n2])
  end
  maxdim -= dˡ

  # Returns `NL` such that `norm(ψ.AL[n1] * NL) ≈ 0`
  NL = nullspace(ψ.AL[n1], lⁿ¹; atol=atol)
  NR = nullspace(ψ.AR[n2], rⁿ¹; atol=atol)

  nL = uniqueinds(NL, ψ.AL[n1])
  nR = uniqueinds(NR, ψ.AR[n2])

  ψH2 = noprime(ψ.AL[n1] * H[(n1, n2)] * ψ.C[n1] * ψ.AR[n2])
  ψHN2 = ψH2 * NL * NR

  U, S, V = svd(ψHN2, nL; maxdim=maxdim, cutoff=cutoff, kwargs...)
  if dim(S) == 0
    return (ψ.AL[n1], ψ.AL[n2]), ψ.C[n1], (ψ.AR[n1], ψ.AR[n2])
  end
  @show S[end, end]
  NL *= dag(U)
  NR *= dag(V)

  ALⁿ¹, l = ITensors.directsum(
    ψ.AL[n1], dag(NL), uniqueinds(ψ.AL[n1], NL), uniqueinds(NL, ψ.AL[n1]); tags=("Left",)
  )
  ARⁿ², r = ITensors.directsum(
    ψ.AR[n2], dag(NR), uniqueinds(ψ.AR[n2], NR), uniqueinds(NR, ψ.AR[n2]); tags=("Right",)
  )

  C = ITensor(dag(l)..., dag(r)...)
  ψCⁿ¹ = permute(ψ.C[n1], lⁿ¹..., rⁿ¹...)
  for I in eachindex(ψ.C[n1])
    v = ψCⁿ¹[I]
    if !iszero(v)
      C[I] = ψCⁿ¹[I]
    end
  end

  # Also expand the dimension of the neighboring MPS tensors
  ALⁿ² = ITensor(dag(l)..., uniqueinds(ψ.AL[n2], ψ.AL[n1])...)

  il = only(uniqueinds(ψ.AL[n2], ALⁿ²))
  ĩl = only(uniqueinds(ALⁿ², ψ.AL[n2]))
  for IV in eachindval(inds(ψ.AL[n2]))
    ĨV = replaceind_indval(IV, il => ĩl)
    v = ψ.AL[n2][IV...]
    if !iszero(v)
      ALⁿ²[ĨV...] = v
    end
  end

  ARⁿ¹ = ITensor(dag(r)..., uniqueinds(ψ.AR[n1], ψ.AR[n2])...)

  ir = only(uniqueinds(ψ.AR[n1], ARⁿ¹))
  ĩr = only(uniqueinds(ARⁿ¹, ψ.AR[n1]))
  for IV in eachindval(inds(ψ.AR[n1]))
    ĨV = replaceind_indval(IV, ir => ĩr)
    v = ψ.AR[n1][IV...]
    if !iszero(v)
      ARⁿ¹[ĨV...] = v
    end
  end

  CL = combiner(l; tags=tags(only(lⁿ¹)))
  CR = combiner(r; tags=tags(only(rⁿ¹)))
  ALⁿ¹ *= CL
  ALⁿ² *= dag(CL)
  ARⁿ² *= CR
  ARⁿ¹ *= dag(CR)
  C = (C * dag(CL)) * dag(CR)

  # TODO: delete or only print when verbose
  ## ψ₂ = ψ.AL[n1] * ψ.C[n1] * ψ.AR[n2]
  ## ψ̃₂ = ALⁿ¹ * C * ARⁿ²
  ## local_energy(ψ, H) = (noprime(ψ * H) * dag(ψ))[]

  return (ALⁿ¹, ALⁿ²), C, (ARⁿ¹, ARⁿ²)
end

function subspace_expansion(ψ, H; kwargs...)
  ψ = copy(ψ)
  N = nsites(ψ)
  AL = ψ.AL
  C = ψ.C
  AR = ψ.AR
  for n in 1:N
    n1, n2 = n, n + 1
    ALⁿ¹², Cⁿ¹, ARⁿ¹² = subspace_expansion(ψ, H, (n1, n2); kwargs...)
    ALⁿ¹, ALⁿ² = ALⁿ¹²
    ARⁿ¹, ARⁿ² = ARⁿ¹²
    AL[n1] = ALⁿ¹
    AL[n2] = ALⁿ²
    C[n1] = Cⁿ¹
    AR[n1] = ARⁿ¹
    AR[n2] = ARⁿ²
    ψ = InfiniteCanonicalMPS(AL, C, AR)
  end
  return ψ
end


text:
# module Bukdu

export Assoc

"""
    Assoc
"""
struct Assoc
    __bukdu_assoc::Vector{Pair{String,String}}
end
Assoc(pair::Pair) = Assoc([pair])
Assoc(ps...) = Assoc(collect(ps))

function Base.getindex(assoc::Assoc, key::String)::String
    for (k, v) in assoc.__bukdu_assoc
        k == key && return v
    end
    string()
end

function Base.getindex(assoc::Assoc, key::Symbol)::String
    getindex(assoc, String(key))
end

function Base.setindex!(assoc::Assoc, value::String, key::Symbol)
    setindex!(assoc, value, String(key))
end

function Base.setindex!(assoc::Assoc, value::String, key::String)
    ind = something(findfirst(isequal(key), keys(assoc)), 0)
    if ind > 0
        assoc.__bukdu_assoc[ind] = Pair(key, value)
    else
        push!(assoc, Pair(key, value))
    end
end

function Base.getproperty(assoc::Assoc, prop::Symbol)
    if :__bukdu_assoc == prop
        getfield(assoc, prop)
    else
        getindex(assoc, String(prop))
    end
end

function Base.length(assoc::Assoc)
    length(assoc.__bukdu_assoc)
end

function Base.isempty(assoc::Assoc)::Bool
    isempty(assoc.__bukdu_assoc)
end

function Base.empty!(assoc::Assoc)
    empty!(assoc.__bukdu_assoc)
end

function Base.haskey(assoc::Assoc, key::String)::Bool
    key in keys(assoc)
end

function Base.haskey(assoc::Assoc, key::Symbol)::Bool
    String(key) in keys(assoc)
end

function Base.keys(assoc::Assoc)::Vector{String}
    getindex.(assoc.__bukdu_assoc, 1)
end

function Base.values(assoc::Assoc)::Vector{String}
    getindex.(assoc.__bukdu_assoc, 2)
end

function Base.pairs(assoc::Assoc)::Base.Iterators.Pairs
    pairs(assoc.__bukdu_assoc)
end

function Base.push!(assoc::Assoc, kv::Pair{String,String})
    push!(assoc.__bukdu_assoc, kv)
end

function Base.:(==)(left::Assoc, right::Assoc)::Bool
    left.__bukdu_assoc == right.__bukdu_assoc
end

function Base.get(assoc::Assoc, key::Symbol, value::Any)
    if haskey(assoc, key)
        v = assoc[key]
        value isa String ? v : parse(typeof(value), v)
    else
        value
    end
end

function Base.merge(ps::T...)::T where {T <: Vector{Pair{String,String}}}
    assoc = Assoc()
    for p in ps
        for (k, v) in p
            assoc[k] = v
        end
    end
    assoc.__bukdu_assoc
end

function Base.show(io::IO, mime::MIME"text/plain", assoc::Assoc)
    body = join(assoc.__bukdu_assoc, ", ")
    println(io, Assoc, '(', body, ')')
end

function Base.iterate(assoc::Assoc, v::Int=1)
    iterate(assoc.__bukdu_assoc, v)
end

Base.lastindex(assoc::Assoc) = lastindex(assoc.__bukdu_assoc)

import JSON2
JSON2.write(io::IO, assoc::Assoc) = JSON2.write(io, assoc.__bukdu_assoc)

# module Bukdu

max_stars_repo_name: wookay/TestBukdu

max_issues_repo_name: wookay/TestBukdu

max_forks_repo_name: wookay/TestBukdu


text:
using ForwardDiff
using Setfield: @set

"""
Find vectors tJ2 and cotJ such that ker(J) = span{tJ1, tJ2} and
ker(J') = span{cotJ} (left nullspace or cokernel):
"""
function find_more_nullspaces(Q, L, rtol, atol, max_steps)
    y = _zeros(Q, size(Q, 1) - 1)
    if y isa SArray
        y = @set y[length(y)] = 1   # y[end] doesn't work
    else
        y[end] = 1
    end

    #=
    if abs(L[end-1, end-1]) < atol
        tJ2 = Q[end-1, :]
        cotJ = (@view Q[1:end-1, 1:end-1]) \ y
        return tJ2, cotJ
    end
    =#

    L2 = popbottomright(L)
    R2 = L2'
    if Q isa StaticArray
        L2 = LowerTriangular(SMatrix{size(L2)...}(L2))
        R2 = UpperTriangular(SMatrix{size(R2)...}(R2))
    end
    y, cotJ = _find_more_nullspaces(L2, R2, y, rtol, atol, max_steps)
    tJ2 = (Q[1:end-1, :])' * y

    if y isa SVector  # TODO: don't
        return SVector(tJ2...), cotJ
    end

    return tJ2, cotJ
end

function _find_more_nullspaces(L2, R2,
                               y :: T,
                               rtol, atol, max_steps,
                               ) where {T <: AbstractVector}
    if abs(det(L2)) < atol
        # TODO: optimize
        ker_L2 = nullspace(Array(L2))
        ker_R2 = nullspace(Array(R2))
        if size(ker_L2, 2) > 0 && size(ker_R2, 2) > 0
            return (T(ker_L2[:, 1]), T(ker_R2[:, 1]))
        end
        # Otherwise, let's fallback to the manual method.
    end
    x0 = _similar(y, length(y))
    x1 = _similar(y, length(y))

    x1 = _normalize!(R2 \ y)
    y = _normalize!(L2 \ x1)
    for _ in 2:max_steps
        (x0, x1) = (x1, x0)
        x1 = _normalize!(R2 \ y)
        y = _normalize!(L2 \ x1)
        if isapprox(x0, x1; rtol=rtol, atol=atol)
            return y::T, x1::T
        end
    end
    error("Failed to find the cokernel.")
end

@with_kw struct SimpleBifurcationSolution{RV, LV, M}
    found::Bool
    tv1::RV
    tv2::RV
    tJ1::RV
    tJ2::RV
    cotJ::LV
    hess::M
end

function solve_simple_bifurcation!(cache, opts,
                                   u0::TV, tJ::TV,
                                   L, Q,
                                   ) where {TV <: AbstractVector}
    rtol = opts.rtol
    atol = opts.atol

    tJ1 = tJ
    tJ2, cotJ = find_more_nullspaces(Q, L, rtol, atol, opts.max_misc_steps)
    tJ2 :: TV

    function g(xi)
        prob_cache = cache.prob_cache
        v = u0 .+ xi[1] * tJ1 .+ xi[2] * tJ2
        return cotJ ⋅ residual(v, prob_cache)
    end
    hess = ForwardDiff.hessian(g, [0.0, 0.0])
    H11, H12, _, H22 = hess

    # Solve (ξ₁, ξ₂)ᵀ H (ξ₁, ξ₂) = 0 where ξ₁=x and ξ₂=1-x
   discriminant = - H11 * H22 + H12^2  # = det(hess)
    if discriminant <= 0
        tv1 = tJ1 .* NaN
        tv2 = tJ2 .* NaN
        return SimpleBifurcationSolution(false, tv1, tv2, tJ1, tJ2, cotJ, hess)
    end

    x1 = (- H12 + H22 + sqrt(discriminant)) / (H11 + H22 - 2 * H12)
    x2 = (- H12 + H22 - sqrt(discriminant)) / (H11 + H22 - 2 * H12)
    if abs(x1) < abs(x2)
        x1, x2 = x2, x1
    end
    tv1 = @. x1 * tJ1 + (1 - x1) * tJ2
    tv2 = @. x2 * tJ1 + (1 - x2) * tJ2
    tv1 /= norm(tv1)
    tv2 /= norm(tv2)

    return SimpleBifurcationSolution(true, tv1, tv2, tJ1, tJ2, cotJ, hess)
end

function new_branches!(cache, opts, sbint::SimpleBifurcationInterval)
    u0, tJ, L, Q = try
        find_simple_bifurcation!(cache, opts, sbint)
    catch err
        # Maybe there is a singularity (around which det(J; tJ)
        # flips) within `sbint` and the root finder hits that?
        if err isa SingularException
            @warn(
                """
                Failed to find bifurcation point within:
                $(sbint)
                """,
                exception = err)
            return []
        end
        rethrow()
    end

    sbsol = solve_simple_bifurcation!(cache, opts, u0, tJ, L, Q)
    @unpack tv1, tv2 = sbsol
    @assert sbsol.found  # TODO: nicer error handling

    # Choose the direction `tv` of the new branch.  Use the one least
    # parallel to the direction `tv0` along the previous curve.
    tv0 = sbint.u1 .- sbint.u0
    tv0 /= norm(tv0)
    if abs(tv1 ⋅ tv0) > abs(tv2 ⋅ tv0)
        tv = tv2
    else
        tv = tv1
    end

    args = (sbint.h, sbint.direction)
    u1, h1 = predictor_corrector_step!(cache, opts, u0, tv, args...)
    simple_bifurcation1 = cache.simple_bifurcation
    u2, h2 = predictor_corrector_step!(cache, opts, u0, -tv, args...)
    simple_bifurcation2 = cache.simple_bifurcation

    # Maybe return BranchResult?
    check(BranchResult(
        tv, tv0, tv1, tv2, u0, u1, u2,
        simple_bifurcation1, simple_bifurcation2,
    ), opts)

    return [
        (u0, u1, sbint.direction, h1),
        (u0, u2, sbint.direction, h2),
    ]
end

struct BranchResult{T}
    tv::T
    tv0::T
    tv1::T
    tv2::T
    u0::T
    u1::T
    u2::T
    simple_bifurcation1::Bool
    simple_bifurcation2::Bool
end

struct BranchError{T} <: Exception
    result::BranchResult{T}
    opts::ContinuationOptions
end

function check(result::BranchResult, opts)
    @unpack tv, tv0, tv1, tv2, u0, u1, u2,
        simple_bifurcation1, simple_bifurcation2 = result

    parallel_to_old = du -> abs(du ⋅ tv0) > abs(du ⋅ tv)

    if simple_bifurcation1 || simple_bifurcation2 ||
            parallel_to_old(u1 .- u0) || parallel_to_old(u2 .- u0) ||
            isapprox(u1, u0; atol=opts.atol, rtol=opts.rtol) ||
            isapprox(u2, u0; atol=opts.atol, rtol=opts.rtol)
        # throw(BranchError(result, opts))
        @warn("Failed to branch.",
              exception=BranchError(result, opts))
    end
end

function Base.showerror(io::IO, err::BranchError)
    @unpack tv, tv0, tv1, tv2, u0, u1, u2,
        simple_bifurcation1, simple_bifurcation2 = err.result
    opts = err.opts

    sep = "\n    * "
    _print = (args...) -> print(io, sep, args...)

    parallel_to_old = du -> abs(du ⋅ tv0) > abs(du ⋅ tv)

    print(io, "Possible errors while branching off from old curve:")

    if simple_bifurcation1 || simple_bifurcation2
        _print("Simple bifurcation during the first corrector.")
        if ! simple_bifurcation1
            _print("Only the second branch had simple bifurcation.")
        elseif ! simple_bifurcation2
            _print("Only the first branch had simple bifurcation.")
        end
    end
    if parallel_to_old(u1 .- u0) || parallel_to_old(u2 .- u0)
        _print("New branch candidates are more parallel to the old curve",
               " than the computed new direction.")
    end
    if (isapprox(u1, u0; atol=opts.atol, rtol=opts.rtol) ||
        isapprox(u2, u0; atol=opts.atol, rtol=opts.rtol))
        _print("New points are approximately equal to the branch point.")
    end
end


text:
@testset "Medium" begin

    # #Test we can call the variable and the function
    # @test begin m = @agent(1, u::Medium); m.declaredSymbols["Medium"] == [:u] end
    # @test begin m = @agent(1, [u,v]::Medium); m.declaredSymbols["Medium"] == [:u,:v] end
    # @test begin m = @agent(1, UpdateMedium = ∂t_u = Δ(u) + δx(0.)*δy(0.)*1.); "UpdateMedium" in keys(m.declaredUpdates) end

    # # Declaration of Mediums
    # for j in ["Newmann",
    #     "Dirichlet",
    #     "Periodic",
    #     "Newmann-Newmann",
    #     "Dirichlet-Dirichlet",
    #     "Periodic-Periodic",
    #     "Newmann-Dirichlet",
    #     "Dirichlet-Newmann"]

    #     @test_nowarn MediumFlat(j,10)
    # end

    # for j in ["Newmann-Periodic",
    #     "Dirichlet-Periodic",
    #     "Periodic-Dirichlet",
    #     "Periodic-Newmann"]

    #     @test_throws ErrorException MediumFlat(j,10)
    # end

    # @test_throws ErrorException MediumFlat("k",10)

    # # Declaration of MediumFlat in spaces
    # for i in 1:3

    #     if i == 1
    #         m = @agent(1, u::Medium)
    #     elseif i == 2
    #         m = @agent(2, u::Medium)
    #     else
    #         m = @agent(3, u::Medium)
    #     end

    #     box = [(:x,-10.,10.),(:y,-10.,10.),(:z,-10.,10.)][1:i]
    #     medium = [MediumFlat("Newmann",10) for k in 1:i]

    #     @test_nowarn SimulationFree(m,box=box,medium=medium)
    #     @test_nowarn SimulationGrid(m,box,1.,medium=medium)
    # end

    # m = @agent(3, u::Medium)
    # @test_throws ErrorException SimulationFree(m,box=[(:x,0,1)],medium=[MediumFlat("Dirichlet",10)])
    # @test_throws ErrorException SimulationFree(m,box=[(:x,0,1)],medium=[MediumFlat("Dirichlet",10),MediumFlat("Dirichlet",10)])

    # # Compilation
    # for i in testplatforms
    #     @test_nowarn begin
    #         m = @agent(3, u::Medium)
    #         s = SimulationFree(m,
    #                         box=[(:x,-10.,10.),(:y,-10.,10.),(:z,-10.,10.)],
    #                         medium=[MediumFlat("Dirichlet",10),MediumFlat("Dirichlet",10),MediumFlat("Dirichlet",10)])
    #         compile(m,s,platform=i)
    #     end
    # end

    # #Community construction
    # @test_nowarn begin
    #     m = @agent(3, u::Medium)
    #     s = SimulationFree(m,
    #                     box=[(:x,-10.,10.),(:y,-10.,10.),(:z,-10.,10.)],
    #                     medium=[MediumFlat("Dirichlet",10),MediumFlat("Dirichlet",10),MediumFlat("Dirichlet",10)])
    #     mc = compile(m,s)

    #     Community(mc,N=10)
    # end

    # #Medium access in Community
    # @test_nowarn begin
    #     m = @agent(3, [u,v]::Medium)
    #     s = SimulationFree(m,
    #                     box=[(:x,-10.,10.),(:y,-10.,10.),(:z,-10.,10.)],
    #                     medium=[MediumFlat("Dirichlet",10),MediumFlat("Dirichlet",10),MediumFlat("Dirichlet",10)])
    #     mc = compile(m,s)

    #     com = Community(mc,N=10)

    #     b = com.u #Read
    #     com.v .= 1 #Assign scalar
    #     com.v .= zeros(10,10,10) #Assign vector
    #     com.v[1,2,3] = 1; #Assign single value
    # end

    # m = @agent(3, u::Medium)
    # s = SimulationFree(m,
    #                 box=[(:x,-10.,10.),(:y,-10.,10.),(:z,-10.,10.)],
    #                 medium=[MediumFlat("Dirichlet",10),MediumFlat("Dirichlet",10),MediumFlat("Dirichlet",10)])
    # mc = compile(m,s)
    # com = Community(mc,N=10)
    # com.u .= 0
    # comt = mc.evolve(com,dt=0.1,tMax=1)
    # @test all(comt[end].u .== zeros(10,10,10))

    # # Evolve Community Save Working
    # for i in testplatforms
    #     for j in ["RAM","CSV"]
    #         for medium in [[MediumFlat("Dirichlet",10),MediumFlat("Dirichlet",10),MediumFlat("Dirichlet",10)],
    #                 [MediumFlat("Periodic",10),MediumFlat("Dirichlet",10),MediumFlat("Dirichlet",10)],
    #                 [MediumFlat("Dirichlet",10),MediumFlat("Periodic",10),MediumFlat("Dirichlet",10)],
    #                 [MediumFlat("Dirichlet",10),MediumFlat("Dirichlet",10),MediumFlat("Periodic",10)],
    #                 [MediumFlat("Periodic",10),MediumFlat("Periodic",10),MediumFlat("Dirichlet",10)],
    #                 [MediumFlat("Periodic",10),MediumFlat("Dirichlet",10),MediumFlat("Periodic",10)],
    #                 [MediumFlat("Dirichlet",10),MediumFlat("Periodic",10),MediumFlat("Periodic",10)],
    #                 [MediumFlat("Periodic",10),MediumFlat("Periodic",10),MediumFlat("Periodic",10)]
    #                         ]

    #             m = @agent(3, [u,v]::Medium)
    #             s = SimulationFree(m,
    #                             box=[(:x,-10.,10.),(:y,-10.,10.),(:z,-10.,10.)],
    #                             medium=medium)
    #             mc = compile(m,s,platform=i,save=j)
    #             com = Community(mc,N=10)
    #             @test_nowarn mc.evolve(com,dt=0.1,tMax=10.)
    #         end

    #         dir = readdir("./")
    #         for i in dir
    #             if occursin(".csv",i)
    #                 rm(i)
    #             end
    #         end
    #     end
    # end

    # # Test appropiate Update working
    # for i in testplatforms
    #     m = @agent(3, [u,v]::Medium, UpdateMedium = ∂t_u = 1)
    #     s = SimulationFree(m,
    #                     box=[(:x,-1.,1.),(:y,-1.,1.),(:z,-1.,1.)],
    #                     medium=[MediumFlat("Newmann",10),MediumFlat("Newmann",10),MediumFlat("Newmann",10)])
    #     mc = compile(m,s,platform=i)
    #     com = Community(mc,N=10)
    #     com.u .= 0
    #     com.v .= 1
    #     @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
    #     @test begin
    #         comt = mc.evolve(com,dt=0.1,tMax=.2)

    #         all(comt[end].u .≈ 0.2)
    #     end
    #     @test begin
    #         comt = mc.evolve(com,dt=0.1,tMax=.2)

    #         all(comt[end].v .≈ 1.)
    #     end
    # end

    # for i in testplatforms #Test single deltas working
    #     m = nothing
    #     for j in 1:3
    #         if j == 1
    #             m = @agent(3, [u,v]::Medium, UpdateMedium = ∂t_u = δx(0.)*1)
    #         elseif j == 2
    #             m = @agent(3, [u,v]::Medium, UpdateMedium = ∂t_u = δy(0.)*1)
    #         elseif j == 3
    #             m = @agent(3, [u,v]::Medium, UpdateMedium = ∂t_u = δz(0.)*1)
    #         end

    #         s = SimulationFree(m,
    #                         box=[(:x,-1.,1.),(:y,-2.,2.),(:z,-3.,3.)],
    #                         medium=[MediumFlat("Newmann",10),MediumFlat("Newmann",100),MediumFlat("Newmann",1000)])
    #         mc = compile(m,s,platform=i)
    #         com = Community(mc,N=10)
    #         com.u .= 0
    #         com.v .= 1
    #         comt = mc.evolve(com,dt=0.1,tMax=.2)
    #         @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
    #         if j == 1
    #             @test sum(comt[end].u[5,:,:] .≈ 0.2) == 100000
    #         elseif j == 2
    #             @test sum(comt[end].u[:,50,:] .≈ 0.2) == 10000
    #         elseif j == 3
    #             @test sum(comt[end].u[:,:,500] .≈ 0.2) == 1000
    #         end
    #     end
    # end

    # for i in testplatforms #Test joint deltas working
    #     m = @agent(3, [u,v]::Medium, UpdateMedium = ∂t_u = δx(0.)*δy(0.)*δz(0.)*1)
    #     s = SimulationFree(m,
    #                     box=[(:x,-1.,1.),(:y,-1.,1.),(:z,-1.,1.)],
    #                     medium=[MediumFlat("Newmann",10),MediumFlat("Newmann",10),MediumFlat("Newmann",10)])
    #     mc = compile(m,s,platform=i)
    #     com = Community(mc,N=10)
    #     com.u .= 0
    #     com.v .= 1
    #     @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
    #     @test begin
    #         comt = mc.evolve(com,dt=0.1,tMax=.2)

    #         sum(comt[end].u .≈ 0.2) == 1
    #     end
    # end

    # for i in testplatforms #Check diffussion is conservative Periodic

    #     #1D
    #     m = @agent(1, [u,v]::Medium, UpdateMedium = ∂t_u = Δx(u))
    #     s = SimulationFree(m,
    #                     box=[(:x,-20.,20.)],
    #                     medium=[MediumFlat("Periodic",100)])
    #     mc = compile(m,s,platform=i)
    #     com = Community(mc,N=10)
    #     com.u .= pdf.(Normal(0,0.5),range(-20,20,length=100))
    #     # println(mc.program)
    #     @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
    #     @test begin
    #         comt = mc.evolve(com,dt=0.01,tMax=1000,dtSave=100)

    #         sum(comt[end].u) .≈ sum(comt[1].u)
    #     end
    #     @test begin
    #         comt = mc.evolve(com,dt=0.01,tMax=1000,dtSave=100)

    #         maximum(comt[end].u) .≈ minimum(comt[end].u)
    #     end

    #     #2D
    #     m = @agent(2, [u,v]::Medium, UpdateMedium = ∂t_u = Δx(u) + Δy(u))
    #     s = SimulationFree(m,
    #                     box=[(:x,-20.,20.),(:y,-20.,20.)],
    #                     medium=[MediumFlat("Periodic",100),MediumFlat("Periodic",100)])
    #     mc = compile(m,s,platform=i)
    #     com = Community(mc,N=10)
    #     com.u .= 1
    #     # println(mc.program)
    #     @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
    #     @test begin
    #         comt = mc.evolve(com,dt=0.01,tMax=10,dtSave=1)
    #         (sum(comt[end].u) .- sum(comt[1].u))/10000 < 0.01
    #     end

    #     #3D
    #     m = @agent(3, [u,v]::Medium, UpdateMedium = ∂t_u = Δx(u) + Δy(u) + Δz(u))
    #     s = SimulationFree(m,
    #                     box=[(:x,-20.,20.),(:y,-20.,20.),(:z,-20.,20.)],
    #                     medium=[MediumFlat("Periodic",100),MediumFlat("Periodic",100),MediumFlat("Periodic",100)])
    #     mc = compile(m,s,platform=i)
    #     com = Community(mc,N=10)
    #     com.u .= 1
    #     # println(mc.program)
    #     @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
    #     @test begin
    #         comt = mc.evolve(com,dt=0.01,tMax=10,dtSave=1)

    #         (sum(comt[end].u) .- sum(comt[1].u))/1000000 < 0.01
    #     end

    # end

    # for i in testplatforms #Check diffussion disappears Dirichlet

    #     m = @agent(1, [u,v]::Medium, UpdateMedium = ∂t_u = Δx(u))
    #     s = SimulationFree(m,
    #                     box=[(:x,-20.,20.)],
    #                     medium=[MediumFlat("Dirichlet",100)])
    #     mc = compile(m,s,platform=i)
    #     com = Community(mc,N=10)
    #     com.u .= pdf.(Normal(0,0.5),range(-20,20,length=100))
    #     # println(mc.program)
    #     @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
    #     @test begin
    #         comt = mc.evolve(com,dt=0.01,tMax=2000,dtSave=100)

    #         sum(comt[end].u) < 0.001
    #     end

    # end

    # for i in testplatforms #Check diffussion conservative when Newmann

    #     m = @agent(1, [u,v]::Medium, UpdateMedium = ∂t_u = Δx(u))
    #     s = SimulationFree(m,
    #                     box=[(:x,-20.,20.)],
    #                     medium=[MediumFlat("Newmann",100)])
    #     mc = compile(m,s,platform=i)
    #     com = Community(mc,N=10)
    #     com.u .= pdf.(Normal(0,0.5),range(-20,20,length=100))
    #     # println(mc.program)
    #     @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
    #     @test begin
    #         comt = mc.evolve(com,dt=0.01,tMax=2000,dtSave=100)

    #         sum(comt[end].u)-sum(comt[1].u) < 0.1
    #     end

    # end

    # for i in testplatforms #Check steady state with source

    #     m = @agent(1, [u,v]::Medium, UpdateMedium = ∂t_u = Δx(u) + δx(0.))
    #     s = SimulationFree(m,
    #                     box=[(:x,-20.,20.)],
    #                     medium=[MediumFlat("Dirichlet",100)])
    #     mc = compile(m,s,platform=i)
    #     com = Community(mc,N=10)
    #     com.u .= 0
    #     # println(mc.program)
    #     @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
    #     @test begin
    #         comt = mc.evolve(com,dt=0.01,tMax=4000,dtSave=200)

    #         all(abs.(comt[end].u .- comt[10].u) .< 0.001)
    #     end

    # end

    # for i in testplatforms #Check coumpling agents with medium

    #     m = @agent(1,
    #             p::Local,
    #             [u,v]::Medium,
    #             UpdateMedium = begin
    #                 ∂t_u = Δx(u)
    #                 ∂t_v = Δx(v) + δx(0.)
    #             end,
    #             UpdateMediumInteraction = u += p*dt
    #             )
    #     s = SimulationFree(m,
    #                     box=[(:x,-20.,20.)],
    #                     medium=[MediumFlat("Dirichlet",100)])
    #     mc = compile(m,s,platform=i)
    #     com = Community(mc,N=1)
    #     com.u .= 0.
    #     com.v .= 0.
    #     com.x .= 0.
    #     com.p .= 1.
    #     # println(mc.program)
    #     @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
    #     @test begin
    #         comt = mc.evolve(com,dt=0.01,tMax=4000,dtSave=200)

    #         all(abs.(comt[end].u .- comt[end].v) .< 0.1)
    #     end

    # end

    # for i in testplatforms #Check agents reading from medium

    #     m = @agent(1,
    #             p::Local,
    #             u::Medium,
    #             k::Global,

    #             Equation = begin
    #                 d_p = (u-k*p)*dt
    #             end
    #             )
    #     s = SimulationFree(m,
    #                     box=[(:x,0.,100.)],
    #                     medium=[MediumFlat("Dirichlet",100)])
    #     mc = compile(m,s,platform=i)
    #     com = Community(mc,N=100)
    #     com.x .= [Float64(i) for i in 0:1:99]
    #     com.u .= [Float64(i) for i in 0:99]
    #     com.k = 1.
    #     com.p .= 0.
    #     # println(mc.program)
    #     @test_nowarn mc.evolve(com,dt=0.1,tMax=.2)
    #     @test begin
    #         comt = mc.evolve(com,dt=0.01,tMax=4000,dtSave=200)
    #         all(comt[end].p .<= [Float64(i) for i in 0:1:99])
    #     end

    # end

end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
using JAC: HalfInt
using JAC.Basics: twice

@testset "HalfInt" begin
    @testset "Construction" begin
        @test HalfInt(3, 1) === HalfInt(6, 2)
        @test HalfInt(3) === HalfInt(3, 1)
        @test HalfInt(3.0, 1//1) === HalfInt(3)
        @test HalfInt(5, 2.0) === HalfInt(5, 2)
        @test HalfInt(0.5) === HalfInt(1, 2)
        @test HalfInt(1.0) === HalfInt(1, 1)
        @test HalfInt(-1//2) === HalfInt(-1, 2)
        @test HalfInt(AngularJ64(-3, 2)) === HalfInt(-3, 2)
        @test HalfInt(AngularJ64(4)) === HalfInt(4)
        @test HalfInt(AngularM64(-3, 2)) === HalfInt(-3, 2)
        @test HalfInt(AngularM64(4)) === HalfInt(4)
        @test_throws DomainError HalfInt(1, 3)
        @test_throws DomainError HalfInt(5.0, 2.5)
        @test_throws DomainError HalfInt(1//3)
        @test_throws InexactError HalfInt(0.8)
    end

    @testset "Conversion" begin
        # To real types
        for T in (:AbstractFloat, :Float16, :Float32, :Float64, :Rational, :(Rational{Int}), :(Rational{Int8}))
            @eval @test @inferred($T(HalfInt(5//2))) === $T(5//2)
        end
        @test @inferred(Number(HalfInt(-2))) === HalfInt(-2)
        @test @inferred(Real(HalfInt(3/2))) === HalfInt(3/2)
        @test BigFloat(HalfInt(4)) isa BigFloat
        @test BigFloat(HalfInt(4)) == BigFloat(4)
        @test BigFloat(HalfInt(-5//2)) == BigFloat(-2.5)
        # To integer types
        for T in (:Integer, :Signed, :Unsigned,
                  :Int8, :Int16, :Int32, :Int64, :Int128,
                  :UInt8, :UInt16, :UInt32, :UInt64, :UInt128)
            @eval @test @inferred($T(HalfInt(3))) === $T(3)
            @eval @test_throws InexactError $T(HalfInt(3/2))
        end
        @test BigInt(HalfInt(3)) isa BigInt
        @test BigInt(HalfInt(3)) == BigInt(3)
        @test_throws InexactError BigInt(HalfInt(3/2))
        @test Bool(HalfInt(0)) === false
        @test Bool(HalfInt(1)) === true
        @test_throws InexactError Bool(HalfInt(3/2))
        # To angular momentum types
        @test AngularJ64(HalfInt(1)) === AngularJ64(1, 1)
        @test AngularJ64(HalfInt(-5/2)) === AngularJ64(-5, 2)
        @test AngularM64(HalfInt(1)) === AngularM64(1, 1)
        @test AngularM64(HalfInt(-5/2)) === AngularM64(-5, 2)
        # Other functions
        @test @inferred(complex(HalfInt(-13//2))) === HalfInt(-13/2) + HalfInt(0)*im
        @test @inferred(real(HalfInt(-13//2))) === HalfInt(-13/2)
        @test @inferred(imag(HalfInt(-13//2))) === HalfInt(0)
        @test @inferred(float(HalfInt(-3//2))) === -1.5
    end

    @testset "Comparison" begin
        for T in (:Int8, :Int16, :Int32, :Int64, :Int128, :BigInt,
                  :UInt8, :UInt16, :UInt32, :UInt64, :UInt128)
            @eval @test HalfInt(3)   == $T(3)
            @eval @test HalfInt(5/2) != $T(3)
            @eval @test HalfInt(5/2) ≤ $T(3)
            @eval @test HalfInt(3)   ≤ $T(3)
            @eval @test HalfInt(5/2) < $T(3)
            @eval @test HalfInt(5/2) ≥ $T(2)
            @eval @test HalfInt(2)   ≥ $T(2)
            @eval @test HalfInt(5/2) > $T(2)
        end
        for T in (:HalfInt, :Float16, :Float32, :Float64, :BigFloat, :Rational)
            @eval @test HalfInt(5/2) == $T(5//2)
            @eval @test HalfInt(5/2) != $T(3)
            @eval @test HalfInt(5/2) ≤ $T(5//2)
            @eval @test HalfInt(5/2) ≤ $T(3)
            @eval @test HalfInt(5/2) < $T(3)
            @eval @test HalfInt(5/2) ≥ $T(5//2)
            @eval @test HalfInt(5/2) ≥ $T(2)
            @eval @test HalfInt(5/2) > $T(2)
        end
    end

    @testset "Hashing" begin
        @test hash(HalfInt(1)) === hash(1)
        @test hash(HalfInt(5//2)) === hash(5//2)
    end

    @testset "Arithmetic" begin
        # Unary plus, minus
        @test @inferred(+HalfInt(3/2)) === HalfInt(3/2)
        @test @inferred(-HalfInt(3/2)) === HalfInt(-3/2)
        # Addition, subtraction
        @test HalfInt(3/2) + HalfInt(5/2) === HalfInt(4)
        @test HalfInt(3/2) - HalfInt(5/2) === HalfInt(-1)
        @test HalfInt(5/2) + 2 === HalfInt(9/2)
        @test HalfInt(5/2) - 2 === HalfInt(1/2)
        @test 2 + HalfInt(5/2) === HalfInt(9/2)
        @test 2 - HalfInt(5/2) === HalfInt(-1/2)
        @test HalfInt(2) + 1.5 === 3.5
        @test HalfInt(2) - 1.5 === 0.5
        @test 1.5 + HalfInt(2) === 3.5
        @test 1.5 - HalfInt(2) === -0.5
        @test HalfInt(-1/2) + 4//3 === 5//6
        @test HalfInt(-1/2) - 4//3 === -11//6
        @test 4//3 + HalfInt(-1/2) === 5//6
        @test 4//3 - HalfInt(-1/2) === 11//6
        # Multiplication
        @test HalfInt(3/2)*HalfInt(7/2) === 5.25
        @test HalfInt(3/2)*5 === HalfInt(15/2)
        @test HalfInt(3/2)*5.0 === 7.5
        @test HalfInt(3/2)*(5//2) === 15//4
        @test 3*HalfInt(-1/2) === HalfInt(-3/2)
        @test 3.0*HalfInt(-1/2) === -1.5
        @test (3//2)*HalfInt(-1/2) === -3//4
        # Division
        @test HalfInt(5/2)/HalfInt(3/2) === 5/3
        @test HalfInt(5/2)/2.5 === 1.0
        @test HalfInt(5/2)/2 === 5/4
        @test HalfInt(5/2)/(2//1) === 5//4
        @test 2/HalfInt(3/2) === 4/3
        @test 2.0/HalfInt(3/2) === 4/3
        @test (2//1)/HalfInt(3/2) === 4//3
        # ^
        @test HalfInt(5/2)^2 == 6.25
        @test HalfInt(2)^3 == 8.0
        @test HalfInt(5/2)^2.0 ≈ 6.25
        @test HalfInt(2)^3.0 ≈ 8.0
        @test (-1)^HalfInt(2) ≈ 1
        @test 2^HalfInt(5/2) ≈ √32
        @test (-1.0)^HalfInt(2) ≈ 1
        @test (2.0)^HalfInt(5/2) ≈ √32
        @test HalfInt(5/2)^HalfInt(2) ≈ 6.25
        @test HalfInt(2)^HalfInt(3) ≈ 8.0
        # //
        @test HalfInt(3/2)//2 === 3//4
        @test 1//HalfInt(-1/2) === -2//1
        @test HalfInt(3/2)//HalfInt(5/2) === 3//5
        # floor
        @test floor(HalfInt(3/2)) === HalfInt(1)
        @test floor(HalfInt(0)) === HalfInt(0)
        @test floor(HalfInt(-1)) === HalfInt(-1)
        @test floor(HalfInt(-3/2)) === HalfInt(-2)
        # isinteger
        @test isinteger(HalfInt(2))
        @test !isinteger(HalfInt(3/2))
        # lcm
        @test lcm(1, HalfInt(1/2)) === HalfInt(1)
        @test lcm(1, HalfInt(3/2)) === HalfInt(3)
        @test lcm(2, HalfInt(3/2)) === HalfInt(6)
        @test lcm(HalfInt(1), HalfInt(1/2)) === HalfInt(1)
        @test lcm(HalfInt(1), HalfInt(3/2)) === HalfInt(3)
        @test lcm(HalfInt(2), HalfInt(3/2)) === HalfInt(6)
        @test lcm(HalfInt(3/2), HalfInt(3/2)) === HalfInt(3/2)
        @test lcm(HalfInt(5/2), HalfInt(3/2)) === HalfInt(15/2)
        # gcd
        @test gcd(1, HalfInt(3/2)) === HalfInt(1/2)
        @test gcd(2, HalfInt(3/2)) === HalfInt(1/2)
        @test gcd(3, HalfInt(3/2)) === HalfInt(3/2)
        @test gcd(3, HalfInt(9/2)) === HalfInt(3/2)
        @test gcd(HalfInt(1), HalfInt(3/2)) === HalfInt(1/2)
        @test gcd(HalfInt(2), HalfInt(3/2)) === HalfInt(1/2)
        @test gcd(HalfInt(3), HalfInt(3/2)) === HalfInt(3/2)
        @test gcd(HalfInt(3), HalfInt(9/2)) === HalfInt(3/2)
        @test gcd(HalfInt(3/2), HalfInt(3/2)) === HalfInt(3/2)
        @test gcd(HalfInt(-5/2), HalfInt(3/2)) === HalfInt(1/2)
        # gcdx
        @test ((d,u,v) = gcdx(1, HalfInt(3/2)); u*1+HalfInt(3/2)*v == d === HalfInt(1/2))
        @test ((d,u,v) = gcdx(2, HalfInt(3/2)); u*2+HalfInt(3/2)*v == d === HalfInt(1/2))
        @test ((d,u,v) = gcdx(3, HalfInt(3/2)); u*3+HalfInt(3/2)*v == d === HalfInt(3/2))
        @test ((d,u,v) = gcdx(3, HalfInt(9/2)); u*3+HalfInt(9/2)*v == d === HalfInt(3/2))
        @test ((d,u,v) = gcdx(HalfInt(1), HalfInt(3/2)); u*1+HalfInt(3/2)*v == d === HalfInt(1/2))
        @test ((d,u,v) = gcdx(HalfInt(2), HalfInt(3/2)); u*2+HalfInt(3/2)*v == d === HalfInt(1/2))
        @test ((d,u,v) = gcdx(HalfInt(3), HalfInt(3/2)); u*3+HalfInt(3/2)*v == d === HalfInt(3/2))
        @test ((d,u,v) = gcdx(HalfInt(3), HalfInt(9/2)); u*3+HalfInt(9/2)*v == d === HalfInt(3/2))
        @test ((d,u,v) = gcdx(HalfInt(3/2), HalfInt(3/2)); u*HalfInt(3/2)+HalfInt(3/2)*v == d === HalfInt(3/2))
        @test ((d,u,v) = gcdx(HalfInt(-5/2), HalfInt(3/2)); u*HalfInt(-5/2)+HalfInt(3/2)*v == d === HalfInt(1/2))
        # rem
        @test rem(HalfInt(5/2), 2) === HalfInt(1/2)
        @test rem(HalfInt(5/2), HalfInt(1/2)) === HalfInt(0)
        @test rem(HalfInt(5/2), HalfInt(7/2)) === HalfInt(5/2)
        @test rem(4, HalfInt(3/2)) === HalfInt(1)
        # sign
        @test sign(HalfInt(5/2)) == 1
        @test sign(HalfInt(1)) == 1
        @test sign(HalfInt(0)) == 0
        @test sign(HalfInt(-1/2)) == -1
    end

    @testset "Ranges" begin
        # UnitRange
        @test HalfInt(1/2):HalfInt(5) isa UnitRange{HalfInt}
        @test HalfInt(1/2):5 isa UnitRange{HalfInt}
        @test HalfInt(1/2):HalfInt(5) === HalfInt(1/2):HalfInt(9/2)
        @test length(HalfInt(1/2):HalfInt(0)) == 0
        @test length(HalfInt(1/2):HalfInt(1/2)) == 1
        @test length(HalfInt(1/2):HalfInt(1)) == 1
        @test length(HalfInt(-1/2):HalfInt(1)) == 2
        @test 2.5 ∉ HalfInt(1):HalfInt(5)
        @test 2.5 ∈ HalfInt(1/2):HalfInt(5)
        @test 5//2 ∉ HalfInt(1):HalfInt(5)
        @test 5//2 ∈ HalfInt(1/2):HalfInt(5)
        @test 2 ∈ HalfInt(1):HalfInt(5)
        @test 2 ∉ HalfInt(1/2):HalfInt(5)
        # StepRange
        @test HalfInt(2):HalfInt(1/2):HalfInt(5) isa StepRange{HalfInt}
        @test HalfInt(2):HalfInt(2):HalfInt(9/2) isa StepRange{HalfInt}
        @test HalfInt(2):2:HalfInt(9/2) isa StepRange{HalfInt}
        @test 2:2:HalfInt(9/2) isa StepRange{HalfInt}
        @test HalfInt(2):HalfInt(2):HalfInt(9/2) === HalfInt(2):HalfInt(2):HalfInt(4)
        @test HalfInt(2):HalfInt(-1):HalfInt(1/2) === HalfInt(2):HalfInt(-1):HalfInt(1)
        @test length(HalfInt(2):HalfInt(1/2):HalfInt(5)) == 7
        @test length(HalfInt(2):HalfInt(2):HalfInt(9/2)) == 2
        @test 2.5 ∈ HalfInt(1/2):HalfInt(2):HalfInt(7)
        @test 2.5 ∉ HalfInt(1/2):HalfInt(5/2):HalfInt(7)
        @test 2.5 ∉ HalfInt(1/2):HalfInt(-1):HalfInt(-1/2)
        @test 5//2 ∈ HalfInt(1/2):HalfInt(2):HalfInt(7)
        @test 5//2 ∈ HalfInt(1/2):2:HalfInt(7)
        @test 5//2 ∈ HalfInt(1/2):2:7
        @test 5//2 ∉ HalfInt(1/2):HalfInt(5/2):HalfInt(7)
        @test 5//2 ∉ HalfInt(1/2):HalfInt(-1):HalfInt(-1/2)
        @test 5//2 ∉ HalfInt(1/2):-1:HalfInt(-1/2)
        @test 2 ∈ HalfInt(0):HalfInt(1):HalfInt(5/2)
        @test 2 ∉ HalfInt(1/2):HalfInt(1):HalfInt(5/2)
        # intersect
        @test @inferred((HalfInt(1):HalfInt(3)) ∩ (HalfInt(1/2):HalfInt(3))) isa UnitRange{HalfInt}
        @test @inferred((1:3) ∩ (HalfInt(1/2):HalfInt(3))) isa UnitRange{HalfInt}
        @test @inferred((HalfInt(1/2):HalfInt(3)) ∩ (1:3)) isa UnitRange{HalfInt}
        @test @inferred((HalfInt(3/2):HalfInt(1):HalfInt(5)) ∩ (HalfInt(1):HalfInt(3))) isa StepRange{HalfInt}
        @test @inferred((HalfInt(3/2):HalfInt(1):HalfInt(5)) ∩ (1:3)) isa StepRange{HalfInt}
        @test @inferred((1:3) ∩ (HalfInt(3/2):HalfInt(1):HalfInt(5))) isa StepRange{HalfInt}
        @test @inferred((1:3:7) ∩ (HalfInt(3/2):HalfInt(5))) isa StepRange{HalfInt}
        @test @inferred((HalfInt(3/2):HalfInt(5)) ∩ (1:3:7)) isa StepRange{HalfInt}
        @test @inferred((HalfInt(3/2):HalfInt(1/2):HalfInt(5)) ∩ (HalfInt(-1):HalfInt(3):HalfInt(5))) isa StepRange{HalfInt}
        @test @inferred((-1:3:5) ∩ (HalfInt(3/2):HalfInt(1/2):HalfInt(5))) isa StepRange{HalfInt}
        @test @inferred((HalfInt(3/2):HalfInt(1/2):HalfInt(5)) ∩ (-1:3:5)) isa StepRange{HalfInt}
        @test isempty((HalfInt(1):HalfInt(3)) ∩ (HalfInt(1/2):HalfInt(3)))
        @test (HalfInt(1/2):HalfInt(2)) ∩ (HalfInt(3/2):HalfInt(3)) == HalfInt(3/2):HalfInt(3/2)
        @test isempty((HalfInt(1/2):HalfInt(2)) ∩ (1:2))
        @test isempty((1:2) ∩ (HalfInt(1/2):HalfInt(2)))
        @test (HalfInt(2):HalfInt(5)) ∩ (1:3) == HalfInt(2):HalfInt(3)
        @test (1:3) ∩ (HalfInt(2):HalfInt(5)) == HalfInt(2):HalfInt(3)
        @test (HalfInt(2):HalfInt(1):HalfInt(5)) ∩ (1:3) == HalfInt(2):HalfInt(1):HalfInt(3)
        @test (1:3) ∩ (HalfInt(2):HalfInt(1):HalfInt(5)) == HalfInt(2):HalfInt(1):HalfInt(3)
        @test isempty((HalfInt(3/2):HalfInt(1):HalfInt(5)) ∩ (1:3))
        @test isempty((1:3) ∩ (HalfInt(3/2):HalfInt(1):HalfInt(5)))
        @test (HalfInt(2):HalfInt(1/2):HalfInt(5)) ∩ (1:3) == HalfInt(2):HalfInt(1):HalfInt(3)
        @test (1:3) ∩ (HalfInt(2):HalfInt(1/2):HalfInt(5)) == HalfInt(2):HalfInt(1):HalfInt(3)
        @test (HalfInt(3/2):HalfInt(1/2):HalfInt(5)) ∩ (1:3) == HalfInt(2):HalfInt(1):HalfInt(3)
        @test (1:3) ∩ (HalfInt(3/2):HalfInt(1/2):HalfInt(5)) == HalfInt(2):HalfInt(1):HalfInt(3)
        @test (HalfInt(2):HalfInt(1):HalfInt(5)) ∩ (HalfInt(1):HalfInt(3)) == HalfInt(2):HalfInt(1):HalfInt(3)
        @test (HalfInt(1):HalfInt(3)) ∩ (HalfInt(2):HalfInt(1):HalfInt(5)) == HalfInt(2):HalfInt(1):HalfInt(3)
        @test isempty((HalfInt(3/2):HalfInt(1):HalfInt(5)) ∩ (HalfInt(1):HalfInt(3)))
        @test isempty((HalfInt(1):HalfInt(3)) ∩ (HalfInt(3/2):HalfInt(1):HalfInt(5)))
        @test (HalfInt(2):HalfInt(1/2):HalfInt(5)) ∩ (HalfInt(1):HalfInt(3)) == HalfInt(2):HalfInt(1):HalfInt(3)
        @test (HalfInt(1):HalfInt(3)) ∩ (HalfInt(2):HalfInt(1/2):HalfInt(5)) == HalfInt(2):HalfInt(1):HalfInt(3)
        @test (HalfInt(3/2):HalfInt(1/2):HalfInt(5)) ∩ (HalfInt(1):HalfInt(3)) == HalfInt(2):HalfInt(1):HalfInt(3)
        @test (HalfInt(1):HalfInt(3)) ∩ (HalfInt(3/2):HalfInt(1/2):HalfInt(5)) == HalfInt(2):HalfInt(1):HalfInt(3)
        @test (HalfInt(5/2):HalfInt(3/2):HalfInt(7)) ∩ (HalfInt(1):HalfInt(8)) == HalfInt(4):HalfInt(3):HalfInt(7)
        @test (HalfInt(1):HalfInt(8)) ∩ (HalfInt(5/2):HalfInt(3/2):HalfInt(7)) == HalfInt(4):HalfInt(3):HalfInt(7)
        @test (1:3:7) ∩ (HalfInt(3):HalfInt(5)) == HalfInt(4):HalfInt(3):HalfInt(4)
        @test (HalfInt(3):HalfInt(5)) ∩ (1:3:7) == HalfInt(4):HalfInt(3):HalfInt(4)
        @test isempty((1:3:7) ∩ (HalfInt(3/2):HalfInt(5)))
        @test isempty((HalfInt(3/2):HalfInt(5)) ∩ (1:3:7))
        @test (HalfInt(3/2):HalfInt(1/2):HalfInt(5)) ∩ (-1:3:5) == HalfInt(2):HalfInt(3):HalfInt(5)
        @test (-1:3:5) ∩ (HalfInt(3/2):HalfInt(1/2):HalfInt(5)) == HalfInt(2):HalfInt(3):HalfInt(5)
        @test (HalfInt(2):HalfInt(1):HalfInt(5)) ∩ (-1:3:5) == HalfInt(2):HalfInt(3):HalfInt(5)
        @test (-1:3:5) ∩ (HalfInt(2):HalfInt(1):HalfInt(5)) == HalfInt(2):HalfInt(3):HalfInt(5)
        @test isempty((HalfInt(3/2):HalfInt(1):HalfInt(5)) ∩ (-1:3:5))
        @test isempty((-1:3:5) ∩ (HalfInt(3/2):HalfInt(1):HalfInt(5)))
        @test isempty((HalfInt(3/2):HalfInt(1):HalfInt(5)) ∩ (HalfInt(-1):HalfInt(3):HalfInt(5)))
        @test isempty((HalfInt(-1):HalfInt(3):HalfInt(5)) ∩ (HalfInt(3/2):HalfInt(1):HalfInt(5)))
        @test (HalfInt(3/2):HalfInt(1):HalfInt(5)) ∩ (HalfInt(-1):HalfInt(3/2):HalfInt(5)) == HalfInt(7/2):HalfInt(3):HalfInt(7/2)
    end

    @testset "string" begin
        @test string(HalfInt(5)) == "5"
        @test string(HalfInt(-1)) == "-1"
        @test string(HalfInt(1/2)) == "1/2"
        @test string(HalfInt(-7/2)) == "-7/2"
    end

    @testset "twice" begin
        @test @inferred(twice(HalfInt(3/2))) === 3
        @test @inferred(twice(-2)) === -4
        @test @inferred(twice(1.6)) === 3.2
        @test @inferred(twice(2//3)) === 4//3
    end
end


text:
@testset "FeynmanAthermal tests" begin

println("FeynmanAthermal Tests - check athermal feynmanvw(α) -> v,w literature values.")

# Results from Feynman & Hibbs, Emended Edition, p 319.
# α v w E
Schultz=[
3.00 3.44 2.55 -3.1333
5.00 4.02 2.13 -5.4401
7.00 5.81 1.60 -8.1127
9.00 9.85 1.28 -11.486
11.0 15.5 1.15 -15.710
]

rows(M::Matrix) = map(x->reshape(getindex(M, x, :), :, size(M)[2]), 1:size(M)[1])

for r in rows(Schultz)
    α,vSchultz,wSchultz,ESchultz=r # Unpack row of results
    v,w=feynmanvw(α) # performs the optimisation
    E=F(v,w,α) # energy at the optimised parameters

    println("α=$α v=$v w=$w E=$E  | Schultz: v=$vSchultz w=$wSchultz E=$ESchultz")

    @test v ≈ vSchultz atol=0.1 # Strangely these need more tolerance than the Energies
    @test w ≈ wSchultz atol=0.1
    @test E ≈ ESchultz atol=0.001
end

end


text:
"`alwaysfalse(ags...; kwargs...) = false`"
alwaysfalse(ags...; kwargs...) = false

export WithinRange, NeighborNumber, SearchType
export searchstructure
export search, isearch, inrange, knn
export bulksearch, bulkisearch

"""
Supertype of all possible search types of the Neighborhood.jl common API.
"""
abstract type SearchType end

"""
    searchstructure(S, data, metric; kwargs...) → ss
Create a search structure `ss` of type `S` (e.g. `KDTree, BKTree, VPTree` etc.) based on the
given `data` and `metric`. The data types and supported metric types are package-specific,
but typical choices are subtypes of `<:Metric` from Distances.jl.
Some common metrics are re-exported by Neighborhood.jl.
"""
function searchstructure(::Type{S}, data::D, metric::M; kwargs...) where
                         {S, D, M}
    error("Given type $(S) has not implemented the Neighborhood.jl public API "*
          "for data type $(D) and metric type $(M).")
end

"""
    WithinRange(r::Real) <: SearchType
Search type representing all neighbors with distance `≤ r` from the query
(according to the search structure's metric).
"""
struct WithinRange{R} <: SearchType; r::R; end

"""
    NeighborNumber(k::Int) <: SearchType
Search type representing the `k` nearest neighbors of the query (or approximate
neighbors, depending on the search structure).
"""
struct NeighborNumber <: SearchType; k::Int; end

"""
    datatype(::Type{S}) :: Type

Get the type of the data points (arguments to the metric function) of a search
struture of type `S`.
"""
datatype(::Type) = Any
datatype(ss) = datatype(typeof(ss))

"""
    getmetric(ss)

Get the metric function used by the search structure `ss`.
"""
function getmetric end

"""
    search(ss, query, t::SearchType [, skip]; kwargs... ) → idxs, ds
Perform a neighbor search in the search structure `ss` for the given
`query` with search type `t` (see [`SearchType`](@ref)). Return the indices of the neighbors (in the original data)
and the distances from the query.

Optional `skip` function takes as input the index of the found neighbor
(in the original data) `skip(i)` and returns `true` if this neighbor should be skipped.

Package-specific keywords are possible.
"""
function search(ss::S, query::Q, t::T; kwargs...) where {S, Q, T<: SearchType}
    error("Given type $(S) has not implemented the Neighborhood.jl public API "*
          "for data type $(Q) and search type $(T).")
end
function search(ss::S, query::Q, t::T, skip; kwargs...) where {S, Q, T<: SearchType}
    error("Given type $(S) has not implemented the Neighborhood.jl public API "*
          "for data type $(Q), search type $(T) and skip function.")
end

"""
    isearch(args...; kwargs... ) → idxs
Same as [`search`](@ref) but only return the neighbor indices.
"""
isearch(args...; kwargs...) = search(args...; kwargs...)[1]

"""
    inrange(ss, query, r::Real [, skip]; kwargs...)
[`search`](@ref) for `WithinRange(r)` search type.
"""
inrange(a, b, r, args...; kwargs...) = search(a, b, WithinRange(r), args...; kwargs...)

"""
    knn(ss, query, k::Int [, skip]; kwargs...)
[`search`](@ref) for `NeighborNumber(k)` search type.
"""
knn(a, b, k::Integer, args...; kwargs...) =
search(a, b, NeighborNumber(k), args...; kwargs...)

###########################################################################################
# Bulk
###########################################################################################
"""
    bulksearch(ss, queries, t::SearchType [, skip]; kwargs... ) → vec_of_idxs, vec_of_ds
Same as [`search`](@ref) but many searches are done for many input query points.

In this case `skip` takes two arguments `skip(i, j)` where now `j` is simply
the index of the query that we are currently searching for (`j` is the index in
`queries` and goes from 1 to `length(queries)`).
"""
function bulksearch(ss, queries, t; kwargs...)
    i1, d1 = search(ss, queries[1], t; kwargs...)
    idxs, ds = [i1], [d1]
    sizehint!(idxs, length(queries))
    sizehint!(ds, length(queries))
    for j in 2:length(queries)
        i, d = search(ss, queries[j], t; kwargs...)
        push!(idxs, i); push!(ds, d)
    end
    return idxs, ds
end
function bulksearch(ss, queries, t, skip; kwargs...)
    i1, d1 = search(ss, queries[1], t, i -> skip(i, 1); kwargs...)
    idxs, ds = [i1], [d1]
    sizehint!(idxs, length(queries))
    sizehint!(ds, length(queries))
    for j in 2:length(queries)
        sk = k -> skip(j, k)
        i, d = search(ss, queries[j], t, sk; kwargs...)
        push!(idxs, i); push!(ds, d)
    end
    return idxs, ds
end

"""
    bulkisearch(ss, queries, t::SearchType [, skip]; kwargs... ) → vec_of_idxs
Same as [`bulksearch`](@ref) but return only the indices.
"""
bulkisearch(args...; kwargs...) = bulksearch(args...; kwargs...)[1]


text:
function checkIC(a::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Tri3},b::ReproducingKernel{SNode,𝒑,𝑠,𝜙,:Tri3}) where {𝒑,𝑠,𝜙}
    x₁ = a.𝓒[1].x;y₁ = a.𝓒[1].y
    x₂ = a.𝓒[2].x;y₂ = a.𝓒[2].y
    x₃ = a.𝓒[3].x;y₃ = a.𝓒[3].y
    n₁₁ = y₃-y₂;n₂₁ = y₁-y₃;n₃₁ = y₂-y₁
    n₁₂ = x₂-x₃;n₂₂ = x₃-x₁;n₃₂ = x₁-x₂
    nᵈ = length(a.𝓒)
    nᵖ = get𝑛𝒒(a)
    # nᵖ = get𝑛𝒑(a)
    f₁ = zeros(nᵈ,nᵖ)
    f₂ = zeros(nᵈ,nᵖ)
    for ξ in a.𝓖
        N,B₁,B₂ = get∇𝝭(a,ξ)
        𝒙 = get𝒙(a,ξ)
        𝒒 = get𝒒(a,ξ)
        p = get𝒑(a,𝒙)
        𝑤 = get𝑤(a,ξ)
        for i in 1:nᵈ
            for j in 1:nᵖ
                f₁[i,j] += B₁[i]*p[j]*𝑤
                f₂[i,j] += B₂[i]*p[j]*𝑤
            end
        end
    end
    for ξ in b.𝓖
        N = get𝝭(b,ξ)
        𝒒,∂𝒒∂ξ,∂𝒒∂η = get∇𝒒(b,ξ)
        𝒙 = get𝒙(b,ξ)
        p,∂p∂x,∂p∂y = get∇𝒑(b,𝒙)
        wᵇ = ξ.wᵇ
        w = ξ.w
        𝑤 = get𝑤(b,ξ)
        nᵇ₁ = 0.0;nᵇ₂ = 0.0
        ξ.ξ == 0.0 ? (nᵇ₁ += n₁₁;nᵇ₂ += n₁₂) : nothing
        ξ.η == 0.0 ? (nᵇ₁ += n₂₁;nᵇ₂ += n₂₂) : nothing
        ξ.ξ+ξ.η ≈ 1.0 ? (nᵇ₁ += n₃₁;nᵇ₂ += n₃₂) : nothing
        for j in 1:nᵖ
            W₁ = p[j]*nᵇ₁*wᵇ - ∂p∂x[j]*𝑤
            W₂ = p[j]*nᵇ₂*wᵇ - ∂p∂y[j]*𝑤
            for i in 1:nᵈ
                f₁[i,j] -= N[i]*W₁
                f₂[i,j] -= N[i]*W₂
            end
        end
    end
    return f₁,f₂
end

function checkCC(a::ReproducingKernel)
    nᵖ = get𝑛𝒑(a)
    nⁱ = length(a.𝓖)
    f = zeros(nⁱ,nᵖ)
    f₁ = zeros(nⁱ,nᵖ)
    f₂ = zeros(nⁱ,nᵖ)
    for i in 1:length(a.𝓖)
        ξ = a.𝓖[i]
        𝒙 = get𝒙(a,ξ)
        p,∂p∂x,∂p∂y = get∇𝒑(a,𝒙)
        N,B₁,B₂ = get∇𝝭(a,ξ)
        for j in 1:nᵖ
            for k in 1:length(a.𝓒)
                𝒙̄ = (a.𝓒[k].x,a.𝓒[k].y,a.𝓒[k].z)
                p̄ = get𝒑(a,𝒙̄)
                f[i,j]  += N[k]*p̄[j]
                f₁[i,j] += B₁[k]*p̄[j]
                f₂[i,j] += B₂[k]*p̄[j]
            end
            f[i,j]  -= p[j]
            f₁[i,j] -= ∂p∂x[j]
            f₂[i,j] -= ∂p∂y[j]
        end
    end
    return f,f₁,f₂
end

function checkConsistency(as::Vector{T}) where T<:ReproducingKernel
    nᵖ = get𝑛𝒑(as[1])
    f = zeros(nᵖ)
    𝒑ʰ = zeros(nᵖ)
    for a in as
        checkConsistency(a,f,𝒑ʰ)
    end
    return f.^0.5
end

function checkConsistency(a::T,f::Vector{Float64},𝒑ʰ::Vector{Float64}) where T<:ReproducingKernel
    for ξ in (a.𝓖)
        𝒙 = get𝒙(a,ξ)
        𝝭 = get𝝭(a,ξ)
        𝑤 = get𝑤(a,ξ)
        𝒑 = get𝒑(a,𝒙)
        fill!(𝒑ʰ,0.0)
        for (i,node) in enumerate(a.𝓒)
            𝒙ᵢ = (node.x,node.y,node.z)
            𝒑ᵢ = get𝒑(a,𝒙ᵢ)
            𝒑ʰ .+= 𝝭[i].*𝒑ᵢ
        end
        f .+= (𝒑 .- 𝒑ʰ).^2 .* 𝑤
    end
end

function checkConsistency(as::Vector{T},get𝝭_::F,get𝒑_::H) where {T<:ReproducingKernel,F<:Function,H<:Function}
    nᵖ = get𝑛𝒑(as[1])
    n = length(get𝒑_(as[1],(0.0,0.0,0.0)))
    f = zeros(nᵖ,n)
    𝒑 = zeros(nᵖ,n)
    𝒑ʰ = zeros(nᵖ,n)
    for a in as
        checkConsistency(a,get𝝭_,get𝒑_,f,𝒑,𝒑ʰ)
    end
    return f.^0.5
end

function checkConsistency(a::T,get𝝭_::F,get𝒑_::H,f::Matrix{Float64},𝒑::Matrix{Float64},𝒑ʰ::Matrix{Float64}) where {T<:ReproducingKernel,F<:Function,H<:Function}
    for ξ in (a.𝓖)
        𝒙 = get𝒙(a,ξ)
        𝑤 = get𝑤(a,ξ)
        𝝭_ = get𝝭_(a,ξ)
        𝒑_ = get𝒑_(a,𝒙)
        for (j,𝒑__) in enumerate(𝒑_)
            𝒑[:,j] .= 𝒑__
        end
        fill!(𝒑ʰ,0.0)
        for (i,node) in enumerate(a.𝓒)
            𝒙ᵢ = (node.x,node.y,node.z)
            𝒑ᵢ = get𝒑(a,𝒙ᵢ)
            for (j,𝝭) in enumerate(𝝭_)
                𝒑ʰ[:,j] .+= 𝝭[i].*𝒑ᵢ
            end
        end
        f .+= (𝒑 .- 𝒑ʰ).^2 .* 𝑤
    end
end

max_stars_repo_name: jmwjc/ApproxOperator

max_issues_repo_name: jmwjc/ApproxOperator

max_forks_repo_name: jmwjc/ApproxOperator


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

using BAT
using Compat.Test
using Distributions

@testset "const_density" begin
    cd = @inferred ConstDensity(
        BAT.HyperRectBounds{Float64}(HyperRectVolume([-1., 0.5], [2.,1]),
        [BAT.hard_bounds, BAT.hard_bounds]),
        one)

    @testset "ConstDensity" begin
        @test typeof(cd) <: ConstDensity{HyperRectBounds{Float64}, Float64}
        @test cd.log_value == zero(Float64)

        cd = @inferred ConstDensity(HyperRectBounds([-1., 0.5], [2.,1],
            hard_bounds), normalize)
        @test cd.log_value ≈ -0.4054651081081644

        pbounds = @inferred param_bounds(cd)
        @test pbounds.vol.lo ≈ [-1., 0.5]
        @test nparams(cd) == 2
    end

    @testset "unsafe_density_logval" begin
        @test BAT.unsafe_density_logval(cd, [1.,2.]) ≈ -0.4054651081081644
        res = ones(2)
        BAT.unsafe_density_logval!(res, cd, [1. 2.; 3. 4.],
            ExecContext())
        @test res ≈ -0.4054651081081644 * ones(2)
    end

    @testset "ExecCapabilities" begin
        ec = @inferred BAT.exec_capabilities(BAT.unsafe_density_logval, cd)
        @test ec.nthreads == 0
        @test ec.threadsafe == true
        @test ec.nprocs == 0
        @test ec.remotesafe == true

        ec = @inferred BAT.exec_capabilities(BAT.unsafe_density_logval!, ones(2) ,cd)
        @test ec.nthreads == 0
        @test ec.threadsafe == true
        @test ec.nprocs == 0
        @test ec.remotesafe == true
    end

    @testset "sampler" begin
        s = @inferred sampler(cd)
        @test typeof(s) <: SpatialVolume
        @test s.lo ≈ [-1., 0.5]
        @test s.hi ≈ [2., 1.]
    end
end


text:
struct Account
    dropletlimit::Integer
    floatingiplimit::Integer
    email::String
    uuid::String
    emailverified::Bool
    status::String
    statusmessage::String

    function Account(data::Dict{String})
        new(
            data["droplet_limit"],
            data["floating_ip_limit"],
            data["email"],
            data["uuid"],
            data["email_verified"],
            data["status"],
            data["status_message"]
        )
    end
end

function show(io::IO, a::Account)
    print(io, "Account ($(a.email))")
end

# Get user information
function getaccount!(client::AbstractClient)
    uri = joinpath(ENDPOINT, "account")
    Account(getdata!(client, uri))
end


text:
using DASSL
using Base.Test
using FactCheck

facts("Testing maxorder") do

    F(t,y,dy)=(dy+y.^2)
    Fy(t,y,dy)=diagm(2y)
    Fdy(t,y,dy)=eye(length(y))

    sol(t)=1./(1 .+ t)
    tspan=[0.0, 1.0]

    atol = 1.0e-5
    rtol = 1.0e-3

    # test the maxorder option
    for order=1:6
        # scalar version
        (tn,yn,dyn)=DASSL.dasslSolve(F, sol(0.0), tspan, maxorder = order)
        aerror = maximum(abs(yn-sol(tn)))
        rerror = maximum(abs(yn-sol(tn))/abs(sol(tn)))
        nsteps = length(tn)

        @fact aerror => less_than(2*nsteps*atol)
        @fact rerror => less_than(2*nsteps*rtol)

        # vector version
        (tnV,ynV,dynV)=DASSL.dasslSolve(F,[sol(0.0)], tspan, maxorder = order)

        @fact  vcat(ynV...) => yn
        @fact vcat(dynV...) => dyn

        # analytical jacobian version (vector)
        (tna,yna,dyna)=dasslSolve(F, [sol(0.0)], tspan, maxorder = order, Fy = Fy, Fdy = Fdy)
        aerror = maximum(abs(map(first,yn)-sol(tn)))
        rerror = maximum(abs(map(first,yn)-sol(tn))/abs(sol(tn)))
        nsteps = length(tn)

        @fact aerror => less_than(2*nsteps*atol)
        @fact rerror => less_than(2*nsteps*rtol)
    end
end

facts("Testing minimal error tolerances") do
    eps=1e-6
    # van der Pol equation
    Fvdp(t,y,dy)=([dy[1]+y[2],
                   eps*dy[2]-y[1]+y[2]*(1.0-y[1]^2)])
    y0 = [2., 0]
    tn = (0.,[2.,0.],[0.,-2/eps])

    tol = 1e0
    while true
        tol/=10
        sol=dasslIterator(Fvdp, y0, 1.0; reltol=tol, abstol=tol)

        try
            (tn,yn,dyn)=consume(sol)
        catch
            break
        end

    end
    @fact tol => roughly(1e-15)
end


text:
using Distributed
using JLD2

@time @everywhere include("../BioClim/src/required.jl")
cd("/Users/kiristern/Documents/GitHub/SDM")

@time begin
    # Load data from CSV files
    df = CSV.read("../data/pred_prey/Vulpes_vulpes_sp.csv", header=true, delim="\t")
    # Prepare data (select columns, arrange values)
    df = prepare_gbif_data(df)
    # Separate species
    taxa_occ = [df[df.species .== u,:] for u in unique(df.species)]

    # Define coordinates range
    lon_range = (-180.0, -46.0)
    lat_range = (-90.0, 90.0)
end

## Get the worldclim data
@time wc_vars = pmap(x -> worldclim(x)[lon_range, lat_range], 1:19);

## Make predictions for all species
@time all_predictions = pmap(x -> species_bclim(x, wc_vars), warblers_occ);

## Export predictions
@save "../data/all_predictions_vulpes_vulpes_world.jld2" all_predictions

max_stars_repo_name: kiristern/SDM

max_issues_repo_name: kiristern/SDM

max_forks_repo_name: kiristern/SDM


text:
"""
    filter_axes(f, a) -> Tuple

Return the axes of `a`, removing those for which `f` is false. The function `f`
is passed one argument.

```jldoctest
julia> axs = (Index{:a}(1:10), Index{:b}(1:10), Index(1:10));

julia> filter_axes(x -> dimnames(x) == :a, axs)
(Index{a}(1:10 => Base.OneTo(10)),)

julia> filter_axes(x -> length(x) == 10, axs)
(Index{a}(1:10 => Base.OneTo(10)), Index{b}(1:10 => Base.OneTo(10)), Index(1:10 => Base.OneTo(10)))
```
"""
filter_axes(f::Function, x::AbstractArray) = filter_axes(f, axes(x))
filter_axes(f::Function, x::Tuple) = _filter_axes(f, x)
function _filter_axes(f, x::Tuple)
    if f(first(x))
        return (first(x), _filter_axes(f, tail(x))...)
    else
        return _filter_axes(f, tail(x))
    end
end
_filter_axes(f, ::Tuple{}) = ()


text:
export mem, defcache
# # Memoization

"Default cache"
defcache() = Dict{defID(), Any}()

@inline tagmem(ω, cache = defcache()) = tag(ω, (memoize = cache))

# Don't ache traits
cache(f) = true

@inline function memapl(f, ω)
  if cache(f)
    result = get(ω.tags.cache, f.id, 0)
    if result === 0
      res = recurse(f, ω)
      ω.tags.cache[f.id] = res
      res
    else
      ω.tags.cache[f.id]::(Core.Compiler).return_type(f, typeof((ω,)))
    end
  else
    recurse(f, ω)
  end
end

"""
Memoize a variable.

`mem(x)` returns a 'memoized' version of `x` such that the result of calls to
`x(ω)` are cached.

```julia
using LinearAlgebra, BenchmarkTools
x = normal(0, 1, (1000, 1000))
h(x) = (println("call!"); svd(x).S)
y = lift(h)(x)
vars = randtuple((y, y*10, y*20))
@benchmark vars(ω) setup = (ω = rand(defΩ()))
@benchmark vars(ω) setup = (ω = tagmem(rand(defΩ())))
```
"""
@inline mem(f, cache = defcache()) = ω -> tagmem(ω, V)


text:
module LAP_julia

using TimerOutputs

export
    # visualisation funcitons
    showflow,
    imgshowflow,
    imgshow,
    warp_imgshowflow,
    showsparseflow,
    addpoints,
    imgoverlay,
    imgoverlay_v2,
    viewimg,
    # useful helpers
    mean,
    classic_alg,
    create_sparse_flow_from_sparse,
    create_sparse_flow_from_full,
    # data generation
    gen_quad_flow,
    gen_tiled_flow,
    gen_uniform_flow,
    gen_homo_flow,
    gen_chess,
    gen_init,
    gen_lena,
    gen_anhir,
    gen_spaghetti,
    load_anhir_image_pair,
    gen_anhir,
    resize_to_diag_size,
    # Main LAP funcitons
    lap,
    sparse_lap,
    sparse_lap_win_sum1,
    sparse_pflap,
    single_lap,
    pflap,
    sparse_pflap,
    single_lap_at_points,
    # Experimental LAP functions
    sparse_pflap_psnr,
    # Types
    Image,
    Flow,
    # Point location
    find_edge_points,
    # Interpolation
    warp_img,
    interpolate_flow,
    interpolate_flow_quad,
    # experimenting
    test_registration_alg,
    assess_flow_quality,
    assess_source_reg_quality,
    time_reg_alg,
    sparse_pflap_multi

"""
    Image{T} = Matrix{T} where T

Image is a `Matrix` with elements that are `Real`.
"""
const Image{T} = Matrix{T} where T

"""
    Flow{T} = Matrix{Complex{T}} where T <: Real

Flow is a `Matrix` with elements that are `Complex`.
"""
const Flow{T} = Matrix{Complex{T}} where T <: Real

include("helpers.jl")
include("inpaint.jl")
include("interpolation.jl")
include("visualise.jl")
include("data_gen.jl")
include("gradient_points.jl")
include("lap_algs/helpers.jl")
include("lap_algs/multi.jl")
include("lap_algs/single.jl")
include("lap_algs/experimental.jl")
include("experimenting/speedtest.jl")
include("experimenting/qualitytest.jl")
include("experimenting/visualise.jl")
include("birl.jl")

loaded() = println("LAP_julia succesfully loaded!")

loaded()

end # module


text:
# initialize outermost container Dash app
include("app.jl");

using Pkg
Pkg.develop(path="./dash-user-guide-components")

using Dash, DashCoreComponents, DashHtmlComponents, DashUserGuideComponents, Match

# Load Chapter, Example, Header, Section, Syntax components
map(include, filter(x->occursin(r".jl$", x), readdir("dash_docs/reusable_components/", join=true)));

# Load chapter container Dash apps
include("dash_docs/chapters/whats_dash/introduction.jl");
include("dash_docs/chapters/installation/index.jl");
include("dash_docs/chapters/getting_started/index.jl");
include("dash_docs/chapters/basic_callbacks/index.jl");
include("dash_docs/chapters/graph_crossfiltering/index.jl");
include("dash_docs/chapters/sharing_data/index.jl");
include("dash_docs/chapters/faq_gotchas/index.jl");
include("dash_docs/chapters/deployment/index.jl");

for example in chapters_callbacks.examples
    example.callback!(app)
end

for example in chapters_interactive_graphing.examples
    example.callback!(app)
end

for example in chapters_sharing_data.examples
    example.callback!(app)
end

header = html_div(
    children = (
        html_div(
            style = Dict("height" => "95%"),
            className = "container-width",
            children = (
                html_div(
                    children = (
                        html_span("📣 Welcome! Dash for Julia is New! Read "),
                        html_a("the community announcement", className = "link", href="https://community.plotly.com/t/welcome-to-dash-julia/46056"),
                        html_span(" regarding the current status of the project.")
                    ),
                    style = Dict("background-color" => "#80CFBE", "text-align" => "center", "color" => "#000000")),
                html_a(
                    html_img(
                        style = Dict("height" => "100%", "padding-left" => "80px"),
                        src = "https://dash.plotly.com/assets/images/logo-plotly.png"
                        ),
                        href = "https://plotly.com/products/dash",
                        className = "logo-link"
                        ),
                        html_div(
                            children = (
                                html_a("Announcements", className = "links--announcements", href = "https://community.plotly.com/tag/announcements"),
                                html_a("Gallery", className = "links--gallery", href = "https://dash-gallery.plotly.host"),
                                html_a("Show & Tell", className = "links--show-and-tell", href = "https://community.plotly.com/tag/show-and-tell"),
                                html_a("Community Forum", className = "links--community-forum", href = "https://community.plotly.com/c/dash/julia/20"),
                                html_iframe(src="https://ghbtns.com/github-btn.html?user=plotly&repo=dash.jl&type=star&count=true&size=small", style = Dict("border" => "none", "height" => "30px", "verticalAlign" => "middle", "marginTop" => "9px", "width" => "120px")),
                                html_a("dash enterprise demo", className = "link--demo-button", href = "https://plotly.com/get-demo/?utm_source=docs&utm_medium=banner&utm_campaign=sept&utm_content=demo", style = Dict("backgroundColor" => "#f4564e", "borderRadius" => "1.22rem", "color" => "white", "cursor" => "pointer", "display" => "inline-block", "fontStyle" => "italic", "fontWeight" => "700", "lineHeight" => "1.2", "letter-spacing" => "1.33px", "outline" => "none", "padding" => ".55rem 1.22rem", "marginRight" => "5px", "textAlign" => "center", "verticalAlign" => "middle", "textDecoration" => "none", "textTransform" => "uppercase", "transition" => "background-color .2s ease-in-out"))
                                ),
                            className = "links"
                        )
                    )
        )
    ),
    className = "header"
);

app.layout = html_div() do
    html_div(id = "wait-for-layout"),
    dcc_location(id = "url", refresh=false),
    header,
    html_div(
        className = "content-wrapper",
        children = (
            html_div(
                (
                    html_div(id = "backlinks-top", className = "backlinks"),
                    html_div(
                        html_div(id = "chapter", className = "content"), # the children of this component is the layout of a dash app, based on URL
                        className = "content-container",
                        style = Dict("margin" => "70px")
                    ),
                    html_div(id = "backlinks-bottom", className = "backlinks")
                ),
                className = "rhs-content container-width"
            ),
            dugc_pagemenu(id = "pagemenu")
        )
    )
end;

callback!(app,
    Output("chapter", "children"),
    Output("pagemenu", "dummy2"),
    Input("url", "pathname")) do pathname
       get_content(pathname) = @match pathname begin
            "/introduction" => chapters_whats_dash.app.layout
            "/installation" => chapters_installation.app.layout
            "/getting-started" => chapters_getting_started.app.layout
            "/basic-callbacks" => chapters_callbacks.app.layout
            "/interactive-graphing" => chapters_interactive_graphing.app.layout
            "/sharing-data-between-callbacks" => chapters_sharing_data.app.layout
            "/deployment" => chapters_deployment.app.layout
            "/faqs" => chapters_faq_gotchas.app.layout
            _ => html_div() do
                html_br(),
                html_h1("Dash for Julia User Guide"),
                Section(
                    "What's Dash?",
                    (
                        Chapter(
                            "Introduction",
                            "/introduction",
                            "A quick paragraph about Dash and a link to the talk at Plotcon that started it all."
                        ),
                        Chapter(
                            "Announcement Essay",
                            "https://medium.com/plotly/announcing-dash-for-julia-f017c90c6ef1",
                            "Our extended essay on Dash. An extended discussion of Dash's architecture and our motivation behind the project."
                        ),
                        Chapter(
                            "Dash App Gallery",
                            "https://dash.plotly.com/gallery",
                            "A glimpse into what's possible with Dash."
                        ),
                        Chapter(
                            "Dash Club",
                            "https://plot.us12.list-manage.com/subscribe?u=28d7f8f0685d044fb51f0d4ee&id=0c1cb734d7",
                            "A fortnightly email newsletter by chriddyp, the creator of Dash."
                        )
                    )
                ),
                Section(
                "Dash Tutorial",
                (
                    Chapter(
                        "Part 1. Installation",
                        "/installation",
                        "How to install and upgrade Dash libraries with the Pkg package manager."
                    ),
                    Chapter(
                        "Part 2. The Dash Layout",
                        "/getting-started",
                        "The Dash `layout` describes what your app will look like and is composed of a set of declarative Dash components."
                    ),
                    Chapter(
                        "Part 3. Basic Callbacks",
                        "/basic-callbacks",
                        "Dash apps are made interactive through Dash Callbacks:
                        Julia functions that are automatically called whenever an input component's property changes. Callbacks can be chained,
                        allowing one update in the UI to trigger several updates across the app."
                    ),
                    Chapter(
                        "Part 4. Interactive Graphing and Crossfiltering",
                        "/interactive-graphing",
                        "Bind interactivity to the Dash `Graph` component whenever you hover, click, or
                        select points on your chart."
                    ),
                    Chapter(
                        "Part 5. Sharing Data Between Callbacks",
                        "/sharing-data-between-callbacks",
                        "`global` variables will break your Dash apps. However, there are other ways
                        to share data between callbacks. This chapter is useful for callbacks that
                        run expensive data processing tasks or process large data."
                    ),
                    Chapter(
                        "Part 6. FAQs and Gotchas",
                        "/faqs",
                        "If you have read through the rest of the tutorial and still have questions
                        or are encountering unexpected behaviour, this chapter may be useful."
                        )
                    )
                ),
                Section(
                    "Production",
                    (
                        Chapter(
                            "Deployment",
                            "/deployment",
                            ""
                        )
                    )
                )
            end
        end
    return get_content(pathname), ""
end;

callback!(
    ClientsideFunction("clientside", "pagemenu"),
    app,
    Output("pagemenu", "dummy"),
    Input("chapter", "children")
)

port = parse(Int64, ENV["PORT"])

print("Binding to PORT $(port)...")

run_server(app, "0.0.0.0", port)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
export check_grad, nparams
export gradient_numeric

using FixedPointNumbers: Fixed

isvar(x) = nparams(x) != 0

nparams(model) = nparams(NiLangCore.struct2namedtuple(model))
nparams(x::AbstractArray{<:AbstractFloat}) = length(x)
nparams(x::AbstractArray{<:GVar}) = length(x)
nparams(x::AbstractArray) = sum(nparams, x)
nparams(x::Fixed) = 1
function nparams(x::Union{Tuple,NamedTuple})
    res = 0
    for xi in x
        res += nparams(xi)
    end
    res
end
nparams(x::AbstractFloat) = 1
nparams(x::GVar) = 1

function tset(vfunc::Function, tp::Tuple, iloss)
    map(i->i===iloss ? vfunc(tp[i]) : tp[i], (1:length(tp)...,))
end
function tset(value, tp::Tuple, iloss)
    map(i->i===iloss ? value : tp[i], (1:length(tp)...,))
end

function update_var(args, iarg, i::Int, val)
    args[iarg][i] += val
    args
end

function update_var(args, iarg, ::Nothing, val)
    tset(x->chfield(x, value, value(x) + val), args, iarg)
end

function ng_single(::Type{T}, f, args, kwargs, iarg, i, iloss, δ) where T
    args = update_var(args, iarg, i, T(δ/2))
    @instr f(args...; kwargs...)
    pos = value(args[iloss])
    @instr (~f)(args...; kwargs...)
    args = update_var(args, iarg, i, -T(δ))
    @instr f(args...; kwargs...)
    neg = value(args[iloss])
    @instr (~f)(args...; kwargs...)
    args = update_var(args, iarg, i, T(δ/2))
    (pos - neg)/δ
end

function ng_single(::Type{T}, f, args, kwargs, iarg, i, iloss, δ) where T<:Complex
    res = zero(T)
    for dd = [δ, im*δ]
        args = update_var(args, iarg, i, dd/2)
        @instr f(args...; kwargs...)
        pos = value(args[iloss])
        @instr (~f)(args...; kwargs...)
        args = update_var(args, iarg, i, -dd)
        @instr f(args...; kwargs...)
        neg = value(args[iloss])
        @instr (~f)(args...; kwargs...)
        args = update_var(args, iarg, i, dd/2)
        if dd == δ
            res += (pos - neg)/δ
        else
            res += im*(pos - neg)/δ
        end
    end
    res
end

function ng(f, args, iarg; iloss::Int, δ=1e-5, kwargs...)
    x = args[iarg]
    T = eltype(x)
    if x isa AbstractArray
        res = zero(x)
        for i = 1:length(x)
            res[i] = ng_single(T, f, args, kwargs, iarg, i, iloss, δ)
        end
        return res
    else
        ng_single(T, f, args, kwargs, iarg, nothing, iloss, δ)
    end
end

"""
    gradient_numeric(f, args...; iloss, kwargs...)

Numeric differentiating f(args..., kwargs...).
"""
function gradient_numeric(f, args; iloss::Int, kwargs...)
    map(1:length(args)) do iarg
        if isvar(args[iarg])
            ng(f, args, iarg; iloss=iloss, kwargs...)
        else
            0
        end
    end
end

"""
    check_grad(f, args; atol::Real=1e-8, verbose::Bool=false, iloss::Int, kwargs...)

Return true if the gradient of `f(args..., kwargs...)` is reversible.
"""
function check_grad(f, args; atol::Real=1e-4, verbose::Bool=false, iloss::Int, kwargs...)
    vars = ((iarg for iarg in 1:length(args) if isvar(args[iarg]))...,)
    initial_vars = deepcopy(vars)
    ngs = gradient_numeric(f, args; kwargs..., iloss=iloss)
    gs = gradient(f, args; kwargs..., iloss=iloss)
    verbose && @show ngs
    verbose && @show gs
    if !all(isapprox.(ngs, gs, atol=atol))
        verbose && println("gradient not match: $ngs v.s. $gs")
        return false
    end

    if !world_similar(initial_vars, vars, atol=atol, verbose=verbose)
        verbose && println("world changed during obtaining gradient.")
        return false
    end
    return true
end