text:
"""
    get_until_correct(prompt::String, dtype::Type)

Prompts user until a valid input is given
"""
function get_until_correct(prompt::String, dtype::Type; parser::Function = parse)
    out = nothing
    try
        print(prompt)
        out = readline()
        out = parse(dtype,out)
    catch
        println("Invalid entry. Try again...")
        out = get_until_correct(prompt,dtype)
    end
    return out
end
export get_until_correct

function parse_array(a::String, dtype::Type; token::String = ",")
    out_str = strip.(split(a,token))
    out = parse.(dtype,out_str)
    return out
end

"""
    request_particle_type()

Prompts user to enter a particle type. Returns a ParticleID.

"""
function request_particleid()
    particle_menu = request("Select a particle:",RadioMenu(particle_keys))
    particle = particle_dict[particle_keys[particle_menu]]
    return particle
end

"""
    request_energy()

Prompts user to enter energy. Returns Quantity{Float64,Energy}
"""
function request_energy()
    eunits = ["MeV","GeV","TeV"]
    eunits_menu = request("Select energy units:",RadioMenu(eunits))
    eunit = eunits[eunits_menu]
    eprompt = "Enter energy value ($units):"
    eval = get_until_correct(eprompt,Float64)
    eunit = uparse(eunit)
    energy = eval*eunit
    return energy
end

"""
    request_position()

Prompts user to enter a position. Returns Array{Quantity{Float64,Length},1}
"""
function request_position()
    posprompt = "Enter initial position in kpc as `XX.XX, YY.YY, ZZ.ZZ`: "
    pos = get_until_correct(posprompt,Float64,parse_array)
    return pos
end

"""
    request_direction()

Prompts user to enter a direction. Returns Array{Float64,1}
"""
function request_direction()
    dirprompt = "Enter initial direction as `XX.XX, YY.YY, ZZ.ZZ`: "
    dir = get_until_correct(dirprompt,Float64,parse_array)
    return dir
end

"""
    request_particle_ic()

Prompts user to enter initial conditions for a single particle
"""
function request_particle_ic()
    particle = request_particleid()
    energy   = request_energy()
    position = request_position()
    direction = request_direction()
    dtmenu = ["Automatic", "Manual"]
    dtselect = request("Time stepping: ",RadioMenu(dtmenu))
    dt = nothing
    if dtselect == 1
        dtprompt = "Enter dt (s): "
        dt       = get_until_correct(dtprompt,Float64)
    else
        timeprompt = "Enter duration (s): "
        time       = get_until_correct(timeprompt, Float64)
    end

end

"""
    particle_prompt()

Gets input from user to construct initial particle tree
"""
function particle_prompt()

end

"""
    clear!(scene::Scene)

Removes everything from scene except axes
"""
function clear!(scene::Scene)
    if length(scene.plots) > 1
        delete!(scene,scene.plots[end])
        clear!(scene)
    elseif length(scene.plots) == 1
        delete!(scene,scene.plots[end])
    end
end

export clear!

"""
    scaleShift(a::T) where {T<:AbstractVector}

Shifts values of array so it is centered on its mean and range is normalized to 1
"""
function scaleshift(a::T) where {T<:AbstractVector}
    aRange = max(a...) - min(a...)
    aMean  = mean(a)
    out = a .- aMean
    if aRange == 0.0
        return out
    else
        return aMean, aRange, out ./ aRange
    end
end

export scaleshift


text:
square(x) = x .* x

struct ReadNoAttention
end

function (l::ReadNoAttention)(x, xhat, hdec;
                              Fx=nothing, Fy=nothing, gamma=nothing)
    vcat(x, xhat)
end

WriteNoAttention = Linear

struct QNet
    mu_layer
    logsigma_layer
end

function (l::QNet)(henc)
    mu = l.mu_layer(henc)
    logsigma = l.logsigma_layer(henc)
    sigma = exp.(logsigma)
    if SEED == -1
        noise = randn!(similar(mu))
    else
        Knet.seed!(SEED)
        noise = randn!(similar(mu))
    end
    sampled = mu .+ noise .* sigma
    return (sampled, mu, logsigma, sigma)
end

function QNet(input_dim::Int, output_dim::Int, atype=_atype, init=xavier)
    mu_layer = Linear(input_dim, output_dim; atype=atype, init=init)
    logsigma_layer = Linear(input_dim, output_dim; atype=atype, init=init)
    return QNet(mu_layer, logsigma_layer)
end

function sample_noise(q::QNet, batchsize::Int; generation=false)
    zdim = size(value(q.mu_layer.w), !generation ? 2 : 1)
    if SEED == -1
        z = randn(zdim, batchsize)
    else
        Knet.seed!(SEED)
        z = randn(zdim, batchsize)
    end
    atype = typeof(value(q.mu_layer.w))
    return convert(atype, z)
end

#
# Attention-related modules
#

mutable struct AttentionWindow
    linear
    A
    B
    N
end

function (l::AttentionWindow)(hdec)
    ps = l.linear(hdec)
    hsize, batchsize = size(ps)
    gx, gy, logsigma2, logdelta, loggamma = map(
        i->reshape(ps[i, :], 1, batchsize) , 1:hsize)
    gx = div(l.A+1, 2) .* (gx .+ 1)
    gy = div(l.B+1, 2) .* (gy .+ 1)
    delta = div(max(l.A, l.B) - 1, (l.N - 1)) .* exp.(logdelta)
    sigma2 = exp.(logsigma2)
    gamma = exp.(loggamma)

    Fx, Fy = filterbank(gx, gy, sigma2, delta, l.A, l.B, l.N)
    return Fx, Fy, gamma
end

function AttentionWindow(input_dim::Int, A::Int, B::Int, N::Int; atype=_atype)
    linear = Linear(input_dim, 5; atype=atype)
    return AttentionWindow(linear, A, B, N)
end

function filterbank(gx, gy, sigma2, delta, A, B, N)
    atype = typeof(value(gx)) <: KnetArray ? KnetArray : Array
    etype = eltype(gx)
    batchsize = size(gx, 2)

    rng = reshape(atype{etype}(0:N-1), N, 1)
    mu_x = compute_mu(gx, rng, delta, N)
    mu_y = compute_mu(gy, rng, delta, N)

    a = reshape(atype{etype}(0:A-1), A, 1, 1)
    b = reshape(atype{etype}(0:B-1), B, 1, 1)

    mu_x = reshape(mu_x, 1, size(mu_x)...)
    mu_y = reshape(mu_y, 1, size(mu_y)...)
    sig2 = reshape(sigma2, 1, 1, length(sigma2))

    Fx = filterbank_matrices(a, mu_x, sig2)
    Fy = filterbank_matrices(b, mu_y, sig2)

    return Fx, Fy
end

function filterbank_matrices(xs, mu, sigma2; epsilon=1e-9)
    y = xs .- mu
    y = y ./ 2sigma2
    y = exp.(-y .* y)
    y = y ./ (sum(y, dims=1) .+ epsilon)
end

function compute_mu(g, rng, delta, N)
    tmp = (rng .- 0.5f0 .- div(N, 2)) .* delta
    mu = tmp .+ g
end

