text:
safe_parseint32(s) = (s != nothing) ? Base.parse(Int32, s) : nothing
safe_parseint64(s) = (s != nothing) ? Base.parse(Int64, s) : nothing
safe_parseint(s) = (s != nothing) ? Base.parse(Int, s) : nothing
safe_parsebool(s) = (s != nothing) ? ((lowercase(s) == "true") ? true : false) : nothing
safe_parseb64(s) = (s != nothing) ? base64decode(s) : nothing

function safe_parse_as(as::Type, s::Union{String, Void})
    if (as == String) || (s == nothing)
        return s
    elseif (as == Int)
        return Base.parse(Int, s)
    elseif (as == Int32)
        return int32(s)
    elseif (as == Float64)
        return float64(s)
    elseif (as == Int64)
        return int64(s)
    elseif (as == Bool)
        return (lowercase(s) == "true" ? true : false)
    elseif (as == DateTime)
        return DateTime(s[1:end-1])
    end
    error("Unsupported parse type $as")
end
export safe_parse_as

macro parse_vector(typ, vect)
     quote
        jl_vect = $(esc(typ))[]
        if ($(esc(vect)) != nothing)
            for pd in $(esc(vect))
            	push!(jl_vect, $(typ)(pd))
            end
        end
        jl_vect
     end
end
export @parse_vector

function parse_vector_as(as_type::Type, typ_str::String, vect)
    jl_vect = as_type[]
    if (vect == nothing) return jl_vect end
    for pd in vect
        val = LightXML.content(LightXML.find_element(pd, typ_str))
        val = safe_parse_as(as_type, val)
        if (val == nothing) error("Invalid $(typ_str) for pd vector") end
        push!(jl_vect, val)
    end

    jl_vect
end
export parse_vector_as

function parse_calendar_time(pd, elem::String, format::String)
    datestr = LightXML.content(LightXML.find_element(pd, elem))
    DateTime(datestr[1:end-1], format)
end
parse_calendar_time(pd, elem::String) = parse_calendar_time(pd, elem, "yyyy-MM-DD'T'HH:mm:ss")

xml(o::Any) = string(o)
xml(tag::String, value::Any) = (value != nothing) ? "<$(tag)>" * xml(value) * "</$(tag)>" : ""

function xml(tag::String, children::Array; xmlns="", xsi_type="")
    if (xsi_type != "")
        open_tag = "<$(tag)$(xmlns) xsi:type=\"$(xsi_type)\">"
    else
        open_tag = "<$(tag)$(xmlns)>"
    end

    children_xml = ""
    for child in children
        if isa(child, Tuple)
            (child_tag, child_value) = child
            children_xml = children_xml * xml(child_tag, child_value)
        else
            children_xml = children_xml * xml(child)
        end
    end

    return open_tag * children_xml * "</$(tag)>"
end

function aws_string(dt::DateTime)
    y,m,d = Dates.yearmonthday(dt)
    h,mi,s = Dates.hour(dt),Dates.minute(dt),Dates.second(dt)
    yy = y < 0 ? @sprintf("%05i",y) : lpad(y,4,"0")
    mm = lpad(m,2,"0")
    dd = lpad(d,2,"0")
    hh = lpad(h,2,"0")
    mii = lpad(mi,2,"0")
    ss = lpad(s,2,"0")
    return "$yy-$mm-$(dd)T$hh:$mii:$ss"
end

aws_string(v::Bool) = v ? "true" : "false"
aws_string(v::Vector{UInt8}) = join(map(x->@sprintf("%%%02d", x), v), "")
aws_string(v::Any) = string(v)
export aws_string

#ISO8601
function get_utc_timestamp(addsecs=0;basic=false)
    dt = Dates.unix2datetime(Dates.datetime2unix(now(Dates.UTC)) + addsecs)
    dstr = aws_string(dt)
    if basic
        dstr = replace(dstr, Set(":-"), "")
    end
    return string(dstr, "Z")
end
export get_utc_timestamp


text:
# One stop shop to install SimpleWorld in DEVELOPMENT mode. Can be
# run from the command line:
# julia install.jl

using Pkg
LOC = "https://github.com/scheinerman/SimpleWorld.jl.git"
Pkg.develop(PackageSpec(url = LOC))

ext: jl


text:
module GeometricFlux
using Requires

using Core.Intrinsics: llvmcall
using Base.Threads
using Statistics: mean
using Flux
using Flux: param, glorot_uniform, leakyrelu, GRUCell
using Flux: @treelike
using SparseArrays: SparseMatrixCSC
using LinearAlgebra: I, issymmetric, diagm, eigmax

import Base.Threads: atomictypes, llvmtypes, inttype, ArithmeticTypes, FloatTypes,
       atomic_cas!,
       atomic_xchg!,
       atomic_add!, atomic_sub!,
       atomic_and!, atomic_nand!, atomic_or!, atomic_xor!,
       atomic_max!, atomic_min!

import Base.Sys: ARCH, WORD_SIZE

export

    # layers/msgpass
    MessagePassing,
    neighboring,

    # layers/conv
    GCNConv,
    ChebConv,
    GraphConv,
    GATConv,
    GatedGraphConv,
    EdgeConv,
    message,
    update,
    propagate,

    # layer/pool
    sumpool,
    subpool,
    prodpool,
    divpool,
    maxpool,
    minpool,
    meanpool,
    pool,

    # models
    GAE,
    VGAE,
    InnerProductDecoder,
    VariationalEncoder,

    # linalg
    degree_matrix,
    laplacian_matrix,
    normalized_laplacian,
    neighbors,

    # scatter
    scatter_add!,
    scatter_sub!,
    scatter_max!,
    scatter_min!,
    scatter_mul!,
    scatter_div!,
    scatter_mean!,
    scatter!,

    # graph/utils
    adjlist

include("scatter.jl")
include("layers/msgpass.jl")
include("layers/conv.jl")
include("layers/pool.jl")
include("models.jl")
include("linalg.jl")

function __init__()
    @require LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d" begin
        include("graph/simplegraphs.jl")
    end
    @require SimpleWeightedGraphs = "47aef6b3-ad0c-573a-a1e2-d07658019622" begin
        include("graph/weightedgraphs.jl")
        include("graph/utils.jl")
    end
    @require MetaGraphs = "626554b9-1ddb-594c-aa3c-2596fe9399a5" include("metagraphs.jl")
end

end


text:
# __precompile__()
module EchoJulia

using Reexport

export transform

@reexport using SimradRaw
@reexport using SimradEK60
@reexport using EchogramPyPlot
#@reexport using EchogramPlots
@reexport using EchogramImages
@reexport using EchoviewEvr
@reexport using EchoviewEcs
@reexport using SimradEK60TestData

function between(p, starttime, endtime)

    if starttime == nothing
        starttime = p.filetime
    end

    if endtime == nothing
        endtime = p.filetime
    end

    starttime <= p.filetime <= endtime

end

function transform(datagrams::Vector{SimradRaw.Datagram}; calibration=nothing, starttime=nothing, endtime=nothing)
    ps = collect(pings(datagrams))
    transform(ps, calibration=calibration, starttime=starttime, endtime=endtime)
end

function transform(ps::Vector{SimradEK60.EK60Ping}; calibration=nothing, starttime=nothing, endtime=nothing)

    transducers= []
    if calibration != nothing
        transducers = calibration
    end

    ps = [p for p in ps if between(p, starttime, endtime)]

    frequencies = unique([p.frequency for p in ps])

    dict = Dict()

    index = 1
    for f in frequencies
        fr = trunc(Int,f / 1000)
        psf = [p for p in ps if p.frequency == f]

        gain = nothing
        equivalentbeamangle = nothing
        soundvelocity =  nothing
        absorptioncoefficient = nothing
        transmitpower = nothing
        pulselength = nothing
        sacorrection = nothing

        if length(transducers) >= index
            transducer = transducers[index]
            gain = get(transducer, "Ek60TransducerGain", nothing)
            equivalentbeamangle = get(transducer,"TwoWayBeamAngle", nothing)
            soundvelocity=  get(transducer,"SoundSpeed", nothing)
            absorptioncoefficient= get(transducer,"AbsorptionCoefficient", nothing)
            transmitpower= get(transducer,"TransmittedPower", nothing)
            pulselength= get(transducer,"TransmittedPulseLength", nothing)
            sacorrection= get(transducer,"EK60SaCorrection", nothing)
            if pulselength != nothing
                pulselength = pulselength / 1000
            end
        end

        dict["Sv$fr"] = SimradEK60.Sv(psf,
                                      gain=gain,
                                      equivalentbeamangle=equivalentbeamangle,
                                      soundvelocity=soundvelocity,
                                      absorptioncoefficient=absorptioncoefficient,
                                      transmitpower=transmitpower,
                                      pulselength=pulselength,
                                      sacorrection=sacorrection);

        dict["ntheta$fr"] = alongshipangle(psf)
        dict["nphi$fr"] = athwartshipangle(psf)
        dict["r$fr"] = R(psf, soundvelocity=soundvelocity)
        dict["t$fr"] = SimradEK60.filetime(psf) # timestamps

        if gain != nothing
            dict["gain$fr"] = gain
        else
            dict["gain$fr"] = [x.gain for x in psf]
        end
        if equivalentbeamangle != nothing
            dict["equivalentbeamangle$fr"]= equivalentbeamangle
        else
            dict["equivalentbeamangle$fr"]= [x.equivalentbeamangle for x in psf]
        end
        if soundvelocity != nothing
            dict["soundvelocity$fr"]= soundvelocity
        else
            dict["soundvelocity$fr"]=  [x.soundvelocity for x in psf]
        end
        if absorptioncoefficient != nothing
            dict["absorptioncoefficient$fr"] = absorptioncoefficient
        else
            dict["absorptioncoefficient$fr"] = [x.absorptioncoefficient for x in psf]
        end
        if transmitpower != nothing
            dict["transmitpower$fr"] = transmitpower
        else
            dict["transmitpower$fr"] = [x.transmitpower for x in psf]
        end
        if pulselength != nothing
            dict["pulselength$fr"] = pulselength
        else
            dict["pulselength$fr"] = [x.pulselength for x in psf]
        end
        if sacorrection != nothing
            dict["sacorrection$fr"] = sacorrection
        else
            dict["sacorrection$fr"] =[x.sacorrection for x in psf]
        end

        dict["sampleinterval$fr"] =[x.sampleinterval for x in psf]

        index +=1
    end

    dict["DESCR"] = "Converted by EchoJulia."

    return dict
end

end # module


text:
#!/usr/bin/env julia

# Create image data VTK file.
# This corresponds to a rectilinear grid with uniform spacing in each direction.

using WriteVTK
const FloatType = Float32
const vtk_filename_noext = "imagedata"

function main()
    # Define grid.
    Ni, Nj, Nk = 20, 30, 40

    # Create some scalar and vectorial data.
    p = zeros(FloatType, Ni, Nj, Nk)
    q = zeros(FloatType, Ni, Nj, Nk)
    vec = zeros(FloatType, 3, Ni, Nj, Nk)

    for k = 1:Nk, j = 1:Nj, i = 1:Ni
        p[i, j, k] = i*i + k
        q[i, j, k] = k*sqrt(j)
        vec[1, i, j, k] = i
        vec[2, i, j, k] = j
        vec[3, i, j, k] = k
    end

    # Create some scalar data at grid cells.
    # Note that in structured grids, the cells are the hexahedra formed between
    # grid points.
    cdata = zeros(FloatType, Ni-1, Nj-1, Nk-1)
    for k = 1:Nk-1, j = 1:Nj-1, i = 1:Ni-1
        cdata[i, j, k] = 2i + 3k * sin(3*pi * (j-1) / (Nj-2))
    end

    # These are all optional:
    extent = map(r -> r .+ 43, (1:Ni, 1:Nj, 1:Nk))
    origin = (1.2, 4.3, -3.1)
    spacing = (0.1, 0.5, 1.2)

    # Initialise new vti file (image data).
    @time outfiles = vtk_grid(vtk_filename_noext, Ni, Nj, Nk, extent=extent,
                              origin=origin, spacing=spacing) do vtk
        # Add data.
        vtk["p_values"] = p
        vtk["q_values"] = q
        vtk["myVector"] = vec
        vtk["myCellData"] = cdata

        # Field data
        vtk["time"] = 4.2
        vtk["field_tuple"] =
            ([2.0, 3.0, 4.0], [1.0, 2.0, 3.0])  # 2 components, 3 "tuples"
        vtk["field_matrix", VTKFieldData()] =
            (1:4) * (2:5)'  # 1 component, 16 "tuples"
        vtk["a_string"] = "blabla"
        vtk["some_strings"] = ["blabla", "ab©dé", "ñeee"]
        vtk["string_tuple"] = ("blabla", "ab©dé", "ñeee")
    end

    # Test 2D dataset
    @time outfiles_2D = vtk_grid(vtk_filename_noext * "_2D", Ni, Nj,
                                 # extent=extent[1:2],
                                 origin=origin[1:2], spacing=spacing[1:2]) do vtk
        vtk["p_values"] = p[:, :, 1]
        vtk["myVector"] = vec[:, :, :, 1]
        vtk["myCellData"] = cdata[:, :, 1]
    end

    append!(outfiles, outfiles_2D)

    # Test specifying coordinates using LinRange
    for D in (2, 3)  # test 2D and 3D datasets
        suffix = D == 2 ? "_2D" : ""
        let xyz = (LinRange(0., 5., Ni),
                   LinRange(1., 3., Nj),
                   LinRange(2., 6., Nk))
            coords = xyz[1:D]
            @time saved = vtk_grid(vtk_filename_noext * "_LinRange" * suffix,
                                   coords, compress=true) do vtk
                vtk["myVector"] = D == 3 ? vec : view(vec, 1:2, :, :, 1)
                vtk["myCellData"] = D == 3 ? cdata : view(cdata, :, :, 1)
            end
            append!(outfiles, saved)
        end

        # Similar using StepRangeLen
        let xyz = (0:0.08:5, 1.2:0.05:2., 1.2:0.1:2.)
            coords = xyz[1:D]
            @time saved = vtk_grid(
                identity, vtk_filename_noext * "_StepRangeLen" * suffix, coords)
            append!(outfiles, saved)
        end
    end

    println("Saved:   ", join(outfiles, "  "))

    return outfiles::Vector{String}
