text:
using DynamicalBilliards, AbstractPlotting, MakieLayout
export interactive_billiard, interactive_billiard_bmap, billiard_video

"""
    interactive_billiard(bd::Billiard [, x, y, φ] [, ω=nothing]; kwargs...)
    interactive_billiard(bd::Billiard, ps::Vector{<:AbstractParticle}; kwargs...)

Launch an interactive application that evolves particles in a dynamical billiard `bd`, using
[DynamicalBilliards.jl](https://juliadynamics.github.io/DynamicalBilliards.jl/dev/).
You can either specify exactly the particles that will be used `ps` or provide
some initial conditions `x,y,φ,ω`, which by default are random in the billiard.

The particles are evolved in real time instead of being pre-calculated,
so the application can be left to run for infinite time.

See also [`interactive_billiard_bmap`](@ref) and [`billiard_video`](@ref).

## Interaction
Push "play" to start evolving particles in the billiard, and "reset" to restore them
to their (latest) initial condition. The "particles" hides or shows the particles.
The "speed" slider controls the animation speed (in fact, it controls how often
are the plots updated).

Clicking and dragging inside the billiard plot shows a line. When the line is selected,
new particles are created that have the direction of this line, as well as its starting
position, using the function `particlebeam` from `DynamicalBilliards`.

## Further keywords
* `N = 100` : if exact particles are not given, `N` are created. Otherwise it is `length(ps)`.
* `dx = 0.01` : width of the particle beam.
* `dt = 0.001` : time resolution of the animation.
* `tail = 1000` : length of the tail of the particles (multiplies `dt`).
* `colors = JULIADYNAMICS_COLORS` : If a symbol (colormap name) each particle gets
  a color from the map. If Vector of length `N`, each particle gets a color form the vector.
  If Vector with length < `N`, linear interpolation across contained colors is done.
* `fade = true` : Whether to add fadeout to the particle tail.
* `sleept = nothing` : If the slowest speed of the animation is already too fast,
  give a small number to `sleept`.
* `plot_particles = true` : If false, the particles are not plotted (as balls and arrows).
  This makes the application faster (you cannot show them again with the button).
"""
interactive_billiard(bd::Billiard, ω::Union{Nothing, Real} = nothing; kwargs...) =
interactive_billiard(bd::Billiard, randominside_xyφ(bd)..., ω; kwargs...)

function interactive_billiard(bd::Billiard, x::Real, y::Real, φ::Real, ω = nothing;
    kwargs...)
    N = get(kwargs, :N, 100)
    dx = get(kwargs, :dx, 0.01)
    ps = particlebeam(x, y, φ, N, dx, ω, Float32)
    interactive_billiard(bd::Billiard, ps; kwargs...)
end

function interactive_billiard(bd::Billiard, ps::Vector{<:AbstractParticle};
        dt = 0.001, tail = 1000, dx = 0.01, colors = JULIADYNAMICS_COLORS,
        plot_particles = true, α = 1.0, N = 100, res = (800, 800),
        intervals = nothing, sleept = nothing, fade = true,
        backgroundcolor = DEFAULT_BG,
        vr = _estimate_vr(bd),
        add_controls = true,
    )

    if eltype(bd) ≠ Float32 || eltype(ps[1]) ≠ Float32
        error("Only Float32 number type is possible for the billiard applications. "*
        "Please initialize billiards and particles by explicitly passing Float32 numbers "*
        "in all numeric fields (e.g. `bd = billiard_mushroom(1f0, 0.2f0, 1f0, 0f0)`)")
    end
    N = length(ps)
    p0s = deepcopy(ps) # deep is necessary because vector of mutables
    ω0 = ismagnetic(ps[1]) ? ps[1].ω : nothing

    # Initialized inside process
    cs = (!(colors isa Vector) || length(colors) ≠ N) ? colors_from_map(colors, α, N) : colors
    scene, layout = layoutscene(resolution = res, backgroundcolor = backgroundcolor)
    ax = layout[1, 1] = LAxis(scene, backgroundcolor = backgroundcolor)
    tight_ticklabel_spacing!(ax)
    ax.autolimitaspect = 1
    allparobs = [ParObs(p, bd, tail) for p in ps]
    bdplot!(ax, bd)

    # Plot tails
    for (i, p) in enumerate(allparobs)
        x = to_color(cs[i])
        if fade
            x = [RGBAf0(x.r, x.g, x.b, i/tail) for i in 1:tail]
        end
        lines!(ax, p.tail, color = x)
    end

    # Plot particles
    if plot_particles
        balls = Observable([Point2f0(p.p.pos) for p in allparobs])
        vels = Observable([vr * Point2f0(p.p.vel) for p in allparobs])
        particle_plots = (
            scatter!(ax, balls; color = darken_color.(cs), marker = MARKER, markersize = 12AbstractPlotting.px),
            arrows!(ax, balls, vels; arrowcolor = darken_color.(cs), linecolor = darken_color.(cs),
                normalize = false, arrowsize = 0.02AbstractPlotting.px,
                linewidth  = 4,
            )
        )
    end

    # Controls:
    if add_controls
        resetbutton = LButton(scene, label = "reset",
            buttoncolor = RGBf0(0.8, 0.8, 0.8),
            height = 40, width = 80
        )
        runbutton = LButton(scene, label = Observable("run"),
            buttoncolor = Observable(RGBf0(0.8, 0.8, 0.8)),
            buttoncolor_hover = Observable(RGBf0(0.7, 0.7, 0.9)),
            buttoncolor_active = Observable(RGBf0(0.6, 0.6, 1.0)),
            labelcolor = Observable((RGBf0(0,0,0))),
            labelcolor_active = Observable((RGBf0(1,1,1))),
            height = 40, width = 70,
        )
        nslider = LSlider(scene, range = 0:50, startvalue=0)
        controls = [resetbutton, runbutton, LText(scene, "speed:"), nslider]
        if plot_particles
            particlebutton = LButton(scene, label = "particles",
                buttoncolor = RGBf0(0.8, 0.8, 0.8),
                height = 40, width = 100
            )
            pushfirst!(controls, particlebutton)
        end

        # Functionality that CREATES the play/stop button
        # TODO: will be deleted once MakieLayout has proper togglable button
        on(runbutton.clicks) do n
            t = runbutton.label[] == "run" ? "stop" : "run"
            runbutton.label[] = t
            for (s1, s2) in ((:buttoncolor, :buttoncolor_active), (:labelcolor, :labelcolor_active))
            getproperty(runbutton, s1)[], getproperty(runbutton, s2)[] =
                getproperty(runbutton, s2)[], getproperty(runbutton, s1)[]
            end
            runbutton.labelcolor_hover[] = runbutton.labelcolor[]
        end

        # Create the "control panel"
        layout[2, 1] = grid!(hcat(
            controls...,
        ), tellwidth = false, tellheight = true)

        # Play/stop
        on(runbutton.clicks) do nclicks
            @async while runbutton.label[] == "stop"
                for i in 1:N
                    p = ps[i]
                    parobs = allparobs[i]
                    animstep!(parobs, bd, dt, true, intervals)
                    if plot_particles
                        balls[][i] = parobs.p.pos
                        vels[][i] = vr * parobs.p.vel
                    end
                end
                if plot_particles
                    balls[] = balls[]
                    vels[] = vels[]
                end
                for _ in 1:nslider.value[]
                    for i in 1:N
                        p = ps[i]
                        parobs = allparobs[i]
                        animstep!(parobs, bd, dt, false, intervals)
                    end
                end
                if sleept == nothing
                    yield()
                else
                    sleep(sleept)
                end
                isopen(scene) || break
            end
        end

        # Resetting
        on(resetbutton.clicks) do nclicks
            # update parobs
            for i in 1:N
                parobs = allparobs[i]
                rebind_partobs!(parobs, p0s[i], bd)
                if plot_particles
                    balls[][i] = parobs.p.pos
                    vels[][i] = vr * parobs.p.vel
                end
            end
            if plot_particles
                balls[] = balls[]
                vels[] = vels[]
            end
            yield()
        end

        # Show/hide particles
        if plot_particles
            on(particlebutton.clicks) do nclicks
                particle_plots[1].visible[] = !particle_plots[1].visible[]
                for i in 1:2
                particle_plots[2].plots[i].visible[] = !particle_plots[2].plots[i].visible[]
                end
            end
        end

        # Selecting new particles
        sline = select_line(ax.scene)
        on(sline) do val
            pos = val[1]
            dir = val[2] - val[1]
            φ = atan(dir[2], dir[1])
            p0s .= particlebeam(pos..., φ, N, dx, ω0, Float32)
            resetbutton.clicks[] += 1
        end

    end # adding controls blocks

    display(scene)
    if add_controls
        return scene, layout, allparobs, resetbutton, p0s, sline
    elseif plot_particles
        return scene, layout, allparobs, balls, vels, vr
    else
        return scene, layout, allparobs, nothing, nothing, nothing
    end
end

"""
    billiard(file, bd::Billiard [, x, y, φ] [, ω=nothing]; kwargs...)
    billiard(file, bd::Billiard, ps::Vector{<:AbstractParticle}; kwargs...)

Perform the same animation like in [`interactive_billiard`](@ref), but there is no
interaction; the result is saved directly as a video in `file` (no buttons are shown).

## Keywords
* `N, dt, tail, dx, colors, plot_particles, fade`: same as `interactive_billiard`, but
  with a bit "denser" defaults. `plot_particles` is `false` by default here.
* `speed = 4`: Animation "speed" (how many `dt` steps are taken before a frame is recorded)
* `frames = 1000`: amount of frames to record.
* `framerate = 60`: exported framerate.
* `backgroundcolor = RGBf0(0.99, 0.99, 0.99)`.
* `res = nothing`: resolution of the frames. If nothing, a resolution matching the
  the billiard aspect ratio is estimated. Otherwise pass a 2-tuple.
"""
billiard_video(file::String, bd::Billiard, ω::Union{Nothing, Real} = nothing; kwargs...) =
billiard_video(file::String, bd::Billiard, randominside_xyφ(bd)..., ω; kwargs...)

function billiard_video(file::String, bd::Billiard, x::Real, y::Real, φ::Real, ω = nothing;
    kwargs...)
    N = get(kwargs, :N, 500)
    dx = get(kwargs, :dx, 0.01)
    ps = particlebeam(x, y, φ, N, dx, ω, Float32)
    billiard_video(file, bd, ps; kwargs...)
end

function billiard_video(file::String, bd::Billiard, ps::Vector{<:AbstractParticle};
        plot_particles = false, res = nothing, dt = 0.001,
        speed = 4, frames = 1000, framerate = 60, kwargs...
    )

    if res == nothing
        xmin, ymin, xmax, ymax = DynamicalBilliards.cellsize(bd)
        aspect = (xmax - xmin)/(ymax-ymin)
        res = (round(Int, aspect*800), 800)
    end

    scene, layout, allparobs, balls, vels, vr = interactive_billiard(bd, ps;
        res = res, plot_particles=plot_particles, kwargs..., add_controls =false
    )
    N = length(ps)

    record(scene, file, 1:frames; framerate = framerate) do j
        for i in 1:N
            p = ps[i]
            parobs = allparobs[i]
            animstep!(parobs, bd, dt, true)
            if plot_particles
                balls[][i] = parobs.p.pos
                vels[][i] = vr * parobs.p.vel
            end
        end
        if plot_particles
            balls[] = balls[]
            vels[] = vels[]
        end
        for _ in 1:speed
            for i in 1:N
                p = ps[i]
                parobs = allparobs[i]
                animstep!(parobs, bd, dt, false)
            end
        end
    end
    return
end

"""
    interactive_billiard_bmap(bd::Billiard, ω=nothing; kwargs...)
Launch an interactive application whose left part is [`interactive_billiard`](@ref)
and whose write part is an interactive boundary map of the billiard (see "Phase spaces"
in DynamicalBilliards.jl).

A particle evolved in the real billiard is also shown on the boundary map.
All interaction of the billiard works as before, but there is also interaction in
the boundary map: clicking on it will generate a particle whose boundary map is
the clicked point.

The mean collision time "m.c.t." of the particle is shown as well.

## Keywords
* `newcolor = randomcolor` A function which takes as input `(pos, vel, ξ, sφ)` and
  outputs a color (for the scatter points in the boundary map).
* `ms = 12` markersize (in pixels).
* `dt, tail, sleept, fade` : propagated to `interactive_billiard`.
"""
function interactive_billiard_bmap(bd::Billiard, ω=nothing;
    newcolor = randomcolor, ms = 12, lock = true,
    kwargs...)

    intervals = arcintervals(bd)
    scene, layout, allparobs, resetbutton, p0s, sline = interactive_billiard(
        bd, ω; kwargs..., N = 1, intervals = intervals, res = (1600, 800)
    )
    parobs = allparobs[1] # only one exists.

    sublayout = GridLayout()
    cleanbutton = LButton(scene, label = "clean", tellwidth = false)
    sublayout[2, 1] = cleanbutton
    mct = Observable("m.c.t. = 0.0")
    mcttext = LText(scene, mct, haligh = :left, tellwidth = false)
    sublayout[2, 2] = mcttext
    bmapax = sublayout[1,:] = LAxis(scene)
    bmapax.xlabel = "arclength ξ"
    bmapax.ylabel = "normal angle sin(θ)"
    bmapax.targetlimits[] = BBox(intervals[1], intervals[end], -1, 1)

    current_color = Observable(newcolor(parobs.p.pos, parobs.p.vel, parobs.ξsin[]...))
    scatter_points = Observable(Point2f0[])
    scatter_colors = Observable(RGBAf0[])

    scatter!(bmapax.scene, scatter_points; color = scatter_colors,
        marker = MARKER, markersize = ms*AbstractPlotting.px
    )

    # Make obstacle axis, add info about where each obstacle terminates
    if length(intervals) > 2 # at least 2 obstacles
        add_obstacle_axis!(scene, sublayout, intervals, bmapax, lock)
    end

    # Obtain new color when selecting line in main plot
    on(sline) do val
        current_color[] = newcolor(parobs.p.pos, parobs.p.vel, parobs.ξsin[]...)
    end
    # Whenever boundary map is updated, plot the update
    on(parobs.ξsin) do v
        pushupdate!(scatter_points, v)
        pushupdate!(scatter_colors, current_color[])
        τ = parobs.T/parobs.n
        mct[] = "m.c.t. = $(rpad(string(τ), 25))"
    end

    # Selecting a point on the boundary map:
    spoint = select_point(bmapax.scene)
    on(spoint) do ξsin
        ξ, sφ = ξsin
        sφ = clamp(sφ, -1, 1)
        ξsin = Point2f0(ξ, sφ)
        if !(intervals[1] ≤ ξ ≤ intervals[end])
            error("selected point is outside the boundary map")
        end
        pos, vel = from_bcoords(ξsin..., bd, intervals)
        current_color[] = newcolor(pos, vel, ξsin...)
        p0s[1] = ω ≠ nothing ? MagneticParticle(pos, vel, ω) : Particle(pos, vel)
        propagate!(p0s[1], 10eps(Float32)) # ensure no bad behavior with boundary
        rebind_partobs!(parobs, p0s[1], bd, ξsin)
        resetbutton.clicks[] += 1 # reset main billiard
    end

    # Clean button functionality
    on(cleanbutton.clicks) do nclicks
        scatter_points[] = Point2f0[]
        scatter_colors[] = RGBAf0[]
    end

    # Lock zooming
    if lock
        bmapax.xpanlock = true
        bmapax.ypanlock = true
        bmapax.xzoomlock = true
        bmapax.yzoomlock = true
    end

    layout[:, 2] = sublayout
    display(scene)
    return scene, layout, parobs
