text:
"""
	BINS
"""
# ( min, max, # bins - 1 )
function bins( inputs::Tuple{<:Real,<:Real,<:Real} )
    min, max, divs = convert.( Float32, inputs )
    return bins( min, max, (max-min)/divs == 0 ? 1.0 : (max-min)/divs  )
end

# ( min, max, step )
bins( min, max, step ) = bins( Float32(min), Float32(max), Float32(step)  )
bins( min::Float32, max::Float32, step::Float32 ) = collect( min:step:max )

# range
bins( range::AbstractRange ) = Float32.( collect( range ) )

"""
	HISTOGRAMS
"""
function histogram( image, numbins::Integer )
	min, max = extrema( image );
	return histogram( image, Float32.((min,max,numbins)) )
end

function histogram( image, bin_data::Union{Tuple{<:Real,<:Real,<:Real}, AbstractRange }, T=Int32 )

    bin_intervals = bins( bin_data );
    histogram = zeros( T, length(bin_intervals) - 1 );
    outs = 0; # keeps count of pixels that don't fall into any interval

    @inbounds for pixel in image
        idx = binaryBinSearch( bin_intervals, pixel )

        if idx == 0
            outs += 1;
        else
            histogram[idx] += T(1);
        end
    end

    return histogram, outs, bin_intervals
end

function normhistogram( image, numbins::Integer )
	min, max = extrema( image );
	return normhistogram( image, Float32.((min,max,numbins)) )
end

function normhistogram( image, bin_data::Union{Tuple{<:Real,<:Real,<:Real}, AbstractRange } )

	h, o, b = histogram( image, bin_data, Float32 )

	N = length( image );
	@inbounds @simd for e in 1:length(h)
		h[e] /= N
	end

    return h, o, b
end

# Cumulative density function from histogram

function cdf( image, bin_data::Union{Tuple{T1,T2,T3}, AbstractRange{T4}} ) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real}

    bin_intervals = bins( bin_data );
    cdf = zeros( Float32, length(bin_intervals) - 1 );
    onef = Float32(1.0);
    outs = 0;

    @inbounds for pixel in image
        idx = binaryBinSearch( bin_intervals, pixel )
        if idx == 0
            outs += 1;
        else
            cdf[idx] += onef;
        end
    end

	N = length( image )
	cdf[1] = cdf[1]/N
	@inbounds for e in 2:length(cdf)
		cdf[e] = cdf[e]/N + cdf[e-1]
	end

    return cdf, outs, bin_intervals
end

# Histogram image equalization

function eqim( img, bin_data::Union{Tuple{T1,T2,T3}, AbstractRange{T4}} ) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real}

	bs   = bins( bin_data );
	_cdf, _, _ = cdf( img, bin_data );
    eqim = copy( img );

    for e in 1:length(img)
        pixel = img[e]
        for idx in 1:length(bs)-1
            if pixel > bs[idx] && pixel <= bs[idx+1]
                eqim[e] = floor( eltype(img), length(bs)*_cdf[idx] );
            end
        end

    end
    return eqim
end

function eqim( img::Array{<:AbstractFloat,2}, bin_data::Union{Tuple{T1,T2,T3}, AbstractRange{T4}} ) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real}

	bs   = bins( bin_data );
	_cdf, _, _ = cdf( img, bin_data );
    eqim = copy( img );

    for e in 1:length(img)
        pixel = img[e]
        for idx in 1:length(bs)-1
            if pixel > bs[idx] && pixel <= bs[idx+1]
                eqim[e] = floor( length(bs)*_cdf[idx] );
            end
        end

    end
    return eqim
end

# Otsu thresholding
# Does it only work on integers...?

function o_varB( k, sumP, sumIP, meanT )
	return ( meanT*sumP[k] - sumIP[k] )^2/( sumP[k]*( 1 - sumP[k] ) )
end

function otsu( p::Array{Float32,1} )

	L = length( p )
	sumP  = copy( p )
	sumIP = copy( p )
	@simd for i in 2:L
		 @inbounds sumP[i] += sumP[ i - 1 ]
	end
	@simd for i in 2:L
		 @inbounds sumIP[i] = i*p[i] + sumIP[ i - 1 ]
	end
	meanT = sumIP[end]

	maxK = 2
	maxvarB = o_varB( maxK, sumP, sumIP, meanT )
	for k in 3:L-1
		varB = o_varB( k, sumP, sumIP, meanT )
		if ( varB > maxvarB )
			maxK = k
			maxvarB = varB
		end
	end

	return maxK
end

function otsu( img, bin_data::Union{Tuple{T1,T2,T3}, AbstractRange{T4}} ) where {T1<:Real,T2<:Real,T3<:Real,T4<:Real}

	nh, _, _ = normhistogram( img, bin_data )
	ot = otsu( nh );

	println( "Otsu threshold = ", ot );
	mask = zeros( Bool, size(img) );

	for e in 1:length(img)
		mask[e] = img[e] > ot
	end

	return mask
end

function otsu( img::Array{<:Real,N} ) where {N}

	mn, mx = extrema( img );
	bin_data = mn:(mx-mn)/100:mx;
	nh, _, _ = normhistogram( img, bin_data )
	ot = otsu( nh );

	println( "Otsu threshold = ", ot );
	mask = zeros( UInt8, size(img) );

	for e in 1:length(img)
		mask[e] = reinterpret( UInt8, img[e] > ot )
	end

	return mask
end

function otsuArr( p::Array{Float32,1} )

	arr   = copy( p )

	L = length( p )
	sumP  = copy( p )
	sumIP = copy( p )
	@simd for i in 2:L
		 @inbounds sumP[i] += sumP[ i - 1 ]
	end
	@simd for i in 2:L
		 @inbounds sumIP[i] = i*p[i] + sumIP[ i - 1 ]
	end
	meanT = sumIP[end]

	for k in 1:L
		arr[k] = o_varB( k, sumP, sumIP, meanT )
	end

	return arr
end


text:
#############
# Promotion #
#############

# Promotion between two tensors promote the eltype and promotes
# symmetric tensors to tensors

function Base.promote_rule{dim , A <: Number, B <: Number, order, M}(::Type{SymmetricTensor{order, dim, A, M}},
                                                                     ::Type{SymmetricTensor{order, dim, B, M}})
    SymmetricTensor{order, dim, promote_type(A, B), M}
end

function Base.promote_rule{dim , A <: Number, B <: Number, order, M}(::Type{Tensor{order, dim, A, M}},
                                                                     ::Type{Tensor{order, dim, B, M}})
    Tensor{order, dim, promote_type(A, B), M}
end

function Base.promote_rule{dim , A <: Number, B <: Number, order, M1, M2}(::Type{SymmetricTensor{order, dim, A, M1}},
                                                                          ::Type{Tensor{order, dim, B, M2}})
    Tensor{order, dim, promote_type(A, B), M2}
end

function Base.promote_rule{dim , A <: Number, B <: Number, order, M1, M2}(::Type{Tensor{order, dim, A, M1}},
                                                                          ::Type{SymmetricTensor{order, dim, B, M2}})
    Tensor{order, dim, promote_type(A, B), M1}
end

###############
# Conversions #
###############

# Identity conversions
@inline Base.convert{order, dim, T}(::Type{Tensor{order, dim, T}}, t::Tensor{order, dim, T}) = t
@inline Base.convert{order, dim, T, M}(::Type{Tensor{order, dim, T, M}}, t::Tensor{order, dim, T, M}) = t
@inline Base.convert{order, dim, T}(::Type{SymmetricTensor{order, dim, T}}, t::SymmetricTensor{order, dim, T}) = t
@inline Base.convert{order, dim, T, M}(::Type{SymmetricTensor{order, dim, T, M}}, t::SymmetricTensor{order, dim, T, M}) = t

# Change element type
@inline Base.convert{order, dim, T1}(::Type{Tensor{order, dim, T1}}, t::Tensor{order, dim}) = Tensor{order, dim}(convert(SVector{n_components(Tensor{order, dim}), T1}, tovector(t)))
@inline Base.convert{order, dim, T1}(::Type{SymmetricTensor{order, dim, T1}}, t::SymmetricTensor{order, dim}) = SymmetricTensor{order, dim}(convert(SVector{n_components(SymmetricTensor{order, dim}), T1}, tovector(t)))

# Peel off the M but define these so that convert(typeof(...), ...) works
@inline Base.convert{order, dim, T1, M}(::Type{Tensor{order, dim, T1, M}}, t::Tensor{order, dim}) = convert(Tensor{order, dim, T1}, t)
@inline Base.convert{order, dim, T1, M}(::Type{SymmetricTensor{order, dim, T1, M}}, t::SymmetricTensor{order, dim}) = convert(SymmetricTensor{order, dim, T1}, t)
@inline Base.convert{order, dim, T1, M}(::Type{Tensor{order, dim, T1, M}}, t::SymmetricTensor{order, dim}) = convert(Tensor{order, dim, T1}, t)
@inline Base.convert{order, dim, T1, M}(::Type{SymmetricTensor{order, dim, T1, M}}, t::Tensor{order, dim}) = convert(SymmetricTensor{order, dim, T1}, t)

@inline Base.convert{order, dim, T}(::Type{Tensor{order, dim}}, t::SymmetricTensor{order, dim, T}) = convert(Tensor{order, dim, T}, t)
@inline Base.convert{order, dim, T}(::Type{SymmetricTensor{order, dim}}, t::Tensor{order, dim, T}) = convert(SymmetricTensor{order, dim, T}, t)
@inline Base.convert{order, dim, T}(::Type{Tensor}, t::SymmetricTensor{order, dim, T}) = convert(Tensor{order, dim, T}, t)
@inline Base.convert{order, dim, T}(::Type{SymmetricTensor}, t::Tensor{order, dim, T}) = convert(SymmetricTensor{order, dim, T}, t)

# SymmetricTensor -> Tensor
# We unroll the creation by calling the compute_index function
@generated function Base.convert{order, dim, T1, T2}(::Type{Tensor{order, dim, T1}}, t::SymmetricTensor{order, dim, T2})
    exp = Expr(:tuple)
    # Compute (row, col) from linear index
    if order == 2
        for j in 1:dim, i in 1:dim
            push!(exp.args, :(T1(data[$(compute_index(SymmetricTensor{order, dim}, i, j))])))
        end
    else
        for l in 1:dim, k in 1:dim, j in 1:dim, i in 1:dim
            push!(exp.args, :(T1(data[$(compute_index(SymmetricTensor{order, dim}, i, j, k, l))])))
        end
    end
    return quote
        $(Expr(:meta, :inline))
        data = get_data(t)
        v = $exp
        Tensor{order, dim}(v)
    end
end

# Tensor -> SymmetricTensor
function Base.convert{dim, order, T1}(::Type{SymmetricTensor{order, dim, T1}}, t::Tensor{order, dim})
    if issymmetric(t)
        return convert(SymmetricTensor{order, dim, T1}, symmetric(t))
    else
        throw(InexactError())
    end
end

ext: jl


text:
function _check_hdf5_shapes(shapes)
    nobs = map(last, filter(!isempty, shapes))

    return all(==(first(nobs)), nobs[2:end])
end

"""
    HDF5Dataset(file::AbstractString, paths)
    HDF5Dataset(fid::HDF5.File, paths::Union{HDF5.Dataset, Vector{HDF5.Dataset}})
    HDF5Dataset(fid::HDF5.File, paths::Union{AbstractString, Vector{<:AbstractString}})
    HDF5Dataset(fid::HDF5.File, paths::Union{HDF5.Dataset, Vector{HDF5.Dataset}}, shapes)

Wrap several HDF5 datasets (`paths`) as a single dataset container.
Each dataset `p` in `paths` should be accessible as `fid[p]`.
Calling `getobs` on a `HDF5Dataset` returns a tuple with each element corresponding
to the observation from each dataset in `paths`.
See [`close(::HDF5Dataset)`](@ref) for closing the underlying HDF5 file pointer.

For array datasets, the last dimension is assumed to be the observation dimension.
For scalar datasets, the stored value is returned by `getobs` for any index.
"""
struct HDF5Dataset{T<:Union{HDF5.Dataset, Vector{HDF5.Dataset}}} <: AbstractDataContainer
    fid::HDF5.File
    paths::T
    shapes::Vector{Tuple}

    function HDF5Dataset(fid::HDF5.File, paths::T, shapes::Vector) where T<:Union{HDF5.Dataset, Vector{HDF5.Dataset}}
        _check_hdf5_shapes(shapes) ||
            throw(ArgumentError("Cannot create HDF5Dataset for datasets with mismatched number of observations."))

        new{T}(fid, paths, shapes)
    end
end

HDF5Dataset(fid::HDF5.File, path::HDF5.Dataset) = HDF5Dataset(fid, path, [size(path)])
HDF5Dataset(fid::HDF5.File, paths::Vector{HDF5.Dataset}) =
    HDF5Dataset(fid, paths, map(size, paths))
HDF5Dataset(fid::HDF5.File, path::AbstractString) = HDF5Dataset(fid, fid[path])
HDF5Dataset(fid::HDF5.File, paths::Vector{<:AbstractString}) =
    HDF5Dataset(fid, map(p -> fid[p], paths))
HDF5Dataset(file::AbstractString, paths) = HDF5Dataset(h5open(file, "r"), paths)

_getobs_hdf5(dataset::HDF5.Dataset, ::Tuple{}, i) = read(dataset)
function _getobs_hdf5(dataset::HDF5.Dataset, shape, i)
    I = map(s -> 1:s, shape[1:(end - 1)])

    return dataset[I..., i]
