text:
module PSApp

# needed for use of JuliaPaintedItem
ENV["QSG_RENDER_LOOP"] = "basic"

using QML
if ! (:loadqml in names(QML))
    const loadqml = QML.load
    @warn "You seem to be using an old version of QML; problems are likely."
end

using Observables

using ..PseudospectraView: _verbosity, _pbackend
using ..PseudospectraView: _doomed
verbosity() = _verbosity[]

# the App will store return values here:
const results = Dict{Symbol,Any}()

# which entry in QML stacks to use for graphs
live_pane() = _pbackend[] == :pyplot ? 2 : 1

using Plots

using CxxWrap

if _pbackend[] == :pyplot
    # Python must not manage a gui
    ENV["MPLBACKEND"] = "Agg"
    pyplot()
elseif _pbackend[] == :gr
    ENV["GKSwstype"] = "use_default"
    ENV["GKS_WSTYPE"] = 381
    ENV["GKS_QT_VERSION"] = 5
    gr(show=false)
else
    error("backend must be :pyplot or :gr")
end

(verbosity() > 1) && println("Plots loaded, backend is $(typeof(Plots.backend())).")

using Pseudospectra
# shorthand
PSA = Pseudospectra

PseudospectraPlots = Pseudospectra.PseudospectraPlots
PlotsGUIState = Pseudospectra.PseudospectraPlots.PlotsGUIState
function drawcmd end # forward decl
const gs = PlotsGUIState(nothing,0,drawcmd,separate_subplots=(_pbackend[] == :gr))

using PseudospectraView

PSAData = PseudospectraView.PSAData

using Formatting

using Printf: @sprintf
using SparseArrays: issparse
# for schur!
using LinearAlgebra

################################################################
# Global data
# Since the App runs in global scope, and dependencies make this godawful
# slow anyway, we may as well use global vars where it's convenient.
# A few more are created by big function calls below.

maindisplay = nothing
mainpltobj = nothing
const mainsize = [0,0]

otherdisplay = nothing
pltobj2 = nothing
const othersize = [0,0]

# the instance of PSA opts maintained by the GUI
guiopts = Dict{Symbol,Any}()

ps_data = nothing

firstcall = true

@enum DirVsIter dvi_unset dvi_direct dvi_iter
dvi = dvi_unset

const timer = QTimer()

const mcounter = Ref(0)

# flags indicating whether items changed by GUI require action
need_recomp = false
need_redraw = false

# In order to force updates of the GUI window, we use tasks
# and a QTimer. For now we use produce/consume logic, and
# do everything synchronously, in an effort to ensure GUI updates.
compute_chnl = nothing
const running = Ref(false)

const axname2idx = Dict(:xMin=>1,:xMax=>2,:yMin=>3,:yMax=>4)

tmpdata = nothing

# Until we can get a plot location from the pointer and GUI canvas,
# we use a dialog which saves its result here:
zselected = NaN+0.0im

################################################################
# Global state
# Observables for data exchange between GUI and Julia

const o_banner_msg = Observable("Welcome to PseudospectraView")
const o_info_msg = Observable("")
const computation = Observable{Any}(nothing)
const computation_key = Observable(Cint(0))
const o_progress = Observable(0.0)
const o_gridpts = Observable("")
const o_xmin = Observable("")
const o_xmax = Observable("")
const o_ymin = Observable("")
const o_ymax = Observable("")
const o_firstlev = Observable("")
const o_lastlev = Observable("")
const o_levstep = Observable("")
const o_fov = Observable(Cint(0))
const o_mousez = Observable(Cint(0))
const o_imaxis = Observable(Cint(0))
const o_unitcircle = Observable(Cint(0))
const o_arpack_nev = Observable("6")
const o_arpack_ncv = Observable("8")
const o_arpack_which = Observable("LM")
const o_ticks = Observable(Cint(0)) # Number of times the timer has ticked
const o_surf_ok = Observable(_pbackend[] == :pyplot ? Cint(1) : Cint(0))
const o_savefig_ok = Observable(_pbackend[] == :pyplot ? Cint(1) : Cint(0))
const o_timed = Observable(Cint(_doomed[] ? 1 : 0))

propmap = JuliaPropertyMap(
    "computationkey" => computation_key,
    "gridpts" => o_gridpts,
    "xmin" => o_xmin,
    "xmax" => o_xmax,
    "ymin" => o_ymin,
    "ymax" => o_ymax,
    "ticks" => o_ticks,
    "banner_msg" => o_banner_msg,
    "info_msg" => o_info_msg,
    "firstlev" => o_firstlev,
    "lastlev" => o_lastlev,
    "levstep" => o_levstep,
    "fov" => o_fov,
    "mousez" => o_mousez,
    "imag_axis" => o_imaxis,
    "unit_circle" => o_unitcircle,
    "arpack_nev" => o_arpack_nev,
    "arpack_ncv" => o_arpack_ncv,
    "arpack_which" => o_arpack_which,
    "surface_ok" => o_surf_ok,
    "savefig_ok" => o_savefig_ok,
    "timed" => o_timed,
)

"""
attempt to shut down the Qt app cleanly
"""
function panic(msg=nothing)
    if msg === nothing
        @warn "Unspecified fatal error"
    else
        @warn "Fatal error: $msg"
    end
    @emit dieDieDie()
end

mutable struct Monitor
    status::Int # 0:initialized 1:running 2:interrupted 3:done 4:error
    key::Int32
    function Monitor()
        me = new(0,computation_key[])
        on(computation_key) do key
            me.key = key
            me.status = 0
        end
        return me
    end
end

isfinished(c::Monitor) = c.status >= 3

struct Runner
    channel::Channel
    ready::Bool
end

# This is used for debugging the task/channel stuff.
# There should be exactly zero chance for libuv to confuse me.
function cprint(s::String)
    verbosity() > 2 || return
    n = length(s)+1
    fid = 2
    ccall(:write,Clong,(Cint,Cstring,Clong),fid,s,n)
end

function compute_chan(channel::Channel)
    verbosity() > 1 && println("starting computation")
    c = Monitor()
    mcounter[] = 0
    while !isfinished(c)
        run!(c)
        if isready(channel)
            cprint("blocking in cc\n")
        end
        cprint("putting $(c.status) in cc\n")
        put!(channel, c.status)
        cprint("continuing in cc\n")
    end
    (verbosity() > 1) && println("paint calls: $(mcounter[])")
    cprint("cc done\n")
    put!(channel, 4)
    cprint("cc exiting\n")
    nothing
end

const pvals = [0.0, 0.5, 0.99, 1.0, -1.0]
function step(r::Runner)
    if !isready(r.channel)
        cprint("blocking in step\n")
    end
    s = take!(r.channel)
    cprint("taking $s in step\n")
    o_progress[] = pvals[s]
    @emit updateMainPlot()
end

on(o_ticks) do t
    cprint("tick\n")
    if t > 0
        step(computation[])
    end
end

function run!(c::Monitor)
    c.status = 1
    running[] = true
    (verbosity() > 0) && println("calculating...")
    printinfo("calculating...")
    PSA.driver!(ps_data,guiopts,gs,
                myprintln=println # should be printinfo, but we're not getting updates yet
                )
    setedittext(ps_data.zoom_list[ps_data.zoom_pos])
    running[] = false
    (verbosity() > 0) && println("calculation done.")
    printinfo("done.")
    c.status = 3
end

function setup_computation(key)
    computation_key[] = key
    o_progress[] = 0.0
    computation[] = Runner(Channel(compute_chan), false)
end

on(o_progress) do p
    if p >= 1
        # this acts as finalizer for Runner
        println("ticks: ", o_ticks[])
        QML.stop(timer)
        o_ticks[] = Cint(0)
        computation[] = nothing
    end
end

on(o_gridpts) do s
    local n
    n = tryparse(Int, s)
    if n === nothing
        @warn "Invalid entry for npts"
        return
    end
    global need_recomp
    global need_redraw
    if ! (ps_data === nothing)
        curzoom = ps_data.zoom_list[ps_data.zoom_pos]
        old_n = curzoom.npts
        if n >= 5
            if n != old_n
                curzoom.npts = n
                need_recomp = true
                need_redraw = true
                @emit enableGo(Int32(1))
            end
        else
            # this should be obviated by QML logic.
            @warn "npts must be at least 5; got $n"
            o_gridpts[] = max(curzoom.npts, 5)
        end
    end
end

function update_ax(k, s::AbstractString)
    global need_recomp
    global need_redraw
    local x
    x = tryparse(Float64,s)
    if x === nothing
        @warn "Invalid entry for $k"
        return
    end
    (ps_data === nothing) && return
    curzoom = ps_data.zoom_list[ps_data.zoom_pos]
    isempty(curzoom.ax) && (curzoom.ax = fill(NaN,(4,)))
    curzoom.ax[axname2idx[k]] = x
    need_recomp = true
    if !any(isnan.(curzoom.ax)) # goable
        @emit enableGo(Int32(1))
    end
end

on(o_xmin) do s
    update_ax(:xMin, s)
end
on(o_xmax) do s
    update_ax(:xMax, s)
end
on(o_ymin) do s
    update_ax(:yMin, s)
end
on(o_ymax) do s
    update_ax(:yMax, s)
end

function update_levels(k, s)
    global need_redraw
    x = tryparse(Float64, s)
    if x === nothing
        @warn "Invalid entry for contour levels $k"
        return
    end
    curzoom = ps_data.zoom_list[ps_data.zoom_pos]
    setfield!(curzoom.levels, k, x)
    curzoom.autolev = false
    need_redraw = true
    @emit enableGo(Int32(1))
end
on(o_firstlev) do s
    update_levels(:first, s)
end
on(o_lastlev) do s
    update_levels(:last, s)
end
on(o_levstep) do s
    update_levels(:step, s)
end

on(o_arpack_nev) do s
    global need_recomp
    ps_dict = ps_data.ps_dict
    if haskey(ps_dict,:arpack_opts)
        n = tryparse(Int, s)
        if n === nothing
            @warn "Invalid entry for ARPACK nev"
            return
        end
        if (n < 1) || (n > size(ps_data.input_matrix,1)-2)
            @warn "Invalid entry for ARPACK nev"
        end
        if ps_dict[:arpack_opts].nev != n
            ps_dict[:arpack_opts].nev = n
            ps_dict[:proj_valid] = false
            need_recomp = true
            @emit enableGo(Int32(1))
        end
    end
end

on(o_arpack_ncv) do s
    global need_recomp
    ps_dict = ps_data.ps_dict
    if haskey(ps_dict,:arpack_opts)
        n = tryparse(Int, s)
        if n === nothing
            @warn "Invalid entry for ARPACK ncv"
            return
        end
        if (n < 1) || (n > size(ps_data.input_matrix,1)-2)
            @warn "Invalid entry for ARPACK ncv"
        end
        if ps_dict[:arpack_opts].ncv != n
            ps_dict[:arpack_opts].ncv = n
            ps_dict[:proj_valid] = false
            need_recomp = true
            @emit enableGo(Int32(1))
        end
    end
end

on(o_arpack_which) do s
    (verbosity() > 1) && println("attempting to set which to $s")
    global need_recomp
    ps_dict = ps_data.ps_dict
    if s ∈ ["LM", "SM", "LR", "SR", "LI", "SI"]
        if haskey(ps_dict,:arpack_opts)
            if ps_dict[:arpack_opts].which != Symbol(s)
                ps_dict[:arpack_opts].which = Symbol(s)
                ps_dict[:proj_valid] = false
                need_recomp = true
                @emit enableGo(Int32(1))
            end
        end
    else
        @warn "Invalid entry for ARPACK \"which\""
    end
end

on(o_fov) do s
    tf = Bool(s)
    global need_redraw
    if get(guiopts,:showfov,false) != tf
        need_redraw = true
        @emit enableGo(Int32(1))
    end
    guiopts[:showfov] = tf
end

on(o_imaxis) do s
    tf = Bool(s)
    global need_redraw
    if get(guiopts,:showimagax,false) != tf
        need_redraw = true
        @emit enableGo(Int32(1))
    end
    guiopts[:showimagax] = tf
end

on(o_unitcircle) do s
    tf = Bool(s)
    global need_redraw
    if get(guiopts,:showunitcircle,false) != tf
        need_redraw = true
        @emit enableGo(Int32(1))
    end
    guiopts[:showunitcircle] = tf
end

################################################################
# QMLFunctions: how the GUI talks to Julia
# Apparently qmlfunctions need to be defined in global scope
# and QML.exec() needs to run in global scope to see them.

# These are called by the App code, apparently only on setup or resize
# Display is generally forced via drawcmd().

# Somebody doesn't manage dimensions correctly
# This is an attempt at a work-around.
const pyplot_fudge = 0.9

function init_backend(width::Float64, height::Float64, idx)
    if width < 5 || height < 5
        return
    end
    if idx == 1
        mainsize .= round.(Int,[width,height])
    else
        othersize .= round.(Int,[width,height])
    end
    if isa(Plots.backend(),Plots.PyPlotBackend)
        pyplot(size=(round(Int64,pyplot_fudge*width),
                     round(Int64,height)))
    elseif isa(Plots.backend(),Plots.GRBackend)
        gr(size=(round(Int64,width),round(Int64,height)))
        Plots.GR.inline()
    else
        @error "unsupported Plots backend"
    end
    return
end

# This is supposed to help with z-picking. It does, sometimes.
function set_paint_size(width::Float64, height::Float64, idx)
    if width < 5 || height < 5
        return
    end
    if idx == 1
        mainsize .= round.(Int,[width,height])
    else
        othersize .= round.(Int,[width,height])
    end
end

function mainplot(d::JuliaDisplay)
    global maindisplay
    maindisplay = d
    if mainpltobj != nothing
        display(d, mainpltobj)
    end
    return
end

function otherplot(d::JuliaDisplay)
    global otherdisplay
    otherdisplay = d
    if pltobj2 != nothing
        display(d, pltobj2)
    end
    return
end

qmlfunction("mainplot", mainplot)
qmlfunction("init_backend", init_backend)
qmlfunction("otherplot", otherplot)
qmlfunction("size_paint", set_paint_size)

# Functions called by various GUI controls

# setoptbydlg() is called when a GUI dialog accepted a value from user.
# Named for principal use case, but can also invoke a compute/plot routine.
function setoptbydlg(optkey::AbstractString,val::AbstractString)
    k = Symbol(optkey)
    global need_recomp
