text:
"""

DexArray is a container for GPU arrays that supports most of the
AbstractArray interface.  Important differences from the alternative
CudaArray are: (1) a custom memory manager that minimizes the number
of calls to the slow cudaMalloc by reusing already allocated but
garbage collected GPU pointers.  (2) a custom getindex that handles
ranges such as `a[5:10]` as views (with memory shared with the
original array) instead of copies.  DexArrays can be created by
specifying the element type and dimensions or by conversion from
regular Arrays and they can be converted back to regular Arrays (which
involve copying to and from the GPU memory):

    a = DexArray(Float32,2,3)
    b = DexArray(zeros(2,3))
    c = Array(b)

"""
type DexArray{T,N}
    ptr::DexPtr
    dims::NTuple{N,Int}
end

# Aliases:
typealias DexMatrix{T} DexArray{T,2}
typealias DexVector{T} DexArray{T,1}
typealias DexVecOrMat{T} Union{DexVector{T}, DexMatrix{T}}

# Constructors:
DexArray{T,N}(::Type{T}, dims::NTuple{N,Int})=DexArray{T,N}(DexPtr(sizeof(T)*prod(dims)), dims)
DexArray(T::Type, dims::Int...)=DexArray(T,dims)
DexArray(T::Type, d::Integer...)=DexArray(T,convert(Tuple{Vararg{Int}}, d))

# Conversions:
import Base: convert, reshape, vec, unsafe_convert, pointer
# DexArray <- DexArray
convert{T,N}(::Type{DexArray}, x::DexArray{T,N}) = x
convert{T,N}(::Type{DexArray{T}}, x::DexArray{T,N}) = x
convert{T,N}(::Type{DexArray{T,N}}, x::DexArray{T,N}) = x
convert{T,N,S}(::Type{DexArray{T}}, x::DexArray{S,N}) = convert(DexArray{T,N}, x)
convert{T,N,S}(::Type{DexArray{T,N}}, x::DexArray{S,N}) = convert(DexArray{T,N},unsafe_copy!(Array(S, size(x)), 1, x, 1, length(x)))
reshape{T}(a::DexArray{T},dims::Dims)=(if dims==size(a); a; elseif prod(dims)!=length(a); throw(DimensionMismatch()); else; DexArray{T,length(dims)}(a.ptr,dims); end)
reshape(a::DexArray, dims::Int...) = reshape(a, dims)
vec(a::DexArray) = reshape(a, length(a))
# DexArray <- AbstractArray
convert{T,N}(::Type{DexArray}, x::AbstractArray{T,N}) = convert(DexArray{T,N}, x)
convert{T,N,S}(::Type{DexArray{T}}, x::AbstractArray{S,N}) = convert(DexArray{T,N}, x)
convert{T,N,S}(::Type{DexArray{T,N}}, x::AbstractArray{S,N}) = unsafe_copy!(DexArray(T, size(x)), 1, convert(Array{T,N},x), 1, length(x))
# Array <- DexArray
convert{T,N}(::Type{Array}, x::DexArray{T,N}) = convert(Array{T,N}, x)
convert{T,N,S}(::Type{Array{T}}, x::DexArray{S,N}) = convert(Array{T,N}, x)
convert{T,N,S}(::Type{Array{T,N}}, x::DexArray{S,N}) = convert(Array{T,N},unsafe_copy!(Array(S, size(x)), 1, x, 1, length(x)))
# Ptr <- DexArray
unsafe_convert{T}(::Type{Ptr{T}}, a::DexArray) = unsafe_convert(Ptr{T}, pointer(a))
pointer{T}(a::DexArray{T})=convert(Ptr{T}, a.ptr.ptr)
pointer{T}(a::DexArray{T},i)=convert(Ptr{T}, a.ptr.ptr + (i-1)*sizeof(T))

# AbstractArray interface
import Base: eachindex, elsize, eltype, endof, fill!, first, length, linearindexing, ndims, ones, similar, size, stride, zeros
eachindex(a::DexArray) = (1:length(a))
elsize{T}(::DexArray{T}) = sizeof(T)
eltype{T}(::DexArray{T})=T
eltype{T}(::Type{DexArray{T}}) = T
eltype{T,n}(::Type{DexArray{T,n}}) = T
endof(a::DexArray) = length(a)
fill!{T}(a::DexArray{T},x)=(dexfill!(a,T(x),1,length(a));a)
first(a::DexArray) = a[1]
length(a::DexArray)=prod(size(a))
linearindexing(::DexArray)=LinearFast()
ndims{T,N}(a::DexArray{T,N})=N
ones{T}(a::DexArray{T})=fill!(similar(a),one(T))
similar(a::DexArray, T, dims::Dims)      = DexArray(T, dims)
similar(a::DexArray, T, dims::Int...)    = similar(a, T, dims)
similar(a::DexArray, T)                  = similar(a, T, size(a))
similar{T}(a::DexArray{T})               = similar(a, T, size(a))
similar{T}(a::DexArray{T}, dims::Dims)   = similar(a, T, dims)
similar{T}(a::DexArray{T}, dims::Int...) = similar(a, T, dims)
size(a::DexArray)=a.dims
size{T,N}(a::DexArray{T,N},i::Integer)=(if i>N; 1; else; size(a)[i]; end)
stride{T,N}(a::DexArray{T,N},i::Integer)=(if i>N; length(a); else; s=1; for n=1:(i-1); s*=size(a,n); end; s; end)
strides{T,N}(a::DexArray{T,N})=ntuple(n->stride(a,n), N)
zeros{T}(a::DexArray{T})=fill!(similar(a),zero(T))

# Indexing:
import Base: getindex, setindex!

# We will implement indexing ranges as views not copies, if possible (when contiguous).
# For contiguous memory without stride all but the last >1 dimension must be full

# The original getindex(a,i:j...) for AbstractArray copies:
# function _getindex(l::LinearIndexing, A::AbstractArray, I::Union{Real, AbstractArray, Colon}...)
# in abstractarray.jl:487,multidimensional.jl:184.

# which getindex ops does array implement?
# getindex(A::Array, i1::Real)
# getindex(A::Array, i1::Real, i2::Real, I::Real...)
# getindex(A::Array, I::UnitRange{Int})
# getindex(A::Array, c::Colon)
# getindex{T<:Real}(A::Array, I::Range{T})

# Julia #14770
# If I is shorter than ndims(A) but longer than 1 the remaining indices assumed =1
# Also extra 1's at the end of I are ignored

# These two are not sufficient in spite of what the documentation says:
# display goes into an infinite loop!
# getindex{T}(A::DexArray{T}, i::Int)=unsafe_copy!(T[0], 1, A, i, 1)[1]
# setindex!{T}(A::DexArray{T}, v, i::Int)=unsafe_copy!(A, i, T[v], 1, 1)

# First deal with the easy cases: integer indices, a Colon or a UnitRange.

function getindex{T}(A::DexArray{T}, I::Real)
    J = Int(I)
    1 <= J <= length(A) || throw(BoundsError(A,J))
    unsafe_copy!(T[0], 1, A, J, 1)[1]
end

function setindex!{T}(A::DexArray{T}, v, I::Real)
    J = Int(I)
    1 <= J <= length(A) || throw(BoundsError(A,J))
    unsafe_copy!(A, J, T[v], 1, 1)
end

function getindex{T}(A::DexArray{T}, I::Real...)
    J = Base.to_indexes(I...)
    @inbounds for j=1:length(J)
        1 <= J[j] <= size(A,j) || throw(BoundsError(A,J))
    end
    i = sub2ind(size(A), J...)
    unsafe_copy!(T[0], 1, A, i, 1)[1]
end

function setindex!{T}(A::DexArray{T}, v, I::Real...)
    J = Base.to_indexes(I...)
    @inbounds for j=1:length(J)
        1 <= J[j] <= size(A,j) || throw(BoundsError(A,J))
    end
    i = sub2ind(size(A), J...)
    unsafe_copy!(A, i, T[v], 1, 1)
end

function getindex{T}(A::DexArray{T}, I::UnitRange)
    1 <= first(I) <= last(I) <= length(A) || throw(BoundsError(A,I))
    off = 1+(first(I)-1)*sizeof(T)
    len = length(I)*sizeof(T)
    ptr = DexPtr(A.ptr, off, len)
    DexArray{T,1}(ptr, (length(I),))
end

function setindex!{T}(A::DexArray{T}, v, I::UnitRange)
    1 <= first(I) <= last(I) <= length(A) || throw(BoundsError(A,I))
    if isa(v,Number)
        dexfill!(A,T(v),first(I),length(I))
    elseif (isa(v,DexArray) || isa(v,Array))
        length(v)==length(I) || throw(DimensionMismatch())
        eltype(v)==T || (v = convert(Array{T},v))
        unsafe_copy!(A,first(I),v,1,length(I))
    else
        throw(MethodError(setindex!, A, v, I))
    end
end

function getindex(A::DexArray, I::Colon)
    reshape(A,length(A))
end

function setindex!{T}(A::DexArray{T}, v, I::Colon)
    if isa(v,Number)
        dexfill!(A,T(v),1,length(A))
    elseif (isa(v,DexArray) || isa(v,Array))
        length(v)==length(A) || throw(DimensionMismatch())
        eltype(v)==T || (v = convert(Array{T},v))
        unsafe_copy!(A,1,v,1,length(A))
    else
        throw(MethodError(setindex!, A, v, I))
    end
end

# TODO: the following getindex, setindex! work for 1 and 2 dimensions only, write general versions.

function getindex{T,N}(A::DexArray{T,N}, I::Union{Real, UnitRange, Colon}...)
    (nelts,nrows,ncols,firstindex,astep) = indexparams(A,I...)
    B1 = isa(I[1],Colon) ? size(A,1) : length(I[1])
    B2 = isa(I[2],Colon) ? size(A,2) : length(I[2])
    if ncols == 1
        off = 1+(firstindex-1)*sizeof(T)
        len = nrows*sizeof(T)
        ptr = DexPtr(A.ptr, off, len)
        DexArray{T,2}(ptr, (B1,B2))
    else
        B = similar(A, (B1,B2))
        nrows *= sizeof(T); astep *= sizeof(T)
        ccall((:xcopy,libcudex),Void,(Cint,Cint,Cptr,Cint,Cptr,Cint),
              nrows, ncols, pointer(A,firstindex), astep, B, nrows)
        return B
    end
end

function setindex!{T,N}(A::DexArray{T,N}, B, I::Union{Real, UnitRange, Colon}...)
    (nelts,nrows,ncols,firstindex,astep) = indexparams(A,I...)
    aptr0 = pointer(A, firstindex)
    if isa(B,Number)
        B = T(B)
        if ncols == 1
            if T <: Float32;    ccall((:fill_32,libcudex),Void,(Cint,Cfloat, Ptr{Cfloat}), nelts,B,aptr0)
            elseif T<: Float64; ccall((:fill_64,libcudex),Void,(Cint,Cdouble,Ptr{Cdouble}),nelts,B,aptr0)
            else error("$T not supported"); end
        else
            if T <: Float32;    ccall((:xfill_32,libcudex),Void,(Cint,Cint,Cfloat, Ptr{Cfloat}, Cint),nrows,ncols,B,aptr0,astep)
            elseif T<: Float64; ccall((:xfill_64,libcudex),Void,(Cint,Cint,Cdouble,Ptr{Cdouble},Cint),nrows,ncols,B,aptr0,astep)
            else error("$T not supported"); end
        end
    else
        length(B) == nelts || throw(DimensionMismatch())
        B = convert(DexArray{T},B)
        if ncols == 1
            @cuda(cudart,cudaMemcpyAsync,(Cptr,Cptr,Csize_t,UInt32,Cptr),
                  aptr0, B, nelts*sizeof(T), cudadir(A,B), C_NULL)
        else
            nrows *= sizeof(T); astep *= sizeof(T)
            ccall((:xcopy,libcudex),Void,(Cint,Cint,Cptr,Cint,Cptr,Cint), nrows, ncols, B, nrows, aptr0, astep)
        end
    end
    return A
end

function indexparams{T,N}(A::DexArray{T,N}, I::Union{Real, UnitRange, Colon}...)
    N > 2 && error("setindex for ndims > 2 not implemented yet")
    skipped = false
    nrows = nelts = 1
    subs1 = ones(Int,2)
    astep = length(A)
    for i=1:length(I)
        Ii = I[i]
        Ai = size(A,i)
        if isa(Ii, Colon)
            Li = Ai
            subs1[i] = 1
        elseif isa(Ii, Real)
            1 <= Ii <= Ai || throw(DimensionMismatch())
            Li = 1
            subs1[i] = Int(Ii)
        else
            1 <= first(Ii) <= last(Ii) <= Ai || throw(DimensionMismatch())
            Li = length(Ii)
            subs1[i] = first(Ii)
        end
        nelts *= Li
        if !skipped
            nrows *= Li
            if Li < Ai
                skipped = true
                astep = stride(A,i+1)
            end
        end
    end
    ncols = div(nelts, nrows)
    firstindex = sub2ind(size(A),subs1...)
    return (nelts,nrows,ncols,firstindex,astep)
end

# Concatenation:
import Base: hcat, vcat, cat

