text:
# ---
# title: 874. Walking Robot Simulation
# id: problem874
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Greedy
# link: <https://leetcode.com/problems/walking-robot-simulation/description/>
# hidden: true
# ---
#
# A robot on an infinite grid starts at point (0, 0) and faces north.  The robot
# can receive one of three possible types of commands:
#
#   * `-2`: turn left 90 degrees
#   * `-1`: turn right 90 degrees
#   * `1 <= x <= 9`: move forward `x` units
#
# Some of the grid squares are obstacles.
#
# The `i`-th obstacle is at grid point `(obstacles[i][0], obstacles[i][1])`
#
# If the robot would try to move onto them, the robot stays on the previous grid
# square instead (but still continues following the rest of the route.)
#
# Return the **square** of the maximum Euclidean distance that the robot will be
# from the origin.
#
#
#
# **Example 1:**
#
#
#
#     Input: commands = [4,-1,3], obstacles = []
#     Output: 25
#     Explanation: robot will go to (3, 4)
#
#
# **Example 2:**
#
#
#
#     Input: commands = [4,-1,4,-2,4], obstacles = [[2,4]]
#     Output: 65
#     **Explanation** : robot will be stuck at (1, 4) before turning left and going to (1, 8)
#
#
#
#
# **Note:**
#
#   1. `0 <= commands.length <= 10000`
#   2. `0 <= obstacles.length <= 10000`
#   3. `-30000 <= obstacle[i][0] <= 30000`
#   4. `-30000 <= obstacle[i][1] <= 30000`
#   5. The answer is guaranteed to be less than `2 ^ 31`.
#
#
## @lc code=start
using LeetCode

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


text:
using Revise
using MccDaqHats
using Infiltrator
includet(joinpath(@__DIR__, "scan_utils.jl"))

# Constants
const DEVICE_COUNT = 2
const MASTER = 0
const CURSOR_SAVE = "\x1b[s"
const CURSOR_RESTORE = "\x1b[u"
const CURSOR_BACK_2 = "\x1b[2D"
const ERASE_TO_END_OF_LINE = "\x1b[0K"

"""
multi_hat_synchronous_scan()
        MCC 172 Functions Demonstrated:
        mcc172_trigger_config
        mcc172_a_in_clock_config_write
        mcc172_a_in_clock_config_read
        mcc172_a_in_scan_start
        mcc172_a_in_scan_status
        mcc172_a_in_scan_read
        mcc172_a_in_scan_stop

    Purpose:
        Get synchronous data from multiple MCC 172 devices.

    Description:
        This example demonstrates acquiring data synchronously from multiple
        MCC 172 devices.  This is done using the shared clock and trigger
        options.  An external trigger source must be provided to the TRIG
        terminal on the master MCC 172 device.  The clock and trigger on the
        master device are configured for SOURCE_MASTER and the remaining devices
        are configured for SOURCE_SLAVE.
Enahancements:
Deinterleave data using deinterleave()
Save data use one of the libraries in github.com/juliaIO, use similar format to matlab format
Check if trigger will work for synchronizing
"""
function multi_hat_synchronous_scan()
    hats = HatInfo[]

    # Define the channel list for each HAT device
    chans = Vector{Integer}[]
    push!(chans, [0, 1])
    push!(chans, [0, 1])

    # Define the options for each HAT device
    options = [Set([:OPTS_EXTTRIGGER]);       # first hat
                    Set([:OPTS_EXTTRIGGER])]  # second hat
    samples_per_channel = 10000
    sample_rate = 10240.0  # Samples per second
    trigger_mode = :TRIG_RISING_EDGE

    try
        # Get an instance of the selected hat device.
        hats = select_hat_devices(:MCC_172, DEVICE_COUNT)

        # Validate the selected channels.
        for (i, hat) in enumerate(hats)
            validate_channels(chans[i], mcc172_info().NUM_AI_CHANNELS)
        end

        # Turn on IEPE supply?
        iepe_enable = get_iepe()

        for (i, hat) in enumerate(hats)
            mcc172_open(hat.address)
            for channel in chans[i]
                # Configure IEPE.
                mcc172_iepe_config_write(hat.address, channel, iepe_enable)
            end

            if hat.address != MASTER
                # Configure the slave clocks.
                mcc172_a_in_clock_config_write(hat.address, :SOURCE_SLAVE, sample_rate)
                # Configure the trigger.
                mcc172_trigger_config(hat.address, :SOURCE_SLAVE, trigger_mode)
            end
        end

        # Configure the master clock and start the sync.
        mcc172_a_in_clock_config_write(MASTER, :SOURCE_MASTER, sample_rate)
        synced = false
        actual_rate = 0
        while !synced
            _source_type, actual_rate, synced = mcc172_a_in_clock_config_read(MASTER)
            if !synced
                sleep(0.005)
            end
        end

        # Configure the master trigger.
        mcc172_trigger_config(MASTER, :SOURCE_MASTER, trigger_mode)

        println("MCC 172 multiple HAT example using external clock and external trigger options")
        println("    Functions demonstrated:")
        println("         mcc172_trigger_mode")
        println("         mcc172_a_in_clock_config_write")
        println("         mcc172_a_in_clock_config_read")
        println("         mcc172_a_in_scan_start")
        println("         mcc172_a_in_scan_read")
        println("         mcc172_a_in_scan_stop")
        println("         mcc172_a_in_scan_cleanup")
        println("    IEPE power: $(iepe_enable ? "on" : "off")")
        println("    Samples per channel:", samples_per_channel)
        println("    Requested Sample Rate: $(round(sample_rate, digits=3))")
        println("    Actual Sample Rate: $(round(actual_rate, digits=3))")
        println("    Trigger type: $trigger_mode")

        for (i, hat) in enumerate(hats)
            println("    HAT: $i")
            println("      Address: $(hat.address)")
            println("      Channels: $(join(chans, ", "))")
            # options_str = enum_mask_to_string(OptionFlags, options[i])
            println("      Options: $options")
        end

        # determine if internal (gpio) or external trigger
        trigger_dialog()

        # Start the scan.
        for (i, hat) in enumerate(hats)
            # @show(chans[i])
            chan_mask = chan_list_to_mask(chans[i])
            mcc172_a_in_scan_start(hat.address, chan_mask, UInt32(samples_per_channel), options[i])
        end

        println("\nWaiting for trigger ... Press Ctrl-C to stop scan\n")
        trigger_dialog()

        try
            # Monitor the trigger status on the master device.
            wait_for_trigger(MASTER)
            # Read and display data for all devices until scan completes
            # or overrun is detected.
            read_and_display_data(hats, chans)

        catch e
            if isa(e, InterruptException)  #KeyboardInterrupt "^C"
                # Clear the "^C" from the display.
                println("$CURSOR_BACK_2 $ERASE_TO_END_OF_LINE \nAborted\n")
            else
                println("\n $e")
            end
        end

    finally
        for (i, hat) in enumerate(hats)
            mcc172_a_in_scan_stop(hat.address)
            mcc172_a_in_scan_cleanup(hat.address)
        end
    end
end

"""
    function read_and_display_data(hats, chans::Vector{Vector{Integer}})

Reads data from the specified channels on the specified DAQ HAT devices
and updates the data on the terminal display.  The reads are executed in a
loop that continues until either the scan completes or an overrun error
is detected.

Args:
hats (list[mcc172]): A list of mcc172 HAT device objects.
chans (list[int][int]): A 2D list to specify the channel list for each
mcc172 HAT device.

Returns:
None
"""
function read_and_display_data(hats::Vector{HatInfo}, chans::Vector{Vector{Integer}})
    # assume two channels per HAT
    num_chan = 2
    samples_to_read = 1000
    timeout = 5  # Seconds
    samples_per_chan_read = zeros(Integer, DEVICE_COUNT)
    total_samples_per_chan = zeros(Integer, DEVICE_COUNT)
    is_running = true

    # Since the read_request_size is set to a specific value, a_in_scan_read()
    # will block until that many samples are available or the timeout is
    # exceeded.

    # Create blank lines where the data will be displayed
    for _ in 1:(DEVICE_COUNT * 4 + 1)
        println()
    end

    # Move the cursor up to the start of the data display.
    print("\x1b[$(DEVICE_COUNT * 4 + 1)A")
    print(CURSOR_SAVE)

    while true
        data = Matrix{Float64}(undef, samples_to_read*num_chan, DEVICE_COUNT)
        # Read the data from each HAT device.
        for (i, hat) in enumerate(hats)
            if num_chan != length(chans[i])
                error("Expecting $num_chan channels, got $(length(chans[i])) for hat $i")
            end
            resultcode, statuscode, result, samplesread =
                mcc172_a_in_scan_read(hat.address, Int32(samples_to_read), num_chan, timeout)
            data[:,i] = result
            status = mcc172_status_decode(statuscode)
            is_running &= status.running  # use to break out of while loop
            samples_per_chan_read[i] = length(result) / num_chan
            total_samples_per_chan[i] += samples_per_chan_read[i]
            # @show(length(result), samples_per_chan_read[i], total_samples_per_chan[i])

            if status.bufferoverrun
                print("\nError: Buffer overrun")
                break
            elseif status.hardwareoverrun
                print("\nError: Hardware overrun")
                break
            end
        end

        print(CURSOR_RESTORE)

        # Display the data for each HAT device
        for (i, hat) in enumerate(hats)
            print("HAT $i")

            # Print the header row for the data table.
            print("  Samples Read    Scan Count")
            for chan in chans[i]
                print("     Channel $chan")
            end
            println()

            # Display the sample count information.
            print("          $(samples_per_chan_read[i])            $(total_samples_per_chan[i])     ")

            # Display the data for all selected channels
            #for chan_idx in range(len(chans[i])):
            #    sample_idx = ((samples_per_chan_read[i] * len(chans[i]))
            #                  - len(chans[i]) + chan_idx)
            #    println("[:>12.5f} V".format(data[i][sample_idx]), end="")

            # Display the RMS voltage for each channel.
            # can use deinterleave to separate interleaved data.
            if samples_per_chan_read[i] > 0
                for channel in chans[i]
                    value = calc_rms(data[:,i], channel, length(chans[i]), samples_per_chan_read[i])
                    print("  $(round(value, digits=5)) Vrms ")
                end
                # stdout.flush()
            end
            print("\n")
        end

        if !is_running
            break
        end
    end
end

ext: jl

max_stars_repo_name: Spectrum-Tec/MccDaqHats

max_issues_repo_name: Spectrum-Tec/MccDaqHats

max_forks_repo_name: Spectrum-Tec/MccDaqHats


text:
# Loading of Matthias' Peninsula data
#
# big ones: 49, 197, 281, 402, 422, 428

function LoadPara(gid::PeninsulaGlacier; pl_kwargs...)
    glacier = getid(gid)
    pl_kwargs = deepcopy(pl_kwargs)
    pl = LoadPara(; gid=gid, data_root=joinpath(LoadPara().data_root,
                                          "Peninsula"), pl_kwargs...)
    for (k,v) in (ThicknessData=>10.0,
                  BdotData=>0.2,
                  DhdtData=>0.2,
                  IVData=>1.0)
        pl.dataset_opts[k][:sigma]=v
    end
    updated_pm_defaults = Dict(
        :bandsize => 30,
        :window_width_smooth=>100, # can make a noisy tau
        :cap_at_floatation=>false,
        :dist_exp => 0.5,
        :h_at_top => -1,
    )

    # update pm but keeping any explicitly set (via kwargs) options
    merge!((a,b)->a, pl.dataset_opts[ParaData][:pm], updated_pm_defaults)

    # Set loaders & files
    id = getid(gid)
    id5 = VAWTools.int2str5(id)
    pl = deepcopy(pl)

    dl = pl.dataset_LOADERS
    fl = pl.dataset_files
    # Set the LOADERS
    dl[OutlineData] = PeninsulaLoader
    if isfile(joinpath(pl.data_root, "outline/outcrops/$id5.xyn"))
        fl[OutlineData] = [joinpath(pl.data_root, "outline/$id5.xyn"),
                           joinpath(pl.data_root, "outline/outcrops/$id5.xyn")]
    else
        fl[OutlineData] = [joinpath(pl.data_root, "outline/$id5.xyn")]
    end
    dl[ParaData] = PeninsulaLoader

    dl[DEMData] = PeninsulaLoader
    fl[DEMData] = [joinpath(pl.data_root, "dem/APDEM100m.bin")]

    dl[IVData] = PeninsulaLoader
    fl[IVData] = [joinpath(pl.data_root, "velocity/speed_measured.grid")]

    dl[BdotData] = PeninsulaLoader
    fl[BdotData] = [joinpath(pl.data_root, "RACMO/downscale_RACMO_nn.bin")]

    dl[DhdtData] = PeninsulaLoader

    dl[ThicknessData] = PeninsulaLoader
    f = joinpath(pl.data_root, "icebridge/single_glacier/$id5.dat")
    if isfile(f)
        fl[ThicknessData] = [f]
    end

    dl[FslData] = SyntheticBenchLoader
    dl[TempData] = SyntheticBenchLoader
    dl[TerminusFluxData] = SyntheticBenchLoader

    return pl
end

function load!(ds::DataSet{OutlineData, PeninsulaLoader}, gid)
    iscalving = true
    isicecap = false
    proj = ""

    # outlines & outcrops
    outline, splits = VAWTools.concat_poly(read_xyn(ds.files[1]))

    if length(ds.files)==2
        outcrops, splits_oc = VAWTools.concat_poly(read_xyn(ds.files[2], fix=true))

        return Outline(gid, outline, splits, outcrops, splits_oc),
          iscalving::Bool, isicecap::Bool, proj::String
    else
        return Outline(gid, outline, splits),
          iscalving::Bool, isicecap::Bool, proj::String
    end
end

function load!(ds::DataSet{DEMData, PeninsulaLoader}, outline::Outline, pp::Phys
               )::Tuple{Gridded, Matrix{Bool}, Matrix{Bool}, Matrix{F}, F, Tuple, String}
    load!(convert(DataSet{DEMData, :ASCIIGridLoader}, ds), outline, pp)
end
function generate_DEM_masks(gid::PeninsulaGlacier, dem, outline)
    if hasoutcrops(outline)
        glaciermask, landmask = maskit(dem, outline, outline.outcrops)
    else
        glaciermask, _ = maskit(dem, outline)
        landmask = convert(Matrix{Bool}, falses(size(glaciermask)))
    end
    return glaciermask, landmask
end