end

function add_obstacle_axis!(scene, sublayout, intervals, bmapax, lock)
    ticklabels = ["$(round(ξ, sigdigits=4))" for ξ in intervals[2:end-1]]
    bmapax.xticks = (Float32[intervals[2:end-1]...], ticklabels)
    for (i, ξ) in enumerate(intervals[2:end-1])
        lines!(bmapax.scene, [Point2f0(ξ, -1), Point2f0(ξ, 1)], linestyle = :dash, color = :black)
    end
    obstacle_ticklabels = String[]
    obstacle_ticks = Float32[]
    for (i, ξ) in enumerate(intervals[1:end-1])
        push!(obstacle_ticks, ξ + (intervals[i+1] - ξ)/2)
        push!(obstacle_ticklabels, string(i))
    end

    obax = sublayout[1,:] = LAxis(scene)
    obax.xticks = (obstacle_ticks, obstacle_ticklabels)
    obax.xaxisposition = :top
    obax.xticklabelalign = (:center, :bottom)
    obax.xlabel = "obstacle index"
    obax.xgridvisible = false
    hideydecorations!(obax)
    hidespines!(obax)
    xlims!(obax, 0, intervals[end])
    if lock
        obax.xpanlock = true
        obax.ypanlock = true
        obax.xzoomlock = true
        obax.yzoomlock = true
    end
    return obax
end

export billiard_bmap_plot

"""
    billiard_bmap_plot(bd::Billiard, ps::Vector{<:AbstractParticle}; kwargs...)
Return a static scene which has particles plotted on both the real billiard as well
the boundary map, each with its own color (keyword `colors`), and the same color is used
for the corresponding scatter points in the boundary map.

All keyword arguments are the same as [`interactive_billiard_bmap`](@ref), besides
the interaction part of course.
"""
function billiard_bmap_plot(bd::Billiard, ps::Vector{<:AbstractParticle};
        ms = 8, plot_particles=true, colors = JULIADYNAMICS_COLORS,
        dt = 0.001, steps = round(Int, 10/dt), kwargs...
    )

    N = length(ps)
    intervals = arcintervals(bd)
    scene, layout, allparobs, balls, vels, vr = interactive_billiard(bd, ps;
        kwargs..., dt = dt, add_controls =false, plot_particles=plot_particles,
        intervals = intervals, res = (1600, 800), colors = colors
    )
    cs = (!(colors isa Vector) || length(colors) ≠ N) ? colors_from_map(colors, 1.0, N) : colors

    sublayout = GridLayout()
    bmapax = sublayout[1,1] = LAxis(scene)
    bmapax.xlabel = "arclength ξ"
    bmapax.ylabel = "normal angle sin(θ)"
    ylims!(bmapax, -1, 1)
    xlims!(bmapax, intervals[1], intervals[end])
    bmapax.xticklabelsize = 28
    bmapax.yticklabelsize = 28
    bmapax.xlabelsize = 36
    bmapax.ylabelsize = 36
    bmapax.ylabelpadding = 20
    if length(intervals) > 2 # at least 2 obstacles
        obax = add_obstacle_axis!(scene, sublayout, intervals, bmapax, false)
        obax.xticklabelsize = 28
        obax.yticklabelsize = 28
        obax.xlabelsize = 36
        obax.ylabelsize = 36
        obax.ylabelpadding = 20
    end

    # create listeners that update boundary map points
    all_bmap_scatters = [Point2f0[] for i in 1:N]
    for i in 1:N
        parobs = allparobs[i]
        vector = all_bmap_scatters[i]
        on(parobs.ξsin) do val
            push!(vector, val)
        end
    end

    # evolve particles until necessary time
    for j in 1:steps-1
        for i in 1:N
            p = ps[i]
            parobs = allparobs[i]
            animstep!(parobs, bd, dt, false, intervals)
            if plot_particles
                balls[][i] = parobs.p.pos
                vels[][i] = vr * parobs.p.vel
            end
        end
    end
    # last step (that actually updates plot)
    for i in 1:N
        p = ps[i]
        parobs = allparobs[i]
        animstep!(parobs, bd, dt, true)
    end
    if plot_particles
        balls[] = balls[]
        vels[] = vels[]
    end

    # Scatterplots:
    for i in 1:N
        c = cs[i]
        vector = all_bmap_scatters[i]
        scatter!(bmapax, vector; color = cs[i],
            marker = MARKER, markersize = ms*AbstractPlotting.px
        )
    end
    ylims!(bmapax, -1, 1)
    xlims!(bmapax, intervals[1], intervals[end])
    layout[:, 2] = sublayout
    return scene, bmapax
end


text:
function train_layer_1_rbf_pde(neurons,vars,obs,n_per_part,inds_all,x,cc,xb,yb,y_source)

    # a_all=zeros(neurons,1000vars)
    a_all=Vector{Vector{Float64}}()
    ooops=0
    for i=1:neurons
        if n_per_part[i+1]-n_per_part[i]!=0
            try
                x1=copy(x[inds_all[n_per_part[i]+1:n_per_part[i+1]],:])
                obs1=size(x1,1)
                phi_x = calc_phi_deriv_deriv(x1,x1,obs1,obs1,cc,vars,1)
                phi_y = calc_phi_deriv_deriv(x1,x1,obs1,obs1,cc,vars,2)
                phi = calc_phi(x1,xb,obs1,size(xb,1),cc,vars)
                atmp=[(phi_x+phi_y);phi]\[y_source[inds_all[n_per_part[i]+1:n_per_part[i+1]]];yb]
                # if sum(isnan.(atmp))==0 && sum(isinf.(atmp))==0 a_all[i,1:obs1].=atmp else ooops+=1 end
                if sum(isnan.(atmp))==0 && sum(isinf.(atmp))==0 push!(a_all,atmp) else ooops+=1 end
            catch ex
                ooops+=1; println(i," ",ex); break
            end
        end
    end
    println("ooops=",ooops)

    layer1=zeros(obs+size(yb,1),neurons)
    for i=1:neurons
        x1=copy(x[inds_all[n_per_part[i]+1:n_per_part[i+1]],:])
        obs1=size(x1,1)
        phi_x = calc_phi_deriv_deriv(x1,x,obs1,obs,cc,vars,1)
        phi_y = calc_phi_deriv_deriv(x1,x,obs1,obs,cc,vars,2)
        phi = calc_phi(x1,xb,obs1,size(xb,1),cc,vars)
        # layer1[:,i].=[(phi_x+phi_y);phi]*a_all[i,1:obs1]
        layer1[:,i].=[(phi_x+phi_y);phi]*a_all[i]
    end

    xxxx=[layer1 ones(obs+size(yb,1))]
    yyyy=[y_source;yb]
    a_layer1=xxxx\yyyy

    return a_all,a_layer1,layer1

end

ext: jl


text:
# https://www.sciencedirect.com/science/article/pii/S0377042705007958
# https://link.springer.com/article/10.1007/BF02592055
# https://epubs.siam.org/doi/abs/10.1137/0917003
# https://epubs.siam.org/doi/10.1137/0911026
# https://pdfs.semanticscholar.org/b321/3084f663260076dcb92f2fa6031b362dc5bc.pdf
# https://www.sciencedirect.com/science/article/abs/pii/0098135483800027
#
abstract type ForcingSequence end

struct FixedForceTerm{T} <: ForcingSequence
    η::T
end
η(fft::FixedForceTerm, info) = fft.η

# Truncated-Newton algorithms for large-scale unconstrained optimization
# https://link.springer.com/article/10.1007/BF02592055
struct DemboSteihaug <: ForcingSequence
end
η(fft::DemboSteihaug, info) = min(1/(info.k+2), info.ρFz)

# See page 19 https://epubs.siam.org/doi/abs/10.1137/0917003
# Choice 1 (2.2) - it appears that the extra F'(x)s evaluation would
# not be worth it. It is at least as small as (2.1). If it's signifi-
# cantly smaller we risk oversolving.
# η₀ ∈ [0, 1)
struct EisenstatWalkerA{T} <: ForcingSequence
    s::T
    ηmax::T
end
EisenstatWalkerA() = EisenstatWalkerA(0.1, 0.9999)
function η(fft::EisenstatWalkerA, info)

    η = (info.ρFz - info.residual_old)/info.ρFx
    β = info.η_old^(1+sqrt(5)/2)
    if β ≥ fft.s
        η = max(η, β)
    end
    min(η, fft.ηmax)
end

# See page 20 https://epubs.siam.org/doi/abs/10.1137/0917003
# Choice 2 (2.6) seems to result in less oversolving.
# γ ∈ [0, 1] - but recommend γ ≥ 0.7
# α ∈ (1, 2]
# η₀ ∈ [0, 1)
# Default to EisenstatA-values. The paper [[[]]] suggests that γ ≥ 0.7 and
# ω ≥ (1+sqrt(5))/2
struct EisenstatWalkerB{T} <: ForcingSequence
    s::T
    α::T
    γ::T
    ηmax::T
end
EisenstatWalkerB() = EisenstatWalkerB(0.1, 1+sqrt(5)/2, 1.0, 100.0)

struct BrownSaad{T}

end
function η(fft::EisenstatWalkerB, info)
    T = typeof(info.ρFz)
    γ, α = fft.γ, fft.α

    ηx = info.η_old

    ρ = info.ρFz/info.ρFx
    η = γ*ρ^α
    β = γ*ηx^α
    if β ≥ fft.s
        η = max(η, β)
    end
    min(η, fft.ηmax)
end

# INB here:
# https://pdfs.semanticscholar.org/b321/3084f663260076dcb92f2fa6031b362dc5bc.pdf
# Originally Walker
# S. C. EiSENSTAT AND H. F. Walker, Globally convergent inexact Newton
# methods, SIAM J. Optim., 4 (1994), pp. 393-422.
"""
    InexactNewton(; force_seq, eta0, maxiter)

Constructs a method type for the Inexact Newton's method with Linesearch.

"""
struct InexactNewton{LinSolve, ForcingType<:ForcingSequence, Tη}
    linsolve::LinSolve
    force_seq::ForcingType
    η₀::Tη
    maxiter::Int
end
#InexactNewton(; force_seq=FixedForceTinexacerm(1e-4), eta0 = 1e-4, maxiter=300)=InexactNewton(force_seq, eta0, maxiter)

InexactNewton(; linsolve, force_seq=DemboSteihaug(), eta0 = 1e-4, maxiter=300)=InexactNewton(linsolve, force_seq, eta0, maxiter)
# map from method to forcing sequence
η(fft::InexactNewton, info) = η(fft.force_seq, info)

function solve(problem::NEqProblem, x, method::InexactNewton, options::NEqOptions)
    if !(mstyle(problem) === InPlace())
        throw(ErrorException("solve() not defined for OutOfPlace() with InexactNewton"))
    end
    t0 = time()

    F = problem.R.F
    JvGen = problem.R.Jv

    Tx = eltype(x)
    xp, Fx = copy(x), copy(x)
    Fx = F(Fx, x)
    JvOp = JvGen(x)

    ρ2F0 = norm(Fx, 2)
    ρFz = ρ2F0
    ρF0 = norm(Fx, Inf)
    ρs = norm(x, 2)
    z = copy(x)
    stoptol = Tx(options.f_reltol)*ρFz + Tx(options.f_abstol)

    force_info = (k = 1, ρFz=ρFz, ρFx=nothing, η_old=nothing)

    iter = 0

    while iter < options.maxiter
        iter += 1
        x .= z
        # Refactor this
        if iter == 1 && !isa(method.force_seq, FixedForceTerm)
            ηₖ = method.η₀
        else
            ηₖ = η(method, force_info)
        end
        JvOp = JvGen(x)

        xp .= 0
        xp, res = method.linsolve(xp, JvOp, Fx, ηₖ)

        ρs = norm(xp)
        t = 1e-4
        # use line searc functions here
        btk_conv = false
        ρFx = force_info.ρFz

        # Backtracking
        #
        it = 0
        while !btk_conv
            it += 1
            z = retract(problem, z, x, xp, -1)

            Fx = problem.R.F(Fx, z)
            btk_conv = norm(Fx, 2) ≤ (1-t*(1-ηₖ))*ρFx || it > 20
        end
        if norm(Fx, 2) < stoptol
            break
        end
        η_old = ηₖ
        ρFz = norm(Fx, 2)
        force_info = (k = iter, ρFz=ρFz, ρFx=ρFx, η_old=η_old, residual_old=res)
    end
    return ConvergenceInfo(method, (solution=x, best_residual=Fx, ρF0=ρF0, ρ2F0=ρ2F0, ρs=ρs, iter=iter, time=time()-t0), options)
