text:
## Overridden behaviours
"""
    WrapperLogType

Abstract wrapper type to override the default backend to which a type would be
logged in TensorBoard.

Every concrete implementation should implement the method `content` for the
concrete type, returning the content.
"""
abstract type WrapperLogType end

"""
    content(x::WrapperLogType)

Returns the wrapped data inside of `x`.
"""
content(x::WrapperLogType) = @error "You should define `content($(typeof(x)))`"

# When we hit a wrapped type, we dispatch based on it's type.
preprocess(name, val::WrapperLogType, data) = push!(data, name=>val)

# When logging to text or console, unwrap the LogType to fix formatting.
Base.show(io::IO, mime::AbstractString, x::WrapperLogType) =
    Base.show(io, mime, content(x))

########## For things going to LogAudio ##############################
"""
    TBAudios(data, samplerate)

Forces elements of Array `data` to be serialized as Audio to TensorBoard.
"""
struct TBAudios <:WrapperLogType
    data::AbstractArray
    samplerate::Real
end
content(x::TBAudios) = x.data
function preprocess(name, val::TBAudios, data)
    for (n, audio) in enumerate(val.data)
        preprocess(name*"/$n", TBAudio(audio, val.samplerate), data)
    end
end

"""
    TBAudio(data, samplerate)

Forces `data` to be serialized as Audio to TensorBoard.
"""
struct TBAudio <: WrapperLogType
    data::AbstractArray
    samplerate::Real
end
content(x::TBAudio) = x.data
preprocess(name, val::TBAudio, data) = push!(data, name=>val)
summary_impl(name, val::TBAudio) = audio_summary(name, val.data, val.samplerate)

########## For things going to LogImage ########################
"""
    TBImages(data, format)

Forces elements of Array `data` to be serialized as an Image to TensorBoard.
"""
struct TBImages <:WrapperLogType
    data::AbstractArray
    format::ImageFormat
end

"""
    TBImages(image)

Forces elements of Array `data` to be serialized as images to TensorBoard.
`showable("image/png", image)` must be true.
"""
TBImages(image) = TBImages(image, PngImage)

content(x::TBImages) = x.data
function preprocess(name, val::TBImages, data)
    for (n, img) in enumerate(val.data)
        preprocess(name*"/$n", TBImage(img, val.format), data)
    end
end

"""
    TBImage(data, format)

Forces `data` to be serialized as an Image to TensorBoard.
"""
struct TBImage{T} <: WrapperLogType
    data::T
    format::ImageFormat
end
"""
    TBImage(image)

Renders `image` to PNG and logs it as an image to TensorBoard.
`showable("image/png", image)` must be true.
"""
TBImage(image) = TBImage(image, PngImage)

content(x::TBImage) = x.data

function Base.convert(T::Type{PngImage}, img::TBImage)
    if img.format == PNG
        if !showable("image/png", img.data)
            @error "cannot log as an image object" obj=img.data
            return
        end
        pb=PipeBuffer()
        show(pb, "image/png", img.data)
        return PngImage(pb)
    else
        @error "Unknown format for datatype." format=img.format datatype=typeof(img.data)
        return nothing
    end
end

function Base.convert(T::Type{PngImage}, img::TBImage{<:AbstractArray})
    imgArray = img.data
    format   = img.format

    #unpack RGB, RGBA value to channels using channelview
    imgArray = channelview(imgArray)
    #if data contains integer numbers, scale them to 0-255
    if isa(first(imgArray), Integer)
        imgArray = (imgArray./255)
    end

    #convert all values to `Float64` for uniformity and scale all values to 0-1
    imgArray = Float64.(imgArray)
    imgArray = (imgArray./(max(maximum(imgArray), 1)))
    imgArray = convert_to_CHW(imgArray, format)

    channels, height, width = size(imgArray)
    color = ColorType_from_nchannels(channels)

    #if it is a single channel Array, convert it to HW
    if color == Gray
        imgArray = imgArray[1, :, :]
    end
    #convert Array to PNG pass it to image_summary
    img = colorview(color, imgArray)
    return convert(PngImage, img)
end

function preprocess(name, val::TBImage, data)
    imgArray = val.data
    format = val.format
    obsdim = obs_dim(format)

    # obs_dim of PNG format is 0
    if iszero(obsdim)
        preprocess(name, convert(PngImage, val), data)
    else
        imgArray = channelview(imgArray)
        @assert ndims(imgArray) == expected_ndims(format)

        format = strip_obs[format]
        for n in 1:size(imgArray, obsdim)
            nth_img = selectdim(imgArray, obsdim, n)
            preprocess(name*"/$n", convert(PngImage, TBImage(nth_img, format)), data)
        end
    end
    return data
end

########## For things going to LogText ##############################
"""
    TBText(data)

Forces `data` to be serialized as text to TensorBoard.
"""
struct TBText <: WrapperLogType data; end
content(x::TBText) = x.data
summary_impl(name, val::TBText) = text_summary(name, val.data)

"""
    RealArrayWrapperLogType{T,N} <: WrapperLogType

A wrapper type that accepts only arrays and converts complex-valued arrays
into two real value ones
"""
abstract type RealArrayWrapperLogType{T,N} <: WrapperLogType end

########## For things going to LogHistograms ########################
"""
    TBHistogram(data)

Forces `data` to be serialized as an histogram to TensorBoard.
"""
struct TBHistogram{T,N} <: RealArrayWrapperLogType{T,N}
    data::AbstractArray{T,N};
end
content(x::TBHistogram) = x.data
preprocess(name, val::TBHistogram{T,N}, data) where {T<:Complex,N} =
    push!(data, name*"/re"=>TBHistogram(real.(content(val))), name*"/im"=>TBHistogram(imag.(content(val))))
function summary_impl(name, val::TBHistogram)
    hvals = fit(Histogram, val.data)
    return histogram_summary(name, collect(hvals.edges[1]), hvals.weights)
end

"""
    TBVector(data)

Forces `data` to be serialized as a vector in the histogram backend of
TensorBoard.
"""
struct TBVector{T,N} <: RealArrayWrapperLogType{T,N} data::AbstractArray{T,N}; end
content(x::TBVector) = x.data
preprocess(name, val::TBVector{T,N}, data) where {T<:Complex,N} =
    push!(data, name*"/re"=>TBVector(real.(content(val))), name*"/im"=>TBVector(imag.(content(val))))
summary_impl(name, val::TBVector) = histogram_summary(name, collect(0:length(val.data)),
                                                            val.data)


text:
T = 250
σ = 0.2
t = range(0, 20, length = T)
s = sin.(t)
y = s + randn(T) * σ

h = (t[end] - t[1])/T
H = [1.0, 0.0, 0.0]
A = [1.0 h   0.5*h^2;
     0.0 1.0 h;
     0.0 0.0 1.0]

m = Model(y, H, A, Matrix(0.01*I, 3, 3), σ^2, [0.0, 0.0, 0.0], Matrix(1000.0I, 3, 3))
kf = kalmanfilter(m)
ks = kalmansmoother(m, kf.priors, kf.posteriors)

using Plots

scatter(y, label = "Measurements", legend = :bottom, ylim = [-3, 2],
        markerstrokewidth = 0, markercolor = "black", foreground_color_legend = nothing)
plot!(s, label = "Noiseless signal", linewidth = 3, color = "grey")
plot!(kf.priors.μ, label = "Prediction", linewidth = 2, color = "blue")
plot!(kf.posteriors.μ, label = "Filter", linewidth = 2, color = "red")
plot!(ks.μ, label = "Smoother", linewidth = 2, color = "green")

scatter(y, label = "Measurements", legend = :bottom, ylim = [-3, 2],
        markerstrokewidth = 0, markercolor = "black", foreground_color_legend = nothing)
plot!(s, label = "Noiseless signal", linewidth = 3, color = "grey")
plot!(ks.x, label = ["Position" "Velocity" "Acceleration"], linewidth = 2, color = ["blue" "red" "green"])

ext: jl


text:
module OpenGene

##############################################################################
##
## Dependencies
##
##############################################################################

using DataFrames

# make it compatible for different version of Julia
include("compat.jl")

# all code in Base should have no code dependency on others
include("Base/Base.jl")

# Settings
include("Settings/Settings.jl")

# IO code should only depend on Base
include("IO/IO.jl")

# module OpenGene.Algorithm
include("Algorithm/Algorithm.jl")

# module OpenGene.Reference
include("Reference/Reference.jl")

include("utils/print.jl")

end # module


text:
using Gadfly
using RDatasets

plot(dataset("lattice", "singer"), x="VoicePart", y="Height", Geom.beeswarm)


text:
module Models

using DataFrames
using CSV
using Statistics
using LinearAlgebra
using Distributions
using Random
using ImageFiltering
import GeoPhyInv.Smooth
import GeoPhyInv.Interpolation

"""
Store reference model parameters
"""
mutable struct Seismic_ref{T<:Real}
	vp::T
	vs::T
	ρ::T
	K::T
	μ::T
	KI::T
	μI::T
	ρI::T
end

"""
Data type fo represent a seismic model.
A contrast function for a model m is given by ``χ(m) = \frac{m}{m0}-1``.

# Fields

* `vp0::Vector{Float64}` : [vpmin, vpmax]
* `vs0::Vector{Float64}` : [vsmin, vsmax]
* `ρ0::Vector{Float64}` : [ρmin, ρmax]
* `χvp::Array{Float64,2}` : two-dimensional contrast model (χ) for vp, for e.g., zeros(length(mgrid[1]), length(mgrid[2]))
* `χvs::Array{Float64}` : two-dimensional contrast model (χ) for vs, for e.g., zeros(length(mgrid[1]), length(mgrid[2]))
* `χρ::Array{Float64}` : two-dimensional contrast model (χ) for density, for e.g., zeros(length(mgrid[1]), length(mgrid[2]))
* `mgrid` : array of ranges to determine the dimensions of models
"""
mutable struct Seismic
	vp0::Vector{Float64}
	vs0::Vector{Float64}
	ρ0::Vector{Float64}
	χvp::Array{Float64,2}
	χvs::Array{Float64,2}
	χρ::Array{Float64,2}
	mgrid::Vector{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}
	# all derived fields from previous
	K0::Vector{Float64}
	μ0::Vector{Float64}
	KI0::Vector{Float64}
	μI0::Vector{Float64}
	ρI0::Vector{Float64}
	ref::Seismic_ref{Float64}
	"adding conditions that are to be false while construction"
	Seismic(vp0,vs0,ρ0,χvp,χvs,χρ,mgrid,K0,μ0,KI0,μI0,ρI0,ref) =
		any([
		     any(vp0.<0.0),
		     any(vs0.<0.0),
		     any(ρ0.<0.0),
		     #all([all(vp0 .≠ 0.0), any(χ(χvp,vp0,-1) .< vp0[1])]), # check vp bounds
		     #all([all(vp0 .≠ 0.0), any(χ(χvp,vp0,-1) .> vp0[2])]), # check vp bounds
		     #all([all(vs0 .≠ 0.0), any(χ(χvs,vs0,-1) .< vs0[1])]), # check vs bounds
		     #all([all(vs0 .≠ 0.0), any(χ(χvs,vs0,-1) .> vs0[2])]), # check vs bounds
		     #all([all(ρ0 .≠ 0.0), any(χ(χρ,ρ0,-1) .< ρ0[1])]), # check ρ bounds
		     #all([all(ρ0 .≠ 0.0), any(χ(χρ,ρ0,-1) .> ρ0[2])]), # check ρ bounds
		     size(χvp) != (length(mgrid[1]), length(mgrid[2])), # dimension check
		     size(χvs) != (length(mgrid[1]), length(mgrid[2])), # dimension check
		     size(χρ) != (length(mgrid[1]), length(mgrid[2])) # dimension check
		    ]) ?
		       error("error in Seismic construction") : new(vp0,vs0,ρ0,χvp,χvs,χρ,mgrid,K0,μ0,KI0,μI0,ρI0,ref)
end

# check whether a seismic model is bounded
function isbounded(mod::Seismic)
	@warn "need to be implemented"
	result=true
end

# method to create Seismic object without using derived fields
function Seismic(vp0, vs0, ρ0, χvp, χvs, χρ,
		 mgrid::Vector{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}})
	mod=Seismic(vp0, vs0, ρ0, χvp, χvs, χρ, mgrid,
	    zeros(2), zeros(2), zeros(2), zeros(2), zeros(2),
	    Seismic_ref(zeros(8)...))
	update_derived!(mod)
	return mod
end

# update the derived fields of the Seisimc model
function update_derived!(mod::Seismic)
	"note that: mean(K0) ≠ 1/mean(KI0)"
	mod.K0 = mod.vp0 .* mod.vp0 .* mod.ρ0;
	mod.KI0 = reverse(inv.(mod.K0),dims=1);
	ρ0 = mod.ρ0;
	mod.ρI0 = reverse(inv.(ρ0),dims=1);
	mod.μ0 = mod.vs0 .* mod.vs0 .* mod.ρ0;
	mod.μI0 = reverse(inv.(mod.μ0),dims=1);
	mod.ref = Seismic_ref(Statistics.mean.([mod.vp0,mod.vs0,mod.ρ0,mod.K0,mod.μ0,mod.KI0,mod.μI0,mod.ρI0])...)
end

"""
Print information about `Seismic`
"""
function Base.print(mod::Seismic, name::String="")
	println("\tSeismic Model:\t",name)
	println("\t> number of samples:\t","x\t",length(mod.mgrid[2]),"\tz\t",length(mod.mgrid[1]))
	println("\t> sampling intervals:\t","x\t",step(mod.mgrid[2]),"\tz\t",step(mod.mgrid[1]))
	println("\t> vp:\t","min\t",minimum(Seismic_get(mod,:vp)),"\tmax\t",maximum(Seismic_get(mod,:vp)))
	println("\t> vp bounds:\t","min\t",mod.vp0[1],"\tmax\t",mod.vp0[2])
	println("\t> ρ:\t","min\t",minimum(Seismic_get(mod,:ρ)),"\tmax\t",maximum(Seismic_get(mod,:ρ)))
	println("\t> ρ bounds:\t","min\t",mod.ρ0[1],"\tmax\t",mod.ρ0[2])
end