function load!(ds::DataSet{IVData, PeninsulaLoader}, dem, glaciermask, target_proj,
               outline, cropbox, pp::Phys)::Tuple{SomeData, Matrix{Bool}}
    load!(convert(DataSet{IVData, :ASCIIGridLoader}, ds),
          dem, glaciermask, target_proj, outline, cropbox, pp)
end

function load!(ds::DataSet{ThicknessData, PeninsulaLoader}, dem, glaciermask, target_proj, outline, cropbox,
               pp::Phys)::SomeData
    @unpack files, opts, daterange = ds
    sigma = opts[:sigma]
    if isfile(files[1])
        out = get_cache_ram!(ds, :thick, ()->readdlm(files[1]))
    else
        return SomeData(nothing)
    end

    if size(out,2)==5
        out = out[out[:,5].==0,:] # drop wall reflection points
    end
    tr = Traj(out[:,1], out[:,2], out[:,4], out[:,3])
    split_traj!(tr, 1e3)
    return SomeData(vals=tr, sigma=sigma)
end

function load!(ds::DataSet{BdotData, PeninsulaLoader}, dem, glaciermask, target_proj, outline, cropbox,
               pp::Phys)
    bdot = load_asciigrid!(ds, :bdot, cropbox)

    # fill in any gaps inside the glacier
    fillgaps!(bdot, glaciermask)

    # scale from mm w.e./a to m ice / a
    scale_m_ice!(bdot, pp, 1/1000)

    return SomeData(vals=bdot, sigma=ds.opts[:sigma])
end

"This is a synthtetic dh/dt as in H&F 2014"
function load!(ds::DataSet{DhdtData, PeninsulaLoader},
               bdot, terminus_flux,
               dem, glaciermask, target_proj, outline, cropbox,
               iscalving, pp::Phys)
    opts = ds.opts
    cutoff = get(opts, :cutoff, 700)
    base = get(opts, :base, -1)
    sigma = opts[:sigma]

    # TODO: larger dh/dt for glaciers feeding into Larson A & B
    out = (base .- dem.v .* base./cutoff) # m ice/a
    out[dem.v.>cutoff] = 0
    return SomeData(vals=Gridded(dem.x, dem.y, out),
                    sigma=1.0)
end

function load!(ds::DataSet{FslData, PeninsulaLoader}, dem, glaciermask, target_proj, outline, cropbox,
               pp::Phys)
    error("FSL data not provided in Peninsula.  Use different loader, probably syntetic")
end
function load!(ds::DataSet{TempData, PeninsulaLoader}, dem, glaciermask, target_proj, outline, cropbox,
               pp::Phys)
    error("Ice-temp data not provided in Peninsula.  Use different loader, probably syntetic")
end
function load!(ds::DataSet{TerminusFluxData, PeninsulaLoader}, dem, glaciermask, iscalving,
               pp::Phys)
    error("Calving flux not provided in Peninsula.  Use different loader, probably synthetic")
end

##############
# helpers


text:
export zipped

import Base: show
import DataStructures: Deque

"""
    zipped(sources...)

Combines multiple Observables to create an Observable whose values are calculated from the values,
in order, of each of its input Observables.

# Arguments
- `sources`: input sources

# Examples
```jldoctest
using Rocket

source = zipped(of(1), from(2:5))

subscribe!(source, logger())
;

# output

[LogActor] Data: (1, 2)
[LogActor] Completed
```

```jldoctest
using Rocket

source = zipped(from(1:3), from(1:5))

subscribe!(source, logger())
;

# output

[LogActor] Data: (1, 1)
[LogActor] Data: (2, 2)
[LogActor] Data: (3, 3)
[LogActor] Completed
```

```jldoctest
using Rocket

source = zipped(completed(), of(0.0))

subscribe!(source, logger())
;

# output

[LogActor] Completed
```

See also: [`Subscribable`](@ref), [`subscribe!`](@ref)
"""
zipped()                                 = error("zipped operator expects at least one inner observable on input")
zipped(args...)                          = zipped(args)
zipped(sources::S) where { S <: Tuple }  = ZipObservable{combined_type(sources), S}(sources)
zipped(sources::V) where { V <: Vector } = ZipObservable{combined_type(sources), V}(sources)

##

struct ZipInnerActor{L, W, I} <: Actor{L}
    wrapper :: W
end

Base.show(io::IO, inner::ZipInnerActor{L, W, I}) where { L, W, I } = print(io, "ZipInnerActor($L, $I)")

on_next!(actor::ZipInnerActor{L, W, I}, data::L) where { L, W, I } = next_received!(actor.wrapper, data, Val{I}())
on_error!(actor::ZipInnerActor{L, W, I}, err)    where { L, W, I } = error_received!(actor.wrapper, err, Val{I}())
on_complete!(actor::ZipInnerActor{L, W, I})      where { L, W, I } = complete_received!(actor.wrapper, Val{I}())

##

struct ZipActorWrapper{S, A}
    actor  :: A
    values :: Deque{S}

    nsize   :: Int
    cstatus :: BitArray{1}
    vstatus :: Deque{BitArray{1}}

    subscriptions :: Vector{Teardown}

    ZipActorWrapper{S, A}(nsize::Int, actor::A) where { S, A } = begin
        values  = Deque{S}()
        cstatus = falses(nsize)
        vstatus = Deque{BitArray{1}}()

        push!(values, S())
        push!(vstatus, falses(nsize))

        subscriptions = fill!(Vector{Teardown}(undef, nsize), voidTeardown)
        return new(actor, values, nsize, cstatus, vstatus, subscriptions)
    end
end

cstatus(wrapper::ZipActorWrapper, index)       = wrapper.cstatus[index]
first_vstatus(wrapper::ZipActorWrapper, index) = first(wrapper.vstatus)[index]
last_vstatus(wrapper::ZipActorWrapper, index)  = last(wrapper.vstatus)[index]

dispose(wrapper::ZipActorWrapper) = begin fill!(wrapper.cstatus, true); foreach(s -> unsubscribe!(s), wrapper.subscriptions) end

function next_received!(wrapper::ZipActorWrapper{S}, data, index::Val{I}) where { S, I }
    cindex, cvstatus, cstorage = nothing, nothing, nothing # c - current

    for (counter, (vstatus, storage)) in enumerate(zip(wrapper.vstatus, wrapper.values))
        if vstatus[I] === false
            cindex   = counter
            cvstatus = vstatus
            cstorage = storage
            break
        end
    end

    if cindex === nothing || cvstatus === nothing || cstorage === nothing
        cindex   = length(wrapper.values) + 1
        cvstatus = falses(wrapper.nsize)
        cstorage = S()
        push!(wrapper.values, cstorage)
        push!(wrapper.vstatus, cvstatus)
    end

    setstorage!(cstorage, data, index)
    cvstatus[I] = true
    if cindex === 1 && all(cvstatus) && !all(wrapper.cstatus)
        next!(wrapper.actor, snapshot(cstorage))
        popfirst!(wrapper.values)
        popfirst!(wrapper.vstatus)

        if isempty(wrapper.values) || isempty(wrapper.vstatus)
            if any(wrapper.cstatus)
                dispose(wrapper)
                complete!(wrapper.actor)
            else
                push!(wrapper.values, S())
                push!(wrapper.vstatus, falses(wrapper.nsize))
            end
        end

    end
end

function error_received!(wrapper::ZipActorWrapper, err, index::Val{I}) where I
    if !wrapper.cstatus[I]
        dispose(wrapper)
        error!(wrapper.actor, err)
    end
end

function complete_received!(wrapper::ZipActorWrapper, ::Val{I}) where I
    if !all(wrapper.cstatus)
        wrapper.cstatus[I] = true
        if all(wrapper.cstatus) || first_vstatus(wrapper, I) === false
            dispose(wrapper)
            complete!(wrapper.actor)
        end
    end
end

##

##

struct ZipObservable{T, C} <: Subscribable{T}
    sources  :: C
end

function on_subscribe!(observable::ZipObservable{T, C}, actor::A) where { T, C, A }
    S       = typeof(getmstorage(T))
    nsize   = length(observable.sources)
    wrapper = ZipActorWrapper{S, A}(nsize, actor)

    for (index, source) in enumerate(observable.sources)
        wrapper.subscriptions[index] = subscribe!(source, ZipInnerActor{eltype(source), typeof(wrapper), index}(wrapper))
        if all(wrapper.cstatus) || (cstatus(wrapper, index) === true && first_vstatus(wrapper, index) === false)
            dispose(wrapper)
            break
        end
    end

    return ZipSubscription(wrapper)
end

##

struct ZipSubscription{W} <: Teardown
    wrapper :: W
end

as_teardown(::Type{ <: ZipSubscription }) = UnsubscribableTeardownLogic()

function on_unsubscribe!(subscription::ZipSubscription)
    dispose(subscription.wrapper)
    return nothing
end

Base.show(io::IO, ::ZipObservable{D}) where D = print(io, "ZipObservable($D)")
Base.show(io::IO, ::ZipSubscription)          = print(io, "ZipSubscription()")


text:
# A single script to contain most of this package.
# For use within a Pluto notebook
using Plots, SymPy, Roots, ForwardDiff, LinearAlgebra, SpecialFunctions, QuadGK
using Printf, PlutoUI

begin
    Base.adjoint(f::Function) = x -> ForwardDiff.derivative(f, float(x))
    const e = exp(1)
    fisheye(f)=x->atan(f(tan(x)))
    tangent(f,c) = x -> f(c) + f'(c) * (x-c)
    secant(f, a, b) = x -> f(a) + (f(b) - f(a)) / (b-a) * (x - a)

        """

`lim(f, c, n, dir="+")`: means to generate numeric table of values of `f` as `h` gets close to `c`.

Example:
```
f(x) = sin(x) / x
lim(f, 0)
```
    """
    function lim(f::Function, c::Real; n::Int=6, dir="+")
	hs = [(1/10)^i for i in 1:n] # close to 0
	if dir == "+"
	    xs = c .+ hs
	else
	    xs = c .- hs
	end
	ys = map(f, xs)
	[xs ys]
    end

    "`newton(f,x; kwargs)` run Newton's method. Use `with_terminal` to display trace"
    function newton(f, x; kwargs...)
    	p = Roots.ZeroProblem((f,f'),x)
	i = init(p, Roots.Newton(); kwargs...)
    	println(
            @sprintf(
	    	"%s = % 18.16f,\t %s = % 18.16f",
	    	"x0", float(x), "f(x0)", float(f(x)))
	)
	α = x
	for (i,xᵢ) ∈ enumerate(i)
	    α = xᵢ
    	    println(
                @sprintf(
                    "%s = % 18.16f,\t %s = % 18.16f",
                    "x$i", float(xᵢ), "f(x$i)", float(f(xᵢ)))
            )
	end
	α
    end

    """
   sign_chart(f, a, b; atol=1e-4)

Create a sign chart for `f` over `(a,b)`. Returns a tuple with an identified zero or vertical asymptote and the corresponding sign change. The tolerance is used to disambiguate numerically found values.

# Example

```
julia> sign_chart(x -> x/(x-1)^2, -5, 5)
2-element Vector{NamedTuple{(:∞0, :sign_change), Tuple{Float64, String}}}:
 (∞0 = 0.0, sign_change = "- → +")
 (∞0 = 1.0000000000000002, sign_change = "+ → +")
```

"""
function sign_chart(f, a, b; atol=1e-6)
    pm(x) = x < 0 ? "-" : x > 0 ? "+" : "0"
    summarize(f,cp,d) = (DNE_0_∞=cp, sign_change=pm(f(cp-d)) * " → " * pm(f(cp+d)))

    if Roots._is_f_approx_0(f(a),a, eps(), eps()) ||
        Roots._is_f_approx_0(f(b), b, eps(), eps())
        return "Sorry, the endpoints must not be zeros for the function"
    end

    zs = find_zeros(f, a, b)
    pts = vcat(a, zs, b)
    for (u,v) ∈ zip(pts[1:end-1], pts[2:end])
        zs′ = find_zeros(x -> 1/f(x), u, v)
        for z′ ∈ zs′
            flag = false
            for z ∈ zs
                if isapprox(z′, z, atol=atol)
                    flag = true
                    break
                end
            end
            !flag && push!(zs, z′)
        end
    end

    if isempty(zs)
	fc = f(a + (b-a)/2)
	return "No sign change, always " * (fc > 0 ? "positive" : iszero(fc) ? "zero" : "negative")
    end

    sort!(zs)
    m,M = extrema(zs)
    d = min((m-a)/2, (b-M)/2)
    if length(zs) > 1
        d′ = minimum(diff(zs))/2
        d = min(d, d′ )
    end
    summarize.(f, zs, d)
end

    "`riemann(f,a,b,n,[method])` compute Riemann summs"
    function riemann(f::Function, a::Real, b::Real, n::Int; method="right")
	if method == "right"
            meth = (f,l,r) -> f(r) * (r-l)
        elseif method == "left"
            meth= (f,l,r) -> f(l) * (r-l)
        elseif method == "trapezoid"
            meth = (f,l,r) -> (1/2) * (f(l) + f(r)) * (r-l)
        elseif method == "simpsons"
            meth = (f,l,r) -> (1/6) * (f(l) + 4*(f((l+r)/2)) + f(r)) * (r-l)
        end

        xs = range(a, b, length=n+1)
        as = [meth(f, l, r) for (l,r) in zip(xs[1:end-1], xs[2:end])]
        sum(as)
    end

    # Visualization
    "`bisection(f,a,b; n=0)` visualize `n` steps of the bisection method"
    function bisection(f, a, b; n=0)
	p=plot(f, a , b, legend=false)
        a′,b′ = a,b
	show_segment=true
	for i in 1:n
	    c = (a + b)/2
	    if f(c) == 0
		scatter!([c], [f(c)])
		show_segment=false
		break
	    end
	    if f(c) * f(a) < 0
		a,b = a,c
	    else
		a,b = c,b
	    end
	end
	if show_segment
            plot!([a′,b′], [0,0], linewidth=1, color=:gray)
	    plot!([a,b], [0,0], linewidth=5)
	end
	p
    end

    "`trimplot(f,a,b, c= 20;...)` trim large values of `f(x)`, then plot"
    function trimplot(f, a, b, c=20; color=:black, legend=false, kwargs...)
        F = (a,b) -> begin
            fa, fb = f(a), f(b)
            M = max(fa, fb)
            m = min(fa, fb)
            m < -c && return false
            M > c && return false
            true
        end
        xs = range(a, b, length=251)
        cols = find_colors(F, xs, (color, :transparent, :red))
        Plots.plot(xs, f.(xs), colors=cols, legend=legend, kwargs...)
    end

    "`newton_vis` visualize Newton's method"
    function newton_vis(f, x0, a=Inf,b=-Inf; steps=5, kwargs...)
        xs = Float64[x0]
        for i in 1:steps
            push!(xs, xs[end] - f(xs[end]) / f'(xs[end]))
        end

        m,M = extrema(xs)
        m = min(m, a)
        M = max(M, b)

        p = Plots.plot(f, m, M; linewidth=3, legend=false, kwargs...)
        Plots.plot!(p, zero)
        for i in 1:steps
            Plots.plot!(p, [xs[i],xs[i],xs[i+1]], [0,f(xs[i]), 0])
            Plots.scatter!(p, xs[i:i],[0])
        end
        Plots.scatter!(p, [xs[steps+1]], [0])
        p
    end

    # for plotif. This identifies a vector of colors
    function identify_colors(g, xs, colors=(:red, :blue, :black))
        F = (a,b) -> begin
            ga,gb=g(a),g(b)
            ga * gb < 0 && return nothing
            ga >= 0 && return true
            return false
        end
        find_colors(F, xs, colors)
    end

    # F(a,b) returns true, false, or nothing
    function find_colors(F, xs, colors=(:red, :blue, :black))
        n = length(xs)
        cols = repeat([colors[1]], n)
        for i in 1:n-1
            a,b = xs[i], xs[i+1]
            val = F(a,b)
            if val == nothing
                cols[i] = colors[3]
            elseif val
                cols[i] = colors[1]
            else
                cols[i] = colors[2]
            end
        end
        cols[end] = cols[end-1]
        cols
    end

    "`plotif(f, g, a, b; kwargs...)` Plot `f`, using red to indicate where `g(x)>0`."
    function plotif(f, g, a, b; kwargs...)
    	xs = range(a, b, length=251)
	cols = identify_colors(g, xs)
    	Plots.plot(xs, f; color=cols, legend=false, kwargs...)
    end

    nothing
end


text:
#####
##### Weighted Essentially Non-Oscillatory (WENO) scheme of order 5
#####

struct WENO5 <: AbstractUpwindBiasedAdvectionScheme end

@inline boundary_buffer(::WENO5) = 2

@inline symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, ::WENO5, c) = symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, centered_fourth_order, c)
@inline symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, ::WENO5, c) = symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, centered_fourth_order, c)
@inline symmetric_interpolate_zᵃᵃᶠ(i, j, k, grid, ::WENO5, c) = symmetric_interpolate_zᵃᵃᶠ(i, j, k, grid, centered_fourth_order, c)