end

main()


text:
function utotal(p, data :: Data)

	@unpack N, cutoff, side, l = data

	nt = Threads.nthreads()

	U = zeros(nt)

	npt = div_thread(N, nt)
	first_atom, next_atom = initial_linklist(p, data)

	Threads.@threads for id in 1:nt

		@inbounds @simd for i in npt[id, 1]:npt[id, 2]

			icell = trunc(Int64, p[i][1]/cutoff) + 1
			jcell = trunc(Int64, p[i][2]/cutoff) + 1

			for ic in icell-1:icell+1, jc in jcell-1:jcell+1

				iw, jw = wrap_cell(ic, jc, l)

				j = first_atom[iw, jw]

				while j > 0

					if j > i

						r = rpbc(p[i], p[j], side)

						if r < cutoff
							U[id] += upair(r, data)
						end
					end

					j = next_atom[j]
				end
			end
		end
	end

	return sum(U)
end


text:
@testset "DeviceModel Tests" begin
    @test_throws ArgumentError DeviceModel(ThermalGen, ThermalStandardUnitCommitment)
    @test_throws ArgumentError DeviceModel(ThermalStandard, PSI.AbstractDeviceFormulation)
end

@testset "OperationsProblem Tests" begin
    sys = build_system("c_sys5")
    for p in [true, false]
        t = OperationsProblem(TestOpProblem, CopperPlatePowerModel, sys, use_parameters = p)
        @test PSI.model_has_parameters(t.psi_container) == p
    end
end


text:
using ..metrics

import OrdinaryDiffEq, DiffEqCallbacks, StochasticDiffEq

function recast! end

"""
    integrate(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128},
            state::T, dstate::T, fout::Function; kwargs...)

Integrate using OrdinaryDiffEq
"""
function integrate(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128},
            state::T, dstate::T, fout::Function;
            alg::OrdinaryDiffEq.OrdinaryDiffEqAlgorithm = OrdinaryDiffEq.DP5(),
            steady_state = false, tol = 1e-3, save_everystep = false,
            callback = nothing, kwargs...) where T

    function df_(dx::Vector{Complex128}, x::Vector{Complex128}, p, t)
        recast!(x, state)
        recast!(dx, dstate)
        df(t, state, dstate)
        recast!(dstate, dx)
    end
    function fout_(x::Vector{Complex128}, t::Float64, integrator)
        recast!(x, state)
        fout(t, state)
    end

    out_type = pure_inference(fout, Tuple{eltype(tspan),typeof(state)})

    out = DiffEqCallbacks.SavedValues(Float64,out_type)

    scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan,
                                         save_everystep=save_everystep,
                                         save_start = false)

    prob = OrdinaryDiffEq.ODEProblem{true}(df_, x0,(tspan[1],tspan[end]))

    if steady_state
        affect! = function (integrator)
            !save_everystep && scb.affect!(integrator,true)
            OrdinaryDiffEq.terminate!(integrator)
        end
        _cb = OrdinaryDiffEq.DiscreteCallback(
                                SteadyStateCondtion(copy(state),tol,state),
                                affect!;
                                save_positions = (false,false))
        cb = OrdinaryDiffEq.CallbackSet(_cb,scb)
    else
        cb = scb
    end

    full_cb = OrdinaryDiffEq.CallbackSet(callback,cb)

    sol = OrdinaryDiffEq.solve(
                prob,
                alg;
                reltol = 1.0e-6,
                abstol = 1.0e-8,
                save_everystep = false, save_start = false,
                save_end = false,
                callback=full_cb, kwargs...)
    out.t,out.saveval
end

function integrate(tspan::Vector{Float64}, df::Function, x0::Vector{Complex128},
            state::T, dstate::T, ::Void; kwargs...) where T
    function fout(t::Float64, state::T)
        copy(state)
    end
    integrate(tspan, df, x0, state, dstate, fout; kwargs...)
end

struct SteadyStateCondtion{T,T2,T3}
    rho0::T
    tol::T2
    state::T3
end
function (c::SteadyStateCondtion)(rho,t,integrator)
    timeevolution.recast!(rho,c.state)
    dt = integrator.dt
    drho = metrics.tracedistance(c.rho0, c.state)
    c.rho0.data[:] = c.state.data
    drho/dt < c.tol
end

"""
    integrate_stoch(tspan::Vector{Float64}, df::Function, dg::Vector{Function}, x0::Vector{Complex128},
            state::T, dstate::T, fout::Function; kwargs...)

Integrate using StochasticDiffEq
"""
function integrate_stoch(tspan::Vector{Float64}, df::Function, dg::Function, x0::Vector{Complex128},
            state::T, dstate::T, fout::Function, n::Int;
            save_everystep = false, callback=nothing,
            kwargs...) where T

    function df_(dx::Vector{Complex128}, x::Vector{Complex128}, p, t)
        recast!(x, state)
        recast!(dx, dstate)
        df(t, state, dstate)
        recast!(dstate, dx)
    end

    function fout_(x::Vector{Complex128}, t::Float64, integrator)
        recast!(x, state)
        fout(t, state)
    end

    out_type = pure_inference(fout, Tuple{eltype(tspan),typeof(state)})

    out = DiffEqCallbacks.SavedValues(Float64,out_type)

    scb = DiffEqCallbacks.SavingCallback(fout_,out,saveat=tspan,
                                         save_everystep=save_everystep,
                                         save_start = false)

    full_cb = OrdinaryDiffEq.CallbackSet(callback, scb)
    if n == 1
        integrate_stoch_scalar(tspan, df_, dg, x0, state, dstate, fout_, out, full_cb; kwargs...)
    else
        integrate_stoch_ndiag(tspan, df_, dg, x0, state, dstate, fout_, out, n, full_cb; kwargs...)
    end
end

"""
Dispatch for scalar noise problems.
"""
function integrate_stoch_scalar(tspan::Vector{Float64}, df_::Function, dg::Function,
            x0::Vector{Complex128}, state::T, dstate::T, fout_::Function,
            out::DiffEqCallbacks.SavedValues, callback;
            alg::StochasticDiffEq.StochasticDiffEqAlgorithm = StochasticDiffEq.RKMil(interpretation=:Stratonovich),
            noise=StochasticDiffEq.RealWienerProcess(0.0, randn()),
            kwargs...) where T

    function dg_(dx::Vector{Complex128}, x::Vector{Complex128}, p, t)
        recast!(x, state)
        recast!(dx, dstate)
        dg(t, state, dstate, 1)
        recast!(dstate, dx)
    end

    prob = StochasticDiffEq.SDEProblem{true}(df_, dg_, x0,(tspan[1],tspan[end]),
                noise=noise)

    sol = StochasticDiffEq.solve(
                prob,
                alg;
                reltol = 1.0e-6,
                abstol = 1.0e-8,
                save_everystep = false, save_start = false,
                save_end = false,
                callback=callback, kwargs...)

    out.t,out.saveval
end

"""
Dispatch for non-diagonal noise.
"""
function integrate_stoch_ndiag(tspan::Vector{Float64}, df_::Function, dg::Function,
            x0::Vector{Complex128}, state::T, dstate::T, fout_::Function,
            out::DiffEqCallbacks.SavedValues, n::Int, callback;
            alg::StochasticDiffEq.StochasticDiffEqAlgorithm = StochasticDiffEq.EulerHeun(),
            noise_rate_prototype=Array{Complex128}(length(state), n),
            noise=StochasticDiffEq.RealWienerProcess(0.0, randn(n)),
            kwargs...) where T

    function dg_(dx::Array{Complex128, 2}, x::Vector{Complex128}, p, t)
        recast!(x, state)
        @inbounds for i=1:n
            dx_i = @view dx[:, i]
            recast!(dx_i, dstate)
            dg(t, state, dstate, i)
            recast!(dstate, dx_i)
        end
    end

    prob = StochasticDiffEq.SDEProblem{true}(df_, dg_, x0,(tspan[1],tspan[end]),
                noise=noise,
                noise_rate_prototype=noise_rate_prototype)

    sol = StochasticDiffEq.solve(
                prob,
                alg;
                reltol = 1.0e-6,
                abstol = 1.0e-8,
                save_everystep = false, save_start = false,
                save_end = false,
                callback=callback, kwargs...)

    out.t,out.saveval
end

"""
    integrate_stoch

Define fout if it was omitted.
"""
function integrate_stoch(tspan::Vector{Float64}, df::Function, dg::Function, x0::Vector{Complex128},
    state::T, dstate::T, ::Void, n::Int; kwargs...) where T
    function fout(t::Float64, state::T)
        copy(state)
    end
    integrate_stoch(tspan, df, dg, x0, state, dstate, fout, n; kwargs...)
end

Base.@pure pure_inference(fout,T) = Core.Inference.return_type(fout, T)


text:
#################################
"""
    gauges = loadgauge(dirname::AbstractString; eta0::Float64=0.0, labelhead::AbstractString="Gauge ", loadeta::Bool=true, loadvel::Bool=false)

gauge*.txt reader
"""
function loadgauge(dirname::AbstractString, gaugeid::AbstractVector{Int64}=0:0;
                   eta0=0.0, labelhead::AbstractString="Gauge ", loadeta::Bool=true, loadvel::Bool=false)
    # check args
    if !isdir(dirname); error("$dirname is not found or directory"); end
    files = readdir(dirname)
    filter!(x->occursin(r"^gauge\d+\.txt$", x), files)
    if isempty(files); println("No gauge file"); return empty([], VisClaw.Gauge) end;
    nf = length(files)

    gaugeid == 0:0 && (gaugeid = 1:nf)

    nfload = length(gaugeid)

    # preallocate
    gauges = Vector{VisClaw.Gauge}(undef,nf)
    for k in gaugeid
        filename=joinpath(dirname,files[k])

        # read header
        f = open(filename,"r")
        header1 = readline(f)
        close(f)
        id = parse(Int64,header1[13:17])
        loc = [parse(Float64,header1[30:46]), parse(Float64,header1[48:64])]

        # read time-series of vars in the certain colmns
        ###dataorg = readdlm(filename, skipstart=2)
        dataorg = readdlm(filename, skipstart=3) # v5.7.0
        AMRlevel = convert.(Int64,dataorg[:,1])
        time = convert.(Float64,dataorg[:,2])
        D = convert.(Float64,dataorg[:,3])

        nt = length(time)
        if loadvel
            u = convert.(Float64,dataorg[:,4])./D
            v = convert.(Float64,dataorg[:,5])./D
        else
            u = v = empty([0.0])
        end
        if loadeta
            eta = convert.(Float64,dataorg[:,6])
            eta[D.<=1e-3] .= 0.0
            eta = eta.-eta0
        else
            eta = empty([0.0])
        end

        # label
        label = labelhead*@sprintf("%d",id)

        # instance
        gauges[k] = VisClaw.Gauge(label,id,nt,loc,AMRlevel,time,eta,u,v)
    end

    return gauges
end
#################################

loadgauge(dirname::AbstractString, gaugeid::Integer; eta0=0.0, labelhead::AbstractString="Gauge ", loadeta::Bool=true, loadvel::Bool=false) =
loadgauge(dirname, gaugeid:gaugeid; eta0=eta0, labelhead=labelhead, loadeta=loadeta, loadvel=loadvel)


text:
@cdef :PyUnicode_DecodeUTF8 PyPtr (Ptr{Cchar}, Py_ssize_t, Ptr{Cvoid})
@cdef :PyUnicode_AsUTF8String PyPtr (PyPtr,)
@cdef :PyUnicode_InternInPlace Cvoid (Ptr{PyPtr},)

const PyUnicode_Type__ref = Ref(PyNULL)
PyUnicode_Type() = pyglobal(PyUnicode_Type__ref, :PyUnicode_Type)

PyUnicode_Check(o) = Py_TypeCheckFast(o, Py_TPFLAGS_UNICODE_SUBCLASS)
PyUnicode_CheckExact(o) = Py_TypeCheckExact(o, PyUnicode_Type())

PyUnicode_From(s::Union{Vector{Cuchar},Vector{Cchar},String,SubString{String}}) =
    PyUnicode_DecodeUTF8(pointer(s), sizeof(s), C_NULL)

PyUnicode_AsString(o) = begin
    b = PyUnicode_AsUTF8String(o)
    isnull(b) && return ""
    r = PyBytes_AsString(b)
    Py_DecRef(b)
    r
end

PyUnicode_AsVector(o, ::Type{T} = UInt8) where {T} = begin
    b = PyUnicode_AsUTF8String(o)
    isnull(b) && return T[]
    r = PyBytes_AsVector(b, T)
    Py_DecRef(b)
    r
end

PyUnicode_TryConvertRule_string(o, ::Type{String}) = begin
    r = PyUnicode_AsString(o)
    isempty(r) && PyErr_IsSet() && return -1
    putresult(r)
end

PyUnicode_TryConvertRule_vector(o, ::Type{Vector{X}}) where {X} = begin
    r = PyUnicode_AsVector(o, X)
    isempty(r) && PyErr_IsSet() && return -1
    putresult(r)
end

PyUnicode_TryConvertRule_symbol(o, ::Type{Symbol}) = begin
    r = PyUnicode_AsString(o)
    isempty(r) && PyErr_IsSet() && return -1
    putresult(Symbol(r))
end

PyUnicode_TryConvertRule_char(o, ::Type{Char}) = begin
    r = PyUnicode_AsString(o)
    isempty(r) && PyErr_IsSet() && return -1
    length(r) == 1 || return 0
    putresult(first(r))
end


text:
export NpyHeader, NpyDescr, readnpy, readnpyheader, readnpydata, writenpy

using ..PyObjects, ..PyExprIO, ..PickleIO

const NPY_MAGIC = b"\x93NUMPY"

const BIGENDIAN = ENDIAN_BOM == 0x01020304

"""
    NpyDescr

The descr portion of a .npy header.
"""
struct NpyDescr
    eltype :: DataType
    bigendian :: Bool
end

"""
    NpyHeader

The header portion of a .npy file.
"""
struct NpyHeader
    raw :: String
    version_major :: Int
    version_minor :: Int
    descr :: NpyDescr
    fortran_order :: Bool
    shape :: Vector{Int}
end