"""
Return medium property bounds based on maximum and minimum values of the array and frac.
The bounds cannot be less than zero
"""
function bounds(mod::Array{Float64}, frac::Float64=0.1)
	any(mod .< 0.0) && error("model values less than zero")
	bounds=zeros(2)
	bound=frac*Statistics.mean(mod)
	bounds[1] = ((minimum(mod) - bound)<0.0) ? 0.0 : (minimum(mod) - bound)
	bounds[2] = maximum(mod)+bound
	return bounds
end

function adjust_bounds!(mod::Seismic,mod0::Seismic)
	adjust_bounds!(mod,mod0.vp0,mod0.vs0,mod0.ρ0)
end

function adjust_bounds!(mod::Seismic,vp0::Vector{T},vs0::Vector{T},ρ0::Vector{T}) where {T<:Real}
	copyto!(mod.vp0, Float64.(vp0))
	copyto!(mod.vs0, Float64.(vs0))
	copyto!(mod.ρ0, Float64.(ρ0))
	update_derived!(mod)
end

"""
Adjust the bounds and hence the reference values.
Since the reference values are adjust the χ fields should also be changed
"""
function adjust_bounds!(mod::Seismic, frac::Float64=0.1)
	for f in [:χvp, :χρ, :χvs]
		f0 = Symbol((replace("$(f)", "χ" => "")),0)
		fm = Symbol((replace("$(f)", "χ" => "")))
		m = Seismic_get(mod, fm)
		m0 = bounds(m, frac)
		setfield!(mod, f0, m0)
		setfield!(mod, f, χ(m, mean(m0), 1))
	end
	update_derived!(mod)
	return mod
end

"""
Return `Seismic` with zeros everywhere;
this method is used for preallocation.

# Arguments
* `mgrid` : used for sizes of χ fields
"""
function Seismic_zeros(mgrid::Vector{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}})
	return Seismic(fill(0.0,2), fill(0.0,2), fill(0.0,2),
		zeros(length(mgrid[1]), length(mgrid[2])), zeros(length(mgrid[1]), length(mgrid[2])),
		zeros(length(mgrid[1]), length(mgrid[2])), deepcopy(mgrid))
end

function Base.fill!(mod::Seismic, k::Float64=0.0)
	fill!(mod.χvp,k)
	fill!(mod.χρ,k)
	fill!(mod.χvs,k)
	return mod
end

"""
Return true if a `Seismic` model is just allocated with zeros.
"""
function Base.iszero(mod::Seismic)
	return any((mod.vp0 .* mod.ρ0) .== 0.0) ? true : false
end

"""
Return a similar model to the input model, used for allocation.
"""
function Base.similar(mod::Seismic)
	modout=deepcopy(mod)
	fill!(modout, 0.0)
	return modout
end

"Compare if two `Seismic` models are equal"
function Base.isequal(mod1::T, mod2::T) where {T<:Union{Seismic,Seismic_ref}}
	fnames = fieldnames(T)
	vec=[(isequal(getfield(mod1, name),getfield(mod2, name))) for name in fnames]
	return all(vec)
end

"""
Return if two `Seismic` models have same dimensions and bounds.
"""
function Base.isapprox(mod1::Seismic, mod2::Seismic)
	vec=([(mod1.vp0==mod2.vp0), (mod1.vs0==mod2.vs0), (mod1.ρ0==mod2.ρ0),
       		(size(mod1.χvp)==size(mod2.χvp)),
       		(size(mod1.χvs)==size(mod2.χvs)),
       		(size(mod1.χρ)==size(mod2.χρ)),
		isequal(mod1.mgrid, mod2.mgrid),
		])
	return all(vec)
end

"""
Copy for `Seismic` models. The models should have same bounds and sizes.
Doesn't allocate any memory.
"""
function Base.copyto!(modo::Seismic, mod::Seismic)
	if(isapprox(modo, mod))
		for f in [:χvp, :χρ, :χvs]
			modof=getfield(modo,f)
			modf=getfield(mod,f)
			copyto!(modof,modf)
		end
		return modo
	else
		error("attempt to copy dissimilar models")
	end
end

function Seismic_get(mod::Seismic, attrib::Symbol)
	# allocate
	rout=zeros(length(mod.mgrid[1]), length(mod.mgrid[2]))
	Seismic_get!(rout, mod, [attrib])
	return rout
end

"""
Get other dependent model parameters of a seismic model
that are not present in `Seismic`.

* `:ρI` : inverse of density
* `:Zp` : P-wave impedance
"""
function Seismic_get!(x, mod::Seismic, attribvec::Vector{Symbol})
	ρ=mod.χρ; χ!(ρ, mod.ref.ρ,-1) # undo it later
	vp=mod.χvp; χ!(vp, mod.ref.vp,-1) # undo it later
	vs=mod.χvs; χ!(vs, mod.ref.vs,-1) # undo it later
	nznx=length(ρ)
	(length(x)≠(count(attribvec.≠ :null)*nznx)) &&  error("size x")
	i0=0;
	for attrib in attribvec
		if(attrib ≠:null)
		if(attrib == :ρI)
			@inbounds for i in 1:nznx; x[i0+i]=inv(ρ[i]); end
		elseif(attrib == :ρ)
			@inbounds for i in 1:nznx; x[i0+i]=ρ[i]; end
		elseif(attrib == :vp)
			@inbounds for i in 1:nznx; x[i0+i]=vp[i]; end
		elseif(attrib == :vs)
			@inbounds for i in 1:nznx; x[i0+i]=vs[i]; end
		elseif(attrib == :χρI)
			@inbounds for i in 1:nznx; x[i0+i]=χ(inv(ρ[i]),mod.ref.ρI,1); end
		elseif(attrib == :χρ)
			@inbounds for i in 1:nznx; x[i0+i]=χ(ρ[i],mod.ref.ρ,1); end
		elseif(attrib == :χvp)
			@inbounds for i in 1:nznx; x[i0+i]=χ(vp[i],mod.ref.vp,1); end
		elseif(attrib == :χK)
			@inbounds for i in 1:nznx; x[i0+i]=χ(vp[i]*vp[i]*ρ[i],mod.ref.K,1); end
		elseif(attrib == :χμ)
			@inbounds for i in 1:nznx; x[i0+i]=χ(vs[i]*vs[i]*ρ[i],mod.ref.μ,1); end
		elseif(attrib == :χKI)
			@inbounds for i in 1:nznx; x[i0+i]=χ(inv(vp[i]*vp[i]*ρ[i]),mod.ref.KI,1); end
		elseif(attrib == :KI)
			@inbounds for i in 1:nznx; x[i0+i]=inv(vp[i]*vp[i]*ρ[i]); end
		elseif(attrib == :K)
			@inbounds for i in 1:nznx; x[i0+i]=vp[i]*vp[i]*ρ[i]; end
		elseif(attrib == :Zp)
			@inbounds for i in 1:nznx; x[i0+i]=vp[i]*ρ[i]; end
		else
			error("invalid attrib")
		end
		i0+=nznx
		end
	end
	χ!(ρ, mod.ref.ρ,1)
	χ!(vp, mod.ref.vp,1)
	χ!(vs, mod.ref.vs,1)
	return x
end

include("reparameterize.jl")

include("chainrule.jl")

"""
Return dimensionless contrast model parameter using the
reference value.

# Arguments
* `mod::Array{Float64}` : subsurface parameter
* `mod0::Vector{Float64}` : reference value is mean of this vector
* `flag::Int64=1` :
"""
function χ(mod::AbstractArray{T}, mod0::T, flag::Int64=1) where {T<:Real}
	mod_out=copy(mod)
	χ!(mod_out, mod0, flag)
	return mod_out
end
function χ!(mod::AbstractArray{T}, mod0::T, flag::Int64=1) where {T<:Real}
	m0 = mod0
	if(flag == 1)
		if(iszero(m0))
			return nothing
		else
			@inbounds for i in eachindex(mod)
				mod[i] = ((mod[i] - m0) * inv(m0))
			end
			return mod
		end
	elseif(flag == -1)
		@inbounds for i in eachindex(mod)
			mod[i] = mod[i] * m0 + m0
		end
		return mod
	end
end

function χ(m::T, m0::T, flag::Int64=1) where {T<:Real}
	if(flag==1)
		return (m-m0) * inv(m0)
	else(flag==-1)
		return (m*m0+m0)
	end
end

"""
Gradients
Return contrast model parameter using the
reference value.
"""
function χg(mod::AbstractArray{T}, mod0::T, flag::Int64=1) where {T<:Real}
	mod_out=copy(mod)
	χg!(mod_out, mod0, flag)
	return mod_out
end
function χg!(mod::AbstractArray{T}, mod0::T, flag::Int64=1) where {T<:Real}
	if(flag == 1)
		rmul!(mod, mod0)
	elseif(flag == -1)
		m0 = inv(mod0)
		rmul!(mod, m0)
	end
	return mod
end

function χg(m::T, m0::T, flag::Int64=1) where {T<:Real}
	if(flag==1)
		return m*m0
	else(flag==-1)
		return m*inv(m0)
	end
end

"""
In-place method to add features to the input model.

# Arguments
* `mod::Seismic` : model that is modified

# Keyword Arguments

* `point_loc::Vector{Float64}=[0., 0.,]` : approx location of point pert
* `point_pert::Float64=0.0` : perturbation at the point scatterer
* `ellip_loc::Vector{Float64}=nothing` : location of center of perturbation, [z, x]
* `ellip_rad::Float64=0.0` : size of elliptic perturbation
* `ellip_pert::Float64=0.1` : perturbation inside the ellipse
* `ellip_α=0.0` : rotate the ellipse
* `rect_loc::Array{Float64}=nothing` : rectangle location, [zmin, xmin, zmax, xmax]
* `rect_pert::Float64=0.1` : perturbation in a rectangle
* `constant_pert::Float64=0.0` : constant perturbation
* `randn_pert::Float64=0.0` : percentage of reference values for additive random noise
* `fields::Vector{Symbol}=[:χvp,:χρ,:χvs]` : which fields are to be modified?
* `onlyin` : `mod` is modified only when field values are in these ranges
"""
function Seismic_addon!(mod::Seismic;
		       point_loc=[0., 0.,],
		       point_pert::Float64=0.0,
		       ellip_loc=[0., 0.,],
		       ellip_rad=0.0,
		       ellip_pert::Float64=0.0,
		       ellip_α=0.0,
		       rect_loc=[0., 0., 0., 0.,],
		       rect_pert::Float64=0.0,
		       constant_pert::Float64=0.0,
		       randn_perc::Real=0.0,
		       fields::Vector{Symbol}=[:χvp,:χρ,:χvs],
		       onlyin::Vector{Vector{Float64}}=[[typemin(Float64), typemax(Float64)] for i in fields]
		       )
	rect_loc=convert.(Float64,rect_loc);
	ellip_loc=convert.(Float64,ellip_loc);

	temp = zeros(length(mod.mgrid[1]), length(mod.mgrid[2]))

	ipointlocx = Interpolation.indminn(mod.mgrid[2], Float64(point_loc[2]), 1)[1]
	ipointlocz = Interpolation.indminn(mod.mgrid[1], Float64(point_loc[1]), 1)[1]
	temp[ipointlocz, ipointlocx] += point_pert

	if(!(ellip_pert == 0.0))
		α=ellip_α*pi/180.
		# circle or ellipse
		rads= (length(ellip_rad)==1) ? [ellip_rad[1],ellip_rad[1]] : [ellip_rad[1],ellip_rad[2]]

		temp += [(((((mod.mgrid[2][ix]-ellip_loc[2])*cos(α)+(mod.mgrid[1][iz]-ellip_loc[1])*sin(α))^2*inv(rads[1]^2) +
	      ((-mod.mgrid[1][iz]+ellip_loc[1])*cos(α)+(mod.mgrid[2][ix]-ellip_loc[2])*sin(α))^2*inv(rads[2]^2)) <= 1.) ? ellip_pert : 0.0)  for
	   		iz in 1:length(mod.mgrid[1]), ix in 1:length(mod.mgrid[2]) ]
	end
	if(!(rect_pert == 0.0))
		temp += [
			(((mod.mgrid[2][ix]-rect_loc[4]) * (mod.mgrid[2][ix]-rect_loc[2]) < 0.0) &
			((mod.mgrid[1][iz]-rect_loc[3]) * (mod.mgrid[1][iz]-rect_loc[1]) < 0.0)) ?
			rect_pert : 0.0  for
			iz in 1:length(mod.mgrid[1]), ix in 1:length(mod.mgrid[2]) ]
	end
	if(!(constant_pert == 0.0))
		temp .+= constant_pert
	end

	for (iff,field) in enumerate(fields)
		m=getfield(mod, field)
		for i in eachindex(m)
			if(((m[i]-onlyin[iff][1])*(m[i]-onlyin[iff][2]))<0.0)
				m[i] += temp[i]
			end
		end
	end

	# random noise
	for field in fields
		f0 = Symbol((replace("$(field)", "χ" => "")),0)
		m=getfield(mod, field)
		m0=getfield(mod, f0)
		for i in eachindex(m)
			m[i] += randn() * randn_perc * 1e-2
		end

	end
	return mod
end

"""
Apply smoothing to `Seismic` using a Gaussian filter of zwidth and xwidth

# Arguments

* `mod::Seismic` : argument that is modified
* `zperc::Real` : smoothing percentage in z-direction
* `xperc::Real=zperc` : smoothing percentage in x-direction

# Keyword Arguments

* `zmin::Real=mod.mgrid[1][1]` :
* `zmax::Real=mod.mgrid[1][end]` :
* `xmin::Real=mod.mgrid[2][1]` :
* `xmax::Real=mod.mgrid[2][end]` :
* `fields` : fields of seismic model that are to be smooth
"""
function Seismic_smooth(mod::Seismic, zperc::Real, xperc::Real=zperc;
		 zmin::Real=mod.mgrid[1][1], zmax::Real=mod.mgrid[1][end],
		 xmin::Real=mod.mgrid[2][1], xmax::Real=mod.mgrid[2][end],
		 fields=[:χvp, :χρ, :χvs]
			)
	xwidth = Float64(xperc) * 0.01 * abs(mod.mgrid[2][end]-mod.mgrid[2][1])
	zwidth = Float64(zperc) * 0.01 * abs(mod.mgrid[1][end]-mod.mgrid[1][1])
	xnwin=Int(div(xwidth,step(mod.mgrid[2])*2.))
	znwin=Int(div(zwidth,step(mod.mgrid[1])*2.))

	izmin = Interpolation.indminn(mod.mgrid[1], Float64(zmin), 1)[1]
	izmax = Interpolation.indminn(mod.mgrid[1], Float64(zmax), 1)[1]
	ixmin = Interpolation.indminn(mod.mgrid[2], Float64(xmin), 1)[1]
	ixmax = Interpolation.indminn(mod.mgrid[2], Float64(xmax), 1)[1]

	modg=deepcopy(mod)
	for (i,iff) in enumerate(fields)
		m=view(getfield(mod, iff),izmin:izmax,ixmin:ixmax)
		mg=view(getfield(modg, iff),izmin:izmax,ixmin:ixmax)
		imfilter!(mg, m, Kernel.gaussian([znwin,xnwin]));
	end
	return modg
