text:
abstract type Recorder end

struct RecStateNow{PS<:PopulationState}
  ps::PS
  isdone::Ref{Bool}
  krec::Int64 # record every k timesteps
  nrecords::Int64 # max recorded steps
  k_start::Int64 # starts to record only after this
  idx_save::Vector{Int64} # which neurons to save. empty -> ALL neurons
  times::Vector{Float64}
  state_now::Matrix{Float64}
  function RecStateNow(ps::PS,everyk::Integer,dt::Float64,Tend::Float64;
      idx_save::Vector{Int64}=Int64[],Tstart::Float64=0.0) where PS
    k_start = floor(Int64,Tstart/dt)
    nrecs = floor(Int64,(Tend-Tstart)/(everyk*dt))
    isdone = Ref(false)
    if nrecs < 1
      @warn "something wrong in the recorder... I assume this is a test ?"
      nrecs=1
      isdone = Ref(true)
    end
    times = fill(NaN,nrecs)
    nneus = isempty(idx_save) ? nneurons(ps) : length(idx_save)
    states = fill(NaN,nneus,nrecs)
    return new{PS}(ps,isdone,everyk,nrecs,k_start,idx_save,times,states)
  end
end

function reset!(rec::RecStateNow)
  fill!(rec.times,NaN)
  fill!(rec.state_now,NaN)
  rec.isdone[]=false
  return nothing
end

function (rec::RecStateNow)(t::Float64,k::Integer,::AbstractNetwork)
  # I assume k starts from 1, which corresponds to t=0
  kless,_rem=divrem((k-1-rec.k_start),rec.krec)
  kless += 1 # vector index must start from 1
  if (_rem != 0) || kless<=0 || rec.isdone[]
    return nothing
  end
  if kless > rec.nrecords
    rec.isdone[]=true
    return nothing
  end
  if isempty(rec.idx_save)
    rec.state_now[:,kless] .= rec.ps.state_now
  else
    rec.state_now[:,kless] .= rec.ps.state_now[rec.idx_save]
  end
  rec.times[kless] = t
  return nothing
end

# counts the number of spikes
# to simplify things I do it for all neurons
struct RecCountSpikes{PS<:PopulationState}
  ps::PS
  k_warmup::Int64 # starts to record only after this
  spikecount::Vector{Int64} # how many spikes for each neuron?
  function RecCountSpikes(ps::PS,dt::Float64;
      idx_save::Vector{Int64}=Int64[],t_warmup::Float64=0.0) where PS
    nneus = isempty(idx_save) ? nneurons(ps) : length(idx_save)
    k_warmup = floor(Int64,t_warmup/dt)
    spikecount = fill(0,nneus)
    return new{PS}(ps,k_warmup,spikecount)
  end
end

function reset!(rec::RecCountSpikes)
  fill!(rec.spikecount,0)
  return nothing
end

function (rec::RecCountSpikes)(t::Float64,k::Integer,ntw::AbstractNetwork)
  if k<rec.k_warmup
    return nothing
  else
    for i in findall(rec.ps.isfiring)
      rec.spikecount[i] += 1
    end
    return nothing
  end
end

# utiliy function
function get_mean_rates(rec::RecCountSpikes,dt::Float64,Ttot::Float64)
  ΔT = Ttot-(dt*rec.k_warmup)
  return rec.spikecount ./ ΔT
end

struct RecSpikes{PS<:PopulationState}
  ps::PS
  isdone::Ref{Bool}
  nrecords::Int64 # max recorded steps
  k_now::Ref{Int64} # current spike counter
  Tstart::Float64 # starts to record only after this time
  Tend::Float64  # max recorded time
  idx_save::Vector{Int64} # which neurons to save. empty -> ALL neurons
  spiketimes::Vector{Float32} # resolution 1E-7 , still below dt
  spikeneurons::Vector{UInt64} # max is about 60_000 neurons, which is plenty
  function RecSpikes(ps::PS, expected_rate::Float64,Tend::Float64;
      idx_save::Vector{Int64}=Int64[],Tstart::Float64=0.0,
      nrecmax::Int64 = 10_000_000) where PS
    @assert (isempty(idx_save)  || issorted(idx_save)) "Plase store indexes to save in sorted order"
    nneus = isempty(idx_save) ? nneurons(ps) : length(idx_save)
    if nneus == nneurons(ps)  # if 1:Ntot, no need to keep it
      idx_save = Int64[]
    end
    nrecs = ceil(Integer,nneus * expected_rate * (Tend-Tstart))
    isdone = Ref(false)
    if nrecs < 1
      @warn "something wrong in the recorder... I assume this is a test ?"
      nrecs=1
      isdone = Ref(true)
    end
    @assert nrecs <= nrecmax "Saving data might require too much memory!"
    spiketimes = fill(NaN32,nrecs)
    spikeneurons = fill(UInt16(0),nrecs) # no neuron should be 0
    return new{PS}(ps,isdone,nrecs,Ref(1),Tstart,Tend,idx_save,
      spiketimes,spikeneurons)
  end
end
function reset!(rec::RecSpikes)
  rec.k_now[]=1
  rec.isdone[]=false
  fill!(rec.spikeneurons,UInt16(0))
  fill!(rec.spiketimes,NaN32)
  return nothing
end

# easy to save wrapper
struct RecSpikesContent
  Tstart::Float64 # starts to record only after this time
  Tend::Float64  # max recorded time
  idx_save::Vector{Int64} # which neurons to save. empty -> ALL neurons
  spiketimes::Vector{Float32} # resolution 1E-7 , still below dt
  spikeneurons::Vector{UInt64} # max is about 60_000 neurons, which is plenty
  function RecSpikesContent(r::RecSpikes)
    to_keep = findall(isfinite,r.spiketimes)
    new(r.Tstart,r.Tend,r.idx_save,r.spiketimes[to_keep],r.spikeneurons[to_keep])
  end
end
function get_content(rec::RecSpikes)
  return RecSpikesContent(rec)
end

# little helper function
function _in_sorted(v::Vector{I},el::I) where I
  return !isempty(searchsorted(v,el))
end

function (rec::RecSpikes)(t::Float64,::Integer,::AbstractNetwork)
  if t < rec.Tstart || rec.isdone[]
    return nothing
  end
  if (rec.k_now[] > rec.nrecords) || (t > rec.Tend )
    rec.isdone[]=true
    return nothing
  end
  read_all_neus = isempty(rec.idx_save)
  for (neu,isfiring) in enumerate(rec.ps.isfiring)
    if isfiring && (read_all_neus || _in_sorted(rec.idx_save,neu))
      know = rec.k_now[]
      rec.spiketimes[know] = Float32(t)
      rec.spikeneurons[know] = UInt16(neu)
      rec.k_now[] = know+1
      if know+1 > rec.nrecords # must check for each spike
        rec.isdone[]=true
        return nothing
      end
    end
  end
  return nothing
end

function get_spiketimes_spikeneurons(rec::Union{RecSpikes,RecSpikesContent})
  _idx = findall(isfinite,rec.spiketimes)
  return Float64.(rec.spiketimes[_idx]),Int64.(rec.spikeneurons[_idx])
end

function get_spiketimes_dictionary(rec::Union{RecSpikes,RecSpikesContent})
  spk_t,spk_neu = get_spiketimes_spikeneurons(rec)
  ret=Dict{Int64,Vector{Float64}}()
  for neu in unique(spk_neu)
    idx = spk_neu.==neu
    ret[neu] = spk_t[idx]
  end
  return ret
end

function get_mean_rates(rec::Union{RecSpikes,RecSpikesContent};
    (Tend::Float64)=0.0)
  Tend = min(Tend,rec.Tend)
  ΔT = Tend-rec.Tstart
  dict = get_spiketimes_dictionary(rec)
  ret = Dict{Int64,Float64}()
  for (neu,spks) in pairs(dict)
    ret[neu] = length(spks)/ΔT
  end
  return ret
end

# adds a pulse at time of spiking
# voltage_traces : rows are neurons, columns are timesteps
function add_fake_spikes!(v_spike::R,
    voltage_times::Vector{R},
    voltage_traces::Matrix{R},
    spiketimes::Vector{R},
    spikeneurons::Vector{I},
    idx_save::Vector{I}) where {R,I}
  nneus,nvtimes = size(voltage_traces)
  if isempty(idx_save)
    idx_save = collect(1:nneus)
  end
  for i in 1:nneus
    idx_neu = idx_save[i]
    idx_tspike = findall(==(idx_neu),spikeneurons)
    for k in idx_tspike
      tspike = spiketimes[k]
      idx_vtime = searchsortedfirst(voltage_times,tspike)
      if idx_vtime <= nvtimes
        voltage_traces[i,idx_vtime] = v_spike
      end
    end
  end
  return nothing
end

function add_fake_spikes!(v_spike::Float64,rtrace::RecStateNow,rspk::Union{RecSpikes,RecSpikesContent})
  spiketimes,spikeneurons = get_spiketimes_spikeneurons(rspk)
  return add_fake_spikes!(v_spike,rtrace.times,rtrace.state_now,
    spiketimes,spikeneurons,rtrace.idx_save)
end

##

function binned_spikecount(dt::Float64,rspk::Union{RecSpikes,RecSpikesContent};
    Nneurons::Int64=-1,Tend::Float64=-1.0)
  return binned_spikecount(dt,get_spiketimes_spikeneurons(rspk)...;
    Nneurons=Nneurons,Tend=max(Tend,rspk.Tend),Tstart=rspk.Tstart)
end
function binned_spikecount(dt::Float64,spktimes::Vector{Float64},
    spkneurons::Vector{Int64};Nneurons::Int64=-1,
      Tend::Float64=-1.0,Tstart::Float64=0.0)
  Tend = max(Tend, maximum(spktimes)+dt)
  Nneus = max(Nneurons,maximum(spkneurons))
  tbins = Tstart:dt:Tend
  ntimes = length(tbins)-1
  binnedcount = fill(0,(Nneus,ntimes))
  for (t,neu) in zip(spktimes,spkneurons)
    if (tbins[1] < t <= tbins[end]) # just in case
      tidx = searchsortedfirst(tbins,t)-1
      binnedcount[neu,tidx]+=1
    end
  end
  binsc=midpoints(tbins)
  return binsc,binnedcount
end

# convolve spiketimes with negative exponential
# for SINGLE / GROUP of neurons (not all)
function get_spike_exp_convolution(dt::R,tau_exp::R,
    spiktimes::Vector{R},spkneurons::Vector{I},
    Nneurons::I,Tstart::R,Tend::R ;
    idx_neurons::Union{Vector{I},Nothing}=nothing) where {R,I}
  binnedt,binned_spikes = binned_spikecount(dt,spiktimes,spkneurons;
    Nneurons=Nneurons,Tstart=Tstart,Tend=Tend)
  @assert tau_exp > dt "convolution is too narrow!"
  convvect = @. (exp(-(0:dt:(10tau_exp))/tau_exp))/tau_exp
  idx_neurons = something(idx_neurons,(1:Nneurons))
  bin_select = sum(view(binned_spikes,idx_neurons,:);dims=1)[:] ./ length(idx_neurons)
  ret = conv(convvect,bin_select)
  return binnedt,ret[1:length(binnedt)]
end

function get_spike_exp_convolution(dt::R,tau_exp::R,
    rspk::Union{RecSpikes,RecSpikesContent},Nneurons::I;
    idx_neurons::Union{Vector{I},Nothing}=nothing) where {R,I}
  return get_spike_exp_convolution(dt,tau_exp,
    get_spiketimes_spikeneurons(rspk)...,Nneurons,
    rspk.Tstart,rspk.Tend;
    idx_neurons=idx_neurons)
end

# simular to the above, but averages over selected neurons
# and returns a value in Hz
function get_psth(idxs_neu::AbstractVector{<:Integer},
    dt::Float64,rspk::Union{RecSpikes,RecSpikesContent};
    Nneurons::Int64=-1,Tend::Float64=-1.0)
  return get_psth(idxs_neu,dt,get_spiketimes_spikeneurons(rspk)...;
    Nneurons=Nneurons,Tend=max(Tend,rspk.Tend),Tstart=rspk.Tstart)
end
function get_psth(idxs_neu::AbstractVector{<:Integer},dt::Float64,spktimes::Vector{Float64},
  spkneurons::Vector{Int64};Nneurons::Int64=-1,Tend::Float64=-1.0,
  Tstart::Float64=0.0)
 tmid,counts = binned_spikecount(dt,spktimes,spkneurons;
  Nneurons=Nneurons,Tend=Tend,Tstart=Tstart)
 ret2 = mean(view(counts,idxs_neu,:);dims=1)[:]
 return tmid,ret2 ./ dt
end

function raster_png(dt::Float64,rspk::Union{RecSpikes,RecSpikesContent};
    Nneurons::Int64=-1,
    Tend::Float64=0.0,spike_height::Int64=5,
    reorder::Vector{Int64}=Int64[])
  spkt,spkn = get_spiketimes_spikeneurons(rspk)
  if iszero(Tend)
    Tend = max(rspk.Tend,maximum(spkt)+dt)
  end
  if isempty(spkn) && (Nneurons==-1)
    error("""
     No spikes recorded !
    Impossible to determine the image size.
    Please set the parameters `Nneurons` and `Tend`
    """)
  end
  Nneurons = let maxn = isempty(spkn) ? Nneurons : maximum(spkn)
     max(Nneurons,maxn)
  end
  tbins = rspk.Tstart-eps(rspk.Tstart):dt:Tend
  ntimes = length(tbins)-1
  rasterbin = falses(Nneurons,ntimes)
  for (neu,t) in zip(spkn,spkt)
    if t >= tbins[1] # just in case
      tidx = searchsortedfirst(tbins,t)-1
      rasterbin[neu,tidx]=true
    end
  end
  if !isempty(reorder)
    rasterbin .= rasterbin[reorder,:]
  end
  ret = @. RGB(Gray(!rasterbin))
  if spike_height > 1
    ret = repeat(ret;inner=(spike_height,1))
  end
  return ret
end

# record weights
# might be useful when testing plasticity at small scale,
# otherwise use RecWeightsFull, below