"""
    readnpyheader(io::IO) :: NpyHeader

Read the .npy header.
"""
function readnpyheader(io::IO)
    # magic string
    for b in NPY_MAGIC
        if read(io, UInt8) != b
            error("Not a .npy file (incorrect magic string)")
        end
    end

    # version
    version_major = read(io, UInt8)
    version_minor = read(io, UInt8)

    if 1 ≤ version_major ≤ 3
        # header size
        if version_major < 2
            header_size = Int(ltoh(read(io, UInt16)))
        else
            header_size = Int(ltoh(read(io, UInt32)))
        end

        # read header
        raw = String(read(io, header_size))
        sizeof(raw) == header_size || error("reached end of file while reading header")
        version_major < 2 && !isascii(raw) && error("header contains non-ASCII characters")

        hdr = readpyexpr(IOBuffer(raw))::PyDict
        descr = nothing
        fortran_order = nothing
        shape = nothing
        for (k, v) in hdr.items
            k::String
            if k == "descr"
                descr = parse_descr(v::String)
            elseif k == "fortran_order"
                fortran_order = v::Bool
            elseif k == "shape"
                shape = collect(Int, (v::PyTuple).values)
            else
                error("bad header key, got $(repr(k))")
            end
        end
        descr === nothing && error("missing header key $(repr("descr"))")
        fortran_order === nothing && error("missing header key $(repr("fortran_order"))")
        shape === nothing && error("missing header key $(repr("shape"))")

        # done
        NpyHeader(strip(raw), version_major, version_minor, descr, fortran_order, shape)

    else
        error("npy version $(version_major).$(version_minor) not supported")
    end
end

function parse_descr(str)
    e = str[1]
    e in "<>|" || error("bad descr, got $(repr(str))")
    t = str[2:end]
    if t == "b1"
        T = Bool
    elseif t == "i1"
        T = Int8
    elseif t == "i2"
        T = Int16
    elseif t == "i4"
        T = Int32
    elseif t == "i8"
        T = Int64
    elseif t == "u1"
        T = UInt8
    elseif t == "u2"
        T = UInt16
    elseif t == "u4"
        T = UInt32
    elseif t == "u8"
        T = UInt64
    elseif t == "f2"
        T = Float16
    elseif t == "f4"
        T = Float32
    elseif t == "f8"
        T = Float64
    elseif t == "c4"
        T = ComplexF16
    elseif t == "c8"
        T = ComplexF32
    elseif t == "c16"
        T = ComplexF64
    elseif t == "O"
        T = Any
    else
        error("unsupported descr: $(repr(str))")
    end
    if e == '|' && T != Any && sizeof(T) > 1
        error("invalid descr, got $(repr(str))")
    end
    return NpyDescr(T, e == '>')
end

"""
    readnpydata(io::IO, h::NpyHeader; transpose=false, simplify=false)

Read the data portion of a .npy file.

If `transpose=true` then the array is transposed.

If the data is pickled, the `simplify` argument has the same meaning as for
[`PythonIO.readpkl`](@ref). Otherwise it has no effect.
"""
function readnpydata(io::IO, h::NpyHeader; transpose::Bool=false, simplify=false)
    fo = h.fortran_order
    dt = h.descr
    sz = h.shape
    if dt.eltype == Any
        x = readpkl(io, simplify=false)
        @assert fo == false
        @assert x isa PyFuncCall
        @assert x.func isa PyGlobal
        @assert x.func.mod == "numpy.core.multiarray"
        @assert x.func.attr == "_reconstruct"
        @assert x.state isa PyTuple
        @assert length(x.state.values) ≥ 5
        @assert x.state.values[2] isa PyTuple
        @assert x.state.values[2].values == sz
        @assert x.state.values[5] isa PyList
        x = x.state.values[5].values::Vector{Any}
        @assert length(x) == prod(sz)
        if simplify === nothing
            x = map(identity, x)::Vector
        else
            x = PyObjects.map_simplify(x)::Vector
            x = PyObjects.map_walk(simplify, x)::Vector
        end
        x = reshape(x, Tuple(reverse(sz)))::Array
        if !transpose
            x = _reversedims(x)::Array
        end
    else
        if transpose
            fo = !fo
            sz = reverse(sz)
        end
        x = _readarray(io, dt.eltype, Tuple(fo ? sz : reverse(sz)))::Array
        if BIGENDIAN != dt.bigendian && sizeof(st.eltype) > 1
            x = map(bswap, x)::Array
        end
        if !fo
            x = _reversedims(x)::Array
        end
    end
    return x
end

_reversedims(x::AbstractArray{T,N}) where {T,N} = N<2 ? x : permutedims(x, ntuple(i->(N+1-i), N))

function _readarray(io::IO, ::Type{T}, size::NTuple{N,Int}) where {N,T}
    @assert isbitstype(T)
    x = Array{T,N}(undef, size)
    read!(io, x)
    return x
end

"""
    readnpy(io_or_filename; transpose=false, simplify=false)

Read a .npy array from an IO stream or filename.

If `transpose=true` then the array is transposed.

If the data is pickled, the `simplify` argument has the same meaning as for
[`PythonIO.readpkl`](@ref). Otherwise it has no effect.
"""
function readnpy(io::IO; transpose::Bool=false, simplify=false)
    readnpydata(io, readnpyheader(io); transpose, simplify)
end
function readnpy(fn::AbstractString; kw...)
    open(io->readnpy(io; kw...), fn)
end

"""
    writenpy(io_or_filename, x::AbstractArray)

Write a .npy array to an IO stream or filename.
"""
function writenpy(io::IO, x::AbstractArray{T,N}) where {T,N}
    if T == Bool
        t = "b1"
    elseif T == Int8
        t = "i1"
    elseif T == Int16
        t = "i2"
    elseif T == Int32
        t = "i4"
    elseif T == Int64
        t = "i8"
    elseif T == UInt8
        t = "u1"
    elseif T == UInt16
        t = "u2"
    elseif T == UInt32
        t = "u4"
    elseif T == UInt64
        t = "u8"
    elseif T == Float16
        t = "f2"
    elseif T == Float32
        t = "f4"
    elseif T == Float64
        t = "f8"
    elseif T == ComplexF16
        t = "c4"
    elseif T == ComplexF32
        t = "c8"
    elseif T == ComplexF64
        t = "c16"
    else
        error("unsupported eltype: $T")
    end
    e = sizeof(T) == 1 ? '|' : BIGENDIAN ? '>' : '<'
    header = "{'descr': '$e$t', 'fortran_order': True, 'shape': $(Tuple{Vararg{Int}}(size(x)))}"
    if sizeof(header) + 100 > typemax(UInt16)
        version = 3
    else
        version = 1
    end
    curlen = sizeof(NPY_MAGIC) + 2 + (version==1 ? 2 : 4) + sizeof(header) + 1
    npad = mod(-curlen, 64)
    @assert 0 ≤ npad < 64
    @assert mod(curlen + npad, 64) == 0
    write(io, NPY_MAGIC)
    write(io, UInt8(version))
    write(io, UInt8(0))
    if version == 1
        write(io, UInt16(sizeof(header) + npad + 1))
    else
        write(io, UInt32(sizeof(header) + npad + 1))
    end
    write(io, header)
    write(io, ' '^npad)
    write(io, '\n')
    write(io, x)
    return
end
function writenpy(fn::AbstractString, x::AbstractArray)
    open(io->writenpy(io, x), fn, "w")
end


text:
function skip_indent(x)
    if x.typ === CSTParser.LITERAL && x.val == ""
        return true
    elseif x.typ === NEWLINE || x.typ === NOTCODE
        return true
    end
    false
end

function print_leaf(io, x, s)
    if x.typ === NOTCODE
        print_notcode(io, x, s)
    elseif x.typ === INLINECOMMENT
        print_inlinecomment(io, x, s)
    else
        write(io, x.val)
    end
end

function print_tree(io::IOBuffer, x::PTree, s::State)
    if is_leaf(x)
        print_leaf(io, x, s)
        return
    end

    ws = repeat(" ", x.indent)
    for (i, n) in enumerate(x.nodes)
        if is_leaf(n)
            print_leaf(io, n, s)
        else
            print_tree(io, n, s)
        end

        if n.typ === NEWLINE && i < length(x.nodes)
            if is_closer(x.nodes[i+1])
                write(io, repeat(" ", x.nodes[i+1].indent))
            elseif x.nodes[i+1].typ === CSTParser.Block
                write(io, repeat(" ", x.nodes[i+1].indent))
            elseif x.nodes[i+1].typ === CSTParser.Begin
                write(io, repeat(" ", x.nodes[i+1].indent))
            elseif !skip_indent(x.nodes[i+1])
                write(io, ws)
            end
        end
    end
end

@inline function print_notcode(io, x, s)
    prev_nl = true
    for l = x.startline:x.endline
        v = getline(s.doc, l)
        v == "" && continue
        if l == x.endline && v[end] == '\n'
            v = v[1:prevind(v, end)]
        end
        write(io, v)
        prev_nl = v == "\n" ? true : false
    end
end

@inline function print_inlinecomment(io, x, s)
    v = get(s.doc.comments, x.startline, "")
    isempty(v) && return
    v = getline(s.doc, x.startline)
    idx = findlast(c -> c == '#', v)
    idx === nothing && return
    idx = findlast(c -> !isspace(c), v[1:prevind(v, idx)])
    v = v[end] == '\n' ? v[nextind(v, idx):prevind(v, end)] : v[nextind(v, idx):end]
    write(io, v)
end


text:
import Multirate
import Multirate: NaiveResamplers
using Base.Test

function Base.isapprox( x1::Vector, x2::Vector )
    Nx1 = length( x1 )
    Nx2 = length( x2 )

    if Nx1 != Nx2
        @printf("x1 & x2 are different lengths vectors\n")
        return false
    end

    for i = 1:Nx1
        if !isapprox( x1[i], x2[i] )
            @printf( "Something went wrong at index %d\n", i )
            return false
        end
    end

    return true
end

Tx           = Float32
Th           = Float32
numFilters   = 32
x            = Tx[1.0:101]
resampleRate = 33/32

cutoffFreq      = 0.45
transitionWidth = 0.05
(hLen, β)       = Multirate.kaiserlength( transitionWidth, samplerate = numFilters )
hLen            = iceil(  hLen/numFilters  ) .* numFilters
h               = Multirate.firdes( hLen, cutoffFreq, DSP.kaiser, samplerate = numFilters, beta = β ) .* numFilters
h               = convert( Vector{Th}, h)

@time yNaive = NaiveResamplers.naivefilt( h, x, resampleRate, numFilters )
@time yArb   = Multirate.filt( h, x, resampleRate )

commonLen = min( length(yNaive), length(yArb) )
isapprox( yNaive[1:commonLen], yArb[1:commonLen] ) ||  display( [ [1:commonLen] yNaive[1:commonLen] yArb[1:commonLen] yNaive[1:commonLen].-yArb[1:commonLen]] )


text:
"""
	im2targets(img)

Transforms image `im` into a data representation and returns a tuple `(X,Y)` where
`X` and `Y` represent the coordinates and pixels values of the image respectively.
The image can be any of the supported `Images.jl` image representations. The rows
(e.g. variables) of `X` will correspond to the `Array` dimensions of `X`.

Read the `Images.jl` documentation for more on images.

#Examples
```
julia> A = rand(2,2)
2×2 Array{Float64,2}:
 0.825531  0.400416
 0.798304  0.187369

julia> myimg = Images.Gray.(A);

julia> X,Y = im2targets(myimg) # in X, the first line corresponds to the first dimension (e.g. y axis)
([1.0 2.0 1.0 2.0; 1.0 1.0 2.0 2.0], ColorTypes.Gray{Float64}[Gray{Float64}(0.825531), Gray{Float64}(0.798304), Gray{Float64}(0.400416), Gray{Float64}(0.187369)])
```
"""
# Function to transform from/to an image to/from a datacell where the pixel values are the targets
# and the data is represented by the pixel coordinates
function im2targets(img::I where I<:AbstractArray{T,N}) where {T,N}
	# Pre-allocate
	X = zeros(N, length(img))
	Y = Vector{T}(length(img))

	for i in eachindex(img)
		idxs = ind2sub(img, i)
		for j in eachindex(idxs)
			X[j,i] = idxs[j]
		end
		Y[i] = img[i]
	end
	return (X,Y)
end

"""
	targets2im(timg [,dim=Val{2}])

Converts a `DataCell` obtained with `im2targets` back to an image. The data variable `timg` can be either
a tuple of the form `(coordinates, pixels)` or a datacell with the coordinates as data and pixel values as
targets. `dim` is a `Val{T}` and specifies how many dimensions the output image will have. For color images,
if the channel views are used, `dim` has to be explicitly changed to `Val{3}`.

Read the `Images.jl` documentation for more on images.

#Examples
```
julia> A=rand(2,2) # grayscale image
2×2 Array{Float64,2}:
 0.516513  0.649544
 0.390564  0.455701

julia> img = Images.Gray.(A); D = im2targets(A)
([1.0 2.0 1.0 2.0; 1.0 1.0 2.0 2.0], [0.516513, 0.390564, 0.649544, 0.455701])

julia> rimg = targets2im(D, Val{2}())
2×2 Array{Float64,2}:
 0.516513  0.649544
 0.390564  0.455701

julia> D2 = datacell(D)
DataCell, 4 obs, 2 vars, 1 target(s)/obs, 4 distinct values: "0.39056363119663073"(1),"0.4557007500598864"(1),"0.5165134793129167"(1),"0.6495441098194403"(1)

julia> rimg = j4pr.targets2im(D2, Val{2}())
2×2 Array{Float64,2}:
 0.516513  0.649544
 0.390564  0.455701

julia> A=rand(3,2,2)
3×2×2 Array{Float64,3}:
[:, :, 1] =
 0.412733  0.543216
 0.603818  0.521376
 0.738603  0.295367

[:, :, 2] =
 0.701765   0.735347
 0.0223506  0.751951
 0.93287    0.121406

julia> img = j4pr.targets2im(j4pr.im2targets(img)) # FAILS
ERROR: ...

julia> img = j4pr.targets2im(j4pr.im2targets(img),Val{3}())
3×2×2 Array{Float64,3}:
[:, :, 1] =
 0.158023  0.944338
 0.52807   0.40711
 0.638118  0.853228

[:, :, 2] =
 0.0279565  0.820002
 0.896604   0.966971
 0.504744   0.0373515
```
"""
function targets2im(timg::I where I<:Tuple{AbstractArray{T,N}, AbstractVector{V}}, dim::Val{X}=Val{2}()) where {T,N,V,X}
	# Pre-allocate
	M = size(timg[1],1)
	dims = zeros(M)
	for i = 1:M
		dims[i] = maximum(timg[1][i,:]) # maximum(X[i,:])
	end
	img = Array{V,X}((Int.(dims)...))
	img[:] = timg[2][:]
	return img::Array{V,X}