end
Base.getindex(dataset::HDF5Dataset{HDF5.Dataset}, i) =
    _getobs_hdf5(dataset.paths, only(dataset.shapes), i)
Base.getindex(dataset::HDF5Dataset{<:Vector}, i) =
    Tuple(map((p, s) -> _getobs_hdf5(p, s, i), dataset.paths, dataset.shapes))
Base.length(dataset::HDF5Dataset) = last(first(filter(!isempty, dataset.shapes)))

"""
    close(dataset::HDF5Dataset)

Close the underlying HDF5 file pointer for `dataset`.
"""
Base.close(dataset::HDF5Dataset) = close(dataset.fid)


text:
module FESpacesInterfacesTests

using FillArrays
using Test
using Gridap.Arrays
using Gridap.Fields
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.FESpaces
using Gridap.CellData
using LinearAlgebra

order = 1
domain =(0,1,0,1)
partition = (3,3)
model = CartesianDiscreteModel(domain,partition)
reffe = ReferenceFE(lagrangian,Float64,order)
V = FESpace(model,reffe,dirichlet_tags=["tag_1","tag_6"])
test_fe_space(V)

vh = FEFunction(V,rand(num_free_dofs(V)))
@test isa(vh,FEFunction)
test_fe_function(vh)

dv = get_fe_basis(V)
du = get_trial_fe_basis(V)

Ω = Triangulation(model)
∂Ω = BoundaryTriangulation(model)
trian_Γ = SkeletonTriangulation(model)
x_Γ = get_cell_points(trian_Γ)

@test isa(dv.minus(x_Γ)[1],ArrayBlock)
@test isa(du.plus(x_Γ)[1],ArrayBlock)
@test isa(∇(dv).plus(x_Γ)[1],ArrayBlock)
@test isa(∇(du).minus(x_Γ)[1],ArrayBlock)

glue = get_glue(∂Ω,Val(num_cell_dims(model)))
cellids = glue.tface_to_mface
cell_vals = get_cell_dof_values(vh,∂Ω)
@test cell_vals == get_cell_dof_values(vh)[cellids]

cell_vals = get_cell_dof_values(vh,trian_Γ)
@test isa(cell_vals[1],ArrayBlock)

zh = zero(V)
@test isa(zh,FEFunction)
@test has_constraints(V) == false

cellids = [1,3,5,2]
Ω1 = Triangulation(model,cellids)
@test get_cell_dof_ids(V,Ω1) == [[-1, 1, 4, 5], [2, 3, 6, 7], [5, 6, 9, 10], [1, 2, 5, 6]]

@test isa(get_cell_dof_ids(V,trian_Γ)[1],ArrayBlock)
@test get_cell_dof_ids(V,trian_Γ)[1][1] == [-1, 1, 4, 5]
@test get_cell_dof_ids(V,trian_Γ)[1][2] == [4, 5, 8, 9]

cell_constr = get_cell_constraints(V)
@test cell_constr == [Matrix(I,4,4) for cell in 1:num_cells(model)]

cell_constr = get_cell_constraints(V,Ω1)
@test cell_constr == [Matrix(I,4,4) for cell in 1:length(cellids)]

@test get_cell_isconstrained(V,Ω1) == Fill(false,length(cellids))

cell_constr = get_cell_constraints(V,trian_Γ)
@test isa(cell_constr[1],ArrayBlock)

du = get_trial_fe_basis(V)
du_data = get_data(du)
@test size(du_data[1]) == (1,4)

@test get_cell_isconstrained(V,trian_Γ) == Fill(false,num_cells(trian_Γ))

cellmat = [rand(4,4) for cell in 1:num_cells(model)]
cellvec = [rand(4) for cell in 1:num_cells(model)]
cellmatvec = pair_arrays(cellmat,cellvec)

@test cellmat === attach_constraints_rows(V,cellmat,Ω)
@test cellmat === attach_constraints_cols(V,cellmat,Ω)
test_fe_space(V,cellmatvec,cellmat,cellvec,Ω)

end # module


text:
# cache such that creating svgs from LaTeX don't need to be created every time
# this is also used for test cases such that `tex2svg` doesn't need to be installed on Github Objects
include("latexsvgfile.jl")
latex(text::LaTeXString) = latex(text, O)
latex(text::LaTeXString, pos::Point, valign::Symbol, halign::Symbol) =
    latex(text, pos, :stroke, valign = valign, halign = halign)
latex(text::LaTeXString, x, y) = latex(text, Point(x, y), :stroke)
latex(text::LaTeXString, x, y, valign::Symbol, halign::Symbol) =
    latex(text, Point(x, y), :stroke, valign = valign, halign = halign)
latex(text::LaTeXString, pos::Point; valign = :top, halign = :left) =
    latex(text, pos, :stroke, valign = valign, halign = halign)

"""
    latex(text::LaTeXString, pos::Point, object::Symbol; valign = :top, halign = :left)

Add the latex string `text` to the top left corner of the LaTeX path.
Can be added to `Luxor.jl` graphics via [`Video`](@ref).

**NOTES:**
- **This only works if `tex2svg` is installed.**
    It can be installed using the following command (you may have to prefix this command with `sudo` depending on your installation):

        npm install -g mathjax-node-cli

- **The `latex` method must be called from within an [`Object`](@ref).**

# Arguments
- `text::LaTeXString`: a LaTeX string to render.
- `pos::Point`: position of the upper left corner of the latex text. Default: `O`
    - can be written as `x, y` instead of `Point(x, y)`
- `object::Symbol`: graphics objects defined by `Luxor.jl`. Default `:stroke`.
Available objects:
  - `:stroke` - Draws the latex string on the canvas. For more info check `Luxor.strokepath`
  - `:path` - Creates the path of the latex string but does not render it to the canvas.

# Keywords:
  - `valign::Symbol=:top`: vertical alignment with respect to the specified `pos` parameter.
      - Options available are `:top`, `:middle`, `:bottom`
  - `halign::Symbol=:left`: horizontal alignment with respect to the specified `pos` parameter.
      - Options available are `:left`, `:center/:centre`, `:right`

# Throws
- `IOError`: mathjax-node-cli is not installed

# Warning
Shows a warning if either of the alignment options are unrecognised.

# Example

```julia
using Javis
using LaTeXStrings

function ground(args...)
    background("white")
    sethue("black")
end

function draw_latex(video, object, frame)
    fontsize(50)
    x = 100
    y = 120
    latex(L"\\sqrt{5}", x, y)
end

demo = Video(500, 500)
Background(1:2, ground)
Object(draw_latex)
render(demo; pathname = "latex.gif")
```

"""
function latex(
    text::LaTeXString,
    pos::Point,
    draw_object::Symbol;
    valign = :top,
    halign = :left,
)
    object = CURRENT_OBJECT[1]
    opts = object.opts
    t = get(opts, :draw_text_t, 1.0)
    return animate_latex(text, pos, t, valign, halign, draw_object)
end

function animate_latex(text, pos::Point, t, valign::Symbol, halign::Symbol, object)
    svg = get_latex_svg(text)
    object == :stroke && (object = :fill)

    w, h = svgwh(svg)
    halignment = findfirst(isequal(halign), (:left, :center, :right, :centre))

    if halignment === nothing
        @warn "Unknown horizontal alignment option: $(halign). Defaulting to left alignment"
        halignment = 1
    elseif halignment == 4
        halignment = 2
    end

    textpointx = pos.x - (0, w / 2, w)[halignment]
    valignment = findfirst(isequal(valign), (:top, :middle, :bottom))

    if valignment === nothing
        @warn "Unknown vertical alignment option: $(valign). Defaulting to top alignment"
        valignment = 1
    end

    textpointy = pos.y - (0, h / 2, h)[valignment]
    upperleft = Point(textpointx, textpointy)

    if t >= 1
        translate(upperleft)
        pathsvg(svg)
        do_action(object)
        translate(-upperleft)
        return
    end

    translate(upperleft)
    pathsvg(svg)
    do_action(:clip)
    r = t * sqrt(w^2 + h^2)
    circle(O, r, :fill)
    translate(-upperleft)
end
"""
    strip_eq(text::LaTeXString)

Strips `\$\$` from `text.s` if present and returns the resulting string.

# Arguments
- `text::LaTeXString`: a LaTeX string
"""
function strip_eq(text::LaTeXString)
    ts = text.s
    if ts[1] == '$'
        ts = ts[2:(end - 1)]
    end
    ts
end

# \todo update LaTeXSVG cache to use output of strip_eq as the key. See https://github.com/JuliaAnimators/Javis.jl/pull/307#issuecomment-749616375
function get_latex_svg(text::LaTeXString)
    # check if it's cached
    if haskey(LaTeXSVG, text)
        svg = LaTeXSVG[text]
    else
        ts = replace(strip_eq(text), "\n" => " ")
        command = if Sys.iswindows()
            `cmd /C tex2svg $ts`
        else
            `tex2svg $ts`
        end
        try
            svg = read(command, String)
        catch e
            @warn "Using LaTeX needs the program `tex2svg` which might not be installed"
            @info "It can be installed using `npm install -g mathjax-node-cli`"
            throw(e)
        end
        LaTeXSVG[text] = svg
    end
    return svg
end


text:
"""
    LayerProvider

A `LayerProvider` is an abstract type used to dispatch the correct call of
`SimpleSDMPredictor` to a specific dataset. A dataset is specified by a
`LayerProvider` and a `LayerDataset`, as well as optionally one or multiple
layers, and future climate information, resolution, or dates.
"""
abstract type LayerProvider end

"""
    LayerDataset

A `LayerDataset` is a specific set of rasters provided by a `LayerProvider`. For
a number of dataset types that are very broad (`LandCover`,
`HabitatHeterogeneity`), the precise mapping of layers is documented in their
`SimpleSDMPredictor` method.
"""
abstract type LayerDataset end

"""
    WorldClim

[WorldClim](https://worldclim.org/) offers bioclimatic data both historical and
future, under `CMIP6` scenarios.

This provider currently offers `BioClim` data, both historical and future under
`CMIP6`, and the `Elevation` raster.
"""
struct WorldClim <: LayerProvider end

"""
    CHELSA

[CHELSA](https://chelsa-climate.org/) offers high resolution climatologies. This
provider currently offers `BioClim` data, both historical and future under
`CMIP5` *and* `CMIP6`.

Note that CHELSA offers a subset of all possible CMIP6 combinations, which is
supposed to be the most informative.
"""
struct CHELSA <: LayerProvider end

"""
    EarthEnv

Data from the [EarthEnv](https://www.earthenv.org//) project, all released under
a CC-BY-NC license to Tuanmu & Jetz. This provider currently offers `LandCover`
and `HabitatHeterogeneity` rasters.
"""
struct EarthEnv <: LayerProvider end

"""
    BioClim

A list of 19 bioclimatic variables derived from the monthly temperature and
precipitation data. This dataset is provided by `WorldClim` and `CHELSA`, both
of which offer future versions under `CMIP5` and `CMIP6` models.

| Variable | Description                                                |
| ------   | ------                                                     |
| 1        | Annual Mean Temperature                                    |
| 2        | Mean Diurnal Range (Mean of monthly (max temp - min temp)) |
| 3        | Isothermality (BIO2/BIO7) (* 100)                          |
| 4        | Temperature Seasonality (standard deviation *100)          |
| 5        | Max Temperature of Warmest Month                           |
| 6        | Min Temperature of Coldest Month                           |
| 7        | Temperature Annual Range (BIO5-BIO6)                       |
| 8        | Mean Temperature of Wettest Quarter                        |
| 9        | Mean Temperature of Driest Quarter                         |
| 10       | Mean Temperature of Warmest Quarter                        |
| 11       | Mean Temperature of Coldest Quarter                        |
| 12       | Annual Precipitation                                       |
| 13       | Precipitation of Wettest Month                             |
| 14       | Precipitation of Driest Month                              |
| 15       | Precipitation Seasonality (Coefficient of Variation)       |
| 16       | Precipitation of Wettest Quarter                           |
| 17       | Precipitation of Driest Quarter                            |
| 18       | Precipitation of Warmest Quarter                           |
| 19       | Precipitation of Coldest Quarter                           |
"""
struct BioClim <: LayerDataset end

"""
    LandCover

Information on land cover, currently only provided by `EarthEnv`.
"""
struct LandCover <: LayerDataset end

"""
    Elevation

General type for a DEM, currently available through `WorldClim`
"""
struct Elevation <: LayerDataset end

"""
    HabitatHeterogeneity

Information on [habitat heterogeneity](https://www.earthenv.org/texture),
currently only provided by `EarthEnv`.
"""
struct HabitatHeterogeneity <: LayerDataset end

"""
    Topography

Information on habitat topography, currently provided by `EarthEnv`.
"""
struct Topography <: LayerDataset end

"""
    SharedSocioeconomicPathway

Enumeration of the four SSPs, which can be listed with
`instances(SharedSocioeconomicPathway)`. These are meant to be used with `CMIP6`
models.
"""
@enum SharedSocioeconomicPathway SSP126 SSP245 SSP370 SSP585

"""
    RepresentativeConcentrationPathway

Enumeration of the four RCPs, which can be listed with
`instances(RepresentativeConcentrationPathway)`. These are meant to be used with
`CMIP5` models.
"""
@enum RepresentativeConcentrationPathway RCP26 RCP45 RCP60 RCP85

"""
    CMIP6

Enumeration of the models from CMIP6, which can be listed with
`instances(CMIP6)`.
"""
@enum CMIP6 BCCCSM2MR CNRMCMS61 CNRMESM21 CanESM5 GFDLESM4 IPSLCM6ALR MIROCES2L MIROC6 MRIESM2