struct RecWeightsApplyFunction
  weights::SparseMatrixCSC{Float64,Int64}
  isdone::Ref{Bool}
  krec::Int64 # record every k timesteps
  nrecords::Int64 # max recorded steps
  k_start::Int64 # starts to record only after this
  fun_things::Function # stuff to do, takes SparseMatrix, returns vector
  out_size::Int64 # size of function output vector
  times::Vector{Float64}
  fun_output::Matrix{Float64}
  function RecWeightsApplyFunction(conn::Connection,thefun::Function,
      everyk::Integer,dt::Float64,Tmax::Float64;
      Tstart::Float64=0.0)
    if iszero(n_plasticity_rules(conn))
      @warn "The connection has no plasticity, there is no need of tracking it!"
    end
    out_sample = thefun(conn.weights)
    typeassert(out_sample,Vector{Float64})
    out_size = length(out_sample)
    k_start = floor(Int64,Tstart/dt)
    nrecs = floor(Int64,(Tmax-Tstart)/(everyk*dt))
    isdone = Ref(false)
    if nrecs < 1
      @warn "something wrong in the recorder... I assume this is a test ?"
      nrecs=1
      isdone = Ref(true)
    end
    times = fill(NaN,nrecs)
    fun_output = fill(NaN,out_size,nrecs)
    return new(conn.weights,isdone,everyk,nrecs,k_start,
      thefun,out_size,
      times,fun_output)
  end
end
Base.length(rec::RecWeightsApplyFunction) = length(rec.times)

# easy to save wrapper
struct RecWeightsApplyFunctionContent
  times::Vector{Float64}
  fun_output::Matrix{Float64}
  function RecWeightsApplyFunctionContent(r::RecWeightsApplyFunction)
    new(r.times,r.fun_output)
  end
end
function get_content(rec::RecWeightsApplyFunction)
  return RecWeightsApplyFunctionContent(rec)
end

function reset!(rec::RecWeightsApplyFunction)
  fill!(rec.times,NaN)
  fill!(rec.fun_output,NaN)
  rec.isdone[]=false
  return nothing
end

function (rec::RecWeightsApplyFunction)(t::Float64,k::Integer,::AbstractNetwork)
  kless,_rem=divrem((k-1-rec.k_start),rec.krec)
  kless += 1 # vector index must start from 1
  if (_rem != 0) || kless<=0 || rec.isdone[]
    return nothing
  end
  if kless > rec.nrecords
    rec.isdone[]=true
    return nothing
  end
  rec.fun_output[:,kless] .= rec.fun_things(rec.weights)
  rec.times[kless] = t
  return nothing
end

# function rebuild_weights(rec::RecWeights)
#   times = rec.times
#   ntimes = length(times)
#   wout = Vector{SparseMatrixCSC{Float64,Int64}}(undef,ntimes)
#   for kt in 1:ntimes
#     _w = zeros(Float64,size(rec.weights))
#     for (k,ij) in enumerate(rec.idx_save)
#       _w[ij] = rec.weights_now[k,kt]
#     end
#     wout[kt] = sparse(_w)
#   end
#   return times,wout
# end

# records the full weight matrix
# this makes life much easier where there is structural plasticity.
# BUT please use it with a large k (e.g. 2 Hz )

struct RecWeightsFull
  weights::SparseMatrixCSC{Float64,Int64}
  isdone::Ref{Bool}
  krec::Int64 # record every k timesteps
  nrecords::Int64 # max recorded steps
  k_start::Int64 # starts to record only after this
  times::Vector{Float64}
  weights_now:: Vector{SparseMatrixCSC{Float64,Int64}}
  function RecWeightsFull(conn::Connection,
      everyk::Integer,dt::Float64,Tend::Float64;
      Tstart::Float64=0.0)
    if iszero(n_plasticity_rules(conn))
      @warn "The connection has no plasticity, there is no need of tracking it!"
    end
    k_start = floor(Int64,Tstart/dt)
    nrecs = floor(Int64,(Tend-Tstart)/(everyk*dt))
    isdone = Ref(false)
    if nrecs < 1
      @warn "something wrong in the recorder... I assume this is a test ?"
      nrecs=1
      isdone = Ref(true)
    end
    times = fill(NaN,nrecs)
    weights_now = Vector{SparseMatrixCSC{Float64,Int64}}(undef,nrecs)
    return new(conn.weights,isdone,everyk,nrecs,k_start,times,weights_now)
  end
end
Base.length(rec::RecWeightsFull) = length(rec.times)

# easy to save wrapper
struct RecWeightsFullContent
  times::Vector{Float64}
  weights_now::Vector{SparseMatrixCSC{Float64,Int64}}
  function RecWeightsFullContent(r::RecWeightsFull)
    new(r.times,r.weights_now)
  end
end
function get_content(rec::RecWeightsFull)
  return RecWeightsFullContent(rec)
end

function reset!(rec::RecWeightsFull)
  fill!(rec.times,NaN)
  fill!(rec.weights_now,sparse(zeros(1,1)))
  rec.isdone[]=false
  return nothing
end

function (rec::RecWeightsFull)(t::Float64,k::Integer,::AbstractNetwork)
  # I assume k starts from 1, which corresponds to t=0
  kless,_rem=divrem((k-rec.k_start-1),rec.krec)
  kless += 1 # vector index must start from 1
  if (_rem != 0) || kless<=0 || rec.isdone[]
    return nothing
  end
  if kless > rec.nrecords
    rec.isdone[]=true
    return nothing
  end
  rec.weights_now[kless] = copy(rec.weights)
  rec.times[kless] = t
  return nothing
end

# utility functions for differences in weights

_wdiff_plain(w1::R,w2::R) where R<:Real = w1 - w2
_wdiff_plain(wmat1::Matrix{R},wmat2::Matrix{R}) where R<:Real = _wdiff_plain.(wmat1,wmat2)
function _wdiff_rel(w1::R,w2::R) where R<:Real
  wavg = 0.5(w1+w2)
  (wavg == 0.0) && (return 0.0)
  return (w1-w2)/wavg
end
_wdiff_rel(wmat1::Matrix{R},wmat2::Matrix{R}) where R<:Real = _wdiff_rel.(wmat1,wmat2)

function _plus_minus_rescale(mat::Array{R}) where R
  down,up = extrema(mat)
  down=abs(down)
  ret = copy(mat)
  for (i,x) in enumerate(ret)
    if x > 0
        ret[i] = x/up
    elseif x < 0
        ret[i] = x/down
    end
  end
  return ret
end
function _wdiff_pm(wmat1::Matrix{R},wmat2::Matrix{R}) where R
  return _plus_minus_rescale(wmat1.-wmat2)
end

# +1 if both elements exist , 0 if both absent, -1 if mismatch
function _wdiff_structure(wmat1::Matrix{R},wmat2::Matrix{R}) where R
  ret = Matrix{Int64}(undef,size(wmat1)...)
  for i in eachindex(ret)
    if (wmat1[i] != 0.0) && (wmat2[i] != 0.0)
      ret[i] = 1
    elseif (wmat1[i] != 0.0) || (wmat2[i] != 0.0)
      ret[i] = -1
    else
      ret[i] = 0
    end
  end
  return ret
end

function compare_weigth_matrices(weights1::AbstractMatrix,
    weights2::AbstractMatrix,compare_funs...)
  if isempty(compare_funs)
    compare_funs=(_wdiff_plain,_wdiff_rel,_wdiff_pm)
  end
  w1 = Matrix(weights1)
  w2 = Matrix(weights2)
  @assert size(w1) == size(w2)
  ret = map(compare_funs) do _f
    _f(w1,w2)
  end
  return ret
end


text:
import Base.==

immutable Length{U, T} <: Measure
    value::T
end
Length{T}(unit::Symbol, x::T) = Length{unit, T}(x)

typealias AbsoluteLength Length{:mm, Float64}

Base.convert{u, T1 <: Number, T2 <: Number}(::Type{Length{u, T1}}, x::Length{u, T2}) =
    Length{u, T1}(x.value)

=={u}(x::Length{u}, y::Length{u}) = x.value == y.value
Base.isequal{u}(x::Length{u}, y::Length{u}) = isequal(x.value, y.value)
Base.hash{u}(x::Length{u}) = hash(x.value, hash(u))

# Operations
# ----------

# TODO: figure out if we actually need this:
# Add{P <: Length, Q<:Length}(x::P, y::Q) = Add{P, Q}(x, y)
# Add(x::Measure, y::Measure) = Add{Measure, Measure}(x, y)

Neg{T <: Length}(x::T) = T(-x.value)
Div{u}(a::Length{u}, b::Number) = Length(u, a.value / b)
Div{u}(a::Length{u}, b::Length{u}) = a.value / b.value
Mul{u}(a::Length{u}, b::Number) = Length(u, a.value * b)
Max{u}(a::Length{u}, b::Length{u}) = Length(u, max(a.value, b.value))
Min{u}(a::Length{u}, b::Length{u}) = Length(u, min(a.value, b.value))

@compat Base.:+{u}(a::Length{u}, b::Length{u}) = Length(u, a.value + b.value)
@compat Base.:-{u}(a::Length{u}, b::Length{u}) = Length(u, a.value - b.value)

iszero(x::Length) = x.value == 0.0

Base.abs{T <: Length}(a::T) = T(abs(a.value))
Base.isless{u}(a::Length{u}, b::Length{u}) = a.value < b.value

# Constants
# ---------

const mm   = AbsoluteLength(1.0)
const cm   = AbsoluteLength(10.0)
const inch = AbsoluteLength(25.4)
const pt   = inch/72.0

const w    = Length{:w, Float64}(1.0)
const h    = Length{:h, Float64}(1.0)
const d    = Length{:d, Float64}(1.0)

ext: jl


text:
using FractalTools
using Makie

# Construct a tesselation
pts = [0., 0.], [1., 0], [0.5, 1.]
tridln = TriDelaunay(pts...)

# Plot the delaunay.
tridelaunayplot(tridln)


text:
# Copyright 2016 Martin Holters
# See accompanying license file.

import Base.isless
import Base.isempty

type KDTree{Tcv<:AbstractVector,Tp<:AbstractMatrix}
    cut_dim::Vector{Int}
    cut_val::Tcv
    ps_idx::Vector{Int}
    ps::Tp
end

function KDTree(p::AbstractMatrix, Np=size(p,2))
    function calc_cut_idx(min_idx, max_idx)
        N = max_idx - min_idx + 1
        N2 = 2^floor(Int, log2(N-1))
        if 3*div(N2, 2) ≤ N
            return min_idx+N2-1
        else
            return min_idx+N-div(N2, 2)-1
        end
    end

    if Np == 0
        return KDTree{Vector{eltype(p)},typeof(p)}([], [], [], p)
    end

    min_idx = zeros(Int, Np-1)
    max_idx = zeros(Int, Np-1)
    cut_idx = zeros(Int, Np-1)
    cut_dim = zeros(Int, Np-1)
    cut_val = zeros(eltype(p), Np-1)

    if Np == 1
        return KDTree{typeof(cut_val),typeof(p)}(cut_dim, cut_val, [1], p)
    end

    dim = indmax(var(p[:,1:Np],2))
    p_idx = sortperm(vec(p[dim,:]))

    min_idx[1] = 1
    max_idx[1] = Np
    cut_idx[1] = calc_cut_idx(min_idx[1], max_idx[1])
    cut_dim[1] = dim
    cut_val[1] = mean(p[dim, p_idx[cut_idx[1]:cut_idx[1]+1]])

    for n in 2:Np-1
        parent_n = div(n, 2)
        if mod(n, 2) == 0
            min_idx[n] = min_idx[parent_n]
            max_idx[n] = cut_idx[parent_n]
        else
            min_idx[n] = cut_idx[parent_n]+1
            max_idx[n] = max_idx[parent_n]
        end
        dim = indmax(var(p[:,p_idx[min_idx[n]:max_idx[n]]],2))
        idx = sortperm(vec(p[dim,p_idx[min_idx[n]:max_idx[n]]]))
        p_idx[min_idx[n]:max_idx[n]] = p_idx[idx + min_idx[n] - 1]
        cut_idx[n] = calc_cut_idx(min_idx[n], max_idx[n])
        cut_dim[n] = dim
        cut_val[n] = mean(p[dim, p_idx[cut_idx[n]:cut_idx[n]+1]])
    end

    p_idx_final = zeros(Int, Np)
    for n in 1:Np
        parent_n = div(n+Np-1, 2);
        if mod(n+Np, 2) == 1
            p_idx_final[n] = p_idx[min_idx[parent_n]]
        else
            p_idx_final[n] = p_idx[max_idx[parent_n]]
        end
    end

    return KDTree{typeof(cut_val),typeof(p)}(cut_dim, cut_val, p_idx_final, p)
end

immutable AltEntry{T}
    idx::Int
    delta::Vector{T}
    delta_norm::T
end

isless(e1::AltEntry, e2::AltEntry) = isless(e1.delta_norm, e2.delta_norm)

type Alts{T}
    entries::Vector{AltEntry{T}}
    best_dist::T
    best_pidx::Int
end

Alts{T}(p::Vector{T}) = Alts([AltEntry(1, zeros(p), zero(T))], typemax(T), 0)

function siftup!(alts::Alts, i)
    entries = alts.entries
    if i > length(entries)
        throw(BoundsError())
    end
    parent = div(i, 2)
    @inbounds while i > 1 && entries[i] < entries[parent]
        entries[i], entries[parent] = entries[parent], entries[i]
        i = parent
        parent = div(i, 2)
    end
end

function siftdown!(alts::Alts, i)
    if i < 1
        throw(BoundsError())
    end
    entries = alts.entries
    N = length(entries)
    @inbounds while true
        min = i
        if 2i ≤ N && entries[2i] < entries[min]
            min = 2i
        end
        if 2i+1 ≤ N && entries[2i+1] < entries[min]
            min = 2i+1
        end
        if min == i
            break
        end
        entries[i], entries[min] = entries[min], entries[i]
        i = min
    end
end

isempty(alts::Alts) = isempty(alts.entries)
peek(alts::Alts) = alts.entries[1]

function dequeue!(alts::Alts)
    e = alts.entries[1]
    alts.entries[1] = alts.entries[end]
    deleteat!(alts.entries, length(alts.entries))
    siftdown!(alts, 1)
    return e
end

function enqueue!(alts::Alts, entry::AltEntry)
    if entry.delta_norm < alts.best_dist
        push!(alts.entries, entry)
        siftup!(alts, length(alts.entries))
    end
end

