text:
using .FFTW
laplacianunwrap(ϕ) = laplacianunwrap!(copy(ϕ))
function laplacianunwrap!(ϕ::AbstractArray)
    FFTW.set_num_threads(Threads.nthreads())
    ϕ .+= 2π .* k(ϕ)
end

# Schofield and Zhu 2003, https://doi.org/10.1364/OL.28.001194
k(ϕw) = 1 / 2π .* ∇⁻²(∇²_nw(ϕw) - ∇²(ϕw))  # (1)

∇²(x) = -(2π)^ndims(x) / length(x) .* idct(pqterm(size(x)) .* dct(x))  # (2)

∇⁻²(x) = -length(x) / (2π)^ndims(x) .* idct(dct(x) ./ pqterm(size(x)))  # (3)

∇²_nw(ϕw) = cos.(ϕw) .* ∇²(sin.(ϕw)) .- sin.(ϕw) .* ∇²(cos.(ϕw))  # (in text)

pqterm(sz::NTuple{1}) = (1:sz[1]).^2  # 1D case
pqterm(sz::NTuple{2}) = [p^2 + q^2 for p in 1:sz[1], q in 1:sz[2]]  # 2D case
pqterm(sz::NTuple{3}) = [p^2 + q^2 + t^2 for p in 1:sz[1], q in 1:sz[2], t in 1:sz[3]]  # 3D case
pqterm(sz::NTuple{4}) = [p^2 + q^2 + t^2 + r^2 for p in 1:sz[1], q in 1:sz[2], t in 1:sz[3], r in 1:sz[4]]  # 4D case

"""
    laplacianunwrap(ϕ::AbstractArray)

Performs laplacian unwrapping on the input phase. (1D - 4D)
The phase has to be scaled to radians.
The implementation is close to the original publication: Schofield and Zhu 2003, https://doi.org/10.1364/OL.28.001194.
It is not the fastest implementation of laplacian unwrapping (doesn't use discrete laplacian).
"""
laplacianunwrap, laplacianunwrap!


text:
#############################################
###MatrixEquivalent of the trace operation###
#############################################
"""
    buildTraceMatrixRight(maxLen; typed=Float64)

Create the matrix that takes the trace on the right-most site of a DM of length maxLen
typed is the type of the data store in the matrix: for rational arythmetic, choose SymPy.Sym
"""
function buildTraceMatrixRight(maxLen; typed=Float64)
    trmat=spzeros(typed, 2^(2*maxLen), 2^(2*maxLen+2))
    for x in 1:2^maxLen
        for y=1:2^maxLen
            trmat[x+(y-1)*2^maxLen, x+(y-1)*2^(maxLen+1)]=1
            trmat[x+(y-1)*2^maxLen, (x+2^maxLen)+(y-1+2^maxLen)*2^(maxLen+1)]=1
        end
    end
    return trmat
end
"""
    buildTraceMatrixLeft(maxLen; typed=Float64)

create the matrix that takes the trace on the left-most site of a DM of length maxLen
typed is the type of the data store in the matrix: for rational arythmetic, choose SymPy.Sym
"""
function buildTraceMatrixLeft(maxLen; typed=Float64)
    trmat=spzeros(typed, 2^(2*maxLen), 2^(2*maxLen+2))
    for x in 1:2^maxLen
        for y=1:2^maxLen
            trmat[x+(y-1)*2^maxLen, 2*(x-1)+2*(y-1)*2^(maxLen+1)+1]=1
            trmat[x+(y-1)*2^maxLen, 1+2*(x-1)+(2*(y-1)+1)*2^(maxLen+1)+1]=1
        end
    end
    return trmat
end

"""
    buildTraceMatrix(maxLen)

create the matrix that takes the trace of a DM of length maxLen
"""
function buildTraceMatrix(maxLen)
    trmat=spzeros(1, 2^(2*maxLen+2))
    for x in 1:2^(maxLen+1)
        trmat[x+2^(maxLen+1)*(x-1)]=1
    end
    return trmat
end

"""
    initializeTraceMatrices(maxLen)

Build a complete table of Traces
Traces[i, j] takes (i-1) the trace on the left, and (j-1) the trace on the right
"""
function initializeTraceMatrices(maxLen)
    TrLeft=buildTraceMatrixLeft(2);
    #global const pseudoIdentity=sparse(triming.(MatJoin*invMatJoin));
    Traces=Array{typeof(TrLeft)}(maxLen, maxLen)
    Traces[1, 1]=speye(2^(2*maxLen), 2^(2*maxLen));
    for j=2:maxLen
        Traces[j, 1]=transpose(buildTraceMatrixLeft(maxLen-j+1)*transpose(Traces[j-1, 1]))
    end
    for j=1:maxLen-1
        for k=2:maxLen+1-j
            Traces[j, k]=transpose(buildTraceMatrixRight(maxLen-j-k+2)*transpose(Traces[j, k-1]))
        end
    end
    return Traces
end
##########################################
###Building and reducing density matrix###
##########################################
"""
    buildReducedDensityMatrix(rho, site, Traces)

Compute the single-site density matrix of rho on site site.
Convention site=0 is the leftmost site.
Traces is the set of Trace matrices previously computed
"""
@inline function buildReducedDensityMatrix(rho, site, Traces)
    return transpose(At_mul_B(rho, Traces[1+site, end-site]))
end

"""
    buildReducedDensityMatrix(rho, site, nb, Traces)

Compute the reduced density matrix of rho on sites site:site+nb-1.
Convention site=0 is the leftmost site.
Traces is the set of Trace matrices previously computed
"""
@inline function buildReducedDensityMatrix(rho, site, nb, Traces)    #Local densityMatrix for sites site:site+nb-1
    len=getLength(length(rho));
    if len==nb
        return rho
    end
    if nb==1
        return buildReducedDensityMatrix(rho, site, Traces)
    end
    return transpose(At_mul_B(rho, Traces[1+site, len-nb-site+1]))
end

"""
    reduceDensityMatrixLeft(rho, Traces)

Takes the left trace of rho
Traces is the set of Trace matrices previously computed
"""
@inline function reduceDensityMatrixLeft(rho, Traces)
    return transpose(At_mul_B(rho, Traces[2, 1]))
end

"""
    reduceDensityMatrixLeft!(redrho, rho, Traces)

Takes the left trace of rho to redrho
Traces is the set of Trace matrices previously computed
"""
@inline function reduceDensityMatrixLeft!(redrho, rho, Traces)    #Local densityMatrix for site site
    redrho.=transpose(At_mul_B(rho, Traces[2, 1]))
end

"""
    reduceDensityMatrixRight(rho, Traces)

Takes the right trace of rho
Traces is the set of Trace matrices previously computed
"""
@inline function reduceDensityMatrixRight(rho, Traces)
    return transpose(At_mul_B(rho, Traces[1, 2]))
end

"""
    reduceDensityMatrixRight!(redrho, rho, Traces)

Takes the right trace of rho to redrho
Traces is the set of Trace matrices previously computed
"""
@inline function reduceDensityMatrixRight!(redrho, rho, Traces)   #Local densityMatrix for site site
    redrho.=transpose(At_mul_B(rho, Traces[1, 2]))
end

####################
###Reimplementation of kron###
####################
"""
    leftFusionIdentity(rho)

Compute I2⊗rho in vector form
Note that it would correspond to kron(rho, I2) for matrices
"""
function leftFusionIdentity(rho)
        len=getLength(length(rho));
        stepp=2^len
        extrho=zeros(Complex128, 2^(2*len+2))
        offset=0; offset2=0
        for j=1:stepp
            for k=1:stepp
                @inbounds extrho[offset2+2*k-1]=0.5*rho[offset+k]
                @inbounds extrho[offset2+2*stepp+2*k]=0.5*rho[offset+k]
            end
            offset+=stepp
            offset2+=4*stepp
        end
        return extrho
end

"""
    leftFusionIdentity(extrho, rho)

Compute extrho=I2⊗rho in vector form
Note that it would correspond to kron(rho, I2) for matrices
"""
function leftFusionIdentity!(extrho, rho)
        len=getLength(length(rho));
        @simd for i =1:length(extrho)
            @inbounds extrho[i]=zero(Complex128)
        end
        stepp=2^len
        offset=0; offset2=0
        for j=1:stepp
            for k=1:stepp
                @inbounds extrho[offset2+2*k-1]=0.5*rho[offset+k]
                @inbounds extrho[offset2+2*stepp+2*k]=0.5*rho[offset+k]
            end
            offset+=stepp
            offset2+=4*stepp
        end
end

"""
    leftFusion(smallrho, rho)

Compute smallrho⊗rho in vector form
Note that it would correspond to kron(rho, smallrho) for matrices
"""
function leftFusion(smallrho, rho)
        len=getLength(length(rho));
        extrho=Array{Complex128}(2^(2*len+2))
        stepp=2^len
        offset=0; offset2=0
        for j=1:stepp
            for k=1:stepp
                @inbounds extrho[offset2+2*k-1]=smallrho[1]*rho[offset+k]
                @inbounds extrho[offset2+2*k]=smallrho[2]*rho[offset+k]
                @inbounds extrho[offset2+2*stepp+2*k-1]=smallrho[3]*rho[offset+k]
                @inbounds extrho[offset2+2*stepp+2*k]=smallrho[4]*rho[offset+k]
            end
            offset+=stepp
            offset2+=4*stepp
        end
        return extrho
end

"""
    leftFusion!(extrho, smallrho, rho)

Compute extrho=smallrho⊗rho in vector form
Note that it would correspond to kron(rho, smallrho) for matrices
"""
function leftFusion!(extrho, smallrho, rho)
        len=getLength(length(rho));
        stepp=2^len
        offset=0; offset2=0
        for j=1:stepp
            for k=1:stepp
                @inbounds extrho[offset2+2*k-1]=smallrho[1]*rho[offset+k]
                @inbounds extrho[offset2+2*k]=smallrho[2]*rho[offset+k]
                @inbounds extrho[offset2+2*stepp+2*k-1]=smallrho[3]*rho[offset+k]
                @inbounds extrho[offset2+2*stepp+2*k]=smallrho[4]*rho[offset+k]
            end
            offset+=stepp
            offset2+=4*stepp
        end
end

"""
    rightFusionIdentity(rho)

Compute rho⊗I2 in vector form
Note that it would correspond to kron(I2, rho) for matrices
"""
function rightFusionIdentity(rho)
        len=getLength(length(rho));
        stepp=2^len; temp=2^(2*len+1)+stepp
        extrho=zeros(Complex128, 2^(2*len+2))
        offset=0; offset2=0
        for j=1:stepp
            for k=1:stepp
                @inbounds extrho[offset2+k]=0.5*rho[offset+k]
                @inbounds extrho[offset2+k+temp]=0.5*rho[offset+k]
            end
            offset+=stepp
            offset2+=2*stepp
        end
        return extrho
end

"""
    rightFusionIdentity!(extrho, rho)

Compute extrho=rho⊗I2 in vector form
Note that it would correspond to kron(I2, rho) for matrices
"""
function rightFusionIdentity!(extrho, rho)
        len=getLength(length(rho));
        stepp=2^len; temp=2^(2*len+1)+stepp
        @simd for i=1:length(extrho)
            @inbounds extrho[i]=zero(Complex128)
        end
        offset=0; offset2=0
        for j=1:stepp
            for k=1:stepp
                @inbounds extrho[offset2+k]=0.5*rho[offset+k]
                @inbounds extrho[offset2+k+temp]=0.5*rho[offset+k]
            end
            offset+=stepp
            offset2+=2*stepp
        end
end

"""
    rightFusion(rho, smallrho)

Compute rho⊗smallrho in vector form
Note that it would correspond to kron(smallrho, rho) for matrices
"""
function rightFusion(rho, smallrho)
    len=getLength(length(rho));
    stepp=2^len; temp=2^(2*len+1)
    extrho=Array{Complex128}(2^(2*len+2))
    offset=0; offset2=0
    for j=1:stepp
        for k=1:stepp
            @inbounds extrho[offset2+k]=smallrho[1]*rho[offset+k]
            @inbounds extrho[offset2+k+stepp]=smallrho[2]*rho[offset+k]
            @inbounds extrho[offset2+k+temp]=smallrho[3]*rho[offset+k]
            @inbounds extrho[offset2+k+temp+stepp]=smallrho[4]*rho[offset+k]
        end
        offset+=stepp
        offset2+=2*stepp
    end
    return extrho
end

"""
    rightFusion!(extrho, rho, smallrho);

Compute extrho=rho⊗smallrho in vector form
Note that it would correspond to kron(smallrho, rho) for matrices
"""
function rightFusion!(extrho, rho, smallrho)
    len=getLength(length(rho));
    stepp=2^len; temp=2^(2*len+1)
    offset=0; offset2=0
    for j=1:stepp
        for k=1:stepp
            @inbounds extrho[offset2+k]=smallrho[1]*rho[offset+k]
            @inbounds extrho[offset2+k+stepp]=smallrho[2]*rho[offset+k]
            @inbounds extrho[offset2+k+temp]=smallrho[3]*rho[offset+k]
            @inbounds extrho[offset2+k+temp+stepp]=smallrho[4]*rho[offset+k]
        end
        offset+=stepp
        offset2+=2*stepp
    end
end

########################################################################################
###Build from an operator acting on a DM the corresponding superoperator for a vector###
########################################################################################
"""
    buildSuperOperatorLeft(operator)

Return the sparse matrix representing operator acting on the left on a DM written as a vector
"""
function buildSuperOperatorLeft(operator)
    return kron(speye(typeof(operator[1]), size(operator)[1]), operator)''
end

"""
    buildSuperOperatorRight(operator)

Return the sparse matrix representing operator acting on the right on a DM written as a vector
"""
function buildSuperOperatorRight(operator)
    return kron(transpose(operator), speye(typeof(operator[1]), size(operator)[1]))''
end

"""
    buildSuperCommutator(operator)

Return the sparse matrix representing [operator, ρ] for ρ written as a vector
"""
function buildSuperCommutator(operator)
    return (buildSuperOperatorLeft(operator)-buildSuperOperatorRight(operator))''
end

ext: jl

max_stars_repo_name: LoicHerviou/MLQuantumMarginal

max_issues_repo_name: LoicHerviou/MLQuantumMarginal

max_forks_repo_name: LoicHerviou/MLQuantumMarginal


text:
using PyCall, PyPlot