end

"""
Return a truncated seismic model using input bounds.
Note that there is no interpolation going on here, but only truncation, so
the input bounds cannot be strictly imposed.

# Arguments
* `mod::Seismic` : model that is truncated

# Keyword Arguments

* `zmin::Float64=mod.mgrid[1][1]` :
* `zmax::Float64=mod.mgrid[1][end]` :
* `xmin::Float64=mod.mgrid[2][1]` :
* `xmax::Float64=mod.mgrid[2][end]` :
"""
function Seismic_trun(mod::Seismic;
			 zmin::Float64=mod.mgrid[1][1], zmax::Float64=mod.mgrid[1][end],
			 xmin::Float64=mod.mgrid[2][1], xmax::Float64=mod.mgrid[2][end],
			 )

	izmin = Interpolation.indminn(mod.mgrid[1], zmin, 1)[1]
	izmax = Interpolation.indminn(mod.mgrid[1], zmax, 1)[1]
	ixmin = Interpolation.indminn(mod.mgrid[2], xmin, 1)[1]
	ixmax = Interpolation.indminn(mod.mgrid[2], xmax, 1)[1]

	x = mod.mgrid[2][ixmin:ixmax]
	z = mod.mgrid[1][izmin:izmax]

	# allocate model
	mod_trun = Seismic_zeros([z,x])
	adjust_bounds!(mod_trun, mod)

	for f in [:χvp, :χvs, :χρ]
		f0 = Symbol((replace("$(f)", "χ" => "")),0)
		setfield!(mod_trun, f, getfield(mod, f)[izmin:izmax, ixmin:ixmax])
		setfield!(mod_trun, f0, getfield(mod, f0)) # adjust bounds later after truncation if necessary
	end
	return mod_trun

end

"""
Extend a seismic model into PML layers
"""
function Seismic_pml_pad_trun(mod::Seismic, nlayer_rand, npml)

	mg=mod.mgrid
	mgex=[
        range(mg[1][1] - npml*step(mg[1]),
	     stop=mg[1][end] + npml*step(mg[1]), length=length(mg[1])+2*npml),
        range(mg[2][1] - npml*step(mg[2]),
	     stop=mg[2][end] + npml*step(mg[2]), length=length(mg[2])+2*npml)
	]

	modex=Seismic_zeros(mgex)
	adjust_bounds!(modex,mod)
	Seismic_pml_pad_trun!(modex, mod, nlayer_rand)

	return modex
end

"only padding implemented"
function Seismic_pml_pad_trun!(modex::Seismic, mod::Seismic, nlayer_rand)
	vp=mod.χvp; χ!(vp,mod.ref.vp,-1)
	vs=mod.χvs; χ!(vs,mod.ref.vs,-1)
	ρ=mod.χρ; χ!(ρ,mod.ref.ρ,-1)

	vpex=modex.χvp; vsex=modex.χvs; ρex=modex.χρ;
	pml_pad_trun!(vpex,vp,mod.vp0,nlayer_rand,50.0);
	pml_pad_trun!(vsex,vs,[1,2],nlayer_rand,0.0);
	pml_pad_trun!(ρex,ρ,mod.ρ0,nlayer_rand,0.0)
	χ!(vpex,modex.ref.vp,1); χ!(vsex,modex.ref.vs,1); χ!(ρex,modex.ref.ρ,1)
	χ!(vp,mod.ref.vp,1); χ!(vs,mod.ref.vs,1); χ!(ρ,mod.ref.ρ,1)
end

function pml_pad_trun(mod::Array{Float64,2}, np::Int64, flag::Int64=1)
	if(isequal(flag,1))
		nz = size(mod,1); nx = size(mod,2)
		modo = zeros(nz + 2*np, nx + 2*np)
		pml_pad_trun!(modo,mod,1)
	elseif(isequal(flag,-1))
		nz = size(mod,1); nx = size(mod,2)
		modo = zeros(nz - 2*np, nx - 2*np)
		pml_pad_trun!(mod,modo,-1)
	end
	return modo
end

"""
PML Extend a model on all four sides
"""
function pml_pad_trun!(modex::Array{Float64,2}, mod::Array{Float64,2}, bounds, nlayer_rand, var_fact=1.)
	np=(size(modex,1)-size(mod,1) == size(modex,2)-size(mod,2)) ?
			div(size(modex,2)-size(mod,2),2) : error("modex size")
	nz = size(mod,1); nx = size(mod,2)
	#=
	# first fill with random numbers
	Random.seed!(1)
	#for i in eachindex(modex)
	#	modex[i]=rand(Distributions.Uniform(bounds...))
	#end
	# experiments with random boundary conditions
	sample = v->rand(Distributions.Truncated(Distributions.Normal(mean(bounds),v),bounds[1],bounds[2]))
	for ix in 1+np:np+nx
		for iz in 1:np
			v=var_fact*(np-iz)/np+1e-5
			modex[iz,ix] = sample(v)
		end
		for iz in nz+1+np:size(modex,1)
			v=var_fact*(iz-np-nz-1)/np+1e-5
			modex[iz,ix] = sample(v)
		end
	end
	for iz in 1+np:np+nz
		for ix in 1:np
			v=var_fact*(np-ix)/np+1e-5
			modex[iz,ix] = sample(v)
		end
		for ix in nx+1+np:size(modex,2)
			v=var_fact*(ix-np-nz-1)/np+1e-5
			modex[iz,ix] = sample(v)
		end
	end

	for ix in 1:np
		for iz in 1:np
			v=var_fact*sqrt((np-iz)^2+(np-ix)^2)/np+1e-5
			modex[iz,ix] = sample(v)
		end
		for iz in nz+1+np:size(modex,1)
			v=var_fact*sqrt((iz-np-nz-1)^2+(np-ix)^2)/np+1e-5
			modex[iz,ix] = sample(v)
		end
	end
	for ix in np+nx+1:size(modex,1)
		for iz in 1:np
			v=var_fact*sqrt((np-iz)^2+(ix-np-nx-1)^2)/np+1e-5
			modex[iz,ix] = sample(v)
		end
		for iz in nz+1+np:size(modex,1)
			v=var_fact*sqrt((iz-np-nz-1)^2+(ix-np-nx-1)^2)/np+1e-5
			modex[iz,ix] = sample(v)
		end
	end

	=#
	for ix in 1:nx
		for iz in 1:nz
			modex[np+iz,np+ix]=mod[iz,ix]
		end
	end

	# flood borders
	for ix in 1:size(modex,2)
		for iz in 1:np-nlayer_rand
			modex[iz,ix] = modex[np+1-nlayer_rand,ix]
		end
		for iz in nz+1+np+nlayer_rand:size(modex,1)
			modex[iz,ix] = modex[nz+np+nlayer_rand,ix];
		end
	end
	for iz in 1:size(modex,1)
		for ix in 1:np-nlayer_rand
			modex[iz,ix] = modex[iz,np+1-nlayer_rand]
		end
		for ix in nx+1+np+nlayer_rand:size(modex,2)
			modex[iz,ix] = modex[iz,nx+np+nlayer_rand]
		end
	end
	return modex
end

"""
function to resample in the model domain

# Arguments
* `mod::Seismic` : model
* `modi::Seismic` : model after interpolation
"""
function interp_spray!(mod::Seismic, modi::Seismic, attrib::Symbol, Battrib::Symbol=:B2; pa=nothing)
	if(pa===nothing)
		pa=Interpolation.Kernel([mod.mgrid[2], mod.mgrid[1]], [modi.mgrid[2], modi.mgrid[1]], Battrib)
	end

	"loop over fields in `Seismic`"
	Interpolation.interp_spray!(mod.χvp, modi.χvp, pa, attrib)
	Interpolation.interp_spray!(mod.χvs, modi.χvs, pa, attrib)
	Interpolation.interp_spray!(mod.χρ, modi.χρ,  pa, attrib)

end

#
#macro inter

function save(mod::Seismic, folder; N=100)
	!(isdir(folder)) && error("invalid directory")
	error("need to be updated")

	nx=length(mod.mgrid[2])
	nz=length(mod.mgrid[1])
	n=max(nx,nz);
	fact=(n>N) ? round(Int,n/N) : 1
	#mgrid=resamp(mod.mgrid, step(mod.mgrid[2])*fact, step(mod.mgrid[1])*fact)
	x=mgrid[2]
	z=mgrid[1]
	nx=length(x)
	nz=length(z)
	modo=Seismic_zeros(mgrid)
	adjust_bounds!(modo, mod)
	interp_spray!(mod, modo, :interp)

	for m in [:vp, :ρ, :Zp]
		# save original gf
		file=joinpath(folder, string("im", m, ".csv"))
		CSV.write(file,DataFrame(hcat(repeat(z,outer=nx),
					      repeat(x,inner=nz),vec(Seismic_get(modo, m)))))
	end
end

end # module


text:
using Test
using SafeTestsets

@safetestset "Reference shapes tests" begin include("hyperrectangle_test.jl") end

@safetestset "HyperRectangle tests" begin include("hyperrectangle_test.jl") end

@safetestset "Domain tests" begin include("domain_test.jl") end

@safetestset "Element tests" begin include("parametricentity_test.jl") end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# -*- coding: utf-8 -*-
using Catlab
using MacroTools
import MacroTools: postwalk, striplines
using ModelingToolkit
import ModelingToolkit: Constant
using Test
using Petri
using SemanticModels.ModelTools.PetriModels

using Catlab.WiringDiagrams
using Catlab.Doctrines
import Catlab.Doctrines.⊗
import Catlab.Graphics: to_graphviz, to_tikz
import Catlab.Graphics.Graphviz: run_graphviz
⊗(a::WiringDiagram, b::WiringDiagram) = otimes(a,b)
import Base: ∘
∘(a::WiringDiagram, b::WiringDiagram) = compose(b, a)
⊚(a,b) = b ∘ a

# Generators
S, E, I, R, D, G= Ob(FreeSymmetricMonoidalCategory, :S, :E, :I, :R, :D, :G)

# +
infecting = Hom(:infection, S ⊗ I, I⊗I)

ss = mcopy(Ports([S]), 2)

function build(dict, y)
    push!(dict, codom(y)=>compose(typeof(dom(y)) == FreeSymmetricMonoidalCategory.Ob{:generator} ? dict[dom(y)] : otimes(map(x->dict[x], dom(y).args)), y))
end

dict = Dict()
push!(dict, I=>id(I))

rain = Hom(Symbol("rain"), I, R)
sprinkler = Hom(Symbol("sprinkler"), R, S)
grasswet = Hom(Symbol("grasswet"), S⊗R, G)
whatever = Hom(Symbol("whatever"), G⊗S, D)
set = [rain, sprinkler, grasswet, whatever]

map((x)->build(dict, x), set)

to_graphviz(dict[D], labels=true)

#to_graphviz(WiringDiagram()⊚(ss⊗WiringDiagram(Hom(:id, I, I))), labels=true)
#=
s1 = WiringDiagram(Hom(:id, S, S))

s_split = s1 ⊚ ss

to_graphviz(s_split, labels=true)
=#

# +
inf  = WiringDiagram(infecting)
expo = WiringDiagram(Hom(:exposure, S⊗I, E⊗I))
rec  = WiringDiagram(Hom(:recovery, I,   R))
wan  = WiringDiagram(Hom(:waning,   R,   S))

rain = Hom(:rain, I, R)
sprinkler = Hom(:sprinkler, R, S)
grasswet = Hom(Symbol("grasswet"), S⊗R, G)

set = [rain, sprinkler, grasswet]

foldl((x, y)->println(y), set)
# -

si    = WiringDiagram(Hom(:infection,   S⊗I, I⊗I))
se    = WiringDiagram(Hom(:exposure,    S⊗I, E⊗I))
prog  = WiringDiagram(Hom(:progression, E,   I))
fatal = WiringDiagram(Hom(:die,  I, D))
rip   = WiringDiagram(Hom(:rest, D, D))

sir    = si    ⊚ (rec   ⊗ rec)
seir   = se    ⊚ (prog  ⊗ rec)
seirs  = seir  ⊚ (wan   ⊗ wan)
seird  = seir  ⊚ (fatal ⊗ WiringDiagram(Hom(:id, R, R)))
seirds = seird ⊚ (rip   ⊗ wan)

model(PetriModel, sir)

models = [sir, seir, seirs, seird, seirds]
nets   = model.(PetriModel, models)

modelnames = ["sir",
         "seir",
         "seirs",
         "seird",
         "seirds",
         ]
map(zip(nets, modelnames)) do (net, name)
    Δ = reverse(net.model.Δ)
    println("Model: $name\n  Connections:")
    println("    $Δ\n")
    Δ
end

function writesvg(f::Union{IO,AbstractString}, d::WiringDiagram)
    write(f, to_graphviz(d, labels=true)
          |>g->run_graphviz(g, format="svg")
          )
end

writesvg(x) = writesvg(x...)

try
    mkpath("img")
    writesvg.(zip(["img/$f.svg"
                   for f in modelnames],
                  models))
catch e
    error("Could not write images, does ./img exist?\n $e")
end
println("Model Drawings are in ./img")


text:
using Random
using GPR

include("utils/utils.jl")
include("utils/predictdynamics.jl")
include("utils/data/datasets.jl")
include("utils/data/simulations.jl")
include("utils/data/transformations.jl")

include("parallel/core.jl")
include("parallel/dataframes.jl")

include("maximal_coordinates/P1param.jl")
include("maximal_coordinates/P2param.jl")
include("maximal_coordinates/CPparam.jl")
include("maximal_coordinates/FBparam.jl")
include("minimal_coordinates/P1param.jl")
include("minimal_coordinates/P2param.jl")
include("minimal_coordinates/CPparam.jl")
include("minimal_coordinates/FBparam.jl")

