text:
# https://github.com/mbaz/Gaston.

# --------------------------------------------
# These functions are called by Plots
# --------------------------------------------

# Create the window/figure for this backend.
function _create_backend_figure(plt::Plot{GastonBackend})
    state_handle = Gaston.nexthandle() # for now all the figures will be kept
    plt.o = Gaston.newfigure(state_handle)
end

function _before_layout_calcs(plt::Plot{GastonBackend})
    # Initialize all the subplots first
    plt.o.subplots = Gaston.SubPlot[]

    n1 = n2 = 0
    if length(plt.inset_subplots) > 0
        n1, sps = gaston_get_subplots(0, plt.inset_subplots, plt.layout)
        gaston_init_subplots(plt, sps)
    end

    if length(plt.subplots) > 0
        n2, sps = gaston_get_subplots(0, plt.subplots, plt.layout)
    end

    if (n = n1 + n2) != length(plt.subplots)
        @error "Gaston: $n != $(length(plt.subplots))"
    end

    plt.o.layout = gaston_init_subplots(plt, sps)

    # Then add the series (curves in gaston)
    for series ∈ plt.series_list
        gaston_add_series(plt, series)
    end

    for sp ∈ plt.subplots
        sp === nothing && continue
        for ann in sp[:annotations]
            x, y, val = locate_annotation(sp, ann...)
            sp.o.axesconf *= "\nset label '$(val.str)' at $x,$y $(gaston_font(val.font))"
        end
    end
    nothing
end

function _update_min_padding!(sp::Subplot{GastonBackend})
    sp.minpad = 0mm, 0mm, 0mm, 0mm
    nothing
end

function _update_plot_object(plt::Plot{GastonBackend})
    # respect the layout ratio
    dat = gaston_multiplot_pos_size(plt.layout, (0, 0, 1, 1))
    gaston_multiplot_pos_size!(dat)
    nothing
end

for (mime, term) ∈ (
    "application/eps"        => "epscairo",
    "image/eps"              => "epslatex",
    "application/pdf"        => "pdfcairo",
    "application/postscript" => "postscript",
    "image/png"              => "png",
    "image/svg+xml"          => "svg",
    "text/latex"             => "tikz",
    "application/x-tex"      => "epslatex",
    "text/plain"             => "dumb",
)
    @eval function _show(io::IO, ::MIME{Symbol($mime)}, plt::Plot{GastonBackend})
        term = String($term); tmpfile = "$(Gaston.tempname()).$term"

        Gaston.save(term=term, output=tmpfile, handle=plt.o.handle, saveopts=gaston_saveopts(plt))
        while !isfile(tmpfile) end  # avoid race condition with read in next line
        write(io, read(tmpfile))
        rm(tmpfile, force=true)
        nothing
    end
end

_display(plt::Plot{GastonBackend}) = display(plt.o)

# --------------------------------------------
# These functions are gaston specific
# --------------------------------------------

function gaston_saveopts(plt::Plot{GastonBackend})
    saveopts = String["size $(join(plt.attr[:size], ","))"]

    push!(saveopts, gaston_font(plottitlefont(plt), rot=false, align=false, color=false, scale=1))

    push!(saveopts, "background $(gaston_color(plt.attr[:background_color]))")

    # push!(saveopts, "title '$(plt.attr[:window_title])'")

    # Scale all plot elements to match Plots.jl DPI standard
    scaling = plt.attr[:dpi] / Plots.DPI
    push!(saveopts, "fontscale $scaling lw $scaling dl $scaling ps $scaling")

    return join(saveopts, " ")
end

function gaston_get_subplots(n, plt_subplots, layout)
    nr, nc = size(layout)
    sps = Array{Any}(nothing, nr, nc)
    for r ∈ 1:nr, c ∈ 1:nc  # NOTE: col major
        l = layout[r, c]
        if l isa GridLayout
            n, sub = gaston_get_subplots(n, plt_subplots, l)
            sps[r, c] = size(sub) == (1, 1) ? only(sub) : sub
        else
            sps[r, c] = get(l.attr, :blank, false) ? nothing : plt_subplots[n += 1]
        end
    end
    return n, sps
end

function gaston_init_subplots(plt, sps)
    sz = nr, nc = size(sps)
    for c ∈ 1:nc, r ∈ 1:nr  # NOTE: row major
        sp = sps[r, c]
        if sp isa Subplot || sp === nothing
            gaston_init_subplot(plt, sp)
        else
            gaston_init_subplots(plt, sp)
            sz = max.(sz, size(sp))
        end
    end
    return sz
end

function gaston_init_subplot(plt::Plot{GastonBackend}, sp::Union{Nothing,Subplot{GastonBackend}})
    if sp === nothing
        push!(plt.o.subplots, sp)
    else
        dims = RecipesPipeline.is3d(sp) || sp.attr[:projection] == "3d" || needs_any_3d_axes(sp) ? 3 : 2
        any_label = false
        for series ∈ series_list(sp)
            if dims == 2 && series[:seriestype] ∈ (:heatmap, :contour)
                dims = 3  # we need heatmap/contour to use splot, not plot
            end
            any_label |= should_add_to_legend(series)
        end
        sp.o = Gaston.Plot(
            dims=dims, curves=[], axesconf=gaston_parse_axes_args(plt, sp, dims, any_label),
        )
        push!(plt.o.subplots, sp.o)
    end
    nothing
end

function gaston_multiplot_pos_size(layout, parent_xy_wh)
    nr, nc = size(layout)
    dat = Array{Any}(nothing, nr, nc)
    for r ∈ 1:nr, c ∈ 1:nc
        l = layout[r, c]
        # width and height (pct) are multiplicative (parent)
        w = layout.widths[c].value * parent_xy_wh[3]
        h = layout.heights[r].value * parent_xy_wh[4]
        if isa(l, EmptyLayout)
            dat[r, c] = (c - 1) * w, (r - 1) * h, w, h, nothing
        else
            # previous position (origin)
            prev_r = r > 1 ? dat[r - 1, c] : nothing
            prev_c = c > 1 ? dat[r, c - 1] : nothing
            prev_r isa Array && (prev_r = prev_r[end, end])
            prev_c isa Array && (prev_c = prev_c[end, end])
            x = prev_c !== nothing ? prev_c[1] + prev_c[3] : parent_xy_wh[1]
            y = prev_r !== nothing ? prev_r[2] + prev_r[4] : parent_xy_wh[2]
            if l isa GridLayout
                sub = gaston_multiplot_pos_size(l, (x, y, w, h))
                dat[r, c] = size(sub) == (1, 1) ? only(sub) : sub
            else
                dat[r, c] = x, y, w, h, l
            end
        end
    end
    return dat
end

function gaston_multiplot_pos_size!(dat)
    nr, nc = size(dat)
    for r ∈ 1:nr, c ∈ 1:nc
        xy_wh_sp = dat[r, c]
        if xy_wh_sp isa Array
            gaston_multiplot_pos_size!(xy_wh_sp)
        elseif xy_wh_sp isa Tuple
            x, y, w, h, sp = xy_wh_sp
            sp === nothing && continue
            sp.o === nothing && continue
            # gnuplot screen coordinates: bottom left at 0,0 and top right at 1,1
            sp.o.axesconf = "set origin $x, $(1 - y - h)\nset size $w, $h\n" * sp.o.axesconf
        end
    end
    nothing
end

function gaston_add_series(plt::Plot{GastonBackend}, series::Series)
    sp = series[:subplot]; gsp = sp.o
    x, y, z = series[:x], series[:y], series[:z]
    st = series[:seriestype]

    curves = []
    if gsp.dims == 2 && z === nothing
        for (n, seg) ∈ enumerate(series_segments(series, st; check=true))
            i, rng = seg.attr_index, seg.range
            for sc ∈ gaston_seriesconf!(sp, series, i, n == 1)
                push!(curves, Gaston.Curve(x[rng], y[rng], nothing, nothing, sc))
            end
        end
    else
        if z isa Surface
            z = z.surf
            if st == :image
                z = reverse(Float32.(Gray.(z)), dims=1)  # flip y axis
                nr, nc = size(z)
                if (ly = length(y)) == 2 && ly != nr
                    y = collect(range(y[1], y[2], length=nr))
                end
                if (lx = length(x)) == 2 && lx != nc
                    x = collect(range(x[1], x[2], length=nc))
                end
            elseif st == :heatmap
                length(x) == size(z, 2) + 1 && (x = @view x[1:end-1])
                length(y) == size(z, 1) + 1 && (y = @view y[1:end-1])
            end
        end
        if st == :mesh3d
            x, y, z = mesh3d_triangles(x, y, z, series[:connections])
        end
        for sc ∈ gaston_seriesconf!(sp, series, 1, true)
            push!(curves, Gaston.Curve(x, y, z, nothing, sc))
        end
    end

    for c ∈ curves
        append = length(gsp.curves) > 0; push!(gsp.curves, c)
        Gaston.write_data(c, gsp.dims, gsp.datafile, append=append)
    end
    nothing
end

function gaston_seriesconf!(sp::Subplot{GastonBackend}, series::Series, i::Int, add_to_legend::Bool)
    #=
    gnuplot abbreviations (see gnuplot/src/set.c)
    ---------------------------------------------
    dl: dashlength
    dt: dashtype
    fc: fillcolor
    fs: fillstyle
    lc: linecolor
    lp: linespoints
    ls: linestyle
    lt: linetype
    lw: linewidth
    pi: pointinterval
    pn: pointnumber
    ps: pointscale
    pt: pointtype
    tc: textcolor
    w: with
    =#
    gsp = sp.o; st = series[:seriestype]; extra = []
    add_to_legend &= should_add_to_legend(series)
    curveconf = String[add_to_legend ? "title '$(series[:label])'" : "notitle"]

    clims = get_clims(sp, series)
    if st ∈ (:scatter, :scatter3d)
        lc, dt, lw = gaston_lc_ls_lw(series, clims, i)
        pt, ps, mc = gaston_mk_ms_mc(series, clims, i)
        push!(curveconf, "w points pt $pt ps $ps lc $mc")
    elseif st ∈ (:path, :straightline, :path3d)
        lc, dt, lw = gaston_lc_ls_lw(series, clims, i)
        if series[:markershape] == :none  # simplepath
            push!(curveconf, "w lines lc $lc dt $dt lw $lw")
        else
            pt, ps, mc = gaston_mk_ms_mc(series, clims, i)
            push!(curveconf, "w lp lc $mc dt $dt lw $lw pt $pt ps $ps")
        end
    elseif st == :shape
        fc = gaston_color(get_fillcolor(series, i), get_fillalpha(series, i))
        lc, _ = gaston_lc_ls_lw(series, clims, i)
        push!(curveconf, "w filledcurves fc $fc fs solid border lc $lc")
    elseif st ∈ (:steppre, :stepmid, :steppost)
        step = if st == :steppre
            "fsteps"
        elseif st == :stepmid
            "histeps"
        elseif st == :steppost
            "steps"
        end
        push!(curveconf, "w $step")
        lc, dt, lw = gaston_lc_ls_lw(series, clims, i)
        push!(extra, "w points lc $lc dt $dt lw $lw notitle")
    elseif st == :image
        palette = gaston_palette(series[:seriescolor])
        gsp.axesconf *= "\nset palette model RGB defined $palette"
        push!(curveconf, "w image pixels")
    elseif st ∈ (:contour, :contour3d)
        push!(curveconf, "w lines")
        st == :contour && (gsp.axesconf *= "\nset view map\nunset surface")  # 2D
        levels = join(map(string, collect(contour_levels(series, clims))), ", ")
        gsp.axesconf *= "\nset contour base\nset cntrparam levels discrete $levels"
    elseif st ∈ (:surface, :heatmap)
        push!(curveconf, "w pm3d")
        palette = gaston_palette(series[:seriescolor])
        gsp.axesconf *= "\nset palette model RGB defined $palette"
        st == :heatmap && (gsp.axesconf *= "\nset view map")
    elseif st ∈ (:wireframe, :mesh3d)
        lc, dt, lw = gaston_lc_ls_lw(series, clims, i)
        push!(curveconf, "w lines lc $lc dt $dt lw $lw")
    elseif st == :quiver
        push!(curveconf, "w vectors filled")
    else
        @warn "Gaston: $st is not implemented yet"
    end

    return [join(curveconf, " "), extra...]
end

function gaston_parse_axes_args(
    plt::Plot{GastonBackend}, sp::Subplot{GastonBackend}, dims::Int, any_label::Bool
)
    # axesconf = String["set margins 2, 2, 2, 2"]  # left, right, bottom, top
    axesconf = String[]

    polar = ispolar(sp) && dims == 2  # cannot splot in polar coordinates

    for letter ∈ (:x, :y, :z)
        (letter == :z && dims == 2) && continue
        axis = sp.attr[Symbol(letter, :axis)]
        # label names
        push!(axesconf, "set $(letter)label '$(axis[:guide])' $(gaston_font(guidefont(axis)))")
        mirror = axis[:mirror] ? "mirror" : "nomirror"

        if axis[:scale] == :identity
            logscale, base = "nologscale", ""
        elseif axis[:scale] == :log10
            logscale, base = "logscale", "10"
        elseif axis[:scale] == :log2
            logscale, base = "logscale", "2"
        elseif axis[:scale] == :ln
            logscale, base = "logscale", "e"
        end
        push!(axesconf, "set $logscale $letter $base")

        # handle ticks
        if polar
            push!(axesconf, "set size square\nunset $(letter)tics")
        else
            push!(axesconf, "set $(letter)tics $(mirror) $(axis[:tick_direction]) $(gaston_font(tickfont(axis)))")

            # major tick locations
            if axis[:ticks] != :native
                if axis[:flip]
                    hi, lo = axis_limits(sp, letter)
                else
                    lo, hi = axis_limits(sp, letter)
                end
                push!(axesconf, "set $(letter)range [$lo:$hi]")

                ticks = get_ticks(sp, axis)
                gaston_set_ticks!(axesconf, ticks, letter, "", "")

                if axis[:minorticks] != :native
                    minor_ticks = get_minor_ticks(sp, axis, ticks)
                    gaston_set_ticks!(axesconf, minor_ticks, letter, "m", "add")
                end
            end
        end

        if axis[:grid]
            push!(axesconf, "set grid " * (polar ? "polar" : "$(letter)tics"))
            axis[:minorgrid] && push!(axesconf, "set grid " * (polar ? "polar" : "m$(letter)tics"))
        end

        ratio = get_aspect_ratio(sp)
        if ratio != :none
            ratio == :equal && (ratio = -1)
            push!(axesconf, "set size ratio $ratio")
        end
    end
    gaston_set_legend!(axesconf, sp, any_label)

    if hascolorbar(sp)
        push!(axesconf, "set cbtics $(gaston_font(colorbartitlefont(sp)))")
    end

    if sp[:title] !== nothing
        push!(axesconf, "set title '$(sp[:title])' $(gaston_font(titlefont(sp)))")
    end

    if polar
        push!(axesconf, "unset border\nset polar\nset border polar")
        tmin, tmax = axis_limits(sp, :x, false, false)
        rmin, rmax = axis_limits(sp, :y, false, false)
        rticks = get_ticks(sp, :y)
        if (ttype = ticksType(rticks)) == :ticks
            gaston_ticks = string.(rticks)
        elseif ttype == :ticks_and_labels
            gaston_ticks = String["'$l' $t" for (t, l) ∈ zip(rticks...)]
        end
        push!(axesconf, "set rtics ( " * join(gaston_ticks, ", ") * " ) $(gaston_font(tickfont(sp.attr[:yaxis])))")
        push!(axesconf, "set trange [$(min(0, tmin)):$(max(2π, tmax))]")
        push!(axesconf, "set rrange [$rmin:$rmax]")
        push!(axesconf, "set ttics 0,30 format \"%g\".GPVAL_DEGREE_SIGN $(gaston_font(tickfont(sp.attr[:xaxis])))")
        push!(axesconf, "set mttics 3")
    end

    return join(axesconf, "\n")