## Example 1
py"""
import numpy as np
def func(x, y):
    return x*(1-x)*np.cos(4*np.pi*x) * np.sin(4*np.pi*y**2)**2
"""

np = pyimport("numpy")
scipy = pyimport("scipy")
interpolate = pyimport("scipy.interpolate")
griddata = interpolate.griddata

s = pybuiltin(:slice)
XY = get(np.mgrid, (s(0,1,100im), s(0,1,200im)))
grid_x, grid_y = XY[1,:,:], XY[2,:,:]

points = np.random.rand(1000, 2)
values = py"func"(points[:,1], points[:,2])

grid_z0 = griddata(points, values, (grid_x, grid_y), method="nearest")
grid_z1 = griddata(points, values, (grid_x, grid_y), method="linear")
grid_z2 = griddata(points, values, (grid_x, grid_y), method="cubic")

plt.subplot(221)
plt.imshow(py"func"(grid_x, grid_y), extent=(0,1,0,1), origin="lower")
plt.plot(points[:,1], points[:,2], "k.", ms=1)
plt.title("Original")
plt.subplot(222)
plt.imshow(grid_z0, extent=(0,1,0,1), origin="lower")
plt.title("Nearest")
plt.subplot(223)
plt.imshow(grid_z1, extent=(0,1,0,1), origin="lower")
plt.title("Linear")
plt.subplot(224)
plt.imshow(grid_z2, extent=(0,1,0,1), origin="lower")
plt.title("Cubic")
plt.gcf().set_size_inches(6, 6)
plt.show()

## Example 2
py"""
def func(x, y):
    return x + y
"""

s = pybuiltin(:slice)
XY = get(np.mgrid, (s(0,1,5im),s(0,1,5im)))
grid_x, grid_y = XY[1,:,:], XY[2,:,:]
points = np.random.rand(25, 2)
values = py"func"(points[:, 1], points[:, 2])

griddata(points, values, (grid_x, grid_y), method="nearest")
griddata(points, values, (grid_x, grid_y), method="cubic")

using PyPlot
plt.plot(points[:,1], points[:,2], "k.", ms=1)

## test 1
py"""
def func(x, y, z):
    return x + y + z
"""
np = pyimport("numpy")
scipy = pyimport("scipy")
interpolate = pyimport("scipy.interpolate")
griddata = interpolate.griddata

s = pybuiltin(:slice)
XYZ = get(np.mgrid, (s(0,1,10im), s(0,1,20im), s(0,1,10im)))
grid_x, grid_y, grid_z = XYZ[1,:,:,:], XYZ[2,:,:,:], XYZ[3,:,:,:]

points = np.random.rand(1000, 3)
values = py"func"(points[:,1], points[:,2], points[:,3])

grid_z0 = griddata(points, values, (grid_x, grid_y, grid_z), method="nearest")
grid_z1 = griddata(points, values, (grid_x, grid_y, grid_z), method="linear")
grid_z2 = griddata(points, values, (grid_x, grid_y, grid_z), method="cubic")


text:
struct Params
  order::Vector{Any}
  params::IdSet{Any}
  Params() = new([], IdSet())
end

@forward Params.order Base.iterate, Base.length

function Base.push!(ps::Params, x)
  if !(x in ps.params)
    push!(ps.order, x)
    push!(ps.params, x)
  end
  return ps
end

Base.push!(ps::Params, x...) = (foreach(x -> push!(ps, x), x); ps)

Params(xs) = push!(Params(), xs...)

function Base.show(io::IO, ps::Params)
  print(io, "Params([")
  join(io, ps.order, ", ")
  print(io, "])")
end

struct Grads
  grads::IdDict{Any,Any}
end

Base.show(io::IO, ps::Grads) = println(io, "Grads(...)")

Grads() = Grads(IdDict())

@forward Grads.grads Base.setindex!, Base.haskey, Base.length, Base.iterate

Grads(ps::Params) = Grads(IdDict(tracker(p) => init_grad(data(p)) for p in ps))

Base.getindex(g::Grads, x::Tracked) = g.grads[x]

function Base.getindex(g::Grads, x)
  istracked(x) || error("Object not tracked: $x")
  g[tracker(x)]
end

accum!(g::Grads, x, Δ) = g[x] = haskey(g, x) ? g[x] .+ Δ : Δ


text:
# NOTE: Adapted from the implementation found in Julia package
# `DifferentialEquations.jl`

using Parameters

@with_kw struct Tsit5Tableau
    c₁::Float64 = 0.161
    c₂::Float64 = 0.327
    c₃::Float64 = 0.9
    c₄::Float64 = 0.9800255409045097
    c₅::Float64 = 1
    c₆::Float64 = 1
    a₂₁::Float64 = 0.161
    a₃₁::Float64 = -0.008480655492356989
    a₃₂::Float64 = 0.335480655492357
    a₄₁::Float64 = 2.8971530571054935
    a₄₂::Float64 = -6.359448489975075
    a₄₃::Float64 = 4.3622954328695815
    a₅₁::Float64 = 5.325864828439257
    a₅₂::Float64 = -11.748883564062828
    a₅₃::Float64 = 7.4955393428898365
    a₅₄::Float64 = -0.09249506636175525
    a₆₁::Float64 = 5.86145544294642
    a₆₂::Float64 = -12.92096931784711
    a₆₃::Float64 = 8.159367898576159
    a₆₄::Float64 = -0.071584973281401
    a₆₅::Float64 = -0.028269050394068383
    a₇₁::Float64 = 0.09646076681806523
    a₇₂::Float64 = 0.01
    a₇₃::Float64 = 0.4798896504144996
    a₇₄::Float64 = 1.379008574103742
    a₇₅::Float64 = -3.290069515436081
    a₇₆::Float64 = 2.324710524099774
    b̃₁::Float64 = -0.00178001105222577714
    b̃₂::Float64 = -0.0008164344596567469
    b̃₃::Float64 = 0.007880878010261995
    b̃₄::Float64 = -0.1447110071732629
    b̃₅::Float64 = 0.5823571654525552
    b̃₆::Float64 = -0.45808210592918697
    b̃₇::Float64 = 0.015151515151515152
end

function tsit5(f, t, y, dt, P, tableau)
    (@unpack c₁,c₂,c₃,c₄,c₅,c₆,a₂₁,a₃₁,a₃₂,a₄₁,a₄₂,a₄₃,a₅₁,a₅₂,a₅₃,a₅₄,a₆₁,a₆₂,
             a₆₃,a₆₄,a₆₅,a₇₁,a₇₂,a₇₃,a₇₄,a₇₅,a₇₆ = tableau)
    k1 = f(t, y, P)
    k2 = f(t + c₁*dt, y + dt*a₂₁*k1, P)
    k3 = f(t + c₂*dt, y + dt*(a₃₁*k1 + a₃₂*k2), P)
    k4 = f(t + c₃*dt, y + dt*(a₄₁*k1 + a₄₂*k2 + a₄₃*k3), P)
    k5 = f(t + c₄*dt, y + dt*(a₅₁*k1 + a₅₂*k2 + a₅₃*k3 + a₅₄*k4), P)
    k6 = f(t + c₅*dt, y + dt*(a₆₁*k1 + a₆₂*k2 + a₆₃*k3 + a₆₄*k4 + a₆₅*k5), P)
    y + dt*(a₇₁*k1 + a₇₂*k2 + a₇₃*k3 + a₇₄*k4 + a₇₅*k5 + a₇₆*k6)
end

struct Tsit5 <: ODESolverType end

#=
        FOUR components
=#
function update(::Tsit5, fs, t, A, B, C, D, dt, P, tableau)
    (@unpack c₁,c₂,c₃,c₄,c₅,c₆,a₂₁,a₃₁,a₃₂,a₄₁,a₄₂,a₄₃,a₅₁,a₅₂,a₅₃,a₅₄,a₆₁,a₆₂,
             a₆₃,a₆₄,a₆₅,a₇₁,a₇₂,a₇₃,a₇₄,a₇₅,a₇₆ = tableau)
    kA1 = update(fs[1], t, A, B, C, D, P)
    kB1 = update(fs[2], t, A, B, C, D, P)
    kC1 = update(fs[3], t, A, B, C, D, P)
    kD1 = update(fs[4], t, A, B, C, D, P)

    Ai = A + dt*a₂₁*kA1
    Bi = B + dt*a₂₁*kB1
    Ci = C + dt*a₂₁*kC1
    Di = D + dt*a₂₁*kD1

    kA2 = update(fs[1], t + c₁*dt, Ai, Bi, Ci, Di, P)
    kB2 = update(fs[2], t + c₁*dt, Ai, Bi, Ci, Di, P)
    kC2 = update(fs[3], t + c₁*dt, Ai, Bi, Ci, Di, P)
    kD2 = update(fs[4], t + c₁*dt, Ai, Bi, Ci, Di, P)

    Ai = A + dt*(a₃₁*kA1 + a₃₂*kA2)
    Bi = B + dt*(a₃₁*kB1 + a₃₂*kB2)
    Ci = C + dt*(a₃₁*kC1 + a₃₂*kC2)
    Di = D + dt*(a₃₁*kD1 + a₃₂*kD2)

    kA3 = update(fs[1], t + c₂*dt, Ai, Bi, Ci, Di, P)
    kB3 = update(fs[2], t + c₂*dt, Ai, Bi, Ci, Di, P)
    kC3 = update(fs[3], t + c₂*dt, Ai, Bi, Ci, Di, P)
    kD3 = update(fs[4], t + c₂*dt, Ai, Bi, Ci, Di, P)

    Ai = A + dt*(a₄₁*kA1 + a₄₂*kA2 + a₄₃*kA3)
    Bi = B + dt*(a₄₁*kB1 + a₄₂*kB2 + a₄₃*kB3)
    Ci = C + dt*(a₄₁*kC1 + a₄₂*kC2 + a₄₃*kC3)
    Di = D + dt*(a₄₁*kD1 + a₄₂*kD2 + a₄₃*kD3)

    kA4 = update(fs[1], t + c₃*dt, Ai, Bi, Ci, Di, P)
    kB4 = update(fs[2], t + c₃*dt, Ai, Bi, Ci, Di, P)
    kC4 = update(fs[3], t + c₃*dt, Ai, Bi, Ci, Di, P)
    kD4 = update(fs[4], t + c₃*dt, Ai, Bi, Ci, Di, P)

    Ai = A + dt*(a₅₁*kA1 + a₅₂*kA2 + a₅₃*kA3 + a₅₄*kA4)
    Bi = B + dt*(a₅₁*kB1 + a₅₂*kB2 + a₅₃*kB3 + a₅₄*kB4)
    Ci = C + dt*(a₅₁*kC1 + a₅₂*kC2 + a₅₃*kC3 + a₅₄*kC4)
    Di = D + dt*(a₅₁*kD1 + a₅₂*kD2 + a₅₃*kD3 + a₅₄*kD4)

    kA5 = update(fs[1], t + c₄*dt, Ai, Bi, Ci, Di, P)
    kB5 = update(fs[2], t + c₄*dt, Ai, Bi, Ci, Di, P)
    kC5 = update(fs[3], t + c₄*dt, Ai, Bi, Ci, Di, P)
    kD5 = update(fs[4], t + c₄*dt, Ai, Bi, Ci, Di, P)

    Ai = A + dt*(a₆₁*kA1 + a₆₂*kA2 + a₆₃*kA3 + a₆₄*kA4 + a₆₅*kA5)
    Bi = B + dt*(a₆₁*kB1 + a₆₂*kB2 + a₆₃*kB3 + a₆₄*kB4 + a₆₅*kB5)
    Ci = C + dt*(a₆₁*kC1 + a₆₂*kC2 + a₆₃*kC3 + a₆₄*kC4 + a₆₅*kC5)
    Di = D + dt*(a₆₁*kD1 + a₆₂*kD2 + a₆₃*kD3 + a₆₄*kD4 + a₆₅*kD5)

    kA6 = update(fs[1], t + c₅*dt, Ai, Bi, Ci, Di, P)
    kB6 = update(fs[2], t + c₅*dt, Ai, Bi, Ci, Di, P)
    kC6 = update(fs[3], t + c₅*dt, Ai, Bi, Ci, Di, P)
    kD6 = update(fs[4], t + c₅*dt, Ai, Bi, Ci, Di, P)

    (A + dt*(a₇₁*kA1 + a₇₂*kA2 + a₇₃*kA3 + a₇₄*kA4 + a₇₅*kA5 + a₇₆*kA6),
     B + dt*(a₇₁*kB1 + a₇₂*kB2 + a₇₃*kB3 + a₇₄*kB4 + a₇₅*kB5 + a₇₆*kB6),
     C + dt*(a₇₁*kC1 + a₇₂*kC2 + a₇₃*kC3 + a₇₄*kC4 + a₇₅*kC5 + a₇₆*kC6),
     D + dt*(a₇₁*kD1 + a₇₂*kD2 + a₇₃*kD3 + a₇₄*kD4 + a₇₅*kD5 + a₇₆*kD6))
end

#=
        THREE components