function update_best_dist!(alts, dist, p_idx)
    if dist < alts.best_dist
        alts.best_dist = dist
        alts.best_pidx = p_idx
        i = length(alts.entries)
        @inbounds while i > 0
            if alts.entries[i].delta_norm .≥ alts.best_dist
                alts.entries[i] = alts.entries[end]
                last_idx = length(alts.entries)
                deleteat!(alts.entries, last_idx)
                if last_idx ≠ i
                    if i==1 || alts.entries[i] > alts.entries[div(i,2)]
                        siftdown!(alts, i)
                    else
                        siftup!(alts, i)
                    end
                end
            end
            i -= 1
        end
    end
end

indnearest(tree::KDTree, p::AbstractVector, alt = Alts(p)) =
    indnearest(tree, p, typemax(Int), alt)

function indnearest(tree::KDTree, p::AbstractVector, max_leaves::Int,
                    alt = Alts(p))
    if length(p) ≠ size(tree.ps, 1)
        throw(DimensionMismatch())
    end
    l = 0
    p_idx = 0
    while l < max_leaves && ~isempty(alt)
        best_alt = dequeue!(alt)
        idx = best_alt.idx
        delta = best_alt.delta
        delta_norm = best_alt.delta_norm

        while idx ≤ length(tree.cut_dim)
            dim = tree.cut_dim[idx]
            new_alt_delta_norm = delta_norm - delta[dim]^2 + (p[dim] - tree.cut_val[idx])^2
            if new_alt_delta_norm < alt.best_dist
                new_alt_delta = copy(delta)
                new_alt_delta[dim] = p[dim] - tree.cut_val[idx]
                new_alt_idx = p[dim] ≤ tree.cut_val[idx] ? 2idx+1 : 2idx;
                enqueue!(alt, AltEntry(new_alt_idx, new_alt_delta, new_alt_delta_norm))
            end
            if p[dim] ≤ tree.cut_val[idx]
                idx *= 2
            else
                idx = 2idx + 1
            end
        end
        idx -= length(tree.cut_dim)

        p_idx = tree.ps_idx[idx]
        if p_idx < 1 || p_idx > size(tree.ps, 2)
            throw(BoundsError)
        end
        dist = 0.
        @simd for i in 1:length(p)
            @inbounds dist += (p[i] - tree.ps[i, p_idx])^2
        end
        update_best_dist!(alt, dist, p_idx)

        l += 1
    end

    return alt.best_pidx, alt
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libidn_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libidn")
JLLWrappers.@generate_main_file("libidn", UUID("f0f5855d-14eb-5155-b6cb-9f6c251f73bb"))
end  # module libidn_jll


text:
type Shot
    nr   :: Int64
    isz  :: Int64
    isx  :: Int64
    nz   :: Int64
    nx   :: Int64
    ext  :: Int64
    iflag:: Int64
    irz  :: Array{Int64,1}
    irx  :: Array{Int64,1}
    ot   :: Float64
    dt   :: Float64
    nt   :: Int64
    Vx   :: Array{Float64,2}
    Vz   :: Array{Float64,2}
    Txx  :: Array{Float64,2}
    Tzz  :: Array{Float64,2}
    Txz  :: Array{Float64,2}
end

function InitShot(isz::Int64, isx::Int64, nz::Int64, nx::Int64, ext::Int64, iflag::Int64, irz::Array{Int64,1}, irx::Array{Int64,1}, ot::Float64, dt::Float64, nt::Int64)
    nr = length(irz)
    shot = Shot(nr, isz, isx, nz, nx, ext, iflag, irz, irx, ot, dt, nt, zeros(nt,nr), zeros(nt,nr), zeros(nt,nr), zeros(nt,nr), zeros(nt,nr))
    return shot
end

type ShotV
    nr   :: Int64
    isz  :: Int64
    isx  :: Int64
    nz   :: Int64
    nx   :: Int64
    ext  :: Int64
    iflag:: Int64
    irz  :: Array{Int64,1}
    irx  :: Array{Int64,1}
    ot   :: Float64
    dt   :: Float64
    nt   :: Int64
    Vx   :: Array{Float64,2}
    Vz   :: Array{Float64,2}
end

function InitShotV(isz::Int64, isx::Int64, nz::Int64, nx::Int64, ext::Int64, iflag::Int64, irz::Array{Int64,1}, irx::Array{Int64,1}, ot::Float64, dt::Float64, nt::Int64)
    nr = length(irz)
    shotv = ShotV(nr, isz, isx, nz, nx, ext, iflag, irz, irx, ot, dt, nt, zeros(nt,nr), zeros(nt,nr))
    return shotv
end

function WriteShot(pathout::String, shot::Shot)
    fid = open(pathout, "w")
    nr  = length(shot.irz)
    write(fid, Int32(nr),      Int32(shot.isz), Int32(shot.isx))
    write(fid, Int32(shot.nz), Int32(shot.nx) , Int32(shot.ext), Int32(shot.iflag))
    write(fid, convert(Array{Int32},shot.irz) , convert(Array{Int32},shot.irx))
    write(fid, Float32(shot.ot), Float32(shot.dt), Int32(shot.nt))
    write(fid, convert(Array{Float32},vec(shot.Vx)), convert(Array{Float32},vec(shot.Vz)), convert(Array{Float32},vec(shot.Txx)), convert(Array{Float32},vec(shot.Tzz)), convert(Array{Float32},vec(shot.Txz)))
    close(fid)
    return nothing
end

function WriteShotV(path::String, shot::ShotV)
    fid = open(path, "w")
    nr  = length(shot.irz)
    write(fid, Int32(nr),      Int32(shot.isz), Int32(shot.isx))
    write(fid, Int32(shot.nz), Int32(shot.nx) , Int32(shot.ext), Int32(shot.iflag))
    write(fid, convert(Array{Int32},shot.irz) , convert(Array{Int32},shot.irx))
    write(fid, Float32(shot.ot), Float32(shot.dt), Int32(shot.nt))
    write(fid, convert(Array{Float32},vec(shot.Vx)), convert(Array{Float32},vec(shot.Vz)))
    close(fid)
    return nothing
end

function ReadShot(pathin::String)
    fid = open(pathin, "r")
    nr  = Int64(read(fid,Int32)); isz = Int64(read(fid,Int32)); isx = Int64(read(fid,Int32));
    nz  = Int64(read(fid,Int32)); nx  = Int64(read(fid,Int32)); ext = Int64(read(fid,Int32)); iflag = Int64(read(fid,Int32));
    irz = convert(Array{Int64},read(fid,Int32,nr)); irx = convert(Array{Float64},read(fid,Int32,nr));
    ot  = Float64(read(fid,Float32))  ; dt  = Float64(read(fid,Float32))  ; nt = Int64(read(fid,Int32));
    Vx  = convert(Array{Float64}, reshape(read(fid,Float32,nt*nr),nt,nr))
    Vz  = convert(Array{Float64}, reshape(read(fid,Float32,nt*nr),nt,nr))
    Txx = convert(Array{Float64}, reshape(read(fid,Float32,nt*nr),nt,nr))
    Tzz = convert(Array{Float64}, reshape(read(fid,Float32,nt*nr),nt,nr))
    Vxz = convert(Array{Float64}, reshape(read(fid,Float32,nt*nr),nt,nr))
    close(fid)
    shot = Shot(nr, isz, isx, nz, nx, ext, iflag, irz, irx, ot, dt, nt, Vx, Vz, Txx, Tzz, Txz)
    return shot
end

function ReadShotV(pathin::String)
    fid = open(pathin, "r")
    nr  = Int64(read(fid,Int32)); isz = Int64(read(fid,Int32)); isx = Int64(read(fid,Int32));
    nz  = Int64(read(fid,Int32)); nx  = Int64(read(fid,Int32)); ext = Int64(read(fid,Int32)); iflag = Int64(read(fid,Int32));
    irz = convert(Array{Int64},read(fid,Int32,nr)); irx = convert(Array{Int64},read(fid,Int32,nr));
    ot  = Float64(read(fid,Float32))  ; dt  = Float64(read(fid,Float32))  ; nt = Int64(read(fid,Int32));
    Vx  = convert(Array{Float64}, reshape(read(fid,Float32,nt*nr),nt,nr))
    Vz  = convert(Array{Float64}, reshape(read(fid,Float32,nt*nr),nt,nr))
    close(fid)
    shot = ShotV(nr, isz, isx, nz, nx, ext, iflag, irz, irx, ot, dt, nt, Vx, Vz)
    return shot
end

function InfoShotV(pathin::String)
    fid = open(pathin, "r")
    nr  = Int64(read(fid,Int32)); isz = Int64(read(fid,Int32)); isx = Int64(read(fid,Int32));
    nz  = Int64(read(fid,Int32)); nx  = Int64(read(fid,Int32)); ext = Int64(read(fid,Int32)); iflag = Int64(read(fid,Int32));
    irz = convert(Array{Int64},read(fid,Int32,nr)); irx = convert(Array{Int64},read(fid,Int32,nr));
    ot  = Float64(read(fid,Float32))  ; dt  = Float64(read(fid,Float32))  ; nt = Int64(read(fid,Int32));
    return nr, nz, nx, ext, iflag, irz, irx, dt, nt
end

function Spt2Shot!(shot::Shot, spt::SnapShot)
    it = spt.it;
    nz=shot.nz; nx=shot.nx; ext=shot.ext; iflag=shot.iflag;
    if iflag == 1
       Nz = nz + 2*ext
       upper = ext
    elseif iflag == 2
       Nz = nz +   ext
       upper = 0
    end
    nr = length(shot.irz)
    for ir = 1: nr
        iz = shot.irz[ir]
        ix = shot.irx[ir]
        ind = (ext+ix-1)*Nz + (iz+upper)
        shot.Vx[it,ir]  = spt.Vxx[ind] + spt.Vxz[ind]
        shot.Vz[it,ir]  = spt.Vzx[ind] + spt.Vzz[ind]
        shot.Txx[it,ir] = spt.Txxx[ind] + spt.Txxz[ind]
        shot.Tzz[it,ir] = spt.Tzzx[ind] + spt.Tzzz[ind]
        shot.Txz[it,ir] = spt.Txzx[ind] + spt.Txzz[ind]
    end
    return nothing
end

function Spt2ShotV!(shotv::ShotV, spt::SnapShot)
    it = spt.it;
    nz=shotv.nz; nx=shotv.nx; ext=shotv.ext; iflag=shotv.iflag;
    if iflag == 1
       Nz = nz + 2*ext
       upper = ext
    elseif iflag == 2
       Nz = nz +   ext
       upper = 0
    end
    nr = length(shotv.irz)
    for ir = 1: nr
        iz = shotv.irz[ir]
        ix = shotv.irx[ir]
        ind = (ext+ix-1)*Nz + (iz+upper)
        shotv.Vx[it,ir]  = spt.Vxx[ind] + spt.Vxz[ind]
        shotv.Vz[it,ir]  = spt.Vzx[ind] + spt.Vzz[ind]
    end
    return nothing
end

function AddShot2Spt!(spt::SnapShot, shot::Shot)
    it = spt.it
    if spt.iflag == 1
       Nz = spt.nz + 2*spt.ext
       upper = spt.ext
    elseif spt.iflag == 2
       Nz = spt.nz +   spt.ext
       upper = 0
    end
    nr = length(shot.irz)
    for ir = 1: nr
        iz = shot.irz[ir]
        ix = shot.irx[ir]
        ind = (spt.ext+ix-1)*Nz + (iz+upper)
        spt.Vxx[ind] = spt.Vxx[ind] + shot.Vx[it,ir]
        spt.Vxz[ind] = spt.Vxz[ind] + shot.Vx[it,ir]
        spt.Vzx[ind] = spt.Vzx[ind] + shot.Vz[it,ir]
        spt.Vzz[ind] = spt.Vzz[ind] + shot.Vz[it,ir]
        spt.Txxx[ind] = spt.Txxx[ind] + shot.Txx[it,ir]
        spt.Txxz[ind] = spt.Txxz[ind] + shot.Txx[it,ir]
        spt.Tzzx[ind] = spt.Tzzx[ind] + shot.Tzz[it,ir]
        spt.Tzzz[ind] = spt.Tzzz[ind] + shot.Tzz[it,ir]
        spt.Txzx[ind] = spt.Txzx[ind] + shot.Txz[it,ir]
        spt.Txzz[ind] = spt.Txzz[ind] + shot.Txz[it,ir]
    end
    return nothing
end

function AddShotV2Spt!(spt::SnapShot, shotv::ShotV)
    it = spt.it
    if spt.iflag == 1
       Nz = spt.nz + 2*spt.ext
       upper = spt.ext
    elseif spt.iflag == 2
       Nz = spt.nz +   spt.ext
       upper = 0
    end
    nr = length(shotv.irz)
    for ir = 1: nr
        iz = shotv.irz[ir]
        ix = shotv.irx[ir]
        ind = (spt.ext+ix-1)*Nz + (iz+upper)
        spt.Vxx[ind] = spt.Vxx[ind] + shotv.Vx[it,ir]
        spt.Vxz[ind] = spt.Vxz[ind] + shotv.Vx[it,ir]
        spt.Vzx[ind] = spt.Vzx[ind] + shotv.Vz[it,ir]
        spt.Vzz[ind] = spt.Vzz[ind] + shotv.Vz[it,ir]
    end
    return nothing
end

function Dif2Records(records1, records2) #save result of records2 - records1 to records1
  if records1[1].nz != records2[1].nz || records1[1].nx != records2[1].nx || records1[1].ext != records2[1].ext
     error("size does not match")
  end
  if records1[1].nt != records2[1].nt
     error("time length does not match")
  end
  nrec = length(records1)
  for irec = 1: nrec
      records1[irec].Vx =  records2[irec].Vx -  records1[irec].Vx
      records1[irec].Vz =  records2[irec].Vz -  records1[irec].Vz
      records1[irec].Txx = records2[irec].Txx - records1[irec].Txx
      records1[irec].Tzz = records2[irec].Tzz - records1[irec].Tzz
      records1[irec].Txz = records2[irec].Txz - records1[irec].Txz
  end
  return records1
end

function Dif2RecordsNew(records1, records2) #save result of records1 - records2 to records2
  if records1[1].nz != records2[1].nz || records1[1].nx != records2[1].nx || records1[1].ext != records2[1].ext
     error("size does not match")
  end
  if records1[1].nt != records2[1].nt
     error("time length does not match")
  end
  nrec = length(records1)
  for irec = 1: nrec
      records2[irec].Vx =  records1[irec].Vx -  records2[irec].Vx
      records2[irec].Vz =  records1[irec].Vz -  records2[irec].Vz
      records2[irec].Txx = records1[irec].Txx - records2[irec].Txx
      records2[irec].Tzz = records1[irec].Tzz - records2[irec].Tzz
      records2[irec].Txz = records1[irec].Txz - records2[irec].Txz
  end
  return records2