function expand_config(config::Dict, EXPERIMENT_ID::String, nsamples::Integer, dfs::Tuple{DataFrame, DataFrame})
    EXPERIMENT_ID *= string(nsamples)
    nprocessed, kstep_mse, params = 0, [], []
    config = Dict("EXPERIMENT_ID"=>EXPERIMENT_ID,
                  "nruns"=>config["nruns"],
                  "Δtsim"=>config["Δtsim"],
                  "datasets"=>dfs,
                  "trainsamples"=>nsamples,
                  "testsamples"=>config["testsamples"],
                  "simsteps"=>config["simsteps"],
                  "nprocessed"=>nprocessed,
                  "kstep_mse"=>kstep_mse,
                  "params"=>params,  # Optimal hyperparameters in sim. Vector of tested parameters in search
                  "mechanismlock"=>ReentrantLock(),
                  "paramlock"=>ReentrantLock(),
                  "checkpointlock"=>ReentrantLock(),
                  "resultlock"=>ReentrantLock())
    return config
end

config = Dict("nruns" => 100, "Δtsim" => 0.0001, "testsamples" => 100, "simsteps" => 20)

dfkeys = ["FB"]  # "P1", "P2", "CP", "FB"
@info "Main Thread: loading datasets $dfkeys"
datasets = Dict(key => loaddatasets(key) for key in dfkeys)
@info "Main Thead: finished loading datasets"

for nsamples in [1024, 2048]  # [2, 4, 8, 16, 32, 64, 128, 256, 512]
    #parallelsearch(experimentP1Max, expand_config(config, "P1_MAX", nsamples, datasets["P1"]))
    #parallelsearch(experimentP2Max, expand_config(config, "P2_MAX", nsamples, datasets["P2"]))
    #parallelsearch(experimentCPMax, expand_config(config, "CP_MAX", nsamples, datasets["CP"]))
    parallelsearch(experimentFBMax, expand_config(config, "FB_MAX", nsamples, datasets["FB"]))

    #parallelsearch(experimentP1Min, expand_config(config, "P1_MIN", nsamples, datasets["P1"]))
    #parallelsearch(experimentP2Min, expand_config(config, "P2_MIN", nsamples, datasets["P2"]))
    #parallelsearch(experimentCPMin, expand_config(config, "CP_MIN", nsamples, datasets["CP"]))
    parallelsearch(experimentFBMin, expand_config(config, "FB_MIN", nsamples, datasets["FB"]))
end


text:
using CellMLModelRepository
using CSV, DataFrames
using Test

df = cellml_metadata()
@test df isa DataFrame

# test suite stuff
cellml_models()
fns = readdir(joinpath(@__DIR__, "../data/cellml_models/"); join=true)[1:10] # only do 10
@show fns
@test isfile(fns[1])

using Pkg

Pkg.add(url="https://github.com/anandijain/BioXMLSciMLTest.jl")
using BioXMLSciMLTest

files_to_sciml(fns; pmap=true)
df = results_df()
@test df isa DataFrame


text:
abstract AbstractInvariant{T<:Number,D<:AbstractDomain}

immutable NullInvariant <: AbstractBank
end

immutable SumInvariant{T<:Number,D<:AbstractDomain} <: AbstractInvariant{T,D}
    domain::D
    pathkey::PathKey
    signaltype::Type{T}
end
SumInvariant(Ux::ScatteredBlob) =
    SumInvariant(Ux.nodes[collect(keys(Ux.nodes))[1]])

function SumInvariant{T<:Number}(node::AbstractNode{T})
    pathkey = node.ranges[1].first
    K = sum([r.first == pathkey for r in node.ranges])
    domain = SpatialDomain(K)
    signaltype = T
    SumInvariant(domain, pathkey, signaltype)
end

function SumInvariant{T<:Number}(node::AbstractFourierNode{T})
    pathkey = node.ranges[1].first
    K = sum([r.first == pathkey for r in node.ranges])
    domain = FourierDomain(K)
    signaltype = T
    SumInvariant(domain, pathkey, signaltype)
end

function (::Type{SumInvariant{T,SpatialDomain{K}}}){T,K,N}(
        nodes::DataStructures.SortedDict{Path,
            Node{T,N},Base.Order.ForwardOrdering})
    sumnodes =
        DataStructures.SortedDict{Path,Node{T,N},Base.Order.ForwardOrdering}()
    for (path, nodevalue) in nodes
        dims = find([r.first == invariant.pathkey for r in nodevalue.ranges])
        sumnodes[path] = Node{T,N}(sum(nodevalue.data, dims), nodevalue.ranges)
    end
    return sumnodes
end

(invariant::SumInvariant)(Ux::ScatteredBlob) = invariant(Ux.nodes)

function transform!{T<:Number}(
        destination::SubArray,
        invariant::SumInvariant{T,SpatialDomain},
        node::AbstractNode{T},
        dim::Int)
    destination[:] = sum(node.data, dim)
end


text:
@testset "PCM16 WAV file reading" begin
    buf = load_wav(reference_wav)
    @test samplerate(buf) == srate
    @test nchannels(buf) == 2
    @test nframes(buf) == 100
    @test isapprox(domain(buf), collect(0:99)/(srate))
    @test mse(buf, reference_buf) < 1e-10
    @test eltype(buf) == PCM16Sample
end

@testset "PCM32 WAV file reading" begin
    buf = load_wav(reference_wav_pcm24)
    @test samplerate(buf) == srate
    @test nchannels(buf) == 2
    @test nframes(buf) == 100
    @test isapprox(domain(buf), collect(0:99)/(srate))
    @test mse(buf, reference_buf) < 1e-10
    @test eltype(buf) == PCM32Sample
end

@testset "Float32 WAV file reading" begin
    buf = load_wav(reference_wav_float)
    @test samplerate(buf) == srate
    @test nchannels(buf) == 2
    @test nframes(buf) == 100
    @test isapprox(domain(buf), collect(0:99)/(srate))
    @test mse(buf, reference_buf) < 1e-10
    @test eltype(buf) == Float32
end

@testset "Float64 WAV file reading" begin
    buf = load_wav(reference_wav_double)
    @test samplerate(buf) == srate
    @test nchannels(buf) == 2
    @test nframes(buf) == 100
    @test isapprox(domain(buf), collect(0:99)/(srate))
    @test mse(buf, reference_buf) < 1e-10
    @test eltype(buf) == Float64
end

@testset "FLAC file reading" begin
    buf = load_flac(reference_flac)
    @test samplerate(buf) == srate
    @test nchannels(buf) == 2
    @test nframes(buf) == 100
    @test isapprox(domain(buf), collect(0:99)/srate)
    @test mse(buf, reference_buf) < 1e-10
    @test eltype(buf) == PCM16Sample
end

@testset "OGG file reading" begin
    buf = load_ogg(reference_ogg)
    @test samplerate(buf) == srate
    @test nchannels(buf) == 2
    @test nframes(buf) == 100
    @test isapprox(domain(buf), collect(0:99)/srate)
    # lossy compression, so relax the accuracy a bit
    @test mse(buf, reference_buf) < 1e-5
end

@testset "Read errors" begin
    @test_throws ErrorException load_wav("doesnotexist.wav")
end


text:
struct FunctionExpr
    fn::Symbol
    fargs::Array{T,1} where T
    output::Symbol
    whstmt::Array{T,1} where T
    argnames::Tuple{Vararg{Symbol}}
    template::Union{Expr, Symbol}
end

function f_parser(Ex::Union{Expr,Symbol}; head=:call)::Expr

    if typeof(Ex) != Expr
        return Expr(:(::), :(nop(nothing)), :Any)
    elseif Ex.head == head
        return Ex
    else
        return f_parser(Ex.args[1]; head=head)
    end

end

function f_header(e::FunctionExpr)::Expr

    if (e.output == :Any) && isempty(e.whstmt)
        let call = :($(e.fn)($(e.fargs...)))
            :(function $call end)
        end
    elseif isempty(e.whstmt)
        let call = :($(e.fn)($(e.fargs...))::$(e.output))
            :(function $call end)
        end
    else
        let call = :($(e.fn)($(e.fargs...))::$(e.output)), cpwhere = copy(e.whstmt)
            while !isempty(cpwhere)
                annotation = pop!(cpwhere)
                call = :($call where $annotation)
            end
            :(function $call end)
        end
    end

end

function f_closure(e::FunctionExpr)::Expr

    let closure = f_header(e)

        push!(closure.args, e.template)
        closure

    end

end

function f_expr(f::Expr)::Expr

    whstmt = let node = f.args[1], annotations = []
        while node.head==:where
            typeinfo = f_parser(node; head=:where).args[2]
            push!(annotations, typeinfo)
            node = node.args[1]
        end
        annotations
    end

    type_annotations = if isempty(whstmt)
        Dict()
    else
        Dict(ex.args[1] => (ex.head, ex.args[2]) for ex in whstmt)
    end

    fn, fargs = let root = f_parser(f.args[1]; head=:call)
        root.args[1], root.args[2:end]
    end

    argnames = let nameonly = map(fargs) do arg
            typeof(arg) == Expr ? arg.args[1] : arg
        end
        tuple(nameonly...)
    end

    OutType = f_parser(f.args[1]; head=:(::)).args[2]

    InTypes = map(fargs) do arg
        let name = typeof(arg)==Symbol ? :Any : f_parser(arg, head=:(::)).args[2]
            if isempty(type_annotations) || !haskey(type_annotations, name)
                name
            else
                op, tp = get(type_annotations, name, (:(<:), :Any))
                Expr(op, tp)
            end
        end
    end

   template = let block =
        let root = f.args[2]
            if typeof(root.args[1])!=Expr
                root.args[2]
            else
                root.args[1]
            end
        end

        let f_block = block.head == :return ? block.args[1] : block
            :(
                let args = tuple($(argnames...))
                    if haskey(cache, args)
                        cache[args]
                    else
                        get!(
                            cache, args,
                            $f_block
                        )
                    end
                end
            )
        end

    end

    let fexpr = FunctionExpr(fn, fargs, OutType, copy(whstmt), argnames, template)
        let f_decstmt = f_header(fexpr), f_block = f_closure(fexpr)
            push!(
                f_decstmt.args,
                :(

                    let cache = Dict{Tuple{$(InTypes...)}, $OutType}()

                        let $fn = $f_block

                            $fn($(argnames...))

                        end

                    end
                )
            )

            return f_decstmt
        end
    end

end


text:
#==
Functions related to identifying resonant modes (eigenangles) within the earth-ionosphere
waveguide.

These functions use Budden's model, which represents the ground and ionosphere as sharply
reflecting boundaries. The middle of the waveguide is filled with a fictitious medium with a
refractive index that mimicks the propagation of the radio wave through free space over
curved earth.
==#

@doc raw"""
    PhysicalModeEquation{W<:HomogeneousWaveguide} <: ModeEquation

Parameters for solving the physical mode equation ``\det(Rg*R - I)``.

Fields:

    - ea::EigenAngle
    - frequency::Frequency
    - waveguide::W
"""
struct PhysicalModeEquation{W<:HomogeneousWaveguide} <: ModeEquation
    ea::EigenAngle
    frequency::Frequency
    waveguide::W
end

"""
    PhysicalModeEquation(f::Frequency, w::HomogeneousWaveguide)

Create a `PhysicalModeEquation` struct with `ea = complex(0.0)`.

See also: [`setea`](@ref)
"""
PhysicalModeEquation(f::Frequency, w::HomogeneousWaveguide) =
    PhysicalModeEquation(EigenAngle(complex(0.0)), f, w)

"""
    setea(ea, modeequation)

Return `modeequation` with eigenangle `ea`.

`ea` will be converted to an `EigenAngle` if necessary.
"""
setea(ea, modeequation::PhysicalModeEquation) =
    PhysicalModeEquation(EigenAngle(ea), modeequation.frequency, modeequation.waveguide)

##########
# Reflection coefficients
##########

"""
    wmatrix(ea::EigenAngle, T)

Compute the four submatrix elements of `W` used in the equation ``dR/dz`` from the
ionosphere with `T` matrix returned as a tuple `(W₁₁, W₂₁, W₁₂, W₂₂)`.

Following Budden's formalism for the reflection matrix of a plane wave obliquely incident on
the ionosphere [Budden1955a], the wave below the ionosphere can be resolved into
upgoing and downgoing waves of elliptical polarization, each of whose components are
themselves resolved into a component with the electric field in the plane of propagation and
a component perpendicular to the plane of propagation. The total field can be written in
matrix form as ``e = Lf`` where ``L`` is a 4×4 matrix that simply selects and specifies the
incident angle of the components and ``f`` is a column matrix of the complex amplitudes of
the component waves. By inversion, ``f = L⁻¹e`` and its derivative with respect to height
``z`` is ``f′ = -iL⁻¹TLf = -½iWf``. Then ``W = 2L⁻¹TL`` describes the change in amplitude of
the upgoing and downgoing component waves.

``W`` is also known as ``S`` in many texts.

# References

[Budden1955a]: K. G. Budden, “The numerical solution of differential equations governing
    reflexion of long radio waves from the ionosphere,” Proc. R. Soc. Lond. A, vol. 227,
    no. 1171, pp. 516–537, Feb. 1955.
"""
function wmatrix(ea::EigenAngle, T)
    C, Cinv = ea.cosθ, ea.secθ

    # Precompute
    T12Cinv = T[1,2]*Cinv
    T14Cinv = T[1,4]*Cinv
    T32Cinv = T[3,2]*Cinv
    T34Cinv = T[3,4]*Cinv
    CT41 = C*T[4,1]

    #==
    T = SMatrix{4,4}(T11, 0, T31, T41,
                     T12, 0, T32, T42,
                     0,   1,   0,   0,
                     T14, 0, T34, T44)
    L = SMatrix{4,4}(C,  0,  0, 1,
                     0, -1, -C, 0,
                    -C,  0,  0, 1,
                     0, -1,  C, 0)
    # 2*inv(L) = Linv2
    L2inv = SMatrix{4,4}(Cinv,  0, -Cinv,    0,
                         0,    -1,     0,   -1,
                         0, -Cinv,     0, Cinv,
                         1,     0,     1,    0)
    W = Linv2*T*L

    # W = 2*(L\T)*L

    ---

    W = | W11 | W12 |
        | W21 | W22 |

    W11 = | a11+a11r | -b11 |
          | -c11     | d11  |

    W12 = | a21+a21r | -b11 |
          | c12      | d12  |

    W21 = | a21-a21r | b21  |
          | c11      | -d12 |

    W22 = | a11-a11r | b21  |
          | -c12     | -d11 |
    ==#

    a11 = T[1,1] + T[4,4]
    a11r = T14Cinv + CT41
    a21 = T[4,4] - T[1,1]
    a21r = T14Cinv - CT41

    b11 = T12Cinv + T[4,2]
    b21 = T12Cinv - T[4,2]

    c11 = T[3,1] + T34Cinv
    c12 = T[3,1] - T34Cinv

    d11 = C + T32Cinv
    d12 = T32Cinv - C

    # Form the four 2x2 submatrices of `S`
    W11 = SMatrix{2,2}(a11+a11r, -c11, -b11, d11)
    W12 = SMatrix{2,2}(a21+a21r, c12, -b11, d12)
    W21 = SMatrix{2,2}(a21-a21r, c11, b21, -d12)
    W22 = SMatrix{2,2}(a11-a11r, -c12, b21, -d11)

    return W11, W21, W12, W22