end


text:
using Test, Flux
using DiffEqDevTools
using DiffEqBase, NeuralNetDiffEq

# Homogeneous on scalar
f = (du, u, p, t) -> -du+6*u
tspan = (0.0f0, 1.0f0)
u0 = 1.0f0
du0 = 0.0f0
dt = 1/5f0
f_analytic = t -> [1/5 * exp(-3*t) *(3*exp(5*t)+2.0)]
prob = SecondOrderODEProblem(ODEFunction{false}(f, analytic=f_analytic), u0, du0, tspan)
chain = Chain(Dense(1,5,σ),Dense(5,1))
opt = ADAM(1, (0.9, 0.95))
sol = solve(prob, NeuralNetDiffEq.NNODE2(chain,opt), dt=dt, verbose = true, maxiters = 200)
err = sqrt(sum((sol.u - [f_analytic(t)[1] for t in tspan[1]:dt:tspan[2]]).^2)[1])
@test err< 0.5

#Example 2: Test on vector
f = (du, u, p, t) -> -2*du-5*u+5*t^2+12
tspan = (0.0f0, 1.0f0)
u0 = [1.0f0]
du0 = [0.0f0]
dt = 1/5f0
f_analytic = t -> [1/50 * exp(-t) * (2*exp(t) * (25*t^2 - 20*t + 58) - 13*sin(2*t) - 66*cos(2*t))]
prob = SecondOrderODEProblem(ODEFunction{false}(f, analytic=f_analytic), u0, du0, tspan)
chain = Chain(Dense(1,5,σ),Dense(5,1))
opt = ADAM(0.1, (0.9, 0.95))
sol = solve(prob, NeuralNetDiffEq.NNODE2(chain,opt), dt=dt, verbose = true, maxiters = 300)
err = sqrt(sum((sol.u - [f_analytic(t)[1] for t in tspan[1]:dt:tspan[2]]).^2)[1])
@test err < 0.5

#Example 3: Test on system of 2nd order ODEs
function func(ddu, du, u, p, t)
	ddu[1] = (-p[3]*u[1]-p[4]*(u[1]-u[2]))/p[1]
	ddu[2] = (-p[4]*(u[2]-u[1]))/p[2]
end
function func(du, u, p, t)
	[(-p[3]*u[1]-p[4]*(u[1]-u[2]))/p[1], (-p[4]*(u[2]-u[1]))/p[2]]
end
p = [1.0f0, 1.0f0, 6.0f0, 4.0f0]
u0 = [1.0f0, 0.0f0]
du0 = [2.0f0, 0.0f0]
prob = SecondOrderODEProblem{false}(func, [u0], [du0], (0.0f0, 3.0f0), p)
dt = 1/20f0
opt = ADAM(1e-01)
chain = Chain(Dense(1,5,σ),Dense(5,2))
sol = solve(prob, NeuralNetDiffEq.NNODE2(chain, opt), dt=dt, verbose=true, abstol=1e-10, maxiters=50)


text:
using InteractionNetworkModels, Plots, Distributions, JLD2
using StatsPlots, Plots.Measures, StatsBase

cd("Z:/simulations/")
readdir()
chains, times = load(
    "iex_sim_50_match_dist_multi.jld2",
    "chains", "times"
)

plot([x.suff_stat_trace for x in chains])
plot([length.(x.S_sample) for x in chains])

plot(chains, E)

posterior = chains[1].posterior
data = posterior.data
d = posterior.dist

Ê, ind = sample_frechet_mean(data, d)
sum(x->d(x,Ê), data)


text:
using Test
using ExtensibleEffects
using DataTypesBasic
using TypeClasses
using Suppressor
splitln(str) = split(strip(str), "\n")
import ExtensibleEffects: autorun

result = @syntax_eff begin
    a = NoEffect(1)
    b = NoEffect(a + 5.0)
    @pure a, b
end
@test result == (1, 6.0)

max_stars_repo_name: schlichtanders/ExtensibleEffects

max_issues_repo_name: schlichtanders/ExtensibleEffects

max_forks_repo_name: schlichtanders/ExtensibleEffects


text:
function waterfall(;x::AbstractVector = Int[], y::AbstractVector = Int[])
    if length(x) != length(y)
        error("X and Y need to be the same length")
    end

    #Calculate waterfall
    _start = zeros(Int, length(x))
    _end = cumsum(y) #End values for bar

    #Do lag calculation for bar starting value
    for i in 2:length(_start)
        _start[i] = _end[i-1]
    end

    #Calculate group based on positive or negative
    _group = [Int(val < 0) for val in y]

    #Make a copy so x isn't mutated
    _x = deepcopy(x)

    #Add Summary to Arrays
    push!(_start, 0)
    push!(_end, _end[end])
    push!(_group, _group[end] + 1)
    push!(_x, "Total")

    v = barplot(x = _x, y = _end, y2 = _start, group = _group)

    #Clean-up presentation
    v.width = 110 * length(x) #Make the graph variable width based on number of levels
    v.marks[1].properties.enter.width = VegaValueRef(scale = "x", band = true, mult= 0.5)
    v.marks[1].properties.enter.x  = VegaValueRef(scale = "x", field = "x", offset = 31)
    v.scales[1].domain.sort = false

    #Hide unneeded labels
    xlab!(v)
    ylab!(v)
    legend!(v, show = false)

    colorscheme!(v; palette = ["#4daf4a", "#e41a1c", "#377eb8"])

    return v
end

max_stars_repo_name: jpthor/vega


text:
struct ADAM
    η::Float64
    β::Tuple{Float64,Float64}
end
struct ADAMState{A}
    mv::A
    vt::A
    βp::Base.RefValue{Tuple{Float64,Float64}}