end

function NormOfShot(shpt::Shot; itype="l2")
  tmp = 0.0
  nr  = length(shot.irz)
  if itype == "l2"
     tmp = tmp + dot(vec(shot.Vx), vec(shot.Vx )) + dot(vec(shot.Vz ),vec(shot.Vz ))
               + dot(vec(shot.Txx),vec(shot.Txx)) + dot(vec(shot.Tzz),vec(shot.Tzz)) + dot(vec(shot.Txz),vec(shot.Txz))
     tmp = sqrt(tmp)
  elseif itype == "l1"
     tmp = tmp + sum(abs(shot.Vx )) + sum(abs(shot.Vz ))
               + sum(abs(shot.Txx)) + sum(abs(shot.Tzz)) + sum(abs(shot.Txz))
  end
  return tmp
end

function imaging!(dm::Array{Float64,1}, du::Array{Float64,1}, spt::SnapShot, path::String)
    nz = spt.nz; nx = spt.nx; ext = spt.ext; iflag = spt.iflag;
    it = spt.it
    if spt.iflag == 1
       Nz = spt.nz + 2*spt.ext
       upper = spt.ext
    elseif spt.iflag == 2
       Nz = spt.nz +   spt.ext
       upper = 0
    end
    Nx = spt.nx + 2*spt.ext
    (vxx, vxz, vzx, vzz) = ReadPv(path, it)
    Txx = (spt.Txxx+spt.Txxz) * 1/2
    Tzz = (spt.Tzzx+spt.Tzzz) * 1/2
    Txz = (spt.Txzx+spt.Txzz) * 1/2
    tmp = Txx+Tzz
    mysum2!(dm, (vxx+vzz).*tmp)
    mysum2!(du, 2*vxx.*Txx)
    mysum2!(du, 2*vzz.*Tzz)
    mysum2!(du, (vxz+vzx).*Txz)
    return nothing
end

# function imaging(dm::Array{Float64,1}, du::Array{Float64,1}, spt::SnapShot, path::String)
#     nz = spt.nz; nx = spt.nx; ext = spt.ext; iflag = spt.iflag;
#     it = spt.it
#     if spt.iflag == 1
#        Nz = spt.nz + 2*spt.ext
#        upper = spt.ext
#     elseif spt.iflag == 2
#        Nz = spt.nz +   spt.ext
#        upper = 0
#     end
#     Nx = spt.nx + 2*spt.ext
#     (vxx, vxz, vzx, vzz) = ReadPv(path, it)
#     Txx = reshape(spt.Txxx+spt.Txxz, Nz, Nx) * 1/2
#     Tzz = reshape(spt.Tzzx+spt.Tzzz, Nz, Nx) * 1/2
#     Txz = reshape(spt.Txzx+spt.Txzz, Nz, Nx) * 1/2
#     tmp = Txx[upper+1:upper+nz, ext+1:ext+nx] + Tzz[upper+1:upper+nz, ext+1:ext+nx]
#     dm  = dm + (vxx+vzz) .* tmp
#     du  = du + 2*vxx.*Txx[upper+1:upper+nz, ext+1:ext+nx]
#     du  = du + 2*vzz.*Tzz[upper+1:upper+nz, ext+1:ext+nx]
#     du  = du + (vxz+vzx) .* Txz[upper+1:upper+nz, ext+1:ext+nx]
#     return dm, du
# end

function SptsPv2Born(path_spt, path_pv)
    (nz, nx, ext, iflag, dt, nt) = InfoPv(path_pv)
    if iflag == 1
       Nz = nz + 2*ext
       upper = ext
    elseif iflag == 2
       Nz = nz +   ext
       upper = 0
    end
    Nx = nx + 2*ext
    dm = zeros(nz, nx); du = zeros(nz, nx);
    for it = 1 : nt
        (vxx, vxz, vzx, vzz) = ReadPv(path_pv, it)
        spt = ReadSnapShot(path_spt, it)
        Txx = reshape(spt.Txxx+spt.Txxz, Nz, Nx) * 1/2
        Tzz = reshape(spt.Tzzx+spt.Tzzz, Nz, Nx) * 1/2
        Txz = reshape(spt.Txzx+spt.Txzz, Nz, Nx) * 1/2
        tmp = Txx[upper+1:upper+nz, ext+1:ext+nx] + Tzz[upper+1:upper+nz, ext+1:ext+nx]
        dm  = dm + (vxx+vzz) .* tmp
        du  = du + 2*vxx.*Txx[upper+1:upper+nz, ext+1:ext+nx]
        du  = du + 2*vzz.*Tzz[upper+1:upper+nz, ext+1:ext+nx]
        du  = du + (vxz+vzx) .* Txz[upper+1:upper+nz, ext+1:ext+nx]
    end
    return dm, du
end

function removeDirect(shotv::ShotV, v::Float64, src::Source, dx::Float64, dt::Float64)
    nr = length(shotv.irz)
    nt = size(shotv.Vx, 1 )
    isz = src.isz; isx = src.isx;
    lsrc= length(src.Txx)
    damp= hanning(21)[2:11]
    if shotv.irz[1] != isz
       error("source and receivers are not on same level")
    end
    for ir = 1 : nr
        dp  = ones(nt)
        l  = abs((shotv.irx[ir]-isx)*dx)
        it = ceil(Int64,l/(v*dt))+1+lsrc
        dp[1:it] = zeros(it)
        dp[it+1:it+10] = damp
        shotv.Vx[:,ir] = shotv.Vx[:,ir] .* dp
        shotv.Vz[:,ir] = shotv.Vz[:,ir] .* dp
    end
    return nothing
end

function binWfd2shotv(irz::Array{Int64,1}, irx::Array{Int64,1}, path::String)
    fid = open(path, "r")
    nr  = length(irz)
    nz = Int64(read(fid, Int32))
    nx = Int64(read(fid, Int32))
    ext= Int64(read(fid, Int32))
    iflag= Int64(read(fid, Int32))
    dt = Float64(read(fid, Float32))
    lwfd = sizeof(Float32)*nz*nx*5
    pre  = sizeof(Float32)*5
    nt = round(Int64, filesize(fid)/lwfd)
    shotv = InitShotV(0, 0, nz, nx, ext, iflag, irz, irx, 0.0, dt, nt)
    for it = 1 : nt
        pos = pre + (it-1)*lwfd
        seek(fid, pos)
        vx = reshape(read(fid, Float32, nz*nx), nz, nx)
        vz = reshape(read(fid, Float32, nz*nx), nz, nx)
        for ir = 1 : nr
            shotv.Vx[it,ir] = vx[irz[ir], irx[ir]]
            shotv.Vz[it,ir] = vz[irz[ir], irx[ir]]
        end
    end
    return shotv
end

max_stars_repo_name: WneleiGao/ElasticWave


text:
import Base: map, eachindex
import Distances: pairwise
export transform, pick_dims, periodic, scale, Indexer

"""
    ITMean{Tμ<:MeanFunction, Tf} <: MeanFunction

"InputTransformationMean": A `MeanFunction` `μ_it` defined by applying the function `f` to
the inputs to another `MeanFunction` `μ`. Concretely, `mean(μ_it, X) = mean(μ, f(X))`.
"""
struct ITMean{Tμ<:MeanFunction, Tf} <: MeanFunction
    μ::Tμ
    f::Tf
end
(μ::ITMean)(x) = μ.μ(μ.f(x))
ITMean(μ::MeanFunction, ::typeof(identity)) = μ
length(μ::ITMean) = length(μ.μ)
eachindex(μ::ITMean) = eachindex(μ.μ)
_map(μ::ITMean, X::AV) = map(μ.μ, map(μ.f, X))

"""
    ITKernel{Tk<:Kernel, Tf} <: Kernel

"InputTransformationKernel": An `ITKernel` `kit` is the kernel defined by applying a
transform `f` to the argument to a kernel `k`. Concretely: `k(x, x′) = k(f(x), f(x′))`.
"""
struct ITKernel{Tk<:Kernel, Tf} <: Kernel
    k::Tk
    f::Tf
end
(k::ITKernel)(x) = k.k(k.f(x))
(k::ITKernel)(x, x′) = k.k(k.f(x), k.f(x′))
ITKernel(k::Kernel, ::typeof(identity)) = k
length(k::ITKernel) = length(k.k)
isstationary(k::ITKernel) = isstationary(k.k)
eachindex(k::ITKernel) = eachindex(k.k)

"""
    LhsITCross{Tk<:CrossKernel, Tf} <: CrossKernel

"LhsInputTransformationCrossKernel": The kernel `k′` given by `k′(x, x′) = k(f(x), x′)`.
"""
struct LhsITCross{Tk<:CrossKernel, Tf} <: CrossKernel
    k::Tk
    f::Tf
end
(k::LhsITCross)(x, x′) = k.k(k.f(x), x′)
size(k::LhsITCross, n::Int) = size(k.k, n)

"""
    RhsITCross{Tk<:CrossKernel, Tf} <: CrossKernel

"RhsInputTransformationCrossKernel": The kernel `k′` given by `k′(x, x′) = k(x, f(x′))`.
"""
struct RhsITCross{Tk<:CrossKernel, Tf} <: CrossKernel
    k::Tk
    f::Tf
end
(k::RhsITCross)(x, x′) = k.k(x, k.f(x′))
size(k::RhsITCross, n::Int) = size(k.k, n)

"""
    ITCross{Tk<:Kernel, Tf, Tf′} <: CrossKernel

"InputTransformationCrossKernel": the kernel `k′` given by `k′(x, x′) = k(f(x), f′(x′))`.
"""
struct ITCross{Tk<:Kernel, Tf, Tf′} <: CrossKernel
    k::Tk
    f::Tf
    f′::Tf′
end
(k::ITCross)(x, x′) = k.k(k.f(x), k.f′(x′))
size(k::ITCross, n::Int) = size(k.k, n)

_map(k::ITKernel, X::AV) = map(k.k, map(k.f, X))
_map(k::ITKernel, X::AV, X′::AV) = map(k.k, map(k.f, X), map(k.f, X′))
_pairwise(k::ITKernel, X::AV) = pairwise(k.k, map(k.f, X))
_pairwise(k::ITKernel, X::AV, X′::AV) = pairwise(k.k, map(k.f, X), map(k.f, X′))

"""
    transform(f::Union{MeanFunction, Kernel}, ϕ)

Applies the input-transform `ϕ` to `f`.
"""
transform(μ::MeanFunction, ϕ) = ITMean(μ, ϕ)

transform(μ::MeanFunction, ::typeof(identity)) = μ
transform(μ::ZeroMean, ϕ) = μ
transform(μ::ZeroMean, ::typeof(identity)) = μ

transform(k::Kernel, ϕ) = ITKernel(k, ϕ)
transform(k::Kernel, ::typeof(identity)) = k
transform(k::ZeroKernel, ϕ) = k
transform(k::ZeroKernel, ::typeof(identity)) = k

transform(k::CrossKernel, ϕ, ::Val{1}) = LhsITCross(k, ϕ)
transform(k::CrossKernel, ϕ, ::Val{2}) = RhsITCross(k, ϕ)

transform(k::ZeroKernel, ϕ, ::Val{1}) = k
transform(k::ZeroKernel, ϕ, ::Val{2}) = k

transform(k::CrossKernel, ϕ, ϕ′) = ITCross(k, ϕ, ϕ′)

"""
    scale(f::Union{MeanFunction, Kernel}, l::Real)

Multiply each element of the input by `l`.
"""
scale(f::Union{MeanFunction, Kernel}, l::Real) = transform(f, Scale(l))

struct Scale{T<:Real}
    l::T
end
(s::Scale)(x) = s.l * x
map(s::Scale, x::AVM) = s.l .* x

"""
    pick_dims(x::Union{MeanFunction, Kernel}, I)

Returns either an `ITMean` or `ITKernel` which uses the columns of the input matrix `X`
specified by `I`.
"""
pick_dims(μ::MeanFunction, I) = ITMean(μ, X::AV->X[I])
pick_dims(k::Kernel, I) = ITKernel(k, X::AV->X[I])

"""
    periodic(k::Union{MeanFunction, Kernel}, θ::Real)

Make `k` periodic with period `f`.
"""
periodic(μ::MeanFunction, f::Real) = ITMean(μ, Periodic(f))
periodic(k::Kernel, f::Real) = ITKernel(k, Periodic(f))
periodic(k::EQ, p::Real) = PerEQ(p)
periodic(f::Union{MeanFunction, Kernel}) = periodic(f, 1.0)

struct Periodic{Tf<:Real}
    f::Tf
end
(p::Periodic)(t::Real) = [cos((2π * p.f) * t), sin((2π * p.f) * t)]
function map(p::Periodic, t::AV)
    return ColsAreObs(vcat(
        map(x->cos((2π * p.f) * x), t)',
        map(x->sin((2π * p.f) * x), t)',
    ))
end

"""
    Indexer

Pulls out the `n`th index of whatever object is passed.
"""
struct Indexer
    n::Int
end
(indexer::Indexer)(x) = x[indexer.n]
map(indexer::Indexer, x::ColsAreObs) = view(x.X, indexer.n, :)


text:
@model AHMCGaussian(y, N) = begin
    mu ~ Normal(0, 1)
    sigma ~ Truncated(Cauchy(0, 5), 0, Inf)
    y ~ MvNormal(Fill(mu, N), sigma)  # use `Fill(mu, N)` is memory freiendly
end

AHMCconfig = Turing.NUTS(1000, .85)

CmdStanGaussian = "
data {
  int<lower=0> N;
  vector[N] y;
}
parameters {
  real mu;
  real<lower=0> sigma;
}
model {
  mu ~ normal(0,1);
  sigma ~ cauchy(0,5);
  y ~ normal(mu,sigma);
}
"

CmdStanConfig = Stanmodel(
  name = "CmdStanGaussian", model=CmdStanGaussian, nchains=1,
  Sample(
    num_samples=1000, num_warmup=1000, adapt=CmdStan.Adapt(delta=0.8),
    save_warmup=true
  )
)

struct GaussianProb{TY <: AbstractVector}
    "Observations."
    y::TY
end

function (problem::GaussianProb)(θ)
    @unpack y = problem   # extract the data
    @unpack mu, sigma = θ
    N = length(y)
    logpdf(MvNormal(Fill(mu, N), sigma), y) + logpdf(Normal(0, 1), mu) +
      logpdf(Truncated(Cauchy(0, 5), 0, Inf), sigma)