#    println("setting $optkey by dialog to $val")
    ao = get(ps_data.ps_dict,:arpack_opts,nothing)
    curzoom = ps_data.zoom_list[ps_data.zoom_pos]
    ao_current = true
    if (optkey[1:6] == "Arpack") && ao == nothing
        @warn "no arpack opts in ps_data"
        return
    end
    local n,x

    if k == :ArpackMaxiter
        n = tryparse(Int,val)
        if n === nothing
            @warn "Invalid entry for ARPACK maxiter"
            return
        end
        if ao.maxiter != n
            ao.maxiter = n
            need_recomp = true
            ao_current = false
        end
    elseif k == :ArpackNcv
        n = tryparse(Int,val)
        if n === nothing
            @warn "Invalid entry for ARPACK ncv"
            return
        end
        if ao.ncv != n
            ao.ncv = n
            need_recomp = true
            ao_current = false
        end
    elseif k == :ArpackTol
        x = tryparse(Float64,val)
        if x === nothing
            @warn "Invalid entry for ARPACK tol"
            return
        end
        if !isapprox(ao.tol,x)
            ao.tol = x
            need_recomp = true
            ao_current = false
        end
    elseif k == :AbscissaEps
        x = tryparse(Float64,val)
        if x ===  nothing
            @warn "Invalid entry for ϵ"
            return
        end
        if x < 0
            @warn "Invalid entry for ϵ"
        end
        abscissacb(x)
    elseif k == :RadiusEps
        x = tryparse(Float64,val)
        if x === nothing
            @warn "Invalid entry for ϵ"
            return
        end
        if x < 0
            @warn "Invalid entry for ϵ"
        end
        radiuscb(x)
    elseif k == :ProjLev
        x = tryparse(Float64,val)
        if x === nothing
            @warn "Invalid entry for Projection Level"
            return
        end
        if !isapprox(curzoom.proj_lev, x)
            curzoom.proj_lev = x
            curzoom.computed = false
            ps_data.ps_dict[:proj_valid] = false
            need_recomp = true
        end
    else
        @warn "setoptbydlg: Unexpected key $k"
    end
    # TODO:  (upstream also has v0, v0Cmd)
    if !ao_current
        setarpackgui(ao)
    end
    return
end

function smartlevels()
    global need_redraw
    if ps_data == nothing
        @warn "smartlev invoked before population of ps_data"
        return
    end
    curzoom = ps_data.zoom_list[ps_data.zoom_pos]
    # FIXME: needs much more logic for zooming
    t_levels,err = PSA.recalc_levels(curzoom.Z,curzoom.ax)
    if err == -1
        printwarning("Range is too small: no contours to plot. Refine grid or zoom out.")
        return
    elseif err == -2
        printwarning("Matrix is too non-normal: resolvent norm → ∞. Zoom out!")
        return
    end
    lev = PSA.LevelDesc(t_levels)
    curzoom.levels = lev
    curzoom.autolev = true
    @emit setLevels(lev.first,lev.last,lev.step)
    need_redraw = true
    @emit enableGo(Int32(1))
    nothing
end

function go()
    global need_recomp
    global need_redraw
    if ps_data == nothing
        @warn "go invoked before population"
        return
    end
    @emit showPlot(Int32(0), Int32(live_pane()))
    if need_recomp
        (verbosity() > 1) && println("opening channel for computation task")
        @emit startTimer(Int32(10))
        setup_computation(0)
        need_recomp = false
        need_redraw = false
    elseif need_redraw
        printinfo("redrawing")
        (verbosity() > 1) && println("redrawing")
        Pseudospectra.redrawcontour(gs,ps_data,guiopts)
        need_redraw = false
    end
#    ax = (mainpltobj != nothing) ? PSA.getxylims(mainpltobj) : zeros(0)
    setedittext(ps_data.zoom_list[ps_data.zoom_pos])
    printinfo("ready")
    @emit enableGo(Int32(0))
end

# assorted callbacks

function setselectedz(zrs::AbstractString,zis::AbstractString)
    global zselected
    zr,zi = NaN,NaN
    zr = tryparse(Float64,zrs)
    zi = tryparse(Float64,zis)
    if zr === nothing || zi === nothing
        printwarning("invalid z")
        zr,zi = NaN,NaN
    end
    zselected = zr + zi*im
#    notify(dialogcv)
    nothing # returns to QML-land, which doesn't like Complex
end

function loadmtx(varname::AbstractString,myexpr::AbstractString)
    #  any existing plots are obsolete, so clear
    global mainplotobj, plotobj2
    mainplotobj = nothing
    plotobj2 = nothing
    @emit showPlot(Int32(0), Int32(0))
    @emit showPlot(Int32(1), Int32(0))

    # this should really be done by the GUI, but this is easier:
    if isempty(varname)
        varname = "A"
    end

    global firstcall
    if firstcall
        printinfo("Please be patient, plotting code needs to be compiled...")
    end

    try
        Core.eval(PSAData, Meta.parse("$(varname) = $(myexpr)"))
        Core.eval(PseudospectraView.PSApp,
             Meta.parse("ps_data = Pseudospectra.new_matrix(PSAData.$(varname),guiopts)"))
        (verbosity() > 1) && println("Matrix has been digested.")
    catch JE
        printwarning("New Matrix failed. See console.")
        @warn "loadmtx failed, exception was $JE"
        return nothing
    end
    global guiopts
    printbanner("Working on $(varname) = $(myexpr)")
    PSAData.clearopts() # delete obsolete matrix-specific stuff
    guiopts = PSA.fillopts(gs,PSAData.getopts())
    refresh()
    firstcall = false
    (verbosity() > 0) && println("matrix loaded")
    nothing
end

function savefigx(filename::AbstractString, mykey)
    if mykey == "portrait"
        try
            Plots.savefig(gs.mainph, filename)
        catch JE
            println("savefig failed; exception:")
            println(JE)
        end
    else
        @warn "unrecognized key"
    end
end

function savedata(varname::AbstractString,mykey)
    # FIXME: test for valid variable name
    if mykey == "portrait"
        tmpdata = ps_data.zoom_list[ps_data.zoom_pos]
    elseif mykey == "ps_data"
        tmpdata = ps_data
    elseif mykey == "zpsradius"
        tmpdata = ps_data.ps_dict[:zpsradius]
    elseif mykey == "zpsabscissa"
        tmpdata = ps_data.ps_dict[:zpsabscissa]
    else
        @warn "unrecognized save key $mykey"
        return
    end
    try
        results[Symbol(varname)] = deepcopy(tmpdata)
    catch JE
        @warn "savedata failed, exception was $JE"
    end
    tmpdata = nothing
end

function use_eigs(x::Int32)
    global dvi
    direct = (x == 1)
    ps_dict = ps_data.ps_dict
    ps_dict[:direct] = direct
    m,n = size(ps_data.input_matrix)
    if direct
        dvi = dvi_direct
        if !isempty(get(ps_dict,:schur_mtx,zeros(0)))
            # revert to input matrix if dense square
            ps_data.matrix = ps_dict[:schur_mtx]
            ps_dict[:ews] = ps_dict[:orig_ews]
            ps_dict[:ew_estimates] = false
        elseif (m==n) && issparse(ps_data.input_matrix)
            # if sparse square, forget ARPACK compression
            ps_data.matrix = ps_data.input_matrix
        end
        # if there's a unitary matrix from Schur decomposition, use that
        # otherwise use the one give by the user, if any
        if !isempty(get(ps_dict,:schur_unitary_mtx,zeros(0)))
            ps_data.unitary_mtx = (ps_data.input_unitary_mtx
                                     * ps_dict[:schur_unitary_mtx])
        else
            ps_data.unitary_mtx = ps_data.input_unitary_mtx
        end
        ps_dict[:proj_valid] = false
        # if we're reverting to a square matrix, forget ARPACK projection
        ss = size(ps_data.matrix)
        if ss[1] == ss[2]
            ps_dict[:isHessenberg] = false
        end
        if (ps_data.zoom_list[ps_data.zoom_pos].computed
            && !isempty(get(ps_dict,:schur_mtx,zeros(0))))
            # CHECKME: what's the reasoning here?
            @emit toggleFoV(Int32(1))
        end
    else
        if dvi != dvi_iter
            ps_dict[:proj_valid] = false
        end
        dvi = dvi_iter
        if !haskey(ps_dict,:arpack_opts)
            ps_dict[:arpack_opts] = PSA.ArpackOptions{eltype(ps_data.matrix)}()
            @warn "setting default ARPACK options (from GUI)"
        end
        @emit toggleFoV(Int32(0))
    end
    global need_recomp
    need_recomp = true
    @emit enableGo(Int32(1))
end

function zmoused(leftright::Int32, xpxl::Float64, ypxl::Float64, yref::Float64)
    if isa(Plots.backend(),Plots.PyPlotBackend)
        fig = gs.mainph.o
        ax_plt = fig.axes[1]
        inv = ax_plt.transData.inverted()
        xscaled = xpxl * pyplot_fudge
        xy = inv.transform((xscaled,yref+1-ypxl))
    elseif isa(Plots.backend(),Plots.GRBackend)
        fig = Plots.GR.gcf()
        nc, nr = mainsize # fig[:size]
        # println("GR thinks size is $nc $nr")
        fac = (nc > nr) ? (1.0 / nc) : (1.0 / nr)
        xndc = fac * xpxl
        yndc = fac * (yref+1-ypxl)
        xy = Plots.GR.ndctowc(xndc, yndc)
    else
        printwarning("zooming or point selection not implemented for $(typeof(Plots.backend()))")
        return nothing
    end
    (verbosity() > 1) && println("selected $xy from yvals $yref $ypxl")
    z = xy[1]+1im*xy[2]
    if true
        zooming_in = (leftright == 0)
        mainzoom(zooming_in, z)
    else
    end
end

function mainzoom(zooming_in::Bool, z)
    ax = PseudospectraPlots.getxylims(gs.mainph)
    if zooming_in
        res = zoomin!(ps_data,z,ax)
    else
        res = zoomout!(ps_data,z,ax,include_fov=get(guiopts,:showfov,false))
    end
    if res < 0
        printwarning("unable to zoom based on requested point")
    elseif res == 0
        # redraw existing portrait
        printinfo("redrawing...")
        Pseudospectra.redrawcontour(gs, ps_data, guiopts)
        printinfo("done")
    else
        # compute and draw new portrait
        printinfo("recomputing...")
        driver!(ps_data, guiopts, gs)
        printinfo("done")
    end
    setedittext(ps_data.zoom_list[ps_data.zoom_pos])
    nothing
end

qmlfunction("loadmtx", loadmtx)
qmlfunction("savedata", savedata)
qmlfunction("savefigx", savefigx)
qmlfunction("go", go)
qmlfunction("use_eigs", use_eigs)
qmlfunction("setoptbydlg", setoptbydlg)
qmlfunction("setselectedz", setselectedz)
qmlfunction("zmoused", zmoused)

#=
# This looks clever but deadlocks
function zdlg()
    @emit getZ("Select z for pseudomode")
end

function guizdlg(gs; what="a z-value")
    global dialogcv
    global zselected
    t = Task(zdlg)
    schedule(t)
    wait(dialogcv)
    println("dialog selection: $zselected")
    zselected
end
=#
function getguiz(gs; what="a z-value")
    zselected
end

function abscissacb(epsilon)
    f,z = PSA.psa_abscissa(ps_data.input_matrix,epsilon)
    println("abscissa(ϵ = $epsilon): $f at ",z)
    msg = "abscissa(ϵ = $epsilon) = $f at $(z[1])"
    printinfo(msg)
    ps_data.ps_dict[:zpsabscissa] = (ϵ = epsilon, x = f, z = z)
    @emit saveVariable("Save point(s) at pseudospectral abscissa","zpsabscissa")
end

function radiuscb(epsilon)
    f,z = PSA.psa_radius(ps_data.input_matrix,epsilon)
    println("radius(ϵ = $epsilon): $f at ",z)
    msg = "radius(ϵ = $epsilon) = $f at $(z[1])"
    printinfo(msg)
    ps_data.ps_dict[:zpsradius] = (ϵ = epsilon, r = f, z = z)
    @emit saveVariable("Save point(s) at pseudospectral radius","zpsradius")
end

function ploteigv()
    if isnan(zselected)
        printwarning("Select z first.")
    else
        # TODO: wrap in gui stupefy/waken
        PSA.modeplot(ps_data, 0; gs=gs, zgetter=getguiz)
        @emit showPlot(Int32(1), Int32(live_pane()))
    end
end
function plotpmode()
    if isnan(zselected)
        printwarning("Select z first.")
    else
        # TODO: wrap in gui stupefy/waken
        PSA.modeplot(ps_data, 1; gs=gs, zgetter=getguiz)
        @emit showPlot(Int32(1), Int32(live_pane()))
    end
end

function mtxpwrcb(nmax)
    flag = false
    try
        PSA.mtxpowersplot(gs,ps_data,nmax,gradual=false)
    catch JE
        println("Matrix power exception:")
        println(JE)
        flag = true
    end
    if flag
        printwarning("Matrix power comp/plot failed. See console.")
    else
        @emit showPlot(Int32(1), Int32(live_pane()))
    end
end

function mtxexpcb(dt,nmax)
    flag = false
    try
        PSA.mtxexpsplot(gs,ps_data,dt,nmax, gradual=false)
    catch JE
        println("Matrix exponential exception:")
        println(JE)
        flag = true
    end
    if flag
        printwarning("Matrix exp comp/plot failed. See console.")
    else
        @emit showPlot(Int32(1), Int32(live_pane()))
    end
end

function plotps3d()
    PSA.surfplot(gs, ps_data, guiopts)
    @emit showPlot(Int32(1), Int32(live_pane()))
end

qmlfunction("ploteigv", ploteigv)
qmlfunction("plotpmode", plotpmode)
qmlfunction("mtxpwrcb", mtxpwrcb)
qmlfunction("mtxexpcb", mtxexpcb)
qmlfunction("plotps3d", plotps3d)
qmlfunction("smartlevels", smartlevels)