end

"""
    dwmatrix(ea::EigenAngle, T, dT)

Compute the four submatrix elements of ``dW/dθ`` returned as the tuple
`(dW₁₁, dW₂₁, dW₁₂, dW₂₂)` from the ionosphere with `T` matrix and its derivative with
respect to ``θ``, `dT`.
"""
function dwmatrix(ea::EigenAngle, T, dT)
    C, S, C², Cinv = ea.cosθ, ea.sinθ, ea.cos²θ, ea.secθ
    C²inv = Cinv^2

    dC = -S
    dCinv = S*C²inv

    dt12Cinv = dT[1,2]*Cinv + T[1,2]*dCinv
    dt14Cinv = dT[1,4]*Cinv + T[1,4]*dCinv
    dt32Cinv = dT[3,2]*Cinv + T[3,2]*dCinv
    dt34Cinv = dT[3,4]*Cinv + T[3,4]*dCinv
    dt41C = dC*T[4,1]

    ds11a = dT[1,1] + dT[4,4]
    dd11a = dT[1,1] - dT[4,4]
    ds11b = dt14Cinv + dt41C
    dd11b = dt14Cinv - dt41C
    ds12 = dt12Cinv
    dd12 = dt12Cinv
    ds21 = dt34Cinv
    dd21 = -dt34Cinv
    ds22 = dC + dt32Cinv
    dd22 = dC - dt32Cinv

    # Form the four 2x2 submatrices of `dW`
    dW11 = SMatrix{2,2}(ds11a+ds11b, -ds21, -ds12, ds22)
    dW12 = SMatrix{2,2}(-dd11a+dd11b, dd21, -ds12, -dd22)
    dW21 = SMatrix{2,2}(-dd11a-dd11b, ds21, dd12, dd22)
    dW22 = SMatrix{2,2}(ds11a-ds11b, -dd21, dd12, -ds22)

    return dW11, dW21, dW12, dW22
end

"""
    dRdz(R, modeequation, z, susceptibilityfcn=z->susceptibility(z, modeequation; params=LMPParams()))

Compute the differential of the reflection matrix `R`, ``dR/dz``, at height `z`.
`susceptibilityfcn` is a function returning the ionosphere susceptibility at height `z`.

Following the Budden formalism for the reflection of an (obliquely) incident plane wave from
a horizontally stratified ionosphere [Budden1955a], the differential of the
reflection matrix `R` with height `z` can be described by
```math
dR/dz = k/(2i)⋅(W₂₁ + W₂₂R - RW₁₁ - RW₁₂R)
```
Integrating ``dR/dz`` downwards through the ionosphere gives the reflection matrix ``R`` for
the ionosphere as if it were a sharp boundary at the stopping level with free space below.

# References

[Budden1955a]: K. G. Budden, “The numerical solution of differential equations governing
    reflexion of long radio waves from the ionosphere,” Proc. R. Soc. Lond. A, vol. 227,
    no. 1171, pp. 516–537, Feb. 1955.
"""
function dRdz(R, modeequation, z, susceptibilityfcn=z->susceptibility(z, modeequation; params=LMPParams()))
    @unpack ea, frequency = modeequation

    k = frequency.k

    M = susceptibilityfcn(z)
    T = tmatrix(ea, M)
    W11, W21, W12, W22 = wmatrix(ea, T)

    # the factor k/(2i) isn't explicitly in [Budden1955a] because of his change of variable
    # ``s = kz``
    return k/2im*(W21 + W22*R - R*W11 - R*W12*R)
end

"""
    dRdθdz(RdRdθ, p, z)

Compute the differential ``dR/dθ/dz`` at height `z` returned as an `SMatrix{4,2}` with
``dR/dz`` in the first 2 rows and ``dR/dθ/dz`` in the bottom 2 rows.

`p` is a tuple containing instances `(PhysicalModeEquation(), LMPParams())`.
"""
function dRdθdz(RdRdθ, p, z)
    modeequation, params = p
    @unpack ea, frequency = modeequation

    k = frequency.k

    M = susceptibility(z, modeequation; params=params)
    T = tmatrix(ea, M)
    dT = dtmatrix(ea, M)
    W11, W21, W12, W22 = wmatrix(ea, T)
    dW11, dW21, dW12, dW22 = dwmatrix(ea, T, dT)

    R = RdRdθ[SVector(1,2),:]
    dRdθ = RdRdθ[SVector(3,4),:]

    dz = k/2im*(W21 + W22*R - R*W11 - R*W12*R)
    dθdz = k/2im*(dW21 + dW22*R + W22*dRdθ - (dRdθ*W11 + R*dW11) -
                    (dRdθ*W12*R + R*dW12*R + R*W12*dRdθ))

    return vcat(dz, dθdz)
end

"""
    integratedreflection(modeequation::PhysicalModeEquation;
        params=LMPParams(), susceptibilityfcn=z->susceptibility(z, modeequation; params=params))

Integrate ``dR/dz`` downward through the ionosphere described by `modeequation` from
`params.topheight`, returning the ionosphere reflection coefficient `R` at the ground.
`susceptibilityfcn` is a function returning the ionosphere susceptibility tensor as a
function of altitude `z` in meters.

`params.integrationparams` are passed to `DifferentialEquations.jl`.
"""
function integratedreflection(modeequation::PhysicalModeEquation;
    params=LMPParams(), susceptibilityfcn=z->susceptibility(z, modeequation; params=params))

    @unpack topheight, integrationparams = params
    @unpack tolerance, solver, dt, force_dtmin, maxiters = integrationparams

    Mtop = susceptibility(topheight, modeequation; params=params)
    Rtop = bookerreflection(modeequation.ea, Mtop)

    prob = ODEProblem{false}((R,p,z)->dRdz(R,p,z,susceptibilityfcn), Rtop, (topheight, BOTTOMHEIGHT), modeequation)

    # WARNING: When save_on=false, don't try interpolating the solution!
    sol = solve(prob, solver; abstol=tolerance, reltol=tolerance,
                force_dtmin=force_dtmin, dt=dt, maxiters=maxiters,
                save_on=false, save_start=false, save_end=true)

    R = sol[end]

    return R
end

"""
    integratedreflection(modeequation::PhysicalModeEquation, ::Dθ; params=LMPParams())

Compute ``R`` and ``dR/dθ`` as an `SMatrix{4,2}` with ``R`` in rows (1, 2) and ``dR/dθ`` in
rows (3, 4).

The `params.integrationparams.tolerance` is hardcoded to `1e-10` in this version of the
function.
"""
function integratedreflection(modeequation::PhysicalModeEquation, ::Dθ; params=LMPParams())
    @unpack topheight, integrationparams = params
    @unpack solver, dt, force_dtmin, maxiters = integrationparams

    # Tolerance is overridden for this `::Dθ` form.
    # Using an identical accuracy appears to result in relatively less accurate solutions
    # compared to the non-Dθ form.
    tolerance = 1e-10

    Mtop = susceptibility(topheight, modeequation; params=params)
    Rtop, dRdθtop = bookerreflection(modeequation.ea, Mtop, Dθ())
    RdRdθtop = vcat(Rtop, dRdθtop)

    prob = ODEProblem{false}(dRdθdz, RdRdθtop, (topheight, BOTTOMHEIGHT),
                             (modeequation, params))

    # WARNING: When save_on=false, don't try interpolating the solution!
    sol = solve(prob, solver; abstol=tolerance, reltol=tolerance,
                force_dtmin=force_dtmin, dt=dt, maxiters=maxiters,
                save_on=false, save_start=false, save_end=true)

    RdR = sol[end]

    return RdR
end

##########
# Ground reflection coefficient matrix
##########

"""
    fresnelreflection(ea::EigenAngle, ground::Ground, frequency::Frequency)
    fresnelreflection(m::PhysicalModeEquation)

Compute the Fresnel reflection coefficient matrix for the ground-freespace interface at the
ground.
"""
fresnelreflection

function fresnelreflection(ea::EigenAngle, ground::Ground, frequency::Frequency)
    C, S² = ea.cosθ, ea.sin²θ
    ω = frequency.ω

    Ng² = complex(ground.ϵᵣ, -ground.σ/(ω*E0))

    CNg² = C*Ng²
    sqrtNg²mS² = sqrt(Ng² - S²)

    Rg11 = (CNg² - sqrtNg²mS²)/(CNg² + sqrtNg²mS²)
    Rg22 = (C - sqrtNg²mS²)/(C + sqrtNg²mS²)

    Rg = SDiagonal(Rg11, Rg22)

    return Rg
end

fresnelreflection(m::PhysicalModeEquation) =
    fresnelreflection(m.ea, m.waveguide.ground, m.frequency)

"""
    fresnelreflection(ea::EigenAngle, ground::Ground, frequency::Frequency, ::Dθ)
    fresnelreflection(m::PhysicalModeEquation, ::Dθ)

Compute the Fresnel reflection coefficient matrix for the ground as well as its derivative
with respect to ``θ`` returned as the tuple `(Rg, dRg)`.
"""
function fresnelreflection(ea::EigenAngle, ground::Ground, frequency::Frequency, ::Dθ)
    C, S, S² = ea.cosθ, ea.sinθ, ea.sin²θ
    S2 = 2*S
    ω = frequency.ω

    Ng² = complex(ground.ϵᵣ, -ground.σ/(ω*E0))

    CNg² = C*Ng²
    sqrtNg²mS² = sqrt(Ng² - S²)

    Rg11 = (CNg² - sqrtNg²mS²)/(CNg² + sqrtNg²mS²)
    Rg22 = (C - sqrtNg²mS²)/(C + sqrtNg²mS²)

    Rg = SDiagonal(Rg11, Rg22)

    dRg11 = (S2*Ng²*(1 - Ng²))/(sqrtNg²mS²*(CNg² + sqrtNg²mS²)^2)
    dRg22 = (S2*(C - sqrtNg²mS²))/(sqrtNg²mS²*(sqrtNg²mS² + C))

    dRg = SDiagonal(dRg11, dRg22)

    return Rg, dRg
end

fresnelreflection(m::PhysicalModeEquation, ::Dθ) =
    fresnelreflection(m.ea, m.waveguide.ground, m.frequency, Dθ())

##########
# Identify EigenAngles
##########

"""
    modalequation(R, Rg)

Compute the determinental mode equation ``det(Rg R - I)`` given reflection coefficients `R`
and `Rg`.

A propagating waveguide mode requires that a wave, having reflected from the ionosphere and
then the ground, must be identical with the original upgoing wave. This criteria is met at
roots of the mode equation [Budden1962].

# References

[Budden1962]: K. G. Budden and N. F. Mott, “The influence of the earth’s magnetic field on
    radio propagation by wave-guide modes,” Proceedings of the Royal Society of London.
    Series A. Mathematical and Physical Sciences, vol. 265, no. 1323, pp. 538–553,
    Feb. 1962.
"""
function modalequation(R, Rg)
    return det(Rg*R - I)
end

"""
    dmodalequation(R, dR, Rg, dRg)

Compute the derivative of the determinantal mode equation with respect to ``θ``.
"""
function dmodalequation(R, dR, Rg, dRg)
    # See e.g. https://folk.ntnu.no/hanche/notes/diffdet/diffdet.pdf
    A = Rg*R - I
    dA = dRg*R + Rg*dR
    return det(A)*tr(A\dA)
end

"""
    solvemodalequation(modeequation::PhysicalModeEquation;
        params=LMPParams(), susceptibilityfcn=z->susceptibility(z, modeequation; params=params))

Compute the ionosphere and ground reflection coefficients and return the value of the
determinental modal equation associated with `modeequation`. `susceptibilityfcn` is a
function that returns the ionosphere susceptibility as a function of altitude `z` in meters.

See also: [`solvedmodalequation`](@ref)
"""
function solvemodalequation(modeequation::PhysicalModeEquation;
    params=LMPParams(), susceptibilityfcn=z->susceptibility(z, modeequation; params=params))

    R = integratedreflection(modeequation; params=params, susceptibilityfcn=susceptibilityfcn)
    Rg = fresnelreflection(modeequation)

    f = modalequation(R, Rg)
    return f
end

"""
    solvemodalequation(θ, modeequation::PhysicalModeEquation;
        params=LMPParams(), susceptibilityfcn=z->susceptibility(z, modeequation; params=params))

Set `θ` for `modeequation` and then solve the modal equation.
"""
function solvemodalequation(θ, modeequation::PhysicalModeEquation;
    params=LMPParams(), susceptibilityfcn=z->susceptibility(z, modeequation; params=params))

    # Convenience function for `grpf`
    modeequation = setea(EigenAngle(θ), modeequation)
    solvemodalequation(modeequation; params=params, susceptibilityfcn=susceptibilityfcn)
end

"""
    solvedmodalequation(modeequation::PhysicalModeEquation; params=LMPParams())

Compute the derivative of the modal equation with respect to ``θ`` returned as the tuple
`(dF, R, Rg)` for the ionosphere and ground reflection coefficients.
"""
function solvedmodalequation(modeequation::PhysicalModeEquation; params=LMPParams())
    RdR = integratedreflection(modeequation, Dθ(); params=params)
    R = RdR[SVector(1,2),:]
    dR = RdR[SVector(3,4),:]

    Rg, dRg = fresnelreflection(modeequation, Dθ())

    dF = dmodalequation(R, dR, Rg, dRg)
    return dF, R, Rg