end

targets2im(timg::T where T<:DataCell, dim::Val{X}=Val{2}()) where {X} = targets2im((getx!(timg), gety!(timg)),dim)

### Note: A function to transform from/to an image to/from a datacell where the pixel values are the data
# is not needed at this point. Since Images.jl uses Arrays as data container, transforming into a
# DataCell is pretty trivial. The reverse transform (e.g. visualizing some array) does and should not require
# j4pr at all.
#
# For example, to transform an N-d array into a DataCell with 1 sample:
# julia> A=rand(2,2)
# 2×2 Array{Float64,2}:
#  0.472233  0.936257
#  0.947118  0.737812
#
# julia> img = Images.Gray.(A);
#
# julia> typeof(img) # still an array
# Array{ColorTypes.Gray{Float64},2}
#
# julia> j4pr.datacell(j4pr.mat(img[:])) # convert to column matrix
# DataCell, 1 obs, 4 vars, 0 target(s)/obs
#


text:
abstract type DispersionFormula end

getindex(d::DispersionFormula, i) = getindex(d.coefs, i)

struct Sellmeier{N} <: DispersionFormula
    coefs::NTuple{N,Float64}
end

function (c::Sellmeier{N})(λ) where {N}
    rhs = c[1]
    for i = 2:2:N
        rhs += c[i]*λ^2 / (λ^2 - c[i+1]^2)
    end
    return sqrt(rhs + 1)
end

struct Sellmeier2{N} <: DispersionFormula
    coefs::NTuple{N,Float64}
end

function (c::Sellmeier2{N})(λ) where {N}
    rhs = c[1]
    for i = 2:2:N
        rhs += c[i]*λ^2 / (λ^2 - c[i+1])
    end
    return sqrt(rhs + 1)
end

struct Polynomial{N} <: DispersionFormula
    coefs::NTuple{N,Float64}
end

function (c::Polynomial{N})(λ) where {N}
    rhs = c[1]
    for i = 2:2:N
        rhs += c[i]*λ^c[i+1]
    end
    return sqrt(rhs)
end

struct RIInfo{N} <: DispersionFormula
    coefs::NTuple{N,Float64}
end

function (c::RIInfo{N})(λ) where {N}
    rhs = c[1]
    for i = 2:4:min(N, 9)
        rhs += (c[i]*λ^c[i+1]) / (λ^2 - c[i+2]^c[i+3])
    end
    for i = 10:2:N
        rhs += c[i]*λ^c[i+1]
    end
    return sqrt(rhs)
end

struct Cauchy{N} <: DispersionFormula
    coefs::NTuple{N,Float64}
end

function (c::Cauchy{N})(λ) where {N}
    rhs = c[1]
    for i = 2:2:N
        rhs += c[i]*λ^c[i+1]
    end
    return rhs
end

struct Gases{N} <: DispersionFormula
    coefs::NTuple{N,Float64}
end

function (c::Gases{N})(λ) where {N}
    rhs = c[1]
    for i = 2:2:N
        rhs += c[i] / (c[i+1] - 1/λ^2)
    end
    return rhs + 1
end

struct Herzberger{N} <: DispersionFormula
    coefs::NTuple{N,Float64}
end

function (c::Herzberger{N})(λ) where {N}
    rhs = c[1]
    rhs += c[2] / (λ^2 - 0.028)
    rhs += c[3] * (1/(λ^2 - 0.028))^2
    for i = 4:N
        pow = 2*(i - 3)
        rhs += c[i]*λ^pow
    end
    return rhs
end

struct Retro{N} <: DispersionFormula
    coefs::NTuple{N,Float64}
end

function (c::Retro{N})(λ) where {N}
    rhs = c[1] + c[2]*λ^2 / (λ^2 - c[3]) + c[4]*λ^2
    return sqrt((-2rhs - 1) / (rhs - 1))
end

struct Exotic{N} <: DispersionFormula
    coefs::NTuple{N,Float64}
end

function (c::Exotic{N})(λ) where {N}
    rhs = c[1] + c[2]/(λ^2 - c[3]) + c[4]*(λ - c[5]) / ((λ - c[5])^2 + c[6])
    return sqrt(rhs)
end

abstract type Tabulated <: DispersionFormula end

const ITP_TYPE = typeof(LinearInterpolation(sort(rand(2)), zeros(2)))

struct TabulatedNK <: Tabulated
    n::ITP_TYPE
    k::ITP_TYPE
end

function TabulatedNK(raw::Matrix{Float64})
    λ = raw[:, 1]
    n = raw[:, 2]
    k = raw[:, 3]
    TabulatedNK(LinearInterpolation(λ, n), LinearInterpolation(λ, k))
end

struct TabulatedN <: Tabulated
    n::ITP_TYPE
end

function TabulatedN(raw::Matrix{Float64})
    λ = raw[:, 1]
    n = raw[:, 2]
    TabulatedN(LinearInterpolation(λ, n))
end

struct TabulatedK <: Tabulated
    k::ITP_TYPE
end

function TabulatedK(raw::Matrix{Float64})
    λ = raw[:, 1]
    k = raw[:, 2]
    TabulatedK(LinearInterpolation(λ, k))
end


text:
# Rohan McLure, Australian National University (2019)
# Distributed implementation of the Transpose Parallel Research Kernel
using Distributed

function main()
    # Parameters
    local iterations, order
    if length(ARGS) != 2 || nprocs() == 1
        ("Usage: julia -p <# workers> transpose.jl <# iterations> <# order>")
        exit(1)
    else
        argv = map(x -> parse(Int, x), ARGS)
        iterations, order = argv
    end

    if order % nworkers() != 0
        println("Number of workers must divide order of matrix")
        exit(1)
    end

    # Rank setup
    rank = nworkers()

    # Dimension setup
    std_dim = Int(order / rank)
    RC = RemoteChannel{Channel{Matrix{Float64}}}

    # Setup channels
    in_channels = Vector{Vector{Union{RC, Nothing}}}(undef, rank)
    out_channels = Vector{Vector{Union{RC, Nothing}}}(undef, rank)
    for i in 1:rank
        in_channels[i] = Vector{Union{RC, Nothing}}(undef, rank)
        out_channels[i] = Vector{Union{RC, Nothing}}(undef, rank)
    end

    channel_lambda = () -> Channel{Matrix{Float64}}(0)

    for i in 1:rank
        for j in 1:rank
            if i == j
                in_channels[i][j] = nothing
            end
            in_channels[i][j] = RemoteChannel(channel_lambda, workers()[i])
        end
    end

    for i in 1:rank
        for j in 1:rank
            out_channels[i][j] = in_channels[j][i]
        end
    end

    dims = (order, std_dim)

    futures = Vector{Future}(undef, rank)
    @sync for k in 1 : rank
        in_c = in_channels[k]
        out_c = out_channels[k]
        futures[k] = @spawnat workers()[k] work(dims, in_c, out_c, iterations)
    end

    epsilon = 1.e-8
    results = map(fetch, futures)
    abserr = sum(map(x -> x[1], results))
    time = maximum(map(x -> x[2], results))
    nbytes = 2 * order^2 * 8

    if abserr < epsilon
        println("Solution validates")
    else
        println("error ",abserr, " exceeds threshold ",epsilon)
        println("ERROR: solution did not validate")
    end

    avgtime = time/iterations
    println("Rate (MB/s): ",1.e-6*nbytes/avgtime, " Avg time (s): ", avgtime)
end

# Local transpose operation
@inline @everywhere function do_transpose!(A, B, order)
    for jt in 1 : 32 : order
        for it in 1 : 32 : order
            for j in jt : min(jt+32-1, order)
                for i in it : min(it+32-1, order)
                    @inbounds B[i,j] += A[j,i]
                    @inbounds A[j,i] += 1.0
                end
            end
        end
    end
end

# Same function as do_transpose!, but just reassigns output
@inline @everywhere function write_block!(A, B, order)
    for jt in 1 : 32 : order
        for it in 1 : 32 : order
            for j in jt : min(jt+32-1, order)
                for i in it : min(it+32-1, order)
                    @inbounds B[i,j] = A[j,i]
                    @inbounds A[j,i] += 1.0
                end
            end
        end
    end
end

@everywhere function work(
    dims::Tuple{Int,Int}, # First component will be order
    in_channels::Vector{Union{RemoteChannel{Channel{Matrix{Float64}}}, Nothing}},
    out_channels::Vector{Union{RemoteChannel{Channel{Matrix{Float64}}}, Nothing}},
    iterations::Int
)
    # Set up local process data
    id = myid()
    n = nworkers()
    order, width_local = dims

    # Local 'column blocks' of the dist. matrix
    A = zeros(dims...)
    B = zeros(dims...)
    for i in 1:order
        for j in 1:width_local
            jl = j + (id-2)*width_local
            A[i,j] = order * (jl - 1) + (i - 1)
        end
    end

    # Views into the diagonal block
    local_frame_A = view(A, (id-2)*width_local+1:(id-1)*width_local, :)
    local_frame_B = view(B, (id-2)*width_local+1:(id-1)*width_local, :)

    # Arrays used for communication blocks
    local block_in
    block_out = Matrix{Float64}(undef, width_local, width_local)

    local t0, t1
    for k in 0 : iterations
        if k == 1
            t0 = time_ns()
        end

        do_transpose!(local_frame_A, local_frame_B, width_local)

        for phase in 1 : n-1
            send_idx = (id-2 + phase) % n + 1
            recv_idx = (id-2 + n - phase) % n + 1

            inbox = in_channels[recv_idx]
            outbox = out_channels[send_idx]

            frame = view(A, (send_idx-1)*width_local+1:send_idx*width_local, :)
            write_block!(frame, block_out, width_local)

            @sync begin
                @async begin
                    frame = view(B, (recv_idx-1)*width_local+1:recv_idx*width_local, :)
                    block_in = take!(inbox)
                    frame .+= block_in
                end
                @async begin
                    put!(outbox, block_out)
                end
            end
        end
    end
    t1 = time_ns()

    # Validate the local solution - sending the error back to master
    addit = (0.5*iterations) * (iterations+1)
    abserr = 0.0
    for i in 1:order
        for j in 1:width_local
            jl = j + (id-2)*width_local
            temp = (order * (i-1) + (jl-1)) * (iterations+1)
            @inbounds abserr = abserr + abs(B[i,j] - (temp+addit))
        end
    end
    return (abserr, (t1-t0)*1.e-9)
end

main()


text:
module EmpiricalModeDecomposition

using Interpolations
using Dierckx
using IterTools
using Random
using Statistics
using Base.Iterators

import Base: iterate, IteratorSize, eltype

export EMDIterable
export sift, sift!
export emd, eemd, ceemd, maketestdata

include("testdata.jl")
include("utils.jl")
include("sift.jl")

"""
Iterator for the Empirical Mode Decomposition.
The time series values are an AbstractVector of type `T`,
and the time positions are an AbstractVector of type `U`.
"""
struct EMDIterable{T<:AbstractVector, U<:AbstractVector}
    yvec::T
    xvec::U
end

eltype(::Type{EMDIterable{T,U}}) where {T,U} = T

function iterate(iter::EMDIterable, imf_prev=(iter.yvec,false))
    if imf_prev[2]
        return nothing
    elseif sum(abs, imf_prev[1]) > 0.0 && !ismonotonic(imf_prev[1])
        imf = sift(imf_prev[1], iter.xvec; reltol=0.1)
        return imf, (imf_prev[1]-imf, false)
    else
        return imf_prev[1], (imf_prev[1], true)
    end
end

Base.IteratorSize(::Type{EMDIterable{U,V}}) where {U,V} = Base.SizeUnknown()

"""
    emd(measurements, xvec, numimfs=6)

Return the Intrinsic Mode Functions (IMF) and the residual of the Empirical Mode
Decomposition (EMD) of the `measurements` given on time steps given in `xvec`.
"""
function emd(measurements, xvec, num_imfs=6)
    collect(take(EMDIterable(measurements, xvec), num_imfs))
end

"""
    eemd(measurements, xvec, numtrails=100)

Return the Intrinsic Mode Functions and the residual of the ensemble Empirical Mode
Decomposition (EMD) of the `measurements` given on time steps `xvec`.
"""
function eemd(measurements, xvec, numtrails=100, num_imfs=6)
    random = randn(length(measurements))
    imfs_mean = emd(measurements .+ random, xvec, num_imfs)
    num_imfs = length(imfs_mean)
    for i in 1:numtrails
        randn!(random)
        imfs = emd(measurements .+ random, xvec, num_imfs)
        if length(imfs) < length(imfs_mean)
            imfs_mean[length(imfs)] = sum(imfs_mean[length(imfs):end])
            imfs_mean = imfs_mean[1:length(imfs)]
        elseif length(imfs) > length(imfs_mean)
            imfs[length(imfs_mean)] = sum(imfs[length(imfs_mean):end])
            imfs = imfs[1:length(imfs_mean)]
        end
        imfs_mean .+= imfs
    end

    imfs_mean ./= numtrails
end

"""
Iterator for the Complete Empirical Mode Decomposition.
The time series values are an AbstractVector of type T,
and the time positions are an AbstractVector of type U.
"""
struct CEEMDIterable{U<:AbstractVector,V<:AbstractVector,T<:AbstractVector}
    "Vector with the values of the time series"
    yvec::U
    "Vector with the positions of the time series"
    xvec::V
    "Ensemble of noise which is a Vector of Vectors of the same size as yvec."
    noise_ens::T
end

"""
Intermediate results of the CEEMD Iteration
"""
struct CEEMDState
    yvec
    iter_ens
    imf_state_ens
    finished::Bool