function filter_image(image, Fxt, Fy, gamma, A, B, N)
    batchsize = size(image, 2)
    img = reshape(image, A, B, batchsize)
    glimpse = bmm(bmm(Fxt, img), Fy)
    glimpse = reshape(glimpse, N, N, batchsize)
    out = glimpse .* reshape(gamma, 1, 1, batchsize)
    out = reshape(out, div(length(out), batchsize), batchsize)
end

mutable struct ReadAttention
    window
end

function (l::ReadAttention)(x, xhat, hdec)
    Fx, Fy, gamma = l.window(hdec)
    A, B, N = l.window.A, l.window.B, l.window.N
    Fxt = permutedims(Fx, (2, 1, 3))
    xnew = filter_image(x, Fxt, Fy, gamma, A, B, N)
    xhatnew = filter_image(xhat, Fxt, Fy, gamma, A, B, N)
    return vcat(xnew, xhatnew)
end

mutable struct WriteAttention
    window
    linear
end

function (l::WriteAttention)(hdec)
    A, B, N = l.window.A, l.window.B, l.window.N
    batchsize = size(hdec, 2)
    w = l.linear(hdec)
    w = reshape(w, N, N, batchsize)
    Fx, Fy, gamma = l.window(hdec)
    Fyt = permutedims(Fy, (2, 1, 3))
    wr = bmm(bmm(Fx, w), Fyt)
    wr = reshape(wr, A*B, batchsize)
    return wr ./ gamma
end

function WriteAttention(window, decoder_dim, N; atype=_atype)
    linear = Linear(decoder_dim, N*N; atype=atype)
    return WriteAttention(window, linear)
end

ext: jl

max_stars_repo_name: ilkerkesen/draw

max_issues_repo_name: ilkerkesen/draw

max_forks_repo_name: ilkerkesen/draw


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

"""
    partition_space(
		samples::DensitySampleVector,
		npartitions::Integer,
		algorithm::KDTreePartitioning
	)

*BAT-internal, not part of stable public API.*

The function generates a space partition tree with the number of partitions
given by `npartitions`, using `KDTreePartitioning` algorithm. Exploration samples
are given by `samples`. The output contains `SpacePartTree` and the values of
the cost function.
"""
function partition_space(samples::DensitySampleVector, npartitions::Integer, algorithm::KDTreePartitioning)

	# Maybe there is a better way to get n_params:
	n_params = size(flatview(unshaped.(samples.v)))[1]

	# Check whether dimensions for partition are specified. Use all if not.
	if algorithm.partition_dims == :auto
		partition_dims = collect(Base.OneTo(n_params))
	else
		if !prod([i in Base.OneTo(n_params) for i in algorithm.partition_dims])
			@warn "Partition parameters are not correct"
		end
		partition_dims = sort(intersect(algorithm.partition_dims, collect(Base.OneTo(n_params)))) # to be safe
	end

	flat_scaled_data, μ, δ = scale_samples(samples)
	bounds = repeat([0.0 1.0],size(flat_scaled_data.samples)[1])
	partition_tree = def_init_node(flat_scaled_data, bounds)
	cost_values = Float64[]
	for i in 1:npartitions-1
		# @info "KDTreePartitioning: Increasing tree depth (depth = $i)"
		initialize_partitioning!(partition_tree, flat_scaled_data, partition_dims)
		ind, sum_cost = det_part_node(partition_tree)
		append!(cost_values, sum_cost)
		generate_node!(partition_tree, flat_scaled_data, ind)
	end

	rescale_tree!(partition_tree, μ, δ)

	return partition_tree, cost_values
end

function scale_samples(samples::DensitySampleVector)
	flat_samples = collect(flatview(unshaped.(samples.v)))
    μ = minimum(flat_samples, dims=2)
    δ = maximum(flat_samples, dims=2).-minimum(flat_samples, dims=2)
    return (samples = (flat_samples .- μ) ./ δ, weights = samples.weight, loglik = samples.logd), μ, δ
end

function def_init_node(data::T, bounds::Array{F}) where {T<:NamedTuple, F<:AbstractFloat}
    cost = evaluate_total_cost(data)
    return SpacePartTree(true, bounds, missing, missing, missing, missing, cost, missing)
end

function initialize_partitioning!(tree::SpacePartTree, data::T, axes::Array{I,1}) where {T<:NamedTuple, I<:Integer}
    if tree.terminated_leaf == true
        if ismissing(tree.cost_part)
            cost, (div_axis, cut_pos) = find_min_along_axes(mask_data(data, tree.bounds), axes)
            tree.cost_part = cost
            tree.cut_coordinate = cut_pos
            tree.cut_axis = div_axis
            return tree
        else
            return tree
        end
    else
        initialize_partitioning!(tree.left_child, data, axes)
        initialize_partitioning!(tree.right_child, data, axes)
    end
end

function find_min_along_axes(data::T, axes::Array{I,1}) where {T<:NamedTuple, I<:Integer}

    min_function = Float64[]
    min_positions = Float64[]

    for i in axes
        _, _, min_x, min_y = find_min_along_axis(data, i)
        push!(min_function, min_y)
        push!(min_positions, min_x)
    end
    min_glob_axis = argmin(min_function)
    divide_along = axes[min_glob_axis]
    cut_position = min_positions[min_glob_axis]
    return minimum(min_function), (divide_along, cut_position)
end

function find_min_along_axis(data::T, axis::Integer; n_δ=50) where {T<:NamedTuple}

    grid = range(minimum(data.samples[axis,:]),stop=maximum(data.samples[axis,:]), length=n_δ)
    cost_values = Float64[]

    for cut in collect(grid)
        cost_value = evaluate_param_cost(data, cut, axis)
        push!(cost_values, cost_value)
    end
    nan_mask = .!isnan.(cost_values)
    min_value = minimum(cost_values[nan_mask])

	min_indicies = findall(x->x==min_value, cost_values[nan_mask])
	mid_index = 1 + (length(min_indicies)-1)÷2
	mid_index = min_indicies[mid_index]

    # min_bin = argmin(cost_values[nan_mask])
    min_bin_value = grid[nan_mask][mid_index]

    return (grid, cost_values, min_bin_value, min_value)
end

function evaluate_param_cost(data::T, cut_position::F, axis::I) where {T<:NamedTuple, I<:Integer, F<:AbstractFloat}
    size_s = size(data.samples)
	n_samples_critical = 3 #if less than 3 samples then return Inf
    if size_s[2] > n_samples_critical
        mask = data.samples[axis, :] .< cut_position
		a1 = evaluate_total_cost((samples = data.samples[:,mask], weights = data.weights[mask], loglik = data.loglik[mask]))
		a2 = evaluate_total_cost((samples = data.samples[:,.!mask], weights = data.weights[.!mask], loglik = data.loglik[.!mask]))
        return  a1+a2
    else
        return Inf
    end
end

function mask_data(data::T, bounds::Array{F}) where {T<:NamedTuple, F<:AbstractFloat}
    vect_mask = [prod(bounds[:,2] .> data.samples .> bounds[:,1], dims=1)...]
    return (samples = data.samples[:,vect_mask], weights = data.weights[vect_mask], loglik = data.loglik[vect_mask] )
end