end

"""
    solvedmodalequation(θ, modeequation::PhysicalModeEquation; params=LMPParams())

Set `θ` for `modeequation` and then solve the derivative of the mode equation with respect
to `θ`.
"""
function solvedmodalequation(θ, modeequation::PhysicalModeEquation; params=LMPParams())
    modeequation = setea(EigenAngle(θ), modeequation)
    solvedmodalequation(modeequation; params=params)
end

"""
    findmodes(modeequation::ModeEquation, mesh=nothing; params=LMPParams())

Find `EigenAngle`s associated with `modeequation.waveguide` within the domain of
`mesh`.

`mesh` should be an array of complex numbers that make up the original grid over which
the GRPF algorithm searches for roots of `modeequation`. If `mesh === nothing`,
it is computed with [`defaultmesh`](@ref).

There is a check for redundant modes that requires modes to be separated by at least
1 orders of magnitude greater than `grpfparams.tolerance` in real and/or imaginary
component. For example, if `grpfparams.tolerance = 1e-5`, then either the real or imaginary
component of each mode must be separated by at least 1e-4 from every other mode.
"""
function findmodes(modeequation::ModeEquation, mesh=nothing; params=LMPParams())
    @unpack approxsusceptibility, grpfparams = params

    if isnothing(mesh)
        mesh = defaultmesh(modeequation.frequency)
    end

    # WARNING: If tolerance of mode finder is much less than the R integration tolerance, it
    # is possible multiple identical modes will be identified. Checks for valid and
    # redundant modes help ensure valid eigenangles are returned from this function.

    if approxsusceptibility
        susceptibilityfcn = susceptibilityspline(modeequation; params=params)
    else
        susceptibilityfcn = z -> susceptibility(z, modeequation; params=params)
    end

    roots, _ = grpf(θ->solvemodalequation(θ, modeequation;
        params=params, susceptibilityfcn=susceptibilityfcn), mesh, grpfparams)

    # Scale tolerance for filtering
    # if tolerance is 1e-8, this rounds to 7 decimal places
    ndigits = round(Int, abs(log10(grpfparams.tolerance)+1), RoundDown)

    # Remove any redundant modes
    sort!(roots; by=reim, rev=true)
    unique!(z->round(z; digits=ndigits), roots)

    return EigenAngle.(roots)
end

#==
Mesh grids for `GRPF`
==#

"""
    defaultmesh(frequency; rmin=deg2rad(30.0), imin=deg2rad(-10.0),
        Δr_coarse=deg2rad(0.5), Δr_fine=deg2rad(0.1),
        rtransition=deg2rad(75.0), itransition=deg2rad(-1.5))

Generate vector of complex coordinates to be used by GRPF in the search for
waveguide modes.

`rmin` is the lower bound of the real axis and `imin` is the lower bound of the imaginary
axis.

At frequencies above 12 kHz the mesh spacing in the upper right corner of the domain
with real values above `rtransition` and imaginary values above `itransition` is
`Δr_fine` and is `Δr_coarse` everywhere else.

At frequencies below 12 kHz the mesh spacing is always `Δr_coarse`.

The lower right diagonal of the lower right quadrant of the complex plane is excluded
from the mesh.

See also: [`findmodes`](@ref)
"""
function defaultmesh(frequency;
    rmin=deg2rad(30.0), imin=deg2rad(-10.0),
    Δr_coarse=deg2rad(0.5), Δr_fine=deg2rad(0.1),
    rtransition=deg2rad(75.0), itransition=deg2rad(-1.5))

    # TODO: get a better idea of frequency transition
    if frequency > 12000
        zbl_coarse = complex(rmin, imin)
        ztr_coarse = complex(deg2rad(89.9), 0.0)

        mesh = trianglemesh(zbl_coarse, ztr_coarse, Δr_coarse)

        filter!(z->(real(z) < rtransition || imag(z) < itransition), mesh)

        zbl_fine = complex(rtransition, itransition)
        ztr_fine = complex(deg2rad(89.9), 0.0)

        append!(mesh, trianglemesh(zbl_fine, ztr_fine, Δr_fine))
    else
        zbl = complex(rmin, imin)
        ztr = complex(deg2rad(89.9), 0.0)

        mesh = trianglemesh(zbl, ztr, Δr_coarse)
    end

    return mesh
end
defaultmesh(f::Frequency; kw...) = defaultmesh(f.f; kw...)

"""
    trianglemesh(zbl, ztr, Δr)

Generate initial mesh node coordinates for a right triangle domain from complex
coordinate `zbl` in the bottom left and `ztr` in the top right with initial mesh step `Δr`.

`zbl` and `ztr` are located on the complex plane at the locations marked in the
diagram below:

```
       im
       |
-re ---|------ ztr
       |     /
       |    /
       |   /
       |  /
       zbl
```
"""
function trianglemesh(zbl, ztr, Δr)
    rzbl, izbl = reim(zbl)
    rztr, iztr = reim(ztr)

    X = rztr - rzbl
    Y = iztr - izbl

    n = ceil(Int, Y/Δr)
    dy = Y/n

    ## dx = sqrt(Δr² - (dy/2)²), solved for equilateral triangle
    m = ceil(Int, X/sqrt(Δr^2 - dy^2/4))
    dx = X/m
    half_dx = dx/2

    slope = 1  # 45° angle

    T = promote_type(ComplexF64, typeof(zbl), typeof(ztr), typeof(Δr))
    mesh = Vector{T}()

    shift = false  # we will displace every other line by dx/2
    for j = 0:n
        y = izbl + dy*j

        for i = 0:m
            x = rzbl + dx*i

            if shift && i == 0
                continue  # otherwise, we shift out of left bound
            elseif shift
                x -= half_dx
            end

            if i == m
                shift = !shift
            end

            ## NEW: check if `x, y` is in upper left triangle
            if y >= slope*x - π/2
                push!(mesh, complex(x, y))
            end
        end
    end

    return mesh
end


text:
using HorizonSideRobots
include("../basics.jl")

md"""
    Посчитать число всех горизонтальных прямолинейных перегородок (вертикальных - нет)
"""

"""
Считает кол-во горизонтальных перегородок в строке, двигаясь в направлении {side}
"""
function count_borders_inline!(r, side)
    counter = 0
    is_under_border = false
    while !isborder(r, side)
        if !is_under_border && isborder(r, Nord)
            counter+=1
            is_under_border = true
        elseif !isborder(r, Nord)
            is_under_border = false
        end
        move!(r, side)
    end
    return counter
end

"""Решение задачи 20"""
function solve_20!(r)
    path = move_to_corner!(r, West, Sud)

    counter = 0
    while !isborder(r, Nord)
        counter += count_borders_inline!(r, Ost)
        carefull_move!(r, Nord)
        if !isborder(r, Nord)
            counter += count_borders_inline!(r, West)
            carefull_move!(r, Nord)
        end
    end

    move_to_corner!(r, West, Sud)
    move_back_from_corner!(r, West, Sud, reverse(path))
    return counter
end


text:
# Most of this file is shamelessly copied from LaTeXStrings.jl. Here's its license:

#=
Copyright © 2013-2014 by Steven G. Johnson.

Permission is hereby granted, free of charge, to any person obtaining a copy of this software and associated documentation files (the "Software"), to deal in the Software without restriction, including without limitation the rights to use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of the Software, and to permit persons to whom the Software is furnished to do so, subject to the following conditions:

The above copyright notice and this permission notice shall be included in all copies or substantial portions of the Software.

THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
=#

@doc raw"""
    Lsvg"..."

Renders the input as SVG with the current preamble and LaTeX engine.

Similar to `raw"..."` strings, you do not have to escape special characters like `\` and `$`, which makes it easier to write LaTeX.

You can interpolate variables using the `%$` syntax, which is the same as the `L"..."` string macro in `LaTeXStrings.jl`. For instance, if there is a variable `x = 10`, then `Lsvg"number %$x"` renders the string `"number 10"`. However, unlike the `L"..."` macro, `Lsvg"..."` does not automatically wrap your input in a pair of dollar signs, since it is very much possible that you'll want to use some other LaTeX environments.
"""
macro Lsvg_str(s::String)
    return esc(_Lsvg(s, String(__source__.file)))
end

function _Lsvg(s::String, filename)
    i = firstindex(s)
    buf = IOBuffer(maxsize=ncodeunits(s))
    ex = Expr(:call, GlobalRef(LatexSVG, :latexsvg))
    while i <= ncodeunits(s)
        c = @inbounds s[i]
        i = nextind(s, i)
        if c === '%' && i <= ncodeunits(s)
            c = @inbounds s[i]
            if c === '$'
                position(buf) > 0 && push!(ex.args, String(take!(buf)))
                atom, i = Meta.parseatom(s, nextind(s, i), filename=filename)
                Meta.isexpr(atom, :incomplete) && error(atom.args[1])
                atom !== nothing && push!(ex.args, atom)
                continue
            else
                print(buf, '%')
            end
        else
            print(buf, c)
        end
    end
    position(buf) > 0 && push!(ex.args, String(take!(buf)))
    return ex
end


text:
abstract AbstractRemoteSegments

remotesegments(p::AbstractPolicy,m::AbstractModel,s::AbstractSampler,nproposals::Int) = _remotesegments(p,m,s,nproposals)

numsegments(s::AbstractRemoteSegments) = s.numsegments
numproposalspersegment(s::AbstractRemoteSegments) = s.numproposalspersegment
numtotalproposals(s::AbstractRemoteSegments) = s.numsegments*s.numproposalspersegment

@inline _numjobsegments(::Type{Val{:procs}}) = nprocs()
@inline _numjobsegments(::Type{Val{:workers}}) = nworkers()
@inline _numjobsegments(::Type{Val{:none}}) = 0
@inline _numjobsegments(::Type{Val{:one}}) = 1
@inline _numjobsegments(::Type{Val{:two}}) = 2

@inline _processnumbers(::Type{Val{:procs}}) = procs()
@inline _processnumbers(::Type{Val{:workers}}) = workers()
@inline _processnumbers(::Type{Val{:none}}) = []
@inline _processnumbers(::Type{Val{:one}}) = workers()
@inline _processnumbers(::Type{Val{:two}}) = workers()

@inline _numproposalspersegment(nproposals::Int,njobsegments::Int) = ceil(Int,nproposals/njobsegments)


text:
filenames = ["f1.txt", "f2.txt"]
for filename in filenames
   txt = read(filename, String)
   open(filename, "w") do f
      write(f, replace(txt, "Goodbye London!" => "Hello New York!"))
   end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Base.Threads

include("update.jl")

function rw_map(xs, nsteps=10000)
    for i in 1:nsteps
        map!(update, xs, xs)
    end
end

function rw_loop(xs, nsteps=10000)
    num = length(xs)
    for i in 1:nsteps
        @inbounds for j in 1:num
            xs[j] += update(xs[j])
        end
    end
end

function rw_thread(xs, nsteps=10000)
    num = length(xs)
    for i in 1:nsteps
        @inbounds @threads for j in 1:num
            xs[j] += update(xs[j])
        end
    end
end

max_stars_repo_name: JuliaTagBot/JuliaBook-Samples

max_issues_repo_name: JuliaTagBot/JuliaBook-Samples

max_forks_repo_name: JuliaTagBot/JuliaBook-Samples


text:
# originates from sisue #458

# using Revise

using Test
using IncrementalInference

@testset "Test clique factors, #458 Example 1" begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addVariable!(fg, :x2, ContinuousScalar)
addVariable!(fg, :x3, ContinuousScalar)
addVariable!(fg, :x4, ContinuousScalar)

addVariable!(fg, :l0, ContinuousScalar)
addVariable!(fg, :l1, ContinuousScalar)

lc = LinearConditional(Normal())
lp = Prior(Normal())
addFactor!(fg, [:x0;:x1], lc, autoinit=false)
addFactor!(fg, [:x1;:x2], lc, autoinit=false)
addFactor!(fg, [:x2;:x3], lc, autoinit=false)
addFactor!(fg, [:x3;:x4], lc, autoinit=false)

addFactor!(fg, [:x0;:l0], lc, autoinit=false)
addFactor!(fg, [:x2;:l0], lc, autoinit=false)

addFactor!(fg, [:x0;:l1], lc, autoinit=false)
addFactor!(fg, [:x2;:l1], lc, autoinit=false)

addFactor!(fg, [:x0;], lp, autoinit=false)
addFactor!(fg, [:l0;], lp, autoinit=false)

# build tree
vo = Symbol[:x2, :x0, :l0, :x3, :x1, :l1, :x4]
tree = resetBuildTreeFromOrder!(fg, vo)
# drawTree(tree, show=true)

# check that frontal variables only show up once
frontals = getCliqFrontalVarIds.((x->getCliq(tree, x)).([:x0; :l0; :x4]))

@test intersect(frontals[1], frontals[2]) |> length == 0
@test intersect(frontals[2], frontals[3]) |> length == 0
@test intersect(frontals[1], frontals[3]) |> length == 0

# check that all variables exist as frontals
lsvars = ls(fg)
@test intersect(union(frontals...), lsvars) |> length == lsvars |> length

## Now check if factors in cliques are okay

C3 = getCliq(tree, :x0)
C3_fg = buildCliqSubgraph(fg, tree, C3 )
# drawGraph(C3_fg, show=true)

C3_fcts = [:x0l0f1;:x0l1f1;:x0x1f1;:x0f1]

@test intersect(ls(C3_fg), [:x0; :x1; :l0; :l1]) |> length == 4
@test intersect(lsf(C3_fg), C3_fcts) |> length == length(C3_fcts)
@test intersect(getCliqFactorIdsAll(C3), C3_fcts) |> length == length(C3_fcts)

C2 = getCliq(tree, :l0)
C2_fg = buildCliqSubgraph(fg, tree, C2)
# drawGraph(C2_fg, show=true)

C2_fcts = [:x1x2f1; :x2x3f1; :x2l0f1; :x2l1f1; :l0f1]

@test intersect(ls(C2_fg), [:x3; :x2; :x1; :l0; :l1]) |> length == 5
@test intersect(lsf(C2_fg), C2_fcts) |> length == length(C2_fcts)
@test intersect(getCliqFactorIdsAll(C2), C2_fcts) |> length == length(C2_fcts)