=#
function update(::Tsit5, fs, t, A, B, C, dt, P, tableau)
    (@unpack c₁,c₂,c₃,c₄,c₅,c₆,a₂₁,a₃₁,a₃₂,a₄₁,a₄₂,a₄₃,a₅₁,a₅₂,a₅₃,a₅₄,a₆₁,a₆₂,
             a₆₃,a₆₄,a₆₅,a₇₁,a₇₂,a₇₃,a₇₄,a₇₅,a₇₆ = tableau)
    kA1 = update(fs[1], t, A, B, C, P)
    kB1 = update(fs[2], t, A, B, C, P)
    kC1 = update(fs[3], t, A, B, C, P)

    Ai = A + dt*a₂₁*kA1
    Bi = B + dt*a₂₁*kB1
    Ci = C + dt*a₂₁*kC1

    kA2 = update(fs[1], t + c₁*dt, Ai, Bi, Ci, P)
    kB2 = update(fs[2], t + c₁*dt, Ai, Bi, Ci, P)
    kC2 = update(fs[3], t + c₁*dt, Ai, Bi, Ci, P)

    Ai = A + dt*(a₃₁*kA1 + a₃₂*kA2)
    Bi = B + dt*(a₃₁*kB1 + a₃₂*kB2)
    Ci = C + dt*(a₃₁*kC1 + a₃₂*kC2)

    kA3 = update(fs[1], t + c₂*dt, Ai, Bi, Ci, P)
    kB3 = update(fs[2], t + c₂*dt, Ai, Bi, Ci, P)
    kC3 = update(fs[3], t + c₂*dt, Ai, Bi, Ci, P)

    Ai = A + dt*(a₄₁*kA1 + a₄₂*kA2 + a₄₃*kA3)
    Bi = B + dt*(a₄₁*kB1 + a₄₂*kB2 + a₄₃*kB3)
    Ci = C + dt*(a₄₁*kC1 + a₄₂*kC2 + a₄₃*kC3)

    kA4 = update(fs[1], t + c₃*dt, Ai, Bi, Ci, P)
    kB4 = update(fs[2], t + c₃*dt, Ai, Bi, Ci, P)
    kC4 = update(fs[3], t + c₃*dt, Ai, Bi, Ci, P)

    Ai = A + dt*(a₅₁*kA1 + a₅₂*kA2 + a₅₃*kA3 + a₅₄*kA4)
    Bi = B + dt*(a₅₁*kB1 + a₅₂*kB2 + a₅₃*kB3 + a₅₄*kB4)
    Ci = C + dt*(a₅₁*kC1 + a₅₂*kC2 + a₅₃*kC3 + a₅₄*kC4)

    kA5 = update(fs[1], t + c₄*dt, Ai, Bi, Ci, P)
    kB5 = update(fs[2], t + c₄*dt, Ai, Bi, Ci, P)
    kC5 = update(fs[3], t + c₄*dt, Ai, Bi, Ci, P)

    Ai = A + dt*(a₆₁*kA1 + a₆₂*kA2 + a₆₃*kA3 + a₆₄*kA4 + a₆₅*kA5)
    Bi = B + dt*(a₆₁*kB1 + a₆₂*kB2 + a₆₃*kB3 + a₆₄*kB4 + a₆₅*kB5)
    Ci = C + dt*(a₆₁*kC1 + a₆₂*kC2 + a₆₃*kC3 + a₆₄*kC4 + a₆₅*kC5)

    kA6 = update(fs[1], t + c₅*dt, Ai, Bi, Ci, P)
    kB6 = update(fs[2], t + c₅*dt, Ai, Bi, Ci, P)
    kC6 = update(fs[3], t + c₅*dt, Ai, Bi, Ci, P)

    (A + dt*(a₇₁*kA1 + a₇₂*kA2 + a₇₃*kA3 + a₇₄*kA4 + a₇₅*kA5 + a₇₆*kA6),
     B + dt*(a₇₁*kB1 + a₇₂*kB2 + a₇₃*kB3 + a₇₄*kB4 + a₇₅*kB5 + a₇₆*kB6),
     C + dt*(a₇₁*kC1 + a₇₂*kC2 + a₇₃*kC3 + a₇₄*kC4 + a₇₅*kC5 + a₇₆*kC6))
end

#=
        TWO components
=#
function update(::Tsit5, fs, t, A, B, dt, P, tableau)
    (@unpack c₁,c₂,c₃,c₄,c₅,c₆,a₂₁,a₃₁,a₃₂,a₄₁,a₄₂,a₄₃,a₅₁,a₅₂,a₅₃,a₅₄,a₆₁,a₆₂,
             a₆₃,a₆₄,a₆₅,a₇₁,a₇₂,a₇₃,a₇₄,a₇₅,a₇₆ = tableau)
    kA1 = update(fs[1], t, A, B, P)
    kB1 = update(fs[2], t, A, B, P)

    Ai = A + dt*a₂₁*kA1
    Bi = B + dt*a₂₁*kB1

    kA2 = update(fs[1], t + c₁*dt, Ai, Bi, P)
    kB2 = update(fs[2], t + c₁*dt, Ai, Bi, P)

    Ai = A + dt*(a₃₁*kA1 + a₃₂*kA2)
    Bi = B + dt*(a₃₁*kB1 + a₃₂*kB2)

    kA3 = update(fs[1], t + c₂*dt, Ai, Bi, P)
    kB3 = update(fs[2], t + c₂*dt, Ai, Bi, P)

    Ai = A + dt*(a₄₁*kA1 + a₄₂*kA2 + a₄₃*kA3)
    Bi = B + dt*(a₄₁*kB1 + a₄₂*kB2 + a₄₃*kB3)

    kA4 = update(fs[1], t + c₃*dt, Ai, Bi, P)
    kB4 = update(fs[2], t + c₃*dt, Ai, Bi, P)

    Ai = A + dt*(a₅₁*kA1 + a₅₂*kA2 + a₅₃*kA3 + a₅₄*kA4)
    Bi = B + dt*(a₅₁*kB1 + a₅₂*kB2 + a₅₃*kB3 + a₅₄*kB4)

    kA5 = update(fs[1], t + c₄*dt, Ai, Bi, P)
    kB5 = update(fs[2], t + c₄*dt, Ai, Bi, P)

    Ai = A + dt*(a₆₁*kA1 + a₆₂*kA2 + a₆₃*kA3 + a₆₄*kA4 + a₆₅*kA5)
    Bi = B + dt*(a₆₁*kB1 + a₆₂*kB2 + a₆₃*kB3 + a₆₄*kB4 + a₆₅*kB5)

    kA6 = update(fs[1], t + c₅*dt, Ai, Bi, P)
    kB6 = update(fs[2], t + c₅*dt, Ai, Bi, P)

    (A + dt*(a₇₁*kA1 + a₇₂*kA2 + a₇₃*kA3 + a₇₄*kA4 + a₇₅*kA5 + a₇₆*kA6),
     B + dt*(a₇₁*kB1 + a₇₂*kB2 + a₇₃*kB3 + a₇₄*kB4 + a₇₅*kB5 + a₇₆*kB6))
end

"""
    createTableau(::Tsit5)

Tableau of coefficients for the `Tsit5` ODE solver
"""
createTableau(::Tsit5) = Tsit5Tableau()


text:
module NRRD

using Mmap, Printf
# Packages needed to return the possible range of element types
using FixedPointNumbers, Colors, ColorVectorSpace, StaticArrays, Quaternions
# Other packages
using AxisArrays, ImageAxes, Unitful, MappedArrays
using FileIO
import Libz
import FixedPointNumbers

using Colors: AbstractGray
using AxisArrays: HasAxes

string2type = Dict(
    "signed char" => Int8,
    "int8" => Int8,
    "int8_t" => Int8,
    "uchar" => UInt8,
    "unsigned char" => UInt8,
    "uint8" => UInt8,
    "uint8_t" => UInt8,
    "short" => Int16,
    "short int" => Int16,
    "signed short" => Int16,
    "signed short int" => Int16,
    "int16" => Int16,
    "int16_t" => Int16,
    "ushort" => UInt16,
    "unsigned short" => UInt16,
    "unsigned short int" => UInt16,
    "uint16" => UInt16,
    "uint16_t" => UInt16,
    "int" => Int32,
    "signed int" => Int32,
    "int32" => Int32,
    "int32_t" => Int32,
    "uint" => UInt32,
    "unsigned int" => UInt32,
    "uint32" => UInt32,
    "uint32_t" => UInt32,
    "longlong" => Int64,
    "long long" => Int64,
    "long long int" => Int64,
    "signed long long" => Int64,
    "signed long long int" => Int64,
    "int64" => Int64,
    "int64_t" => Int64,
    "ulonglong" => UInt64,
    "unsigned long long" => UInt64,
    "unsigned long long int" => UInt64,
    "uint64" => UInt64,
    "uint64_t" => UInt64,
    "float16" => Float16,
    "float" => Float32,
    "double" => Float64
)

# The opposite of string2type
type2string(::Type{Float16}) = "float16"
type2string(::Type{Float32}) = "float"
type2string(::Type{Float64}) = "double"
type2string(::Type{T}) where {T<:Integer} = lowercase(string(T.name.name))
type2string(::Type{Normed{T,f}}) where {T<:Unsigned,f} = type2string(T)
type2string(::Type{T}) where {T} = type2string(eltype(T), T)
type2string(::Type{T}, ::Type{T}) where {T} = error("type $T unrecognized")
type2string(::Type{T1}, ::Type{T2}) where {T1,T2} = type2string(T1)

space2axes = Dict(
    "right-anterior-superior" => (3,(:R,:A,:S)),
    "ras" => (3,(:R,:A,:S)),
    "left-anterior-superior" => (3,(:L,:A,:S)),
    "las" => (3,(:L,:A,:S)),
    "left-posterior-superior" => (3,(:L,:P,:S)),
    "lps" => (3,(:L,:P,:S)),
    "right-anterior-superior-time" => (4,(:R,:A,:S,:time)),
    "rast" => (4,(:R,:A,:S,:time)),
    "left-anterior-superior-time" => (4,(:L,:A,:S,:time)),
    "last" => (4,(:L,:A,:S,:time)),
    "left-posterior-superior-time" => (4,(:L,:P,:S,:time)),
    "lpst" => (4,(:L,:P,:S,:time)),
    "scanner-xyz" => (3,(:scannerx,:scannery,:scannerz)),
    "scanner-xyz-time" => (4,(:scannerx,:scannery,:scannerz,:time)),
    "3d-right-handed" => (3,(:xrcs,:yrcs,:zrcs)),
    "3d-left-handed" => (3,(:xlcs,:ylcs,:zlcs)),
    "3d-right-handed-time" => (4,(:xrcs,:yrcs,:zrcs,:time)),
    "3d-left-handed-time" => (4,(:xlcs,:ylcs,:zlcs,:time)),
)

const axes2space = Dict(
    (:R,:A,:S) => "right-anterior-superior",
    (:L,:A,:S) => "left-anterior-superior",
    (:L,:P,:S) => "left-posterior-superior",
    (:R,:A,:S,:time) => "right-anterior-superior-time",
    (:L,:A,:S,:time) => "left-anterior-superior-time",
    (:L,:P,:S,:time) => "left-posterior-superior-time",
    (:scannerx,:scannery,:scannerz) => "scanner-xyz",
    (:scannerx,:scannery,:scannerz,:time) => "scanner-xyz-time",
    (:xrcs,:yrcs,:zrcs) => "3d-right-handed",
    (:xlcs,:ylcs,:zlcs) => "3d-left-handed",
    (:xrcs,:yrcs,:zrcs,:time) => "3d-right-handed-time",
    (:xlcs,:ylcs,:zlcs,:time) => "3d-left-handed-time")

# We put these in a dict so that we don't eval untrusted
# strings. Please submit PRs to add to this list if you need
# additional unit support.
unit_string_dict = Dict("" => 1, "m" => u"m", "mm" => u"mm", "s" => u"s",
                        "um" => u"μm", "μm" => u"μm", "microns" => u"μm",
                        "pixel" => 1)

struct QString end                 # string with quotes around it: "mm"
VTuple{T} = Tuple{Vararg{T}}  # space-delimited tuple: 80 150
struct PTuple{T} end               # parenthesis-delimited tuple: (80,150)
struct StringPTuple{T} end         # string or PTuple{T}

struct IntFloat end

# This should list anything that DOESN'T parse to a string
const parse_type = Dict(
    # basic
    "dimension"=>Int,
    "block size"=>Int,
    "blocksize"=>Int,
    "min"=>Float64,
    "max"=>Float64,
    "old min"=>Float64,
    "oldmin"=>Float64,
    "old max"=>Float64,
    "oldmax"=>Float64,
    "line skip"=>Int,
    "lineskip"=>Int,
    "byte skip"=>Int,
    "byteskip"=>Int,
    # orientation
    "space dimension"=>Int,
    "space units"=>VTuple{QString},
    "space origin"=>PTuple{IntFloat},
    "space directions"=>VTuple{StringPTuple{Float64}},
    "measurement frame"=>VTuple{PTuple{Float64}},
    # per-axis
    "sizes"=>VTuple{Int},
    "spacings"=>VTuple{IntFloat},
    "thicknesses"=>VTuple{Float64},
    "axis mins"=>VTuple{IntFloat},
    "axismins"=>VTuple{IntFloat},
    "axis maxs"=>VTuple{IntFloat},
    "axismaxs"=>VTuple{IntFloat},
    "centers"=>VTuple{String},
    "centerings"=>VTuple{String},
    "labels"=>VTuple{QString},
    "units"=>VTuple{QString},
    "kinds"=>VTuple{String},
)

const fieldorder = ["content", "type", "dimension", "space", "space dimension",
                    "sizes", "spacings", "space directions", "kinds",
                    "centers", "centerings", "thickness",
                    "axis mins", "axismins", "axis maxs", "axismaxs",
                    "labels", "units",
                    "min", "max", "old min", "oldmin", "old max", "oldmax",
                    "block size", "blocksize", "endian", "encoding",
                    "space units", "space origin", "measurement frame",
                    "line skip", "lineskip", "byte skip", "byteskip",
                    "sample units", "sampleunits",
                    "data file", "datafile"]

const per_axis = [# orientation-related fields
                  "space directions",
                  # other per-axis fields
                  "sizes", "spacings", "thicknesses", "axis mins", "axismins",
                  "axis maxs", "axismaxs", "centers", "centerings", "labels",
                  "units", "kinds"]

const per_spacedim = ["space units", "space origin", "measurement frame"]

# version >= n is required if it has any fields in version_reqs[n]
const version_reqs = ([],
                      [],   # key/value tested separately,
                      ["kinds"],
                      ["thicknesses", "sample units", "space",
                       "space dimension", "space directions",
                       "space origin", "space units",
                       "data file", "datafile"],
                      ["measurement frame"],
)

function myendian()
    if ENDIAN_BOM == 0x04030201
        return "little"
    elseif ENDIAN_BOM == 0x01020304
        return "big"
    end
end

# Don't extend FileIO.load
# Set mode to "r+" if you want to be able to modify values in the
# image and have them update in the disk file
function load(f::File{format"NRRD"}, args...; mode="r", mmap=:auto)
    open(f, mode) do io
        skipmagic(io)
        load(io, args...; mode=mode, mmap=mmap)
    end
end