end
Base.IteratorSize(::Type{CEEMDIterable{U,V,T}}) where {U,V,T} = Base.SizeUnknown()

function Base.iterate(iter::CEEMDIterable)
    imf0 = mean([sift(iter.yvec+noise, iter.xvec; reltol=0.1) for noise in iter.noise_ens])
    iter_ens = EMDIterable.(iter.noise_ens, [iter.xvec])
    state_ens = Tuple[]
    imf_state_ens = iterate.(iter_ens)
    imf0, CEEMDState(iter.yvec-imf0, iter_ens, imf_state_ens, false)
end

function Base.iterate(iter::CEEMDIterable, state::CEEMDState)

    vstop = var(iter.yvec)*1e-10

    if state.finished

        return nothing

    elseif sum(abs,state.yvec)>vstop && !ismonotonic(state.yvec)

        # TODO: Format for clarity?
        imf = vec(median(hcat([sift(state.yvec+noise[1], iter.xvec; reltol=0.1) for noise in state.imf_state_ens]...),dims = 2))

        for iens in 1:length(state.iter_ens)
            r = iterate(state.iter_ens[iens], state.imf_state_ens[iens][2])
            if isnothing(r)
                fill!(state.imf_state_ens[iens][1], 0)
            else
                imf_noise, ensstate = r
                state.imf_state_ens[iens] = imf_noise, ensstate
            end
        end

        newstate = CEEMDState(state.yvec-imf, state.iter_ens, state.imf_state_ens, false)
        return imf,newstate
    else
        @show state.yvec
        return state.yvec,CEEMDState(state.yvec, state.iter_ens, state.imf_state_ens, true)
    end
end

"""
    ceemd(measurements, xvec; num_imfs=6)

Compute the Complete Empirical Mode Decomposition (CEMD) of the time series with values
`measurements` and time steps `xvec`.
`num_imfs` is the number of Intrinsic Mode Functions.
Returns a list of num_imfs + 1 vectors of the same size as measurements.
"""
function ceemd(measurements, xvec; num_imfs=6, numtrails=100, β=0.04,
    noise_ens = [β*std(measurements) .* randn(length(xvec)) for i in 1:numtrails])
    imfs = collect(take(CEEMDIterable(measurements, xvec, noise_ens), num_imfs))
    @show size.(imfs)
    residual = measurements - sum(imfs)
    push!(imfs, residual)
    return imfs
end

function tone_masking(ys, xs, tone)
    ys_plus = ys .+ tone
    ys_minus = ys .- tone
    phi_plus = sift(ys_plus)
    phi_minus = sift(ys_minus)
    @. 0.5 * (phi_plus + phi_minus)
end

function iaestimation(imf, xs)
    r = abs.(imf)
    maxs = Int[]
    mins = Int[]
    localmaxmin!(imf, maxs, mins)
    interpolate(xs[maxs], imf[maxs], xs, DierckXInterp())
end

function iterAMremoval(imf, xvec)
    g = imf
    n=1
    b = zero(imf)
    while any(b .!=1) && n<=3
        b = iaestimation(imf, xvec)
        g = g ./ b
        n+=1
    end
    g
end

#
# function ceemd(measurements, xvec, num_imfs=6, numtrails=100, β=0.02, noise_ens = [β .* randn(length(xvec)) for i in 1:numtrails])
#     imfs = typeof(measurements)[]
#     #noise_imfs = [EMD(noise) for noise in noise_ens]
#     ycur = sift(measurements)
#     for i in 1:numtrails
#         ycur += sift(measurements .+ noise_ens[i],xvec)
#     end
#     ycur ./= numtrails
#     k=0
#     while length(imfs)<num_imfs && sum(abs, ycur) > 0.0 && !ismonotonic(ycur)
#         println(sum(ycur))
#         global k+=1
# #        plot(ycur)
#         y = sift(ycur, )
#         y = @parallel (+) for i in 1:numtrails
#             summand = ycur + 0.02*EmpiricalModeDecomposition.emd(randn(length(ycur)), xvec, k)[end]
#             EmpiricalModeDecomposition.sift(summand, xvec)
#         end
#         y ./=numtrails
#         println(maximum(y))
#         ycur = ycur-y
#         push!(imfs, ycur)
#     end
#     imfs
# end

end # module


text:
module BoundingBoxes

export @boundingbox

macro boundingbox(ex)
    bound_name = ex.args[1]
    axes = ex.args[2:end]
    # construct Bounding Box type
    types = [Expr(:(::),Symbol(axe*"_max"),:T) for axe in axes]
    append!(types,[Expr(:(::),Symbol(axe*"_min"),:T) for axe in axes])
    curly = Expr(:curly, bound_name, Expr(:(<:), :T, :Number))
    bound = Expr(:type, true, curly, Expr(:block, types...))

    # construct equality
    eqs = [Expr(:call, :(==),
                 dots("a", "_max", axe),
                 dots("b", "_max", axe)) for axe in axes]
    append!(eqs, [Expr(:call, :(==),
                 dots("a", "_min", axe),
                 dots("b", "_min", axe)) for axe in axes])
    eq_return = Expr(:return, andlist(eqs))
    eq_block = Expr(:block, eq_return)
    eq_args = :(==(a::$(bound_name), b::$(bound_name)))
    equality = Expr(:function, eq_args, eq_block)

    # construct empty bounds
    empty_max = [Expr(:call, :typemin, :T) for i = 1:length(axes)]
    empty_min = [Expr(:call, :typemax, :T) for i = 1:length(axes)]
    empty_bounds = :($(bound_name)(T) = $(bound_name){T}($(empty_max...), $(empty_min...)))

    # promotion handling function
    promote_bounds = :($(bound_name)(x...) = $(bound_name)(promote(x...)...))

    # construct bounds updating
    update_args = :(update!(a::$(bound_name), b::AbstractArray))
    update_body = [Expr(:(=), dots("a", "_max", axes[i]),
                        Expr(:call,
                            :(max), dots("a", "_max", axes[i]),
                            Expr(:ref, :b, i))
                        ) for i = 1:length(axes)]
    append!(update_body, [Expr(:(=), dots("a", "_min", axes[i]),
                            Expr(:call,
                                    :(min), dots("a", "_min", axes[i]),
                                    Expr(:ref, :b, i))
                            ) for i = 1:length(axes)])
    update_block = Expr(:block, update_body..., :(return nothing))
    update_bounds = Expr(:function, update_args, update_block)

    # construct is (a in b)
    in_eqs = [Expr(:call, :(>),
                 dots("b", "_max", axe),
                 dots("a", "_max", axe)) for axe in axes]
    append!(in_eqs, [Expr(:call, :(<),
                 dots("b", "_min", axe),
                 dots("a", "_min", axe)) for axe in axes])
    in_return = Expr(:return, andlist(in_eqs))
    in_block = Expr(:block, in_return)
    in_args = :(Base.in(a::$(bound_name), b::$(bound_name)))
    in_bounds = Expr(:function, in_args, in_block)

    quote
        $(esc(bound)) # create bounding box type
        $(esc(equality)) # create equality method (==)
        $(esc(promote_bounds)) # create arg promoting method
        $(esc(update_bounds)) # create method for updating against AbstractArray
        $(esc(in_bounds)) # create in method
        $(esc(empty_bounds)) # create empty bounds (Bounds(T))
    end
end

# member access
function dots(parent, elt, axe)
    return Expr(:(.), Symbol(parent), QuoteNode(Symbol(axe*elt)))
end

# create a list of chained "&&" Exprs
function andlist(exprs)
    last_expr = exprs[1]
    expr = Expr(:block)
    for i = 2:length(exprs)
        expr = Expr(:(&&), exprs[i], last_expr)
        last_expr = expr
    end
    return expr
end

end # module


text:
"""
    Time(duration, clock)

An instant in `Time` is represented as the `duration` since the "epoch" of the
`clock`, where the epoch is the instant the `clock` reads zero.

Times can be subtracted to get a relative duration, if a relationship between
their clocks is known.
"""
struct Time{C <: Clock, D <: Duration}
    duration::D
    clock::C # This might be singleton, or it might have complex data
end

function Base.show(io::IO, t::Time)
    print(io, t.duration)
    print(io, " (since ")
    print(io, t.clock)
    print(io, ")")
end

# Mathematical operations
@inline Base.:+(d::Duration, t::Time) = Time(d + t.duration, t.clock)
@inline Base.:+(t::Time, d::Duration) = Time(t.duration + d, t.clock)
@inline Base.:-(t::Time, d::Duration) = Time(t.duration - d, t.clock)

@inline Base.:-(t1::Time{C1}, t2::Time{C2}) where {C1<:TAIClock, C2<:TAIClock} = t1.duration - t2.duration + (t1.clock - t2.clock)

"""
    Time(time, clock)

Convert the input `time` to use a new `clock`, if a relationship is known.
"""
function Time(time::Time{C1}, clock::C2) where {C1 <: TAIClock, C2 <: TAIClock}
    return Time(time.duration + (time.clock - clock), clock)
end


text:
"""
    FDEProblem

General type for all kinds of problems in FractionalDiffEq.jl.
"""
abstract type FDEProblem end

"""
    MultiTermsFODEProblem(parameters, orders, rightfun, u0, T)

    MultiTermsFODEProblem(parameters, orders, rightfun, rparameters, rorders, u0, T)

Define a multi-terms fractional ordinary differential equation.
"""
struct MultiTermsFODEProblem <: FDEProblem
    parameters::AbstractArray
    orders::AbstractArray
    rightfun::Union{Function, Number}
    rparameters::Union{AbstractArray, Nothing}
    rorders::Union{AbstractArray, Nothing}
    u0::AbstractArray
    tspan::Union{Tuple{Number, Number}, Number}
end

#=MultiTermsFODEProblem constructor, dispatch for closedform problem=#
MultiTermsFODEProblem(parameters, orders, rightfun, u0, T) = MultiTermsFODEProblem(parameters, orders, rightfun, nothing, nothing, u0, T)
#MultiTermsFODEProblem(parameters, orders, rightfun, u0, t0, T) = MultiTermsFODEProblem(parameters, orders, rightfun, nothing, nothing, u0, T)

"""

    SingleTermFODEProblem(f, α, u0, tspan)

Define a single term fractional ordinary differential equation, there are only one fractional differential operator in this problem.
"""
struct SingleTermFODEProblem <: FDEProblem
    f::Function
    α::Float64
    u0
    tspan::Union{Tuple, Number}
end

"""
Base type for all of the FractionalDiffEq algorithms
"""
abstract type FractionalDiffEqAlgorithm end

"""
    FPDEProblem(α, β, T, M, N)

Define fractional order partial differential equations problem
"""
struct FPDEProblem <: FDEProblem
    α
    β
    T
    M
    N
end

"""
    FDDEProblem(f, ϕ, α, τ, tspan)

Construct a fractional delayed differential equation problem.
"""
struct FDDEProblem <: FDEProblem
    f::Function
    ϕ::Union{Number, Function}
    α::Union{Number, Function}
    τ::Union{Number, AbstractArray}
    tspan::Union{Number, Tuple}
end

#=FDDEProblem constructor=#
FDDEProblem(f, ϕ, α, τ, T) = FDDEProblem(f, ϕ, α, τ, T, nothing)

"""
    FDDESystem(f, ϕ, α, τ, T)

Construct system of fractional delay differential equations problem.
"""
struct FDDESystem <: FDEProblem
    f::Function
    ϕ::AbstractArray
    α
    τ::Number
    T
end

"""
    FDDEMatrixProblem(α, τ, A, B, f, x0, tspan)

Construct a fractional matrix differential equation with delay. The general type is:

```math
D_{t_0}^\\alpha\\textbf{x}(t)=\\textbf{A}(t)\\textbf{x}(t)+\\textbf{B}(t)\\textbf{x}(t-\\tau)+\\textbf{f}(t)
```
"""
struct FDDEMatrixProblem <: FDEProblem
    α::Float64
    τ::Float64
    A::Matrix
    B::Matrix
    f::Vector
    x0::Function
    tspan::Tuple{Float64, Float64}
end

"""
    FODESystem(f, α, u0, tspan)

Define system of fractional differential equations.
"""
struct FODESystem <: FDEProblem
    f::Function
    α::AbstractArray
    u0::AbstractArray
    tspan::Union{Tuple, Number}
end

FODESystem(f, α, u0, T) = FODESystem(f, α, u0, 0, T)# Start from t=0

"""
    DODEProblem(parameters, orders, interval, tspan, rightfun)

Define distributed order differential equation problem.
"""
struct DODEProblem <: FDEProblem
    parameters::AbstractArray
    orders::AbstractArray
    interval::Tuple
    rightfun::Function
    u0
    tspan
end

"""
    FractionalDifferenceProblem(f, α, x0)

Define fractional difference equation problems.
"""
struct FractionalDifferenceProblem <: FDEProblem
    fun::Function
    α
    u0
end

"""
    FractionalDifferenceSystem(f, α, u0)

Define Fractional Difference System with the general constructure:

```math
{_{a-m}^G\\nabla_k^\\alpha x(k+1)}=f(x(k))\\
```

With given initial condition ``x(i)``.
"""
struct FractionalDifferenceSystem <: FDEProblem
    fun::Function
    α
    u0
end

"""
    FIEProblem(parameters, orders, rightfun, tspan)

Define fractional integral equation problems.
"""
struct FIEProblem <: FDEProblem
    parameters::AbstractArray
    orders::AbstractArray
    rightfun::Union{Function, Number}
    tspan
end

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

abstract type AbstractFDESolution end

struct FODESolution <: AbstractFDESolution
    t::AbstractArray
    u::AbstractArray
end

struct FDifferenceSolution <: AbstractFDESolution
    t::AbstractArray
    u::AbstractArray
end

struct FIESolution <: AbstractFDESolution
    t::AbstractArray
    u::AbstractArray
end

struct DODESolution <: AbstractFDESolution
    t::AbstractArray
    u::AbstractArray
end