@inline symmetric_interpolate_xᶜᵃᵃ(i, j, k, grid, ::WENO5, u) = symmetric_interpolate_xᶜᵃᵃ(i, j, k, grid, centered_fourth_order, u)
@inline symmetric_interpolate_yᵃᶜᵃ(i, j, k, grid, ::WENO5, v) = symmetric_interpolate_yᵃᶜᵃ(i, j, k, grid, centered_fourth_order, v)
@inline symmetric_interpolate_zᵃᵃᶜ(i, j, k, grid, ::WENO5, w) = symmetric_interpolate_zᵃᵃᶜ(i, j, k, grid, centered_fourth_order, w)

#####
##### ENO interpolants of size 3
#####

@inline left_biased_px₀(i, j, k, ψ) = @inbounds + 1/3 * ψ[i-1, j, k] + 5/6 * ψ[i,   j, k] -  1/6 * ψ[i+1, j, k]
@inline left_biased_px₁(i, j, k, ψ) = @inbounds - 1/6 * ψ[i-2, j, k] + 5/6 * ψ[i-1, j, k] +  1/3 * ψ[i,   j, k]
@inline left_biased_px₂(i, j, k, ψ) = @inbounds + 1/3 * ψ[i-3, j, k] - 7/6 * ψ[i-2, j, k] + 11/6 * ψ[i-1, j, k]

@inline left_biased_py₀(i, j, k, ψ) = @inbounds + 1/3 * ψ[i, j-1, k] + 5/6 * ψ[i, j,   k] -  1/6 * ψ[i, j+1, k]
@inline left_biased_py₁(i, j, k, ψ) = @inbounds - 1/6 * ψ[i, j-2, k] + 5/6 * ψ[i, j-1, k] +  1/3 * ψ[i, j  , k]
@inline left_biased_py₂(i, j, k, ψ) = @inbounds + 1/3 * ψ[i, j-3, k] - 7/6 * ψ[i, j-2, k] + 11/6 * ψ[i, j-1, k]

@inline left_biased_pz₀(i, j, k, ψ) = @inbounds + 1/3 * ψ[i, j, k-1] + 5/6 * ψ[i, j,   k] -  1/6 * ψ[i, j, k+1]
@inline left_biased_pz₁(i, j, k, ψ) = @inbounds - 1/6 * ψ[i, j, k-2] + 5/6 * ψ[i, j, k-1] +  1/3 * ψ[i, j,   k]
@inline left_biased_pz₂(i, j, k, ψ) = @inbounds + 1/3 * ψ[i, j, k-3] - 7/6 * ψ[i, j, k-2] + 11/6 * ψ[i, j, k-1]

@inline right_biased_px₀(i, j, k, ψ) = @inbounds + 11/6 * ψ[i,   j, k] - 7/6 * ψ[i+1, j, k] + 1/3 * ψ[i+2, j, k]
@inline right_biased_px₁(i, j, k, ψ) = @inbounds +  1/3 * ψ[i-1, j, k] + 5/6 * ψ[i,   j, k] - 1/6 * ψ[i+1, j, k]
@inline right_biased_px₂(i, j, k, ψ) = @inbounds -  1/6 * ψ[i-2, j, k] + 5/6 * ψ[i-1, j, k] + 1/3 * ψ[i,   j, k]

@inline right_biased_py₀(i, j, k, ψ) = @inbounds + 11/6 * ψ[i,   j, k] - 7/6 * ψ[i, j+1, k] + 1/3 * ψ[i, j+2, k]
@inline right_biased_py₁(i, j, k, ψ) = @inbounds +  1/3 * ψ[i, j-1, k] + 5/6 * ψ[i,   j, k] - 1/6 * ψ[i, j+1, k]
@inline right_biased_py₂(i, j, k, ψ) = @inbounds -  1/6 * ψ[i, j-2, k] + 5/6 * ψ[i, j-1, k] + 1/3 * ψ[i,   j, k]

@inline right_biased_pz₀(i, j, k, ψ) = @inbounds + 11/6 * ψ[i, j,   k] - 7/6 * ψ[i, j, k+1] + 1/3 * ψ[i, j, k+2]
@inline right_biased_pz₁(i, j, k, ψ) = @inbounds +  1/3 * ψ[i, j, k-1] + 5/6 * ψ[i, j,   k] - 1/6 * ψ[i, j, k+1]
@inline right_biased_pz₂(i, j, k, ψ) = @inbounds -  1/6 * ψ[i, j, k-2] + 5/6 * ψ[i, j, k-1] + 1/3 * ψ[i, j,   k]

#####
##### Jiang & Shu (1996) WENO smoothness indicators. See also Equation 2.63 in Shu (1998).
#####

@inline left_biased_βx₀(i, j, k, ψ) = @inbounds 13/12 * (ψ[i-1, j, k] - 2ψ[i,   j, k] + ψ[i+1, j, k])^2 + 1/4 * (3ψ[i-1, j, k] - 4ψ[i,   j, k] +  ψ[i+1, j, k])^2
@inline left_biased_βx₁(i, j, k, ψ) = @inbounds 13/12 * (ψ[i-2, j, k] - 2ψ[i-1, j, k] + ψ[i,   j, k])^2 + 1/4 * ( ψ[i-2, j, k]                 -  ψ[i,   j, k])^2
@inline left_biased_βx₂(i, j, k, ψ) = @inbounds 13/12 * (ψ[i-3, j, k] - 2ψ[i-2, j, k] + ψ[i-1, j, k])^2 + 1/4 * ( ψ[i-3, j, k] - 4ψ[i-2, j, k] + 3ψ[i-1, j, k])^2

@inline left_biased_βy₀(i, j, k, ψ) = @inbounds 13/12 * (ψ[i, j-1, k] - 2ψ[i, j,   k] + ψ[i, j+1, k])^2 + 1/4 * (3ψ[i, j-1, k] - 4ψ[i,   j, k] +  ψ[i, j+1, k])^2
@inline left_biased_βy₁(i, j, k, ψ) = @inbounds 13/12 * (ψ[i, j-2, k] - 2ψ[i, j-1, k] + ψ[i, j,   k])^2 + 1/4 * ( ψ[i, j-2, k]                 -  ψ[i,   j, k])^2
@inline left_biased_βy₂(i, j, k, ψ) = @inbounds 13/12 * (ψ[i, j-3, k] - 2ψ[i, j-2, k] + ψ[i, j-1, k])^2 + 1/4 * ( ψ[i, j-3, k] - 4ψ[i, j-2, k] + 3ψ[i, j-1, k])^2

@inline left_biased_βz₀(i, j, k, ψ) = @inbounds 13/12 * (ψ[i, j, k-1] - 2ψ[i, j,   k] + ψ[i, j, k+1])^2 + 1/4 * (3ψ[i, j, k-1] - 4ψ[i, j,   k] +  ψ[i, j, k+1])^2
@inline left_biased_βz₁(i, j, k, ψ) = @inbounds 13/12 * (ψ[i, j, k-2] - 2ψ[i, j, k-1] + ψ[i, j,   k])^2 + 1/4 * ( ψ[i, j, k-2]                 -  ψ[i, j,   k])^2
@inline left_biased_βz₂(i, j, k, ψ) = @inbounds 13/12 * (ψ[i, j, k-3] - 2ψ[i, j, k-2] + ψ[i, j, k-1])^2 + 1/4 * ( ψ[i, j, k-3] - 4ψ[i, j, k-2] + 3ψ[i, j, k-1])^2

# Right-biased smoothness indicators are a reflection or "symmetric modification" of the left-biased smoothness
# indicators around grid point `i-1/2`.

@inline right_biased_βx₀(i, j, k, ψ) = @inbounds 13/12 * (ψ[i,   j, k] - 2ψ[i+1, j, k] + ψ[i+2, j, k])^2 + 1/4 * ( ψ[i,   j, k] - 4ψ[i+1, j, k] + 3ψ[i+2, j, k])^2
@inline right_biased_βx₁(i, j, k, ψ) = @inbounds 13/12 * (ψ[i-1, j, k] - 2ψ[i,   j, k] + ψ[i+1, j, k])^2 + 1/4 * ( ψ[i-1, j, k]                 -  ψ[i+1, j, k])^2
@inline right_biased_βx₂(i, j, k, ψ) = @inbounds 13/12 * (ψ[i-2, j, k] - 2ψ[i-1, j, k] + ψ[i,   j, k])^2 + 1/4 * (3ψ[i-2, j, k] - 4ψ[i-1, j, k] +  ψ[i,   j, k])^2

@inline right_biased_βy₀(i, j, k, ψ) = @inbounds 13/12 * (ψ[i, j,   k] - 2ψ[i, j+1, k] + ψ[i, j+2, k])^2 + 1/4 * ( ψ[i,   j, k] - 4ψ[i, j+1, k] + 3ψ[i, j+2, k])^2
@inline right_biased_βy₁(i, j, k, ψ) = @inbounds 13/12 * (ψ[i, j-1, k] - 2ψ[i, j,   k] + ψ[i, j+1, k])^2 + 1/4 * ( ψ[i, j-1, k]                 -  ψ[i, j+1, k])^2
@inline right_biased_βy₂(i, j, k, ψ) = @inbounds 13/12 * (ψ[i, j-2, k] - 2ψ[i, j-1, k] + ψ[i, j,   k])^2 + 1/4 * (3ψ[i, j-2, k] - 4ψ[i, j-1, k] +  ψ[i,   j, k])^2

@inline right_biased_βz₀(i, j, k, ψ) = @inbounds 13/12 * (ψ[i, j,   k] - 2ψ[i, j, k+1] + ψ[i, j, k+2])^2 + 1/4 * ( ψ[i, j,   k] - 4ψ[i, j, k+1] + 3ψ[i, j, k+2])^2
@inline right_biased_βz₁(i, j, k, ψ) = @inbounds 13/12 * (ψ[i, j, k-1] - 2ψ[i, j,   k] + ψ[i, j, k+1])^2 + 1/4 * ( ψ[i, j, k-1]                 -  ψ[i, j, k+1])^2
@inline right_biased_βz₂(i, j, k, ψ) = @inbounds 13/12 * (ψ[i, j, k-2] - 2ψ[i, j, k-1] + ψ[i, j,   k])^2 + 1/4 * (3ψ[i, j, k-2] - 4ψ[i, j, k-1] +  ψ[i, j,   k])^2

#####
##### WENO-5 optimal weights
#####

const C3₀ = 3/10
const C3₁ = 3/5
const C3₂ = 1/10

#####
##### WENO-5 raw weights
#####

# Note: these constants may need to be changed for smooth solutions and/or fine grid.
const ε = 1e-6
const ƞ = 2  # WENO exponent

@inline left_biased_αx₀(i, j, k, ψ) = C3₀ / (left_biased_βx₀(i, j, k, ψ) + ε)^ƞ
@inline left_biased_αx₁(i, j, k, ψ) = C3₁ / (left_biased_βx₁(i, j, k, ψ) + ε)^ƞ
@inline left_biased_αx₂(i, j, k, ψ) = C3₂ / (left_biased_βx₂(i, j, k, ψ) + ε)^ƞ

@inline left_biased_αy₀(i, j, k, ψ) = C3₀ / (left_biased_βy₀(i, j, k, ψ) + ε)^ƞ
@inline left_biased_αy₁(i, j, k, ψ) = C3₁ / (left_biased_βy₁(i, j, k, ψ) + ε)^ƞ
@inline left_biased_αy₂(i, j, k, ψ) = C3₂ / (left_biased_βy₂(i, j, k, ψ) + ε)^ƞ

@inline left_biased_αz₀(i, j, k, ψ) = C3₀ / (left_biased_βz₀(i, j, k, ψ) + ε)^ƞ
@inline left_biased_αz₁(i, j, k, ψ) = C3₁ / (left_biased_βz₁(i, j, k, ψ) + ε)^ƞ
@inline left_biased_αz₂(i, j, k, ψ) = C3₂ / (left_biased_βz₂(i, j, k, ψ) + ε)^ƞ