function load(io::Stream{format"NRRD"}, Tuser::Type=Any; mode="r", mmap=:auto)
    # Assemble all the information about the array we're about to
    # read: element type, size, and the "meaning" of axes
    version, header, keyvals, comments = parse_header(io)
    Traw, need_bswap = raw_eltype(header)
    szraw = (header["sizes"]...,)  # "sizes" may change in outer_eltype!, grab it now
    T, nd, perm = outer_eltype!(header, Traw, Tuser)
    axs = get_axes(header, nd)
    sz = get_size(axs)

    # Read the data
    iodata = find_datafile(io, header; mode=mode)
    compressed = in(header["encoding"], ("gzip", "gz"))
    if compressed
        iodata = Libz.ZlibInflateInputStream(iodata)
    end

    can_mmap = header["encoding"] == "raw"

    if mmap == true && (!can_mmap)
        error("Cannot use memory-mapped for reading a non-raw or bswapped file")
    end

    # Use memory-mapping for large files
    do_mmap = can_mmap && (prod(szraw) > 10^8) && (mmap == :auto)
    do_mmap |= can_mmap && (mmap == true)

    if !compressed
        szraw, sz = checked_size(Traw, szraw, sz, iodata)
    end

    if do_mmap
        # Recent Julia versions are picky about alignment
        pos = position(iodata)
        if pos % sizeof(Traw) != 0
            szcor = (szraw[1]*sizeof(Traw), szraw[2:end]...)
            A = reinterpret(Traw, Mmap.mmap(iodata, Array{UInt8,length(szraw)}, szcor, pos;
                          grow=false))
        else
            A = Mmap.mmap(iodata, Array{Traw,length(szraw)}, szraw, pos;
                          grow=false)
        end
        if need_bswap
            f = mode == "r+" ? (bswap, bswap) : bswap
            A = mappedarray(f, A)
        end
    elseif header["encoding"] == "raw" || in(header["encoding"], ("gzip", "gz"))
        A = read!(iodata, Array{Traw}(undef, szraw...))
        if need_bswap
            A = [bswap(a) for a in A]
        end
    else
        error("\"", header["encoding"], "\" encoding not supported.")
    end

    if perm == ()
        if T != eltype(A)
            A = need_bswap ? A = mappedarray(x->T(x), A) : reshape(reinterpret(T, A), sz)
        end
    else
        A = permuteddimsview(A, perm)
        if T<:Color
            A = colorview(T, A)
        end
    end

    isa(axs, Dims) ? A : AxisArray(A, axs)
end

function save(f::File{format"NRRD"}, img::AbstractArray; kwargs...)
    open(f, "w") do io
        write(io, magic(format"NRRD"))
        save(io, img; kwargs...)
    end
end

function save(io::Stream{format"NRRD"}, img::AbstractArray{T}; props::Dict = Dict{String,Any}(), keyvals=nothing, comments=nothing, kwargs...) where T
    axs = axisinfo(img)
    header = headerinfo(T, axs; kwargs...)
    header_eltype!(header, T)
    # copy fields from props to override those in header
    for (k, v) in props
        header[k] = v
    end
    v = version(header, !(keyvals==nothing || isempty(keyvals)))
    write_header(io, v, header, keyvals, comments)
    datafilename = get(props, "datafile", "")
    if isempty(datafilename)
        datafilename = get(props, "data file", "")
    end
    if isempty(datafilename)
        nrrd_write(io, img)
    else
        println(io.io, "data file: ", datafilename)
        if !get(props, "headeronly", false)
            open(datafilename, "w") do file
                nrrd_write(file, img)
            end
        end
    end
end

axisinfo(img) = axisinfo(HasAxes(img), img)
axisinfo(::HasAxes{true}, img) = AxisArrays.axes(img)
axisinfo(::HasAxes, img) = size(img)

### Interpreting header settings

"""
    arraytype!(header, version) -> T, axs, perm, need_bswap

Analyze the `header` dictionary to extract the element-type `T`, size
or axes information `axs`, the permutation `perm` (if any) that julia
should use for "wrapping" the read data, and a boolean `need_bswap`
indicating whether the data need to be byte-swapped (to account for
differences in endianness). `T` includes any color information (in
which case a dimension of the array will be "consumed"). `axs` will be
a Dims-tuple in simple cases, or an `Axes` tuple (from AxisArrays.jl)
if dimensions are labeled or have their spatial information
(pixelspacing, spacedirections, etc) specified. `perm` is the
permutation needed to move the color data to the first dimension, or
an empty tuple if no permutation is required.

This function may modify the `header` dictionary (the reason for the !
in the name), so make a copy first if necessary.
"""
function arraytype!(header, version)
    Traw, need_bswap = raw_eltype(header)
    T, nd, perm = outer_eltype!(header, Traw)
    axs = get_axes(header, nd)
    T, axs, perm, need_bswap
end

"""
    arraytype(filename)

Parse NRRD header and calls `arraytype!(header, version)`. See
`arraytype!` for information about the return values.
"""
function arraytype(filename)
    version, header, keyvals, comments = parse_header(filename)
    arraytype!(header, version)
end

function headerinfo(T, axs)
    header = Dict{String,Any}()
    Traw = raw_eltype(T)
    header["type"] = type2string(Traw)
    header["endian"] = myendian()
    header["encoding"] = "raw"
    if T <: Gray
        val = gray(oneunit(T))
        val = isa(val, FixedPoint) ? reinterpret(val) : val
        header["sample units"] = string("gray ", val)
    elseif T <: Union{RGB,RGBA}
        val = red(oneunit(T))
        val = isa(val, FixedPoint) ? reinterpret(val) : val
        valfmt = isa(val, Integer) ? 'd' : 'f'
        if T <: RGB
            colstr = "rgb"
            valfmtstr = "(%$valfmt,%$valfmt,%$valfmt)"
            vals = (val,val,val)
        else
            colstr = "rgba"
            valfmtstr = "(%$valfmt,%$valfmt,%$valfmt,%$valfmt)"
            vals = (val,val,val,val)
        end
        fmtstr = "%s $valfmtstr"
        x = (colstr, vals...)
        header["sample units"] = @eval @sprintf($fmtstr, $(x...))
    end
    # Do the axes information
    header["dimension"] = length(axs)
    if isa(axs, Base.Indices)
        axs = map(length, axs)
    end
    specifyorientation = false
    if isa(axs, Dims)
        header["sizes"] = [axs...]
    else
        # axs is an Axis-tuple
        header["sizes"] = [map(length, axs)...]
        axnames = map(ax->axisnames(ax)[1], axs)
        isspace = map(s->!startswith(string(s), "time"), axnames)
        if haskey(axes2space, axnames)
            header["space"] = axes2space[axnames]
            specifyorientation = true
        end
        header["kinds"] = [isspc ? "domain" : "time" for isspc in isspace]
        if !all(isdefaultname, axnames)
            header["labels"] = [string(s) for s in axnames]
        end
        rng = map(ax->axisvalues(ax)[1], axs)
        stepval = map(step, rng)
        unitstr = map(x->isa(x, Quantity) ? string(unit(x)) : "", stepval)
        spacing = map(x->isa(x, Quantity) ? ustrip(x) : x,        stepval)
        if !all(x->x=="", unitstr)
            header["units"] = [unitstr...]
        end
        if !all(x->x==1, spacing)
            if specifyorientation
                header["space directions"] = [ntuple(d->d==i ? spacing[i] : 0, length(spacing)) for i = 1:length(spacing)]
            else
                header["spacings"] = [Float64.(spacing)...]
            end
        end
        origin = map(x->isa(x, Quantity) ? ustrip(x) : x, map(first, rng))
        if specifyorientation && any(x->x!=0, origin)
            header["space origin"] = [origin[[isspace...]]...]
        end
    end
    # Adjust the axes for color
    if T <: Colorant && !(T <: AbstractGray)
        header["dimension"] = length(axs)+1
        pushfirst!(header["sizes"], length(T))
        if haskey(header, "spacings")
            pushfirst!(header["spacings"], NaN)
        end
        if haskey(header, "labels")
            pushfirst!(header["labels"], lowercase(string(T.name.name)))
        end
        if haskey(header, "units")
            pushfirst!(header["units"], "")
        end
        if !haskey(header, "kinds")
            header["kinds"] = ["domain" for d = 1:length(axs)]
        end
        if T <: RGB
            colkind = "RGB-color"
        elseif T <: HSV
            colkind = "HSV-color"
        elseif T <: XYZ
            colkind = "XYZ-color"
        elseif T <: RGBA
            colkind = "RGBA-color"
        else
            colkind = string(length(T), "-color")
        end
        pushfirst!(header["kinds"], colkind)
    end
    header
end

function version(header, has_keyvalue::Bool=false)
    for n = length(version_reqs):-1:1
        vr = version_reqs[n]
        for f in vr
            if haskey(header, f)
                return n
            end
        end
    end
    has_keyvalue ? 2 : 1
end

function isdefaultname(s::AbstractString)
    startswith(s, "dim_") || startswith(s, "space") ||
        startswith(s, "time") || startswith(s, "domain")
end
isdefaultname(s::Symbol) = isdefaultname(string(s))

"""
    raw_eltype(header) -> Traw, need_bswap
    raw_eltype(::Type{T}) -> Traw

Get the "basic" element type of the data, e.g., `UInt16` or
`Float32`.

This function does not try to determine whether the image is color
(`Traw` does not contain any color information), nor does it try to
interpret `Traw` as a `Normed` type.

See also: outer_eltype!, fixedtype.
"""
function raw_eltype(header)
    Traw = string2type[header["type"]]
    need_bswap = haskey(header, "endian") && header["endian"] != myendian() && sizeof(Traw) > 1
    Traw, need_bswap
end

raw_eltype(::Type{C}) where {C<:Colorant} = raw_eltype(eltype(C))
raw_eltype(::Type{T}) where {T<:FixedPoint}   = FixedPointNumbers.rawtype(T)
raw_eltype(::Type{T}) where {T} = T

"""
    fixedtype(Traw, header) -> Tu

Attempt to interpret type `Traw` in terms of FixedPoint numbers. The
interpretation depends on whether `header` has a "sample units" field
of the form:

    sample units: <colorspace> <whitepoint>

There must be a space between `colorspace` and `whitepoint`, and
`colorspace` must be one of "gray", "rgb", "rgba", "hsv", "xyz", or
their uppercase variants.  (Other than "gray", all of these are
supported "kinds" values. "gray" does not typically correspond to an
axis, which is why it isn't encoded in "kinds".) The presence of any
of these words indicates that the data represent an image rather than
some other kind of array.

Any other choices are ignored, as "sample units" can also be an
arbitrary string.

# Examples:

    # conventional uint8 grayscale
    sample units: gray 255

    # a 14-bit grayscale camera (numbers can be represented in hex format)
    sample units: gray 0x3fff

    # RGB encoded with float or double
    sample units: rgb (1.0,1.0,1.0)

    # RGB encoded with float or double, but using the scaling of uint8
    sample units: rgb (255.0,255.0,255.0)

    # conventional XYZ
    sample units: xyz (95.047,100.000,108.883)

    # HSV, hue measured in degrees
    sample units: hsv (360, 0, 1)

    # HSV, hue normalized to [0, 1]
    sample units: hsv (1, 0, 1)

If `Traw` cannot be interpreted as `Normed`, `Tu = Traw`.
"""
function fixedtype(::Type{Traw}, header) where Traw<:Unsigned
    # Note that "max" is not useful in this context.
    # See https://sourceforge.net/p/teem/bugs/14/
    if haskey(header, "sample units")
        su = header["sample units"]
    elseif haskey(header, "sampleunits")
        su = header["sampleunits"]
    else
        return Traw
    end
    idxsplit = findfirst(isequal(' '), su)
    idxsplit === nothing && return Traw
    cm, rest = lowercase(su[1:idxsplit-1]), strip(su[idxsplit+1:end])
    if cm ∈ ("gray", "rgb", "rgba", "xyz", "hsv")
        # It looks like a whitepoint definition
        if cm == "gray"
            val = numberparse(rest)
            return Gray{fixedtype_max(Traw, val)}
        elseif cm ∈ ("rgb", "rgba")
            s = nrrd_parse(PTuple{String}, rest)
            val = map(numberparse, s)
            if all(v->v==val[1], val)
                return fixedtype_max(Traw, val[1])
            end
        end
    end
    Traw
end
fixedtype(::Type{Traw}, header) where {Traw} = Traw

function fixedtype_max(::Type{Traw}, mx) where Traw<:Unsigned
    fmx = log2(mx+1)
    if round(fmx) == fmx
        return Normed{Traw,round(Int,fmx)}
    end
    Traw
end

header_eltype!(header, ::Type) = header
function header_eltype!(header, ::Type{T}) where T<:FixedPoint
    header["sample units"] = string("gray ", reinterpret(one(T)))
    header
end
function header_eltype!(header, ::Type{C}) where C<:Colorant
    _header_eltype!(header, C, eltype(C))
    header
end
function _header_eltype!(header, ::Type{C}, ::Type{T}) where {C<:AbstractGray,T<:FixedPoint}
    header["sample units"] = string("gray ", reinterpret(one(T)))
end
function _header_eltype!(header, ::Type{C}, ::Type{T}) where {C<:AbstractGray,T}
    header["sample units"] = string("gray ", one(T))
end
function _header_eltype!(header, ::Type{C}, ::Type{T}) where {C<:AbstractRGB,T<:FixedPoint}
    o = reinterpret(one(T))
    header["sample units"] = "rgb ($o,$o,$o)"
end
function _header_eltype!(header, ::Type{C}, ::Type{T}) where {C<:AbstractRGB,T}
    o = one(T)
    header["sample units"] = "rgb ($o,$o,$o)"
end
function _header_eltype!(header, ::Type{C}, ::Type) where C<:XYZ
    header["sample units"] = "xyz (95.047,100.000,108.883)"
end
function _header_eltype!(header, ::Type{C}, ::Type) where C<:HSV
    header["sample units"] = "hsv (360, 0, 1)"
end

"""
    colorant_eltype(C, T) -> Tc

Return a valid "inner" element type `Tc` for colorant type `C`. When
`T` != `Tc`, values must be "converted" before they can be interpreted
as type `C`.
"""
colorant_eltype(::Type{C}, ::Type{T}) where {C<:Colorant, T<:AbstractFloat} = C{T}
colorant_eltype(::Type{C}, ::Type{T}) where {C<:Colorant, T} = C{Float32}

"""
    UnknownColor{T,N}

An unknown Color. This type gets returned when one of the "kind"
settings is "3-color" or "4-color".
"""
struct UnknownColor{T,N} <: Color{T,N}
    col::NTuple{N,T}
end