end

function gaston_set_ticks!(axesconf, ticks, letter, maj_min, add)
    ticks == :auto && return
    if ticks ∈ (:none, nothing, false)
        push!(axesconf, "unset $(maj_min)$(letter)tics")
        return
    end

    gaston_ticks = String[]
    if (ttype = ticksType(ticks)) == :ticks
        tick_locs = @view ticks[:]
        for i ∈ eachindex(tick_locs)
            tick = if maj_min == "m"
                "'' $(tick_locs[i]) 1"  # see gnuplot manual 'Mxtics'
            else
                "$(tick_locs[i])"
            end
            push!(gaston_ticks, tick)
        end
    elseif ttype == :ticks_and_labels
        tick_locs = @view ticks[1][:]
        tick_labels = @view ticks[2][:]
        for i ∈ eachindex(tick_locs)
            lab = gaston_enclose_tick_string(tick_labels[i])
            push!(gaston_ticks, "'$lab' $(tick_locs[i])")
        end
    else
        gaston_ticks = nothing
        @error "Gaston: invalid input for $(maj_min)$(letter)ticks: $ticks"
    end
    if gaston_ticks !== nothing
        push!(axesconf, "set $(letter)tics $add (" * join(gaston_ticks, ", ") * ")")
    end
    nothing
end

function gaston_set_legend!(axesconf, sp, any_label)
    leg = sp[:legend]
    if sp[:legend] ∉ (:none, :inline) && any_label
        leg == :best && (leg = :topright)

        push!(axesconf, "set key " * (occursin("outer", string(leg)) ? "outside" : "inside"))
        for position ∈ ("top", "bottom", "left", "right")
            occursin(position, string(leg)) && push!(axesconf, "set key $position")
        end
        push!(axesconf, "set key $(gaston_font(legendfont(sp), rot=false, align=false))")
        if sp[:legendtitle] !== nothing
            # NOTE: cannot use legendtitlefont(sp) as it will override legendfont
            push!(axesconf, "set key title '$(sp[:legendtitle])'")
        end
        push!(axesconf, "set key box lw 1 opaque")
        push!(axesconf, "set border back")
    else
        push!(axesconf, "set key off")
    end
    nothing
end

# --------------------------------------------
# Helpers
# --------------------------------------------

gaston_halign(k) = (left=:left, hcenter=:center, right=:right)[k]
gaston_valign(k) = (top=:top, vcenter=:center, bottom=:bottom)[k]

gaston_alpha(alpha) = alpha === nothing ? 0 : alpha

gaston_lc_ls_lw(series::Series, clims, i::Int) = (
    gaston_color(get_linecolor(series, clims, i), get_linealpha(series, i)),
    gaston_linestyle(get_linestyle(series, i)),
    get_linewidth(series, i),
)

gaston_mk_ms_mc(series::Series, clims, i::Int) = (
    gaston_marker(_cycle(series[:markershape], i), get_markeralpha(series, i)),
    _cycle(series[:markersize], i) * 1.3 / 5,
    gaston_color(get_markercolor(series, clims, i), get_markeralpha(series, i)),
)

function gaston_font(f; rot=true, align=true, color=true, scale=1)
    font = String["font '$(f.family),$(round(Int, scale * f.pointsize))'"]
    align && push!(font, "$(gaston_halign(f.halign))")
    rot && push!(font, "rotate by $(f.rotation)")
    color && push!(font, "textcolor $(gaston_color(f.color))")
    return join(font, " ")
end

function gaston_palette(gradient)
    palette = String[]; n = -1
    for rgba ∈ gradient  # FIXME: naive conversion, inefficient ?
        push!(palette, "$(n += 1) $(rgba.r) $(rgba.g) $(rgba.b)")
    end
    return '(' * join(palette, ", ") * ')'
end

function gaston_marker(marker, alpha)
    # NOTE: :rtriangle, :ltriangle, :hexagon, :heptagon, :octagon seems unsupported by gnuplot
    filled = gaston_alpha(alpha) == 0
    marker == :none && return -1
    marker == :pixel && return 0
    marker ∈ (:+, :cross) && return 1
    marker ∈ (:x, :xcross) && return 2
    marker == :star5 && return 3
    marker == :rect && return filled ? 5 : 4
    marker == :circle && return filled ? 7 : 6
    marker == :utriangle && return filled ? 9 : 8
    marker == :dtriangle && return filled ? 11 : 10
    marker == :diamond && return filled ? 13 : 12
    marker == :pentagon && return filled ? 15 : 14

    @warn "Gaston: unsupported marker $marker"
    return 1
end

function gaston_color(col, alpha=0)
    col = single_color(col)  # in case of gradients
    col = alphacolor(col, gaston_alpha(alpha))  # add a default alpha if non existent
    return "rgb '#$(hex(col, :aarrggbb))'"
end

function gaston_linestyle(style)
    style == :solid && return "1"
    style == :dash && return "2"
    style == :dot && return "3"
    style == :dashdot && return "4"
    style == :dashdotdot && return "5"
end

function gaston_enclose_tick_string(tick_string)
    findfirst("^", tick_string) === nothing && return tick_string
    base, power = split(tick_string, "^")
    return "$base^{$power}"
end


text:
module MyApp

using Example
using HelloWorldC_jll
using Artifacts
using Distributed
using Random

const myrand = rand()

const outputo = begin
    o = Base.JLOptions().outputo
    o == C_NULL ? "ok" : unsafe_string(o)
end

fooifier_path() = joinpath(artifact"fooifier", "bin", "fooifier" * (Sys.iswindows() ? ".exe" : ""))

function julia_main()
    try
        real_main()
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function is_crayons_loaded()
    Base.PkgId(Base.UUID("a8cc5b0e-0ffa-5ad4-8c14-923d3ee1735f"), "Crayons") in keys(Base.loaded_modules)
end

function real_main()
    @show ARGS
    @show Base.PROGRAM_FILE
    @show DEPOT_PATH
    @show LOAD_PATH
    @show pwd()
    @show Base.active_project()
    @show Sys.BINDIR

    @show Base.JLOptions().opt_level
    @show Base.JLOptions().nthreads
    @show Base.JLOptions().check_bounds

    display(Base.loaded_modules)
    println()

    println("Running a jll package:")
    HelloWorldC_jll.hello_world() do x
        println("HelloWorld artifact at $(realpath(x))")
        run(`$x`)
    end
    println()

    @show is_crayons_loaded()

    println("Running the artifact")
    res = read(`$(fooifier_path()) 5 10`, String)
    println("The result of 2*5^2 - 10 == $res")

    @show unsafe_string(Base.JLOptions().image_file)
    @show Example.domath(5)
    @show sin(0.0)

    println("outputo: $outputo")

    # Check that the RNG is seeded during precompilation
    println("myrand: ", myrand == 0.0 ? "fail" : "ok")
    rand() # Check that RNG state is ok
    if nworkers() != 4
        addprocs(4)
        @eval @everywhere using MyApp
    end

    n = @distributed (+) for i = 1:20000000
        1
    end
    println("n = $n")
    @eval @everywhere using Example
    @everywhere println(Example.domath(3))
    return
end

if abspath(PROGRAM_FILE) == @__FILE__
    real_main()
end

end # module

ext: jl


text:
"""
armijosc(xt, xc, ft, fc, d, residm, ItRules, derivative_is_old)
Line search for Newton's method. Read the notebook or print book
for the explanation. This is an internal function and I did not
design it to be hackable by the novice.
"""
function armijosc(xt, xc, ft, fc, d, residm, ItRules, derivative_is_old)
    idid = true
    alpha = 1.e-4
    iarm = -1
    lambda = 1.0
    lam0 = 0.0
    lamc = lambda
    lamm = lamc
    f = ItRules.f
    #    fp=ItRules.fp
    #    dx=ItRules.dx
    ResidC = residm
    armmax = ItRules.armmax
    armfix = ItRules.armfix
    #
    if derivative_is_old
        armmax = 0
    end
    #
    #   Take the full step and, if happy, go home.
    #
    (xt, ft, residt) = UpdateIteration(xt, xc, ft, lambda, d, ItRules)
    armfail = residt > (1 - alpha * lambda) * residm
    iarm += 1
    #
    #
    # At this point I've taken a full step. I'll enter the loop only if
    # that full step has failed.
    #
    ffc = residt^2
    ff0 = residm^2
    ffm = ffc
    while armfail && iarm < armmax
        #
        #   At this point the full step has failed. Now it's time to be
        #   serious about the line search.
        #
        lambda = update_lambda(iarm, armfix, lambda, lamc, ff0, ffc, ffm)
        (xt, ft, residt) = UpdateIteration(xt, xc, ft, lambda, d, ItRules)
        ffm = ffc
        ffc = residt^2
        iarm += 1
        armfail = residt > (1 - alpha * lambda) * residm
    end
    if iarm >= armmax
        idid = false
    end
    return (ax = xt, afc = ft, resnorm = residt, aiarm = iarm, idid = idid)
end

function update_lambda(iarm, armfix, lambda, lamc, ff0, ffc, ffm)
    if iarm == 0 || armfix == true
        lambda = lambda * 0.5
    else
        lamm = lamc
        lamc = lambda
        lambda = parab3p(lamc, lamm, ff0, ffc, ffm)
    end
    return lambda
end

"""
parab3p(lambdac, lambdam, ff0, ffc, ffm)

Three point parabolic line search.

input:\n
       lambdac = current steplength
       lambdam = previous steplength
       ff0 = value of || F(x_c) ||^2
       ffc = value of || F(x_c + lambdac d) ||^2
       ffm = value of || F(x_c + lambdam d) ||^2

output:\n
       lambdap = new value of lambda

internal parameters:\n
       sigma0 = .1, sigma1=.5, safeguarding bounds for the linesearch

You get here if cutting the steplength in half doesn't get you
sufficient decrease. Now you have three points and can build a parabolic
model. I do not like cubic models because they either need four points
or a derivative.

So let's think about how this works. I cheat a bit and check the model
for negative curvature, which I don't want to see.

 The polynomial is

 p(lambda) = ff0 + (c1 lambda + c2 lambda^2)/d1

 d1 = (lambdac - lambdam)*lambdac*lambdam < 0
 So if c2 > 0 we have negative curvature and default to
      lambdap = sigma0 * lambda
 The logic is that negative curvature is telling us that
 the polynomial model is not helping much, so it looks better
 to take the smallest possible step. This is not what I did in the
 matlab code because I did it wrong. I have sinced fixed it.

 So (Students, listen up!) if c2 < 0 then all we gotta do is minimize
 (c1 lambda + c2 lambda^2)/d1 over [.1* lambdac, .5*lambdac]
 This means to MAXIMIZE c1 lambda + c2 lambda^2 becase d1 < 0.
 So I find the zero of the derivative and check the endpoints.

"""
function parab3p(lambdac, lambdam, ff0, ffc, ffm)
    #
    # internal parameters
    #
    sigma0 = 0.1
    sigma1 = 0.5
    #
    c2 = lambdam * (ffc - ff0) - lambdac * (ffm - ff0)
    if c2 >= 0
        #
        # Sanity check for negative curvature
        #
        lambdap = sigma0 * lambdac
    else
        #
        # It's a convex parabola, so use calculus!
        #
        c1 = lambdac * lambdac * (ffm - ff0) - lambdam * lambdam * (ffc - ff0)
        lambdap = -c1 * 0.5 / c2
        #
        lambdaup = sigma1 * lambdac
        lambdadown = sigma0 * lambdac
        lambdap = max(lambdadown, min(lambdaup, lambdap))
    end
end


text:
# ---
# title: 572. Subtree of Another Tree
# id: problem572
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree
# link: <https://leetcode.com/problems/subtree-of-another-tree/description/>
# hidden: true
# ---
#
# Given two **non-empty** binary trees **s** and **t** , check whether tree
# **t** has exactly the same structure and node values with a subtree of **s**.
# A subtree of **s** is a tree consists of a node in **s** and all of this
# node's descendants. The tree **s** could also be considered as a subtree of
# itself.
#
# **Example 1:**
# Given tree s:
#
#
#
#          3
#         / \
#        4   5
#       / \
#      1   2
#
#
# Given tree t:
#
#
#
#        4
#       / \
#      1   2
#
#
# Return **true** , because t has the same structure and node values with a
# subtree of s.
#
#
#
# **Example 2:**
# Given tree s:
#
#
#
#          3
#         / \
#        4   5
#       / \
#      1   2
#         /
#        0
#
#
# Given tree t:
#
#
#
#        4
#       / \
#      1   2
#
#
# Return **false**.
#
#
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
export pmult, imult
export psparse, ParallelSparseMatrix
export SparseMatrixCSR, sparse_csr

type ParallelSparseMatrix{TF}
  F::TF
  refs::Vector{Any}
  procs::Vector{Int}
end

function psparse(F, procs)
  ParallelSparseMatrix(
    F,
    Any[@spawnat(i, fetch(F)) for i in procs],
    procs)
end

import Base.At_mul_B
import Base.isempty
import Base.size
import Base.eltype
import Base.Ac_mul_B

At_mul_B(A::ParallelSparseMatrix, U::AbstractMatrix) = pmult(size(A.F,2), A.refs, U, A.procs, imult)
*(A::ParallelSparseMatrix, U::AbstractMatrix)        = pmult(size(A.F,1), A.refs, U, A.procs, mult)
isempty(A::ParallelSparseMatrix) = isempty(A.F)
size(A::ParallelSparseMatrix)    = size(A.F)
size(A::ParallelSparseMatrix, i::Int) = size(A.F, i::Int)
eltype(A::ParallelSparseMatrix)       = eltype(A.F)
Ac_mul_B(A::ParallelSparseMatrix, B::ParallelSparseMatrix) = Ac_mul_B(A.F, B.F)
At_mul_B(A::ParallelSparseMatrix, B::ParallelSparseMatrix) = At_mul_B(A.F, B.F)

###### CSR matrix ######

type SparseMatrixCSR{Tv,Ti}
  csc::SparseMatrixCSC{Tv,Ti}
end