# Benchmarks in μs for hcat and vcat: a=rand(1000,1000) v=rand(1000), t=v'
#		cpu	gpu	g->c->g	vkernel
# hcat(a,a)	2350	225	16160
# hcat(a,v)	1230	115	6490
# hcat(v,a)	1220	120	6490
# hcat(v,v)	3.53	12.53	48.49
# vcat(a,a)	2630	10980	16590	665
# vcat(a,t)	1350	10860	6550	338
# vcat(t,a)	1360	10850	6570	338
# vcat(v,v)	2.13	12.33	45.40	13.58

# setindex! methods called by hcat/vcat:
# hcat(v,v): I = (Colon(),1:1) I = (Colon(),2:2)
# vcat(v,v): uses single index
# hcat(m,m): I = (Colon(),1:5) I = (Colon(),6:10)
# vcat(m,m): I = (1:3,Colon()) I = (4:6,Colon())

function hcat{T}(a::DexVecOrMat{T}, b::DexVecOrMat{T})
    size(a,1)==size(b,1) || throw(DimensionMismatch())
    c1 = size(a,1)
    c2 = size(a,2) + size(b,2)
    c = DexArray(T, (c1,c2))
    c[:,1:size(a,2)] = a
    c[:,1+size(a,2):end] = b
    return c
end

function vcat{T}(a::DexVector{T}, b::DexVector{T})
    c = DexArray(T, length(a)+length(b))
    c[1:length(a)] = a
    c[1+length(a):end] = b
    return c
end

function vcat{T}(a::DexVecOrMat{T}, b::DexVecOrMat{T})
    size(a,2)==size(b,2) || throw(DimensionMismatch())
    c1 = size(a,1) + size(b,1)
    c2 = size(a,2)
    c = DexArray(T, (c1,c2))
    c[1:size(a,1),:] = a
    c[1+size(a,1):end,:] = b
    return c
end

function cat{T}(d, a::DexVecOrMat{T}, b::DexVecOrMat{T})
    if     d==1; vcat(a,b)
    elseif d==2; hcat(a,b)
    else error("cat($d) not implemented.")
    end
end

# Utilities:

# Generalizing low level copy using linear indexing to/from gpu arrays:
# copy!{T}(dest::Array{T}, doffs::Integer, src::Array{T}, soffs::Integer, n::Integer)
# Note that this is an unsafe operation, no argument or bounds checking performed.
# Defined in Base:
# unsafe_copy!{T}(dest::Ptr{T}, src::Ptr{T}, n) at array.jl:73
# unsafe_copy!{T}(dest::Array{T,N}, doffs, src::Array{T,N}, soffs, n) at array.jl:79

import Base: unsafe_copy!, copy

function unsafe_copy!{T}(dest::Union{DexArray{T},Array{T}}, doffs, src::Union{DexArray{T},Array{T}}, soffs, n; stream=C_NULL)
    @cuda(cudart,cudaMemcpyAsync,(Cptr,Cptr,Csize_t,UInt32,Cptr),
          pointer(dest,doffs), pointer(src,soffs), n*sizeof(T), cudadir(dest,src), stream)
    return dest
end

function cudadir(a,b)
    deva = isa(a,DexArray) && a.ptr.dev >= 0
    devb = isa(b,DexArray) && b.ptr.dev >= 0
    if !deva && !devb; return 0
    elseif deva && !devb; return 1
    elseif !deva && devb; return 2
    elseif deva && devb;  return 3
    end
end

copy(a::DexArray)=unsafe_copy!(similar(a),1,a,1,length(a))

# Efficient fill:
for S in (32,64); T = Symbol("Float$S"); F = "fill_$S"
    @eval function dexfill!(a::DexArray{$T},v::$T,off,len)
        ccall(($F,$libcudex),Void,(Cint,$T,Ptr{$T}),len,v,pointer(a,off))
    end
end

# # AutoGrad functions:
# import AutoGrad: zeroslike, sum_outgrads, OneHot, unary_nd, indexed_function, isequivalent
# zeroslike(a::DexArray)=zeros(a)
# sum_outgrads{T}(a::DexArray{T},b::DexArray{T})=(a+b)
# sum_outgrads(a::DexArray,b::OneHot)=setindex!(a,sum_outgrads(getindex(a,b.index...),b.value),b.index...)
# unary_nd(f, x::DexArray, eps) = reshape(eltype(x)[unary_nd(indexed_function(f, x, i), x[i], eps) for i in 1:length(x)], size(x))
# isequivalent(x::Union{DexArray,AbstractArray}, y::Union{DexArray,AbstractArray}; o...)=(length(x)==length(y) && all(i->isequivalent(x[i],y[i];o...), 1:length(x)))

# Hack for printing without copying the whole DexArray and without inheriting AbstractArray:
import Base: display, summary
type DexDisplay{T,N} <: AbstractArray{T,N}; a::DexArray{T,N}; end
getindex(a::DexDisplay, i...) = getindex(a.a, i...)
size(a::DexDisplay) = size(a.a)
summary(a::DexDisplay) = summary(a.a)
summary(a::DexArray) = string(Base.dims2string(size(a)), " ", typeof(a))
display(a::DexArray) = display(DexDisplay(a))
# AutoGrad._dbg(a::DexArray) = "K$(join([AutoGrad.id2(a),size(a)...],'_'))"

# curand functions:

import Base: rand!
rand!(a::DexArray{Float32})=(@cuda(curand,curandGenerateUniform,(Cptr,Ptr{Cfloat},Csize_t),rng(),a,length(a)); a)
rand!(a::DexArray{Float64})=(@cuda(curand,curandGenerateUniformDouble,(Cptr,Ptr{Cdouble},Csize_t),rng(),a,length(a)); a)

let RNG=0
global rng
function rng(init=false)
    if RNG==0 || init
        ptr = Cptr[0]
        # CURAND_RNG_PSEUDO_DEFAULT = 100, ///< Default pseudorandom generator
        @cuda(curand,curandCreateGenerator,(Cptr,Cint),ptr,100)
        RNG = ptr[1]
    end
    return RNG
end
end


text:
using DiscontinuousGherkins
using DiscontinuousGherkins: jacobi_gauss_lobatto, vandermonde_1d, Dmatrix1D

NODETOL = 1e-10

N = 8  # Order of polymomials used for approximation
K = 10 # Number of elements

domain = (0, 2)
grid = DiscontinuousGalerkin1DGrid(domain, K)

Np = N+1
Nfp = 1
Nfaces = 2

r = jacobi_gauss_lobatto(0, 0, N)

V = vandermonde_1d(N, r)
V⁻¹ = inv(V)

Dr = Dmatrix1D(N, r, V)