end

# Define problem with data and inits.
function sampleDHMC(obs, N, nsamples, autodiff)
  p = GaussianProb(obs);
  p((mu=0.0, sigma=1.0))

  # Write a function to return properly dimensioned transformation.
  trans = as((mu=as(Real, -25, 25), sigma=asℝ₊), )

  # Use Flux for the gradient.

  P = TransformedLogDensity(trans, p)
  ∇P = ADgradient(autodiff, P)

  # Sample from the posterior.
  results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, nsamples; reporter = NoProgressReport())

  # Undo the transformation to obtain the posterior from the chain.
  posterior = transform.(trans, results.chain)

  chns = nptochain(results, posterior)
  return chns
end

simulateGaussian(; μ=0, σ=1, Nd, kwargs...) = (y=rand(Normal(μ, σ), Nd), N=Nd)


text:
export Dense

contract_t = 0.0

struct Dense{T} <: TensorStorage
  data::Vector{T}
  Dense{T}(data::Vector) where {T} = new{T}(convert(Vector{T},data))
  Dense{T}(size::Integer) where {T} = new{T}(Vector{T}(undef,size))
  Dense{T}(x::Number,size::Integer) where {T} = new{T}(fill(convert(T,x),size))
  Dense{T}() where {T} = new{T}(Vector{T}())
end

data(D::Dense) = D.data
length(D::Dense) = length(data(D))
eltype(D::Dense) = eltype(data(D))
getindex(D::Dense,i::Int) = data(D)[i]
#TODO: this should do proper promotions of the storage data
#e.g. ComplexF64*Dense{Float64} -> Dense{ComplexF64}
*(D::Dense{T},x::S) where {T,S<:Number} = Dense{promote_type(T,S)}(x*data(D))
*(x::Number,D::Dense) = D*x

convert(::Type{Dense{T}},D::Dense) where {T} = Dense{T}(data(D))

# convert to complex
storage_complex(D::Dense{T}) where {T} = Dense{complex(T)}(complex(data(D)))

copy(D::Dense{T}) where {T} = Dense{T}(copy(data(D)))

outer(D1::Dense{T},D2::Dense{S}) where {T, S <:Number} = Dense{promote_type(T,S)}(vec(data(D1)*transpose(data(D2))))

storage_convert(::Type{Array},D::Dense,is::IndexSet) = reshape(data(D),dims(is))

storage_fill!(D::Dense,x::Number) = fill!(data(D),x)

function storage_getindex(Tstore::Dense{T},
                          Tis::IndexSet,
                          vals::Union{Int, AbstractVector{Int}}...) where {T}
  return getindex(reshape(data(Tstore),dims(Tis)),vals...)
end

function storage_setindex!(Tstore::Dense,
                           Tis::IndexSet,
                           x::Union{<:Number, AbstractArray{<:Number}},
                           vals::Union{Int, AbstractVector{Int}}...)
  return setindex!(reshape(data(Tstore),dims(Tis)),x,vals...)
end

# For use in _add!
using Base.Cartesian: @nexprs,
                      @ntuple,
                      @nloops

#
# A generalized permutedims!(P,B,perm) that also allows
# a function to be applied elementwise (defaults to addition of
# the elements of P and the permuted elements of B)
#
# Based off of the permutedims! implementation in Julia's base:
# https://github.com/JuliaLang/julia/blob/91151ab871c7e7d6689d1cfa793c12062d37d6b6/base/multidimensional.jl#L1355
#
@generated function _add!(P::Array{T,N},
                          B::Array{S,N},
                          perm,
                          f = (x,y)->x+y) where {T,S,N}
  quote
    Base.checkdims_perm(P, B, perm)

    #calculates all the strides
    native_strides = Base.size_to_strides(1, size(B)...)
    strides_1 = 0
    @nexprs $N d->(strides_{d+1} = native_strides[perm[d]])

    #Creates offset, because indexing starts at 1
    offset = 1 - sum(@ntuple $N d->strides_{d+1})

    ind = 1
    @nexprs 1 d->(counts_{$N+1} = strides_{$N+1}) # a trick to set counts_($N+1)
    @nloops($N, i, P,
            d->(counts_d = strides_d), # PRE
            d->(counts_{d+1} += strides_{d+1}), # POST
            begin # BODY
                sumc = sum(@ntuple $N d->counts_{d+1})
                @inbounds P[ind] = f(P[ind],B[sumc+offset])
                ind += 1
            end)

    return P
  end
end

function add!(Bstore::Dense,
              Bis::IndexSet,
              Astore::Dense,
              Ais::IndexSet,
              x::Number = 1)
  p = calculate_permutation(Bis,Ais)
  Adata = data(Astore)
  Bdata = data(Bstore)
  if is_trivial_permutation(p)
    if x == 1
      Bdata .+= Adata
    else
      Bdata .= Bdata .+ x .* Adata
    end
  else
    reshapeBdata = reshape(Bdata,dims(Bis))
    reshapeAdata = reshape(Adata,dims(Ais))
    if x == 1
      _add!(reshapeBdata,reshapeAdata,p)
    else
      _add!(reshapeBdata,reshapeAdata,p,(a,b)->a+x*b)
    end
  end
end

function storage_add!(Bstore::Dense{BT},
                      Bis::IndexSet,
                      Astore::Dense{AT},
                      Ais::IndexSet,
                      x::Number = 1) where {BT,AT}
  NT = promote_type(AT,BT)
  if NT == BT
    add!(Bstore,Bis,Astore,Ais, x)
    return Bstore
  end
  Nstore = convert(Dense{NT},Bstore)
  add!(Nstore,Bis,Astore,Ais, x)
  return Nstore
end

function storage_copyto!(Bstore::Dense,
                         Bis::IndexSet,
                         Astore::Dense,
                         Ais::IndexSet,
                         x::Number = 1)
  p = calculate_permutation(Bis,Ais)
  Adata = data(Astore)
  Bdata = data(Bstore)
  if is_trivial_permutation(p)
    if x == 1
      Bdata .= Adata
    else
      Bdata .= x .* Adata
    end
  else
    reshapeBdata = reshape(Bdata,dims(Bis))
    reshapeAdata = reshape(Adata,dims(Ais))
    if x == 1
      _add!(reshapeBdata,reshapeAdata,p,(a,b)->b)
    else
      _add!(reshapeBdata,reshapeAdata,p,(a,b)->x*b)
    end
  end
end

function storage_mult!(Astore::Dense,
                       x::Number)
  Adata = data(Astore)
  rmul!(Adata, x)
end

function storage_mult(Astore::Dense,
                      x::Number)
  Bstore = copy(Astore)
  storage_mult!(Bstore, x)
  return Bstore
end

# For Real storage and complex scalar, promotion
# of the storage is needed
function storage_mult(Astore::Dense{T},
                      x::S) where {T<:Real,S<:Complex}
  Bstore = convert(Dense{promote_type(S,T)},Astore)
  storage_mult!(Bstore, x)
  return Bstore
end

# TODO: make this a special version of storage_add!()
# Make sure the permutation is optimized
function storage_permute!(Bstore::Dense,
                          Bis::IndexSet,
                          Astore::Dense,
                          Ais::IndexSet)
  p = calculate_permutation(Bis,Ais)
  Adata = data(Astore)
  Bdata = data(Bstore)
  if is_trivial_permutation(p)
    Bdata .= Adata
  else
    reshapeBdata = reshape(Bdata,dims(Bis))
    permutedims!(reshapeBdata,reshape(Adata,dims(Ais)),p)
  end
end

function storage_dag(Astore::Dense,Ais::IndexSet)
  return dag(Ais),storage_conj(Astore)
end

function storage_scalar(D::Dense)
  if length(D)==1
    return D[1]
  else
    throw(ErrorException("Cannot convert Dense -> Number for length of data greater than 1"))
  end
end

function is_outer(l1::Vector{Int},l2::Vector{Int})
  for l1i in l1
    if l1i < 0
      return false
    end
  end
  for l2i in l2
    if l2i < 0
      return false
    end
  end
  return true
end

# TODO: make this storage_contract!(), where C is pre-allocated.
#       This will allow for in-place multiplication
# TODO: optimize the contraction logic so C doesn't get permuted?
function storage_contract(Astore::TensorStorage,
                          Ais::IndexSet,
                          Bstore::TensorStorage,
                          Bis::IndexSet)
  if length(Ais)==0
    Cis = Bis
    Cstore = storage_scalar(Astore)*Bstore
  elseif length(Bis)==0
    Cis = Ais
    Cstore = storage_scalar(Bstore)*Astore
  else
    #TODO: check for special case when Ais and Bis are disjoint sets
    #I think we should do this analysis outside of storage_contract, at the ITensor level
    #(since it is universal for any storage type and just analyzes in indices)
    (Alabels,Blabels) = compute_contraction_labels(Ais,Bis)
    if is_outer(Alabels,Blabels)
      Cis = IndexSet(Ais,Bis)
      Cstore = outer(Astore,Bstore)
    else
      (Cis,Clabels) = contract_inds(Ais,Alabels,Bis,Blabels)
      Cstore = contract(Cis,Clabels,Astore,Ais,Alabels,Bstore,Bis,Blabels)
    end
  end
  return (Cis,Cstore)
end

function storage_svd(Astore::Dense{T},
                     Lis::IndexSet,
                     Ris::IndexSet;
                     kwargs...) where {T}
  maxdim::Int = get(kwargs,:maxdim,min(dim(Lis),dim(Ris)))
  mindim::Int = get(kwargs,:mindim,1)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  absoluteCutoff::Bool = get(kwargs,:absoluteCutoff,false)
  doRelCutoff::Bool = get(kwargs,:doRelCutoff,true)
  utags::String = get(kwargs,:utags,"Link,u")
  vtags::String = get(kwargs,:vtags,"Link,v")
  fastSVD::Bool = get(kwargs,:fastSVD,false)

  if fastSVD
    MU,MS,MV = svd(reshape(data(Astore),dim(Lis),dim(Ris)))
  else
    MU,MS,MV = recursiveSVD(reshape(data(Astore),dim(Lis),dim(Ris)))
  end
  MV = conj!(MV)

  P = MS.^2
  #@printf "  Truncating with maxdim=%d cutoff=%.3E\n" maxdim cutoff
  truncate!(P;mindim=mindim,
              maxdim=maxdim,
              cutoff=cutoff,
              absoluteCutoff=absoluteCutoff,
              doRelCutoff=doRelCutoff)
  dS = length(P)
  if dS < length(MS)
    MU = MU[:,1:dS]
    resize!(MS,dS)
    MV = MV[:,1:dS]
  end

  u = Index(dS,utags)
  v = settags(u,vtags)
  Uis,Ustore = IndexSet(Lis...,u),Dense{T}(vec(MU))
  #TODO: make a diag storage
  Sis,Sstore = IndexSet(u,v),Dense{Float64}(vec(Matrix(Diagonal(MS))))
  Vis,Vstore = IndexSet(Ris...,v),Dense{T}(Vector{T}(vec(MV)))

  return (Uis,Ustore,Sis,Sstore,Vis,Vstore)
end

function storage_eigen(Astore::Dense{T},
                       Lis::IndexSet,
                       Ris::IndexSet;
                       kwargs...) where {T}
  maxdim::Int = get(kwargs,:maxdim,min(dim(Lis),dim(Ris)))
  mindim::Int = get(kwargs,:mindim,1)
  cutoff::Float64 = get(kwargs,:cutoff,0.0)
  absoluteCutoff::Bool = get(kwargs,:absoluteCutoff,false)
  doRelCutoff::Bool = get(kwargs,:doRelCutoff,true)
  tags::TagSet = get(kwargs,:lefttags,"Link,u")
  lefttags::TagSet = get(kwargs,:lefttags,tags)
  righttags::TagSet = get(kwargs,:righttags,prime(lefttags))

  dim_left = dim(Lis)
  dim_right = dim(Ris)
  MD,MU = eigen(Hermitian(reshape(data(Astore),dim_left,dim_right)))

  # Sort by largest to smallest eigenvalues
  p = sortperm(MD; rev = true)
  MD = MD[p]
  MU = MU[:,p]

  #@printf "  Truncating with maxdim=%d cutoff=%.3E\n" maxdim cutoff
  truncate!(MD;maxdim=maxdim,
              cutoff=cutoff,
              absoluteCutoff=absoluteCutoff,
              doRelCutoff=doRelCutoff)
  dD = length(MD)
  if dD < size(MU,2)
    MU = MU[:,1:dD]
  end

  #TODO: include truncation parameters as keyword arguments
  u = Index(dD,lefttags)
  v = settags(u,righttags)
  Uis,Ustore = IndexSet(Lis...,u),Dense{T}(vec(MU))
  #TODO: make a diag storage
  Dis,Dstore = IndexSet(u,v),Dense{T}(vec(Matrix(Diagonal(MD))))
  return (Uis,Ustore,Dis,Dstore)
end

function polar(A::Matrix)
  U,S,V = svd(A) # calls LinearAlgebra.svd()
  return U*V',V*Diagonal(S)*V'
end

function storage_qr(Astore::Dense{T},
                    Lis::IndexSet,
                    Ris::IndexSet;
                    kwargs...) where {T}
  tags::TagSet = get(kwargs,:tags,"Link,u")
  dim_left = dim(Lis)
  dim_right = dim(Ris)
  MQ,MP = qr(reshape(data(Astore),dim_left,dim_right))
  dim_middle = min(dim_left,dim_right)
  u = Index(dim_middle,tags)
  #Must call Matrix() on MQ since the QR decomposition outputs a sparse
  #form of the decomposition
  Qis,Qstore = IndexSet(Lis...,u),Dense{T}(vec(Matrix(MQ)))
  Pis,Pstore = IndexSet(u,Ris...),Dense{T}(vec(Matrix(MP)))
  return (Qis,Qstore,Pis,Pstore)
end

function storage_polar(Astore::Dense{T},
                       Lis::IndexSet,
                       Ris::IndexSet) where {T}
  dim_left = dim(Lis)
  dim_right = dim(Ris)
  MQ,MP = polar(reshape(data(Astore),dim_left,dim_right))
  dim_middle = min(dim_left,dim_right)
  Uis = prime(Ris)
  Qis,Qstore = IndexSet(Lis...,Uis...),Dense{T}(vec(MQ))
  Pis,Pstore = IndexSet(Uis...,Ris...),Dense{T}(vec(MP))
  return (Qis,Qstore,Pis,Pstore)