@inline right_biased_αx₀(i, j, k, ψ) = C3₂ / (right_biased_βx₀(i, j, k, ψ) + ε)^ƞ
@inline right_biased_αx₁(i, j, k, ψ) = C3₁ / (right_biased_βx₁(i, j, k, ψ) + ε)^ƞ
@inline right_biased_αx₂(i, j, k, ψ) = C3₀ / (right_biased_βx₂(i, j, k, ψ) + ε)^ƞ

@inline right_biased_αy₀(i, j, k, ψ) = C3₂ / (right_biased_βy₀(i, j, k, ψ) + ε)^ƞ
@inline right_biased_αy₁(i, j, k, ψ) = C3₁ / (right_biased_βy₁(i, j, k, ψ) + ε)^ƞ
@inline right_biased_αy₂(i, j, k, ψ) = C3₀ / (right_biased_βy₂(i, j, k, ψ) + ε)^ƞ

@inline right_biased_αz₀(i, j, k, ψ) = C3₂ / (right_biased_βz₀(i, j, k, ψ) + ε)^ƞ
@inline right_biased_αz₁(i, j, k, ψ) = C3₁ / (right_biased_βz₁(i, j, k, ψ) + ε)^ƞ
@inline right_biased_αz₂(i, j, k, ψ) = C3₀ / (right_biased_βz₂(i, j, k, ψ) + ε)^ƞ

#####
##### WENO-5 normalized weights
#####

@inline function left_biased_weno5_weights_x(i, j, k, ψ)
    α₀ = left_biased_αx₀(i, j, k, ψ)
    α₁ = left_biased_αx₁(i, j, k, ψ)
    α₂ = left_biased_αx₂(i, j, k, ψ)

    Σα = α₀ + α₁ + α₂
    w₀ = α₀ / Σα
    w₁ = α₁ / Σα
    w₂ = α₂ / Σα

    return w₀, w₁, w₂
end

@inline function left_biased_weno5_weights_y(i, j, k, ψ)
    α₀ = left_biased_αy₀(i, j, k, ψ)
    α₁ = left_biased_αy₁(i, j, k, ψ)
    α₂ = left_biased_αy₂(i, j, k, ψ)

    Σα = α₀ + α₁ + α₂
    w₀ = α₀ / Σα
    w₁ = α₁ / Σα
    w₂ = α₂ / Σα

    return w₀, w₁, w₂
end

@inline function left_biased_weno5_weights_z(i, j, k, ψ)
    α₀ = left_biased_αz₀(i, j, k, ψ)
    α₁ = left_biased_αz₁(i, j, k, ψ)
    α₂ = left_biased_αz₂(i, j, k, ψ)

    Σα = α₀ + α₁ + α₂
    w₀ = α₀ / Σα
    w₁ = α₁ / Σα
    w₂ = α₂ / Σα

    return w₀, w₁, w₂
end

@inline function right_biased_weno5_weights_x(i, j, k, ψ)
    α₀ = right_biased_αx₀(i, j, k, ψ)
    α₁ = right_biased_αx₁(i, j, k, ψ)
    α₂ = right_biased_αx₂(i, j, k, ψ)

    Σα = α₀ + α₁ + α₂
    w₀ = α₀ / Σα
    w₁ = α₁ / Σα
    w₂ = α₂ / Σα

    return w₀, w₁, w₂
end

@inline function right_biased_weno5_weights_y(i, j, k, ψ)
    α₀ = right_biased_αy₀(i, j, k, ψ)
    α₁ = right_biased_αy₁(i, j, k, ψ)
    α₂ = right_biased_αy₂(i, j, k, ψ)

    Σα = α₀ + α₁ + α₂
    w₀ = α₀ / Σα
    w₁ = α₁ / Σα
    w₂ = α₂ / Σα

    return w₀, w₁, w₂
end

@inline function right_biased_weno5_weights_z(i, j, k, ψ)
    α₀ = right_biased_αz₀(i, j, k, ψ)
    α₁ = right_biased_αz₁(i, j, k, ψ)
    α₂ = right_biased_αz₂(i, j, k, ψ)

    Σα = α₀ + α₁ + α₂
    w₀ = α₀ / Σα
    w₁ = α₁ / Σα
    w₂ = α₂ / Σα

    return w₀, w₁, w₂
end

#####
##### WENO-5 reconstruction
#####

@inline function left_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, ::WENO5, ψ)
    w₀, w₁, w₂ = left_biased_weno5_weights_x(i, j, k, ψ)
    return w₀ * left_biased_px₀(i, j, k, ψ) + w₁ * left_biased_px₁(i, j, k, ψ) + w₂ * left_biased_px₂(i, j, k, ψ)
end

@inline function left_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, ::WENO5, ψ)
    w₀, w₁, w₂ = left_biased_weno5_weights_y(i, j, k, ψ)
    return w₀ * left_biased_py₀(i, j, k, ψ) + w₁ * left_biased_py₁(i, j, k, ψ) + w₂ * left_biased_py₂(i, j, k, ψ)
end

@inline function left_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, ::WENO5, ψ)
    w₀, w₁, w₂ = left_biased_weno5_weights_z(i, j, k, ψ)
    return w₀ * left_biased_pz₀(i, j, k, ψ) + w₁ * left_biased_pz₁(i, j, k, ψ) + w₂ * left_biased_pz₂(i, j, k, ψ)
end

@inline function right_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, ::WENO5, ψ)
    w₀, w₁, w₂ = right_biased_weno5_weights_x(i, j, k, ψ)
    return w₀ * right_biased_px₀(i, j, k, ψ) + w₁ * right_biased_px₁(i, j, k, ψ) + w₂ * right_biased_px₂(i, j, k, ψ)
end

@inline function right_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, ::WENO5, ψ)
    w₀, w₁, w₂ = right_biased_weno5_weights_y(i, j, k, ψ)
    return w₀ * right_biased_py₀(i, j, k, ψ) + w₁ * right_biased_py₁(i, j, k, ψ) + w₂ * right_biased_py₂(i, j, k, ψ)
end

@inline function right_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, ::WENO5, ψ)
    w₀, w₁, w₂ = right_biased_weno5_weights_z(i, j, k, ψ)
    return w₀ * right_biased_pz₀(i, j, k, ψ) + w₁ * right_biased_pz₁(i, j, k, ψ) + w₂ * right_biased_pz₂(i, j, k, ψ)
end

@inline left_biased_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme::WENO5, ψ) = left_biased_interpolate_xᶠᵃᵃ(i+1, j, k, grid, scheme, ψ)
@inline left_biased_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme::WENO5, ψ) = left_biased_interpolate_yᵃᶠᵃ(i, j+1, k, grid, scheme, ψ)
@inline left_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme::WENO5, ψ) = left_biased_interpolate_zᵃᵃᶠ(i, j, k+1, grid, scheme, ψ)

@inline right_biased_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme::WENO5, ψ) = right_biased_interpolate_xᶠᵃᵃ(i+1, j, k, grid, scheme, ψ)
@inline right_biased_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme::WENO5, ψ) = right_biased_interpolate_yᵃᶠᵃ(i, j+1, k, grid, scheme, ψ)
@inline right_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme::WENO5, ψ) = right_biased_interpolate_zᵃᵃᶠ(i, j, k+1, grid, scheme, ψ)


text:
module MultiFieldFESpacesTests

using BlockArrays
using FillArrays
using Gridap.Arrays
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.Fields
using Gridap.Integration
using Gridap.ReferenceFEs
using Gridap.CellData
using Test

using Gridap.MultiField
using Gridap.Arrays: BlockArrayCooMap

order = 2

domain = (0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)

trian = get_triangulation(model)
degree = order
quad = CellQuadrature(trian,degree)

V = TestFESpace(model,ReferenceFE(lagrangian,Float64,order);conformity=:H1)
Q = TestFESpace(model,ReferenceFE(lagrangian,Float64,order-1),conformity=:L2)

U = TrialFESpace(V)
P = TrialFESpace(Q)

multi_field_style = ConsecutiveMultiFieldStyle()

Y = MultiFieldFESpace(Vector{Float64},[V,Q],multi_field_style)
X = MultiFieldFESpace(Vector{Float64},[U,P],multi_field_style)

@test num_free_dofs(X) == num_free_dofs(U) + num_free_dofs(P)
@test num_free_dofs(X) == num_free_dofs(Y)

dy = get_cell_shapefuns(Y)
dv, dq = dy

dx = get_cell_shapefuns_trial(X)
du, dp = dx

cellmat = integrate(dv*du,quad)
cellvec = integrate(dv*2,quad)
cellids = get_cell_to_bgcell(trian)
cellmatvec = pair_arrays(cellmat,cellvec)
@test isa(cellmat, LazyArray{<:Fill{<:BlockArrayCooMap}})
@test is_nonzero_block(cellmat[1],1,1)
@test is_zero_block(cellmat[1],1,2)
@test isa(cellvec, LazyArray{<:Fill{<:BlockArrayCooMap}})
@test is_nonzero_block(cellvec[1],1)
@test is_zero_block(cellvec[1],2)

matvecdata = (cellmatvec,cellids,cellids)
matdata = (cellmat,cellids,cellids)
vecdata = (cellvec,cellids)

free_values = rand(num_free_dofs(X))
xh = FEFunction(X,free_values)
test_fe_function(xh)
@test isa(xh,FEFunction)
uh, ph = xh
@test isa(uh,FEFunction)
@test isa(ph,FEFunction)

cell_isconstr = get_cell_isconstrained(X)
@test cell_isconstr == Fill(false,num_cells(model))

cell_constr = get_cell_constraints(X)
@test isa(cell_constr,LazyArray{<:Fill{<:BlockArrayCooMap}})

cell_dof_ids = get_cell_dof_ids(X)
@test isa(cell_dof_ids,LazyArray{<:Fill{<:BlockArrayCooMap}})

cf = CellField(X,get_cell_dof_ids(X))
@test isa(cf,MultiFieldCellField)

test_fe_space(X,matvecdata,matdata,vecdata)
test_fe_space(Y,matvecdata,matdata,vecdata)

#using Gridap.Visualization
#writevtk(trian,"trian";nsubcells=30,cellfields=["uh" => uh, "ph"=> ph])

f(x) = sin(4*pi*(x[1]-x[2]^2))+1
fh = interpolate([f,f],X)
fh = interpolate_everywhere([f,f],X)
fh = interpolate_dirichlet([f,f],X)

end # module


text:
# TODO integration underway with Manifolds.jl, see RoME #244, also see IIF #467 regarding consolidation effort.

export projectCartesian

"""
$(TYPEDEF)

XY Euclidean manifold variable node softtype.
"""
@defVariable Point2 2 (:Euclid, :Euclid)

"""
$(TYPEDEF)

XYZ Euclidean manifold variable node softtype.

Example
-------
```julia
p3 = Point3()
```
"""
@defVariable Point3 3 (:Euclid,:Euclid,:Euclid)

"""
$(TYPEDEF)

Pose2 is a SE(2) mechanization of two Euclidean translations and one Circular rotation, used for general 2D SLAM.
"""
@defVariable Pose2 3 (:Euclid,:Euclid,:Circular)

"""
$(TYPEDEF)

Pose3 is currently a Euler angle mechanization of three Euclidean translations and three Circular rotation.

Future:
------
- Work in progress on AMP3D for proper non-Euler angle on-manifold operations.
- TODO the AMP upgrade is aimed at resolving 3D to Quat/SE3/SP3 -- current Euler angles will be replaced
"""
@defVariable Pose3 6 (:Euclid,:Euclid,:Euclid,:Circular,:Circular,:Circular)

"""
$(TYPEDEF)

Dynamic point in 2D space with velocity components: `x, y, dx/dt, dy/dt`

"""
@defVariable DynPoint2 4 (:Euclid,:Euclid,:Euclid,:Euclid)

"""
$(TYPEDEF)

Dynamic pose variable with velocity components: `x, y, theta, dx/dt, dy/dt`
"""
@defVariable DynPose2 5 (:Euclid,:Euclid,:Circular,:Euclid,:Euclid)

"""
$SIGNATURES

Function to project only XY data onto Cartesian plane for 2D plotting.
"""
projectCartesian(pose::Union{Point2,Point3, Pose2, Pose3, DynPoint2, DynPose2},
                 x::Vector{Float64}) = [x[1]; x[2]; 0]

#


text:
module Examples

# Utilities
include("reparametrization.jl")

# Dynamical systems
include("pitchfork.jl")
include("transcritical.jl")
include("bautin.jl")
include("cusp.jl")
include("calcium.jl")
include("predator_prey.jl")
include("bazykin_85.jl")
include("bazykin_khibnik_81.jl")
include("duffing_van_der_pol.jl")
include("morris_lecar.jl")

example_modules() = [
    Pitchfork,
    Transcritical,
    Bautin,
    Cusp,
    Calcium,
    PredatorPrey,
    Bazykin85,
    BazykinKhibnik81,
    DuffingVanDerPol,
    MorrisLecar,
]

examples() = [nameof(ex) => ex for ex in example_modules()]

end  # module

using .Examples: examples
using .Examples.Reparametrization: reparametrize


text:
# ------------------------------------------------------------
# 1D change of basis matrix

function Freq2Wave(samples::StridedVector, wavename::AbstractString, J::Int, B::Float64=NaN; uniform::Bool=true, args...)
	vm = van_moment(wavename)
	( Nint = 2^J ) >= 2*vm-1 || throw(AssertionError("Scale it not large enough for this wavelet"))

	M = length(samples)
	if Nint >= M
		warn("The scale is high compared to the number of samples")
	end

	# Weights for non-uniform samples
	uni_samples = isuniform(samples)
	if uni_samples
		sample_dist = samples[2] - samples[1]
		W = Nullable{ Vector{Complex{Float64}} }()
	else
		isnan(B) && error("Samples are not uniform; supply bandwidth")
		Nint <= 2*B || warn("The scale is high compared to the bandwidth")
		W = sqrt(weights(samples, B))
		W = Nullable(complex( W ))
	end

	# Fourier transform of the internal scaling function
	internal = FourScalingFunc( samples, wavename, J; args... )

	# The number of internal wavelets in the reconstruction
	if hasboundary(wavename)
		Nint -= 2*vm
		Nint <= 0 && error("Too few wavelets: Boundary functions overlap")
	end

    if uni_samples && isinteger(1/sample_dist) && uniform
	# For suitable uniform samples, the ordinary FFT can be used instead of NFFT
		p = FFTPlan(samples, J, Nint)
	else
		# NFFTPlans: Frequencies must be in the torus [-1/2, 1/2)
		xi = samples*2.0^(-J)
		frac!(xi)
		p = NFFTPlan(xi, Nint)
	end

	# Wavelets w/o boundary
	if !hasboundary(wavename)
		return Freq2NoBoundaryWave1D(internal, W, J, wavename, p)
	else
		left = FourScalingFunc( samples, wavename, 'L', J; args... )
		right = FourScalingFunc( samples, wavename, 'R', J; args... )

		return Freq2BoundaryWave1D(internal, W, J, wavename, p, left, right)
	end