C1 = getCliq(tree, :x4)
C1_fg = buildCliqSubgraph(fg, tree, C1)
# drawGraph(C1_fg, show=true)

C1_fcts = [:x3x4f1;]

@test intersect(ls(C1_fg), [:x3; :x4; :x1; :l1]) |> length == 4
@test intersect(lsf(C1_fg), C1_fcts) |> length == length(C1_fcts)
@test intersect(getCliqFactorIdsAll(C1), C1_fcts) |> length == length(C1_fcts)

# check that all factors are counted
allCliqFcts = union(C1_fcts, C2_fcts, C3_fcts)

@test length(intersect(lsf(fg), allCliqFcts)) == length(allCliqFcts)

end

@testset "Test clique factors, #458 Example 2" begin

@warn "Test for Example 2 from 458 must still be coded."

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)
addVariable!(fg, :x2, ContinuousScalar)
addVariable!(fg, :x3, ContinuousScalar)

addVariable!(fg, :lm0, ContinuousScalar)
addVariable!(fg, :lm3, ContinuousScalar)

lc = LinearConditional(Normal())
lp = Prior(Normal())
addFactor!(fg, [:x0;:x1], lc, autoinit=false)
addFactor!(fg, [:x1;:x2], lc, autoinit=false)
addFactor!(fg, [:x2;:x3], lc, autoinit=false)

addFactor!(fg, [:x0;:lm0], lc, autoinit=false)
addFactor!(fg, [:x1;:lm0], lc, autoinit=false)
addFactor!(fg, [:x2;:lm3], lc, autoinit=false)
addFactor!(fg, [:x3;:lm3], lc, autoinit=false)

# particular order from 458
vo = Symbol[:x0, :x2, :x1, :lm3, :lm0, :x3]
tree = resetBuildTreeFromOrder!(fg, vo)
# drawTree(tree, show=true)

sfg_1 = buildCliqSubgraph(fg,tree,:x3)
sfg_2 = buildCliqSubgraph(fg,tree,:x2)
sfg_3 = buildCliqSubgraph(fg,tree,:x0)
# drawGraph(sfg_1, show=true)
# drawGraph(sfg_2, show=true)
# drawGraph(sfg_3, show=true)

## WIP
C1_fcts = [:x1lm0f1; :x3lm3f1]
C2_fcts = [:x1x2f1; :x2x3f1; :x2lm3f1]
C3_fcts = [:x0x1f1; :x0lm0f1]

# check all factors are accounted for
@test union(C1_fcts, C2_fcts, C3_fcts) |> length == lsf(fg) |> length
@test intersect(C1_fcts, C2_fcts) |> length == 0
@test intersect(C1_fcts, C3_fcts) |> length == 0
@test intersect(C2_fcts, C3_fcts) |> length == 0

# clique 1
@test intersect(getCliqFactorIdsAll(tree, :x3), C1_fcts) |> length == 2
@test intersect(lsf(sfg_1), C1_fcts) |> length == 2

# clique 2
@test intersect( getCliqFactorIdsAll(tree, :x2), C2_fcts) |> length == 3
@test intersect( lsf(sfg_2), C2_fcts) |> length == 3

# clique 3
@test intersect( getCliqFactorIdsAll(tree, :x0), C3_fcts) |> length == 2
@test intersect( lsf(sfg_3), C3_fcts) |> length == 2

end

#


text:
type Set{T}
    dict::Dict{T,Nothing}

    Set() = new(Dict{T,Nothing}())
    Set(itr) = union!(new(Dict{T,Nothing}()), itr)
end
Set() = Set{Any}()
Set(itr) = Set{eltype(itr)}(itr)

show(io::IO, s::Set) = (show(io, typeof(s)); show_comma_array(io, s,"({","})"))

isempty(s::Set) = isempty(s.dict)
length(s::Set)  = length(s.dict)
eltype{T}(s::Set{T}) = T

in(x, s::Set) = haskey(s.dict, x)

push!(s::Set, x) = (s.dict[x] = nothing; s)
pop!(s::Set, x) = (pop!(s.dict, x); x)
pop!(s::Set, x, deflt) = pop!(s.dict, x, deflt) == deflt ? deflt : x
delete!(s::Set, x) = (delete!(s.dict, x); s)

union!(s::Set, xs) = (for x=xs; push!(s,x); end; s)
setdiff!(s::Set, xs) = (for x=xs; delete!(s,x); end; s)

similar{T}(s::Set{T}) = Set{T}()
copy(s::Set) = union!(similar(s), s)

sizehint(s::Set, newsz) = sizehint(s.dict, newsz)
empty!{T}(s::Set{T}) = (empty!(s.dict); s)

start(s::Set)       = start(s.dict)
done(s::Set, state) = done(s.dict, state)
# NOTE: manually optimized to take advantage of Dict representation
next(s::Set, i)     = (s.dict.keys[i], skip_deleted(s.dict,i+1))

# TODO: simplify me?
pop!(s::Set) = (val = s.dict.keys[start(s.dict)]; delete!(s.dict, val); val)

join_eltype() = None
join_eltype(v1, vs...) = typejoin(eltype(v1), join_eltype(vs...))

union() = Set()
union(s::Set) = copy(s)
function union(s::Set, sets::Set...)
    u = Set{join_eltype(s, sets...)}()
    union!(u,s)
    for t in sets
        union!(u,t)
    end
    return u
end

intersect(s::Set) = copy(s)
function intersect(s::Set, sets::Set...)
    i = copy(s)
    for x in s
        for t in sets
            if !in(x,t)
                delete!(i,x)
                break
            end
        end
    end
    return i
end

function setdiff(a::Set, b::Set)
    d = copy(a)
    for x in b
        delete!(d, x)
    end
    d
end

isequal(l::Set, r::Set) = (length(l) == length(r)) && (l <= r)
<(l::Set, r::Set) = (length(l) < length(r)) && (l <= r)
<=(l::Set, r::Set) = issubset(l, r)

function issubset(l, r)
    for elt in l
        if !in(elt, r)
            return false
        end
    end
    return true
end

function unique(C)
    out = Array(eltype(C),0)
    seen = Set{eltype(C)}()
    for x in C
        if !in(x, seen)
            push!(seen, x)
            push!(out, x)
        end
    end
    out
end

function filter!(f::Function, s::Set)
    for x in s
        if !f(x)
            delete!(s, x)
        end
    end
    return s
end
filter(f::Function, s::Set) = filter!(f, copy(s))

hash(s::Set) = hash(sort(s.dict.keys[s.dict.slots .!= 0]))

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
# syntax: proto3
using ProtoBuf
import ProtoBuf.meta

