text:
module VAWTools
using Parameters
using Printf, SparseArrays, Statistics

# General tools, maybe applicable more widely.
export read_agr, write_agr, read_xyn, inpoly, AGridded, Gridded, Gridded1d, Traj,
    smooth_vector, absslope, gradient3by3,
    downsample, split_traj!, boxcar, boxcar_matrix, apply_boxcar_matrix, bin_grid, piecewiselinear, split_poly,
    transform_proj, max_smooth, min_smooth

include("smoothing-functions.jl")

import Base: ==, size, length, step, +, -, *

# const spatialorder = "yx"  # order of imagestorage, https://github.com/timholy/Images.jl#storage-order-and-changing-the-representation-of-images

abstract type AGridded{T} end
"Check whether an error-field is defined."
haserror(g::AGridded) = size(g.err)!=(0,0)
hasproj(g::AGridded) = g.proj!=""
size(g::AGridded) = size(g.v)
length(g::AGridded) = length(g.v)
step(g::AGridded) = step(g.x)

function Base.show(io::IO, g::AGridded)
    show(io, typeof(g))
    print(io, " ")
    # if haserror(g)
    #     println(io, " holding values and errors on grid:")
    # else
    #     println(io, " holding values on grid:")
    # end
    show(io, g.x);
    if isdefined(g, :y)
        print(io, " "); show(io, g.y)
    end
end

"""
Gridded holds 2D data (and its error) on a regular grid.

- rows of v correspond to x-dir, columns to ydir
- be sure to be clear whether x corresponds to cell centers or
  boundaries.

TODO:
- hold NA value?
- default of err, see https://github.com/JuliaLang/julia/issues/23926
"""
@with_kw struct Gridded{T} <: AGridded{T}
    x::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
    y::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
    v::Matrix{T}  # values
    @assert size(v)==(length(x), length(y))
    err::Matrix{T}=Matrix{eltype(v)}(0,0)  # error of values
    midpoint::Bool=true  # if true, the (x,y) is cell midpoint, otherwise lower-left corner
    proj::String="" # Proj4 string
end
Gridded(x,y,v::Matrix{T}) where T = Gridded{T}(x,y,v,Matrix{T}(undef,0,0),true,"")
Gridded(x,y,v::Matrix{T},err) where T = Gridded{T}(x,y,v,err,true,"")

"""
    downsample(g::Gridded, step, start=1, average=true, averagemask=all_points)

Makes a `Gridded` smaller by downsampling.  If averaging (by default),
the window is chosen such that there is no overlap (for odd `step`) or
one cell overlap (for even `step`).

"""
function downsample(g::Gridded, step, start=1, average=true, averagemask=UniformArray{Bool,2}(true); stop=size(g.v))
    if start isa Number
        rx = start:step:stop[1]
        ry = start:step:stop[2]
    else
        rx = start[1]:step:stop[1]
        ry = start[2]:step:stop[2]
    end
    nx,ny = length(rx),length(ry)
    @assert nx>1
    @assert ny>1
    x = g.x[rx]
    y = g.y[ry]
    midpoint = g.midpoint
    # Window: no overlap on odd step sizes
    window = average ? Int(floor((step)/2)) : 1
    vnew = _downsample_boxcar(g.v,rx,ry,window,averagemask)

    errnew = haserror(g) ? _downsample_boxcar(g.err,rx,ry,window,averagemask) : g.err
    Gridded(x,y,vnew,errnew,g.midpoint,g.proj)
end
function _downsample_boxcar(v,rx,ry,window,averagemask)
    # make type which allows for averaging
    T = typeof(one(eltype(v))/2)
    vnew = zeros(T,length(rx), length(ry))
    if window==0
        for (j,rj) in enumerate(ry)
            for (i,ri) in enumerate(rx)
                vnew[i,j]=v[ri,rj]
            end
        end
        return vnew
    end
    # Follows boxcar filter below
    for (j,rj) in enumerate(ry)
        for (i,ri) in enumerate(rx)
            # keep NaNs or points outside averagemask as they are:
            if isnan(v[ri,rj]) || !averagemask[ri,rj]
                vnew[i,j] = v[ri,rj]
                continue
            end
            # other points make an average:
            v0 = zero(T)
            n = 0
            for jj=max(1,rj-window):min(size(v,2),rj+window)
                for ii=max(1,ri-window):min(size(v,1),ri+window)
                    tmp = v[ii,jj]
                    if !isnan(tmp) && averagemask[ii,jj]
                        v0 += tmp
                        n+=1
                    end
                end
            end
            vnew[i,j] = v0/n
        end
    end
    return vnew
end

import Interpolations
"""
    upsample(g::Gridded, x, y)

Upsample to match x, y.
"""
function upsample(g::Gridded, x, y)
    gi  = Interpolations.interpolate((g.x, g.y), g.v, Interpolations.Gridded(Interpolations.Linear()) )
    if haserror(g)
        gie  = Interpolations.interpolate((g.x, g.y), g.err, Interpolations.Gridded(Interpolations.Linear()) )
        return Gridded(x, y, [gi[xx,yy] for xx=x, yy=y],
                [gie[xx,yy] for xx=x, yy=y])
    else
        return Gridded(x, y, [gi[xx,yy] for xx=x, yy=y])
    end
end

function (+)(g1::Gridded, g2::Gridded)
    @assert g1.midpoint==g2.midpoint
    @assert g1.proj==g2.proj
    Gridded(g1.x, g1.y, g1.v.+g2.v, g1.err.+g2.err, g1.midpoint, g1.proj)
end
(*)(n::Number, g::Gridded) = Gridded(g.x, g.y, g.v*n, g.err*n, g.midpoint, g.proj)
(*)(g::Gridded, n::Number) = n*g

"""
    fillgaps(g::Gridded, di=1, fillval=NaN)

Fills gaps on a gridded array by averaging over a distance `di` (number of cells cells).

Note,
- if gridded-dataset is surrounded by fillvalues then there might be
  odd edge effects up to `di` number of cells.
- the test against `fillval` is done with `isequal` (to allow NaNs)
"""
function fillgaps(g::Gridded, di=1, fillval=NaN)
    sz1,sz2 = size(g.v)
    gg = deepcopy(g)
    vv = g.v
    @inbounds for (j,y)=enumerate(g.y), (i,x)=enumerate(g.x)
        # fill in FILL-gaps inside the glacier:
        if isequal(g.v[i,j], fillval)
            val = 0.0
            n = 0
            for jj=max(1,j-di):min(sz2,j+di)
                for ii=max(1,i-di):min(sz1,i+di)
                    tmp = vv[ii,jj]
                    if !isnan(tmp)
                        val += tmp
                        n+=1
                    end
                end
            end
            if n>0
                gg.v[i,j] = val/n
            else
                # leave unfilled hole
            end
        end
    end
    return gg
end

"""
Holds 1D fields, e.g. elevation band data.

Be sure to be clear whether x corresponds to cell centers or boundaries.
"""
@with_kw struct Gridded1d{T} <:  AGridded{T}
    x::StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}
    v::Vector{T} # values
    @assert size(v)==(length(x), )
    err::Vector{T}=Vector{eltype(v)}(undef, 0) # error of values
    @assert size(v)==size(err) || size(err)==(0,)
    midpoint::Bool=true  # if true, the (x) is cell midpoint, otherwise lower corner
    proj::String="" # Proj4 string
end
Gridded1d(x,v::Vector{T}) where T = Gridded1d{T}(x,v,T[],true,"")
Gridded1d(x,v::Vector{T},err) where T = Gridded1d{T}(x,v,err,true,"")

# function (+)(g1::Gridded1d, g2::Gridded1d)
#     @assert g1.midpoint==g2.midpoint
#     @assert g1.proj==g2.proj
#     Gridded1d(g1.x, g1.v+g2.v, g1.err+g2.err, g1.midpoint, g1.proj)
# end
(*)(n::Number, g::Gridded1d) = Gridded1d(g.x, g.v*n, g.err*n, g.midpoint, g.proj)
(*)(g::Gridded1d, n::Number) = n*g

"""
Holds values on a trajectory or generally unstructured data.

Fields:
- x,y -- coordinates
- v::Vector{T} -- values
- err::Vector{T} -- errors (or second set of values)
- splits::Vector{UnitRange{Int}} -- If the trajectory consists of several sub-trajectories
- proj::String a string interpretable by Proj4

There are constructors to leave off v, err, splits and proj.

TODO: having the `splits` is probably worse than just using a list of trajectories.
"""
@with_kw mutable struct Traj{T}
    x::Vector{Float64}
    y::Vector{Float64}
    @assert length(x)==length(y)
    v::Vector{T}=Void[]
    @assert length(v)==length(y) || length(v)==0
    err::Vector{T}=Vector{eltype(v)}(undef, 0)
    @assert length(err)==length(y) || length(err)==0
    splits::Vector{UnitRange{Int}}=UnitRange{Int}[1:length(x)]
    proj::String="" # Proj4 string
end
Traj(x,y) = Traj{Nothing}(x,y,Nothing[],Nothing[],UnitRange{Int}[1:length(x)],"")
Traj(x,y,v::AV) where AV<:AbstractVector{T} where T = Traj{T}(x,y,v,T[],UnitRange{Int}[1:length(x)],"")
Traj(x,y,v::AV,err) where AV<:AbstractVector{T} where T = Traj{T}(x,y,v,err,UnitRange{Int}[1:length(x)],"")
Traj(x,y,v::AV,err,splits) where AV<:AbstractVector{T} where T = Traj{T}(x,y,v,err,splits,"")
function Traj(t::Traj, inds::Vector{Int})
    x,y,v,err,splits = similar(t.x,0), similar(t.y,0), similar(t.v,0), similar(t.err,0), similar(t.splits,0)
    s1 = length(x)
    for i in inds
        append!(x, t.x[t.splits[i]])
        append!(y, t.y[t.splits[i]])
        if length(t.v)>0
            append!(v, t.v[t.splits[i]])
        end
        if length(t.err)>0
            append!(err, t.err[t.splits[i]])
        end
        push!(splits, s1+1:length(x))
        s1 = length(x)
    end
    Traj(x,y,v,err,splits,t.proj)
end

hasvalues(g::Traj) = length(g.v)!=0
haserror(g::Traj) = length(g.err)!=0
hasproj(g::Traj) = g.proj!=""
length(g::Traj) = length(g.x)

function distances(t::Traj)
    dist = zeros(length(t)-1)
    @inbounds for i=1:length(dist)
        dist[i] = sqrt( (t.x[i+1]-t.x[i])^2 + (t.y[i+1]-t.y[i])^2)
    end
    dist
end

"""
    downsample(t::Traj{T}, window::Float64, average=true) where T
    downsample(t::Traj{T}, step::Int)

Downsample a
"""
function downsample(t::Traj, step::Int)
    @assert length(t.splits)==
    return Traj(t.x[1:step:end], t.y[1:step:end], t.v[1:step:end], t.err[1:step:end], [1:length(t.err[1:step:end])], t.proj)
end
function downsample(t::Traj{T}, window::Float64, average=true) where T
    dist = distances(t)
    x, y = Float64[], Float64[]
    v = T[]
    err = T[]
    splits = similar(t.splits, 0)
    ttt = time()
    for s in t.splits
        i = s[1]
        s1 = length(x)+1
        while i<s[end]
            ##
            d = 0.0
            n = 1
            val = t.v[i]
            eval = t.err[i]
            if average
                # before point
                for j=i-1:-1:s[1]
                    d += dist[j]
                    if d<window/2
                        n += 1
                        val += t.v[j]
                        eval += t.err[j]
                    else
                        break
                    end
                end
                # after point
                d = 0.0
                for j=i+1:s[end]
                    d += dist[j-1]
                    if d<window/2
                        n += 1
                        val += t.v[j]
                        eval += t.err[j]
                    else
                        break
                    end
                end
            end
            push!(x, t.x[i])
            push!(y, t.y[i])
            push!(v, val/n)
            push!(err, eval/n)

            ## Find next i
            d = 0.0
            for j=i+1:s[end] # returns also s[end]
                d += dist[j-1]
                if d>=window/2 || j==s[end]
                    i = j
                    break
                end
            end
        end
        push!(splits, s1:length(x))
    end
    return Traj(x, y, v, err, splits, t.proj)
end

function mindistance(point, x, y, inds=1:length(x))
    dist = Inf
    ind = 0
    @inbounds for i=inds
        d = (x[i]-point[1])^2 + (y[i]-point[2])^2
        if d<dist
            dist = d
            ind = i
        end
    end
    sqrt(dist), ind
end

"""
     sort_traj(tr)

Brute force sorting of a trajectory.  Ignores splits.
"""
function sort_traj(tr)
    x,y = tr.x, tr.y
    inds = collect(1:length(tr))
    todo = trues(inds)
    out = Int[]
    j = findmin(tr.x)[2]
    push!(out, j)
    todo[out[1]] = false

    @inbounds for i=2:length(tr)
        j = mindistance((tr.x[j], tr.y[j]), x, y, inds[todo])[2]
        push!(out, j)
        todo[j] = false
    end
    return Traj(tr.x[out], tr.y[out], tr.v[out], tr.err[out], [1:length(tr)], tr.proj)
end