function det_part_node(tree::SpacePartTree)
    costs_arrays::Array{Float64} = []
    cost_delta::Array{Float64} = []
    collect_part_costs!(tree, costs_arrays, cost_delta)
    ind = argmax(cost_delta)
    return costs_arrays[ind], sum(costs_arrays)
end

function collect_part_costs!(tree::SpacePartTree, cost_values::Array{Float64,1}, cost_delta::Array{Float64,1})
    if tree.terminated_leaf == true
        append!(cost_values, tree.cost)
        append!(cost_delta, tree.cost - tree.cost_part)
    else
        collect_part_costs!(tree.left_child,cost_values, cost_delta)
        collect_part_costs!(tree.right_child,cost_values, cost_delta)
    end
end

function generate_node!(tree::SpacePartTree, data::T, cut_val::F) where {T<:NamedTuple, F<:AbstractFloat}
    if tree.terminated_leaf == true
        if tree.cost == cut_val

            new_bounds_left = copy(tree.bounds)
            new_bounds_right = copy(tree.bounds)

            new_bounds_left[tree.cut_axis, 2] = tree.cut_coordinate
            new_bounds_right[tree.cut_axis, 1] = tree.cut_coordinate

            tree.left_child = def_init_node(mask_data(data, new_bounds_left), new_bounds_left)
            tree.right_child = def_init_node(mask_data(data, new_bounds_right), new_bounds_right)
            tree.terminated_leaf = false

            return tree
        else
            return tree
        end

    else
        generate_node!(tree.left_child, data, cut_val)
        generate_node!(tree.right_child, data, cut_val)
    end
end

function evaluate_total_cost(data::T) where {T<:NamedTuple}
    sise_s = size(data.samples)
    if sise_s[2] > 3
        μ = mean(data.samples, weights(data.weights), dims=2)
        return sum(data.weights.*sum((data.samples .- μ).^2, dims=1))
    else
        return Inf
    end
end

function cost_f_2(data::T) where {T<:NamedTuple}
	# Alternative Cost Function #2
    sise_s = size(data.samples)
    if sise_s[2] > 3
		μ = mean(exp.(data.loglik), weights(data.weights))
        return sum(data.weights.*(exp.(data.loglik) .- μ).^2)
    else
        return Inf
    end
end

function cost_f_3(data::T) where {T<:NamedTuple}
	# Alternative Cost Function #3
    sise_s = size(data.samples)
    if sise_s[2] > 3
        μ = data.samples[:,argmax(data.weights)]
        return   sum( data.weights.*sum((data.samples .- μ).^2, dims=1))
    else
        return Inf
    end
end

function rescale_tree!(tree::SpacePartTree, μ::Array{F}, δ::Array{F}) where {F<:AbstractFloat}
    if tree.terminated_leaf == true
        tree.bounds = tree.bounds.*δ .+ μ
    else
        tree.bounds = tree.bounds.*δ .+ μ
        rescale_tree!(tree.left_child, μ, δ)
        rescale_tree!(tree.right_child, μ, δ)
    end
end

function get_tree_par_bounds!(tree::SpacePartTree, bounds_part::Array{Array{F},1}) where {F<:AbstractFloat}
    if tree.terminated_leaf == true
        push!(bounds_part, tree.bounds)
    else
        get_tree_par_bounds!(tree.left_child, bounds_part)
        get_tree_par_bounds!(tree.right_child, bounds_part)
    end
end

function get_tree_par_bounds(tree::SpacePartTree)
    param_bounds::Array{Array{Float64},1} = []
	return get_tree_par_bounds!(tree, param_bounds)
end

function extend_tree_bounds!(tree::SpacePartTree, lo::Array{F,1}, hi::Array{F,1}) where {F<:AbstractFloat}

	subspaces_rect_bounds = get_tree_par_bounds(tree)

	lo_tree_bounds = [minimum(hcat([tree_bound[:,1] for tree_bound in subspaces_rect_bounds]...), dims=2)...]
    hi_tree_bounds = [maximum(hcat([tree_bound[:,2] for tree_bound in subspaces_rect_bounds]...), dims=2)...]

	to_be_extended = [Pair.(lo_tree_bounds, lo); Pair.(hi_tree_bounds, hi)]

	extend_tree_bounds!(tree, to_be_extended)
end

function extend_tree_bounds!(tree::SpacePartTree, to_be_extended::Array{P}) where {P<:Pair}
    if tree.terminated_leaf == true
		tree.bounds = replace(tree.bounds, to_be_extended...)
    else
        extend_tree_bounds!(tree.left_child, to_be_extended)
        extend_tree_bounds!(tree.right_child, to_be_extended)
    end
end


text:
###############################################################################
## Coefficient Tableaus
###############################################################################

# Butcher Tableaus, or more generally coefficient tables
# see Hairer & Wanner 1992, p. 134, 166

abstract Tableau{Name, S, T<:Real}
# Name is the name of the tableau/method (a symbol)
# S is the number of stages (an int)
# T is the type of the coefficients
#
# TODO: have a type parameter which specifies adaptive vs non-adaptive
#
# For all types of tableaus it assumes fields:
# order::(Int...) # order of the method(s)
#
# For Runge-Kutta methods it assumes fields:
# a::Matrix{T}  # SxS matrix
# b::Matrix{T}  # 1 or 2 x S matrix (fixed step/ adaptive)
# c::Vector{T}  # S
#
# For a tableau:
#  c1  | a_11   ....   a_1s
#  .   | a_21 .          .
#  .   | a_31     .      .
#  .   | ....         .  .
#  c_s | a_s1  ....... a_ss
# -----+--------------------
#      | b_1     ...   b_s   this is the one used for stepping
#      | b'_1    ...   b'_s  this is the one used for error-checking

Base.eltype{N,S,T}(b::Tableau{N,S,T}) = T
order(b::Tableau) = b.order
# Subtypes need to define a convert method to convert to a different
# eltype with signature:
Base.convert{Tnew<:Real}(::Type{Tnew}, tab::Tableau) = error("Define convert method for concrete Tableau types")

###########################################
# Tableaus for explicit Runge-Kutta methods
###########################################

# m3: these tableaus should go into rk.jl as they belong to R-K methods

immutable TableauRKExplicit{Name, S, T} <: Tableau{Name, S, T}
    order::(@compat(Tuple{Vararg{Int}})) # the order of the methods
    a::Matrix{T}
    # one or several row vectors.  First row is used for the step,
    # second for error calc.
    b::Matrix{T}
    c::Vector{T}
    function TableauRKExplicit(order,a,b,c)
        @assert isa(S,Integer)
        @assert isa(Name,Symbol)
        @assert c[1]==0
        @assert istril(a)
        @assert S==length(c)==size(a,1)==size(a,2)==size(b,2)
        @assert size(b,1)==length(order)
        @assert norm(sum(a,2)-c'',Inf)<T(1e-10) # consistency.
        new(order,a,b,c)
    end
end

function TableauRKExplicit{T}(name::Symbol, order::(@compat(Tuple{Vararg{Int}})),
                              a::Matrix{T}, b::Matrix{T}, c::Vector{T})
    TableauRKExplicit{name,length(c),T}(order, a, b, c)
end