end
ADAMState(

function update!(θ::AbstractArray{T}, s::ADAMState, Δ::AbstractArray{T}, a::ADAM) where {T}
    @unpack η, (β₁, β₂) = a
    @unpack mv, vt, βp = s
    βp₁, βp₂ = βp
    omβ₁ = T(1 - β₁); βp₁⁻¹ = T(inv(1 - βp₁))
    omβ₂ = T(1 - β₂); βp₂⁻ʰ = T(inv(sqrt(1 - βp₂)))
    @avx for i ∈ eachindex(Δ)
        Δₒ = Δ[i]
        mtₒ = mt[i]
        mtₙ = β₁ * mtₒ + omβ₁ * Δₒ
        vtₒ = vt[i]
        vtₙ = β₂ * vtₒ + omβ₂ * Δₒ * Δₒ
        Δₙ = mtₙ * βp₁⁻¹ / (sqrt(vtₙ)*βp₂⁻ʰ + ϵ)
        mt[i] = mtₙ
        vt[i] = vtₙ
        θ[i] -= Δₙ * η
    end
    s.βp[] = (βp₁*β₁, βp₂*β₂); nothing
end
update!(::Nothing, ::Nothing, ::Nothing, ::ADAM) = nothing


text:
"""
Type: LAStopping

Methods: `start!`, `stop!`, `update_and_start!`, `update_and_stop!`, `fill_in!`, `reinit!`, `status`,
`linear_system_check`, `normal_equation_check`

Specialization of GenericStopping. Stopping structure for linear algebra
solving either

``Ax = b``

or

```math
min\\_{x} \\tfrac{1}{2}\\|Ax - b\\|^2
```

Attributes:
- `pb`         : a problem using, for instance, either `LLSModel` (designed for linear least square problem, see https://github.com/JuliaSmoothOptimizers/LLSModels.jl ) or `LinearSystem`.
- `current_state`      : The information relative to the problem, see `GenericState`.
- (opt) `meta` : Metadata relative to stopping criteria, see `StoppingMeta`.
- (opt) `main_stp` : Stopping of the main loop in case we consider a Stopping
                          of a subproblem.
                          If not a subproblem, then `VoidStopping`.
- (opt) `listofstates` : ListofStates designed to store the history of States.
- (opt) `stopping_user_struct` : Contains a structure designed by the user.

Constructors:
- `LAStopping(pb, meta::AbstractStoppingMeta, stop_remote::AbstractStopRemoteControl, state::AbstractState; main_stp::AbstractStopping=VoidStopping(), list::AbstractListofStates = VoidListofStates(), user_struct::AbstractDict = Dict(), zero_start::Bool = false)`
     The default constructor.
- `LAStopping(pb, meta::AbstractStoppingMeta, state::AbstractState; main_stp::AbstractStopping=VoidStopping(), list::AbstractListofStates = VoidListofStates(), user_struct::AbstractDict = Dict(), zero_start::Bool = false, kwargs...)`
     The one passing the `kwargs` to the `stop_remote`.
- `LAStopping(pb, state::AbstractState; stop_remote::AbstractStopRemoteControl = StopRemoteControl(), main_stp::AbstractStopping=VoidStopping(), list::AbstractListofStates = VoidListofStates(), user_struct::AbstractDict = Dict(), zero_start::Bool = false, kwargs...)`
     The one passing the `kwargs` to the `meta`.
- `LAStopping(:: Union{AbstractLinearOperator, AbstractMatrix}, :: AbstractVector; sparse::Bool = true, n_listofstates::Int = 0, kwargs...)`
     The one setting up a default problem (`sparse ? LLSModel(A, b) : LinearSystem(A, b)`), a default `GenericState` using x, and initializing the list of states if `n_listofstates>0`.
- `LAStopping(:: Union{AbstractLinearOperator, AbstractMatrix}, :: AbstractVector, :: AbstractState; sparse::Bool = true, kwargs...)`
     The one setting up a default problem (`sparse ? LLSModel(A, b) : LinearSystem(A, b)`).

Notes:
- No specific State targeted
- State don't necessarily keep track of evals
- Evals are checked only for `pb.A` being a LinearOperator
- `zero_start` is true if 0 is the initial guess (not check automatically)
- `LLSModel` counter follow `NLSCounters` (see `init_max_counters_NLS` in `NLPStoppingmod.jl`)
- By default, `meta.max_cntrs` is initialized with an NLSCounters

See also `GenericStopping`, `NLPStopping`, `linear_system_check`, `normal_equation_check`
 """
mutable struct LAStopping{Pb, M, SRC, T, MStp, LoS} <: AbstractStopping{Pb, M, SRC, T, MStp, LoS}

  # problem
  pb::Pb
  # Common parameters
  meta::M
  stop_remote::SRC
  # current state of the problem
  current_state::T
  # Stopping of the main problem, or nothing
  main_stp::MStp
  # History of states
  listofstates::LoS
  # User-specific structure
  stopping_user_struct::AbstractDict

  #zero is initial point
  zero_start::Bool
end

function LAStopping(
  pb::Pb,
  meta::M,
  stop_remote::SRC,
  current_state::T;
  main_stp::AbstractStopping = VoidStopping(),
  list::AbstractListofStates = VoidListofStates(),
  user_struct::AbstractDict = Dict(),
  zero_start::Bool = false,
) where {Pb <: Any, M <: AbstractStoppingMeta, SRC <: AbstractStopRemoteControl, T <: AbstractState}
  return LAStopping(pb, meta, stop_remote, current_state, main_stp, list, user_struct, zero_start)
end

function LAStopping(
  pb::Pb,
  meta::M,
  current_state::T;
  main_stp::AbstractStopping = VoidStopping(),
  list::AbstractListofStates = VoidListofStates(),
  user_struct::AbstractDict = Dict(),
  zero_start::Bool = false,
  kwargs...,
) where {Pb <: Any, M <: AbstractStoppingMeta, T <: AbstractState}
  stop_remote = StopRemoteControl(; kwargs...) #main_stp == VoidStopping() ? StopRemoteControl() : cheap_stop_remote_control()

  return LAStopping(pb, meta, stop_remote, current_state, main_stp, list, user_struct, zero_start)
end

function LAStopping(
  pb::Pb,
  current_state::T;
  stop_remote::AbstractStopRemoteControl = StopRemoteControl(), #main_stp == VoidStopping() ? StopRemoteControl() : cheap_stop_remote_control()
  main_stp::AbstractStopping = VoidStopping(),
  list::AbstractListofStates = VoidListofStates(),
  user_struct::AbstractDict = Dict(),
  zero_start::Bool = false,
  kwargs...,
) where {Pb <: Any, T <: AbstractState}
  if :max_cntrs in keys(kwargs)
    mcntrs = kwargs[:max_cntrs]
  elseif Pb <: LLSModel
    mcntrs = init_max_counters_NLS()
  else
    mcntrs = init_max_counters_linear_operators()
  end

  if :optimality_check in keys(kwargs)
    oc = kwargs[:optimality_check]
  else
    oc = linear_system_check
  end

  meta = StoppingMeta(; max_cntrs = mcntrs, optimality_check = oc, kwargs...)

  return LAStopping(pb, meta, stop_remote, current_state, main_stp, list, user_struct, zero_start)
end

function LAStopping(
  A::TA,
  b::Tb;
  x::Tb = zeros(eltype(Tb), size(A, 2)),
  sparse::Bool = true,
  n_listofstates::Int = 0,
  kwargs...,
) where {TA <: Any, Tb <: AbstractVector}
  pb = sparse ? LLSModel(A, b) : LinearSystem(A, b)
  state = GenericState(x)

  mcntrs = sparse ? init_max_counters_NLS() : init_max_counters_linear_operators()

  if n_listofstates > 0 && :list ∉ keys(kwargs)
    list = ListofStates(n_listofstates, Val{typeof(state)}())
    return LAStopping(pb, state, max_cntrs = mcntrs, list = list; kwargs...)
  end

  return LAStopping(pb, state, max_cntrs = mcntrs; kwargs...)
end

function LAStopping(
  A::TA,
  b::Tb,
  state::S;
  sparse::Bool = true,
  kwargs...,
) where {TA <: Any, Tb <: AbstractVector, S <: AbstractState}
  pb = sparse ? LLSModel(A, b) : LinearSystem(A, b)

  mcntrs = sparse ? init_max_counters_NLS() : init_max_counters_linear_operators()

  return LAStopping(pb, state, max_cntrs = mcntrs; kwargs...)
end

"""
Type: LACounters
"""
mutable struct LACounters{T <: Int}
  nprod::T
  ntprod::T
  nctprod::T
  sum::T

  function LACounters(nprod::T, ntprod::T, nctprod::T, sum::T) where {T <: Int}
    return new{T}(nprod, ntprod, nctprod, sum)
  end
end

function LACounters(; nprod::Int64 = 0, ntprod::Int64 = 0, nctprod::Int64 = 0, sum::Int64 = 0)
  return LACounters(nprod, ntprod, nctprod, sum)
end

"""
init\\_max\\_counters\\_linear\\_operators: counters for LinearOperator

`init_max_counters_linear_operators(; allevals :: T = 20000, nprod = allevals, ntprod = allevals, nctprod = allevals, sum = 11 * allevals)`
"""
function init_max_counters_linear_operators(;
  allevals::T = 20000,
  nprod::T = allevals,
  ntprod::T = allevals,
  nctprod::T = allevals,
  sum::T = allevals * 11,
) where {T <: Int}
  cntrs =
    Dict{Symbol, T}([(:nprod, nprod), (:ntprod, ntprod), (:nctprod, nctprod), (:neval_sum, sum)])

  return cntrs
end

"""
LinearSystem: Minimal structure to store linear algebra problems

`LinearSystem(:: Union{AbstractLinearOperator, AbstractMatrix}, :: AbstractVector)`

Note:
Another option is to convert the `LinearSystem` as an `LLSModel`.
"""
mutable struct LinearSystem{
  TA <: Union{AbstractLinearOperator, AbstractMatrix},
  Tb <: AbstractVector,
}
  A::TA
  b::Tb

  counters::LACounters

  function LinearSystem(
    A::TA,
    b::Tb;
    counters::LACounters = LACounters(),
    kwargs...,
  ) where {TA <: Union{AbstractLinearOperator, AbstractMatrix}, Tb <: AbstractVector}
    return new{TA, Tb}(A, b, counters)
  end
end

function LAStopping(
  A::TA,
  b::Tb;
  x::Tb = zeros(eltype(Tb), size(A, 2)),
  kwargs...,
) where {TA <: AbstractLinearOperator, Tb <: AbstractVector}
  return LAStopping(A, b, GenericState(x), kwargs...)
end

function LAStopping(
  A::TA,
  b::Tb,
  state::AbstractState;
  kwargs...,
) where {TA <: AbstractLinearOperator, Tb <: AbstractVector}
  return LAStopping(
    LinearSystem(A, b),
    state,
    max_cntrs = init_max_counters_linear_operators(),
    kwargs...,
  )
end

"""
 \\_resources\\_check!: check if the optimization algorithm has
exhausted the resources. This is the Linear Algebra specialized version.

 Note:
 - function does _not_ keep track of the evals in the state
 - check `:nprod`, `:ntprod`, and `:nctprod` in the `LinearOperator` entries
 """
function _resources_check!(stp::LAStopping, x::T) where {T}

  #GenericState has no field evals.
  #_smart_update!(stp.current_state, evals = cntrs)

  # check all the entries in the counter
  # global user limit diagnostic
  stp.meta.resources = _counters_loop!(stp.pb.counters, stp.meta.max_cntrs)

  return stp.meta.resources
end

function _counters_loop!(cntrs::LACounters{T}, max_cntrs::Dict{Symbol, T}) where {T}
  sum, max_f = 0, false

  for f in (:nprod, :ntprod, :nctprod)
    ff = getfield(cntrs, f)
    max_f = max_f || (ff > max_cntrs[f])
    sum += ff
  end

  return max_f || (sum > max_cntrs[:neval_sum])
end

function _counters_loop!(cntrs::NLSCounters, max_cntrs::Dict{Symbol, T}) where {T}
  sum, max_f = 0, false

  for f in intersect(fieldnames(NLSCounters), keys(max_cntrs))
    max_f = f != :counters ? (max_f || (getfield(cntrs, f) > max_cntrs[f])) : max_f
  end
  for f in intersect(fieldnames(Counters), keys(max_cntrs))
    max_f = max_f || (getfield(cntrs.counters, f) > max_cntrs[f])
  end

  return max_f || (sum > max_cntrs[:neval_sum])
end

"""
linear\\_system\\_check: return ||Ax-b||_p

`linear_system_check(:: Union{LinearSystem, LLSModel}, :: AbstractState; pnorm :: Real = Inf, kwargs...)`

Note:
- Returns the p-norm of state.res
- state.res is filled in if nothing.
"""
function linear_system_check(pb::LinearSystem, state::AbstractState; pnorm::Real = Inf, kwargs...)
  pb.counters.nprod += 1
  if state.res == _init_field(typeof(state.res))
    update!(state, res = pb.A * state.x - pb.b)
  end

  return norm(state.res, pnorm)
end

function linear_system_check(pb::LLSModel, state::AbstractState; pnorm::Real = Inf, kwargs...)
  if state.res == _init_field(typeof(state.res))
    Axmb = if xtype(state) <: SparseVector
      sparse(residual(pb, state.x))
    else
      residual(pb, state.x)
    end
    update!(state, res = Axmb)
  end

  return norm(state.res, pnorm)
end

"""
normal\\_equation\\_check: return ||A'Ax-A'b||_p

`normal_equation_check(:: Union{LinearSystem, LLSModel}, :: AbstractState; pnorm :: Real = Inf, kwargs...)`

Note: pb must have A and b entries
"""
function normal_equation_check(pb::LinearSystem, state::AbstractState; pnorm::Real = Inf, kwargs...)
  pb.counters.nprod += 1
  pb.counters.ntprod += 1
  return norm(pb.A' * (pb.A * state.x) - pb.A' * pb.b, pnorm)
end

function normal_equation_check(pb::LLSModel, state::AbstractState; pnorm::Real = Inf, kwargs...)
  nres = jtprod_residual(pb, state.x, residual(pb, state.x))
  return norm(nres, pnorm)
end


text:
"""
    ReluVal(max_iter::Int64, tree_search::Symbol)

ReluVal combines symbolic reachability analysis with iterative interval refinement to minimize over-approximation of the reachable set.

# Problem requirement
1. Network: any depth, ReLU activation
2. Input: hyperrectangle
3. Output: hpolytope

# Return
`CounterExampleResult` or `ReachabilityResult`

# Method
Symbolic reachability analysis and iterative interval refinement (search).
- `max_iter` default `10`.
- `tree_search` default `:DFS` - depth first search.

# Property
Sound but not complete.

# Reference
[S. Wang, K. Pei, J. Whitehouse, J. Yang, and S. Jana,
"Formal Security Analysis of Neural Networks Using Symbolic Intervals,"
*CoRR*, vol. abs/1804.10829, 2018. arXiv: 1804.10829.](https://arxiv.org/abs/1804.10829)

[https://github.com/tcwangshiqi-columbia/ReluVal](https://github.com/tcwangshiqi-columbia/ReluVal)
"""
@with_kw struct ReluVal
    max_iter::Int64     = 10
    tree_search::Symbol = :DFS # only :DFS/:BFS allowed? If so, we should assert this.
end

struct SymbolicInterval
    Low::Matrix{Float64}
    Up::Matrix{Float64}
    interval::Hyperrectangle
end

SymbolicInterval(x::Matrix{Float64}, y::Matrix{Float64}) = SymbolicInterval(x, y, Hyperrectangle([0],[0]))

# Data to be passed during forward_layer
struct SymbolicIntervalMask
    sym::SymbolicInterval
    LΛ::Vector{Vector{Int64}}
    UΛ::Vector{Vector{Int64}}
end

function solve(solver::ReluVal, problem::Problem)
    reach = forward_network(solver, problem.network, problem.input)
    result = check_inclusion(reach.sym, problem.output, problem.network)
    result.status == :Unknown || return result
    reach_list = SymbolicIntervalMask[reach]
    for i in 2:solver.max_iter
        length(reach_list) > 0 || return BasicResult(:SAT)
        reach = pick_out!(reach_list, solver.tree_search)
        LG, UG = get_gradient(problem.network, reach.LΛ, reach.UΛ)
        feature = get_smear_index(problem.network, reach.sym.interval, LG, UG)
        intervals = split_interval(reach.sym.interval, feature)
        for interval in intervals
            reach = forward_network(solver, problem.network, interval)
            result = check_inclusion(reach.sym, problem.output, problem.network)
            result.status == :UNSAT && return result
            result.status == :SAT || (push!(reach_list, reach))
        end
    end
    return BasicResult(:Unknown)
end

function pick_out!(reach_list, tree_search)
    if tree_search == :BFS
        reach = reach_list[1]
        deleteat!(reach_list, 1)
    else
        n = length(reach_list)
        reach = reach_list[n]
        deleteat!(reach_list, n)
    end
    return reach
end

function symbol_to_concrete(reach::SymbolicInterval)
    n_output = size(reach.Low, 1)
    upper = zeros(n_output)
    lower = zeros(n_output)
    for i in 1:n_output
        lower[i] = lower_bound(reach.Low[i, :], reach.interval)
        upper[i] = upper_bound(reach.Low[i, :], reach.interval)
    end
    return Hyperrectangle(low = lower, high = upper)
end

function check_inclusion(reach::SymbolicInterval, output::AbstractPolytope, nnet::Network)
    reachable = symbol_to_concrete(reach)
    issubset(reachable, output) && return BasicResult(:SAT)
    is_intersection_empty(reachable, output) && return BasicResult(:UNSAT)

    # Sample the middle point
    middle_point = (high(reach.interval) + low(reach.interval))./2
    y = compute_output(nnet, middle_point)
    ∈(y, output) || return CounterExampleResult(:UNSAT, middle_point)

    return BasicResult(:Unknown)
end

function forward_layer(solver::ReluVal, layer::Layer, input::Union{SymbolicIntervalMask, Hyperrectangle})
    return forward_act(forward_linear(input, layer))
end

# Symbolic forward_linear for the first layer
function forward_linear(input::Hyperrectangle, layer::Layer)
    (W, b) = (layer.weights, layer.bias)
    sym = SymbolicInterval(hcat(W, b), hcat(W, b), input)
    LΛ = Vector{Vector{Int64}}(undef, 0)
    UΛ = Vector{Vector{Int64}}(undef, 0)
    return SymbolicIntervalMask(sym, LΛ, UΛ)
end

# Symbolic forward_linear
function forward_linear(input::SymbolicIntervalMask, layer::Layer)
    (W, b) = (layer.weights, layer.bias)
    output_Up = max.(W, 0) * input.sym.Up + min.(W, 0) * input.sym.Low
    output_Low = max.(W, 0) * input.sym.Low + min.(W, 0) * input.sym.Up
    output_Up[:, end] += b
    output_Low[:, end] += b
    sym = SymbolicInterval(output_Low, output_Up, input.sym.interval)
    return SymbolicIntervalMask(sym, input.LΛ, input.UΛ)
end

# Symbolic forward_act
function forward_act(input::SymbolicIntervalMask)
    n_node, n_input = size(input.sym.Up)
    output_Low, output_Up = input.sym.Low[:, :], input.sym.Up[:, :]
    mask_lower, mask_upper = zeros(Int, n_node), ones(Int, n_node)
    for i in 1:n_node
        if upper_bound(input.sym.Up[i, :], input.sym.interval) <= 0.0
            # Update to zero
            mask_lower[i], mask_upper[i] = 0, 0
            output_Up[i, :] = zeros(n_input)
            output_Low[i, :] = zeros(n_input)
        elseif lower_bound(input.sym.Low[i, :], input.sym.interval) >= 0
            # Keep dependency
            mask_lower[i], mask_upper[i] = 1, 1
        else
            # Concretization
            mask_lower[i], mask_upper[i] = 0, 1
            output_Low[i, :] = zeros(n_input)
            if lower_bound(input.sym.Up[i, :], input.sym.interval) < 0
                output_Up[i, :] = zeros(n_input)
                output_Up[i, end] = upper_bound(input.sym.Up[i, :], input.sym.interval)
            end
        end
    end
    sym = SymbolicInterval(output_Low, output_Up, input.sym.interval)
    LΛ = push!(input.LΛ, mask_lower)
    UΛ = push!(input.UΛ, mask_upper)
    return SymbolicIntervalMask(sym, LΛ, UΛ)
end

# Return the splited intervals
function get_smear_index(nnet::Network, input::Hyperrectangle, LG::Matrix, UG::Matrix)
    largest_smear = - Inf
    feature = 0
    r = input.radius
    for i in 1:dim(input)
        smear = sum(max.(abs.(UG[:, i]), abs.(LG[:, i]))) * r[i]
        if smear > largest_smear
            largest_smear = smear
            feature = i
        end
    end
    return feature
end

# Get upper bound in concretization
function upper_bound(map::Vector{Float64}, input::Hyperrectangle)
    bound = map[dim(input)+1]
    input_upper = high(input)
    input_lower = low(input)
    for i in 1:dim(input)
        bound += ifelse( map[i]>0, map[i]*input_upper[i], map[i]*input_lower[i])
    end
    return bound
end

# Get lower bound in concretization
function lower_bound(map::Vector{Float64}, input::Hyperrectangle)
    bound = map[dim(input)+1]
    input_upper = high(input)
    input_lower = low(input)
    for i in 1:dim(input)
        bound += ifelse(map[i]>0, map[i]*input_lower[i], map[i]*input_upper[i])
    end
    return bound
end

# Concrete forward_linear
# function forward_linear_concrete(input::Hyperrectangle, W::Matrix{Float64}, b::Vector{Float64})
#     n_output = size(W, 1)
#     output_upper = zeros(n_output)
#     output_lower = zeros(n_output)
#     for j in 1:n_output
#         output_upper[j] = upper_bound(W[j, :], input) + b[j]
#         output_lower[j] = lower_bound(W[j, :], input) + b[j]
#     end
#     output = Hyperrectangle(low = output_lower, high = output_upper)
#     return output
# end

# Concrete forward_act
# function forward_act(input::Hyperrectangle)
#     input_upper = high(input)
#     input_lower = low(input)
#     output_upper = zeros(dim(input))
#     output_lower = zeros(dim(input))
#     mask_upper = ones(Int, dim(input))
#     mask_lower = zeros(Int, dim(input))
#     for i in 1:dim(input)
#         if input_upper[i] <= 0
#             mask_upper[i] = mask_lower[i] = 0
#             output_upper[i] = output_lower[i] = 0.0
#         elseif input_lower[i] >= 0
#             mask_upper[i] = mask_lower[i] = 1
#         else
#             output_lower[i] = 0.0
#         end
#     end
#     return (output_lower, output_upper, mask_lower, mask_upper)
# end

# This function is similar to forward_linear
# function backward_linear(Low::Matrix{Float64}, Up::Matrix{Float64}, W::Matrix{Float64})
#     n_output, n_input = size(W)
#     n_symbol = size(Low, 2) - 1

#     output_Low = zeros(n_output, n_symbol + 1)
#     output_Up = zeros(n_output, n_symbol + 1)
#     for k in 1:n_symbol + 1
#         for j in 1:n_output
#             for i in 1:n_input
#                 output_Up[j, k] += ifelse(W[j, i]>0, W[j, i] * Up[i, k], W[j, i] * Low[i, k])
#                 output_Low[j, k] += ifelse(W[j, i]>0, W[j, i] * Low[i, k], W[j, i] * Up[i, k])
#             end
#         end
#     end
#     return (output_Low, output_Up)
# end


text:
"""
    NormalInverseChisq(μ, σ2, κ, ν)

A Normal-χ^-2 distribution is a conjugate prior for a Normal distribution with
unknown mean and variance.  It has parameters:

* μ: expected mean
* σ2 > 0: expected variance
* κ ≥ 0: mean confidence
* ν ≥ 0: variance confidence

The parameters have a natural interpretation when used as a prior for a Normal
distribution with unknown mean and variance: μ and σ2 are the expected mean and
variance, while κ and ν are the respective degrees of confidence (expressed in
"pseudocounts").  When interpretable parameters are important, this makes it a
slightly more convenient parametrization of the conjugate prior.

Equivalent to a `NormalInverseGamma` distribution with parameters:

* m0 = μ
* v0 = 1/κ
* shape = ν/2
* scale = νσ2/2

Based on Murphy "Conjugate Bayesian analysis of the Gaussian distribution".
"""
struct NormalInverseChisq{T<:Real} <: ContinuousUnivariateDistribution
    μ::T
    σ2::T
    κ::T
    ν::T

    function NormalInverseChisq{T}(μ::T, σ2::T, κ::T, ν::T) where T<:Real
        if ν < 0 || κ < 0 || σ2 ≤ 0
            throw(ArgumentError("Variance and confidence (κ and ν) must all be positive"))
        end
        new{T}(μ, σ2, κ, ν)
    end
end

NormalInverseChisq() = NormalInverseChisq{Float64}(0.0, 1.0, 0.0, 0.0)

function NormalInverseChisq(μ::Real, σ2::Real, κ::Real, ν::Real)
    T = promote_type(typeof(μ), typeof(σ2), typeof(κ), typeof(ν))
    NormalInverseChisq{T}(T(μ), T(σ2), T(κ), T(ν))
end

Base.convert(::Type{NormalInverseGamma}, d::NormalInverseChisq) =
    NormalInverseGamma(d.μ, 1/d.κ, d.ν/2, d.ν*d.σ2/2)

Base.convert(::Type{NormalInverseChisq}, d::NormalInverseGamma) =
    NormalInverseChisq(d.mu, d.scale/d.shape, 1/d.v0, d.shape*2)

insupport(::Type{NormalInverseChisq}, μ::T, σ2::T) where T<:Real =
    isfinite(μ) && zero(σ2) <= σ2 < Inf

params(d::NormalInverseChisq) = d.μ, d.σ2, d.κ, d.ν

function pdf(d::NormalInverseChisq, μ::T, σ2::T) where T<:Real
    Zinv = sqrt(d.κ / 2pi) / gamma(d.ν*0.5) * (d.ν * d.σ2 / 2)^(d.ν*0.5)
    Zinv * σ2^(-(d.ν+3)*0.5) * exp( (d.ν*d.σ2 + d.κ*(d.μ - μ)^2) / (-2 * σ2))
end

function logpdf(d::NormalInverseChisq, μ::T, σ2::T) where T<:Real
    logZinv = (log(d.κ) - log(2pi))*0.5 - lgamma(d.ν*0.5) + (log(d.ν) + log(d.σ2) - log(2)) * (d.ν/2)
    logZinv + log(σ2)*(-(d.ν+3)*0.5) + (d.ν*d.σ2 + d.κ*(d.μ - μ)^2) / (-2 * σ2)
end

function mean(d::NormalInverseChisq)
    μ = d.μ
    σ2 = d.ν/(d.ν-2)*d.σ2
    return μ, σ2
end

function mode(d::NormalInverseChisq)
    μ = d.μ
    σ2 = d.ν*d.σ2/(d.ν + 3)
    return μ, σ2
end

rand(d::NormalInverseChisq) = rand(convert(NormalInverseGamma, d))


text:
module MINPACK

using Distances
using Printf, LinearAlgebra

export fsolve

const _dl_ext = @static Sys.isapple() ? "dylib" : (Sys.iswindows() ? "dll" : "so")

const cminpack = joinpath(
    dirname(dirname(@__FILE__)), "deps", "libcminpack.$(_dl_ext)"
)

struct IterationState
    iteration::Int
    fnorm::Float64
    xnorm::Float64
    step_time::Float64
end

function Base.show(io::IO, t::IterationState)
    @printf io "%6d   %14e   %14e   %14f\n" t.iteration t.fnorm t.xnorm t.step_time
end

mutable struct AlgoTrace
    f_calls::Int
    g_calls::Int
    show_trace::Bool
    tracing::Bool
    maxit::Int
    x_old::Vector{Float64}
    trace::Vector{IterationState}
    prev_iflag::Int  # allows us to track when we stop computing deriv
    start_time::Float64
    last_feval_time::Float64
    tot_time::Float64

    function AlgoTrace(x_init::Vector{Float64}, verbose::Bool=false, tracing::Bool=false,
                       maxit::Int=typemax(Int))
        if verbose
            tracing = true
        end
        x_old = tracing ? copy(x_init) : Array{Float64}(undef, 0)
        states = Array{IterationState}(undef,0)
        new(0, 0, verbose, tracing, maxit, x_old, states, 1, time(), time(), NaN)
    end
end

Base.unsafe_convert(::Type{Ptr{Cvoid}}, o::AlgoTrace) = o

function Base.show(io::IO, trace::AlgoTrace)
    @printf io "Iter     f(x) inf-norm    Step 2-norm      Step time\n"
    @printf io "------   --------------   --------------   --------------\n"
    for state in trace.trace
        show(io, state)
    end
end

function Base.push!(trace::AlgoTrace, x::Vector{Float64}, fvec::Vector{Float64},
                    iflag::Cint)
    if trace.tracing
        if iflag == 2  # computing derivative
            if trace.prev_iflag == 1
                # only increment if just starting to compute deriv
                trace.g_calls += 1
            end
        elseif iflag == 1  # computing function
            trace.f_calls += 1
            x_step = sqeuclidean(trace.x_old, x)
            f_norm = maximum(abs, fvec)
            now = time()
            elapsed = now - trace.last_feval_time
            trace.last_feval_time = now
            ss = IterationState(trace.f_calls, f_norm, x_step, elapsed)
            trace.show_trace && show(ss)
            push!(trace.trace, ss)
            copyto!(trace.x_old, x)
        end
        trace.prev_iflag = iflag
    end
end

"""
    SolverResults

Holds results of attempting to solve system of equations.

Fields should be accessed directly and are:

- `algo::String`: the name of the lagorithm used
- `initial_x::Vector{Float64}`: the starting point
- `x::Vector{Float64}`: The final value of `x`. When converged evaluating the objective at `x` should give zeros
- `f::Vector{Float64}`: The function value at `x`
- `return_code::Int`: The return code from MINPACK
- `converged::Bool`: Whether or not the algorithm converged
- `msg::String`: The message from MINPACK describing outcome
- `trace::AlgoTrace`: If tracing was enabled, a detailed trace of all iterations
"""
struct SolverResults
    algo::String
    initial_x::Vector{Float64}
    x::Vector{Float64}
    f::Vector{Float64}
    return_code::Int
    converged::Bool
    msg::String
    trace::AlgoTrace
end

# NOTE: this method was adapted from NLsolve.jl
function Base.show(io::IO, s::SolverResults)
    @printf io "Results of Nonlinear Solver Algorithm\n"
    @printf io " * Algorithm: %s\n" s.algo
    @printf io " * Starting Point: %s\n" string(s.initial_x)
    @printf io " * Zero: %s\n" string(s.x)
    @printf io " * Inf-norm of residuals: %f\n" norm(s.f, Inf)
    @printf io " * Convergence: %s\n" s.converged
    @printf io " * Message: %s\n" s.msg
    @printf io " * Total time: %f seconds\n" s.trace.tot_time
    @printf io " * Function Calls: %d\n" s.trace.f_calls
    @printf io " * Jacobian Calls (df/dx): %d" s.trace.g_calls
end

function fsolve(f!::Function, x0::Vector{Float64}, m::Int=length(x0); tol::Float64=1e-8,
                show_trace::Bool=false, tracing::Bool=false, method::Symbol=:hybr,
                iterations::Int=typemax(Int), kwargs...)
    if method == :hybr
        return hybrd1(f!, x0, tol, show_trace, tracing, iterations)
    elseif method == :lm
        return lmdif1(f!, x0, m, tol, show_trace, tracing, iterations)
    elseif method == :lmdif
        return lmdif(f!, x0, m, tol, show_trace, tracing, iterations; kwargs...)
    elseif method == :hybrd
        return hybrd(f!, x0, tol, show_trace, tracing, iterations; kwargs...)
    else
        error("unknown method $(method)")
    end
end

function fsolve(f!::Function, g!::Function, x0::Vector{Float64}, m::Int=length(x0);
                tol::Float64=1e-8, show_trace::Bool=false, tracing::Bool=false,
                method::Symbol=:hybr, iterations::Int=typemax(Int),
                kwargs...)
    if method == :hybr
        return hybrj(f!, g!, x0, tol, show_trace, tracing, iterations; kwargs...)
    elseif method == :lm
        return lmder(f!, g!, x0, m, tol, show_trace, tracing, iterations; kwargs...)
    else
        error("unknown method $(method)")
    end
end

include("wrappers.jl")

end  # module


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
dot.jl defines methods for the dot(x, y) function
"""

# --- Exports

import LinearAlgebra.dot
export dot

# --- Method definitions

"""
    dot(M::AbstractMultivector, N::AbstractMultivector)::AbstractMultivector
    ⋅(M::AbstractMultivector, N::AbstractMultivector)::AbstractMultivector

Compute the dot product of `M` and `N`.
"""
dot(M::AbstractMultivector, N::AbstractMultivector; left=true) =
    left ? contract_left(M, N) : nothing  # TODO

# ------ Specializations involving an AbstractMultivector instance

dot(M::AbstractMultivector, x::Real; left=true) =
    left ? contract_left(M, x) : nothing  # TODO
dot(x::Real, M::AbstractMultivector; left=true) =
    left ? contract_left(x, M) : nothing  # TODO

dot(M::AbstractMultivector, v::Vector{<:Real}; left=true) =
    left ? contract_left(M, v) : nothing  # TODO
dot(v::Vector{<:Real}, M::AbstractMultivector; left=true) =
    left ? contract_left(v, M) : nothing  # TODO


text:
module LasIO

using Base.Meta
using Dates
using FileIO
using FixedPointNumbers
using ColorTypes
using GeometryTypes # for conversion

export

    # Types
    LasHeader,
    LasVariableLengthRecord,
    LasPoint,
    LasPoint_0_5,
    LasPoint_6_10,
    LasPoint0,
    LasPoint1,
    LasPoint2,
    LasPoint3,
    LasPoint4,
    LasPoint5,
    LasPoint6,
    LasPoint7,
    LasPoint8,
    LasPoint9,
    LasPoint10,
    PointVector,

    # Functions on LasHeader
    update!,
    pointformat,
    boundingbox,

    # Functions on LasPoint
    return_number,
    number_of_returns,
    scanner_channel,
    scan_direction,
    edge_of_flight_line,
    classification,
    synthetic,
    key_point,
    withheld,
    overlap,
    xcoord,
    ycoord,
    zcoord,
    intensity,
    scan_angle,
    user_data,
    pt_src_id,
    gps_time,
    raw_classification,
    flag_byte,
    flag_byte_1,
    flag_byte_2,

    # extended from ColorTypes
    red,
    green,
    blue,
    RGB

include("vlrs.jl")
include("header.jl")
include("meta.jl")
include("point.jl")
include("util.jl")
include("fileio.jl")
include("srs.jl")

function __init__()
    # these should eventually go in
    # https://github.com/JuliaIO/FileIO.jl/blob/master/src/registry.jl
    add_format(format"LAS", "LASF", ".las", [:LasIO])
    add_format(format"LAZ", (), ".laz", [:LasIO])
    add_loader(format"LAZ", :LasIO)
    add_loader(format"LAZ", :LasIO)
end

end # module


text:
module JUDI4Cloud

import Base.+

using AzureClusterlessHPC, Reexport
@reexport using JUDI

export init_culsterless

_njpi = 1
_default_container = "mloubout/judi-cpu:1.4.3"

_judi_defaults = Dict("_POOL_ID"                => "JudiPool",
                    "_POOL_VM_SIZE"           => "Standard_E8s_v3",
                    "_VERBOSE"                => "0",
                    "_NODE_OS_OFFER"          => "ubuntu-server-container",
                    "_NODE_OS_PUBLISHER"      => "microsoft-azure-batch",
                    "_CONTAINER"              => _default_container,
                    "_NODE_COUNT_PER_POOL"    => "4",
                    "_NUM_RETRYS"             => "1",
                    "_POOL_COUNT"             => "1",
                    "_PYTHONPATH"             => "/usr/local/lib/python3.8/dist-packages",
                    "_JULIA_DEPOT_PATH"       => "/root/.julia",
                    "_OMP_NUM_THREADS"        => "4",
                    "_NODE_OS_SKU"            => "20-04-lts")

"""
Define auto scale formula to avoid idle pools
"""
auto_scale_formula(x) = """
\$TargetDedicatedNodes = $x
// Get pending tasks for the past 15 minutes.
\$samples = \$ActiveTasks.GetSamplePercent(TimeInterval_Minute * 15);
// If we have fewer than 70 percent data points, we use the last sample point, otherwise we use the maximum of last sample point and the history average.
\$tasks = \$samples < 70 ? max(0, \$ActiveTasks.GetSample(1)) :
max( \$ActiveTasks.GetSample(1), avg(\$ActiveTasks.GetSample(TimeInterval_Minute * 15)));
// If number of pending tasks is not 0, set targetVM to pending tasks, otherwise 25% of current dedicated.
\$targetVMs = \$tasks > 0 ? \$tasks : max(0, \$TargetDedicatedNodes / 4);
// The pool size is capped at NWORKERS, if target VM value is more than that, set it to NWORKERS.
cappedPoolSize = $x;
// Always start the pool at full size and keep it there for 10 minutes
\$TargetDedicatedNodes = max(0, min(\$targetVMs, cappedPoolSize));"""

len_vm(s::String) = 1
len_vm(s::Array{String, 1}) = len(s)
len_vm(s) = throw(ArgumentError("`vm_size` must be a String Array{String, 1}"))

"""
    init_culsterless(nworkers; kw...)

Initialze the serverless Azure Batch pool with `nworkers` instances/nodes

Parameters
* `nworkers`: Number of instances ine the Azure batch pool.
* `credentials`: Path to credentials file. See AzureClusterlessHPC for more details.

* `pool_name`: Name of the pool
* `vm_size`: Type of virtual machine (vm) to use for the pool. You can provide a list of vm sizes in which case it will creeate one pool per type of vm.
Check your azure batch quotas to make sure this is available for you.
* `verbose`: whether to turn on (1) or of (0) AzureClusterlessHPC verbosity.
* `nthreads`: Sets `OMP_NUM_THREADS=nthreads` on each node.
* `auto_scale`: Whether to enable autoscale. If enabled, the pool will start with no node and automatically scale based on the number of tasks.
* `n_julia_per_instance`: Number of julia worker per node (Default 1). If >1, julia will start a mini distributed setup on each node based on the number of source.
"""
function init_culsterless(nworkers=2; credentials=nothing, vm_size="Standard_E8s_v3",
                          pool_name="JudiPool", verbose=0, nthreads=4,
                          auto_scale=true, n_julia_per_instance=1,
                          container=_default_container, kw...)
    isnothing(credentials) && throw(InputError("`credentials` must be a valid path to Azure json Credentials"))
    # Check input
    npool = len_vm(vm_size)
    blob_name = lowercase("$(pool_name)tmp")
    # Update verbosity and parameters
    @eval(AzureClusterlessHPC, global __verbose__ =  Bool($verbose))
    global AzureClusterlessHPC.__params__["_NODE_COUNT_PER_POOL"] = "$(nworkers)"
    global AzureClusterlessHPC.__params__["_POOL_ID"] = "$(pool_name)"
    global AzureClusterlessHPC.__params__["_POOL_COUNT"] = "$(npool)"
    global AzureClusterlessHPC.__params__["_POOL_VM_SIZE"] = vm_size
    global AzureClusterlessHPC.__params__["_OMP_NUM_THREADS"] = "$(nthreads)"
    global AzureClusterlessHPC.__params__["_VERBOSE"] = "$(verbose)"
    global AzureClusterlessHPC.__params__["_BLOB_CONTAINER"] = blob_name
    global AzureClusterlessHPC.__params__["_CONTAINER"] = container
    # reinit everything
    isfile(credentials) || throw(FileNotFoundError(credentials))
    creds = AzureClusterlessHPC.JSON.parsefile(credentials)
    @eval(AzureClusterlessHPC, global __container__ = $blob_name)
    @eval(AzureClusterlessHPC, global __credentials__ = [$creds])
    @eval(AzureClusterlessHPC, global __resources__ = [[] for i=1:length(__credentials__)])
    @eval(AzureClusterlessHPC, global __clients__ = create_clients(__credentials__, batch=true, blob=true))

    # Create pool with idle autoscale. This will be much more efficient with a defined image rather than docker.
    if auto_scale
        create_pool(;enable_auto_scale=auto_scale,
                    auto_scale_formula=auto_scale_formula(nworkers),
                    auto_scale_evaluation_interval_minutes=5)
    else
        create_pool()
    end

    # Export JUDI on azure
    eval(macroexpand(JUDI4Cloud, quote @batchdef using Distributed, JUDI end))
    global _njpi = isnothing(n_julia_per_instance) ? 1 : n_julia_per_instance
    # Define number of local julia worker on each node in batch
    eval(macroexpand(JUDI4Cloud, quote @batchdef _nproc_loc = $_njpi end))
    include(joinpath(@__DIR__, "batch_defs.jl"))
    include(joinpath(@__DIR__, "modeling.jl"))
end

"""
    finalize_culsterless()

Finalize the clusterless job and deletes all resources (pool, tmp container, jobs)
"""
function finalize_culsterless()
    delete_all_jobs()
    delete_container()
    try delete_pool() catch; nothing end
end

function __init__()
    merge!(AzureClusterlessHPC.__params__, _judi_defaults)
    atexit(finalize_culsterless)
end

end # module


text:
###################################################
## network/tools/congestion.jl
## Create networks with particular congestion patterns
###################################################

"""
	Creates a square grid where the center is twice as congested as the rest
	Args:
		- size: length of edge of square
"""
function centerCongestion(size::Int)
	n = squareNetwork(size)
	timings = roadTypeRouting(n)
	nodeList = div(size * (size+1), 4) + squareNodes(div(size, 2), size)
	for edge in findEdges(n, nodeList)
		timings.times[src(edge), dst(edge)] *= 2
	end
	timings = NetworkTimings(n, timings.times)
	return n, timings
end

"""
	Creates a square grid where two regions (in the corners) are heavily congested (one twice as slow as the rest, the other four times as slow)
	Args:
		- size: length of edge of square
"""
function twoCongestions(size::Int)
	n = squareNetwork(size)
	timings = roadTypeRouting(n)
	nodeList1 = squareNodes(div(size, 2) - 1, size) + size + 1
	nodeList2 = nodeList1 + (div(size, 2) - 1) * (size + 1)
	for edge in findEdges(n, nodeList1)
		timings.times[src(edge), dst(edge)] *= 2
	end
	for edge in findEdges(n, nodeList2)
		timings.times[src(edge), dst(edge)] *= 4
	end
	timings = NetworkTimings(n, timings.times)
	return n, timings
end

"""
	Creates a square grid where a peripheral road is twice as slow as the rest.
	Args:
		- `size`: length of square side
"""
function squareCongestion(size::Int)
	n = squareNetwork(size)
	timings = roadTypeRouting(n)
	nodeList = div(size, 8) * (size + 1) + setdiff(squareNodes(size - div(size, 8) * 2, size), (size + 1) + squareNodes(size - div(size, 8) * 2 - 2, size))
	for edge in findEdges(n, nodeList)
		timings.times[src(edge), dst(edge)] *= 2
	end
	timings = NetworkTimings(n, timings.times)
	return n, timings
end

"""
	Creates a square grid of edge length 'size' with 4 zones in increasing congestion
"""
function gradientCongestion(size::Int)
	n = squareNetwork(size)
	g = n.graph
	times = maxSpeedTimes(n)
	for u in vertices(g)
		for v in outneighbors(g, u)
			if v > u
				times[u,v] *= div(u-1, size^2/4) + 1
			end
		end
		for v in in_neighbors(g, u)
			if v > u
				times[v,u] *= div(u-1, size^2/4) + 1
			end
		end
	end
	timings = NetworkTimings(n, times)
	return n, timings
end

"""
	Helper function to construct squares
"""
function squareNodes(edge::Int, size::Int)
	nodeList = Int[]
	for i=1:edge
		append!(nodeList, collect(1:edge) + ((i-1) * size))
	end
	return nodeList
end


text:
##
cd(@__DIR__())
using DelimitedFiles

##
path = "input.txt"
vals = readdlm(path, ' ', Int)

##
function find_triple(vals, res = 2020)
    for x in vals, y in vals, z in vals
        x + y + z == res && return x*y*z
    end
end

##
find_triple(vals)


text:
module ExtractMacroTests

using ExtractMacro
using Compat
if VERSION < v"0.7.0-DEV.1995"
    using Base.Test
else
    using Test
end

mutable struct X
    i::Int
    v::Vector{Int}
end

function test1C()
    x = X(1, [2, 3, 4])

    @extract x : i v

    @test i == 1
    @test v == [2, 3, 4]

    v[3] = 5

    @test x.v == [2, 3, 5]

    x.i = 3

    @extract x : a=v[i]

    @test a == 2

    @extract x : j=i a=v[j]

    @test j == 3
    @test a == 5

    i = 2
    @extract x : a=(v[1]+esc(i))

    @test a == 4

    @extract x : s=sum(v)+2i-1

    @test s == 15
end

function test1X()
    x = X(1, [2, 3, 4])

    @extract x

    @extract x i v

    @test i == 1
    @test v == [2, 3, 4]

    v[3] = 5

    @test x.v == [2, 3, 5]

    x.i = 3

    @extract x a=v[i]

    @test a == 2

    @extract x j=i a=v[j]

    @test j == 3
    @test a == 5

    i = 2
    @extract x a=(v[1]+esc(i))

    @test a == 4

    @extract x s=sum(v)+2i-1

    @test s == 15
end

mutable struct Y
    i::Int
    v::Vector{Int}
    z::Vector{Int}
end

function test2C()
    y = Y(1, [2,3,4], [5,6,7])

    @extract y : t=[i for i in v]
    @test t == [i for i in y.v]

    @extract y : t=[2k .+ v' .- j .- i for (k,j) in zip(v,z)]
    @test t == [2k .+ y.v' .- j .- y.i for (k,j) in zip(y.v,y.z)]

    @extract y : t=[2z .+ v' .- j .- i for (z,j) in zip(v,z)]
    @test t == [2z .+ y.v' .- j .- y.i for (z,j) in zip(y.v,y.z)]
end

function test2X()
    y = Y(1, [2,3,4], [5,6,7])

    @extract y t=[i for i in v]
    @test t == [i for i in y.v]

    @extract y t=[2k .+ v' .- j .- i for (k,j) in zip(v,z)]
    @test t == [2k .+ y.v' .- j .- y.i for (k,j) in zip(y.v,y.z)]

    @extract y t=[2z .+ v' .- j .- i for (z,j) in zip(v,z)]
    @test t == [2z .+ y.v' .- j .- y.i for (z,j) in zip(y.v,y.z)]
end

mutable struct Z
    x::X
end

function test3C()
    z = Z(X(1, [2,3,4]))

    @extract z : t=x.i
    @test t == z.x.i

    @extract z.x : t=i
    @test t == z.x.i
end

function test3X()
    z = Z(X(1, [2,3,4]))

    @extract z t=x.i
    @test t == z.x.i

    @extract z.x t=i
    @test t == z.x.i
end

macro test_extract_failure(ex...)
    VERSION ≥ v"0.7-DEV.357" && (ex = [nothing, ex...])
    ex = Expr(:call, :macroexpand, @__MODULE__, Expr(:quote, Expr(:macrocall, Symbol("@extract"), ex...)))
    if VERSION ≥ v"0.7-DEV.1676"
        quote
            @test_throws LoadError $ex
        end
    else
        quote
            ex = $ex
            @test Meta.isexpr(ex, :error)
            @test isa(ex.args[1], ErrorException)
        end
    end
end

function test4C()
    @test_extract_failure x : abs(v)
    @test_extract_failure x : v'
    @test_extract_failure x : v[1]
end

function test4X()
    @test_extract_failure x abs(v)
    @test_extract_failure x v'
    @test_extract_failure x v[1]
end

test1C()
test1X()
test2C()
test2X()
test3C()
test3X()
test4C()
test4X()

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 102adb52-8b83-11eb-1876-a52236a70e1f
using Images,MosaicViews

# ╔═╡ 4f027e0c-8b88-11eb-2f3f-237e69a18582
begin
	using TestImages, ImageContrastAdjustment, HistogramThresholding

img = testimage("cameraman")
end

# ╔═╡ 5bd1efe0-8bf2-11eb-3f4c-6d5a43febc6d
begin
edges, counts = build_histogram(img,256)
t = find_threshold(Otsu(), counts[1:end], edges)
# Apply the threshold to img and visualize the result.
img_final = zeros(Gray{Bool}, axes(img))
for i in CartesianIndices(img)
        img_final[i] = img[i] < t ? 0 : 1
end
img_final
end

# ╔═╡ Cell order:
# ╠═102adb52-8b83-11eb-1876-a52236a70e1f
# ╠═4f027e0c-8b88-11eb-2f3f-237e69a18582
# ╠═5bd1efe0-8bf2-11eb-3f4c-6d5a43febc6d

max_stars_repo_name: ashwani-rathee/Demos

max_issues_repo_name: ashwani-rathee/Demos

max_forks_repo_name: ashwani-rathee/Demos


text:
using LinearAlgebra

"""
This file contains a list of cost functions that are ready to be used in optimisations
I'm going to try and keep some standard inputs that make sense

Needs citation from the Schuster lab paper
"""

"""
Density matrix fidelity
"""
function C1(KT, KN)
    D = size(KT)[1]
    1 - abs2(tr(KT' * KN) / D)
    # 1 - real(tr(KT' * KN))/D
end

"""
Target-state infidelity
"""
function C2(ψT, ψN)
    1 - abs2(ψT' * ψN)
end

"""
Control amplitudes
"""
function C3(u)
    sum(abs2.(u))
end

"""
Control variations
"""
function C4(u)
    (K, J) = size(u)
    sum(abs2.(diff(u, dims = 2)))
end

"""
Occupation of forbidden states
"""
function C5(ψF, ψj)
    # since ψj is the intermediate states we can assume its an array
    sum(abs2(ψF' * ψ) for ψ in ψj)
end

"""
Evolution time (target gate)
"""
function C6(KT, KJ, N, D)
    1 - 1 / N * sum(abs2(tr(KT' * Kj) / D) for Kj in KJ)
end

"""
Evolution time (target state)
"""
function C7(ψT, ψJ, N)
    1 - 1 / N * sum(abs2(tr(ψT * ψj)) for ψj in ψJ)
end

"""
Stores information, relative weights and functions to compute the penalty functions
"""
struct PenaltyFunctionals
    weights::Any # relative weights of each
    functions::Any # array of penalty functionals
end

"""
Cost functions for the GRAPE algorithm
"""

"""
Unitary synthesis figure of merit, from Khaneja et. al. paper (proper citation needed)
"""
function fom_func(system_type::UnitaryGate, t, k, state, costate, props)::Float64
    @views rael(tr(costate[t, k]' * state[t, k]) * tr(state[t, k]' * costate[t, k]))
end

function fom_func(
    system_type::Union{StateTransfer,CoherenceTransfer},
    t,
    k,
    state,
    costate,
    props,
)::Float64
    # recall that target is always the last entry of L
    # and that we have in U[end] the propagated forward target state
    @views real(C1(costate[t, k], state[t, k]))
end

function fom_func(system_type::UnitaryGate, t, state, costate, props)::Float64
    @views real(tr(state[t]' * costate[t]) * tr(state[t]' * costate[t]))
end

function fom_func(
    system_type::Union{StateTransfer,CoherenceTransfer},
    t,
    state,
    costate,
    props,
)::Float64
    @views real(C1(costate[t], state[t]))
end


text:
using LinearAlgebra
using SparseArrays
using Plots
using ForwardDiff

"This routine solves Poisson's equation using piecewise linear finite element methods.
The matrix is assembled locally, and boundary nodes are not removed."

m = 100 # number of elements

α = 0.0
β = 0.0
f(x) = sin(2*pi*x)
κ(x) = 1.0 + 2.0*(x > 0)

# define spatial grid
x = LinRange(-1,1,m+2) # x_0, x_1, ..., x_m, x_{m+1} = x_0
xint = x[1:end-1]
h = x[2]-x[1]

# x = x + randn(size(x))/(2*m)
# x = @. x + (1+x)*(1-x)/3

# construct local FEM matrix
map_point(x,a,b) = a + (b-a) * x # maps x ∈ [0,1] to interval [a,b]
A = spzeros(m+2,m+2)
b = zeros(m+2)
A_local = [1 -1;-1 1]
for e = 1:m+1
    h_e = x[e+1] - x[e]
    ids = e:e+1

    # use 1-point quadrature to approximate integrals
    x_e = map_point(0, x[e], x[e+1])
    w_e = 2.0

    # accumulate local contributions
    A[ids,ids] += A_local * κ(x_e) / h_e
    b[ids] += h_e/2 * w_e * f(x_e) * [.5;.5]
end

# modify for boundary conditions
function impose_Dirichlet_BC!(A,b,i,val)
    b .-= A[:,i]*val
    b[i] = val
    A[i,:] .= 0
    A[:,i] .= 0
    A[i,i] = 1.0
end
function impose_Neumann_BC!(A,b,i,val)
    b[i] += val
end

impose_Dirichlet_BC!(A,b,1,α)
impose_Dirichlet_BC!(A,b,m+2,β)

u = A\b

# uex = @. (1+x)*(1-x)/2
plot(x,u,legend=false,mark=:dot,ms=2)


text:
using Base: Int64, Float64, NamedTuple
using Printf
using Parameters        # helps setting default parameters in structures
using SpecialFunctions: erfc

# LaMEM_geometry
#
# These are routines that help to create LaMEM input geometries, such as slabs with a given angle
#

export  AddBox!, AddSphere!, AddEllipsoid!, AddCylinder!,
        ConstantTemp, LinearTemp, HalfspaceCoolingTemp, SpreadingRateTemp,
        ConstantPhase, LithosphericPhases,
        Compute_ThermalStructure, Compute_Phase

"""
    AddBox!(Phase, Temp, Grid::LaMEM_grid; xlim=Tuple{2}, [ylim=Tuple{2}], zlim=Tuple{2},
            Origin=nothing, StrikeAngle=0, DipAngle=0,
            phase = ConstantPhase(1),
            T=nothing )

Adds a box with phase & temperature structure to a 3D model setup.  This simplifies creating model geometries in geodynamic models

Parameters
====
- Phase - Phase array (consistent with Grid)
- Temp  - Temperature array (consistent with Grid)
- Grid - LaMEM grid structure (usually obtained with ReadLaMEM_InputFile)
- xlim - left/right coordinates of box
- ylim - front/back coordinates of box [optional; if not specified we use the whole box]
- zlim - bottom/top coordinates of box
- Origin - the origin, used to rotate the box around. Default is the left-front-top corner
- StrikeAngle - strike angle of slab
- DipAngle - dip angle of slab
- phase - specifies the phase of the box. See `ConstantPhase()`,`LithosphericPhases()`
- T - specifies the temperature of the box. See `ConstantTemp()`,`LinearTemp()`,`HalfspaceCoolingTemp()`,`SpreadingRateTemp()`

Examples
========

Example 1) Box with constant phase and temperature & a dip angle of 10 degrees:
```julia
julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> AddBox!(Phases,Temp,Grid, xlim=(0,500), zlim=(-50,0), phase=ConstantPhase(3), DipAngle=10, T=ConstantTemp(1000))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> Write_Paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"
```