"""
    issorted_traj(tr, inds=1:min(length(tr), 500))

Checks whether a trajectory is sorted.  Only checks inds.
"""
function issorted_traj(tr, inds=1:min(length(tr), 500))
    x,y = tr.x[inds], tr.y[inds]
    todo = trues(inds)
    out = Int[]
    j = findmin(tr.x[inds])[2]
    push!(out, j)
    todo[out[1]] = false

    @inbounds for i=inds[2:end]
        j = mindistance((tr.x[j], tr.y[j]), x, y, inds[todo])[2]
        push!(out, j)
        todo[j] = false
    end
    return Base.issorted(out)
end

"""
    split_traj!{T<:Traj}(t::T, dist)

Splits trajectory into several between points farther apart than `dist`.
"""
function split_traj!(t::T, dist) where T<:Traj
    ii = 1
    empty!(t.splits)
    for i=1:length(t.x)-1
        if (t.x[i]-t.x[i+1])^2+(t.y[i]-t.y[i+1])^2 > dist^2
            push!(t.splits,ii:i)
            ii = i+1
        end
    end
    push!(t.splits, ii:length(t.x))
    nothing
end

function Base.convert(::Type{Traj{T2}}, g::Traj{T1}) where {T2,T1}
    if typeof(g)==Traj{T2}
        return g
    end
    x = convert(Vector{T2}, g.x)
    y = convert(Vector{T2}, g.y)
    return Traj(x, y, convert(Vector{T2}, g.v), convert(Vector{T2}, g.err), g.splits, proj)
end

"""
Holds a rectangular grid, mirrors the ASCII-grid .agr, .grid, .asc, .bin files.

Note that the indexing into the value matrix ins awkward, see Ref
below.  In general use Gridded.

Ref: https://en.wikipedia.org/wiki/Esri_grid
"""
struct AGR{T} # Ascii GRid
    v::Matrix{T}    # values: orientation is awkward, as in the AGR file, see https://en.wikipedia.org/wiki/Esri_grid
    nc::Int64       # NCOLS TODO: remove those
    nr::Int64       # NROWS
    xll::Float64    # XLLCORNER
    yll::Float64    # YLLCORNER
    dx::Float64     # CELLSIZE
    NA::T           # NODATA
    hasutm::Bool    # Matthias sometimes abuses the NODATA_value field as UTM-zone field
    extra_header::Vector{Float32} # the .bin files have space for
                                  # extra information in the header
    function AGR{T}(va, nc, nr, xll, yll, dx, NA,
                    hasutm=false, extra_header=zeros(Float32,6)) where T
        @assert size(va)==(nr,nc)
        @assert dx>=0
        new{T}(va, nc, nr, xll, yll, dx, NA, hasutm, extra_header)
    end
end
AGR(va::Matrix{T}, nc, nr, xll, yll, dx, NA, hasutm=false, extra_header=zeros(Float32,6)) where {T} = AGR{T}(va, nc, nr, xll, yll, dx, NA, hasutm, extra_header)

"""
Transform Gridded to AGR

Optional:
- NA_g: no-data value of g:Gridded, defaults to NaN
- NA_agr: no-data value of AGR output, defaults to NaN
- utmzone: if 0<utmzone<61, then above two are ignored and instead a file with UTM_ZONE is written.
- write_err: write g.err instead of g.v

"""
function AGR(g::Gridded{T}; NA_g=convert(T,NaN), NA_agr=convert(T,NaN), write_err=false, utmzone=0) where T
    v = write_err ? g.err : g.v # whether to write the values or errors into the ascii grid
    if 0<utmzone<61
        hasutm = true
    else
        hasutm = false
        if !isequal(NA_g,NA_agr)
            v = deepcopy(v)
            for i=eachindex(v)
                if isequal(v[i],NA_g)
                    v[i] =  NA_agr
                end
            end
        end
    end

    v = rotl90(v)
    nc,nr = size(g)
    dx = convert(T, step(g.x))
    if g.midpoint
        xll = g.x[1]-dx/2
        yll = g.y[1]-dx/2
    else
        xll = g.x[1]
        yll = g.y[1]
    end
    if hasutm
        AGR(v, nc, nr, xll, yll, dx, utmzone, hasutm)
    else
        AGR(v, nc, nr, xll, yll, dx, NA_agr)
    end
end

# xs(g::AGR) = linspace(g.xll+g.dx/2, nr*g.dx, nr)
# ys(g::AGR) = linspace(g.yll+g.dx/2, nc*g.dx, nc)

size(g::AGR) = size(g.v)
length(g::AGR) = length(g.v)
function ==(g1::AGR,g2::AGR)
    g1.v==g2.v &&
    g1.nc==g2.nc &&
    g1.nr==g2.nr &&
    g1.xll==g2.xll &&
    g1.yll==g2.yll &&
    g1.dx==g2.dx &&
    g1.NA==g2.NA &&
    g1.extra_header==g2.extra_header
end

# function my_rotr90(A::AbstractMatrix)
#     ind1, ind2 = indices(A)
#     B = similar(A, (ind2,ind1))
#     m = first(ind1)+last(ind1)
#     @inbounds for i=ind1, j=indices(A,2) #added inbounds
#         B[j,m-i] = A[i,j]
#     end
#     return B
# end

"""
Construct Gridded from AGR.  Replace AGR-no-value NA with NaN.  By
default assumes AGR-no-value==NaN.  Moves (x,y)-coords to midpoints
"""
function Gridded(agr::AGR{T}; NA=convert(T,NaN)) where T
    #    v = my_rotr90(agr.v)
    v = rotr90(agr.v)
    if agr.hasutm
        proj = "+proj=utm +zone=$(Int(agr.NA)) +datum=WSG84"
    else
        # only swap NA if it has a FILL value:
        if !isequal(NA, agr.NA)
            for i=eachindex(v)
                if isequal(v[i], agr.NA); v[i] = NA end
            end
        end
        proj = ""
    end
    Gridded(range(agr.xll+agr.dx/2, step=agr.dx, length=agr.nc),
            range(agr.yll+agr.dx/2, step=agr.dx, length=agr.nr),
            v,
            zeros(T,0,0),
            true,
            proj)
end

# File readers
##############

"""
Check if .bin file nor not.  With IO types this may not work 100%.
"""
isbin_file(fn::AbstractString) = splitext(fn)[2]==".bin"
function isbin_file(fn::IO)
    if :name in fieldnames(typeof(fn))
        f = fn.name
    else # this should catch TranscodingStreams.jl
        f = fn.stream.name
    end
    f = strip(f, ['<','>']) # often looks like "<file asdf.ext>"
    f = strip(f)
    isbin_file(f)
end

"""
    read_agr(fl, T=Float32; NA=convert(T,NaN))

Read a Ascii grid https://en.wikipedia.org/wiki/Esri_grid.  Can also
read from .gz compressed files, if `import CodecZlib`.

In:
- fn -- file name

Optional keywords:
- T -- type of output array. Defaults to Float32.
- NA -- replace the fill value with this. Defaults to use NaN.

Out: Gridded
"""
function read_agr(fl, T=Float32; NA=convert(T,NaN))
    agr = _read_agr(fl, T)
    @unpack v, hasutm = agr
    NA_old = agr.NA
    # replace missing value with something else
    if !hasutm && NA_old!=NA
        _refill!(v, NA_old, NA)
    end
    Gridded(agr)::Gridded{T}
end

function _read_agr(fl::AbstractString, T=Float32)
    if !isfile(fl)
        error("File $fl does not exist!")
    end
    if endswith(fl, ".gz")
        @eval import CodecZlib
        io = CodecZlib.GzipDecompressorStream(open(fl))
    else
        io = open(fl, "r")
    end
    out = _read_agr(io, T)
    close(io)
    return out
end

function _read_agr(io::IO, T=Float32)
    if isbin_file(io)
        toT = (io,T) -> convert(T, read(io, Float32))
    else
        toT = (io,T) -> parse(T, split(readline(io))[2])
    end

    local va::Matrix{T}
    # read header
    nc = toT(io, Int)
    nr = toT(io, Int)
    xll = toT(io,Float64)
    yll = toT(io,Float64)
    dx = toT(io,Float64)
    # Matthias sometimes abuses the NODATA_value field as UTM-zone field
    # in non-binary grids:
    if isbin_file(io)
        fill = toT(io,T)
        hasutm = false
    else
        prop, val = split(readline(io))
        if lowercase(prop)=="nodata_value"
            hasutm = false
        elseif lowercase(prop)=="utm_zone"
            hasutm = true
        else
            error("Unrecognized last header field: $prop")
        end
        fill = parse(T,val)
    end

    if isbin_file(io)
        # read extra header
        extra_header = Array{Float32}(undef, 6)
        read!(io, extra_header)
        # read values
        va = Array{Float32}(undef, nc, nr)
        va = permutedims(read!(io, va))
        if eltype(va)!=T
            error("Not implemented yet")
        end
        if !eof(io)
            @warn("End-of-file was not reached!")
        end
    else
        va = Array{T}(undef, nr, nc)
        # no extra header for ascii .agr
        extra_header = zeros(Float32, 6)
        # read values
        tmp = split(read(io, String))
        if length(tmp)!=nc*nr
            error("Something's wrong with the file/stream $io.  It should contain $(nc*nr) values but has $(length(tmp))")
        end
        for i=1:nr, j=1:nc
            va[i,j] = parse(T, tmp[(i-1)*nc + j])
        end
    end
    # make a AGR data structure
    AGR(va, nc, nr, xll, yll, dx, fill, hasutm, extra_header)
end
# A function barrier is needed, thus this helper function:
function _refill!(a::Array, oldfill, newfill)
    @inbounds for i in eachindex(a)
        if a[i]==oldfill
            a[i] = newfill
        end
    end
    newfill
end

"""
    get_utm_asciigrid(fl_or_io)

Get the UTM zone of an ASCII-grid file (a unofficial format change used by Matthias).

If there is no zone return "".
"""
function get_utm_asciigrid(fl)
    open(fl) do io
        get_utm_asciigrid(io)
    end
end
function get_utm_asciigrid(io::IO)
    T = Float32
    if isbin_file(io)
        return ""
        # error("Ascii-bin files cannot have UTM")
    else
        toT = (io,T) -> parse(T, split(readline(io))[2])
    end

    # read header
    _ = toT(io, Int)
    _ = toT(io, Int)
    _ = toT(io,T)
    _ = toT(io,T)
    _ = toT(io,T)
    # Matthias sometimes abuses the NODATA_value field as UTM-zone field
    # in non-binary grids:
    prop, val = split(readline(io))
    if lowercase(prop)=="nodata_value"
        return ""
        # error("No field UTM_ZONE found")
    end
    utm = parse(Int,val)
    return "+proj=utm +zone=$(utm) +datum=WGS84"
end

"""
    write_agr(g::Gridded{T_}, fn::AbstractString; T=T_, NA_g=convert(T_,NaN), NA_agr=convert(T,NaN), write_err=false) where T_

Write Ascii grid.

Output format is determined by the extension:
- .bin -- binary
- .agr -- ASCII
- .grid or .asc are renamed to .agr

Always writes .bin in Float32

Optional:
- NA_g: no-data value of input g::Gridded, defaults to NaN
- NA_agr: no-data value of AGR-write output, defaults to NaN
- write_err: write g.err instead of g.v
- T: choose output type: Int or Float (only for ASCII files)
"""
function write_agr(g::Gridded{T_}, fn::AbstractString; T=T_, NA_g=convert(T_,NaN), NA_agr=convert(T,NaN), write_err=false) where T_
    write_agr(AGR(g, NA_g=NA_g, NA_agr=NA_agr, write_err=write_err), fn, T=T)
end

# setting NA will transform the NA value to that
function write_agr(g::AGR{T_}, fn::AbstractString; NA=nothing, T=T_) where T_
    if !isbin_file(fn)
        ext = splitext(fn)[2]
        if ext==".grid" || ext==".asc"
            # println("Changing extension to .agr")
            # ext = ".agr"
        elseif ext!=".agr"
            error("unsupported extension")
        end
    end
    if isbin_file(fn)
        convInt(x) = convert(Float32, x)
        convT(x) = convert(Float32, x)
        convMT = convT
    else # conversion to ASCII
        convInt = (x) -> @sprintf("%d\n", x)
        if T<:AbstractFloat
            convT = (x) ->  @sprintf("%f\n", x)
            convMT = (x) ->  @sprintf("%f  ", x)
        elseif T<:Integer
            convT = (x) ->  @sprintf("%d\n", x)
            convMT = (x) ->  @sprintf("%d  ", x)
        end
    end
    open(fn, "w") do io
        # write header
        isbin_file(fn) || write(io, "ncols         ")
        write(io, convInt(g.nc))
        isbin_file(fn) || write(io, "nrows         ")
        write(io, convInt(g.nr))
        isbin_file(fn) || write(io, "xllcorner     ")
        write(io, convT(g.xll))
        isbin_file(fn) || write(io, "yllcorner     ")
        write(io, convT(g.yll))
        isbin_file(fn) || write(io, "cellsize      ")
        write(io, convT(g.dx))
        if !isbin_file(fn)
            if g.hasutm
                write(io, "UTM_ZONE      ")
            else
                write(io, "NODATA_value  ")
            end
        end
        # replace NA value if necessary
        if NA==nothing
            fill = g.NA
            va = g.v
        else
            fill = NA
            va = deepcopy(g.v)
            for i=1:length(va)
                if va[i]==g.NA
                    va[i] = fill
                end
            end
        end
        if g.hasutm
            write(io, convInt(fill))
        else
            write(io, convT(fill))
        end
        # write extra header if .bin:
        isbin_file(fn) && write(io, g.extra_header)
        # write values
        for i=1:g.nr, j=1:g.nc
            v = va[i,j]
            write(io, convMT(v))
        end
        return nothing
    end