function TableauRKExplicit(name::Symbol, order::(@compat(Tuple{Vararg{Int}})), T::Type,
                           a::Matrix, b::Matrix, c::Vector)
    TableauRKExplicit{name,length(c),T}(order, convert(Matrix{T},a),
                                        convert(Matrix{T},b), convert(Vector{T},c) )
end

conv_field{T,N}(D,a::Array{T,N}) = convert(Array{D,N}, a)

lengthks(tab::TableauRKExplicit) = length(tab.c)

function Base.convert{Tnew<:Real,Name,S,T}(::Type{Tnew}, tab::TableauRKExplicit{Name,S,T})
    # Converts the tableau coefficients to the new type Tnew
    newflds = ()
    @compat for n in fieldnames(tab)
        fld = getfield(tab,n)
        if eltype(fld)==T
            newflds = tuple(newflds..., conv_field(Tnew, fld))
        else
            newflds = tuple(newflds..., fld)
        end
    end
    TableauRKExplicit{Name,S,Tnew}(newflds...) # TODO: could this be done more generically in a type-stable way?
end

isexplicit(b::TableauRKExplicit) = istril(b.a) # Test whether it's an explicit method
isadaptive(b::TableauRKExplicit) = size(b.b, 1)==2

# First same as last.  Means ks[:,end]=ks_nextstep[:,1], c.f. H&W p.167
isFSAL(btab::TableauRKExplicit) = btab.a[end,:]==btab.b[1,:] && btab.c[end]==1 # the latter is not needed really

## Tableaus for explicit RK methods
# Fixed step:
const bt_feuler = TableauRKExplicit(:feuler,(1,), Rational{Int64},
                                    zeros(Int,1,1),
                                    [1]',
                                    [0]
                                    )