end


text:
# Notes (E - everywhere/every worker, L - local/this worker)
# [E] @everywhere using MonteCarlo (this is blocking)
# [E] => creates constants
# [L] create MC, updates, scheduler
# [L] => scheduler generates get_conf
# [L] => notify connected workers (write to constants)
# [L] yield to recieve notfications
# [L] run simulation, yield once per sweep

# Look through these to find a partner
const connected_ids = Vector{Int}(undef, 0)

# process local
add_worker!(id) = begin in(id, connected_ids) || push!(connected_ids, id); true end
remove_worker!(id) = begin filter!(x -> x != id, connected_ids); true end
# this should be replaced when a parallel update is created
get_conf(args...) = false

# multi-process communication
# don't care if they return, we're just saying "call me maybe?"
connect(targets...; wait = false) = connect(collect(targets), wait = wait)
function connect(targets::Vector; wait = false)
    if wait
        @sync for target in targets
            @async remotecall_wait(add_worker!, target, myid())
        end
    else
        for target in targets
            remotecall(add_worker!, target, myid())
        end
    end
    nothing
end
# here we do care, because otherwise remote might be requesting from us when we're gone
disconnect(targets...; wait = false) = disconnect(collect(targets), wait = wait)
function disconnect(targets; wait=false)
    if wait
        @sync for target in targets
            @async remotecall_wait(remove_worker!, target, myid())
        end
    else
        for target in targets
            remotecall(remove_worker!, target, myid())
        end
    end
    nothing
end
pull_conf_from_remote(trg) = remotecall_fetch(get_conf, trg)

const weight_probability = Channel{Tuple{Int, Float64, Float64}}(1)
function put_weight_prob!(id, weight, probability)
    put!(weight_probability, (id, weight, probability))
    nothing
end

# 2 process barrier
const barrier_counter = Ref(1)
const waiting_list = Set{Tuple{Int, Int}}()

function wait_for(id, counter)
    push!(waiting_list, (id, counter))
    filter!(x -> x[2] >= barrier_counter[], waiting_list)
    nothing
end

function wait_for_remote(ids::Int...; timeout=Inf)
    t0 = time()
    for trg in ids
        remotecall(wait_for, trg, myid(), barrier_counter[])
    end

    while !all(trg -> (trg, barrier_counter[]) in waiting_list, ids)
        yield()
        if time() - t0 > timeout
            barrier_counter[] += 1
            return false
        end
    end

    barrier_counter[] += 1
    return true
end

function generate_communication_functions(local_conf)
    @eval MonteCarlo begin
        get_conf() = $local_conf
    end
end

################################################################################
### Updates
################################################################################

abstract type AbstractParallelUpdate <: AbstractGlobalUpdate end
init!(mc, ::AbstractParallelUpdate) = generate_communication_functions(conf(field(mc)))

"""
    ReplicaExchange([mc, model], target[, timeout = 600.0])

Represents a replica exchange update with a given `target` worker.

Note that this update is blocking, i.e. a simulation will wait on its partner to
perform the update. You can define a `timeout` in seconds to avoid the simulation
getting stuck for too long.
"""
struct ReplicaExchange <: AbstractParallelUpdate
    target::Int
    timeout::Float64
end
ReplicaExchange(target) = ReplicaExchange(target, 600.0)
ReplicaExchange(mc, model, target, timeout=600.0) = ReplicaExchange(target, timeout)
name(::ReplicaExchange) = "ReplicaExchange"

@bm function update(u::ReplicaExchange, mc, model, field)
    tc = temp_conf(field)
    # Need to sync at the start here because else th weights might be based on different confs
    # barrier
    if !wait_for_remote(u.target, timeout = u.timeout)
        return 0
    end

    # swap conf
    conf = pull_conf_from_remote(u.target)
    tc .= conf

    # compute weight
    detratio, ΔE_boson, passthrough = propose_global_from_conf(mc, model, tc)
    local_weight = exp(- ΔE_boson) * detratio
    local_prob = rand()

    # swap weight, probability
    remotecall(put_weight_prob!, u.target, myid(), local_weight, local_prob)
    while !isready(weight_probability)
        yield()
    end
    remote, remote_weight, remote_prob = take!(weight_probability)
    @assert remote == u.target

    # accept/deny
    # We need to pick one of the random values and do so consistently across
    # processes.
    w = local_weight * remote_weight
    if ifelse(myid() < u.target, local_prob, remote_prob) < w
        accept_global!(mc, model, tc, passthrough)
        return 1
    end

    return 0
end

"""
    ReplicaPull([mc, model])

This update will pull a configuration from a connected worker and attempt a
global update with it. The target worker will be cycled as the simulation
progresses.

To connect workers for this update, the remote worker has to call
`connect(this_worker)`.

This update only blocks locally.
"""
mutable struct ReplicaPull <: AbstractParallelUpdate
    cycle_idx::Int
end
ReplicaPull() = ReplicaPull(1)
ReplicaPull(mc::MonteCarloFlavor, model::Model) = ReplicaPull(1)
name(::ReplicaPull) = "ReplicaPull"
Base.:(==)(a::ReplicaPull, b::ReplicaPull) = a.cycle_idx == b.cycle_idx

@bm function update(u::ReplicaPull, mc, model, field)
    tc = temp_conf(field)
    # cycle first to make sure the idx is in bounds
    @sync if !isempty(connected_ids)
        idx = mod1(u.cycle_idx, length(connected_ids))
        conf = pull_conf_from_remote(connected_ids[idx])
        tc .= conf
        return global_update(mc, model, tc)
    end
    return 0
end


text:
model = Particle()
withloss = WithLossDo(x-> nothing)
withtraininglosses = WithTrainingLossesDo(x->nothing)
step = Step(1)

@test_throws IC.ERR_TRAIN IterationControl.train!(model)
@test_throws IC.err_train(model) IterationControl.train!(model, 1)

# lifting train!:
IC.train!(model::Particle, n) = train!(model, n)

@test_throws(IC.ERR_NEEDS_LOSS,
             IC.train!(model, step, withloss, NumberLimit(1)))

# lifting loss!:
IterationControl.loss(m::Particle) = loss(m)

@test_logs IC.train!(model, step, withloss, NumberLimit(1), verbosity=0)

@test_throws(IC.ERR_NEEDS_TRAINING_LOSSES,
             IC.train!(model, step, PQ(), NumberLimit(1)))

@test_throws(IC.ERR_NEEDS_TRAINING_LOSSES,
             IC.train!(model, step, withtraininglosses, NumberLimit(1)))

# lifting training_losses:
IterationControl.training_losses(m::Particle) = training_losses(m)

@test_logs IC.train!(model, Step(2), PQ(), NumberLimit(1), verbosity=0)
@test_logs IC.train!(model, Step(2), withtraininglosses,
                     NumberLimit(1), verbosity=0)

@test_throws(IC.err_ingest(model),
             IC.train!(model, Data(1:2), NumberLimit(1)))

#lifting ingest!:
IC.ingest!(model::Particle, datum) = ingest!(model, datum)

@test_logs IC.train!(model, Step(2), Data(1:1), NumberLimit(1), verbosity=0)


text:
# Chemfiles.jl, a modern library for chemistry file reading and writing
# Copyright (C) Guillaume Fraux and contributors -- BSD license

export add_atom!, remove_atom!
export has_velocities, add_velocities!, positions, velocities
export set_cell!, set_topology!, set_step!, guess_bonds!
export distance, dihedral, out_of_plane

__ptr(frame::Frame) = __ptr(frame.__handle)
__const_ptr(frame::Frame) = __const_ptr(frame.__handle)

"""
Create a new empty ``Frame``.
"""
function Frame()
    ptr = @__check_ptr(lib.chfl_frame())
    return Frame(CxxPointer(ptr, is_const=false))
end

"""
Get the number of atoms in the ``frame``.
"""
function Base.size(frame::Frame)
    count = Ref{UInt64}(0)
    __check(lib.chfl_frame_atoms_count(__const_ptr(frame), count))
    return count[]
end

"""
Get the number of atoms in the ``frame``.
"""
function Base.length(frame::Frame)
    size(frame)
end

"""
Resize the positions and the velocities in the ``frame``, to make space for
``natoms`` atoms. This function may invalidate any pointer to the positions or
the velocities if the new size is bigger than the old one. In all the cases,
previous data is conserved. This function conserve the presence or absence of
velocities.
"""
function Base.resize!(frame::Frame, natoms::Integer)
    __check(lib.chfl_frame_resize(__ptr(frame), UInt64(natoms)))
end

"""
Get the positions in a ``Frame`` as an array. The positions are readable and
writable from this array. If the frame is resized (by writing to it, or calling
``resize!``), the array is invalidated.
"""
function positions(frame::Frame)
    ptr = Ref{Ptr{Float64}}()
    natoms = Ref{UInt64}(0)
    __check(lib.chfl_frame_positions(__ptr(frame), ptr, natoms))
    return unsafe_wrap(Array{Float64, 2}, ptr[], (3, Int(natoms[])); own=false)
end

"""
Get the velocities in a ``Frame`` as an array. The velocities are readable and
writable from this array. If the frame is resized (by writing to it, or calling
``resize!``), the array is invalidated.

If the frame do not have velocity, this function will error. You can use
``add_velocities!`` to add velocities to a frame before calling this function.
"""
function velocities(frame::Frame)
    ptr = Ref{Ptr{Float64}}()
    natoms = Ref{UInt64}(0)
    __check(lib.chfl_frame_velocities(__ptr(frame), ptr, natoms))
    return unsafe_wrap(Array{Float64, 2}, ptr[], (3, Int(natoms[])); own=false)
end

"""
Add velocities to this ``frame``. The storage is initialized with the result of
``size(frame)`` as number of atoms. If the frame already have velocities, this
does nothing.
"""
function add_velocities!(frame::Frame)
    __check(lib.chfl_frame_add_velocities(__ptr(frame)))
    return nothing
end

"""
Check if a ``frame`` contains velocity data or not.
"""
function has_velocities(frame::Frame)
    result = Ref{UInt8}(0)
    __check(lib.chfl_frame_has_velocities(__const_ptr(frame), result))
    return convert(Bool, result[])
end

"""
Set the ``cell`` associated with a ``frame``.
"""
function set_cell!(frame::Frame, cell::UnitCell)
    __check(lib.chfl_frame_set_cell(__ptr(frame), __const_ptr(cell)))
    return nothing
end

"""
Set the ``topology`` associated with a ``frame``.
"""
function set_topology!(frame::Frame, topology::Topology)
    __check(lib.chfl_frame_set_topology(__ptr(frame), __const_ptr(topology)))
    return nothing
end

"""
Get the ``frame`` step, *i.e.* the frame number in the trajectory.
"""
function Base.step(frame::Frame)
    result = Ref{UInt64}(0)
    __check(lib.chfl_frame_step(__const_ptr(frame), result))
    return result[]
end

"""
Set the ``frame`` step to ``step``.
"""
function set_step!(frame::Frame, step::Integer)
    __check(lib.chfl_frame_set_step(__ptr(frame), UInt64(step)))
    return nothing
end

"""
Guess the bonds, angles and dihedrals in the ``frame`` using a distance criteria.
"""
function guess_bonds!(frame::Frame)
    __check(lib.chfl_frame_guess_bonds(__ptr(frame)))
    return nothing
end

"""
Calculate the distance between two atoms.
"""
function distance(frame::Frame, i::Integer, j::Integer)
    result = Ref{Float64}(0)
    __check(lib.chfl_frame_distance(__const_ptr(frame), UInt64(i), UInt64(j), result))
    return result[]
end

"""
Calculate the angle made by three atoms.
"""
function Base.angle(frame::Frame, i::Integer, j::Integer, k::Integer)
    result = Ref{Float64}(0)
    __check(lib.chfl_frame_angle(__const_ptr(frame), UInt64(i), UInt64(j), UInt64(k), result))
    return result[]
end

"""
Calculate the dihedral (torsional) angle made by four unbranched atoms.
"""
function dihedral(frame::Frame, i::Integer, j::Integer, k::Integer, m::Integer)
    result = Ref{Float64}(0)
    __check(lib.chfl_frame_dihedral(__const_ptr(frame), UInt64(i), UInt64(j), UInt64(k), UInt64(m), result))
    return result[]
end

"""
Calculate the out-of-plane (improper) angle made by four atoms.
"""
function out_of_plane(frame::Frame, i::Integer, j::Integer, k::Integer, m::Integer)
    result = Ref{Float64}(0)
    __check(lib.chfl_frame_out_of_plane(__const_ptr(frame), UInt64(i), UInt64(j), UInt64(k), UInt64(m), result))
    return result[]
end

"""
Add an ``atom`` and the corresponding ``position`` and ``velocity`` data to a
``frame``.
"""
function add_atom!(frame::Frame, atom::Atom, position::Vector{Float64}, velocity::Vector{Float64} = Float64[0.0,0.0,0.0])
    __check(lib.chfl_frame_add_atom(__ptr(frame), __const_ptr(atom), position, velocity))
    return nothing
end

"""
Remove the ``atom`` at ``index`` from the ``frame``.

This modify all the ``atoms`` indexes after ``index``, and invalidate any
array obtained using ``positions`` or ``velocities``.
"""
function remove_atom!(frame::Frame, index::Integer)
    __check(lib.chfl_frame_remove(__ptr(frame), UInt64(index)))
    return nothing
end

"""
Set a named property for the given ``Frame``.
"""
function set_property!(frame::Frame, name::String, value)
    property = Property(value)
    __check(lib.chfl_frame_set_property(
        __ptr(frame), pointer(name), __const_ptr(property)
    ))
    return nothing
end

"""
Get a named property for the given atom.
"""
function property(frame::Frame, name::String)
    ptr = lib.chfl_frame_get_property(__const_ptr(frame), pointer(name))
    return extract(Property(CxxPointer(ptr, is_const=false)))
end

"""
Get the number of properties associated with a frame.
"""
function properties_count(frame::Frame)
    count = Ref{UInt64}(0)
    __check(lib.chfl_frame_properties_count(__const_ptr(frame), count))
    return count[]
end

"""
Get the names of all properties associated with a frame.
"""
function list_properties(frame::Frame)
    count = properties_count(frame)
    names = Array{Ptr{UInt8}}(undef, count)
    __check(lib.chfl_frame_list_properties(__const_ptr(frame), pointer(names), count))
    return map(unsafe_string, names)
end

