text:
module PyPlot

using PyCall
import PyCall: PyObject, pygui
import Base: convert, ==, isequal, hash, writemime, getindex, setindex!, haskey, keys, show, mimewritable
export Figure, plt, matplotlib, pygui, withfig

###########################################################################
# file formats supported by Agg backend, from MIME types
const aggformats = [ "application/eps" => "eps", "image/eps" => "eps",
                     "application/pdf" => "pdf",
                     "image/png" => "png",
                     "application/postscript" => "ps",
                     "image/svg+xml" => "svg" ]

function isdisplayok()
    for mime in keys(aggformats)
        if displayable(mime)
            return true
        end
    end
    false
end

###########################################################################
# We allow the user to turn on or off the Python gui interactively via
# pygui(true/false).  This is done by loading pyplot with a GUI backend
# if possible, then switching to a Julia-display backend (if available),
# hooking into pyplot via a monkey-patched draw_if_interactive.

const isjulia_display = Bool[isdisplayok()]
const matplotlib = pyimport("matplotlib")
const version = try
    convert(VersionNumber, matplotlib[:__version__])
catch
    v"0.0" # fallback
end

pymodule_exists(s::String) = try
    pyimport(s)
    true
catch
    false
end

const (backend, gui) = begin
    const gui2matplotlib = [ :wx=>"WXAgg", :gtk=>"GTKAgg", :qt=>"Qt4Agg" ]
    try
        # We will get an exception when we import pyplot below
        # (on Unix) if an X server is not available, even though
        # pygui_works and matplotlib.use(backend) succeed, at
        # which point it will be too late to switch backends.  So,
        # throw exception (drop to catch block below) if DISPLAY is not set.
        # [Might be more reliable to test success(`xdpyinfo`), but only
        #  if xdpyinfo is installed.]
        @unix_only (@osx ? nothing : ENV["DISPLAY"])

        local gui::Symbol = :none
        if PyCall.gui == :default
            # try to ensure that GUI both exists and has a matplotlib backend
            for g in (@linux? (:gtk, :qt, :wx) : (:qt, :wx, :gtk))
                if PyCall.pygui_works(g)
                    # must call matplotlib.use *before* loading backends module
                    matplotlib[:use](gui2matplotlib[g])
                    if g == :qt && !PyCall.pyexists("PyQt4")
                        PyDict(matplotlib["rcParams"])["backend.qt4"]="PySide"
                    end
                    if pymodule_exists(string("matplotlib.backends.backend_",
                                              lowercase(gui2matplotlib[g])))
                        gui = g
                        break
                    end
                end
            end
            if gui == :none
                error("no gui found") # go to catch clause below
            end
        else
            gui = pygui()
            matplotlib[:use](gui2matplotlib[gui])
        end
        if !isjulia_display[1]
            pygui_start(gui)
        end
        matplotlib[:interactive](true)
        (gui2matplotlib[gui], gui)
    catch
        if !isjulia_display[1]
            warn("No working GUI backend found for matplotlib.")
            isjulia_display[1] = true
        end
        pygui(:default)
        matplotlib[:use]("Agg") # GUI not available
        matplotlib[:interactive](isdisplayok())
        ("Agg", :none)
    end
end

const pltm = pyimport("matplotlib.pyplot") # raw Python module

function pygui(b::Bool)
    if !b != isjulia_display[1]
        if backend != "Agg"
            pltm[:switch_backend](b ? backend : "Agg")
            monkeypatch()
            if b
                pygui_start(gui) # make sure event loop is started
            end
        elseif b
            error("No working GUI backend found for matplotlib.")
        end
        isjulia_display[1] = !b
    end
    return b
end

###########################################################################
# Wrapper around matplotlib Figure, supporting graphics I/O and pretty display

type Figure
    o::PyObject
end

PyObject(f::Figure) = f.o
convert(::Type{Figure}, o::PyObject) = Figure(o)
==(f::Figure, g::Figure) = f.o == g.o
isequal(f::Figure, g::Figure) = f.o == g.o # Julia 0.2 compatibility
==(f::Figure, g::PyObject) = f.o == g
==(f::PyObject, g::Figure) = f == g.o
hash(f::Figure) = hash(f.o)

getindex(f::Figure, x) = getindex(f.o, x)
setindex!(f::Figure, v, x) = setindex!(f.o, v, x)
haskey(f::Figure, x) = haskey(f.o, x)
keys(f::Figure) = keys(f.o)

pytype_mapping(pltm["Figure"], Figure)

for (mime,fmt) in aggformats
    @eval function writemime(io::IO, m::MIME{symbol($mime)}, f::Figure)
        if !haskey(pycall(f.o["canvas"]["get_supported_filetypes"], PyDict),
                   $fmt)
            throw(MethodError(writemime, (io, m, f)))
        end
        f.o["canvas"][:print_figure](io, format=$fmt, bbox_inches="tight")
    end
    if fmt != "svg"
        @eval mimewritable(::MIME{symbol($mime)}, f::Figure) = !isempty(f) && haskey(pycall(f.o["canvas"]["get_supported_filetypes"], PyDict), $fmt)
    end
end

# disable SVG output by default, since displaying large SVGs (large datasets)
# in IJulia is slow, and browser SVG display is buggy.  (Similar to IPython.)
const SVG = [false]
mimewritable(::MIME"image/svg+xml", f::Figure) = SVG[1] && !isempty(f) && haskey(pycall(f.o["canvas"]["get_supported_filetypes"], PyDict), "svg")
svg() = SVG[1]
svg(b::Bool) = (SVG[1] = b)

###########################################################################
# Monkey-patch pylab to call redisplay after each drawing command
# (which calls draw_if_interactive) for Julia displays.

const Gcf = pyimport("matplotlib._pylab_helpers")["Gcf"]
const drew_something = [false]
const orig_draw = pltm["draw_if_interactive"]

Base.isempty(f::Figure) = isempty(pycall(f["get_axes"], PyVector))

# monkey-patch draw_if_interactive to queue the figure for drawing in IJulia
function draw_if_interactive()
    if isjulia_display[1]
        if pycall(matplotlib["is_interactive"], Bool)
            manager = Gcf[:get_active]()
            if manager != nothing
                fig = Figure(manager["canvas"]["figure"])
                redisplay(fig)
                drew_something[1] = true
            end
        end
    else
        pycall(orig_draw, PyObject)
    end
    nothing
end

# The logic of display/redisplay/close is a bit complicated.  In
# IJulia, we want to automatically display any figures that were
# queued (via redisplay in draw_if_interactive) at the end of the cell
# execution, and then close them.  However, we don't want to
# close/display *all* figures, as there may be some other figures that
# the user is keeping track of in some other way, e.g. for interactive
# widgets.  So, we only want to close figures in the
# IJulia.displayqueue.  Furthermore, if the user explicitly calls
# display() on a figure or does so implicitly by returning a Figure
# object from the cell, we still want to eventually close the figure
# even though display removes it from the displayqueue, so we need to
# keep track of it in a separate closequeue in that case.

# queue of figures that need closing despite being removed from displayqueue
const closequeue = Figure[]
function pushclose(f::Figure)
    if !in(f, closequeue)
        push!(closequeue, f)
    end
    return f
end

function display_figs()
    if drew_something[1] && isjulia_display[1]
        for manager in Gcf[:get_all_fig_managers]()
            display(pushclose(Figure(manager["canvas"]["figure"])))
        end
        drew_something[1] = false # reset until next drawing command
    end
    nothing
end

function close_queued_figs()
    if isjulia_display[1] && (drew_something[1] || !isempty(closequeue))
        for f in Main.IJulia.displayqueue
            if isa(f, Figure)
                pltm[:close](f[:number])
            end
        end
        for f in closequeue
            pltm[:close](f[:number])
        end
        empty!(closequeue)
        drew_something[1] = false # reset until next drawing command
    end
end

# hook to force new IJulia cells to create new figure objects
const gcf_isnew = [false] # true to force next gcf() to be new figure
force_new_fig() = gcf_isnew[1] = true

# monkey-patch gcf() and figure() so that we can force the creation
# of new figures in new IJulia cells (e.g. after @manipulate commands
# that leave the figure from the previous cell open).
const orig_gcf = pltm["gcf"]
const orig_figure = pltm["figure"]
function figure(args...; kws...)
    gcf_isnew[1] = false
    pycall(orig_figure, PyAny, args...; kws...)
end
function gcf()
    if isjulia_display[1] && gcf_isnew[1]
        return figure()
    else
        return pycall(orig_gcf, PyAny)
    end
end

function monkeypatch()
    pltm["draw_if_interactive"] = draw_if_interactive
    pltm["show"] = display_figs
    pltm["gcf"] = gcf
    pltm["figure"] = figure
end

if isdefined(Main, :IJulia) && Main.IJulia.inited
    Main.IJulia.push_preexecute_hook(force_new_fig)
    Main.IJulia.push_postexecute_hook(close_queued_figs)
    Main.IJulia.push_posterror_hook(close_queued_figs)
end

if isjulia_display[1]
    if backend != "Agg"
        pltm[:switch_backend]("Agg")
    end
    monkeypatch()
end

###########################################################################
# Base.Help.FUNCTION_DICT is undocumented, but it is better than nothing
# until Julia gets a documented docstring-like facility.

function addhelp(f::String, o::PyObject)
    try
        Base.Help.init_help()
        if haskey(o, "__doc__")
            if !haskey(Base.Help.FUNCTION_DICT, f)
                Base.Help.FUNCTION_DICT[f] = Any[]
            end
            push!(Base.Help.FUNCTION_DICT[f], convert(String, o["__doc__"]))
        end
    end
end
addhelp(f::Symbol, o::PyObject) = addhelp(string(f), o)
addhelp(f, o::PyObject, key::String) = haskey(o, key) && addhelp(f, o[key])

###########################################################################

const plt = pywrap(pltm)

# export documented pyplot API (http://matplotlib.org/api/pyplot_api.html)
export acorr,annotate,arrow,autoscale,autumn,axes,axhline,axhspan,axis,axvline,axvspan,bar,barbs,barh,bone,box,boxplot,broken_barh,cla,clabel,clf,clim,cohere,colorbar,colors,contour,contourf,cool,copper,csd,delaxes,disconnect,draw,errorbar,eventplot,figimage,figlegend,figtext,figure,fill_between,fill_betweenx,findobj,flag,gca,gcf,gci,get_current_fig_manager,get_figlabels,get_fignums,get_plot_commands,ginput,gray,grid,hexbin,hist2d,hlines,hold,hot,hsv,imread,imsave,imshow,ioff,ion,ishold,jet,legend,locator_params,loglog,margins,matshow,minorticks_off,minorticks_on,over,pause,pcolor,pcolormesh,pie,pink,plot,plot_date,plotfile,polar,prism,psd,quiver,quiverkey,rc,rc_context,rcdefaults,rgrids,savefig,sca,scatter,sci,semilogx,semilogy,set_cmap,setp,show,specgram,spectral,spring,spy,stackplot,stem,step,streamplot,subplot,subplot2grid,subplot_tool,subplots,subplots_adjust,summer,suptitle,switch_backend,table,text,thetagrids,tick_params,ticklabel_format,tight_layout,title,tricontour,tricontourf,tripcolor,triplot,twinx,twiny,vlines,waitforbuttonpress,winter,xkcd,xlabel,xlim,xscale,xticks,ylabel,ylim,yscale,yticks