"""
# Usage

    solve(prob::SingleTermFODEProblem, h, PECE())

Predict-Evaluate-Correct-Evaluate algorithm.

For more details, please refer to [Predictor-Corrector algorithms](https://en.wikipedia.org/wiki/Predictor%E2%80%93corrector_method)

This PECE algorithm is taken from Diethelm's paper.

### References

```tex
@article{
title={A predictor-corrector approach for the numerical solution of fractional differential equations},
author={Diethelm, Kai and Ford, Neville J. and Freed, Alan D.}
doi={https://doi.org/10.1023/A:1016592219341}
}
```
"""
struct PECE <: FractionalDiffEqAlgorithm end
#TODO: Use Richardson extrapolation to refine the PECE algorithms
#TODO: Rename as ABM

"""
    solve(problem::SingleTermFODEProblem, h, PECE())

Generalproblem solving API for solving FDE problems.
"""
function solve(FODE::SingleTermFODEProblem, h::Float64, ::PECE)
    @unpack f, α, u0, tspan = FODE
    t0 = tspan[1]; T = tspan[2]
    N::Int64 = round(Int, (T-t0)/h)
    y = zeros(N+1) # Initialize the solution array

    @fastmath @inbounds @simd for n ∈ 0:N
        # Handling initial value
        temp=0
        for k=0:ceil(Int, α)-1
            temp += u0[k+1]*(t0+(n+1)*h)^k/factorial(k)
        end
        y[n+1] = temp + h^α/gamma(α+2)*(f(t0+(n+1)*h, predictor(f, y, α, n, h, u0, t0)) + right(f, y, α, n, h))
    end

    tspan = collect(t0:h:T)
    return FODESolution(tspan, y)
end

function right(f, y, α, n, h::Float64)
    temp = zero(Float64)

    @fastmath @inbounds @simd for j = 0:n
        temp += A(j, n, α)*f(j*h, y[j+1])
    end

    return temp
end

function predictor(f, y, α::Float64, n::Integer, h::Float64, u0, t0)
    predict = 0
    leftsum = 0

    l = ceil(Int, α)

    # Handling initial value
    for k=0:l-1
        leftsum += u0[k+1]*(t0+(n+1)*h)^k/factorial(k)
    end

    @turbo for j ∈ 0:n
        predict += B(j, n, α)*f(j*h, y[j+1])
    end

    return leftsum + h^α/α*predict
end

function A(j, n, α)
    if j == 0
        return n^(α+1) - (n-α)*(n+1)^α
    elseif 1 ≤ j ≤ n
        return (n-j+2)^(α+1) + (n-j)^(α+1) - 2*(n-j+1)^(α+1)
    elseif j == n+1
        return 1
    end
end

B(j, n, α) = ((n + 1 - j)^α - (n - j)^α) # Moved the h^α/α to the end of predictor: return leftsum + h^α/α*predict


text:
using Thunks
using Documenter

DocMeta.setdocmeta!(Thunks, :DocTestSetup, :(using Thunks); recursive=true)

makedocs(;
    modules=[Thunks],
    authors="Tyler Benster <thunks.jl@tylerbenster.com> and contributors",
    repo="https://github.com/tbenst/Thunks.jl/blob/{commit}{path}#{line}",
    sitename="Thunks.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://tbenst.github.io/Thunks.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/tbenst/Thunks.jl",
)


text:
module SortMyLife

export sortmylife

using Cassette
Cassette.@context SortCtx

struct Task
  x::String
end
Base.show(io::IO, t::Task) = print(io, t.x)

function Base.isless(x::Task, y::Task)
  println("Which should you do first?")
  println("a) $(x.x)")
  println("b) $(y.x)")
  @label readinput
  line = readline(stdin)
  if line == "a"
    false
  elseif line == "b"
    true
  else
    println("Did not receive a or b, try again")
    @goto readinput
  end
end

function Cassette.overdub(ctx::SortCtx, ::typeof(Base.isless), x, y)
  ctx.metadata[] += 1
  x < y
end

function sortmylife(tasks::Vector{Task})
  count1 = Ref(0)
  ctx = SortCtx(metadata = count1)
  res = Cassette.overdub(ctx, sort, unique(tasks))
  println("I sorted your life out and it only took $(ctx.metadata.x) comparisons!")
  res
end

"""
papertasks = ["Add section on control flow",
	        		"Add inverse graphics plots",
              "Add code examples",
              "Discuss scaling in section 3",
              "Relook at glucose example",
              "add Kernel comparison plot",
              "Add approximate bayesian computation benchmarks",
              "Add theorem about composing primitives",
              "Revise Draft",
              "Add a third example.. Program learning, Physics?"]

Usage: sortmylife(papertasks)
"""
sortmylife(tasks::AbstractVector) = sortmylife(Task.(tasks))

"""
tasks = \"\"\"
. Send Mum a birthday present
. Email people from hacker news
. Put tomorrow variational talk into calednar
. Setup Papaxray
. Pay Ilans thing
. Pay Taras thing
. Email people from hacker news
. Email back people who emailed me from hacker news
. Investigate Incorporate Basis.
. Think about what kind of structure I want Basis to have
. Plan San Francisco trip
. Try to formalize causal entailment / causal sufficiency
. Add HP definition of explanation to overleaf paper
. Read situation calculus paper
. Make a larger plan for submission
. Try to develop a prototype
. Fix one Omega issue
. Explore Omega redesign.
. Read GAIKAS spectactular Empire
\"\"\"

sortmylife(tasks)
"""
sortmylife(::String; delim = "\n") = strip.(split(tasks, delim))

function readtasks()
  tasks = Task[]
  while true
    println("Please add a task to add or press enter to be done.")
    line = readline(stdin)
    if line == ""
      break
    else
      println("Adding task: $line")
      push!(tasks, Task(line))
    end
  end
  tasks
end

sortmylife() = sortmylife(readtasks())

end


text:
using DiffEqFlux, GalacticOptim, OrdinaryDiffEq, RecursiveArrayTools

u0 = Float32[0.; 2.]
du0 = Float32[0.; 0.]
tspan = (0.0f0, 1.0f0)
t = range(tspan[1], tspan[2], length=20)

model = FastChain(FastDense(2, 50, tanh), FastDense(50, 2))
p = initial_params(model)
ff(du,u,p,t) = model(u,p)
prob = SecondOrderODEProblem{false}(ff, du0, u0, tspan, p)

function predict(p)
    Array(solve(prob, Tsit5(), p=p, saveat=t, sensealg = InterpolatingAdjoint(autojacvec=ZygoteVJP())))
end

correct_pos = Float32.(transpose(hcat(collect(0:0.05:1)[2:end], collect(2:-0.05:1)[2:end])))

function loss_n_ode(p)
    pred = predict(p)
    sum(abs2, correct_pos .- pred[1:2, :]), pred
end

data = Iterators.repeated((), 1000)
opt = ADAM(0.01)

l1 = loss_n_ode(p)

cb = function (p,l,pred)
    @show l
    l < 0.01 && Flux.stop()
end

optfunc = GalacticOptim.OptimizationFunction((x, p) -> loss_n_ode(x), GalacticOptim.AutoZygote())
optprob = GalacticOptim.OptimizationProblem(optfunc, p)
res = GalacticOptim.solve(optprob, opt, cb=cb, maxiters = 100)
l2 = loss_n_ode(res.minimizer)
@test l2 < l1

function predict(p)
    Array(solve(prob, Tsit5(), p=p, saveat=t, sensealg = QuadratureAdjoint(autojacvec=ZygoteVJP())))
end

correct_pos = Float32.(transpose(hcat(collect(0:0.05:1)[2:end], collect(2:-0.05:1)[2:end])))

function loss_n_ode(p)
    pred = predict(p)
    sum(abs2, correct_pos .- pred[1:2, :]), pred
end

data = Iterators.repeated((), 1000)
opt = ADAM(0.01)

loss_n_ode(p)

cb = function (p,l,pred)
    @show l
    l < 0.01 && Flux.stop()
end
optfunc = GalacticOptim.OptimizationFunction((x, p) -> loss_n_ode(x), GalacticOptim.AutoZygote())
optprob = GalacticOptim.OptimizationProblem(optfunc, p)
res = GalacticOptim.solve(optprob, opt, cb=cb, maxiters = 100)
l2 = loss_n_ode(res.minimizer)
@test l2 < l1

function predict(p)
    Array(solve(prob, Tsit5(), p=p, saveat=t, sensealg = BacksolveAdjoint(autojacvec=ZygoteVJP())))
end

correct_pos = Float32.(transpose(hcat(collect(0:0.05:1)[2:end], collect(2:-0.05:1)[2:end])))

function loss_n_ode(p)
    pred = predict(p)
    sum(abs2, correct_pos .- pred[1:2, :]), pred
end

data = Iterators.repeated((), 1000)
opt = ADAM(0.01)

loss_n_ode(p)

cb = function (p,l,pred)
    @show l
    l < 0.01 && Flux.stop()
end

optfunc = GalacticOptim.OptimizationFunction((x, p) -> loss_n_ode(x), GalacticOptim.AutoZygote())
optprob = GalacticOptim.OptimizationProblem(optfunc, p)
res = GalacticOptim.solve(optprob, opt, cb=cb, maxiters = 100)
l2 = loss_n_ode(res.minimizer)
@test l2 < l1


text:
# using strings in Julia
# Julia command line

julia> a = "praveen"

julia> a[1]
'p': ASCII/Unicode U+0070 (category Ll: Letter, lowercase)

# index starts at 1!

# for loop in Julia
julia> for i in a
           println(i)
       end

p
r
a
v
e
e
n

# no indentations unlike Python ;)

# collect(<string>)
julia> collect(a)
7-element Array{Char,1}:
 'p'
 'r'
 'a'
 'v'
 'e'
 'e'
 'n'

# isvalid.()
julia> isvalid.(collect(a))
7-element BitArray{1}:
 true
 true
 true
 true
 true
 true
 true


text:
# Equation abstraction and timesteppers for OceanTurb.jl

import LinearAlgebra: Tridiagonal
import Base: @propagate_inbounds

"""
    Timestepper(stepper, args...)

Generalized timestepper constructor.
"""
function Timestepper(stepper, args...)
    fullsteppername = Symbol(stepper, :Timestepper)
    return eval(Expr(:call, fullsteppername, args...))
end

function iterate!(model, Δt, Nt)
    for step = 1:Nt
        iterate!(model, Δt)
    end
    return nothing
end

"""
    iterate!(model; Δt, Nt)

Step `model` forward in time for `Nt` steps with step size Δt.
"""
iterate!(model; Δt, Nt) = iterate!(model, Δt, Nt)

"Update the clock after one iteration."
function tick!(clock, Δt)
    clock.time += Δt
    clock.iter += 1
    return nothing
end

#
# Explicit and implicit kernels
#

"Evaluate the right-hand-side of ∂ϕ∂t for an explicit time-stepping method."
function calc_explicit_rhs!(rhs, eqn, solution, m)

    # Function barrier for better performance and forced specialization
    # on Kϕ and Rϕ
    function kernel!(rhs, ϕ, M::MF, L::LF, K::KF, R::RF, m) where {MF, LF, KF, RF}
        for i in eachindex(ϕ)
            @inbounds rhs[i] = (-L(m, i) * ϕ[i] - ∂zM(M(m, i+1), M(m, i), ϕ, i)
                                    + ∇K∇ϕ(K(m, i+1), K(m, i), ϕ, i) + R(m, i))
        end
    end

    ntuple(Val(length(solution))) do j
        Base.@_inline_meta
           ϕ = solution[j]
        rhsϕ = rhs[j]
          Mϕ = eqn.M[j]
          Lϕ = eqn.L[j]
          Kϕ = eqn.K[j]
          Rϕ = eqn.R[j]

        kernel!(rhsϕ, ϕ, Mϕ, Lϕ, Kϕ, Rϕ, m)
    end
    return nothing
end

"Evaluate the right-hand-side of ∂ϕ∂t for an implicit time-stepping method."
function calc_implicit_rhs!(rhs, eqn, solution, m)

    N = m.grid.N

    ntuple(Val(length(solution))) do j
        Base.@_inline_meta
           ϕ = solution[j]
        rhsϕ = rhs[j]
          Mϕ = eqn.M[j]
          Kϕ = eqn.K[j]
          Rϕ = eqn.R[j]

        for i in eachindex(rhsϕ)
            @inbounds rhsϕ[i] = Rϕ(m, i)
        end

        @inbounds rhsϕ[N] = (-∂zM(Mϕ(m, N+1), Mϕ(m, N), ϕ, N)
                                + ∇K∇ϕ(Kϕ(m, N+1), -zero(eltype(ϕ)), ϕ, N) + Rϕ(m, N))

        # Advective divergence included in lhs
        @inbounds rhsϕ[1] = ∇K∇ϕ(-zero(eltype(ϕ)), Kϕ(m, 1), ϕ, 1) + Rϕ(m, 1)

    end
    return nothing
end

@inline K_op(m, K, iᶠ, iᶜ) = K(m, iᶠ) / Δc(m.grid, iᶠ) / Δf(m.grid, iᶜ)
@inline M_op(m, M, iᴹ, iᶜ) = M(m, iᴹ) / Δc(m.grid, iᶜ)

"Build backward Euler operator for diffusive problems."
function calc_diffusive_lhs!(lhs, eqn, solution, Δt::T, m) where T
    ntuple(Val(length(solution))) do j
        Base.@_inline_meta
         ϕ = solution[j]
        Mϕ = eqn.M[j]
        Lϕ = eqn.L[j]
        Kϕ = eqn.K[j]
         L = lhs[j]

        for i in interiorindices(ϕ)
            @inbounds begin
                L.du[i]   = Δt * (M_op(m, Mϕ, i+1, i+1) - K_op(m, Kϕ, i+1, i))
                L.d[i]    = one(T) + Δt * (Lϕ(m, i) + K_op(m, Kϕ, i+1, i) + K_op(m, Kϕ, i, i) - M_op(m, Mϕ, i, i+1))
                L.dl[i-1] = -Δt * K_op(m, Kϕ, i, i)
            end
        end

        # Bottom row
        @inbounds L.du[1] = Δt * (M_op(m, Mϕ, 2, 2) - K_op(m, Kϕ, 2, 1))
        @inbounds L.d[1] = one(T) + Δt*(Lϕ(m, 1) + K_op(m, Kϕ, 2, 1) - M_op(m, Mϕ, 1, 2))

        # Top row
        N = length(ϕ)
        @inbounds L.dl[end] = -Δt * K_op(m, Kϕ, N, N)
        @inbounds L.d[end]  = one(T) + Δt * (Lϕ(m, N) + K_op(m, Kϕ, N, N) - M_op(m, Mϕ, N, N+1))
    end

    return nothing