# these will be wrapped below
function paint_main(p::CxxPtr{QPainter}, item::CxxPtr{JuliaPaintedItem})
    (_pbackend[] == :gr) || return
    ENV["GKS_CONID"] = split(repr(p.cpp_object), "@")[2]
    dev = device(p[])[]
    r = effectiveDevicePixelRatio(window(item[])[])
    w, h = width(dev) / r, height(dev) / r
    if running[]
        c = mcounter[] + 1
        mcounter[] = c
    end
    if (gs.mainph === nothing)
        (verbosity() > 1) && println("null paint_main called")
        return
    end
    if mainsize[1] > 0
        plot!(gs.mainph, size=(w,h)) # Tuple(mainsize))
        display(gs.mainph)
    end
    return
end
function paint_second(p::CxxPtr{QPainter}, item::CxxPtr{JuliaPaintedItem})
    (_pbackend[] == :gr) || return
    ENV["GKS_CONID"] = split(repr(p.cpp_object), "@")[2]
    dev = device(p[])[]
    r = effectiveDevicePixelRatio(window(item[])[])
    w, h = width(dev) / r, height(dev) / r
    if (gs.ph2 === nothing)
        (verbosity() > 1) && println("null paint_second called")
        return
    end
    if othersize[1] > 0
        plot!(gs.ph2, size=(w,h))
        # for a while, we needed this:
        # px = plot(gs.ph2[1],layout=(1,1))
        # px = plot(gs.ph2[1])
        px = gs.ph2
        display(px)
    end
    return
end
function paint_third(p::CxxPtr{QPainter}, item::CxxPtr{JuliaPaintedItem})
    (_pbackend[] == :gr) || return
    ENV["GKS_CONID"] = split(repr(p.cpp_object), "@")[2]
    dev = device(p[])[]
    r = effectiveDevicePixelRatio(window(item[])[])
    w, h = width(dev) / r, height(dev) / r
    if (gs.ph3 === nothing)
        (verbosity() > 1) && println("null paint_third called")
        return
    end
    if othersize[1] > 0
        plot!(gs.ph3, size=(w,h))
        # px = plot(gs.ph3[1],layout=(1,1))
        # px = plot(gs.ph3[1])
        px = gs.ph3
        display(px)
    end
    return
end

################################################################
# functions to be invoked from Julialand

# drawing command, invoked as a callback inside PSA
function drawcmd(gs::PlotsGUIState,ph,id)
    global mainpltobj
    global pltobj2
    if id == 1
        if _pbackend[] == :pyplot
            mainpltobj = ph
            if maindisplay != nothing
                display(maindisplay, ph)
            end
        else
            if running[]
                # There is often a call before anything is ready.
                # AFAICT it is not a problem.
                if !(computation[] === nothing)
                    if isready(computation[].channel)
                        cprint("blocking in drawcmd\n")
                    end
                    cprint("putting 2 in drawcmd\n")
                    put!(computation[].channel, 2)
                    cprint("continuing in drawcmd\n")
                end
            else
                @emit updateMainPlot()
            end
        end
    else
        if _pbackend[] == :pyplot
            pltobj2 = ph
            if otherdisplay != nothing
                display(otherdisplay, ph)
            end
        else
            if id == 2
                gs.ph2 = ph
                @emit updatePlot2()
            elseif id == 3
                gs.ph3 = ph
                @emit updatePlot3()
            else
                @warn "unimplemented plot destination id=$id"
            end
        end
    end
end

"""
format a pair of reals so that they are distinguishable (if distinct)
and short (if possible)
"""
function formatpair(x1::Real,x2::Real)
    dx = abs(x2-x1)
    if !isfinite(dx) || (dx == 0) # punt
        s1 = @sprintf("%g",x1)
        s2 = @sprintf("%g",x2)
    else
        xscale = max(abs(x1),abs(x2))
        ndig = 2+ceil(Int,log10(xscale/dx))
        # this gave world-age problems (revisit if interpolation yields bad fmt)
        # "%." * @sprintf("%dg",ndig)
        fmt = "%.$(ndig)g"
        s1 = sprintf1(fmt,x1)
        s2 = sprintf1(fmt,x2)
    end
    return [s1,s2]
end

# This is a remedy for the annoying behavior of
# Qt toString(): 1/10 -> 0.100000000001 etc.
formataxislims(ax) = vcat(formatpair(ax[1],ax[2]), formatpair(ax[3],ax[4]))

# fills in the main text fields
function setedittext(zoom)
    ax = zoom.ax
    if !isempty(ax)
        @emit setAxisLims(formataxislims(ax)...)
    else
        @emit setAxisLims("NaN","NaN","NaN","NaN")
    end
    if zoom.levels.isunif
        @emit setLevels(zoom.levels.first,zoom.levels.last,zoom.levels.step)
    else
        # println("Not setting levels fields; zoom is")
        # println(zoom)
        @emit setLevels(NaN,NaN,NaN)
    end
    @emit setGrid(Int32(zoom.npts))
    if dvi != dvi_unset
        @emit toggleDirect(dvi == dvi_direct ? Int32(1) : Int32(0))
    end
    global need_recomp
    global need_redraw
    need_recomp = false
    need_redraw = false
end

# Warning: must be coordinated w/ ListModel in QML file
const which2idx = Dict{Symbol,Int}(:LM => 0, :LR => 1, :SR => 2,
                                   :LI => 3, :SI => 4)
function setarpackgui(ao)
    @emit setArpackOpts(Int32(ao.nev),
                        Int32(ao.ncv),
                        Int32(ao.maxiter),
                        ao.tol,
                        Int32(which2idx[ao.which]))
end

"""
place some text in the main App banner
"""
function printbanner(msg::AbstractString)
    o_banner_msg[] = msg
end

"""
place some text into the informational text area
"""
function printinfo(msg::AbstractString)
    o_info_msg[] = msg
    # TODO: set an info color or icon
end

"""
place a warning into the informational text area
"""
function printwarning(msg::AbstractString)
    o_info_msg[] = msg
    # TODO: set a warning color or icon
end

"""
called after a new matrix is available in ps_data
"""
function refresh()
    global mainplotobj,plotobj2,need_recomp
    mainplotobj,plotobj2 = nothing,nothing
    # Fixme: clear canvasses
    m,n = size(ps_data.input_matrix)
    zoom = ps_data.zoom_list[ps_data.zoom_pos]
    ps_dict = ps_data.ps_dict
    ax = zoom.ax
    if ps_dict[:direct]
        # FIXME: need a function to check whether ew's will be available
        if ((issparse(ps_data.input_matrix)
             && get(ps_dict,:sparse_direct,false))
            || (m != n)
            || length(methods(schur!,(typeof(ps_data.input_matrix),))) == 0
            )
            # rectangular, sparse-direct, or just plain weird
            @emit toggleAutoAx(Int32(-1))
            if isempty(ax)
                # flag need to get ax
                printinfo("Axis limits must be specified to proceed")
                # FIXME: should wait until they are ready
                @emit enableGo(Int32(0))
                need_recomp = true
            end
        else
            # should be ready for this
            @emit toggleDirect(Int32(1))
            # WONTFIXME: this almost certainly won't display unless
            # we wrap it in a QTimer-monitored task. I don't care.
            if isempty(ax) && !isempty(get(ps_dict,:ews,[]))
                PSA.origplot!(ps_data,guiopts,gs)
            else
                PSA.driver!(ps_data,guiopts,gs)
            end
            @emit showPlot(Int32(0), Int32(live_pane()))
            @emit enableGo(Int32(0))
        end
    else # iterative
        @emit toggleDirect(Int32(0))
        if !haskey(ps_dict,:arpack_opts)
            panic("ps_data needs ARPACK opts")
        end
        setarpackgui(ps_dict[:arpack_opts])
        # If the user specified (some) ARPACK options, we might be ready to go,
        # but esp. w/o a stop button, premature start is more likely and
        # more troublesome.

        # FIXME: implement a "clear" command, so we don't get junk while starting up
        begin
            printinfo("Set/check ARPACK options before proceeding")
            need_recomp = true
            # Since we provide defaults, one may proceed immediately.
            @emit enableGo(Int32(1))
        end
    end
    setedittext(ps_data.zoom_list[ps_data.zoom_pos])
end

################################################################
# This is the main App execution sequence.
# QML seems to require that key data structures be in global scope.

# Load up and check data first so if PSA throws an exception
# we are not inchoate.

if !isempty(PSAData.getdefaultmtx())
    (verbosity() > 0) && println("Using existing matrix from PSAData")
    # guiopts is empty until now
    merge!(guiopts,PSA.fillopts(gs,PSAData.opts))
    ps_data = PSA.new_matrix(PSAData.getdefaultmtx(),guiopts)
    (verbosity() > 0) && println("Matrix has been digested.")
    # forget opts that should not apply to newly ingested data
    # viz. when new_matrix() is invoked from GUI.
    PSAData.clearopts(false)
end

# Now start up the engine

qml_file = joinpath(dirname(@__FILE__), "psagui.qml")

loadqml(qml_file, timer=timer, observables=propmap,
     paint_main_wrapped = @safe_cfunction(paint_main, Cvoid,
                                          (CxxPtr{QPainter}, CxxPtr{JuliaPaintedItem})),
     paint_second_wrapped = @safe_cfunction(paint_second, Cvoid,
                                          (CxxPtr{QPainter}, CxxPtr{JuliaPaintedItem})),
     paint_third_wrapped = @safe_cfunction(paint_third, Cvoid,
                                          (CxxPtr{QPainter}, CxxPtr{JuliaPaintedItem}))
     )

if !isempty(PSAData.getdefaultmtx())
    refresh()
end

function runme()
    @emit setInterval(Cint(100))
    # Run the application
    QML.exec()
    return results
end

end


text:
using DSP: conv
using LinearAlgebra

Funcs3D = (
    (x, y, z) -> cat(x, y, z,dims=4),
    (x, y, z) -> cat(
        sqrt(3) * x .* z,
        sqrt(3) * x .* y,
        y .^ 2 - (x .^ 2 + z .^ 2) / 2,
        sqrt(3) * y .* z,
        (z .^ 2 - x .^ 2) * sqrt(3) / 2,
    dims=4),
)

 Funcs2D = (
    (x, y) -> [x, y],
    #     # (x,y)->        [
    #     #             sqrt(3) * x*z,
    #     #             sqrt(3) * x*y,
    #     #             y^ 2 - (x^ 2 + z^ 2) / 2,
    #     #             sqrt(3) *y*z,
    #     #             (z^ 2 - x^ 2) * sqrt(3) / 2,
    #     #         ],
)

function make_spherical_harmonics(dx, size, origin, rank_max; dims = 3)
    if rank_max < 0
        return [], []
    end
    if dims == 3
        rvecs =
            [
                [x, y, z] .- origin for x = 1:size[1], y = 1:size[2],
                z = 1:size[3]
            ] * dx
        Funcs = Funcs3D
    elseif dims == 2
        rvecs = [[x, y] .- origin for x = 1:size[1], y = 1:size[2]] * dx
        Funcs = Funcs2D
    end
    r = norm.(rvecs)
    if rank_max > 0
        rhat = rvecs ./ broadcast(x -> x ≈ 0.0 ? 1.0 : x, r)
        rhat = [getindex.(rhat, i) for i = 1:dims]
    end
    sh = [Funcs[i](rhat...) for i = 1:rank_max]
    r, sh
end

ext: jl


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

using Markdown
using InteractiveUtils

# ╔═╡ bbd5a62e-fbcd-11ea-2bbd-e7be1b736f9f
using Plots

# ╔═╡ afa5cbc0-fbcc-11ea-0dd2-35821062dd18
c = 3

# ╔═╡ 59ddff8e-fbcd-11ea-1124-d308b43c24d9
acr =10

# ╔═╡ d8f55b7e-fbcc-11ea-1881-15dbeb706e20
xvars = [float(c) - i for i in 0:acr]

# ╔═╡ 1c685b10-fbcd-11ea-1ecd-75e8c03e7f87
begin
	function f(x)
		return cos(x)
	end
end

# ╔═╡ 29568700-fbcf-11ea-3096-7f18898aeba0
begin
	function slope(xl,xr,yl,yr,x)
		m = (yr-yl)/(xr-xl)
		b = yl - m*xl
		return m*x + b
	end
end

# ╔═╡ 358ecd40-fbcd-11ea-00fe-755c43fd372e
yvals = [f(xvars[i]) for i in 1:acr]

# ╔═╡ 55799400-fbcd-11ea-05dc-85c5e87f6f63
secs = [(yvals[i]-f(c))/(xvars[i]-c) for i in 1:acr]

# ╔═╡ ae19b810-fbcd-11ea-2ee7-a9f08af0abaf
begin
	plt = plot((x)->f(x),-2c,2c, title ="Limits of secant lines at x=$c",lw=3,legend=false)
	xlabel!("x")
	ylabel!("f(x) ")
	for i in 1:acr
		plot!((x)->slope(xvars[i],c,yvals[i],f(c),x),xvars[i]-1,c+1)

	end
	plt
	savefig("cos(x)-secants.png")
end

