text:
#= ForwardDiff and Interpolations removed as dependencies from this version.
"""
    draw_streamlines(pt, f_xy; physwidth = 10.0m, physheight = 4.0m, centered = true,
        h = 0.1s, nsteps = 10, probability = 0.0001)

where
    pt         point on screen corresponding to local origo
    centered   if true, physwidth and physheight extend symmetrically from origo
               if false, physwidth and physheight extend on the positive x and y axes.
                The streamlines are drawn to the right (+x) and above (+y) pt.
    physwidth, physheight Together with 'scale_to_pt', determines input to f_xy which correspond to pixels.
    f_xy is a function of quantities x and y where the range is also a quantity
    h is a step-size quantity, default time. If the integration variable is time,
        it defines the length of time betweeen points. It can be negative.
    nsteps is the number of steps to trace, forwards and backwards from each
        randomly selected point
    probability is the number of streamlines divided by length(xs)*length(ys)

Draws a few random streamlines, starting at a random number of points picked from (xs, ys).
Streamlines are found by starting at the randomly picked centre of their extent.
They may extend outside of the rectangle defined by xs and ys.
A better visualization might be based on where particles enter.
"""
function draw_streamlines(pt, f_xy; physwidth = 10.0m, physheight = 4.0m, centered = true,
                          h = 0.1s, nsteps = 10, probability = 0.0001)
    gsave()
    xs, ys = x_y_iterators_at_pixels(;physwidth, physheight, centered)

    cx = xs[1]::Quantity{Float64}
    cy = ys[1]
    @assert cx isa Quantity{Float64}
    @assert cy isa Quantity{Float64}
    # prepare a mutable but fast working buffer forward, and another backward.
    vxf = similar(SVector{nsteps}(fill(cx, nsteps)))
    vyf = similar(vxf)
    vxb = similar(vxf)
    vyb = similar(vxf)
    for (i::Int64, cx::Quantity{Float64}) in enumerate(xs), (j::Int64, cy::Quantity{Float64}) in enumerate(ys)
        vxf[1], vyf[1], vxb[1], vyb[1] = cx, cy, cx, cy
        if rand() < probability
            # Find the streamline path
            rk4_steps!(f_xy, vxf, vyf, h)
            rk4_steps!(f_xy, vxb, vyb, -h)
            # Reorder to most recent position -> oldest position
            vx = vcat(reverse(vxf), vxb)
            vy = vcat(reverse(vyf), vyb)
            # Draw it
            trace_diminishing(pt, vx, vy)
        end
    end
    grestore()
    nothing
end

"""
    sawtooth(x, maxabs)

Return x for x in the interval [-maxabs, maxabs]

Repeating between these values.
"""
sawtooth(x, maxabs) =   rem(x, 2maxabs, RoundNearest)

"""
    line_integral_convolution_complex(f_xy, n_xy, x_mid, y_mid, f_s, f_0)

Find the phase and magnitude for one pixel by projecting noise on 20 points on the streamline passing through it.
Combine Runge-Kutta 4th order and convolution.

Arguments:
f_xy    Function of x and y for visualization
n_xy    Noise function of x and y. The spectrum should be adapted for f_xy.
x_mid   Coordinate for which to evaluate. Streamlines forward and back are found for convolution.
y_mid   Coordinate for which to evaluate
f_s     Sampling frequency
f_0     Frequency of interest

The resulting complex number ŝ yields:
    phase = arg(ŝ)
    magnitude = hypot(ŝ)
"""
function line_integral_convolution_complex(f_xy, n_xy, x_mid, y_mid, f_s, f_0)
    x = x_mid
    y = y_mid
    # Unitless, normalized circular frequency of interest
    ω_0n = 2π ∙ f_0 / f_s
    # Frequency of interest position on the complex unit circle
    z = exp(-im * ω_0n)
    # z-transform aggregator, complex, start at mid point
    ŝ = n_xy(x, y) * z^0
    # Step length (time), negative for backwards.
    h = -1 / f_s
    for nn in -1:-1:-9
        x, y = rk4_step(f_xy, h, x, y)                         # This takes 0.37 of total time in the function
                                                               # 1.130 μs (1 allocation: 32 bytes)
        ŝ += (n_xy(x, y)) * z^-nn * (cos(π * nn / 20)^2)       # This takes 144/1539 = 0.09 of total time in the function
                                                               #  449.239 ns (7 allocations: 144 bytes)
    end
    # Jump back to start
    x = x_mid
    y = y_mid
    # Switch to walking forwards
    for nn in 1:10
        x, y = rk4_step(f_xy, -h, x, y)                       # This takes 637 / 1539 = 0.41 of total time in the function
        ŝ += n_xy(x, y) * z^-nn * (cos(π * nn / 20)^2)        # This take 163 / 1539 = 0.11 of total time in the function
    end
    ŝ
end

=#

"""
    noise_between_wavelengths(λ_min, λ_max, x)

Generate one value of 1D OpenSimplex noise. Output has  rougly linear amplitude spectrum
between wavelengths. There is also some noise at longer wavelengths.

Wave lengths ´λ_min´, ´λ_max´ and position ´x´ are typically quantities of dimension length.

Standard deviation for random ´x´ is 0.103.
"""
function noise_between_wavelengths(λ_min, λ_max, x)
    @assert λ_min < λ_max "λ_min < λ_max"
    octaves = round(Int, log2(λ_max / λ_min), RoundUp)
    noise( 2x / λ_max, detail = octaves, persistence = 0.7)
end

"""
    noise_between_wavelengths(λ_min, λ_max, x, y)

Generate one value of 2D OpenSimplex noise. Output has  rougly linear amplitude spectrum
between wavelengths. There is also some noise at longer wavelengths.

Wave lengths ´λ_min´, ´λ_max´ and positions ´x´, ´y´ are typically quantities of dimension length.

Standard deviation for random (´x´, ´y´) is 0.103.
"""
function noise_between_wavelengths(λ_min, λ_max, x, y)
    @assert λ_min <= λ_max "λ_min < λ_max"
    octaves = if λ_min != 0 * λ_min
        min(6, max(1, round(Int, log2(λ_max / λ_min), RoundUp)))
    else
        6
    end
    noise( 2x / λ_max, 2y/ λ_max, detail = octaves, persistence = 0.7)
end

"""
    noise_between_wavelengths(λ_min, λ_max, xs, normalize = true) -> Array{T,1} where T <: Number

Generate a vector of 1D OpenSimplex noise. Output has  rougly linear amplitude spectrum
between wavelengths. There is also some noise at longer wavelengths.

Wave lengths ´λ_min´, ´λ_max´ and position 'xs´ are typically quantities of dimension length.
Position xs is an iterator or vector.

If ´normalize = true´, as is default, the standard deviation for a long ´xs´ is  0.1604.

If ´normalize = false´, the standard deviation for random ´xs´ is 0.103.
"""
function noise_between_wavelengths(λ_min, λ_max, xs::T; normalize = true) where T <: Union{AbstractRange, Vector}
    no = [ noise_between_wavelengths(λ_min, λ_max, x) for x in xs]
    normalize ? normalize_datarange(no) : no
end
"""
    noise_between_wavelengths(λ_min, λ_max, xs, ys; normalize = true) -> Array{T,2} where T <: Number

Generate a matrix of 2D OpenSimplex noise. Output has  rougly linear amplitude spectrum
between wavelengths. There is also some noise at longer wavelengths.

Wave lengths ´λ_min´, ´λ_max´ and position iterators ´xs´, ´ys´ are typically quantities of dimension length.
Positions ´xs´ and ´ys´ are iterators or vectors.

If ´normalize = true´, as is default, the standard deviation for a large matrix is 0.1604.

If ´normalize = false´, the standard deviation for random ´xs´ is 0.103.

NOTE that the output matrix follows the image manipulation convention: xs correspond to output rows.
"""
function noise_between_wavelengths(λ_min, λ_max, xs::T, ys::T; normalize = true) where T <: Union{AbstractRange, Vector}
    no = [ noise_between_wavelengths(λ_min, λ_max, x, y) for y in ys, x in xs]
    normalize ? normalize_datarange(no) : no
end
#=
"""
    noise_for_lic(f_xy, xs, ys)

Create a simplex-based noise function to be used for visualizing f_xy.

Input: f_xy: (Q, Q)  → (Q, Q) or CQ: Function taking a tuple of quantities, outputs quantities
    xs and ys are iterators of Q, i.e. pixel cooridinates

Output: f: (Q, Q)  → R Function taking a tuple of quantities, outputs a real number

The range is [0.0, 1.0]

The spectrum amplitudes are adapted to f_xy in order to produce larger noise amplitude for longer wavelength.
This corresponds to larger velocities.
"""
function noise_for_lic(f_xy, xs, ys)
    # Velocity scale from, to
    v_min, v_max = lenient_min_max(f_xy, xs, ys)
    # Number of waves over one streamline
    waves_per_cycle = 2
    # One streamline is traced over a duration of
    Δt = 2.0s
    # Noise spectrum wavelengths
    λ_min, λ_max = Δt∙(v_min, v_max) / waves_per_cycle
    # Simplex noise matrix with linear spectrum - x in rows, y in columns
    mat_no = noise_between_wavelengths(λ_min, λ_max, xs, ys);
    # Make a function that interpolates between noise pixels:
    matrix_to_function(mat_no)
end

"""
    noise_for_lic(f_xy, matrix::BitMatrix; centered = true)

Create a simplex-based noise function to be used for visualizing f_xy.

Input:  f_xy: (Q, Q)  → (Q, Q) or CQ: Function taking a tuple of quantities, outputs quantities
        matrix:    Used for finding the size of output. The relation between pixel indices and spatial
            positions is implicitly given by ´scale_to_pt(m)´
        centered   keyword argument. If true, values at the centre of the matrix corresponds to (x,y) = (0, 0)m

Output: f: (Q, Q)  → R Function taking a tuple of quantities, outputs a real number

The range is [0.0, 1.0]

The spectrum amplitudes are adapted to f_xy in order to produce larger noise amplitude for longer wavelength.
This corresponds to larger velocities.
"""
function noise_for_lic(f_xy, matrix::BitMatrix; centered = true)
    ny, nx = size(matrix)
    physwidth = nx * scale_pt_to_unit(m)
    physheight =  ny * scale_pt_to_unit(m)
    xmin = -centered * physwidth / 2
    xmax = (1 - centered / 2) * physwidth
    ymin = -centered * physheight / 2
    ymax = (1 - centered / 2) * physheight
    xs = range(xmin, xmax, length = nx)
    ys = range(ymin, ymax, length = ny)
    noise_for_lic(f_xy, xs, ys)
end

"""
    lic_matrix_current(scene, framenumber; waves_per_cycle = 2)

When called in the context of a Scene, the complex convolution matrix is input via
the scene definition, see SceneOpts, Scene and Movie.

Extract one instance in time from the phase + amplitude info in the complex input matrix.

Output: R² Matrix of floating point numbers, representing the animation at this frame
    in this scene context, intended to be displayed as an image.

    If the maximum magnitude in scene.opts.phase_magnitude_matrix, output elements
    are in the range -1.0, 1.0.
"""
function lic_matrix_current(scene, framenumber)
    @assert scene.opts isa LicSceneOpts
    @assert scene.opts.phase_magnitude_matrix isa Array{Complex{Float64},2} string(scene.opts)
    frames_per_cycle = scene.opts.cycle_duration ∙ scene.opts.framerate
    # Where we are in the repeating cycle, [0, 1 - frame_duration]
    normalized_time = framenumber / frames_per_cycle
    # Number of waves over one streamline in the noise image
    CR² = scene.opts.phase_magnitude_matrix
    lic_matrix_current(CR², normalized_time; scene.opts.waves_per_cycle)
end