for f in (:acorr,:annotate,:arrow,:autoscale,:autumn,:axes,:axhline,:axhspan,:axis,:axvline,:axvspan,:bar,:barbs,:barh,:bone,:box,:boxplot,:broken_barh,:cla,:clabel,:clf,:clim,:cohere,:colorbar,:colors,:contour,:contourf,:cool,:copper,:csd,:delaxes,:disconnect,:draw,:errorbar,:eventplot,:figimage,:figlegend,:figtext,:fill_between,:fill_betweenx,:findobj,:flag,:gca,:gci,:get_current_fig_manager,:get_figlabels,:get_fignums,:get_plot_commands,:ginput,:gray,:grid,:hexbin,:hist2d,:hlines,:hold,:hot,:hsv,:imread,:imsave,:imshow,:ioff,:ion,:ishold,:jet,:legend,:locator_params,:loglog,:margins,:matshow,:minorticks_off,:minorticks_on,:over,:pause,:pcolor,:pcolormesh,:pie,:pink,:plot,:plot_date,:plotfile,:polar,:prism,:psd,:quiver,:quiverkey,:rc,:rc_context,:rcdefaults,:rgrids,:savefig,:sca,:scatter,:sci,:semilogx,:semilogy,:set_cmap,:setp,:specgram,:spectral,:spring,:spy,:stackplot,:stem,:streamplot,:subplot,:subplot2grid,:subplot_tool,:subplots,:subplots_adjust,:summer,:suptitle,:switch_backend,:table,:text,:thetagrids,:tick_params,:ticklabel_format,:tight_layout,:title,:tricontour,:tricontourf,:tripcolor,:triplot,:twinx,:twiny,:vlines,:waitforbuttonpress,:winter,:xkcd,:xlabel,:xlim,:xscale,:xticks,:ylabel,:ylim,:yscale,:yticks)
    py_f = symbol(string("py_", f))
    sf = string(f)
    if haskey(pltm, sf)
        @eval begin
            const $py_f = pltm[$sf]
            $f(args...; kws...) = pycall($py_f, PyAny, args...; kws...)
        end
        addhelp(f, pltm[sf])
    else # using a different (older?) version of matplotlib
        @eval $f(args...; kws...) = error("matplotlib ", version,
                                          " does not have pyplot.", $sf)
    end
end

addhelp("figure", orig_figure)
addhelp("gcf", orig_gcf)

# The following pyplot functions must be handled specially since they
# overlap with standard Julia functions:
#          close, connect, fill, hist, xcorr

import Base: close, connect, fill, step

show() = display_figs()

const py_step = pltm["step"]
step(x, y; kws...) = pycall(py_step, PyAny, x, y; kws...)
addhelp("PyPlot.step", py_step)

const py_close = pltm["close"]
close(f::Union(Figure,String,Symbol,Integer)) = pycall(py_close, PyAny, f)
close() = pycall(py_close, PyAny)
addhelp("PyPlot.close", py_close)

const py_connect = pltm["connect"]
connect(s::Union(String,Symbol), f::Function) = pycall(py_connect, PyAny, s, f)
addhelp("PyPlot.connect", py_connect)

const py_fill = pltm["fill"]
fill(x::AbstractArray,y::AbstractArray, args...; kws...) =
    pycall(py_fill, PyAny, x, y, args...; kws...)
addhelp("PyPlot.fill", py_fill)

# no way to use method dispatch for hist or xcorr, since their
# argument signatures look too much like Julia's

include("colormaps.jl")

###########################################################################
# Support array of string labels in bar chart

function bar{T<:String}(x::AbstractVector{T}, y; kws...)
    xi = 1:length(x)
    if !any(kw -> kw[1] == :align, kws)
        push!(kws, (:align, "center"))
    end
    p = bar(xi, y; kws...)
    ax = any(kw -> kw[1] == :orientation && lowercase(kw[2]) == "horizontal",
             kws) ? gca()["yaxis"] : gca()["xaxis"]
    ax[:set_ticks](xi)
    ax[:set_ticklabels](x)
    return p
end

bar{T<:Symbol}(x::AbstractVector{T}, y; kws...) =
    bar(map(string, x), y; kws...)

###########################################################################
# Include mplot3d for 3d plotting.

export art3d, Axes3D, surf, mesh, bar3d, bar3D, contour3D, contourf3D, plot3D, plot_surface, plot_trisurf, plot_wireframe, scatter3D, text2D, text3D, zlabel, zlim, zscale, zticks

const mplot3d = pyimport("mpl_toolkits.mplot3d")
const axes3d = pyimport("mpl_toolkits.mplot3d.axes3d")

const art3d = pywrap(pyimport("mpl_toolkits.mplot3d.art3d"))
const Axes3D = axes3d[:Axes3D]

for f in (:bar3d, :contour3D, :contourf3D, :plot3D, :plot_surface,
          :plot_trisurf, :plot_wireframe, :scatter3D, :text2D, :text3D)
    fs = string(f)
    @eval function $f(args...; kws...)
        ax = gca(projection="3d")
        pycall(ax[$fs], PyAny, args...; kws...)
    end
    addhelp(fs, axes3d["Axes3D"], fs)
end
const bar3D = bar3d # correct for annoying mplot3d inconsistency
addhelp("bar3D", axes3d["Axes3D"], "bar3d")

# it's annoying to have xlabel etc. but not zlabel
for f in (:zlabel, :zlim, :zscale, :zticks)
    fs = string("set_", f)
    @eval function $f(args...; kws...)
        ax = gca(projection="3d")
        pycall(ax[$fs], PyAny, args...; kws...)
    end
    addhelp(f, axes3d["Axes3D"], fs)
end

# export Matlab-like names