end

@doc """
	isuniform(Freq2Wave)

Is the change of basis matrix based on uniform samples.
"""->
function isuniform(T::Freq2Wave)
	isnull(T.weights)
end

@doc """
	hasboundary(T::Freq2Wave) -> Bool

Does the wavelet in the change of basis matrix `T` have boundary correction.
"""->
function hasboundary(T::Freq2Wave)
	isdefined(T, :left)
end

# ------------------------------------------------------------
# Basic operations for 1D Freq2Wave

function Base.collect(T::Freq2NoBoundaryWave1D)
	M, N = size(T)

	F = Array{Complex{Float64}}(M, N)
	offset = div(N, 2) + 1
	for n in 1:N
		for m in 1:M
			@inbounds F[m,n] = T.internal[m]*cis( -twoπ*(n-offset)*T.NFFT.x[m] )
		end
	end

	if !isuniform(T)
		scale!(get(T.weights), F)
	end

	return F
end

function Base.collect(T::Freq2BoundaryWave1D)
	M, N = size(T)
	F = Array{Complex{Float64}}(M, N)

	# Left boundary
	p = van_moment(T)
	F[:,1:p] = T.left

	# Internal function
	offset = div(N, 2) + 1
	for n in (p+1):(N-p)
		for m in 1:M
			@inbounds F[m,n] = T.internal[m]*cis( -twoπ*(n-offset)*T.NFFT.x[m] )
		end
	end

	# Right boundary
	F[:,N-p+1:N] = T.right

	if !isuniform(T)
		scale!(get(T.weights), F)
	end

	return F
end

# ------------------------------------------------------------
# 2D change of basis matrix

@doc """
	dim(Freq2wave)

Return the dimension of the `T`.
"""->
dim(::Freq2Wave1D) = 1
dim(::Freq2Wave2D) = 2

@doc """
	wscale(Freq2Wave)

Return the scale of the wavelet coefficients.
"""->
function wscale(T::Freq2Wave)
	T.J
end

@doc """
	wsize(T::Freq2Wave)

The number of reconstructed wavelet coefficients (in each dimension).

- When `D` == 1, the output is (Int,)
- When `D` == 2, the output is (Int,Int)
"""->
function wsize(T::Freq2BoundaryWave2D)
	N = 2^wscale(T)
	return (N, N)
end

wsize(T::Freq2BoundaryWave1D) = (2^wscale(T),)
wsize(T::Freq2NoBoundaryWave2D) = T.NFFT.N
wsize(T::Freq2NoBoundaryWave1D) = T.NFFT.N

function UnifFourScalingFunc(samples::StridedMatrix{Float64}, wavename::AbstractString, J::Integer; args...)
	# Test if samples are on a grid and in the correct order
	usamplesx = unique(view(samples,:,1))
	usamplesy = unique(view(samples,:,2))
	Mx = length(usamplesx)
	My = length(usamplesy)
	grid_samples = grid( (Mx, My), usamplesx[2]-usamplesx[1] )

	samples == grid_samples ||
	throw(AssertionError("Samples are not on a grid or not ordered properly"))

	# Internal scaling function
	internalx = FourScalingFunc( usamplesx, wavename, J )
	internalx = kron(internalx, ones(My))
	internaly = FourScalingFunc( usamplesy, wavename, J )
	internaly = repmat(internaly, Mx)
	internal = Dict( :x => internalx, :y => internaly )

	# Boundary scaling functions
	if !hasboundary(wavename)
		return internal
	else
		vm = van_moment(wavename)

		leftx = FourScalingFunc( usamplesx, wavename, 'L', J; args... )
		lefty = FourScalingFunc( usamplesy, wavename, 'L', J; args... )
		left = Dict( :x => repmat(leftx, My), :y => kron(lefty, ones(Mx)) )

		rightx = FourScalingFunc( usamplesx, wavename, 'R', J; args... )
		righty = FourScalingFunc( usamplesy, wavename, 'R', J; args... )
		right = Dict( :x => repmat(rightx, My), :y => kron(righty, ones(Mx)) )

		return internal, left, right
	end
end

function NotUnifFourScalingFunc(samples::StridedMatrix{Float64}, wavename::AbstractString, J::Integer, B::Float64; args...)
	# Fourier transform of the internal scaling function
	int = FourScalingFunc( samples, wavename, J )
	internal = Dict( :x => int[:,1], :y => int[:,2] )

	# Boundary scaling functions
	if !hasboundary(wavename)
		return internal
	else
		samplesx = view(samples, :, 1)
		samplesy = view(samples, :, 2)

		left = Dict{ Symbol, Matrix{Complex{Float64}} }()
		left[:x] = FourScalingFunc( samplesx, wavename, 'L', J; args... )
		left[:y] = FourScalingFunc( samplesy, wavename, 'L', J; args... )

		right = Dict{ Symbol, Matrix{Complex{Float64}} }()
		right[:x] = FourScalingFunc( samplesx, wavename, 'R', J; args... )
		right[:y] = FourScalingFunc( samplesy, wavename, 'R', J; args... )

		return internal, left, right
	end
end

function Freq2Wave(samples::StridedMatrix{Float64}, wavename::AbstractString, J::Int, B::Float64=NaN; uniform::Bool=true, args...)
	vm = van_moment(wavename)
	( Nint = 2^J ) >= 2*vm-1 || throw(AssertionError("Scale it not large enough for this wavelet"))
	M = size(samples, 1)
	size(samples,2) == 2 || throw(DimensionMismatch("Samples must have two columns"))

	if Nint >= M
		warn("The scale is high compared to the number of samples")
	end

	uni_samples = isuniform(samples)
	if uni_samples
		W = Nullable{ Vector{Complex{Float64}} }()

		sample_dist = samples[2,2] - samples[1,2]
		scaling_funcs = UnifFourScalingFunc(samples, wavename, J; args...)
	else
		# Weights for non-uniform samples
		isnan(B) && error("Samples are not uniform; supply bandwidth")
		Nint <= 2*B || warn("The scale is high compared to the bandwidth")
		W = sqrt(weights(samples, B))
		W = Nullable(complex( W ))

		scaling_funcs = NotUnifFourScalingFunc(samples, wavename, J, B; args...)
	end

	if hasboundary(wavename)
		Nint -= 2*vm
		Nint <= 0 && error("Too few wavelets: Boundary functions overlap")

		internal = scaling_funcs[1]
		left = scaling_funcs[2]
		right = scaling_funcs[3]
	else
		internal = scaling_funcs
	end

    xi = samples'
	if uni_samples && isinteger(1/sample_dist) && uniform
        # For suitable uniform samples, the ordinary FFT can be used instead of NFFT
        p = FFTPlan(xi, J, (Nint,Nint))
	else
        # NFFTPlans: Frequencies must be in the torus [-1/2, 1/2)^2
        scale!(xi, 2.0^(-J))
        frac!(xi)
        p = NFFTPlan(xi, (Nint,Nint))
    end

	# Wavelets w/o boundary
	if hasboundary(wavename)
		return Freq2BoundaryWave2D(internal, W, J, wavename, p, left, right)
	else
		return Freq2NoBoundaryWave2D(internal, W, J, wavename, p)
	end
end

function Base.A_mul_B!(y::DenseVector{Complex{Float64}}, T::Freq2NoBoundaryWave1D, x::DenseVector{Complex{Float64}})
	size(T,1) == length(y) || throw(DimensionMismatch())
	size(T,2) == length(x) || throw(DimensionMismatch())

	myfft!(y, T.NFFT, x)
	had!(y, T.internal)

	isuniform(T) || had!(y, get(T.weights))

	return y
end

function Base.A_mul_B!(y::DenseVector{Complex{Float64}}, T::Freq2NoBoundaryWave2D, X::DenseMatrix{Complex{Float64}})
	(M = size(T,1)) == length(y) || throw(DimensionMismatch())
	wsize(T) == size(X) || throw(DimensionMismatch())

	myfft!(y, T.NFFT, X)
	for m in 1:M
		@inbounds y[m] *= T.internal[:x][m] * T.internal[:y][m]
	end

	isuniform(T) || had!(y, get(T.weights))

	return y
end

function Base.A_mul_B!(y::DenseVector{Complex{Float64}}, T::Freq2BoundaryWave1D, x::DenseVector{Complex{Float64}})
	(M = size(T,1)) == length(y) || throw(DimensionMismatch())
	size(T,2) == length(x) || throw(DimensionMismatch())

	xleft, xint, xright = split(x, van_moment(T))

	# Internal scaling function
	myfft!(y, T.NFFT, xint)
	had!(y, T.internal)

	# Contribution from the boundaries
	BLAS.gemv!('N', ComplexOne, T.left, xleft, ComplexOne, y)
	BLAS.gemv!('N', ComplexOne, T.right, xright, ComplexOne, y)

	isuniform(T) || had!(y, get(T.weights))

	return y
end

function Base.A_mul_B!(y::DenseVector{Complex{Float64}}, T::Freq2BoundaryWave2D, X::DenseMatrix{Complex{Float64}})
	(M = size(T,1)) == length(y) || throw(DimensionMismatch())
	(N = wsize(T)) == size(X) || throw(DimensionMismatch())

	# Internal scaling functions
	vm = van_moment(T)
	S = split(X, vm)
	myfft!(y, T.NFFT, S.II)
	for m in 1:M
		@inbounds y[m] *= T.internal[:x][m] * T.internal[:y][m]
	end

	onesp = ones(Complex{Float64}, vm)

	# ------------------------------------------------------------
	# Corners

	# LL
	A_mul_B!( T.tmpMulVec, T.left[:x], S.LL )
	had!( T.tmpMulVec, T.left[:y] )
	# y += sum(T.tmpMulVec,2) :
	BLAS.gemv!('N', ComplexOne, T.tmpMulVec, onesp, ComplexOne, y)

	# RL
	A_mul_B!( T.tmpMulVec, T.right[:x], S.RL )
	had!( T.tmpMulVec, T.left[:y] )
	BLAS.gemv!('N', ComplexOne, T.tmpMulVec, onesp, ComplexOne, y)

	# RR
	A_mul_B!( T.tmpMulVec, T.right[:x], S.RR )
	had!( T.tmpMulVec, T.right[:y] )
	BLAS.gemv!('N', ComplexOne, T.tmpMulVec, onesp, ComplexOne, y)

	# LR
	A_mul_B!( T.tmpMulVec, T.left[:x], S.LR )
	had!( T.tmpMulVec, T.right[:y] )
	BLAS.gemv!('N', ComplexOne, T.tmpMulVec, onesp, ComplexOne, y)

	# ------------------------------------------------------------
	# Sides

	# IL
	NFFT.nfft!(T.NFFTx, S.IL, T.tmpMulVec)
	scale!(T.internal[:x], T.tmpMulVec)
	had!(T.tmpMulVec, T.left[:y])
	BLAS.gemv!('N', ComplexOne, T.tmpMulVec, onesp, ComplexOne, y)

	# IR
	NFFT.nfft!(T.NFFTx, S.IR, T.tmpMulVec)
	scale!(T.internal[:x], T.tmpMulVec)
	had!(T.tmpMulVec, T.right[:y])
	BLAS.gemv!('N', ComplexOne, T.tmpMulVec, onesp, ComplexOne, y)

	# LI
	NFFT.nfft!(T.NFFTy, S.LI, T.tmpMulVecT)
	transpose!(T.tmpMulVec, T.tmpMulVecT)
	scale!(T.internal[:y], T.tmpMulVec)
	had!(T.tmpMulVec, T.left[:x])
	BLAS.gemv!('N', ComplexOne, T.tmpMulVec, onesp, ComplexOne, y)

	# RI
	NFFT.nfft!(T.NFFTy, S.RI, T.tmpMulVecT)
	transpose!(T.tmpMulVec, T.tmpMulVecT)
	scale!(T.internal[:y], T.tmpMulVec)
	had!(T.tmpMulVec, T.right[:x])
	BLAS.gemv!('N', ComplexOne, T.tmpMulVec, onesp, ComplexOne, y)

	isuniform(T) || had!(y, get(T.weights))

	return y
end

function Base.A_mul_B!(y::DenseVector{Complex{Float64}}, T::Freq2Wave2D, x::DenseVector{Complex{Float64}})
	size(T,2) == length(x) || throw(DimensionMismatch())

	X = reshape(x, wsize(T))
	A_mul_B!(y, T, X)

	return y
end

@compat function Base.:*(T::Freq2Wave, x::DenseArray)
	if !isa(x, Array{Complex{Float64}})
		x = map(Complex{Float64}, x)
	end

	y = Array{Complex{Float64}}( size(T,1) )
	A_mul_B!(y, T, x)

	return y
end

function Base.Ac_mul_B!(z::DenseVector{Complex{Float64}}, T::Freq2NoBoundaryWave1D, v::DenseVector{Complex{Float64}})
	(M = size(T,1)) == length(v) || throw(DimensionMismatch())
	size(T,2) == length(z) || throw(DimensionMismatch())

	hadc!(T.tmpMulVec, v, T.internal)
	isuniform(T) || had!(T.tmpMulVec, get(T.weights))

	myfft_adjoint!(z, T.NFFT, T.tmpMulVec)

	return z
end

function Base.Ac_mul_B!(Z::DenseMatrix{Complex{Float64}}, T::Freq2NoBoundaryWave2D, v::DenseVector{Complex{Float64}})
	(M = size(T,1)) == length(v) || throw(DimensionMismatch())
	wsize(T) == size(Z) || throw(DimensionMismatch())

	for m in 1:M
		@inbounds T.tmpMulVec[m] = v[m] * conj(T.internal[:x][m]) * conj(T.internal[:y][m])
	end
	isuniform(T) || had!(T.tmpMulVec, get(T.weights))

	myfft_adjoint!(Z, T.NFFT, T.tmpMulVec)

	return Z
end

function Base.Ac_mul_B!(z::DenseVector{Complex{Float64}}, T::Freq2BoundaryWave1D, v::DenseVector{Complex{Float64}})
	(M = size(T,1)) == length(v) || throw(DimensionMismatch())
	size(T,2) == length(z) || throw(DimensionMismatch())

	zleft, zint, zright = split(z, van_moment(T))

	# Boundary contributions are first as they don't use T.internal
	copy!(T.tmpMulVec, v)
	isuniform(T) || had!(T.tmpMulVec, get(T.weights))
	Ac_mul_B!( zleft, T.left, T.tmpMulVec )
	Ac_mul_B!( zright, T.right, T.tmpMulVec )

	# Internal scaling function
	hadc!(T.tmpMulVec, T.internal)
	myfft_adjoint!(zint, T.NFFT, T.tmpMulVec)

	return z