"""
    lic_matrix_current(CR²:Matrix{Complex{Float64}}, normalized_time;
                       waves_per_cycle = 2, zeroval = 0.0)

Input: CR²    Matrix of complex floating point numbers
                  phase and amplitude information. Used for rendering line-integral-convolution
                  showing a static vector field.
       normalized_time       a number in the range [0, 1], where 0 and 1 gives the same result, so pick one when looping frames.
       waves_per_cycle       Keyword argument, depending on the assumptions made when calculating CR².
       zeroval               Keyword argument, the value which replaces 0 in CR², could be used to emphasize streamlines

Output: R² Matrix of floating point numbers, representing the animation at this normalized time, intended to be displayed as an image.
Output elements are in the range -1.0, 1.0 (possible exception if 'zeroval' is outside.)
"""
function lic_matrix_current(CR²::Matrix{Complex{Float64}}, normalized_time;
                           waves_per_cycle = 2, zeroval = 0.0)
    @assert 0 <= normalized_time <= 1
    θ = 2π ∙ waves_per_cycle ∙ normalized_time
    mi, ma = lenient_min_max(CR²)
    map(CR²) do pixel_complex
        θ_pixel = sawtooth(θ + angle(pixel_complex), π )
        magn = hypot(pixel_complex) / ma
        magn > 0 ? magn∙cos(θ_pixel) : zeroval
    end
end

"""
    streamlines_matrix!(streamlinepixels::BitMatrix, v_xy; centered = true, targetdensity = 0.42)

Input:
    v_xy: (Q, Q)  → (Q, Q)    Function taking a tuple of quantities, outputs tuple of quantities
    streamlinepixels: A matrix of booleans or real numbers where 0 indicates unoccupied pixel. The matrix
        follows the image manipulation convention: xs correspond to output rows.
    centered          Bool      Keyword argument. If true, values at the centre of the matrix corresponds to (x,y) = (0, 0)m
    targetdensity     Real      Keyword argument. Fraction of pixels covered by streamlines.

Input matrix streamlinepixels is modified in place.
"""
function streamlines_matrix!(streamlinepixels::BitMatrix, v_xy; centered = true, targetdensity = 0.42)
    ny, nx = size(streamlinepixels)
    physwidth = nx * scale_pt_to_unit(m)
    physheight = ny * scale_pt_to_unit(m)
    xmin = -centered * physwidth / 2
    xmax = (1 - centered / 2) * physwidth
    ymin = -centered * physheight / 2
    ymax = (1 - centered / 2) * physheight
    xs = range(xmin, xmax, length = nx)
    ys = range(ymin, ymax, length = ny)
    # Find an approiate step length (time), negative for backwards.
    v_min, v_max = lenient_min_max(v_xy, xs, ys)
    maxsteplen = 3 * physwidth / nx
    h = maxsteplen / v_max
    #
    min_px_step = 1.0 * v_min / v_max
    bb = (first(xs), last(xs),  first(ys), last(ys))
    count = 0
    while true
        count +=1
        x = rand(xs)
        y = rand(ys)
        add_streamline!(v_xy, streamlinepixels, bb, h, min_px_step, x, y)
        # Check pixel density
        if mod(count, 10) == 0
            sum(streamlinepixels) / (nx*ny) >= targetdensity && break
        end
        # Check against a reasoable upper limit, 2500 * 4000 pixels
        count > 10^7 && break
    end
    streamlinepixels
end

"""
    streamlines_matrix(v_xy; physwidth = 10.0m, physheight = 4.0m, centered = true, targetdensity = 0.42)

Input:
    v_xy: (Q, Q)  → (Q, Q)      Function taking a tuple of quantities, outputs tuple of quantities
    centered          Bool      Keyword argument. If true, values at the centre of the matrix corresponds to (x,y) = (0, 0)m
    targetdensity     Real      Keyword argument. Fraction of pixels covered by streamlines.
    physwidth         Quantity  Together with 'scale_to_pt', determines input to v_xy which correspond to pixels.
    physheight        Quantity
    centered          Bool

Output:
    streamlinepixels: A matrix of booleans where false indicates unoccupied pixel. The matrix
        follows the image manipulation convention: xs correspond to output rows.
"""
function streamlines_matrix(v_xy; physwidth = 10.0m, physheight = 4.0m, centered = true, targetdensity = 0.42)
    xs, ys = x_y_iterators_at_pixels(;physwidth, physheight, centered)
    velocity_matrix = [v_xy(x, y) for y in ys, x in xs]
    streamlinepixels = falses(size(velocity_matrix))
    streamlines_matrix!(streamlinepixels, v_xy)
    streamlinepixels
end

"""
    view_neighbourhood_float_index(matrix, fi, fj)

A 3x3 view around integer indices. Reduced view if at edge of matrix.
"""
function view_neighbourhood_float_index(matrix, fi, fj)
    i = round(Int, fi)
    j = round(Int, fj)
    imin = max(1, i - 1)
    imax = min(size(matrix)[1], i + 1)
    jmin = max(1, j - 1)
    jmax = min(size(matrix)[2], j + 1)
    view(matrix, imin:imax, jmin:jmax)
end

"""
    float_index(maxindex::Int, mini::T, maxi::T, q::T) where T

Linear floating index of q given that index 1 corresponds to 'mini' and 'maxindex' corresponds to 'maxi'.
"""
float_index(maxindex::Int, mini::T, maxi::T, q::T) where T =  1 + (maxindex - 1) * (q - mini) / (maxi - mini)

"""
    streamlineindices(v_xy, streamlinepixels, bb, h_abs, min_px_step, x, y)
    -> Set{Tuple{Int64, Int64}}

Return indices which would be occupied for a streamline through (x, y).
The streamline is traced forwards and backwards until
    - boundary is reached
    - flow stops / starts
    - a previously occupied pixel is hit by the streamline
    - a starting or ending point of a runge-kutta step hits within two pixels of an occupied pixel
    - the number of steps (forward or back) exceed both dimensions of streamlinepixels

Input:
    v_xy: (Q, Q)  → (Q, Q) or CQ: Function taking a tuple of quantities, outputs tuple of quantities
    streamlinepixels: A matrix of booleans or real numbers where 0 indicates unoccupied pixel. The matrix
        follows the image manipulation convention: xs correspond to output rows.
    bb: Tuple{Q, Q, Q, Q} = (minx, maxx, miny, maxy)
    h_abs:  Q, step length in time, for the Runge-Kutta 4th order method. Will be run with +h_abs and -h_abs.
    min_px_step: Exit criterion, a float number. If the process of the streamline falls below this value over one time step,
      the streamline is considered to end (or start) at that point. Measured in pixels.
    x:Q
    y:Q

where Q is a quantity.
"""
function streamlineindices(v_xy, streamlinepixels, bb, h_abs, min_px_step, x_mid, y_mid)
    ny, nx = size(streamlinepixels)
    maxsteps = 3 * max(ny, nx)
    countsteps = 0
    minx, maxx, miny, maxy = bb

    indices = Set{Tuple{Int64, Int64}}()

    x, y = x_mid, y_mid

    # Check neighborhood at start point for previous occupants

    startview = view_neighbourhood_float_index(streamlinepixels,
        float_index(ny, maxy, miny, y),
        float_index(nx, minx, maxx, x))
    if sum(startview) > 0
        return indices
    end
    # First trace forwards from x (positive time step h), then backwards
    for h_bothways in (h_abs, -h_abs)
        while true

            ro = float_index(ny, maxy, miny, y)
            ro < 1.0 && break

            ro > ny && break

            co = float_index(nx, minx, maxx, x)
            co < 1.0 && break
            co > nx && break

            # Find the end of this straight segment
            x1, y1 = rk4_step(v_xy, h_bothways, x, y)

            ro1 = float_index(ny, maxy, miny, y1)
            ro1 < 1.0 && break
            ro1 > ny && break

            co1 = float_index(nx, minx, maxx, x1)
            co1 < 1.0 && break
            co1 > nx && break

            # Is progress shorter than the given minimum step?
            if hypot(ro1 - ro, co1 - co) < min_px_step
                push!(indices, circle_algorithm(ny, nx, ro1, co1)...)
                break
            end
            # Is neighborhood at end of segment previously occupied?
            if sum(view_neighbourhood_float_index(streamlinepixels, ro1, co1)) > 0
                break
            end

            # Integer indices in this linear segment
            lineindices = bresenhams_line_algorithm(ny, nx, ro, co, ro1, co1)

            # Is any point touching previous occupants?
            pixeloccupied = sum(lineindices) do (i, j)
                    streamlinepixels[i, j] == true
                end > 0
            pixeloccupied && break

            # It's ok, store this straight line
            push!(indices, lineindices...)
            countsteps +=1

            # Avoid endless loop: Is the streamline too long (perhaps spiralling, circling?)
            if countsteps > maxsteps
                push!(indices, crossing_line_algorithm(ny, nx, ro1, co1)...)
                break
            end

            # Prepare next iteration
            x, y = x1, y1
        end
        countsteps = 0
        x, y = x_mid, y_mid
    end
    indices
end

"""
    set_true!(matrix, indices::Set{Tuple{Int64, Int64}})

Mofifies the matrix elements at indices to one / true
"""
function set_true!(matrix, indices::Set{Tuple{Int64, Int64}})
    occupied = one(eltype(matrix))
    for (i, j) in indices
        matrix[i, j] = occupied
    end
    nothing
end

"""
    add_streamline!(v_xy, occupied_matrix, bb, h_abs, min_px_step, x, y)

Modify occupied_matrix: set elements entered by a streamline through (x, y) to true or to 1 depending on type.

The streamline is traced forwards and backwards until
    - boundary is reached
    - flow stops / starts
    - a previously occupied pixel is hit by the streamline
    - a starting or ending point of a runge-kutta step hits within two pixels of an occupied pixel
    - the number of steps  (forward or back) exceed both dimensions of streamlinepixels

Input:
    v_xy: (Q, Q)  → (Q, Q) or CQ: Function taking a tuple of quantities, outputs tuple of quantities
    streamlinepixels: A matrix of booleans or real numbers where 0 indicates unoccupied pixel. The matrix
        follows the image manipulation convention: xs correspond to output rows.
    bb: Tuple{Q, Q, Q, Q} = (minx, maxx, miny, maxy), coordinates for 'v_xy'
    h_abs:  Q, step length in time, for the Runge-Kutta 4th order method. Will be run with +h_abs and -h_abs.
    min_px_step: Exit criterion, a float number. If the process of the streamline falls below this value over one time step,
        the streamline is considered to end (or start) at that point. Measured in pixels.
    x:Q  Position coordinate
    y:Q  Position coordinate

where Q is a quantity.
"""
function add_streamline!(v_xy, streamlinepixels, bb, h_abs, min_px_step, x, y)
    indices = streamlineindices(v_xy, streamlinepixels, bb, h_abs, min_px_step, x, y)
    set_true!(streamlinepixels, indices)
    nothing
end

"""
    convolution_matrix(matrix::Array{Quantity{Complex{Float64}, D, U}, 2};
                       centered = true; targetdensity = 0.42) where {D, U}

Input:
    matrix:           CQ²       Matrix of complex quantities, representing 2d vectors. Typically velocities.
    centered          Bool      Keyword argument. If true, values at the centre of the matrix corresponds to (x,y) = (0, 0)m
    targetdensity     Real      Keyword argument. Fraction of pixels covered by streamlines.

    The relation between pixel indices and spatial positions is implicitly given by ´scale_to_pt(m)´ together with 'centered'

Output:
    CR²

where
    Q   = Quantity
    CQ² = Complex quantity matrix
    CR² = Matrix of complex floating point numbers

Streamlines are randomly added until the density of pixels reaches targetdensity'
Output phase and amplitude information for 'streamlinepixels'.
Used for rendering line-integral-convolution showing a static vector field.
"""
#
function convolution_matrix(matrix::Array{Quantity{Complex{Float64}, D, U}, 2};
                           centered = true, targetdensity = 0.42) where {D, U}
    v_xy = matrix_to_function(matrix)
    streamlinepixels = falses(size(matrix))
    streamlines_matrix!(streamlinepixels::BitMatrix, v_xy; centered, targetdensity)
    convolution_matrix(v_xy, streamlinepixels; centered)