"""
    outer_eltype!(header, Traw) -> T, nd, perm

Extract the julia array `eltype` `T`, the number of dimensions `nd`
**excluding** color/complex/vector/matrix element data, and any
permutation needed to put the eltype dimension first. Any dimensions
in the header corresponding to color (or if "kind" is set to one of
the vector types) will be "consumed" upon exit. `Traw` is the
element type as determined by `raw_eltype`.

See also: raw_eltype.
"""
function outer_eltype!(header, Traw, Tuser=Any)
    nd = header["dimension"]
    sz = header["sizes"]
    length(sz) == nd || error("parsing of sizes: $(header["sizes"]) is inconsistent with $nd dimensions")
    perm = ()
    T = Traw
    if !(T <: Tuser)
        if Tuser <: Union{AbstractRGB,ColorTypes.TransparentRGB}
            if T <: Unsigned
                T = Normed{T,8*sizeof(T)}
            end
            T = ccolor(Tuser, base_colorant_type(Tuser){T})
            length(T) == sz[1] || error("first dimension of size $(sz[1]), expected $(length(T))")
            # We've just handled the first dimension, so delete the per-axis data for it
            nd -= 1
            for fn in per_axis
                if haskey(header, fn)
                    deleteat!(header[fn], 1)
                end
            end
        else
            T = Tuser
        end
    end
    if haskey(header, "kinds")
        kinds = header["kinds"]
        length(kinds) == nd || error("parsing of kinds: $(header["kinds"]) is inconsistent with $nd dimensions")
        for i = 1:nd
            k = kinds[i]
            if k == "RGB-color"
                chksize(sz[i], 3)
                Tu = fixedtype(Traw, header)
                Tu = Tu == UInt8 ? N0f8 : (Tu == UInt16 ? N0f16 : Tu)
                T = RGB{Tu}
            elseif k == "HSV-color"
                chksize(sz[i], 3)
                T = colorant_eltype(HSV, fixedtype(Traw, header))
            elseif k == "XYZ-color"
                chksize(sz[i], 3)
                T = colorant_eltype(XYZ, fixedtype(Traw, header))
            elseif k == "RGBA-color"
                chksize(sz[i], 4)
                Tu = fixedtype(Traw, header)
                T = RGBA{Tu}
            elseif k == "3-color"
                chksize(sz[i], 3)
                T = UnknownColor{Traw,3}
            elseif k == "4-color"
                chksize(sz[i], 4)
                T = UnknownColor{Traw,4}
            elseif k == "complex"
                chksize(sz[i], 2)
                T = Complex{Traw}
            elseif k == "quaternion"
                checksize(sz[i], 2)
                T = Quaternion{Traw}
            elseif k == "2-vector"
                chksize(sz[i], 2)
                T = SVector{2,Traw}
            elseif k == "3-vector" || k == "3-gradient" || k == "3-normal"
                chksize(sz[i], 3)
                T = SVector{3,Traw}
            elseif k == "4-vector"
                chksize(sz[i], 4)
                T = SVector{4,Traw}
            elseif k == "2D-matrix"
                chksize(sz[i], 4)
                T = SMatrix{2,2,Traw}
            elseif k == "3D-matrix"
                chksize(sz[i], 9)
                T = SMatrix{3,3,Traw}
            end
            if T != Traw
                # We've handled this dimension, adjust the dimensionality
                if i > 1
                    perm = (i, setdiff(1:nd, i)...)
                end
                nd -= 1
                for fn in per_axis
                    if haskey(header, fn)
                        deleteat!(header[fn], i)
                    end
                end
                break
            end
        end
    end
    if T == Traw && nd == header["dimension"]
        T = fixedtype(Traw, header)
    end
    T, nd, perm
end

function get_axes(header, nd)
    # Validate the per-axis fields...
    for f in per_axis
        if haskey(header, f)
            length(header[f]) == nd || error("expected $nd (remaining) dimensions in field $f, got $(header[f])")
        end
    end
    # ...and the orientation fields
    if haskey(header, "space")
        sd = space2axes[lowercase(header["space"])][1]
    else
        sd = get(header, "space dimension", nd)
    end
    for f in per_spacedim
        if haskey(header, f)
            length(header[f]) == sd || error("expected $sd dimensions in field $f, got $(header[f])")
        end
    end
    ## Test whether we have fields that require an AxisArray
    axes_fields = setdiff(per_axis, ["sizes", "kinds", "centers", "centerings"])
    need_axes = haskey(header, "space")
    for f in axes_fields
        need_axes |= haskey(header, f)
    end
    # also check for a time axis
    istime = falses(nd)
    isspace = falses(nd)
    if haskey(header, "kinds")
        kinds = header["kinds"]
        istime = map(x->x=="time", kinds)
        isspace = map(x->x=="space" || x=="domain", kinds)
        need_axes |= any(istime)
    elseif haskey(header, "space dimension")
        length(isspace) == sd || error("confused, have \"space dimension\" $sd but got $nd dimensions: $header")
        fill!(isspace, true)
    elseif haskey(header, "space")
        nd == sd || error("confused, have \"space\" but can't tell which axes are spatial: $header")
        fill!(isspace, true)
    end
    sz = header["sizes"]
    if !need_axes
        return (sz...,)
    end
    ## Axis names
    axnames = [string("dim_", d) for d = 1:nd]
    if haskey(header, "labels")
        axnames = header["labels"]
        if any(istime)
            for idx in findall(istime)
                labelt = axnames[idx]
                if !istimeaxis(Axis{Symbol(labelt)})
                    @warn("label $labelt is not defined as a time axis, define it with `@traitimpl TimeAxis{Axis{:$labelt}}` (see ImageAxes for more information)")
                end
            end
        end
    elseif haskey(header, "space")
        spcnames = map(string, space2axes[lowercase(header["space"])][2])
        copy_space!(axnames, spcnames, isspace)
    elseif haskey(header, "kinds")
        # Give axes default names based on their kind: space1, space2, etc.
        kindcounter = Dict{String,Int}()
        for (i,k) in enumerate(header["kinds"])
            if !haskey(kindcounter, k)
                kindcounter[k] = 0
            end
            n = kindcounter[k]+1
            kindcounter[k] = n
            # append numeric label if more than 1
            axnames[i] = n > 1 ? string(k, '_', n) : k
        end
        # ...but if there was more than 1, go back and append "1" at
        # the end of the first one
        for (k, n) in kindcounter
            if n > 1
                idx = findall(x->x==k, axnames)
                axnames[idx] .= string(k, '_', 1)
            end
        end
    end
    ## Ranges
    # Can a "space directions" field be encoded as a "spacings" field?
    if haskey(header, "space directions")
        can_convert, spc = spacings(header["space directions"])
        if all(isnan, spc)
            # they were all "none"
            delete!(header, "space directions")
            can_convert = false
        end
        if can_convert
            !haskey(header, "spacings") || error("cannot have both \"space directions\" and \"spacings\"")
            header["spacings"] = spc
            delete!(header, "space directions")
        end
    end
    if haskey(header, "axis mins") || haskey(header, "axismins")
        axmin = haskey(header, "axis mins") ? header["axis mins"] : header["axismins"]
        axmax = haskey(header, "axis maxs") ? header["axis maxs"] : header["axismaxs"]
        rng = Any[startstoplen(axmin[d], axmax[d], sz[d]) for d = 1:nd]
    else
        startval = Any[zeros(Int, nd)...]
        stepval = Any[ones(Int, nd)...]
        if haskey(header, "space origin")
            so = header["space origin"]
            copy_space!(startval, header["space origin"], isspace)
        end
        if haskey(header, "spacings")
            copy_not_nan!(stepval, header["spacings"])
        end
        rng = Any[startsteplen(startval[d], stepval[d], sz[d]) for d = 1:nd]
    end
    # These are more laborious than they should be because of issues
    # with rounding and ranges in julia-0.5.
    if haskey(header, "units")
        ua = [unit_string_dict[x] for x in header["units"]]
        rng = Any[startstepstop(first(r)*u, step(r)*u, last(r)*u) for (r,u) in zip(rng,ua)]
    elseif haskey(header, "space units")
        us = [get(unit_string_dict, x, 1) for x in header["space units"]]
        ua = fill!(Array{Any}(undef, nd), 1)
        copy_space!(ua, us, isspace)
        for d = 1:nd
            r, u = rng[d], ua[d]
            newr = startstepstop(first(r)*u, step(r)*u, last(r)*u)
            rng[d] = newr
        end
    end
    l = map(safelength, rng)
    l == sz || error("expect axis lengths to match array sizes, got $((rng...,)) of length $l and $sz")
    # Return the axes
    ntuple(d->Axis{Symbol(axnames[d])}(rng[d]), nd)
end

get_size(sz::Dims) = sz
get_size(axs) = map(safelength, axs)

### Parsing

"""
    parse_header(io) -> version, header, keyvals, comments

Parse the NRRD header (the top of the .nrrd or the separate .nhdr
file). `io` should be positioned just after the initial "NRRD" in the
file. This reads up to and including the first blank line of the file,
so at the end (if this is a properly-formatted NRRD file) `io` is
positioned at the first byte of the data (if present).

Outputs:
- `version` is a 4-character string, e.g., "0002", giving the NRRD version of the header.
- `header` is a `Dict{String,Any}` of `field=>setting` pairs (the settings are parsed, they are not necessarily strings)
- `keyvals` is a `Dict{String,String}` containing `key=>value` pairs (NRRD0002 or higher, lines like key:=value; many NRRD files do not contain any of these)
- `comments` is an array containing lines of the header that began with `#` (but with the `#` and leading whitespace stripped out)

See also: write_header.
"""
function parse_header(io)
    version = ascii(String(read!(io, Vector{UInt8}(undef, 4))))
    skipchars(isspace, io)
    header = Dict{String, Any}()
    keyvals = Dict{String, String}()
    comments = String[]
    # Read until we encounter a blank line, which is the separator
    # between the header and data
    line = strip(readline(io))
    while !isempty(line)
        if line[1] != '#'
            idx = findfirst(isequal(':'), line)
            idx === nothing && error("no colon found in $line")
            key, value = line[1:idx-1], line[idx+1:end]
            if !isempty(value) && value[1] == '='
                # This is a NRRD key/value pair, insert into keyvals
                keyvals[key] = value[2:end]
            else
                lkey = lowercase(key)
                T = get(parse_type, lkey, String)
                header[lkey] = nrrd_parse(T, strip(value))
            end
        else
            cmt = strip(lstrip(line, ['#', ' ']))
            if !isempty(cmt)
                push!(comments, cmt)
            end
        end
        line = strip(readline(io))
    end
    version, header, keyvals, comments
end

parse_header(s::Stream{format"NRRD"}) = parse_header(stream(s))

function parse_header(filename::AbstractString)
    f = File{format"NRRD"}(filename)
    open(f) do io
        skipmagic(io)
        parse_header(io)
    end
end

"""
    write_header(io, version, header, [keyvals, [comments]])

Write an NRRD header, as the top of the .nrrd or the separate .nhdr
file. `io` should be positioned just after the initial "NRRD" in the
file. This writes the header and a blank line, so that at the end `io`
is positioned at the first byte of the data (if present).

Note that if you're writing a header for a "detached" data file
(separate .nhdr and .raw files), `header` should contain a "data file"
or "datafile" field storing the name of the .raw file.

Inputs:
- `version` is a 4-character string, e.g., "0002", giving the NRRD version of the header, or a integer corresponding to a recognized NRRD header version number
- `header` is a `Dict{String,Any}` of `field=>setting` pairs (as returned by `parse_header`)
- `keyvals` is a `Dict{String,String}` containing `key=>value` pairs (NRRD0002 or higher, lines like key:=value; many NRRD files do not contain any of these)
- `comments` is an array containing lines of the header that began with `#` (but with the `#` and leading whitespace stripped out)

See also: parse_header.
"""
function write_header(io::IO, version, header, keyvals=nothing, comments=nothing)
    writeversionstr(io, version)

    for fn in fieldorder
        if haskey(header, fn)
            print(io, fn, ": ")
            T = get(parse_type, fn, String)
            nrrd_format(io, T, header[fn])
            println(io, "")
        end
    end
    if keyvals != nothing
        for (k,v) in keyvals
            println(io, k, ":=", v)
        end
    end
    if comments != nothing
        for c in comments
            println(io, "# ", c)
        end
    end
    println(io)
    nothing
end
write_header(s::Stream{format"NRRD"}, args...) = write_header(stream(s), args...)

nrrd_parse(::Type{T}, str) where {T} = parse(T, str)  # fallback

function nrrd_parse(::Type{IntFloat}, str)
    x = parse(Float64, str)
    x == round(x) ? Int(x) : x
end

nrrd_parse(::Type{String}, str) = str

function nrrd_parse(::Type{QString}, str)
    str[1] == '"' && str[end] == '"' || error("$str does not start and end with double-quote")
    str[2:end-1]
end

function nrrd_parse(::Type{VTuple{T}}, s::AbstractString) where T
    ss = split(s)  # r"[ ,;]")
    v = Vector{alloctype(T)}(undef, length(ss))
    for i = 1:length(ss)
        v[i] = nrrd_parse(T, ss[i])
    end
    return v
end

function nrrd_parse(::Type{VTuple{QString}}, s::AbstractString)
    str = nrrd_parse(QString, s)  # to strip the first and last ""
    split(str, r"\" +\"")
end

function nrrd_parse(::Type{PTuple{T}}, s::AbstractString) where T
    s[1] == '(' && s[end] == ')' || error("$s should begin and end parentheses")
    str = s[2:end-1]
    ss = split(str, ',', keepempty=false)
    v = Vector{alloctype(T)}(undef, length(ss))
    for i = 1:length(ss)
        v[i] = nrrd_parse(T, ss[i])
    end
    return v
end

function nrrd_parse(::Type{StringPTuple{T}}, s::AbstractString) where T
    if s[1] == '(' && s[end] == ')'
        return nrrd_parse(PTuple{T}, s)
    end
    s
end

writeversionstr(io, v::Integer) = @printf(io, "%04d\n", v)
function writeversionstr(io, v::AbstractString)
    length(v) == 4 || error("version string must have 4 characters, got $v")
    println(io, v)
end

nrrd_format(io, ::Type{T}, x) where {T} = print(io, x)

function nrrd_format(io, ::Type{QString}, str)
    print(io, '"', str, '"')
end

function nrrd_format(io, ::Type{VTuple{T}}, container) where T
    isfirst = true
    for x in container
        if isfirst
            isfirst = false
        else
            print(io, ' ')
        end
        nrrd_format(io, T, x)
    end
    nothing
end