end

function Base.Ac_mul_B!(Z::DenseMatrix{Complex{Float64}}, T::Freq2BoundaryWave2D, v::DenseVector{Complex{Float64}})
	(M = size(T,1)) == length(v) || throw(DimensionMismatch())
	(N = wsize(T)) == size(Z) || throw(DimensionMismatch())

	# As in Ac_mul_B! for Freq2BoundaryWave1D
	if isuniform(T)
		copy!(T.weigthedVec, v)
	else
		had!(T.weigthedVec, v, get(T.weights))
	end

	vm = van_moment(T)
	S = split(Z, vm)

	# Internal coefficients
	for m in 1:M
		@inbounds T.tmpMulcVec[m] = T.weigthedVec[m] * conj(T.internal[:x][m]) * conj(T.internal[:y][m])
	end
	myfft_adjoint!(S.II, T.NFFT, T.tmpMulcVec)

	# ------------------------------------------------------------
	# Left blocks: All use 'L' for the y coordinate

	# LL
	conj!(T.tmpMulVec, T.left[:y])
	scale!(T.weigthedVec, T.tmpMulVec)
	Ac_mul_B!( S.LL, T.left[:x], T.tmpMulVec )

	# RL: Reuse T.tmpMulVec
	Ac_mul_B!( S.RL, T.right[:x], T.tmpMulVec )

	# IL
	conj!(T.tmpMulcVec, T.internal[:x])
	scale!(T.tmpMulcVec, T.tmpMulVec)
	NFFT.nfft_adjoint!( T.NFFTx, T.tmpMulVec, S.IL )

	# ------------------------------------------------------------
	# Middle blocks: All use 'I' for the y coordinate

	# LI
	conj!(T.tmpMulVec, T.left[:x])
	scale!(T.weigthedVec, T.tmpMulVec)

	conj!(T.tmpMulcVec, T.internal[:y])
	scale!(T.tmpMulcVec, T.tmpMulVec)
	transpose!(T.tmpMulVecT, T.tmpMulVec)
	NFFT.nfft_adjoint!( T.NFFTy, T.tmpMulVecT, S.LI )

	# RI: Reuse T.tmpMulcVec
	conj!(T.tmpMulVec, T.right[:x])
	scale!(T.weigthedVec, T.tmpMulVec)
	scale!(T.tmpMulcVec, T.tmpMulVec)
	transpose!(T.tmpMulVecT, T.tmpMulVec)
	NFFT.nfft_adjoint!( T.NFFTy, T.tmpMulVecT, S.RI )

	# ------------------------------------------------------------
	# Right blocks: All use 'R' for the y coordinate

	# LR
	conj!(T.tmpMulVec, T.right[:y])
	scale!(T.weigthedVec, T.tmpMulVec)
	Ac_mul_B!( S.LR, T.left[:x], T.tmpMulVec )

	# RR: Reuse T.tmpMulVec
	Ac_mul_B!( S.RR, T.right[:x], T.tmpMulVec )

	# IR
	conj!(T.tmpMulcVec, T.internal[:x])
	scale!(T.tmpMulcVec, T.tmpMulVec)
	NFFT.nfft_adjoint!( T.NFFTx, T.tmpMulVec, S.IR )

	return Z
end

function Base.Ac_mul_B!(z::DenseVector{Complex{Float64}}, T::Freq2Wave2D, v::DenseVector{Complex{Float64}})
	size(T,2) == length(z) || throw(DimensionMismatch())

	Z = reshape(z, wsize(T))
	Ac_mul_B!(Z, T, v)

	return z
end

function Base.Ac_mul_B(T::Freq2Wave, v::AbstractVector)
	if !isa(v, Array{Complex{Float64}})
		v = map(Complex{Float64}, v)
	end

	z = Array{Complex{Float64}}(size(T,2))
	Ac_mul_B!(z, T, v)

	return z
end

@compat function Base.:\(T::Freq2Wave, Y::AbstractMatrix)
	length(Y) == (M = size(T,1)) || throw(DimensionMismatch())

	y = flatten_view(Y)
	x = T \ y
end

@compat function Base.:\(T::Freq2Wave, y::AbstractVector)
	if !isa(y, Array{Complex{Float64}})
		y = map(Complex{Float64}, y)
	end

	# Non-uniform samples: Scale observations
	if !isuniform(T)
		y .*= get(T.weights)
	end

	x0 = zeros(Complex{Float64}, wsize(T))
	x = cgnr(T, y, x0)
	#x, h = lsqr(T, y)
end

@doc """
	collect(Freq2Wave) -> Matrix

Return the full change of basis matrix.

In 2D, the reconstruction grid is sorted by the `y` coordinate, i.e., the order is
(1,1),
(2,1),
(3,1),
(1,2),
(2,2),
(3,2)
etc
"""->
function Base.collect(T::Freq2NoBoundaryWave2D)
	M = size(T,1)
	Nx, Ny = wsize(T)
	F = Array{Complex{Float64}}(M, Nx*Ny)

	phi = T.internal[:x] .* T.internal[:y]

	offsetx = div(Nx, 2) + 1
	offsety = div(Ny, 2) + 1
	row_idx = 0
	for ny in 1:Ny, nx in 1:Nx
		row_idx += 1
		for m in 1:M
			@inbounds F[m,row_idx] = phi[m]*cis( -twoπ*((nx-offsetx)*T.NFFT.x[1,m] + (ny-offsety)*T.NFFT.x[2,m]) )
		end
	end

	if !isuniform(T)
		scale!(get(T.weights), F)
	end

	return F
end

function Base.collect(T::Freq2BoundaryWave2D)
	M = size(T,1)
	Nx, Ny = wsize(T)
	F = Array{Complex{Float64}}(M, Nx*Ny)

	phix = Array{Complex{Float64}}(M)
	phiy = similar(phix)

	row_idx = 0
	for ny in 1:Ny
		unsafe_FourScaling!(phiy, T, ny, :y)
		for nx in 1:Nx
			unsafe_FourScaling!(phix, T, nx, :x)
			had!(phix, phiy)
			F[:,row_idx+=1] = phix
		end
	end

	if !isuniform(T)
		scale!(get(T.weights), F)
	end

	return F
end

@doc """
	unsafe_FourScaling!(phi, T::Freq2BoundaryWave2D, n::Int, d::Symbol)

Replace `phi` with the `n`'th "column" from dimension `d` of `T`.
"""->
function unsafe_FourScaling!(phi::Vector{Complex{Float64}}, T::Freq2BoundaryWave2D, n::Integer, d::Symbol)
	M = length(phi)
	N = 2^wscale(T)
	p = van_moment(T)

	if p < n <= N-p
		offset = div(N, 2) + 1
        D = (d == :x ? 1 : 2)
		for m in 1:M
			@inbounds phi[m] = T.internal[d][m]*cis( -twoπ*(n-offset)*T.NFFT.x[D,m] )
		end
	elseif 1 <= n <= p
		unsafe_copy!( phi, 1, T.left[d], (n-1)*M+1, M )
	else
		unsafe_copy!( phi, 1, T.right[d], (n-N+p-1)*M+1, M )
	end
end

# ------------------------------------------------------------
# Common

function Base.size(T::Freq2Wave)
	( size(T,1), size(T,2) )
end

function Base.size(T::Freq2Wave, d::Integer)
	if d == 1
        prod(T.NFFT.M)
	elseif d == 2
		prod( wsize(T) )
	else
		throw(AssertionError())
	end
end

van_moment(T::Freq2Wave1D) = hasboundary(T) ? size(T.left,2) : van_moment(T.wavename)
van_moment(T::Freq2Wave2D) = hasboundary(T) ? size(T.left[:x],2)::Int64 : van_moment(T.wavename)

function Base.eltype(::Freq2Wave)
	return Complex{Float64}
end

function Base.show(io::IO, T::Freq2Wave)
	D = dim(T)
	println(io, D, "D change of basis matrix")

	isuniform(T) ?  U = " " : U = " non-"
	M = size(T,1)
	println(io, "From: ", M, U, "uniform frequency samples")

	D == 1 ? N = size(T,2) : N = wsize(T)
	print(io, "To: ", N, " ", T.wavename, " wavelets")
end


text:
using BioMedQuery.Entrez

#convert to a dictionary
efetch_dict = eparse_from_file("efetch_sample.xml")

articles = []

# Decide type of article based on structrure of efetch
if haskey(efetch_dict, "PubmedArticle")
    articles = efetch_dict["PubmedArticle"]
else
    error("Not a PubMed search")
end

for (i, xml_article) in enumerate(articles)

    article = PubMedArticle(xml_article)
    #PMIDs should not be null
    @test !isnull(article.pmid)
    mesh_heading_list = MeshHeadingList(xml_article)

    for heading in mesh_heading_list
        @test !isnull(heading.descriptor_name)
    end

    if i==1
        @test article.title.value == "Five Tips to Building a Successful Sleep Practice."
        @test article.pmid.value == 27483622
        @test article.authors[1][:LastName].value == "Poss"
        @test mesh_heading_list[1].descriptor_name.value == "bruxism"
        @test length(mesh_heading_list)== 6
    end
end


text:
### BEGIN: cpu_times
"""
    cpu_times(cpu_id=nothing)

Return system CPU times as a named tuple.

If `cpu_id::Integer` is provided, the result is for that specific CPU (logical) core.

# Outputs

Every attribute represents the seconds the CPU has spent in the given mode. The attributes
availability varies depending on the platform:

- `user`: time spent by normal processes executing in user mode; on Linux this also includes guest
  time
- `system`: time spent by processes executing in kernel mode
- `idle`: time spent doing nothing

Platform-specific fields:

- `nice` (UNIX): time spent by niced (prioritized) processes executing in user mode; on Linux this
  also includes guest_nice time
- `iowait` (Linux): time spent waiting for I/O to complete. This is not accounted in idle time
  counter.
- `irq` (Linux, BSD): time spent for servicing hardware interrupts
- `softirq` (Linux): time spent for servicing software interrupts
- `steal` (Linux 2.6.11+): time spent by other operating systems running in a virtualized
  environment
- `guest` (Linux 2.6.24+): time spent running a virtual CPU for guest operating systems under the
  control of the Linux kernel
- `guest_nice` (Linux 3.2.0+): time spent running a niced guest (virtual CPU for guest operating
  systems under the control of the Linux kernel)
- `interrupt` (Windows): time spent for servicing hardware interrupts ( similar to “irq” on UNIX)
- `dpc` (Windows): time spent servicing deferred procedure calls (DPCs); DPCs are interrupts that
  run at a lower priority than standard interrupts.

See also: For detailed explaination of each fields, please refer to `/proc/stat` section in `man
proc`.

# Examples

```julia
julia> using PSUtil

julia> PSUtil.cpu_times()
(user = 5243197, nice = 24218, system = 1518950, idle = 621296540, iowait = 114376, irq = 0, softirq = 491843, steal = 0, guest = 0, guest_nice = 0)
```

To get cpu times stat per thread:

```julia
julia> PSUtil.cpu_times.(0:Sys.CPU_THREADS-1);
```
"""
cpu_times

function cpu_times(::Linux, cpu_id::Union{Nothing, Integer}=nothing)
    if isnothing(cpu_id)
        record = readline("/proc/stat")
    else
        n = Sys.CPU_THREADS
        cpu_id <= n || throw(ArgumentError("cpu_id $cpu_id exceeds the number of CPU cores $n"))
        full_records = readlines("/proc/stat") # Enhancement: This get read multiple times if `cpu_times.(0:4)`
        record = full_records[findfirst(startswith("cpu$cpu_id"), full_records)]
    end

    k = (:user, :nice, :system, :idle, :iowait, :irq, :softirq, :steal, :guest, :guest_nice)
    v = map(x->parse(Int, x), split(record)[2:end])
    return (; zip(k, v)...)
end

### END: cpu_times

### BEGIN: cpu_percent

"""
    cpu_percent(cpu_id=nothing; interval=1) -> percentage
    cpu_percent(last_idle, last_total, cpu_id=nothing)

Return a float representing the current system-wide CPU utilization as a percentage.

# Inputs

To query usage for specific CPU core, you can optionally pass `cpu_id::Integer`.

If `lastidle` and `lasttotal` are provided, it calculates the CPU utilization since when the numbers
are recorded. Otherwise, it will calculate the CPU utilization from now to `interval` seconds later,
which is a blocking routine.

# Examples

```julia
julia> using PSUtil

julia> PSUtil.cpu_percent() # blocks the current thread by `interval` seconds
1.2432656444260282
```

`last_idle` and `last_total` can be book kept so that you can formulate a daemon script, for
example, the following script reports the CPU utilization every 5 seconds:

```julia
last_idle, last_total = 0, 0
while true
    utilization, last_idle, last_total = cpu_percent(last_idle, last_total)
    println(utilization)
    sleep(5)
end
```

To query cpu percentage for all cpus asynchronously (so that the process only get blocked by
`interval` seconds instead of `n * interval` seconds):

```julia
results = map(t->t.result,
    @sync map(0:PSUtil.cpu_count()-1) do i
        t = @async PSUtil.cpu_percent(i; interval=1)
    end
)
```
"""
cpu_percent

function cpu_percent(sys::AbstractSystem, cpu_id=nothing; interval=1)
    _, idle, total = cpu_percent(sys, 0, 0, cpu_id)
    sleep(interval)
    utilization, _, _ = cpu_percent(sys, idle, total)
    return utilization
end

function cpu_percent(::AbstractSystem, last_idle, last_total, cpu_id::Union{Nothing, Integer}=nothing; interval=nothing)
    isnothing(interval) || @warn "keyword `interval` is not used in non-block version."

    records = cpu_times(cpu_id)
    idle, total = records.idle, sum(records)
    Δidle, Δtotal = idle - last_idle, total - last_total

    utlization = 100 * (1 - Δidle/Δtotal)
    return utlization, idle, total
end

### END: cpu_percent

### BEGIN: cpu_count