Example 2) Box with halfspace cooling profile
```julia
julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> AddBox!(Phases,Temp,Grid, xlim=(0,500), zlim=(-50,0), phase=ConstantPhase(3), DipAngle=10, T=ConstantTemp(1000))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> Write_Paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"
```
"""
function AddBox!(Phase, Temp, Grid::LaMEM_grid;                 # required input
                xlim=Tuple{2}, ylim=nothing, zlim=Tuple{2},     # limits of the box
                Origin=nothing, StrikeAngle=0, DipAngle=0,      # origin & dip/strike
                phase = ConstantPhase(1),                       # Sets the phase number(s) in the box
                T=nothing )                                     # Sets the thermal structure (various fucntions are available)

    # Limits of block
    if ylim==nothing
        ylim = (Grid.coord_y...,)
    end

    if Origin==nothing
        Origin = (xlim[1], ylim[1], zlim[2])  # upper-left corner
    end

    # Perform rotation of 3D coordinates:
    Xrot = Grid.X .- Origin[1];
    Yrot = Grid.Y .- Origin[2];
    Zrot = Grid.Z .- Origin[3];

    Rot3D!(Xrot,Yrot,Zrot, StrikeAngle, DipAngle)

    # Set phase number & thermal structure in the full domain
    ztop = zlim[2] - Origin[3]
    zbot = zlim[1] - Origin[3]
    ind = findall(  (Xrot .>= (xlim[1] - Origin[1])) .& (Xrot .<= (xlim[2] - Origin[1])) .&
                    (Yrot .>= (ylim[1] - Origin[2])) .& (Yrot .<= (ylim[2] - Origin[2])) .&
                    (Zrot .>= zbot) .& (Zrot .<= ztop)  )

    # Compute thermal structure accordingly. See routines below for different options
    if T != nothing
        Temp[ind] = Compute_ThermalStructure(Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], T)
    end

    # Set the phase. Different routines are available for that - see below.
    Phase[ind] = Compute_Phase(Phase[ind], Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], phase)

    return nothing
end

"""
    AddSphere!(Phase, Temp, Grid::LaMEM_grid; cen=Tuple{3}, radius=Tuple{1},
            phase = ConstantPhase(1).
            T=nothing )