function nrrd_format(io, ::Type{PTuple{T}}, container) where T
    print(io, '(')
    isfirst = true
    for x in container
        if isfirst
            isfirst = false
        else
            print(io, ',')
        end
        nrrd_format(io, T, x)
    end
    print(io, ')')
end

nrrd_format(io, ::Type{StringPTuple{T}}, s::AbstractString) where {T} = nrrd_format(io, String, s)
nrrd_format(io, ::Type{StringPTuple{T}}, container) where {T} = nrrd_format(io, PTuple{T}, container)

alloctype(::Type{T}) where {T} = T
alloctype(::Type{StringPTuple{T}}) where {T} = Any # Union{String,Vector{T}}
alloctype(::Type{IntFloat}) = Union{Int,Float64}

### Utilities

function chksize(sz, sztarget)
    sz == sztarget || error("dimension size should be $sztarget, got $sz")
    nothing
end

function startstoplen(min::Integer, max::Integer, len::Integer)
    if max-min == len-1
        return Int(min):Int(max)
    end
    stepval = (max-min)/(len-1)
    if stepval ≈ round(stepval)
        return Int(min):round(Int,stepval):Int(max)
    end
    startstopsteplen(min, max, stepval, len)
end

function startstoplen(min, max, len)
    stepval = (max-min)/(len-1)
    startstopsteplen(min, max, stepval, len)
end

function startsteplen(min, stepval, len)
    max = min+stepval*(len-1)
    startstoplen(min, max, len)
end

function startstepstop(min, stepval, max)
    n = (max-min)/stepval + 1
    len = floor(Int, n+10*eps(n))
    startstopsteplen(min, max, stepval, len)
end

function startstopsteplen(min::T, max::T, stepval::T, len) where T<:Integer
    stepval == 1 ? (min:max) : (min:stepval:max)
end
function startstopsteplen(min::T, max::T, stepval::T, len) where T<:AbstractFloat
    imin, imax, istepval = round(Int, min), round(Int, max), round(Int, stepval)
    if min ≈ imin && max ≈ imax && stepval ≈ istepval
        return startstopsteplen(imin, imax, istepval, len)
    end
    r = range(min, stop=max, length=len)
    abs(step(r)/stepval - 1) < 0.01 || error("wanted length $len from $min to $max, got $r of length $(length(r))")
    r
end
function startstopsteplen(min::T, max::T, stepval::T, len) where T
    r = range(min, stop=max, length=len)
    abs(step(r)/stepval - 1) < 0.01 || error("wanted length $len from $min to $max, got $r of length $(length(r))")
    r
end
startstopsteplen(min, max, stepval, len) = startstopsteplen(promote(min, max, stepval)..., len)

function copy_space!(dest, src, isspace)
    length(src) == sum(isspace) || throw(DimensionMismatch("src length $(length(src)) disagrees with isspace $isspace"))
    length(dest) == length(isspace) || throw(DimensionMismatch("dest length $(length(dest)) does not agree with isspace of length $(length(isspace))"))
    d = 1
    for ds = 1:length(src)
        while !isspace[d]
            d += 1
        end
        dest[d] = src[ds]
        d += 1
    end
    dest
end

function copy_not_nan!(dest, src)
    nd = length(dest)
    length(src) == nd || throw(DimensionMismatch("dest length $nd does not match source length $(length(src))"))
    for d = 1:nd
        s = src[d]
        if !isa(s, AbstractFloat) || !isnan(s)
            dest[d] = s
        end
    end
    dest
end

function spacings(spacedirections)
    spacings = zeros(length(spacedirections))
    can_convert = true
    d = 1
    for (i, v) in enumerate(spacedirections)
        if v == "none"
            spacings[i] = NaN
        else
            n = 0
            for (j,x) in enumerate(v)
                if x != 0
                    if j == d   # make certain they are in order
                        spacings[i] = x
                        d += 1
                    else
                        can_convert = false
                    end
                    n += 1
                end
            end
            can_convert &= n < 2
        end
    end
    can_convert, spacings
end

### File-related functions

# Adjust the array size if the file is not big enough for reading to succeed
function checked_size(Traw, szraw, sz, iodata)
    cpos = position(iodata)
    datalen = div(filesize(stat(iodata)) - cpos, sizeof(Traw))
    if datalen < prod(szraw)
        if szraw != sz
            # If we've dropped a dimension due to "absorbing"
            # color, etc, don't try to figure out how to correct
            # it, just punt to the user
            error("data length $datalen too small for size $sz array of $T")
        end
        # If the data are smaller than the header suggests, read in as
        # much "complete" data as are available
        sznew = [szraw...]
        strds = [1;cumprod(sznew)]
        k = length(sznew)
        sznew[k] = div(datalen, strds[k])
        while sznew[k] == 0 && k > 1
            pop!(sznew)
            k -= 1
            sznew[k] = div(datalen, strds[k])
        end
        tsznew = (sznew...,)
        @warn("header indicates an array size $szraw, but the file size is consistent with at most $tsznew")
        szraw = sz = tsznew
    end
    szraw, sz
end

function stream2name(s::IO)
    name = s.name
    if !startswith(name, "<file ")
        error("io name ", name, " doesn't fit expected pattern")
    end
    name[7:end-1]
end

function find_datafile(iodata::IOStream, header; mode="r", path=nothing)
    if haskey(header, "data file") || haskey(header, "datafile")
        # Separate header and data files
        fdata = haskey(header, "data file") ? header["data file"] : header["datafile"]
        if path == nothing
            # Check current directory first
            if isfile(fdata)
                iodata = open(fdata)
            else
                path = dirname(stream2name(iodata))
                iodata = open(joinpath(path, fdata), mode)
            end
        else
            iodata = open(joinpath(path, fdata), mode)
        end
    end
    iodata
end

function find_datafile(s::Stream{format"NRRD"}, header; mode="r")
    find_datafile(stream(s), header; mode=mode)
end

nrrd_write(io, A::AxisArray) = nrrd_write(io, A.data)
nrrd_write(io, A::AbstractArray) = nrrd_write_elty(io, A)
nrrd_write_elty(io, A::AbstractArray{C}) where {C<:Colorant} = nrrd_write_elty(io, channelview(A))
nrrd_write_elty(io, A::AbstractArray{T}) where {T<:FixedPoint} = nrrd_write_elty(io, rawview(A))
nrrd_write_elty(io, A::AbstractArray) = write(io, A)

function numberparse(str)
    lstr = lowercase(str)
    startswith(lstr, "0x") && return parse(UInt, lstr)   # hex
    val = parse(Float64, lstr)
end

# Avoids roundoff error in Base's implementation of length
function safelength(r::StepRange)
    n = ((last(r)-first(r))/step(r)) + 1
    floor(Int, n+10*eps(n))
end
safelength(r) = length(r)

end


text:
"""
Reads pcap data from an array of bytes.
"""
mutable struct PcapBufferReader <: PcapReader
    data::Vector{UInt8}
    raw_header::Vector{UInt8}
    header::PcapHeader
    offset::Int64
    mark::Int64
    usec_mul::Int64
    bswapped::Bool

    @doc """
        PcapBufferReader(data::Vector{UInt8})

    Create reader over `data`. Will read and process pcap header,
    and yield records through `read(::PcapBufferReader)`.
    """
    function PcapBufferReader(data::Vector{UInt8})
        length(data) < sizeof(PcapHeader) && throw(EOFError())
        rh = data[1:sizeof(PcapHeader)]
        h = unsafe_load(Ptr{PcapHeader}(pointer(data)))
        h, bswapped, nanotime = process_header(h)
        new(data, rh, h, sizeof(h), -1, nanotime ? 1 : 1000, bswapped)
    end
end

"""
    PcapBufferReader(path::AbstractString)

Memory map file in `path` and create PcapBufferReader over its content.
"""
function PcapBufferReader(path::AbstractString)
    io = open(path)
    data = Mmap.mmap(io)
    PcapBufferReader(data)
end

function Base.close(x::PcapBufferReader)
    x.data = UInt8[]
    x.offset = 0
    nothing
end

Base.length(x::PcapBufferReader) = length(x.data)
Base.position(x::PcapBufferReader) = x.offset; nothing
Base.seek(x::PcapBufferReader, pos) = x.offset = pos; nothing

function Base.mark(x::PcapBufferReader)
    x.mark = x.offset
    x.mark
end

function Base.unmark(x::PcapBufferReader)
    if x.mark >= 0
        x.mark = -1
        true
    else
        false
    end
end

Base.ismarked(x::PcapBufferReader) = x.mark >= 0

function Base.reset(x::PcapBufferReader)
    !ismarked(x) && error("PcapBufferReader not marked")
    x.offset = x.mark
    x.mark = -1
    x.offset
end

Base.eof(x::PcapBufferReader) = (length(x) - x.offset) < sizeof(RecordHeader)

"""
    read(x::PcapBufferReader) -> PcapRecord

Read one record from pcap data.
Throws `EOFError` if no more data available.
"""
@inline function Base.read(x::PcapBufferReader)
    eof(x) && throw(EOFError())
    record_offset = x.offset
    p = pointer(x.data) + record_offset
    GC.@preserve x begin
        h = unsafe_load(Ptr{RecordHeader}(p))
    end
    if x.bswapped
        h = bswap(h)
    end
    t1 = (h.ts_sec + x.header.thiszone) * 1_000_000_000
    t2 = Int64(h.ts_usec) * x.usec_mul
    t = UnixTime(Dates.UTInstant(Nanosecond(t1 + t2)))
    x.offset += sizeof(RecordHeader) + h.incl_len
    x.offset > length(x) && error("Insufficient data in pcap record")
    PcapRecord(h, t, x.data, record_offset)
end


text:
export Shift, Reduce, RightArc, LeftArc, Break
export cost, is_valid, transition, valid_moves

using ..DependencyParsing

module Moves
  const SHIFT = 'S'
  const REDUCE = 'D'
  const RIGHT = 'R'
  const LEFT = 'L'
  const BREAK = 'B'
end

struct Shift <: Move end;
struct Reduce <: Move end;
struct RightArc <: Move end;
struct LeftArc <: Move end;
struct Break <: Move end;

#=
  Cost: push_cost
=#
function cost(state::GoldState, move::Shift, label::String)
  cost = 0
  if head_state(state, state.config.buffer[begin]) == HEAD_IN_STACK &&
    state.gold_tree.nodes[state.config.buffer[begin]].head_id != first(state.config.stack)
    cost += 1
  end

  cost += dependents_in_stack_count(state, state.config.buffer[begin])

  cost
end

#=
  Validity:
  * If stack is empty
  * At least two words in sentence
  * Word has not been shifted before
=#
function is_valid(config::Configuration, move::Shift)
  if stack_depth(config) == 0
    return true
  elseif buffer_length(config) < 2
    return false
  elseif is_unshiftable(config, config.buffer[begin])
    return false
  else
    return true
  end
end

# Move the first word of the buffer onto the stack and mark it as "shifted"
function transition(config::Configuration, label::String, move::Shift)
  config_push!(config)
end

#=
Cost:
  * If B[0] is the start of a sentence, cost is 0
  * Arcs between stack and buffer
  * If arc has no head, we're saving arcs between S[0] and S[1:], so decrement
      cost by those arcs.
=#
function cost(state::GoldState, move::Reduce, label::String)
  cost = 0

  cost += dependents_in_buffer_count(state, first(state.config.stack))
  if head_state(state, first(state.config.stack)) == HEAD_IN_BUFFER
    cost += 1
  end

  if !has_head(state.config.tree, first(state.config.stack))
    if head_state(state, first(state.config.stack)) == HEAD_IN_STACK
      cost -= 1
    end

    cost -= dependents_in_stack_count(state, first(state.config.stack))
  end

  cost
end

#=
  Validity:
    * Stack not empty
=#
function is_valid(config::Configuration, move::Reduce)
  stack_depth(config) == 0 && return false

  true
end

# Pop from the stack. If it has no head and the stack isn't empty, place it back on the buffer.
function transition(config::Configuration, label::String, move::Reduce)
  if has_head(config.tree, first(config.stack)) || stack_depth(config) == 1
    config_pop!(config)
  else
    unshift!(config)
  end
end

# push_cost + (not shifted[b0] and Arc(B[1:], B[0])) - Arc(S[0], B[0], label): NOT VALID COMMENT
function cost(state::GoldState, move::RightArc, label::String)
  cost = 0
  stack = state.config.stack |> collect
  s0 = stack[begin]
  b0 = state.config.buffer[begin]
  b0_gold_head = state.gold_tree.nodes[b0].head_id

  cost += dependents_in_stack_count(state, b0)
  if head_state(state, b0) == HEAD_IN_STACK
    cost += 1
  end

  if b0_gold_head == s0
    cost -= 1
    if !label_correct(state, b0, label)
      cost += 1
    end
  elseif head_state(state, b0) == HEAD_IN_BUFFER && !state.config.unshiftable[b0]
    cost += 1
  end

  cost
end
#=
  Validity:
  * len(S) >= 1
  * len(B) >= 1
  * not is_sent_start(B[0]) ?
=#
function is_valid(config::Configuration, move::RightArc)
  if stack_depth(config) < 1
    return false
  elseif buffer_length(config) < 1
    return false
  else
    return true
  end

end

# Add an arc from S[0] to B[0]. Push B[0].
function transition(config::Configuration, label::String, move::RightArc)
  dependent = config.buffer[begin]
  head = first(config.stack)
  add_arc(config, head, dependent, label)
  config_push!(config)
end

# pop_cost - Arc(B[0], S[0], label) + (Arc(S[1], S[0]) if H(S[0]) else Arcs(S, S[0])): NOT VALID COMMENT
function cost(state::GoldState, move::LeftArc, label::String)
  cost = 0
  stack = state.config.stack |> collect
  s0 = stack[begin]
  b0 = state.config.buffer[begin]
  s0_gold_head = s0 == 0 ? EMPTY_NODE : state.gold_tree.nodes[s0].head_id

  cost += dependents_in_buffer_count(state, s0)
  if head_state(state, s0) == HEAD_IN_BUFFER
    cost += 1
  end

  if has_head(state.config.tree, s0)
    if stack_depth(state.config) >= 2 && s0_gold_head == stack[begin + 1]
      cost += 1
    end
  else
    if head_state(state, s0) == HEAD_IN_STACK
      cost += 1
    end
    cost += dependents_in_stack_count(state, s0)
  end

  if s0_gold_head == b0
    cost -= 1
    if !label_correct(state, s0, label)
      cost += 1
    end
  end

  cost