end

#
# Timestepping methods
#

function update!(bcs, eqn, solution, m)
    N = m.grid.N
    eqn.update!(m)
    ntuple(Val(length(solution))) do j
        Base.@_inline_meta
        ϕ = solution[j]
        Kϕ = eqn.K[j]
        bcsϕ = bcs[j]
        fill_bottom_ghost_cell!(bcsϕ.bottom, ϕ, Kϕ(m, 1), m)
        fill_top_ghost_cell!(bcsϕ.top, ϕ, Kϕ(m, N+1), m, N)
    end
    return nothing
end

#
# ForwardEuler timestepper
#

"""
    ForwardEulerTimestepper(eqn, solution)

Construct a `FowardEulerTimestepper`.
"""
struct ForwardEulerTimestepper{E, R} <: Timestepper
    eqn :: E
    rhs :: R
    function ForwardEulerTimestepper(eqn, solution, args...)
        rhs = deepcopy(solution)
        for fld in rhs
            set!(fld, 0)
        end
        new{typeof(eqn), typeof(rhs)}(eqn, rhs)
    end
end

"Take one forward Euler step."
function forward_euler_step!(rhs, solution, Δt)
    for j in 1:length(solution)
        ϕ = solution[j]
        rhsϕ = rhs[j]

        for i in eachindex(ϕ)
            @inbounds ϕ[i] += Δt * rhsϕ[i]
        end
    end
    return nothing
end

# Forward Euler timestepping
function iterate!(m::AbstractModel{TS}, Δt) where TS <: ForwardEulerTimestepper
    update!(m.bcs, m.timestepper.eqn, m.solution, m)
    calc_explicit_rhs!(m.timestepper.rhs, m.timestepper.eqn, m.solution, m)
    forward_euler_step!(m.timestepper.rhs, m.solution, Δt)
    tick!(m.clock, Δt)
    return nothing
end

#
# BackwardEuler timestepper
#

struct BackwardEulerTimestepper{E, R, L} <: Timestepper
    eqn :: E
    rhs :: R
    lhs :: L
    function BackwardEulerTimestepper(eqn, solution, lhs)
        rhs = deepcopy(solution)
        for fld in rhs
            set!(fld, 0)
        end
        new{typeof(eqn), typeof(rhs), typeof(lhs)}(eqn, rhs, lhs)
    end
end

function Tridiagonal(fld::AbstractField)
    T = eltype(fld)
    A = arraytype(fld)
    N = length(fld)
    return Tridiagonal{T, A}(zeros(N-1), zeros(N), zeros(N-1))
end

"""
    build_lhs(solution)

Build Tridiagonal matrices corresponding to the LHS
for each of `solution`'s equations.
"""
function build_lhs(solution)
    return ntuple(Val(length(solution))) do i
        Tridiagonal(solution[i])
    end
end

"Update solution by inverting Tridiagonal lhs matrix."
function backward_euler_step!(rhs, lhs, solution, Δt)
    ntuple(Val(length(solution))) do j
        Base.@_inline_meta
        lhsϕ = lhs[j]
        rhsϕ = data(rhs[j])
        ϕ = data(solution[j])

        for i in eachindex(ϕ)
            @inbounds rhsϕ[i] = ϕ[i] + Δt*rhsϕ[i]
        end

        t_ldiv!(ϕ, lhsϕ, rhsϕ)
    end
    return nothing
end

"Step forward `m` by `Δt` with the backward Euler method."
function iterate!(m::AbstractModel{TS}, Δt) where TS <: BackwardEulerTimestepper
    update!(m.bcs, m.timestepper.eqn, m.solution, m)
    calc_implicit_rhs!(m.timestepper.rhs, m.timestepper.eqn, m.solution, m)
    calc_diffusive_lhs!(m.timestepper.lhs, m.timestepper.eqn, m.solution, Δt, m)
    backward_euler_step!(m.timestepper.rhs, m.timestepper.lhs, m.solution, Δt)
    tick!(m.clock, Δt)
    return nothing
end


text:
@testset "preprocessors" begin
    obs1 = (state = [1, 2, 3],)
    p = CloneStatePreprocessor()
    obs2 = p(obs1)

    @test get_state(obs1) !== get_state(obs2)
    @test get_state(obs1) == get_state(obs2)
end


text:
using PowerSystems
using PowerSimulations
const PSI = PowerSimulations
IS = PowerSystems.IS

using Dates
using DataFrames

using Cbc # mip solver
solver = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 1, "ratioGap" => 0.5)
using Ipopt # solver that supports duals
ipopt_solver = optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)

base_dir = PowerSystems.download(PowerSystems.TestData; branch = "master");
pm_data = PowerSystems.PowerModelsData(joinpath(base_dir, "matpower", "case5_re_uc.m"))

FORECASTS_DIR = joinpath(base_dir, "forecasts", "5bus_ts", "7day")

tsp_da = IS.read_time_series_file_metadata(
    joinpath(FORECASTS_DIR, "timeseries_pointers_da_7day.json"),
)
tsp_rt = IS.read_time_series_file_metadata(
    joinpath(FORECASTS_DIR, "timeseries_pointers_rt_7day.json"),
)
tsp_agc = IS.read_time_series_file_metadata(
    joinpath(FORECASTS_DIR, "timeseries_pointers_agc_7day.json"),
)

sys_DA = System(pm_data)
reserves = [
    VariableReserve{ReserveUp}("REG1", true, 5.0, 0.1),
    VariableReserve{ReserveUp}("REG2", true, 5.0, 0.06),
    VariableReserve{ReserveUp}("REG3", true, 5.0, 0.03),
    VariableReserve{ReserveUp}("REG4", true, 5.0, 0.02),
]
contributing_devices = get_components(Generator, sys_DA)
for r in reserves
    add_service!(sys_DA, r, contributing_devices)
end

add_time_series!(sys_DA, tsp_da)
transform_single_time_series!(sys_DA, 48, Hour(24))

sys_RT = System(pm_data)
add_time_series!(sys_RT, tsp_rt)
transform_single_time_series!(sys_RT, 12, Hour(1))

sys_AGC = System(pm_data)
add_time_series!(sys_AGC, tsp_agc)

template_uc = template_unit_commitment()
devices = Dict(
    :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
    :Ren => DeviceModel(RenewableDispatch, RenewableFullDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    :HydroROR => DeviceModel(HydroDispatch, HydroDispatchRunOfRiver),
    :RenFx => DeviceModel(RenewableFix, FixedOutput),
)
template_ed = template_economic_dispatch(devices = devices)

stages_definition = Dict(
    "UC" => Stage(
        GenericOpProblem,
        template_uc,
        sys_DA,
        solver,
        balance_slack_variables = true,
    ),
    "ED" => Stage(
        GenericOpProblem,
        template_ed,
        sys_RT,
        ipopt_solver,
        constraint_duals = [:CopperPlateBalance],
    ),
)

feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24))

feedforward = Dict(
    ("ED", :devices, :Generators) => SemiContinuousFF(
        binary_source_stage = PSI.ON,
        affected_variables = [PSI.ACTIVE_POWER],
    ),
)

cache = Dict("UC" => [TimeStatusChange(ThermalStandard, PSI.ON)])

order = Dict(1 => "UC", 2 => "ED")
horizons = Dict("UC" => 48, "ED" => 12)
intervals = Dict("UC" => (Hour(24), Consecutive()), "ED" => (Hour(1), Consecutive()))

DA_RT_sequence = SimulationSequence(
    step_resolution = Hour(24),
    order = order,
    horizons = horizons,
    intervals = intervals,
    ini_cond_chronology = InterStageChronology(),
    feedforward_chronologies = feedforward_chronologies,
    feedforward = feedforward,
)

file_path = mkpath(joinpath(".","5-bus-simulation"))
sim = Simulation(
    name = "5bus-test",
    steps = 1,
    stages = stages_definition,
    stages_sequence = DA_RT_sequence,
    simulation_folder = file_path,
)

build!(sim)

execute!(sim)

# Results

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
module VegaGraphs

using VegaLite
using DataStructures
using LightGraphs, GraphPlot, SimpleWeightedGraphs, LinearAlgebra, DataFrames

export graphplot
include("graph_functions.jl")
include("plotting.jl")

"""
graphplot(df, layout = GraphPlot.spring_layout; kwargs...)

Input:
* g is SimpleGraph of SimpleWeightedGraph.
* layout is a GraphPlot layout function (e.g GraphPlot.circular_layout)

This function produces an standard graph plot.
The possible arguments to be used are the following:
* node_label       = true
* node_labelsize   = 10
* tooltip          = true
* node_size        = 500
* node_color       = "#9ecae9"
* node_sizefield   = nothing
* node_colorfield  = nothing
* node_colorscheme = "blues"
* node_opacity     = 1.0
* edge_opacity     = 0.5
* edge_weightfield = nothing
* width            = 600
* height           = 400
"""
# function graphplot(g::union{simpleweightedgraph,simplegraph},
function graphplot(g::AbstractGraph,
                layout = GraphPlot.spring_layout;node_labels=nothing,kwargs...)
    if node_labels == nothing
        node_labels = collect(vertices(g))
    end
    node_x, node_y = layout(g)
    graph_edges = get_edges(g,node_x,node_y);
    graph_nodes = DataFrame(keywords = node_labels,node_x=node_x,node_y=node_y);

    return vl_graph_plot(graph_nodes,graph_edges;kwargs...)
end

end


text:
export AFM

struct Interface
		  uϵ::Symbol
		  uσ::Symbol
		  from_ϵσ::Function
		  to_ϵσ::Function
	   end

function Base.getindex(d::RheoTimeData, i::Interface)
		  return( NamedTuple{(i.uϵ, i.uσ),Tuple{Vector{RheoFloat},Vector{RheoFloat}}}( i.from_ϵσ(d.ϵ,d.σ) )   )
	   end

function RheoTimeData(interface::Interface ; t::Vector{T3} = RheoFloat[], comment="Created from generic constructor", savelog = true, log = savelog ? RheoLogItem(comment) : nothing, kwargs...)  where {T1<:Real, T2<:Real, T3<:Real}

   if interface.uϵ in keys(kwargs)
	  uϵ = convert(Vector{RheoFloat}, kwargs[interface.uϵ])
   else
	  uϵ = RheoFloat[]
   end

   if interface.uσ in keys(kwargs)
	  uσ = convert(Vector{RheoFloat}, kwargs[interface.uσ])
   else
	  uσ = RheoFloat[]
   end

  typecheck = check_time_data_consistency(t,uϵ,uσ)
  RheoTimeData(convert(Vector{RheoFloat},σ), convert(Vector{RheoFloat},ϵ), convert(Vector{RheoFloat},t),
  log == nothing ? nothing : [ RheoLogItem(log.action,merge(log.info, (type=typecheck, interface = interface)))]     )

end

function importcsv(filepath::String, interface::Interface; t_col::IntOrNone = nothing, delimiter = ',', comment = "Imported from csv file", savelog = true, kwargs...)
   cols = kwargs.data   # should contain column numbers for the strain and stress equivalent data
   uϵ_col_sym = Symbol(string(interface.uϵ) * "_col")
   uσ_col_sym = Symbol(string(interface.uσ) * "_col")

   if uϵ_col_sym in keys(cols)
	  uϵ_col = cols[uϵ_col_sym]
   else
	  uϵ_col = nothing
   end

   if uσ_col_sym in keys(cols)
	  uσ_col = cols[uσ_col_sym]
   else
	  uσ_col = nothing
   end

   data = importcsv(filepath, t_col = t_col, σ_col = uσ_col, ϵ_col = uϵ_col, delimiter = delimiter, comment = comment, savelog = savelog)

   # need to make sure that missing columns are properly transfered.
   ϵ,σ = interface.to_ϵσ(data.ϵ, data.σ)

	 log = if savelog
			 info = (comment=comment, folder=pwd(), stats=(t_min=data.t[1],t_max=data.t[end], n_sample=size(data.t[:])))
			 kwds = NamedTuple{(:t_col, uϵ_col_sym, uσ_col_sym),Tuple{IntOrNone,IntOrNone,IntOrNone}}( (t_col, uϵ_col, uσ_col) )
			 RheoLogItem( (type=:source, funct=:importcsv, params=(filepath=filepath, interface=interface), keywords=kwds), info )
		   else
			 RheoLogItem(type = nothing, info = "nothing")
		   end

   return RheoTimeData(σ, ϵ, data.t, [log])

end

function AFM(R::Real)
   Interface(:d, :f, (ϵ,σ)->(d = (R*(3. /4.).^(2. /3)) .* ϵ.^(2. /3), f = R.^2 .* σ ), (d,f)->(ϵ = 4. /(3. *R^1.5) .* d.^1.5, σ = f./R^2) )
end

function Tweezers(R::Real)
   Interface(:d, :f, (ϵ,σ)->(d = (R) * ϵ, f = R^2 * σ ), (d,f)->(ϵ = 1 .*(d/R)^1.5, σ = f/R^2) )
end


text:
# cachsim output
export linear_solve_cachesim
# import ..Femshop: init_cachesimout, add_cachesim_array, cachesim_load, cachesim_store, cachesim_load_range, cachesim_store_range