Adds a sphere with phase & temperature structure to a 3D model setup.  This simplifies creating model geometries in geodynamic models

Parameters
====
- Phase - Phase array (consistent with Grid)
- Temp  - Temperature array (consistent with Grid)
- Grid - LaMEM grid structure (usually obtained with ReadLaMEM_InputFile)
- cen - center coordinates of sphere
- radius - radius of sphere
- phase - specifies the phase of the box. See `ConstantPhase()`,`LithosphericPhases()`
- T - specifies the temperature of the box. See `ConstantTemp()`,`LinearTemp()`,`HalfspaceCoolingTemp()`,`SpreadingRateTemp()`

Example
========

Sphere with constant phase and temperature:
```julia
julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> AddSphere!(Phases,Temp,Grid, cen=(0,0,-1), radius=0.5, phase=ConstantPhase(2), T=ConstantTemp(800))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> Write_Paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"
```
"""
function AddSphere!(Phase, Temp, Grid::LaMEM_grid;      # required input
    cen=Tuple{3}, radius=Tuple{1},                         # center and radius of the sphere
    phase = ConstantPhase(1),                           # Sets the phase number(s) in the sphere
    T=nothing )                                         # Sets the thermal structure (various fucntions are available)

    # Set phase number & thermal structure in the full domain
    ind = findall(((Grid.X .- cen[1]).^2 + (Grid.Y .- cen[2]).^2 + (Grid.Z .- cen[3]).^2).^0.5 .< radius)

    # Compute thermal structure accordingly. See routines below for different options
    if T != nothing
        Temp[ind] = Compute_ThermalStructure(Temp[ind], Grid.X[ind], Grid.Y[ind], Grid.Z[ind], T)
    end

    # Set the phase. Different routines are available for that - see below.
    Phase[ind] = Compute_Phase(Phase[ind], Temp[ind], Grid.X[ind], Grid.Y[ind], Grid.Z[ind], phase)

    return nothing
end

"""
    AddEllipsoid!(Phase, Temp, Grid::LaMEM_grid; cen=Tuple{3}, axes=Tuple{3},
            Origin=nothing, StrikeAngle=0, DipAngle=0,
            phase = ConstantPhase(1).
            T=nothing )