sparse_csr(csc::SparseMatrixCSC) = SparseMatrixCSR(csc')
sparse_csr(rows, cols, vals) = SparseMatrixCSR(sparse(cols, rows, vals))

At_mul_B(A::SparseMatrixCSR, u::AbstractVector) = A.csc * u
Ac_mul_B(A::SparseMatrixCSR, u::AbstractVector) = A.csc * u
*(A::SparseMatrixCSR, u::AbstractVector) = At_mul_B(A.csc, u)
At_mul_B(A::SparseMatrixCSR, B::SparseMatrixCSR) = A_mul_Bt(A.csc, B.csc)
isempty(A::SparseMatrixCSR) = isempty(A.csc)

eltype(A::SparseMatrixCSR) = eltype(A.csc)
function size(A::SparseMatrixCSR)
  m,n = size(A.csc)
  return n,m
end
size(A::SparseMatrixCSR,d) = (d>2 ? 1 : size(A)[d])

###### parallel multiplication ######

function imult(Fref, u)
  return At_mul_B(fetch(Fref), u)
end

function mult(Fref, u)
  #Flocal = fetch(Fref)
  #if size(Flocal,2) != size(u,1)
  #  error(@sprintf("#columns of F(%d) has to equal number of rows of U(%d).", size(Flocal,2), size(u,1)) )
  #end
  return fetch(Fref) * u
end

## setup:
## feat = genes.F
## Frefs = map(i -> @spawnat(i, fetch(feat)), workers())
function pmult(nrows::Integer, Frefs, U, procs, mfun)
    np = length(procs)  # determine the number of processes available
    n  = size(U,2)
    results = zeros(nrows, n)
    i = 1
    # function to produce the next work item from the queue.
    # in this case it's just an index.
    nextidx() = (idx=i; i+=1; idx)
    @sync begin
        for p in 1:length(procs)
            if procs[p] != myid() || np == 1
                @async begin
                    while true
                        idx = nextidx()
                        if idx > n
                            break
                        end
                        results[:,idx] = remotecall_fetch(mfun, procs[p], Frefs[p], U[:,idx])
                    end
                end
            end
        end
    end
    results
end


text:
## names.jl retrieve and set dimension names
## (c) 2013 David A. van Leeuwen

## This code is licensed under the MIT license
## See the file LICENSE.md in this distribution

import Base.names

## `names` is somewhat loaded, as some of the semantics were transferred tp `fieldnames`
names(dict::AbstractDict) = collect(keys(dict))
names(n::NamedArray) = [names(dict) for dict in n.dicts]
names(n::NamedArray, d::Integer) = names(n.dicts[d])
defaultnames(a::AbstractArray) = [defaultnames(a, d) for d in 1:ndims(a)]
defaultnames(a::AbstractArray, d::Integer) = defaultnames(size(a, d))

@deprecate allnames(n::NamedArray) names(n::NamedArray)
@deprecate allnames(a::AbstractArray) defaultnames(a::AbstractArray)

## dimnames gives array, for tuple use n.dimnames
dimnames(n::NamedArray) = [n.dimnames...]
dimnames(n::NamedArray, d::Integer) = n.dimnames[d]
dimnames(a::AbstractArray) = [defaultdimnames(a)...]
dimnames(a::AbstractArray, d::Integer) = defaultdimname(d)

## string versions of the above
strnames(dict::AbstractDict) = [isa(name, String) ? name : sprint(show, name, context=:compact => true) for name in names(dict)]
strnames(n::NamedArray) = [strnames(d) for d in n.dicts]
strnames(n::NamedArray, d::Integer) = strnames(n.dicts[d])
strdimnames(n::NamedArray) = [string(dn) for dn in n.dimnames]
strdimnames(n::NamedArray, d::Integer) = string(n.dimnames[d])

## seting names, dimnames
function setnames!(n::NamedArray{T,N}, v::Vector{KT}, d::Integer) where {T,N,KT}
    size(n.array, d) == length(v) || throw(DimensionMismatch("inconsistent vector length"))
    keytype(n.dicts[d]) == KT || throw(TypeError(:setnames!, "second argument", keytype(n.dicts[d]), KT))
    ## n.dicts is a tuple, so we need to replace it as a whole...
    vdicts = OrderedDict{Any,Int}[]
    for i = 1:length(n.dicts)
        if i==d
            push!(vdicts, OrderedDict(zip(v, 1:length(v))))
        else
            push!(vdicts, n.dicts[i])
        end
    end
    n.dicts = tuple(vdicts...)::NTuple{N}
end

function setnames!(n::NamedArray, v, d::Integer, i::Integer)
    1 <= d <= ndims(n) || throw(BoundsError("dimension"))
    1 <= i <= size(n, d) || throw(BoundsError("index"))
    isa(v, keytype(n.dicts[d])) || throw(TypeError(:setnames!, "second argument", keytype(n.dicts[d]), typeof(v)))
    filter!(pair -> pair[2]!=i, n.dicts[d]) # remove old name
    n.dicts[d][v] = i
end

function setdimnames!(n::NamedArray{T,N}, dn::NTuple{N,Any}) where {T,N}
    n.dimnames = dn
end

setdimnames!(n::NamedArray, dn::Vector) = setdimnames!(n, tuple(dn...))

function setdimnames!(n::NamedArray{T,N}, v, d::Integer) where {T,N}
    1 <= d <= N || throw(BoundsError(size(n), d))
    vdimnames = Array{Any}(undef, N)
    for i=1:N
        if i==d
            vdimnames[i] = v
        else
            vdimnames[i] = n.dimnames[i]
        end
    end
    n.dimnames = tuple(vdimnames...)
end


text:
@testset "maximumlikelihood_gev.jl" begin
    n = 100

    x₁ = Variable("x₁", randn(n))
    x₂ = Variable("x₂", randn(n) / 3)
    x₃ = Variable("x₃", randn(n) / 10)

    μ = 1.0 .+ x₁.value
    ϕ = -.05 .+ x₂.value
    ξ = x₃.value

    σ = exp.(ϕ)
    θ = [1.0; 1.0; -0.05; 1.0; 0.0; 1.0]

    pd = GeneralizedExtremeValue.(μ, σ, ξ)
    y = rand.(pd)

    @testset "gevfit(y; locationcov, logscalecov, shapecov)" begin
        # model building with non-stationary location, logscale and shape
        fm = Extremes.gevfit(y,
            locationcov = [x₁],
            logscalecov = [x₂],
            shapecov = [x₃])

        # data is y
        @test fm.model.data.value ≈ y

        # location is x₁
        @test length(fm.model.location.covariate) == 1
        @test fm.model.location.covariate[1].value ≈ x₁.value

        # logscale is x₂
        @test length(fm.model.logscale.covariate) == 1
        @test fm.model.logscale.covariate[1].value ≈ x₂.value

        # shape is x₃
        @test length(fm.model.shape.covariate) == 1
        @test fm.model.shape.covariate[1].value ≈ x₃.value

    end

    @testset "gevfit(y, initialvalues; locationcov, logscalecov, shapecov)" begin
        # Using initialvalues does not throw nor warn
        @test_logs fm = Extremes.gevfit(y, zeros(Float64, 6),
            locationcov = [x₁],
            logscalecov = [x₂],
            shapecov = [x₃])

    end

    @testset "gevfit(df, datacol; locationcovid, logscalecovid, shapecovid)" begin
        # model building with non-stationary location, logscale and shape
        df = DataFrame(y = y, x1 = x₁.value, x2 = x₂.value, x3 = x₃.value)

        fm = Extremes.gevfit(df, :y, locationcovid = [:x1], logscalecovid = [:x2], shapecovid = [:x3])

        # data is y
        @test fm.model.data.value ≈ y

        # location is x₁
        @test length(fm.model.location.covariate) == 1
        @test fm.model.location.covariate[1].value ≈ x₁.value

        # logscale is x₂
        @test length(fm.model.logscale.covariate) == 1
        @test fm.model.logscale.covariate[1].value ≈ x₂.value

        # shape is x₃
        @test length(fm.model.shape.covariate) == 1
        @test fm.model.shape.covariate[1].value ≈ x₃.value

    end

    @testset "gevfit(df, datacol, initialvalues; locationcovid, logscalecovid, shapecovid)" begin
        # Using initialvalues does not throw nor warn
        df = DataFrame(y = y, x1 = x₁.value, x2 = x₂.value, x3 = x₃.value)

        @test_logs Extremes.gevfit(df, :y, zeros(Float64, 6), locationcovid = [:x1], logscalecovid = [:x2], shapecovid = [:x3])

    end

    @testset "gevfit(model, initialvalues)" begin
        # non-stationary location, logscale and shape
        model = Extremes.BlockMaxima(Variable("y", y),
            locationcov = [x₁],
            logscalecov = [x₂],
            shapecov = [x₃])

        fm = Extremes.gevfit(model, zeros(Float64, 6))

        # data is y
        @test fm.model.data.value ≈ y

        # location is x₁
        @test length(fm.model.location.covariate) == 1
        @test fm.model.location.covariate[1].value ≈ x₁.value

        # logscale is x₂
        @test length(fm.model.logscale.covariate) == 1
        @test fm.model.logscale.covariate[1].value ≈ x₂.value

        # shape is x₃
        @test length(fm.model.shape.covariate) == 1
        @test fm.model.shape.covariate[1].value ≈ x₃.value

    end

end


text:
module RollCallDataIO
    include("roll_call_encoding.jl")
    include("roll_call.jl")
    include("zero_one_roll_call.jl")
    include("sparse_roll_call.jl")
    include("ord_file.jl")
end

max_stars_repo_name: johnmyleswhite/RollCallDataIO

max_issues_repo_name: johnmyleswhite/RollCallDataIO

max_forks_repo_name: johnmyleswhite/RollCallDataIO


text:
#---------------------------------------------------------------------
#*	only liquid phase
#*--------------------------------------------------------------------
type tank_liq
	tank_liq()=begin
		new(
			tank_basic(),
			[
				:(_base_1.Outletm.v = 0),
				:(_base_1.Outletm.h = _base_1.PP.LiquidEnthalpy(_base_1.Outletm.T,_base_1.Outletm.P,_base_1.Outletm.z)),
				:(_base_1.Tank.V = _base_1.Mt*_base_1.PP.LiquidVolume(_base_1.Outletm.T,_base_1.Outletm.P,_base_1.Outletm.z)),
				:(E = _base_1.Mt*_base_1.Outletm.h - _base_1.Outletm.P*_base_1.Tank.V),
			],
			[
				"Vapourisation fraction","Liquid Enthalpy","Volume constraint","Total internal energy",
			],
		)
	end
	_base_1::tank_basic
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	attributes::Dict{Symbol,Any}
end
export tank_liq
function setEquationFlow(in::tank_liq)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
end
function atributes(in::tank_liq,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Brief]="Model of a generic liquid-phase tank"
	drive!(fields,_)
	return fields
end
tank_liq(_::Dict{Symbol,Any})=begin
	newModel=tank_liq()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
using Plots
using Distributions

μfast = 0
σfast = 1
dfast = LogNormal(μfast, σfast)
x = rand(dfast, 1000)
histogram(x, dpi=600, xlabel="Page Load Time (s)", ylabel="Number of Loads", title="Fast User Experience", legend = false)
png("images/goodpageload.png")

dslow = LogNormal(3, 0.4)
x = rand(dslow, 1000)
histogram(x, dpi=600, xlabel="Page Load Time (s)", ylabel="Number of Loads", title="Slow User Experience", legend = false)
png("images/badpageload.png")

users = zeros(100)
usercount = zeros(100)

for i in 1:100
    while true
        if i < 90
            users[i] += rand(dfast)
        else
            users[i] += rand(dslow)
        end
        usercount[i] += 1.0
        if users[i] > 3600
            break
        end
    end
end

pageloads = rand(dfast, Int64(sum(usercount[1:90])))
append!(pageloads, rand(dslow, Int64(sum(usercount[91:100]))))
pageloads = log.(pageloads)
histogram(pageloads, label = "Simulated Page Loads", xlabel="Logarithm of Page Load Times", ylabel="Number of Page Loads", dpi=600)
vline!([μfast], label="Metric 1", linewidth=2)
vline!([log(mean(exp.(pageloads)))], label="Metric 2", linewidth=2)
vline!([log(median(users./usercount))], label="Metric 3", linewidth=2)
vline!([log(mean(users./usercount))], label="Metric 4", linewidth=2)
png("images/examplesimulation.png")

for k in 75:99
    users = zeros(100)
    usercount = zeros(100)

    for i in 1:100
        while true
            if i < k
                users[i] += rand(dfast)
            else
                users[i] += rand(dslow)
            end
            usercount[i] += 1.0
            if users[i] > 3600
                break
            end
        end
    end

    pageloads = rand(dfast, Int64(sum(usercount[1:k])))
    append!(pageloads, rand(dslow, Int64(sum(usercount[k+1:100]))))
    pageloads = log.(pageloads)
    if k==75
        scatter([k], [exp(μfast)], label="Metric 1", linewidth=2, color="red", dpi=600, xlabel="Percentage of Users with fast experience", ylabel="Metric Calculated Page Load Time (s)")
        scatter!([k], [mean(exp.(pageloads))], label="Metric 2", linewidth=2, color="black")
        scatter!([k], [median(users./usercount)], label="Metric 3", linewidth=2, color="green")
        scatter!([k], [mean(users./usercount)], label="Metric 4", linewidth=2, color="blue")
    else
        scatter!([k], [exp(μfast)], linewidth=2, color="red", label=false)
        scatter!([k], [mean(exp.(pageloads))], linewidth=2, color="black", label=false)
        scatter!([k], [median(users./usercount)], linewidth=2, color="green", label=false)
        scatter!([k], [mean(users./usercount)], linewidth=2, color="blue", label=false)
    end
end
png("images/varyingusers.png")


text:
struct Octonion{T<:Real} <: Number
  s::T
  v1::T
  v2::T
  v3::T
  v4::T
  v5::T
  v6::T
  v7::T
  norm::Bool
end

Octonion{T}(x::Real) where {T<:Real} = Octonion(convert(T, x))
Octonion{T}(x::Complex) where {T<:Real} = Octonion(convert(Complex{T}, x))
Octonion{T}(q::Quaternion) where {T<:Real} = Octonion(convert(Quaternion{T}, q))
Octonion{T}(o::Octonion) where {T<:Real} =
  Octonion{T}(o.s, o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7, o.norm)

Octonion(s::Real, v1::Real, v2::Real, v3::Real, v4::Real, v5::Real, v6::Real, v7::Real, n::Bool = false) =
  Octonion(promote(s, v1, v2, v3, v4, v5, v6, v7)..., n)
Octonion(x::Real) = Octonion(x, zero(x), zero(x), zero(x), zero(x), zero(x), zero(x), zero(x), abs(x) == one(x))
Octonion(z::Complex) = Octonion(z.re, z.im, zero(z.re), zero(z.re), zero(z.re), zero(z.re), zero(z.re), zero(z.re), abs(z) == one(z.re))
Octonion(q::Quaternion) = Octonion(q.s, q.v1, q.v2, q.v3, zero(q.s), zero(q.s), zero(q.s), zero(q.s), q.norm)
Octonion(s::Real, a::Vector) = Octonion(s, a[1], a[2], a[3], a[4], a[5], a[6], a[7])
Octonion(a::Vector) = Octonion(0, a[1], a[2], a[3], a[4], a[5], a[6], a[7])

const OctonionF16 = Octonion{Float16}
const OctonionF32 = Octonion{Float32}
const OctonionF64 = Octonion{Float64}

promote_rule(::Type{Octonion{T}}, ::Type{S}) where {T <: Real, S <: Real} = Octonion{promote_type(T, S)}
promote_rule(::Type{Octonion{T}}, ::Type{Complex{S}}) where {T <: Real, S <: Real} = Octonion{promote_type(T, S)}
promote_rule(::Type{Octonion{T}}, ::Type{Quaternion{S}}) where {T <: Real, S <: Real} = Octonion{promote_type(T, S)}
promote_rule(::Type{Octonion{T}}, ::Type{Octonion{S}}) where {T <: Real, S <: Real} = Octonion{promote_type(T, S)}

octo(p, v1, v2, v3, v4, v5, v6, v7) = Octonion(p, v1, v2, v3, v4, v5, v6, v7)
octo(p, v1, v2, v3, v4, v5, v6, v7, n) = Octonion(p, v1, v2, v3, v4, v5, v6, v7, n)
octo(x) = Octonion(x)
octo(s, a) = Octonion(s, a)

function show(io::IO, o::Octonion)
  pm(x) = x < 0 ? " - $(-x)" : " + $x"
  print(io, o.s, pm(o.v1), "im", pm(o.v2), "jm", pm(o.v3), "km", pm(o.v4), "ilm", pm(o.v5), "jlm", pm(o.v6), "klm", pm(o.v7), "lm")
end

real(o::Octonion) = o.s
imag_part(o::Octonion) = (o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7)
@deprecate imag(o::Octonion) collect(imag_part(o)) false

(/)(o::Octonion, x::Real) = Octonion(o.s / x, o.v1 / x, o.v2 / x, o.v3 / x, o.v4 / x, o.v5 / x, o.v6 / x, o.v7 / x)
(*)(o::Octonion, x::Real) = Octonion(o.s * x, o.v1 * x, o.v2 * x, o.v3 * x, o.v4 * x, o.v5 * x, o.v6 * x, o.v7 * x)
(*)(x::Real, o::Octonion) = o * x

conj(o::Octonion) = Octonion(o.s, -o.v1, -o.v2, -o.v3, -o.v4, -o.v5, -o.v6, -o.v7, o.norm)
abs(o::Octonion) = sqrt(abs2(o))
float(q::Octonion{T}) where T = convert(Octonion{float(T)}, q)
abs_imag(o::Octonion) = sqrt((o.v4 * o.v4 + (o.v2 * o.v2 + o.v6 * o.v6)) + ((o.v1 * o.v1 + o.v5 * o.v5) + (o.v3 * o.v3 + o.v7 * o.v7))) # ordered to match abs2
abs2(o::Octonion) = ((o.s * o.s + o.v4 * o.v4) + (o.v2 * o.v2 + o.v6 * o.v6)) + ((o.v1 * o.v1  + o.v5 * o.v5) + (o.v3 * o.v3 + o.v7 * o.v7))
inv(o::Octonion) = o.norm ? conj(o) : conj(o) / abs2(o)

isreal(o::Octonion) = iszero(o.v1) & iszero(o.v2) & iszero(o.v3) & iszero(o.v4) & iszero(o.v5) & iszero(o.v6) & iszero(o.v7)
isfinite(o::Octonion) = o.norm | (isfinite(real(o)) & isfinite(o.v1) & isfinite(o.v2) & isfinite(o.v3) & isfinite(o.v4) & isfinite(o.v5) & isfinite(o.v6) & isfinite(o.v7))
iszero(o::Octonion) = ~o.norm & iszero(real(o)) & iszero(o.v1) & iszero(o.v2) & iszero(o.v3) & iszero(o.v4) & iszero(o.v5) & iszero(o.v6) & iszero(o.v7)
isnan(o::Octonion) = isnan(real(o)) | isnan(o.v1) | isnan(o.v2) | isnan(o.v3) | isnan(o.v4) | isnan(o.v5) | isnan(o.v6) | isnan(o.v7)
isinf(o::Octonion) = ~o.norm & (isinf(real(o)) | isinf(o.v1) | isinf(o.v2) | isinf(o.v3) | isinf(o.v4) | isinf(o.v5) | isinf(o.v6) | isinf(o.v7))

function normalize(o::Octonion)
  if (o.norm)
    return o
  end
  o = o / abs(o)
  Octonion(o.s, o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7, true)
end

function normalizea(o::Octonion)
  if (o.norm)
    return (o, one(o.s))
  end
  a = abs(o)
  o = o / a
  (Octonion(o.s, o.v1, o.v2, o.v3, o.v4, o.v5, o.v6, o.v7, true), a)
end

(-)(o::Octonion) = Octonion(-o.s, -o.v1, -o.v2, -o.v3, -o.v4, -o.v5, -o.v6, -o.v7, o.norm)

(+)(o::Octonion, w::Octonion) = Octonion(o.s + w.s,
                                            o.v1 + w.v1,
                                            o.v2 + w.v2,
                                            o.v3 + w.v3,
                                            o.v4 + w.v4,
                                            o.v5 + w.v5,
                                            o.v6 + w.v6,
                                            o.v7 + w.v7)

(-)(o::Octonion, w::Octonion) = Octonion(o.s - w.s,
                                            o.v1 - w.v1,
                                            o.v2 - w.v2,
                                            o.v3 - w.v3,
                                            o.v4 - w.v4,
                                            o.v5 - w.v5,
                                            o.v6 - w.v6,
                                            o.v7 - w.v7)

function (*)(o::Octonion, w::Octonion)
    s  = ((o.s * w.s - o.v4 * w.v4) - (o.v2 * w.v2 + o.v6 * w.v6)) - ((o.v1 * w.v1 + o.v5 * w.v5) + (o.v3 * w.v3 + o.v7 * w.v7))
    v1 = ((o.s * w.v1 + o.v1 * w.s) + (o.v6 * w.v5 - o.v5 * w.v6)) + ((o.v2 * w.v3 - o.v3 * w.v2) + (o.v7 * w.v4 - o.v4 * w.v7))
    v2 = ((o.s * w.v2 + o.v2 * w.s) + (o.v4 * w.v6 - o.v6 * w.v4)) + ((o.v3 * w.v1 - o.v1 * w.v3) + (o.v7 * w.v5 - o.v5 * w.v7))
    v3 = ((o.s * w.v3 + o.v3 * w.s) + (o.v5 * w.v4 - o.v4 * w.v5)) + ((o.v1 * w.v2 - o.v2 * w.v1) + (o.v7 * w.v6 - o.v6 * w.v7))
    v4 = ((o.s * w.v4 + o.v4 * w.s) + (o.v3 * w.v5 - o.v5 * w.v3)) + ((o.v1 * w.v7 - o.v7 * w.v1) + (o.v6 * w.v2 - o.v2 * w.v6))
    v5 = ((o.s * w.v5 + o.v5 * w.s) + (o.v2 * w.v7 - o.v7 * w.v2)) + ((o.v1 * w.v6 - o.v6 * w.v1) + (o.v4 * w.v3 - o.v3 * w.v4))
    v6 = ((o.s * w.v6 + o.v6 * w.s) + (o.v3 * w.v7 - o.v7 * w.v3)) + ((o.v2 * w.v4 - o.v4 * w.v2) + (o.v5 * w.v1 - o.v1 * w.v5))
    v7 = ((o.s * w.v7 + o.v7 * w.s) + (o.v5 * w.v2 - o.v2 * w.v5)) + ((o.v4 * w.v1 - o.v1 * w.v4) + (o.v6 * w.v3 - o.v3 * w.v6))
    return Octonion(s, v1, v2, v3, v4, v5, v6, v7, o.norm & w.norm)
end

(/)(o::Octonion, w::Octonion) = o * inv(w)

(==)(q::Octonion, w::Octonion) = (q.s == w.s) & (q.v1 == w.v1) & (q.v2 == w.v2) & (q.v3 == w.v3) &
                                 (q.v4 == w.v4) & (q.v5 == w.v5) & (q.v6 == w.v6) & (q.v7 == w.v7) # ignore .norm field

function exp(o::Octonion)
  s = o.s
  se = exp(s)
  scale = se
  th = abs_imag(o)
  if th > 0
    scale *= sin(th) / th
  end
  Octonion(se * cos(th),
            scale * o.v1,
            scale * o.v2,
            scale * o.v3,
            scale * o.v4,
            scale * o.v5,
            scale * o.v6,
            scale * o.v7,
            abs(s) == 0)
end

function log(o::Octonion)
  o, a = normalizea(o)
  s = o.s
  M = abs_imag(o)
  th = atan(M, s)
  if M > 0
    M = th / M
    return Octonion(log(a),
                     o.v1 * M,
                     o.v2 * M,
                     o.v3 * M,
                     o.v4 * M,
                     o.v5 * M,
                     o.v6 * M,
                     o.v7 * M)
  else
    return Octonion(complex(log(a), ifelse(iszero(a), zero(th), th)))
  end
end

(^)(o::Octonion, w::Octonion) = exp(w * log(o))

function sqrt(o::Octonion)
  exp(0.5 * log(o))
end

octorand() = octo(randn(), randn(), randn(), randn(), randn(), randn(), randn(), randn())

function rand(rng::AbstractRNG, ::Random.SamplerType{Octonion{T}}) where {T<:Real}
  Octonion{T}(rand(rng, T), rand(rng, T), rand(rng, T), rand(rng, T),
              rand(rng, T), rand(rng, T), rand(rng, T), rand(rng, T), false)
end

function randn(rng::AbstractRNG, ::Type{Octonion{T}}) where {T<:AbstractFloat}
  Octonion{T}(
      randn(rng, T) * INV_SQRT_EIGHT,
      randn(rng, T) * INV_SQRT_EIGHT,
      randn(rng, T) * INV_SQRT_EIGHT,
      randn(rng, T) * INV_SQRT_EIGHT,
      randn(rng, T) * INV_SQRT_EIGHT,
      randn(rng, T) * INV_SQRT_EIGHT,
      randn(rng, T) * INV_SQRT_EIGHT,
      randn(rng, T) * INV_SQRT_EIGHT,
      false,
  )
end


text:
"""
    EigenAngle

Plane-wave propagation direction in angle `θ` from the vertical in radians.

Common trigonometric function values of this angle are calculated to increase performance.

# Fields

- `θ::ComplexF64`
- `cosθ::ComplexF64`
- `sinθ::ComplexF64`
- `secθ::ComplexF64`
- `cos²θ::ComplexF64`
- `sin²θ::ComplexF64`

`Real` `θ` will be automatically converted to `Complex`.

!!! note

    Technically this is an angle of incidence from the vertical, and not necessarily an
    _eigen_angle unless it is found to be associated with a propagated mode in the waveguide.
"""
struct EigenAngle
    θ::ComplexF64  # radians, because df/dθ are in radians
    cosθ::ComplexF64
    sinθ::ComplexF64
    secθ::ComplexF64  # == 1/cosθ
    cos²θ::ComplexF64
    sin²θ::ComplexF64
end

EigenAngle(ea::EigenAngle) = ea
EigenAngle(θ::Real) = EigenAngle(complex(θ))
function EigenAngle(θ::Complex)
    rθ, iθ = reim(θ)
    ((abs(rθ) > 2π) || (abs(iθ) > 2π)) && @warn "θ > 2π. Make sure θ is in radians."
    S, C = sincos(θ)
    Cinv = inv(C)  # == sec(θ)
    C² = C^2
    S² = 1 - C²
    EigenAngle(θ, C, S, Cinv, C², S²)
end

"""
    isless(x::EigenAngle, y::EigenAngle)

Calls `Base.isless` for complex `EigenAngle` by `real` and then `imag` component.

By defining `isless`, calling `sort` on `EigenAngle`s sorts by real component first and
imaginary component second.
"""
function Base.isless(x::EigenAngle, y::EigenAngle)
    return isless((real(x.θ), imag(x.θ)), (real(y.θ), imag(y.θ)))
end

function Base.show(io::IO, e::EigenAngle)
    compact = get(io, :compact, false)

    if compact
        print(io, e.θ)
    else
        print(io, e.θ)
    end
    return nothing
end

function Base.show(io::IO, ::MIME"text/plain", e::EigenAngle)
    println(io, "EigenAngle: ")
    show(io, e)
end

"""
    isdetached(ea::EigenAngle, frequency::Frequency; params=LMPParams())

Return `true` if this is likely an earth detached (whispering gallery) mode according to the
criteria in [Pappert1981] eq. 1 with the additional criteria that the frequency be
above 100 kHz.

# References

[Pappert1981]: R. A. Pappert, “LF daytime earth ionosphere waveguide calculations,” Naval
    Ocean Systems Center, San Diego, CA, NOSC/TR-647, Jan. 1981.
    [Online]. Available: https://apps.dtic.mil/docs/citations/ADA096098.
"""
function isdetached(ea::EigenAngle, frequency::Frequency; params=LMPParams())
    C, C² = ea.cosθ, ea.cos²θ
    f, k = frequency.f, frequency.k
    @unpack earthradius, curvatureheight = params
    α = 2/earthradius

    return f > 100e3 && real(2im/3*(k/α)*(C²- α*curvatureheight)^(3/2)) > 12.4
end

"""
    referencetoground(ea::EigenAngle; params=LMPParams())

Reference `ea` from `params.curvatureheight` to ground (``z = 0``).
"""
function referencetoground(ea::EigenAngle; params=LMPParams())
    # see, e.g. PS71 pg 11
    @unpack earthradius, curvatureheight = params
    return EigenAngle(asin(ea.sinθ/sqrt(1 - 2/earthradius*curvatureheight)))
end

"""
    referencetoground(sinθ; params=LMPParams())

Reference ``sin(θ)`` from `params.curvatureheight` to the ground and return ``sin(θ₀)`` at
the ground.
"""
function referencetoground(sinθ; params=LMPParams())
    @unpack earthradius, curvatureheight = params
    return sinθ/sqrt(1 - 2/earthradius*curvatureheight)
end

"""
    attenuation(ea, frequency::Frequency; params=LMPParams())

Compute attenuation of eigenangle `ea` at the ground.

`ea` will be converted to an `EigenAngle` if needed.

This function internally references `ea` to the ground.
"""
function attenuation(ea, frequency::Frequency; params=LMPParams())
    ea = EigenAngle(ea)
    S₀ = referencetoground(ea.sinθ)
    neper2dB = 20log10(exp(1))  # 1 Np ≈ 8.685 dB
    return -neper2dB*frequency.k*imag(S₀)*1e6
end

"""
    phasevelocity(ea)

Compute the relative phase velocity ``v/c`` associated with the eigenangle `θ`.

`ea` will be converted to an `EigenAngle` if needed.

This function internally references `ea` to the ground.
"""
function phasevelocity(ea)
    ea = EigenAngle(ea)
    S₀ = referencetoground(ea.sinθ)
    return 1/real(S₀)
end


text:
module test_lin
using ApproXD
using Test
# using SparseArrays
# using LinearAlgebra

@testset "constructor for 1D Lininterp on a single function" begin

	vs = rand(4)
	gs = Array{Float64,1}[]
	push!(gs, range(2.0, stop = 3, length = 4))

	l = Lininterp(vs,gs)

	@test ApproXD.getDims(l) == [4]
	@test l.n == 1
	@test l.nfunc == 1
	@test l.ifunc == [1]
	@test l.infs == zeros(1)
	@test l.sups == zeros(1)
	@test l.hits == 0
	@test l.miss == 0
	@test l.vals[1] == vs
	@test l.grids == gs
	@test ApproXD.getCache(l) == [1]

	# errors
	pop!(gs)   # missing grid
	@test_throws ArgumentError Lininterp(vs,gs)
	push!(gs, range(-1,stop =3.0, length=5)) 	# wrong length
	@test_throws ArgumentError Lininterp(vs,gs)
	pop!(gs)
	push!(gs, [1.1,0,2,3]) 	# not sorted
	@test_throws ArgumentError Lininterp(vs,gs)

end

@testset "constructor for 2D Lininterp on a single function" begin

	vs = rand(3,4)
	gs = Array{Float64,1}[]
	push!(gs, range(1.0, stop = 3, length = 3))
	push!(gs, range(2.0, stop = 3, length = 4))

	l = Lininterp(vs,gs)

	@test ApproXD.getDims(l) == [3,4]
	@test l.n == 2
	@test l.nfunc == 1
	@test l.ifunc == [1]
	@test l.infs == zeros(2)
	@test l.sups == zeros(2)
	@test l.hits == 0
	@test l.miss == 0
	@test l.vals[1] == vs
	@test l.grids == gs
	@test ApproXD.getCache(l) == [1,1]

	# errors
	pop!(gs)   # missing grid
	@test_throws ArgumentError Lininterp(vs,gs)
	push!(gs, range(-1,stop =3.0, length=5)) 	# wrong length
	@test_throws ArgumentError Lininterp(vs,gs)
	pop!(gs)
	push!(gs, [1.1,0,2,3]) 	# not sorted
	@test_throws ArgumentError Lininterp(vs,gs)

end

@testset "constructor for 2D Lininterp on two functions" begin

	vs1 = rand(3,4)
	vs2 = rand(3,4)
	gs = Array{Float64,1}[]
	push!(gs, range(1.0, stop = 3, length = 3))
	push!(gs, range(2.0, stop = 3, length = 4))

	l = Lininterp(vs1,vs2,gs)

	@test ApproXD.getDims(l) == [3,4]
	@test l.n == 2
	@test l.nfunc == 2
	@test l.ifunc == [1,2]
	@test l.infs == zeros(2)
	@test l.sups == zeros(2)
	@test l.hits == 0
	@test l.miss == 0
	@test l.vals[1] == vs1
	@test l.vals[2] == vs2
	@test l.grids == gs
	@test ApproXD.getCache(l) == [1,1]

	# errors
	pop!(gs)   # missing grid
	@test_throws UndefVarError Lininterp(vs,gs)
	push!(gs, range(-1,stop =3.0, length=5)) 	# wrong length
	@test_throws UndefVarError Lininterp(vs,gs)
	pop!(gs)
	push!(gs, [1.1,0,2,3]) 	# not sorted
	@test_throws UndefVarError Lininterp(vs,gs)

end

@testset "constructor for 3D Lininterp" begin

	vs = rand(3,4,5)
	gs = Array{Float64,1}[]
	push!(gs, range(1.0, stop = 3, length = 3))
	push!(gs, range(2.0, stop = 3, length = 4))
	push!(gs, range(-1,stop =3.0, length=5))

	l = Lininterp(vs,gs)

	@test ApproXD.getDims(l) == [3,4,5]
	@test l.n == 3
	@test l.ifunc == [1]
	@test l.infs == zeros(3)
	@test l.sups == zeros(3)
	@test l.hits == 0
	@test l.miss == 0
	@test l.vals[1] == vs
	@test l.grids == gs
	@test ApproXD.getCache(l) == [1,1,1]

	# errors
	pop!(gs)   # missing grid
	@test_throws ArgumentError Lininterp(vs,gs)
	push!(gs, range(-1, stop = 3.0, length = 4)) 	# wrong length
	@test_throws ArgumentError Lininterp(vs,gs)
	pop!(gs)
	push!(gs, [1.1,0,2,3,4]) 	# not sorted
	@test_throws ArgumentError Lininterp(vs,gs)

end

@testset "constructor for 4D Lininterp" begin

	vs = rand(3,4,5,7)
	gs = Array{Float64,1}[]
	push!(gs, range(1.0, stop = 3, length = 3))
	push!(gs, range(2.0, stop = 3, length = 4))
	push!(gs, range(-1,stop =3.0, length=5))
	push!(gs, range(-0.1,stop = 4.8,length = 7))

	l = Lininterp(vs,gs)

	@test ApproXD.getDims(l) == [3,4,5,7]
	@test l.n == 4
	@test l.infs == zeros(4)
	@test l.sups == zeros(4)
	@test l.hits == 0
	@test l.miss == 0
	@test l.vals[1] == vs
	@test l.grids == gs
	@test ApproXD.getCache(l) == [1,1,1,1]

	# errors
	pop!(gs)   # missing grid
	@test_throws ArgumentError Lininterp(vs,gs)
	push!(gs, range(-1, stop = 3.0, length = 4)) 	# wrong length
	@test_throws ArgumentError Lininterp(vs,gs)
	pop!(gs)
	push!(gs, [1.1,0,2,3,4,5,6]) 	# not sorted
	@test_throws ArgumentError Lininterp(vs,gs)

end

@testset "testing Cache functions" begin

	vs = rand(3,4,5)
	gs = Array{Float64,1}[]
	push!(gs, range(1.0, stop = 3, length = 3))
	push!(gs, range(2.0, stop = 3, length = 4))
	push!(gs, range(-1,stop =3.0, length=5))

	l = Lininterp(vs,gs)
	myc = [2,1,3]
	l.cache = myc	# manually set a cache

	for i in 1:length(gs)
		@test ApproXD.getCache(l,i) == myc[i]
		@test ApproXD.getCachedVal(l,i) == gs[i][myc[i]]
		@test ApproXD.getNextCachedVal(l,i) == gs[i][myc[i]+1]
	end

end

@testset "testing findBracket!(x) for 1D" begin

	vs = rand(5)
	gs = Array{Float64,1}[]
	push!(gs, range(-1,stop =3.0, length=5))

	l = Lininterp(vs,gs)
	@test l.hascache == false

	x = [2.8] 	# == brackets 4
	ApproXD.findBracket!(l,x)

	# brackets is index of lower bound of each bracket
	brackets = [4]

	@test ApproXD.getCache(l) == brackets
	@test ApproXD.hitmiss(l) == (0,1)
	@test l.infs[1] == gs[1][brackets[1]]
	@test l.sups[1] == gs[1][brackets[1]+1]

	@test l.hascache == true
	ApproXD.findBracket!(l,x)
	@test ApproXD.getCache(l) == brackets
	@test ApproXD.hitmiss(l) == (1,1)

	# value x1 and x2 out of bounds
	x = [3.8] 	# == brackets 4
	ApproXD.findBracket!(l,x)

	brackets = [4]
	@test ApproXD.getCache(l) == brackets
	@test ApproXD.hitmiss(l) == (1,2)
	@test l.infs[1] == gs[1][end-1]
	@test l.sups[1] == gs[1][end]

	# side-effect: x is forced into bounds!
	# may or may not be desirable.
	@test x[1] == l.sups[1]
end

@testset "testing findBracket!(x) for 3D" begin

	vs = rand(3,4,5)
	gs = Array{Float64,1}[]
	push!(gs, range(1.0, stop = 3, length = 3))
	push!(gs, range(2.0,stop = 6.0,length = 4))
	push!(gs, range(-1,stop =3.0, length=5))

	l = Lininterp(vs,gs)
	@test l.hascache == false

	x = [1.1,3.5,2.8] 	# == brackets 1,2,4
	ApproXD.findBracket!(l,x)

	# brackets is index of lower bound of each bracket
	brackets = [1,2,4]

	@test ApproXD.getCache(l) == brackets
	@test ApproXD.hitmiss(l) == (0,3)
	for i in 1:3
		@test l.infs[i] == gs[i][brackets[i]]
		@test l.sups[i] == gs[i][brackets[i]+1]
	end

	@test l.hascache == true
	ApproXD.findBracket!(l,x)
	@test ApproXD.getCache(l) == brackets
	@test ApproXD.hitmiss(l) == (3,3)

	# value x1 and x2 out of bounds
	x = [0.9,7.5,2.8] 	# == brackets 1,3,4
	ApproXD.findBracket!(l,x)

	brackets = [1,3,4]
	@test ApproXD.getCache(l) == brackets
	@test ApproXD.hitmiss(l) == (4,5)
	@test l.infs[1] == gs[1][1]
	@test l.infs[2] == gs[2][end-1]
	@test l.infs[3] == gs[3][brackets[3]]
	@test l.sups[1] == gs[1][2]
	@test l.sups[2] == gs[2][end]
	@test l.sups[3] == gs[3][brackets[3]+1]

	# side-effect: x is forced into bounds!
	# may or may not be desirable.
	@test x[1] == l.infs[1]
	@test x[2] == l.sups[2]
end

@testset "testing findBracket!(x) for 4D" begin

	vs = rand(3,4,5,7)
	gs = Array{Float64,1}[]
	push!(gs, range(1.0, stop = 3, length = 3))
	push!(gs, range(2.0,stop = 6.0,length = 4))
	push!(gs, range(-1,stop =3.0, length=5))
	push!(gs, range(-2,stop = 5.0,length = 7))

	l = Lininterp(vs,gs)
	@test l.hascache == false

	x = [1.1,3.5,2.8,3.9] 	# == brackets 1,2,4,6
	ApproXD.findBracket!(l,x)

	# brackets is index of lower bound of each bracket
	brackets = [1,2,4,6]

	@test ApproXD.getCache(l) == brackets
	@test ApproXD.hitmiss(l) == (0,4)
	for i in 1:4
		@test l.infs[i] == gs[i][brackets[i]]
		@test l.sups[i] == gs[i][brackets[i]+1]
	end

	@test l.hascache == true
	ApproXD.findBracket!(l,x)
	@test ApproXD.getCache(l) == brackets
	@test ApproXD.hitmiss(l) == (4,4)

	# value x1 and x2 out of bounds
	x = [0.9,7.5,2.8,3.9] 	# == brackets 1,3,4,6
	ApproXD.findBracket!(l,x)

	brackets = [1,3,4,6]
	@test ApproXD.getCache(l) == brackets
	@test ApproXD.hitmiss(l) == (6,6)
	@test l.infs[1] == gs[1][1]
	@test l.infs[2] == gs[2][end-1]
	@test l.infs[3] == gs[3][brackets[3]]
	@test l.infs[4] == gs[4][brackets[4]]
	@test l.sups[1] == gs[1][2]
	@test l.sups[2] == gs[2][end]
	@test l.sups[3] == gs[3][brackets[3]+1]
	@test l.sups[4] == gs[4][brackets[4]+1]

	# side-effect: x is forced into bounds!
	# may or may not be desirable.
	@test x[1] == l.infs[1]
	@test x[2] == l.sups[2]
end

@testset "testing find3DVertices!(l) for 1 function" begin

	vs = rand(3,4,5)
	gs = Array{Float64,1}[]
	push!(gs, range(1.0, stop = 3, length = 3))
	push!(gs, range(2.0,stop = 6.0,length = 4))
	push!(gs, range(-1,stop =3.0, length=5))

	l = Lininterp(vs,gs)
	x = [1.1,3.5,2.8] 	# == brackets 1,2,4
	ApproXD.findBracket!(l,x)
	@test ApproXD.hitmiss(l) == (0,3)
	ApproXD.find3DVertices!(l)

	@test length(l.vertex[1]) == 8
	@test l.vertex[1][1] == vs[1,2,4]  # v000
	@test l.vertex[1][2] == vs[2,2,4]  # v100
	@test l.vertex[1][3] == vs[1,3,4]  # v010
	@test l.vertex[1][4] == vs[1,2,5]  # v001
	@test l.vertex[1][5] == vs[2,3,4]  # v110
	@test l.vertex[1][6] == vs[2,2,5]  # v101
	@test l.vertex[1][7] == vs[1,3,5]  # v011
	@test l.vertex[1][8] == vs[2,3,5]  # v111

	x = [3,6,3.0] 	# == brackets 2,3,4
	ApproXD.findBracket!(l,x)
	@test ApproXD.hitmiss(l) == (0,6)
	ApproXD.find3DVertices!(l)

	@test l.vertex[1][1] == vs[2,3,4]  # v000
	@test l.vertex[1][2] == vs[3,3,4]  # v100
	@test l.vertex[1][3] == vs[2,4,4]  # v010
	@test l.vertex[1][4] == vs[2,3,5]  # v001
	@test l.vertex[1][5] == vs[3,4,4]  # v110
	@test l.vertex[1][6] == vs[3,3,5]  # v101
	@test l.vertex[1][7] == vs[2,4,5]  # v011
	@test l.vertex[1][8] == vs[3,4,5]  # v111

	# out of bounds: no hits by convention
	x = [-1,6.5,-3.0] 	# == brackets 1,3,1
	ApproXD.findBracket!(l,x)
	@test ApproXD.hitmiss(l) == (0,9)
	ApproXD.find3DVertices!(l)

	@test l.vertex[1][1] == vs[1,3,1]  # v000
	@test l.vertex[1][2] == vs[2,3,1]  # v100
	@test l.vertex[1][3] == vs[1,4,1]  # v010
	@test l.vertex[1][4] == vs[1,3,2]  # v001
	@test l.vertex[1][5] == vs[2,4,1]  # v110
	@test l.vertex[1][6] == vs[2,3,2]  # v101
	@test l.vertex[1][7] == vs[1,4,2]  # v011
	@test l.vertex[1][8] == vs[2,4,2]  # v111
end

@testset "testing find3DVertices!(l) for 2 functions" begin

	vs1 = rand(3,4,5)
	vs2 = rand(3,4,5)
	gs = Array{Float64,1}[]
	push!(gs, range(1.0, stop = 3, length = 3))
	push!(gs, range(2.0,stop = 6.0,length = 4))
	push!(gs, range(-1,stop =3.0, length=5))

	l = Lininterp(vs1,vs2,gs)
	x = [1.1,3.5,2.8] 	# == brackets 1,2,4
	ApproXD.findBracket!(l,x)
	@test ApproXD.hitmiss(l) == (0,3)
	ApproXD.find3DVertices!(l)

	@test length(l.vertex[1]) == 8
	@test l.vertex[1][1] == vs1[1,2,4]  # v000
	@test l.vertex[1][2] == vs1[2,2,4]  # v100
	@test l.vertex[1][3] == vs1[1,3,4]  # v010
	@test l.vertex[1][4] == vs1[1,2,5]  # v001
	@test l.vertex[1][5] == vs1[2,3,4]  # v110
	@test l.vertex[1][6] == vs1[2,2,5]  # v101
	@test l.vertex[1][7] == vs1[1,3,5]  # v011
	@test l.vertex[1][8] == vs1[2,3,5]  # v111

	@test l.vertex[2][1] == vs2[1,2,4]  # v000
	@test l.vertex[2][2] == vs2[2,2,4]  # v100
	@test l.vertex[2][3] == vs2[1,3,4]  # v010
	@test l.vertex[2][4] == vs2[1,2,5]  # v001
	@test l.vertex[2][5] == vs2[2,3,4]  # v110
	@test l.vertex[2][6] == vs2[2,2,5]  # v101
	@test l.vertex[2][7] == vs2[1,3,5]  # v011
	@test l.vertex[2][8] == vs2[2,3,5]  # v111

	x = [3,6,3.0] 	# == brackets 2,3,4
	ApproXD.findBracket!(l,x)
	@test ApproXD.hitmiss(l) == (0,6)
	ApproXD.find3DVertices!(l)

	@test l.vertex[1][1] == vs1[2,3,4]  # v000
	@test l.vertex[1][2] == vs1[3,3,4]  # v100
	@test l.vertex[1][3] == vs1[2,4,4]  # v010
	@test l.vertex[1][4] == vs1[2,3,5]  # v001
	@test l.vertex[1][5] == vs1[3,4,4]  # v110
	@test l.vertex[1][6] == vs1[3,3,5]  # v101
	@test l.vertex[1][7] == vs1[2,4,5]  # v011
	@test l.vertex[1][8] == vs1[3,4,5]  # v111

	@test l.vertex[2][1] == vs2[2,3,4]  # v000
	@test l.vertex[2][2] == vs2[3,3,4]  # v100
	@test l.vertex[2][3] == vs2[2,4,4]  # v010
	@test l.vertex[2][4] == vs2[2,3,5]  # v001
	@test l.vertex[2][5] == vs2[3,4,4]  # v110
	@test l.vertex[2][6] == vs2[3,3,5]  # v101
	@test l.vertex[2][7] == vs2[2,4,5]  # v011
	@test l.vertex[2][8] == vs2[3,4,5]  # v111

	# out of bounds: no hits by convention
	x = [-1,6.5,-3.0] 	# == brackets 1,3,1
	ApproXD.findBracket!(l,x)
	@test ApproXD.hitmiss(l) == (0,9)
	ApproXD.find3DVertices!(l)

	@test l.vertex[1][1] == vs1[1,3,1]  # v000
	@test l.vertex[1][2] == vs1[2,3,1]  # v100
	@test l.vertex[1][3] == vs1[1,4,1]  # v010
	@test l.vertex[1][4] == vs1[1,3,2]  # v001
	@test l.vertex[1][5] == vs1[2,4,1]  # v110
	@test l.vertex[1][6] == vs1[2,3,2]  # v101
	@test l.vertex[1][7] == vs1[1,4,2]  # v011
	@test l.vertex[1][8] == vs1[2,4,2]  # v111

	@test l.vertex[2][1] == vs2[1,3,1]  # v000
	@test l.vertex[2][2] == vs2[2,3,1]  # v100
	@test l.vertex[2][3] == vs2[1,4,1]  # v010
	@test l.vertex[2][4] == vs2[1,3,2]  # v001
	@test l.vertex[2][5] == vs2[2,4,1]  # v110
	@test l.vertex[2][6] == vs2[2,3,2]  # v101
	@test l.vertex[2][7] == vs2[1,4,2]  # v011
	@test l.vertex[2][8] == vs2[2,4,2]  # v111
end

@testset "testing find4DVertices!(l) single function" begin

	vs = rand(3,4,5,7)
	gs = Array{Float64,1}[]
	push!(gs, range(1.0, stop = 3, length = 3))
	push!(gs, range(2.0,stop = 6.0,length = 4))
	push!(gs, range(-1,stop =3.0, length=5))
	push!(gs, range(-2,stop = 5.0,length = 7))

	l = Lininterp(vs,gs)
	x = [1.1,3.5,2.8,3.9] 	# == brackets 1,2,4,6
	ApproXD.findBracket!(l,x)
	@test ApproXD.hitmiss(l) == (0,4)
	ApproXD.find4DVertices!(l)

	@test length(l.vertex[1]) == 2^l.n
	@test l.vertex[1][1]  == vs[1,2,4,6]  # v0000
	@test l.vertex[1][2]  == vs[2,2,4,6]  # v1000
	@test l.vertex[1][3]  == vs[1,3,4,6]  # v0100
	@test l.vertex[1][4]  == vs[1,2,5,6]  # v0010
	@test l.vertex[1][5]  == vs[1,2,4,7]  # v0001
	@test l.vertex[1][6]  == vs[2,3,4,6]  # v1100
	@test l.vertex[1][7]  == vs[2,2,5,6]  # v1010
	@test l.vertex[1][8]  == vs[2,2,4,7]  # v1001
	@test l.vertex[1][9]  == vs[1,3,5,6]  # v0110
	@test l.vertex[1][10] == vs[1,3,4,7]  # v0101
	@test l.vertex[1][11] == vs[1,2,5,7]  # v0011
	@test l.vertex[1][12] == vs[2,3,5,6]  # v1110
	@test l.vertex[1][13] == vs[2,3,4,7]  # v1101
	@test l.vertex[1][14] == vs[2,2,5,7]  # v1011
	@test l.vertex[1][15] == vs[1,3,5,7]  # v0111
	@test l.vertex[1][16] == vs[2,3,5,7]  # v1111

	x = [3,6,3.0,1.5] 	# == brackets 2,3,4,4
	ApproXD.findBracket!(l,x)
	@test ApproXD.hitmiss(l) == (0,8)
	ApproXD.find4DVertices!(l)

	@test l.vertex[1][1]  == vs[2,3,4,4]  # v0000
	@test l.vertex[1][2]  == vs[3,3,4,4]  # v1000
	@test l.vertex[1][3]  == vs[2,4,4,4]  # v0100
	@test l.vertex[1][4]  == vs[2,3,5,4]  # v0010
	@test l.vertex[1][5]  == vs[2,3,4,5]  # v0001
	@test l.vertex[1][6]  == vs[3,4,4,4]  # v1100
	@test l.vertex[1][7]  == vs[3,3,5,4]  # v1010
	@test l.vertex[1][8]  == vs[3,3,4,5]  # v1001
	@test l.vertex[1][9]  == vs[2,4,5,4]  # v0110
	@test l.vertex[1][10] == vs[2,4,4,5]  # v0101
	@test l.vertex[1][11] == vs[2,3,5,5]  # v0011
	@test l.vertex[1][12] == vs[3,4,5,4]  # v1110
	@test l.vertex[1][13] == vs[3,4,4,5]  # v1101
	@test l.vertex[1][14] == vs[3,3,5,5]  # v1011
	@test l.vertex[1][15] == vs[2,4,5,5]  # v0111
	@test l.vertex[1][16] == vs[3,4,5,5]  # v1111

	# out of bounds
	x = [-1,6.5,-3.0,1.5] 	# == brackets 1,3,1,4
	ApproXD.findBracket!(l,x)
	@test ApproXD.hitmiss(l) == (1,11)
	ApproXD.find4DVertices!(l)

	@test l.vertex[1][1]  == vs[1,3,1,4]  # v0000
	@test l.vertex[1][2]  == vs[2,3,1,4]  # v1000
	@test l.vertex[1][3]  == vs[1,4,1,4]  # v0100
	@test l.vertex[1][4]  == vs[1,3,2,4]  # v0010
	@test l.vertex[1][5]  == vs[1,3,1,5]  # v0001
	@test l.vertex[1][6]  == vs[2,4,1,4]  # v1100
	@test l.vertex[1][7]  == vs[2,3,2,4]  # v1010
	@test l.vertex[1][8]  == vs[2,3,1,5]  # v1001
	@test l.vertex[1][9]  == vs[1,4,2,4]  # v0110
	@test l.vertex[1][10] == vs[1,4,1,5]  # v0101
	@test l.vertex[1][11] == vs[1,3,2,5]  # v0011
	@test l.vertex[1][12] == vs[2,4,2,4]  # v1110
	@test l.vertex[1][13] == vs[2,4,1,5]  # v1101
	@test l.vertex[1][14] == vs[2,3,2,5]  # v1011
	@test l.vertex[1][15] == vs[1,4,2,5]  # v0111
	@test l.vertex[1][16] == vs[2,4,2,5]  # v1111
end

@testset "testing find4DVertices!(l) 2 functions" begin

	vs1= rand(3,4,5,7)
	vs2= rand(3,4,5,7)
	gs = Array{Float64,1}[]
	push!(gs, range(1.0, stop = 3, length = 3))
	push!(gs, range(2.0,stop = 6.0,length = 4))
	push!(gs, range(-1,stop =3.0, length=5))
	push!(gs, range(-2,stop = 5.0,length = 7))

	l = Lininterp(vs1,vs2,gs)
	x = [1.1,3.5,2.8,3.9] 	# == brackets 1,2,4,6
	ApproXD.findBracket!(l,x)
	@test ApproXD.hitmiss(l) == (0,4)
	ApproXD.find4DVertices!(l)

	@test length(l.vertex[1]) == 2^l.n
	@test l.vertex[1][1]  == vs1[1,2,4,6]  # v0000
	@test l.vertex[1][2]  == vs1[2,2,4,6]  # v1000
	@test l.vertex[1][3]  == vs1[1,3,4,6]  # v0100
	@test l.vertex[1][4]  == vs1[1,2,5,6]  # v0010
	@test l.vertex[1][5]  == vs1[1,2,4,7]  # v0001
	@test l.vertex[1][6]  == vs1[2,3,4,6]  # v1100
	@test l.vertex[1][7]  == vs1[2,2,5,6]  # v1010
	@test l.vertex[1][8]  == vs1[2,2,4,7]  # v1001
	@test l.vertex[1][9]  == vs1[1,3,5,6]  # v0110
	@test l.vertex[1][10] == vs1[1,3,4,7]  # v0101
	@test l.vertex[1][11] == vs1[1,2,5,7]  # v0011
	@test l.vertex[1][12] == vs1[2,3,5,6]  # v1110
	@test l.vertex[1][13] == vs1[2,3,4,7]  # v1101
	@test l.vertex[1][14] == vs1[2,2,5,7]  # v1011
	@test l.vertex[1][15] == vs1[1,3,5,7]  # v0111
	@test l.vertex[1][16] == vs1[2,3,5,7]  # v1111

	@test l.vertex[2][1]  == vs2[1,2,4,6]  # v0000
	@test l.vertex[2][2]  == vs2[2,2,4,6]  # v1000
	@test l.vertex[2][3]  == vs2[1,3,4,6]  # v0100
	@test l.vertex[2][4]  == vs2[1,2,5,6]  # v0010
	@test l.vertex[2][5]  == vs2[1,2,4,7]  # v0001
	@test l.vertex[2][6]  == vs2[2,3,4,6]  # v1100
	@test l.vertex[2][7]  == vs2[2,2,5,6]  # v1010
	@test l.vertex[2][8]  == vs2[2,2,4,7]  # v1001
	@test l.vertex[2][9]  == vs2[1,3,5,6]  # v0110
	@test l.vertex[2][10] == vs2[1,3,4,7]  # v0101
	@test l.vertex[2][11] == vs2[1,2,5,7]  # v0011
	@test l.vertex[2][12] == vs2[2,3,5,6]  # v1110
	@test l.vertex[2][13] == vs2[2,3,4,7]  # v1101
	@test l.vertex[2][14] == vs2[2,2,5,7]  # v1011
	@test l.vertex[2][15] == vs2[1,3,5,7]  # v0111
	@test l.vertex[2][16] == vs2[2,3,5,7]  # v1111

	x = [3,6,3.0,1.5] 	# == brackets 2,3,4,4
	ApproXD.findBracket!(l,x)
	@test ApproXD.hitmiss(l) == (0,8)
	ApproXD.find4DVertices!(l)

	@test l.vertex[1][1]  == vs1[2,3,4,4]  # v0000
	@test l.vertex[1][2]  == vs1[3,3,4,4]  # v1000
	@test l.vertex[1][3]  == vs1[2,4,4,4]  # v0100
	@test l.vertex[1][4]  == vs1[2,3,5,4]  # v0010
	@test l.vertex[1][5]  == vs1[2,3,4,5]  # v0001
	@test l.vertex[1][6]  == vs1[3,4,4,4]  # v1100
	@test l.vertex[1][7]  == vs1[3,3,5,4]  # v1010
	@test l.vertex[1][8]  == vs1[3,3,4,5]  # v1001
	@test l.vertex[1][9]  == vs1[2,4,5,4]  # v0110
	@test l.vertex[1][10] == vs1[2,4,4,5]  # v0101
	@test l.vertex[1][11] == vs1[2,3,5,5]  # v0011
	@test l.vertex[1][12] == vs1[3,4,5,4]  # v1110
	@test l.vertex[1][13] == vs1[3,4,4,5]  # v1101
	@test l.vertex[1][14] == vs1[3,3,5,5]  # v1011
	@test l.vertex[1][15] == vs1[2,4,5,5]  # v0111
	@test l.vertex[1][16] == vs1[3,4,5,5]  # v1111

	@test l.vertex[2][1]  == vs2[2,3,4,4]  # v0000
	@test l.vertex[2][2]  == vs2[3,3,4,4]  # v1000
	@test l.vertex[2][3]  == vs2[2,4,4,4]  # v0100
	@test l.vertex[2][4]  == vs2[2,3,5,4]  # v0010
	@test l.vertex[2][5]  == vs2[2,3,4,5]  # v0001
	@test l.vertex[2][6]  == vs2[3,4,4,4]  # v1100
	@test l.vertex[2][7]  == vs2[3,3,5,4]  # v1010
	@test l.vertex[2][8]  == vs2[3,3,4,5]  # v1001
	@test l.vertex[2][9]  == vs2[2,4,5,4]  # v0110
	@test l.vertex[2][10] == vs2[2,4,4,5]  # v0101
	@test l.vertex[2][11] == vs2[2,3,5,5]  # v0011
	@test l.vertex[2][12] == vs2[3,4,5,4]  # v1110
	@test l.vertex[2][13] == vs2[3,4,4,5]  # v1101
	@test l.vertex[2][14] == vs2[3,3,5,5]  # v1011
	@test l.vertex[2][15] == vs2[2,4,5,5]  # v0111
	@test l.vertex[2][16] == vs2[3,4,5,5]  # v1111

	# out of bounds
	x = [-1,6.5,-3.0,1.5] 	# == brackets 1,3,1,4
	ApproXD.findBracket!(l,x)
	@test ApproXD.hitmiss(l) == (1,11)
	ApproXD.find4DVertices!(l)

	@test l.vertex[1][1]  == vs1[1,3,1,4]  # v0000
	@test l.vertex[1][2]  == vs1[2,3,1,4]  # v1000
	@test l.vertex[1][3]  == vs1[1,4,1,4]  # v0100
	@test l.vertex[1][4]  == vs1[1,3,2,4]  # v0010
	@test l.vertex[1][5]  == vs1[1,3,1,5]  # v0001
	@test l.vertex[1][6]  == vs1[2,4,1,4]  # v1100
	@test l.vertex[1][7]  == vs1[2,3,2,4]  # v1010
	@test l.vertex[1][8]  == vs1[2,3,1,5]  # v1001
	@test l.vertex[1][9]  == vs1[1,4,2,4]  # v0110
	@test l.vertex[1][10] == vs1[1,4,1,5]  # v0101
	@test l.vertex[1][11] == vs1[1,3,2,5]  # v0011
	@test l.vertex[1][12] == vs1[2,4,2,4]  # v1110
	@test l.vertex[1][13] == vs1[2,4,1,5]  # v1101
	@test l.vertex[1][14] == vs1[2,3,2,5]  # v1011
	@test l.vertex[1][15] == vs1[1,4,2,5]  # v0111
	@test l.vertex[1][16] == vs1[2,4,2,5]  # v1111

	@test l.vertex[2][1]  == vs2[1,3,1,4]  # v0000
	@test l.vertex[2][2]  == vs2[2,3,1,4]  # v1000
	@test l.vertex[2][3]  == vs2[1,4,1,4]  # v0100
	@test l.vertex[2][4]  == vs2[1,3,2,4]  # v0010
	@test l.vertex[2][5]  == vs2[1,3,1,5]  # v0001
	@test l.vertex[2][6]  == vs2[2,4,1,4]  # v1100
	@test l.vertex[2][7]  == vs2[2,3,2,4]  # v1010
	@test l.vertex[2][8]  == vs2[2,3,1,5]  # v1001
	@test l.vertex[2][9]  == vs2[1,4,2,4]  # v0110
	@test l.vertex[2][10] == vs2[1,4,1,5]  # v0101
	@test l.vertex[2][11] == vs2[1,3,2,5]  # v0011
	@test l.vertex[2][12] == vs2[2,4,2,4]  # v1110
	@test l.vertex[2][13] == vs2[2,4,1,5]  # v1101
	@test l.vertex[2][14] == vs2[2,3,2,5]  # v1011
	@test l.vertex[2][15] == vs2[1,4,2,5]  # v0111
	@test l.vertex[2][16] == vs2[2,4,2,5]  # v1111
end

@testset "testing getValue 1D" begin

	lbs = [1.0]
	ubs = [3.0]

	gs = Array{Float64,1}[]
	push!(gs, range(lbs[1],stop = ubs[1],length = 4))

	myfun(i1) = 0.75*i1

	vs = Float64[ myfun(i) for i in gs[1]]

	l = Lininterp(vs,gs)

	# check value on bounds
	@test getValue(l,lbs)[1] == vs[1]
	@test getValue(l,ubs)[1] == vs[4]
	@test getValue(l,[5.0])[1] == vs[4]
	@test getValue(l,[1.0])[1]== vs[1]

	# check values out of bounds
	@test getValue(l,[-1.0])[1] == vs[1]
	@test getValue(l,[200.0])[1] == vs[4]
	println(l)

	# close to bounds
	y=2.9
	@test isapprox(getValue(l,[y])[1] - myfun(y) ,0.0,atol=1e-6)
	@test isapprox(getValue(l,[y])[1] - myfun(y) ,0.0,atol=1e-6)

	# check at random vals in interval
	lbs = [1.0]
	ubs = [33.0]
	gs = Array{Float64,1}[]
	push!(gs, range(lbs[1],stop = ubs[1], length = 40))

	vs = Float64[ myfun(i) for i in gs[1]]

	l = Lininterp(vs,gs)
	for i in 1:30
		x = rand() * (ubs[1]-lbs[1]) + lbs[1]
		@test isapprox(getValue(l,[x])[1] ,myfun(x),atol=1e-6)
	end
	println(l)
end

@testset "testing getValue! 1D" begin

	myfun(i1) = 0.75*i1

	lbs = [1.0]
	ubs = [33.0]
	gs = Array{Float64,1}[]
	push!(gs, range(lbs[1],stop = ubs[1], length = 40))

	vs = Float64[ myfun(i) for i in gs[1]]

	l = Lininterp(vs,gs)
	y = [0.0]
	# check at increasing vals in interval
	x = range(lbs[1],stop = ubs[1],length = 300)
	for i in 1:length(x)
		getValue!(y,l,[x[i]],[1])
		@test isapprox(y[1] ,myfun(x[i]),atol=1e-8)
	end
	println(l)
	resetCache!(l)
	# check at random vals in interval
	for i in 1:30
		x = rand() * (ubs[1]-lbs[1]) + lbs[1]
		getValue!(y,l,[x],[1])
		@test isapprox(y[1] ,myfun(x),atol=1e-8)
	end
	println(l)
end

@testset "testing getValue! 1D on irregularly spaced grid" begin

	lb = 1.0
	ub = 3.0
	n = 40

	xg = zeros(n)
	xg[1] = log(lb + 1)
	xg[n] = log(ub + 1)
	xg    = range(xg[1],stop = xg[n],length = n)
	xg    = exp.(xg) .- 1

	gs = Array{Float64,1}[]
	push!(gs, xg)

	# I'm using a linear function here: the more non-linear your function,
	# the worse a linear approximation will perform - this is obvious.
	myfun(i1) = 0.75*i1

	vs = Float64[ myfun(i) for i in gs[1]]

	l = Lininterp(vs,gs)
	y = [0.0]
	# check at increasing vals in interval
	x = range(lb,stop = ub,length = 300)
	for i in 1:length(x)
		getValue!(y,l,[x[i]],[1])
		@test isapprox(y[1] ,myfun(x[i]),atol=1e-8)
		# println(y[1] - myfun(x[i]))
	end
	println(l)
	resetCache!(l)
	# check at random vals in interval
	for i in 1:30
		x = rand() * (ub-lb) + lb
		getValue!(y,l,[x],[1])
		@test isapprox(y[1] ,myfun(x),atol=1e-8)
	end
	println(l)
end

@testset "testing getValue2D" begin

	lbs = [1.0,2.0]
	ubs = [3.0,5.0]

	gs = Array{Float64,1}[]
	push!(gs, range(lbs[1],stop = ubs[1],length = 3))
	push!(gs, range(lbs[2],stop = ubs[2],length = 4))

	myfun(i1,i2) = 0.5*i1 + 2*i2

	vs = Float64[ myfun(i,j) for i in gs[1], j in gs[2]]

	l = Lininterp(vs,gs)

	# check value on bounds
	@test getValue(l,lbs)[1] == vs[1,1]
	@test getValue(l,ubs)[1] == vs[3,4]
	@test getValue(l,[1.0,5.0])[1] == vs[1,4]
	@test getValue(l,[1.0,5.0])[1]== vs[1,4]

	# check values out of bounds
	@test getValue(l,[-1.0,2.0])[1] == vs[1,1]
	@test getValue(l,[1.0,200.0])[1] == vs[1,4]
	println(l)

	# close to bounds
	x=1.99
	y=4.9
	@test isapprox(getValue(l,[x,y])[1] - myfun(x,y) ,0.0,atol=1e-6)
	@test isapprox(getValue(l,[x,y])[1] - myfun(x,y) ,0.0,atol=1e-6)

	x=1.4861407584066377
	y=3.5646251730324234
	@test isapprox(getValue(l,[x,y])[1] - myfun(x,y) ,0.0,atol=1e-6)

	x = 2.53
	y = 2.58
	@test isapprox(getValue(l,[x,y])[1] - myfun(x,y) ,0.0,atol=1e-6)

	# check at random vals in interval
	lbs = [1.0,2.0]
	ubs = [33.0,5.0]
	gs = Array{Float64,1}[]
	push!(gs, range(lbs[1],stop = ubs[1],length = 30))
	push!(gs, range(lbs[2],stop = ubs[2],length = 40))

	vs = Float64[ myfun(i,j) for i in gs[1], j in gs[2]]

	l = Lininterp(vs,gs)
	for i in 1:30
		x = rand() * (ubs[1]-lbs[1]) + lbs[1]
		y = rand() * (ubs[2]-lbs[2]) + lbs[2]
		@test isapprox(getValue(l,[x,y])[1] ,myfun(x,y),atol=1e-6)
	end
	println(l)
end

@testset "testing getValue3D" begin

	lbs = [1.0,2.0,-1]
	ubs = [3.0,5.0,3]

	gs = Array{Float64,1}[]
	push!(gs, range(lbs[1],stop = ubs[1],length = 3))
	push!(gs, range(lbs[2],stop = ubs[2],length = 4))
	push!(gs, range(lbs[3],stop = ubs[3],length = 5))

	myfun(i1,i2,i3) = 0.5*i1 + 2*i2 + 3*i3

	vs = Float64[ myfun(i,j,k) for i in gs[1], j in gs[2], k in gs[3] ]

	l = Lininterp(vs,gs)

	# check value on bounds
	@test getValue(l,lbs)[1] == vs[1,1,1]
	@test getValue(l,ubs)[1] == vs[3,4,5]
	@test getValue(l,[1.0,5.0,-1])[1] == vs[1,4,1]
	@test getValue(l,[1.0,5.0,3])[1] == vs[1,4,5]

	# check values out of bounds
	@test getValue(l,[-1.0,2.0,-1])[1] == vs[1,1,1]
	@test getValue(l,[1.0,200.0,-1])[1] == vs[1,4,1]
	println(l)
	# ApproXD.resetCache!(l)

	# close to bounds
	x=1.99
	y=4.9
	z=2.9
	@test isapprox(getValue(l,[x,y,z])[1] - myfun(x,y,z) ,0.0,atol=1e-6)
	@test isapprox(getValue(l,[x,y,z])[1] - myfun(x,y,z) ,0.0,atol=1e-6)

	x=1.4861407584066377
	y=3.5646251730324234
	z=2.7832610606532944
	@test isapprox(getValue(l,[x,y,z])[1] - myfun(x,y,z) ,0.0,atol=1e-6)

	x = 2.53
	y = 2.58
	z = -0.87
	@test isapprox(getValue(l,[x,y,z])[1] - myfun(x,y,z) ,0.0,atol=1e-6)

	# check at random vals in interval
	lbs = [1.0,2.0,-10]
	ubs = [33.0,5.0,3]
	gs = Array{Float64,1}[]
	push!(gs, range(lbs[1],stop = ubs[1],length = 30))
	push!(gs, range(lbs[2],stop = ubs[2],length = 40))
	push!(gs, range(lbs[3],stop = ubs[3],length = 50))

	vs = Float64[ myfun(i,j,k) for i in gs[1], j in gs[2], k in gs[3] ]

	l = Lininterp(vs,gs)
	for i in 1:30
		x = rand() * (ubs[1]-lbs[1]) + lbs[1]
		y = rand() * (ubs[2]-lbs[2]) + lbs[2]
		z = rand() * (ubs[3]-lbs[3]) + lbs[3]
		@test isapprox(getValue(l,[x,y,z])[1] ,myfun(x,y,z),atol=1e-6)
	end
	println(l)
end

@testset "testing getValue3D on 2 functions" begin

	lbs = [1.0,2.0,-1]
	ubs = [3.0,5.0,3]

	gs = Array{Float64,1}[]
	push!(gs, range(lbs[1],stop = ubs[1],length = 3))
	push!(gs, range(lbs[2],stop = ubs[2],length = 4))
	push!(gs, range(lbs[3],stop = ubs[3],length = 5))

	myfun1(i1,i2,i3) = 0.5*i1 + 2*i2 + 3*i3
	myfun2(i1,i2,i3) = 3*i1 + pi*i2 + 0.1*i3

	vs1 = Float64[ myfun1(i,j,k) for i in gs[1], j in gs[2], k in gs[3] ]
	vs2 = Float64[ myfun2(i,j,k) for i in gs[1], j in gs[2], k in gs[3] ]

	l = Lininterp(vs1,vs2,gs)

	# check value on bounds
	@test getValue(l,lbs) .- [vs1[1,1,1],vs2[1,1,1]] == zeros(2)
	@test getValue(l,ubs) .- [vs1[3,4,5],vs2[3,4,5]] == zeros(2)
	@test getValue(l,[1.0,5.0,-1])[1] == vs1[1,4,1]
	@test getValue(l,[1.0,5.0,-1])[2] == vs2[1,4,1]
	@test getValue(l,[1.0,5.0,3])[1] == vs1[1,4,5]
	@test getValue(l,[1.0,5.0,3])[2] == vs2[1,4,5]

	# check values out of bounds
	@test getValue(l,[-1.0,2.0,-1])[1] == vs1[1,1,1]
	@test getValue(l,[-1.0,2.0,-1])[2] == vs2[1,1,1]
	@test getValue(l,[1.0,200.0,-1])[1] == vs1[1,4,1]
	@test getValue(l,[1.0,200.0,-1])[2] == vs2[1,4,1]
	println(l)
	# ApproXD.resetCache!(l)

	# close to bounds
	x=1.99
	y=4.9
	z=2.9
	@test isapprox(getValue(l,[x,y,z])[1] - myfun1(x,y,z) ,0.0,atol=1e-6)
	@test isapprox(getValue(l,[x,y,z])[2] - myfun2(x,y,z) ,0.0,atol=1e-6)
	@test isapprox(getValue(l,[x,y,z])[1] - myfun1(x,y,z) ,0.0,atol=1e-6)
	@test isapprox(getValue(l,[x,y,z])[2] - myfun2(x,y,z) ,0.0,atol=1e-6)

	x=1.4861407584066377
	y=3.5646251730324234
	z=2.7832610606532944
	@test isapprox(getValue(l,[x,y,z])[1] - myfun1(x,y,z) ,0.0,atol=1e-6)
	@test isapprox(getValue(l,[x,y,z])[2] - myfun2(x,y,z) ,0.0,atol=1e-6)

	x = 2.53
	y = 2.58
	z = -0.87
	@test isapprox(getValue(l,[x,y,z])[1] - myfun1(x,y,z) ,0.0,atol=1e-6)
	@test isapprox(getValue(l,[x,y,z])[2] - myfun2(x,y,z) ,0.0,atol=1e-6)

	# check at random vals in interval
	lbs = [1.0,2.0,-10]
	ubs = [33.0,5.0,3]
	gs = Array{Float64,1}[]
	push!(gs, range(lbs[1],stop = ubs[1],length = 30))
	push!(gs, range(lbs[2],stop = ubs[2],length = 40))
	push!(gs, range(lbs[3],stop = ubs[3],length = 50))

	vs1 = Float64[ myfun1(i,j,k) for i in gs[1], j in gs[2], k in gs[3] ]
	vs2 = Float64[ myfun2(i,j,k) for i in gs[1], j in gs[2], k in gs[3] ]

	l = Lininterp(vs1,vs2,gs)
	for i in 1:30
		x = rand() * (ubs[1]-lbs[1]) + lbs[1]
		y = rand() * (ubs[2]-lbs[2]) + lbs[2]
		z = rand() * (ubs[3]-lbs[3]) + lbs[3]
		@test isapprox(getValue(l,[x,y,z])[1] ,myfun1(x,y,z),atol=1e-6)
		@test isapprox(getValue(l,[x,y,z])[2] ,myfun2(x,y,z),atol=1e-6)
	end
	println(l)
end

@testset "testing getValue! 3D on 2 functions with ifunc switch" begin

	lbs = [1.0,2.0,-1]
	ubs = [3.0,5.0,3]

	gs = Array{Float64,1}[]
	push!(gs, range(lbs[1],stop = ubs[1],length = 3))
	push!(gs, range(lbs[2],stop = ubs[2],length = 4))
	push!(gs, range(lbs[3],stop = ubs[3],length = 5))

	myfun1(i1,i2,i3) = 0.5*i1 + 2*i2 + 3*i3
	myfun2(i1,i2,i3) = 3*i1 + pi*i2 + 0.1*i3

	vs1 = Float64[ myfun1(i,j,k) for i in gs[1], j in gs[2], k in gs[3] ]
	vs2 = Float64[ myfun2(i,j,k) for i in gs[1], j in gs[2], k in gs[3] ]

	l = Lininterp(vs1,vs2,gs)

	@test_throws MethodError getValue!(l,lbs,[1,2,3])

	# check value on bounds
	y = [0.0]
	getValue!(y,l,lbs,[1])
	@test y == [vs1[1,1,1]]

	getValue!(y,l,lbs,[2])
	@test y == [vs2[1,1,1]]
	getValue!(y,l,ubs,[1])
	@test y == [vs1[3,4,5]]
	getValue!(y,l,ubs,[2])
	@test y == [vs2[3,4,5]]

	getValue!(y,l,[1.0,5.0,-1],[1])
	@test y == [vs1[1,4,1]]
	getValue!(y,l,[1.0,5.0,-1],[2])
	@test y == [vs2[1,4,1]]
	getValue!(y,l,[1.0,5.0,3],[1])
	@test y  == [vs1[1,4,5]]
	getValue!(y,l,[1.0,5.0,3],[2])
	@test y  == [vs2[1,4,5]]

	# check values out of bounds
	getValue!(y,l,[-1.0,2.0,-1],[1])
	@test y == [vs1[1,1,1]]
	getValue!(y,l,[-1.0,2.0,-1],[2])
	@test y == [vs2[1,1,1]]
	y2 = [0.0,0.0]
	getValue!(y2,l,[1.0,200.0,-1],[1,2])
	@test y2 == [vs1[1,4,1],vs2[1,4,1]]
	getValue!(y,l,[1.0,200.0,-1],[2])
	@test y == [vs2[1,4,1]]
	println(l)
	# ApproXD.resetCache!(l)

	# close to bounds
	x=1.99
	y=4.9
	z=2.9
	@test isapprox(getValue(l,[x,y,z])[1] - myfun1(x,y,z) ,0.0,atol=1e-6)
	@test isapprox(getValue(l,[x,y,z])[2] - myfun2(x,y,z) ,0.0,atol=1e-6)
	@test isapprox(getValue(l,[x,y,z])[1] - myfun1(x,y,z) ,0.0,atol=1e-6)
	@test isapprox(getValue(l,[x,y,z])[2] - myfun2(x,y,z) ,0.0,atol=1e-6)

	x=1.4861407584066377
	y=3.5646251730324234
	z=2.7832610606532944
	@test isapprox(getValue(l,[x,y,z])[1] - myfun1(x,y,z) ,0.0,atol=1e-6)
	@test isapprox(getValue(l,[x,y,z])[2] - myfun2(x,y,z) ,0.0,atol=1e-6)

	x = 2.53
	y = 2.58
	z = -0.87
	@test isapprox(getValue(l,[x,y,z])[1] - myfun1(x,y,z) ,0.0,atol=1e-6)
	@test isapprox(getValue(l,[x,y,z])[2] - myfun2(x,y,z) ,0.0,atol=1e-6)

	# check at random vals in interval
	lbs = [1.0,2.0,-10]
	ubs = [33.0,5.0,3]
	gs = Array{Float64,1}[]
	push!(gs, range(lbs[1],stop = ubs[1], length = 30))
	push!(gs, range(lbs[2],stop = ubs[2], length = 40))
	push!(gs, range(lbs[3],stop = ubs[3], length = 50))

	vs1 = Float64[ myfun1(i,j,k) for i in gs[1], j in gs[2], k in gs[3] ]
	vs2 = Float64[ myfun2(i,j,k) for i in gs[1], j in gs[2], k in gs[3] ]

	l = Lininterp(vs1,vs2,gs)
	yout = [0.0,0.0]
	for i in 1:30
		x = rand() * (ubs[1]-lbs[1]) + lbs[1]
		y = rand() * (ubs[2]-lbs[2]) + lbs[2]
		z = rand() * (ubs[3]-lbs[3]) + lbs[3]
		getValue!(yout,l,[x,y,z],[1,2])
		@test isapprox(yout[1] ,myfun1(x,y,z),atol=1e-6)
		@test isapprox(yout[2] ,myfun2(x,y,z),atol=1e-6)
	end
	println(l)
end

@testset "testing getValue hit/miss for 3D" begin

	lbs = [1.0,2.0,-1]
	ubs = [3.0,5.0,3]

	gs = Array{Float64,1}[]
	push!(gs, range(lbs[1],stop = ubs[1], length = 3))
	push!(gs, range(lbs[2],stop = ubs[2], length = 4))
	push!(gs, range(lbs[3],stop = ubs[3], length = 5))

	myfun(i1,i2,i3) = i1 + 2*i2 + 3*i3

	vs = Float64[ myfun(i,j,k) for i in gs[1], j in gs[2], k in gs[3] ]

	l = Lininterp(vs,gs)
	@test hitmiss(l) == (0,0)

	x=1.99
	y=4.9
	z=2.9
	@test isapprox(getValue(l,[x,y,z])[1] - myfun(x,y,z) ,0.0,atol=1e-6)
	@test hitmiss(l) == (0,3)
	@test isapprox(getValue(l,[x,y,z])[1] - myfun(x,y,z) ,0.0,atol=1e-6)
	@test hitmiss(l) == (3,3)
	x=2.99
	@test isapprox(getValue(l,[x,y,z])[1] - myfun(x,y,z) ,0.0,atol=1e-6)
	@test hitmiss(l) == (5,4)
	@test isapprox(getValue(l,[x,y,z])[1] - myfun(x,y,z) ,0.0,atol=1e-6)
	@test hitmiss(l) == (8,4)
	y=2.11
	@test isapprox(getValue(l,[x,y,z])[1] - myfun(x,y,z) ,0.0,atol=1e-6)
	@test hitmiss(l) == (10,5)

end

@testset "testing getValue4D on 3 functions" begin

	lbs = [1.0,2.0,-1,4]
	ubs = [3.0,5.0,3,18.0]

	gs = Array{Float64,1}[]
	push!(gs, range(lbs[1],stop = ubs[1], length = 3))
	push!(gs, range(lbs[2],stop = ubs[2], length = 4))
	push!(gs, range(lbs[3],stop = ubs[3], length = 5))
	push!(gs, range(lbs[4],stop = ubs[4], length = 9))

	myfun1(i1,i2,i3,i4) = 0.5*i1 + 2*i2 + 3*i3 + i4/0.98
	myfun2(i1,i2,i3,i4) = i1 + 0.2*i2 + 3*i3 + i4
	myfun3(i1,i2,i3,i4) = 0.7*i1 + 1.5*i2 + 2*i3 + i4/0.2

	vs1 = Float64[ myfun1(i,j,k,m) for i in gs[1], j in gs[2], k in gs[3], m in gs[4] ]
	vs2 = Float64[ myfun2(i,j,k,m) for i in gs[1], j in gs[2], k in gs[3], m in gs[4] ]
	vs3 = Float64[ myfun3(i,j,k,m) for i in gs[1], j in gs[2], k in gs[3], m in gs[4] ]

		vs = Array{Float64}[]
		push!(vs,vs1)
		push!(vs,vs2)
		push!(vs,vs3)

		l = Lininterp(vs,gs)

	@testset "returning all functions stored on L" begin

		# check value on bounds
		v = getValue(l,lbs)
		@test v[1] == vs1[1,1,1,1]
		@test v[2] == vs2[1,1,1,1]
		@test v[3] == vs3[1,1,1,1]

		v = getValue(l,ubs)
		@test v[1] == vs1[3,4,5,9]
		@test v[2] == vs2[3,4,5,9]
		@test v[3] == vs3[3,4,5,9]

		v = getValue(l,[1.0,5.0,-1,4])
		@test v[1] == vs1[1,4,1,1]
		@test v[2] == vs2[1,4,1,1]
		@test v[3] == vs3[1,4,1,1]

		v = getValue(l,[1.0,5.0,3,18])
		@test v[1] == vs1[1,4,5,9]
		@test v[2] == vs2[1,4,5,9]
		@test v[3] == vs3[1,4,5,9]

		# check values out of bounds
		v = getValue(l,[-1.0,2.0,-1,4])
		@test v[1] == vs1[1,1,1,1]
		@test v[2] == vs2[1,1,1,1]
		@test v[3] == vs3[1,1,1,1]

		v = getValue(l,[1.0,200.0,-1,18])
		@test v[1] == vs1[1,4,1,9]
		@test v[2] == vs2[1,4,1,9]
		@test v[3] == vs3[1,4,1,9]
		println(l)
		# ApproXD.resetCache!(l)

		# close to bounds
		x=1.99
		y=4.9
		z=2.9
		w=17.95
		v = getValue(l,[x,y,z,w])
		@test isapprox(v[1] - myfun1(x,y,z,w) ,0.0,atol=1e-6)
		@test isapprox(v[2] - myfun2(x,y,z,w) ,0.0,atol=1e-6)
		@test isapprox(v[3] - myfun3(x,y,z,w) ,0.0,atol=1e-6)

		x=1.4861407584066377
		y=3.5646251730324234
		z=2.7832610606532944
		w=10.05
		v = getValue(l,[x,y,z,w])
		@test isapprox(v[1] - myfun1(x,y,z,w) ,0.0,atol=1e-6)
		@test isapprox(v[2] - myfun2(x,y,z,w) ,0.0,atol=1e-6)
		@test isapprox(v[3] - myfun3(x,y,z,w) ,0.0,atol=1e-6)

		x = 2.53
		y = 2.58
		z = -0.87
		w = 5.77
		v = getValue(l,[x,y,z,w])
		@test isapprox(v[1] - myfun1(x,y,z,w) ,0.0,atol=1e-6)
		@test isapprox(v[2] - myfun2(x,y,z,w) ,0.0,atol=1e-6)
		@test isapprox(v[3] - myfun3(x,y,z,w) ,0.0,atol=1e-6)

		# check at random vals in interval

		for i in 1:100
			x = rand() * (ubs[1]-lbs[1]) + lbs[1]
			y = rand() * (ubs[2]-lbs[2]) + lbs[2]
			z = rand() * (ubs[3]-lbs[3]) + lbs[3]
			w = rand() * (ubs[4]-lbs[4]) + lbs[4]
			v = getValue(l,[x,y,z,w])
			@test isapprox(v[1] - myfun1(x,y,z,w) ,0.0,atol=1e-6)
			@test isapprox(v[2] - myfun2(x,y,z,w) ,0.0,atol=1e-6)
			@test isapprox(v[3] - myfun3(x,y,z,w) ,0.0,atol=1e-6)
		end
	end

	@testset "returning selected functions stored on L" begin
		f = [myfun1;myfun2;myfun3]
		for i in 1:10
			x = rand() * (ubs[1]-lbs[1]) + lbs[1]
			y = rand() * (ubs[2]-lbs[2]) + lbs[2]
			z = rand() * (ubs[3]-lbs[3]) + lbs[3]
			w = rand() * (ubs[4]-lbs[4]) + lbs[4]
			gets = ApproXD.sample(1:l.nfunc,rand(1:l.nfunc))
			v = zeros(length(gets))
			getValue!(v,l,[x,y,z,w],gets)
			for ig in 1:length(gets)
				@test isapprox(v[ig] - f[gets[ig]](x,y,z,w) ,0.0,atol=1e-6)
			end
		end
	end

end

end


text:
# ---
# name: algorihms.jl
# author: Jhelum Chakravorty, Jayakumar Subramanian, Aditya Mahajan
# date: 1 Feb, 2017
# license: BSD 3 clause
# ---
#
# Introduction
# ============
#
# The code below performs the following calculations. Consider a first-order
# auto-regressive process $\{E_t\}_{t \ge 0}$ which evolves as follows:
# $$ E_{t+1} = a E_t + W_t, $$
# where $a$ is a real number and $\{W_t\}_{t \ge 0}$ is a i.i.d.\ process with
# symmetric and unimodal distribution. In our simulations, we restrict
# attention to $a = 1$ and $W_t$ distributed according to zero-mean Gaussian
# distribution with unit variance. These values are hard-coded in the
# `nextState` function defined below.

@inline function nextState(E)
  a = 1.0
  σ = 1.0
  a*E + σ*randn()
end

# Whenever $|E_t| < k$, where $k$ is a threshold that can be tuned, we incur a
# distortion $d(E_t)$. We hard-code $d(e) = e^2$ in the `distortion` function
# defined below.

@inline distortion(E) = E^2

# Let $S_t$ denote the state of the channel. It is assumed that
# $\{S_t\}_{t \ge 0}$ is an i.i.d. Bernoulli process with $\Pr(S_t = 0) =
# p_d$.
#
# Sampling functions
# ==================
#
# Let $τ$ denote the stopping time when $\{|E_t| \ge k\}$ and $\{S_t = 1\}$.
# The code below performs the following sample-path calculations.
# $$
# L = \sum_{t = 0}^{τ - 1} β^t d(E_t), \quad
# M = \sum_{t = 0}^{τ -1} β^t,         \quad
# K = \sum_{t = 0}^{τ} β^t U_t,
# $$
# where $β$ is the discount factor and $U_t = \IND\{|E_t| \ge k\}$.
#
# Since the calculation is stochastic, there is a positive probability that
# $τ$ is a big number, which can slow down the calculations. So, we set a
# bound `maxIterations` on the maximum size of $τ$. We set the default value
# of `maxIterations` as `10_000`. This number may need to be increased when
# computing the performance for large `threshold`.

@inline function sample(threshold, discount, dropProb; maxIterations = 10_000)

    E, L, M, K = 0.0, 0.0, 0.0, 0.0

    count = 0
    scale = 1.0

    while count <= maxIterations
        channel_on = rand() > dropProb
        transmit   = !(-threshold < E < threshold)
        success    = transmit && channel_on

        if !success
          L += scale * distortion(E)
          K += scale * transmit
          M += scale
        else
          K += scale * transmit
          # Note that we could have written
          #     k += scale
          # because `transmit` is true in this branch of the code.
          break
        end

        scale *= discount
        E = nextState(E)
        count += 1
    end

    (L, M, K)
end

# Mini-batch averaging
# --------------------
#
# The value obtained by one sample path is usually noisy. So, we smoothen it
# out by averaging over a mini-batch. The default size of the mini-batch is
# `1000` iterations.

@inline function sample_average(threshold, discount, dropProb; iterations::Int=1000)

    ell, emm, kay = 0.0, 0.0, 0.0

    for i in 1:iterations
        L, M, K = sample(threshold, discount, dropProb)
        ell += L
        emm += M
        kay += K
    end
    ell /= iterations
    emm /= iterations
    kay /= iterations
    (ell, emm, kay)
end

# Stochastic approximation
# ========================
#
# Stochastic gradient descent for costly communication
# ----------------------------------------------------
#
# It has been shown in the paper that a threshold $k$ is optimal for
# communication cost $λ$ if
# $$ \partial_k C(k) = 0, \quad \hbox{where } C(k) = D(k) + λ N(k), \quad
#  D(k) = \frac{L(k)}{M(k)},  \quad
#  N(k) = \frac{K(k)}{M(k)}. $$
#
# The function `sa_costly` computes the optimal threshold for given values
# of `cost`, `discount`, and `dropProb`. It uses Kiefer Wolfowitz algorithm.
# In particular, the gradient is calculated using finite differences:
# $$ \nabla L \approx \frac{1}{2c} [ L(k + c) - L(k - c) ]. $$
# By default, `c` is set to `0.1`. If we were writing this code for higher
# dimensions, we would replace Kiefer-Wolfowitz with the simultaneous
# perturbation (SPSA) algorithm, which is more sample efficient for higher
# dimensions.
#
# The stochastic approximation iteration starts from an initial guess (the
# parameter `initial`). Its default value is `1.0`. This initialization could
# be useful if we have a reasonable guess for optimal threshold (e.g., the
# exact solution obtained by Fredholm integral equations for the case when
# `dropProb = 0`; see TAC 2017 paper for details).
#
# It is not possible to detect convergence of stochastic approximation
# algorithms. So we run the algorithm for a fixed number of iterations (the
# parameter `iterations`, whose default value is `1_000`).
#
# Stochastic approximation algorithms are sensitive to the choice of learning
# rates. We use ADAM to adapt the learning rates according to the sample
# path. The parameters `decay1`, `decay2`, `alpha`, and `epsilon` can be used
# to tune ADAM. In our experience, these should be left to their default
# values.
#
# It is not possible to detect convergence of stochastic approximation
# algorithms. So we run the algorithm for a fixed number of iterations (the
# parameter `iterations`, whose default value is `1_000`).
#
# Sometimes it is useful to visualize the estimates (of the threshold) as the
# algorithm is running. To do so, set `debug` to `true`, which will print the
# current estimate of the threshold after every 100 iterations.
#
# The output of the function is a trace of the estimates of the
# threshold (therefore, it is a 1D array fo size `iterations`).

@fastmath function sa_costly(cost, discount, dropProb ;
    iterations :: Int     = 1_000,
       initial :: Float64 = 1.0,
        decay1 :: Float64 = 0.9,
	    decay2 :: Float64 = 0.999,
	   epsilon :: Float64 = 1e-8,
        alpha  :: Float64 = 0.01,
             c :: Float64 = 0.1,
	     debug :: Bool    = false,
    )

    threshold = initial
    trace     = zeros(iterations)

    moment1 = 0.0
    moment2 = 0.0

    weight1 = decay1
    weight2 = decay2

    @inbounds for k in 1:iterations
        threshold_plus  = threshold + c
        threshold_minus = threshold - c

        L_plus , M_plus,  K_plus  = sample_average(threshold_plus,  discount, dropProb)
        L_minus, M_minus, K_minus = sample_average(threshold_minus, discount, dropProb)

        C_plus  = (L_plus  + cost*K_plus )/M_plus
        C_minus = (L_minus + cost*K_minus)/M_minus

        gradient = (C_plus - C_minus)/2c

        moment1 = decay1 * moment1 + (1 - decay1) * gradient
        moment2 = decay2 * moment2 + (1 - decay2) * gradient^2

        corrected1 = moment1/(1 - weight1)
        corrected2 = moment2/(1 - weight2)

        weight1 *= decay1
        weight2 *= decay2

        threshold_delta = corrected1 / ( sqrt(corrected2) + epsilon)

        threshold -= alpha * threshold_delta

        # The above analysis is valid for positive values of threshold, but
        # the above line can make threshold negative. Ideally, if the
        # threshold is negative, we should set it to zero. But if the
        # threshold is set to zero, the value of M in the next iteration will
        # also be zero; therefore, the C at the next instant will be infinity.
        # So, we set the minimum value of threshold to be 2c.

        threshold  = max(threshold, 2c)

        if debug && mod(k,100) == 0
          @printf("#:%8d, threshold=%0.6f\n", k, threshold)
        end

        trace[k] = threshold
    end

    return trace
end

# Stochastic gradient descent for constrained communication
# ----------------------------------------------------------
#
# It has been shown in the paper that a threshold $k$ is optimal for a
# communication rate $α$ if
# $$ α M(k) - K(k) = 0. $$
#
# The function `sa_constrained` computes the optimal threshold for a given
# value of `rate`, `discount`, and `dropProb`. It uses Robbins-Monro
# algorithm. This algorithm converges fairly fast, so we do not use
# mini-batches to smoothen out the samples and set the learning rate at
# iteration $n$ to be $\mono{alpha}/n$, where `alpha` is set to `1`, by
# default.
#
# Sometimes it is useful to visualize the estimates (of the threshold) as the
# algorithm is running. To do so, set `debug` to `true`, which will print the
# current estimate of the threshold after every 100 iterations.
#
# The output of the function is a trace of the estimates of the
# threshold (therefore, it is a 1D array fo size `iterations`).

@fastmath function sa_constrained(rate, discount, dropProb ;
    iterations :: Int     = 1_000,
       initial :: Float64 = 1.0,
        alpha  :: Float64 = 1.0,
	     debug :: Bool    = false,
    )

    threshold = initial
    trace     = zeros(iterations)

    @inbounds for k in 1:iterations
        _, M, K = sample_average(threshold, discount, dropProb; iterations=1)

        rl         = alpha/k
        threshold -= rl*(rate * M - K)

        # As in the case of costly communication, the above line can make the
        # threshold negative, so we set its minimum value to be a small
        # positive number (0.05). Note that hard-coding this value means that
        # the code will not work correctly for really small values of the
        # communication rate, where we expect the optimal threshold to be
        # less than 0.05.

        threshold  = max(threshold, 0.05)

        if debug && mod(k,100) == 0
          @printf("#:%8d, threshold=%0.6f\n", k, threshold)
        end

        trace[k] = threshold
    end

    return trace
end

max_stars_repo_name: adityam/remote-estimation-packet-drops

max_issues_repo_name: adityam/remote-estimation-packet-drops

max_forks_repo_name: adityam/remote-estimation-packet-drops