end

"""
    convolution_matrix(matrix::Array{Quantity{Complex{Float64}, D, U}, 2},
                       streamlinepixels::BitMatrix) where {D, U};
                       centered = true)

Input:
    matrix:           CQ²       Matrix of complex quantities, representing 2d vectors. Typically velocities.
    streamlinepixels: Bool²     Matrix of bits, where 'true' indicates that a line-integral-convolution should be found for this pixel.
    centered          Bool      Keyword argument. If true, values at the centre of the matrix corresponds to (x,y) = (0, 0)m

    Both matrix arguments must have the same number of rows and columns.
    The relation between pixel indices and spatial positions is implicitly given by ´scale_to_pt(m)´ together with 'centered'

Output:
    CR²

where
    Q   = Quantity
    CQ² = Complex quantity matrix
    CR² = Matrix of complex floating point numbers

Output phase and amplitude information for 'streamlinepixels'. Used for rendering line-integral-convolution showing a static vector field.
"""
function convolution_matrix(matrix::Array{Quantity{Complex{Float64}, D, U}, 2},
                            streamlinepixels::BitMatrix;
                            centered = true) where {D, U}
    @assert size(matrix) == size(streamlinepixels)
    # Make a continuous function, interpolated between pixels.
    fxy = matrix_to_function(matrix)
    convolution_matrix(fxy, streamlinepixels, centered = centered)
end

"""
    convolution_matrix(f_xy, streamlinepixels::BitMatrix; centered = true)

Input:
    f_xy: (Q, Q)  → (Q, Q)   Function taking two coordinates, outputs a tuple of quantities. The range of function f_xy typically represents velocities
    streamlinepixels:        Bool², matrix where 'true' indicates that a line-integral-convolution should be calculate for this pixel.
    centered                 Keyword argument. If true, values at the centre of the matrix corresponds to (x,y) = (0, 0)m

    The relation between pixel indices and spatial positions is implicitly given by ´scale_to_pt(m)´ and 'centered'

Output:
    CR²

where
    Q   = Quantity
    CR² = Matrix of complex floating point numbers

Output phase and amplitude information for 'streamlinepixels'. Used for rendering line-integral-convolution showing a static vector field.
"""
function convolution_matrix(f_xy, streamlinepixels::BitMatrix; centered = true)
    # Make a continuous function, interpolated between pixels. This is supposed to
    # be faster than the original function (which it may not be)
    fxy = if f_xy isa Extrapolation
        f_xy
    else
        ny, nx = size(streamlinepixels)
        physwidth = nx * scale_pt_to_unit(m)
        physheight = ny * scale_pt_to_unit(m)
        function_to_interpolated_function(f_xy; physwidth, physheight)
    end
    # Prepare a noise function:
    # Simplex based continuous noise function
    nxy = noise_for_lic(f_xy, streamlinepixels)
    # Convolute forward and back for every pixel
    convolution_matrix(f_xy, nxy, streamlinepixels, centered = centered)
end

"""
    convolution_matrix(f_xy, n_xy, streamlinepixels::BitMatrix; centered = true)

Input:
    f_xy: (Q, Q)  → (Q, Q)   Function taking two coordinates, outputs a tuple of quantities
    n_xy: (Q, Q)  → R        Noise function of x and y. The spectrum should be adapted for f_xy
    streamlinepixels:        Bool², matrix where 'true' indicates that a line-integral-convolution should be calculate for this pixel.
    centered                 keyword argument. If true, values at the centre of the matrix corresponds to (x,y) = (0, 0)m

    The relation between pixel indices and spatial positions is implicitly given by ´scale_to_pt(m)´

Output:
    CR²

where
    Q = Quantity
    R = Float64
    CR² = Matrix of complex floating point numbers

Output phase and amplitude information for 'streamlinepixels'. Used for rendering line-integral-convolution showing a static vector field.
"""
function convolution_matrix(f_xy, n_xy::Extrapolation, streamlinepixels::BitMatrix; centered = true)
    # used to be f_xy::Extrapolation
    ny, nx = size(streamlinepixels)
    physwidth = nx * scale_pt_to_unit(m)
    physheight = ny * scale_pt_to_unit(m)
    xmin = -centered * physwidth / 2
    xmax = (1 - centered / 2) * physwidth
    ymin = -centered * physheight / 2
    ymax = (1 - centered / 2) * physheight
    xs = range(xmin, xmax, length = nx)
    ys = range(ymax, ymin, length = ny)
    # Image processing convention: a column has a horizontal line of pixels
    CR² = zeros(Complex{Float64}, ny, nx)
    @assert eltype(xs) <: Quantity{Float64}
    @assert eltype(ys) <: Quantity{Float64}
    @assert n_xy(xs[1], ys[1]) isa Float64
    # Length of a streamline section in time, including forward and back
    Δt = 2.0s
    # Number of waves over one streamline
    wave_per_streamline = 2
    # Frequency of interest
    freq_0 = wave_per_streamline / Δt
    # Number of sample points per set (1 to 9 from tracking backward, 10 (start) to including 20 forward)
    n = 20
    # Sampling frequency.
    freq_s = n / Δt
    for (ro, y) in zip(1:ny, ys), (co, x) in zip(1:nx, xs)
        if streamlinepixels[ro, co]
            # Find the phase and magnitude for one pixel by projecting noise on 20 points on the streamline passing through it
            pv  = line_integral_convolution_complex(f_xy, n_xy, x, y, freq_s, freq_0)
            # Find the original indexes and update the image matrix
            CR²[ro, co] = pv
        end
    end
    CR²
end

"""
LicSceneOpts(cycle_duration, O;
            framerate = 30.0s⁻¹,
            phase_magnitude_matrix = Matrix{Complex{Float64}}(undef, 0, 0),
            data = missing,
            legend = missing,
            waves_per_cycle = 2,
            luminosity_variation = 0.67 )

This info may be passed through function calls to frame drawing functions.
LicSceneOpts optionally contains information needed to process frames based on line-integral-convolution,
ref. lic_matrix_current.

Input:
    cycle_duration::Time{Float64}
    framerate::typeof(30.0s⁻¹)
    O::Point                                    The origo in the local context of this scene.
    phase_magnitude_matrix::Matrix{ComplexF64}  See below
    data                                        See below
    legend                                      If 'data' exists, also supply a ColorLegend.
    waves_per_cycle                             Default is two, should be according to how 'phase_magnitude_matrix' was found.
    luminosity_variation                        Default is 0.67, possible range is [0,1], where 1.0 is full luminosity variation with convolution.

The phase at 'framenumber' and 'phase_magnitude_pixel' (from '..matrix') is found thus:
    θ = 2π ∙ waves_per_cycle ∙ framenumber / frames_per_cycle
    θ_pixel = sawtooth(θ + angle(pixel_complex), π )

Along with the magnitude, the phase information and 'framenumber' this is used to find the
current luminosity at a pixel.

If LicSceneOpts.data contains a matrix of the same size as 'phase_magnitude_matrix', that
can be used to find chroma and hue for the pixel. The conversion function is given by LicSceneOpts.legend.
"""
struct LicSceneOpts
    cycle_duration::Time{Float64}
    framerate::typeof(30.0s⁻¹)
    O::Point
    phase_magnitude_matrix::Matrix{ComplexF64}
    data
    legend
    waves_per_cycle
    luminosity_variation
end
# Constructor
function LicSceneOpts(cycle_duration, O;
            framerate = 30.0s⁻¹,
            phase_magnitude_matrix = Matrix{Complex{Float64}}(undef, 0, 0),
            data = missing,
            legend = missing,
            waves_per_cycle = 2,
            luminosity_variation = 0.67)
    LicSceneOpts(cycle_duration, framerate, O,
                 phase_magnitude_matrix, data, legend, waves_per_cycle, luminosity_variation)
end

"""
    color_matrix_current(scene, framenumber)

Input: scene     Defines options. See Movie, Scene, SceneOpts.

Output: Color² Matrix of colors representing the animation at this frame.
"""
function color_matrix_current(scene, framenumber)
    @assert scene isa Scene
    @assert scene.opts isa LicSceneOpts
    curmat = lic_matrix_current(scene, framenumber)
    speedmat = scene.opts.data
    color_matrix_mix(curmat, speedmat, scene.opts.legend, scene.opts.luminosity_variation)
end

"""
    color_matrix_mix(lumin_mat, quant_scalar_mat, legend, luminosity_variation)

Input:
    lumin_mat              Matrix{Float64} - determines luminosity at a pixel. Values in the range -1.0 to 1.0.
    quant_scalar_mat       Matrix of quantities, which are converted to colors using 'legend'
    legend                 <: ColorLegend
    luminosity_variation   Possible range is [0,1], recommended is 0.67. 1.0:full luminosity variation with convolution.

Output: Color² Matrix of colors representing the animation at this frame.
"""
function color_matrix_mix(lumin_mat, quant_scalar_mat, legend, luminosity_variation)
    map(legend.(quant_scalar_mat), lumin_mat) do col, lu
        # lu will be in the range -1.0 to 1.0
        lum = 50 * (lu * luminosity_variation + 1.0)
        # lum will be in the range 0 to 100.
        color_with_lumin(col, lum)
    end
end
=#


text:
# <hide>
function printsum(a)
    println(summary(a), ": ", repr(a))
end
# </hide>

# dicts can be initialised directly:
a1 = {1=>"one", 2=>"two"}
printsum(a1) #> Dict{Any,Any}: {2=>"two",1=>"one"}

# then added to:
a1[3]="three"
printsum(a1) #> Dict{Any,Any}: {2=>"two",3=>"three",1=>"one"}
# (note dicts cannot be assumed to keep their original order)

# dicts may also be created with the type explicitly set
a2 = Dict{Int64, String}()
a2[0]="zero"

# dicts, like arrays, may also be created from [comprehensions](http://julia.readthedocs.org/en/latest/manual/arrays/#comprehensions)
a3 = {i => @sprintf("%d", i) for i = 1:10}
printsum(a3)#> Dict{Any,Any}: {5=>"5",4=>"4",6=>"6",7=>"7",2=>"2",10=>"10",9=>"9",8=>"8",3=>"3",1=>"1"}

# as you would expect, Julia comes with all the normal helper functions
# for dicts, e.g., (haskey)[http://docs.julialang.org/en/latest/stdlib/base/#Base.haskey]
println(haskey(a1,1)) #> true

# which is equivalent to
println(1 in keys(a1)) #> true
# where [keys](http://docs.julialang.org/en/latest/stdlib/base/#Base.keys) creates an iterator over the keys of the dictionary

# similar to keys, (values)[http://docs.julialang.org/en/latest/stdlib/base/#Base.values] get iterators over the dict's values:
printsum(values(a1)) #> ValueIterator{Dict{Any,Any}}: ValueIterator{Dict{Any,Any}}({2=>"two",3=>"three",1=>"one"})

# use [collect](http://docs.julialang.org/en/latest/stdlib/base/#Base.collect) to get an array:
printsum(collect(values(a1)))#> 3-element Array{Any,1}: {"two","three","one"}

ext: jl

max_stars_repo_name: floswald/JuliaByExample

max_issues_repo_name: floswald/JuliaByExample

max_forks_repo_name: floswald/JuliaByExample