Adds an Ellipsoid with phase & temperature structure to a 3D model setup.  This simplifies creating model geometries in geodynamic models

Parameters
====
- Phase - Phase array (consistent with Grid)
- Temp  - Temperature array (consistent with Grid)
- Grid - LaMEM grid structure (usually obtained with ReadLaMEM_InputFile)
- cen - center coordinates of sphere
- axes - semi-axes of ellipsoid in X,Y,Z
- Origin - the origin, used to rotate the box around. Default is the left-front-top corner
- StrikeAngle - strike angle of slab
- DipAngle - dip angle of slab
- phase - specifies the phase of the box. See `ConstantPhase()`,`LithosphericPhases()`
- T - specifies the temperature of the box. See `ConstantTemp()`,`LinearTemp()`,`HalfspaceCoolingTemp()`,`SpreadingRateTemp()`

Example
========

Ellipsoid with constant phase and temperature, rotated 90 degrees and tilted by 45 degrees:
```julia
julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> AddEllipsoid!(Phases,Temp,Grid, cen=(-1,-1,-1), axes=(0.2,0.1,0.5), StrikeAngle=90, DipAngle=45, phase=ConstantPhase(3), T=ConstantTemp(600))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> Write_Paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"
```
"""
function AddEllipsoid!(Phase, Temp, Grid::LaMEM_grid;      # required input
    cen=Tuple{3}, axes=Tuple{3},                           # center and semi-axes of the ellpsoid
    Origin=nothing, StrikeAngle=0, DipAngle=0,             # origin & dip/strike
    phase = ConstantPhase(1),                              # Sets the phase number(s) in the box
    T=nothing )                                            # Sets the thermal structure (various fucntions are available)

    if Origin==nothing
        Origin = cen  # center
    end

    # Perform rotation of 3D coordinates:
    Xrot = Grid.X .- Origin[1];
    Yrot = Grid.Y .- Origin[2];
    Zrot = Grid.Z .- Origin[3];

    Rot3D!(Xrot,Yrot,Zrot, StrikeAngle, DipAngle)

    # Set phase number & thermal structure in the full domain
    x2     = axes[1]^2
    y2     = axes[2]^2
    z2     = axes[3]^2
    cenRot = cen .- Origin
    ind = findall((((Xrot .- cenRot[1]).^2)./x2 + ((Yrot .- cenRot[2]).^2)./y2 +
                   ((Zrot .- cenRot[3]).^2)./z2) .^0.5 .<= 1)

    # Compute thermal structure accordingly. See routines below for different options
    if T != nothing
        Temp[ind] = Compute_ThermalStructure(Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], T)
    end

    # Set the phase. Different routines are available for that - see below.
    Phase[ind] = Compute_Phase(Phase[ind], Temp[ind], Xrot[ind], Yrot[ind], Zrot[ind], phase)

    return nothing
end

"""
    AddCylinder!(Phase, Temp, Grid::LaMEM_grid; base=Tuple{3}, cap=Tuple{3}, radius=Tuple{1},
            phase = ConstantPhase(1).
            T=nothing )

Adds a cylinder with phase & temperature structure to a 3D model setup.  This simplifies creating model geometries in geodynamic models

Parameters
====
- Phase - Phase array (consistent with Grid)
- Temp  - Temperature array (consistent with Grid)
- Grid - LaMEM grid structure (usually obtained with ReadLaMEM_InputFile)
- base - center coordinate of bottom of cylinder
- cap - center coordinate of top of cylinder
- radius - radius of the cylinder
- phase - specifies the phase of the box. See `ConstantPhase()`,`LithosphericPhases()`
- T - specifies the temperature of the box. See `ConstantTemp()`,`LinearTemp()`,`HalfspaceCoolingTemp()`,`SpreadingRateTemp()`

Example
========