"""
Add an additional bond to the ``Frame``'s ``Topology``.
"""
function add_bond!(frame::Frame, i::Integer, j::Integer, order=nothing)
    if order == nothing
        __check(lib.chfl_frame_add_bond(__ptr(frame), UInt64(i), UInt64(j)))
    else
        # Check that the order is a valid BondOrder
        order = BondOrder(Integer(order))
        __check(lib.chfl_frame_bond_with_order(
            __ptr(frame), UInt64(i), UInt64(j), lib.chfl_bond_order(order)
        ))
    end
    return nothing
end

"""
Remove a bond from the ``Frame``'s ``Topology``.
"""
function remove_bond!(frame::Frame, i::Integer, j::Integer)
    __check(lib.chfl_frame_remove_bond(__ptr(frame), UInt64(i), UInt64(j)))
    return nothing
end

"""
Add a residue to the ``Frame``'s ``Topology``.
"""
function add_residue!(frame::Frame, residue::Residue)
    __check(lib.chfl_frame_add_residue(__ptr(frame), __const_ptr(residue)))
    return nothing
end

"""
Make a deep copy of a ``Frame``.
"""
function Base.deepcopy(frame::Frame)
    ptr = lib.chfl_frame_copy(__const_ptr(frame))
    return Frame(CxxPointer(ptr, is_const=false))
end

# Iteration support
Base.iterate(frame::Frame, state=0) = state >= size(frame) ? nothing : (Atom(frame, state), state+1)
Base.eltype(frame::Frame) = Atom


text:
# Populate a SeisData structure with random channels
"""
    randSeisData([, c=0.2::Float64, s=0.6::Float64, nx=L::Integer])

Generate 8 to 24 channels of random seismic data as a SeisData object.

    randSeisData(N, [, c=0.2::Float64, s=0.6::Float64, nx=L::Integer])

Generate `N` channels of random seismic data as a SeisData object.

### Notes and Keyword Behavior
* `nx=N` sets the length of each channel to `N` samples.
* 100*s is the minimum percentage of channels with guaranteed seismic data
* 100*c is the maximum percentage of channels with irregularly-sampled data
* `s` takes precedence over `c`; they are not renormalized. Thus, for example:
  + randSeisData(12, c=1.0, s=1.0) works like randSeisData(12, c=0.0, s=1.0)
  + randSeisData(c=1.0) works like randSeisData(s=0.0)
* For `N` channels, `s*N` is rounded up; `c*N` is rounded down. Thus:
  + randSeisData(10, c=0.28) and randSeisData(10, c=0.2) are equivalent
  + randSeisData(10, s=0.28) and randSeisData(10, s=0.3) are equivalent
"""
function randSeisData(N::Int64; c::Float64=0.2, s::Float64=0.6, nx::Integer=0)
  S = SeisData(N)

  # Evaluate probabilities, with s taking precedence
  if s != 0.6
    c = max(0.0, 1.0-s)
  end
  if c != 0.2
    s = max(0.0, 1.0-c)
  end

  # determine number of channels of each
  n_seis = max(min(ceil(Int, s*S.n), S.n), 0)
  n_irr = max(min(floor(Int, c*S.n), S.n-n_seis), 0)

  # fill in and shuffle data_spec
  data_spec = zeros(UInt8, S.n)
  data_spec[1:n_seis] .= 0x01
  data_spec[n_seis+1:n_seis+n_irr] .= 0x02
  shuffle!(data_spec)

  # populate all channels according to data_spec
  for i = 1:S.n
    if data_spec[i] == 0x01
      S[i] = randSeisChannel(s=true, nx=nx)
    elseif data_spec[i] == 0x02
      S[i] = randSeisChannel(c=true, nx=nx)
    else
      S[i] = randSeisChannel(nx=nx)
    end
  end
  return S
end
randSeisData(; c::Float64=0.2, s::Float64=0.6, nx::Int64=0) = randSeisData(rand(8:24), c=c, s=s, nx=nx)


text:
"""
Support of vectors of random variables, i.e vectors of Distributions
"""
module DistributionVectors

export AbstractDistributionVector, SimpleDistributionVector, ParamDistributionVector,
    vectuptotupvec,
    AutoCorrelationFunction, cormatrix_for_acf

using Distributions, Missings, Statistics
using MissingStrategies
using StatsBase, Random
using FillArrays, RecursiveArrayTools, BandedMatrices, LinearAlgebra, OffsetArrays

import Base: size, length, IndexStyle, similar, getindex, setindex!, sum, ismissing
import Random: rand, rand!
import Statistics: mean
import StatsBase: params, coef

include("distributionvector.jl")
include("autocor.jl")
include("sumdistributionvector.jl")

# specific distributions
include("univariates.jl")

end # module


text:
#####
##### Prognostic Variable types
#####

struct VolumetricLiquidFraction <: AbstractPrognosticVariable end
struct VolumetricInternalEnergy <: AbstractPrognosticVariable end
struct VolumetricIceFraction <: AbstractPrognosticVariable end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""MutatingWebhook describes an admission webhook and the resources and operations it applies to.

    IoK8sApiAdmissionregistrationV1MutatingWebhook(;
        admissionReviewVersions=nothing,
        clientConfig=nothing,
        failurePolicy=nothing,
        matchPolicy=nothing,
        name=nothing,
        namespaceSelector=nothing,
        objectSelector=nothing,
        reinvocationPolicy=nothing,
        rules=nothing,
        sideEffects=nothing,
        timeoutSeconds=nothing,
    )

    - admissionReviewVersions::Vector{String} : AdmissionReviewVersions is an ordered list of preferred `AdmissionReview` versions the Webhook expects. API server will try to use first version in the list which it supports. If none of the versions specified in this list supported by API server, validation will fail for this object. If a persisted webhook configuration specifies allowed versions and does not include any versions known to the API Server, calls to the webhook will fail and be subject to the failure policy.
    - clientConfig::IoK8sApiAdmissionregistrationV1WebhookClientConfig : ClientConfig defines how to communicate with the hook. Required
    - failurePolicy::String : FailurePolicy defines how unrecognized errors from the admission endpoint are handled - allowed values are Ignore or Fail. Defaults to Fail.
    - matchPolicy::String : matchPolicy defines how the \"rules\" list is used to match incoming requests. Allowed values are \"Exact\" or \"Equivalent\".  - Exact: match a request only if it exactly matches a specified rule. For example, if deployments can be modified via apps/v1, apps/v1beta1, and extensions/v1beta1, but \"rules\" only included `apiGroups:[\"apps\"], apiVersions:[\"v1\"], resources: [\"deployments\"]`, a request to apps/v1beta1 or extensions/v1beta1 would not be sent to the webhook.  - Equivalent: match a request if modifies a resource listed in rules, even via another API group or version. For example, if deployments can be modified via apps/v1, apps/v1beta1, and extensions/v1beta1, and \"rules\" only included `apiGroups:[\"apps\"], apiVersions:[\"v1\"], resources: [\"deployments\"]`, a request to apps/v1beta1 or extensions/v1beta1 would be converted to apps/v1 and sent to the webhook.  Defaults to \"Equivalent\"
    - name::String : The name of the admission webhook. Name should be fully qualified, e.g., imagepolicy.kubernetes.io, where \"imagepolicy\" is the name of the webhook, and kubernetes.io is the name of the organization. Required.
    - namespaceSelector::IoK8sApimachineryPkgApisMetaV1LabelSelector : NamespaceSelector decides whether to run the webhook on an object based on whether the namespace for that object matches the selector. If the object itself is a namespace, the matching is performed on object.metadata.labels. If the object is another cluster scoped resource, it never skips the webhook.  For example, to run the webhook on any objects whose namespace is not associated with \"runlevel\" of \"0\" or \"1\";  you will set the selector as follows: \"namespaceSelector\": {   \"matchExpressions\": [     {       \"key\": \"runlevel\",       \"operator\": \"NotIn\",       \"values\": [         \"0\",         \"1\"       ]     }   ] }  If instead you want to only run the webhook on any objects whose namespace is associated with the \"environment\" of \"prod\" or \"staging\"; you will set the selector as follows: \"namespaceSelector\": {   \"matchExpressions\": [     {       \"key\": \"environment\",       \"operator\": \"In\",       \"values\": [         \"prod\",         \"staging\"       ]     }   ] }  See https://kubernetes.io/docs/concepts/overview/working-with-objects/labels/ for more examples of label selectors.  Default to the empty LabelSelector, which matches everything.
    - objectSelector::IoK8sApimachineryPkgApisMetaV1LabelSelector : ObjectSelector decides whether to run the webhook based on if the object has matching labels. objectSelector is evaluated against both the oldObject and newObject that would be sent to the webhook, and is considered to match if either object matches the selector. A null object (oldObject in the case of create, or newObject in the case of delete) or an object that cannot have labels (like a DeploymentRollback or a PodProxyOptions object) is not considered to match. Use the object selector only if the webhook is opt-in, because end users may skip the admission webhook by setting the labels. Default to the empty LabelSelector, which matches everything.
    - reinvocationPolicy::String : reinvocationPolicy indicates whether this webhook should be called multiple times as part of a single admission evaluation. Allowed values are \"Never\" and \"IfNeeded\".  Never: the webhook will not be called more than once in a single admission evaluation.  IfNeeded: the webhook will be called at least one additional time as part of the admission evaluation if the object being admitted is modified by other admission plugins after the initial webhook call. Webhooks that specify this option *must* be idempotent, able to process objects they previously admitted. Note: * the number of additional invocations is not guaranteed to be exactly one. * if additional invocations result in further modifications to the object, webhooks are not guaranteed to be invoked again. * webhooks that use this option may be reordered to minimize the number of additional invocations. * to validate an object after all mutations are guaranteed complete, use a validating admission webhook instead.  Defaults to \"Never\".
    - rules::Vector{IoK8sApiAdmissionregistrationV1RuleWithOperations} : Rules describes what operations on what resources/subresources the webhook cares about. The webhook cares about an operation if it matches _any_ Rule. However, in order to prevent ValidatingAdmissionWebhooks and MutatingAdmissionWebhooks from putting the cluster in a state which cannot be recovered from without completely disabling the plugin, ValidatingAdmissionWebhooks and MutatingAdmissionWebhooks are never called on admission requests for ValidatingWebhookConfiguration and MutatingWebhookConfiguration objects.
    - sideEffects::String : SideEffects states whether this webhook has side effects. Acceptable values are: None, NoneOnDryRun (webhooks created via v1beta1 may also specify Some or Unknown). Webhooks with side effects MUST implement a reconciliation system, since a request may be rejected by a future step in the admission change and the side effects therefore need to be undone. Requests with the dryRun attribute will be auto-rejected if they match a webhook with sideEffects == Unknown or Some.
    - timeoutSeconds::Int32 : TimeoutSeconds specifies the timeout for this webhook. After the timeout passes, the webhook call will be ignored or the API call will fail based on the failure policy. The timeout value must be between 1 and 30 seconds. Default to 10 seconds.
"""
mutable struct IoK8sApiAdmissionregistrationV1MutatingWebhook <: SwaggerModel
    admissionReviewVersions::Any # spec type: Union{ Nothing, Vector{String} } # spec name: admissionReviewVersions
    clientConfig::Any # spec type: Union{ Nothing, IoK8sApiAdmissionregistrationV1WebhookClientConfig } # spec name: clientConfig
    failurePolicy::Any # spec type: Union{ Nothing, String } # spec name: failurePolicy
    matchPolicy::Any # spec type: Union{ Nothing, String } # spec name: matchPolicy
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    namespaceSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: namespaceSelector
    objectSelector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: objectSelector
    reinvocationPolicy::Any # spec type: Union{ Nothing, String } # spec name: reinvocationPolicy
    rules::Any # spec type: Union{ Nothing, Vector{IoK8sApiAdmissionregistrationV1RuleWithOperations} } # spec name: rules
    sideEffects::Any # spec type: Union{ Nothing, String } # spec name: sideEffects
    timeoutSeconds::Any # spec type: Union{ Nothing, Int32 } # spec name: timeoutSeconds

    function IoK8sApiAdmissionregistrationV1MutatingWebhook(;admissionReviewVersions=nothing, clientConfig=nothing, failurePolicy=nothing, matchPolicy=nothing, name=nothing, namespaceSelector=nothing, objectSelector=nothing, reinvocationPolicy=nothing, rules=nothing, sideEffects=nothing, timeoutSeconds=nothing)
        o = new()
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("admissionReviewVersions"), admissionReviewVersions)
        setfield!(o, Symbol("admissionReviewVersions"), admissionReviewVersions)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("clientConfig"), clientConfig)
        setfield!(o, Symbol("clientConfig"), clientConfig)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("failurePolicy"), failurePolicy)
        setfield!(o, Symbol("failurePolicy"), failurePolicy)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("matchPolicy"), matchPolicy)
        setfield!(o, Symbol("matchPolicy"), matchPolicy)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("namespaceSelector"), namespaceSelector)
        setfield!(o, Symbol("namespaceSelector"), namespaceSelector)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("objectSelector"), objectSelector)
        setfield!(o, Symbol("objectSelector"), objectSelector)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("reinvocationPolicy"), reinvocationPolicy)
        setfield!(o, Symbol("reinvocationPolicy"), reinvocationPolicy)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("rules"), rules)
        setfield!(o, Symbol("rules"), rules)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("sideEffects"), sideEffects)
        setfield!(o, Symbol("sideEffects"), sideEffects)
        validate_property(IoK8sApiAdmissionregistrationV1MutatingWebhook, Symbol("timeoutSeconds"), timeoutSeconds)
        setfield!(o, Symbol("timeoutSeconds"), timeoutSeconds)
        o
    end
end # type IoK8sApiAdmissionregistrationV1MutatingWebhook

const _property_map_IoK8sApiAdmissionregistrationV1MutatingWebhook = Dict{Symbol,Symbol}(Symbol("admissionReviewVersions")=>Symbol("admissionReviewVersions"), Symbol("clientConfig")=>Symbol("clientConfig"), Symbol("failurePolicy")=>Symbol("failurePolicy"), Symbol("matchPolicy")=>Symbol("matchPolicy"), Symbol("name")=>Symbol("name"), Symbol("namespaceSelector")=>Symbol("namespaceSelector"), Symbol("objectSelector")=>Symbol("objectSelector"), Symbol("reinvocationPolicy")=>Symbol("reinvocationPolicy"), Symbol("rules")=>Symbol("rules"), Symbol("sideEffects")=>Symbol("sideEffects"), Symbol("timeoutSeconds")=>Symbol("timeoutSeconds"))
const _property_types_IoK8sApiAdmissionregistrationV1MutatingWebhook = Dict{Symbol,String}(Symbol("admissionReviewVersions")=>"Vector{String}", Symbol("clientConfig")=>"IoK8sApiAdmissionregistrationV1WebhookClientConfig", Symbol("failurePolicy")=>"String", Symbol("matchPolicy")=>"String", Symbol("name")=>"String", Symbol("namespaceSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("objectSelector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("reinvocationPolicy")=>"String", Symbol("rules")=>"Vector{IoK8sApiAdmissionregistrationV1RuleWithOperations}", Symbol("sideEffects")=>"String", Symbol("timeoutSeconds")=>"Int32")
Base.propertynames(::Type{ IoK8sApiAdmissionregistrationV1MutatingWebhook }) = collect(keys(_property_map_IoK8sApiAdmissionregistrationV1MutatingWebhook))
Swagger.property_type(::Type{ IoK8sApiAdmissionregistrationV1MutatingWebhook }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiAdmissionregistrationV1MutatingWebhook[name]))}
Swagger.field_name(::Type{ IoK8sApiAdmissionregistrationV1MutatingWebhook }, property_name::Symbol) =  _property_map_IoK8sApiAdmissionregistrationV1MutatingWebhook[property_name]

function check_required(o::IoK8sApiAdmissionregistrationV1MutatingWebhook)
    (getproperty(o, Symbol("admissionReviewVersions")) === nothing) && (return false)
    (getproperty(o, Symbol("clientConfig")) === nothing) && (return false)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("sideEffects")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiAdmissionregistrationV1MutatingWebhook }, name::Symbol, val)
end


text:
mutable struct DisplayParameters
    format::Symbol
    decorations::Bool
    sigfigs::Int
end

const display_params = DisplayParameters(:standard, false, 6)

doc"""
    setformat(;kw)