text:
@testset "Series: $TSeries" for (TSeries, T) in [(ArbSeries, Arb), (AcbSeries, Acb)]
    @testset "Length and degree" begin
        for (p, l) in [
            (zero(TSeries), 1),
            (one(TSeries), 1),
            (TSeries(T[0, 1, 0]), 3),
            (TSeries(T[0, 1], degree = 4), 5),
        ]
            @test firstindex(p) == 0
            @test size(p) == ()
            @test size(p, 1) == 1
            @test length(p) == Arblib.degree(p) + 1 == lastindex(p) + 1 == l
        end
    end

    @testset "Get and set coefficients" begin
        @test eltype(TSeries()) == T

        p = TSeries(T[i for i = 0:10])

        @test firstindex(p) == 0
        @test lastindex(p) == 10
        @test lastindex(TSeries(degree = 5)) == 5

        @test Arblib.coeffs(p) == p[:] == 0:10

        @test all(p[i] == i for i = 0:10)

        p[3] = T(7)
        @test p[3] == 7
        p[4] = 8
        @test p[4] == 8
        p[5] = π
        @test isequal(p[5], T(π))
        if T == Arb
            p[6] = ArbRefVector(T[9])[1]
        else
            p[6] = AcbRefVector(T[9])[1]
        end
        @test p[6] == 9

        q = TSeries(prec = 512)
        q[0] = π
        @test Arblib.rel_accuracy_bits(q[0]) > 500

        @test_throws BoundsError p[-1]
        @test_throws BoundsError p[11]

        p = TSeries(T[i for i = 0:10])
        @test all(Arblib.ref(p, i) == p[i] for i = 0:10)

        Arblib.set!(Arblib.ref(p, 0), 5)
        @test p[0] == 5

        @test_throws BoundsError Arblib.ref(p, -1)
        @test_throws BoundsError Arblib.ref(p, 11)
    end

    @testset "Constructors" begin
        @test TSeries() == TSeries(T[0]) == zero(TSeries) == zero(TSeries())
        @test TSeries(degree = 3) ==
              TSeries(T(0), degree = 3) ==
              TSeries(T[0], degree = 3) ==
              zero(TSeries(degree = 3))
        @test TSeries(T[1]) == one(TSeries) == one(TSeries())
        @test Arblib.isx(TSeries(T[0, 1]))
        @test TSeries(T[1, 2, 0]) != TSeries(T[1, 2])
        @test TSeries([5.0]) == TSeries([5]) == TSeries(T[5])
        @test TSeries(TSeries([1, 2])) == TSeries(ArbPoly([1, 2])) == TSeries([1, 2])

        @test precision(TSeries(degree = 1, prec = 64)) == 64
        @test precision(TSeries(0, degree = 1, prec = 64)) == 64
        @test precision(TSeries([0], prec = 64)) == 64
        @test precision(zero(TSeries(degree = 1, prec = 64))) == 64
        @test precision(one(TSeries(degree = 1, prec = 64))) == 64
        @test precision(TSeries(TSeries([1, 2], prec = 64))) == 64
        @test precision(TSeries(ArbPoly([1, 2], prec = 64))) == 64

        if TSeries == AcbSeries
            @test TSeries(AcbPoly([1, 2])) == TSeries(ArbSeries([1, 2])) == TSeries([1, 2])
            @test Arblib.degree(TSeries(ArbSeries(degree = 4))) == 4
            @test precision(TSeries(AcbPoly(prec = 64))) == 64
            @test precision(TSeries(ArbSeries(prec = 64))) == 64
        end
    end

    @testset "Arithmetic" begin
        p = TSeries([1, 2, 3])
        q = TSeries([2, 3, 4])

        @test p + q == TSeries([3, 5, 7])
        @test p - q == TSeries([-1, -1, -1])
        @test p * q == TSeries([2, 7, 16])
        @test q / p == TSeries([2, -1, 0])

        @test -p == TSeries([-1, -2, -3])

        @test inv(p) == TSeries([1, -2, 1])

        let p = setprecision(p, 80), q = setprecision(q, 90)
            @test precision(p + q) == 90
            @test precision(p - q) == 90
            @test precision(p * q) == 90
            @test precision(p / q) == 90
            @test precision(-p) == 80
            @test precision(inv(p)) == 80
        end

        if TSeries == AcbSeries
            q = ArbSeries([2, 3, 4])
            @test p + q == q + p == TSeries([3, 5, 7])
            @test p - q == TSeries([-1, -1, -1])
            @test q - p == TSeries([1, 1, 1])
            @test p * q == q * p == TSeries([2, 7, 16])
            @test p / ArbSeries(1, degree = 2) == p
            @test q / AcbSeries(1, degree = 2) == AcbSeries(q)
        end
    end

    @testset "Scalar arithmetic" begin
        p = TSeries([1, 2, 3])

        @test p + T(2) ==
              T(2) + p ==
              p + 2 ==
              2 + p ==
              p + 2.0 ==
              2.0 + p ==
              TSeries([3, 2, 3])
        @test p - T(2) == p - 2 == p - 2.0 == TSeries([-1, 2, 3])
        @test T(2) - p == 2 - p == 2.0 - p == TSeries([1, -2, -3])
        @test p * T(2) ==
              T(2) * p ==
              p * 2 ==
              2 * p ==
              p * 2.0 ==
              2.0 * p ==
              TSeries([2, 4, 6])
        @test p / T(2) == p / 2 == p / 2.0 == TSeries([0.5, 1, 1.5])
        @test T(2) / p == 2 / p == 2.0 / p == TSeries([2, -4, 2])

        # Test with zero polynomial
        @test zero(TSeries) + 1 == TSeries(1)
        @test zero(TSeries) - 1 == TSeries(-1)
        @test 1 - zero(TSeries) == TSeries(1)
        @test 1 * zero(TSeries) == TSeries()
        @test zero(TSeries) / 1 == TSeries()

        # Test that the normalisation works
        @test iszero(TSeries(-1) + 1)
        @test iszero(TSeries(1) - 1)
        @test iszero(1 - TSeries(1))

        let p = setprecision(p, 80)
            @test precision(p + T(2)) ==
                  precision(T(2) + p) ==
                  precision(2 + p) ==
                  precision(p + 2) ==
                  80
            @test precision(p - T(2)) == precision(p - 2) == 80
            @test precision(T(2) - p) == precision(2 - p) == 80
            @test precision(p * T(2)) ==
                  precision(T(2) * p) ==
                  precision(2 * p) ==
                  precision(p * 2) ==
                  80
            @test precision(p / T(2)) == precision(p / 2) == 80
            @test precision(T(2) / p) == precision(2 / p) == 80
        end

        if TSeries == ArbSeries
            @test p + im == p + Acb(im) == AcbSeries([1 + im, 2, 3])
            @test p - im == p - Acb(im) == AcbSeries([1 - im, 2, 3])
            @test im - p == Acb(im) - p == AcbSeries([im - 1, -2, -3])
            @test p * im == p * Acb(im) == AcbSeries([im, 2im, 3im])
            @test p / im == p / Acb(im) == AcbSeries([-im, -2im, -3im])
            @test im / p == Acb(im) / p == AcbSeries([im, -2im, im])

            let p = setprecision(p, 80)
                @test precision(p + im) == 80
                @test precision(p - im) == 80
                @test precision(im - p) == 80
                @test precision(p * im) == 80
                @test precision(p / im) == 80
            end
        end

        if TSeries == AcbSeries
            @test p + Arb(2) == Arb(2) + p == TSeries([3, 2, 3])
            @test p - Arb(2) == TSeries([-1, 2, 3])
            @test Arb(2) - p == TSeries([1, -2, -3])
            @test p * Arb(2) == Arb(2) * p == TSeries([2, 4, 6])
            @test p / Arb(2) == TSeries([0.5, 1, 1.5])
            @test Arb(2) / p == TSeries([2, -4, 2])
        end
    end

    @testset "Composition" begin
        p = TSeries([1, 2])
        q = TSeries([0, 3])

        @test Arblib.taylor_shift(p, T(2)) == Arblib.taylor_shift(p, 2) == TSeries([5, 2])
        @test Arblib.compose(p, q) == TSeries([1, 6])
        @test Arblib.revert(TSeries([0, 2])) == TSeries([0, 0.5])

        @test precision(Arblib.taylor_shift(setprecision(p, 80), T(2))) == 80
        @test precision(Arblib.compose(setprecision(p, 80), setprecision(q, 90))) == 90
        @test precision(Arblib.revert(TSeries([0, 2], prec = 80))) == 80

        @test_throws ArgumentError Arblib.compose(p, TSeries([1, 1]))
        @test_throws ArgumentError Arblib.revert(TSeries([0]))
        @test_throws ArgumentError Arblib.revert(TSeries([1, 1]))
        @test_throws ArgumentError Arblib.revert(TSeries([0, 0, 1]))
    end

    @testset "Evaluation" begin
        p = TSeries([1, 2, 3])

        @test p(Arb(2)) == p(2) == p(2.0) == 17
        @test p(Acb(1, 2)) == p(1 + 2im) == p(1.0 + 2.0im) == -6 + 16im

        @test Arblib.evaluate2(p, Arb(2)) == Arblib.evaluate2(p, 2) == (17, 14)
        @test Arblib.evaluate2(p, Acb(1, 2)) ==
              Arblib.evaluate2(p, 1 + 2im) ==
              (-6 + 16im, 8 + 12im)

        @test precision(p(T())) == precision(p(T(prec = 2precision(p)))) == precision(p)
        @test precision(TSeries(prec = 80)(T())) == 80
    end

    @testset "Differentiation and integration" begin
        p = TSeries([1, 2, 3])

        @test Arblib.derivative(p) == TSeries([2, 6])
        @test Arblib.integral(p) == TSeries([0, 1, 1, 1])

        @test precision(Arblib.derivative(p)) == precision(p)
        @test precision(Arblib.integral(p)) == precision(p)
        @test precision(Arblib.derivative(TSeries(prec = 80))) == 80
        @test precision(Arblib.integral(TSeries(prec = 80))) == 80

        @test Arblib.degree(Arblib.derivative(p)) == 1
        @test Arblib.degree(Arblib.integral(p)) == 3
    end

    @testset "Power methods" begin
        p = TSeries([1, 2, 3])
        q = TSeries([2, 3, 0])

        @test p^q == TSeries([1, 4, 16])

        @test p^T(2) ==
              p^Int(2) ==
              p^UInt(2) ==
              p^(2 // 1) ==
              p^2.0 ==
              p^2 ==
              TSeries([1, 4, 10])
        @test p^T(-1) == p^Int(-1) == p^(-1 // 1) == p^-1.0 == p^-1 == TSeries([1, -2, 1])

        @test p^(2 + 2im) == AcbSeries([1, 4 + 4im, 2 + 18im])

        @test 2^TSeries([1, 0]) == TSeries([2, 0])
        @test (2 + im)^TSeries([1, 0]) == AcbSeries([2 + im, 0])

        @test precision(setprecision(p, 80)^setprecision(q, 90)) == 90
        @test precision(setprecision(p, 80)^T(2)) == 80
    end

    @testset "Series methods" begin
        x = T(0.8)
        p = TSeries([x, 1])

        for f in [
            sqrt,
            log,
            log1p,
            exp,
            sin,
            cos,
            tan,
            atan,
            sinh,
            cosh,
            sinpi,
            cospi,
            Arblib.cotpi,
            Arblib.rsqrt,
        ]
            res = f(p)
            @test res[0] ≈ f(x)
            @test !iszero(res[1])
            iszero(res[1]) && @show f
        end

        @test all(isequal.(sincos(p), (sin(p), cos(p))))
        @test all(Arblib.coeffs(Arblib.sincospi(p)[1]) .≈ Arblib.coeffs(sinpi(p)))
        @test all(Arblib.coeffs(Arblib.sincospi(p)[2]) .≈ Arblib.coeffs(cospi(p)))
        @test all(isequal.(Arblib.sinhcosh(p), (sinh(p), cosh(p))))

        if TSeries == ArbSeries
            for f in [asin, acos, sinc]
                res = f(p)
                @test res[0] ≈ f(x)
                @test !iszero(res[1])
            end
        end
    end
end


text:
# Portions translated from SLICOT-Reference distribution
# Copyright (c) 2002-2020 NICONET e.V.
function run_mb03zd(datfile, io=stdout)
    ZERO = 0.0e0
    ONE = 1.0e0
    NIN = 5
    NOUT = 6
    NMAX = 200
    LDG = NMAX
    LDRES = 2*NMAX
    LDS = NMAX
    LDT = NMAX
    LDU1 = NMAX
    LDU2 = NMAX
    LDUS = 2*NMAX
    LDUU = 2*NMAX
    LDV1 = NMAX
    LDV2 = NMAX
    LDWORK = 3*NMAX*NMAX + 7*NMAX
    SELECT = Array{BlasBool,1}(undef, NMAX)
    SCALE = Array{Float64,1}(undef, NMAX)
    S = Array{Float64,2}(undef, LDS,NMAX)
    T = Array{Float64,2}(undef, LDT,NMAX)
    G = Array{Float64,2}(undef, LDG,NMAX)
    U1 = Array{Float64,2}(undef, LDU1,NMAX)
    U2 = Array{Float64,2}(undef, LDU2,NMAX)
    V1 = Array{Float64,2}(undef, LDV1,NMAX)
    V2 = Array{Float64,2}(undef, LDV2,NMAX)
    WR = Array{Float64,1}(undef, NMAX)
    WI = Array{Float64,1}(undef, NMAX)
    US = Array{Float64,2}(undef, LDUS,2*NMAX)
    UU = Array{Float64,2}(undef, LDUU,2*NMAX)
    IWORK = Array{BlasInt,1}(undef, 2*NMAX)
    LWORK = Array{BlasBool,1}(undef, 2*NMAX)
    DWORK = Array{Float64,1}(undef, LDWORK)
    RES = Array{Float64,2}(undef, LDRES,NMAX)
    f = open(datfile,"r")
    readline(f)
    vs = split(readline(f))
    N = parse(BlasInt, vs[1])
    ILO = parse(BlasInt, vs[2])
    WHICH =  vs[3][1]
    METH =  vs[4][1]
    STAB =  vs[5][1]
    BALANC =  vs[6][1]
    ORTBAL =  vs[7][1]
    if ( N<=0 || N>NMAX )
        @error "illegal N=$N"
    end
    if LSAME(WHICH,'S')
        vs = String[]
        _isz = N
        while length(vs) < _isz
            append!(vs, replace.(split(readline(f)),'D'=>'E'))
        end
       SELECT[1:_isz] .= parsex.(Bool, vs[1:_isz])
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       S[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       T[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    if LSAME(WHICH,'A') && LSAME(METH,'L')
        vs = String[]
        _isz,_jsz = (N,N)
        while length(vs) < _isz*_jsz
            append!(vs, replace.(split(readline(f)),'D'=>'E'))
        end
        for i in 1:_isz
            _i0 = (i-1)*_jsz
            G[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
        end
    end

    if LSAME(BALANC,'P') || LSAME(BALANC,'S') || LSAME(BALANC,'B')
        vs = String[]
        _isz = N
        while length(vs) < _isz
            append!(vs, replace.(split(readline(f)),'D'=>'E'))
        end
       SCALE[1:_isz] .= parsex.(Float64, vs[1:_isz])
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       U1[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       U2[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       V1[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end

    vs = String[]
    _isz,_jsz = (N,N)
    while length(vs) < _isz*_jsz
        append!(vs, replace.(split(readline(f)),'D'=>'E'))
    end
    for i in 1:_isz
       _i0 = (i-1)*_jsz
       V2[i,1:_jsz] .= parsex.(Float64, vs[_i0+1:_i0+_jsz])
    end
    close(f)
# interp call 1

    M, INFO = SLICOT.mb03zd!(WHICH, METH, STAB, BALANC, ORTBAL, SELECT, N, 2*N, ILO, SCALE, S, T, G, U1, U2, V1, V2, WR, WI, US, UU, IWORK)
    println(io, "M = $M")
    @test INFO == 0
    INFO == 0 || return

    println(io, "eigvals:")
    show(io, "text/plain", WR[1:N] + im *WI[1:N])
    println(io)

    if ( LSAME( STAB, 'S' )||LSAME( STAB, 'B' ) )
# interp output 1
        println(io, "US:")
        _nc = M
        _nr = 2*N
        show(io, "text/plain", US[1:_nr,1:_nc])
        println(io)

        if ( LSAME( ORTBAL, 'B' )||LSAME( BALANC, 'N' )|| LSAME( BALANC, 'P' ) )
# interp call 2

            # FOREIGN.dgemm!( 'Transpose', 'No Transpose', M, M, 2*N, ONE, US, LDUS, US, LDUS, ZERO, RES, LDRES )
            BLAS.gemm!('T','N',ONE,US[1:2N,1:M],US[1:2N,1:M],ZERO,view(RES,1:M,1:M))
            resid = norm(RES[1:M,1:M] - I)
            @test resid < 1e-12
            println(io, "orth. resid norm of US: $resid")
        end # if
# interp call 3

        # FOREIGN.dgemm!( 'Transpose', 'No Transpose', M, M, N, ONE, US, LDUS, US(N+1,1), LDUS, ZERO, RES, LDRES )
        BLAS.gemm!('T','N',ONE,US[1:N,1:M],US[N+1:2N,1:M],ZERO,view(RES,1:M,1:M))

# interp call 4

        # FOREIGN.dgemm!( 'Transpose', 'No Transpose', M, M, N, -ONE, US(N+1,1), LDUS, US, LDUS, ONE, RES, LDRES )
        BLAS.gemm!('T','N',-ONE,US[N+1:2N,1:M],US[1:N,1:M],ONE,view(RES,1:M,1:M))

        resid = norm(RES[1:M,1:M])
        @test resid < 1e-12
        println(io, "sympl. resid norm of US: $resid")
    end
    if ( LSAME( STAB, 'U' )||LSAME( STAB, 'B' ) )
        # interp output 2
        println(io, "UU:")
        _nc = M
        _nr = 2*N
        show(io, "text/plain", UU[1:_nr,1:_nc])
        println(io)

        if ( LSAME( ORTBAL, 'B' )||LSAME( BALANC, 'N' )|| LSAME( BALANC, 'P' ) )
# interp call 5

            #FOREIGN.dgemm!( 'Transpose', 'No Transpose', M, M, 2*N, ONE, UU, LDUU, UU, LDUU, ZERO, RES, LDRES )
            BLAS.gemm!('T','N',ONE,UU[1:2N,1:M],UU[1:2N,1:M],ZERO,view(RES,1:M,1:M))
            resid = norm(RES[1:M,1:M] - I)
            @test resid < 1e-12
            println(io, "orth. resid norm of UU: $resid")
        end # if
# interp call 6

        # FOREIGN.dgemm!( 'Transpose', 'No Transpose', M, M, N, ONE, UU, LDUU, UU(N+1,1), LDUU, ZERO, RES, LDRES )
        BLAS.gemm!('T','N',ONE,UU[1:N,1:M],UU[N+1:2N,1:M],ZERO,view(RES,1:M,1:M))

# interp call 7

        #FOREIGN.dgemm!( 'Transpose', 'No Transpose', M, M, N, -ONE, UU(N+1,1), LDUU, UU, LDUU, ONE, RES, LDRES )
        BLAS.gemm!('T','N',-ONE,UU[N+1:2N,1:M],UU[1:N,1:M],ONE,view(RES,1:M,1:M))

        resid = norm(RES[1:M,1:M])
        @test resid < 1e-12
        println(io, "sympl. resid norm of UU: $resid")

    end # if
end # run_mb03zd()


text:
struct SSPRK22Cache{uType,rateType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType # superfluous, only needed for callbacks...
  fsalfirst::rateType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end

u_cache(c::SSPRK22Cache) = ()
du_cache(c::SSPRK22Cache) = (c.k,c.fsalfirst)

struct SSPRK22ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::SSPRK22,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  SSPRK22Cache(u,uprev,k,tmp,fsalfirst,alg.stage_limiter!,alg.step_limiter!)
end

alg_cache(alg::SSPRK22,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = SSPRK22ConstantCache()

struct SSPRK33Cache{uType,rateType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType # superfluous, only needed for callbacks...
  fsalfirst::rateType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end

u_cache(c::SSPRK33Cache) = ()
du_cache(c::SSPRK33Cache) = (c.k,c.fsalfirst)

struct SSPRK33ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::SSPRK33,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  SSPRK33Cache(u,uprev,k,tmp,fsalfirst,alg.stage_limiter!,alg.step_limiter!)
end

alg_cache(alg::SSPRK33,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = SSPRK33ConstantCache()

struct SSPRK53Cache{uType,rateType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  fsalfirst::rateType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  α30::Float64
  α32::Float64
  α40::Float64
  α43::Float64
  α52::Float64
  α54::Float64
  β10::Float64
  β21::Float64
  β32::Float64
  β43::Float64
  β54::Float64
  c1::Float64
  c2::Float64
  c3::Float64
  c4::Float64
end

u_cache(c::SSPRK53Cache) = (c.tmp,)
du_cache(c::SSPRK53Cache) = (c.k,c.fsalfirst)

struct SSPRK53ConstantCache <: OrdinaryDiffEqConstantCache
  α30::Float64
  α32::Float64
  α40::Float64
  α43::Float64
  α52::Float64
  α54::Float64
  β10::Float64
  β21::Float64
  β32::Float64
  β43::Float64
  β54::Float64
  c1::Float64
  c2::Float64
  c3::Float64
  c4::Float64
end

function alg_cache(alg::SSPRK53,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  α30 = 0.355909775063327
  α32 = 0.644090224936674
  α40 = 0.367933791638137
  α43 = 0.632066208361863
  α52 = 0.237593836598569
  α54 = 0.762406163401431
  β10 = 0.377268915331368
  β21 = 0.377268915331368
  β32 = 0.242995220537396
  β43 = 0.238458932846290
  β54 = 0.287632146308408
  c1 = 0.377268915331368
  c2 = 0.754537830662736
  c3 = 0.728985661612188
  c4 = 0.699226135931670
  SSPRK53Cache(u,uprev,k,tmp,fsalfirst,alg.stage_limiter!,alg.step_limiter!,
                α30,α32,α40,α43,α52,α54,β10,β21,β32,β43,β54,c1,c2,c3,c4)
end

function alg_cache(alg::SSPRK53,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  α30 = 0.355909775063327
  α32 = 0.644090224936674
  α40 = 0.367933791638137
  α43 = 0.632066208361863
  α52 = 0.237593836598569
  α54 = 0.762406163401431
  β10 = 0.377268915331368
  β21 = 0.377268915331368
  β32 = 0.242995220537396
  β43 = 0.238458932846290
  β54 = 0.287632146308408
  c1 = 0.377268915331368
  c2 = 0.754537830662736
  c3 = 0.728985661612188
  c4 = 0.699226135931670
  SSPRK53ConstantCache(α30,α32,α40,α43,α52,α54,β10,β21,β32,β43,β54,c1,c2,c3,c4)
end

struct SSPRK63Cache{uType,rateType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  u₂::uType
  fsalfirst::rateType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  α40::Float64
  α41::Float64
  α43::Float64
  α62::Float64
  α65::Float64
  β10::Float64
  β21::Float64
  β32::Float64
  β43::Float64
  β54::Float64
  β65::Float64
  c1::Float64
  c2::Float64
  c3::Float64
  c4::Float64
  c5::Float64
end

u_cache(c::SSPRK63Cache) = (c.tmp,c.u₂)
du_cache(c::SSPRK63Cache) = (c.k,c.fsalfirst)

struct SSPRK63ConstantCache <: OrdinaryDiffEqConstantCache
  α40::Float64
  α41::Float64
  α43::Float64
  α62::Float64
  α65::Float64
  β10::Float64
  β21::Float64
  β32::Float64
  β43::Float64
  β54::Float64
  β65::Float64
  c1::Float64
  c2::Float64
  c3::Float64
  c4::Float64
  c5::Float64
end

function alg_cache(alg::SSPRK63,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  u₂ = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  α40 = 0.476769811285196
  α41 = 0.098511733286064
  α43 = 0.424718455428740
  α62 = 0.155221702560091
  α65 = 0.844778297439909
  β10 = 0.284220721334261
  β21 = 0.284220721334261
  β32 = 0.284220721334261
  β43 = 0.120713785765930
  β54 = 0.284220721334261
  β65 = 0.240103497065900
  c1 = 0.284220721334261
  c2 = 0.568441442668522
  c3 = 0.852662164002783
  c4 = 0.510854218958172
  c5 = 0.795074940292433
  SSPRK63Cache(u,uprev,k,tmp,u₂,fsalfirst,alg.stage_limiter!,alg.step_limiter!,
                α40,α41,α43,α62,α65,β10,β21,β32,β43,β54,β65,c1,c2,c3,c4,c5)
end

function alg_cache(alg::SSPRK63,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  α40 = 0.476769811285196
  α41 = 0.098511733286064
  α43 = 0.424718455428740
  α62 = 0.155221702560091
  α65 = 0.844778297439909
  β10 = 0.284220721334261
  β21 = 0.284220721334261
  β32 = 0.284220721334261
  β43 = 0.120713785765930
  β54 = 0.284220721334261
  β65 = 0.240103497065900
  c1 = 0.284220721334261
  c2 = 0.568441442668522
  c3 = 0.852662164002783
  c4 = 0.510854218958172
  c5 = 0.795074940292433
  SSPRK63ConstantCache(α40,α41,α43,α62,α65,β10,β21,β32,β43,β54,β65,c1,c2,c3,c4,c5)
end

struct SSPRK73Cache{uType,rateType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  u₁::uType
  fsalfirst::rateType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  α40::Float64
  α43::Float64
  α50::Float64
  α51::Float64
  α54::Float64
  α73::Float64
  α76::Float64
  β10::Float64
  β21::Float64
  β32::Float64
  β43::Float64
  β54::Float64
  β65::Float64
  β76::Float64
  c1::Float64
  c2::Float64
  c3::Float64
  c4::Float64
  c5::Float64
  c6::Float64
end

u_cache(c::SSPRK73Cache) = (c.tmp,c.u₁)
du_cache(c::SSPRK73Cache) = (c.k,c.fsalfirst)

struct SSPRK73ConstantCache <: OrdinaryDiffEqConstantCache
  α40::Float64
  α43::Float64
  α50::Float64
  α51::Float64
  α54::Float64
  α73::Float64
  α76::Float64
  β10::Float64
  β21::Float64
  β32::Float64
  β43::Float64
  β54::Float64
  β65::Float64
  β76::Float64
  c1::Float64
  c2::Float64
  c3::Float64
  c4::Float64
  c5::Float64
  c6::Float64
end

function alg_cache(alg::SSPRK73,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  u₁ = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  α40 = 0.184962588071072
  α43 = 0.815037411928928
  α50 = 0.180718656570380
  α51 = 0.314831034403793
  α54 = 0.504450309025826
  α73 = 0.120199000000000
  α76 = 0.879801000000000
  β10 = 0.233213863663009
  β21 = 0.233213863663009
  β32 = 0.233213863663009
  β43 = 0.190078023865845
  β54 = 0.117644805593912
  β65 = 0.233213863663009
  β76 = 0.205181790464579
  c1 = 0.233213863663009
  c2 = 0.466427727326018
  c3 = 0.699641590989027
  c4 = 0.760312095463379
  c5 = 0.574607439040817
  c6 = 0.807821302703826
  SSPRK73Cache(u,uprev,k,tmp,u₁,fsalfirst,alg.stage_limiter!,alg.step_limiter!,
                α40,α43,α50,α51,α54,α73,α76,β10,β21,β32,β43,β54,β65,β76,c1,c2,c3,c4,c5,c6)
end

function alg_cache(alg::SSPRK73,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  α40 = 0.184962588071072
  α43 = 0.815037411928928
  α50 = 0.180718656570380
  α51 = 0.314831034403793
  α54 = 0.504450309025826
  α73 = 0.120199000000000
  α76 = 0.879801000000000
  β10 = 0.233213863663009
  β21 = 0.233213863663009
  β32 = 0.233213863663009
  β43 = 0.190078023865845
  β54 = 0.117644805593912
  β65 = 0.233213863663009
  β76 = 0.205181790464579
  c1 = 0.233213863663009
  c2 = 0.466427727326018
  c3 = 0.699641590989027
  c4 = 0.760312095463379
  c5 = 0.574607439040817
  c6 = 0.807821302703826
  SSPRK73ConstantCache(α40,α43,α50,α51,α54,α73,α76,β10,β21,β32,β43,β54,β65,β76,c1,c2,c3,c4,c5,c6)
end

struct SSPRK83Cache{uType,rateType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType
  u₂::uType
  u₃::uType
  fsalfirst::rateType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  α50::Float64
  α51::Float64
  α54::Float64
  α61::Float64
  α65::Float64
  α72::Float64
  α73::Float64
  α76::Float64
  β10::Float64
  β21::Float64
  β32::Float64
  β43::Float64
  β54::Float64
  β65::Float64
  β76::Float64
  β87::Float64
  c1::Float64
  c2::Float64
  c3::Float64
  c4::Float64
  c5::Float64
  c6::Float64
  c7::Float64
end

u_cache(c::SSPRK83Cache) = (c.tmp,c.u₂,c.u₃)
du_cache(c::SSPRK83Cache) = (c.k,c.fsalfirst)

struct SSPRK83ConstantCache <: OrdinaryDiffEqConstantCache
  α50::Float64
  α51::Float64
  α54::Float64
  α61::Float64
  α65::Float64
  α72::Float64
  α73::Float64
  α76::Float64
  β10::Float64
  β21::Float64
  β32::Float64
  β43::Float64
  β54::Float64
  β65::Float64
  β76::Float64
  β87::Float64
  c1::Float64
  c2::Float64
  c3::Float64
  c4::Float64
  c5::Float64
  c6::Float64
  c7::Float64
end

function alg_cache(alg::SSPRK83,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  u₂ = similar(u)
  u₃ = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  α50 = 0.421366967085359
  α51 = 0.005949401107575
  α54 = 0.572683631807067
  α61 = 0.004254010666365
  α65 = 0.995745989333635
  α72 = 0.104380143093325
  α73 = 0.243265240906726
  α76 = 0.652354615999950
  β10 = 0.195804015330143
  β21 = 0.195804015330143
  β32 = 0.195804015330143
  β43 = 0.195804015330143
  β54 = 0.112133754621673
  β65 = 0.194971062960412
  β76 = 0.127733653231944
  β87 = 0.195804015330143
  c1 = 0.195804015330143
  c2 = 0.391608030660286
  c3 = 0.587412045990429
  c4 = 0.783216061320572
  c5 = 0.561833689734037
  c6 = 0.755247658555329
  c7 = 0.804195984669857
  SSPRK83Cache(u,uprev,k,tmp,u₂,u₃,fsalfirst,alg.stage_limiter!,alg.step_limiter!,
                α50,α51,α54,α61,α65,α72,α73,α76,β10,β21,β32,β43,β54,β65,β76,β87,c1,c2,c3,c4,c5,c6,c7)
end

function alg_cache(alg::SSPRK83,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  α50 = 0.421366967085359
  α51 = 0.005949401107575
  α54 = 0.572683631807067
  α61 = 0.004254010666365
  α65 = 0.995745989333635
  α72 = 0.104380143093325
  α73 = 0.243265240906726
  α76 = 0.652354615999950
  β10 = 0.195804015330143
  β21 = 0.195804015330143
  β32 = 0.195804015330143
  β43 = 0.195804015330143
  β54 = 0.112133754621673
  β65 = 0.194971062960412
  β76 = 0.127733653231944
  β87 = 0.195804015330143
  c1 = 0.195804015330143
  c2 = 0.391608030660286
  c3 = 0.587412045990429
  c4 = 0.783216061320572
  c5 = 0.561833689734037
  c6 = 0.755247658555329
  c7 = 0.804195984669857
  SSPRK83ConstantCache(α50,α51,α54,α61,α65,α72,α73,α76,β10,β21,β32,β43,β54,β65,β76,β87,c1,c2,c3,c4,c5,c6,c7)
end

struct SSPRK432Cache{uType,rateType,uArrayType,uEltypeNoUnits,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType # superfluous, only needed for callbacks
  fsalfirst::rateType
  utilde::uArrayType
  atmp::uEltypeNoUnits
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end

u_cache(c::SSPRK432Cache) = (c.utilde,c.atmp)
du_cache(c::SSPRK432Cache) = (c.k,c.fsalfirst)

struct SSPRK432ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::SSPRK432,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  utilde = similar(u,axes(u))
  atmp = similar(u,uEltypeNoUnits,axes(u))
  SSPRK432Cache(u,uprev,k,tmp,fsalfirst,utilde,atmp,alg.stage_limiter!,alg.step_limiter!)
end

alg_cache(alg::SSPRK432,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = SSPRK432ConstantCache()

struct SSPRK932Cache{uType,rateType,uArrayType,uEltypeNoUnits,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  tmp::uType # superfluous, only needed for callbacks
  fsalfirst::rateType
  utilde::uArrayType
  atmp::uEltypeNoUnits
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end

u_cache(c::SSPRK932Cache) = (c.utilde,c.atmp)
du_cache(c::SSPRK932Cache) = (c.k,c.fsalfirst)

struct SSPRK932ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::SSPRK932,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  utilde = similar(u,axes(u))
  atmp = similar(u,uEltypeNoUnits,axes(u))
  SSPRK932Cache(u,uprev,k,tmp,fsalfirst,utilde,atmp,alg.stage_limiter!,alg.step_limiter!)
end

alg_cache(alg::SSPRK932,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = SSPRK932ConstantCache()

struct SSPRK54Cache{uType,rateType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  k₃::rateType
  u₂::uType
  u₃::uType
  tmp::uType # should be u₄, but tmp is needed for callbacks
  fsalfirst::rateType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
  β10::Float64
  α20::Float64
  α21::Float64
  β21::Float64
  α30::Float64
  α32::Float64
  β32::Float64
  α40::Float64
  α43::Float64
  β43::Float64
  α52::Float64
  α53::Float64
  β53::Float64
  α54::Float64
  β54::Float64
  c1::Float64
  c2::Float64
  c3::Float64
  c4::Float64
end

u_cache(c::SSPRK54Cache) = (c.u₂,c.u₃)
du_cache(c::SSPRK54Cache) = (c.k,c.fsalfirst,c.k₃)

struct SSPRK54ConstantCache <: OrdinaryDiffEqConstantCache
  β10::Float64
  α20::Float64
  α21::Float64
  β21::Float64
  α30::Float64
  α32::Float64
  β32::Float64
  α40::Float64
  α43::Float64
  β43::Float64
  α52::Float64
  α53::Float64
  β53::Float64
  α54::Float64
  β54::Float64
  c1::Float64
  c2::Float64
  c3::Float64
  c4::Float64
end

function alg_cache(alg::SSPRK54,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  u₂ = similar(u)
  u₃ = similar(u)
  tmp = similar(u)
  k = zero(rate_prototype)
  k₃ = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  β10 = 0.391752226571890
  α20 = 0.444370493651235
  α21 = 0.555629506348765
  β21 = 0.368410593050371
  α30 = 0.620101851488403
  α32 = 0.379898148511597
  β32 = 0.251891774271694
  α40 = 0.178079954393132
  α43 = 0.821920045606868
  β43 = 0.544974750228521
  α52 = 0.517231671970585
  α53 = 0.096059710526147
  β53 = 0.063692468666290
  α54 = 0.386708617503269
  β54 = 0.226007483236906
  c1 = 0.391752226571890
  c2 = 0.586079689311540
  c3 = 0.474542363121400
  c4 = 0.935010630967653
  SSPRK54Cache(u,uprev,k,k₃,u₂,u₃,tmp,fsalfirst,alg.stage_limiter!,alg.step_limiter!,
                β10,α20,α21,β21,α30,α32,β32,α40,α43,β43,α52,α53,β53,α54,β54,c1,c2,c3,c4)
end

function alg_cache(alg::SSPRK54,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}})
  β10 = 0.391752226571890
  α20 = 0.444370493651235
  α21 = 0.555629506348765
  β21 = 0.368410593050371
  α30 = 0.620101851488403
  α32 = 0.379898148511597
  β32 = 0.251891774271694
  α40 = 0.178079954393132
  α43 = 0.821920045606868
  β43 = 0.544974750228521
  α52 = 0.517231671970585
  α53 = 0.096059710526147
  β53 = 0.063692468666290
  α54 = 0.386708617503269
  β54 = 0.226007483236906
  c1 = 0.391752226571890
  c2 = 0.586079689311540
  c3 = 0.474542363121400
  c4 = 0.935010630967653
  SSPRK54ConstantCache(β10,α20,α21,β21,α30,α32,β32,α40,α43,β43,α52,α53,β53,α54,β54,c1,c2,c3,c4)
end

struct SSPRK104Cache{uType,rateType,StageLimiter,StepLimiter} <: OrdinaryDiffEqMutableCache
  u::uType
  uprev::uType
  k::rateType
  k₄::rateType
  tmp::uType
  fsalfirst::rateType
  stage_limiter!::StageLimiter
  step_limiter!::StepLimiter
end

u_cache(c::SSPRK104Cache) = ()
du_cache(c::SSPRK104Cache) = (c.k,c.fsalfirst,c.k₄)

struct SSPRK104ConstantCache <: OrdinaryDiffEqConstantCache end

function alg_cache(alg::SSPRK104,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{true}})
  tmp = similar(u)
  k = zero(rate_prototype)
  k₄ = zero(rate_prototype)
  fsalfirst = zero(rate_prototype)
  SSPRK104Cache(u,uprev,k,k₄,tmp,fsalfirst,alg.stage_limiter!,alg.step_limiter!)
end

alg_cache(alg::SSPRK104,u,rate_prototype,uEltypeNoUnits,uBottomEltypeNoUnits,tTypeNoUnits,uprev,uprev2,f,t,dt,reltol,p,calck,::Type{Val{false}}) = SSPRK104ConstantCache()


text:
@testset "utils.jl" begin
    @testset "addlogprob!" begin
        @model function testmodel()
            global lp_before = getlogp(__varinfo__)
            @addlogprob!(42)
            global lp_after = getlogp(__varinfo__)
        end

        model = testmodel()
        varinfo = VarInfo(model)
        @test iszero(lp_before)
        @test getlogp(varinfo) == lp_after == 42
    end

    @testset "getargs_dottilde" begin
        # Some things that are not expressions.
        @test getargs_dottilde(:x) === nothing
        @test getargs_dottilde(1.0) === nothing
        @test getargs_dottilde([1.0, 2.0, 4.0]) === nothing

        # Some expressions.
        @test getargs_dottilde(:(x ~ Normal(μ, σ))) === nothing
        @test getargs_dottilde(:((.~)(x, Normal(μ, σ)))) == (:x, :(Normal(μ, σ)))
        @test getargs_dottilde(:((~).(x, Normal(μ, σ)))) == (:x, :(Normal(μ, σ)))
        @test getargs_dottilde(:(x .~ Normal(μ, σ))) == (:x, :(Normal(μ, σ)))
        @test getargs_dottilde(:(@. x ~ Normal(μ, σ))) === nothing
        @test getargs_dottilde(:(@. x ~ Normal(μ, $(Expr(:$, :(sqrt(v))))))) === nothing
        @test getargs_dottilde(:(@~ Normal.(μ, σ))) === nothing
    end

    @testset "getargs_tilde" begin
        # Some things that are not expressions.
        @test getargs_tilde(:x) === nothing
        @test getargs_tilde(1.0) === nothing
        @test getargs_tilde([1.0, 2.0, 4.0]) === nothing

        # Some expressions.
        @test getargs_tilde(:(x ~ Normal(μ, σ))) == (:x, :(Normal(μ, σ)))
        @test getargs_tilde(:((.~)(x, Normal(μ, σ)))) === nothing
        @test getargs_tilde(:((~).(x, Normal(μ, σ)))) === nothing
        @test getargs_tilde(:(@. x ~ Normal(μ, σ))) === nothing
        @test getargs_tilde(:(@. x ~ Normal(μ, $(Expr(:$, :(sqrt(v))))))) === nothing
        @test getargs_tilde(:(@~ Normal.(μ, σ))) === nothing
    end
end


text:
include("../IsolationForest.jl")

"""
    IForestDetector()

Determine the anomaly score of a sample based on their average path lengths on trees in a forest, see [1].

Parameters
----------
$(SCORE_UNSUPERVISED("IForest"))

References
----------
[1] Liu, Fei Tony, Kai Ming Ting, and Zhi-Hua Zhou. “Isolation Forest.” In 2008 Eighth IEEE International Conference on
Data Mining, 413–22. Pisa, Italy: IEEE, 2008.
"""
OD.@detector mutable struct IForestDetector <: UnsupervisedDetector
    num_trees::Int = 100
    sub_sampling_size::Int = 256
    normalize::Bool = false
end

struct IForestModel <: DetectorModel
    forest::IsolationForest.Forest
end

to_sample_set(X) = eachcol(X)
to_feature_set(X) = map(collect, eachrow(X))
make_scorer(forest, normalize::Bool) = normalize ?
    x -> IsolationForest.score_sample_against_forest_normalized(forest, x) :
    x -> IsolationForest.score_sample_against_forest(forest, x)

function OD.fit(detector::IForestDetector, X::Data; verbosity)::Fit
    num_features = size(X, 1)
    feature_values = to_feature_set(X)
    sample_values = to_sample_set(X)
    forest = IsolationForest.Forest(detector.num_trees, detector.sub_sampling_size, num_features, feature_values)
    score = make_scorer(forest, detector.normalize)
    return IForestModel(forest), score.(sample_values)
end

function OD.transform(detector::IForestDetector, model::IForestModel, X::Data)::Scores
    sample_values = to_sample_set(X)
    score = make_scorer(model.forest, detector.normalize)
    return score.(sample_values)
end


text:
#####
##### Interface for players
#####

"""
    AbstractPlayer{Game}

Abstract type for a player of `Game`.
"""
abstract type AbstractPlayer{Game} end

GameType(::AbstractPlayer{Game}) where Game = Game

"""
    think(::AbstractPlayer, game)

Return a probability distribution over actions as a `(actions, π)` pair.
"""
function think end

"""
    reset_player!(::AbstractPlayer)

Reset the internal memory of a player (e.g. the MCTS tree).
The default implementation does nothing.
"""
function reset_player!(::AbstractPlayer)
  return
end

"""
    player_temperature(::AbstractPlayer, game, turn_number)

Return the player temperature, given the number of actions that have
been played before by both players in the current game.

A default implementation is provided that always returns 1.
"""
function player_temperature(::AbstractPlayer, game, turn)
  return 1.0
end

"""
    select_move(player::AbstractPlayer, game, turn_number)

Return a single action. A default implementation is provided that samples
an action according to the distribution computed by [`think`](@ref), with a
temperature given by [`player_temperature`](@ref).
"""
function select_move(player::AbstractPlayer, game, turn_number)
  actions, π = think(player, game)
  τ = player_temperature(player, game, turn_number)
  π = apply_temperature(π, τ)
  return actions[Util.rand_categorical(π)]
end

#####
##### Random Player
#####

"""
    RandomPlayer{Game} <: AbstractPlayer{Game}

A player that picks actions uniformly at random.
"""
struct RandomPlayer{Game} <: AbstractPlayer{Game} end

function think(player::RandomPlayer, game)
  actions = GI.available_actions(game)
  n = length(actions)
  π = ones(n) ./ length(actions)
  return actions, π
end

#####
##### Epsilon-greedy player
#####

"""
    EpsilonGreedyPlayer{Game, Player} <: AbstractPlayer{Game}

A wrapper on a player that makes it choose a random move
with a fixed ``ϵ`` probability.
"""
struct EpsilonGreedyPlayer{G, P} <: AbstractPlayer{G}
  player :: P
  ϵ :: Float64
  function EpsilonGreedyPlayer(p::AbstractPlayer{G}, ϵ) where G
    return new{G, typeof(p)}(p, ϵ)
  end
end

function think(p::EpsilonGreedyPlayer, game)
  actions, π = think(p.player, game)
  n = length(actions)
  η = ones(n) ./ n
  return actions, (1 - p.ϵ) * π + p.ϵ * η
end

function reset!(p::EpsilonGreedyPlayer)
  reset!(p.player)
end

function player_temperature(p::EpsilonGreedyPlayer, game, turn)
  return player_temperature(p.player, game, turn)
end

#####
##### Player with a custom temperature
#####

"""
    PlayerWithTemperature{Game, Player} <: AbstractPlayer{Game}

A wrapper on a player that enables overwriting the temperature schedule.
"""
struct PlayerWithTemperature{G, P} <: AbstractPlayer{G}
  player :: P
  temperature :: AbstractSchedule{Float64}
  function PlayerWithTemperature(p::AbstractPlayer{G}, τ) where G
    return new{G, typeof(p)}(p, τ)
  end
end

function think(p::PlayerWithTemperature, game)
  return think(p.player, game)
end

function reset!(p::PlayerWithTemperature)
  reset!(p.player)
end

function player_temperature(p::PlayerWithTemperature, game, turn)
  return p.temperature[turn]
end

#####
##### An MCTS-based player
#####

"""
    MctsPlayer{Game, MctsEnv} <: AbstractPlayer{Game}

A player that selects actions using MCTS.

# Constructors

    MctsPlayer(mcts::MCTS.Env; τ, niters, timeout=nothing)

Construct a player from an MCTS environment. When computing each move:

- if `timeout` is provided,
  MCTS simulations are executed for `timeout` seconds by groups of `niters`
- otherwise, `niters` MCTS simulations are run

The temperature parameter `τ` can be either a real number or a
[`AbstractSchedule`](@ref).

    MctsPlayer(oracle::MCTS.Oracle, params::MctsParams; timeout=nothing)

Construct an MCTS player from an oracle and an [`MctsParams`](@ref) structure.
"""
struct MctsPlayer{G, M} <: AbstractPlayer{G}
  mcts :: M
  niters :: Int
  timeout :: Union{Float64, Nothing}
  τ :: AbstractSchedule{Float64} # Temperature
  function MctsPlayer(mcts::MCTS.Env{G}; τ, niters, timeout=nothing) where G
    @assert niters > 0
    @assert isnothing(timeout) || timeout > 0
    new{G, typeof(mcts)}(mcts, niters, timeout, τ)
  end
end

# Alternative constructor
function MctsPlayer(
    oracle::MCTS.Oracle{G}, params::MctsParams; timeout=nothing) where G
  mcts = MCTS.Env{G}(oracle,
    gamma=params.gamma,
    cpuct=params.cpuct,
    noise_ϵ=params.dirichlet_noise_ϵ,
    noise_α=params.dirichlet_noise_α,
    prior_temperature=params.prior_temperature)
  return MctsPlayer(mcts,
    niters=params.num_iters_per_turn,
    τ=params.temperature,
    timeout=timeout)
end

# MCTS with random oracle
function RandomMctsPlayer(::Type{G}, params::MctsParams) where G
  oracle = MCTS.RandomOracle{G}()
  mcts = MCTS.Env{G}(oracle,
    cpuct=params.cpuct,
    gamma=params.gamma,
    noise_ϵ=params.dirichlet_noise_ϵ,
    noise_α=params.dirichlet_noise_α)
  return MctsPlayer(mcts,
    niters=params.num_iters_per_turn,
    τ=params.temperature)
end

function think(p::MctsPlayer, game)
  if isnothing(p.timeout) # Fixed number of MCTS simulations
    MCTS.explore!(p.mcts, game, p.niters)
  else # Run simulations until timeout
    start = time()
    while time() - start < p.timeout
      MCTS.explore!(p.mcts, game, p.niters)
    end
  end
  return MCTS.policy(p.mcts, game)
end

function player_temperature(p::MctsPlayer, game, turn)
  return p.τ[turn]
end

function reset_player!(player::MctsPlayer)
  MCTS.reset!(player.mcts)
end

#####
##### Network player
#####

"""
    NetworkPlayer{Game, Net} <: AbstractPlayer{Game}

A player that uses the policy output by a neural network directly,
instead of relying on MCTS. The given neural network must be in test mode.
"""
struct NetworkPlayer{G, N} <: AbstractPlayer{G}
  network :: N
   function NetworkPlayer(nn::MCTS.Oracle{G}) where G
    return new{G, typeof(nn)}(nn)
  end
end

function think(p::NetworkPlayer, game)
  actions = GI.available_actions(game)
  state = GI.current_state(game)
  π, _ = p.network(state)
  return actions, π
end

#####
##### Merging two players into one
#####

"""
    TwoPlayers{Game} <: AbstractPlayer{Game}

If `white` and `black` are two [`AbstractPlayer`](@ref), then
`TwoPlayers(white, black)` is a player that behaves as `white` when `white`
is to play and as `black` when `black` is to play.
"""
struct TwoPlayers{G, W, B} <: AbstractPlayer{G}
  white :: W
  black :: B
  function TwoPlayers(white, black)
    G = GameType(white)
    @assert G == GameType(black)
    return new{G, typeof(white), typeof(black)}(white, black)
  end
end

flipped_colors(p::TwoPlayers) = TwoPlayers(p.black, p.white)

function think(p::TwoPlayers, game)
  if GI.white_playing(game)
    return think(p.white, game)
  else
    return think(p.black, game)
  end
end

function select_move(p::TwoPlayers, game, turn)
  if GI.white_playing(game)
    return select_move(p.white, game, turn)
  else
    return select_move(p.black, game, turn)
  end
end

function reset!(p::TwoPlayers)
  reset!(p.white)
  reset!(p.black)
end

function player_temperature(p::TwoPlayers, game, turn)
  if GI.white_playing(game)
    return player_temperature(p.white, game, turn)
  else
    return player_temperature(p.black, game, turn)
  end
end

#####
##### Players can play against each other
#####

"""
    play_game(player; flip_probability=0.) :: Trace

Simulate a game by an [`AbstractPlayer`](@ref) and return a trace.

- For two-player games, please use [`TwoPlayers`](@ref).
- If the `flip_probability` argument is set to ``p``, the board
  is _flipped_ randomly at every turn with probability ``p``,
  using [`GI.apply_random_symmetry`](@ref).
"""
function play_game(player; flip_probability=0.)
  Game = GameType(player)
  game = Game()
  trace = Trace{Game}(GI.current_state(game))
  while true
    if GI.game_terminated(game)
      return trace
    end
    if !iszero(flip_probability) && rand() < flip_probability
      game = GI.apply_random_symmetry(game)
    end
    actions, π_target = think(player, game)
    τ = player_temperature(player, game, length(trace))
    π_sample = apply_temperature(π_target, τ)
    a = actions[Util.rand_categorical(π_sample)]
    GI.play!(game, a)
    push!(trace, π_target, GI.white_reward(game), GI.current_state(game))
  end
end

#####
##### Utilities to manage inference servers
#####

function fill_and_evaluate(net, batch; batch_size, fill)
  n = length(batch)
  @assert n > 0
  if !fill
    return Network.evaluate_batch(net, batch)
  else
    nmissing = batch_size - n
    @assert nmissing >= 0
    if nmissing == 0
      return Network.evaluate_batch(net, batch)
    else
      append!(batch, [batch[1] for _ in 1:nmissing])
      return Network.evaluate_batch(net, batch)[1:n]
    end
  end
end

# Start a server that processes inference requests for TWO networks
# Expected queries must have fields `state` and `netid`.
function inference_server(
    net1::AbstractNetwork,
    net2::AbstractNetwork,
    nworkers;
    fill_batches)
  return Batchifier.launch_server(nworkers) do batch
    n = length(batch)
    mask1 = findall(b -> b.netid == 1, batch)
    mask2 = findall(b -> b.netid == 2, batch)
    @assert length(mask1) + length(mask2) == n
    state(x) = x.state
    batch1 = state.(batch[mask1])
    batch2 = state.(batch[mask2])
    if isempty(mask2) # there are only queries from net1
      return fill_and_evaluate(net1, batch1; batch_size=n, fill=fill_batches)
    elseif isempty(mask1) # there are only queries from net2
      return fill_and_evaluate(net2, batch2; batch_size=n, fill=fill_batches)
    else # both networks sent queries
      res1 = fill_and_evaluate(net1, batch1; batch_size=n, fill=fill_batches)
      res2 = fill_and_evaluate(net2, batch2; batch_size=n, fill=fill_batches)
      @assert typeof(res1) == typeof(res2)
      res = similar(res1, n)
      res[mask1] = res1
      res[mask2] = res2
      return res
    end
  end
end

# Start an inference server for one agent
function inference_server(net::AbstractNetwork, n; fill_batches)
  return Batchifier.launch_server(n) do batch
    fill_and_evaluate(net, batch; batch_size=n, fill=fill_batches)
  end
end

ret_oracle(x) = () -> x
do_nothing!() = nothing
send_done!(reqc) = () -> Batchifier.client_done!(reqc)
zipthunk(f1, f2) = () -> (f1(), f2())

# Two neural network oracles
function batchify_oracles(os::Tuple{AbstractNetwork, AbstractNetwork}, fill, n)
  reqc = inference_server(os[1], os[2], n, fill_batches=fill)
  G = GameType(os[1])
  make1() = Batchifier.BatchedOracle{G}(st -> (state=st, netid=1), reqc)
  make2() = Batchifier.BatchedOracle{G}(st -> (state=st, netid=2), reqc)
  return zipthunk(make1, make2), send_done!(reqc)
end

function batchify_oracles(os::Tuple{<:Any, AbstractNetwork}, fill, n)
  reqc = inference_server(os[2], n, fill_batches=fill)
  make2() = Batchifier.BatchedOracle{GameType(os[2])}(reqc)
  return zipthunk(ret_oracle(os[1]), make2), send_done!(reqc)
end

function batchify_oracles(os::Tuple{AbstractNetwork, <:Any}, fill, n)
  reqc = inference_server(os[1], n, fill_batches=fill)
  make1() = Batchifier.BatchedOracle{GameType(os[1])}(reqc)
  return zipthunk(make1, ret_oracle(os[2])), send_done!(reqc)
end

function batchify_oracles(os::Tuple{<:Any, <:Any}, fill, n)
  return zipthunk(ret_oracle(os[1]), ret_oracle(os[2])), do_nothing!
end

function batchify_oracles(o::AbstractNetwork, fill, n)
  reqc = inference_server(o, n, fill_batches=fill)
  make() = Batchifier.BatchedOracle{GameType(o)}(reqc)
  return make, send_done!(reqc)
end

function batchify_oracles(o::Any, fill, n)
  return ret_oracle(o), do_nothing!
end

#####
##### Distributed simulator
#####

"""
    Simulator(make_player, oracles, measure)

A distributed simulator that encapsulates the details of running simulations
across multiple threads and multiple machines.

# Arguments

    - `oracles`: the oracles used by the player. This argument can be either
      `nothing`, a single oracle or a pair of oracles.
    - `make_player`: a function that takes as an argument the `oracles` field
      above and nuild a player from it.
    - `measure(trace, colors_flipped, player)`: the function that is used to
      take measurements after each game simulation.
"""
struct Simulator{MakePlayer, Oracles, Measure}
  make_player :: MakePlayer
  oracles :: Oracles
  measure :: Measure
end

"""
    record_trace

A measurement function to be passed to a [`Simulator`](@ref) that produces
named tuples with two fields: `trace::Trace` and `colors_flipped::Bool`.
"""
record_trace(t, cf, p) = (trace=t, colors_flipped=cf)

"""
    @enum ColorPolicy ALTERNATE_COLORS BASELINE_WHITE CONTENDER_WHITE

Policy for attributing colors in a duel between a baseline and a contender.
"""
@enum ColorPolicy ALTERNATE_COLORS BASELINE_WHITE CONTENDER_WHITE

"""
    simulate(::Simulator; <keyword arguments>)

Play a series of games using a given [`Simulator`](@ref).

# Keyword Arguments

  - `num_games`: number of games to play
  - `num_workers`: number of workers tasks to spawn
  - `handler`: called every time a game is played with the simulated trace
  - `reset_every`: if set, players are reset every `reset_every` games
  - `color_policy`: either `nothing` or a [`ColorPolicy`](@ref)
  - `flip_probability=0.`: see [`play_game`](@ref)

# Return

Return a vector of objects returned by `simulator.measure`.
"""
function simulate(
    simulator::Simulator;
    num_games,
    num_workers,
    handler,
    reset_every=nothing,
    fill_batches=true,
    flip_probability=0.,
    color_policy=nothing)

  lock = ReentrantLock() # only used to surround the calls to `handler`
  spawn_oracles, done =
    batchify_oracles(simulator.oracles, fill_batches, num_workers)
  return Util.mapreduce(1:num_games, num_workers, (x, y) -> [x;y], []) do
    oracles = spawn_oracles()
    player = simulator.make_player(oracles)
    worker_sim_id = 0
    # For each worker
    function simulate_game(sim_id)
      worker_sim_id += 1
      # Switch players' colors if necessary
      if !isnothing(color_policy)
        @assert isa(player, TwoPlayers)
        colors_flipped =
          (color_policy == BASELINE_WHITE) ||
          (color_policy == ALTERNATE_COLORS && sim_id % 2 == 1)
        # "_pf" stands for "possibly flipped"
        player_pf = colors_flipped ? flipped_colors(player) : player
      else
        colors_flipped = false
        player_pf = player
      end
      # Play the game and generate a report
      trace = play_game(player_pf, flip_probability=flip_probability)
      report = simulator.measure(trace, colors_flipped, player)
      # Reset the player periodically
      if !isnothing(reset_every) && worker_sim_id % reset_every == 0
        reset_player!(player)
      end
      # Call the handler
      Base.lock(lock)
      handler(trace)
      Base.unlock(lock)
      return report
    end
    return (process=simulate_game, terminate=done)
  end
end

function compute_redundancy(states)
  unique = Set(states)
  return 1. - length(unique) / length(states)
end

# samples is a vector of named tuples with fields `trace` and `colors_flipped`
function rewards_and_redundancy(samples; gamma)
  rewards = map(samples) do s
    wr = total_reward(s.trace, gamma)
    return s.colors_flipped ? -wr : wr
  end
  states = [st for s in samples for st in s.trace.states]
  redundancy = compute_redundancy(states)
  return rewards, redundancy
end

#####
##### Simple utilities to play an interactive game
#####

"""
    Human{Game} <: AbstractPlayer{Game}

Human player that queries the standard input for actions.

Does not implement [`think`](@ref) but instead implements
[`select_move`](@ref) directly.
"""
struct Human{Game} <: AbstractPlayer{Game} end

struct Quit <: Exception end

function select_move(::Human, game, turn)
  a = nothing
  while isnothing(a) || a ∉ GI.available_actions(game)
    print("> ")
    str = readline()
    print("\n")
    isempty(str) && throw(Quit())
    a = GI.parse_action(game, str)
  end
  return a
end

"""
    interactive!(game, white, black)

Launch an interactive session for `game::AbstractGame` between players
`white` and `black`. Both players have type `AbstractPlayer` and one of them
is typically [`Human`](@ref).
"""
function interactive!(game, player)
  try
  GI.render(game)
  turn = 0
  while !GI.game_terminated(game)
    action = select_move(player, game, turn)
    GI.play!(game, action)
    GI.render(game)
    turn += 1
  end
  catch e
    isa(e, Quit) || rethrow(e)
    return
  end
end

interactive!(game, white, black) = interactive!(game, TwoPlayers(white, black))

interactive!(game::G) where G = interactive!(game, Human{G}(), Human{G}())