function surf(Z::AbstractMatrix; kws...)
    plot_surface([1:size(Z,1)]*ones(1,size(Z,2)),
                 ones(size(Z,1))*[1:size(Z,2)]', Z; kws...)
end

function surf(X, Y, Z::AbstractMatrix, args...; kws...)
    plot_surface(X, Y, Z, args...; kws...)
end

function surf(X, Y, Z::AbstractVector, args...; kws...)
    plot_trisurf(X, Y, Z, args...; kws...)
end

mesh(args...; kws...) = plot_wireframe(args...; kws...)

function mesh(Z::AbstractMatrix; kws...)
    plot_wireframe([1:size(Z,1)]*ones(1,size(Z,2)),
                   ones(size(Z,1))*[1:size(Z,2)]', Z; kws...)
end

addhelp(:surf, axes3d["Axes3D"], "plot_surface")
addhelp(:mesh, axes3d["Axes3D"], "plot_wireframe")

###########################################################################
# Allow plots with 2 independent variables (contour, surf, ...)
# to accept either 2 1d arrays or a row vector and a 1d array,
# to simplify construction of such plots via broadcasting operations.
# (Matplotlib is inconsistent about this.)

for f in (:contour, :contourf)
    @eval function $f(X::AbstractMatrix, Y::AbstractVector, args...; kws...)
        if size(X,1) == 1 || size(X,2) == 1
            $f(reshape(X, length(X)), Y, args...; kws...)
        else
            throw(ArgumentError("if 2nd arg is column vector, 1st arg must be row or column vector"))
        end
    end
end

for f in (:surf,:mesh,:plot_surface,:plot_wireframe,:contour3D,:contourf3D)
    @eval begin
        function $f(X::AbstractVector, Y::AbstractVector, Z::AbstractMatrix, args...; kws...)
            m, n = length(X), length(Y)
            $f(repmat(reshape(X,1,m),n,1), repmat(Y,1,m), Z, args...; kws...)
        end
        function $f(X::AbstractMatrix, Y::AbstractVector, Z::AbstractMatrix, args...; kws...)
            if size(X,1) != 1 && size(X,2) != 1
                throw(ArgumentError("if 2nd arg is column vector, 1st arg must be row or column vector"))
            end
            m, n = length(X), length(Y)
            $f(repmat(reshape(X,1,m),n,1), repmat(Y,1,m), Z, args...; kws...)
        end
    end
end

# Already work: barbs, pcolor, pcolormesh, quiver

# Matplotlib pcolor* functions accept 1d arrays but not ranges
for f in (:pcolor, :pcolormesh)
    @eval begin
        $f(X::Range, Y::Range, args...; kws...) = $f([X...], [Y...], args...; kws...)
        $f(X::Range, Y::AbstractArray, args...; kws...) = $f([X...], Y, args...; kws...)
        $f(X::AbstractArray, Y::Range, args...; kws...) = $f(X, [Y...], args...; kws...)
    end
end

###########################################################################
# a more pure functional style, that returns the figure but does *not*
# have any display side-effects.  Mainly for use with @manipulate (Interact.jl)

function withfig(actions::Function, f::Figure; clear=true)
    ax_save = gca()
    figure(f[:number])
    finalizer(f, close)
    try
        if clear && !isempty(f)
            clf()
        end
        actions()
    catch
        rethrow()
    finally
        try
            sca(ax_save) # may fail if axes were overwritten
        end
        Main.IJulia.undisplay(f)
    end
    return f
end

###########################################################################

using LaTeXStrings
export LaTeXString, latexstring, @L_str, @L_mstr

###########################################################################

end # module PyPlot


text:
import BenchPerf
import Random
using BenchmarkTools

include("utils.jl")

function sum_gather(xs, indices)
    s = zero(eltype(xs))
    for i in indices
        s += @inbounds xs[i]
    end
    return s
end

function threaded_sum_gather(xss, iss, nrepeats)
    sums = zeros(eltype(eltype(xss)), length(xss))
    Threads.@threads for j in eachindex(xss, iss)
        s = sums[j]
        for _ in 1:nrepeats
            s += sum_gather(xss[j], iss[j])
        end
        sums[j] = s
    end
    return sum(sums)
end

group = BenchmarkGroup()
for e in 5:22
    n = 2^e
    nrepeats = nrepeats_from_n(n)
    group["n=$n"] = bench = @benchmarkable(
        threaded_sum_gather(xss, iss, $nrepeats),
        setup = begin
            Random.seed!(43)
            # Iteration over nthreads can be done in a single process.  But it's
            # done via BenchmarkConfigSweeps just for the demo.
            nt = Threads.nthreads()
            n = $n
            xss = [zeros(Int, n) for _ in 1:nt]
            iss = [rand(1:n, n) for _ in 1:nt]
        end,
        # Larger `samples` (than the default: 10_000) to bound the benchmark by
        # the time limit.
        samples = 1_000_000,
    )
    # Not tuning it like ../random_increments/benchmarks.jl since `nrepeats`
    # already takes care of it. To verify this, uncomment the following line,
    # run the benchmarks, and then check `unique(df.evals)`.
    #=
    if e < 15
        tune!(bench)
    end
    =#
end

include("../quick.jl")
maybe_quick!(group)

SUITE = BenchPerf.wrap(group; detailed = 1)

ext: jl


text:
#!/usr/bin/env julia --project
using ClimateMachine
ClimateMachine.init(parse_clargs = true)

using ClimateMachine.MPIStateArrays: euclidean_distance, norm
using ClimateMachine.GenericCallbacks
using ClimateMachine.ODESolvers
using ClimateMachine.Mesh.Filters
using ClimateMachine.VariableTemplates
using ClimateMachine.Mesh.Grids: polynomialorder
using ClimateMachine.Ocean.HydrostaticBoussinesq
using ClimateMachine.Ocean.OceanProblems

using Test

using CLIMAParameters
using CLIMAParameters.Planet: grav
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

function config_simple_box(FT, N, resolution, dimensions; BC = nothing)
    if BC == nothing
        problem = SimpleBox{FT}(dimensions...)
    else
        problem = SimpleBox{FT}(dimensions...; BC = BC)
    end

    model = HydrostaticBoussinesqModel{FT}(
        param_set,
        problem;
        cʰ = FT(1),
        αᵀ = FT(0),
        κʰ = FT(0),
        κᶻ = FT(0),
        fₒ = FT(0),
        β = FT(0),
    )

    config = ClimateMachine.OceanBoxGCMConfiguration(
        "hydrostatic_spindown",
        N,
        resolution,
        param_set,
        model;
        periodicity = (true, true, false),
        boundary = ((0, 0), (0, 0), (1, 2)),
    )

    return config
end

function run_hydrostatic_test(; imex::Bool = false, BC = nothing, refDat = ())
    FT = Float64

    # DG polynomial order
    N = Int(4)

    # Domain resolution and size
    Nˣ = Int(5)
    Nʸ = Int(5)
    Nᶻ = Int(8)
    resolution = (Nˣ, Nʸ, Nᶻ)

    Lˣ = 1e6    # m
    Lʸ = 1e6    # m
    H = 400   # m
    dimensions = (Lˣ, Lʸ, H)

    timestart = FT(0)   # s
    timeout = FT(6 * 3600)  # s
    timeend = FT(86400) # s
    dt = FT(120)    # s

    if imex
        solver_type =
            ClimateMachine.IMEXSolverType(implicit_model = LinearHBModel)
    else
        solver_type = ClimateMachine.ExplicitSolverType(
            solver_method = LSRK144NiegemannDiehlBusch,
        )
    end

    driver = config_simple_box(FT, N, resolution, dimensions; BC = BC)

    grid = driver.grid
    vert_filter = CutoffFilter(grid, polynomialorder(grid) - 1)
    exp_filter = ExponentialFilter(grid, 1, 8)
    modeldata = (vert_filter = vert_filter, exp_filter = exp_filter)

    solver = ClimateMachine.SolverConfiguration(
        timestart,
        timeend,
        driver,
        init_on_cpu = true,
        ode_dt = dt,
        ode_solver_type = solver_type,
        modeldata = modeldata,
    )

    output_interval = ceil(Int64, timeout / solver.dt)

    ClimateMachine.Settings.vtk = "never"
    # ClimateMachine.Settings.vtk = "$(output_interval)steps"

    ClimateMachine.Settings.diagnostics = "never"
    # ClimateMachine.Settings.diagnostics = "$(output_interval)steps"

    cb = ClimateMachine.StateCheck.sccreate(
        [(solver.Q, "state"), (solver.dg.state_auxiliary, "aux")],
        output_interval;
        prec = 12,
    )

    result = ClimateMachine.invoke!(solver; user_callbacks = [cb])

    Q_exact = ClimateMachine.DGMethods.init_ode_state(
        solver.dg,
        timeend,
        solver.init_args...;
        init_on_cpu = solver.init_on_cpu,
    )

    error = euclidean_distance(solver.Q, Q_exact) / norm(Q_exact)

    println("error = ", error)
    @test isapprox(error, FT(0.0); atol = 0.005)

    ## Check results against reference
    ClimateMachine.StateCheck.scprintref(cb)
    if length(refDat) > 0
        @test ClimateMachine.StateCheck.scdocheck(cb, refDat)
    end
end

@testset "$(@__FILE__)" begin

    include("../refvals/3D_hydrostatic_spindown_refvals.jl")

    run_hydrostatic_test(imex = false, refDat = refVals.explicit) # error = 0.0011289879366523504
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
# Minimal implementation of 3-layer MLP that solves Xor problem.

X = [[0, 0], [0, 1], [1, 0], [1, 1]]
Y = (x -> [xor(x...)]).(X)
η = 0.1
W = rand(2, 3)
V = rand(1, 3)
W ∗ x = W * [x; 1]
for _ in 1:1000, (x, y) in zip(X, Y)
    z = W ∗ x
    h = tanh.(z)
    o = V ∗ h
    ∇o = o - y
    ∇h = ∇o * V[:, 1:(end - 1)]
    ∇V = ∇o * [h; 1]'
    ∇W = (@. ∇h * (1.0 - h' ^ 2))' * [x; 1]'
    V .-= η * ∇V
    W .-= η * ∇W
end
predict(x) = round.(Int, V ∗ tanh.(W ∗ x))
@assert predict.(X) == Y

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
"""
    basicblock(inplanes, outplanes, downsample = false)

Create a basic residual block
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: a list of the number of output feature maps for each convolution
               within the residual block
- `downsample`: set to `true` to downsample the input
"""
basicblock(inplanes, outplanes, downsample = false) = downsample ?
  Chain(conv_bn((3, 3), inplanes, outplanes[1]; stride = 2, pad = 1, usebias = false)...,
        conv_bn((3, 3), outplanes[1], outplanes[2]; stride = 1, pad = 1, usebias = false)...) :
  Chain(conv_bn((3, 3), inplanes, outplanes[1]; stride = 1, pad = 1, usebias = false)...,
        conv_bn((3, 3), outplanes[1], outplanes[2]; stride = 1, pad = 1, usebias = false)...)

"""
    bottleneck(inplanes, outplanes, downsample = false)

Create a bottleneck residual block
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: a list of the number of output feature maps for each convolution
               within the residual block
- `downsample`: set to `true` to downsample the input
"""
bottleneck(inplanes, outplanes, downsample = false) = downsample ?
  Chain(conv_bn((1, 1), inplanes, outplanes[1]; stride = 2, usebias = false)...,
        conv_bn((3, 3), outplanes[1], outplanes[2]; stride = 1, pad = 1, usebias = false)...,
        conv_bn((1, 1), outplanes[2], outplanes[3]; stride = 1, usebias = false)...) :
  Chain(conv_bn((1, 1), inplanes, outplanes[1]; stride = 1, usebias = false)...,
        conv_bn((3, 3), outplanes[1], outplanes[2]; stride = 1, pad = 1, usebias = false)...,
        conv_bn((1, 1), outplanes[2], outplanes[3]; stride = 1, usebias = false)...)

"""
    skip_projection(inplanes, outplanes, downsample = false)

Create a skip projection
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: a list of the number of output feature maps
- `downsample`: set to `true` to downsample the input
"""
skip_projection(inplanes, outplanes, downsample = false) = downsample ?
  Chain(conv_bn((1, 1), inplanes, outplanes; stride = 2, usebias = false)...) :
  Chain(conv_bn((1, 1), inplanes, outplanes; stride = 1, usebias = false)...)

# array -> PaddedView(0, array, outplanes) for zero padding arrays
"""
    skip_identity(inplanes, outplanes)

Create a identity projection
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments:
- `inplanes`: the number of input feature maps
- `outplanes`: a list of the number of output feature maps
"""
function skip_identity(inplanes, outplanes)
  if outplanes[end] > inplanes
    return Chain(MaxPool((1, 1), stride = 2),
                 y -> cat(y, zeros(eltype(y),
                                   size(y, 1),
                                   size(y, 2),
                                   outplanes[end] - inplanes, size(y, 4)); dims = 3))
  else
    return identity
  end
end

"""
    resnet(; block, shortcut_config, channel_config, block_config, nclasses = 1000)

Create a ResNet model
([reference](https://arxiv.org/abs/1512.03385v1)).

# Arguments
- `block`: a function with input `(inplanes, outplanes, downsample=false)` that returns
           a new residual block (see [`Metalhead.basicblock`](#) and [`Metalhead.bottleneck`](#))
- `shortcut_config`: the type of shortcut style (either `:A`, `:B`, or `:C`)
- `channel_config`: the growth rate of the output feature maps within a residual block
- `block_config`: a list of the number of residual blocks at each stage
- `nclasses`: the number of output classes
"""
function resnet(; block, shortcut_config, channel_config, block_config, nclasses = 1000)
  inplanes = 64
  baseplanes = 64
  layers = []
  append!(layers, conv_bn((7, 7), 3, inplanes; stride = 2, pad = (3, 3)))
  push!(layers, MaxPool((3, 3), stride = (2, 2), pad = (1, 1)))
  for (i, nrepeats) in enumerate(block_config)
    outplanes = baseplanes .* channel_config
    if shortcut_config == :A
      push!(layers, Parallel(+, block(inplanes, outplanes, i != 1),
                                skip_identity(inplanes, outplanes)))
    elseif shortcut_config == :B || shortcut_config == :C
      push!(layers, Parallel(+, block(inplanes, outplanes, i != 1),
                                skip_projection(inplanes, outplanes[end], i != 1)))
    end
    inplanes = outplanes[end]
    for _ in 2:nrepeats
      if shortcut_config == :A || shortcut_config == :B
        push!(layers, Parallel(+, block(inplanes, outplanes, false),
                                  skip_identity(inplanes, outplanes[end])))
      elseif shortcut_config == :C
        push!(layers, Parallel(+, block(inplanes, outplanes, false),
                                  skip_projection(inplanes, outplanes, false)))
      end
      inplanes = outplanes[end]
    end
    baseplanes *= 2
  end

  return Chain(Chain(layers..., AdaptiveMeanPool((1, 1))),
               Chain(flatten, Dense(inplanes, nclasses)))
end

const resnet_config =
  Dict(:resnet18 => ([1, 1], [2, 2, 2, 2], :A),
       :resnet34 => ([1, 1], [3, 4, 6, 3], :A),
       :resnet50 => ([1, 1, 4], [3, 4, 6, 3], :B),
       :resnet101 => ([1, 1, 4], [3, 4, 23, 3], :B),
       :resnet152 => ([1, 1, 4], [3, 8, 36, 3], :B))

"""
    ResNet(channel_config, block_config, shortcut_config; block, nclasses = 1000)

Create a `ResNet` model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`resnet`](#).

# Arguments
- `channel_config`: the growth rate of the output feature maps within a residual block
- `block_config`: a list of the number of residual blocks at each stage
- `shortcut_config`: the type of shortcut style (either `:A`, `:B`, or `:C`)
- `block`: a function with input `(inplanes, outplanes, downsample=false)` that returns
           a new residual block (see [`Metalhead.basicblock`](#) and [`Metalhead.bottleneck`](#))
- `nclasses`: the number of output classes
"""
struct ResNet{T}
  layers::T
end

function ResNet(channel_config, block_config, shortcut_config; block, nclasses = 1000)
  layers = resnet(block = block,
                  shortcut_config = shortcut_config,
                  channel_config = channel_config,
                  block_config = block_config,
                  nclasses = nclasses)

  ResNet(layers)
end

@functor ResNet

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

"""
    ResNet18(; pretrain = false, nclasses = 1000)

Create a ResNet-18 model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`Metalhead.ResNet`](#).

# Arguments
- `nclasses`: the number of output classes

!!! warning
    `ResNet18` does not currently support pretrained weights.
"""
function ResNet18(; pretrain = false, nclasses = 1000)
  model = ResNet(resnet_config[:resnet18]...; block = basicblock, nclasses = nclasses)

  pretrain && pretrain_error("ResNet18")
  return model
end

"""
    ResNet34(; pretrain = false, nclasses = 1000)

Create a ResNet-34 model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`Metalhead.ResNet`](#).

# Arguments
- `pretrain`: set to `true` to load pre-trained weights for ImageNet
- `nclasses`: the number of output classes

!!! warning
    `ResNet34` does not currently support pretrained weights.
"""
function ResNet34(; pretrain = false, nclasses = 1000)
  model = ResNet(resnet_config[:resnet34]...; block = basicblock, nclasses = nclasses)

  pretrain && pretrain_error("ResNet34")
  return model
end

"""
    ResNet50(; pretrain = false, nclasses = 1000)

Create a ResNet-50 model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`Metalhead.ResNet`](#).

# Arguments
- `pretrain`: set to `true` to load pre-trained weights for ImageNet
- `nclasses`: the number of output classes
"""
function ResNet50(; pretrain = false, nclasses = 1000)
  model = ResNet(resnet_config[:resnet50]...; block = bottleneck, nclasses = nclasses)

  pretrain && Flux.loadparams!(model.layers, weights("resnet50"))
  return model
end

"""
    ResNet101(; pretrain = false, nclasses = 1000)

Create a ResNet-101 model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`Metalhead.ResNet`](#).

# Arguments
- `pretrain`: set to `true` to load pre-trained weights for ImageNet
- `nclasses`: the number of output classes

!!! warning
    `ResNet101` does not currently support pretrained weights.
"""
function ResNet101(; pretrain = false, nclasses = 1000)
  model = ResNet(resnet_config[:resnet101]...; block = bottleneck, nclasses = nclasses)

  pretrain && pretrain_error("ResNet101")
  return model
end

"""
    ResNet152(; pretrain = false, nclasses = 1000)

Create a ResNet-152 model
([reference](https://arxiv.org/abs/1512.03385v1)).
See also [`Metalhead.ResNet`](#).

# Arguments
- `pretrain`: set to `true` to load pre-trained weights for ImageNet
- `nclasses`: the number of output classes

!!! warning
    `ResNet152` does not currently support pretrained weights.
"""
function ResNet152(; pretrain = false, nclasses = 1000)
  model = ResNet(resnet_config[:resnet152]...; block = bottleneck, nclasses = nclasses)

  pretrain && pretrain_error("ResNet152")
  return model
end


text:
"""
The `UnmeasuredElementRule` mechanism provides a method to implement rules for adding unmeasured elements to
the fitting process.  Examples include element-by-stoichiometry or element-by-difference.

`UnmeasuredElementRule` implementations must include methods for:

    NeXLUncertainties.compute(uer::<:UnmeasuredElementRule, inp::Dict{Element,Float64})::Dict{Element,Float64}
    isunmeasured(uer::<:UnmeasuredElementRule, elm::Element)::Bool
"""
abstract type UnmeasuredElementRule end

"""
The NullUnmeasuredRule adds no additional elements in the iteration process.
"""
struct NullUnmeasuredRule <: UnmeasuredElementRule end

"""
    NeXLUncertainties.compute(::NullUnmeasuredRule, inp::Dict{Element,Float64})::Dict{Element,Float64}

A null UnmeasuredElementRule.  Just returns the inputs.
"""
NeXLUncertainties.compute(::NullUnmeasuredRule, inp::Dict{Element,<:AbstractFloat})::Dict{Element,Float64} = inp

isunmeasured(::NullUnmeasuredRule, elm::Element) = false

"""
The ElementByStoichiometry computes the oxygen mass fraction using the standard stoichiometric rules.
"""
struct ElementByStoichiometry <: UnmeasuredElementRule
    element::Element
    valences
    ElementByStoichiometry(elm::Element, val) = new(elm, val)
end

"""
    OByStoichiometry(valences = NeXLCore.valences)

Returns an ElementByStoichiometry for computing O using the valences provided as a tuple or Array indexed by z.
"""
OByStoichiometry(valences = NeXLCore.defaultValences) = ElementByStoichiometry(n"O", valences)

"""
    NeXLUncertainties.compute(stoic::ElementByStoichiometry, inp::Dict{Element,<:AbstractFloat})::Dict{Element,Float64}

Computes an element using stoichiometric rules.
"""
function NeXLUncertainties.compute(stoic::ElementByStoichiometry, inp::Dict{Element,<:AbstractFloat})::Dict{Element,Float64}
    inp[stoic.element] = obystoichiometry(inp, valences = stoic.valences)
    inp
end

function isunmeasured(stoic::ElementByStoichiometry, elm::Element)::Bool
    isequal(elm, stoic.element)
end

function Base.show(io::IO, stoic::ElementByStoichiometry)
    print(io,"ElementByStoichiometry[$(stoic.element)]")
end

struct ElementByFiat <: UnmeasuredElementRule
    element::Element
    massfraction::Float64
end

"""
    NeXLUncertainties.compute(stoic::ElementByFiat, inp::Dict{Element,<:AbstractFloat})::Dict{Element,Float64}

Simply assigns a quantity for an element.
"""
function NeXLUncertainties.compute(fiat::ElementByFiat, inp::Dict{Element,<:AbstractFloat})::Dict{Element,Float64}
    inp[fiat.element] = fiat.massfraction
    inp
end

function isunmeasured(fiat::ElementByFiat, elm::Element)::Bool
    isequal(elm, fiat.element)
end

function Base.show(io::IO, fiat::ElementByFiat)
    print(io,"ElementByFiat[$(fiat.element), $(fiat.massfraction)]")
end
"""
    MultiUnmeasuredElementRule

Apply multiple UnmeasuredElementRule items sequentially.
"""
struct MultiUnmeasuredElementRule <: UnmeasuredElementRule
    rules::Vector{<:UnmeasuredElementRule}
end

function NeXLUncertainties.compute(muer::MultiUnmeasuredElementRule, inp::Dict{Element,<:AbstractFloat})::Dict{Element,Float64}
    res = inp
    for rule in muer.rules
        res = compute(rule, res)
    end
    return res
end

function isunmeasured(muer::MultiUnmeasuredElementRule, elm::Element)::Bool
    any(rule->isunmeasured(rule, elm), muer.rules)
end

function Base.show(io::IO, muer::MultiUnmeasuredElementRule)
    print(io,"MultiUnmeasuredElementRule[$(join(muer.rules, ", "))]")
end

max_stars_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_issues_repo_name: NicholasWMRitchie/NeXLMatrixCorrection

max_forks_repo_name: NicholasWMRitchie/NeXLMatrixCorrection


text:
import Base: print, replace!, show

# Indentation for output.
const INDENT_STEP = 4

abstract type TreeElement end

mutable struct Sample <: TreeElement
    fields::Array{String}
end

mutable struct Clade <: TreeElement
    name::String
    fields::Array{String}
    samples::Array{Sample}
    subclades::Array{Clade}
end

function Clade(fields::Array{<:AbstractString})
    name = strip(fields[1])
    Clade(name, fields, [], [])
end

function readtree(filename::AbstractString)
    # Read tree from file in text format.
    lines = String[]
    open(filename) do f
        for line in eachline(f)
            if !startswith(line, "//") && !all(isspace, line)
                push!(lines, line)
            end
        end
    end
    # Build tree.
    if !isempty(lines)
        root = parsetree!(lines)
    else
        error("Empty file.")
    end
end

function readtrees(filenames::Array{<:AbstractString})
    if isempty(filenames)
        error("No filenames provided.")
    end
    root = readtree(filenames[1])
    if length(filenames) > 1
        for filename in filenames[2:end]
            tree = readtree(filename)
            root = mergetrees!(root, tree)
        end
    end
    return root
end

function parsetree!(lines::Array{String})
    # Create top node.
    line = popfirst!(lines)
    indent = indentation(line)
    fields = split(strip(line), ",")
    node = Clade(fields)

    # Parse samples and subclades
    while !isempty(lines)
        if indentation(lines[1]) <= indent
            break
        end
        content = strip(lines[1])
        if startswith(content, "id:")
            fields = split(content, ",")
            push!(node.samples, Sample(fields))
            popfirst!(lines)
        else
            subclade = parsetree!(lines)
            push!(node.subclades, subclade)
        end
    end
    return node
end

function Base.show(io::IO, node::Clade)
    print(io, node, 0)
end

function print(io::IO, node::Clade, indent::Int)
    subindent = indent + INDENT_STEP
    printspaces(io, indent)
    printfields(io, node.fields)
    println(io)
    for sample in node.samples
        printspaces(io, subindent)
        printfields(io, sample.fields)
        println(io)
    end
    for subclade in node.subclades
        print(io, subclade, subindent)
    end
end

function writetree(filename::String, tree::Clade)
    buffer = IOBuffer()
    print(buffer, tree, 0)
    open(filename, "w") do f
        write(f, String(take!(buffer)))
    end
end

function mergetrees!(tree::Clade, subtree::Clade)
    for i in eachindex(tree.subclades)
        if tree.subclades[i].name == subtree.name
            tree.subclades[i] = subtree
            break
        else
            mergetrees!(tree.subclades[i], subtree)
        end
    end
    return tree
end

function replace!(tree::Clade, pat_f::Pair)
    for i in eachindex(tree.fields)
        tree.fields[i] = replace(tree.fields[i], pat_f)
    end
    for sample in tree.samples
        for i in eachindex(sample.fields)
            sample.fields[i] = replace(sample.fields[i], pat_f)
        end
    end
    for subclade in tree.subclades
        replace!(subclade, pat_f)
    end
end

function deletefields!(tree::Clade, keys::Array{<:AbstractString})
    tree.fields = deletefields(tree.fields, keys)
    for sample in tree.samples
        sample.fields = deletefields(sample.fields, keys)
    end
    for subclade in tree.subclades
        deletefields!(subclade, keys)
    end
end

function deletefields(fields:: Array{<:AbstractString}, keys::Array{<:AbstractString})
    newfields = []
    for field in fields
        isclean = true
        for key in keys
            if occursin(key, field)
                isclean = false
                break
            end
        end
        if isclean
            push!(newfields, field)
        end
    end
    return newfields
end

# Needs Phyloage to be installed.
# https://github.com/yogischogi/phyloage
function calculateTMRCAs(tree::Clade; cal::Float64=140.0, offset::Float64=60.0)
    caltree = deepcopy(tree)
    replace!(caltree, "snps:" => "STR-Count:")

    # Call Phyloage with temporary files to calculate results.
    infile, inio = mktemp(pwd())
    outfile, outio = mktemp(pwd())
    close(inio)
    close(outio)
    writetree(infile, caltree)
    run(`phyloage -treein=$infile -treeout=$outfile -offset=$offset -cal=$cal`)
    resulttree = readtree(outfile)
    rm(infile)
    rm(outfile)

    # Remove clutter from result tree.
    deletefields!(resulttree, ["STR-Count:", "STRs Downstream:", "formed:"])
    return resulttree
end

# Count white spaces.
function indentation(line::String)
    indent = 0
    for char in line
        isspace(char) ? indent += 1 : break
    end
    return indent
end

function printspaces(io::IO, count::Int)
    for i = 1:count
        print(io, " ")
    end
end

function printfields(io::IO, fields::Array{<:AbstractString})
    if isempty(fields)
        return
    end
    for i in 1:length(fields)-1
        print(io, fields[i])
        print(io, ",")
    end
    print(io, fields[end])
end


text:
using WeakRefStrings, Test, Random, Parsers
using DataAPI: refarray, refvalue

include("poslenstrings.jl")
include("inlinestrings.jl")

@testset "WeakRefString{UInt8}" begin
    data = codeunits("hey there sailor")

    str = WeakRefStrings.WeakRefString(pointer(data), 3)

    @test str.len == 3
    @test string(str) == "hey"
    @test String(str) == "hey"
    @test str[1] === 'h'
    @test str[2] === 'e'
    @test str[3] === 'y'

    io = IOBuffer()
    show(io, str)
    @test String(take!(io)) == "\"hey\""

    show(io, typeof(str))
    @test String(take!(io)) == "WeakRefString{UInt8}"
end

@testset "WeakRefString{UInt16}" begin
    # simulate UTF16 data
    data = [0x0068, 0x0065, 0x0079]

    str = WeakRefStrings.WeakRefString(pointer(data), 3)
    @test typeof(str) == WeakRefStrings.WeakRefString{UInt16}
    @test String(str) == "hey"
    @test str.len == 3
end

@testset "WeakRefString{UInt32}" begin
    # simulate UTF32 data with trailing null byte
    data = [0x00000068, 0x00000065, 0x00000079, 0x00000000]

    str = WeakRefStrings.WeakRefString(pointer(data), 4)
    @test typeof(str) == WeakRefStrings.WeakRefString{UInt32}
    @test String(str) == "hey"
    @test str.len == 4
end

@testset "StringVector" begin
    s = "Julia is a name without special letters such as æ, ø, and å. Such letters require more than a single byte when encoded in UTF8"
    @testset "split on $splits" for splits in (['.'], [',', '.', ' '])
        sa     = split(s, splits)
        svinit = StringVector{WeakRefString{UInt8}}(sa)
        @testset "version" for sv in (copy(svinit),
                                   copyto!(similar(svinit), svinit),
                                   convert(StringVector{String}, sa),
                                   StringVector{WeakRefString{UInt8}}(sa))
            @test sa == sv
            @test sort(sa) == sort(sv)
            @test sortperm(sa) == sortperm(sv)
            @test sort!(copy(sa)) == sort(sv)
            @test copy(sv) == sv

            @testset "setindex with WeakRefString" begin
                # important to start with end because of special branch when
                # lengths are empty and setting last element
                tmp = sv[end]
                sv[end] = "🍕"
                @test sv[end] == "🍕"

                sv[end] = tmp
                sv[1]   = sv[1]
                @test sa == sv
            end

            @testset "setindex with String" begin
                sv[1]   = String(sa[1])
                sv[end] = String(sa[end])
                @test sa == sv
            end

            @testset "setindex with missing" begin
                sv1 = StringVector{String}(["foo", "bar"])
                @test_throws MethodError setindex!(sv1, missing, 1)
                sv1 = StringVector{Union{String, Missing}}(["foo", "bar"])
                sv1[1] = missing
                @test sv1[1] === missing
            end

            push!(sv, "🍕") == push!(copy(sa), "🍕")
            @test length(empty!(sv)) == 0
        end
    end

    @testset "resize+index" begin
        sv = StringVector(["JuliaDB", "TextParse"])
        @test resize!(sv, 1)[1] == "JuliaDB"
    end

    @testset "filter regex" begin
        sv = StringVector(["TextParse", "TextParse", "JuliaDB", "TextParse", "TextParse", "TextParse", "TextParse", "JuliaDB", "JuliaDB"])
        sv[end] = "Dagger"
        sv[1] = "Dagger"
        @test length(filter(x->occursin(r"JuliaDB", x), sv)) == 2
    end

    @testset "test WeakRefString element type constructor" begin
        @test eltype(StringVector{WeakRefString}(undef, 1)) <: WeakRefString
    end

    @testset "StringVector constructors" begin
        @test length(StringVector()) == 0
        @test length(StringVector{WeakRefString{UInt8}}()) == 0
        a = StringVector{WeakRefString{UInt8}}(undef, 10)
        @test length(a) == 10
        a = StringVector(undef, 10)
        @test length(a) == 10
        data = Vector{UInt8}("hey there sailor")
        a = StringVector{String}(data, 3)
        a[1] = WeakRefString(pointer(data), 3)
        a[2] = WeakRefString(pointer(data), 3)
        a[3] = WeakRefString(pointer(data), 3)
        @test length(a) == 3
        @test a[1] == "hey"
        @test length(a.buffer) == 16
    end

    @testset "test permute!" begin
        sv = StringVector(["TextParse", "TextParse", "JuliaDB", "TextParse", "TextParse", "TextParse", "TextParse", "JuliaDB", "JuliaDB"])
        permute!(sv, reverse!([1:length(sv);]))
        @test sv[1] == "JuliaDB"
        @test sv[end] == "TextParse"
    end

    @testset "vcat" begin
        sv1 = StringVector{String}(unsafe_wrap(Array, pointer(randstring(1024)), 1024), UInt64[1:10:1000;], ones(UInt32,100)*9);
        sv2 = StringVector{String}(unsafe_wrap(Array, pointer(randstring(1024)), 1024), UInt64[1:10:1000;], ones(UInt32,100)*9);
        sv3 = vcat(sv1, sv2)
        @test length(sv3) == length(sv1) + length(sv2)
    end

    @testset "append!" begin
        sv1 = StringVector{WeakRefString{UInt8}}(["foo", "bar"])
        sv2 = StringVector{String}(["baz", "qux"])

        append!(sv1, sv2)
        @test sv1 == ["foo", "bar", "baz", "qux"]

        append!(sv2, ["yep", "nope"])
        @test sv2 == ["baz", "qux", "yep", "nope"]
    end

    @testset "deleteat!" begin
        sv1 = StringVector{WeakRefString{UInt8}}(["foo", "bar"])
        @test deleteat!(copy(sv1), 1) == ["bar"]
        @test deleteat!(copy(sv1), [1,2]) == []
        @test deleteat!(copy(sv1), 1:2) == []

        sv2 = StringVector{String}(["baz", "qux"])
        @test deleteat!(copy(sv2), 1) == ["qux"]
        @test deleteat!(copy(sv2), [1,2]) == []
        @test deleteat!(copy(sv2), 1:2) == []
    end

    @testset "resize!" begin
        sv1 = StringVector{WeakRefString{UInt8}}(["foo", "bar"])
        resize!(sv1, 0)
        @test length(sv1) == 0
        resize!(sv1, 10)
        @test length(sv1) == 10
    end

    @testset "push!" begin
        sv1 = StringVector{String}(["foo", "bar"])
        push!(sv1, "hey")
        @test length(sv1) == 3
        @test sv1[end] == "hey"
        @test_throws MethodError push!(sv1, missing)
        sv1 = StringVector{Union{String, Missing}}(["foo", "bar"])
        push!(sv1, missing)
        @test length(sv1) == 3
        @test sv1[end] === missing
    end

    @testset "insert!" begin
        sv1 = StringVector{String}(["foo", "bar"])
        insert!(sv1, 2, "hey")
        @test length(sv1) == 3
        @test sv1[2] == "hey"
        @test sv1[3] == "bar"
        @test_throws MethodError insert!(sv1, 2, missing)
        sv1 = StringVector{Union{String, Missing}}(["foo", "bar"])
        insert!(sv1, 2, missing)
        @test sv1[2] === missing
        @test length(sv1) == 3
    end

    @testset "splice!" begin
        sv1 = StringVector{String}(["foo", "bar"])
        v = splice!(sv1, 2)
        @test v == "bar"
        @test length(sv1) == 1
        push!(sv1, "bar")
        v = splice!(sv1, 1:2)
        @test v[1] == "foo"
        @test v[2] == "bar"
        @test length(sv1) == 0
        sv1 = StringVector{String}(["foo", "bar"])
        v = splice!(sv1, 1:2, ["bar", "foo"])
        @test sv1 == reverse(v)
        sv1 = StringVector{String}(["foo", "bar"])
        v = splice!(sv1, 2, ["foo", "foo", "bar"])
        @test length(sv1) == 4
        @test sv1[end] == "bar"
    end

    @testset "prepend!/pushfirst!/pop!/popfirst!" begin
        sv1 = StringVector{String}(["foo", "bar"])
        prepend!(sv1, ["hey", "there"])
        @test length(sv1) == 4
        @test sv1[1] == "hey"
        @test sv1[2] == "there"
        pushfirst!(sv1, "sailor")
        @test length(sv1) == 5
        @test sv1[1] == "sailor"
        v = popfirst!(sv1)
        @test v == "sailor"
        @test length(sv1) == 4
        v = pop!(sv1)
        @test v == "bar"
        @test length(sv1) == 3
    end

    @testset "DataAPI" begin
        a = StringVector(["a", "b", "c"])
        v = refarray(a)
        @test all(v .== a)
        @test eltype(v) == WeakRefString{UInt8}
        for i in 1:3
            @test isequal(refvalue(a, v[i]), a[i])
        end

        b = StringVector(["a", "b", missing])
        w = refarray(b)
        @test w[1] == "a"
        @test w[2] == "b"
        @test ismissing(w[3])
        @test eltype(w) == Union{WeakRefString{UInt8}, Missing}
        for i in 1:3
            @test isequal(refvalue(b, w[i]), b[i])
        end
    end
end


text:
function tourl(path)
    if Sys.iswindows()
        # THere might be a nicer way?
        # Anyways, this seems to be needed on windows
        if !startswith(path, "http")
            path = "file:///" * replace(path, "\\" => "/")
        end
    end
    return repr(path)
end

# NOTE: `save_media` is the function you want to overload
# if you want to create a Gallery with custom types.
# Simply overloading the function should do the trick
# and MakieGallery will take care of the rest.

function save_media(entry, x::Scene, path::String)
    path = joinpath(path, "image.jpg")
    save(FileIO.File(DataFormat{:JPEG}, path), x) # work around FileIO bug for now
    [path]
end

function save_media(entry, x::String, path::String)
    out = joinpath(path, basename(x))
    if out != x
        mv(x, out, force = true)
    end
    [out]
end

function save_media(entry, x::AbstractPlotting.Stepper, path::String)
    # return a list of all file names
    images = filter(x-> endswith(x, ".jpeg") || endswith(x, ".jpg"), readdir(x.folder))
    return map(images) do img
        p = joinpath(x.folder, img)
        out = joinpath(path, basename(p))
        mv(p, out, force = true)
        out
    end
end

function save_media(entry, results::AbstractVector, path::String)
    paths = String[]
    for (i, res) in enumerate(results)
        # Only save supported results
        if res isa Union{Scene, String}
            img = joinpath(path, "image$i.jpg")
            save(FileIO.File(DataFormat{:JPEG}, img), res) # work around FileIO
            push!(paths, img)
        end
    end
    paths
end

function save_media(example, events::RecordEvents, path::String)
    # the path is fixed at record time to be stored relative to the example
    epath = event_path(example, "")
    isfile(epath) || error("Can't find events for example: $(example.unique_name). Please run `record_example_events()`")
    # the current path of RecordEvents is where we now actually want to store the video
    video_path = joinpath(path, "video.mp4")
    record(events.scene, video_path) do io
        replay_events(events.scene, epath) do
            recordframe!(io)
        end
    end
    return [video_path]
end

"""
    embed_image(path::AbstractString)
Returns the html to embed an image
"""
function embed_image(path::AbstractString, alt = "")
    if splitext(path)[2] == "pdf"
        return """

        """
    end
    """

    """
end

"""
    embed_video(path::AbstractString)

Generates a html formatted string for embedding video into Documenter Markdown files
(since `Documenter.jl` doesn't support directly embedding mp4's using ![]() syntax).
"""
function embed_video(path::AbstractString, alt = "")
    """
    <video controls autoplay loop muted>
      <source src=$(tourl(path)) type="video/mp4">
      Your browser does not support mp4. Please use a modern browser like Chrome or Firefox.
    </video>
    """
end

"""
Embeds the most common media types as html
"""
function embed_media(path::String, alt = "")
    file, ext = splitext(path)
    if ext in (".png", ".jpg", ".jpeg", ".JPEG", ".JPG", ".gif", ".pdf", ".svg")
        return embed_image(path, alt)
    elseif ext == ".mp4"
        return embed_video(path, alt)
    else
        error("Unknown media extension: $ext with path: $path")
    end
end

"""
Embeds a vector of media files as HTML
"""
function embed_media(io::IO, paths::AbstractVector{<: AbstractString}, caption = "")
    for (i, path) in enumerate(paths)
        occursin("thumb", path) && continue
        println(io, """

                $(embed_media(path, caption))

        """)
    end
end

"""
Replaces raw html code nodes with an actual RawHTML node.
"""
function preprocess!(md)
    map!(md.content, md.content) do c
        if c isa Markdown.Code && c.language == "@raw html"
            Documenter.Documents.RawHTML(c.code)
        else
            c
        end
    end
    md
end

function md2html(md_path; stylesheets = Vector{String}[])
    open(joinpath(dirname(md_path), "index.html"), "w") do io
        md = preprocess!(Markdown.parse_file(md_path))
        hio = IOBuffer(write=true, read=true)
        for sheetref in stylesheets
            println(hio, "<link rel=\"stylesheet\" href=\"$(sheetref)\">")
        end
        println(io, """
        <!doctype html>

            <meta charset="UTF-8">
            <link rel="stylesheet" href="https://cdn.jsdelivr.net/gh/JuliaDocs/Documenter.jl/assets/html/documenter.css">
            $(String(take!(hio)))

        """)
        println.(Ref(io), string.(HTMLWriter.mdconvert(md)))
        println(io, """

            """
        )
    end
end

"""
    save_markdown(mdpath::String, example::CellEntry, media::Vector{String})

Creates a Markdown representation from an example at `mdpath`.
`media` is a vector of media output files belonging to the example
"""
function save_markdown(mdpath::String, example::CellEntry, media::Vector{String})
    # src = example2source(
    #     example,
    #     scope_start = "",
    #     scope_end = "",
    #     indent = "",
    #     outputfile = (entry, ending)-> string("output", ending)
    # )
    src = example2source(
        example,
        scope_start = "",
        scope_end = "",
        indent = "",
        outputfile = (entry, ending)-> string("output", ending)
    )
    open(mdpath, "w") do io
        println(io, "## ", example.title, "\n")
        println(io, "```julia\n$src\n```")
        println(io, "```@raw html\n")
        embed_media(io, media)
        println(io, "```")
    end
end

"""
    save_highlighted_markdown(path::String, example::CellEntry, media::Vector{String};
                                highlighter = Highlights.Themes.DefaultTheme)
Creates a Markdown representation from an example at `mdpath`.
`media` is a vector of media output files belonging to the example
"""
function save_highlighted_markdown(
                                path::String, example::CellEntry, media::Vector{String},
                                highlighter = Highlights.Themes.DefaultTheme;
                                print_toplevel = true
                            )

    src = example2source(
        example,
        scope_start = "",
        scope_end = "",
        indent = "",
        outputfile = (entry, ending)-> string("output", ending),
        print_toplevel = print_toplevel,
        print_backends = true
    )
    hio = IOBuffer(read = true, write = true)
    highlight(hio, MIME("text/html"), src, Highlights.Lexers.JuliaLexer, highlighter)
    html = String(take!(hio))
    open(path, "w") do io
        println(io, "## ", example.title, "\n")
        println(io, "```@raw html\n$html\n```")
        println(io, "```@raw html\n")
        embed_media(io, media)
        println(io, "```")
    end

end

const last_evaled = Ref{Int}()

"""
    set_last_evaled!(database_idx::Int)
Set's last evaled for resume
"""
function set_last_evaled!(database_idx::Int)
    last_evaled[] = database_idx - 1
    database_idx
end

"""
    set_last_evaled!(database_idx::Int)
Set's last evaled for resume
"""
function set_last_evaled!(unique_name::Symbol)
    idx = findfirst(e-> e.unique_name == unique_name, database)
    if idx === nothing
        error("Unique name $unique_name not found in database")
    end
    last_evaled[] = idx - 1 # minus one, because record_example will start at idx + 1
end

backend_reset_theme!(; kwargs...) = AbstractPlotting.set_theme!(; kwargs...)
# for Plots: default(; kwargs...)

"""
    record_examples(folder = ""; resolution = (500, 500), resume = false)

Records all examples in the database. If error happen, you can fix them and
start record with `resume = true`, to start at the last example that errored.
"""
function record_examples(
        folder = "";
        resolution = (500, 500), resume::Union{Bool, Integer} = false,
        generate_thumbnail = false, display = false,
        display_output_toplevel = true
    )

    function output_path(entry, ending)
        joinpath(folder, "tmp", string(entry.unique_name, ending))
    end
    ispath(folder) || mkpath(folder)
    ispath(joinpath(folder, "tmp")) || mkdir(joinpath(folder, "tmp"))
    result = []
    start = if resume isa Int
        start = resume
    elseif (resume isa Bool) && resume && last_evaled[] >= 0 && last_evaled[] < length(database)
        last_evaled[] + 1
    else
        1
    end
    @info("starting from index $start")
    backend_reset_theme!(resolution = resolution)
    AbstractPlotting.inline!(true)
    @testset "Full Gallery recording" begin
        eval_examples(outputfile = output_path, start = start) do example, value
            uname = example.unique_name
            @testset "$(example.title)" begin
                try
                    subfolder = joinpath(folder, string(uname))
                    outfolder = joinpath(subfolder, "media")
                    ispath(outfolder) || mkpath(outfolder)
                    save_media(example, value, outfolder)
                    push!(result, subfolder)
                    set_last_evaled!(uname)
                    backend_reset_theme!(resolution = resolution) # reset before next example
                    AbstractPlotting.inline!(true)
                    @test true
                    if generate_thumbnail && !isfile(outfolder) && ispath(outfolder)
                        sample = joinpath(outfolder, first(readdir(outfolder)))
                        MakieGallery.generate_thumbnail(sample, joinpath(outfolder, "thumb.jpg"))
                    end
                catch e
                    @warn "Error thrown when evaluating $(example.title)" exception=CapturedException(e, Base.catch_backtrace())
                    @test false
                end
            end
        end
    end
    rm(joinpath(folder, "tmp"), recursive = true, force = true)
    gallery_from_recordings(folder, joinpath(folder, "index.html"); print_toplevel = display_output_toplevel)
    result
end

"""
    gallery_from_recordings(
        folder::String,
        html_out::String = abspath(joinpath(pathof(MakieGallery), "..", "..", "index.html"));
        tags = [
            string.(AbstractPlotting.atomic_function_symbols)...,
            "interaction",
            "record",
            "statsmakie",
            "vbox",
            "layout",
            "legend",
            "colorlegend",
            "vectorfield",
            "poly",
            "camera",
            "recipe",
            "theme",
            "annotations"
        ],
        hltheme = Highlights.Themes.DefaultTheme
    )
Creates a Gallery in `html_out` from already recorded examples in `folder`.
"""
function gallery_from_recordings(
        folder::String,
        html_out::String = abspath(joinpath(pathof(MakieGallery), "..", "..", "index.html"));
        tags = [
            string.(AbstractPlotting.atomic_function_symbols)...,
            "interaction",
            "record",
            "statsmakie",
            "vbox",
            "layout",
            "legend",
            "colorlegend",
            "vectorfield",
            "poly",
            "camera",
            "recipe",
            "theme",
            "annotations",
            "layout",
            "grid",
            "geomakie"
        ],
        hltheme = Highlights.Themes.DefaultTheme,
        print_toplevel = true
    )

    items = map(MakieGallery.database) do example
        base_path = joinpath(folder, string(example.unique_name))
        media_path = joinpath(base_path, "media")
        media = master_url.(abspath(folder), joinpath.(abspath(media_path), readdir(media_path)))
        mdpath = joinpath(base_path, "index.md")
        save_highlighted_markdown(mdpath, example, media, hltheme; print_toplevel = print_toplevel)
        md2html(mdpath; stylesheets = [relpath(joinpath(dirname(html_out), "syntaxtheme.css"), base_path)])
        MediaItem(base_path, example)
    end
    open(html_out, "w") do io
        println(io, create_page(items, tags))
    end
    open(joinpath(dirname(html_out), "syntaxtheme.css"), "w") do io
        stylesheet(io, MIME("text/css"), hltheme)
    end
end

function rescale_image(path::AbstractString, target_path::AbstractString, sz::Int = 200)
    !isfile(path) && error("Input argument must be a file!")
    img = FileIO.load(path)
    # calculate new image size `newsz`
    (height, width) = size(img)
    (scale_height, scale_width) = sz ./ (height, width)
    scale = min(scale_height, scale_width)
    newsz = round.(Int, (height, width) .* scale)

    # filter image + resize image
    gaussfactor = 0.4
    σ = map((o,n) -> gaussfactor*o/n, size(img), newsz)
    kern = KernelFactors.gaussian(σ)   # from ImageFiltering
    imgf = ImageFiltering.imfilter(img, kern, NA())
    newimg = ImageTransformations.imresize(imgf, newsz)
    # save image
    FileIO.save(target_path, newimg)
end

"""
    generate_thumbnail(path::AbstractString, target_path, thumb_size::Int = 200)

Generates a (proportionally-scaled) thumbnail with maximum side dimension `sz`.
`sz` must be an integer, and the default value is 200 pixels.
"""
function generate_thumbnail(path, thumb_path, thumb_size = 128)
    if any(ext-> endswith(path, ext), (".png", ".jpeg", ".jpg"))
        rescale_image(path, thumb_path, thumb_size)
    elseif any(ext-> endswith(path, ext), (".gif", ".mp4", ".webm"))
        seektime = get_video_duration(path) / 2
        FFMPEG.ffmpeg_exe(`-loglevel quiet -ss $seektime -i $path -vframes 1 -vf "scale=$(thumb_size):-2" -y -f image2 $thumb_path`)
    else
        @warn("Unsupported return file format in $path")
    end
end

function generate_thumbnails(media_root)
    for folder in readdir(media_root)
        media = joinpath(media_root, folder, "media")
        if !isfile(media) && ispath(media)
            isempty(readdir(media)) && error("Media $(media) doesn't contain anything")
            sample = joinpath(media, first(readdir(media)))
            generate_thumbnail(sample, joinpath(media, "thumb.jpg"))
        end
    end
end

"""
Embedds all produced media in one big html file
"""
function generate_preview(media_root, path = joinpath(@__DIR__, "preview.html"))
    open(path, "w") do io
        for folder in readdir(media_root)
            media = joinpath(media_root, folder, "media")
            if !isfile(media) && ispath(media)
                medias = joinpath.(media, readdir(media))
                println(io, "
 $folder
")
                MakieGallery.embed_media(io, medias)
            end
        end
    end
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Write

import ...LibGit2, ..Cache, ..Read, ...Pkg.PkgError
importall ...LibGit2

function prefetch(pkg::AbstractString, sha1::AbstractString)
    isempty(Cache.prefetch(pkg, Read.url(pkg), sha1)) && return
    throw(PkgError("$pkg: couldn't find commit $(sha1[1:10])"))
end

function fetch(repo::GitRepo, pkg::AbstractString, sha1::AbstractString)
    cache = Cache.path(pkg)
    LibGit2.fetch(repo, remoteurl=cache, refspecs=["+refs/*:refs/remotes/cache/*"])
    LibGit2.need_update(repo)
    LibGit2.iscommit(sha1, repo) && return
    f = with(GitRepo, cache) do repo
         LibGit2.iscommit(sha1, repo)
    end ? "fetch" : "prefetch"
    url = Read.issue_url(pkg)
    if isempty(url)
        throw(PkgError("$pkg: $f failed to get commit $(sha1[1:10]), please file a bug report with the package author."))
    else
        throw(PkgError("$pkg: $f failed to get commit $(sha1[1:10]), please file an issue at $url"))
    end
end

function checkout(repo::GitRepo, pkg::AbstractString, sha1::AbstractString)
    LibGit2.set_remote_url(repo, "origin", Cache.normalize_url(Read.url(pkg)))
    LibGit2.checkout!(repo, sha1)
end

function install(pkg::AbstractString, sha1::AbstractString)
    prefetch(pkg, sha1)
    repo = if isdir(".trash/$pkg")
        mv(".trash/$pkg", "./$pkg") #TODO check for newer version in cache before moving
        GitRepo(pkg)
    else
        LibGit2.clone(Cache.path(pkg), pkg)
    end
    try
        fetch(repo, pkg, sha1)
        checkout(repo, pkg, sha1)
    finally
        close(repo)
    end
end

function update(pkg::AbstractString, sha1::AbstractString)
    prefetch(pkg, sha1)
    with(GitRepo, pkg) do repo
        fetch(repo, pkg, sha1)
        checkout(repo, pkg, sha1)
    end
end

function remove(pkg::AbstractString)
    isdir(".trash") || mkdir(".trash")
    ispath(".trash/$pkg") && rm(".trash/$pkg", recursive=true)
    mv(pkg, ".trash/$pkg")
end

end # module

max_stars_repo_name: caodi/julia

max_issues_repo_name: caodi/julia

max_forks_repo_name: caodi/julia


text:
module Jultran
export compute_κs

# look up table for gass' hitran codes, molecular weights, partition function
gases = {
    :h2o => [1, 18, T -> T^1.5],
    :co2 => [2, 44, T -> T],
    :ch4 => [6, 16, T -> T^1.5]
}

# contants
N = 6.022136736e23
h = 6.626075540e-34
c = 2.99792458e8
k = 1.38065812e-23

# compute absorption coefficients (kg/m^2 per cm^-1) directly from HITRAN lines
function compute_κs(gas, νs::Range, T, p)
    num_widths = 1e3 # number of multiples of line width (γ) before truncation

    gas_array = gases[gas]
    file_name = "$(Pkg.dir())/Jultran/data/$(@sprintf "%02d" gas_array[1])_hit04.par"
    ν_grid_min, ν_grid_step, ν_grid_max = νs[1], step(νs), νs[end]
    length_νs = length(νs)
    κs = zeros(Float64, length_νs)

    Q = gas_array[3] # partition function
    T_fact = (.1N/gas_array[2]) * (Q(296.)/Q(T))
    for line in eachline(open(file_name))
        # check if line is nonempty and is first isotopologue
        if (length(line) > 0) && (line[3] == '1')
            ν = float(line[4:15]) # wavenumber
            # air-broadened line width, with temperature dependence:
            γ = float(line[36:40]) * (p/1.013e5) * (296/T)^float(line[56:59])

            ν_min, ν_max  = ν + [-1,1] * num_widths * γ
            if (ν_min <= ν_grid_max) && (ν_max >= ν_grid_min)
                i_min = max(1, 1 + iceil((ν_min-ν_grid_min)/ν_grid_step))
                i_max = min(length_νs, 1 + ifloor((ν_max-ν_grid_min)/ν_grid_step))

                # intensity (S), converted to (m^2/kg per cm^-1),
                # and multiplied by two factors expressing T-dependence,
                # and the partition function:
                f = -100h * c/k
                S = float(line[16:25]) * T_fact * exp(f*float(line[46:55])*(1/T-1/296)) * expm1(f*ν/T) / expm1(f*ν/296)

                κs[i_min:i_max] += S * γ / (π * ((ν-νs[i_min:i_max]).^2 + γ^2))
            end
        end
    end

    κs
end

# iw = 1+int(length_νs*(ν-ν_grid_min)/(ν_grid_max-ν_grid_min))
# nsum = int(num_widths*γ/ν_grid_step)
# i_min = max(1,iw-nsum)
# i_max = min(length_νs,iw+nsum)
#
# if (i_max > 1) && (i_min < length_νs)
#     νs = ((i_min-iw):(i_max-iw))*ν_grid_step
#     println(i_min)
#     f = -100h * c/k
#     S = float(line[16:25]) * T_fact * exp(f*float(line[46:55])*(1/T-1/296)) * expm1(f*ν/T) / expm1(f*ν/296)
#     κs[i_min:i_max] += S * γ / (π * (νs.^2 + γ^2))
# end

end # module


text:
# ---
# title: 100. Same Tree
# id: problem100
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/same-tree/description/>
# hidden: true
# ---
#
# Given two binary trees, write a function to check if they are the same or not.
#
# Two binary trees are considered the same if they are structurally identical
# and the nodes have the same value.
#
# **Example 1:**
#
#
#
#     Input:     1         1
#               / \       / \
#              2   3     2   3
#
#             [1,2,3],   [1,2,3]
#
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input:     1         1
#               /           \
#              2             2
#
#             [1,2],     [1,null,2]
#
#     Output: false
#
#
# **Example 3:**
#
#
#
#     Input:     1         1
#               / \       / \
#              2   1     1   2
#
#             [1,2,1],   [1,1,2]
#
#     Output: false
#
#
#
## @lc code=start
using LeetCode

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


text:
"""
A basic cost function, where the computed cost is the size
(number of children) of the current expression.
"""
astsize(n) = 1
astsize(n::Expr) = 1 + length(n.args) - (iscall(n) ? 1 : 0)

"""
A basic cost function, where the computed cost is the size
(number of children) of the current expression, times -1.
Strives to get the largest expression
"""
astsize_inv(n) = -astsize(n)
astsize_inv(n::Expr) = -1 * astsize(n)

const CostData = Dict{Int64, Tuple{Any, Number}}

"""
An [`AbstractAnalysis`](@ref) that computes the cost of expression nodes
and chooses the node with the smallest cost for each E-Class.
"""
struct ExtractionAnalysis <: AbstractAnalysis
    egraph::EGraph
    costfun::Function
    data::CostData
end

ExtractionAnalysis(g::EGraph, costfun::Function) =
    ExtractionAnalysis(g, costfun, CostData())

make(a::ExtractionAnalysis, n) = (n, a.costfun(n))

function make(an::ExtractionAnalysis, n::Expr)
    ncost = an.costfun(n)

    for child_eclass ∈ getfunargs(n)
        !haskey(an, child_eclass) && return (n, Inf)
        if haskey(an, child_eclass) && an[child_eclass] != nothing
            ncost += last(an[child_eclass])
        end
    end

    return (n, ncost)
end

function join(analysis::ExtractionAnalysis, from, to)
    last(from) <= last(to) ? from : to
end

modify!(analysis::ExtractionAnalysis, id::Int64) = nothing

Base.setindex!(an::ExtractionAnalysis, value, id::Int64) =
    setindex!(an.data, value, id)
Base.getindex(an::ExtractionAnalysis, id::Int64) = an.data[id]
Base.haskey(an::ExtractionAnalysis, id::Int64) = haskey(an.data, id)
Base.delete!(an::ExtractionAnalysis, id::Int64) = delete!(an.data, id)
islazy(an::ExtractionAnalysis) = true

function rec_extract(G::EGraph, an::ExtractionAnalysis, id::Int64)
    (cn, ck) = an[id]
    (!(cn isa Expr) || ck == Inf) && return cn

    expr = copy(cn)
    setfunargs!(expr, getfunargs(expr) .|> a -> rec_extract(G, an, a.id))
    return expr
end

"""
Given an [`ExtractionAnalysis`](@ref), extract the expression
with the smallest computed cost from an [`EGraph`](@ref)
"""
function extract!(G::EGraph, extran::ExtractionAnalysis)
    islazy(extran) && analyze!(G, extran, G.root)
    !(extran ∈ G.analyses) && error("Extraction analysis is not associated to EGraph")
    rec_extract(G, extran, G.root)
end

macro extract(expr, theory, costfun)
    quote
        let g = EGraph($expr)
            saturate!(g, $theory)
            extran = addanalysis!(g, ExtractionAnalysis, $costfun)
            ex = extract!(g, extran)
            (g, ex)
        end
    end |> esc
end


text:
using Test, YaoBlocks
import YaoBlocks: print_annotation

g1 = kron(2, 1 => (im * X))
g2 = control(2, 1, 2 => (im * X))
g3 = chain(im * X, X, Y)

io = IOBuffer()
print_annotation(io, g1, g1, g1[1])
@test String(take!(io)) == "1=>[+im] "

print_annotation(io, g2, g2, content(g2))
@test String(take!(io)) == "(2,) [+im] "


text:
using JTetra
using Base.Test

println("==== Test persistency of linear operators ====")

# read tetrahedral mesh
mesh = getTetraMesh("points.txt","tetras.txt")

# scalar and tensor cell properties for mass matrices
a = ones(mesh.nc)
A = [ones(mesh.nc, 3) zeros(mesh.nc, 3)]

# build all matrices
Grad = getNodalGradientMatrix(mesh)
Curl = getCurlMatrix(mesh)
Div  = getDivergenceMatrix(mesh)

L = getEdgeLength(mesh)
F = getFaceArea(mesh)
V = getCellVolume(mesh)

Mn  = getNodalMassMatrix(mesh, a)
Mea = getEdgeMassMatrix(mesh, a)
MeA = getEdgeMassMatrix(mesh, A)
Mfa = getFaceMassMatrix(mesh, a)
MfA = getFaceMassMatrix(mesh, A)

@test Mea == MeA
@test Mfa == MfA

# now all fields should be non-empty
@test !isempty(mesh.Grad)
@test !isempty(mesh.Curl)
@test !isempty(mesh.Div)
@test !isempty(mesh.L)
@test !isempty(mesh.F)
@test !isempty(mesh.V)

clear!(mesh)

# now all three matrices should be empty
@test isempty(mesh.Grad)
@test isempty(mesh.Curl)
@test isempty(mesh.Div)

print("Passed!\n")

println("=== Test getEdgeIntegralOfPolygonalChain ===")

# define one transmitter
trx = [
  -5.0 0.0 0.0
   5.0 0.0 0.0]

Src = getEdgeIntegralOfPolygonalChain(mesh, trx, normalize = false)

# define a receiver coil - closed polygon
rcv1 = [
  -5.0 -5.0 0.0
  -5.0  5.0 0.0
   5.0  5.0 0.0
   5.0 -5.0 0.0
  -5.0 -5.0 0.0]

Obs1 = getEdgeIntegralOfPolygonalChain(mesh, rcv1, normalize = true)

# define another receiver - open polygon
rcv2 = [
  -5.0 -5.0 0.0
  -5.0  5.0 0.0]

Obs2 = getEdgeIntegralOfPolygonalChain(mesh, rcv2, normalize = true)

print("Passed!\n")

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
module Domain

using Printf, LinearAlgebra

include("torus.jl")

end # module Domain

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore


text:
function _fix_outlayer!(v, th)
    last = first(v)
    for i in eachindex(v)
        err = abs(v[i] - last)
        ref = abs((v[i] + last)/ 2) * th
        if (err > ref)
            v[i] = last
        end
        last = v[i]
    end
    v
end

_normalize!(v) = v .= v ./ maximum(v)

function _extract_hist_vecs(h::Dict)
    xs = collect(keys(h))
    sort!(xs)
    ys = [h[x] for x in xs]
    # _fix_outlayer!(ys, 1.0)
    n_ave_conv!(ys, 25)
    _normalize!(ys)
    return xs, ys
end

## ------------------------------------------------------
function _plot_marginal!(plt, P, Vv; pltargs...)
    v_av = sum(P .* Vv)
    vline!(plt, [v_av]; pltargs..., label = "", ls = :dot)

    Pv = hist(Vv, P)
    xs, ys = _extract_hist_vecs(Pv)
    m, M = minimum(ys), maximum(ys)
    ylim = (m, M) .+ ((M - m) .* 0.1 .* [-1, 1])
    plot!(plt, xs, ys; label = "", ylim,
        pltargs...,
    )

    plt
end

## ------------------------------------------------------
function plot_flux_marginals3D(simid, D, ϵs, cg, is_glclim_stst)
    MEmode = "ME_AVZ_EXPECTED_AVUG_BOUNDED"

    # Plots
    dyn_pz, dyn_pug, dyn_puo = plot(), plot(), plot()
    me_pz, me_pug, me_puo = plot(), plot(), plot()
    ps = [dyn_pz, dyn_pug, dyn_puo, me_pz, me_pug, me_puo]
    Vz, Vug, Vuo = nothing, nothing, nothing
    for (ϵi, ϵ) in enumerate(ϵs)

        simparams = (;D, ϵ, cg)

        # Check and load
        status = Dyn.get_status(simid, simparams)
        if (status != STST_SIM_STATUS)
            @warn("Wrong status ", status)
            continue
        end

        datfile = Dyn.maxent_file(simid, MEmode, simparams)
        if !isfile(datfile)
            @warn("File missing ", datfile)
            continue
        end

        if !is_glclim_stst[simparams]
            @warn("Not glc-lim STST ", status)
            continue
        end

        dat = Dyn.load_maxent(simid, MEmode, simparams)
        PME = dat.PME

        S = Dyn.load_simdat(simid, simparams)
        PX = S.P
        V = S.V
        Vz = Vi(V, :z)
        Vug = Vi(V, :ug)
        Vuo = Vi(V, :uo)

        z_avPME = sum(PME .* Vz)
        if (abs(z_avPME - D) / D) > 0.05
            @warn("D and z not converged", status, z_avPME, D)
            continue
        end

        max_ϵs_ = maximum(ϵs)
        ls0 = 3

        for (P, pz, pug, puo) in [
                (PX, dyn_pz, dyn_pug, dyn_puo),
                (PME, me_pz, me_pug, me_puo)
            ]
            comparams = (;
                ylim = [0.0, 1.2],
                ylabel = _textbf("non-normalized pdf"),
                lw =  ls0 + (2 * ls0) * (ϵ / max_ϵs_),
                color = :black,
                alpha = 0.7
            )
            _plot_marginal!(pz, P, Vz;
                xlabel = _textbf("z"),
                comparams...
            )

            _plot_marginal!(pug, P, Vug;
                xlabel = _textbf("u_g"),
                comparams...
            )

            _plot_marginal!(puo, P, Vuo;
                xlabel = _textbf("u_o"),
                comparams...
            )
        end
    end

    return ps

end


text:
# Finite Horizon algorithm evaluating problem with use of other POMDP algorithms
# User has to implement FiniteHorizon function (stage_states, stage_actions and stage_stateindex) in such a way that the function returns current epoch and the following one (the one that has already been evaluated)
# For example see FiniteHorizonPOMDP/test/instances/1DFiniteHorizonGridWorld.jl
#
# Currently supporting: ValueIterationSolver
#
# Possible future improvements: Change arrays from row oriented to column oriented
#                               Create fixhorizon(m::Union{MDP,POMDP}, T::Int) utility

# Policy struct created according to one used in DiscreteValueIteration
mutable struct FiniteHorizonValuePolicy{Q<:AbstractArray, U<:AbstractMatrix, P<:AbstractMatrix, A, M<:MDP} <: Policy
    qmat::Q
    util::U
    policy::P
    action_map::Vector{A}
    include_Q::Bool
    mdp::M
end

# Policy constructor
function FiniteHorizonValuePolicy(mdp::MDP)
    return FiniteHorizonValuePolicy(zeros(horizon(mdp), no_states(mdp), length(actions(mdp))), zeros(horizon(mdp), no_states(mdp)), ones(Int64, horizon(mdp), no_states(mdp)), ordered_actions(mdp), true, mdp)
end

function no_states(mdp::MDP)
    length(stage_states(mdp,1))
end

function action(policy::FiniteHorizonValuePolicy, s::S) where S
    sidx = stage_stateindex(policy.mdp, s, epoch(s))
    aidx = policy.policy'[sidx]
    return policy.action_map[aidx]
end

# Method stores record for each evaluated epoch to FiniteHorizonPolicy and returns it
function addepochrecord(fhpolicy::FiniteHorizonValuePolicy, qmat, util, policy)
    global fhepoch
    fhpolicy.qmat[fhepoch, :, :] = qmat
    fhpolicy.util[fhepoch, :] = util
    fhpolicy.policy[fhepoch, :] = policy
    return fhpolicy
end

# Global variable for storing number of epoch in order to pass it to functions from outer solvers
# Is there a better way to achieve this?
fhepoch = -1

# Is it possible to change name of this function to solve?
# The problem is that I am not able to use for example DiscreteValueIteration.solve() as I do not know the Solver in advance

# MDP given horizon 5 assumes that agent can move 4 times
function solve(mdp::MDP; verbose::Bool=false, new_VI::Bool=true)
    fhpolicy = FiniteHorizonValuePolicy(mdp)
    util = fill(0., no_states(mdp)) # XXX not all MDPs have a no_states field. Suggest using length(states(mdp))

    for epoch=horizon(mdp)-1:-1:1
        # Store number of epoch to global variable in order to work properly
        # Is there a better way to achieve this?
        global fhepoch = epoch

        if verbose
            println("EPOCH: $epoch")
        end

        stage_q, util, pol = valueiterationsolver(mdp, epoch, util)

        fhpolicy = addepochrecord(fhpolicy, stage_q, util, pol)

        if verbose
            println("POLICY\n")
            println("QMAT")
            println(stage_q)
            println("util")
            println(util)
            println("policy")
            println(policy)
            # println("action_map")
            # println(policy.action_map)
            # println("include_Q")
            # println(policy.include_Q)
            # println("mdp")
            # println(policy.mdp)
            println("\n\n\n")
        end
    end

    return fhpolicy
end


text:
struct ColumnBlockComponentCostModel{Tv, α_Col, β_Col} <: AbstractConnectivityModel
    α_col::α_Col
    β_col::β_Col
end

function ColumnBlockComponentCostModel{Tv}(; α_col=false, β_col=false) where {Tv}
    ColumnBlockComponentCostModel{Tv}(α_col, β_col)
end

@deprecate ColumnBlockComponentCostModel{Tv}(w_max, α_col, β_col) where {Tv, Ti} ConstrainedCost(ColumnBlockComponentCostModel{Tv}(α_col, β_col), VertexCount(), w_max)
function ColumnBlockComponentCostModel{Tv}(α_col::α_Col, β_col::β_Col) where {Tv, α_Col, β_Col}
    return ColumnBlockComponentCostModel{Tv, α_Col, β_Col}(α_col, β_col)
end

@inline cost_type(::Type{<:ColumnBlockComponentCostModel{Tv}}) where {Tv} = Tv

(mdl::ColumnBlockComponentCostModel{Tv, α_Col, β_Col})(n_vertices, n_pins, n_nets) where {Tv, α_Col, β_Col} = block_component(mdl.α_col, n_vertices) + n_nets * block_component(mdl.β_col, n_vertices)

struct BlockComponentCostModel{Tv, R, α_Row, α_Col, β_Row<:Tuple{Vararg{Any, R}}, β_Col<:Tuple{Vararg{Any, R}}}
    α_row::α_Row
    α_col::α_Col
    β_row::β_Row
    β_col::β_Col
end

function BlockComponentCostModel{Tv}(; α_row=false, α_col=false, β_row = (), β_col = ()) where {Tv}
    BlockComponentCostModel{Tv}(α_row, α_col, β_row, β_col)
end

function Base.permutedims(cst::BlockComponentCostModel{Tv}) where {Tv}
    return BlockComponentCostModel{Tv}(cst.α_col, cst.α_row, cst.β_col, cst.β_row)
end

@deprecate BlockComponentCostModel{Tv}(w_max, U, α_row, α_col, β_row, β_col) where {Tv, Ti} ConstrainedCost(BlockComponentCostModel{Tv}(α_row, α_col, β_row, β_col), VertexCount(), w_max)
function BlockComponentCostModel{Tv}(α_row::α_Row, α_col::α_Col, β_row::β_Row, β_col::β_Col) where {Tv, R, α_Row, α_Col, β_Row<:Tuple{Vararg{Any, R}}, β_Col<:Tuple{Vararg{Any, R}}}
    return BlockComponentCostModel{Tv, R, α_Row, α_Col, β_Row, β_Col}(α_row, α_col, β_row, β_col)
end

@inline cost_type(::Type{<:BlockComponentCostModel{Tv}}) where {Tv} = Tv

@propagate_inbounds block_component(f::F, w) where {F} = f(w)
@propagate_inbounds block_component(f::Number, w) = f
@propagate_inbounds block_component(f::Tuple, w) = f[w]
@propagate_inbounds block_component(f::AbstractArray, w) = f[w]

mutable struct BlockComponentCostStepOracle{Tv, Ti, Tc, R, Mdl<:BlockComponentCostModel{Tc, R}} <: AbstractOracleCost{Mdl}
    A::SparseMatrixCSC{Tv, Ti}
    Π_asg::MapPartition{Ti}
    Π_spl::SplitPartition{Ti}
    mdl::Mdl
    hst::Vector{Ti}
    Δ::Matrix{Tc}
    d::Vector{Tc}
    j::Ti
    j′::Ti
end

function oracle_stripe(hint::AbstractHint, mdl::BlockComponentCostModel{Tc, R}, A::SparseMatrixCSC{Tv, Ti}, Π; kwargs...) where {Tv, Ti, Tc, R}
    @inbounds begin
        m, n = size(A)
        K = length(Π)
        return BlockComponentCostStepOracle(A, convert(MapPartition, Π), convert(SplitPartition, Π), mdl, ones(Ti, K), undefs(Tc, R, n + 1), zeros(Tc, R), Ti(1), Ti(1))
    end
end

oracle_model(ocl::BlockComponentCostStepOracle) = ocl.mdl

@inline function (ocl::BlockComponentCostStepOracle{Tv, Ti, Tc, R})(j::Ti, j′::Ti, k...) where {Tv, Ti, Tc, R}
    @inbounds begin
        ocl_j = ocl.j
        ocl_j′ = ocl.j′
        A = ocl.A
        Π_asg = ocl.Π_asg.asg
        Π_spl = ocl.Π_spl.spl
        f = ocl.mdl
        d = ocl.d
        Δ = ocl.Δ
        pos = A.colptr
        idx = A.rowval
        hst = ocl.hst

        if j′ < ocl_j′
            ocl_j = Ti(1)
            ocl_j′ = Ti(1)
            zero!(d)
            one!(ocl.hst)
        end
        while ocl_j′ < j′
            q = pos[ocl_j′]
            q′ = pos[ocl_j′ + 1]

            for r = 1:R
                Δ[r, ocl_j′ + 1] = zero(Tc)
            end

            for _q = q:q′ - 1
                i = idx[_q]
                k = Π_asg[i]
                j₀ = hst[k] - 1
                u = Π_spl[k + 1] - Π_spl[k]
                if j₀  < ocl_j′
                    for r = 1:R
                        Δ[r, j₀ + 1] -= block_component(f.β_row[r], u)
                    end
                    for r = 1:R
                        Δ[r, ocl_j′ + 1] += block_component(f.β_row[r], u)
                    end
                end
                if j₀ < ocl_j
                    for r = 1:R
                        d[r] += block_component(f.β_row[r], u)
                    end
                end
                hst[k] = ocl_j′ + 1
            end

            ocl_j′ += 1
        end
        while j < ocl_j
            ocl_j -= 1
            for r = 1:R
                d[r] += Δ[r, ocl_j + 1]
            end
        end
        while j > ocl_j
            for r = 1:R
                d[r] -= Δ[r, ocl_j + 1]
            end
            ocl_j += 1
        end

        w = j′ - j
        c = block_component(f.α_col, w)
        for r = 1:R
            c += d[r] * block_component(f.β_col[r], w)
        end

        ocl.j = ocl_j
        ocl.j′ = ocl_j′
        return c
    end
end

function row_component_value(Π, mdl::BlockComponentCostModel)
    Π = convert(SplitPartition, Π)
    c_α = zero(cost_type(mdl))
    for k = 1:Π.K
        u = Π.spl[k + 1] - Π.spl[k]
        c_α += block_component(mdl.α_row, u)
    end
    return c_α
end

#=
function pack_stripe(A::SparseMatrixCSC{Tv, Ti}, method::DynamicTotalChunker{F}, args...; kwargs...) where {F<:BlockComponentCostModel, Tv, Ti}
    return pack_stripe(A, DynamicTotalChunker(ConstrainedCost(method.f, FeasibleCost(), Feasible())), args..., kwargs...)
end

function pack_stripe(A::SparseMatrixCSC{Tv, Ti}, method::DynamicTotalChunker{<:ConstrainedCost{F}}, Π; kwargs...) where {Tc, R, F<:BlockComponentCostModel{Tc, R}, Tv, Ti}
    @inbounds begin
        m, n = size(A)
        A_pos = A.colptr
        A_idx = A.rowval

        f = method.f.f
        w = oracle_stripe(StepHint(), method.f.w, A, Π) #TODO reverse outer loop so step oracles can have stationary j′
        w_max = method.f.w_max
        Π_asg = convert(MapPartition, Π).asg
        Π_spl = convert(DomainPartition, Π).spl
        K = length(Π)
        hst = fill(n + 1, K)
        Δ = zeros(Tc, R, n + 1)
        d = zeros(Tc, R)
        cst = Vector{cost_type(f)}(undef, n + 1) # cst[j] is the best cost of a partition from j to n
        spl = Vector{Int}(undef, n + 1)

        for r = 1:R
            Δ[r, n + 1] = zero(Tc)
        end
        cst[n + 1] = zero(Tc)
        for j = n:-1:1
            for q = A_pos[j] : A_pos[j + 1] - 1
                i = A_idx[q]
                k = Π_asg[i]
                u = Π_spl[k + 1] - Π_spl[k]
                if hst[k] > j
                    for r = 1:R
                        Δ[r, hst[k]] -= block_component(f.β_row[r], u)
                    end
                    for r = 1:R
                        Δ[r, j] += block_component(f.β_row[r], u)
                    end
                end
                hst[k] = j
            end
            zero!(d)
            for r = 1:R
                d[r] += Δ[r, j]
            end
            best_c = cst[j + 1]
            for r = 1:R
                best_c += d[r] * block_component(f.β_col[r], 1)
            end
            best_c += block_component(f.α_col, 1)
            best_j′ = j + 1
            for j′ = j + 2: n + 1
                if w(j, j′) > w_max
                    break
                end
                for r = 1:R
                    d[r] += Δ[r, j′ - 1]
                end
                c = cst[j′]
                for r = 1:R
                    c += d[r] * block_component(f.β_col[r], j′ - j)
                end
                c += block_component(f.α_col, j′ - j)
                if c < best_c
                    best_c = c
                    best_j′ = j′
                end
            end
            cst[j] = best_c
            spl[j] = best_j′
        end

        K = 0
        j = 1
        while j != n + 1
            j′ = spl[j]
            K += 1
            spl[K] = j
            j = j′
        end
        spl[K + 1] = j
        resize!(spl, K + 1)
        return SplitPartition{Ti}(K, spl)
    end
end
=#