end

"""
Read .xyn or .xyzn files which can contain several, joined polygons.

x              y           label
-2031744.122   833011.310  21
...

x              y           z   label
-2031744.122   833011.310  23.4  21
...

where the label marks the beginning or end of line (21 --start, 22 --
inside, 23 -- end).

There can be embedded polygons, with label sequence 23,20,21.  The
line with 20 is then dropped.

The polygon must be closed, i.e. first element equals last.

Set `hasz` if it also contains a z-coordinate, i.e. a .xyzn file.

Return:

- a list of [x,y(,z)] arrays, one for each polygon.
"""
function read_xyn(fn; hasz=false, fix=false)
    if !isfile(fn)
        error("File $fn cannot be found.")
    end
    x,y,z,l = open(fn, "r") do io
        x = Float64[]
        y = Float64[]
        z = Float64[]
        l = Int[]
        for ls in readlines(io)
            if hasz
                x_,y_,z_,l_ = split(ls)
            else
                x_,y_,l_ = split(ls)
            end
            push!(x, parse(Float64, x_))
            push!(y, parse(Float64, y_))
            if hasz
                push!(z, parse(Float64, z_))
            end
            push!(l, parse(Int, l_))
        end
        return x,y,z,l
    end
    # Split up into individual polygons
    n = length(x)
    out = Matrix{Float64}[]
    is = 1 # start index of one poly
    iend = -1  # end index of one poly
    i = 1
    drop20 = false
    while i<=n
        if l[i]!=21
            error("Malformed file $fn: Expected 21 on line $i.")
        end
        # start a new polygon
        i += 1
        while true  # go through one poly:
            if l[i]==23 # end of a poly
                if i!=n && l[i+1]==20 # drop the 20-line
                    drop20 = true
                end
                break
            elseif l[i]==22
                i+=1
            else
                error("shouldn't get here")
            end
        end
        iend = i
        if length(is:iend)<3
            error("Polygon needs at least three vertices")
        end
        # fill:
        if hasz
            push!(out, (hcat(x[is:iend],y[is:iend],z[is:iend]))' )
        else
            push!(out, (hcat(x[is:iend],y[is:iend]))' )
        end
        if x[is]!=x[iend] || y[is]!=y[iend]
            if fix
                out[end] = hcat(out[end], out[end][:,1])
            else
                error("Malformed file $fn: polygon $is:$iend not closed")
            end
        end

        is = iend+1
        i+=1
        if drop20
            is+=1
            i+=1
            drop20 = false
        end
    end
    return out
end
function read_xyz(fn)
    if !isfile(fn)
        error("File $fn cannot be found.")
    end
    x,y,z = open(fn, "r") do io
        x = Float64[]
        y = Float64[]
        z = Float64[]
        for ls in readlines(io)
            x_,y_,z_ = split(ls)
            push!(x, parse(Float64, x_))
            push!(y, parse(Float64, y_))
            push!(z, parse(Float64, z_))
        end
        return x,y,z
    end
    # # turn into grid
    tmp = diff(x)
    dx = median(tmp[tmp.>0])
    tmp = diff(y)
    dy = median(tmp[tmp.>0])
    xrange = extrema(x)
    yrange = extrema(y)
    xx = xrange[1]:dx:xrange[2]
    yy = yrange[1]:dy:yrange[2]
    zz = zeros(length(xx), length(yy)).*NaN
    for (X,Y,Z) in zip(x,y,z)
        i,j = findfirst(isequal(X), xx), findfirst(isequal(Y), yy)
        zz[i,j] = Z
    end
    return Gridded(xx,yy,zz)
end

"""
    concat_poly(mpoly::Vector)

Concatenates a split poly, as returned from read_xyn, into one. Also
returns indices where to split apart again.  The concatenated polygon
is fully connected, with an edge going back to the first point.  This
allows to use `inpoly`, at least if the inner polygons have different
orientation to the outer. Also note that the input and output polygons
are closed, i.e. last point == first point (this can be fixed by setting
the option `close_poly=true`)

Return:
- bigpoly -- with size==(2,n)
- splits -- ith poly has indices splits[i]:splits[i+1]-1
"""
function concat_poly(mpoly::Vector; close_poly=false)
    T = eltype(mpoly[1])
    # check that they are all closed
    if close_poly
        mpoly = deepcopy(mpoly)
    end
    for i=1:length(mpoly)
        if mpoly[i][:,1]!=mpoly[i][:,end]
            if close_poly
                mpoly[i] = hcat(mpoly[i], mpoly[i][:,1])
            else
                error("All input polys need to be closed.")
            end
        end
    end
    # total size is sum of sizes plus one extra point for all but the
    # first poly.
    totsize = mapreduce(x->size(x,2), +, mpoly) + length(mpoly) -1
    bigpoly = Array{T}(undef, size(mpoly[1],1), totsize)
    splits = Int[]
    is = 1
    for i=1:length(mpoly)
        push!(splits, is)
        bigpoly[:,is:is+size(mpoly[i],2)-1] = mpoly[i]
        if i==1
            is = is+size(mpoly[i],2)
        else
            # add extra point to connect back to mpoly[1][:,1]
            is = is+size(mpoly[i],2) + 1
            bigpoly[:,is-1] = mpoly[1][:,1]
        end
    end
    push!(splits, totsize+1)
    return bigpoly, splits
end

"Finds splitting points in a big-poly"
function find_poly_splits(bigpoly::Matrix)
    splits = Int[]
    p1 = bigpoly[:,1]
    for i=1:size(bigpoly,2)
        if p1==bigpoly[:,i]
            push!(splits,i+1)
        end
    end
    splits
end

"Split up concatenated polygon."
function split_poly(bigpoly::Matrix{T}, splits) where T
    out = Matrix{T}[]
    if length(splits)==0
        return out
    else
        push!(out, bigpoly[:,splits[1]:splits[1+1]-1])
    end
    for i=2:length(splits)-1
        # remove the extra point joining to the first poly
        push!(out, bigpoly[:,splits[i]:splits[i+1]-2])
    end
    out
end

import Proj4

import ArchGDAL
"""
    read_geotiff(fn::AbstractString, T=Float32; bandnr=1, NA=convert(T,NaN))

Reads a geotiff raster and put output
into a Gridded instance (including the Proj4 projection string).

Link: https://github.com/yeesian/ArchGDAL.jl/issues/68
"""
function read_geotiff(filepath::AbstractString, T=Float32; bandnr=1, NA=convert(T,NaN))
    AG = ArchGDAL
    out = AG.registerdrivers() do
        AG.read(filepath) do dataset
            band = AG.getband(dataset, bandnr)
            w, h = AG.width(band), AG.height(band)
            # scale, off = AG.getscale(band), AG.getoffset(band)
            na = AG.getnodatavalue(band)
            mat = convert(Matrix{T}, AG.read(band))
            mat[mat.==na] = NA
            gt = AG.getgeotransform(dataset)
            dx, dy = gt[2], -gt[end]
            x0 = gt[1] + dx/2
            x1 = x0 + (w-1) * dx
            y1 = gt[4] - dy/2
            y0 = y1 - (h-1)*dy

            proj4 = strip(AG.toPROJ4(AG.importWKT(AG.getproj(dataset))))
            Gridded(x0:dx:x1, y0:dy:y1, mat[:,end:-1:1], Matrix{T}(0,0), true, proj4)
        end
    end
end

"""
    write_geotiff(g::Vector{Gridded{T}}, filepath::AbstractString;
                       size_=size(g[1]), # if the geotiff size is different to gridded size
                       xinds = 1:length(g[1].x), # to place it
                       yinds = 1:length(g[1].y), #
                       nodataval=[NaN, -3.697314e28][1] # https://maurow.bitbucket.io/notes/packbits-geotiffs.html
                       ) where T

Write a vector of Gridded into the bands of a geotiff.
"""
function write_geotiff(g::Vector{Gridded{T}}, filepath::AbstractString;
                       size_=size(g[1]), # if the geotiff size is different to gridded size
                       xinds = 1:length(g[1].x), # to place it
                       yinds = 1:length(g[1].y), #
                       nodataval=[NaN, -3.697314e28][1], # https://maurow.bitbucket.io/notes/packbits-geotiffs.html
                       colornames=["gray", ["undefined" for i=2:length(g)]...]) where T
    dx = step(g[1].x)
    dy = step(g[1].y)
    AG = ArchGDAL
    AG.registerdrivers() do
        AG.create(
            filepath,
            AG.getdriver("GTiff"),
            width = size_[1],
            height = size_[2],
            nbands = length(g),
            dtype = T
        ) do raster
            AG.setproj!(raster, AG.toWKT( AG.importPROJ4(g[1].proj)))
            # https://www.gdal.org/gdal_datamodel.html
            top = g[1].y[end] + dy/2
            left = g[1].x[1] - dx/2
            AG.setgeotransform!(raster, [left, dx, 0, top, 0, -dy])
            for (i,gg) in enumerate(g)
                @assert gg.x==g[1].x
                @assert gg.y==g[1].y
                AG.write!(
                    raster,
                    gg.v[:,end:-1:1],
                    i, # update band i
                    yinds, # along (window) xcoords
                    xinds # along (window) ycoords
                )
                AG.setnodatavalue!(AG.getband(raster,i), nodataval)
                # AG.setname!(AG.getband(raster,i), colornames[i])  # does not work
            end
        end
    end
    nothing
end

"""
Use gdal_merge.py

https://github.com/yeesian/ArchGDAL.jl/issues/39
"""
function update_geotiff(g::Gridded{T}, filepath::AbstractString) where T
    error("not implemented")

    AG = ArchGDAL
    # figure out extent
    w, h, na, gt, proj4 = AG.registerdrivers() do
        AG.read(filepath) do dataset
            band = AG.getband(dataset, 1)
            w, h = AG.width(band), AG.height(band)
            na = AG.getnodatavalue(band)
            gt = AG.getgeotransform(dataset)
            proj4 = strip(AG.toPROJ4(AG.importWKT(AG.getproj(dataset))))
            (w, h, na, gt, proj4)
        end
    end

    AG.registerdrivers() do
        AG.update(filepath) do raster
            AG.write!(
                raster,
                g.v[:,end:-1:1],
                1, # update band 1
                xcoords, # along (window) xcoords
                ycoords # along (window) ycoords
            )
        end
    end
    nothing
end

function compress_geotiff(fl, method=PACKBITS)
    AG = ArchGDAL
    fl1 = splitext(fl)[1]*"-comp.tif"
    nodata = AG.registerdrivers() do
        AG.read(fl) do dataset
            band = AG.getband(dataset, 1)
            AG.getnodatavalue(band)
        end
    end
    run(`gdalwarp -overwrite -srcnodata $nodata -dstnodata -3.697314e28  -ot float32 -co COMPRESS=$method $fl $fl1`)
end

## Shapefiles
import Shapefile, GeoInterface

"""
    read_shapefile_poly(filename)

Reads a shapefile and returns a list of the polygons it contains.
"""
function read_shapefile_poly(fln)
    shp = open(fln) do fd
        read(fd, Shapefile.Handle)
    end
    outt = []
    for ii =1:length(shp.shapes)
        poly = shp.shapes[1]
        np = length(poly.points)
        pp = poly.parts
        push!(pp, np)

        out = zeros(2,np+length(pp)-2)
        extra_j = 0
        for i = 2:length(pp)
            for j = (pp[i-1]+1:pp[i])+extra_j
                out[:,j] = GeoInterface.coordinates(poly.points[j-extra_j])
            end
            if i!=length(pp)
                out[:,pp[i]+extra_j+1] = out[:,1]
                extra_j += 1
            end
        end
        push!(outt, out)
    end
    return outt
end

## Misc helpers
###############
# """
# Cuts a smaller array out of a bigger one.
# """
#function cut_array(ar, )

# Does not work.  Probably better call out to gdalsrsinfo, see above.
# """
# Convert Well-Know-Text to Proj4.  Only works if there is a ESRI or
# ESPG number embedded in the WKT.
# """
# function wkt2proj4(wkt)
#     auth = search(wkt, "AUTHORITY")
#     if length(auth)==0
#         error("No AUTHORITY entry found in well-known-text")
#     end
#     endpos = search(wkt[auth[end]+1:end], "]")[1]
#     li = eval(parse(wkt[auth[end]+1:auth[end]+endpos]))
#     nr = parse(Int, li[2])
#     if li[1]=="EPSG"
#         return Proj4.epsg[nr]
#     elseif li[1]=="EPSG"
#         return Proj4.esri[nr]
#     else
#         error("AUTHORITY entry mal-formed: $li")
#     end
# end

"""
    transform_proj(xy or xyz, from, to)

Transform between projections.  Uses Proj4.jl.  Example:

    longlat = "+proj=longlat +datum=WGS84"
    utm56 = "+proj=utm +zone=56 +south +datum=WGS84"
    transform_proj([6e3,197e3], utm56, longlat)

Input:
- xy or xyz vector, matrix or trajectory of input points (x,y) or (x,y,z)
- `from` and `to` are either strings for Proj4.Projections (more preformant).

See VAWTools.Projections for pre-defined projections
"""
function transform_proj(xyz, from, to)
    from = Proj4.Projection(from)
    to = Proj4.Projection(to)
    transform_proj(xyz, from, to)
end
function transform_proj(xyz, from::Proj4.Projection, to::Proj4.Projection)
    Proj4.transform(from, to, xyz)
end
function transform_proj(tr::Traj, from, to)
    xy = transform_proj(hcat(tr.x, tr.y), from, to) # TODO: avoid temporary
    if haserror(tr)
        Traj(xy[:,1], xy[:,2], tr.v, tr.err, tr.splits)
    else
        Traj(xy[:,1], xy[:,2], tr.v, tr.splits)
    end
end

module Projections
using Proj4
"Swiss grid 1903"
const swiss1903 = Proj4.Projection("+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs")
#"Longitude latitude"
const longlat = Proj4.Projection("+proj=longlat +datum=WGS84")

# UTM zones: utm1n, utm1s, etc
# https://github.com/JuliaGeo/Proj4.jl
for z=1:60
    zs = Symbol("utm$(z)s")
    zn = Symbol("utm$(z)n")
    strs = "+proj=utm +zone=$z +south +datum=WGS84 +units=m +no_defs"
    strn = "+proj=utm +zone=$z +north +datum=WGS84 +units=m +no_defs"
    @eval $zs = Projection($strs)
    @eval $zn = Projection($strn)
end
end

######
# Distances
######
"""
    dist(x1,y1,x2,y2)

Cartesian distance
"""
dist(x1,y1,x2,y2) = sqrt((x1-x2)^2 + (y1-y2)^2)

"""
    dist_longlat(lon1,lat1,lon2,lat2,R=6373)

Distance on a sphere (default earth size).
"""
function dist_longlat(lon1,lat1,lon2,lat2,R=6373) # the radius of the Earth
    # http://andrew.hedges.name/experiments/haversine/

    dlon = deg2rad(lon2 - lon1)
    dlat = deg2rad(lat2 - lat1)
    lat1 = deg2rad(lat1)
    lat2 = deg2rad(lat2)
    a = (sin(dlat/2))^2 + cos(lat1) * cos(lat2) * (sin(dlon/2))^2
    c = 2 * atan2( sqrt(a), sqrt(1-a) )
    return R * c
end

########
"Convert int to string and pad with 0 to get to length 5"
int2str5(i) = @sprintf "%05d" i
"Convert int to string and pad with 0 to get to length 2bed"
int2str2(i) = @sprintf "%02d" i

"""
Represents an array all filled with one value.

Pretty incomplete implementation.  Probably needs something similar as
https://github.com/JuliaLang/julia/blob/86bf95fe0a76e4750d41f569bfcb6fa1fb1805e4/base/linalg/uniformscaling.jl.

TODO:
- arithmetic
"""
struct UniformArray{T,N} <: AbstractArray{T,N}
    val::T
end
Base.size(::UniformArray{T,N}) where {T,N} = ntuple(x->typemax(Int), Val(N))
Base.getindex(A::UniformArray, i::Int) = A.val
Base.IndexStyle(::Type{U}) where {U<:UniformArray} = IndexLinear()
Base.iterate(::UniformArray, args...) = error("Cannot iterate over UniformArray")
# Does not work https://github.com/JuliaLang/julia/issues/18004
#Base.show{T,N}(io::IO, u::UniformArray{T,N}) = print(io, "UniformArray{$T,$N} with value $(u.val)")
Base.show(io::IO, ::MIME"text/plain", u::UniformArray{T,N}) where {T,N} = print(io, "UniformArray{$T,$N} with value $(u.val)")

#############
# Polygons
#############
"""
    windnr(p, poly::Matrix)

Determines the winding number of a point and a polygon, i.e. how many
times a polygon winds around the point.

It follows Dan Sunday: http://geomalgorithms.com/a03-_inclusion.html.
"""
function windnr(p, poly::AbstractMatrix)
    @assert size(poly,1)==2
    @assert size(poly,2)>1
    @assert length(p)==2
    @assert poly[:,1]==poly[:,end]
    # Loop over edges
    px = p[1]
    py = p[2]
    wn = 0 # winding number, inside if odd
    len = length(poly)-2
    @inbounds for k=1:2:len #size(poly,2)-1 # @fastmath makes it worse
        # Coordinates of edge endpoints
        ex1,ey1,ex2,ey2 = poly[k], poly[k+1], poly[k+2], poly[k+3]

        # Check edges intersecting a horizontal ray:
        # rules http://geomalgorithms.com/a03-_inclusion.html#Edge-Crossing-Rules
        orient = leftorright(px,py,ex1,ey1,ex2,ey2)
        if up(ey1,ey2)
            # an upward edge includes its starting endpoint, and
            # excludes its final endpoint;
            !(py>=ey1 && py<ey2) && continue
            if orient==-1; wn-=1 end # p strictly left of oriented e
        elseif down(ey1,ey2)
            # a downward edge excludes its starting endpoint, and
            # includes its final endpoint;
            !(py<ey1 && py>=ey2) && continue
            if orient==1; wn+=1 end # p strictly right of oriented e
        else # Horizontal edges are excluded.
            continue
        end
        # NOTE: Rule "the edge-ray intersection point must be strictly
        # right of the point P."  implies that points on a right-side
        # boundary edge being outside, and ones on a left-side edge
        # being inside.
    end
    return wn
end
# edge goes up
up(ey1,ey2)   = ey1<ey2
# edge goes down
down(ey1,ey2) = ey1>ey2
function leftorright(px,py,ex1,ey1,ex2,ey2)
    # returns:
    # -1 if on left of line
    #  0 if on line
    #  1 if on right of line
    vx,vy = ex2-ex1, ey2-ey1
    px,py =  px-ex1,  py-ey1
    sign(px*vy-py*vx)
end

"""
Determines if a point is inside a polygon.

- p -- point (x,y) or [x,y]
- poly -- polygon vertices [x1 x2 ... xn x1
                            y1 y2 ... yn y1]
          (a closed poly)

Returns true if point has an odd winding number.  This should label
points as exterior which are inside outcrops.  See test for a test.

Note that points located on a right-side boundary edge are outside,
and ones on a left-side edge are inside.
"""
inpoly(p, poly::AbstractMatrix) = isodd(windnr(p,poly))

# Other inpoly algo here:
# https://github.com/helenchg/PolygonClipping.jl/blob/1fc74ab797c6585795283749b7c4cc9cb2000243/src/PolygonClipping.jl#L139

"""
    absslope(g::Gridded,weights::AbstractMatrix=UniformArray{Bool,2}(true),retnan=false)
    absslope(x::Range,y::Range,v,weights::AbstractMatrix=UniformArray{Bool,2}(true),retnan=true)

Absolute value of slope angle (by finite differences)

- uses a 3x3 points stencil.
- on outermost points uses one-sided stencil

In:
- gridded elevation set
Out:
- slope angle (rad)

Note:
- slope and angle are very similar up to about ~0.4 or 20deg
"""
absslope(g::Gridded, weights::AbstractMatrix=UniformArray{Bool,2}(true), retnan=true) =
    absslope(g.x,g.y,g.v,weights,retnan)
function absslope(x::AbstractRange, y::AbstractRange, v, weights::AbstractMatrix=UniformArray{Bool,2}(true), retnan=true)
    dvx,dvy = gradient3by3(x,y,v,weights,retnan)
    for i in eachindex(dvx)
        dvx[i] = atan(sqrt(dvx[i]^2+dvy[i]^2))
    end
    dvx
end

## Old, probably faster implementation:
# absslope(g::Gridded) = absslope(g.x,g.y,g.v)
# function absslope(x::Range,y::Range,v)
#     nx, ny = length(x),length(y)
#     dx = step(x)
#     alphas = zeros(Float64, nx, ny)
#     for j=2:ny-1, i=2:nx-1
#         dvx =  ((v[i+1,j+1]+2*v[i+1,j]+v[i+1,j-1])-(v[i-1,j+1]+2*v[i-1,j]+v[i-1,j-1]))/(8*dx)
#         dvy = -((v[i-1,j-1]+2*v[i,j-1]+v[i+1,j-1])-(v[i-1,j+1]+2*v[i,j+1]+v[i+1,j+1]))/(8*dx)
#         alphas[i,j] = atan(sqrt(dvx^2+dvy^2))
#     end
#     # on edge and corners use one-sided f-d
#     for j=[1,ny], i=1:nx
#         if i==1
#             dvx = (v[i+1,j]-v[i,j])/dx
#         elseif i==nx
#             dvx = (v[i,j]-v[i-1,j])/dx
#         else
#             dvx = (v[i+1,j]-v[i-1,j])/(2*dx)
#         end
#         if j==1
#             dvy = (v[i,j+1]-v[i,j])/dx
#         else
#             dvy = (v[i,j]-v[i,j-1])/dx
#         end
#         alphas[i,j] = atan(sqrt(dvx^2+dvy^2))
#     end
#     for j=1:ny, i=[1,nx]
#         if i==1
#             dvx = (v[i+1,j]-v[i,j])/dx
#         else
#             dvx = (v[i,j]-v[i-1,j])/dx
#         end
#         if j==1
#             dvy = (v[i,j+1]-v[i,j])/dx
#         elseif j==ny
#             dvy = (v[i,j]-v[i,j-1])/dx
#         else
#             dvy = (v[i,j+1]-v[i,j-1])/(2*dx)
#         end
#         alphas[i,j] = atan(sqrt(dvx^2+dvy^2))
#     end
#     return alphas
# end

"""
    gradient3by3(g::Gridded)
    gradient3by3(x::Range,y::Range,v)
    gradient3by3(x::Range,y::Range,v,weights::AbstractMatrix,retnan=true)

2D gradient (by finite differences)

- Averaged over 3x3 points.
- no slope is calculated for the outermost points
- with `weights` is given then it can be used for masked areas. If
  `retnan==false` a 0 is returned where no gradient can be made.

In:
- gridded elevation set or x,y,z
Out:
- dvx,dvy

Reference: gradient3by3_fast

Notes:
- maybe this should be done using one of the Images.jl filters instead
"""
gradient3by3(g::Gridded,weights::AbstractMatrix=UniformArray{Bool,2}(true),retnan=true) =
    gradient3by3(g.x,g.y,g.v,weights,retnan)
# Implementation using weights: drops pairs which don't work.
function gradient3by3(x::AbstractRange,y::AbstractRange,v,weights::AbstractMatrix=UniformArray{Bool,2}(true),retnan=true)
    nx, ny = length(x),length(y)
    dx = step(x)
    dvx = zeros(Float64, nx, ny)
    dvy = zeros(Float64, nx, ny)

    # local weights
    lw = zeros(Float64,3,3)
    f = [1,2,1]
    for j=1:ny, i=1:nx
        lw[:] .= 0
        for jj=max(1,j-1):min(ny,j+1), ii=max(1,i-1):min(nx,i+1)
            lw[ii-i+2,jj-j+2] = weights[ii,jj]
        end

        # xdir
        n = 0
        for k=-1:1
            if lw[3,k+2]!=0 && lw[1,k+2]!=0
                # centered difference
                dvx[i,j] += (v[i+1,j+k]-v[i-1,j+k])/(2*dx)*f[k+2]
                n+=1*f[k+2]
            elseif lw[2,k+2]!=0 && lw[1,k+2]!=0
                # right difference
                dvx[i,j] += (v[i,j+k]-v[i-1,j+k])/dx*f[k+2]
                n+=1*f[k+2]
            elseif lw[3,k+2]!=0 && lw[2,k+2]!=0
                # left difference
                dvx[i,j] += (v[i+1,j+k]-v[i,j+k])/dx*f[k+2]
                n+=1*f[k+2]
            end
        end
        if n==0;
            if retnan
                dvx[i,j] = NaN
            else
                dvx[i,j] = 0
            end
        else
            dvx[i,j] /= n
        end
        # ydir
        n = 0
        for k=-1:1
            if lw[k+2,3]!=0 && lw[k+2,1]!=0
                # centered difference
                dvy[i,j] += (v[i+k,j+1]-v[i+k,j-1])/(2*dx)*f[k+2]
                n+=1*f[k+2]
            elseif lw[k+2,3]!=0 && lw[k+2,2]!=0
                # right difference
                dvy[i,j] += (v[i+k,j+1]-v[i+k,j])/dx*f[k+2]
                n+=1*f[k+2]
            elseif lw[k+2,2]!=0 && lw[k+2,1]!=0
                # left difference
                dvy[i,j] += (v[i+k,j]-v[i+k,j-1])/dx*f[k+2]
                n+=1*f[k+2]
            end
        end
        if n==0;
            if retnan
                dvy[i,j] = NaN
            else
                dvy[i,j] = 0
            end
        else
            dvy[i,j] /= n
        end
    end
    return dvx, dvy
end

"Faster than gradient3by3 but does not deal with weights nor the edge."
gradient3by3_fast(g::Gridded) = gradient3by3_fast(g.x,g.y,g.v)
function gradient3by3_fast(x::AbstractRange,y::AbstractRange,v)
    nx, ny = length(x),length(y)
    dx = step(x)
    dvx = zeros(Float64, nx, ny)
    dvy = zeros(Float64, nx, ny)
    for j=2:ny-1, i=2:nx-1
        dvx[i,j] =  ((v[i+1,j+1]+2*v[i+1,j]+v[i+1,j-1]) - (v[i-1,j+1]+2*v[i-1,j]+v[i-1,j-1]))/(8*dx)
        dvy[i,j] = -((v[i-1,j-1]+2*v[i,j-1]+v[i+1,j-1]) - (v[i-1,j+1]+2*v[i,j+1]+v[i+1,j+1]))/(8*dx)
    end
    return dvx, dvy
end

##############
# Bands
##############
include("elevation-bands.jl")

#################
# Filtering
#################

"""
Return the index range of points within the distance dist.
"""
function around(g::Gridded, i, j, dist)
    inds = floor(Int,dist/step(g))
    nx, ny = size(g)
    return max(1,i-inds):min(nx,i+inds), max(1,j-inds):min(ny,j+inds)
end

"""
Return the linear index of points within the distance dist provided
they are within the mask.
"""
function around(g::Gridded, mask, i, j, dist)
    is, js = around(g::Gridded, i, j, dist)
    dims = size(mask)
    l = Int[] # linear index
    for j=js, i=is
        if mask[i,j]
            push!(l, ind2sub(mask, i, j))
        end
    end
    return l
end

"Mean ignoring NaNs"
function meannan(a)
    n = 0
    cum = zero(eltype(a))
    @inbounds for i in eachindex(a)
        if !isnan(a[i])
            n+=1
            cum+=a[i]
        end
    end
    cum/n # == NaN if n==0
end

"Extrema ignoring NaNs"
function extremanan(a)
    mi,ma = zero(eltype(a)),zero(eltype(a))
    @inbounds for i in eachindex(a)
        if !isnan(a[i])
            mi = min(mi, a[i])
            ma = max(ma, a[i])
        end
    end
    (mi,ma)
end

"Maximunm ignoring NaNs"
function maximumnan(a)
    mi,ma = zero(eltype(a)),zero(eltype(a))
    @inbounds for i in eachindex(a)
        if !isnan(a[i])
            ma = max(ma, a[i])
        end
    end
    ma
end

"Std ignoring NaNs"
function stdnan(a)
    n = 0
    inds = Int[]
    @inbounds for i in eachindex(a)
        if !isnan(a[i])
            push!(inds,i)
        end
    end
    std(a[inds])
end

"Mean ignoring fill"
function meanfill(a,fill,verbose=false)
    n = 0
    cum = zero(eltype(a))
    @inbounds for i in eachindex(a)
        if a[i]!=fill
            n+=1
            cum+=a[i]
        end
    end
    if n==0
        verbose && println("Could not fill a gap! Returning fill.")
        return fill
    end
    cum/n # == NaN if n==0
end

"weighted mean"
function mean_weighted(a,w)
    # make an accumulator type closed under addition (needed for Bools):
    n = zero(one(eltype(w)) + one(eltype(w)))
    cum = zero(eltype(a))
    @inbounds for i in eachindex(a)
        n+=w[i]
        cum+=a[i]*w[i]
    end
    if n==0
        error("zero weights")
    end
    cum/n # == NaN if n==0
end

"""
    boxcar(A::AbstractArray, window, [, weights, keepmask])
    boxcar(A::AbstractArray, windows::Tuple, [, weights, keepmask])
    boxcar(A::AbstractArray, window::AbstractArray, [, weights, keepmask])

Boxcar filter.  The two argument call skips NaNs.  The three & four
argument call uses weights and propagates NaNs, it can be a lot faster.

Smoothing occurs over +/-window indices, 0 corresponds to no smoothing.
The window can be specified as:
- integer, for a symmetric window in all dimensions
- a tuple to give lower and upper windows
- a tuple of tuples to give different lower and upper windows for all dimensions
- a array of size(A) for a different, symmetric window at each point.

Weights, if given, will use those relative weights for averaging.  Note that
points which have value==NaN and weight==0 will not poison the result.

No average is calculated for points where keepmask==true, instead
their original value will be kept.

Notes:
- For the weights it may be faster to use non-Bool arrays nor BitArrays,
  say Int8.  Note that NaNs where weight==0 will not poison the result.
- Also works for Vectors.

From http://julialang.org/blog/2016/02/iteration
"""
boxcar(A::AbstractArray, window) = boxcar(A, (window,window))
function boxcar(A::AbstractArray, windows::Tuple)
    window_lower, window_upper = windows
    out = similar(A)
    R = CartesianIndices(size(A))
    I1, Iend = first(R), last(R)
    I_l = CartesianIndex(I1.I.*window_lower)
    I_u = CartesianIndex(I1.I.*window_upper)
    for I in R # @inbounds does not help
        out[I] = NaN
        n, s = 0, zero(eltype(out))
        for J in CartesianIndices(UnitRange.(max(I1, I-I_l).I , min(Iend, I+I_u).I) ) # used to be
        # for J in max(I1, I-I_l):min(Iend, I+I_u) ## in Julia 1.1
            if !isnan(A[J])
                s += A[J]
                n += 1
            end
        end
        out[I] = s/n
    end
    out
end
function boxcar(A::AbstractArray, windows::Tuple{<:AbstractFloat,<:AbstractFloat})
    window_lower, window_upper = map(x->floor(Int,x), windows)
    weight_lower, weight_upper = windows[1]-window_lower, windows[2]-window_upper
    out = similar(A)
    # make an accumulator type closed under addition (needed for Bools):
    R = CartesianIndices(size(A))
    I1, Iend = first(R), last(R)
    I_l = CartesianIndex(I1.I.*window_lower)
    I_u = CartesianIndex(I1.I.*window_upper)
    I_ll = CartesianIndex(I1.I.*(window_lower+1))
    I_uu = CartesianIndex(I1.I.*(window_upper+1))
    for I in R # @inbounds does not help
        out[I] = NaN
        n, s = zero(eltype(out)), zero(eltype(out))
        # lower fractional-cells
        for J in CartesianIndices( UnitRange.(  max(I1, I-I_ll).I , (I-I_l-one(I)).I) ) # CartesianIndices(max(I1, I-I_ll), I-I_l-1)
            if !isnan(A[J])
                s += A[J] * weight_lower
                n += weight_lower
            end
        end
        # normal window
        for J in CartesianIndices(UnitRange.(max(I1, I-I_l).I , min(Iend, I+I_u).I) ) # CartesianIndices(max(I1, I-I_l), min(Iend, I+I_u))
            if !isnan(A[J])
                s += A[J]
                n += 1
            end
        end
        # upper fractional-cells
        for J in CartesianIndices(UnitRange.((I+I_u+one(I)).I , min(Iend, I+I_uu).I) ) # CartesianIndices(I+I_u+1, min(Iend, I+I_uu))
            if !isnan(A[J])
                s += A[J] * weight_upper
                n += weight_upper
            end
        end
        out[I] = s/n
    end
    out
end
function boxcar(A::AbstractArray, windows::Tuple{<:AbstractArray, <:AbstractArray})
    window_lower, window_upper = windows
    out = similar(A)
    R = CartesianIndices(size(A))
    I1, Iend = first(R), last(R)
    for I in R # @inbounds does not help
        out[I] = NaN
        n, s = 0, zero(eltype(out))
        I_l = CartesianIndex(I1.I.*window_lower[I])
        I_u = CartesianIndex(I1.I.*window_upper[I])
        for J in CartesianIndices(UnitRange.(max(I1, I-I_l).I , min(Iend, I+I_u).I) ) # used to be CartesianRange(max(I1, I-I_l), min(Iend, I+I_u) )
            if !isnan(A[J])
                s += A[J]
                n += 1
            end
        end
        out[I] = s/n
    end
    out
end

boxcar(A::AbstractArray, window, weights::AbstractArray, keepmask::AbstractArray=(weights.==0)) =
    boxcar(A, (window,window), weights, keepmask)
function boxcar(A::AbstractArray{T,N}, windows::Tuple,
                     weights::AbstractArray,
                     keepmask::AbstractArray=(weights.==0)) where {T,N}
    @assert size(weights)==size(A)
    window_lower, window_upper = windows
    out = fill(zero(T), size(A))
    # make an accumulator type closed under addition (needed for Bools):
    AT = typeof(one(eltype(weights)) + one(eltype(weights)))
    R = CartesianIndices(size(A))
    I1, Iend = first(R), last(R)
    I_l = CartesianIndex(I1.I.*window_lower)
    I_u = CartesianIndex(I1.I.*window_upper)
    @inbounds for I in R
        if keepmask[I]
            out[I] = A[I]
        else
            n, s = zero(AT), zero(T)
            for J in CartesianIndices(UnitRange.(max(I1, I-I_l).I , min(Iend, I+I_u).I) ) # used to be CartesianRange(max(I1, I-I_l), min(Iend, I+I_u) )
                AJ, w = A[J], weights[J]
                if w!=0
                    s += AJ * convert(T, w)
                    n += w
                end
            end
            if n==0
                #error("At location $I no contributing cells found")
                out[I] = NaN
            else
                out[I] = s/n
            end
        end
    end
    out
end

"""
    boxcar_matrix(T::DataType, window::Integer, weights::AbstractMatrix[, keepmask])

This produces a sparse matrix which can be used to apply the filter:
`bx*hs2d`. Relatively expensive to create but very fast to apply, thus
use when needing the same filter several times.

Notes:
- NaNs where weight==0 do not poison the result.

Apply with

    apply_boxcar_matrix(M, orig)
"""
function boxcar_matrix(::Type{T}, window::Integer,
                               weights::AbstractArray{TT,N},
                               keepmask::AbstractArray=(weights.==0)) where {T,TT,N}
                               # keepmask::AbstractArray=BitArray{N}(ntuple(x->0, Val{N})...))
    # make an accumulator type closed under addition (needed for Ints and Bools):
    Tacc = promote_type(T, eltype(weights))
    nr = size(weights,1)
    nc = size(weights,2)
    is = Int[]
    js = Int[]
    vs = T[]
    sizehint!(is, 2*window*nr*nc)
    sizehint!(js, 2*window*nr*nc)
    sizehint!(vs, 2*window*nr*nc)
    R = CartesianIndices(size(weights))
    I1, Iend = first(R), last(R)
    @inbounds @fastmath for I in R
        i = (I.I[2]-1)*nr + I.I[1] # row of output matrix
        if keepmask[I]
            # do not average at this location, preserve original value
            push!(is, i)
            push!(js, i)
            push!(vs, 1)
            continue
        end
        nrows = 0 # number of contributing cells
        acc = zero(Tacc) # sum of all weights for one cell
        for J in CartesianIndices(UnitRange.(max(I1, I-I1*window).I , min(Iend, I+I1*window).I) ) # CartesianIndices(max(I1, I-I1*window), min(Iend, I+I1*window))
            if weights[J]==0
                continue
            end
            j = (J.I[2]-1)*nr + J.I[1]
            push!(is, i)
            push!(js, j)
            push!(vs, weights[J])
            nrows +=1
            acc += weights[J]
        end
        # divide the current batch by the summed weights:
        lvs = length(vs)
        for n=lvs .- (0:nrows-1)
            vs[n] /= acc
        end
    end
    return sparse(is, js, vs, length(weights), length(weights))
end

"Apply the boxcar filter matrix and reshape result"
apply_boxcar_matrix(M, orig) = reshape(M*reshape(orig, length(orig)), size(orig))

"""
Returns a piecewise linear function through points (xs,ys).
For x<xs[1]->ys[1] and for x>xs[end] -> ys[end]

When length(xs)==1, then it returns a constant function.

Currently extrapolates using the last value. (TODO)

Almost 2x faster than using Interpolations.jl.
"""
function piecewiselinear(xs::AbstractVector, ys::AbstractVector)
    if length(xs)==1
        @assert length(xs)==length(ys)
        return xq -> ys[1]
    end
    rats = diff(ys)./diff(xs)
    if issorted(xs)
        return function (xq)
            # xq<xs[1] && error("cannot extrapolate")
            # xq>xs[end] && error("cannot extrapolate")
            # xq==xs[end] && return ys[end]
            xq<=xs[1] && return ys[1]
            xq>=xs[end] && return ys[end]
            ii = searchsortedlast(xs, xq)
            out = ys[ii] - (xs[ii]-xq)*rats[ii]
            return out
        end
    elseif issorted(xs, Base.Order.Reverse)
        return function (xq)
            # xq>xs[1] && error("cannot extrapolate")
            # xq<xs[end] && error("cannot extrapolate")
            # xq==xs[end] && return ys[end]
            xq>=xs[1] && return ys[1]
            xq<=xs[end] && return ys[end]
            ii = searchsortedlast(xs, xq, Base.Order.Reverse)
            out = ys[ii] - (xs[ii]-xq)*rats[ii]
            return out
        end
    else
        error("xs must be sorted")
    end
end
# """
# Includes extrapolation:
# Returns a piecewise linear function through points (xs,ys).
# For x<xs[1]->ys[1] and for x>xs[end] -> ys[end]

# - assumes xs is sorted

# Almost 2x faster than using Interpolations.jl.
# """
# function piecewiselinear(xs_, ys_)
#     rats = diff(ys_)./diff(xs_)
#     ys = copy(ys_)
#     prepend!(ys,ys[1:1])
#     append!(ys,ys[end:end])

#     xs = copy(xs_)
#     prepend!(xs, xs[1:1])
#     append!(xs, xs[end:end])

#     prepend!(rats, [0.0])
#     append!(rats,  [0.0])

#     function (xq)
#         ii = searchsortedfirst(xs_, xq)
#         return ys[ii] - (xs[ii]-xq)*rats[ii]
#     end
# end

###########
# Smoothing splines
###########
import StatsBase, SmoothingSplines
const SSp = SmoothingSplines

"""

Smoothes a vector of datapoints with a spline with minimum wave-length
components `len`.  Returns a smoothed vector of the same length.

Notes:
- if len=0 no smoothing occurs
- if length(x)==1 or ==2 no smoothing occurs
- better than boxcar if having no edge-effects is important
"""
function smooth_vector(x, y::AbstractVector{T}, len, out=y) where T
    if len==0 && out==y
        return y
    elseif length(y)==1 # no smoothing, horizontal line
        return fill!(similar(y, length(out)), y)
    elseif length(y)==2 # no smoothing, line
        a = (y[2]-y[1])/(x[2]-x[1])
        b = -a*x[2] + y[2]
        return convert(Vector{T}, a.*out + b)
    else
        lambda0 = 0.001 # by trial and error
        lambda = convert(T, lambda0 * len^3) # the ^3 seems correct (trial and error)
        if out==y
            return SSp.predict(StatsBase.fit(SSp.SmoothingSpline, convert(Vector{T}, x), y, lambda))
        else
            return SSp.predict(StatsBase.fit(SSp.SmoothingSpline, convert(Vector{T}, x), y, lambda), out)
        end
    end
end

# using Dierckx
## using the number of nodes does funny stuff to the curve...
# """
# Smoothes a vector using a spline with nknots.  Returns a smoothed
# vector of the same length.

# Notes:
# - if nknots+2>=length(x) then it just returns y
# - better than boxcar if having no edge-effects is important
# """
# function smooth(x, y, nknots)
#     if nknots+2>=length(x)
#         return y
#     else
#         if !issorted(x)
#             perm = sortperm(x)
#             x_ = x[perm]
#             y = y[perm]
#         else
#             x_ = x
#         end
#         xknots = linspace(x_[1],x_[end],nknots+2)[2:end-1]
#         return Spline1D(x_, y, xknots)(x)
#     end
# end

end # module


text:
module Data

using Artifacts
using VsRecorderBase

import SimpleI18n
using SimpleI18n: i18n
import ..VsRecorder
using ..VsRecorder: Stats, Pokemon, Team,
    VsAbstractID,
    PokemonID, ItemID, MoveID, AbilityID,
    STATS_FIELDS, DEX_STATS_FIELDS,
    Nature, PokemonType, EnumIDTypes,
    Gender, GENDER_NULL, GENDER_MALE, GENDER_FEMALE,
    PokemonBattle, PokemonContext,
    datapath
using ..VsI18n

include("utils.jl")

export poke_dex, move_dex,
    item_list, ability_list, poke_list, move_list,
    search_dex
include("dex.jl")

export pokemon_icons, item_icons, gender_icons, type_icons
include("icon_sheets.jl")

export import_team, export_team
include("poke_pastes.jl")

function __init__()
    initialize_dex()
    nothing
end

end

ext: jl


text:
import Base.LinAlg.checksquare
import Base.LinAlg.BlasFloat
import MathProgBase
import Clp

"""
This function computes the incidence matrix for the
matrix A assuming that it is symmetric, and only
uses the upper-triangular portion.
"""
function incidence_matrix_and_weight{T}(A::SparseMatrixCSC{T,Int64})
    n = checksquare(A)
    At = triu(A,1) # find the upper-triangular portion

    ei,ej,ev = findnz(At)
    m = length(ei)
    #@show length([1:m; 1:m]), length([ei; ej]), length([ones(T,m); -ones(T,m)])

    B = sparse([1:m; 1:m], [ei; ej], [ones(T,m); -ones(T,m)], m, n)
    return (B,ev)
end

export incidence_matrix_and_weight

function stmincut{T<:BlasFloat}(A::SparseMatrixCSC{T,Int64}, s::Int64, t::Int64)
    # create the mincut LP
    # minimize_x sum_(edge i,j) |w_{i,j} (x_i - x_j)| s.t. x_s = 1, x_t = 0
    # <->
    # minimize_x,y+,y- w'(y+ + y-)
    #   y+ - y- = B*x
    #   y+ >= 0
    #   y- >= 0
    #   x_s = 1, x_t = 0
    # <->
    # minimize_z [0 w' w']*z
    # s.t.
    #   [B -I +I]*[x; y+; y-] = B*x - y+ + y- = 0
    #   xs = 1
    #
    #   y+ >= 0

    B,w = incidence_matrix_and_weight(A)

    n = size(A,1)
    m = size(B,1)

    @assert 1 <= s <= n
    @assert 1 <= t <= n

    eint = Array(Int64,0)
    et = Array(T,0)

    zn = sparse(eint,eint,et,1,n)
    zm = sparse(eint,eint,et,1,m)
    es = sparse([1], [s], one(T), 1, n)
    et = sparse([1], [t], one(T), 1, n)

    f = [zeros(n); w; w]
    Aeq = [B -speye(m) speye(m);
           es zm zm;
           et zm zm]
    beq = [zeros(m); 1.; 0.]

    sol = MathProgBase.linprog(f, Aeq, '=', beq, 0., 1., Clp.ClpSolver())
    x = sol.sol[1:n]
    cutval = sol.objval

    return x, cutval, sol
end

export stmincut

function cut(A::SparseMatrixCSC{Float64,Int},R::Set{Int64})
    n = size(A,1)
    colptr = A.colptr
    rowval = A.rowval
    nzval = A.nzval
    cutval = 0.
    for z in R
        # z in a vertex in A
        for nzi in colptr[z]:(colptr[z+1] - 1)
            v = rowval[nzi]
            if !(v in R)
                cutval += nzval[nzi]
            end
        end
    end
    return cutval
end

function _create_mqi_reference_cut_graph(B,dr,dout,wcut,wvol)
    C = [spzeros(1,1) wcut*sparse(dr') spzeros(1,1);
        wcut*sparse(dr) wvol*B wvol*sparse(dout);
        spzeros(1,1) wvol*sparse(dout') spzeros(1,1)]
end

function mqi(A,R::Set{Int64};maxiter=10)

    # turn R into a vector
    rvec = collect(R)

    # get the degree vector
    d = sum(A,1) # get a summation
    dr = d[rvec]
    rvol = sum(dr)
    @assert rvol <= sum(d)/2

    # get the initial cut
    cutval = cut(A,R)

    B = A[rvec,rvec]
    dout = dr-vec(sum(B,1))

    nb = size(B,1)

    W = deepcopy(R) # make a copy
    wvol = rvol
    wcut = cutval

    @printf("mqi: cut=%7.2g  vol=%7.2g  r=%.5f graph-edges=%7i\n", wcut, wvol, wcut/wvol, 0)

    # create the st graph
    C = _create_mqi_reference_cut_graph(B,dr,dout,wcut,wvol)
    x,val = stmincut(C,1,size(C,1))
    #display(full(C))

    iter = 1

    # add the MQI loop
    while val < wvol*wcut && iter <= maxiter

        # then there is a smaller conductance set inside
        W = Set{Int64}(rvec[find(x[2:end-1] .> 0.5)]) # TODO handle non-integer solutions
        wcut = cut(A,W)
        wvol = sum(dr[ x[2:end-1] .> 0.5 ])

        @printf("mqi: cut=%7.2g  vol=%8.2g  r=%.5f iter=%3i\n", wcut, wvol, wcut/wvol, iter)

        C = _create_mqi_reference_cut_graph(B,dr,dout,wcut,wvol)
        x,val = stmincut(C,1,size(C,1))

        iter += 1
    end

    @printf("mqi: cut=%7.2g  vol=%8.2g  r=%.5f iter=%3i\n", wcut, wvol, wcut/wvol, iter)

    if iter == maxiter && val < wvol*wcut
        warn("iteration hit limit")
    end

    #S = Set{Int64}(find(x[2:end-1] .> 0.5))
    S = W

    return S, wcut, wvol, x

end

export mqi

function _create_flow_improve_reference_cut_graph(A,dr,drbar,alpha,ep)
    C = [spzeros(1,1) alpha*dr' spzeros(1,1);
        alpha*dr A alpha*ep*sparse(drbar);
        spzeros(1,1) alpha*ep*sparse(drbar') spzeros(1,1)]
end

function flow_improve(A,R::Set{Int64};maxiter=10,kappa=0.)

    # kappa is the parameter that controls the localization
    # kappa = 0. is the original flow improve algorithm
    # kappa > 0. gives the Orecchia & Zhu locally modified
    # version

    n = size(A,1)

    # turn R into a vector
    rvec = collect(R)
    rbarvec = setdiff(1:n,rvec)

    # get the degree vector
    d = sum(A,1) # get a summation
    rvol = sum(d[rvec])
    gvol = sum(d)
    @assert rvol <= gvol/2

    ep = (rvol/(gvol-rvol))*exp(kappa)

    # get the initial cut
    cutval = cut(A,R)

    dr = sparse(rvec,ones(Int64,length(rvec)),d[rvec],n,1)
    drbar = sparse(rbarvec,ones(Int64,length(rbarvec)),d[rbarvec],n,1)

    W = deepcopy(R) # make a copy

    wvol = rvol
    wadjvol = rvol # initially the adjusted volume is the same
    wcut = cutval
    alpha = wcut/wadjvol
    alphaold = alpha

    @printf("flow_improve: cut=%15g  vol=%15g  r=%.5f  alpha=%.5f graph-edges=%7i\n",
                wcut, wvol, wcut/wvol, alpha, gvol)

    # create the st graph
    C = _create_flow_improve_reference_cut_graph(A,dr,drbar,alpha,ep)
    x,val = stmincut(C,1,size(C,1))

    #@printf("Target val: %15g\n", alpha*rvol - alpha*wadjvol + wcut)
    #@printf(" Flow  val: %15g\n", val)

    iter = 1

    #@show val, alpha*wvol + wcut - alpha*wadjvol, wcut, wvol, wadjvol, alpha, ep

    # add the loop
    while val < alpha*rvol - alpha*wadjvol + wcut && iter <= maxiter

        # the cost of a cut in graph C is
        # alpha*rvol - alpha*adjusted-vol + cut
        # so our current score is
        # alpha*rvol - alpha*wadjvol + cut
        # and if the previous computation gave a better score,
        # then we found a better set!

        alphaold = alpha

        # get the cut and compute it's value
        W = Set{Int64}(find(x[2:end-1] .> 0.5)) # TODO handle non-integer solutions
        wcut = cut(A,W)
        wvol = sum(d[ x[2:end-1] .> 0.5 ])
        wadjvol = sum( dr[collect(W)]) - ep*sum(drbar[collect(W)])
        alpha = wcut/wadjvol

        if alphaold <= alpha
            # due to some floating point issues, we might not have detected
            # the stopping condition above, this check will get it
            # essentially, we made no progress, which means the algorithm
            # will stall
            break
        end

        @printf("flow_improve: cut=%15g  vol=%15g  r=%.5f  alpha=%.5f iter=%3i\n", wcut, wvol, wcut/wvol, alpha, iter)

        C = _create_flow_improve_reference_cut_graph(A,dr,drbar,alpha,ep)

        x,val = stmincut(C,1,size(C,1))

        #@printf("finished iteration %i: val=%.16f, cut=%16g, vol=%16g, alpha=%16g\n", iter, val, wcut, wvol, alpha )

        iter += 1
    end

    @printf("flow_improve: cut=%15g  vol=%15g  r=%.5f  alpha=%.5f iter=%3i (done)\n", wcut, wvol, wcut/wvol, alpha, iter)

    if iter == maxiter && alpha < alphaold
        warn("iteration hit limit")
    end

    S = W

    return S, wcut, wvol, x
end

export flow_improve

function connected_cliques(sizes::Vector{Int64},ring::Bool)
    # TODO add xy coords
    # create the index sets for each clique
    n = 0
    nedges = 2*(length(sizes)-1) # the number of connecting edges
    for s in sizes
        nedges += s*(s-1)
        n += s
    end
    if ring
        nedges += 2
    end

    ei = zeros(Int64,nedges)
    ej = zeros(Int64,nedges)

    z = 1
    offset = 0
    for s in sizes
        # create a clique of size s with edges starting at z
        for i=1:s
            for j=1:s
                if i==j
                    continue
                end
                ei[z] = i+offset
                ej[z] = j+offset
                z += 1
            end
        end

        # add the connection
        offset += s
        if offset < n
            # we aren't at the end yet
            ei[z] = offset
            ej[z] = offset + 1
            z += 1
            ei[z] = offset + 1
            ej[z] = offset
            z += 1
        end
    end

    if ring
        ei[z] = n
        ej[z] = 1
        z += 1
        ei[z] = 1
        ej[z] = n
        z += 1
    end

    @assert z == nedges + 1
    A = sparse(ei,ej,1.,n,n)
    return A
end

export connected_cliques


text:
module SpacetimeFields

using NetCDF, DataFrames, Statistics
ncm = NetCDF
import Base.convert, Base.copy, Base.show

export extent, stfield
export copy, llgrid, nc2field, convert, show
export field2ts, array2field

"""
    extent(xmin, xmax, ymin, ymax)

An `extent` object
"""
struct extent
    xmin::Float64
    xmax::Float64
    ymin::Float64
    ymax::Float64
end

"""
    stfield(data, lon, lat, good, time)

A `stfield` object
"""
mutable struct stfield
    data::Array
    lon::Vector
    lat::Vector
    good::BitArray
    time::Vector
end

"""
    copy(x::stfield; layers::AbstractVector=[1])

Creates a copy of an `stfield` object
"""
function copy(x::stfield; layers::AbstractVector=[1])
    m1 = x.data[:,:,layers]
    stfield(m1, x.lon, x.lat, x.good, collect(layers))
end

"""
    llgrid(fname::String; start=[1], count=[-1])

Read a NetCDF file with the given `fname`, and return  a grid of longitude and latitude values
"""
function llgrid(fname::String; start=[1], count=[-1], llnames::Vector{String}=["lon","lat"])
    la = ncread(fname, llnames[1], start, count)
    lb = ncread(fname, llnames[2], start, count)
    ll = vcat([ [i j ] for i in la, j in lb ]...)
    lli = vcat([ [i j ] for i in 1:length(la), j in 1:length(lb) ]...)
    return ll, lli
end

"""
    llgrid(x::stfield)

Return a grid of longitude and latitude values from a `stfield` object
"""
function llgrid(x::stfield)
    ll = vcat([ [i j ] for i in x.lon, j in x.lat ]...)
    lli = vcat([ [i j ] for i in 1:length(x.lon), j in 1:length(x.lat) ]...)
    return ll, lli
end

"""
    nc2field(fname::String, varname::String, ext::extent; mean_dims=[], llnames=["lat","lon"])

Convert NetCDF data to a `stfield` object, after averaging over dimensions `mean_dims` (optional)
"""
function nc2field(fname::String, varname::String, ext::extent; mean_dims::AbstractVector=[], llnames::Vector{String}=String[])

    ncf = ncm.open(fname)
    v = ncf.vars[varname]
    dimnames = getfield.(v.dim, :name)
    isempty(llnames) && (llnames = dimnames[1:2])

    ll,lli = llgrid(fname, llnames=llnames)
    imin = argmin( hypot.(ll[:,1].-ext.xmin, ll[:,2].-ext.ymin ) )
    imax = argmin( hypot.(ll[:,1].-ext.xmax, ll[:,2].-ext.ymax ) )

    a, b = lli[imin,:], lli[imax,:]
    d = abs.(b - a) .+ [1,1]
    sign(b[2]-a[2])==-1 && (a[2] = b[2])

    start = ncm.defaultstart(ncf[varname])
    count = -ncm.defaultstart(ncf[varname])
    start[1:2] = a; count[1:2] = d
    a1 = ncm.readvar(ncf[varname], start=start, count=count)
    !isempty(mean_dims) && (a1 = mean(a1, dims=mean_dims))
    a1 = dropdims(a1, dims=(findall(size(a1).== 1)...,))

    miss_value = ncgetatt(fname, varname,"missing_value")
    good = vec(a1[:,:,1] .≠ miss_value)
    lon = ncm.readvar(ncf[llnames[1]], start=a[1:1], count=d[1:1])
    lat = ncm.readvar(ncf[llnames[2]], start=a[2:2], count=d[2:2])
    time = collect(1:size(a1, 3))
    ncm.close(ncf)
    return stfield(a1, lon, lat, good, time)
end

"""
    array2field(field::AbstractArray, lon, lat, ext::extent)

Convert a 3d-array to a `stfield` object
"""
function array2field(field::AbstractArray, lon::T, lat::T,  ext::extent) where {T<:AbstractVector}

    ll = vcat([ [i j] for i in lon, j in lat ]...)
    lli = vcat([ [i j ] for i in 1:length(lon), j in 1:length(lat) ]...)

    imin = argmin( hypot.(ll[:,1].-ext.xmin, ll[:,2].-ext.ymin ) )
    imax = argmin( hypot.(ll[:,1].-ext.xmax, ll[:,2].-ext.ymax ) )

    a = lli[imin,:];    b = lli[imax,:]
#    dump(a); dump(b)

    data = field[a[1]:b[1], a[2]:b[2], :]
    lon = lon[a[1]:b[1]]
    lat = lat[a[2]:b[2]]
    good = trues(length(data[:,:,1]))
    time = collect(1:size(field, 3))

    return stfield(data, lon, lat, good, time)
end

function field2ts(st::stfield, lon::T, lat::T) where T<:AbstractVector
    ll, lli = llgrid(st)
    gi = [argmin(hypot.(ll[:,1].-x, ll[:,2].-y)) for (x,y) in zip(lon, lat)]
    x, y = lli[gi,:]
    return st.data[x,y,:]
end

"""
    convert(::Type{Matrix}, x::stfield)

Convert `stfield` object to a `Matrix`.
"""
function convert(::Type{Matrix}, x::stfield)
    a, b, n = size(x.data)
    d0 = reshape(x.data, a*b, n)'
    d1 = d0[:,x.good]
    return(d1)
end

"""
    convert(::Type{DataFrame}, x::stfield, time::AbstractVector)

Convert `stfield` object to a `DataFrame`
"""
function convert(::Type{DataFrame}, x::stfield, time::AbstractVector)
    ll = vcat([[i j ] for i in x.lon, j in x.lat ]...)[x.good,:]
    D = vcat( [DataFrame(lon=ll[:,1],lat=ll[:,2], z = x.data[:,:,k][x.good], g="$k") for k in time]...)
    return D
end

"""
    convert(::Type{DataFrame}, x::stfield, time, label::Vector{String})

Convert `stfield` object to a `DataFrame`
"""
function convert(::Type{DataFrame}, x::stfield, time, label::Vector{String})
    ll = vcat([ [i j] for i in x.lon, j in x.lat ]...)[x.good,:]
    D = vcat( [DataFrame(lon=ll[:,1], lat=ll[:,2], z=x.data[:,:,k][x.good], g=label[k]) for k in time]...)
    return D
end

"""
    show(io::IO,x::stfield)

Display info about a `stfield` object
"""
function show(io::IO, x::stfield)
    println(typeof(x))
    print("Data: "); print(typeof(x.data)); println(size(x.data))
    print("Longitude: "); println(x.lon')
    print("Latitude: "); println(x.lat')
    print("Time: "); println(x.time')
end

end # module


text:
# okay, test to see if rejection sampling works

# a very simple rejection sampling algorithm

using Distributions

abstract type Sampler end
abstract type MCMCSampler <: Sampler end
abstract type RejectionSampler <: MCMCSampler end

type StandardRejectionSampler <: RejectionSampler
	envelope:: Function
	tries:: AbstractArray{Integer}
	samples:: AbstractArray{Number}

function rejectionsample(sampler<:RejectionSampler, distribution, num_samples, envelope_fn)

	samples =[]
	num_acceptances =0
	num_rejections = 0

	for i in 1:num_samples
		#somehow figure this out
		u = random_normal(0,1) # sample from a random random normal
		xstar = sample(envelope_fn) # get our thing
		sample_height = u*envelope_fn(xstar)
		distribution_height = distribution(xstar)
		if sample_height <= distribution_height
			push!(samples, sample_height)
			num_acceptances+=1
		end
	else
		num_rejections +=1
	end
	end
end

max_stars_repo_name: Bmillidgework/Misc-Maths

max_issues_repo_name: Bmillidgework/Misc-Maths

max_forks_repo_name: Bmillidgework/Misc-Maths


text:
# Julia wrapper for header: /Users/jgoldfar/Documents/work/projects/Hypre.jl/deps/usr/include/HYPRE.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /Users/jgoldfar/Documents/work/projects/Hypre.jl/deps/usr/include/HYPRE_FEI_includes.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

# Julia wrapper for header: /Users/jgoldfar/Documents/work/projects/Hypre.jl/deps/usr/include/HYPRE_IJ_mv.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function HYPRE_GetError()
    ccall((:HYPRE_GetError,HYPRE_utilities),HYPRE_Int,())
end

function HYPRE_CheckError(hypre_ierr::HYPRE_Int,hypre_error_code::HYPRE_Int)
    ccall((:HYPRE_CheckError,HYPRE_utilities),HYPRE_Int,(HYPRE_Int,HYPRE_Int),hypre_ierr,hypre_error_code)
end

function HYPRE_GetErrorArg()
    ccall((:HYPRE_GetErrorArg,HYPRE_utilities),HYPRE_Int,())
end

function HYPRE_DescribeError(hypre_ierr::HYPRE_Int,descr::Ptr{UInt8})
    ccall((:HYPRE_DescribeError,HYPRE_utilities),Void,(HYPRE_Int,Ptr{UInt8}),hypre_ierr,descr)
end

function HYPRE_ClearAllErrors()
    ccall((:HYPRE_ClearAllErrors,HYPRE_utilities),HYPRE_Int,())
end

function HYPRE_ClearError(hypre_error_code::HYPRE_Int)
    ccall((:HYPRE_ClearError,HYPRE_utilities),HYPRE_Int,(HYPRE_Int,),hypre_error_code)
end

function HYPRE_AssumedPartitionCheck()
    ccall((:HYPRE_AssumedPartitionCheck,HYPRE_utilities),HYPRE_Int,())
end

function HYPRE_IJMatrixCreate(comm::Cint,ilower::HYPRE_Int,iupper::HYPRE_Int,jlower::HYPRE_Int,jupper::HYPRE_Int,matrix::Ptr{HYPRE_IJMatrix})
    ccall((:HYPRE_IJMatrixCreate,HYPRE_IJ_mv),HYPRE_Int,(Cint,HYPRE_Int,HYPRE_Int,HYPRE_Int,HYPRE_Int,Ptr{HYPRE_IJMatrix}),comm,ilower,iupper,jlower,jupper,matrix)
end

function HYPRE_IJMatrixDestroy(matrix::HYPRE_IJMatrix)
    ccall((:HYPRE_IJMatrixDestroy,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,),matrix)
end

function HYPRE_IJMatrixInitialize(matrix::HYPRE_IJMatrix)
    ccall((:HYPRE_IJMatrixInitialize,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,),matrix)
end

function HYPRE_IJMatrixSetValues(matrix::HYPRE_IJMatrix,nrows::HYPRE_Int,ncols::Ptr{HYPRE_Int},rows::Ptr{HYPRE_Int},cols::Ptr{HYPRE_Int},values::Ptr{Cdouble})
    ccall((:HYPRE_IJMatrixSetValues,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,HYPRE_Int,Ptr{HYPRE_Int},Ptr{HYPRE_Int},Ptr{HYPRE_Int},Ptr{Cdouble}),matrix,nrows,ncols,rows,cols,values)
end

function HYPRE_IJMatrixAddToValues(matrix::HYPRE_IJMatrix,nrows::HYPRE_Int,ncols::Ptr{HYPRE_Int},rows::Ptr{HYPRE_Int},cols::Ptr{HYPRE_Int},values::Ptr{Cdouble})
    ccall((:HYPRE_IJMatrixAddToValues,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,HYPRE_Int,Ptr{HYPRE_Int},Ptr{HYPRE_Int},Ptr{HYPRE_Int},Ptr{Cdouble}),matrix,nrows,ncols,rows,cols,values)
end

function HYPRE_IJMatrixAssemble(matrix::HYPRE_IJMatrix)
    ccall((:HYPRE_IJMatrixAssemble,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,),matrix)
end

function HYPRE_IJMatrixGetRowCounts(matrix::HYPRE_IJMatrix,nrows::HYPRE_Int,rows::Ptr{HYPRE_Int},ncols::Ptr{HYPRE_Int})
    ccall((:HYPRE_IJMatrixGetRowCounts,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,HYPRE_Int,Ptr{HYPRE_Int},Ptr{HYPRE_Int}),matrix,nrows,rows,ncols)
end

function HYPRE_IJMatrixGetValues(matrix::HYPRE_IJMatrix,nrows::HYPRE_Int,ncols::Ptr{HYPRE_Int},rows::Ptr{HYPRE_Int},cols::Ptr{HYPRE_Int},values::Ptr{Cdouble})
    ccall((:HYPRE_IJMatrixGetValues,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,HYPRE_Int,Ptr{HYPRE_Int},Ptr{HYPRE_Int},Ptr{HYPRE_Int},Ptr{Cdouble}),matrix,nrows,ncols,rows,cols,values)
end

function HYPRE_IJMatrixSetObjectType(matrix::HYPRE_IJMatrix,_type::HYPRE_Int)
    ccall((:HYPRE_IJMatrixSetObjectType,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,HYPRE_Int),matrix,_type)
end

function HYPRE_IJMatrixGetObjectType(matrix::HYPRE_IJMatrix,_type::Ptr{HYPRE_Int})
    ccall((:HYPRE_IJMatrixGetObjectType,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,Ptr{HYPRE_Int}),matrix,_type)
end

function HYPRE_IJMatrixGetLocalRange(matrix::HYPRE_IJMatrix,ilower::Ptr{HYPRE_Int},iupper::Ptr{HYPRE_Int},jlower::Ptr{HYPRE_Int},jupper::Ptr{HYPRE_Int})
    ccall((:HYPRE_IJMatrixGetLocalRange,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,Ptr{HYPRE_Int},Ptr{HYPRE_Int},Ptr{HYPRE_Int},Ptr{HYPRE_Int}),matrix,ilower,iupper,jlower,jupper)
end

function HYPRE_IJMatrixGetObject(matrix::HYPRE_IJMatrix,object::Ptr{Ptr{Void}})
    ccall((:HYPRE_IJMatrixGetObject,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,Ptr{Ptr{Void}}),matrix,object)
end

function HYPRE_IJMatrixSetRowSizes(matrix::HYPRE_IJMatrix,sizes::Ptr{HYPRE_Int})
    ccall((:HYPRE_IJMatrixSetRowSizes,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,Ptr{HYPRE_Int}),matrix,sizes)
end

function HYPRE_IJMatrixSetDiagOffdSizes(matrix::HYPRE_IJMatrix,diag_sizes::Ptr{HYPRE_Int},offdiag_sizes::Ptr{HYPRE_Int})
    ccall((:HYPRE_IJMatrixSetDiagOffdSizes,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,Ptr{HYPRE_Int},Ptr{HYPRE_Int}),matrix,diag_sizes,offdiag_sizes)
end

function HYPRE_IJMatrixSetMaxOffProcElmts(matrix::HYPRE_IJMatrix,max_off_proc_elmts::HYPRE_Int)
    ccall((:HYPRE_IJMatrixSetMaxOffProcElmts,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,HYPRE_Int),matrix,max_off_proc_elmts)
end

function HYPRE_IJMatrixSetPrintLevel(matrix::HYPRE_IJMatrix,print_level::HYPRE_Int)
    ccall((:HYPRE_IJMatrixSetPrintLevel,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,HYPRE_Int),matrix,print_level)
end

function HYPRE_IJMatrixRead(filename::Ptr{UInt8},comm::Cint,_type::HYPRE_Int,matrix::Ptr{HYPRE_IJMatrix})
    ccall((:HYPRE_IJMatrixRead,HYPRE_IJ_mv),HYPRE_Int,(Ptr{UInt8},Cint,HYPRE_Int,Ptr{HYPRE_IJMatrix}),filename,comm,_type,matrix)
end

function HYPRE_IJMatrixPrint(matrix::HYPRE_IJMatrix,filename::Ptr{UInt8})
    ccall((:HYPRE_IJMatrixPrint,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJMatrix,Ptr{UInt8}),matrix,filename)
end

function HYPRE_IJVectorCreate(comm::Cint,jlower::HYPRE_Int,jupper::HYPRE_Int,vector::Ptr{HYPRE_IJVector})
    ccall((:HYPRE_IJVectorCreate,HYPRE_IJ_mv),HYPRE_Int,(Cint,HYPRE_Int,HYPRE_Int,Ptr{HYPRE_IJVector}),comm,jlower,jupper,vector)
end

function HYPRE_IJVectorDestroy(vector::HYPRE_IJVector)
    ccall((:HYPRE_IJVectorDestroy,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJVector,),vector)
end

function HYPRE_IJVectorInitialize(vector::HYPRE_IJVector)
    ccall((:HYPRE_IJVectorInitialize,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJVector,),vector)
end

function HYPRE_IJVectorSetMaxOffProcElmts(vector::HYPRE_IJVector,max_off_proc_elmts::HYPRE_Int)
    ccall((:HYPRE_IJVectorSetMaxOffProcElmts,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJVector,HYPRE_Int),vector,max_off_proc_elmts)
end

function HYPRE_IJVectorSetValues(vector::HYPRE_IJVector,nvalues::HYPRE_Int,indices::Ptr{HYPRE_Int},values::Ptr{Cdouble})
    ccall((:HYPRE_IJVectorSetValues,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJVector,HYPRE_Int,Ptr{HYPRE_Int},Ptr{Cdouble}),vector,nvalues,indices,values)
end

function HYPRE_IJVectorAddToValues(vector::HYPRE_IJVector,nvalues::HYPRE_Int,indices::Ptr{HYPRE_Int},values::Ptr{Cdouble})
    ccall((:HYPRE_IJVectorAddToValues,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJVector,HYPRE_Int,Ptr{HYPRE_Int},Ptr{Cdouble}),vector,nvalues,indices,values)
end

function HYPRE_IJVectorAssemble(vector::HYPRE_IJVector)
    ccall((:HYPRE_IJVectorAssemble,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJVector,),vector)
end

function HYPRE_IJVectorGetValues(vector::HYPRE_IJVector,nvalues::HYPRE_Int,indices::Ptr{HYPRE_Int},values::Ptr{Cdouble})
    ccall((:HYPRE_IJVectorGetValues,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJVector,HYPRE_Int,Ptr{HYPRE_Int},Ptr{Cdouble}),vector,nvalues,indices,values)
end

function HYPRE_IJVectorSetObjectType(vector::HYPRE_IJVector,_type::HYPRE_Int)
    ccall((:HYPRE_IJVectorSetObjectType,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJVector,HYPRE_Int),vector,_type)
end

function HYPRE_IJVectorGetObjectType(vector::HYPRE_IJVector,_type::Ptr{HYPRE_Int})
    ccall((:HYPRE_IJVectorGetObjectType,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJVector,Ptr{HYPRE_Int}),vector,_type)
end

function HYPRE_IJVectorGetLocalRange(vector::HYPRE_IJVector,jlower::Ptr{HYPRE_Int},jupper::Ptr{HYPRE_Int})
    ccall((:HYPRE_IJVectorGetLocalRange,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJVector,Ptr{HYPRE_Int},Ptr{HYPRE_Int}),vector,jlower,jupper)
end

function HYPRE_IJVectorGetObject(vector::HYPRE_IJVector,object::Ptr{Ptr{Void}})
    ccall((:HYPRE_IJVectorGetObject,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJVector,Ptr{Ptr{Void}}),vector,object)
end

function HYPRE_IJVectorSetPrintLevel(vector::HYPRE_IJVector,print_level::HYPRE_Int)
    ccall((:HYPRE_IJVectorSetPrintLevel,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJVector,HYPRE_Int),vector,print_level)
end

function HYPRE_IJVectorRead(filename::Ptr{UInt8},comm::Cint,_type::HYPRE_Int,vector::Ptr{HYPRE_IJVector})
    ccall((:HYPRE_IJVectorRead,HYPRE_IJ_mv),HYPRE_Int,(Ptr{UInt8},Cint,HYPRE_Int,Ptr{HYPRE_IJVector}),filename,comm,_type,vector)
end

function HYPRE_IJVectorPrint(vector::HYPRE_IJVector,filename::Ptr{UInt8})
    ccall((:HYPRE_IJVectorPrint,HYPRE_IJ_mv),HYPRE_Int,(HYPRE_IJVector,Ptr{UInt8}),vector,filename)
end


text:
function Makie.plot!(
        lscene::LScene, P::Makie.PlotFunc,
        attributes::Makie.Attributes, args...;
        kw_attributes...)
    # We store the show_axis attribute in the LScene
    if haskey(attributes, :show_axis)
        lscene.attributes[:show_axis] = pop!(attributes, :show_axis)
    end

    if haskey(attributes, :limits)
        lscene.attributes[:limits] = pop!(attributes, :limits)
    end

    show_axis = get!(lscene.attributes, :show_axis, true)
    plot = Makie.plot!(lscene.scene, P, attributes, args...; kw_attributes...)

    function get_lims()
        return get(lscene.attributes, :limits) do
            return data_limits(lscene.scene, Makie.isaxis)
        end
    end

    if isnothing(lscene.scene[OldAxis])
        # Add axis on first plot!, if requested
        to_value(show_axis) && Makie.axis3d!(lscene.scene, get_lims())
    else
        # Update limits when plotting new objects
        axis_plot = lscene.scene[OldAxis]
        axis_plot[1] = get_lims()
    end
    # Make sure axis is always in pos 1
    sort!(lscene.scene.plots, by=!Makie.isaxis)
    center!(lscene.scene)
    plot
end

function Makie.plot!(P::Makie.PlotFunc, ls::LScene, args...; kw_attributes...)
    attributes = Makie.Attributes(kw_attributes)
    Makie.plot!(ls, P, attributes, args...)
end

function layoutable(::Type{LScene}, fig_or_scene; bbox = nothing, scenekw = NamedTuple(), kwargs...)

    topscene = get_topscene(fig_or_scene)
    default_attrs = default_attributes(LScene, topscene).attributes
    theme_attrs = subtheme(topscene, :LScene)
    attrs = merge!(merge!(Attributes(kwargs), theme_attrs), default_attrs)
    layoutobservables = LayoutObservables{LScene}(attrs.width, attrs.height, attrs.tellwidth, attrs.tellheight,
        attrs.halign, attrs.valign, attrs.alignmode; suggestedbbox = bbox)
    # pick a camera and draw axis.
    scenekw = merge((clear = false, camera=cam3d!), scenekw)
    scene = Scene(topscene, lift(round_to_IRect2D, layoutobservables.computedbbox); scenekw...)
    return LScene(fig_or_scene, layoutobservables, attrs, Dict{Symbol, Any}(), scene)
end

function Base.delete!(ax::LScene, plot::AbstractPlot)
    delete!(ax.scene, plot)
    ax
end

can_be_current_axis(ls::LScene) = true


text:
export Group, terms, terms!, elements, elements!

const TermsDict = Dict{Symbol,Union{Node,Array{Node,1}}}
"""
    Group([elements [, terms[, units]]])
    Group(<keyword arguments>)

`Group` objects `show` FastHenry commands of their
elements.

**Keyword Arguments**

- `elements`    -- Array of `Element`'s
- `terms`       -- Dict{Symbol,Node} # connections to the Group
- `units`       -- Units used for the Group

Groups may be nested.  Automatic name generation will create a name for
each `Element` in the `Group` and all subgroups as needed when `show` is
called.  Generated names start with an underscore.  Do not use these names
for user named elements.

`getindex`, `setindex!`, and `merge!` operate on `terms`.
`push!`, `pop!`, `pushfirst!`, `popfirst!`, `append!`, `prepend!` operate on `elements`.
"""
mutable struct Group <: Element
  elements :: Array{Element,1}
  terms :: TermsDict
  units :: Units
end
# Group() = Group([],TermsDict(),Units())
Group(e) = Group(e,TermsDict(),Units())
Group(e,d) = Group(e,d,Units())
Group(;elements=[],terms=TermsDict(),units=Units())= Group(elements,terms,units)

units(g::Group) = g.units
function units!(g::Group, u::Units)
  g.units = u
  return nothing
end

elements(g::Group) = g.elements
function elements!(g::Group, e::Array{T,1}) where T<:Element
  g.elements = e
  return nothing
end
"""
    elements(g::Group)
    elements!(g::Group, e::Element)

Get and set `elements` field in a `Group`.
"""
elements, elements!

terms(g::Group) = g.terms
function terms!(g::Group, p)
  g.terms = p
  return nothing
end
"""
    terms(g::Group)
    terms!(g::Group, p::Dict{Symbol,Node})

Get and set the `terms` field in a `Group`.
"""
terms, terms!

Base.getindex(g::Group, key::Symbol) = g.terms[key]
function Base.setindex!(g::Group, n::Union{Node,Array{Node,1}}, key::Symbol)
  g.terms[key] = n
  return nothing
end
function Base.merge!(g::Group, args...)
  for arg in args
    merge!(g.terms, arg.terms)
  end
end
Base.push!(g::Group, args...) = push!(g.elements, args...)
Base.pop!(g::Group) = pop!(g.elements)
Base.pushfirst!(g::Group, args...) = pushfirst!(g.elements, args...)
Base.popfirst!(g::Group) = popfirst!(g.elements)
Base.append!(g1::Group, g2::Group) = append!(g1.elements, g2.elements)
Base.prepend!(g1::Group, g2::Group) = prepend!(g1.elements, g2.elements)


text:
@testset "JDF.jl eachcol" begin
    df = DataFrame([collect(1:100) for i =1:3000])
    df[!, :int_missing] =
        rand(
            [rand(rand([UInt, Int, Float64, Float32, Bool])), missing],
            nrow(df)
        )

    df[!, :missing] .= missing
    df[!, :strs] = [randstring(8) for i in 1:nrow(df)]
    df[!, :stringarray] = StringVector([randstring(8) for i in 1:nrow(df)])

    df[!, :strs_missing] = [rand([missing, randstring(8)]) for i in 1:nrow(df)]
    df[!, :stringarray_missing] = StringVector([rand([missing, randstring(8)]) for i in 1:nrow(df)])
    df[!, :symbol_missing] = [rand([missing, Symbol(randstring(8))]) for i in 1:nrow(df)]
    df[!, :char] = getindex.(df[!, :strs], 1)
    df[!, :char_missing] = allowmissing(df[!, :char])
    df[rand(1:nrow(df), 10), :char_missing] .= missing

    savejdf("a.jdf", df)

    df2 = jdf"a.jdf"

    @test ncol(df2) == 3009
    @test nrow(df2) == 100

    df3 = [a for a in eachcol(df2)]

    df4 = DataFrame(df3)

    @test size(df4) == size(df)
    @test all([isequal(df4[!, n], df[!, n]) for n in 1:ncol(df4)])

    # clean up
    rm("a.jdf", force=true, recursive=true)
end


text:
module HardSphereDynamics

export FixedPlane, MovableBall, HardSphereFluid
export initial_condition!, evolve!
export FreeFlow, ExternalFieldFlow
export ElasticCollision
export AllToAll
export HardSphereSimulation

# stdlib:
using LinearAlgebra

# libraries:
using StaticArrays
using Parameters
using Distributions

using Requires

function __init__()
    @require AbstractPlotting="537997a7-5e4e-5d89-9595-2241ea00577e" begin
        include("visualization/makie_visualization.jl")
    end
end

normsq(v) = sum(abs2, v)

# types:
include("ball.jl")
include("box.jl")
include("hard_sphere_fluid.jl")

# flow dynamics:
include("flow.jl")
include("free_flow.jl")
include("external_field_flow.jl")

# simulation:
include("collisions.jl")
include("event_handler.jl")
include("simulation.jl")

end # module

max_stars_repo_name: dpsanders/HardSphereDynamics

max_forks_repo_name: dpsanders/HardSphereDynamics