"""
    cpu_count(;logical=true)

Return the number of logical CPUs in the system (same as `Sys.CPU_THREADS`).

`logical` cores means the number of physical cores multiplied by the number of threads that can run on
each core (this is known as Hyper Threading). If `logical` is `false` return the number of physical
cores only (Hyper Thread CPUs are excluded) or `nothing` if undetermined.

# Examples

Example on a system having 2 physical hyper-thread CPU cores:

```julia
julia> import psutil
julia> psutil.cpu_count()
4
julia> psutil.cpu_count(logical=false)
2
```

See also: For more detailed CPU information, you can use [`CpuId.jl`](https://github.com/m-j-w/CpuId.jl).
"""

cpu_count(::AbstractSystem; logical=true) = logical ?  CpuId.cputhreads() : CpuId.cpucores()

### END: cpu_count


text:
MLJBase.evaluate!(mach, resampling=MLJBase.CV(; shuffle = true), measure=MLJBase.rms, operation=MLJBase.predict_mean)
particles = predict_particles(mach, X; response = :y)
@test rms_expected(particles, truth.y) == Statistics.mean(rms_distribution(particles, truth.y))
mach_1 = MLJBase.machine(model, X, truth.y)
@test_throws ErrorException predict_particles(mach_1, X; response = :y)


text:
@testset "Line - basic" begin
    c = RCanvas()
    box1 = RBox()
    box2 = RBox()

    @test isempty(box1.outgoing_edges)
    @test isempty(box1.incoming_edges)
    @test isempty(box2.outgoing_edges)
    @test isempty(box2.incoming_edges)

    line = Line(box1, box2)

    @test line.from == box1
    @test line.to == box2
end


text:
export Nodes

# ~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--
# Nodes API interface
# ~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--~--
module Nodes

export
    visit, enter_node, exit_node,
    set_position!, calc_transform!

using ..Geometry
using ..Ranger
using ..Math

using Base.Math:
    deg2rad, rad2deg

# --------------------------------------------------------
# Life cycle events
# --------------------------------------------------------
function enter_node end
function exit_node end

function transition end
function get_replacement end
function take_replacement end
function visit end
function draw end
function interpolate end

function has_parent end
function is_nil end
function has_replacement end
function has_parent end

function update end
function calc_transform! end

function get_children end

function io_event end

# ---------- INCLUDES --------------------------------------------------
include("transition_properties.jl")
include("timing_properties.jl")
include("transform_properties.jl")

include("node_data.jl")
include("node_nil.jl")
include("node_trivial.jl")
include("node_stack.jl")

include("scenes/scenes.jl")
include("scenes/scene_nil.jl")

include("node_manager.jl")

include("filters/filters.jl")

include("space_mappings.jl")
include("detection.jl")

# -----------------------------------------------------------
# Abstract Nodes
# -----------------------------------------------------------
function is_nil(node::Ranger.AbstractNode)
    node.base.id == 0
end

function has_parent(node::Ranger.AbstractNode)
    !is_nil(node.base.parent)
end

using ..Rendering

# visit() traverse "downward" the heirarchy while space-mappings traverse upward.
function visit(node::Ranger.AbstractNode, context::Rendering.RenderContext, interpolation::Float64)
    # println("AbstractNode::visit ", node);
    if !is_visible(node)
        return;
    end

    Rendering.save!(context)

    # Because position and angles are dependent
    # on lerping we perform interpolation first.
    interpolate(node, interpolation)

    aft = calc_transform!(node)

    Rendering.apply!(context, aft)

    children = get_children(node)
    if children ≠ nothing
        draw(node, context)

        for child in children
            visit(child, context, interpolation)
        end
    else
        draw(node, context)
    end

    Rendering.restore!(context)

    # println("End visit ----------------------------------------------")
end

function calc_transform!(node::Ranger.AbstractNode)
    tr = node.transform

    if is_dirty(node)
        # println("calc_transform dirty: ", node.base.name)
        make_translate!(tr.aft, tr.position.x, tr.position.y)

        if tr.rotation ≠ 0.0
            rotate!(tr.aft, tr.rotation)
        end

        if tr.scale.x ≠ 1.0 || tr.scale.y ≠ 1.0
            scale!(tr.aft, tr.scale.x, tr.scale.y)
        end

        Math.invert!(tr.aft, tr.inverse)
    end

    tr.aft
end

function get_transform(node::Ranger.AbstractNode)
    node.transform.aft
end

function interpolate(node::Ranger.AbstractNode, interpolation::Float64)
    # println("AbstractNode::interpolate ", node);
end

function draw(node::Ranger.AbstractNode, context::RenderContext)
    # println("AbstractNode::draw ", node);
end

function enter_node(node::Ranger.AbstractNode, man::NodeManager)
    # println("AbstractNode::enter ", node);
end

function exit_node(node::Ranger.AbstractNode, man::NodeManager)
    # println("AbstractNode::exit ", node);
end

function io_event(node::Ranger.AbstractNode, event::Ranger.AbstractIOEvent)
    # println("AbstractScene::io_event ", event)
end

function is_visible(node::Ranger.AbstractNode)
    node.base.visible
end

function is_dirty(node::Ranger.AbstractNode)
    node.base.dirty
end

# ripple dirty flag
function set_dirty!(node::Ranger.AbstractNode, dirty::Bool)
    node.base.dirty = dirty
end

function ripple_dirty!(node::Ranger.AbstractNode, dirty::Bool)
    children = get_children(node)
    if children ≠ nothing
        for child in children
            ripple_dirty!(child, dirty)
        end
    end

    # println("ripple_dirty: ", node.base.name)
    set_dirty!(node, dirty)
end

function set_position!(node::Ranger.AbstractNode, x::T, y::T) where {T <: AbstractFloat}
    Geometry.set!(node.transform.position, x, y)
    ripple_dirty!(node, true);
end

function set_rotation_in_degrees!(node::Ranger.AbstractNode, rotation::T) where {T <: AbstractFloat}
    node.transform.rotation = deg2rad(rotation);

    ripple_dirty!(node, true);
end

function set_scale!(node::Ranger.AbstractNode, s::T) where {T <: AbstractFloat}
    Geometry.set!(node.transform.scale, s, s);
    ripple_dirty!(node, true);
end

function set_nonuniform_scale!(node::Ranger.AbstractNode, sx::T, sy::T) where {T <: AbstractFloat}
    Geometry.set!(node.transform.scale, sx, sy)
    ripple_dirty!(node, true);
end

# -----------------------------------------------------------
# Abstract Scenes
# -----------------------------------------------------------
function is_nil(node::Ranger.AbstractScene)
    node.base.id == 0
end

function has_replacement(node::Ranger.AbstractScene)
    !is_nil(node.replacement)
end

function has_parent(node::Ranger.AbstractScene)
    !is_nil(node.base.parent)
end

function enter_node(node::Ranger.AbstractScene, man::NodeManager)
    println("AbstractScene::enter ", node);
end

function exit_node(node::Ranger.AbstractScene, man::NodeManager)
    println("AbstractScene::exit ", node);
end

# --------------------------------------------------------
# Grouping
# --------------------------------------------------------
function get_children(node::Ranger.AbstractScene)
    nothing
end

function get_children(node::Ranger.AbstractNode)
    nothing
end

# --------------------------------------------------------
# Polygons
# --------------------------------------------------------
function get_bucket(node::Ranger.AbstractNode)
    node.polygon.mesh.bucket # return AbstractArray{Geometry.Point{T}}
end

# --------------------------------------------------------------------------
# Timing
# --------------------------------------------------------------------------
function update(node::Ranger.AbstractNode, dt::Float64)
    # println("AbstractNode::update : ", node)
end

# --------------------------------------------------------------------------
# Misc debugging
# --------------------------------------------------------------------------
import Base.show

function Base.show(io::IO, node::Ranger.AbstractNode)
    print(io, "|'", node.base.name, "' (", node.base.id, ")|");
end

function print_tree(node::Ranger.AbstractNode)
    println("---------- Tree ---------------")
    print_branch(UInt32(0), node)

    children = get_children(node)
    if children ≠ nothing
        print_sub_tree(children, UInt32(1))
    end
    println("-------------------------------");
end

function print_sub_tree(children::Array{Ranger.AbstractNode,1}, level::UInt32)
    for child in children
        sub_children = get_children(child)
        if sub_children ≠ nothing
            print_branch(level, child)
            print_sub_tree(sub_children, UInt32(level + 1))
        else
            print_branch(level, child)
        end
    end
end

function print_branch(level::UInt32, node::Ranger.AbstractNode)
    for _ in 1:level
        print("  ")
    end
    println(node)
end

end # End Module --------------------------------------------------------

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
using Test
using OMEinsum
using CuArrays

CuArrays.allowscalar(false)

@testset "cuda einsum" begin
    a = [randn(fill(3, i)...) for i=1:4]
    ca = a .|> CuArray
    ixs = ((1,2), (2,3))
    xs = (ca[2], ca[2])
    @test loop_einsum!(EinCode(ixs, (1,3)), xs, zeros(3,3) |> CuArray, OMEinsum.get_size_dict(ixs, xs)) ≈ ca[2]*ca[2]
    for f in [ein"ij,jk->ik", ein"ii->", ein"ijj ->i", ein"ij,ik,il->jkl", ein"ii->i", ein"ijl->i", ein"i->ii", ein"ij,jk,kl->il", ein"ij,ij,ij -> ij"]
        @show f
        cins = map(ix->ca[length(ix)], OMEinsum.getixs(f))
        ins = map(ix->a[length(ix)], OMEinsum.getixs(f))
        @test f(cins...) isa CuArray
        @test f(ins...) ≈ f(cins...)
    end
end

@testset "fallback - getindex IR error" begin
    a = rand(ComplexF64,2,2,2)
    ca = CuArray(a);
    @test Array(ein"npu,por,dom,lmn -> urdl"(ca,ca,ca,ca)) ≈ ein"npu,por,dom,lmn -> urdl"(a,a,a,a)
end

@testset "isbatchmul" begin
    for (ixs, iy) in [(((1,2), (2,3)), (1,3)), (((1,2,3), (2,3)), (1,3)),
                        (((7,1,2,3), (2,4,3,7)), (1,4,3)),
                        (((3,), (3,)), (3,)), (((3,1), (3,)), (3,1))
                        ]
        xs = ([randn(ComplexF64, fill(4,length(ix))...) |> CuArray for ix in ixs]...,)
        @test OMEinsum.batched_contract(ixs[1], xs[1], ixs[2], xs[2], iy) |> Array ≈ loop_einsum(EinCode(ixs, iy), xs, OMEinsum.get_size_dict(ixs, xs))
        @test EinCode(ixs, iy)(xs...) |> Array ≈ loop_einsum(EinCode(ixs, iy), xs, OMEinsum.get_size_dict(ixs, xs))
    end
end


text:
using PatternFolds
using Test

@testset "PatternFolds.jl" begin
    include("vector.jl")
    include("intervals.jl")
end


text:
include("../../utils.jl")
using Cthulhu

module TortureTest

    # Scalability parameters - set here to roughly match the real world usecase
    # though that use case does more work.
    const NFIELDS = 888
    const NSET = 387
    const NREPS = 40

    using Random
    using Base.Meta

    struct Maybe{T}
        x::T
        is_missing::Bool
    end
    Base.convert(::Type{Maybe{T}}, x::Maybe{T}) where {T} = x
    Base.convert(::Type{Maybe{T}}, x::T) where {T} = Maybe{T}(x, true)
    Base.convert(::Type{Maybe{T}}, x) where {T} = Maybe{T}(convert(T, x), true)
    Base.eltype(::Type{Maybe{T}}) where {T} = T
    Maybe(x::T) where {T} = Maybe(x, false)

    @eval struct TestCase
        $([:($(Symbol("x$i"))::Maybe{$(rand(Bool) ? Float64 : Int64)}) for i = 1:NFIELDS]...)
    end

    @eval TestCase($(Expr(:parameters, [Expr(:kw, Symbol("x$i"), :(Maybe(convert(eltype(fieldtype(TestCase, $(quot(Symbol("x$i"))))), $(rand(Bool) ? Float64(i) : i))))) for i = 1:NFIELDS]...))) =
        TestCase($([Symbol("x$i") for i = 1:NFIELDS]...))

    @eval function run_test1(x)
        $([:(TestCase(;$([:($(Symbol("x$i")) = $(rand() < 0.05 ? (:x) : rand(Bool) ? Float64(i) : i)) for i = randperm(NFIELDS)[1:NSET]]...))) for _ = 1:NREPS]...)
    end

    @eval function run_test2(x)
        $([:(TestCase(;$([:($(Symbol("x$i")) = $(rand(Bool) ? Float64(i) : i)) for i = randperm(NFIELDS)[1:NSET]]...))) for _ = 1:NREPS]...)
    end

end

@gctime Cthulhu.mkinterp(TortureTest.run_test1, Tuple{Int})[end].inInference

max_stars_repo_name: d-netto/GCBenchmarks

max_issues_repo_name: d-netto/GCBenchmarks

max_forks_repo_name: d-netto/GCBenchmarks


text:
using Oceananigans.Grids
using Oceananigans.Utils: prettytime, hours, day, days, years

using Statistics
using JLD2
using Printf
using GLMakie

function geographic2cartesian(λ, φ; r=1)
    Nλ = length(λ)
    Nφ = length(φ)

    λ = repeat(reshape(λ, Nλ, 1), 1, Nφ)
    φ = repeat(reshape(φ, 1, Nφ), Nλ, 1)

    λ_azimuthal = λ .+ 180  # Convert to λ ∈ [0°, 360°]
    φ_azimuthal = 90 .- φ   # Convert to φ ∈ [0°, 180°] (0° at north pole)

    x = @. r * cosd(λ_azimuthal) * sind(φ_azimuthal)
    y = @. r * sind(λ_azimuthal) * sind(φ_azimuthal)
    z = @. r * cosd(φ_azimuthal)

    return x, y, z
end