"""
    CMIP5

Enumeration of the models from CMIP5, which can be listed with
`instances(CMIP5)`.
"""
@enum CMIP5 ACCESS10 BNUESM CCSM4 CESM1BGC CESM1CAM5 CMCCCMS CMCCCM CNRMCM5 CSIROMK360 CanESM2 FGOALSG2 FIOESM GFDLCM3 GFDLESM2G GFDLESM2M GISSE2HCC GISSE2H GISSE2RCC GISSE2R HADGEM2AO HADGEM2CC IPSLCM5ALR IPSLCM5AMR MIROCESMCHEM MIROCESM MIROC5 MPIESMLR MPIESMMR MRICGCM3 MRIESM1 NORESM1M BCCCSM11 INMCM4

# Provider paths
_rasterpath(::Type{WorldClim}) = "WorldClim"
_rasterpath(::Type{CHELSA}) = "CHELSA"
_rasterpath(::Type{EarthEnv}) = "EarthEnv"

# Dataset paths
_rasterpath(::Type{BioClim}) = "BioClim"
_rasterpath(::Type{LandCover}) = "LandCover"
_rasterpath(::Type{HabitatHeterogeneity}) = "HabitatHeterogeneity"
_rasterpath(::Type{Elevation}) = "Elevation"
_rasterpath(::Type{Topography}) = "Topography"

# Future paths
_rasterpath(model::CMIP6) = _rasterpath(Val{model})
_rasterpath(model::CMIP5) = _rasterpath(Val{model})
_rasterpath(ssp::SharedSocioeconomicPathway) = _rasterpath(Val{ssp})
_rasterpath(rcp::RepresentativeConcentrationPathway) = _rasterpath(Val{rcp})

# SSP path
_rasterpath(::Type{Val{SSP126}}) = "ssp126"
_rasterpath(::Type{Val{SSP245}}) = "ssp245"
_rasterpath(::Type{Val{SSP370}}) = "ssp370"
_rasterpath(::Type{Val{SSP585}}) = "ssp585"

# RCP path
_rasterpath(::Type{Val{RCP26}}) = "rcp26"
_rasterpath(::Type{Val{RCP45}}) = "rcp45"
_rasterpath(::Type{Val{RCP60}}) = "rcp60"
_rasterpath(::Type{Val{RCP85}}) = "rcp85"

# CMIP6 path
_rasterpath(::Type{Val{BCCCSM2MR}}) = "BCC-CSM2-MR"
_rasterpath(::Type{Val{CNRMCMS61}}) = "CNRM-CM6-1"
_rasterpath(::Type{Val{CNRMESM21}}) = "CNRM-ESM2-1"
_rasterpath(::Type{Val{CanESM5}}) = "CanESM5"
_rasterpath(::Type{Val{GFDLESM4}}) = "GFDL-ESM4"
_rasterpath(::Type{Val{IPSLCM6ALR}}) = "IPSL-CM6A-LR"
_rasterpath(::Type{Val{MIROCES2L}}) = "MIROC-ES2L"
_rasterpath(::Type{Val{MIROC6}}) = "MIROC6"
_rasterpath(::Type{Val{MRIESM2}}) = "MRI-ESM2-0"

# CMIP5 path
_rasterpath(::Type{Val{ACCESS10}}) = "ACCESS1-0"
_rasterpath(::Type{Val{BNUESM}}) = "BNU-ESM"
_rasterpath(::Type{Val{CCSM4}}) = "CCSM4"
_rasterpath(::Type{Val{CESM1BGC}}) = "CESM1-BGC"
_rasterpath(::Type{Val{CESM1CAM5}}) = "CESM1-CAM5"
_rasterpath(::Type{Val{CMCCCMS}}) = "CMCC-CMS"
_rasterpath(::Type{Val{CMCCCM}}) = "CMCC-CM"
_rasterpath(::Type{Val{CNRMCM5}}) = "CNRM-CM5"
_rasterpath(::Type{Val{CSIROMK360}}) = "CSIRO-Mk3-6-0"
_rasterpath(::Type{Val{CanESM2}}) = "CanESM2"
_rasterpath(::Type{Val{FGOALSG2}}) = "FGOALS-g2"
_rasterpath(::Type{Val{FIOESM}}) = "FIO-ESM"
_rasterpath(::Type{Val{GFDLCM3}}) = "GFDL-CM3"
_rasterpath(::Type{Val{GFDLESM2G}}) = "GFDL-ESM2G"
_rasterpath(::Type{Val{GFDLESM2M}}) = "GFDL-ESM2M"
_rasterpath(::Type{Val{GISSE2HCC}}) = "GISS-E2-H-CC"
_rasterpath(::Type{Val{GISSE2H}}) = "GISS-E2-H"
_rasterpath(::Type{Val{GISSE2RCC}}) = "GISS-E2-R-CC"
_rasterpath(::Type{Val{GISSE2R}}) = "GISS-E2-R"
_rasterpath(::Type{Val{HADGEM2AO}}) = "HadGEM2-AO"
_rasterpath(::Type{Val{HADGEM2CC}}) = "HadGEM2-CC"
_rasterpath(::Type{Val{IPSLCM5ALR}}) = "IPSL-CM5A-LR"
_rasterpath(::Type{Val{IPSLCM5AMR}}) = "IPSL-CM5A-MR"
_rasterpath(::Type{Val{MIROCESMCHEM}}) = "MIROC-ESM-CHEM"
_rasterpath(::Type{Val{MIROCESM}}) = "MIROC-ESM"
_rasterpath(::Type{Val{MIROC5}}) = "MIROC5"
_rasterpath(::Type{Val{MPIESMLR}}) = "MPI-ESM-LR"
_rasterpath(::Type{Val{MPIESMMR}}) = "MPI-ESM-MR"
_rasterpath(::Type{Val{MRICGCM3}}) = "MRI-CGCM3"
_rasterpath(::Type{Val{MRIESM1}}) = "MRI-ESM1"
_rasterpath(::Type{Val{NORESM1M}}) = "NorESM1-M"
_rasterpath(::Type{Val{BCCCSM11}}) = "bcc-csm1-1"
_rasterpath(::Type{Val{INMCM4}}) = "inmcm4"


text:
struct AdditiveZeroArray{T,N} <: AbstractArray{T,N} end

Base.:(+)(lhs::AbstractArray{T,N}, rhs::AdditiveZeroArray{T,N}) where {T,N} = lhs
Base.:(+)(lhs::AdditiveZeroArray{T,N}, rhs::AbstractArray{T,N}) where {T,N} = rhs
Base.:(-)(lhs::AbstractArray{T,N}, rhs::AdditiveZeroArray{T,N}) where {T,N} = lhs
Base.:(-)(lhs::AdditiveZeroArray{T,N}, rhs::AbstractArray{T,N}) where {T,N} = -rhs
Base.:(*)(lhs::Number, rhs::AdditiveZeroArray{T,N}) where {T,N} = rhs
Base.:(*)(lhs::AdditiveZeroArray{T,N}, rhs::Number) where {T,N} = lhs


text:
using Oceananigans.Diagnostics: AdvectiveCFL, DiffusiveCFL
using Oceananigans.Simulations: TimeStepWizard
using Oceananigans.Utils: prettytime
using Printf

export SimpleProgressMessenger, SingleLineProgressMessenger, TimedProgressMessenger

get_Δt(Δt) = Δt
get_Δt(wizard::TimeStepWizard) = wizard.Δt

function SimpleProgressMessenger_function(simulation; LES=false, SI_units=true,
                                           initial_wall_time_seconds=1e-9*time_ns())
    model = simulation.model
    Δt = simulation.Δt

    i, t = model.clock.iteration, model.clock.time

    progress = 100 * (t / simulation.stop_time)

    current_wall_time = 1e-9 * time_ns() - initial_wall_time_seconds

    u_max = maximum(abs, model.velocities.u)
    v_max = maximum(abs, model.velocities.v)
    w_max = maximum(abs, model.velocities.w)

    if SI_units
        @info @sprintf("[%06.2f%%] i: % 6d,     time: % 10s,     Δt: % 10s,     wall time: % 8s",
                        progress, i, prettytime(t), prettytime(get_Δt(Δt)), prettytime(current_wall_time),)
        if LES
            ν_max = maximum(abs, model.diffusivity_fields.νₑ)
            @info @sprintf("          └── max(|u⃗|): [%.2e, %.2e, %.2e] m/s,     adv CFL: %.2e,     diff CFL: %.2e,     νₘₐₓ: %.2e m²/s",
                           u_max, v_max, w_max, AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model), ν_max)
        else
            @info @sprintf("          └── max(|u⃗|): [%.2e, %.2e, %.2e] m/s,     adv CFL: %.2e,     diff CFL: %.2e",
                           u_max, v_max, w_max, AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model))
        end

    else
        @info @sprintf("[%06.2f%%] i: % 6d,     time: %13.2f,     Δt: %13.2f,     wall time: % 8s",
                       progress, i, t, get_Δt(Δt), prettytime(current_wall_time))
        if LES
            ν_max = maximum(abs, model.diffusivity_fields.νₑ)
            @info @sprintf("          └── max(|u⃗|): [%.2e, %.2e, %.2e],     adv CFL: %.2e,     diff CFL: %.2e,     νₘₐₓ: %.2e",
                           u_max, v_max, w_max, AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model), ν_max)
        else
            @info @sprintf("          └── max(|u⃗|): [%.2e, %.2e, %.2e],     adv CFL: %.2e,     diff CFL: %.2e",
                           u_max, v_max, w_max, AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model))
        end
    end

    @info ""

    return nothing
end

function SimpleProgressMessenger(; kwargs...)
    return simulation -> SimpleProgressMessenger_function(simulation; kwargs...)
end

function SingleLineProgressMessenger_func(simulation; LES=false, SI_units=true,
                                          initial_wall_time_seconds=1e-9*time_ns())
    model = simulation.model
    Δt = simulation.Δt

    i, t = model.clock.iteration, model.clock.time

    progress = 100 * (t / simulation.stop_time)

    current_wall_time = 1e-9 * time_ns() - initial_wall_time_seconds

    if SI_units
        if LES
            ν_max = maximum(abs, model.diffusivity_fields.νₑ)
            @info @sprintf("[%06.2f%%] i: % 6d,     time: %10s,     Δt: %10s,     wall time: %8s,     adv CFL: %.2e,     diff CFL: %.2e,     νₘₐₓ: %.2e m²/s",
                           progress, i, prettytime(t), prettytime(get_Δt(Δt)), prettytime(current_wall_time),
                           AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model), ν_max)
        else
            @info @sprintf("[%06.2f%%] i: % 6d,     time: %10s,     Δt: %10s,     wall time: %8s,     adv CFL: %.2e,     diff CFL: %.2e",
                           progress, i, prettytime(t), prettytime(get_Δt(Δt)), prettytime(current_wall_time),
                           AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model))
        end

    else
        if LES
            ν_max = maximum(abs, model.diffusivity_fields.νₑ)
            @info @sprintf("[%06.2f%%] i: % 6d,     time: %13.2f,     Δt: %13.2f,     wall time: % 8s,     adv CFL: %.2e,     diff CFL: %.2e,     νₘₐₓ: %.2e",
                           progress, i, t, get_Δt(Δt), prettytime(current_wall_time), AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model), ν_max)
        else
            @info @sprintf("[%06.2f%%] i: % 6d,     time: %13.2f,     Δt: %13.2f,     wall time: % 8s,     adv CFL: %.2e,     diff CFL: %.2e",
                           progress, i, t, get_Δt(Δt), prettytime(current_wall_time), AdvectiveCFL(Δt)(model), DiffusiveCFL(Δt)(model))
        end
    end

    return nothing
end

function SingleLineProgressMessenger(; kwargs...)
    return simulation -> SingleLineProgressMessenger_func(simulation; kwargs...)
end

mutable struct TimedProgressMessenger{T, A, D, L, Δ} <: Function
    wall_time₀ :: T  # Wall time at simulation start
    wall_time⁻ :: T  # Wall time at previous calback
       adv_cfl :: A
       dif_cfl :: D
           LES :: L
            Δt :: Δ
end

TimedProgressMessenger(Δt) =
    TimedProgressMessenger(
                      1e-9 * time_ns(),
                      1e-9 * time_ns(),
                      AdvectiveCFL(Δt),
                      DiffusiveCFL(Δt),
                      false,
                      Δt)

function TimedProgressMessenger(; Δt=0, LES=false,
                       wall_time₀=1e-9*time_ns(), wall_time⁻=1e-9*time_ns(),
                       adv_cfl=nothing, dif_cfl=nothing)

    adv_cfl = adv_cfl==nothing ? AdvectiveCFL(Δt) : adv_cfl
    dif_cfl = dif_cfl==nothing ? DiffusiveCFL(Δt) : dif_cfl
    return TimedProgressMessenger(wall_time₀, wall_time⁻, adv_cfl, dif_cfl, LES, Δt)
end