# ╔═╡ 12bbbe30-fbce-11ea-289f-ff4eeafbdcbf
cd(#= type where you want the images to save on your machine=#)

# ╔═╡ Cell order:
# ╠═bbd5a62e-fbcd-11ea-2bbd-e7be1b736f9f
# ╠═afa5cbc0-fbcc-11ea-0dd2-35821062dd18
# ╠═59ddff8e-fbcd-11ea-1124-d308b43c24d9
# ╠═d8f55b7e-fbcc-11ea-1881-15dbeb706e20
# ╠═1c685b10-fbcd-11ea-1ecd-75e8c03e7f87
# ╠═29568700-fbcf-11ea-3096-7f18898aeba0
# ╠═358ecd40-fbcd-11ea-00fe-755c43fd372e
# ╠═55799400-fbcd-11ea-05dc-85c5e87f6f63
# ╠═ae19b810-fbcd-11ea-2ee7-a9f08af0abaf
# ╠═12bbbe30-fbce-11ea-289f-ff4eeafbdcbf


text:
Pkg.add("Weber",{{version}})
if !isfile("calibrate.jl")
  open("calibrate.jl","w") do s
    print(s,"""
    # call run_calibrate() to select an appropriate attenuation.
    const atten_dB = 30

    # call Pkg.test(\"Weber\"). If the timing test fails, increase
    # moment resolution to avoid warnings.
    const moment_resolution = 0.0015
    """)
  end
end


text:
module GFX

include("gfx_node_description.jl")
include("gfx_rules_prototypes.jl")
include("gfx_update_rules.jl")

end

max_stars_repo_name: wmkouw/gfx-node

max_issues_repo_name: wmkouw/GFX

max_forks_repo_name: wmkouw/GFX


text:
using Test
using MetaArrays
using StaticArrays
using AxisArrays
using Unitful

struct TestMerge
  val::Int
end
MetaArrays.metamerge(x::TestMerge,y::TestMerge) = TestMerge(x.val + y.val)

struct TestIndex <: AbstractArray{Float64,1}
  x::Vector{Float64}
end
Base.similar(val::TestIndex,::Type{S},dims::NTuple{<:Any,Int}) where S =
  TestIndex{typeof(val.x),eltype(val),ndims(val)}(similar(val.x))
Base.size(val::TestIndex) = size(val.x)
Base.getindex(val::TestIndex,i::Int...) = getindex(val.x,i...)
# not realy a useful function, just a test:
Base.getindex(val::TestIndex,i::Float64) = val[floor(Int,i)] + val[ceil(Int,i)]

struct TestMergeFail
  val::Int
end

struct TestArray{A,T,N} <: AbstractArray{T,N}
  val::A
end
Base.similar(val::TestArray,::Type{S},dims::NTuple{<:Any,Int}) where S =
  TestArray{typeof(val.val),eltype(val),ndims(val)}(similar(val.val))
Base.size(val::TestArray) = size(val.val)
Base.getindex(val::TestArray,i::Int...) = getindex(val.val,i...)
Base.setindex!(val::TestArray,v,i::Int...) = setindex!(val.val,v,i...)

testunion(x::MetaUnion{AbstractRange}) = :range
testunion(x) = :notrange

# TODO: add some more robust tests of broadcast machinery
# (greater variety of arguments, with more call variants
# and then check for test coverage

@testset "MetaArrays" begin

  @testset "MetaArray handles standard array operations" begin
    data = collect(1:10)
    x = meta(data,val=1)
    y = meta(collect(1:10),val=1)

    @test (x.^2) == (1:10).^2
    @test x[1] == data[1]
    @test size(x) == size(data)
    @test similar(x) isa MetaArray
    @test x[1:5] == data[1:5]
    @test (y[1:5] .= 1; sum(y[1:5]) == 5)
    @test (y[1] = 2; y[1] == 2)
    @test x .+ (1:10) == data .+ (1:10)
    @test (.-x) isa MetaArray
    @test collect(1:10) .+ x .+ collect(11:20) == (13:3:40)
    @test broadcast(+,collect(1:10),x,collect(11:20)) == (13:3:40)
    @test (1:10) .+ meta(1:10,val=1) .+ (11:20) isa MetaArray
    @test view(x,1:2) isa MetaArray

    x = TestIndex(collect(1:5))
    @test x[1] == 1
    @test x[1.5] == 3
  end

  @testset "MetaArray preserves metadata over array operations" begin
    data = collect(1:10)
    x = meta(data,val=1)

    @test x.val == 1
    @test x[1:5].val == x.val
    @test x[:].val == x.val
    @test (x .+ (1:10)).val == x.val
    @test (x .+= (1:10); x.val == 1)
    @test (broadcast(+,x,1:10).val == x.val)
    @test similar(x).val == x.val
    @test (.-x).val == 1

    x = meta(collect(1:10),val=1)
    y = meta(collect(1:10),val=1)
    @test vcat(x,y) isa MetaArray
    @test hcat(x,y) isa MetaArray
    @test x*y' isa MetaArray

    x = collect(1:10)
    y = meta(collect(1:10),val=1)
    @test hcat(x,y) isa Array
    @test hcat(y,x) isa MetaArray
    @test vcat(x,y) isa Array
    @test vcat(y,x) isa MetaArray
    @test x*y' isa Array
    @test y*x' isa MetaArray
  end

  @testset "MetaArray properly handles strides" begin
    @test_throws MethodError strides(meta(1:10,val=1))
    x = meta(collect(1:10),val=1)
    @test strides(x) == (1,)
    @test strides(x)[1] == stride(x,1)
  end

  @testset "MetaUnion dispatches properly" begin
    @test testunion(1:5) == :range
    @test testunion(meta(1:10,val=1)) == :range
    @test testunion(meta(collect(1:10),val=1)) == :notrange
  end

  @testset "MetaArray properly merges metadata" begin
    x = meta(collect(1:10),val=1)
    y = meta(collect(2:11),string="string")
    @test (x.+y).val == 1
    @test (x.+y).string == "string"
    @test (x .+ (1:10)).val == 1
    @test ((1:10) .+ x).val == 1

    x = MetaArray(Dict(:test => "value"),collect(1:10))
    y = MetaArray(Dict(:test => "value"),collect(2:11))
    @test getmeta((x .+ y))[:test] == "value"

    x = meta(collect(1:10),val=1)
    y = meta(collect(1:10),val=2)
    @test_throws ErrorException x.+y

    x = meta(collect(1:10),val=TestMerge(1))
    y = meta(collect(1:10),val=TestMerge(2))
    @test (x.+y).val == TestMerge(3)
    @test broadcast(+,x,y).val == TestMerge(3)

    x = meta(collect(1:10),val=(joe=2,bob=3))
    y = meta(collect(1:10),val=(bill=4,))
    @test (x.+y).val == (joe=2,bob=3,bill=4)
  end

  f(x) = (x[0u"s" .. 0.5u"s"] .= 0)

  @testset "MetaArray is AxisArray friendly" begin
    x = meta(AxisArray(rand(10,10),Axis{:time}(range(0u"s",1u"s",length=10)),
      Axis{:freq}(1:10)),val=1)
    x[0u"s" .. 0.5u"s"] = fill(0.0,size(x[0u"s" .. 0.5u"s"]))
    @test all(x[1:5] .== 0)
    x[0u"s" .. 0.5u"s"] .= 1
    @test all(x[1:5] .== 1)

    @test axisdim(x,Axis{:time}) == 1
    @test AxisArray(x) isa AxisArray
    @test AxisArrays.axes(x)[1] isa Axis{:time}
    @test AxisArrays.axes(x,1) isa Axis{:time}
    @test AxisArrays.axes(x,Axis{:time}) isa Axis{:time}
    @test axisnames(x) == (:time,:freq)
    @test axisvalues(x)[2] == 1:10
    @test size(similar(x)) == size(x)
  end

  @testset "MetaArray preserves broadcast specialization" begin
    # Range objects define a specialized lazy broadcasting style,
    # we use them to test the presevation of this lazy style
    x = meta(1:10,val=1)
    @test (x .+ 4) isa MetaArray{<:AbstractRange}
    @test (x .+ 4) == [xi+4 for xi in x]
    @test (.-x) isa MetaArray{<:AbstractRange}
    @test (1:10) .+ meta(1:10,val=1) .+ (11:20) isa MetaArray{<:AbstractRange}
    @test .-x == [-xi for xi in x]

    # SVector defines a specialized style for broadcasting, we use it to test
    # the handling of specialized broadcasting styles
    x = meta(SVector(1,2,3),val=1)
    y = meta([1,2,3],val=1)
    @test x.+y == [2,4,6]
    @test (x.+y).val == 1

    x = meta(collect(1:10),val=1)
    y = SVector((1:10)...)
    @test (x .+ y).val == 1
    @test (y .+ x).val == 1

    y = meta(y,val=1)
    x = collect(1:10)
    @test (x .+ y).val == 1
    @test (y .+ x).val == 1
  end

  @testset "MetaArray allows custom metadata type" begin
    x = MetaArray(TestMerge(2),1:10)
    y = MetaArray(TestMerge(3),1:10)
    k = MetaArray(TestMergeFail(1),1:10)
    h = MetaArray(TestMergeFail(1),1:10)
    m = MetaArray(TestMergeFail(2),1:10)

    z = x.+y
    @test x == 1:10
    @test (x.+y) == ((1:10) .+ (1:10))
    @test (x.+y).val == 5
    @test (h.+k).val == 1
    @test_throws ErrorException h.+m
  end

  @testset "Can extract metadata and underlying array" begin
    x = meta(1:10,val=1)

    @test convert(Array,x) isa Array
    @test getcontents(x) isa AbstractRange
    @test getmeta(x).val == 1
    @test getmeta(x) isa NamedTuple
  end

  @testset "Appropriate conversion" begin
    x = meta(1:10,val=1)

    @test convert(AbstractArray{Int},x) === x
    @test convert(AbstractArray{Int,1},x) === x
    @test_throws MethodError convert(AbstractArray{String},x)
    @test_throws MethodError convert(AbstractArray{Int,2},x)
    @test Array(x) == Array(1:10)

    xplus = MetaArray((test=2,),x)
    @test getcontents(xplus) == getcontents(x)
    @test xplus.val == 1
    @test xplus.test == 2

    xplus = meta(x,test=2)
    @test getcontents(xplus) == getcontents(x)
    @test xplus.val == 1
    @test xplus.test == 2
  end

  @testset "Preoper MetaArray display" begin
    expected = "MetaArray of 1:10"
    x = meta(1:10,val=1)
    iobuf = IOBuffer()
    display(TextDisplay(iobuf), x)
    @test String(take!(iobuf)) == expected
  end

  @testset "Conversion to AbstractArray passes through" begin
    x = meta(collect(1:10),val=1)
    @test convert(AbstractArray{Int},x) === x
    @test convert(AbstractArray{Int,1},x) === x
  end

  @testset "Can compute identity of underlying array" begin
    x = meta(TestArray{Array,Int,1}(collect(1:5)),val=1)
    @test similar(x) isa MetaArray{<:TestArray}
    @test zero(x) isa MetaArray{<:TestArray}
    @test one(x) isa MetaArray{<:TestArray}
    @test all(zero(x) .== 0)
    @test all(one(x) .== 1)
  end
end


text:
using TableTraits
using IteratorInterfaceExtensions
using Test

@testset "TableTraits" begin

table_array = [(a=1,), (a=2,)]
any_table_array = Any[(a=1,), (a=2,)]
other_array = [1,2,3]
without_eltype = (i for i in table_array)

@test isiterabletable(table_array)
@test !isiterabletable(other_array)
@test isiterabletable(without_eltype)===missing
@test isiterabletable(any_table_array)===missing
@test !supports_get_columns_copy(table_array)
@test !supports_get_columns_view(table_array)
@test !supports_get_columns_copy_using_missing(table_array)
@test isiterabletable(Union{}[]) == false

end


text:
function spikeTriggeredAvg(stimulus, response, windowSize, filter=[1], cond=1)
	spikeIndices = findall(x-> x == cond, conv(response, filter))
	temp = zeros(Int(length(spikeIndices)), windowSize)
	for i in 1:1:length(spikeIndices)
		if spikeIndices[i][1] <= windowSize
			subArr = stimulus[1:spikeIndices[i][1]]
			prepend!(subArr, zeros(windowSize-spikeIndices[i][1]))
			temp[i,:] = subArr
		else
			temp[i,:] = stimulus[spikeIndices[i][1]-(windowSize-1):spikeIndices[i][1]]
		end
	end

	return mean(temp, dims=1), sqrt.(var(temp, dims=1)./size(temp)[1])
end

function spikeTriggeredAvgNDimFilter(stimulus, response, windowSize, filter=[1], cond=1)
	spikeIndices = findall(x-> x == cond, sum(conv(response, filter), dims=2))
	temp = zeros(Int(length(spikeIndices)), windowSize)
	for i in 1:1:length(spikeIndices)
		if spikeIndices[i][1] <= windowSize
			subArr = stimulus[1:spikeIndices[i][1]]
			prepend!(subArr, zeros(windowSize-spikeIndices[i][1]))
			temp[i,:] = subArr
		else
			temp[i,:] = stimulus[spikeIndices[i][1]-(windowSize-1):spikeIndices[i][1]]
		end
	end

	return mean(temp, dims=1), sqrt.(var(temp, dims=1)./size(temp)[1])
end


text:
module IndexedVectors

export IndexedVector

mutable struct IndexedVector{T, I} <: AbstractVector{T}
    _data::Vector{T}
    _map::Dict{I, Int}
    _indexer::Function
end

"""
    IndexedVector(data, indexer, copydata=false)

A constructor for IndexedVector.

# Arguments

- `data::AbstractVector`: the actual vector data
- `indexer::Function`: the indexing function which take an element of the `data` and return its key
- `copydata::Bool`: flag to specify if copy the `data`
"""

function IndexedVector(data::AbstractVector, indexer::Function, copydata::Bool=false)
    if copydata
        data = copy(data)
    end
    n = length(data)
    _map = Dict(zip(map(indexer, data), 1:n))
    return IndexedVector(data, _map, indexer)
end

Base.getindex(v::IndexedVector, i::Int) = getindex(v._data, i)
Base.size(v::IndexedVector) = size(v._data)

function Base.getindex(v::IndexedVector{T, I}, key::I)::T where {T, I}
    idx = v._map[key]
    return v._data[idx]
end

function Base.sort!(v::IndexedVector; kws...)
    sorted_data = sort!(v._data; kws...)
    return IndexedVector(sorted_data, v._indexer)
end

function Base.sort(v::IndexedVector; kws...)
    sorted_data = sort(v._data; kws...)
    return IndexedVector(sorted_data, v._indexer)
end

end


text:
#-------------------------------------------------------------------------------
# general checks and stuff

function checkazimuth(θ)
    @assert 0 <= θ < π/2 "azimuth angle θ must be ∈ [0,π/2)"
end

function checkstreams(n)
    n < 4 && @warn "careful! using nstream < 4 is likely to be inaccurate!" maxlog=1
end

#handling of inputs for temperature, molar mass, heat capacity
formprofile(P, x::AbstractVector) = AtmosphericProfile(P, x)
formprofile(::Any, x::Real) = (::Any...) -> x
formprofile(::Any, x) = x
formprofiles(P, X...) = map(x -> formprofile(P,x), X)

#------------------------------------------------------------------------------
export opticaldepth

"""
    opticaldepth(P₁, P₂, g, 𝒻T, 𝒻μ, θ, absorbers...; tol=1e-5)

Compute monochromatic [optical depths](https://en.wikipedia.org/wiki/Optical_depth#Atmospheric_sciences) (``\\tau``) between two pressure levels

# Arguments

* `P₁`: first pressure level [Pa]
* `P₂`: second pressure level [Pa]
* `g`: gravitational acceleration [m/s``^2``]
* `𝒻T`: temperature [K] as a function of pressure [Pa], `𝒻T(P)`. This may be any callable object, like [`MoistAdiabat`](@ref), for example.
* `𝒻μ`: mean molar mass as a function of temperature [K] and pressure [Pa], `𝒻μ(T,P)`
* `θ`: angle [radians] of path, must be ∈ [0,π/2), where 0 is straight up/down
* `absorbers`: at least one gas object and any number of [`CIATables`](@ref) and functions in the form σ(ν, T, P)

Returns a vector of optical depths across all wavenumbers stored in gas objects. The `tol` keyword argument adjusts integrator error tolerance.
"""
function opticaldepth(P₁::Real,
                      P₂::Real,
                      g::Real,
                      𝒻T::Q,
                      𝒻μ::R,
                      θ::Real,
                      absorbers...;
                      tol::Float64=1e-5) where {Q,R}
    #parse absorber(s)
    𝒜, ν, nν = unifyabsorbers(absorbers)
    #get pressure limits from high to low
    P₁, P₂ = max(P₁, P₂), min(P₁, P₂)
    #converted pressure coordinates
    ω₁, ω₂ = P2ω(P₁), P2ω(P₂)
    #zenith angle factor
    𝓂 = 1/cos(θ)
    #final checks
    checkpressures(𝒜, P₁, P₂)
    checkazimuth(θ)

    #spawn asynchronous parallel integrations using radau function
    tasks = Vector{Task}(undef, nν)
    for i ∈ eachindex(ν)
        tasks[i] = @spawn 𝓇depth(dτdω, ω₁, ω₂, 𝒜, i, g, 𝓂, 𝒻T, 𝒻μ, tol)
    end
    #fetch and return the results
    [fetch(task) for task ∈ tasks]
end

function opticaldepth(P::AbstractVector{<:Real},
                      g::Real,
                      T::Q,
                      μ::R,
                      θ::Real,
                      absorbers...;
                      nlobatto::Int=4) where {Q,R}
    #parse absorber(s)
    𝒜, ν, nν = unifyabsorbers(absorbers)
    #parse profiles
    𝒻T, 𝒻μ = formprofiles(P, T, μ)
    #let pressure be ascending, shouldn't matter here
    P = sort(P)
    #angle factor
    𝓂 = 1/cos(θ)
    #precalculated constant
    C = 1e-4*𝐍𝐚/g
    #pre-evaluate T and μ at intra-layer lobatto nodes
    T, μ = lobattoevaluations(P, 𝒻T, 𝒻μ, nlobatto)
    #final checks
    checkpressures(𝒜, P[end], P[1])
    checkazimuth(θ)

    #integrate in parallel
    τ = Vector{eltype(T)}(undef,nν)
    @threads for i ∈ eachindex(ν)
        τ[i] = 𝒹depth(P, T, μ, 𝒜, i, C, 𝓂, nlobatto)
    end
    return τ
end

#------------------------------------------------------------------------------
#already exported

"""
    transmittance(P₁, P₂, g, 𝒻T, 𝒻μ, θ, absorbers...; tol=1e-5)

Compute monochromatic [transmittances](https://en.wikipedia.org/wiki/Transmittance). between two pressure levels

Accepts the same arguments as [`opticaldepth`](@ref) and returns a vector of transmittances across all wavenumbers stored in gas objects.
"""
transmittance(args...; kwargs...) = exp.(-opticaldepth(args...; kwargs...))

#------------------------------------------------------------------------------
export outgoing

"""
    outgoing(Pₛ, g, 𝒻T, 𝒻μ, absorbers; Ptop=1.0, nstream=5, tol=1e-5)

Compute outgoing monochromatic radiative fluxes [W/m``^2``/cm``^{-1}``], line-by-line. Integrates the [`schwarzschild`](@ref) equation from `Pₛ` to `Ptop` at each wavenumber in the provided gas object(s) using any number of streams/angles. Total flux [W/m``^2``] can be evaluated with the [`trapz`](@ref) function. This function does not include reflected stellar radiation, which is accounted for in [`topfluxes`](@ref) and [`topimbalance`](@ref)

# Arguments
* `Pₛ`: surface pressure [Pa]
* `g`: gravitational acceleration [m/s``^2``]
* `𝒻T`: temperature [K] as a function of pressure [Pa], `𝒻T(P)`. This may be any callable object, like [`MoistAdiabat`](@ref), for example.
* `𝒻μ`: mean molar mass as a function of temperature [K] and pressure [Pa], `𝒻μ(T,P)`
* `absorbers`: at least one [gas object](gas_objects.md) and any number of [`CIATables`](@ref) and functions in the form σ(ν, T, P)

# Keywords
* `Ptop`: top of atmopshere pressure [Pa]
* `nstream`: the number of atmospheric streams (radiation angles) to calculate and integrate
* `tol`: integrator error tolerance

The keyword argument `nstream` specifies how many independent streams, or beam angles through the atmosphere, to integrate. The keyword argument `tol` is a numerical error tolerance passed to the [`radau`](https://github.com/markmbaum/ScalarRadau.jl) integrator.
"""
function outgoing(Pₛ::Real,
                  g::Real,
                  𝒻T::Q,
                  𝒻μ::R,
                  absorbers...;
                  Ptop::Real=1.0, #top of atmosphere pressure, 1 Pa by default
                  nstream::Int=5, #number of streams to use
                  tol::Real=1e-5) where {Q,R}
    #parse the absorber(s)
    𝒜, ν, nν = unifyabsorbers(absorbers)
    #convert pressure coordinates
    ω₁, ω₂ = P2ω(Pₛ, Ptop)
    #surface temperature
    Tₛ = 𝒻T(Pₛ)
    #final checks
    checkpressures(𝒜, Pₛ, Ptop)
    checkstreams(nstream)

    #integrate in parallel, dynamic schedule
    tasks = Vector{Task}(undef, nν)
    for (i,νᵢ) ∈ enumerate(ν)
        I₀ = planck(νᵢ, Tₛ)
        tasks[i] = @spawn 𝓇streams(dIdω, I₀, ω₁, ω₂, 𝒜, i, g, 𝒻T, 𝒻μ, nstream, tol)
    end
    [fetch(task) for task ∈ tasks]
end

function outgoing(P::AbstractVector{<:Real},
                  g::Real,
                  T::Q,
                  μ::R,
                  absorbers...;
                  nstream::Int=5,
                  nlobatto::Int=3) where {Q,R}
    #parse the absorber(s)
    𝒜, ν, nν = unifyabsorbers(absorbers)
    #parse profiles
    𝒻T, 𝒻μ = formprofiles(P, T, μ)
    #pressure should be descending for outgoing calculation
    P = sort(P, rev=true)
    #surface temperature [K]
    Tₛ = 𝒻T(P[1])
    #constant
    C = 1e-4*𝐍𝐚/g
    #pre-evaluate T and μ at intra-layer lobatto nodes
    T, μ = lobattoevaluations(P, 𝒻T, 𝒻μ, nlobatto)
    #final checks
    checkpressures(𝒜, P[1], P[end])
    checkstreams(nstream)

    #no need for asynchronous when fully-discretized
    olr = zeros(typeof(Tₛ), nν)
    @threads for i ∈ eachindex(ν)
        #initial irradiance
        I₀ = planck(ν[i], Tₛ)
        #monochromatic outgoing flux
        olr[i] = 𝒹streams(P, T, μ, 𝒜, i, I₀, C, nstream, nlobatto)
    end
    return olr
end

#------------------------------------------------------------------------------
export monochromaticfluxes

function monochromaticfluxes!(M⁺::AbstractMatrix,
                              M⁻::AbstractMatrix,
                              τ::AbstractMatrix,
                              core::Radau,
                              P::AbstractVector{<:Real},
                              g::Real,
                              T::Q,
                              μ::R,
                              𝒻S::U,
                              𝒻a::V,
                              absorbers...;
                              θₛ::Real=0.841)::Nothing where {Q,R,U,V}
    #parse absorber(s)
    𝒜, ν, nν = unifyabsorbers(absorbers)
    #parse profiles
    𝒻T, 𝒻μ = formprofiles(P, T, μ)
    #integration instructions
    @unpack nstream, tol = core
    #check ascending order of pressure coordinates
    @assert issorted(P) "pressure coordinates must be in ascending order (sorted)"
    #transformed coordinates
    ω, ι = P2ω.(P), P2ι.(P)
    reverse!(ω)
    #final checks
    checkpressures(𝒜, P[end], P[1])
    checkstreams(nstream)
    checkazimuth(θₛ)
    #make it clear that optical depth is not used
    τ[:] .*= NaN

    #asynchronous, parallel integrations
    tasks = Vector{Task}(undef, nν)
    for j ∈ eachindex(ν)
        Mⱼ⁺ = @view M⁺[:,j]
        Mⱼ⁻ = @view M⁻[:,j]
        tasks[j] = @spawn 𝓇monoflux!(Mⱼ⁺, Mⱼ⁻, P, ω, ι, 𝒜, j, g, 𝒻T, 𝒻μ, 𝒻S, 𝒻a, θₛ, nstream, tol)
    end
    [fetch(task) for task ∈ tasks]
    return nothing
end

function monochromaticfluxes!(M⁺::AbstractMatrix,
                              M⁻::AbstractMatrix,
                              τ::AbstractMatrix,
                              core::Discretized,
                              P::AbstractVector{<:Real},
                              g::Real,
                              T::Q,
                              μ::R,
                              𝒻S::U,
                              𝒻a::V,
                              absorbers...;
                              θₛ::Real=0.841)::Nothing where {Q,R,U,V}
    #parse absorber(s)
    𝒜, ν, _ = unifyabsorbers(absorbers)
    #parse profiles
    𝒻T, 𝒻μ = formprofiles(P, T, μ)
    #integration instructions
    @unpack nstream, nlobatto = core
    #check ascending order of pressure coordinates
    @assert issorted(P) "pressure coordinates must be in ascending order (sorted)"
    #precalculated constant
    C = 1e-4*𝐍𝐚/g
    #temperature and molar mass at lobatto nodes within layers
    T, μ = lobattoevaluations(P, 𝒻T, 𝒻μ, nlobatto)
    #pre-evaluate Planck at pressure coordinates
    B = planckevaluations(P, ν, 𝒻T)
    #final checks
    checkpressures(𝒜, P[end], P[1])
    checkstreams(nstream)
    checkazimuth(θₛ)

    #no need for asynchronous when fully discretized
    @threads for j ∈ eachindex(ν)
        τⱼ = @view τ[:,j]
        𝒹depth!(τⱼ, P, T, μ, 𝒜, j, C, nlobatto)
        Mⱼ⁺ = @view M⁺[:,j]
        Mⱼ⁻ = @view M⁻[:,j]
        Bⱼ = @view B[:,j]
        𝒹monoflux!(Mⱼ⁺, Mⱼ⁻, τⱼ, P, Bⱼ, ν[j], 𝒻S, 𝒻a, θₛ, nstream)
    end
    return nothing
end

function monochromaticfluxes(P::AbstractVector{<:Real},
                             g::Real,
                             T::Q,
                             μ::R,
                             𝒻S::U,
                             𝒻a::V,
                             absorbers...;
                             core::AbstractNumericalCore=Discretized(),
                             θₛ::Real=0.841) where {Q,R,U,V}
    #parse absorber(s)
    𝒜, _, nν = unifyabsorbers(absorbers)
    #parse profile
    𝒻T = formprofile(P, T)
    #output type
    W = typeof(𝒻T(P[1]))
    #big blocks of flux
    np = length(P)
    M⁺ = zeros(W, np, nν) #monochromatic upward fluxes
    M⁻ = zeros(W, np, nν) #monochromatic downward
    τ = zeros(W, np-1, nν)

    #main calculation of all fluxes
    monochromaticfluxes!(M⁺, M⁻, τ, core, P, g, T, μ, 𝒻S, 𝒻a, 𝒜; θₛ=θₛ)

    return M⁺, M⁻
end

#------------------------------------------------------------------------------
export fluxes, netfluxes

function fluxes(P::AbstractVector{<:Real},
                g::Real,
                T::Q,
                μ::R,
                𝒻S::U,
                𝒻a::V,
                absorbers...;
                core::AbstractNumericalCore=Discretized(),
                θₛ::Real=0.841) where {Q,R,U,V} θₛ
    #parse absorber(s)
    𝒜, ν, nν = unifyabsorbers(absorbers)
    #parse profiles
    𝒻T, 𝒻μ = formprofiles(P, T, μ)
    #coordinate count
    np = length(P)
    #output type
    W = typeof(𝒻T(P[1]))

    #big blocks of flux and depth
    M⁺ = zeros(W, np, nν) #monochromatic upward fluxes
    M⁻ = zeros(W, np, nν) #monochromatic downward
    τ = zeros(W, np-1, nν)
    #big blocks of monochromatic fluxes
    monochromaticfluxes!(M⁺, M⁻, τ, core, P, g, 𝒻T, 𝒻μ, 𝒻S, 𝒻a, 𝒜; θₛ=θₛ)
    #integrate across wavenumbers
    F⁺ = similar(M⁺, np)
    F⁻ = similar(M⁻, np)
    ∫F!(F⁺, F⁻, M⁺, M⁻, ν)
    return F⁺, F⁻
end

function netfluxes(P::AbstractVector{<:Real},
                   g::Real,
                   T::Q,
                   μ::R,
                   𝒻S::U,
                   𝒻a::V,
                   absorbers...;
                   kwargs...) where {Q,R,U,V}
    F⁺, F⁻ = fluxes(P, g, T, μ, 𝒻S, 𝒻a, absorbers...; kwargs...)
    return F⁺ .- F⁻
end

#------------------------------------------------------------------------------
export radiate!, radiate

function radiate!(F::FluxPack,
                  core::AbstractNumericalCore,
                  P::AbstractVector{<:Real},
                  g::Real,
                  T::Q,
                  μ::R,
                  𝒻S::U,
                  𝒻a::V,
                  absorbers...;
                  kwargs...)::Nothing where {Q,R,U,V}
    #parse absorber(s)
    𝒜, ν, nν = unifyabsorbers(absorbers)
    #form profile functions
    𝒻T, 𝒻μ = formprofiles(P, T, μ)
    #unpack the model arrays
    @unpack τ, M⁺, M⁻, F⁺, F⁻, Fnet = F
    #check sizes
    np = length(P)
    @assert size(F) == (np, nν) "size of FluxPack does not match number of pressure or wavenumber coordinates"
    #compute monochromatic fluxes in place
    monochromaticfluxes!(M⁺, M⁻, τ, core, P, g, 𝒻T, 𝒻μ, 𝒻S, 𝒻a, 𝒜; kwargs...)
    #the other stuff is easy
    ∫F!(F⁺, F⁻, M⁺, M⁻, ν)
    @. Fnet = F⁺ - F⁻
    #nothing to return
    return nothing
end

function radiate(P::AbstractVector{<:Real},
                 g::Real,
                 T::Q,
                 μ::R,
                 𝒻S::U,
                 𝒻a::V,
                 absorbers...;
                 core::AbstractNumericalCore=Discretized(),
                 θₛ::Real=0.841
                 )::FluxPack where {Q,R,U,V}
    #parse absorber(s)
    𝒜, ν, _ = unifyabsorbers(absorbers)
    #form profile functions
    𝒻T, 𝒻μ = formprofiles(P, T, μ)
    #allocate arrays of the proper type
    F = FluxPack(P, ν, typeof(𝒻T(P[1])))
    #compute
    radiate!(F, core, P, g, T, 𝒻μ, 𝒻S, 𝒻a, 𝒜; θₛ=θₛ)
    return F
end


text:
abstract type MTDSModel end

# Common interface
is_amortized(m::MTDSModel) = error(format("It's not known whether {:s} uses amortized inference. "*
 "Please implement `is_amortized(m)` for this type.", get_strtype_wo_params(m)))
has_x0_encoding(m::MTDSModel) = error(format("It's not known whether {:s} uses an x0 (recurrent "*
 "state) encoder. Please implement `has_x0_encoding(m)` for this type.", get_strtype_wo_params(m)))

forward(m::MTDSModel, x0, z, U; T_steps=size(U, 2)) =
    error(format("Not Implemented Yet: {:s}. Please implement a forward proc (given z).", string(typeof(m))))
encode(m::MTDSModel, Y, U; kl_coeff=1f0, stochastic=true, logstd_prior=nothing, T_enc=size(U, 2)) =
    error(format("Not Implemented Yet: {:s}. Please implement an encode proc.", string(typeof(m))))

function reconstruct(m::MTDSModel, Y, U; T_steps=size(Y, 2), T_enc=T_steps, stochastic=true)
    @argcheck size(Y,2) == size(U,2)
    (smpz, μ_z, σ_z), (smpx0, μ_x0, σ_x0) = encode(m, Y, U; T_enc=T_enc, stochastic=stochastic)
    forward(m, smpx0, smpz, U; T_steps=T_steps)
end

function forward_multiple_z(m::MTDSModel, z, u::AbstractMatrix; maxbatch=200)
    M = size(z, 2)
    ŷs = map(Iterators.partition(1:M, maxbatch)) do cbatch
        u_rep = repeat(u, outer=(1, 1, length(cbatch)));
        m(z[:, cbatch], u_rep)
    end
    cat(ŷs..., dims=3)   # problematic if length(ŷs) is too large, but there's no nice `reduce` version in Julia yet.
end
forward_multiple_z(m::MTDSModel, z, u::AbstractVector; maxbatch=200) = forward_multiple_z(m, z, unsqueeze(u,1); maxbatch=maxbatch)

do_importance_smp(m::MTDSModel, y, u, tT=80, M=200) = error(format("Not implemented yet for model type: {:s}", get_strtype_wo_params(m)))

nllh(m::MTDSModel, args...; kwargs...) = error(format("Not implemented yet for model type: {:s}", get_strtype_wo_params(m)))
nllh_of_recon_batchavg(m::MTDSModel, args...; kwargs...) = error(format("Not implemented yet for model type: {:s}", get_strtype_wo_params(m)))
elbo(m::MTDSModel, args...; kwargs...) = error(format("Not implemented yet for model type: {:s}", get_strtype_wo_params(m)))
elbo_w_kl(m::MTDSModel, args...; kwargs...) = error(format("Not implemented yet for model type: {:s}", get_strtype_wo_params(m)))

check_inputs_outputs(m::MTDSModel, Ys, Us) = error(format("Not implemented yet for model type: {:s}", get_strtype_wo_params(m)))


text:
using DrWatson
@quickactivate
using ArgParse
using GenerativeAD
using GenerativeAD.Models
using BSON
using StatsBase: fit!, predict, sample

using Flux
using MLDataPattern

s = ArgParseSettings()
@add_arg_table! s begin
   "max_seed"
		required = true
		arg_type = Int
		help = "seed"
	"dataset"
		required = true
		arg_type = String
		help = "dataset"
	"anomaly_classes"
		arg_type = Int
		default = 10
		help = "number of anomaly classes"
	"method"
		arg_type = String
		default = "leave-one-out"
		help = "method for data creation -> \"leave-one-out\" or \"leave-one-in\" "
    "contamination"
    	arg_type = Float64
    	help = "contamination rate of training data"
    	default = 0.0
end
parsed_args = parse_args(ARGS, s)
@unpack dataset, max_seed, anomaly_classes, method, contamination = parsed_args

modelname ="adVAE"

function sample_params()
	argnames = (
		:zdim,
		:nf,
		:extra_layers,
		:activation,
		:gamma,
		:lambda,
		:mx,
		:mz,
		:lr,
		:decay,
		:batch_size,
		:iters,
		:check_every,
		:patience,
		:init_seed,
	)
	par_vec = (
		2 .^(3:8), # dim of latent space and number neurons in transformer
		2 .^(2:7), # number of filters
		[0:2 ...], # extra layers
		["relu", "swish", "tanh"], # activation function for dense layers
		[0.0005, 0.001, 0.005], # γ
		[0.005, 0.01, 0.05], # λ
		[1, 1.5], #mx
		[40, 50, 60], #mz
		10f0 .^ (-4:-3), # lr
		0f0:0.1:0.5, # weight decay
		2 .^ (5:6), # batch_size
		[10000],
		[30],
		[10],
		1:Int(1e8),
	)

	return NamedTuple{argnames}(map(x->sample(x,1)[1], par_vec))
end

function fit(data, parameters)
	# define models (Generator, Discriminator)
	advae = Conv_adVAE(;parameters...)

	# define optimiser
	try
		global info, fit_t, _, _, _ = @timed fit!(advae |> gpu, data, parameters)
	catch e
		println("Error caught => $(e).")
		return (fit_t = NaN, model = nothing, history = nothing, n_parameters = NaN), []
	end

	training_info = (
		fit_t = fit_t,
		model = info[2]|>cpu,
		history = info[1], # losses through time
		npars = info[3], # number of parameters
		iters = info[4] # optim iterations of model
		)

	return training_info,
	[(x -> GenerativeAD.Models.anomaly_score(advae|>cpu, x; dims=(1,2,3), L=100)[:], merge(parameters, (L=100, )))]
	# L = samples for one x in anomaly_score computation
end

#_________________________________________________________________________________________________

try_counter = 0
max_tries = 10*max_seed
cont_string = (contamination == 0.0) ? "" : "_contamination-$contamination"
while try_counter < max_tries
	parameters = sample_params()

	for seed in 1:max_seed
		for i in 1:anomaly_classes
			savepath = datadir("experiments/images_$(method)$cont_string/$(modelname)/$(dataset)/ac=$(i)/seed=$(seed)")

			data = GenerativeAD.load_data(dataset, seed=seed, method=method, anomaly_class_ind=i)
			# computing additional parameters
			in_ch = size(data[1][1],3)
			isize = maximum([size(data[1][1],1),size(data[1][1],2)])

			isize = (isize % 16 != 0) ? isize + 16 - isize % 16 : isize
			# update parameter
			parameters = merge(parameters, (isize=isize, in_ch = in_ch, out_ch = 1))
			# here, check if a model with the same parameters was already tested
			if GenerativeAD.check_params(savepath, parameters)

				data = GenerativeAD.Models.preprocess_images(data, parameters)
				#(X_train,_), (X_val, y_val), (X_test, y_test) = data
				training_info, results = fit(data, parameters)
				# saving model separately
				if training_info.model !== nothing
					tagsave(joinpath(savepath, savename("model", parameters, "bson", digits=5)), Dict("model"=>training_info.model), safe = true)
					training_info = merge(training_info, (model = nothing,))
				end
				save_entries = merge(training_info, (modelname = modelname, seed = seed, dataset = dataset, anomaly_class = i,
					contamination=contamination))
				# now loop over all anomaly score funs
				for result in results
					GenerativeAD.experiment(result..., data, savepath; save_entries...)
				end
				global try_counter = max_tries + 1
			else
				@info "Model already present, sampling new hyperparameters..."
				global try_counter += 1
			end
		end
	end
end
(try_counter == max_tries) ? (@info "Reached $(max_tries) tries, giving up.") : nothing


text:
@doc raw"""
`SDE`: Stratonovich Stochastic Differential Equation

Defines a stochastic differential initial value problem
```math
\begin{align*}
\dq (t) &= v(t, q(t)) \, dt + B(t, q(t)) \circ dW , & q(t_{0}) &= q_{0} ,
\end{align*}
```
with drift vector field ``v``, diffusion matrix ``B``,
initial conditions ``q_{0}``, the dynamical variable ``q``
taking values in ``\mathbb{R}^{d}``, and the m-dimensional Wiener process W

### Fields

* `d`:  dimension of dynamical variable ``q`` and the vector field ``v``
* `m`:  dimension of the Wiener process
* `n`:  number of initial conditions
* `ns`: number of sample paths
* `v`:  function computing the deterministic vector field
* `B`:  function computing the d x m diffusion matrix
* `t₀`: initial time
* `q₀`: initial condition for dynamical variable ``q`` (may be a random variable itself)

The functions `v` and `B`, providing the drift vector field and diffusion matrix,
`v(t, q, v)` and `B(t, q, B; col=0)`, where `t` is the current time, `q` is the
current solution vector, and `v` and `B` are the variables which hold the result
of evaluating the vector field ``v`` and the matrix ``B`` on `t` and `q` (if col==0),
or the column col of the matrix B (if col>0).

### Example

```julia
    function v(λ, t, q, v)
        v[1] = λ*q[1]
        v[2] = λ*q[2]
    end

    function B(μ, t, q, B; col=0)
        if col==0 #whole matrix
            B[1,1] = μ*q[1]
            B[2,1] = μ*q[2]
        elseif col==1
            #just first column
        end
    end

    t₀ = 0.
    q₀ = [1., 1.]
    λ  = 2.
    μ  = 1.

    v_sde = (t, q, v) -> v(λ, t, q, v)
    B_sde = (t, q, B) -> B(μ, t, q, B)

    sde = SDE(v_sde, B_sde, t₀, q₀)
```
"""
struct SDE{dType <: Number, tType <: Number, vType <: Function, BType <: Function,
           pType <: Union{Tuple,Nothing}, N} <: Equation{dType, tType}

    d::Int
    m::Int
    n::Int
    ns::Int
    v::vType
    B::BType
    t₀::tType
    q₀::Array{dType, N}           #Initial condition: N=1 - single deterministic, N=2 - single random or multiple deterministic, N=3 - multiple random
    parameters::pType
    periodicity::Vector{dType}

    function SDE(m, n, ns, v::vType, B::BType, t₀::tType, q₀::DenseArray{dType};
                 parameters=nothing, periodicity=zeros(dType,size(q₀,1))) where {
                        dType <: Number, tType <: Number, vType <: Function, BType <: Function}

        N = ndims(q₀)
        d = size(q₀,1)

        if ns==1 && n==1
            # single sample path and single initial condition, therefore N=1
            @assert N == 1
        elseif ns>1 && n==1
            # single initial condition, but may be random (N=2) or deterministic (N=1)
            @assert N ∈ (1,2)
            if N==2
                @assert ns == size(q₀,2)
            end
        elseif ns==1 && n>1
            # multiple  deterministic initial conditions, so N=2
            @assert N == 2
            @assert n == size(q₀,2)
        elseif ns>1 && n>1
            # either multiple random initial conditions (N=3),
            # or multiple deterministic initial conditions (N=2)
            @assert N ∈ (2,3)

            if N==2
                @assert n == size(q₀,2)
            else
                @assert ns == size(q₀,2)
                @assert n == size(q₀,3)
            end
        end

        new{dType,tType,vType,BType,typeof(parameters),N}(d, m, n, ns, v, B, t₀, q₀, parameters, periodicity)
    end
end

function SDE(m::Int, ns::Int, v::Function, B::Function, t₀::Number, q₀::DenseArray{DT,1}; kwargs...) where {DT <: Number}
    # A 1D array q₀ contains a single deterministic initial condition, so n=1, but we still need to specify
    # the number of sample paths ns
    SDE(m, 1, ns, v, B, t₀, q₀; kwargs...)
end

# A 2-dimensional matrix q0 can represent a single random initial condition with ns>1 and n=1,
# or a set of deterministic initial conditions with n>1 (for which we can have both ns=1 and ns>1)
# The function below assumes q0 to represent a single random initial condition (n=1, ns=size(q₀, 2))
function SDE(m::Int, v::Function, B::Function, t₀::Number, q₀::DenseArray{DT,2}; kwargs...) where {DT <: Number}
    SDE(m, 1, size(q₀,2), v, B, t₀, q₀; kwargs...)
end

# On the other hand, the function below assumes q₀ represents multiple deterministic initial conditions
# (n=size(q₀, 2)), but these initial conditions may be run an arbitrary number ns of sample paths, so ns has to be explicitly specified
function SDE(m::Int, ns::Int, v::Function, B::Function, t₀::Number, q₀::DenseArray{DT,2}; kwargs...) where {DT <: Number}
    SDE(m, size(q₀,2), ns, v, B, t₀, q₀; kwargs...)
end

function SDE(m::Int, v::Function, B::Function, t₀::Number, q₀::DenseArray{DT,3}; kwargs...) where {DT <: Number}
    # A 3D array q₀ contains multiple random initial condition, so n=size(q₀,3) and ns=size(q₀,2)
    SDE(m, size(q₀,3), size(q₀,2), v, B, t₀, q₀; kwargs...)
end

function SDE(m::Int, ns::Int, v::Function, B::Function, q₀::DenseArray{DT}; kwargs...) where {DT}
    SDE(m, ns, v, B, zero(DT), q₀; kwargs...)
end

function SDE(m::Int, v::Function, B::Function, q₀::DenseArray{DT}; kwargs...) where {DT}
    SDE(m, v, B, zero(DT), q₀; kwargs...)
end

Base.hash(sde::SDE, h::UInt) = hash(sde.d, hash(sde.m, hash(sde.n, hash(sde.ns,
                               hash(sde.v, hash(sde.B, hash(sde.t₀, hash(sde.q₀,
                               hash(sde.parameters, hash(sde.periodicity, h))))))))))

Base.:(==)(sde1::SDE, sde2::SDE) = (
                                sde1.d == sde2.d
                             && sde1.m == sde2.m
                             && sde1.n == sde2.n
                             && sde1.ns == sde2.ns
                             && sde1.v == sde2.v
                             && sde1.B == sde2.B
                             && sde1.t₀ == sde2.t₀
                             && sde1.q₀ == sde2.q₀
                             && sde1.parameters == sde2.parameters
                             && sde1.periodicity == sde2.periodicity)

function Base.similar(sde::SDE, q₀::DenseArray)
    similar(sde, sde.t₀, q₀)
end

function Base.similar(sde::SDE, t₀::TT, q₀::DenseArray{DT}) where {DT <: Number, TT <: Number}
    @assert sde.d == size(q₀,1)
    SDE(sde.m, sde.ns, sde.v, sde.B, t₀, q₀, periodicity=sde.periodicity)
end

function Base.similar(sde::SDE, q₀::DenseArray, ns::Int)
    similar(sde, sde.t₀, q₀, ns)
end

function Base.similar(sde::SDE, t₀::TT, q₀::DenseArray{DT}, ns::Int) where {DT <: Number, TT <: Number}
    @assert sde.d == size(q₀,1)
    SDE(sde.m, ns, sde.v, sde.B, t₀, q₀, periodicity=sde.periodicity)
end

@inline Base.ndims(sde::SDE) = sde.d

@inline periodicity(equation::SDE) = equation.periodicity


text:
function expm1()

    T11 = 1
    T12 = 0.4
    T21 = 0.08
    T22 = 0.02
    τ12 = 0.05
    τ21 = 0.1

    expA = BlochMcConnellMatrix(2)
    A = BlochSim.BlochMcConnellDynamicsMatrix(2)
    A.A[1].R2 = -1 / T21 - 1 / τ12
    A.A[1].Δω = 2π
    A.A[1].R1 = -1 / T11 - 1 / τ12
    A.A[2].R2 = -1 / T22 - 1 / τ21
    A.A[2].Δω = 2π
    A.A[2].R1 = -1 / T12 - 1 / τ21
    A.E[1].r = 1 / τ12
    A.E[2].r = 1 / τ21

    correct = exp(Matrix(A))
    BlochSim.expm!(expA, A)

    return Matrix(expA) ≈ correct

end

function dfrexp1()

    f = x -> 2x * BlochSim.frexp1(x)
    x = LinRange(-5, 5, 1000)
    df_correct = 2BlochSim.frexp1.(x) .+ 2x .* BlochSim.dfrexp1.(x)
    df_forwarddiff = ForwardDiff.derivative.(f, x)

    return df_forwarddiff ≈ df_correct

end

function dfrexp2()

    f = x -> 2x * BlochSim.frexp2(x)
    x = LinRange(-5, 5, 1000)
    df_correct = 2BlochSim.frexp2.(x) .+ 2x .* BlochSim.dfrexp2.(x)
    df_forwarddiff = ForwardDiff.derivative.(f, x)

    return df_forwarddiff ≈ df_correct

end

function autodiff1()

    Ae = ExcitationMatrix()
    rf = InstantaneousRF(π/2)
    grad = ForwardDiff.gradient([1000.0, 100.0]) do x
        T1, T2 = x
        T = typeof(T1)
        Af = FreePrecessionMatrix{T}()
        Bf = Magnetization{T}()
        BtoM = Magnetization{T}()
        s = Spin(1, T1, T2, 10)
        excite!(Ae, s, rf)
        applydynamics!(s, BtoM, Ae)
        freeprecess!(Af, Bf, s, 10)
        applydynamics!(s, BtoM, Af, Bf)
        abs.(signal(s))
    end
    correct = [0.0, 0.0009048374180359595]

    return grad ≈ correct

end

function autodiff2()

    Ae = ExcitationMatrix()
    rf = InstantaneousRF(π/2)
    grad = ForwardDiff.gradient([1000.0, 100.0]) do x
        T1, T2 = x
        T = typeof(T1)
        Af = BlochMcConnellMatrix{T}(2)
        Bf = MagnetizationMC{T}(2)
        BtoM = MagnetizationMC{T}(2)
        s = SpinMC(1, [0.15, 0.85], [400, T1], [20, T2], [25, 10], [Inf, Inf])
        excite!(Ae, s, rf)
        applydynamics!(s, BtoM, Ae)
        freeprecess!(Af, Bf, s, 10)
        applydynamics!(s, BtoM, Af, Bf)
        abs.(signal(s))
    end
    correct = [0.0, 0.0007660512555728833]

    return grad ≈ correct

end

function autodiff3()

    Ae = ExcitationMatrix()
    rf = InstantaneousRF(π/2)
    grad = ForwardDiff.gradient([1000.0, 100.0]) do x
        T1, T2 = x
        T = typeof(T1)
        Af = BlochMcConnellMatrix{T}(2)
        Bf = MagnetizationMC{T}(2)
        BtoM = MagnetizationMC{T}(2)
        s = SpinMC(1, [0.15, 0.85], [400, T1], [20, T2], [25, 10], [Inf, Inf])
        excite!(Ae, s, rf)
        applydynamics!(s, BtoM, Ae)
        freeprecess!(Af, Bf, s, 2)
        applydynamics!(s, BtoM, Af, Bf)
        abs.(signal(s))
    end
    correct = [0.0, 0.00016657611260161996]

    return grad ≈ correct

end

function autodiff4()

    Ae = ExcitationMatrix()
    rf = InstantaneousRF(π/2)
    grad = ForwardDiff.gradient([1000.0, 100.0]) do x
        T1, T2 = x
        T = typeof(T1)
        Af = BlochMcConnellMatrix{T}(2)
        Bf = MagnetizationMC{T}(2)
        BtoM = MagnetizationMC{T}(2)
        s = SpinMC(1, [0.15, 0.85], [400, T1], [20, T2], [25, 10], [Inf, Inf])
        excite!(Ae, s, rf)
        applydynamics!(s, BtoM, Ae)
        freeprecess!(Af, Bf, s, 1)
        applydynamics!(s, BtoM, Af, Bf)
        abs.(signal(s))
    end
    correct = [0.0, 8.414639502122038e-5]

    return grad ≈ correct

end

function autodiff5()

    Ae = ExcitationMatrix()
    rf = InstantaneousRF(π/2)
    grad = ForwardDiff.gradient([1000.0, 100.0]) do x
        T1, T2 = x
        T = typeof(T1)
        Af = BlochMcConnellMatrix{T}(2)
        Bf = MagnetizationMC{T}(2)
        BtoM = MagnetizationMC{T}(2)
        s = SpinMC(1, [0.15, 0.85], [400, T1], [20, T2], [25, 10], [Inf, Inf])
        excite!(Ae, s, rf)
        applydynamics!(s, BtoM, Ae)
        freeprecess!(Af, Bf, s, 0.1)
        applydynamics!(s, BtoM, Af, Bf)
        abs.(signal(s))
    end
    correct = [0.0, 8.491495820718459e-6]

    return grad ≈ correct

end

function autodiff6()

    Ae = ExcitationMatrix()
    rf = InstantaneousRF(π/2)
    grad = ForwardDiff.gradient([1000.0, 100.0]) do x
        T1, T2 = x
        T = typeof(T1)
        Af = BlochMcConnellMatrix{T}(2)
        Bf = MagnetizationMC{T}(2)
        BtoM = MagnetizationMC{T}(2)
        s = SpinMC(1, [0.15, 0.85], [400, T1], [20, T2], [25, 10], [Inf, Inf])
        excite!(Ae, s, rf)
        applydynamics!(s, BtoM, Ae)
        freeprecess!(Af, Bf, s, 0.01)
        applydynamics!(s, BtoM, Af, Bf)
        abs.(signal(s))
    end
    correct = [0.0, 8.499149957624547e-7]

    return grad ≈ correct

end

@testset "Matrix Exponential" begin

    @testset "expm Accuracy" begin

        @test expm1()

    end

    @testset "frexp Derivatives" begin

        @test dfrexp1()
        @test dfrexp2()

    end

    @testset "AutoDiff" begin

        @test autodiff1()
        @test autodiff2()
        @test autodiff3()
        @test autodiff4()
        @test autodiff5()
        @test autodiff6()

    end

end


text:
using Documenter
import OptimizingIR

makedocs(
    sitename = "OptimizingIR.jl",
    modules = [ OptimizingIR ],
    pages = [ "Home" => "index.md",
              "API Reference" => "api.md" ]
)

deploydocs(
    repo = "github.com/felipenoris/OptimizingIR.jl.git",
    target = "build",
)


text:
include("utils.jl")
#include("fileio.jl")
include("packer.jl")
include("shape/truncated_tet.jl")
include("args.jl")

options = parse_args(ARGS, s)

if options["s"] != nothing
    Random.seed!(options['s'])
end

d_eq, d₀ = parse_dist(options["d"])
r, r₀ = parse_dist(options["r"])
k, k₀ = parse_dist(options["k"])
generator = TTGenerator(r, k, d_eq)

vol_frac_goal = options["vol_frac_goal"]
L = 2. #options["L"]
m = 5 #options["m"]
M = round(Int16, m * L / d₀)
nr_tries = 5 #options["nr_tries"]
Δ = M÷Int8(2) #options["delta"]
spacing = L / M

inclusions = prepare_inclusions(vol_frac_goal, L, TruncatedTriangle, generator)

# Sort triangles w.r.t. volume typically leads to better packing:
sort!(inclusions, by=volume, rev=true)

println("Prepared $(length(inclusions)) inclusions")

#tt = TruncatedTriangle(SA[0., 0., 0.], SA[1. 0. 0.; 0. 1. 0.; 0. 0. 1.], r₀, k₀, d₀)
#lower, upper = bounding_box(tt)
#test = voxelize(inclusions[1], L, M)

grain_ids, overlaps, inclusions_voxels = pack_inclusions!(inclusions, M, L, nr_tries, Δ)

m_phase = phases(grain_ids)
m_euler_angles = euler_angles(grain_ids, inclusions)

make_mcp_bound!(grain_ids, inclusions_voxels, overlaps, voxel_indices, steps::Integer, kBT::Real)

write_dream3d(options["f"], SA[L, L, L], grain_ids, phases, euler_angles)


text:
"""
    `count_case(cases::T) where T <:AbstractVector{Bool}`
    starting from 1 increase the count value after an event signalled on boolvec
    is true
"""
function count_cases(cases::T) where T <:AbstractVector{Bool}
    shifted_cases = vcat([Bool(0)], cases[1:end-1])
    accumulate(+, shifted_cases; init = 1)
end

"""
    `count_bouts!(pokedf::DataFrames.AbstractDataFrame)`
    using a pokes dataframe count reward bouts within each session
"""
function count_bouts!(pokedf::DataFrames.AbstractDataFrame)
    if ((:LastPoke in propertynames(pokedf)) && (:Reward in propertynames(pokedf)))
        return transform!(groupby(pokedf,[:ExpDay,:MouseID]), [:LastPoke, :Reward] => ((l,r) -> count_cases(l .| r)) => :Bout)
    else
        error("Not a pokes dataframe")
    end
end

function process_bouts(pokedf::AbstractDataFrame)
    if !(:Bout in propertynames(pokedf))
        error("No bout information")
    else
        # [:ExpDay, :MouseID, :Day, :Trial, :Poke, :PokeIn, :PokeOut, :Side, :Reward,
        # :ROI_In, :ROI_Out, :Protocol, :Stim, :Poke_within_Trial, :Pre_Interpoke, :Post_Interpoke,
        # :LastPoke, :Poke_Hierarchy, :ReverseTrial, :Correct, :PokeDuration, :Stim_Day, :Phase, :Group,
        # :Treatment, :Injection, :ExpType, :ProtocolDay, :Next_Prew, :Pnextrew, :InstRewRate, :ElapsedTime,
        # :Cumulative_Reward, :AverageRewRate, :TimeFromLeaving, :CumRewTrial, :Bout]

        vals_per_bout = [:Side, :Trial,:Protocol, :Stim,:ROI_In, :ROI_Out]
        vals_per_day = [:Day, :Phase, :Group, :Treatment, :Injection, :ExpType, :ProtocolDay, :Stim_Day];
        bout_table = combine(groupby(pokedf, :Bout)) do dd
            dt = DataFrame(
            Omissions_plus_one = size(dd,1),
            Leave = dd[end,:LastPoke],
            End_Reward = dd[end,:Reward],
            BoutIn = dd[1,:PokeIn],
            BoutOut = dd[end,:PokeOut],
            Bout_duration = (dd[end,:PokeOut]-dd[1,:PokeIn]),
            Next_Prew = dd[end,:Next_Prew],
            Pnextrew = dd[end,:Pnextrew],
            InstRewRate = dd[end,:InstRewRate],
            AverageRewRate = dd[end,:AverageRewRate]
            )

            for v in vals_per_bout
                if v in propertynames(dd)
                    dt[!,v] .= dd[1, v]
                end
            end

            return dt
        end

        for v in vals_per_day
            if v in propertynames(pokedf)
                bout_table[!,v] .= pokedf[1, v]
            end
        end

        return bout_table
    end
end


text:
using DataFrames
using LaTeXStrings

df = DataFrame(A = 'x':'z', B = ["α/β", 1//2, 8])

@test latexify(df, latex=true) ==
Markdown.md"|   A |                      B |
| ---:| ----------------------:|
| $x$ | $\frac{\alpha}{\beta}$ |
| $y$ |          $\frac{1}{2}$ |
| $z$ |                    $8$ |
"

@test_broken latexify(df; env=:array) == replace(
L"\begin{equation}
\left[
\begin{array}{cc}
A & B \\
x & \frac{\alpha}{\beta} \\
y & \frac{1}{2} \\
z & 8 \\
\end{array}
\right]
\end{equation}
", "\r\n"=>"\n")

@test latexify(df; env=:table, latex=true) == replace(
raw"\begin{tabular}{cc}
$A$ & $B$\\
$x$ & $\frac{\alpha}{\beta}$\\
$y$ & $\frac{1}{2}$\\
$z$ & $8$\\
\end{tabular}
", "\r\n"=>"\n")

@test latexify(df; env=:table, latex=false) == replace(
raw"\begin{tabular}{cc}
A & B\\
x & α/β\\
y & 1//2\\
z & 8\\
\end{tabular}
", "\r\n"=>"\n")


text:
using MPI
using StaticArrays
using Logging
using Printf
using LinearAlgebra
using GPUifyLoops
using Test
using CLIMA
using CLIMA.Mesh.Topologies
using CLIMA.Mesh.Grids
using CLIMA.MPIStateArrays
using CLIMA.DGmethods
using CLIMA.DGmethods.NumericalFluxes
using CLIMA.ODESolvers
using CLIMA.GenericCallbacks
using CLIMA.Atmos
using CLIMA.VariableTemplates
using CLIMA.MoistThermodynamics
using CLIMA.PlanetParameters

function run_test1(mpicomm, dim, Ne, N, FT, ArrayType)
  warpfun = (ξ1, ξ2, ξ3) -> begin
    x1 = ξ1 + (ξ1 - 1/2) * cos(2 * π * ξ2 * ξ3) / 4
    x2 = ξ2  + (ξ2 - 1/2) * cos(2 * π * ξ2 * ξ3) / 4
    x3 = ξ3 + ξ1 / 4 + sin(2 * π * ξ1) / 16
    return (x1, x2, x3)
  end

  brickrange = (range(FT(0); length=Ne+1, stop=1),
                range(FT(0); length=Ne+1, stop=1),
                range(FT(0); length=2, stop=1))
  topl = StackedBrickTopology(mpicomm, brickrange,
                              periodicity = (false, false, false))

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                          meshwarp = warpfun,
                                         )

  nrealelem = length(topl.realelems)
  Nq = N + 1
  Nqk = dimensionality(grid) == 2 ? 1 : Nq
  vgeo = grid.vgeo
  host_array = Array ∈ typeof(vgeo).parameters
  localvgeo = host_array ? vgeo : Array(vgeo)

  S = zeros(Nqk)
  S1 = zeros(Nqk)

  for e in 1:nrealelem
    for k in 1:Nqk
      for j in 1:Nq
        for i in 1:Nq
          ijk = i + Nq * ((j-1)+ Nq * (k-1))
          S[k] += localvgeo[ijk,grid.x1id,e] * localvgeo[ijk,grid.MHid,e]
          S1[k] += localvgeo[ijk,grid.MHid,e]
        end
      end
    end
  end

  Stot = zeros(Nqk)
  S1tot = zeros(Nqk)
  Err = 0.0

  for k in 1:Nqk
    Stot[k] = MPI.Reduce(S[k], +, 0, mpicomm)
    S1tot[k] = MPI.Reduce(S1[k], +, 0, mpicomm)
    Err += (0.5 - Stot[k] / S1tot[k])^2
  end
  Err = sqrt(Err / Nqk)
  @test Err < 2e-15
end

function run_test2(mpicomm, dim, Ne, N, FT, ArrayType)
  warpfun = (ξ1, ξ2, ξ3) -> begin
    x1 = sin(2 * π * ξ3)/16 + ξ1
    x2 = ξ2  + (ξ2 - 1/2) * cos(2 * π * ξ2 * ξ3) / 4
    x3 = ξ3 + sin(2π * (ξ1))/20
    return (x1, x2, x3)
  end

  brickrange = (range(FT(0); length=Ne+1, stop=1),
                range(FT(0); length=Ne+1, stop=1),
                range(FT(0); length=2, stop=1))
  topl = StackedBrickTopology(mpicomm, brickrange,
                              periodicity = (false, false, false))

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                          meshwarp = warpfun,
                                         )

  nrealelem = length(topl.realelems)
  Nq = N + 1
  Nqk = dimensionality(grid) == 2 ? 1 : Nq
  vgeo = grid.vgeo
  host_array = Array ∈ typeof(vgeo).parameters
  localvgeo = host_array ? vgeo : Array(vgeo)

  S = zeros(Nqk)
  S1 = zeros(Nqk)
  K = zeros(Nqk)

  for e in 1:nrealelem
    for k in 1:Nqk
      for j in 1:Nq
        for i in 1:Nq
          ijk = i + Nq * ((j-1)+ Nq * (k-1))
          S[k] += localvgeo[ijk,grid.x1id,e] * localvgeo[ijk,grid.MHid,e]
          S1[k] += localvgeo[ijk,grid.MHid,e]
          K[k] = localvgeo[ijk,grid.x3id,e]
        end
      end
    end
  end

  Stot = zeros(Nqk)
  S1tot = zeros(Nqk)
  Err = 0.0

  for k in 1:Nqk
    Stot[k] = MPI.Reduce(S[k], +, 0, mpicomm)
    S1tot[k] = MPI.Reduce(S1[k], +, 0, mpicomm)
    Err += (0.5 + sin(2 * π * K[k]) / 16 - Stot[k] / S1tot[k])^2
  end
  Err = sqrt(Err / Nqk)
  @test 2e-15 > Err
end

function run_test3(mpicomm, dim, Ne, N, FT, ArrayType)
  base_Nhorz = 4
  base_Nvert = 2
  Rinner = 1 // 2
  Router = 1
  expected_result = [-4.5894269717905445e-8  -1.1473566985387151e-8 ;
		                 -2.0621904184281448e-10 -5.155431637149377e-11 ;
	                   -8.72191208145523e-13   -2.1715962361668062e-13]
  for l in 1:3
    Nhorz = 2^(l-1) * base_Nhorz
    Nvert = 2^(l-1) * base_Nvert
    Rrange = grid1d(FT(Rinner), FT(Router); nelem=Nvert)
    topl = StackedCubedSphereTopology(mpicomm, Nhorz, Rrange)
    grid = DiscontinuousSpectralElementGrid(topl,
                                            FloatType = FT,
                                            DeviceArray = ArrayType,
                                            polynomialorder = N,
                                            meshwarp = Topologies.cubedshellwarp,
                                           )

    nrealelem = length(topl.realelems)
    Nq = N + 1
    Nqk = dimensionality(grid) == 2 ? 1 : Nq
    vgeo = grid.vgeo
    host_array = Array ∈ typeof(vgeo).parameters
    localvgeo = host_array ? vgeo : Array(vgeo)

    topology = grid.topology
    nvertelem = topology.stacksize
    nhorzelem = div(nrealelem, nvertelem)

    Surfout = 0
    Surfin = 0

    for ev in 1:nvertelem
      for eh in 1:nhorzelem
        e = ev + (eh - 1) * nvertelem
        for i in 1:Nq
          for j in 1:Nq
            for k in 1:Nqk
              ijk = i + Nq * ((j-1) + Nqk * (k-1))
              if (k == Nqk && ev == nvertelem)
                Surfout += localvgeo[ijk,grid.MHid,e]
              end
              if (k == 1 && ev == 1)
                Surfin += localvgeo[ijk,grid.MHid,e]
              end
            end
          end
        end
      end
    end
    Surfouttot = MPI.Reduce(Surfout, +, 0, MPI.COMM_WORLD)
    Surfintot = MPI.Reduce(Surfin, +, 0, MPI.COMM_WORLD)
    @test (4 * π * Router^2 - Surfouttot) ≈ expected_result[l,1] rtol=1e-3 atol=eps(FT) * 4 * π * Router^2
    @test (4 * π * Rinner^2 - Surfintot)  ≈ expected_result[l,2] rtol=1e-3 atol=eps(FT) * 4 * π * Rinner^2
  end
end

# Test for 2D integral
function run_test4(mpicomm, dim, Ne, N, FT, ArrayType)
  brickrange = ntuple(j->range(FT(0); length=Ne+1, stop=1), 2)
  topl = StackedBrickTopology(mpicomm, brickrange,
                              periodicity=ntuple(j->true, 2))

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
                                         )
  nrealelem = length(topl.realelems)
  Nq = N + 1
  vgeo = grid.vgeo
  host_array = Array ∈ typeof(vgeo).parameters
  localvgeo = host_array ? vgeo : Array(vgeo)

  S = zeros(Nq)

  for e in 1:nrealelem
    for i in 1:Nq
      for j in 1:Nq
        ij = i + Nq * (j - 1)
        S[j] += localvgeo[ij,grid.x1id,e] * localvgeo[ij,grid.MHid,e]
      end
    end
  end

  Err = 0
  for j in 1:Nq
    Err += (S[j] - 0.5)^2
  end

  Err = sqrt(Err / Nq)
  @test Err <= 1e-15
end

function run_test5(mpicomm, dim, Ne, N, FT, ArrayType)
  warpfun = (ξ1, ξ2, ξ3) -> begin
    x1 = cos(π * ξ2)/16 + abs(ξ1)
    x2 = ξ2
    x3 = ξ3
    return (x1, x2, x3)
  end

  brickrange = ntuple(j->range(FT(0); length=Ne+1, stop=1), 2)
  topl = StackedBrickTopology(mpicomm, brickrange,
                              periodicity=ntuple(j->true, 2))

  grid = DiscontinuousSpectralElementGrid(topl,
                                          FloatType = FT,
                                          DeviceArray = ArrayType,
                                          polynomialorder = N,
					  meshwarp = warpfun,
                                         )

  nrealelem = length(topl.realelems)
  Nq = N + 1
  vgeo = grid.vgeo
  host_array = Array ∈ typeof(vgeo).parameters
  localvgeo = host_array ? vgeo : Array(vgeo)

  S = zeros(Nq)
  J = zeros(Nq)

  for e in 1:nrealelem
    for i in 1:Nq
      for j in 1:Nq
        ij = i + Nq * (j - 1)
        S[j] += localvgeo[ij,grid.x1id,e] * localvgeo[ij,grid.MHid,e]
        J[j] = localvgeo[ij,grid.x2id,e]
      end
    end
  end

  Err = 0
  for j in 1:Nq
    Err += (S[j] - 0.5 - cos(π * J[j]) / 16)^2
  end

  Err = sqrt(Err / Nq)
  @test Err < 1e-15
end

let
  CLIMA.init()
  ArrayType = CLIMA.array_type()

  mpicomm = MPI.COMM_WORLD
  ll = uppercase(get(ENV, "JULIA_LOG_LEVEL", "INFO"))
  loglevel = ll == "DEBUG" ? Logging.Debug :
  ll == "WARN"  ? Logging.Warn  :
  ll == "ERROR" ? Logging.Error : Logging.Info
  logger_stream = MPI.Comm_rank(mpicomm) == 0 ? stderr : devnull
  global_logger(ConsoleLogger(logger_stream, loglevel))

  FT = Float64
  dim = 3
  Ne = 1
  polynomialorder = 4

  @info (ArrayType, FT, dim)

  @testset "horizontal_integral" begin
    run_test1(mpicomm, dim, Ne, polynomialorder, FT, ArrayType)
    run_test2(mpicomm, dim, Ne, polynomialorder, FT, ArrayType)
    run_test3(mpicomm, dim, Ne, polynomialorder, FT, ArrayType)
    run_test4(mpicomm, dim, Ne, polynomialorder, FT, ArrayType)
    run_test5(mpicomm, dim, Ne, polynomialorder, FT, ArrayType)
  end
end

max_stars_repo_name: favba/CLIMA

max_issues_repo_name: favba/CLIMA

max_forks_repo_name: favba/CLIMA


text:
"""
# Description
Constructor for the isotropic Squared Exponential kernel (covariance)

k(x,x') = σ²exp(-(x-x')ᵀ(x-x')/2ℓ²)
# Arguments:
* `ll::Float64`: Log of the length scale ℓ
* `lσ::Float64`: Log of the signal standard deviation σ
"""
type SEIso <: Isotropic
    ℓ2::Float64      # Length scale
    σ2::Float64      # Signal std
    SEIso(ll::Float64, lσ::Float64) = new(exp(2*ll), exp(2*lσ))
end

function set_params!(se::SEIso, hyp::Vector{Float64})
    length(hyp) == 2 || throw(ArgumentError("Squared exponential only has two parameters"))
    se.ℓ2, se.σ2 = exp(2.0*hyp)
end

get_params(se::SEIso) = Float64[log(se.ℓ2)/2.0, log(se.σ2)/2.0]
get_param_names(::SEIso) = [:ll, :lσ]
num_params(se::SEIso) = 2

metric(se::SEIso) = SqEuclidean()
cov(se::SEIso, r::Float64) = se.σ2*exp(-0.5*r/se.ℓ2)

function grad_kern(se::SEIso, x::Vector{Float64}, y::Vector{Float64})
    r = distance(se, x, y)
    exp_r = exp(-0.5*r/se.ℓ2)
    g1 = se.σ2*r/se.ℓ2*exp_r # dK_d(log ℓ)
    g2 = 2.0*se.σ2*exp_r     # dK_d(log σ)
    return [g1,g2]
end

function grad_stack!(stack::AbstractArray, se::SEIso, X::Matrix{Float64}, data::IsotropicData)
    nobsv = size(X, 2)
    R = distance(se, X, data)
    exp_R = exp(-0.5*R/se.ℓ2)
    for i in 1:nobsv, j in 1:i
        @inbounds stack[i,j,1] = se.σ2*R[i,j]/se.ℓ2*exp_R[i,j] # dK_dℓ
        @inbounds stack[j,i,1] = stack[i,j,1]
        @inbounds stack[i,j,2] = 2.0*se.σ2*exp_R[i,j]
        @inbounds stack[j,i,2] = stack[i,j,2]
    end

    return stack
end


text:
__precompile__()

module Bio

import BioSequences
import GenomicFeatures
import BioAlignments
import BioStructures
import GeneticVariation
import BioServices
import BioTools

# TODO
#import Phylogenies

const Seq = BioSequences
const Intervals = GenomicFeatures
const Align = BioAlignments
const Structure = BioStructures
const Var = GeneticVariation
const Services = BioServices
const Tools = BioTools

# TODO
#const Phylo = Phylogenies

function __init__()
    @warn "This package has been depreceated, and should not be used for new projects. Please see Bio.jl's repository README for more information."
end

end  # module Bio