`setformat` changes how intervals are displayed using keyword arguments.
The following options are available:

- `format`: interval output format

    - `:standard`: `[1, 2]`
    - `:full`: `Interval(1, 2)`
    - `:midpoint`: 1.5 ± 0.5

- `sigfigs`: number of significant figures to show in `standard` mode; the default is 6

- `decorations` (boolean):  show decorations or not

Example:
```
julia> setformat(:full, decorations=true)
```
"""
function setformat(format = display_params.format;
                    decorations = display_params.decorations, sigfigs::Integer = display_params.sigfigs)

    if format ∉ (:standard, :full, :midpoint)
        throw(ArgumentError("Allowed format option is one of  $display_options."))
    end

    if decorations ∉ (true, false)
        throw(ArgumentError("`decorations` must be `true` or `false`"))
    end

    if sigfigs < 1
        throw(ArgumentError("`sigfigs` must be `>= 1`"))
    end

    # update values in display_params:
    display_params.format = format
    display_params.decorations = decorations
    display_params.sigfigs = sigfigs
end

doc"""
    @format [style::Symbol] [decorations::Bool] [sigfigs::Integer]

The `@format` macro provides a simple interface to control the output format
for intervals. These options are passed to the `setformat` function.

The arguments may be in any order and of type:

- `Symbol`: the output format (`:full`, `:standard` or `:midpoint`)
- `Bool`: whether to display decorations
- `Integer`: the number of significant figures

E.g.
```
julia> x = 0.1..0.3
@[0.0999999, 0.300001]

julia> @format full

julia> x
Interval(0.09999999999999999, 0.30000000000000004)

julia> @format standard 3

julia> x
[0.0999, 0.301]
```
"""
macro format(expr...)

    format = Meta.quot(display_params.format)
    decorations = display_params.decorations
    sigfigs = display_params.sigfigs

    for ex in expr

        if isa(ex, Symbol)
            format = Meta.quot(ex)

        elseif isa(ex, Integer)
            sigfigs = ex

        elseif isa(ex, Bool)
            decorations = ex
        end
    end

    format_code = :(setformat($format, decorations=$decorations, sigfigs=$sigfigs))

    return format_code
end

## Output

# round to given number of signficant digits
# basic structure taken from string(x::BigFloat) in base/mpfr.jl
function round_string(x::BigFloat, digits::Int, r::RoundingMode)

    lng = digits + Int32(8)
    buf = Array{UInt8}(lng + 1)

    lng = ccall((:mpfr_snprintf,:libmpfr), Int32,
    (Ptr{UInt8}, Culong,  Ptr{UInt8}, Int32, Ptr{BigFloat}...),
    buf, lng + 1, "%.$(digits)R*g", Base.MPFR.to_mpfr(r), &x)

    repr = unsafe_string(pointer(buf))

    repr = replace(repr, "nan", "NaN")

    return repr
end

round_string(x::Real, digits::Int, r::RoundingMode) = round_string(big(x), digits, r)

function basic_representation(a::Interval, format=nothing)
    if isempty(a)
        return "∅"
    end

    if format == nothing
        format = display_params.format  # default
    end

    sigfigs = display_params.sigfigs

    local output

    if format == :standard

        aa = round_string(a.lo, sigfigs, RoundDown)
        bb = round_string(a.hi, sigfigs, RoundUp)

        output = "[$aa, $bb]"
        output = replace(output, "inf", "∞")
        output = replace(output, "Inf", "∞")

    elseif format == :full
        output = "Interval($(a.lo), $(a.hi))"

    elseif format == :midpoint
        m = round_string(mid(a), sigfigs, RoundNearest)
        r = round_string(radius(a), sigfigs, RoundUp)
        output = "$m ± $r"
    end

    output
end

function basic_representation(a::Interval{Rational{T}}, format=nothing) where
        T<:Integer

    if isempty(a)
        return "∅"
    end

    if format == nothing
        format = display_params.format  # default
    end

    local output

    if format == :standard
        output = "[$(a.lo), $(a.hi)]"

    elseif format == :full
        output = "Interval($(a.lo), $(a.hi))"

    elseif format == :midpoint
        m = mid(a)
        r = radius(a)
        output = "$m ± $r"
    end

    output
end

function subscriptify(n::Int)
    dig = reverse(digits(n))
    subscript_0 = Int('₀')    # 0x2080
    join( [Char(subscript_0 + i) for i in dig])
end

# fall-back:
representation(a::Interval{T}, format=nothing) where T =
    basic_representation(a, format)

function representation(a::Interval{BigFloat}, format=nothing)

    if format == nothing
        format = display_params.format  # default
    end

    if format == :standard

        return string(basic_representation(a, format), subscriptify(precision(a.lo)))

    else

        return basic_representation(a, format)

    end
end

function representation(a::DecoratedInterval{T}, format=nothing) where T

    if format == nothing
        format = display_params.format  # default
    end

    if format==:full
        return "DecoratedInterval($(representation(interval_part(a), format)), $(decoration(a)))"
    end

    interval = representation(interval_part(a), format)

    if display_params.decorations
        string(interval, "_", decoration(a))
    else
        interval
    end

end

function representation(X::IntervalBox, format=nothing)

    if format == nothing
        format = display_params.format  # default
    end

    if display_params.format == :full
        return string("IntervalBox(", join(X, ", "), ")")

    else
        return join(X, " × ")
    end

end

for T in (Interval, DecoratedInterval)
    @eval show(io::IO, a::$T{S}) where S = print(io, representation(a))
    @eval showall(io::IO, a::$T{S}) where S = print(io, representation(a, :full))
end

T = IntervalBox
@eval show(io::IO, a::$T) = print(io, representation(a))
@eval show(io::IO, ::MIME"text/plain", a::$T) = print(io, representation(a))
@eval showall(io::IO, a::$T) = print(io, representation(a, :full))


text:
#=
 Benchmark two collision checkers:
collision_checker(veh1::Vehicle, veh2::Vehicle) from AutomotivePOMDPs using parallel axis theorem
is_colliding(veh1::Vehicle, veh2::Vehicle) from AutomotiveDrivingModels using Minkowski sums
=#

using AutomotiveDrivingModels
using BenchmarkTools

## Helpers

const ROADWAY = gen_straight_roadway(1, 20.0)

function create_vehicle(x::Float64, y::Float64, θ::Float64 = 0.0; id::Int64=1)
    s = VehicleState(VecSE2(x, y, θ), ROADWAY, 0.0)
    return Entity(s, VehicleDef(), id)
end

const VEH_REF = create_vehicle(0.0, 0.0, 0.0, id=1)

function no_collisions_PA(xspace, yspace, thetaspace)
    for (x,y,th) in zip(xspace, yspace, thetaspace)
        veh2 = create_vehicle(x,y,th,id=2)
        @assert !collision_checker(VEH_REF, veh2) "Error for veh ", veh2
    end
end

function collisions_PA(xspace, yspace, thetaspace)
    for (x,y,th) in zip(xspace, yspace, thetaspace)
        veh2 = create_vehicle(x,y,th,id=2)
        @assert collision_checker(VEH_REF, veh2) "Error for veh ", veh2
    end
end

function no_collisions_MI(xspace, yspace, thetaspace)
    for (x,y,th) in zip(xspace, yspace, thetaspace)
        veh2 = create_vehicle(x,y,th,id=2)
        @assert !is_colliding(VEH_REF, veh2) "Error for veh ", veh2
    end
end

function collisions_MI(xspace, yspace, thetaspace)
    for (x,y,th) in zip(xspace, yspace, thetaspace)
        veh2 = create_vehicle(x,y,th,id=2)
        @assert is_colliding(VEH_REF, veh2) "Error for veh ", veh2
    end
end

function consistency(xspace, yspace, thetaspace)
    for (x,y,th) in zip(xspace, yspace, thetaspace)
        veh2 = create_vehicle(x,y,th,id=2)
        @assert is_colliding(VEH_REF, veh2) == collision_checker(VEH_REF, veh2) "Error for veh ", veh2
    end
end

## Series of test 1: Far field
thetaspace = LinRange(0.0, 2*float(pi), 10)
xspace = LinRange(15.0, 300.0, 10)
yspace = LinRange(10.0, 300.0, 10)

no_collisions_PA(xspace, yspace, thetaspace)
no_collisions_MI(xspace, yspace, thetaspace)
consistency(xspace, yspace, thetaspace)

println(" -- TEST 1 : Far Range Collision Detection -- ")
println(" ")
println("Parallel Axis performance on far range negative collisions: ")
@btime no_collisions_PA($xspace, $yspace, $thetaspace)
println(" ")
println("Minkowski Sum performance on far range negative collisions: ")
@btime no_collisions_MI($xspace, $yspace, $thetaspace)
println(" ")
println(" ------------------------------------------")
println(" ")
## Series of test 2: Superposition
thetaspace = LinRange(0.0, 2*float(pi), 10)
xspace = LinRange(-2.0, 2.0, 10)
yspace = LinRange(1.5, 1.5, 10)

collisions_PA(xspace, yspace, thetaspace)
collisions_MI(xspace, yspace, thetaspace)
consistency(xspace, yspace, thetaspace)

println(" -- TEST 2 : Superposition Collision Detection -- ")
println(" ")
println("Parallel Axis performance on positive collisions: ")
@btime collisions_PA($xspace, $yspace, $thetaspace)
println(" ")
println("Minkowski Sum performance on positive collisions: ")
@btime collisions_MI($xspace, $yspace, $thetaspace)
println(" ")
println(" ------------------------------------------")
println(" ")
## Series of test 3: Close field
thetaspace = [0.0]
yspace = [2.5]
xspace = LinRange(-2.0, 2.0, 10)

no_collisions_PA(xspace, yspace, thetaspace)
no_collisions_MI(xspace, yspace, thetaspace)
consistency(xspace, yspace, thetaspace)

println(" -- TEST 3 : Close Range Collision Detection -- ")
println(" ")
println("Parallel Axis performance on close range negative collisions: ")
@btime no_collisions_PA($xspace, $yspace, $thetaspace)
println(" ")
println("Minkowski Sum performance on close range negative collisions: ")
@btime no_collisions_MI($xspace, $yspace, $thetaspace)
println(" ")
println(" ------------------------------------------")
println(" ")


text:
module GeometricMCMC
  using Base

  require("Distributions")
  using Distributions

  include("types.jl")
  include("stats.jl")
  include("setStep.jl")
  include("mh.jl")
  include("mala.jl")
  include("smmala.jl")
  include("mmala.jl")
  include("hmc.jl")
  include("rmhmc.jl")
  include("zv.jl")

  export
    ## Types
    Model,
    MhOpts,
    MalaOpts,
    SmmalaOpts,
    MmalaOpts,
    HmcOpts,
    RmhmcOpts,

    ## Various stats functions
    logOfNormalCdf, # Function for calculating log of std Normal cdf accurately
    autocov, # Function for calculating the sample autocovariance of a vector
    imse, # Function for calculating Geyer's initial monotone sequence estimator

    ## Auxiliary step functions
    setStep, # Generic step function
    setMalaDriftStep, # Function for adjusting MALA drift step
    setSmmalaDriftStep, # Function for adjusting SMMALA drift step
    setMmalaDriftStep, # Function for adjusting MMALA drift step
    setHmcLeapStep, # Function for adjusting HMC leap step
    setRmhmcLeapStep, # Function for adjusting RMHMC leap step

    ## MCMC functions
    mh, # Metropolis-Hastings (MH)
    mala, # Metropolis adjusted Langevin algorithm (MALA)
    smmala, # Simplified Manifold MALA (SMMALA)
    mmala, # Manifold MALA (MMALA)
    hmc, # Hamiltonian Monte Carlo (HMC)
    rmhmc, # Riemannian manifold HMC (RMHMC)

    ## Functions for calculating zero-variance estimators
    linearZv, # ZV-MCMC using linear polynomial
    quadraticZv # ZV-MCMC using quadratic polynomial
end


text:
module Griddly

    using CxxWrap
    # @wrapmodule(joinpath(@__DIR__,"..","..","..","Debug","bin","JuGriddly"),:define_module_jugriddly)
	@wrapmodule(joinpath(@__DIR__,"..","..","Release","bin","JuGriddly"),:define_module_jugriddly)

    function __init__()
   		@initcxx
    end
	include("docs.jl")
	include("utils.jl")
    ######################################## RenderTool ########################################
    using Makie
	using FixedPointNumbers
	using ImageCore

	include("rendertools.jl")
	export RenderWindow, VideoRecorder, init, render, save_frame, start_video, add_frame!, save_video, MultipleScreen, render_multiple, save_frames

end