function (pm::TimedProgressMessenger)(simulation)
    model = simulation.model
    Δt = simulation.Δt

    i, t = model.clock.iteration, model.clock.time

    progress = 100 * (t / simulation.stop_time)

    current_wall_time = 1e-9 * time_ns() - pm.wall_time₀
    time_since_last_callback = 1e-9 * time_ns() - pm.wall_time⁻
    wall_time_per_step = time_since_last_callback / simulation.iteration_interval
    pm.wall_time⁻ = 1e-9 * time_ns()

    u_max = maximum(abs, model.velocities.u)
    v_max = maximum(abs, model.velocities.v)
    w_max = maximum(abs, model.velocities.w)

    @info @sprintf("[%06.2f%%] iteration: % 6d, time: % 10s, Δt: % 10s, wall time: % 8s (% 8s / time step)",
                    progress, i, prettytime(t), prettytime(get_Δt(Δt)), prettytime(current_wall_time), prettytime(wall_time_per_step))

    if pm.LES
        ν_max = maximum(abs, model.diffusivity_fields.νₑ)
        @info @sprintf("          └── max(|u⃗|): [%.2e, %.2e, %.2e] m/s, CFL: %.2e, νCFL: %.2e, ν_max: %.2e m²/s",
                        u_max, v_max, w_max, pm.adv_cfl(model), pm.dif_cfl(model), ν_max)
    else
        @info @sprintf("          └── max(|u⃗|): [%.2e, %.2e, %.2e] m/s, CFL: %.2e, νCFL: %.2e",
                        u_max, v_max, w_max, pm.adv_cfl(model), pm.dif_cfl(model))
    end

    @info ""

    return nothing
end


text:
# This file contains tests that have become out-of-date with the
# current code.  They no longer run. It may be useful to rewrite them
# and re-include them in the test suite at some point, or perhaps just
# write completely new tests testing the same thing.

import Celeste: detect_sources
import Celeste.ParallelRun: BoundingBox, OptimizedSource,
    one_node_joint_infer, one_node_single_infer, detect_sources
import Celeste.SensitiveFloats: SensitiveFloat
import Celeste.DeterministicVI.ElboMaximize
import Celeste.Coordinates: match_coordinates

"""
load_ea_from_source
Helper function to load elbo args for a particular source
"""
function load_ea_from_source(target_source, target_sources, catalog,
                             images, all_vps)

    patches = Model.get_sky_patches(images, catalog)

    # Get neighbors of the source from which to load elbo args.
    neighbor_ids = Model.find_neighbors(patches, target_source)

    # Create a dictionary to the optimized parameters.
    target_source_variational_params = Dict{Int64, Array{Float64}}()
    for (indx, cur_source) in enumerate(target_sources)
        target_source_variational_params[cur_source] = all_vps[indx]
    end

    patches = Model.get_sky_patches(images, catalog)

    # Get neighbors of the source from which to load elbo args.
    neighbor_ids = Model.find_neighbors(patches, target_source)

    # Load neighbors, patches and variational parameters
    entry = catalog[target_source]
    neighbors = catalog[neighbor_ids]
    ids_local = vcat([target_source], neighbor_ids)
    cat_local = vcat([entry], neighbors)
    patches = patches[ids_local, :]  # limit patches to catalog
    vp = [haskey(target_source_variational_params, x) ?
          target_source_variational_params[x] :
          DeterministicVI.catalog_init_source(catalog[x]) for x in ids_local]

    # Create elbo args
    ElboArgs(images, patches, [1])
end

"""
compute_unconstrained_gradient
"""
function compute_unconstrained_gradient(target_source, target_sources,
                                        catalog, images, all_vps)
    # Load ea
    ea = load_ea_from_source(target_source, target_sources, catalog, images,
                             all_vps)

    # Evaluate in constrained space and then unconstrain
    # (taken from the old maximize_elbo.jl code)
    last_sf::SensitiveFloat{Float64} = SensitiveFloats.SensitiveFloat{Float64}(length(UnconstrainedParams), 1, true, true)
    transform = DeterministicVI.get_mp_transform(vp, ea.active_sources)
    f_res = DeterministicVI.elbo(ea)
    Transform.transform_sensitive_float!(transform, last_sf, f_res, vp, ea.active_sources)
    last_sf.d
end

"""
unconstrained_gradient_near_zero
"""
function unconstrained_gradient_near_zero(target_sources, catalog, images, all_vps)

    # Compute gradient per source
    for (cur_source_indx, source) in enumerate(target_sources)
        gradient = compute_unconstrained_gradient(source, target_sources, catalog, images,all_vps)
        if !isapprox(gradient, zeros(length(gradient)), atol=1)
            return false
        end
    end
    return true
end

"""
compare_vp_params
Helper to check whether two sets of variational parameters from the
result of joint inference are the same.

Return true if vp params are the same, false otherwise
"""
function compare_vp_params(flux_loc, flux_scale)

    length(flux_loc) == length(flux_scale) || return false

    # Check the existence and equivalence of each source's vp in flux_scale
    for i in eachindex(flux_loc)
        a = flux_loc[i].vs
        b = flux_scale[i].vs
        if !(isapprox(a, b))
            println("compare_vp_params: Mismatch - $(a) vs $(b)")
            print("norm(a - b): ", norm(a - b))
            return false
        end
    end
    return true
end

"""
compute_obj_value
computes obj value given set of results from one_node_infer and one_node_joint_infer
"""
function compute_obj_value(images::Vector{<:Image},
                           catalog::Vector{CatalogEntry},
                           box::BoundingBox,
                           results::Vector{OptimizedSource})

    # TODO: This stuff is duplicated from ParallelRun.infer_box.
    # We should refactor infer_box to return the objective value in some way!
    patches = Model.get_sky_patches(images, catalog)
    entry_in_range = entry->((box.ramin < entry.pos[1] < box.ramax) &&
                             (box.decmin < entry.pos[2] < box.decmax))
    target_ids = find(entry_in_range, catalog)

    # There must be a vp for every patch in the call to elbo().
    # So, here we must limit patches to just the targets we optimized
    # and pass [1, 2, 3, ...] as the target indexes.
    ea = ElboArgs(images, patches[target_ids, :],
                  collect(1:length(target_ids)); include_kl=false)
    vp = [r.vs for r in results]
    DeterministicVI.elbo(ea, vp).v[]
end

"""
load_stripe_82_data
"""
function load_stripe_82_data()
    rcf = RunCamcolField(4263, 5, 119)
    cd(datadir)
    run(`make RUN=$(rcf.run) CAMCOL=$(rcf.camcol) FIELD=$(rcf.field)`)
    cd(wd)
    catalog, target_sources = infer_init([rcf], datadir, primary_initialization=false)
    images = SDSSIO.load_field_images([rcf], datadir)
    ([rcf], datadir, target_sources, catalog, images)
end

# TODO - This doesn't pass right now since some light sources are
# not close to zero.
function test_gradient_is_near_zero_on_stripe_82()
    (rcfs, datadir, target_sources, catalog, images) = load_stripe_82_data()

    # Make sure joint infer with few iterations does not pass the gradient near zero check
    joint_few(cnti...) = one_node_joint_infer(cnti...)
    results_few = one_node_infer(rcfs, datadir; infer_callback=joint_few, primary_initialization=false)
    @test !unconstrained_gradient_near_zero(target_sources, catalog, images, [x.vs for x in results_few])

    # Make sure joint infer with many iterations passes the gradient near zero check
    joint_many(cnti...) = one_node_joint_infer(cnti...; n_iters=10)
    results_many = one_node_infer(rcfs, datadir; infer_callback=joint_many, primary_initialization=false)
    @test unconstrained_gradient_near_zero(target_sources, catalog, images, [x.vs for x in results_many])
end

# Tests that the gradient is zero on overlapping sources after
# joint optimization.
# Makes sure that with fewer iterations the gradient is not zero.
function test_gradient_is_near_zero()
    images = SampleData.get_sdss_images(4263, 5, 119)
    catalog = SampleData.get_sdss_catalog(4263, 5, 119)

    # This box has 3 overlapping objects in it.
    box = BoundingBox(0.467582, 0.473275, 0.588383, 0.595095)

    result_few = ParallelRun.infer_box(images, catalog, box;
                                       method=:joint, n_iters=2)

    @test length(result_few) == 3

    # this stuff is copied from infer_box so that we can calculate the gradient
    #catalog, patches = detect_sources(images)
    patches = Model.get_sky_patches(images, catalog)
    entry_in_range = entry->((box.ramin < entry.pos[1] < box.ramax) &&
                             (box.decmin < entry.pos[2] < box.decmax))
    target_sources = find(entry_in_range, catalog)

    @test !unconstrained_gradient_near_zero(target_sources, catalog, images,
                                            [x.vs for x in result_few])

    result_many = ParallelRun.infer_box(images, catalog, box;
                                        method=:joint, n_iters=20)
    @test length(result_many) == 3
    @test unconstrained_gradient_near_zero(target_sources, catalog, images,
                                           [x.vs for x in result_many])
end

"""
test_different_result_with_different_iter()
Using 3 iters instead of 1 iters should result in a different set of parameters.
"""
function test_different_result_with_different_iter()
    # This bounding box has overlapping stars. (neighbor map is not empty)
    box = BoundingBox(154.39, 164.41, 39.11, 39.13)
    field_triplets = [RunCamcolField(3900, 6, 269),]

    infer_few(ctni...) = one_node_joint_infer(ctni...;
                                              n_iters=1,
                                              within_batch_shuffling=false)
    result_iter_1 = one_node_infer(field_triplets, datadir;
                                   infer_callback=infer_few, box=box)

    infer_many(ctni...) = one_node_joint_infer(ctni...;
                                               n_iters=5,
                                               within_batch_shuffling=false)
    result_iter_5 = one_node_infer(field_triplets, datadir;
                                   infer_callback=infer_many, box=box)

    # Make sure that parameters are not the same
    @test !compare_vp_params(result_iter_1, result_iter_5)
end

"""
test_same_result_with_diff_batch_sizes
Varying batch sizes using the cyclades algorithm should not change final objective value.
"""
function test_same_result_with_diff_batch_sizes()
    # This bounding box has overlapping stars. (neighbor map is not empty)
    box = BoundingBox(154.39, 164.41, 39.11, 39.13)
    field_triplets = [RunCamcolField(3900, 6, 269),]

    # With batch size = 7
    infer_few(ctni...) = one_node_joint_infer(ctni...;
                                n_iters=3,
                                batch_size=7,
                                within_batch_shuffling=false)
    result_bs_7 = one_node_infer(field_triplets, datadir;
                                 infer_callback=infer_few,
                                 box=box)

    # With batch size = 39
    infer_many(ctni...) = one_node_joint_infer(ctni...;
                                n_iters=3,
                                batch_size=39,
                                within_batch_shuffling=false)
    result_bs_39 = one_node_infer(field_triplets, datadir;
                                  infer_callback=infer_many,
                                  box=box)

    # Make sure that parameters are exactly the same
    @test compare_vp_params(result_bs_7, result_bs_39)
end

function test_same_one_node_infer_twice()
    # This bounding box has overlapping stars. (neighbor map is not empty)
    box = BoundingBox(154.39, 164.41, 39.11, 39.13)
    field_triplets = [RunCamcolField(3900, 6, 269),]

    infer_few(ctni...) = one_node_joint_infer(ctni...;
                                              n_iters=3,
                                              batch_size=7,
                                              within_batch_shuffling=false)
    result_bs_7_1 = one_node_infer(field_triplets, datadir;
                                   infer_callback=infer_few,
                                   box=box)

    result_bs_7_2 = one_node_infer(field_triplets, datadir;
                                   infer_callback=infer_few,
                                   box=box)

    # Make sure that parameters are exactly the same
    @test compare_vp_params(result_bs_7_1, result_bs_7_2)
end

"""
test infer multi iter obj overlapping.
This makes sure one_node_joint_infer achieves sum objective value less
than single iter on non overlapping sources.
"""
function test_one_node_joint_infer_obj_overlapping()
    # This bounding box has overlapping stars. (neighbor map is not empty)
    box = BoundingBox(154.39, 164.41, 39.11, 39.13)
    field_triplets = [RunCamcolField(3900, 6, 269),]

    # 100 iterations
    tic()
    infer_multi(ctni...) = one_node_joint_infer(ctni...;
                                                n_iters=30,
                                                within_batch_shuffling=true)
    result_multi = one_node_infer(field_triplets, datadir;
                                  infer_callback=infer_multi,
                                  box=box)
    multi_iter_time = toq()
    score_multi = compute_obj_value(result_multi, field_triplets, datadir; box=box)

    # 2 iterations
    tic()
    infer_two(ctni...) = one_node_joint_infer(ctni...;
                                              n_iters=2,
                                              within_batch_shuffling=true)
    result_two = one_node_infer(field_triplets, datadir;
                                infer_callback=infer_two,
                                box=box)
    two_iter_time = toq()
    score_two = compute_obj_value(result_two, field_triplets, datadir; box=box)

    # One node infer (1 iteration, butm ore newton steps)
    tic()
    result_single = one_node_infer(field_triplets, datadir;
                                   infer_callback=one_node_single_infer,
                                   box=box)
    single_iter_time = toq()
    score_single = compute_obj_value(result_single, field_triplets, datadir; box=box)

    println("One node joint infer objective value: $(score_multi)")
    println("One node joint infer 2 iter objective value: $(score_two)")
    println("Single iter objective value: $(score_single)")
    println("One node joint infer time: $(multi_iter_time)")
    println("One node joint infer 2 iter time: $(two_iter_time)")
    println("Single iter time: $(single_iter_time)")
    @test score_multi > score_single
    @test score_multi > score_two
end

# Stripe 82 tests are long running
@testset "stripe82" begin
    # Test gradients near zero (this takes 100 iterations and is a bit slow)
    test_gradient_is_near_zero()

    # Test that we reach a higher objective with more iterations.
    # This is a bit slow.
    test_one_node_joint_infer_obj_overlapping()

    test_same_one_node_infer_twice()
    test_different_result_with_different_iter()
    test_same_result_with_diff_batch_sizes()
end


text:
# NOTE: This file contains machine-generated code by translate_exported_functions.jl.