end

#=
Validity:
  * len(S) >= 1
  * len(B) >= 1
  * not is_sent_start(B[0]) ?
=#
function is_valid(config::Configuration, move::LeftArc)
  if stack_depth(config) < 1
    return false
  elseif buffer_length(config) < 1
    return false
  else
    return true
  end
end

# Add an arc between B[0] and S[0], replacing the previous head of S[0] ifone is set. Pop S[0] from the stack.
function transition(config::Configuration, label::String, move::LeftArc)
  head = config.buffer[begin]
  dependent = config_pop!(config)
  has_head(config.tree, dependent) && del_arc(config.tree, dependent)
  add_arc(config, head, dependent, label)
  set_reshiftable(config, dependent)
end

# Not sure that i need this Move

function cost(config::Configuration, move::Break, system::ArcEagerSystem)
end

#=
Validity:
  * len(buffer) >= 2
  * B[1] == B[0] + 1
  * not is_sent_start(B[1])
  * not cannot_sent_start(B[1])
=#
function is_valid(config::Configuration, move::Break)
  if buffer_length(config) < 2
    return false
  elseif config.buffer[begin] != config.buffer[begin + 1]
    return false
  else
    return true
  end
end

#Mark the second word of the buffer as the start of a sentence.
function transition(config::Configuration, label::String, move::Break)
end

function valid_moves(config::Configuration)
  moves = Move[Shift(), Reduce(), RightArc(), LeftArc()]

  filter(move -> is_valid(config, move), moves)
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
using LinearAlgebra
using Plots
using Plots.PlotMeasures
using LaTeXStrings
using ColorSchemes

export  eigval_img,
        eigval_img_2d,
        berry_img,
        nernst_img,
        T_img,
        μ_img,
        Tμ_img

"""
    eigval_img(d,M,Ham!)

    plot eigen values. You specify the plot area by `k`
"""