mutable struct Point <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Point(; kwargs...)
        obj = new(meta(Point), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Point
const __meta_Point = Ref{ProtoMeta}()
function meta(::Type{Point})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Point)
            __meta_Point[] = target = ProtoMeta(Point)
            allflds = Pair{Symbol,Union{Type,String}}[:latitude => Int32, :longitude => Int32]
            meta(target, Point, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Point[]
    end
end
function Base.getproperty(obj::Point, name::Symbol)
    if name === :latitude
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :longitude
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

mutable struct Rectangle <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Rectangle(; kwargs...)
        obj = new(meta(Rectangle), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Rectangle
const __meta_Rectangle = Ref{ProtoMeta}()
function meta(::Type{Rectangle})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Rectangle)
            __meta_Rectangle[] = target = ProtoMeta(Rectangle)
            allflds = Pair{Symbol,Union{Type,String}}[:lo => Point, :hi => Point]
            meta(target, Rectangle, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Rectangle[]
    end
end
function Base.getproperty(obj::Rectangle, name::Symbol)
    if name === :lo
        return (obj.__protobuf_jl_internal_values[name])::Point
    elseif name === :hi
        return (obj.__protobuf_jl_internal_values[name])::Point
    else
        getfield(obj, name)
    end
end

mutable struct Feature <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function Feature(; kwargs...)
        obj = new(meta(Feature), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct Feature
const __meta_Feature = Ref{ProtoMeta}()
function meta(::Type{Feature})
    ProtoBuf.metalock() do
        if !isassigned(__meta_Feature)
            __meta_Feature[] = target = ProtoMeta(Feature)
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :location => Point]
            meta(target, Feature, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_Feature[]
    end
end
function Base.getproperty(obj::Feature, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :location
        return (obj.__protobuf_jl_internal_values[name])::Point
    else
        getfield(obj, name)
    end
end

mutable struct RouteNote <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RouteNote(; kwargs...)
        obj = new(meta(RouteNote), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct RouteNote
const __meta_RouteNote = Ref{ProtoMeta}()
function meta(::Type{RouteNote})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RouteNote)
            __meta_RouteNote[] = target = ProtoMeta(RouteNote)
            allflds = Pair{Symbol,Union{Type,String}}[:location => Point, :message => AbstractString]
            meta(target, RouteNote, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RouteNote[]
    end
end
function Base.getproperty(obj::RouteNote, name::Symbol)
    if name === :location
        return (obj.__protobuf_jl_internal_values[name])::Point
    elseif name === :message
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    else
        getfield(obj, name)
    end
end

mutable struct RouteSummary <: ProtoType
    __protobuf_jl_internal_meta::ProtoMeta
    __protobuf_jl_internal_values::Dict{Symbol,Any}
    __protobuf_jl_internal_defaultset::Set{Symbol}

    function RouteSummary(; kwargs...)
        obj = new(meta(RouteSummary), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            if fldval !== nothing
                values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
            end
        end
        obj
    end
end # mutable struct RouteSummary
const __meta_RouteSummary = Ref{ProtoMeta}()
function meta(::Type{RouteSummary})
    ProtoBuf.metalock() do
        if !isassigned(__meta_RouteSummary)
            __meta_RouteSummary[] = target = ProtoMeta(RouteSummary)
            allflds = Pair{Symbol,Union{Type,String}}[:point_count => Int32, :feature_count => Int32, :distance => Int32, :elapsed_time => Int32]
            meta(target, RouteSummary, allflds, ProtoBuf.DEF_REQ, ProtoBuf.DEF_FNUM, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_RouteSummary[]
    end
end
function Base.getproperty(obj::RouteSummary, name::Symbol)
    if name === :point_count
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :feature_count
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :distance
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :elapsed_time
        return (obj.__protobuf_jl_internal_values[name])::Int32
    else
        getfield(obj, name)
    end
end

# service methods for RouteGuide
const _RouteGuide_methods = MethodDescriptor[
        MethodDescriptor("GetFeature", 1, Point, Feature),
        MethodDescriptor("ListFeatures", 2, Rectangle, Channel{Feature}),
        MethodDescriptor("RecordRoute", 3, Channel{Point}, RouteSummary),
        MethodDescriptor("RouteChat", 4, Channel{RouteNote}, Channel{RouteNote})
    ] # const _RouteGuide_methods
const _RouteGuide_desc = ServiceDescriptor("routeguide.RouteGuide", 1, _RouteGuide_methods)

RouteGuide(impl::Module) = ProtoService(_RouteGuide_desc, impl)

mutable struct RouteGuideStub <: AbstractProtoServiceStub{false}
    impl::ProtoServiceStub
    RouteGuideStub(channel::ProtoRpcChannel) = new(ProtoServiceStub(_RouteGuide_desc, channel))
end # mutable struct RouteGuideStub

mutable struct RouteGuideBlockingStub <: AbstractProtoServiceStub{true}
    impl::ProtoServiceBlockingStub
    RouteGuideBlockingStub(channel::ProtoRpcChannel) = new(ProtoServiceBlockingStub(_RouteGuide_desc, channel))
end # mutable struct RouteGuideBlockingStub

GetFeature(stub::RouteGuideStub, controller::ProtoRpcController, inp::Point, done::Function) = call_method(stub.impl, _RouteGuide_methods[1], controller, inp, done)
GetFeature(stub::RouteGuideBlockingStub, controller::ProtoRpcController, inp::Point) = call_method(stub.impl, _RouteGuide_methods[1], controller, inp)

ListFeatures(stub::RouteGuideStub, controller::ProtoRpcController, inp::Rectangle, done::Function) = call_method(stub.impl, _RouteGuide_methods[2], controller, inp, done)
ListFeatures(stub::RouteGuideBlockingStub, controller::ProtoRpcController, inp::Rectangle) = call_method(stub.impl, _RouteGuide_methods[2], controller, inp)

RecordRoute(stub::RouteGuideStub, controller::ProtoRpcController, inp::Channel{Point}, done::Function) = call_method(stub.impl, _RouteGuide_methods[3], controller, inp, done)
RecordRoute(stub::RouteGuideBlockingStub, controller::ProtoRpcController, inp::Channel{Point}) = call_method(stub.impl, _RouteGuide_methods[3], controller, inp)

RouteChat(stub::RouteGuideStub, controller::ProtoRpcController, inp::Channel{RouteNote}, done::Function) = call_method(stub.impl, _RouteGuide_methods[4], controller, inp, done)
RouteChat(stub::RouteGuideBlockingStub, controller::ProtoRpcController, inp::Channel{RouteNote}) = call_method(stub.impl, _RouteGuide_methods[4], controller, inp)

export Point, Rectangle, Feature, RouteNote, RouteSummary, RouteGuide, RouteGuideStub, RouteGuideBlockingStub, GetFeature, ListFeatures, RecordRoute, RouteChat


text:
cd(@__DIR__); import Pkg; Pkg.activate(".")
coords = (left = -145.0, right = -50.0, bottom = 20.0, top = 75.0)
Pkg.develop(path="$(homedir())/github/SimpleSDMLayers.jl")
using Revise
using SimpleSDMLayers
using Test
using ArchGDAL
using Plots
using StatsBase

file = joinpath(dirname(pathof(SimpleSDMLayers)), "..", "data", "connectivity.tiff")

# geotiff(SimpleSDMPredictor, file)

dataset = ArchGDAL.read(file)
band = ArchGDAL.getband(dataset, 1)
values = ArchGDAL.read(dataset, 1)
ArchGDAL.getgeotransform(dataset)

landcover_mat = permutedims(values[:,end:-1:1])
landcover_mat = convert(Array{Union{Float32, Nothing}}, landcover_mat)
replace!(landcover_mat, -9999.0 => nothing)

coords = (left = -75.17734, right = -72.36486, bottom = 45.34523, top = 47.38457)

mp = SimpleSDMPredictor(landcover_mat, coords...)

qfunc = ecdf(convert(Vector{Float64}, filter(!isnothing, mp.grid)))

qmap = broadcast(qfunc, mp)

plot(qmap, frame=:grid, c=:cork, clim=(0,1))

file2 = joinpath(dirname(pathof(SimpleSDMLayers)), "..", "data", "connectivity2.tiff")
geotiff(file2, mp)
mp2 = geotiff(SimpleSDMPredictor, file2)
mp2 == mp

mp2.grid == mp.grid
mp2.left == mp.left
mp2.right == mp.right
mp2.bottom == mp.bottom
mp2.top == mp.top

# geotiff
dataset = ArchGDAL.read(file2)
band = ArchGDAL.getband(dataset, 1)
values = ArchGDAL.read(dataset, 1)
ArchGDAL.getgeotransform(dataset)

(mp.top - mp.bottom)/size(mp,1)
(mp.right - mp.left)/size(mp,2)

(mp.top - mp.bottom)/2stride(mp, 2)
(mp.right - mp.left)/2stride(mp, 1)

# find_span
layer = copy(mp)
left, right, bottom, top = (-180.0, 180.0, -90.0, 90.0)
bandnumber = 1
dataset = ArchGDAL.read(file2)
import SimpleSDMLayers: _find_span

n, m, M, pos, side = width, minlon, maxlon, right, :right

stride1 = (M - m) / n
centers = (m + 0.5stride1):stride1:(M-0.5stride1)
length(stride1) # length is wrong, should be 1206 (n)

# Whatever, I don't even know if the bounding coordinates are right in the first place...
# Overwriting the tif file
geotiff(file, mp)


text:
module ExactODEReduction

using DataStructures
using Distributions
using JSON
using OffsetArrays
using Primes
using SparseArrays
using Base.Threads

# Alex: Why uuuusing ??
using AbstractAlgebra
using Nemo

import Nemo
import AbstractAlgebra

# These are extended
import Base: ==, !=, +, -, *, lcm, rand, zero
import Base.Iterators: Stateful
import LinearAlgebra: SingularException
import Nemo: base_ring, gfp_elem, gfp_fmpz_elem, fmpq_mpoly,
            fmpq, terms, monomials, fmpz

import CustomUnitRanges: filename_for_zerorange
include(filename_for_zerorange)

# Utilities
include("utils.jl")
include("typics.jl")
include("gizmos.jl")

# Reading input
include("myeval.jl")
include("parsik.jl")

# Structures for sparse linear algebra
include("structs/sparsik.jl")
include("structs/dok_sparsik.jl")
include("structs/csr_sparsik.jl")
include("structs/densik.jl")
include("structs/subspacik.jl")

# Main functionality
include("wiedemannchik.jl")
include("basis.jl")
include("radical.jl")
include("invariants.jl")

#------------------------------------------------------------------------------

# finds one common invariant subspace of the given matrices
# using the provived default method for finding one invaiant subspace
"""
    invariant_subspace(As; backend_algorithm)

Finds *one* common invariant subspace of the given matrices `As`
using the provived default method `backend_algorithm`. By default, `invariant_subspace_1`
is used.

Each matrix in `As` to be a subtype of `AbstractSparseMatrix`
"""
function invariant_subspace(
    As;
    backend_algorithm=invariant_subspace_1)

    backend_algorithm(As)
end

#------------------------------------------------------------------------------

# finds several common invariant subspaces of the given matrices
# using the provided default method for finding one invaiant subspace
"""
    many_invariant_subspaces(As; backend_algorithm)

Finds *several* common invariant subspaces of the given matrices `As`
using the provived default method `backend_algorithm`. By default, `invariant_subspace_1`
is used.

Each matrix in `As` to be a subtype of `AbstractSparseMatrix`
"""
function many_invariant_subspaces(
    As;
    backend_algorithm=invariant_subspace_1)

    __many_invariant_subspaces(As, backend_algorithm)
end

end


text:
function compute_derivatives(obj::SecondaryODE, arg0::jdouble, arg1::Vector{jdouble}, arg2::Vector{jdouble}, arg3::Vector{jdouble})
    return jcall(obj, "computeDerivatives", Vector{jdouble}, (jdouble, Vector{jdouble}, Vector{jdouble}, Vector{jdouble}), arg0, arg1, arg2, arg3)
end

function get_dimension(obj::SecondaryODE)
    return jcall(obj, "getDimension", jint, ())
end

function init(obj::SecondaryODE, arg0::jdouble, arg1::Vector{jdouble}, arg2::Vector{jdouble}, arg3::jdouble)
    return jcall(obj, "init", void, (jdouble, Vector{jdouble}, Vector{jdouble}, jdouble), arg0, arg1, arg2, arg3)
end


text:
module EUregmort

using CSV, DataFrames, PyCall, PyPlot, Statistics

const shpreader = PyNULL()
const ccrs = PyNULL()

function __init__()
	copy!(shpreader, pyimport("cartopy.io.shapereader"))
	copy!(ccrs, pyimport("cartopy.crs"))
end

export nuts2ids, caprop_regcmp, caprop_regsexplot, fourp, fivep, caprop_mapplot, meanrate
mainpath = normpath(@__DIR__, "..")
datapath = normpath(mainpath, "data")
shpdatapath = normpath(datapath, "NUTS_2013_03M_SH", "data")
ycdr = CSV.read(normpath(datapath, "hlth_cd_ycdr2.csv"); missingstring = ":", rows_for_type_detect = 200)
nuts = CSV.read(normpath(datapath, "NUTS_AT_2013.csv"); rows_for_type_detect = 200)
nids = convert(Array, collect(skipmissing(nuts[:NUTS_ID])))
ageseq = ["Y_LT1"; "Y1-4"; "Y5-9"; "Y10-14"; "Y15-19";
	"Y20-24"; "Y25-29"; "Y30-34"; "Y35-39"; "Y40-44"; "Y45-49";
	"Y50-54"; "Y55-59"; "Y60-64"; "Y65-69"; "Y70-74"; "Y75-79";
	"Y80-84"; "Y85-89"; "Y90-94"; "Y_GE95"]
agesplitter(age) = split(age,  ['Y'; 'M'; '_'; '-'], keepempty = false)

function agealias(age)
	aspl = agesplitter(age)
	if aspl[1] == "TOTAL"
		alias = "alla åldrar"
	elseif aspl[1] == "LT1"
		alias = "0"
	elseif startswith(aspl[1], "LT")
		alias = "0\u2013$(parse(aspl[end][3:end])-1)"
	elseif startswith(aspl[1], "GE")
		alias = "$(aspl[end][3:end])\u2013"
	else
		alias = "$(aspl[1])\u2013$(aspl[end])"
	end

	if (startswith(age, "YM") || startswith(age, "M"))
		return "$alias genomsnitt över åldrar"
	else
		return alias
	end
end

sexaliases = Dict("F" => "kvinnor", "M" => "män", "T" => "alla")
causealiases = Dict("pop" => L"$10^5$ personår", "A-R_V-Y" => "alla orsaker")

function causealias(ca)
	if ca in keys(causealiases)
		return causealiases[ca]
	else
		return ca
	end
end

function nuts2ids(ctry)
	n2ids = map((x)->Symbol(x), filter((x)->(length(x)==4 && startswith(x, ctry)), nids))
	filter!((x) -> x in names(ycdr), n2ids)
end

dfarrmatch(col, arr) = map((x) -> in(x, arr), Vector(col))

function meanrate(sage, eage)
	sexframes = Dict()
	sind = indexin([sage], ageseq)[1]
	eind = indexin([eage], ageseq)[1]
	agelist = ageseq[sind:eind]
	if sage == "Y_LT1"
		if eage == "Y_GE95"
			agestr = "MTOTAL"
		else
			agestr = "YM_LT$(parse(agesplitter(eage)[end])+1)"
		end
	else
		if eage == "Y_GE95"
			agestr = "YM_GE$(agesplitter(sage)[1])"
		else
			agestr = "YM$(agesplitter(sage)[1])-$(agesplitter(eage)[end])"
		end
	end
	for sex in keys(sexaliases)
		ycdr_sub = ycdr[((ycdr[:sex].==sex) & (dfarrmatch(ycdr[:age], agelist))),
			 [4;6:size(ycdr)[2]]]
		ycdr_sub_agg = aggregate(ycdr_sub, :icd10, mean)
		aggregs = rename((x)->Symbol(replace(string(x), "_mean", "")),
			ycdr_sub_agg[:, 2:end])
		aggl = size(ycdr_sub_agg)[1]
		metaframe = DataFrame(unit = fill("RT", aggl), sex = fill(sex, aggl),
			age = fill(agestr, aggl),
			icd10 = ycdr_sub_agg[:icd10], time_geo = fill(2013, aggl))
		sexframes[sex] = hcat(metaframe, aggregs)
	end
	vcat(values(sexframes)...)
end

function ycdr_regcmp(nutsregs, sex, age, icd10, time_geo, inframe)
	inframe_sub = DataFrame(unit = inframe[:unit], sex = inframe[:sex],
		age = inframe[:age], icd10 = inframe[:icd10],
		time_geo = inframe[:time_geo])
	for reg in nutsregs
		inframe_sub[reg] = inframe[reg]
	end
	inframe_long = stack(inframe_sub, nutsregs)
	inframe_long[((inframe_long[:sex].==sex) .& (inframe_long[:age].==age)
		.& (inframe_long[:icd10].==icd10) .& (inframe_long[:time_geo].==time_geo)), :]
end

function caprop_regcmp(nutsregs, sex, age, ca1, ca2, time_geo, inframe)
	ca1frame = ycdr_regcmp(nutsregs, sex, age, ca1, time_geo, inframe)
	if ca2 == "pop"
		capropframe = ca1frame
	else
		ca2frame = ycdr_regcmp(nutsregs, sex, age, ca2, time_geo, inframe)
		capropframe = copy(ca1frame)
		capropframe[:value] = ca1frame[:value]./ca2frame[:value]
	end
	capropframe
end

function caprop_regsexplot(nutsregs, age, ca1; ca2 = "A-R_V-Y", time_geo = 2013,
	inframe = ycdr)
	fframe = caprop_regcmp(nutsregs, "F", age, ca1, ca2, time_geo, inframe)
	mframe = caprop_regcmp(nutsregs, "M", age, ca1, ca2, time_geo, inframe)
	scatter(fframe[:value], mframe[:value], alpha=0.3)
	for i in 1:size(nutsregs)[1]
		text(fframe[:value][i], mframe[:value][i], fframe[:variable][i],
			ha = "center", va = "center")
	end
	grid(1)
	xlabel(sexaliases["F"])
	ylabel(sexaliases["M"])
	title(*("Döda $(time_geo) $(causealias(ca1))/$(causealias(ca2)) ",
		"$(agealias(age))"))
end

perc_round(value) = replace("$(round(value; digits=4))", "." => ",")
fourp(prop) =
	[
	Dict("col" => "lightyellow", "value" => quantile(prop, 1/4));
	Dict("col" => "yellow", "value" => quantile(prop, 2/4));
	Dict("col" => "tomato", "value" => quantile(prop, 3/4));
	Dict("col" => "red", "value" => quantile(prop, 1))
	]
fivep(prop) =
	[
	Dict("col" => "lightyellow", "value" => quantile(prop, 1/5));
	Dict("col" => "yellow", "value" => quantile(prop, 2/5));
	Dict("col" => "orange", "value" => quantile(prop, 3/5));
	Dict("col" => "tomato", "value" => quantile(prop, 4/5));
	Dict("col" => "red", "value" => quantile(prop, 1))
	]

function caprop_mapplot(nutsregs, sex, age, ca1; ca2 = "A-R_V-Y", time_geo = 2013,
		percfunc = fivep,
		shapefname = normpath(shpdatapath, "NUTS_RG_03M_2013_3034.shp"),
		inframe = ycdr)
	pframe = caprop_regcmp(nutsregs, sex, age, ca1, ca2, time_geo, inframe)
	region_shp = shpreader[:Reader](shapefname)
	regstrings = map((x)->string(x), nutsregs)
	proj = ccrs[:LambertConformal](central_longitude = 10, central_latitude = 52,
		standard_parallels = (35,65), false_easting = 4000000,
		false_northing = 2800000, globe = ccrs[:Globe](ellipse = "GRS80"))
	ax = plt[:axes](projection = proj)
	prop = pframe[:value]
	propdict = Dict(zip(regstrings, prop))
	quantiles = percfunc(prop)
	boundlist = []
	fcolor = "blue"
	for region_rec in region_shp[:records]()
		nid = region_rec[:attributes]["NUTS_ID"]
		if nid in regstrings
			boundlist = vcat(boundlist, region_rec[:bounds])
			xmean = mean([boundlist[end][1];
				boundlist[end][3]])
			ymean = mean([boundlist[end][2];
				boundlist[end][4]])
			for quantile in quantiles
				if propdict[nid] <= quantile["value"]
					fcolor = quantile["col"]
					break
				end
			end
			ax[:add_geometries](region_rec[:geometry], proj,
				edgecolor = "black", facecolor = fcolor)
			ax[:annotate](nid, (xmean, ymean), ha = "center")
		end
	end
	xminimum = minimum([bound[1] for bound in boundlist])
	xmaximum = maximum([bound[3] for bound in boundlist])
	yminimum = minimum([bound[2] for bound in boundlist])
	ymaximum = maximum([bound[4] for bound in boundlist])
	ax[:set_xlim](xminimum, xmaximum)
	ax[:set_ylim](yminimum, ymaximum)
	percpatches = []
	perclabels = []
	for (i, quantile) in enumerate(quantiles)
		percpatch = matplotlib[:patches][:Rectangle]((0, 0), 1, 1,
			facecolor = quantile["col"])
		percpatches = vcat(percpatches, percpatch)
		if i == 1
			perclabel = *("\u2265", perc_round(minimum(prop)),
				"\n\u2264", perc_round(quantile["value"]))
		else
			perclabel = *("\u2264", perc_round(quantile["value"]))
		end
		perclabels = vcat(perclabels, perclabel)
	end
	legend(percpatches, perclabels, loc = "upper left",
		framealpha = 0.75, bbox_to_anchor=(1,1))
	title(*("Döda $(time_geo) $(causealias(ca1))/$(causealias(ca2)) ",
		"$(sexaliases[sex]) $(agealias(age))"))
	subplots_adjust(right = 0.8)
end

end # module


text:
bouts = preprocess_bouts(pokes3)
open_html_table(bouts[1:1000,:])