# void hts_set_log_level(enum htsLogLevel level);
@defun hts_set_log_level(level::htsLogLevel)::Cvoid

# enum htsLogLevel hts_get_log_level(void);
@defun hts_get_log_level()::htsLogLevel

# void hts_log(enum htsLogLevel severity, const char *context, const char *format, ...) HTS_FORMAT(HTS_PRINTF_FMT, 3, 4);
# ErrorException("unknown token: '.'")


text:
using LinearAlgebra, LuxurySparse

LinearAlgebra.tr(x::Value) = register(LinearAlgebra.tr, x)
gradient(::typeof(tr), grad, output, x) = (grad * IMatrix{size(x, 1), eltype(grad)}(), )

gradient(::typeof(transpose), grad, output, x::AbstractMatrix) = (transpose(grad), )
function gradient(::typeof(*), grad, output, lhs::AbstractArray{<:Real}, rhs::AbstractArray{<:Real})
    grad * transpose(rhs), transpose(lhs) * grad
end


text:
# Autogenerated wrapper script for rdma_core_jll for aarch64-linux-musl
export lbrdmacm, libibverbs

using libnl_jll
JLLWrappers.@generate_wrapper_header("rdma_core")
JLLWrappers.@declare_library_product(lbrdmacm, "librdmacm.so.1")
JLLWrappers.@declare_library_product(libibverbs, "libibverbs.so.1")
function __init__()
    JLLWrappers.@generate_init_header(libnl_jll)
    JLLWrappers.@init_library_product(
        lbrdmacm,
        "lib/librdmacm.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libibverbs,
        "lib/libibverbs.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test
using Forecast
using TimeSeries

@testset "ccf" begin
    x1 = rand(100);
    x2 = circshift(x1,6);
    @test_throws AssertionError ccf(x1,x2; lag = 100)

    res = ccf(x1, x2; lag = 20, type="cor");
    @test res isa CCF
    @test res.call == "ccf(x1, x2; type=\"cor\", lag=20, alpha=(0.95, 0.99))"
    @test res.type == "cor"
    @test !res.auto
    @test res.lag == 20
    @test res.N == 100
    @test res.alpha == (0.95, 0.99)
    @test res.ci[1] ≈ 0.19900419155027554
    @test res.ci[2] ≈ 0.2615358405397182

    @test res.ccf[20+1+6] ≈ 0.9685785099293228
end


text:
## Metadata For graph creation
@enum VertexLocation LOC_PMEM LOC_DRAM LOC_DRAM_PRE LOC_SOURCE LOC_SINK

# Helpers - generally pretty simple but helpful for abstracting the `isdram` attribute.
ispmem(loc::VertexLocation) = loc == LOC_PMEM
isdram(loc::VertexLocation) = loc == LOC_DRAM || loc == LOC_DRAM_PRE
issource(loc::VertexLocation) = loc == LOC_SOURCE
issink(loc::VertexLocation) = loc == LOC_SINK

@enum EdgeType begin
    EDGE_NONE
    EDGE_SYNC_READ
    EDGE_SYNC_WRITE
    EDGE_ASYNC_READ
    EDGE_ASYNC_WRITE
end

isasync(et::EdgeType) = in(et, (EDGE_ASYNC_READ, EDGE_ASYNC_WRITE))

@enum MoveType MOVE_NONE MOVE_SYNC MOVE_ASYNC

# Metadata to assign to each node in the liveness graph for tensors.
"""
    VertexMetadata

Metadata associated with the vertex of the tensor graph.

$(FIELDS)
"""
struct VertexMetadata
    "The gadget that this vertex belongs to. Used for edge generation."
    gadget::Int
    "The op index that this gadget refers to."
    op::XNode
    "Where the vertex lives."
    location::VertexLocation
    "What type of moves this vertex allows."
    move_type::MoveType
    "`true` is this vertex is a user of the tensor in whose graph this vertex appear."
    isuser::Bool
    vertex_number::Int
end

struct EdgeMetadata
    edgetype::EdgeType
end
isasync(em::EdgeMetadata) = isasync(em.edgetype)

#####
##### Preprocessing
#####

"""
$(SIGNATURES) -> UnitRange{Int64}

Return a range of indices corresponding to nodes during which `xtensor` is live.
"""
liverange(xtensor::XTensor) = (producer(xtensor).index):(consumer(xtensor).index)

function _getgadgets(A::ILPHolder{Asynchronous}, data::FunctionData, t::XTensor)
    range = liverange(t)
    livenodes = (nodes(data)[x] for x in range)
    refs = Vector{NamedTuple{(:node, :move_type),Tuple{XNode,MoveType}}}()

    # Build the referece map
    reference_map = Dict{XNode, XNode}()
    ref = producer(t)

    # To decide if a node should be considered as an aynchronous move point, we check to see
    # if the node is with `bound` distance of a user of the tensor.
    #
    # The intuition here is that moves will probably be located closer to their producers or
    # consumers rather than further.
    #
    # Making `bound` larger increased the search space of the formulation, which may lead to
    # better results at the cost of a larger mode.
    bound = 15
    move_time = sizeof(t) / A.write_bandwidth
    for ind in range
        node = nodes(data)[ind]
        if in(node, users(t))
            push!(refs, (node = node, move_type = MOVE_SYNC))
            ref = node

        # Check if there is a user node within `bound`. If so, make this an async move node.
        elseif hasprofile(node) &&
            any(
                in(users(t)),
                (nodes(data)[i] for i in max(ind-bound, 1):min(ind+bound, length(nodes(data))))
            )

            push!(refs, (node = node, move_type = MOVE_ASYNC))
            ref = node
        end
        reference_map[node] = ref
    end

    return refs, reference_map
end

function _getgadgets(::ILPHolder{Synchronous}, data::FunctionData, t::XTensor)
    range = liverange(t)
    livenodes = (nodes(data)[x] for x in range)

    # Build the referece map
    reference_map = Dict{XNode, XNode}()
    ref = producer(t)
    for ind in range
        node = nodes(data)[ind]
        if in(node, users(t))
            ref = node
        end
        reference_map[node] = ref
    end

    nt = [
        (node = u, move_type = isone(i) ? MOVE_NONE : MOVE_SYNC) for (i,u) in enumerate(users(t))
    ]

    return nt, reference_map
end

function _getgadgets(::ILPHolder{Static}, data::FunctionData, t::XTensor)
    range = liverange(t)
    producer = nodes(data)[first(range)]

    reference_map = Dict{XNode, XNode}()
    for ind in range
        reference_map[nodes(data)[ind]] = producer
    end

    return [(node = producer, move_type = MOVE_NONE)], reference_map
end

function getgadgets(x, data::FunctionData, t::XTensor)
    refs, reference_map = _getgadgets(x, data, t)
    # Update with argument

    # If this tensor `isarg`, then since it is live for the entire function, we need to
    # insert node references for all nodes pointing to the producer.
    if isarg(t)
        ref = producer(t)
        for node in nodes(data)
            # Shortcut here, `get!` will lookup and automatically insert the mapping if it
            # doesn't exist.
            get!(reference_map, node, ref)
        end
    end

    return refs, reference_map
end

# TODO: When can't do anything because exhausted, make ASCII art for what's going on with
# these tensor graphs.
#=
=#

function edge_metadata(src, dst, s, d, src_move_type, tensor::XTensor)
    # Setup the correct move annotations.
    if src_move_type == MOVE_ASYNC
        edge_read_type = EDGE_ASYNC_READ
        edge_write_type = EDGE_ASYNC_WRITE
    else
        edge_read_type = EDGE_SYNC_READ
        edge_write_type = EDGE_SYNC_WRITE
    end

    # Determine if an edge should be added and what kind of edge it is.

    ### LOC_SOURCE node is source
    if (src, dst) == (LOC_SOURCE, LOC_DRAM_PRE)
        isone(d) && return EdgeMetadata(EDGE_NONE)
    elseif (src, dst) == (LOC_SOURCE, LOC_PMEM)
        isone(d) && return EdgeMetadata(EDGE_NONE)
    # If the tensor is an argument, we need to create an edge directly from source to sink.
    # This indicates that the tensor begins in DRAM - so will always be in DRAM
    elseif (src, dst) == (LOC_SOURCE, LOC_SINK)
        isarg(tensor) && return EdgeMetadata(EDGE_NONE)

    ### LOC_DRAM as source

    # If we're in LOC_DRAM, data already exists in PMEM. Thus all edges originating in
    # LOC_DRAM have no metadata. I.E., no read or write type eges.
    elseif (src, dst) == (LOC_DRAM, LOC_DRAM)
        s == d-1 && return EdgeMetadata(EDGE_NONE)
    elseif (src, dst) == (LOC_DRAM, LOC_PMEM)
        s == d-1 && return EdgeMetadata(EDGE_NONE)
    elseif (src, dst) == (LOC_DRAM, LOC_SINK)
        s == d-1 && return EdgeMetadata(EDGE_NONE)

    ### LOC_DRAM_PRE as source

    # Moving from LOC_DRAM_PRE to LOC_PMEM indicates a movement of data from DRAM to PMEM
    # where the data did not exist in PMEM before. Thus, we must record a write.
    elseif (src, dst) == (LOC_DRAM_PRE, LOC_PMEM)
        s == d-1 && return EdgeMetadata(edge_write_type)
    # Movement between DRAM or to the SINK requires no metadata
    elseif (src, dst) == (LOC_DRAM_PRE, LOC_DRAM_PRE)
        s == d-1 && return EdgeMetadata(EDGE_NONE)
    elseif (src, dst) == (LOC_DRAM_PRE, LOC_SINK)
        s == d-1 && return EdgeMetadata(EDGE_NONE)

    ### LOC_PMEM as source

    # Movement edges from PMEM to DRAM. Must not happen on the first `component` because
    # prefetching from PMEM to DRAM at the time of creation makes no sense.
    elseif (src, dst) == (LOC_PMEM, LOC_DRAM)
        (s == d) && !isone(s) && return EdgeMetadata(edge_read_type)
    elseif (src, dst) == (LOC_PMEM, LOC_PMEM)
        (s == d-1) && return EdgeMetadata(EDGE_NONE)
    # Otherwise, no metadata needed
    elseif (src, dst) == (LOC_PMEM, LOC_SINK)
        (s == d-1) && return EdgeMetadata(EDGE_NONE)
    end
    return nothing
end

function preprocess!(S::ILPHolder, data::FunctionData)
     for tensor in tensors(data)

        # Get the users of this node
        # Get two things from getgadgets:
        #
        # 1. A named tuple (node::XNode, move_type::MoveType)
        # 2. A dictionary implementing the `ref` function.
        gadgets, reference_map = getgadgets(S, data, tensor)

        @assert !isempty(gadgets)

        # Graph building time :D
        g = MetaGraph(DiGraph(), EdgeMetadata, VertexMetadata)

        # If this tensor can only be assigned to a single location - don't generate a graph
        # for it - otherwise, populate the nodes in a tensor graph
        fixed_tensor = length(locations(tensor)) == 1
        if !fixed_tensor
            for (count, nt) in enumerate(gadgets)
                # Unpack the XNode
                node = nt.node
                move_type = nt.move_type
                islast = (count == length(gadgets))

                isuser = in(node, users(tensor))

                if count == 1
                    add_vertex!(g, VertexMetadata(0, node, LOC_SOURCE, move_type, isuser, nv(g)+1))
                end
                # Enumerate over locations that this tensor can live.
                #
                # Do it this way because some nodes can only live in DRAM, so iterating
                # then filtering takes care of that
                for location in locations(tensor)
                    if location == DRAM
                        # Add DRAM nodes
                        #
                        # Do not add these nodes if this tensor is an argument because it will
                        # have already originated in PMM
                        if !isarg(tensor)
                            add_vertex!(g,
                                VertexMetadata(count, node, LOC_DRAM_PRE, move_type, isuser, nv(g)+1)
                            )
                        end

                        # only add a DRAM node if there could have been a write to PMEM
                        if count > 1
                            add_vertex!(g,
                                VertexMetadata(count, node, LOC_DRAM, move_type, isuser, nv(g)+1)
                            )
                        end
                    end

                    if location == PMEM
                        @assert !startswith(nGraph.name(tensor), "Constant")
                        # PMEM node
                        add_vertex!(g,
                            VertexMetadata(count, node, LOC_PMEM, move_type, isuser, nv(g)+1)
                        )
                    end
                end
                if islast
                    # Set the gadget number for the sink to one higher than the last count.
                    add_vertex!(g,
                        VertexMetadata(count + 1, node, LOC_SINK, move_type, isuser, nv(g)+1)
                    )
                end
            end
        end

        # Use a quadratic complexity algorithm for doing edge assignment. It's not
        # perfect but it's simple, and as long as the graphs don't get too big should
        # run quickly enough for our purposes.
        for src in vertices(g), dst in vertices(g)
            src == dst && continue

            src_meta = getmeta(g, src)
            dst_meta = getmeta(g, dst)

            metadata = edge_metadata(
                src_meta.location,
                dst_meta.location,
                src_meta.gadget,
                dst_meta.gadget,
                src_meta.move_type,
                tensor,
            )

            isnothing(metadata) && continue

            add_edge!(g, src, dst, metadata)
        end

        # Create the descriptor
        S.descriptors[tensor] = TensorMeta(
            g,
            [g.node for g in gadgets],
            reference_map,
            fixed_tensor
        )
    end
end

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
module InverseIsing

export GBM,
       AbstractGBM,
       fit,
       coef,
       weights,
       infer,
       decode,
       dprint,
       anneal,
       trans

include("core.jl")

end # module

max_stars_repo_name: oilneck/InverseIsing

max_issues_repo_name: oilneck/InverseIsing

max_forks_repo_name: oilneck/InverseIsing


text:
include("test_common.jl")
using Random

new_context()

@testset "errors" begin
    @test_throws NodeError node"x"
    p = node"(async () => { throw('async error') })()"
    @test_throws NodeError wait(p)
    err = try
        node"throw('error')"
    catch e
        e
    end
    @test startswith(sprint(show, err), "NodeError(")
    @test startswith(sprint(showerror, err), "NodeError:")
    @test_throws NodeError begin
        node_throw("error")
        NodeCall.throw_error()
    end
    @test_throws NodeError node"f=>f()"() do
        error("error")
    end
end

@testset "references" begin
    x = rand()
    NodeCall.reference(x)
    ref = NodeCall.reference(x)
    ref2 = NodeCall.reference(x)
    @test pointer(ref) == pointer(ref2)
    NodeCall.dereference(x)
    NodeCall.dereference(pointer(ref))
    @test NodeCall.delete_reference(pointer(ref2))
    @test !NodeCall.delete_reference(C_NULL)
end

@testset "internals" begin
    @test JsValue(2) == 2
    @test JsValue(Bool, node"true"o)
    @test value(value(NodeValueTemp, 1)) == 1
    struct S end
    o = node"{}"o
    @test_throws ErrorException convert(S, o)
    @test NodeCall.is_array(["a", "b"])
    @test NodeCall.get_type(123) == NodeCall.NapiTypes.napi_number
    o2 = node"1"o
    @test startswith(sprint(show, o), "NodeObject:")
    @test startswith(sprint(show, o2), "NodeValueTemp:")
    o = nothing
    o2 = nothing
    GC.gc()
    NodeCall.dispose()
end


text:
@doc """
Amplitude and Phase Estimation of a Sinusoid (APES).
Implemented as described in [Li, Li & Stoica 1998](http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=700967).

This function estimates the complex amplitude of signal `y` at frequencies `f`.
M denotes the filter length, which is set to N/2 by default [[Li & Stoica 1996](http://ieeexplore.ieee.org/xpls/abs_all.jsp?arnumber=506612)].

## References
H. Li, J. Li, and P. Stoica. Performance analysis of forward-backward matched-filterbank spectral esti- mators. Signal Processing, IEEE Transactions on, 46(7):1954–1966, 1998.

J. Li and P. Stoica. An adaptive filtering approach to spectral estimation and sar imaging. Signal Processing, IEEE Transactions on, 44(6):1469–1484, 1996.

Y. Hua, A. Gershman, and Q. Cheng. High-resolution and robust signal processing. Signal Processing and Communications. CRC Press, October 2003.

""" ->
function apes{T <: FloatingPoint}(y::Array{T}, f::Union(AbstractArray{T}, T); M::Int=div(length(y), 2), Q::Union(Array{T, 2}, Nothing)=nothing)

    N = length(y)
    L = N - M + 1  # (eqn 3)

    y_f = zeros(M, L)  # Forward data vectors
    y_b = zeros(M, L)  # Backward data vectors
    R_f = zeros(M, M)  # Forward data covariance matrix
    R_b = zeros(M, M)  # Backward data covariance matrix
    for l = 0:L-1
        y_f[:, l+1] = y[l+1 : 1 : l+M-1+1]               # eqn 3
        y_b[:, l+1] = conj(y[N-l-1+1 : -1 : N-l-M+1])    # eqn 6
        R_f += y_f[:, l+1] * y_f[:, l+1]'                # eqn 12
        R_b += y_b[:, l+1] * y_b[:, l+1]'                # eqn 13
    end
    R_f = R_f ./ L
    R_b = R_b ./ L

    # Forward-backward estimate of covariance matrix (eqn 14)
    R  = (1/2) * (R_f + R_b)

    # Calculate complex amplitude for each frequency
    α = zeros(Complex{T}, size(f))
    for i in 1:length(f)

        ω = 2 * pi * f[i]
        α[i] = apes(ω, M, L, y_f, y_b, R, Q)
    end

    return α
end

function apes{T <: FloatingPoint}(ω::T, M::Int, L::Int, y_f::Array{T, 2}, y_b::Array{T, 2}, R::Array{T, 2}, Q::Union(Nothing, Array{T, 2}))

    a_m = exp(-im * ω .* [0:1:M-1 ; ])
    a_l = exp(-im * ω .* [0:1:L-1 ; ])

    # Normalised Fourier transforms of forward and backward data vectors
    g_f = (1/L) * y_f * conj(a_l)  # eqn 4.3.11 [1]
    g_b = (1/L) * y_b * conj(a_l)
    G   = (1/sqrt(2)) * [g_f g_b]  # eqn 32

    # Estimate of noise covariance matrix
    if Q == nothing
        Q  = R - G * ctranspose(G)  # eqn 31
    end
    Qi = inv(Q)

    # Complex amplitude of sinusoid signal with frequency ω (eqn 34)
    (2 * (a_m' * Qi * g_f) / (a_m' * Qi * a_m))[1]
end


text:
include("utils.jl")

filename = "data/avg_earth_hs_he_12_hp_4_ve_12_vp_4_roefanov.jld2"

fig = Figure(resolution = (1700+600, 1000+400))
add_label = true
println("looking at ", filename)
state_names = []

i = 1
ii = (i-1)%3 + 1 # +1 on why 1 based indexing is wrong
jj = (i-1)÷3 + 1 # +1 on why 1 based indexing is wrong
jl_file = jldopen(filename, "r+")
s_string = "u"
slice_zonal = grab_state(s_string, jl_file)
colorrange, contour_levels, s_string = plot_helper(s_string, slice_zonal)
colorrange = (-40,40) # override
λ, ϕ, r, p_coord = grab_grid(jl_file)
push!(state_names, s_string)
ax1 = fig[jj,ii] = Axis(fig,title = state_names[i], titlesize = 40)
contour_heatmap!(ax1, ϕ, p_coord, slice_zonal,
                 contour_levels, colorrange,
                 add_labels = add_label, random_seed = 1)

i = 2
ii = (i-1)%3 + 1 # +1 on why 1 based indexing is wrong
jj = (i-1)÷3 + 1 # +1 on why 1 based indexing is wrong
jl_file = jldopen(filename, "r+")
s_string = "T"
slice_zonal = grab_state(s_string, jl_file)
colorrange, contour_levels, s_string = plot_helper(s_string, slice_zonal)
colorrange = (180,310)
λ, ϕ, r, p_coord = grab_grid(jl_file)
push!(state_names, s_string)
ax2 = fig[jj,ii] = Axis(fig,title = state_names[i], titlesize = 40)
contour_heatmap!(ax2, ϕ, p_coord, slice_zonal, contour_levels, colorrange, add_labels = add_label, colormap = :thermometer)
hideydecorations!(ax2, grid = false)

i = 3
ii = (i-1)%3 + 1 # +1 on why 1 based indexing is wrong
jj = (i-1)÷3 + 1 # +1 on why 1 based indexing is wrong
jl_file = jldopen(filename, "r+")

s_string = "TT"
slice_zonal, s_string = eddy_variance(s_string, jl_file)
colorrange, contour_levels, s_string = plot_helper(s_string, slice_zonal)
colorrange = (-8,48) # override
λ, ϕ, r, p_coord = grab_grid(jl_file)
push!(state_names, s_string)
ax3 = fig[jj,ii] = Axis(fig,title = state_names[i], titlesize = 40)
contour_heatmap!(ax3, ϕ, p_coord, slice_zonal, contour_levels,
                 colorrange, add_labels = add_label,
                 colormap = :thermometer, random_seed = 12)
hideydecorations!(ax3, grid = false)

i = 4
ii = (i-1)%3 + 1 # +1 on why 1 based indexing is wrong
jj = (i-1)÷3 + 1 # +1 on why 1 based indexing is wrong
jl_file = jldopen(filename, "r+")
s_string = "uv"
slice_zonal, s_string = eddy_variance(s_string, jl_file)
colorrange, contour_levels, s_string = plot_helper(s_string, slice_zonal)
colorrange = (-60, 60)
λ, ϕ, r, p_coord = grab_grid(jl_file)
push!(state_names, s_string)
ax4 = fig[jj,ii] = Axis(fig,title = state_names[i], titlesize = 40)
contour_heatmap!(ax4, ϕ, p_coord, slice_zonal, contour_levels, colorrange, add_labels = add_label)

i = 5
ii = (i-1)%3 + 1 # +1 on why 1 based indexing is wrong
jj = (i-1)÷3 + 1 # +1 on why 1 based indexing is wrong
jl_file = jldopen(filename, "r+")
s_string = "vT"
slice_zonal, s_string = eddy_variance(s_string, jl_file)
colorrange, contour_levels, s_string = plot_helper(s_string, slice_zonal)
colorrange = (-24, 24) # override
λ, ϕ, r, p_coord = grab_grid(jl_file)
push!(state_names, s_string)
ax5 = fig[jj,ii] = Axis(fig,title = state_names[i], titlesize = 40)
contour_heatmap!(ax5, ϕ, p_coord, slice_zonal, contour_levels, colorrange, add_labels = add_label)
hideydecorations!(ax5, grid = false)

i = 6
ii = (i-1)%3 + 1 # +1 on why 1 based indexing is wrong
jj = (i-1)÷3 + 1 # +1 on why 1 based indexing is wrong
jl_file = jldopen(filename, "r+")

s_string = "uu"
slice_zonal1, s_string = eddy_variance(s_string, jl_file)
colorrange, contour_levels = plot_helper(s_string, slice_zonal)

s_string = "vv"
slice_zonal2, s_string = eddy_variance(s_string, jl_file)

slice_zonal = 0.5 .* (slice_zonal1 + slice_zonal2)
colorrange = (0, 360) # modify to extrema

s_string = L"\langle (u' u' + v' v')/2 \rangle"
push!(state_names, s_string)
λ, ϕ, r, p_coord = grab_grid(jl_file)
ax6 = fig[jj,ii] = Axis(fig,title = state_names[i], titlesize = 40)
contour_heatmap!(ax6, ϕ, p_coord, slice_zonal,
                 contour_levels, colorrange, add_labels = add_label,
                 colormap = :thermometer, random_seed = 10)
hideydecorations!(ax6, grid = false)

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
# Copyright (c) 2020-2021, Bank of Canada
# All rights reserved.

@testset "workspace" begin
    # Create an empty workspace
    @test (Workspace(); true)
    # Add keys to the workspace
    @test (work1 = Workspace(); work1.a = 1; true)
    # Create a new workspace
    let work1 = Workspace()
        work1.a = 1
        work1.ts = TSeries(2020Q1,randn(10))
        work1.mvts = MVTSeries(2020Q1,(:a,:b),randn(10,2))
        # propertynames
        @test (isa(propertynames(work1),Tuple))
        # getproperty
        @test (work1.a == 1)
        # setproperty
        @test (work1.a = 2; work1.a == 2)
        # isempty
        @test (isempty(work1) == false)
        @test (isempty(Workspace()))
        # in
        @test (in(:b,work1) == false)
        @test (in(:ts,work1))
        # keys
        @test (collect(keys(work1)) == [:a, :ts, :mvts])
        # values
        @test (isa(values(work1),Base.ValueIterator))
        # iterate
        #
        # show
        let io = IOBuffer()
            @test (show(io, work1); true)
        end
    end
end


text:
module ConstrainedRootSolvers

using DocStringExtensions: METHODLIST, TYPEDEF, TYPEDFIELDS
using UnPack: @unpack

# export public Types
export BisectionMethod, NelderMeadMethod, NewtonBisectionMethod,
       NewtonRaphsonMethod, ReduceStepMethod, ReduceStepMethodND,
       ResidualTolerance, SolutionTolerance, SolutionToleranceND

# export public functions
export find_peak, find_zero

include("method.jl"   )
include("tolerance.jl")
include("find_zero.jl")

include("findpeak/bisection.jl"   )
include("findpeak/neldermead.jl"  )
include("findpeak/reducestep.jl"  )
include("findpeak/reducestepND.jl")

end # module


text:
# Test intersection detection and clustering

module TestIntersections

using OpenStreetMap
using Base.Test
using Compat

MAP_FILENAME = "tech_square.osm"

bounds = Bounds(42.3637, 42.3655, -71.0919, -71.0893)

bounds_ENU = ENU(bounds)

nodesLLA, hwys, builds, feats = getOSMData(MAP_FILENAME)
nodes = ENU(nodesLLA, center(bounds))

# Find intersections in map
intersections = findIntersections(hwys)
@test length(intersections) == 91

# Find Highway Sets
highway_sets = findHighwaySets(hwys)
@test length(highway_sets) == 4

# Cluster intersections
intersection_cluster_mapping = findIntersectionClusters(nodes,intersections,highway_sets,max_dist=15)
intersection_clusters = unique(collect(values(intersection_cluster_mapping)))
@test length(intersection_clusters) == 5

# Replace Nodes in Highways
replaceHighwayNodes!(hwys,intersection_cluster_mapping)
intersections_clustered = findIntersections(hwys)
@test length(intersections_clustered) == 82

# Check how many semi-redundant intersections we were able to remove
removed = length(intersections) - length(intersections_clustered)
@test removed == 9

end # module TestIntersections


text:
module OnlineEstimators

include("parameter_estimators.jl")
include("state_estimators.jl")

end


text:
# This file was generated, do not modify it. # hide
function print_performance(model, data...)
    e = evaluate(model, data...;
                 resampling=CV(rng=StableRNG(1234), nfolds=8),
                 measure=rms,
                 verbosity=0)
    μ = round(e.measurement[1], sigdigits=5)
    ste = round(std(e.per_fold[1])/sqrt(8), digits=5)
    println("$model = $μ ± $(2*ste)")
end;

X, y = @load_boston

print_performance(linear, X, y)
print_performance(knn, X, y)
print_performance(average_two, X, y)


text:
# Alignment functions to create parallel data

using .DTWs

import MelGeneralizedCepstrums: mc2e

# general feature alignment
function align(src::AbstractMatrix, # source feature matrix
               tgt::AbstractMatrix  # target feature matrix
               )
    if size(src, 1) != size(tgt, 1)
        throw(DimentionMismatch("order of feature vector must be equal"))
    end

    # Alignment
    d = DTW(fstep=0, bstep=2) # allow one skip
    path = fit!(d, src, tgt)

    # create aligned tgt
    newtgt = zeros(eltype(src), size(src))
    newtgt[:,path] = tgt[:,1:length(path)]

    # interpolation
    # TODO(ryuichi) better solution
    hole = setdiff(collect(path[1]:path[end]), path)
    for i in hole
        if i > 1 && i < size(src, 2)
            for j=1:size(newtgt, 1)
                @inbounds newtgt[j,i] = (newtgt[j,i-1] + newtgt[j,i+1]) / 2.0
            end
        end
    end

    src, newtgt
end

# alignment function specialized for mel-cepstrum
function align_mcep(src::AbstractMatrix,       # source feature matrix
                    tgt::AbstractMatrix,       # target feature matrix
                    α::AbstractFloat,         # all-pass constant
                    fftlen::Integer;           # fft length
                    threshold::Float64=-14.0,  # threshold to remove silence frames
                    remove_silence::Bool=true
                    )
    src, newtgt = align(src, tgt)

    if remove_silence
        e = mc2e(src, α, fftlen)
        E = log(e)
        src = src[:, E .> threshold]
        newtgt = newtgt[:, E .> threshold]
    end

    src, newtgt
end


text:
using KDTrees, TaxiSimulation, JLD, Base.Dates, DataFrames, LightGraphs

#Ride selector
function good_ride(df::DataFrame, i::Int)
  # testing time of day
  if dayofweek(df[i, :pTime]) <= 5 && 12 <= hour(df[i,:pTime]) <= 14 &&
  dayofweek(df[i, :dTime]) <= 5 && 12 <= hour(df[i,:dTime]) <= 14
    #removing rides too short or too long, in time or space
    if  60*1000 <= (df[i,:dTime] - df[i,:pTime]).value <= 60*60*1000 &&
    (abs(df[i,:pX] - df[i,:dX]) + abs(df[i,:pY] - df[i,:dY])) >= 200
      return true
    end
  end
  return false
end

# Create the KDTree used to project onto the nearest node.
m=Manhattan()

nodesX = Float64[]
nodesY = Float64[]
nodesID = Int[]

## NO PICKUP and DROPOFF on HIGHWAYS
for (i,p) in enumerate(m.positions)
  nodeInHighway = false
  for n in out_neighbors(m.network,i)
    nodeInHighway = nodeInHighway || m.roadType[i,n] == 8
  end
  for n in in_neighbors(m.network,i)
    nodeInHighway = nodeInHighway || m.roadType[n,i] == 8
  end
  if !nodeInHighway
    push!(nodesX, p.x)
    push!(nodesY, p.y)
    push!(nodesID, i)
  end
end
nodes = Array(Float64, (2,length(nodesX)))
nodes[1,:] = nodesX
nodes[2,:] = nodesY

tree = KDTree(nodes)

data = DataFrame(pTime = DateTime[], dTime=DateTime[], pNode=Int[], dNode=Int[])

for j in 1:12
  df = load("data/reduced_trip_$j.jld", "df")
  println("Opening reduced_trip_$j.jld ...")

  #second : adding rows to dataFrame
  for i in 1:nrow(df)
    if i%1000000 == 0
      println(i, " - ", nrow(data))
    end
    if good_ride(df,i)
      pN = nodesID[knn(tree,[Float64(df[i,:pX]),Float64(df[i,:pY])],1)[1][1]]
      dN = nodesID[knn(tree,[Float64(df[i,:dX]),Float64(df[i,:dY])],1)[1][1]]
      if pN != dN
        push!(data, [df[i,:pTime], df[i,:dTime], pN, dN])
      end
    end
  end
end
save("data/traveltimes.jld", "df", data)
println("traveltimes.jld created. $(nrow(data)) rows after filtering")

max_stars_repo_name: semartin/TaxiDataExtraction

max_issues_repo_name: ericytung/TaxiDataTools

max_forks_repo_name: ericytung/TaxiDataTools


text:
using BinaryBuilder

# Collection of sources required to build SHTOOLS
name = "SHTOOLS"
version = v"4.8"
sources = [
    ArchiveSource("https://github.com/SHTOOLS/SHTOOLS/releases/download/v4.8/SHTOOLS-4.8.tar.gz",
                  "c36fc86810017e544abbfb12f8ddf6f101a1ac8b89856a76d7d9801ffc8dac44"),
    DirectorySource("./bundled"),
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/SHTOOLS-*

# Patch source code
# Add missing C interface for MakeGradientDH
atomic_patch -p1 $WORKSPACE/srcdir/patches/add-cMakeGradientDH.patch
# Don't use libtool
atomic_patch -p0 $WORKSPACE/srcdir/patches/no-libtool.patch

# Build and install static libraries
make fortran -j${nproc} F95FLAGS="-fPIC -O3 -std=gnu"
make fortran-mp -j${nproc} F95FLAGS="-fPIC -O3 -std=gnu"
make install PREFIX=${prefix}

# Create shared libraries
gfortran -shared -o ${libdir}/libSHTOOLS.${dlext} -Wl,$(flagon --whole-archive) ${prefix}/lib/libSHTOOLS.a -Wl,$(flagon --no-whole-archive | cut -d' ' -f1) -lfftw3 -lopenblas -lm
gfortran -fopenmp -shared -o ${libdir}/libSHTOOLS-mp.${dlext} -Wl,$(flagon --whole-archive) ${prefix}/lib/libSHTOOLS-mp.a -Wl,$(flagon --no-whole-archive | cut -d' ' -f1) -lfftw3 -lopenblas -lm
"""

platforms = expand_gfortran_versions(supported_platforms())

# The products that we will ensure are always built
products = [
    LibraryProduct("libSHTOOLS", :libSHTOOLS),
    LibraryProduct("libSHTOOLS-mp", :libSHTOOLS_mp),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("FFTW_jll"),
    Dependency("OpenBLAS32_jll"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; preferred_gcc_version=v"5")

max_stars_repo_name: waralex/Yggdrasil

max_issues_repo_name: waralex/Yggdrasil

max_forks_repo_name: waralex/Yggdrasil


text:
using RCall

R"setwd('~/.julia/packages/linearstep/src')"
R"source('adam.R')"

## Conversion tools
function array_r_to_jl(s4arr)
    dims = R"dim($(s4arr))"
    if typeof(dims) <: RObject{NilSxp}
        dims = [rcopy(R"length($(s4arr))")]
    end
    m = zeros(dims...)
    if length(dims) == 2
        for j in 1:dims[2]
            for i in 1:dims[1]
                m[i, j] = R"$(s4arr)[$(i), $(j)]"
            end
        end
    elseif length(dims) == 1
        for i in 1:dims[1]
            m[i] = R"$(s4arr)[$(i)]"
        end
    else
        error("cannot convert because more than 2 dimensions $(dims)")
    end
    m
end

function rtocohortstate(X)
    if typeof(X) <: RObject
        Xj = rcopy(X)
    else
        Xj = X
    end
    X = array_r_to_jl(Xj[:X])
    pit = array_r_to_jl(Xj[:pit])
    V = Xj[:V]
    ageInDays = Float64(Xj[:ageInDays])
    CohortState(X, pit, V, ageInDays, Xj[:alpha1], Xj[:alpha2])
end

function rtopopstate(rXX)
    if typeof(rXX) <: RObject
        rXXj = rcopy(rXX)
    else
        rXXj = rXX
    end
    X = array_r_to_jl(rXXj[:X])
    pit = array_r_to_jl(rXXj[:pit])
    V = rXXj[:V]
    rage = rXXj[:ageInDays]
    if typeof(rage) <: Array
        ageInDays = convert(Array{Float64,2}, reshape(rage, 1, length(rage)))
    else
        @show rage
        ageInDays = array_r_to_jl(rXXj[:ageInDays])
    end
    PopulationState(X, pit, V, ageInDays, rXXj[:alpha1], rXXj[:alpha2])
end

## Setup for testing
alpha = 0.4
reval("alpha <- $(alpha)")
reval("params <- makeParameters_Adam()")
params = make_parameters_adam()

## Check that the parameters match exactly.
param = make_parameters_adam()
rparam = R"makeParameters_Adam()"
rparamj = rcopy(rparam)
typeof(param[:tau])
@test Set(keys(rparamj)) == Set(keys(param))
for (key, value) in param
    if (typeof(value) <: Number)
        rvalue = rparamj[key]
        @test isapprox(value, rvalue)
    elseif (typeof(value) <: Vector)
        rvalue = rparamj[key]
        if length(value) == length(rvalue)
            if !all(isapprox.(value, rvalue))
                @show key
                @show value
                @show rvalue
            end
        else
            @test length(value) == length(rvalue)
        end
    end
end

## calP_SoI

# Both are length 9
nu = [0.50000000, 0.16666667, 0.14285714, 0.12500000, 0.11111111, 0.11111111,
0.11111111, 0.08333333, 0.00000000]
r = [1.0000000, 0.9900498, 0.9607894, 0.9607894, 0.9607894, 0.9607894,
0.9607894, 0.9607894, 0.9753099]
PSoI = calP_SoI(9, 10, nu, r)
@rput(nu)
@rput(r)
rPSoI = R"calP_SoI(list(N=9, tau=10, nu=nu, r=r))"
rPSoIj = array_r_to_jl(rPSoI)
@test PSoI ≈ rPSoIj

rSuS = array_r_to_jl(rparamj[:SuS])
@test param[:SuS] ≈ rSuS
rPSoI = array_r_to_jl(rparamj[:PSoI])
@test param[:PSoI] ≈ rPSoI

## xiF, defined in make_parameters.
reval("empty = emptyX_Adam()")
reval("xiF_0 <- function(x, p.xi = 2, N = 9, C = 10, D = 0) {
  diff(pbeta(c(0:N) / N, (1 - exp(-p.xi * x)) * C, exp(-p.xi * x) * C + D))
}")
rd = R"xiF_0(0.3)"
jd = xiF_inner(0.3)
@test jd ≈ array_r_to_jl(rd)
rd = R"xiF_0(0)"
jd = xiF_inner(0)
@test jd ≈ array_r_to_jl(rd)

rxi = R"params$xiF(empty$pit, empty$V)"
rxij = array_r_to_jl(rxi)
vecX = emptyX_Adam()
jxi = params[:xiF](vecX.pit, vecX.V)
@test jxi ≈ rxij

## calP_Adam
vecX = emptyX_Adam()
calP = calP_Adam(alpha, vecX.pit, vecX.V, params)
reval("vecX = emptyX_Adam()")
rcalP = R"calP_Adam($(alpha), vecX$pit, vecX$V, params)"
rcalPj = rcopy(rcalP)
# there are NaN in rcalPj. Maybe that's a translation problem?
rcalPj[isnan.(rcalPj)] .= 0
@test calP ≈ rcalPj

## emtpyXX_Adam
cXX = emptyXX_Adam(L=2920)
rXX = R"emptyXX_Adam(L=2920)"
rXXj = rtopopstate(rXX)
@test cXX ≈ rXXj
differences(cXX, rXXj)

## PxX_Adam tests
vecX = emptyX_Adam()
PxX_Adam!(vecX, alpha, params)
reval("vecX = emptyX_Adam()")
rvecX = R"PxX_Adam(alpha, vecX, params)"
rvecXj = rtocohortstate(rvecX)
@test vecX ≈ rvecXj

## cohortXX_Adam tests
cXX = cohortXX_Adam(alpha, params)
# The params variable is defined in R at this point.
rXX = R"cohortXX_Adam($(alpha), params)"
rXXj = rtopopstate(rXX)
@test cXX ≈ rXXj
differences(cXX, rXXj)

## cohort2ages_Adam
cXX = cohortXX_Adam(alpha, params)
XX = cohort2ages_Adam(cXX, params)
reval("cXX = cohortXX_Adam($(alpha), params)")
reval("XX = cohort2ages_Adam(cXX, params)")
rXX = R"XX"
rXXj = rtopopstate(rXX)
# Confirm they started with the same values.
@test rtopopstate(R"cXX") ≈ cXX
# This is the new test.
@test_broken XX ≈ rXXj
@test maximum((XX.X .- rXXj.X).^2) < 0.0015
differences(XX, rXXj)

## ar2stableWave_Adam
XX = ar2stableWave_Adam(alpha, params)
rXX = R"ar2stableWave_Adam($(alpha), $(rparam))"
rXXj = rcopy(rXX)

rXXj = rtopopstate(rXX)
@test_broken rXXj ≈ XX
differences(XX, rXXj)
@test maximum((XX.X .- rXXj.X).^2) < 0.0015

pr1 = ar2pr_Adam(alpha, params)
pr2 = R"ar2pr_Adam(alpha, params)"
@test abs(pr1 - rcopy(pr2)) < 1e-5


text:
export fill_pbo!

# https://www.unavco.org/data/strain-seismic/bsm-data/lib/docs/bottle_format.pdf

bottle_chans = Dict{String, Tuple{String,String}}(
  "BatteryVolts" => ("ABV", "V"),
  "CalOffsetCH0G3" => ("AO0", "{unknown}"),
  "CalOffsetCH1G3" => ("AO1", "{unknown}"),
  "CalOffsetCH2G3" => ("AO2", "{unknown}"),
  "CalOffsetCH3G3" => ("AO3", "{unknown}"),
  "CalStepCH0G2" => ("A02", "{unknown}"),
  "CalStepCH0G3" => ("A03", "{unknown}"),
  "CalStepCH1G2" => ("A12", "{unknown}"),
  "CalStepCH1G3" => ("A13", "{unknown}"),
  "CalStepCH2G2" => ("A22", "{unknown}"),
  "CalStepCH2G3" => ("A23", "{unknown}"),
  "CalStepCH3G2" => ("A32", "{unknown}"),
  "CalStepCH3G3" => ("A33", "{unknown}"),
  "DownholeDegC" => ("KD", "Cel"),
  "LoggerDegC" => ("K1", "Cel"),
  "PowerBoxDegC" => ("K2", "Cel"),
  "PressureKPa" => ("DI", "kPa"),
  "Rainfallmm" => ("R0", "mm"),
  "RTSettingCH0" => ("AR0", "{unknown}"),
  "RTSettingCH1" => ("AR1", "{unknown}"),
  "RTSettingCH2" => ("AR2", "{unknown}"),
  "RTSettingCH3" => ("AR3", "{unknown}"),
  "SolarAmps" => ("ASO", "A"),
  "SystemAmps" => ("ASY", "A")
)

bottle_nets = Dict{String, String}(
  "AIRS" => "MC",
  "TEPE" => "GF",
  "BUY1" => "GF",
  "ESN1" => "GF",
  "TRNT" => "MC",
  "HALK" => "GF",
  "B948" => "ARRA",
  "OLV1" => "MC",
  "OLV2" => "MC",
  "GERD" => "MC",
  "SIV1" => "GF",
  "BOZ1" => "GF",
  "B947" => "ARRA"
  )

function check_bads!(x::AbstractArray, nv::T) where T
  # Check for bad samples
  @inbounds for i in x
    if i == nv
      return true
    end
  end
  return false
end

function channel_guess(str::AbstractString, fs::Float64)
  si = fs >= 1.0 ? 14 : 10
  ei = length(str) - (endswith(str, "_20") ? 3 : 0)

  # name, id, units
  str = str[si:ei]
  if length(str) == 3
    units = "m/m"
  else
    (str, units) = get(bottle_chans, str, ("YY", "{unknown}"))
  end

  # form channel string
  if length(str) == 2
    str = (fs > 1.0 ? "B" : fs > 0.1 ? "L" : "R")*str
  end

  return (str, units)
end

function read_bottle!(S::GphysData, fstr::String, nx_new::Int64, nx_add::Int64, memmap::Bool, strict::Bool, v::Integer)
  buf = BUF.buf
  files = ls(fstr)

  for file in files
    io = memmap ? IOBuffer(Mmap.mmap(file)) : open(file, "r")

    # Read header ============================================================
    fastskip(io, 8)
    t0 = round(Int64, fastread(io, Float64)*sμ)
    dt = fastread(io, Float32)
    nx = fastread(io, Int32)
    ty = fastread(io, Int32)
    nv = fastread(io, Int32)
    fastskip(io, 8)
    fs = 1.0/dt
    v > 2 && println("t0 = ", t0, ", fs = ", fs, ", nx = ", nx, ", ty = ", ty, ", nv = ", nv)

    # Read data ==============================================================
    T = ty == 0 ? Int16 : ty == 1 ? Int32 : Float32
    nb = nx*sizeof(T)
    checkbuf_8!(buf, nb)
    fast_readbytes!(io, buf, nb)
    close(io)

    # Try to create an ID from the file name =================================
    # Assumes fname SSSSyyJJJ... (SSSS = station, yy = year, JJJ = Julian day)
    fname = splitdir(file)[2]
    sta = fname[1:4]
    (cha, units) = channel_guess(fname, fs)

    # find relevant entry in station data
    net = get(bottle_nets, sta, "PB")
    id = net * "." * sta * ".." * cha

    # Load into S ============================================================
    i = findid(id, S.id)
    if strict
      i = channel_match(S, i, fs)
    end
    if i == 0

      # Create C.x
      x = Array{Float64,1}(undef, max(nx_new, nx))
      os = 1

      C = SeisChannel()
      setfield!(C, :id, id)
      setfield!(C, :fs, fs)
      setfield!(C, :units, units)
      mk_t!(C, nx, t0)
      setfield!(C, :x, x)
      push!(S, C)
    else
      xi = S.t[i][end,1]
      x = getindex(getfield(S, :x), i)
      check_for_gap!(S, i, t0, nx, v)
      Lx = length(x)
      if xi + nx > Lx
        resize!(x, Lx + max(nx_add, nx))
      end
      os = xi + 1
    end

    # Check for null values
    nv = T(nv)
    y = reinterpret(T, buf)
    b = T == Int16 ? false : check_bads!(y, nv)
    if b
      j = os
      @inbounds for i = 1:nx
        if y[i] == nv
          x[j] = NaN
        else
          x[j] = y[i]
        end
        j += 1
      end
    else
      copyto!(x, os, y, 1, nx)
    end
  end
  trunc_x!(S)
  resize!(buf, 65535)

  return nothing
end

function read_bottle(fstr::String, nx_new::Int64, nx_add::Int64, memmap::Bool, strict::Bool, v::Integer)

  S = SeisData()
  read_bottle!(S, fstr, nx_new, nx_add, memmap, strict, v)
  return S
end

"""
    fill_pbo!(S)

Attempt to fill `:name` and `:loc` fields of S using station names (second field of S.id) cross-referenced against a PBO station info file.
"""
function fill_pbo!(S::GphysData)
  sta_data = readdlm(path * "/Formats/PBO_bsm_coords.txt", ',', comments=false)
  sta_data[:,2] .= strip.(sta_data[:,2])
  sta_data[:,6] .= strip.(sta_data[:,6])
  n_sta = size(sta_data, 1)
  for i = 1:S.n
    sta = split(S.id[i], '.')[2]
    for j = 1:n_sta
      if sta_data[j, 1] == sta
        S.name[i] = String(sta_data[j,2])
        lat = sta_data[j,3]
        lon = sta_data[j,4]
        el = sta_data[j,5]
        S.loc[i] = GeoLoc(lat = lat, lon = lon, el = el)
        break
      end
    end
  end
  return nothing
end


text:
#################### Binary MCMC Model Composition ####################

#################### Types and Constructors ####################

const BMC3Form = Union{Int, Vector{Vector{Int}}}

mutable struct BMC3Tune{F<:BMC3Form} <: SamplerTune
  logf::Union{Function, Missing}
  k::F

  BMC3Tune{F}() where {F<:BMC3Form} = new{F}()

  BMC3Tune{F}(x::Vector, k::F) where {F<:BMC3Form} = new{F}(missing, k)

  BMC3Tune{F}(x::Vector, k::F, logf::Function) where {F<:BMC3Form} =
    new{F}(logf, k)
end

const BMC3IntVariate = SamplerVariate{BMC3Tune{Int}}
const BMC3VecVariate = SamplerVariate{BMC3Tune{Vector{Vector{Int}}}}

BMC3Variate(x::Vector, logf::Function; k::F=1) where {F<:BMC3Form} =
  SamplerVariate{BMC3Tune{F}}(x, k, logf)

function validate(v::BMC3IntVariate)
  n = length(v)
  v.tune.k <= n || throw(ArgumentError("k exceeds variate length $n"))
  validatebinary(v)
end

function validate(v::BMC3VecVariate)
  n = length(v)
  mapreduce(maximum, max, v.tune.k) <= n ||
    throw(ArgumentError("indices in k exceed variate length $n"))
  validatebinary(v)
end

#################### Sampler Constructor ####################

function BMC3(params::ElementOrVector{Symbol}; k::F=1) where {F<:BMC3Form}
  samplerfx = function(model::Model, block::Integer)
    block = SamplingBlock(model, block)
    v = SamplerVariate(block, k)
    sample!(v, x -> logpdf!(block, x))
    relist(block, v)
  end
  Sampler(params, samplerfx, BMC3Tune{F}())
end

#################### Sampling Functions ####################

sample!(v::SamplerVariate{BMC3Tune{F}}) where {F<:BMC3Form} = sample!(v, v.tune.logf)

function sample!(v::SamplerVariate{BMC3Tune{F}}, logf::Function) where {F<:BMC3Form}
  x = v[:]
  idx = randind(v)
  x[idx] = 1.0 .- v[idx]
  if rand() < exp(logf(x) - logf(v.value))
    v[:] = x
  end
  v
end

randind(v::BMC3IntVariate) = sample(1:length(v), v.tune.k, replace=false)
randind(v::BMC3VecVariate) = sample(v.tune.k)


text:
using PrettyLinearAlgebra
using Test
using LinearAlgebra

@testset "PrettyLinearAlgebra.jl" begin
    v1 = [1, 5, 9, 99]
    m1 = [1 5 9; 4 3 1; 33 44 55]

    @test matex(v1) == "\\begin{bmatrix} 1 \\\\ 5 \\\\ 9 \\\\ 99\\end{bmatrix}"
    @test matex(m1) == "\\begin{bmatrix}1 & 5 & 9 \\\\ 4 & 3 & 1 \\\\ 33 & 44 & 55\\end{bmatrix}"
end


text:
# initial approach
function do_z_expected_ug_bounded(method, model_key)

    # orig model
    iterator = Fd.val(:D) |> enumerate |> collect
    @threads for (exp, D) in iterator
        thid = threadid()

        ## -------------------------------------------------------------------
        # Excluded because of unfeasibility
        exp in IGNORE_EXPS && continue

        ## -------------------------------------------------------------------
        # handle cache
        is_cached(;method, exp) && continue

        ## -------------------------------------------------------------------
        # ep
        me_params = lglob(iJR, :maxent, :params)
        @extract me_params: alpha epsconv maxiter damp maxvar minvar

        # approach
        convth = 0.05
        converr = nothing

        # prepare model
        model = iJR.load_model(model_key, exp)
        biomidx = ChU.rxnindex(model, iJR.BIOMASS_IDER)
        M, N = size(model)
        exp_growth = Fd.val(:D, exp)
        growth_ub = ChU.ub(model, iJR.BIOMASS_IDER)
        ChU.ub!(model, biomidx, growth_ub * 1.1) # open a bit helps EP

        lock(WLOCK) do
            nzabs_range = ChU.nzabs_range(model.S)
            @info("Starting... ",
                exp, method,
                size(model),
                nzabs_range,
                thid
            ); println()
        end

        # simulation globals
        epouts = Dict()
        beta_vec = zeros(N)
        epout_seed = nothing
        epout = nothing
        betas = [0.0; 10.0.^(3:0.05:15)]
        nan_beta = first(betas)

        for approach in [:log_approach, :linear_approach]

            lock(WLOCK) do
                @info("Starting", approach,
                    exp, method,
                    length(epouts),
                    thid
                ); println()
            end

            for beta in betas

                nan_beta = beta
                haskey(epouts, beta) && continue

                beta_vec[biomidx] = beta
                try
                    epout = ChEP.maxent_ep(model;
                        beta_vec, alpha, damp, maxiter,
                        epsconv, maxvar, minvar,
                        verbose = false,
                        solution = epout_seed
                    )
                catch err
                    lock(WLOCK) do
                        @info("Log approach errors", exp, beta, err, thid);
                        println()
                    end
                    break
                end

                # stop conditions
                biom_avPME = ChU.av(model, epout, biomidx)
                converr = abs(biom_avPME - exp_growth) / exp_growth

                # Check failing
                if (isnan(biom_avPME) || biom_avPME == 0.0)
                    lock(WLOCK) do
                        @info("Log approach fails", exp, beta, biom_avPME, thid);
                        println()
                    end
                    break
                end

                # info
                lock(WLOCK) do
                    @info("Results", approach,
                        exp, method, beta,
                        exp_growth, growth_ub, biom_avPME,
                        converr, convth,
                        length(epouts),
                        thid
                    ); println()
                end

                # updating
                epout_seed = epouts[beta] = epout

                # convergence
                conv = converr < convth || biom_avPME > exp_growth
                conv && break

            end # for betas

            # setup lineal approach
            isempty(epouts) && error("LOG FAILS TOO SOON")
            last_beta = maximum(keys(epouts))
            step = abs(last_beta - nan_beta) / 100.0
            iszero(step) && break
            betas = range(last_beta, 1.0e15; step)

        end # for approach

        # further conv
        let
            exp_beta = maximum(keys(epouts))

            @info("Further conv",
                exp, method,
                length(epouts),
                exp_beta,
                thid
            ); println()

            beta_vec[biomidx] = exp_beta
            epouts[exp_beta] = ChEP.maxent_ep(model;
                beta_vec, alpha, damp,
                maxiter = max(maxiter * 3, 5000),
                epsconv = epsconv * 0.1,
                maxvar, minvar,
                verbose = false,
                solution = epouts[exp_beta]
            )
        end

        # saving
        lock(WLOCK) do
            datfile = dat_file(;method, exp)
            dat = Dict()
            dat[:model] = model |> ChU.compressed_model
            dat[:epouts] = epouts
            dat[:exp_beta] = maximum(keys(epouts))
            sdat(iJR, dat, datfile)

            @info("Finished", exp, method,
                length(epouts),
                thid
            ); println()
        end

    end # for (exp, D)

end