const bt_midpoint = TableauRKExplicit(:midpoint,(2,), Rational{Int64},
                                      [0  0
                                       1//2  0],
                                      [0, 1]',
                                      [0, 1//2]
                                      )
const bt_heun = TableauRKExplicit(:heun,(2,), Rational{Int64},
                                  [0  0
                                   1  0],
                                  [1//2, 1//2]',
                                  [0, 1])

const bt_rk4 = TableauRKExplicit(:rk4,(4,),Rational{Int64},
                                 [0    0    0 0
                                  1//2 0    0 0
                                  0    1//2 0 0
                                  0    0    1 0],
                                 [1//6, 1//3, 1//3, 1//6]',
                                 [0, 1//2, 1//2, 1])

# Adaptive step:
# Heun Euler https://en.wikipedia.org/wiki/Runge–Kutta_methods
const bt_rk21 = TableauRKExplicit(:heun_euler,(2,1), Rational{Int64},
                                  [0     0
                                   1     0],
                                  [1//2  1//2
                                   1     0],
                                  [0,    1])

# Bogacki–Shampine coefficients
const bt_rk23 = TableauRKExplicit(:bogacki_shampine,(2,3), Rational{Int64},
                                  [0           0      0      0
                                   1/2         0      0      0
                                   0         3/4      0      0
                                   2/9       1/3     4/9     0],
                                  [7/24 1/4 1/3 1/8
                                   2/9 1/3 4/9 0],
                                  [0, 1//2, 3//4, 1]
                         )

# Fehlberg https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta%E2%80%93Fehlberg_method
const bt_rk45 = TableauRKExplicit(:fehlberg,(4,5),Rational{Int64},
                          [  0           0           0            0         0     0
                             1//4        0           0            0         0     0
                             3//32       9//32       0            0         0     0
                          1932//2197 -7200//2197  7296//2197      0         0     0
                           439//216     -8        3680//513    -845//4104   0     0
                            -8//27       2       -3544//2565   1859//4104 -11//40 0 ],
                           [25//216      0        1408//2565   2197//4104  -1//5  0
                            16//135      0        6656//12825 28561//56430 -9//50 2//55],
                            [0,          1//4,       3//8,       12//13,    1,    1//2])

# Dormand-Prince https://en.wikipedia.org/wiki/Dormand%E2%80%93Prince_method
const bt_dopri5 = TableauRKExplicit(:dopri, (5,4), Rational{Int64},
                     [0           0            0                   0            0      0 0
                      1//5        0            0                   0            0      0 0
                      3//40       9//40        0                   0            0      0 0
                      44//45      -56//15      32//9               0            0      0 0
                      19372//6561 -25360//2187 64448//6561 -212//729            0      0 0
                      9017//3168  -355//33     46732//5247   49//176 -5103//18656      0 0
                      35//384     0            500//1113    125//192  -2187//6784 11//84 0],
                     [35//384         0     500//1113       125//192      -2187//6784         11//84      0
                      5179//57600     0     7571//16695     393//640     -92097//339200     187//2100     1//40],
                     [0, 1//5, 3//10, 4//5, 8//9, 1, 1]
                     )

# Fehlberg 7(8) coefficients
# Values from pag. 65, Fehlberg, Erwin. "Classical fifth-, sixth-, seventh-, and eighth-order Runge-Kutta formulas with stepsize control".
# National Aeronautics and Space Administration.
const bt_feh78 = TableauRKExplicit(:feh78, (7,8), Rational{Int64},
                            [     0       0      0       0         0          0        0        0      0       0     0 0 0
                                  2//27   0      0       0         0          0        0        0      0       0     0 0 0
                                  1//36   1//12  0       0         0          0        0        0      0       0     0 0 0
                                  1//24   0      1//8    0         0          0        0        0      0       0     0 0 0
                                  5//12   0    -25//16  25//16     0          0        0        0      0       0     0 0 0
                                  1//20   0      0       1//4      1//5       0        0        0      0       0     0 0 0
                                -25//108  0      0     125//108  -65//27    125//54    0        0      0       0     0 0 0
                                 31//300  0      0       0       61//225    -2//9     13//900   0      0       0     0 0 0
                                  2       0      0     -53//6    704//45   -107//9    67//90    3      0       0     0 0 0
                                -91//108  0      0      23//108 -976//135   311//54  -19//60   17//6  -1//12   0     0 0 0
                               2383//4100 0      0    -341//164 4496//1025 -301//82 2133//4100 45//82 45//164 18//41 0 0 0
                                  3//205  0      0       0        0        -6//41     -3//205  -3//41  3//41   6//41 0 0 0
                              -1777//4100 0      0    -341//164 4496//1025 -289//82 2193//4100 51//82 33//164 12//41 0 1 0],
                              [41//840 0 0 0 0 34//105 9//35 9//35 9//280 9//280 41//840 0 0
                               0 0 0 0 0 34//105 9//35 9//35 9//280 9//280 0     41//840 41//840],
                               [0,    2//27, 1//9, 1//6 , 5//12, 1//2 , 5//6 , 1//6 , 2//3 , 1//3 , 1 , 0, 1]
                            )

# m3:
# - this function is not used anymore!  Is there a reason for this?
# - this should probably stay in this file.
# - update type-variables to
#     ET, EY, EF, CT, CY
#     - I think it is fine to assume that typeof(y0)==typeof(F(t0,y0))
#       i.e. dy has same type as y.
function make_consistent_types(fn, y0, tspan, btab::Tableau)
    # There are a few types involved in a call to a ODE solver which
    # somehow need to be consistent:
    #
    # Et = eltype(tspan)
    # Ey = eltype(y0)
    # Ef = eltype(Tf)
    #
    # There are also the types of the containers, but they are not
    # needed as `similar` is used to make containers.
    # Tt = typeof(tspan)
    # Ty = typeof(y0)              # note, this can be a scalar
    # Tf = typeof(F(tspan(1),y0))  # note, this can be a scalar
    #
    # Returns
    # - Et: eltype of time, needs to be a real "continuous" type, at
    #       the moment a AbstractFloat
    # - Eyf: suitable eltype of y and f(t,y)
    #   --> both of these are set to typeof(y0[1]/(tspan[end]-tspan[1]))
    # - Ty: container type of y0
    # - btab: tableau with entries converted to Et

    # Needed interface:
    # On components: /, -
    # On container: eltype, promote_type
    # On time container: eltype

    Ty = typeof(y0)
    Eyf = typeof(y0[1]/(tspan[end]-tspan[1]))

    Et = eltype(tspan)
    @assert Et<:Real
    if !(Et<:AbstractFloat)
        Et = promote_type(Et, Float64)
    end

    # if all are Floats, make them the same
    if Et<:AbstractFloat &&  Eyf<:AbstractFloat
        Et = promote_type(Et, Eyf)
        Eyf = Et
    end

    !isleaftype(Et) && warn("The eltype(tspan) is not a concrete type!  Change type of tspan for better performance.")
    !isleaftype(Eyf) && warn("The eltype(y0/tspan[1]) is not a concrete type!  Change type of y0 and/or tspan for better performance.")

    btab_ = convert(Et, btab)
    return Et, Eyf, Ty, btab_
end


text:
using MPSKit,TensorKit,Test

let
    ham = nonsym_ising_ham(lambda=4.0);
    gs = FiniteMPS(fill(TensorMap(rand,ComplexF64,ℂ^1*ℂ^2,ℂ^1),10));
    (gs,pars,_) = find_groundstate(gs,ham,Dmrg2(trscheme=truncdim(10)));

    #we are in the groundstate
    #we expect to find a single isolated pole around the gs energy
    polepos = real(sum(expectation_value(gs,ham,pars)));

    vals = (-0.5:0.05:0.5).+polepos
    eta = 0.3im;

    predicted = [1/(v+eta-polepos) for v in vals];

    data = map(vals) do v
        (v,_) = dynamicaldmrg(gs,v+eta,ham,verbose=false)
        v
    end

    @test norm(data-predicted) ≈ 0 atol=1e-8
end


text:
struct PoincareInvariant1stCanonical{ET,DT,TT}
    equ::ET
    Δt::TT
    nloop::Int
    ntime::Int
    nsave::Int
    nt::Int
    I::Vector{DT}
end

function PoincareInvariant1stCanonical(f_equ::Function, f_loop_q::Function, f_loop_p::Function, Δt::TT, d::Int, nloop::Int, ntime::Int, nsave::Int, DT=Float64) where {TT}

    if get_config(:verbosity) > 1
        println()
        println("First Euler-Poincaré Integral Invariant")
        println("=======================================")
        println()
        println(" nloop = ", nloop)
        println(" ntime = ", ntime)
        println(" nsave = ", nsave)
        println(" Δt    = ", Δt)
        println()
    end

    # compute initial conditions
    q₀ = zeros(DT, (d, nloop))
    p₀ = zeros(DT, (d, nloop))

    for i in 1:nloop
        q₀[:,i] .= f_loop_q(i/nloop)
        p₀[:,i] .= f_loop_p(i/nloop)
    end

    # initialise euation
    equ = f_equ(q₀, p₀)

    # create arrays for results
    nt = div(ntime, nsave)

    I = zeros(DT, nt+1)

    PoincareInvariant1stCanonical{typeof(equ),DT,TT}(equ, Δt, nloop, ntime, nsave, nt, I)
end

function evaluate_poincare_invariant(pinv::PoincareInvariant1stCanonical, sol::SolutionPODE)
    v = zeros(size(sol.q.d,1), size(sol.q.d,3))

    for i in 1:size(sol.q.d,2)
        compute_velocity(sol.q.d[:,i,:], v)
        pinv.I[i] = compute_loop_integral(sol.p.d[:,i,:], v)
    end

    pinv.I
end

function CommonFunctions.write_to_hdf5(pinv::PoincareInvariant1stCanonical, sol::Solution, output_file::String)
    # h5open(output_file, isfile(output_file) ? "r+" : "w") do h5
    h5open(output_file, "w") do h5
        write(h5, "t", sol.t.t)
        write(h5, "I", pinv.I)
    end
end


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    SDFileAtom, SmilesAtom,
    setcharge, setstereo, atomnumber, todict

const ATOMTABLE = let
    weightsfile = joinpath(dirname(@__FILE__), "../../assets/const/atomicweights.yaml")
    include_dependency(weightsfile)
    YAML.load(open(weightsfile))
end

const ATOMSYMBOLMAP = let
    symbolfile = joinpath(dirname(@__FILE__), "../../assets/const/symboltonumber.yaml")
    include_dependency(symbolfile)
    YAML.load(open(symbolfile))
end

const ATOM_COVALENT_RADII = let
    radiifile = joinpath(dirname(@__FILE__), "../../assets/const/covalent_radii.csv")
    include_dependency(radiifile)
    tab, headers = readdlm(radiifile, '\t', String; header=true, comments=true)
    radii = Dict{Int,Union{Float32,Dict{String,Float32}}}()
    for i = 1:size(tab, 1)
        an = parse(Int, tab[i, 1])
        container = if !haskey(radii, an)
            if i < size(tab, 1) && tab[i+1,1] == tab[i,1]  # special handling for elements with multiple options
                radii[an] = Dict{String,Float32}()
            else
                nothing
            end
        else
            radii[an]
        end
        ar = parse(Float32, tab[i, 3])
        if container === nothing
            radii[an] = ar
        else
            container[tab[i, 2]] = ar
        end
    end
    radii
end

const ATOM_VANDERWAALS_RADII = let
    radiifile = joinpath(dirname(@__FILE__), "../../assets/const/vanderWaals_radii.csv")
    include_dependency(radiifile)
    tab, headers = readdlm(radiifile, '\t', String; header=true, comments=true)
    radii = Dict{Int,Float32}()
    for i = 1:size(tab, 1)
        an = parse(Int, tab[i, 1])
        ar = parse(Float32, tab[i, 2])
        radii[an] = ar
    end
    radii
end

struct SDFileAtom <: Atom
    symbol::Symbol
    charge::Int
    multiplicity::Int
    mass::Union{Int, Nothing}
    coords::Union{Vector{Float64}, Nothing}
    stereo::Symbol

    function SDFileAtom(sym, chg, multi, mass, coords, stereo)
        new(sym, chg, multi, mass, coords, stereo)
    end
end

SDFileAtom() = SDFileAtom(:C, 0, 1, nothing, nothing)
SDFileAtom(sym) = SDFileAtom(sym, 0, 1, nothing, nothing)
SDFileAtom(sym, chg) = SDFileAtom(sym, chg, 1, nothing, nothing)
SDFileAtom(sym, chg, multi, mass, coords
    ) = SDFileAtom(sym, chg, multi, mass, coords, :unspecified)

SDFileAtom(data::Dict{String,Any}) = SDFileAtom(
    Symbol(data["symbol"]),
    data["charge"],
    data["multiplicity"],
    data["mass"],
    data["coords"],
    Symbol(data["stereo"])
)

function todict(a::SDFileAtom)
    data = Dict{String,Any}()
    for field in fieldnames(SDFileAtom)
        data[string(field)] = getfield(a, field)
    end
    return data
end

setcharge(a::SDFileAtom, chg
    ) = SDFileAtom(a.symbol, chg, a.multiplicity, a.mass, a.coords)

setstereo(a::SDFileAtom, direction) = SDFileAtom(
    a.symbol, a.charge, a.multiplicity, a.mass, a.coords, direction)

struct SmilesAtom <: Atom
    symbol::Symbol
    charge::Int
    multiplicity::Int
    mass::Union{Int, Nothing}
    isaromatic::Union{Bool, Nothing}
    stereo::Symbol

    function SmilesAtom(sym, chg, multi, mass, aromatic, stereo)
        new(sym, chg, multi, mass, aromatic, stereo)
    end
end

SmilesAtom() = SmilesAtom(:C, 0, 1, nothing, nothing, :unspecified)
SmilesAtom(sym) = SmilesAtom(sym, 0, 1, nothing, nothing, :unspecified)
SmilesAtom(sym, chg) = SmilesAtom(sym, chg, 1, nothing, nothing, :unspecified)

SmilesAtom(data::Dict{String,Any}) = SmilesAtom(
    Symbol(data["symbol"]),
    data["charge"],
    data["multiplicity"],
    data["mass"],
    data["isaromatic"],
    Symbol(data["stereo"])
)

function todict(a::SmilesAtom)
    data = Dict{String,Any}()
    for field in fieldnames(SmilesAtom)
        data[string(field)] = getfield(a, field)
    end
    return data
end

setcharge(a::SmilesAtom, chg) = SmilesAtom(
    a.symbol, chg, a.multiplicity, a.mass, a.isaromatic, a.stereo)

setstereo(a::SmilesAtom, direction) = SmilesAtom(
    a.symbol, a.charge, a.multiplicity, a.mass, a.isaromatic, direction)

"""
    atomnumber(atomsymbol::Symbol) -> Int
    atomnumber(atom::Atom) -> Int

Return atom number.
"""
atomnumber(atomsymbol::Symbol) = ATOMSYMBOLMAP[string(atomsymbol)]
atomnumber(atom::Atom) = atomnumber(atom.symbol)

"""
    atomsymbol(n::Int) -> Symbol

Return atom symbol of given atomic number.
"""
atomsymbol(n::Int) = Symbol(ATOMTABLE[n]["Symbol"])


text:
#!/usr/bin/env julia
# vim:set ft=julia ts=4 sw=4 sts=4 autoindent:

# Sanity checking for the perceptron module.
#
# Author:   Pontus Stenetorp    <pontus stenetorp se>
# Version:  2014-05-27

include("../hdr.jl")

require("percep.jl")

using Percep

olapath = "$respath/ola/fruit_data.tsv"

# Labels.
const APPLE     = "apple"
const ORANGE    = "orange"
const LEMON     = "lemon"
const LBLS      = (APPLE, ORANGE, LEMON)

# As per fruit_types.txt
const LBLBYID = [
    1   => APPLE,
    2   => ORANGE,
    3   => APPLE,
    4   => APPLE,
    5   => APPLE,
    6   => ORANGE,
    7   => ORANGE,
    8   => ORANGE,
    9   => LEMON,
    10  => LEMON,
]

const NUMFEATS = 3

type Sample
    label::String
    fvec::Vector{Float64}
end

function readdata(data_f)
    samples = Sample[]
    for line in eachline(data_f)
        idstr, fvecstr = split(rstrip(line, '\n'), '\t', 2)
        lbl = LBLBYID[int(idstr)]
        @assert lbl in LBLS
        fvec = map(float, split(fvecstr, '\t'))
        @assert length(fvec) == NUMFEATS
        push!(samples, Sample(lbl, map(float, split(fvecstr, '\t'))))
    end
    return samples
end

open(olapath, "r") do data_f
    samples = readdata(data_f)

    # Binary classification.
    let
        # Only use the oranges and lemons (they are linearly separable).
        data = collect(filter(s -> s.label in (ORANGE, LEMON), samples))
        percep = AvgPerceptron(NUMFEATS)

        wsum = zeros(size(percep.weights))
        bsum = 0
        count = 0

        errors = true
        while errors
            errors = false
            shuffle!(data)

            for sample in data
                islemon = classify(sample.fvec, percep)
                wsum += percep.weights
                bsum += percep.bias
                count += 1

                if isequal(sample.label, LEMON) && !islemon
                    update!(percep, sample.fvec, true)
                    errors = true
                elseif isequal(sample.label, ORANGE) && islemon
                    update!(percep, sample.fvec, false)
                    errors = true
                end
            end
        end

        @test count == percep.count
        vanilla = Perceptron(percep)
        # Compare with the averages, epsilon to account for rounding errors.
        wavg = wsum / count
        @test_approx_eq wavg vanilla.weights
        bavg = bsum / count
        @test_approx_eq bavg vanilla.bias
        for sample in data
            classify(sample.fvec, vanilla)
        end
    end
end

max_stars_repo_name: ninjin/allen

max_issues_repo_name: ninjin/allen

max_forks_repo_name: ninjin/allen


text:
import Base: resize!
resize!(M::ElasticMatrix, N::Integer) = resize!(M, size(M,1), div(N,size(M,1)))

"""
    hermitianize!(x)

Extract the Hermitian part of `x`.
# Arguments
* `x`: some matrix.
"""
function hermitianize!(x::AbstractMatrix)
    ishermitian(x) && return x
    n = LinearAlgebra.checksquare(x);
    @inbounds @views for j = 1:n, i = 1:j
        x[i,j] = (x[i,j] + conj(x[j,i])) / 2.
        x[j,i] = conj(x[i,j]);
    end
    return x
end

"""
    compute_tstops(circuit, t0=0.0)

Compute all callback and gate switching times of the circuit.
"""
function compute_tstops(circuit::AbstractArray{G}, t0::Float64 = 0.0) where G<:AbstractGate
    # Times at which a PCA step is performed
    tstops  = Float64[t0]
    # Times at which the current gate changes (change of Hamiltonian)
    tswitch = Float64[]

    t = t0
    for gate in circuit
        push!(tstops, [t:gate.dt:t+gate.T;][2:end]...)
        push!(tswitch, t + gate.T)
        t += gate.T
    end

    return tstops[1:end], tswitch[1:end-1]
end

"""
    time_evolution(L0, circuit, J; algo=Tsit5(), ϵ=1e-5, tol=1e-9, fout=missing)

Compute the full time evolution of an initial state `L0` through a `circuit` under
the action of dissipators `J`.
"""
function time_evolution(L0::ElasticMatrix, circuit::AbstractArray, J::AbstractArray; algo=Tsit5(), ϵ=1e-5, tol=1e-9, fout=missing)
    N = size(L0,1)          # Size of the full Hilbert space
    M = size(L0,2)          # Size of the corner
    tstops, tswitch = compute_tstops(circuit)
    i_g = 1                 # Current gate index
    h = circuit[i_g].dt     # Current timestep

    # "kraus" operators (K_i = sqrt(h) * k_i)
    ks = eltype(J) <: AbstractOperator ? map(j -> j.data, J) : J

    # Compute the non-Hermitian term of the Hamiltonian
    nh_term = zero(first(ks))
    for k in ks
        nh_term -= 0.5 .* k'k
    end
    # Initialise the non-Hermitian Hamiltonian
    δk0 = -im.*hamiltonian(circuit[i_g])::typeof(nh_term) + nh_term # Generator of the pseudounitary dynamics

    # results from fout
    res_fout = ismissing(fout) ? nothing : [fout(0.0, L0)]

    # Caches
    KL = ElasticArray{ComplexF64}(undef, N, (length(ks)+1) * M)
    KL.= zero(eltype(KL))
    # Precompute the dissipative Kraus trajectories at t+h
    for i in 1:length(ks), j in 1:M
        KL[:,i*M + j] .= sqrt(h) .* ks[i]*L0[:,j]
    end

    # Pseudounitary time evolution
    function f!(dL, L, p, t)
        mul!(dL, δk0, L)
    end

    # Function called at each callback, performs PCA over the Kraus trajectories
    # generated during the previous time step
    function affect!(integrator)
        if integrator.t in tswitch
            # Switch to next quantum gate
            i_g += 1
            # Adapt timestep
            h = circuit[i_g].dt
            # Adapt the Hamiltonian to match the new gate
            δk0 = -im.*hamiltonian(circuit[i_g])::typeof(δk0) + nh_term
        end
        # Add missing pseudounitary trajectories
        @inbounds @views KL[:,1:size(integrator.u,2)] .= integrator.u
        # Reinitialise the integrator state
        integrator.u .= zero(ComplexF64)

        # Build the low-rank representation of the density matrix
        ρ = KL'KL
        hermitianize!(ρ)
        _λs, cs = eigen(ρ)
        λs = real.(_λs)
        λs[λs .< 0.0] .= 0.0
        # Regularise
        λs ./= sum(λs)

        perm = sortperm(λs; rev=true)
        # New corner dimension
        M = findfirst(trace->1.0-trace < ϵ, accumulate(+,λs[perm]))
        # map to kept columns
        σ = perm[1:M]

        # Resize the integrator state vector and caches to match the new corner dimension
        resize!(integrator, N*M)

        # Initialise the integrator state for the timestep to come
        for (i, σi) in enumerate(σ)
            # Non-allocating equivalent to :
            # integrator.u[:,i] .= sqrt(λs[σi]) .* normalize(droptol!(sparse(KL * cs[:,σi]),tol))
            v = @view integrator.u[:,i]
            mul!(v, KL, view(cs, :, σi), one(eltype(v)), zero(eltype(v)))
            v[abs.(v) .< tol] .= zero(eltype(v))
            normalize!(v)
            v .*= sqrt(λs[σi])
        end

        # Apply fout to the current integrator state
        !ismissing(fout) && push!(res_fout, fout(integrator.t, copy(integrator.u)))

        # Precompute the dissipative Kraus trajectories at t+h
        resize!(KL, N, (length(ks)+1) * M)
        for i in 1:length(ks), j in 1:M
            KL[:,i*M + j] .= sqrt(h) .* ks[i]*integrator.u[:,j]
        end

        integrator
    end

    # Callback that deals with the periodic call of the PCA procedure
    cb = PresetTimeCallback(tstops[2:end],affect!)

    # Build the ODE problem
    pb = ODEProblem(f!, L0, extrema(tstops))

    L = solve(pb, algo, callback=cb, save_on=false, save_end=true).u[2]

    if ismissing(fout)
        return L
    else
        L, res_fout
    end
end

time_evolution(L0::StateVector,    circuit::AbstractArray, J::AbstractArray; kwargs...) = time_evolution(L0.data, circuit, J; kwargs...)
time_evolution(L0::AbstractVector, circuit::AbstractArray, J::AbstractArray; kwargs...) = time_evolution(reshape(L0,length(L0),1), circuit, J; kwargs...)
time_evolution(L0::AbstractMatrix, circuit::AbstractArray, J::AbstractArray; kwargs...) = time_evolution(ElasticMatrix(L0), circuit, J; kwargs...)


text:
###  R E F E R E N C E S
# Scope is to recognise and manage local and global references to positions or ranges

## C E L L   P O S I T I O N S
const colsyntax = "([A-Z]+)"
const rowsyntax = "([1-9][0-9]*)"
const cellsyntax = colsyntax*rowsyntax
const cellregex = Regex("^"*cellsyntax*"\$")

# Recognise a cell in the same sheet
function _LocalPosition(r::AbstractString)  # accept "A1"-syntax, convert to "A",1 syntax
                                            # returns Union{LocalPosition,Nothing}
    m = match(cellregex,r)
    isnothing(m) && return nothing
    colStr = string(m.captures[1])
    rowNum = parse(UInt,m.captures[2])
    LocalPosition(colStr,rowNum)
end
function LocalPosition(r::AbstractString)   # accept "A1"-syntax, convert to "A",1 syntax
                                            # returns LocalPosition or throws error
    s = _LocalPosition(r)
    @assert !isnothing(s) ("That's not a valid local reference: " * r)
    s
end
function colStrAsNum(c::AbstractString) # accept "A", convert to 1
    colsbase26 = map(x -> Int(x)-64,collect(c))
    col = zero(UInt)
    l = length(colsbase26)
    for i in 1:l
        col += colsbase26[i]*26^(l-i)
    end
    col
end
function LocalPosition(colStr::AbstractString,rowNum::Number) # accept "A",1-syntax, convert to 1,1-syntax
    LocalPosition(rowNum,colStrAsNum(colStr))
end

# arithmetic of local positions
Base.isless(p::LocalPosition,q::LocalPosition) = (p.row < q.row) && (p.column < q.column)
Base.:-(p::LocalPosition,q::LocalPosition) = LocalPosition(p.row-q.row, p.column-q.column)
Base.:+(p::LocalPosition,q::LocalPosition) = LocalPosition(p.row+q.row, p.column+q.column)
Base.maximum(pp::Array{LocalPosition}) = LocalPosition(maximum(getfield.(pp,:row)),maximum(getfield.(pp,:column)))
Base.minimum(pp::Array{LocalPosition}) = LocalPosition(minimum(getfield.(pp,:row)),minimum(getfield.(pp,:column)))

# display functions
function colNumAsStr(colNum)
    if colNum > 26
        colNumAsStr(div(colNum,26)) * colNumAsStr(rem(colNum,26))
    else
        string(Char(64+colNum))
    end
end
Base.string(p::LocalPosition) = colNumAsStr(p.column) * string(p.row)
Base.show(io::IO, p::LocalPosition) = Base.print(io, string(p))

## C E L L   R A N G E S
const rangeregex = Regex("([1-9A-Z][0-9A-Z!]*:[1-9A-Z][0-9A-Z]*)")

# cell ranges
const cellrangeregex = Regex("^"*cellsyntax*":"*cellsyntax*"\$")
CellRange(p::AbstractString,q::AbstractString) = CellRange(LocalPosition(p),LocalPosition(q))
Base.:(:)(p::LocalPosition,q::LocalPosition) = CellRange(p,q)
Base.string(r::CellRange) = string(r.start) * ":" * string(r.stop)
Base.show(io::IO, r::CellRange) = Base.print(io, string(r))

function Base.collect(R::CellRange)
    m = min(R.start,R.stop)
    r = m.row
    c = m.column
    M = max(R.start,R.stop)
    delta = M - m
    I = delta.row + 1
    J = delta.column + 1
    pp = Array{LocalPosition}(undef,I,J) # empty array
    for i in 1:I, j in 1:J
        pp[i,j] = LocalPosition(r+i-1,c+j-1)
    end
    pp
end
Base.collect(p::LocalPosition) = [p]
expanded_refs(rr) = vcat(collect.(rr)...)

# larger ranges --> not really there yet
const columnrangeregex = Regex("^"*colsyntax*":"*colsyntax*"\$")
ColumnRange(a::AbstractString,b::AbstractString) = ColumnRange(colStrAsNum(a),colStrAsNum(b))
Base.string(r::ColumnRange) = colNumAsStr(r.start) * ":" * colNumAsStr(r.stop)
Base.show(io::IO, r::ColumnRange) = Base.print(io, string(r))

const rowrangeregex = Regex("^"*rowsyntax*":"*rowsyntax*"\$")
RowRange(a::AbstractString,b::AbstractString) = RowRange(parse(UInt,a),parse(UInt,b))
Base.string(r::RowRange) = string(r.start) * ":" * string(r.stop)
Base.show(io::IO, r::RowRange) = Base.print(io, string(r))

## L O C A L   R E F E R E N C E S
const localsyntax = "([1-9A-Z][0-9A-Z:]*)"

function _LocalRange(s::AbstractString) # returns Union{LocalRangeType,Nothing}
    m = match(cellrangeregex,s)
    isnothing(m) || return CellRange(m.captures[1]*m.captures[2],m.captures[3]*m.captures[4])
    m = match(columnrangeregex,s)
    isnothing(m) || return ColumnRange(m.captures[1],m.captures[2])
    m = match(rowrangeregex,s)
    isnothing(m) || return RowRange(m.captures[1],m.captures[2])
    nothing
end
function LocalRange(s::AbstractString) # returns LocalRangeType or throws error
    r = _LocalRange(s)
    @assert !isnothing(r) ("That's not a valid local range: " * s)
    r
end

function _LocalRef(s::AbstractString) # returns Union{LocalRef,Nothing}
    t = _LocalPosition(s)
    isnothing(t) || return t
    u = _LocalRange(s)
    isnothing(u) || return u
    nothing
end
function LocalRef(s::AbstractString) # returns LocalRef or throws error
    r = _LocalRef(s)
    @assert !isnothing(r) ("That's not a valid local ref: " * s)
    r
end

## G L O B A L   P O S I T I O N
const sheetsyntax = "([a-zA-Z_][0-9A-Za-z]*)"
const sheetregex = Regex("^"*sheetsyntax*"\$")

const globalposregex = Regex("^"*sheetsyntax * "!" * cellsyntax*"\$")
const globalrefregex = Regex("^"*sheetsyntax * "!" * localsyntax*"\$")
const refregex = Regex("^"*sheetsyntax)

# "Sheet1", "A1"-syntax
GlobalPosition(s::AbstractString,p::AbstractString) = GlobalPosition(s, LocalPosition(p))
function _GlobalPosition(s::AbstractString)
    m = match(globalposregex,s)
    isnothing(m) || return GlobalPosition(m.captures[1],LocalPosition(m.captures[2] * m.captures[3]))
    nothing
end
function GlobalPosition(s::AbstractString)
    t = _GlobalPosition(s)
    @assert !isnothing(t) ("That's not a valid global position: " * s)
    t
end
LocalPosition(gp::GlobalPosition) = gp.lref

## G L O B A L   R E F E R E N C E
GlobalRef(s::AbstractString, lref::LocalPosition) = GlobalPosition(s,lref)

function Base.:(:)(p::GlobalPosition,q::GlobalPosition)
    @assert (p.sheet == q.sheet) ("Cannot range between cells of two different sheets!: " * string(p) * " and " * string(q))
    GlobalRange(p.sheet, CellRange(p.lref,q.lref))
end
Base.collect(R::GlobalRange) = GlobalPosition.(R.sheet, collect(R.lref))
GlobalRef(s::AbstractString, lref::LocalRangeType) = GlobalRange(s,lref)

GlobalRef(s::AbstractString, lref::AbstractString) = GlobalRef(s, LocalRef(lref))
function _GlobalRef(s::AbstractString)
    m = match(globalrefregex,s)
    isnothing(m) || return GlobalRef(m.captures[1], LocalRef(m.captures[2]))
    nothing
end
function GlobalRef(s::AbstractString)
    t = _GlobalRef(s)
    @assert !isnothing(t) ("That's not a valid global reference: " * s)
    t
end
function GlobalisedRef(s::AbstractString, default_Sheet::String)
    t = _GlobalRef(s)
    isnothing(t) || return t
    u = LocalRef(s)
    isnothing(u) || return GlobalRef(default_Sheet,u)
    error("That's not a reference: " * s)
end

# display functions
Base.string(r::GlobalRef) = r.sheet * "!" * string(r.lref)
Base.show(io::IO, r::GlobalRef) = Base.print(io, string(r))

## P A R S I N G
# make sure that ranges are correctly grabbed by the Julia AST parser, by adding brackets around them
bracketedRanges(f::AbstractString) = replace(f, rangeregex => s"(\1)")

"""
    extractMatchFromSymbol(symbol)

Identifies and returns the string representing a reference to a position or range in the `symbol`
"""
function extractMatchFromSymbol(s::Symbol)
    t = string(s)
    m = match(globalrefregex,t)
    isnothing(m) || return m.match
    m = match(cellregex,t)
    isnothing(m) || return m.match
    nothing
end

function exprIsRange(e::Expr)
    (e.head == :call) &&
    (length(e.args) >= 3) &&
    (e.args[1] == :(:)) &&
    isa(e.args[2],Symbol) &&
    isa(e.args[3],Symbol)
end
"""
    extractRefsFromExpression(expression)

Identifies and returns the set of strings representing references to positions and ranges in the `expression`
"""
function extractRefsFromExpression(e::Expr)
    rr = String[]

    exprIsRange(e) && return [string(e)]

    m = match(cellregex,string(e.head))
    isnothing(m) || push!(rr,m.match)

    for (index,arg) in enumerate(e.args)
        if isa(arg,Symbol)
            m = match(cellregex,string(arg))
            isnothing(m) || push!(rr,m.match)
        elseif isa(arg,Expr)
            push!(rr,extractRefsFromExpression(arg)...)
        end
    end

    unique(rr)
end