function linear_solve_cachesim(var, bilinear, linear, stepper=nothing)
    # if config.linalg_matrixfree
    #     return solve_matrix_free_sym(var, bilinear, linear, stepper);
    #     #return solve_matrix_free_asym(var, bilinear, linear, stepper);
    # end
    # If more than one variable
    if typeof(var) <: Array
        # multiple variables being solved for simultaneously
        dofs_per_node = 0;
        dofs_per_loop = 0;
        for vi=1:length(var)
            dofs_per_loop += length(var[vi].symvar);
            dofs_per_node += var[vi].total_components;
        end
    else
        # one variable
        dofs_per_loop = length(var.symvar);
        dofs_per_node = var.total_components;
    end
    N1 = size(grid_data.allnodes,2);
    Nn = dofs_per_node * N1;
    Np = refel.Np;
    nel = mesh_data.nel;

    init_cachesimout(N1, refel, nel, dofs_per_node, variables);

    if prob.time_dependent && !(stepper === nothing)
        #TODO time dependent coefficients
        assemble_t = @elapsed((A, b) = assemble_cachesim(var, bilinear, linear, 0, stepper.dt));
        log_entry("Assembly took "*string(assemble_t)*" seconds");

        log_entry("Beginning "*string(stepper.Nsteps)*" time steps.");
        t = 0;
        sol = [];
        start_t = Base.Libc.time();
        for i=1:stepper.Nsteps
            b = assemble_rhs_only_cachesim(var, linear, t, stepper.dt);
            #sol = A\b;
            cachesim_load_range(1);
            cachesim_load_range(2);

            t += stepper.dt;
        end
        end_t = Base.Libc.time();

        log_entry("Cachesim solve took "*string(end_t-start_t)*" seconds");
        return 0;

    else
        assemble_t = @elapsed((A, b) = assemble_cachesim(var, bilinear, linear));
        #sol_t = @elapsed(sol = A\b);
        cachesim_load_range(1);
        cachesim_load_range(2);

        log_entry("Assembly took "*string(assemble_t)*" seconds");

        return 0;
    end
end

function assemble_cachesim(var, bilinear, linear, t=0.0, dt=0.0)
    Np = refel.Np;
    nel = mesh_data.nel;
    N1 = size(grid_data.allnodes,2);
    multivar = typeof(var) <: Array;
    if multivar
        # multiple variables being solved for simultaneously
        dofs_per_node = 0;
        dofs_per_loop = 0;
        for vi=1:length(var)
            dofs_per_loop += length(var[vi].symvar);
            dofs_per_node += var[vi].total_components;
        end
    else
        # one variable
        dofs_per_loop = length(var.symvar);
        dofs_per_node = var.total_components;
    end
    Nn = dofs_per_node * N1;

    b = 0;
    A = 0;
    AI = 0;
    AJ = 0;
    AV = 0;

    allnodes_id = add_cachesim_array(size(grid_data.allnodes),8);
    ak_id = add_cachesim_array(Np*Np*dofs_per_node*dofs_per_node,8);
    aktmp_id = add_cachesim_array(Np*Np*dofs_per_node*dofs_per_node,8);
    bk_id = add_cachesim_array(Np*dofs_per_node,8);
    bktmp_id = add_cachesim_array(Np*dofs_per_node,8);

    #  Elemental loop follows elemental ordering
    for e=elemental_order;
        gis = grid_data.glbvertex[:,e];
        vx = grid_data.allnodes[:,gis];         # coordinates of element's vertices
        glb = grid_data.loc2glb[:,e];                 # global indices of this element's nodes for extracting values from var arrays
        xe = grid_data.allnodes[:,glb[:]];  # coordinates of this element's nodes for evaluating coefficient functions
        cachesim_load_range(allnodes_id, gis[:], 1:refel.dim);
        cachesim_load_range(allnodes_id, glb[:], 1:refel.dim);

        cachesim_load_range(aktmp_id);
        cachesim_load_range(bktmp_id);

        Astart = (e-1)*Np*dofs_per_node*Np*dofs_per_node + 1; # The segment of AI, AJ, AV for this element

        # The linear part. Compute the elemental linear part for each dof
        rhsargs = (var, e, 0, grid_data, geo_factors, refel, t, dt, 0, 0);
        lhsargs = (var, e, 0, grid_data, geo_factors, refel, t, dt, 0, 0);
        if dofs_per_node == 1
            linchunk = linear.func(rhsargs);  # get the elemental linear part
            #b[glb] .+= linchunk;
            cachesim_load_range(2, glb);
            cachesim_load_range(bk_id);
            cachesim_store_range(2, glb);

            bilinchunk = bilinear.func(lhsargs); # the elemental bilinear part
            #A[glb, glb] .+= bilinchunk;         # This will be very inefficient for sparse A
            Arange = Astart:(Astart + Np*Np*dofs_per_node*dofs_per_node);
            cachesim_load_range(1, Arange, 1:3);
            cachesim_load_range(ak_id);
            cachesim_store_range(1, Arange, 1:3);

        elseif typeof(var) == Variable
            # only one variable, but more than one dof
            linchunk = linear.func(rhsargs);
            insert_linear_cachesim!(bk_id, linchunk, glb, 1:dofs_per_node, dofs_per_node);

            bilinchunk = bilinear.func(lhsargs);
            insert_bilinear_cachesim!(ak_id, Astart, bilinchunk, glb, 1:dofs_per_node, dofs_per_node);
        else
            linchunk = linear.func(rhsargs);
            insert_linear_cachesim!(bk_id, linchunk, glb, 1:dofs_per_node, dofs_per_node);

            bilinchunk = bilinear.func(lhsargs);
            insert_bilinear_cachesim!(ak_id, Astart, bilinchunk, glb, 1:dofs_per_node, dofs_per_node);
        end
    end

    cachesim_load_range(1);

    # Boundary conditions are not considered

    return (A, b);
end

# assembles the A and b in Au=b
function assemble_rhs_only_cachesim(var, linear, t=0.0, dt=0.0)
    Np = refel.Np;
    nel = mesh_data.nel;
    N1 = size(grid_data.allnodes,2);
    multivar = typeof(var) <: Array;
    if multivar
        # multiple variables being solved for simultaneously
        dofs_per_node = 0;
        var_to_dofs = [];
        for vi=1:length(var)
            tmp = dofs_per_node;
            dofs_per_node += length(var[vi].symvar);
            push!(var_to_dofs, (tmp+1):dofs_per_node);
        end
    else
        # one variable
        dofs_per_node = length(var.symvar);
    end
    Nn = dofs_per_node * N1;

    #b = zeros(Nn);

    for e=1:nel;
        glb = grid_data.loc2glb[:,e];                 # global indices of this element's nodes for extracting values from var arrays
        xe = grid_data.allnodes[:,glb[:]];  # coordinates of this element's nodes for evaluating coefficient functions

        rhsargs = (var, e, 0, grid_data, geo_factors, refel, t, dt, 0, 0);

        #linchunk = linear.func(args);  # get the elemental linear part
        if dofs_per_node == 1
            linchunk = linear.func(rhsargs);  # get the elemental linear part
            #b[glb] .+= linchunk;

        elseif typeof(var) == Variable
            # only one variable, but more than one dof
            linchunk = linear.func(rhsargs);
            insert_linear_cachesim!(bk_id, linchunk, glb, 1:dofs_per_node, dofs_per_node);

        else
            linchunk = linear.func(rhsargs);
            insert_linear_cachesim!(bk_id, linchunk, glb, 1:dofs_per_node, dofs_per_node);

        end
    end

    return b;
end

# Inset the single dof into the greater construct
function insert_linear_cachesim!(bk_id, bel, glb, dof, Ndofs)
    # group nodal dofs
    for d=1:length(dof)
        ind = glb.*Ndofs .- (Ndofs-dof[d]);
        ind2 = ((d-1)*length(glb)+1):(d*length(glb));

        #b[ind] = b[ind] + bel[ind2];

        cachesim_load_range(2,ind);
        cachesim_load_range(bk_id,ind2);
        cachesim_store_range(2,ind);
    end
end

function insert_bilinear_cachesim!(ak_id, Astart, ael, glb, dof, Ndofs)
    Np = length(glb);
    # group nodal dofs
    for dj=1:length(dof)
        indj = glb.*Ndofs .- (Ndofs-dof[dj]);
        indj2 = ((dj-1)*Np+1):(dj*Np);
        for di=1:length(dof)
            indi = glb.*Ndofs .- (Ndofs-dof[di]);
            indi2 = ((di-1)*Np+1):(di*Np);

            for jj=1:Np
                offset = Astart + (jj-1 + Np*(dj-1))*Np*Ndofs + Np*(di-1) - 1;
                Arange = (offset+1):(offset + Np);
                cachesim_load_range(1, Arange, 1:3);
                cachesim_load_range(ak_id, indi2, indj2);
                cachesim_store_range(1, Arange, 1:3);
            end
        end
    end

end

max_stars_repo_name: paralab/femshop

max_issues_repo_name: paralab/femshop

max_forks_repo_name: paralab/femshop


text:
struct PipelineBindingState
    pipeline::Pipeline
    push_ranges::Vector{PushConstantRange}
end

PipelineBindingState(pipeline) = PipelineBindingState(pipeline, [])

"""
Binding state that must be set in order for
drawing commands to render correctly.
"""
struct BindRequirements
    dependencies::ShaderDependencies
    pipeline_state::PipelineBindingState
    layout::PipelineLayout
    push_data::Any
end

BindRequirements(dependencies, pipeline_state, layout) = BindRequirements(dependencies, pipeline_state, layout, nothing)

"""
Describes the current binding state.
"""
struct BindState
    vertex_buffer::Optional{VertexBuffer}
    index_buffer::Optional{IndexBuffer}
    descriptor_sets::Vector{DescriptorSet}
    pipeline_state::Optional{PipelineBindingState}
    push_data::Any
end

BindState() = BindState(nothing, nothing, [], nothing, nothing)

function Base.bind(cbuffer::CommandBuffer, reqs::BindRequirements, state::BindState = BindState())
    @unpack vertex_buffer, index_buffer, descriptor_sets = reqs.dependencies
    @unpack push_ranges, pipeline = reqs.pipeline_state
    current_push_data = state.push_data
    layout = reqs.layout

    if pipeline ≠ state.pipeline_state.pipeline
        cmd_bind_pipeline(cbuffer, PIPELINE_BIND_POINT_GRAPHICS, pipeline)
        if push_ranges ≠ state.pipeline_state.push_ranges
            # push data gets invalidated
            current_push_data = nothing
        end
    end

    vertex_buffer ≠ state.vertex_buffer && cmd_bind_vertex_buffers(cbuffer, [vertex_buffer], [0])

    if !isnothing(index_buffer) && index_buffer ≠ state.index_buffer
        cmd_bind_index_buffer(cbuffer, index_buffer, 0, INDEX_TYPE_UINT32)
    end

    if !isempty(descriptor_sets)
        dsets = handle.(descriptor_sets)
        dsets ≠ state.descriptor_sets && cmd_bind_descriptor_sets(cbuffer, PIPELINE_BIND_POINT_GRAPHICS, layout, 0, dsets, [])
    end

    if !isnothing(push_data) && push_data ≠ current_push_data
        cmd_push_constants(cbuffer, layout, SHADER_STAGE_VERTEX_BIT, 1, Ref(push_data), sizeof(push_data))
    end

    BindState(vertex_buffer, index_buffer, descriptor_sets, PipelineBindingState(pipeline, push_ranges), push_data)
end


text:
# The BFGS update
# The inverse can either be updated with
#  H = (I-s*y'/y'*s)*H*(I-y*s'/y'*s)/(s'y)^2 - (H*y*s' + s*y'*H)/s'*y
# or
#  H = H + (s'*y + y'*H*y)*(s*s')/(s'*y)^2 - (B*y*s' + s*y'*H)/(s'*y)
# The direct update can be written as
#  B = B = yy'/y's-Bss'B'/s'B*s

struct BFGS{T1} <: QuasiNewton{T1}
   approx::T1
end
BFGS() = BFGS(Inverse())
# function update!(scheme::BFGS, B::Inverse, Δx, y)
#    B.A = B.A + y*y'/dot(Δx, y) - B.A*y*y'*B.A/(y'*B.A*y)
# end

function update!(scheme::BFGS{<:Direct}, B, s, y)
    # We could write this as
    #     B .+= (y*y')/dot(s, y) - (B*s)*(s'*B)/(s'*B*s)
    #     B .+= (y*y')/dot(s, y) - b*b'/dot(s, b)
    # where b = B*s
    # But instead, we split up the calculations. First calculate the denominator
    # in the first term
    σ = dot(s, y)
    ρ = inv(σ) # scalar
    # Then calculate the vector b
    b = B*s # vector temporary
    sBs = dot(s, b)
    # Calculate one vector divided by dot(s, b)
    ρbb = inv(sBs)*b
    # And calculate
    B .+= (ρ*y)*y' .- ρbb*b'
end
function update(scheme::BFGS{<:Direct}, B, s, y)
   # As above, but out of place
   σ = dot(s, y)
   ρ = inv(σ)
   b = B*s
   ρbb = inv(dot(s,b))*b
   B + (ρ*y)*y' - ρbb*b'
end
function update(scheme::BFGS{<:Inverse}, H, s, y)
   σ = dot(s, y)
   ρ = inv(σ)

   C = (I - ρ*s*y')
   H = C*H*C' + ρ*s*s'

   H
end
function update!(scheme::BFGS{<:Inverse}, H, s, y)
   σ = dot(s, y)
   ρ = inv(σ)

   if isfinite(ρ)
      Hy = H*y
      H .= H .+ ((σ+y'*Hy).*ρ^2)*(s*s')
      Hys = Hy*s'
      Hys .= Hys .+ Hys'
      H .= H .- Hys.*ρ
   end
   H
end
# end
# function update!(H, s, y, scheme::BFGS{<:Inverse})
#    sy = dot(s, y)
#    ρ = inv(sy)
#
#    if isfinite(ρ)
#       C = (I - ρ.*s*y')
#       H .= C*H*C' + ρ*s*s'
#    end
#    H
# end

function update!(scheme::BFGS{<:Inverse}, A::UniformScaling, s, y)
   update(scheme, A, s, y)
end
function update!(scheme::BFGS{<:Direct}, A::UniformScaling, s, y)
   update(scheme, A, s, y)
end


text:
# Test these sundials_sundials, sundials_sunlinsolve, sundials_sunlinsolve, sundial_sunmatrix are callable
@test Sundials.SUNLinSolNewEmpty() isa Any


text:
module AStarSearch
export astar, depthfirst, breadthfirst, iterative_deepening

include("common.jl")
include("uninformed_search.jl")
include("a_star.jl")

end # module


text:
using Colors

include("helpers/page.jl")
include("helpers/doc.jl")

intro = md"""
$(title(2, "Widgets"))
$(vskip(1em))

# Functions
$(vskip(1em))

$(
    showdocs([button, slider, checkbox, togglebutton,
    radio, radiogroup, textinput, progress, paper, datepicker, codemirror])
)
"""

function main(window)
    docpage(intro)
end