Cylinder with constant phase and temperature:
```julia
julia> Grid = ReadLaMEM_InputFile("test_files/SaltModels.dat")
LaMEM Grid:
  nel         : (32, 32, 32)
  marker/cell : (3, 3, 3)
  markers     : (96, 96, 96)
  x           ϵ [-3.0 : 3.0]
  y           ϵ [-2.0 : 2.0]
  z           ϵ [-2.0 : 0.0]
julia> Phases = zeros(Int32,   size(Grid.X));
julia> Temp   = zeros(Float64, size(Grid.X));
julia> AddCylinder!(Phases,Temp,Grid, base=(-1,-1,-1.5), cap=(1,1,-0.5), radius=0.25, phase=ConstantPhase(4), T=ConstantTemp(400))
julia> Model3D = ParaviewData(Grid, (Phases=Phases,Temp=Temp)); # Create Cartesian model
julia> Write_Paraview(Model3D,"LaMEM_ModelSetup")           # Save model to paraview
1-element Vector{String}:
 "LaMEM_ModelSetup.vts"
```
"""
function AddCylinder!(Phase, Temp, Grid::LaMEM_grid;    # required input
    base=Tuple{3}, cap=Tuple{3}, radius=Tuple{1},       # center and radius of the sphere
    phase = ConstantPhase(1),                           # Sets the phase number(s) in the sphere
    T=nothing )                                         # Sets the thermal structure (various fucntions are available)

    # axis vector of cylinder
    axVec = cap .- base
    ax2   = (axVec[1]^2 + axVec[2]^2 + axVec[3]^2)

    # distance between grid points and cylinder base
    dx_b  = Grid.X .- base[1]
    dy_b  = Grid.Y .- base[2]
    dz_b  = Grid.Z .- base[3]

    # find normalized parametric coordinate of a point-axis projection
    t     = (axVec[1] .* dx_b .+ axVec[2] .* dy_b .+ axVec[3] .* dz_b) ./ ax2

    # find distance vector between point and axis
    dx    = dx_b .- t.*axVec[1]
    dy    = dy_b .- t.*axVec[2]
    dz    = dz_b .- t.*axVec[3]

    # Set phase number & thermal structure in the full domain
    ind = findall((t .>= 0.0) .& (t .<= 1.0) .& ((dx.^2 + dy.^2 + dz.^2).^0.5 .<= radius))

    # Compute thermal structure accordingly. See routines below for different options
    if T != nothing
        Temp[ind] = Compute_ThermalStructure(Temp[ind], Grid.X[ind], Grid.Y[ind], Grid.Z[ind], T)
    end

    # Set the phase. Different routines are available for that - see below.
    Phase[ind] = Compute_Phase(Phase[ind], Temp[ind], Grid.X[ind], Grid.Y[ind], Grid.Z[ind], phase)

    return nothing
end

# Internal function that rotates the coordinates
function Rot3D!(X,Y,Z, StrikeAngle, DipAngle)

    # rotation matrixes
    roty = [cosd(-DipAngle) 0 sind(-DipAngle) ; 0 1 0 ; -sind(-DipAngle) 0  cosd(-DipAngle)];
    rotz = [cosd(StrikeAngle) -sind(StrikeAngle) 0 ; sind(StrikeAngle) cosd(StrikeAngle) 0 ; 0 0 1]

    for i in eachindex(X)
        CoordVec = [X[i], Y[i], Z[i]]
        CoordRot =  rotz*CoordVec;
        CoordRot =  roty*CoordRot;
        X[i] = CoordRot[1];
        Y[i] = CoordRot[2];
        Z[i] = CoordRot[3];
    end

    return nothing
end

abstract type AbstractThermalStructure end

"""
    ConstantTemp(T=1000)

Sets a constant temperature inside the box

Parameters
===
- T : the value
"""
@with_kw_noshow mutable struct ConstantTemp <: AbstractThermalStructure
    T = 1000
end

function Compute_ThermalStructure(Temp, X, Y, Z, s::ConstantTemp)
    Temp .= s.T
    return Temp
end

"""
    LinearTemp(Ttop=0, Tbot=1000)

Set a linear temperature structure from top to bottom

Parameters
===
- Ttop : the value @ the top
- Tbot : the value @ the bottom

"""
@with_kw_noshow mutable struct LinearTemp <: AbstractThermalStructure
    Ttop = 0
    Tbot = 1350
end

function Compute_ThermalStructure(Temp, X, Y, Z, s::LinearTemp)
    @unpack Ttop, Tbot  = s

    dz   = Z[end]-Z[1];
    dT   = Tbot - Ttop

    Temp = abs.(Z./dz).*dT .+ Ttop
    return Temp
end

"""
    HalfspaceCoolingTemp(Tsurface=0, Tmantle=1350, Age=60, Adiabat=0)

Sets a halfspace temperature structure in plate

Parameters
========
- Tsurface : surface temperature [C]
- Tmantle : mantle temperature [C]
- Age : Thermal Age of plate [Myrs]
- Adiabat : Mantle Adiabat [K/km]

"""
@with_kw_noshow mutable struct HalfspaceCoolingTemp <: AbstractThermalStructure
    Tsurface = 0       # top T
    Tmantle = 1350     # bottom T
    Age  = 60          # thermal age of plate [in Myrs]
    Adiabat = 0        # Adiabatic gradient in K/km
end

function Compute_ThermalStructure(Temp, X, Y, Z, s::HalfspaceCoolingTemp)
    @unpack Tsurface, Tmantle, Age, Adiabat  = s

    kappa       =   1e-6;
    SecYear     =   3600*24*365
    dz          =   Z[end]-Z[1];
    ThermalAge  =   Age*1e6*SecYear;

    MantleAdiabaticT    =   Tmantle .+ Adiabat*abs.(Z);   # Adiabatic temperature of mantle

    for i in eachindex(Temp)
        Temp[i] =   (Tsurface .- Tmantle)*erfc((abs.(Z[i])*1e3)./(2*sqrt(kappa*ThermalAge))) + MantleAdiabaticT[i];
    end
    return Temp
end

"""
    SpreadingRateTemp(Tsurface=0, Tmantle=1350, Adiabat=0, MORside="left",SpreadingVel=3, AgeRidge=0, maxAge=80)

Sets a halfspace temperature structure within the box, combined with a spreading rate (which implies that the plate age varies)

Parameters
========
- Tsurface : surface temperature [C]
- Tmantle : mantle temperature [C]
- Adiabat : Mantle Adiabat [K/km]
- MORside : side of the box where the MOR is located ["left","right","front","back"]
- SpreadingVel : spreading velocity [cm/yr]
- AgeRidge : thermal age of the ridge [Myrs]
- maxAge : maximum thermal Age of plate [Myrs]

"""
@with_kw_noshow mutable struct SpreadingRateTemp <: AbstractThermalStructure
    Tsurface = 0       # top T
    Tmantle = 1350     # bottom T
    Adiabat = 0        # Adiabatic gradient in K/km
    MORside = "left"   # side of box where the MOR is located
    SpreadingVel = 3   # spreading velocity [cm/yr]
    AgeRidge = 0       # Age of the ridge [Myrs]
    maxAge  = 60       # maximum thermal age of plate [Myrs]
end

function Compute_ThermalStructure(Temp, X, Y, Z, s::SpreadingRateTemp)
    @unpack Tsurface, Tmantle, Adiabat, MORside, SpreadingVel, AgeRidge, maxAge  = s

    kappa       =   1e-6;
    SecYear     =   3600*24*365
    dz          =   Z[end]-Z[1];

    MantleAdiabaticT    =   Tmantle .+ Adiabat*abs.(Z);   # Adiabatic temperature of mantle

    if MORside=="left"
        Distance = X .- X[1,1,1];
    elseif MORside=="right"
        Distance = X[end,1,1] .- X;
    elseif MORside=="front"
        Distance = Y .- Y[1,1,1];
    elseif MORside=="back"
        Distance = Y[1,end,1] .- Y;
    else
        error("unknown side")
    end

    for i in eachindex(Temp)
        ThermalAge    =   abs(Distance[i]*1e3*1e2)/SpreadingVel + AgeRidge*1e6;   # Thermal age in years
        if ThermalAge>maxAge*1e6
            ThermalAge = maxAge*1e6
        end

        ThermalAge    =   ThermalAge*SecYear;

        Temp[i] = (Tsurface .- Tmantle)*erfc((abs.(Z[i])*1e3)./(2*sqrt(kappa*ThermalAge))) + MantleAdiabaticT[i];
    end
    return Temp
end

abstract type AbstractPhaseNumber end

"""
    ConstantPhase(phase=1)

Sets a constant phase inside the box

Parameters
===
- phase : the value
"""
@with_kw_noshow mutable struct ConstantPhase <: AbstractPhaseNumber
    phase = 1
end

function Compute_Phase(Phase, Temp, X, Y, Z, s::ConstantPhase)
    Phase .= s.phase
    return Phase
end

"""
    LithosphericPhases(Layers=[10 20 30], Phases=[1 2 3 4], Tlab=nothing )

This allows defining a layered lithosphere. Layering is defined from the top downwards.

Parameters
===
- Layers : the thickness of the layers from top downwards
- Phases : the phases of the layers from top down. Note that this array
- Tlab   : Temperature of the lithosphere asthenosphere boundary. If specified, the phases at locations with T>Tlab is set to Phases[end]

"""
@with_kw_noshow mutable struct LithosphericPhases <: AbstractPhaseNumber
    Layers  = [10., 20., 50.]
    Phases  = [1,   2  , 3,  4]
    Tlab    = nothing
end

function Compute_Phase(Phase, Temp, X, Y, Z, s::LithosphericPhases)
    @unpack Layers, Phases, Tlab  = s

    Phase .= Phases[end]
    Ztop  = 0
    for i=1:length(Layers)
        Zbot = Ztop-Layers[i]
        ind = findall( ( Z .>= Zbot) .&  (Z .<= Ztop) );
        Phase[ind] .= Phases[i]

        Ztop = Zbot
    end

    # set phase to mantle if requested
    if Tlab != nothing
        ind = findall(Temp .> Tlab)
        Phase[ind] .= Phases[end]
    end

    return Phase
end


text:
#= The famous Tower of Hanoi question where there are 3 towers with discs in them. You need
to transfer all the discs from A to C =#

## Function

function towerOfHanoi(discs, source = "A", destination = "C", helper = "B")
    if (discs == 0)
        return
    end
    towerOfHanoi(discs - 1, source, helper, destination)
    println("Move from $(source) to $(destination)")
    towerOfHanoi(discs - 1, helper, destination, source)
end

## Input

println("Enter the number of discs: ")
n = readline()
n = parse(Int64, n)

## Calling the function

towerOfHanoi(n)

#=
Sample Test Case:
Input:
    Enter the number of discs: 3
Output:
    Move from A to C
    Move from A to B
    Move from C to B
    Move from A to C
    Move from B to A
    Move from B to C
    Move from A to C
Time complexity: O(2^N)
=#

max_stars_repo_name: TechSpiritSS/NeoAlgo


text:
#using BinDeps

#if !(is_linux())
#	error("currently, SIT only supports Linux")
#end
#
## get the current path
#currentFilePath = @__FILE__()
#currentDirPath = dirname(currentFilePath)
#
#srcf90dir=joinpath(currentDirPath, "src", "f90")
#makefile=joinpath(currentDirPath,"makefile")
#
#
## use mkdir to create usr/lib
#function mklibdir()
#	usrdir = joinpath(currentDirPath, "usr");
#	if isdir(usrdir) == false
#		mkdir(usrdir);
#	end
#
#	libdir = joinpath(usrdir, "lib");
#	if isdir(libdir) == false
#		mkdir(libdir);
#	end
#end
#
#
#f90cmd=string("ifort -assume byterecl  -fp-model fast=2 -Ofast  -autodouble -fPIC -shared")
#f90cmd=string("gfortran -fopenmp   -O3  -fdefault-real-8 -fdefault-double-8 -fPIC -shared")
#
#libsrcvec=[
#	"precision_mod.f90", "string_routines.f90 precision_mod.f90",
#	"error_messages.f90 precision_mod.f90 string_routines.f90",
#	"unique.f90 error_messages.f90",
#	"sorting.f90",
#	"io.f90 precision_mod.f90 error_messages.f90 unique.f90 sorting.f90 string_routines.f90",
#	"test.f90 precision_mod.f90 string_routines.f90"
#	]
#
#libsrcvec=[
#	"precision_mod.f90", "string_routines.f90",
#	"error_messages.f90",
#	"unique.f90",
#	"sorting.f90",
#	"io.f90",
#	"test.f90"
#	]
#
#libsrcname=broadcast(libsrcvec) do x
#	joinpath(srcf90dir, x)
#end
#
#
#libsovec=[
#	"precision_mod.so",
#	"string_routines.so",
#	"error_messages.so",
#	"unique.so",
#	"sorting.so",
#	"io.so",
#	"test.so"
#	]
#
#libsovec=["lib.so"]
#
#function compile()
#	for ilib in 1:length(libsovec)
#		libsrcname=joinpath(srcf90dir,libsrcvec[ilib])
#		libsoname = joinpath(currentDirPath, "usr", "lib", libsovec[ilib])
#		run(`gfortran -O3  -fdefault-real-8 -fdefault-double-8 -fPIC -shared $libsrcname -o $libsoname`)
#		println(libsrcname, libsoname)
#
#		outputfile = open(joinpath(currentDirPath, "deps.jl"), "w");
#		write(outputfile, string("@checked_lib ", libsovec[ilib], " ", libsoname))
#		close(outputfile)
#	end
#end
#
#
#function writeDeps()
#	libpath = joinpath(currentDirPath, "usr", "lib", "liblbfgsbf.so")
#	outputfile = open(joinpath(currentDirPath, "deps.jl"), "w");
#	write(outputfile, "macro checked_lib(libname, path)\n")
#	write(outputfile, "(dlopen_e(path) == C_NULL) && \(error(\"Unable to load \$libname (\$path)")
#	write(outputfile, "Please re-run Pkg.build(package), and restart Julia.\")\n \)")
#        write(outputfile, "quote const \$(esc(libname)) = \$path end\nend\n")
#	close(outputfile)
#end
#
#
##mklibdir()
##compile()
#writeDeps()