Emat = zeros(Np, Nfaces * Nfp)
Emat[1,  1] = 1
Emat[Np, 2] = 1
lift = V * (V' * Emat)

# build coordinates of all the nodes
VX = grid.node_coordinates
EToV = grid.element_to_node_connectivity
va = EToV[:, 1]'
vb = EToV[:, 2]'
x = ones(N+1, 1) * VX[Int.(va)] .+ 0.5 .* (r .+ 1) * (VX[Int.(vb)] - VX[Int.(va)])

# Compute the metric elements for the local mappings of the 1D elements.
xr = Dr * x
J = xr
rx = 1 ./ J

# Compute masks for edge nodes
Fmask = BitArray(undef, size(x)...)
Fmask .= 0
Fmask[1, 1] = 1
Fmask[end, end] = 1
Fx = x[Fmask]

# % Build surface normals and inverse metric at surface
nx = zeros(Nfp * Nfaces, K)
nx[1, :] .= -1
nx[2, :] .= +1

Fscale = 1 ./ J[Fmask]

## Build global connectivity arrays for 1D grid based on standard EToV input array from grid generator

Nfaces = 2

# Find number of elements and vertices
TotalFaces = Nfaces * K
Nv = K+1

# List of local face to local vertex connections
vn = Int[1, 2]

# Build global face to node sparse array
SpFToV = zeros(TotalFaces, Nv)
sk = 1;
for k in 1:K, face in 1:Nfaces
    SpFToV[sk, EToV[k, vn[face]] |> Int] = 1
end

# Build global face to global face sparse array
SpFToF = SpFToV * SpFToV' - I(TotalFaces)

% Find complete face to face connections
[faces1, faces2] = find(SpFToF==1);

% Convert face global number to element and face numbers
element1 = floor( (faces1-1)/Nfaces ) + 1;
face1 = mod( (faces1-1), Nfaces ) + 1;
element2 = floor( (faces2-1)/Nfaces ) + 1;
face2 = mod( (faces2-1), Nfaces ) + 1;

% Rearrange into Nelements x Nfaces sized arrays
ind = sub2ind([K, Nfaces], element1, face1);
EToE = (1:K)’*ones(1,Nfaces);
EToF = ones(K,1)*(1:Nfaces);
EToE(ind) = element2; EToF(ind) = face2;
return

ext: jl


text:
_isfulldummy(x::AbstractTerm) = false
function _isfulldummy(x::CategoricalTerm)
    return isa(x.contrasts, StatsModels.ContrastsMatrix{StatsModels.FullDummyCoding})
end

function _fulldummycheck(outer::InteractionTerm)
    all(_isfulldummy, outer.terms[1:end-1]) ||
        throw(ArgumentError("Outer interactions in a nesting must consist only " *
                            " of categorical terms with FullDummyCoding, got $outer"))
    return nothing
end

"""
    group / term

Generate predictors for `term` within each level of `group`.  Implemented as
`group + fulldummy(group) & term`.
"""
function Base.:(/)(outer::CategoricalTerm, inner::AbstractTerm)
    return outer + fulldummy(outer) & inner
end

function Base.:(/)(outer::CategoricalTerm, inner::TermTuple)
    fd = fulldummy(outer)
    return mapfoldl(x -> fd & x, +, inner; init=outer)
end

function Base.:(/)(outer::TermTuple, inner::Union{AbstractTerm, TermTuple})
    return outer[1:end-1] + last(outer) / inner
end

function Base.:(/)(outer::InteractionTerm, inner::AbstractTerm)
    # we should only get here via expansion where the interaction term,
    # but who knows what devious things users will try
    _fulldummycheck(outer)
    return outer + outer & inner
end

function Base.:(/)(outer::InteractionTerm, inner::TermTuple)
    # we should only get here via expansion where the interaction term,
    # but who knows what devious things users will try
    _fulldummycheck(outer)
    return mapfoldl(x -> outer & x, +, inner; init=outer)
end

function Base.:(/)(outer::AbstractTerm, inner::AbstractTerm)
    throw(ArgumentError("Nesting terms requires categorical grouping term, got $outer / $inner " *
                        "Manually specify $outer as `CategoricalTerm` in hints/contrasts"))
end

function StatsModels.apply_schema(
    t::FunctionTerm{typeof(/)},
    sch::StatsModels.FullRank,
    Mod::Type{<:RegressionModel},
)
    length(t.args_parsed) == 2 ||
        throw(ArgumentError("malformed nesting term: $t (Exactly two arguments required)"))

    any(x -> isa(x, ConstantTerm), t.args_parsed) && return t

    args = apply_schema.(t.args_parsed, Ref(sch), Mod)

    return first(args) / last(args)
end


text:
########## General statistics ##########

# normalization factor of p-th momeent of a standard normal distribution
normal_moment_factor = p -> 2^(p/2) * gamma((p+1)/2)/sqrt(pi)

"""
Compute the symmetric polynomial weight vector.

# Args
- N: length of time window
- p: order of polynomial, p=0 gives uniform weight
"""
function poly_weight(N::Integer, p::Integer)
    # println("poly_weight $(N)")
    w = ones(N)
    for n=0:p-1
        w = cumsum(w)
    end
    if N>=2
        w[(N÷2):end] = reverse(w[1:N-(N÷2)+1])
    end
    return w/sum(w)
end

"""
Causal polynomial weight for a time window, i.e. the weight is polynomial increasing from the left end to the right end of a time window.

# Note
`p=0` gives uniform weight, `p=1` gives linear weight etc.
"""
causal_weight(N, p) = p == 0 ? poly_weight(N,p) : cumsum(poly_weight(N,p-1))

"""
    cov(X::AbstractVecOrMat, Y::AbstractVecOrMat, w::AbstractWeights)

Reweighted covariance between `X` and `Y`, where each row is an observation.
"""
function cov(X::AbstractVecOrMat, Y::AbstractVecOrMat, w::AbstractWeights)
    # w is always a column vector, by definition of AbstractWeights
    @assert size(X, 1) == size(Y, 1) == length(w)
    # weighted mean
    mX = mean(X, w, 1)
    mY = mean(Y, w, 1)
    return (X .- mX)' * (w .* (Y .- mY)) / sum(w)
end

"""
Robust estimator of variance
"""
function robustvar(X::AbstractArray{T}; dims::Int=1) where {T<:Number}
    lqt = -0.393798799600891  # log of 3/4 quantil of standard normal distribution
    Q = log.((X .- mean(X, dims=dims)).^2)
    return exp.(median(Q, dims=dims) .- 2lqt)
end

"""Principal Component Analysis.

# Args
-X (2d array): each row represents a variable and each column represents an observation
-nc (int): number of components to hold

# Returns
- C, U : coefficients and corresponded principal directions
"""
function pca(X::AbstractMatrix, nc::Int, w::AbstractWeights, center=false)
    U, S, V = svd(cov(X, X, w))
    mX = mean(X, w, 1)
    C = center ? (X .- mX) * U' : X * U'

    return C[:,1:nc], U[:,1:nc], S[1:nc]
end

pca(X::AbstractMatrix, nc::Int, center=false) = pca(X, nc, StatsBase.weights(ones(size(X,1))), center)

"""
    moment_incr(X::AbstractVector{<:Number}, d::Integer, p::Real, w=StatsBase.AbstractWeights[])

Comput the `p`-th moment for the increment vector of `X` with lag `d`.
"""
function moment_incr(X::AbstractVector{<:Number}, d::Integer, p::Real, w=StatsBase.AbstractWeights[])
    return if length(w) > 0
        mean((abs.(X[d+1:end] - X[1:end-d])).^p, w)
    else
        mean((abs.(X[d+1:end] - X[1:end-d])).^p)
    end
end

######## Linear models ########
"""
Multi-linear regression in the column direction.
"""
function multi_linear_regression_colwise(Y::AbstractVecOrMat{<:Real}, X::AbstractVecOrMat{<:Real}, w::AbstractWeights)
    @assert size(Y,1)==size(X,1)==length(w)
    # println("size(Y)=",size(Y))
    # println("size(X)=",size(X))
    μy = mean(Y, w, 1)  # Julia function (version <= 0.7) `mean` does not take keyword argument `dims=1` if weight is passed. This yields a scalar of Y is vector.
    μx = mean(X, w, 1)
    Σyx = cov(Y, X, w)   # This calls user defined cov function. This yields scalar if Y and X are vectors.
    Σxx = cov(X, X, w)

    # Σxx += 1e-5 * Matrix{Float64}(I,size(Σxx))  # perturbation
    # A = Σyx / Σxx  # scalar or matrix, i.e., Σyx * inv(Σxx)

    A = Σyx * pinv(Σxx)  # scalar or matrix, i.e., Σyx * inv(Σxx)

    # A = try
    #     Σyx * pinv(Σxx)  # scalar or matrix, i.e., Σyx * inv(Σxx)
    # catch
    #     zero(Σyx)
    # end

    β = μy - A * μx'  # scalar or vector, transpose has no effect on scalar
    E = Y - (A * X' .+ β)'
    Σ = cov(E, E, w)
    return A, β, E, Σ
end

"""
    linear_regression(Y::AbstractMatrix, X::AbstractMatrix, w::AbstractWeights; dims::Integer=1)

Reweighted (multi-)linear regression of data matrix `Y` versus `X` in the given dimension (dimension of observation).

# Returns
- (A, β): explanatory slope (scalar or matrix) and intercept (scalar or vector)
- E: residual
- Σ: covariance of the residual
"""
function linear_regression(Y::AbstractMatrix, X::AbstractMatrix, w::AbstractWeights; dims::Integer=1)
    if dims==1
        return multi_linear_regression_colwise(Y, X, w)
    else
        return multi_linear_regression_colwise(Y', X', w)
    end
end

linear_regression(Y::AbstractMatrix, X::AbstractMatrix; dims::Integer=1) =
    linear_regression(Y, X, StatsBase.weights(ones(size(Y,dims))); dims=dims)

function linear_regression(Y::AbstractVector, X::AbstractVecOrMat, w::AbstractWeights)
    return multi_linear_regression_colwise(Y, X, w)
end

linear_regression(Y::AbstractVector, X::AbstractVecOrMat) = linear_regression(Y, X, StatsBase.weights(ones(size(Y,1))))

# function linear_prediction(R::Tuple, X0::AbstractVecOrMat, Y0::AbstractVecOrMat=[])
#     A,β = R[1]
#     Yp = A * X0 .+ β
#     Ep = isempty(Y0) ? [] : Y0-Yp  # error of prediction
#     return Yp, Y0, Ep
# end

"""
    conv_regression(X::AbstractVector{<:Real}, Y::AbstractVector{<:Real}, w::Integer; mode::Symbol=:causal, nan::Symbol=:ignore)

Estimate by linear regression the convolution kernel and bias such that `kernel * X + bias = Y`.

# Args
- X: input vector
- Y: output vector
- w: length of the convolution kernel
- mode: mode of convolution, :causal or :anticausal
- nan: how NaN values are handled: :ignore for ignoring or :zero for replacing by zero

# Returns
A, b, E, σ^2: kernel, bias, error and variance of error
"""
function conv_regression(X::AbstractVector{<:Real}, Y::AbstractVector{<:Real}, w::Integer, p::Integer=1; mode::Symbol=:causal, nan::Symbol=:ignore)
    @assert length(X) == length(Y)

    tidx, Xm0 = rolling_vectorize(X, w, 1, p; mode=mode)
    Xm = transpose(Xm0)
    Yr = Y[tidx]

    if nan==:ignore
        idx = findall((any(isnan.(Xm), dims=2)[:] .+ isnan.(Yr)).==0)
        # println((any(isnan.(Xm), dims=2)[:] .+ isnan.(Yr)).==0)
        # println(any(isnan.(Xm), dims=2)[:])
        # println(isnan.(Yr))
        # print(idx)
        Yr = Yr[idx]
        Xm = Xm[idx,:]
    elseif nan==:zero
        Yr[isnan.(Yr)] .= 0
        Xm[isnan.(Xm)] .= 0
    else
        error("Unknown method for NaN values")
    end

    res = linear_regression(Yr, Xm)
    # res = IRLS(Yr, Xm, .5)

    A, b, E, σ2 = reverse(res[1][:]), res[2][end], res[3], sqrt(res[4][end])
    return A, b, E, σ2
end

######## Nonlinear models ########
"""
IRLS
"""
function IRLS(Y::AbstractVecOrMat{<:Real}, X::AbstractVecOrMat{<:Real}, pnorm::Real=2.; maxiter::Integer=10^3, tol::Real=10^-3, vreg::Real=1e-8)
    @assert pnorm > 0

    wfunc = E -> (sqrt.(sum(E.*E, dims=2) .+ vreg).^(pnorm-2))[:]  # function for computing weight vector

    A, β, E, Σ = linear_regression(Y, X)  # initialization
    w0 =  wfunc(E) # weight vector
    # w0 /= sum(w0)
    err = 0.

    for n=1:maxiter
        # println("Iter=$n")
        A, β, E, Σ = linear_regression(Y, X, StatsBase.weights(w0))
        w = wfunc(E)
        err = LinearAlgebra.norm(w - w0)/LinearAlgebra.norm(w0)
        w0 = w / sum(w)
        if  err < tol
            break
        end
    end
    # println(n)
    # println(err)
    Σ = sum(sqrt.(sum(E.*E, dims=2)).^pnorm)
    return A, β, E, Σ, w0
end

###### Online Statistics ######

# function exponential_moving_average!(X::AbstractVector{<:Number}, α::Real)
#     # @assert 0<α<=1
#     for t=2:size(X,1)
#         X[t] = α * X[t] + (1-α) * (isnan(X[t-1]) ? 0. : X[t-1])  # NaN safe
#     end
#     return X
# end

function exponential_moving_average!(X::AbstractVecOrMat{<:Number}, α::Real)
    @assert 0<α<=1
    X[findall(isnan.(X))] .= 0.  # NaN safe

    for t=2:size(X,1)  # first axis is the time
        W = view(X,t,:)
        W .= α * W + (1-α) * view(X,t-1,:)
    end
    return X
end

"""
EMA of a vector or a matrix (in row direction).
"""
function exponential_moving_average(X::AbstractVecOrMat{<:Number}, α::Real, n::Integer=1)
    @assert 0<α<=1
    X1 = copy(X)
    Xm = copy(X1)
    exponential_moving_average!(Xm, α)

    for t=2:n
        X1 += X-Xm  # N-fold EMA
        # or the zero-lag EMA filter: X1 += X-circshift(X, 1)
        Xm = copy(X1)
        exponential_moving_average!(Xm, α)
    end
    return Xm
end

# function simple_moving_average!(X::AbstractVector{<:Number}, n::Integer)
#     @assert n>0
#     for t=2:size(X,1)
#         X[t] = (X[t] + X[t-1] * min(n, t-1) - (t>n ? X[t-n] : 0)) / min(n, t)
#     end
#     return X
# end

function simple_moving_average(X0::AbstractVecOrMat{<:Number}, n::Integer)
    @assert n>0

    X1 = copy(X0)  # X0 must be saved in the causal implementation
    X1[findall(isnan.(X1))] .= 0.  # NaN safe
    X = copy(X1)

    for t=2:n
        X[t,:] .= (X[t,:] + X[t-1,:] * (t-1)) / t
    end
    for t=n+1:size(X,1)
        X[t,:] .= (X[t,:] + X[t-1,:] * n - X1[t-n,:]) / n
    end
    for t=2:size(X,1)
        W = view(X,t,:)
        W .= (W + view(X,t-1,:) * min(n, t-1) .- (t>n ? view(X1,t-n,:) : 0)) / min(n, t)
        # # equivalent to
        # X[t,:] .= (X[t,:] + X[t-1,:] * min(n, t-1) .- (t>n ? X[t-n,:] : 0)) / min(n, t)
    end
    return X
end

max_stars_repo_name: yanncalec/FracFin

max_issues_repo_name: yanncalec/FracFin

max_forks_repo_name: yanncalec/FracFin


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#    Tests of tables with UTF-8 characters.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

@testset "Strings with characters that have variable width" begin
    matrix = ["😄"^10 "😄"^10 "😅"^10; "🧐"^5 "🥺"^5 "😇"^5; "a"^10 "a"^10 "a"^10]
    header = ["😋"^5 "😁"^10 "🤣"^15; "⚡"^15 "👽"^10 "🤩"^5]

    expected = """
┌────────────────────────────────┬──────────────────────┬────────────────────────────────┐
│ 😋😋😋😋😋                     │ 😁😁😁😁😁😁😁😁😁😁 │ 🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣 │
│ ⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡ │ 👽👽👽👽👽👽👽👽👽👽 │                     🤩🤩🤩🤩🤩 │
├────────────────────────────────┼──────────────────────┼────────────────────────────────┤
│ 😄😄😄😄😄😄😄😄😄😄           │ 😄😄😄😄😄😄😄😄😄😄 │           😅😅😅😅😅😅😅😅😅😅 │
│ 🧐🧐🧐🧐🧐                     │      🥺🥺🥺🥺🥺      │                     😇😇😇😇😇 │
│ aaaaaaaaaa                     │      aaaaaaaaaa      │                     aaaaaaaaaa │
└────────────────────────────────┴──────────────────────┴────────────────────────────────┘
"""

    result = pretty_table(String, matrix, header, alignment = [:l, :c, :r])

    @test result == expected

    # Crop
    # ----

    expected = """
┌────────────────────────────────┬──────────────────────┬─────────────────────────────────
│ 😋😋😋😋😋                     │ 😁😁😁😁😁😁😁😁😁😁 │ 🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣 ⋯
│ ⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡ │ 👽👽👽👽👽👽👽👽👽👽 │                     🤩🤩🤩🤩🤩 ⋯
├────────────────────────────────┼──────────────────────┼─────────────────────────────────
│ 😄😄😄😄😄😄😄😄😄😄           │ 😄😄😄😄😄😄😄😄😄😄 │           😅😅😅😅😅😅😅😅😅😅 ⋯
│ 🧐🧐🧐🧐🧐                     │      🥺🥺🥺🥺🥺      │                     😇😇😇😇😇 ⋯
│ aaaaaaaaaa                     │      aaaaaaaaaa      │                     aaaaaaaaaa ⋯
└────────────────────────────────┴──────────────────────┴─────────────────────────────────
"""

    result = pretty_table(String, matrix, header,
                          alignment = [:l, :c, :r],
                          crop = :both,
                          display_size = (0,90))

    @test result == expected

    expected = """
┌────────────────────────────────┬──────────────────────┬────────────────────────────────
│ 😋😋😋😋😋                     │ 😁😁😁😁😁😁😁😁😁😁 │ 🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣🤣  ⋯
│ ⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡⚡ │ 👽👽👽👽👽👽👽👽👽👽 │                     🤩🤩🤩🤩  ⋯
├────────────────────────────────┼──────────────────────┼────────────────────────────────
│ 😄😄😄😄😄😄😄😄😄😄           │ 😄😄😄😄😄😄😄😄😄😄 │           😅😅😅😅😅😅😅😅😅  ⋯
│ 🧐🧐🧐🧐🧐                     │      🥺🥺🥺🥺🥺      │                     😇😇😇😇  ⋯
│ aaaaaaaaaa                     │      aaaaaaaaaa      │                     aaaaaaaaa ⋯
└────────────────────────────────┴──────────────────────┴────────────────────────────────
                                                                         1 column omitted
"""

    result = pretty_table(String, matrix, header,
                          alignment = [:l, :c, :r],
                          crop = :both,
                          display_size = (0,89))

    @test result == expected
end


text:
using ForneyLab
using LinearAlgebra
import ForneyLab: SoftFactor, @ensureVariables, generateId, addNode!, associate!,
                  averageEnergy, Interface, Variable, slug, ProbabilityDistribution,
                  differentialEntropy, unsafeLogMean, unsafeMean, unsafeCov, unsafePrecision, unsafeMeanCov,
                  collectAverageEnergyInbounds, localEdgeToRegion, ultimatePartner, region, Region, isClamped, currentInferenceAlgorithm
import SpecialFunctions: digamma
export NAutoregressiveX, NARX

"""
Description:

    A Nonlinear Autoregressive model with eXogeneous input (NARX)

    y_t = f(y_t-1, …, y_t-M, u_t, u_t-1, …, u_t-N) + e_t

    where M is the auto-regression order of observations and N of inputs.
    Concatenate previous observations: x_t-1 = [y_t-1, …, y_t-M]' and
    previous inputs z_t-1 = [u_t-1, …, u_t-N]'.

    Assume y_t, x_t-1, z_t-1 and u_t are observed and e_t ~ N(0, τ^-1).

    f is assumed to be a linear product of coefficients θ and a basis expansion
    of inputs and outputs ϕ:

        f(...) = θ'*ϕ(y_t-1, …, y_t-n_y, u_t, …, u_t-n_u)

Interfaces:

    1. y (output)
    2. θ (function coefficients)
    3. x (previous observations vector)
    4. u (input)
    5. z (previous inputs vector)
    6. τ (precision)

Construction:

    NAutoregressiveX(y, θ, x, u, z, τ, g=ϕ, id=:some_id)
"""

mutable struct NAutoregressiveX <: SoftFactor
    id::Symbol
    interfaces::Vector{Interface}
    i::Dict{Symbol,Interface}

    g::Function # Scalar function between autoregression coefficients and state variable

    function NAutoregressiveX(y, θ, x, u, z, τ; g::Function, id=generateId(NAutoregressiveX))
        @ensureVariables(y, θ, x, u, z, τ)
        self = new(id, Array{Interface}(undef, 6), Dict{Symbol,Interface}(), g)
        addNode!(currentGraph(), self)
        self.i[:y] = self.interfaces[1] = associate!(Interface(self), y)
        self.i[:θ] = self.interfaces[2] = associate!(Interface(self), θ)
        self.i[:x] = self.interfaces[3] = associate!(Interface(self), x)
        self.i[:u] = self.interfaces[4] = associate!(Interface(self), u)
        self.i[:z] = self.interfaces[5] = associate!(Interface(self), z)
        self.i[:τ] = self.interfaces[6] = associate!(Interface(self), τ)
        return self
    end
end

slug(::Type{NAutoregressiveX}) = "NARX"

function averageEnergy(::Type{NAutoregressiveX},
                       g::Function,
                       marg_y::ProbabilityDistribution{Univariate},
                       marg_θ::ProbabilityDistribution{Multivariate},
                       marg_x::ProbabilityDistribution{Multivariate},
                       marg_u::ProbabilityDistribution{Univariate},
                       marg_z::ProbabilityDistribution{Multivariate},
                       marg_τ::ProbabilityDistribution{Univariate})

    # Extract moments of beliefs
    mθ,Vθ = unsafeMeanCov(marg_θ)
    my = unsafeMean(marg_y)
    mx = unsafeMean(marg_x)
    mu = unsafeMean(marg_u)
    mz = unsafeMean(marg_z)
    mτ = unsafeMean(marg_τ)
    aτ = marg_τ.params[:a]
    bτ = marg_τ.params[:b]

    # Gradient of f evaluated at mθ
	# Jθ = Zygote.gradient(g, mθ, mx, mu, mz)[1]
	Jθ = g([mx; mu; mz])

	# Evaluate f at mθ
	fθ = mθ'*Jθ

    return 1/2*log(2*pi) - 1/2*(digamma(aτ)-log(bτ)) + 1/2*mτ*(my^2 - 2*my*fθ + fθ^2 + Jθ'*Vθ*Jθ)
end

function collectAverageEnergyInbounds(node::NAutoregressiveX)
    inbounds = Any[]

    # Push function to calling signature (g needs to be defined in user scope)
	push!(inbounds, Dict{Symbol, Any}(:g => node.g, :keyword => false))

    local_edge_to_region = localEdgeToRegion(node)

    encountered_regions = Region[] # Keep track of encountered regions
    for node_interface in node.interfaces
        inbound_interface = ultimatePartner(node_interface)
        current_region = region(node_interface.node, node_interface.edge)

        if isClamped(inbound_interface)
            # Hard-code marginal of constant node in schedule
            push!(inbounds, assembleClamp!(copy(inbound_interface.node), ProbabilityDistribution)) # Copy Clamp before assembly to prevent overwriting dist_or_msg field
        elseif !(current_region in encountered_regions)
            # Collect marginal entry from marginal dictionary (if marginal entry is not already accepted)
            target = local_edge_to_region[node_interface.edge]
            current_inference_algorithm = currentInferenceAlgorithm()
            push!(inbounds, current_inference_algorithm.target_to_marginal_entry[target])
        end

        push!(encountered_regions, current_region)
    end

    return inbounds
end

max_stars_repo_name: wmkouw/NARX

max_issues_repo_name: wmkouw/NARX

max_forks_repo_name: wmkouw/NARX


text:
module HeatEquation

using ProgressMeter
using Plots

include("setup.jl")
include("io.jl")
include("core.jl")

export Field, simulate!, initialize, visualize, average_temperature

end # module


text:
### A Pluto.jl notebook ###
# v0.12.21

using Markdown
using InteractiveUtils

# ╔═╡ ebbb8758-126a-11eb-3682-31a80e9f8437
begin
	using Empirikos
	using JuMP
	using RCall
	using LaTeXStrings
	using Plots
	using StatsPlots
	using MosekTools
	using Random
end

# ╔═╡ ff3c85b6-126a-11eb-3fc0-635d27bf1d7e
md"""
# Reproducing the REBayes R vignette

The goal of this notebook is showcase the functionality of the `Empirikos.jl` package by (partially) reproducing results shown in the [REBayes vignette](https://cran.r-project.org/web/packages/REBayes/vignettes/rebayes.pdf).
"""

# ╔═╡ 88f13064-12ba-11eb-2329-b1b7196a3af4
pgfplotsx()

# ╔═╡ 675ba12a-12b8-11eb-0791-77781a1619fa
md"""
## Gaussian Mixture models

### Needles and haystacks
"""

# ╔═╡ 8a5b07c4-12b8-11eb-0884-19f987050ba4
normal_Zs = StandardNormalSample.(rand(MersenneTwister(100), Normal(), 1000) .+ 										  [fill(0,900); fill(2,100)]);

# ╔═╡ cf8c5152-12b8-11eb-220e-b9750bf9cc9f
normal_npmle = fit(NPMLE(DiscretePriorClass(), Mosek.Optimizer), normal_Zs);

# ╔═╡ e7ab7a1e-12b9-11eb-2dc6-0180b9df937d
normal_postmean_targets =  PosteriorMean.(sort(normal_Zs));

# ╔═╡ 1731df1a-12b9-11eb-2aeb-2762823871f4
begin
normal_pl_marginal = plot(marginalize(NormalSample(1.0),
		                  DiscreteNonParametric([0.0;2.0], [0.9;0.1])), 									  components=false, label=nothing, xguide="z",
                          yguide = "f(z)",
						  title = "True mixture")
normal_pl_npmle_prior = plot(support(normal_npmle.prior),
	 	  					 probs(normal_npmle.prior), seriestype=:sticks,
	                         xguide = L"\mu", yguide=L"g(\mu)",
		                     title = "NPMLE prior estimate", label=nothing)
normal_pl_posterior_mean = plot( response.(location.(normal_postmean_targets)),
		                         normal_postmean_targets.(normal_npmle.prior),
	                             label=nothing, title="NPMLE posterior mean",
	                             xguide="z", yguide=L"E[\mu \mid Z=z]")
plot(normal_pl_marginal, normal_pl_npmle_prior, normal_pl_posterior_mean,
	 size=(750,250), layout=(1,3))
end

# ╔═╡ fa7eae00-126a-11eb-1fe6-09bb4826877e
md"""
## Binomial NPMLE

We start by loading the Tacks data by Beckett and Diaconis.
"""

# ╔═╡ 470e7048-126b-11eb-3eb1-7bca6628f555
R"library(REBayes)";

# ╔═╡ ed1ab854-126b-11eb-399f-677c45bf7ecd
tacks_tbl = rcopy(R"tacks");

# ╔═╡ 04590064-126c-11eb-2023-9767e8bead48
first(tacks_tbl,3)

# ╔═╡ 5b6c41cc-126c-11eb-315d-8752bfde9d9a
md"""
Let us convert the table as a vector of Binomial samples
"""

# ╔═╡ 89cd6128-126b-11eb-2164-9d634900847a
tack_Zs = BinomialSample.(tacks_tbl.x, Int64.(tacks_tbl.k))

# ╔═╡ c774a158-126b-11eb-2241-a9b6cf80c681
tacks_npmle = fit(NPMLE(DiscretePriorClass(), Mosek.Optimizer), tack_Zs);

# ╔═╡ 345669c6-126e-11eb-2b12-a976e63e3335
plot(support(tacks_npmle.prior),
	 probs(tacks_npmle.prior),
	 seriestype=:sticks,
	 xguide = "p",
	 yguide = "g(p)",
	 size=(350,250),
	 label=nothing)

# ╔═╡ 5a8ad83e-126e-11eb-2d85-6581c5379a4b
md"""
# Poisson NPMLE
"""

# ╔═╡ b43cdce2-126e-11eb-0da0-fd2ec7d7bbad
norberg = rcopy(R"Norberg"); first(norberg,3)

# ╔═╡ cccfbc70-126e-11eb-0439-4552285a634d
norberg_Zs = Empirikos.PoissonSample.(norberg.Death, norberg.Exposure ./ 344);

# ╔═╡ f6e28740-1273-11eb-1464-ab04ba549059
extrema(response.(norberg_Zs) ./ nuisance_parameter.(norberg_Zs))

# ╔═╡ 3460571e-126f-11eb-176b-a3445db55987
sum(response.(norberg_Zs) .== 0) # how many zeros in the samples?

# ╔═╡ 08c5dc30-1276-11eb-01c1-4baf4dece1d5
mosek_lowreltol = optimizer_with_attributes(Mosek.Optimizer,
					"MSK_DPAR_INTPNT_CO_TOL_REL_GAP" => 10^(-8))

# ╔═╡ e7e4d35a-126f-11eb-34aa-dbb75f86c6a6
norberg_npmle = fit(NPMLE(DiscretePriorClass(),
		                  mosek_lowreltol;
		                  prior_grid_length=1000),
					norberg_Zs);

# ╔═╡ 80fcb26e-12a3-11eb-0977-3f37bdd5759e
norberg_gamma_mle = fit(ParametricMLE(Gamma()), norberg_Zs)

# ╔═╡ 56671378-1271-11eb-1a8f-3f289512a146
begin
pl = plot(plot(support(norberg_npmle.prior),
	 	  probs(norberg_npmle.prior)./sum(probs(norberg_npmle.prior)),
	      seriestype=:sticks,
	      xguide = L"\theta",
	      yguide = L"g(\theta)",
		  legend=:topright,
		  label="NPMLE"),
	 plot(support(norberg_npmle.prior),
	 	  (probs(norberg_npmle.prior) ./ sum(probs(norberg_npmle.prior))).^(1/3),
	      seriestype=:sticks,
	      xguide = L"\theta",
	      yguide = L"\sqrt[3]{g(\theta)}",
		  legend=:topright,
		  label="NPMLE"),
	size=(700,250))
plot!(pl[1], support(norberg_npmle.prior),
		     x->pdf(norberg_gamma_mle, x), color=:purple, label="Γ MLE")
end

# ╔═╡ db81cb56-12a4-11eb-3b5b-910d159781fd
norgberg_postmean_gamma = (PosteriorMean.(norberg_Zs)).(norberg_gamma_mle)

# ╔═╡ 316b4056-12a5-11eb-129f-bf07f3962a39
norgberg_postmean_npmle = (PosteriorMean.(norberg_Zs)).(norberg_npmle.prior)

# ╔═╡ a0f0629e-12a5-11eb-3d7e-c92f58bb13ad
begin
	plot(norgberg_postmean_gamma, norgberg_postmean_npmle, seriestype=:scatter,
		 xguide= "Γ MLE posterior mean",
		 yguide = "NPMLE posterior mean", label=nothing,
		 size=(300,250), aspect_ratio=1)
	plot!(norgberg_postmean_gamma, norgberg_postmean_gamma,
		 linestyle=:dash, label=nothing)
end

# ╔═╡ Cell order:
# ╟─ff3c85b6-126a-11eb-3fc0-635d27bf1d7e
# ╠═ebbb8758-126a-11eb-3682-31a80e9f8437
# ╠═88f13064-12ba-11eb-2329-b1b7196a3af4
# ╟─675ba12a-12b8-11eb-0791-77781a1619fa
# ╠═8a5b07c4-12b8-11eb-0884-19f987050ba4
# ╠═cf8c5152-12b8-11eb-220e-b9750bf9cc9f
# ╠═e7ab7a1e-12b9-11eb-2dc6-0180b9df937d
# ╠═1731df1a-12b9-11eb-2aeb-2762823871f4
# ╟─fa7eae00-126a-11eb-1fe6-09bb4826877e
# ╠═470e7048-126b-11eb-3eb1-7bca6628f555
# ╠═ed1ab854-126b-11eb-399f-677c45bf7ecd
# ╠═04590064-126c-11eb-2023-9767e8bead48
# ╟─5b6c41cc-126c-11eb-315d-8752bfde9d9a
# ╠═89cd6128-126b-11eb-2164-9d634900847a
# ╠═c774a158-126b-11eb-2241-a9b6cf80c681
# ╠═345669c6-126e-11eb-2b12-a976e63e3335
# ╟─5a8ad83e-126e-11eb-2d85-6581c5379a4b
# ╠═b43cdce2-126e-11eb-0da0-fd2ec7d7bbad
# ╠═cccfbc70-126e-11eb-0439-4552285a634d
# ╠═f6e28740-1273-11eb-1464-ab04ba549059
# ╠═3460571e-126f-11eb-176b-a3445db55987
# ╠═08c5dc30-1276-11eb-01c1-4baf4dece1d5
# ╠═e7e4d35a-126f-11eb-34aa-dbb75f86c6a6
# ╠═80fcb26e-12a3-11eb-0977-3f37bdd5759e
# ╠═56671378-1271-11eb-1a8f-3f289512a146
# ╠═db81cb56-12a4-11eb-3b5b-910d159781fd
# ╠═316b4056-12a5-11eb-129f-bf07f3962a39
# ╠═a0f0629e-12a5-11eb-3d7e-c92f58bb13ad


text:
function right!(r)
    r = 0
    while (isborder(r, Ost) && !isborder(r,Nord))==true
        move!(r,Nord)
        r = r + 1
    end
    if !isborder(r,Ost)
        move!(r,Ost)
        while isborder(r,Sud)
            move!(r,Ost)
        end
    end
    while r > 0
        move!(r, Sud)
        r = r - 1
    end
end
function left!(r)
    l = 0
    while (isborder(r, West) && !isborder(r,Nord))==true
        move!(r, Nord)
        l = l + 1
    end
    if !isborder(r,West)
        move!(r,West)
        while isborder(r,Sud)
            move!(r,West)
        end
    end
    while l > 0
        move!(r, Sud)
        l = l - 1
    end
end
function mark_corss!(r::Robot)
    a = 0
    b = 0
    while (isborder(r,Nord) && isborder(r,West))==false
        while !isborder(r , Nord)
            move!(r, Nord)
            a = a + 1
        end
        while !isborder(r, West)
            move!(r,West)
            b = b + 1
        end
    end
    putmarker!(r)
    while (isborder(r,Sud) && isborder(r,Ost))==false
        if isborder(r,West)
            while !isborder(r, Ost)
                move!(r,Ost)
                putmarker!(r)
            end
            right!(r)
            if isborder(r, Ost)
                if !isborder(r, Sud)
                    move!(r, Sud)
                end
            end
        end
        if isborder(r,Ost)
            while !isborder(r, West)
                move!(r,West)
                putmarker!(r)
            end
            left!(r)
            if isborder(r, West)
                if !isborder(r, Sud)
                    move!(r, Sud)
                end
            end
        end
    end
    while !isborder(r,Nord)
        move!(r, Nord)
    end
    while !isborder(r,West)
        move!(r, West)
    end
    while b > 0
        move!(r, Ost)
        b = b - 1
    end
    while a > 0
        if isborder(r,Sud)
            s = 0
            while isborder(r, Sud)
                move!(r, West)
                s = s + 1
            end
            move!(r, Sud)
            a = a - 1
            while isborder(r, Ost)
                move!(r, Sud)
                a = a - 1
            end
            while s > 0
                move!(r, Ost)
                s = s - 1
            end
        else
            move!(r, Sud)
            a = a - 1
    end
end


text:
using DataStructures

# Struct representing a heap abstraction.
abstract type Heap end

# Struct representing the heap of nearest neighbors of a Point
mutable struct NN <: Heap
    # Ordered tree-like data structure, used to represent in a decreasing order the distances of the neighbours of the Point
    valtree::BinaryMaxHeap{Float64}
    # Maximum capacity of the heap
    capacity::Int64
    # Weight of the point, as the sum of the distances between the point and his neighbourhood
    weight::Float64
    function NN(k::T) where T<:Int64
        if k < 0
            error("Invalid cardinality (< 0)")
        end
        new(BinaryMaxHeap{Float64}(),k, 0.)
    end

    function NN()
        new(BinaryMaxHeap{Float64}(),0, 0.)
    end
end

function Base.show(io::IO, obj::NN)
    print(string(obj.valtree))
end

# Data structure representing the solution set
mutable struct MinHeap <: Heap
    # Ordered tree-like data structure, sorted by increasing point weight
    valtree::BinaryMinHeap
    # Minimum capacity of the heap
    capacity::Int64

    # Init MinHeap with min capacity = n
    function MinHeap(n::N) where {N<:Int64}
        if n < 0
            error("Invalid cardinality (< 0)")
        end
        new(BinaryMinHeap{Point}(),n)
    end
    function MinHeap()
        new(BinaryMinHeap{Point}(),0)
    end
end

function Base.show(io::IO, obj::Heap)
    print(string(obj.valtree))
end

# Insert the input distance into NN if the heap is not full, or if this distance is lower than the first element of the tree. The total weight is updated accordingly.
function updateMin(h::NN, d::Float64)
    if (length(h.valtree) < h.capacity)
        push!(h.valtree, d)
        h.weight += d
        return
    elseif d < top(h.valtree)
        h.weight -= pop!(h.valtree)
        push!(h.valtree,d)
        h.weight += d
        return
    end
end

# Insert the input element into MinHeap, if the heap is not full, or if the element’s weight is higher than the first element of the tree.
function updateMax(h::MinHeap, el::Element)
    if (length(h.valtree) < h.capacity)
        push!(h.valtree, el)
        return
    elseif el.nn.weight > top(h.valtree).nn.weight
        pop!(h.valtree)
        push!(h.valtree, el)
        return
    end
end

# get the first element of MinHeap
function getMin(h::MinHeap)
    return isempty(h.valtree) ? 0. : top(h.valtree).nn.weight
end


text:
using SnoopCompile, HomotopyContinuation

inf_timing = @snoopi include(joinpath(@__DIR__, "snoop.jl"))
pc = SnoopCompile.parcel(inf_timing)

dir = mktempdir()
@info "Writing precompile results into $dir"

SnoopCompile.write(dir, pc)

mv(joinpath(dir, "precompile_HomotopyContinuation.jl"),
    joinpath(@__DIR__, "../src/precompile.jl"); force = true)


text:
"""
    MatrixUtilityModule

Module for general utility matrix product functions.
"""
module MatrixUtilityModule

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict

"""
    loc!(loc::FFltMat, ecoords::FFltMat, N::FFltMat)

Compute the location of the quadrature point.

Arguments:
`loc` = matrix of coordinates, overwritten  inside the function
`ecoords` = matrix of the node coordinates for the element.
`N` = matrix of basis function values
"""
function loc!(loc::FFltMat, ecoords::FFltMat, N::FFltMat)
    n = size(N, 1)
    @inbounds for j in 1:size(loc, 2)
        la = 0.0
        @inbounds for k in 1:n
            la += N[k] * ecoords[k, j]
        end
        loc[j] = la
    end
    return loc
end

"""
    jac!(J::FFltMat, ecoords::FFltMat,gradNparams::FFltMat)

Compute the Jacobian matrix at the quadrature point.

Arguments:
`J` = Jacobian matrix, overwritten  inside the function
`ecoords` = matrix of the node coordinates for the element.
`gradNparams` = matrix of basis function gradients
"""
function jac!(J::FFltMat, ecoords::FFltMat, gradNparams::FFltMat)
    n = size(gradNparams, 1)
    @inbounds for j in 1:size(J, 2)
        @inbounds for i in 1:size(J, 1)
            Ja = 0.0
            @inbounds for k in 1:n
                Ja += ecoords[k, i] * gradNparams[k, j]
            end
            J[i, j] = Ja
        end
    end
    return J
end

"""
    locjac!(loc::FFltMat, J::FFltMat, ecoords::FFltMat, N::FFltMat, gradNparams::FFltMat)

Compute location and Jacobian matrix at the quadrature point.

Arguments:
`loc` = matrix of coordinates, overwritten  inside the function
`J` = Jacobian matrix, overwritten  inside the function
`ecoords` = matrix of the node coordinates for the element.
`N` = matrix of basis function values
`gradNparams` = matrix of basis function gradients
"""
function locjac!(loc::FFltMat, J::FFltMat, ecoords::FFltMat, N::FFltMat, gradNparams::FFltMat)
    return loc!(loc, ecoords, N), jac!(J, ecoords, gradNparams)
end

"""
    add_mggt_ut_only!(Ke::FFltMat, gradN::FFltMat, mult::FFlt)

Add the product `gradN*mult*gradNT*(Jac*w[j])` to the elementwise matrix `Ke`.
The argument `mult` is a scalar.
*Only upper triangle* is computed; the lower triangle is not touched.
(Use `complete_lt!` to complete the lower triangle, if needed.)

The matrix `Ke` is assumed to be suitably initialized.

The matrix `Ke` is modified.  The matrix `gradN` is not modified
inside this function.
"""
function add_mggt_ut_only!(Ke::FFltMat, gradN::FFltMat, mult::FFlt)
    Kedim = size(Ke, 1)
    @assert Kedim == size(Ke, 2) # Square matrix?
    nne, mdim = size(gradN)
    @assert nne == Kedim # compatible matrices?
    @inbounds for nx = 1:Kedim # Do: Ce  =  Ce + gradN*((Jac*w[j]))*gradN' ;
        @inbounds for px = 1:mdim
            a = (mult)*gradN[nx, px]
            @inbounds for mx = 1:nx # only the upper triangle
                Ke[mx, nx] +=  gradN[mx, px] * a
            end
        end
    end
    return true
end

"""
    add_gkgt_ut_only!(Ke::FFltMat, gradN::FFltMat, Jac_w::FFlt,
      kappa_bar::FFltMat, kappa_bargradNT::FFltMat)

Add the product `gradN*kappa_bar*gradNT*(Jac*w[j])` to the elementwise matrix `Ke`.
*Only upper triangle* is computed; the lower triangle is not touched.
(Use `complete_lt!` to complete the lower triangle, if needed.)

The matrix `Ke` is assumed to be suitably initialized.

Upon return,  the matrix `Ke` is updated.  The scratch buffer `kappa_bargradNT`
is overwritten during each call of this function. The matrices `gradN` and
`kappa_bar` are not modified inside this function.
"""
function add_gkgt_ut_only!(Ke::FFltMat, gradN::FFltMat, Jac_w::FFlt,
        kappa_bar::FFltMat, kappa_bargradNT::FFltMat)
    @assert size(Ke, 1) == size(Ke, 2)
    Kedim = size(Ke, 1)
    nne, mdim = size(gradN)
    @assert size(kappa_bar) == (mdim, mdim)
    @assert size(kappa_bargradNT) == (mdim, nne)
    # A_mul_Bt!(kappa_bargradNT, kappa_bar, gradN); # intermediate result
    @inbounds for nx = 1:nne
        @inbounds for mx = 1:mdim
	        accum::FFlt  = 0.0
	        @inbounds @simd for px = 1:mdim
	            accum += kappa_bar[mx, px]*gradN[nx, px]
	        end
	        kappa_bargradNT[mx, nx] = Jac_w*accum
        end
    end
    # Ke = Ke + gradN*(kappa_bar*(Jac*w[j]))*gradN'; only upper triangle
    @inbounds for nx = 1:Kedim
        @inbounds for mx = 1:nx # only the upper triangle
            accum::FFlt  = 0.0
            @inbounds @simd for px = 1:mdim
                accum += gradN[mx, px]*kappa_bargradNT[px, nx]
            end
            Ke[mx, nx] += accum
        end
    end
    return true
end

"""
    complete_lt!(Ke::FFltMat)

Complete the lower triangle of the elementwise matrix `Ke`.

The matrix `Ke` is modified  inside this function. The
upper-triangle  entries  are copied  across the diagonal
to the lower triangle.
"""
function complete_lt!(Ke::FFltMat)
    Kedim = size(Ke, 1)
    @assert Kedim == size(Ke, 2)
    @inbounds for nx = 1:Kedim # complete the lower triangle
        @inbounds for mx = nx+1:Kedim
            Ke[mx, nx] = Ke[nx, mx]
        end
    end
    return true
end

"""
    add_btdb_ut_only!(Ke::FFltMat, B::FFltMat, Jac_w::FFlt,
                  D::FFltMat, DB::FFltMat)

Add the product  `(B'*(D*(Jac*w[j]))*B)`, to the elementwise matrix Ke.
*Only upper triangle* is computed; the lower triangle is not touched.
(Use `complete_lt!` to complete the lower triangle, if needed.)

The matrix Ke is assumed to be suitably initialized.

The matrix Ke is modified.  The matrices B and D are not modified
inside this function. The scratch buffer DB is overwritten
during each call of this function.
"""
function add_btdb_ut_only!(Ke::FFltMat, B::FFltMat, Jac_w::FFlt, D::FFltMat, DB::FFltMat)
    @assert size(Ke, 1) == size(Ke, 2)
    @assert size(B, 1) == size(D, 1)
    nstr, Kedim = size(B)
    @assert size(D) == (nstr, nstr)
    @assert size(DB) == (nstr, Kedim)
    # A_mul_B!(DB, D, B)    # intermediate product
    @inbounds for nx = 1:Kedim
        @inbounds for mx = 1:nstr
            accum::FFlt  = 0.0
            @inbounds for px = 1:nstr
                accum += D[mx, px]*B[px, nx]
            end
            DB[mx, nx] = Jac_w*accum
        end
    end
    #  Ke = Ke + (B'*(D*(Jac*w[j]))*B); only the upper triangle
    @inbounds for nx = 1:Kedim
        @inbounds for mx = 1:nx # only the upper triangle
	        accum::FFlt  = 0.0
	        @inbounds for px = 1:nstr
	            accum += B[px, mx]*DB[px, nx]
	        end
	        Ke[mx, nx] += accum
        end
    end
    return true
end

"""
    add_nnt_ut_only!(Ke::FFltMat, Nn::FFltMat, Jac_w_coeff::FFlt)

Add the product  `Nn*(Nn'*(coeff*(Jac*w(j)))`, to the elementwise matrix `Ke`.
*Only upper triangle* is computed; the lower triangle is not touched.

The matrix `Ke` is assumed to be suitably initialized.

The matrix `Ke` is modified.  The matrix `Nn` is not modified
inside this function.
"""
function add_nnt_ut_only!(Ke::FMat{T}, N::FFltMat, Jac_w_coeff::T) where {T<:Number}
    Kedim = size(N, 1)
    @assert Kedim == size(Ke, 2)
    @inbounds for nx = 1:Kedim
        a = (Jac_w_coeff)*N[nx]
        @inbounds for mx = 1:nx # only the upper triangle
            Ke[mx,nx] += N[mx] * a
        end
    end
    return true
end

function add_nnt_ut_only!(Ke::FFltMat, N::FFltVec, Jac_w_coeff::FFlt)
    return add_nnt_ut_only!(Ke, reshape(N, length(N), 1), Jac_w_coeff)
end

"""
    add_btv!(elvec::FFltVec, B::FFltMat, sig::FFltMat, Jac_w_coeff::FFlt)

Add product of the strain-displacement matrix transpose times the stress vector.

Note:  the coefficient `Jac_w_coeff` will be typically  NEGATIVE.

The argument `elvec` needs to be suitably  initialized before the first call
(filled with zeros, for instance),  and it is updated upon return.
The arguments `B`, `sig` are not modified.
"""
function add_btv!(elvec::FFltVec, B::FFltMat, sig::FFltVec, Jac_w_coeff::FFlt)
    @assert size(B, 1) == length(sig)
    @assert size(B, 2) == length(elvec)
    Kedim = length(elvec)
    nsig = length(sig)
    @inbounds for mx = 1:Kedim
        @inbounds for px = 1:nsig
            elvec[mx] += B[px, mx]*Jac_w_coeff*sig[px]
        end
    end
    return true
end

"""
    adjugate3!(B, A)

Compute the adjugate matrix of `A`.
"""
function adjugate3!(B, A)
    # % A = rand(3)
    # % det(A)*inv(A)-B
    # B =0*A;
    B[1,1] =(A[5]*A[9]-A[8]*A[6]);
    B[1,2] =-(A[4]*A[9]-A[7]*A[6]);
    B[1,3] =(A[4]*A[8]-A[7]*A[5]);

    B[2,1] =-(A[2]*A[9]-A[8]*A[3]);
    B[2,2] =(A[1]*A[9]-A[7]*A[3]);
    B[2,3] =-(A[1]*A[8]-A[7]*A[2]);

    B[3,1] =(A[2]*A[6]-A[5]*A[3]);
    B[3,2] =-(A[1]*A[6]-A[4]*A[3]);
    B[3,3] =(A[1]*A[5]-A[4]*A[2]);
    return B
end

"""
    symmetrize!(a)

Make the matrix on input symmetric.
"""
function symmetrize!(a)
	@assert size(a, 1) == size(a, 2)
    @inbounds for c in 1:size(a, 2)
    	@inbounds for r in c:size(a, 1)
    		a[c, r] =  a[r, c] += a[c, r]
    	end
    end
    a .= 0.5 * a
end

end


text:
<%!
    import shutil
    from logging import getLogger

    import onlinejudge_template.generator.julia as julia
    import onlinejudge_template.generator.about as about
    # import onlinejudge_template.generator.hook as hook
%>\
<%
    logger = getLogger(__name__)
    if not shutil.which("yapf"):
        logger.warning("yapf is not installed. If you want to generate well-formatted code, please install it. You can use $ pip3 install yapf")
    else:
        format_config = "{" + ", ".join([
            "BASED_ON_STYLE: google",
            "COLUMN_LIMIT: 9999",
        ]) + "}"
        # hook.register_filter_command(["yapf", "--style", format_config], data=data)
%>\
#!/usr/bin/env julia

${julia.declare_constants(data)}

# def solve(${julia.formal_arguments(data)}) -> ${julia.return_type(data)}:
function solve(${julia.formal_arguments(data, typed=False)})
    pass  # TODO: edit here

end

# generated by ${about.title} ${about.version} (${about.url})
function main()
${julia.read_input(data)}
    ${julia.return_value(data)} = solve(${julia.actual_arguments(data)})
${julia.write_output(data)}

end

main()

max_stars_repo_name: kameyama/template-generator

max_issues_repo_name: kameyama/template-generator

max_forks_repo_name: kameyama/template-generator


text:
using DecisionTreePerformance
using Base.Test

function assert_all_eq(elements)
    head = elements[1]
    for el in elements
        @test head == el
    end
end

function check_all_equal(features::Vector{Float64}, evaluators::Any)
    map(p -> evaluate(p[2], features), evaluators) |> assert_all_eq
end

function test_equal_evaluations(num_trees, num_features, depth, iterations)

    evaluators = random_forest(num_trees, num_features, depth) |> construct_evaluators

    # QuickCheck style test
    for i in 1:iterations
        fv = random_feature_vector(num_features)
        check_all_equal(fv, evaluators)
    end
end

function bench_evaluation(feature_vector, evaluator, iterations)
    @profile begin
        for i in 1:iterations
            evaluate(evaluator, feature_vector)
        end
    end
end

test_equal_evaluations(5, 5, 2, 10)


text:
get_method(f, args::Tuple) = get_method(f, Tuple{args...})
function get_method(f::Type{F}, args::Type{<:Tuple}) where F<:Function
    return get_method(F.instance, args)
end
get_method(f, args::Type{<:Tuple}) = Base.which(f, args)

function get_ast(meth::Method)
    ast = definition(meth)
    return MacroTools.striplines(ast)
end

functiontypeof(m::Method) = parameter_typeof(m.sig)[1]
parameter_typeof(sig::UnionAll) = parameter_typeof(sig.body)
parameter_typeof(sig::DataType) = sig.parameters

function determine_type_params(
        meth::Method, argtypes::Type{<:Tuple},
    )
    ftype = functiontypeof(meth)

    instance_sig = Tuple{ftype, argtypes.parameters...}
    ti, lenv = ccall(
        :jl_type_intersection_with_env,
        Any, (Any, Any),
        instance_sig, meth.sig)
    return lenv
end


text:
#
# Input types
#
abstract type ABCInput end

abstract type ABCRejectionInput <: ABCInput end

RepetitiveTraining(; rt_iterations::Int64=0, rt_extra_training_points::Int64=1, rt_sample_size::Int64=1000) =
    RepetitiveTraining(rt_iterations, rt_extra_training_points, rt_sample_size)

struct DistanceSimulationInput
    reference_summary_statistic::AbstractArray{Float64,1}
    simulator_function::Function
    summary_statistic::Function
    distance_metric::Function
end

"""
	AbstractEmulatorRetraining

Subtypes of this abstract type control the additional retraining procedure that may or may not be carried out before each iteration of emulation-based ABC. Custom strategies may be implemented by creating new subtypes of this type and new [`abc_retrain_emulator`](@ref) methods for them.

The following implementations are shipped:
- [`IncrementalRetraining`](@ref)
- [`PreviousPopulationRetraining`](@ref)
- [`PreviousPopulationThresholdRetraining`](@ref)
- [`NoopRetraining`](@ref)
"""
abstract type AbstractEmulatorRetraining end

"""
	IncrementalRetraining <: AbstractEmulatorRetraining

This emulator retraining strategy samples extra particles from the previous population, and adds them to the set of design points that were used on the previous iteration. The new design points are filtered according to the threshold. The combined set of design points is used to train the new emulator.

# Fields
- `design_points`: number of design points to add on each iteration
- `max_simulations`: maximum number of simulations to perform during re-trainging on each iteration
"""
struct IncrementalRetraining <: AbstractEmulatorRetraining
    design_points::Int
    max_simulations::Int
end

"""
	PreviousPopulationRetraining <: AbstractEmulatorRetraining

This emulator retraining strategy samples extra particles from the previous population, and uses them to re-train the emulator from scratch. No filtering of the new design points is performed. Design points from the previous iteration are discarded.
"""
struct PreviousPopulationRetraining <: AbstractEmulatorRetraining end

"""
	PreviousPopulationThresholdRetraining <: AbstractEmulatorRetraining

This emulator retraining strategy samples extra particles from the previous population, and uses them to re-train the emulator from scratch. Design points from the previous iteration are discarded. This strategy allows to control how many design points are sampled with distance to the reference data below the threshold.

# Fields:
- `n_design_points`: number of design points
- `n_below_threshold`: number of design points below the threshold
- `max_iter`: maximum number of simulations to perform on each re-training iteration
"""
struct PreviousPopulationThresholdRetraining <: AbstractEmulatorRetraining
    n_design_points::Int
    n_below_threshold::Int
    max_iter::Int
end

"""
	NoopRetraining <: AbstractEmulatorRetraining

A sentinel retraining strategy that does not do anything. When used, the emulator is trained only once before
"""
struct NoopRetraining <: AbstractEmulatorRetraining end

"""
    AbstractEmulatorTraining

Subtypes of this abstract type control how the emulator is trained for emulation-based ABC algorithms (rejection and SMC).
At the moment, only [`DefaultEmulatorTraining`](@ref) is shipped. Custom emulator training procedure
can be implemented by creating new subtypes of this type and overriding [`train_emulator`](@ref) for them.

A typical use case would be trying to control the behaviour of [`gp_train`](@ref) more tightly,
or not using it altogeather (e.g. using another optimisation package).
"""
abstract type AbstractEmulatorTraining end

"""
	DefaultEmulatorTraining <: AbstractEmulatorTraining

# Fields
- `kernel::AbstractGPKernel`: the kernel ([`AbstractGPKernel`](@ref)) that will be used with the Gaussian Process ([`GPModel`](@ref)). Defaults to [`SquaredExponentialArdKernel`](@ref).

[`train_emulator`](@ref) method with this argument type calls [`gp_train`](@ref) with default arguments.
"""
struct DefaultEmulatorTraining{K<:AbstractGPKernel} <: AbstractEmulatorTraining
    kernel::K
end
DefaultEmulatorTraining() = DefaultEmulatorTraining(SquaredExponentialArdKernel())

struct EmulatorTrainingInput{ET<:AbstractEmulatorTraining}
    distance_simulation_input::DistanceSimulationInput
    design_points::Int64
    emulator_training::ET
end
EmulatorTrainingInput(dsi::DistanceSimulationInput) = EmulatorTrainingInput(dsi, DefaultEmulatorTraining())
EmulatorTrainingInput(n_design_points, reference_summary_statistic, simulator_function, summary_statistic, distance_metric, et=DefaultEmulatorTraining()) =
    EmulatorTrainingInput(DistanceSimulationInput(
        reference_summary_statistic, simulator_function,
        build_summary_statistic(summary_statistic), distance_metric),
        n_design_points, et)

"""
	AbstractEmulatedParticleSelection

Subtypes of this type control the criteria that determine what particles are included in the posterior for emulation-based ABC. Custom strategies
can be implemented by creating new subtypes of this type and overriding [`abc_select_emulated_particles`](@ref) for them.

Two implementations are shipped:
- [`MeanEmulatedParticleSelection`](@ref)
- [`MeanVarEmulatedParticleSelection`](@ref)
"""
abstract type AbstractEmulatedParticleSelection end

"""
   MeanEmulatedParticleSelection <: AbstractEmulatedParticleSelection

When this strategy is used, the particles for which only the *mean* value returned by
[`gp_regression`](@ref) is below the ABC threshold are included in the posterior.
Variance is not checked.
"""
struct MeanEmulatedParticleSelection <: AbstractEmulatedParticleSelection end

"""
   MeanVarEmulatedParticleSelection <: AbstractEmulatedParticleSelection

When this strategy is used, the particles for which both *mean and standard deviation* returned by
[`gp_regression`](@ref) is below the ABC threshold are included in the posterior.

# Fields
- `variance_threshold_factor`: scaling factor, by which the ABC threshold is multiplied
before checking the standard deviation. Defaults to 1.0.
"""
struct MeanVarEmulatedParticleSelection <: AbstractEmulatedParticleSelection
    variance_threshold_factor::Float64
end
MeanVarEmulatedParticleSelection() = MeanVarEmulatedParticleSelection(1.0)

struct SimulatedABCRejectionInput <: ABCRejectionInput
    n_params::Int64
    n_particles::Int64
    threshold::Float64
    priors::AbstractArray{ContinuousUnivariateDistribution,1}
    distance_simulation_input::DistanceSimulationInput
    max_iter::Int
end

struct EmulatedABCRejectionInput{CUD<:ContinuousUnivariateDistribution, EPS<:AbstractEmulatedParticleSelection} <: ABCRejectionInput
	n_params::Int64
	n_particles::Int64
	threshold::Float64
	priors::AbstractArray{CUD,1}
	batch_size::Int64
    max_iter::Int64
    emulator_training_input::EmulatorTrainingInput
    selection::EPS
end

abstract type ABCSMCInput <: ABCInput end

struct SimulatedABCSMCInput <: ABCSMCInput
    n_params::Int64
    n_particles::Int64
    threshold_schedule::AbstractArray{Float64,1}
    priors::AbstractArray{ContinuousUnivariateDistribution,1}
    distance_simulation_input::DistanceSimulationInput
    max_iter::Int
end

SimulatedABCRejectionInput(smc_input::SimulatedABCSMCInput) =
    SimulatedABCRejectionInput(smc_input.n_params,
        smc_input.n_particles,
        smc_input.threshold_schedule[1],
        smc_input.priors,
        smc_input.distance_simulation_input,
        smc_input.max_iter)

struct EmulatedABCSMCInput{CUD<:ContinuousUnivariateDistribution,
        ER<:AbstractEmulatorRetraining, EPS<:AbstractEmulatedParticleSelection} <: ABCSMCInput
    n_params::Int64
    n_particles::Int64
    threshold_schedule::AbstractArray{Float64,1}
    priors::AbstractArray{CUD,1}
    batch_size::Int64
    max_iter::Int64
    emulator_training_input::EmulatorTrainingInput
    emulator_retraining::ER
    selection::EPS
end

#
# Tracker types
#
abstract type ABCSMCTracker end

mutable struct SimulatedABCSMCTracker <: ABCSMCTracker
    n_params::Int64
    n_accepted::AbstractArray{Int64,1}
    n_tries::AbstractArray{Int64,1}
    threshold_schedule::AbstractArray{Float64,1}
	can_continue::Bool # we used to return this flag from iterateABCSMC!, but it got broken in julia 1.0 - see https://github.com/JuliaLang/julia/issues/29805
    population::AbstractArray{AbstractArray{Float64,2},1}
    distances::AbstractArray{AbstractArray{Float64,1},1}
    weights::AbstractArray{StatsBase.Weights,1}
    priors::AbstractArray{ContinuousUnivariateDistribution,1}
    distance_simulation_input::DistanceSimulationInput
    max_iter::Int
end

mutable struct EmulatedABCSMCTracker{CUD<:ContinuousUnivariateDistribution, ET,
        ER<:AbstractEmulatorRetraining, EPS<:AbstractEmulatedParticleSelection} <: ABCSMCTracker
    n_params::Int64
    n_accepted::AbstractArray{Int64,1}
    n_tries::AbstractArray{Int64,1}
    threshold_schedule::AbstractArray{Float64,1}
	can_continue::Bool # we used to return this flag from iterateABCSMC!, but it got broken in julia 1.0 - see https://github.com/JuliaLang/julia/issues/29805
    population::AbstractArray{AbstractArray{Float64,2},1}
    distances::AbstractArray{AbstractArray{Float64,1},1}
    weights::AbstractArray{StatsBase.Weights,1}
    priors::AbstractArray{CUD,1}
    emulator_training_input::EmulatorTrainingInput
    emulator_retraining::ER
    selection::EPS
    batch_size::Int64
    max_iter::Int64
    emulators::AbstractArray{ET,1}
end

#
# Output types
#
abstract type ABCOutput end

"""
    ABCRejectionOutput

A container for the output of a rejection-ABC computation.

# Fields
- `n_params::Int64`: The number of parameters to be estimated.
- `n_accepted::Int64`: The number of accepted parameter vectors (particles) in the posterior.
- `n_tries::Int64`: The total number of parameter vectors (particles) that were tried.
- `threshold::Float64`: The maximum distance from the summarised model output to summarised observed data for a parameter vector to be included in the posterior.
- `population::AbstractArray{Float64,2}`: The parameter vectors (particles) in the posterior. Size: (`n_accepted`, `n_params`).
- `distances::AbstractArray{Float64,1}`: The distances for each parameter vector (particle) in the posterior to the observed data in summary statistic space. Size: (`n_accepted`).
- `weights::StatsBase.Weights`: The weight of each parameter vector (particle) in the posterior.
"""
abstract type ABCRejectionOutput <: ABCOutput end

struct EmulatedABCRejectionOutput{ET} <: ABCRejectionOutput
    n_params::Int64
    n_accepted::Int64
    n_tries::Int64
    threshold::Float64
    population::AbstractArray{Float64,2}
    distances::AbstractArray{Float64,1}
    weights::StatsBase.Weights
    emulator::ET
end

struct SimulatedABCRejectionOutput <: ABCRejectionOutput
    n_params::Int64
    n_accepted::Int64
    n_tries::Int64
    threshold::Float64
    population::AbstractArray{Float64,2}
    distances::AbstractArray{Float64,1}
    weights::StatsBase.Weights
end

"""
    ABCSMCOutput

A container for the output of a rejection-ABC computation.

# Fields
- `n_params::Int64`: The number of parameters to be estimated.
- `n_accepted::Int64`: The number of accepted parameter vectors (particles) in the posterior.
- `n_tries::Int64`: The total number of parameter vectors (particles) that were tried.
- `threshold_schedule::AbstractArray{Float64}`: A set of maximum distances from the summarised model output to summarised observed data for a parameter vector to be included in the posterior. Each distance will be used in a single run of the ABC-SMC algorithm.
- `population::AbstractArray{Float64,2}`: The parameter vectors (particles) in the posterior. Size: (`n_accepted`, `n_params`).
- `distances::AbstractArray{Float64,1}`: The distances for each parameter vector (particle) in the posterior to the observed data in summary statistic space. Size: (`n_accepted`).
- `weights::StatsBase.Weights`: The weight of each parameter vector (particle) in the posterior.
"""
abstract type ABCSMCOutput <: ABCOutput end

struct SimulatedABCSMCOutput <: ABCSMCOutput
    n_params::Int64
    n_accepted::AbstractArray{Int64,1}
    n_tries::AbstractArray{Int64,1}
    threshold_schedule::AbstractArray{Float64,1}
    population::AbstractArray{AbstractArray{Float64,2},1}
    distances::AbstractArray{AbstractArray{Float64,1},1}
    weights::AbstractArray{StatsBase.Weights,1}
end

struct EmulatedABCSMCOutput <: ABCSMCOutput
    n_params::Int64
    n_accepted::AbstractArray{Int64,1}
    n_tries::AbstractArray{Int64,1}
    threshold_schedule::AbstractArray{Float64,1}
    population::AbstractArray{AbstractArray{Float64,2},1}
    distances::AbstractArray{AbstractArray{Float64,1},1}
    weights::AbstractArray{StatsBase.Weights,1}
    emulators::AbstractArray{GPModel,1}
end

#
# Read/write functions
#
function write_scalar(
        stream::IO,
        x::Real;
        separator="\n",
        )
    write(stream, string(x), separator)
end

function write_vector(
        stream::IO,
        vector::Vector{R};
        separator=" ",
        finalcharacter="\n"
        ) where {
        R<:Real,
        }
    n = length(vector)
    for i in 1:(n - 1)
        write_scalar(stream, vector[i]; separator=separator)
    end
    write_scalar(stream, vector[n]; separator=finalcharacter)
end

function write_matrix(
        stream::IO,
        matrix::Matrix{R};
        separator=" ",
        newline=true,
        ) where {
        R<:Real,
        }
    n = size(matrix, 2)
    for i in 1:(n - 1)
        write_vector(stream, matrix[:, i]; separator=separator, finalcharacter="\n")
    end
    finalcharacter = newline ? "\n" : ""
    write_vector(stream, matrix[:, n]; separator=separator, finalcharacter=finalcharacter)
end

function Base.write(stream::IO, output::ABCRejectionOutput)
    write_scalar(stream, output.n_params)
    write_scalar(stream, output.n_accepted)
    write_scalar(stream, output.n_tries)
    write_scalar(stream, output.threshold)
    write_matrix(stream, output.population)
    write_vector(stream, output.distances)
    write_vector(stream, output.weights.values)

    flush(stream)
end

function Base.write(stream::IO, output::Union{ABCSMCOutput, ABCSMCTracker})
    write_scalar(stream, output.n_params)
    write_vector(stream, output.n_accepted)
    write_vector(stream, output.n_tries)
    write_vector(stream, output.threshold_schedule)
    for i in 1:length(output.threshold_schedule)
        write(stream, "-----\n")
        write_scalar(stream, output.n_accepted[i])
        write_scalar(stream, output.n_tries[i])
        write_scalar(stream, output.threshold_schedule[i])
        write_matrix(stream, output.population[i])
        write_vector(stream, output.distances[i])
        write_vector(stream, output.weights[i].values)
    end

    flush(stream)
end

function write(stream::IO, output::Union{ABCOutput, ABCSMCTracker})
    return Base.write(stream, output)
end

function read_rejection_output(filepath::AbstractString)
    input_file = open(filepath, "r")

    try
        n_params = parse(Int64, readline(input_file))
        n_accepted = parse(Int64, readline(input_file))
        n_tries = parse(Int64, readline(input_file))
        threshold = parse(Float64, readline(input_file))

        population = zeros(n_params, n_accepted)
        for i in 1:n_accepted
            population[:, i] = parse.(Float64, split(readline(input_file)))
        end

        distances = parse.(Float64, split(readline(input_file)))

        wts = parse.(Float64, split(readline(input_file)))
        weights = StatsBase.Weights(wts)

        return ABCRejectionOutput(n_params,
                                  n_accepted,
                                  n_tries,
                                  threshold,
                                  population,
                                  distances,
                                  weights,
                                  )
    finally
        close(input_file)
    end
end

function read_smc_output(filepath::AbstractString)
    input_file = open(filepath, "r")

    try
        n_params = parse(Int64, readline(input_file))
        n_accepted = parse.(Int64, split(readline(input_file)))
        n_tries = parse.(Int64, split(readline(input_file)))
        threshold_schedule = parse.(Float64, split(readline(input_file)))

        @assert length(n_accepted) == length(n_tries) == length(threshold_schedule)

        population = Matrix{Float64}[]
        distances = Vector{Float64}[]
        weights = StatsBase.Weights[]

        for i in 1:length(threshold_schedule)
            separator = readline(input_file)

            @assert parse(Int64, readline(input_file)) == n_accepted[i]
            @assert parse(Int64, readline(input_file)) == n_tries[i]
            @assert parse(Float64, readline(input_file)) == threshold_schedule[i]

            push!(population, zeros(n_params, n_accepted[i]))
            for j in 1:n_accepted[i]
                particle = parse.(Float64, split(readline(input_file)))
                population[i][:, j] = particle
            end

            push!(distances, parse.(Float64, split(readline(input_file))))

            wts = parse.(Float64, split(readline(input_file)))
            push!(weights, StatsBase.Weights(wts))
        end

        return ABCSMCOutput(n_params,
                            n_accepted,
                            n_tries,
                            threshold_schedule,
                            population,
                            distances,
                            weights,
                            )
    finally
        close(input_file)
    end
end

function checkABCInput(input::ABCInput)
	if input.n_params != length(input.priors)
		throw(ArgumentError("There are $(input.n_params) unknown parameters but $(length(input.priors)) priors were provided"))
	end
end


text:
#=
FloydWarshall:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-10
=#

using SuiteSparseGraphBLAS
using SparseArrays

#A is the input matrix for Floyd-Warshall multiple sources multiple destinations, Output: shortest paths
function floydwarshall(A, n)

    #distance matrix = input matrix
    D = A
    print(D)
    result = GBMatrix{Int64}(n, n)

    for i = 1:n
            print("\n\n\n")
            partial = mul(D[:,i], D[i,:], Semirings.MIN_PLUS)
            D = emul(D, partial,  BinaryOps.MIN)
            print(partial)
            result = partial
    end

    return result

end


text:
__precompile__()
module FINUFFT

## Export
export nufft1d1, nufft1d2, nufft1d3
export nufft2d1, nufft2d2, nufft2d3
export nufft3d1, nufft3d2, nufft3d3

export nufft1d1!, nufft1d2!, nufft1d3!
export nufft2d1!, nufft2d2!, nufft2d3!
export nufft3d1!, nufft3d2!, nufft3d3!

export nufft_opts
export nufft_c_opts        # backward-compatibility - remove?
export finufft_plan
export finufft_default_opts
export finufft_makeplan
export finufft_setpts
export finufft_exec
export finufft_destroy
export finufft_exec!
export BIGINT
export finufftReal

## External dependencies
using finufft_jll

const libfinufft = finufft_jll.libfinufft

const BIGINT = Int64        # must match that in FINUFFT include/dataTypes.h

# our allowed real array types...
const finufftReal = Union{Float64,Float32}

## FINUFFT opts struct, must bytewise match that in include/nufft_opts.h
"""
    mutable struct nufft_opts
        modeord            :: Cint
        chkbnds            :: Cint
        debug              :: Cint
        spread_debug       :: Cint
        showwarn           :: Cint
        nthreads           :: Cint
        fftw               :: Cint
        spread_sort        :: Cint
        spread_kerevalmeth :: Cint
        spread_kerpad      :: Cint
        upsampfac          :: Cdouble
        spread_thread      :: Cint
        maxbatchsize       :: Cint
        spread_nthr_atomic :: Cint
        spread_max_sp_size :: Cint
    end

Options struct passed to the FINUFFT library.

# Fields

This is a summary only; see FINUFFT documentation for full descriptions.

    modeord :: Cint
0: CMCL-style increasing mode ordering (neg to pos), or\\
1: FFT-style mode ordering (affects type-1,2 only)

    chkbnds :: Cint
0: don't check if input NU pts in [-3pi,3pi], 1: do

    debug :: Cint
0: silent, 1: text basic timing output

    spread_debug :: Cint
passed to spread_opts, 0 (no text) 1 (some) or 2 (lots)

    showwarn :: Cint
Whether to print warnings to stderr. 0: silent, 1: print warnings

    nthreads :: Cint
How many threads FINUFFT should use, or 0 (use max available in OMP)

    fftw :: Cint
0:`FFTW_ESTIMATE`, or 1:`FFTW_MEASURE` (slow plan but faster FFTs)

    spread_sort :: Cint
passed to spread_opts, 0 (don't sort) 1 (do) or 2 (heuristic)

    spread_kerevalmeth :: Cint
passed to spread_opts, 0: exp(sqrt()), 1: Horner ppval (faster)

    spread_kerpad :: Cint
passed to spread_opts, 0: don't pad to mult of 4, 1: do

    upsampfac :: Cdouble
upsampling ratio sigma: 2.0 (standard), or 1.25 (small FFT), or\\
0.0 (auto).

    spread_thread :: Cint
(for ntrans>1 only)\\
0: auto choice,\\
1: sequential multithreaded,\\
2: parallel singlethreaded spread.

    maxbatchsize :: Cint
(for ntrans>1 only). max blocking size for vectorized, 0 for auto-set

    spread_nthr_atomic :: Cint
if >=0, threads above which spreader OMP critical goes atomic

    spread_max_sp_size :: Cint
if >0, overrides spreader (dir=1 only) max subproblem size
"""
mutable struct nufft_opts{T}
    modeord            :: Cint
    chkbnds            :: Cint
    debug              :: Cint
    spread_debug       :: Cint
    showwarn           :: Cint
    nthreads           :: Cint
    fftw               :: Cint
    spread_sort        :: Cint
    spread_kerevalmeth :: Cint
    spread_kerpad      :: Cint
    upsampfac          :: Cdouble        # regardless of T
    spread_thread      :: Cint
    maxbatchsize       :: Cint
    spread_nthr_atomic :: Cint
    spread_max_sp_size :: Cint
    nufft_opts{T}() where T <: finufftReal = new{T}()
end
# The above must match include/nufft_opts.h in FINUFFT.
# Note that nufft_opts happens to be the same for Float32 vs Float64 types.
# This is as in FINUFFT. It may not always hold in the future.

const nufft_c_opts = nufft_opts        # for backward compatibility - remove?

"""
    p = finufft_default_opts()
    p = finufft_default_opts(dtype=Float32)

Return a [`nufft_opts`](@ref) struct with the default FINUFFT settings. Set up the double precision variant by default.\\
See: <https://finufft.readthedocs.io/en/latest/usage.html#options>
"""
function finufft_default_opts(dtype::DataType=Float64)
    opts = nufft_opts{dtype}()

    if dtype==Float64
        ccall( (:finufft_default_opts, libfinufft),
               Nothing,
               (Ref{nufft_opts},),
               opts
               )

    else
        ccall( (:finufftf_default_opts, libfinufft),
               Nothing,
               (Ref{nufft_opts},),
               opts
               )
    end

    return opts
end

### Error handling

# Following should match error codes in https://github.com/flatironinstitute/finufft/blob/master/include/defs.h
const ERR_EPS_TOO_SMALL        = 1
const ERR_MAXNALLOC            = 2
const ERR_SPREAD_BOX_SMALL     = 3
const ERR_SPREAD_PTS_OUT_RANGE = 4
const ERR_SPREAD_ALLOC         = 5
const ERR_SPREAD_DIR           = 6
const ERR_UPSAMPFAC_TOO_SMALL  = 7
const HORNER_WRONG_BETA        = 8
const ERR_NDATA_NOTVALID       = 9
const ERR_TYPE_NOTVALID        = 10
const ERR_ALLOC                = 11
const ERR_DIM_NOTVALID         = 12
const ERR_SPREAD_THREAD_NOTVALID = 13

struct FINUFFTError <: Exception
    errno::Cint
    msg::String
end
Base.showerror(io::IO, e::FINUFFTError) = print(io, "FINUFFT Error ($(e.errno)): ", e.msg)

function check_ret(ret)
    # Check return value and output corresponding error message to Julia side.
    # This should be kepts up to date with error code interpretation in
    # https://finufft.readthedocs.io/en/latest/error.html
    # which is generated from
    # https://github.com/flatironinstitute/finufft/blob/master/docs/error.rst
    if ret==0     # no error or warning
        return
    elseif ret==ERR_EPS_TOO_SMALL
        msg = "requested tolerance epsilon too small to achieve (warning only)"
    elseif ret==ERR_MAXNALLOC
        msg = "attemped to allocate internal array larger than MAX_NF (defined in defs.h)"
    elseif ret==ERR_SPREAD_BOX_SMALL
        msg = "spreader: fine grid too small compared to spread (kernel) width"
    elseif ret==ERR_SPREAD_PTS_OUT_RANGE
        msg = "spreader: if chkbnds=1, a nonuniform point is out of input range [-3pi,3pi]^d"
    elseif ret==ERR_SPREAD_ALLOC
        msg = "spreader: array allocation error"
    elseif ret==ERR_SPREAD_DIR
        msg = "spreader: illegal direction (should be 1 or 2)"
    elseif ret==ERR_UPSAMPFAC_TOO_SMALL
        msg = "upsampfac too small (should be >1.0)"
    elseif ret==HORNER_WRONG_BETA
        msg = "upsampfac not a value with known Horner eval: currently 2.0 or 1.25 only"
    elseif ret==ERR_NDATA_NOTVALID
        msg = "ntrans not valid in vectorized interface (should be >= 1)"
    elseif ret==ERR_TYPE_NOTVALID
        msg = "invalid transform type, type should be 1, 2, or 3"
    elseif ret==ERR_ALLOC
        msg = "general allocation failure"
    elseif ret==ERR_DIM_NOTVALID
        msg = "invalid dimension, should be 1, 2, or 3"
    elseif ret==ERR_SPREAD_THREAD_NOTVALID
        msg = "spread_thread option not valid"
    else
        msg = "error of type unknown to Julia interface! Check FINUFFT documentation"
    end
    throw(FINUFFTError(ret, msg))
end

# HELPER routines for guru interface...

### validate sizes of inputs for setpts
function valid_setpts(type::Integer,
                      dim::Integer,
                      x::Array{T},
                      y::Array{T}=T[],
                      z::Array{T}=T[],
                      s::Array{T}=T[],
                      t::Array{T}=T[],
                      u::Array{T}=T[]) where T <: finufftReal
    nj = length(x)
    if type==3
        nk = length(s)
    else
        nk = 0
    end

    if dim>1
        @assert nj == length(y)
        if type==3
            @assert nk == length(t)
        end
    end

    if dim>2
        @assert nj == length(z)
        if type==3
            @assert nk == length(u)
        end
    end

    return (nj, nk)
end

### validate number of transforms
function valid_ntr(x::Array{T},
                   c::Array{Complex{T}}) where T <: finufftReal
    ntrans = Cint(length(c) / length(x))
    @assert ntrans*length(x) == length(c)
    return ntrans
end

### infer number of modes from fk array
function get_nmodes_from_fk(dim::Integer,
                            fk::Array{Complex{T}}) where T <: finufftReal
    ndim = ndims(fk)
    @assert dim==1 || dim==2 || dim==3
    @assert ndim==dim || ndim==dim+1
    if ndim==dim
        ntrans = 1
        return (size(fk)...,ntrans)
    else
        return size(fk)
    end
end

### kwargs opt set
function setkwopts!(opts::nufft_opts; kwargs...)
    dtype = Float64
    for (key, value) in kwargs
        if hasproperty(opts, key::Symbol)
            setproperty!(opts, key, value)
        elseif String(key)=="dtype"
            @assert value <: finufftReal
            dtype = value
        else
            @warn "nufft_opts does not have attribute " * String(key)
        end
    end
    return dtype
end

### check kwargs with dtype
function checkkwdtype(dtype::DataType; kwargs...)
    for (key, value) in kwargs
        if String(key)=="dtype"
            @assert  value == dtype
        end
    end
end

# Finally, bring in the main user-facing interfaces...

include("guru.jl")
include("simple.jl")

end # module


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

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

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

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

    function OpDef_ArgDef(; kwargs...)
        obj = new(meta(OpDef_ArgDef), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct OpDef_ArgDef
const __meta_OpDef_ArgDef = Ref{ProtoMeta}()
function meta(::Type{OpDef_ArgDef})
    ProtoBuf.metalock() do
        if !isassigned(__meta_OpDef_ArgDef)
            __meta_OpDef_ArgDef[] = target = ProtoMeta(OpDef_ArgDef)
            fnum = Int[1,2,3,4,5,6,16]
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :description => AbstractString, :_type => Int32, :type_attr => AbstractString, :number_attr => AbstractString, :type_list_attr => AbstractString, :is_ref => Bool]
            meta(target, OpDef_ArgDef, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_OpDef_ArgDef[]
    end
end
function Base.getproperty(obj::OpDef_ArgDef, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :description
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :_type
        return (obj.__protobuf_jl_internal_values[name])::Int32
    elseif name === :type_attr
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :number_attr
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :type_list_attr
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :is_ref
        return (obj.__protobuf_jl_internal_values[name])::Bool
    else
        getfield(obj, name)
    end
end

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

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

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

    function OpDef(; kwargs...)
        obj = new(meta(OpDef), Dict{Symbol,Any}(), Set{Symbol}())
        values = obj.__protobuf_jl_internal_values
        symdict = obj.__protobuf_jl_internal_meta.symdict
        for nv in kwargs
            fldname, fldval = nv
            fldtype = symdict[fldname].jtyp
            (fldname in keys(symdict)) || error(string(typeof(obj), " has no field with name ", fldname))
            values[fldname] = isa(fldval, fldtype) ? fldval : convert(fldtype, fldval)
        end
        obj
    end
end # mutable struct OpDef
const __meta_OpDef = Ref{ProtoMeta}()
function meta(::Type{OpDef})
    ProtoBuf.metalock() do
        if !isassigned(__meta_OpDef)
            __meta_OpDef[] = target = ProtoMeta(OpDef)
            fnum = Int[1,2,3,20,4,8,5,6,18,16,17,19]
            allflds = Pair{Symbol,Union{Type,String}}[:name => AbstractString, :input_arg => Base.Vector{OpDef_ArgDef}, :output_arg => Base.Vector{OpDef_ArgDef}, :control_output => Base.Vector{AbstractString}, :attr => Base.Vector{OpDef_AttrDef}, :deprecation => OpDeprecation, :summary => AbstractString, :description => AbstractString, :is_commutative => Bool, :is_aggregate => Bool, :is_stateful => Bool, :allows_uninitialized_input => Bool]
            meta(target, OpDef, allflds, ProtoBuf.DEF_REQ, fnum, ProtoBuf.DEF_VAL, ProtoBuf.DEF_PACK, ProtoBuf.DEF_WTYPES, ProtoBuf.DEF_ONEOFS, ProtoBuf.DEF_ONEOF_NAMES)
        end
        __meta_OpDef[]
    end
end
function Base.getproperty(obj::OpDef, name::Symbol)
    if name === :name
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :input_arg
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{OpDef_ArgDef}
    elseif name === :output_arg
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{OpDef_ArgDef}
    elseif name === :control_output
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{AbstractString}
    elseif name === :attr
        return (obj.__protobuf_jl_internal_values[name])::Base.Vector{OpDef_AttrDef}
    elseif name === :deprecation
        return (obj.__protobuf_jl_internal_values[name])::OpDeprecation
    elseif name === :summary
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :description
        return (obj.__protobuf_jl_internal_values[name])::AbstractString
    elseif name === :is_commutative
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :is_aggregate
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :is_stateful
        return (obj.__protobuf_jl_internal_values[name])::Bool
    elseif name === :allows_uninitialized_input
        return (obj.__protobuf_jl_internal_values[name])::Bool
    else
        getfield(obj, name)
    end
end

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

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

export OpDef_ArgDef, OpDef_AttrDef, OpDef, OpDeprecation, OpList