function visualize_makie(output_prefix)

    filepath = output_prefix * ".jld2"

    file = jldopen(filepath)

    Nx = file["grid/Nx"]
    Ny = file["grid/Ny"]
    Lλ = file["grid/Lx"]
    Lφ = file["grid/Ly"]
    Lz = file["grid/Lz"]

    # A spherical domain
    grid = LatitudeLongitudeGrid(size = (Nx, Ny, 1),
                                 longitude = (-180, 180),
                                 latitude = (-Lφ/2, Lφ/2),
                                 z = (-Lz, 0))

    iterations = parse.(Int, keys(file["timeseries/t"]))

    λ, ϕ, r = nodes((Face, Face, Center), grid)
    x, y, z = geographic2cartesian(λ, ϕ, r=1.01)

    iter = Node(0)
    ζ′ = @lift file["timeseries/ζ/" * string($iter)][:, :, 1]

    fig = Figure(resolution = (2000, 2000))

    clims = @lift begin
        scale = 1.5e-1
        min_clim = 2e-5
        max_ζ = maximum(abs, file["timeseries/ζ/" * string($iter)][:, :, 1])
        scale * max_ζ < min_clim && (-min_clim, min_clim)
        scale .* (-max_ζ, max_ζ)
    end

    ax = fig[:, :] = LScene(fig) # make plot area wider
    wireframe!(ax, Sphere(Point3f0(0), 1f0), show_axis=false)
    surface!(ax, x, y, z, color=ζ′, colormap=:blues, colorrange=clims)
    rotate_cam!(ax.scene, (π/4, π/6, 0))
    zoom!(ax.scene, (0, 0, 0), 5, false)

    plot_title = @lift @sprintf("Vertical vorticity in decaying, rotating, barotropic turbulence at time = %s",
                                prettytime(file["timeseries/t/" * string($iter)]))
    supertitle = fig[0, :] = Label(fig, plot_title, textsize=50)

    display(fig)

    record(fig, output_prefix * ".mp4", iterations, framerate=12) do i
        @info "Plotting iteration $i of $(iterations[end])..."
        iter[] = i
    end

    return nothing
end

output_prefix = "rotating_freely_decaying_barotropic_turbulence_fly_Nx1080_Ny480"

visualize_makie(output_prefix)


text:
using SICMModelingToolkit
using Test

@testset "SICMModelingToolkit.jl" begin
    # Write your tests here.
end

max_stars_repo_name: zoemcc/SICMModelingToolkit

max_issues_repo_name: zoemcc/SICMModelingToolkit

max_forks_repo_name: zoemcc/SICMModelingToolkit


text:
Dict(
    2 => Symbol[:dot,:floor,:max,:any,:ldexp,:ceil,:normalize,:all,:log,:mod,:cross,:trunc,:round,:min]
    3 => Symbol[:round,:ceil,:trunc,:floor,:fma,:any,:all,:clamp]
    1 => Symbol[:atanh,:tanh,:exp2,:floor,:max,:any,:sin,:tan,:modf,:step,:log2,:acosh,:ceil,:normalize,:length,:atan,:all,:log,:isnan,:trunc,:asinh,:transpose,:abs,:frexp,:sqrt,:round,:sinh,:sign,:acos,:asin,:cosh,:exp,:min,:isinf,:cos]
)


text:
module Test

g(x::Float64, y) = 2x + y
g(x, y::Float64) = x + 2y

function largest_seq(seq::String)
    l = length(seq)
    kmer = seq[1]
    k = 1
    for i = 2:l
        if seq[i] == seq[i-1]
            kmer = string(kmer,seq[i])
        else
            kmer = ""
        end
    end

end

end


text:
module RocketMapOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: map()" begin

    println("Testing: operator map()")

    run_proxyshowcheck("Map", map(Any, d -> d))

    run_testset([
        (
            source      = from(1:5) |> map(Int, d -> d ^ 2),
            values      = @ts([ 1, 4, 9, 16, 25, c ]),
            source_type = Int
        ),
        (
            source      = from(1:5) |> map(Float64, d -> convert(Float64, d)),
            values      = @ts([ 1.0, 2.0, 3.0, 4.0, 5.0, c ]),
            source_type = Float64
        ),
        (
            source      = completed() |> map(Int, d -> d + 1),
            values      = @ts(c),
            source_type = Int
        ),
        (
            source      = faulted(Int, "e") |> map(String, d -> string(d)),
            values      = @ts(e("e")),
            source_type = String
        ),
        (
            source      = never() |> map(Int, d -> 1),
            values      = @ts(),
            source_type = Int
        ),
        (
            source      = from(1:5) |> map(Int, d -> 1.0), # Invalid output: Float64 instead of Int
            values      = @ts(),
            source_type = Int,
            throws      = Exception
        ),
        (
            source      = from(1:5) |> safe() |> map(Int, d -> 1.0), # Invalid output: Float64 instead of Int
            values      = @ts(e),
            source_type = Int
        )
    ])

end

end


text:
######################################################################################################
# Copy from https://bitbucket.org/wmkoolen/tidnabbil/src/master/purex_games_paper/ by Wouter Koolen  #
######################################################################################################
include("binary_search.jl");

struct Gaussian
    σ2;
end

# convenience
Gaussian() = Gaussian(1);

struct Bernoulli
end

struct Exponential
end

struct Poisson
end

# KL divergence
rel_entr(x, y) = x==0 ? 0. : x*log(x/y);
dx_rel_entr(x, y) = x==0 ? 0. : log(x/y);
dy_rel_entr(x, y) = -x/y;

d(expfam::Gaussian,    μ, λ) = (μ-λ)^2/(2*expfam.σ2);
d(expfam::Bernoulli,   μ, λ) = max(0, rel_entr(μ, λ) + rel_entr(1-μ, 1-λ));
d(expfam::Exponential, μ, λ) = λ == 0 ? Inf : μ/λ - log(μ/λ) - 1;
d(expfam::Poisson,     μ, λ) = rel_entr(μ, λ) - μ + λ;

dµ_d(expfam::Gaussian,    μ, λ) = (µ-λ)/expfam.σ2
dµ_d(expfam::Bernoulli,   μ, λ) = dx_rel_entr(μ, λ) - dx_rel_entr(1-μ, 1-λ);
dµ_d(expfam::Exponential, μ, λ) = 1/λ - 1/µ;
dµ_d(expfam::Poisson,     μ, λ) = dx_rel_entr(μ, λ) - 1;

dλ_d(expfam::Gaussian,    μ, λ) = (λ-µ)/expfam.σ2
dλ_d(expfam::Bernoulli,   μ, λ) = dy_rel_entr(μ, λ) - dy_rel_entr(1-μ, 1-λ);
dλ_d(expfam::Exponential, μ, λ) = -µ/λ^2 + 1/λ;
dλ_d(expfam::Poisson,     μ, λ) = dy_rel_entr(μ, λ) + 1;
# solution for λ in \min_λ d(μ, λ) - λ x
# i.e. λ satisfying (λ - μ)ϕ''(λ) == x

invh(expfam::Gaussian,    μ, x) = μ + x*expfam.σ2;
invh(expfam::Bernoulli,   μ, x) = 2μ/(1-x+sqrt((x-1)^2 + 4*x*μ));
invh(expfam::Exponential, μ, x) = x > 0 ? Inf : 2μ/(1 + sqrt(1 - 4*x*μ));
invh(expfam::Poisson,     μ, x) = μ/(1-x);

sample(rng, expfam::Gaussian,    μ) = μ + sqrt(expfam.σ2)*randn(rng);
sample(rng, expfam::Bernoulli,   μ) = rand(rng) ≤ μ;
sample(rng, expfam::Exponential, μ) = randexp(rng)*μ;

# upward and downward confidence intervals
dup(expfam::Gaussian, μ, v) = μ + sqrt(2*expfam.σ2*v);
ddn(expfam::Gaussian, μ, v) = μ - sqrt(2*expfam.σ2*v);

function dup(expfam::Bernoulli, μ, v)
    μ == 1 ? 1. : binary_search(λ -> d(expfam, μ, λ) - v, μ, 1);
end

function ddn(expfam::Bernoulli, μ, v)
    μ == 0 ? 0. : binary_search(λ -> v - d(expfam, μ, λ), 0, μ);
end

function dup(expfam::Exponential, μ, v)
    binary_search(λ -> d(expfam, μ, λ) - v, μ, μ*exp(v+1));
end

function ddn(expfam::Exponential, μ, v)
    binary_search(λ -> v - d(expfam, μ, λ), 0, μ);
end


text:
struct Extractor
    f::Function
    name::Symbol
end

function names(extractors::Array{Extractor,1})
    return map(e -> e.name, extractors)
end

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
# Factorize two different matrices (one real, one complex) and solve for multiple rhs and free memory
@everywhere using MUMPS
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end
using Compat
include("getDivGrad.jl");

A  = getDivGrad(24,20,22);
A2 = getDivGrad(30,31,33);

n  = size(A,1);
n2 = size(A2,1);
A  = A + im*spdiagm(rand(n),0);

nrhs = 10;
rhs = randn(n,nrhs) + im*randn(n,nrhs);
rhs2 = randn(n2,nrhs);

@test nworkers()>=3

println("Factorize complex matrix on worker 1")
tic();
F1 = remotecall_fetch(factorMUMPS,workers()[1], A,1)
toc();
dump(F1)

println("Factorize real matrix on worker 2")
tic();
F2 = remotecall_fetch(factorMUMPS,workers()[2],A2,1)
toc();
dump(F2)

println("Solve complex system on worker 1")
tic();
x = remotecall_fetch(applyMUMPS,workers()[3],F1,rhs)
toc();
err = zeros(nrhs)
for i=1:nrhs
        err[i] =  norm(A*x[:,i]-rhs[:,i]) / norm(rhs[:,i])
end
@test maximum(err) < 1e-14

println("Solve real system on worker 2")
tic();
x2 = remotecall_fetch(applyMUMPS,workers()[2],F2,rhs2)
toc();
err = zeros(nrhs)
for i=1:nrhs
        err[i] =  norm(A2*x2[:,i]-rhs2[:,i]) / norm(rhs2[:,i]);
end
@test maximum(err) < 1e-14

println("Free memory on worker 1")
destroyMUMPS(F1)
# remotecall_fetch(destroyMUMPS,workers()[1],F1);

println("Free memory on worker 2")
destroyMUMPS(F2)
# remotecall_fetch(destroyMUMPS,workers()[2],F2);

println("DONE!")


text:
using QDNN
using Yao
using Test

@testset "QDNN.jl" begin
    # Write your own tests here.
end

@testset "circuit.jl" begin
    ql = encoder_circuit(2,2)
end

@testset "layer.jl" begin
    cir_e = encoder_circuit(3, 5)
    cir_t = transform_circuit(3, 8)
    w = parameters(cir_t)
    H1_Y = [put(3, i=>Y) for i = 1:3]
    b = zeros(3)

    ql = QNNL{3, Float64}(cir_e, cir_t, w, H1_Y, b)
    x = rand(15)
    forward(ql, x)
    back_propagation(ql, x)
end

@testset "model.jl" begin

end


text:
const deps = joinpath(@__DIR__, "..", "deps")
const url = "https://github.com/FluxML/Metalhead.jl/releases/download/Models"

function testimgs end
function valimgs end

function getweights(name)
  mkpath(deps)
  cd(deps) do
    if name == "vgg19.bson"
        isfile(name) || Base.download("$url/$name", name)
    else
        loc = "https://github.com/FluxML/Metalhead.jl/releases/download/v0.1.1"
        isfile(name) || Base.download("$loc/$name", name)
    end
  end
end

function weights(name)
  getweights(name)
  BSON.load(joinpath(deps, name))
end

load_img(im::AbstractMatrix{<:Color}) = im
load_img(str::AbstractString) = load(str)
load_img(val::ValidationImage) = load_img(val.img)

# Resize an image such that its smallest dimension is the given length
function resize_smallest_dimension(im, len)
  reduction_factor = len/minimum(size(im)[1:2])
  new_size = size(im)
  new_size = (
      round(Int, size(im,1)*reduction_factor),
      round(Int, size(im,2)*reduction_factor),
  )
  if reduction_factor < 1.0
    # Images.jl's imresize() needs to first lowpass the image, it won't do it for us
    im = imfilter(im, KernelFactors.gaussian(0.75/reduction_factor), Inner())
  end
  return imresize(im, new_size)
end

# Take the len-by-len square of pixels at the center of image `im`
function center_crop(im, len)
  l2 = div(len,2)
  adjust = len % 2 == 0 ? 1 : 0
  return im[div(end,2)-l2:div(end,2)+l2-adjust,div(end,2)-l2:div(end,2)+l2-adjust]
end

function preprocess(im::AbstractMatrix{<:AbstractRGB})
  # Resize such that smallest edge is 256 pixels long
  im = resize_smallest_dimension(im, 256)

  # Center-crop to 224x224
  im = center_crop(im, 224)

  # Convert to channel view and normalize (these coefficients taken
  # from PyTorch's ImageNet normalization code)
  μ = [0.485, 0.456, 0.406]
  σ = [0.229, 0.224, 0.225]
  im = (channelview(im) .- μ)./σ

  # Convert from CHW (Image.jl's channel ordering) to WHCN (Flux.jl's ordering)
  # and enforce Float32, as that seems important to Flux
  return Float32.(permutedims(im, (3, 2, 1))[:,:,:,:].*255)
end

preprocess(im) = preprocess(load(im))
preprocess(im::AbstractMatrix) = preprocess(RGB.(im))

forward(model, im) = vec(model(preprocess(RGB.(im))))

topk(ps::AbstractVector, k::Int = 5) = topk(1:length(ps), ps, k)
topk(classes, ps::AbstractVector, k::Int = 5) = sort(collect(zip(classes, ps)), by = x -> -x[2])[1:k]

make_fname(s::AbstractMatrix{<:Color}) = ""
make_fname(s::String) = s
make_fname(im::ValidationImage) = "$(im.set) Validation #$(im.idx)"

ground_truth(m, s::Union{AbstractMatrix, AbstractString}, result) = (nothing, 0.)
function ground_truth(m::ClassificationModel{Class}, im::ValidationImage, result) where {Class}
    if typeof(im.ground_truth) == Class
        (im.ground_truth, result[im.ground_truth.class])
    elseif hasmethod(convert, Tuple{Class, typeof(im.ground_truth)})
        (im.ground_truth, result[convert(Class, im.ground_truth).class])
    else
        (im.ground_truth, 0.0)
    end
end

function predict(model::ClassificationModel{Class}, im, k = 5) where {Class}
    fname = make_fname(im)
    img = load_img(im)
    result = forward(model, img)
    PredictionFrame(fname, img,
        Prediction([Class(x)=>y for (x,y) in topk(result, k)]),
        ground_truth(model, im, result)...)
end
classify(model::ClassificationModel, im) = Flux.onecold(forward(model, load_img(im)), labels(model))

function predict(model::ClassificationModel{Class},
        im::Vector{<:Union{AbstractMatrix, String, ValidationImage}}, k=5) where Class
    error("predict is not implicitly vectorized. Use broadcast syntax: predict.(model, imgs)")
end

struct Prediction
    sorted_predictions::Vector{Pair{ObjectClass, Float32}}
end

struct PredictionFrame
    filename::Union{String, Nothing}
    img
    prediction::Prediction
    # Ground truth, if known
    ground_truth::Union{Nothing, ObjectClass}
    ground_truth_confidence::Float32
end