function eigval_img(k :: AbstractArray ,M,Ham!)
    plt = Array{Float64,3}(undef,M,length(k),length(k))
    @inbounds Threads.@threads for x in eachindex(k)
        H1 = zeros(ComplexF64,M,M)
        @inbounds for y in eachindex(k)
            Ham!(k[x],k[y],H1)
            plt[:,x,y] = eigvals!(H1)
        end
    end
    p1 = plot3d(fmt=:png)
    for i in 1:M
        plot3d!(k,k, plt[i,:,:] , st = :surface)
    end
    p2 = plot(k,plt[:,:,length(k)÷2+1]', xlabel = L"k_x")
    p3 = plot(k,plt[:,length(k)÷2+1,:]', xlabel = L"k_y")
    plot(
        p1,p2,p3,layout = @layout([a [b ; c]]),size = (900,500)
    )
end

"""
    eigval_img(d,M,Ham!)

    plot eigen values. You specify the span of `k` with d
"""

function eigval_img(d :: Int,M,Ham!)
    k = -1:1/(1<<d):1
    plt = Array{Float64,3}(undef,M,length(k),length(k))
    @inbounds Threads.@threads for x in eachindex(k)
        H1 = zeros(ComplexF64,M,M)
        @inbounds for y in eachindex(k)
            Ham!(k[x],k[y],H1)
            plt[:,x,y] = eigvals!(H1)
        end
    end
    p1 = plot3d(fmt=:png)
    for i in 1:M
        plot3d!(k,k, plt[i,:,:] , st = :surface)
    end
    p2 = plot(k,plt[:,:,1<<d+1]', xlabel = L"k_x")
    p3 = plot(k,plt[:,1<<d+1,:]', xlabel = L"k_y")
    plot(
        p1,p2,p3,layout = @layout([a [b ; c]]),size = (900,500)
    )
end

function eigval_img_2d(k :: AbstractArray ,M,Ham!)
    plt = Array{Float64,3}(undef,M,length(k),length(k))
    @inbounds Threads.@threads for x in eachindex(k)
        H1 = zeros(ComplexF64,M,M)
        @inbounds for y in eachindex(k)
            Ham!(k[x],k[y],H1)
            plt[:,x,y] = eigvals!(H1)
        end
    end

    plot(
        [heatmap(k,k,plt[n,:,:]') for n in 1:M]...,size = (900,500)
    )
end

function eigval_img_2d(d :: Int ,M,Ham!)
    k = -1:1/(1<<d):1
    plt = Array{Float64,3}(undef,M,length(k),length(k))
    @inbounds Threads.@threads for x in eachindex(k)
        H1 = zeros(ComplexF64,M,M)
        @inbounds for y in eachindex(k)
            Ham!(k[x],k[y],H1)
            plt[:,x,y] = eigvals!(H1)
        end
    end

    plot(
        [heatmap(k,k,plt[n,:,:]') for n in 1:M]...,size = (900,500)
    )
end

function berry_img(d :: Int,M,Ham!)
    k = -1:1/(1<<d):1
    plt = Array{Float64,3}(undef,M,length(k),length(k))
    @inbounds Threads.@threads for x in eachindex(k)
        H1, H2, H3 = zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M)
        @inbounds for y in eachindex(k)
             plt[:,x,y]=Berry(k[x],k[y],Ham!,H1,H2,H3)
        end
    end

    @. plt /= 1<<(2d)

    plot(
        [ heatmap(k,k,plt[i,:,:]', title = "band $(i)" , c = cgrad(:bwr)) for i in 1:M]...,
        xlabel = L"k_x" , ylabel = L"k_y",titlefontsize=8
    )
end

function berry_img(k :: AbstractArray ,M,Ham!)
    plt = Array{Float64,3}(undef,M,length(k),length(k))
    dμ = step(k)^2
    @inbounds Threads.@threads for x in eachindex(k)
        H1, H2, H3 = zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M)
        @inbounds for y in eachindex(k)
             plt[:,x,y]=Berry(k[x],k[y],Ham!,H1,H2,H3)
        end
    end

    plt *= dμ
    plot(
        [ heatmap(k,k,plt[i,:,:]', title = "band $(i) " , c = cgrad(:bwr)) for i in 1:M]...,
        xlabel = L"k_x" , ylabel = L"k_y",titlefontsize=8
    )
end
"""
    nernst_img(d :: Int,T,μ,M,Ham!)
    nernst_img(k :: AbstractArray ,T,μ,M,Ham!)
"""

function nernst_img(d :: Int,T,μ,M,Ham!)
    k = -1:1/(1<<d):1
    plt = Array{Float64}(undef,length(k),length(k))
    @inbounds Threads.@threads for x in eachindex(k)
        H1, H2, H3 = zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M)
        @inbounds for y in eachindex(k)
             plt[x,y]=nernst(k[x],k[y],T,μ,Ham!,H1,H2,H3)
        end
    end
    @. plt /=1<<(2d) # dxdy
    println(sum(plt))
    heatmap(k,k,plt[:,:]',
        xlabel=L"k_x",
        ylabel=L"k_y"
    )
end

function nernst_img(k :: AbstractArray ,T,μ,M,Ham!)
    plt = Array{Float64}(undef,length(k),length(k))
    dμ = step(k)^2
    @inbounds Threads.@threads for x in eachindex(k)
        H1, H2, H3 = zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M)
        @inbounds for y in eachindex(k)
             plt[x,y]=nernst(k[x],k[y],T,μ,Ham!,H1,H2,H3)
        end
    end
    @. plt *= dμ
    println(sum(plt))

    heatmap(k,k,plt[:,:]',
        xlabel=L"k_x",
        ylabel=L"k_y"
    )
end

"""
    T_img(μ,M,Ham! ; rangeT = 50:50:500 , k = -1:1/(1<<4):1)
"""

function T_img(μ,M,Ham! ; rangeT = 50:50:500 , k = -1:1/(1<<4):1)
    plt = zeros(length(rangeT))
    dμ = step(k)^2
    @inbounds for mu in μ
        @. plt = 0.0
        @inbounds Threads.@threads for T in eachindex(rangeT)
            H1, H2, H3 = zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M)
            @inbounds for x in eachindex(k)
                @inbounds for y in eachindex(k)
                    plt[T] += nernst(k[x],k[y],rangeT[T],mu,Ham!,H1,H2,H3)
                end
            end
        end
        plot!(rangeT,plt.*dμ,label= L"\mu = %$(mu)",xlabel = "logT",ylabel="α",marker=true,ms=4,msw=0)
    end
    plot!()
end

"""
    μ_img(T,M,Ham!;rangeμ = -5:0.5:5 ,k = -1:1/(1<<4):1)
"""
function μ_img(T,M,Ham!;rangeμ = -5:0.5:5 ,k = -1:1/(1<<4):1)
    plt = zeros(length(rangeμ))
    dμ = step(k)^2
    @inbounds Threads.@threads for μ in eachindex(rangeμ)
        H1, H2, H3 = zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M)
        @inbounds for x in eachindex(k)
            @inbounds for y in eachindex(k)
                plt[μ] += nernst(k[x],k[y],T,rangeμ[μ],Ham!,H1,H2,H3)
            end
        end
    end
    plot!(rangeμ,plt.*dμ,label= L"T = %$(T)",xlabel = "μ",ylabel="α",marker=true,ms=4,msw=0)
end

"""
    Tμ_img(M,Ham!; rangeT = 100:100:1000 , rangeμ = -5:0.5:5 , k = -1:1/(1<<4):1)
"""

function Tμ_img(M,Ham!; rangeT = 100:100:1000 , rangeμ = -5:0.5:5 , k = -1:1/(1<<4):1)
    plt = zeros(length(rangeT),length(rangeμ))
    dμ = step(k)^2
    @inbounds Threads.@threads for T in eachindex(rangeT)
        H1, H2, H3 = zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M) , zeros(ComplexF64,M,M)
        @inbounds for μ in eachindex(rangeμ)
            @inbounds for x in eachindex(k)
                @inbounds for y in eachindex(k)
                    plt[T,μ] += nernst(k[x],k[y],rangeT[T],rangeμ[μ],Ham!,H1,H2,H3)
                end
            end
        end
    end
    @. plt *= dμ
    heatmap(rangeT,rangeμ,plt',
        c=cgrad(:bwr),
        xlabel = "T",ylabel="μ"
    )
end


text:
export zero_grad!
"""
    zero_grad!(var)

clear gradient storage in the whole comput-graph.
"""
function zero_grad! end

function zero_grad!(x)
    for each in parameters(x)
        zero_grad!(each)
    end
    x
end

function zero_grad!(x::Union{CachedNode, Node})
    for each in args(x)
        zero_grad!(each)
    end

    for each in kwargs(x)
        zero_grad!(each)
    end

    zero_grad!(value(x))
    x
end

function zero_grad!(x::Variable)
    fill!(x.grad, zero(eltype(x.grad)))
    x
end

export register_parameters, parameters

"""
    register_parameters(x::OperatorType) -> iterator

Returns the iterator of all parameters in the instance `x` of `OperatorType`.
Note, here, `OperatorType` does not need to be subtype of [`Operator`](@ref).
"""
function register_parameters end

# TODO: move this to YAAD
register_parameters(op::T) where T = error("$T's parameters are not registered")
register_parameters(op::Function) = ()
register_parameters(op::Operator) = ()
# just an eye candy
parameters(op) = register_parameters(op)


text:
################################################
#
# Module to test Modia/src/symbolic/Tearing.jl
#
# Author: Martin Otter, DLR-SR
#
################################################
module TestTearing

using Base.Test

include(joinpath("..","..","src","symbolic","Tearing.jl"))

@testset "\nTest Tearing" begin

@testset "... Test equation system 1" begin

   G = Any[[1,2,3,4,5], [1,2,3,4,8], [4,1,2], [3,2,1], [2,1,7]]
   es = [2,3,5,1]
   vs = [1,3,7,8]
   td = TraverseDAG(G,8)

   (eSolved,vSolved,eResidue,vTear) = tearEquations!(td,G,es,vs)

   @test vTear    == [3,8]
   @test eSolved  == [2,5]
   @test vSolved  == [1,7]
   @test eResidue == [3,1]
end

@testset "... Test equation system 2" begin

   #=
     eSolved  = [2,5]
     vSolved  = [1,7]
     eResidue = [3,1]
     vTear    = [3,8]
     ignored variables: 2,4,5,6

     known  : 2,3,4,5,6,8
     unknown: 1,7
     [1,2,3,4,8] is solved for 1
     [2,1,7]     is solved for 7
   =#

   #=
      unknowns: x1dd, x2dd, x3dd, x4d, x5
       1: 0 = x1dd - x2dd
      2: 0 = x1dd + x2dd - x3dd + x6ddd
      3: 0 = x1d  + x3dd - x4d
      4: 0 = 2*x1dd + x2dd + x3dd + x4d + x6ddd
      5: 0 = 3*x1dd + 2*x2dd + x5

      Variables
      1: x1dd
      2: x2dd
      3: x3dd
      4: x4d
      5: x5
      6: x1d
      7: x6ddd

      vTear    = [2]
      eSolved  = [1,2,3,5]
      vSolved  = [1,3,4,5]
      eResidue = [4]
   =#

   G = Any[[1,2],
           [1,2,3,7],
           [6,3,4],
           [1,2,3,4,7],
           [1,2,5]]
   es = [1,2,3,4,5]
   vs = [1,2,3,4,5]

   td2 = TraverseDAG(G,7)
   (eSolved,vSolved,eResidue,vTear) = tearEquations!(td2,G,es,vs)

   @test vTear    == [2]
   @test eSolved  == [1,2,3,5]
   @test vSolved  == [1,3,4,5]
   @test eResidue == [4]
end

@testset "... Test equation system 3" begin

   G = Any[ [1,2,4,7],
            [6,3,4],
            [2,5],
            [4,5],
            [1,2,3,7]]
   es = [1,2,3,4,5]
   vs = [1,2,3,4,5]
   td3 = TraverseDAG(G,7)
   (eSolved,vSolved,eResidue,vTear) = tearEquations!(td3,G,es,vs)

   @test vTear    == [5]
   @test eSolved  == [4,3,1,2]
   @test vSolved  == [4,2,1,3]
   @test eResidue == [5]
end

@testset "... Test equation system 4" begin
   #=
      Same as test2, but indices changed

      unknowns: z1, z2, z4, z6, z7
      1: 0 = f1(z2,z6,z7)
       2: 0 = f2(z1,z7)
      3: 0 = f3(z1,z2,z3,z7)
      4: 0 = f4(z6,z3,z4)
      5: 0 = f5(z1,z2,z3,z4,z7)
      6: 0 = f6(z1,z2,z5)
      7: 0 = f7(z4,z5,z6)

      Variables
      1: z1    (x1dd)
      2: z2    (x6ddd)
      3: z3    (x3dd)
      4: z4    (x4d)
      5: z5    (x5)
      6: z6    (x1d)
      7: z7    (x2dd)
   =#

   G = Any[[2,6,7],
           [1,2],
           [1,2,3,7],
           [6,3,4],
           [1,2,3,4,7],
           [1,2,5],
           [4,5,6]]
   es = [2,3,4,5,6]
   vs = [1,3,4,5,7]

   td4 = TraverseDAG(G,7)
   (eSolved,vSolved,eResidue,vTear) = tearEquations!(td4,G,es,vs)

   @test vTear    == [7]
   @test eSolved  == [2,3,4,6]
   @test vSolved  == [1,3,4,5]
   @test eResidue == [5]
end

@testset "... Test equation system 5" begin
   #=
     Check "marked" flag

      unknowns: z1, z2, z3, z4, z5, z6
      1: 0 = f1(z1,z6)
       2: 0 = f2(z2,z1)
      3: 0 = f3(z3,z2)
      4: 0 = f4(z4,z3)
      5: 0 = f5(z5,z4)
      6: 0 = f6(z6,z1)
   =#

   G = Any[[1,6],
           [2,1],
           [3,2],
           [4,3],
           [5,4],
           [6,5]]
   es = [1,2,3,4,5,6]
   vs = [1,2,3,4,5,6]

   td5 = TraverseDAG(G,6)
   (eSolved,vSolved,eResidue,vTear) = tearEquations!(td5,G,es,vs)

   @test vTear    == [6]
   @test eSolved  == [1,2,3,4,5]
   @test vSolved  == [1,2,3,4,5]
   @test eResidue == [6]
end

@testset "... Test equation system 6" begin
   #=
     Check predefined equations

     phi1 = phi2
     w1 = der(phi1)
     w2 = der(phi2)
     der(phi1) = der(phi2)

      unknowns: 1:w1, 2:w2, 3:der(phi1), 4:der(phi2)
      1: 0 = f1(w1,der(phi1))
      2: 0 = f2(w2,der(phi2))
      3: 0 = f3(der(phi1), der(phi2))
   =#

   G = Any[[1,3],
           [2,4],
           [3,4]]
   es = [1,2,3]
   vs = [1,2,3,4]

   td = TraverseDAG(G,4)
   (eSolved,vSolved,eResidue,vTear) = tearEquations!(td,G,es,vs;eSolvedFixed=[3],vSolvedFixed=[3],vTearFixed=[4])

   @test vTear    == [4]
   @test eSolved  == [3,1,2]
   @test vSolved  == [3,1,2]
   @test eResidue == Int[]
end

@testset "... Test equation systems 7a and 7b" begin
   neq = 10
   G = fill(fill(0,0),neq)
   for i in 1:neq
      G[i] = i==1 ? [1,neq] : [i,i-1]
   end
   #G = [[i, (i==1 ? neq : i-1)] for i = 1:neq]
   es = collect(1:neq)
   vs = es
   td5 = TraverseDAG(G,neq)
   (eSolved,vSolved,eResidue,vTear) = tearEquations!(td5,G,es,vs)
   @test vTear == [neq]

   for i in 1:neq
      G[i] = i==neq ? [1,neq] : [neq-i+1,neq-i]
   end
   #println("G = ")
   #display(G)
   #G = [[i, (i==1 ? neq : i-1)] for i = neq:-1:1]
   es = collect(1:neq)
   vs = es
   td5 = TraverseDAG(G,neq)
   (eSolved,vSolved,eResidue,vTear) = tearEquations!(td5,G,es,vs)
   @test vTear == [1]
end

end  # facts

end  # module


text:
module Reflection

export
    # Classes.jl
    findclass, isarray, isinterface, superclass, ClassDescriptor,
    # Modifiers.jl
    ModifiersDescriptor,
    # Methods.jl
    classmethods, classdeclaredmethods,
    isstatic, ispublic, MethodDescriptor,
    # Constructors.jl
    classconstructors, ConstructorDescriptor,
    # Fields.jl
    classfields, FieldDescriptor

include("Classes.jl")
include("Modifiers.jl")
include("Methods.jl")
include("Constructors.jl")
include("Fields.jl")

using .Classes
using .Modifiers
using .Methods
using .Constructors
using .Fields

end


text:
using Pkg
Pkg.activate("..")
using PlayAroundPackage
using Documenter

makedocs(;
    modules=[PlayAroundPackage],
    authors="David Yang",
    repo="https://github.com/yushoteke/PlayAroundPackage.jl/blob/{commit}{path}#L{line}",
    sitename="PlayAroundPackage.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://yushoteke.github.io/PlayAroundPackage.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/yushoteke/PlayAroundPackage.jl",
)

max_stars_repo_name: yushoteke/PlayAroundRepository

max_issues_repo_name: yushoteke/PlayAroundRepository

max_forks_repo_name: yushoteke/PlayAroundRepository


text:
@testset "entropy" begin
    ent = entropy(cars, 7, 4)
    ent2 = entropy(cars, "class", header, 4)
    @test ent >= 0 && ent <= 1
    @test ent2 >= 0 && ent2 <= 1
    @test ent == ent2
end

@testset "information gain" begin
    for i in 1:6
        ig = information_gain(cars, i, 7, 4)
        ig2 = information_gain(cars, header[i], "class", header, 4)
        ig3 = information_gain(cars, header[i], 7, header, 4)
        ig4 = information_gain(cars, i, "class", header, 4)
        @test ig >= 0 && ig <= 1
        @test ig2 >= 0 && ig2 <= 1
        @test ig3 >= 0 && ig3 <= 1
        @test ig4 >= 0 && ig4 <= 1
        @test ig == ig2 == ig3 == ig4
    end
end


text:
using Compat
using MLBase
using Distributions
import Base: length, convert, promote_rule

if VERSION >= v"0.4.0-dev"
    typealias Range1{Int} Range{Int}
    _int(x) = round(Integer, x)
    float(x) = map(Float64, x)
else
    _int(x) = int(x)
end
squeeze(v::Vector, i::Integer) = v

const NO_BEST=(0,0)

immutable Leaf
    majority::Any
    values::Vector
end

@compat immutable Node
    featid::Integer
    featval::Any
    left::Union{Leaf,Node}
    right::Union{Leaf,Node}
end

@compat typealias LeafOrNode Union{Leaf,Node}

immutable Ensemble
    trees::Vector{Node}
end

convert(::Type{Node}, x::Leaf) = Node(0, nothing, x, Leaf(nothing,[nothing]))
promote_rule(::Type{Node}, ::Type{Leaf}) = Node
promote_rule(::Type{Leaf}, ::Type{Node}) = Node

length(ensemble::Ensemble) = length(ensemble.trees)

include("measures.jl")
include("classification.jl")
include("auc.jl")

max_stars_repo_name: paulstey/boosting

max_issues_repo_name: paulstey/boosting

max_forks_repo_name: paulstey/boosting


text:
# ps3000_get_unit_info.jl

export ps3000_get_unit_info, PS3000GetUnitInfoLine

function ps3000_get_unit_info(handle :: Int16,
                              line)
  string = Array(UInt8,32)
  returned_string_length = ccall((:ps3000_get_unit_info, ps3000driver), Int16,
    (Int16, Ptr{UInt8}, Int16, Int16),
    handle, string, length(string), line)
  if returned_string_length == 0
    error("ps3000 Error: Parameter Out Of Range")
  end
  return bytestring(string[1:returned_string_length])
end

function ps3000_get_unit_info(handle :: Int16)
  info = Array(ASCIIString,7)
  for line in Int16(0):Int16(6)
    info[line+1] = ps3000_get_unit_info(handle,line)
  end
  return info
end

baremodule PS3000GetUnitInfoLine
  const PS3000_DRIVER_VERSION = 0
  const PS3000_USB_VERSION = 1
  const PS3000_HARDWARE_VERSION = 2
  const PS3000_VARIANT_INFO = 3
  const PS3000_BATCH_AND_SERIAL = 4
  const PS3000_CAL_DATE = 5
  const PS3000_ERROR_CODE = 6   # see PS3000DriverErrorCodes
end


text:
import Test: @test, @test_throws, @testset
using Hankel
import LinearAlgebra: diagm, mul!, ldiv!
import SpecialFunctions: besseli, besselix, besselj
import HCubature: hquadrature
using Random
using ChainRulesCore
using ChainRulesTestUtils

# Brute-force equivalent of Hankel.dot! - slow but certain to be correct
function slowdot!(out, M, V; dim=1)
    idxlo = CartesianIndices(size(V)[1:dim-1])
    idxhi = CartesianIndices(size(V)[dim+1:end])
    _slowdot!(out, M, V, idxlo, idxhi)
end

function _slowdot!(out, M, V, idxlo, idxhi)
    for lo in idxlo
        for hi in idxhi
            view(out, lo, :, hi) .= M * view(V, lo, :, hi)
        end
    end
end

@testset "multiplication" begin
    M = diagm(0 => [1, 2, 3])
    V = 2 .* ones((3, 3, 2))
    out = similar(V)
    Hankel.dot!(out, M, V)
    @test all(out[1, :, :] .== 2)
    @test all(out[2, :, :] .== 4)
    @test all(out[3, :, :] .== 6)
    Hankel.dot!(out, M, V, dim=2)
    @test all(out[:, 1, :] .== 2)
    @test all(out[:, 2, :] .== 4)
    @test all(out[:, 3, :] .== 6)
    @test_throws DomainError Hankel.dot!(out, M, V, dim=3)

    M = rand(32, 32)
    V = rand(32)
    for N = 1:5
        for n = 1:N
            shape = 16*ones(Int64, N)
            shape[n] = 32
            V = rand(shape...)
            out = similar(V)
            out2 = similar(out)
            Hankel.dot!(out, M, V, dim=n)
            slowdot!(out2, M, V, dim=n)
            @test all(out2 .≈ out)
        end
    end
end

@testset "sphbesselj functions" begin
    @testset "sphbessel_scale" begin
        @test Hankel.sphbesselj_scale(1) == 1
        @test Hankel.sphbesselj_scale(3) == 1
        @test Hankel.sphbesselj_scale(5) == 1
        @test Hankel.sphbesselj_scale(2) == √(π/2)
        @test Hankel.sphbesselj_scale(4) == √(π/2)
        @test Hankel.sphbesselj_scale(6) == √(π/2)
    end

    @testset "sphbesselj" begin
        @test Hankel.sphbesselj(0, 1, 0) ≈ 1
        @test Hankel.sphbesselj(0, 2, 0) ≈ 1
        @test Hankel.sphbesselj(0, 3, 0) ≈ 1 / 2
        @test Hankel.sphbesselj(1, 1, 0) ≈ 0
        @test Hankel.sphbesselj(1, 2, 0) ≈ 0
        @test Hankel.sphbesselj(2, 1, 0) ≈ 0
        @test Hankel.sphbesselj(2, 2, 0) ≈ 0
        @test Hankel.sphbesselj(0, 1, 0.2) == besselj(0, 0.2)
        @test Hankel.sphbesselj(1, 1, 0.2) == besselj(1, 0.2)
        @test Hankel.sphbesselj(2, 1, 0.2) == besselj(2, 0.2)
        @test Hankel.sphbesselj(0, 2, 0.2) ≈ besselj(0.5, 0.2) * √(π / 2 / 0.2)
        @test Hankel.sphbesselj(1, 2, 0.2) ≈ besselj(1.5, 0.2) * √(π / 2 / 0.2)
        @test Hankel.sphbesselj(2, 2, 0.2) ≈ besselj(2.5, 0.2) * √(π / 2 / 0.2)
        @test Hankel.sphbesselj(0, 3, 0.2) == besselj(1, 0.2) / 0.2
        @test Hankel.sphbesselj(1, 3, 0.2) == besselj(2, 0.2) / 0.2
        @test Hankel.sphbesselj(2, 3, 0.2) == besselj(3, 0.2) / 0.2
        @test Hankel.sphbesselj(0, 4, 0.2) ≈ besselj(1.5, 0.2) * √(π / 2 / 0.2^3)
        @test Hankel.sphbesselj(1, 4, 0.2) ≈ besselj(2.5, 0.2) * √(π / 2 / 0.2^3)
        @test Hankel.sphbesselj(2, 4, 0.2) ≈ besselj(3.5, 0.2) * √(π / 2 / 0.2^3)
    end

    @testset "sphbesselj_zero($p, $n, m) for m = 1:200" for p in [1:4; 10; 20; 50; 100; 500], n in 1:4
        zs = Hankel.sphbesselj_zero.(p, n, 1:200)
        # test sequence is monotonic
        @test all(>(0), zs[2:end] .- zs[1:end-1])
        # test are zeros
        @test all(x -> isapprox(x, 0; atol=1e-12), Hankel.sphbesselj.(p, n, zs))
    end
end

include("qdht.jl")
include("chainrules.jl")

max_stars_repo_name: chrisbrahms/Hankel

max_issues_repo_name: chrisbrahms/Hankel

max_forks_repo_name: chrisbrahms/Hankel


text:
function movement!(r,side)
    while !isborder(r,side)
        move!(r,side)
    end
end
function begin!(r)
    for side in (West,Sud)
        while !isborder(r,side)
            move!(r,side)
        end
    end
end
function chess!(r,num_step)
    if num_step % 2 == 0
        putmarker!(r)
    end
end
function chehov!(r,num_step)
    num_step = 0
    for side in (Ost, West)
        while !isborder(r,side)
            chess!(r,num_step)
            num_step += 1
            move!(r,side)
        end
        num_step += 1
        move!(r,Nord)
    end
end
function go!(r)
    begin!(r)
    num_step = 0
    while (!isborder(r,West) || !isborder(r,Ost)) && !isborder(r,Nord)
        chehov!(r,num_step)
    end
    side = West
    while !isborder(r,inverse(side))
        chess!(r,num_step)
        num_step += 1
        move!(r,inverse(side))
    end
    while !isborder(r,side)
        chess!(r,num_step)
        num_step += 1
        move!(r,side)
    end
end
inverse(side::HorizonSide) = HorizonSide(mod(Int(side)+2,4))