text:
#This wraps all the basic ArrayFire routines as is.

#Math functions
import Base: +, -, *, /

# Resolve conflicts
+(x::AFAbstractArray{Bool},y::Bool) = AFArray{Bool}(@cxx +(x.array,y))
+(y::Bool,x::AFAbstractArray{Bool}) = AFArray{Bool}(@cxx +(y,x.array))
-(x::AFAbstractArray{Bool},y::Bool) = AFArray{Bool}(@cxx -(x.array,y))
-(y::Bool,x::AFAbstractArray{Bool}) = AFArray{Bool}(@cxx -(y,x.array))

for (op,cppop) in ((:+,:+),(:(.+),:+),(:-,:-),(:(.-),:-),(:.*,:*),(:./,:/),(:.>>,:>>),(:.<<,:<<))
    @eval function Base.($(quot(op))){T,S}(x::AFAbstractArray{T}, y::AFAbstractArray{S})
        a1 = x.array
        a2 = y.array
        AFArray{af_promote(T,S)}(@cxx ($(cppop))(a1,a2))
    end
    @eval function Base.($(quot(op))){T,S<:Number}(x::AFAbstractArray{T}, y::S)
        a = x.array
        # This is special behavior hardcoded in arrayfire for real floats
        ST = S <: AbstractFloat ? T : S
        AFArray{af_promote(T,ST)}(@cxx ($(cppop))(a, y))
    end
    @eval function Base.($(quot(op))){T,S<:Number}(y::S, x::AFAbstractArray{T})
        a = x.array
        # This is special behavior hardcoded in arrayfire for real floats
        ST = S <: AbstractFloat ? T : S
        AFArray{af_promote(T,ST)}(@cxx ($(cppop))(y, a))
    end
end

#Functions to generate arrays
af_constant(a::Real,b::vcpp"af::dim4", c::Cxx.CppEnum{:af_dtype}) = icxx"af::constant($a,$b,$c);"
af_randu(a::vcpp"af::dim4", b::Cxx.CppEnum{:af_dtype}) = icxx"af::randu($a,$b);"
af_randn(a::vcpp"af::dim4", b::Cxx.CppEnum{:af_dtype}) = icxx"af::randn($a,$b);"
af_identity(a::vcpp"af::dim4", b::Cxx.CppEnum{:af_dtype}) = icxx"af::identity($a,$b);"
af_diag(a::AFAbstractArray, b::Integer) = icxx"af::diag($a,$b);"
af_getSeed() = icxx"af::getSeed();"
af_range(a::vcpp"af::dim4", b::Integer, c::Cxx.CppEnum{:af_dtype}) = icxx"af::range($a,$b,$c);"
af_setSeed(a::Integer) = icxx"af::setSeed($a);"
af_iota(a::vcpp"af::dim4", b::vcpp"af::dim4", c::Cxx.CppEnum{:af_dtype}) = icxx"af::iota($a,$b,$c);"

#Function to modify or reorg arrays
af_tile(a::AFAbstractArray, b::vcpp"af::dim4") = icxx"af::tile($a,$b);"
af_join(a::Integer, b::AFAbstractArray, c::AFAbstractArray) = icxx"af::join($a,$b,$c);"
af_join(a::Integer, b::AFAbstractArray, c::AFAbstractArray, d::AFAbstractArray) = icxx"af::join($a,$b,$c,$d);"
af_join(a::Integer, b::AFAbstractArray, c::AFAbstractArray, d::AFAbstractArray,
            el::AFAbstractArray) = icxx"af::join($a,$b,$c,$d,$el);"
af_moddims(a::AFAbstractArray, b::vcpp"af::dim4") = icxx"af::moddims($a,$b);"
af_reorder(a::AFAbstractArray, b::Integer, c::Integer, d::Integer, el::Integer) = icxx"af::reorder($a,$b,$c,$d,$el);"
af_replace(a::AFAbstractArray, b::AFAbstractArray{Bool}, c::Union{Real,AFAbstractArray}) = icxx"af::replace($a,$b,$c);"
af_select(a,::AFAbstractArray{Bool}, b::AFAbstractArray, c::Union{Real,AFAbstractArray}) = icxx"af::select($a,$b,$c);"
af_shift(a::AFAbstractArray, b::Integer, c::Integer, d::Integer, el::Integer) = icxx"af::shift($a,$b,$c,$d,$el);"

#Some more array methods
af_complex_conj(a::AFAbstractArray) = icxx"$a.H();"
af_as(a::AFAbstractArray, to::Cxx.CppEnum{:af_dtype}) = icxx"(af::array)($a.as($to));"
af_copy(a::AFAbstractArray) = icxx"$a.copy();"
af_bytes(a::AFAbstractArray) = icxx"$a.bytes();"
af_isempty(a::AFAbstractArray) = icxx"$a.isempty();"
af_isscalar(a::AFAbstractArray) = icxx"$a.isscalar();"
af_isvector(a::AFAbstractArray) = icxx"$a.isvector();"
af_isrow(a::AFAbstractArray) = icxx"$a.isrow();"
af_iscolumn(a::AFAbstractArray) = icxx"$a.iscolumn();"
af_iscomplex(a::AFAbstractArray) = icxx"$a.iscomplex();"
af_isreal(a::AFAbstractArray) = icxx"$a.isreal();"
af_isdouble(a::AFAbstractArray) = icxx"$a.isdouble();"
af_issingle(a::AFAbstractArray) = icxx"$a.issingle();"
af_isrealfloating(a::AFAbstractArray) = icxx"$a.isrealfloating();"
af_isfloating(a::AFAbstractArray) = icxx"$a.isfloating();"
af_isinteger(a::AFAbstractArray) = icxx"$a.isinteger();"
af_isbool(a::AFAbstractArray) = icxx"$a.isbool();"

#Helper functions
af_isNaN(a::AFAbstractArray) = icxx"af::isNaN($a);"
af_isinf(a::AFAbstractArray) = icxx"af::isInf($a);"
af_iszero(a::AFAbstractArray) = icxx"af::iszero($a);"

#ArrayFire info functions
af_info() = icxx"af::info();"
af_isDoubleAvailable() = icxx"af::isDoubleAvailable();"
af_sync() = icxx"af::sync();"
af_sync(a) = icxx"af::sync($a);"
af_getDevice() = icxx"af::getDevice();"
af_setDevice(a) = icxx"af::setDevice($a);"

#Loading and Saving Arrays
af_saveArray(a::AFAbstractArray, b::AbstractString, c::AbstractString) = icxx"af::saveArray($a,$b,$c);"
af_readArray(a::AbstractString, b::AbstractString) = icxx"af::readArray($a,$b);"

#Numeric functions
af_abs(a::AFAbstractArray) = icxx"af::abs($a);"
af_max(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = icxx"af::max($a, $b);"
af_min(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = icxx"af::min($a, $b);"
af_arg(a::AFAbstractArray) = icxx"af::arg($a);"
af_ceil(a::AFAbstractArray) = icxx"af::cell($a);"
af_floor(a::AFAbstractArray) = icxx"af::floor($a);"
af_hypot(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::hypot($a, $b);"
af_mod(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::mod($a, $b);"
of_rem(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::rem($a, $b);"
af_round(a::AFAbstractArray) = icxx"af::round($a);"
af_sign(a::AFAbstractArray) = icxx"af::sign($a);"
af_trunc(a::AFAbstractArray) = icxx"af_trunc($a);"

#Logical ops
af_equals(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = icxx"$a == $b;"
af_gt(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = icxx"$a > $b;"
af_ge(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = icxx"$a >= $b;"
af_bitand(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = icxx"$a & $b;"
af_and(a::AFAbstractArray, b::AFAbstractArray) = icxx"$a && $b;"
af_bitor(a::AFAbstractArray, b::AFAbstractArray) = icxx"$a | $b;"
af_bitxor(a::AFAbstractArray, b::AFAbstractArray) = icxx"$a ^ $b;"
af_le(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = icxx"$a <= $b;"
af_lt(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = icxx"$a < $b;"
af_neg(a::AFAbstractArray) = icxx"-$a;"
af_neq(a::AFAbstractArray, b::AFAbstractArray) = icxx"$a != $b;"
af_not(a::AFAbstractArray) = icxx"!$a;"
af_or(a::AFAbstractArray, b::AFAbstractArray) = icxx"$a || $b;"

#Complex functions
af_complex(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::complex($a,$b);"
af_complex(a::AFAbstractArray) = icxx"af::complex($a);"
af_conjg(a::AFAbstractArray) = icxx"af::conjg($a);"
af_imag(a::AFAbstractArray) = icxx"af::imag($a);"
af_real(a::AFAbstractArray) = icxx"af::real($a);"

#Expontential functions
af_pow(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = icxx"af::pow($a, $b);"
af_root(a::Union{Real,AFAbstractArray}, b::Union{Real,AFAbstractArray}) = icxx"af::root($a, $b);"

#Linear algebra
const AF_MAT_NONE = icxx"AF_MAT_NONE;"
const AF_MAT_CTRANS = icxx"AF_MAT_CTRANS;"

af_dot(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::dot($a, $b);"

af_matmul(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::matmul($a, $b);"
af_matmul3(a::AFAbstractArray, b::AFAbstractArray, c::AFAbstractArray) = icxx"af::matmul($a, $b, $c);"
af_matmul4(a::AFAbstractArray, b::AFAbstractArray, c::AFAbstractArray, d::AFAbstractArray) = icxx"af::matmul($a, $b, $c, $d);"
af_matmul_flags(a::AFAbstractArray, b::AFAbstractArray, flag1::Cxx.CppEnum{:af_mat_prop}, flag2::Cxx.CppEnum{:af_mat_prop}) = icxx"af::matmul($a,$b,$flag1,$flag2);"
af_matmulNT(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::matmulNT($a,$b);"
af_matmulTN(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::matmulTN($a,$b);"
af_matmulTT(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::matmulTT($a,$b);"

af_transpose(a::AFAbstractArray) = icxx"af::transpose($a);"
af_ctranspose(a::AFAbstractArray) = icxx"af::transpose($a, true);"
af_transposeInPlace(a::AFAbstractArray) = icxx"af::transposeInPlace($a);"
af_ctransposeInPlace(a::AFAbstractArray) = icxx"af:transposeInPlace($a, true);"

af_solve(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::solve($a,$b);"
af_solveLU(a::AFAbstractArray, b::AFAbstractArray) = icxx"af_solveLU($a,$b);"

af_cholesky(out::vcpp"af::array", a::AFMatrix, flag::Bool) = icxx"af::cholesky($out,$a,$flag);"
af_choleskyInPlace(a::AFMatrix, flag::Bool) = icxx"af::choleskyInPlace($a,$flag);"
af_lu(l::vcpp"af::array", u::vcpp"af::array", p::vcpp"af::array", a::AFMatrix) = icxx"af::lu($l,$u,$p,$a);"
af_qr(q::vcpp"af::array", r::vcpp"af::array", tau::vcpp"af::array", a::AFMatrix) = icxx"af::qr($q,$r,$tau,$a);"
af_svd(u::vcpp"af::array", s::vcpp"af::array", vt::vcpp"af::array", a::AFMatrix) = icxx"af::svd($u,$s, $vt,$a);"
af_upper(a::Union{AFMatrix, vcpp"af::array"}) = icxx"af::upper($a);"
af_lower(a::Union{AFMatrix, vcpp"af::array"}) = icxx"af::lower($a);"

af_det(a::AFAbstractArray) = icxx"af::det<float>($a);"
af_inverse(a::AFAbstractArray) = icxx"af::inverse($a);"
af_norm(a::AFAbstractArray) = icxx"af::norm($a);"
af_rank(a::AFAbstractArray) = icxx"af::rank($a);"

#Signal
af_fft(a::AFVector) = icxx"af::fft($a);"
af_fft2(a::AFMatrix) = icxx"af::fft2($a);"
af_fft3(a::AFAbstractArray) = icxx"af::fft3($a);"
af_ifft(a::AFVector) = icxx"af::ifft($a);"
af_ifft2(a::AFMatrix) = icxx"af::ifft2($a);"
af_ifft3(a::AFAbstractArray) = icxx"af::ifft3($a);"

af_fir(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::fir($a, $b);"
af_iir(a::AFAbstractArray, b::AFAbstractArray, c::AFAbstractArray) = icxx"af::iir($a, $b, $c);"

const AF_CONV_EXPAND = icxx"AF_CONV_EXPAND;"
const AF_CONV_DEFAULT = icxx"AF_CONV_DEFAULT;"

af_convolve(a::AFAbstractArray, b::AFAbstractArray, mode=AF_CONV_DEFAULT) = icxx"af::convolve($a,$b,$mode);"
af_convolve1(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::convolve1($a,$b);"
af_convolve2(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::convolve2($a,$b);"
af_convolve3(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::convolve3($a,$b);"
af_fftConvolve(a::AFAbstractArray, b::AFAbstractArray, flag::Cxx.CppEnum{:af_conv_mode}) = icxx"af::fftConvolve($a,$b,$flag);"
af_fftConvolve2(a::AFAbstractArray, b::AFAbstractArray, flag::Cxx.CppEnum{:af_conv_mode}) = icxx"af::fftConvolve2($a,$b,$flag);"
af_fftConvolve3(a::AFAbstractArray, b::AFAbstractArray, flag::Cxx.CppEnum{:af_conv_mode}) = icxx"af::fftConvolve3($a,$b,$flag);"

af_approx1(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::approx1($a,$b);"
af_approx2(a::AFAbstractArray, b::AFAbstractArray, c::AFAbstractArray) = icxx"af::approx2($a,$b,$c);"

#Statistics
af_mean(a::AFAbstractArray) = icxx"af::mean<float>($a);"
af_mean(a::AFAbstractArray, b::Integer) = icxx"af::mean($a, $b);"
af_mean(a::AFAbstractArray, b::AFAbstractArray, c::Integer) = icxx"af::mean($a, $b, $c);"
af_median(a::AFAbstractArray) = icxx"af::median<float>($a);"
af_median(a::AFAbstractArray, b::Integer) = icxx"af::median($a,$b);"
af_stdev(a::AFAbstractArray) = icxx"af::stdev<float>($a);"
af_stdev(a::AFAbstractArray, b::Integer) = icxx"af::stdev($a,$b);"
af_var(a::AFAbstractArray) = icxx"af::var<float>($a);"
af_var(a::AFAbstractArray, b::Integer) = icxx"af::var($a,$b);"
af_var(a::AFAbstractArray, b::AFAbstractArray, c::Integer) = icxx"af::var($a,$b,$c);"
af_cov(a::AFAbstractArray, b::Integer) = icxx"af::cov($a,$b);"
af_corrcoef(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::corrcoef($a, $b);"

#Vector algorithms
af_sum(a::AFAbstractArray) = icxx"af::sum<float>($a);"
af_sum(a::AFAbstractArray, b::Integer) = icxx"af::sum($a,$b);"
af_product(a::AFAbstractArray) = icxx"af::product<float>($a);"
af_product(a::AFAbstractArray, b::Integer) = icxx"af::product($a,$b);"
af_max(a::AFAbstractArray) = icxx"af::max<float>($a);"
af_max(a::AFAbstractArray, b::Real) = icxx"af::max($a,$b);"
af_max_dim(a::AFAbstractArray, b::Integer) = icxx"""
    int dim = $b;
    af::max($a,dim);
    """
af_min(a::AFAbstractArray) = icxx"af::min<float>($a);"
af_min(a::AFAbstractArray, b::Real) = icxx"af::min($a,$b);"
af_min_dim(a::AFAbstractArray, b::Integer) = icxx"""
    int dim = $b;
    af::min($a,dim);
"""
af_anyTrue(a::AFAbstractArray) = icxx"af::anyTrue<bool>($a);"
af_anyTrue(a::AFAbstractArray, b::Integer) = icxx"af::anyTrue($a,$b);"
af_allTrue(a::AFAbstractArray) = icxx"af::allTrue<bool>($a);"
af_allTrue(a::AFAbstractArray, b::Integer) = icxx"af::allTrue($a,$b);"
af_count(a::AFAbstractArray) = icxx"af::count<int>($a);"
af_count(a::AFAbstractArray, b::Integer) = icxx"af::count($a,$b);"
af_accum(a::AFAbstractArray, b::Integer) = icxx"af::accum($a,$b);"
af_where(a::AFAbstractArray) = icxx"af::where($a);"
af_diff(a::AFAbstractArray, b::Integer) = icxx"af::diff($a,$b);"
af_diff2(a::AFAbstractArray, b::Integer) = icxx"af::diff2($a,$b);"
af_grad(a::AFAbstractArray, b::AFAbstractArray, c::AFAbstractArray) = icxx"af::grad($a,$b,$c);"
af_sort(a::AFAbstractArray, b::Integer, flag::Bool) = icxx"af::sort($a,$b,$flag);"
af_sort(a::AFAbstractArray, b::AFAbstractArray, c::AFAbstractArray, d::Integer,flag::Bool) = icxx"af::sort($a,$b,$c,$d,$flag);"
af_sort(a::AFAbstractArray, b::AFAbstractArray, c::AFAbstractArray, d::AFAbstractArray, el::Integer,flag::Bool) = icxx"af::sort($a,$b,$c,$d,$flag);"
af_setIntersect(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::setIntersect($a,$b);"
af_setUnion(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::setUnion($a,$b);"
af_setUnique(a::AFAbstractArray) = icxx"af::setUnique($a);"
af_flip(a::AFAbstractArray, b::Integer) = icxx"af::flip($a,$b);"
af_flat(a::AFAbstractArray) = icxx"af::flat($a);"

#Image
af_loadImage(a::AbstractString, b::Bool) = icxx"af::loadImage($a,$b);"
af_saveImage(a::AbstractString, b::AFAbstractArray) = icxx"af::saveImage($a,$b);"
af_rotate(a::AFAbstractArray, b::Real) = icxx"af::rotate($a,$b);"
af_scale(a::AFAbstractArray, b::Real, c::Real) = icxx"af::scale($a,$b,$c);"
af_skew(a::AFAbstractArray, b::Real, c::Real) = icxx"af::skew($a,$b,$c);"
af_translate(a::AFAbstractArray, b::Real, c::Real) = icxx"af::translate($a,$b,$c);"
af_transform(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::transform($a,$b);"

af_regions(a::AFAbstractArray) = icxx"$a.as(b8);"

const AF_RGB = icxx"AF_RGB;"
const AF_GRAY= icxx"AF_GRAY;"
const AF_HSV = icxx"AF_HSV;"

af_gray2rgb(a::AFAbstractArray) = icxx"af::gray2rgb($a);"
af_hsv2rgb(a::AFAbstractArray) = icxx"af::hsv2rgb($a);"
af_rgb2gray(a::AFAbstractArray) = icxx"af::rgb2gray($a);"
af_rgb2hsv(a::AFAbstractArray) = icxx"af::rgb2hsv($a);"
af_rgb2ycbcr(a::AFAbstractArray) = icxx"af::rgb2ycbcr($a);"
af_ycbcr2rgb(a::AFAbstractArray) = icxx"af::ycbcr2rgb($a);"
af_colorspace(a::AFAbstractArray, b::Cxx.CppEnum{:af_cspace_t},c::Cxx.CppEnum{:af_cspace_t}) = icxx"af::colorspace($a,$b,$c);"

af_SAT(a::AFAbstractArray) = icxx"af::sat($a);"
af_bilateral(a::AFAbstractArray, b::Real, c::Real) = icxx"af::bilateral($a,$b,$c);"
af_maxfilt(a::AFAbstractArray) = icxx"af::maxfilt($a);"
af_medfilt(a::AFAbstractArray) = icxx"af::medfilt($a);"
af_minfilt(a::AFAbstractArray) = icxx"af::minfilt($a);"
af_sobel(a::AFAbstractArray) = icxx"af::sobel($a);"
af_meanShift(a::AFAbstractArray, b::Real, c::Real, d::Integer) = icxx"af::meanShift($a,$b,$c,$d);"

af_histogram(a::AFAbstractArray, b::Integer) = icxx"af::histogram($a,$b);"
af_histogram(a::AFAbstractArray, b::Integer, c::Real, d::Real) = icxx"af::histogram($a,$b,$c,$d);"
af_histEqual(a::AFAbstractArray, b::AFAbstractArray) = icxx"af::histEqual($a, $b);"

af_dilate(a::AFAbstractArray,b::AFAbstractArray) = icxx"af::dilate($a,$b);"
af_dilate3(a::AFAbstractArray,b::AFAbstractArray) = icxx"af::dilate3($a,$b);"
af_erode(a::AFAbstractArray,b::AFAbstractArray) = icxx"af::erode($a,$b);"
af_erode3(a::AFAbstractArray,b::AFAbstractArray) = icxx"af::erode3($a,$b);"

#Computer Vision
const AF_SAD = icxx"AF_SAD;"
const AF_ZSAD = icxx"AF_ZSAD;"
const AF_LSAD = icxx"AF_LSAD;"
const AF_SSD = icxx"AF_SSD;"
const AF_ZSSD = icxx"AF_ZSSD;"
const AF_LSSD = icxx"AF_LSSD;"
const AF_NCC = icxx"AF_NCC;"
const AF_ZNCC = icxx"AF_ZNCC;"
const AF_SHD = icxx"AF_SHD;"

af_getX(a::AFFeatures) = icxx"$a.getX();"
af_getY(a::AFFeatures) = icxx"$a.getY();"
af_getSize(a::AFFeatures) = icxx"$a.getSize();"
af_getScore(a::AFFeatures) = icxx"$a.getScore();"
af_getOrientation(a::AFFeatures) = icxx"$a.getOrientation();"

af_gaussianKernel(a::Integer, b::Integer) = icxx"af::gaussianKernel($a,$b);"
af_dog(a::AFAbstractArray, b::Integer, c::Integer) = icxx"af::dog($a,$b,$c);"
af_matchTemplate(a::AFAbstractArray, b::AFAbstractArray, c::Cxx.CppEnum{:af_match_type}) = icxx"af::matchTemplate($a,$b,$c);"
af_orb(a::AFAbstractArray, b::AFAbstractArray, c::AFAbstractArray, d::Real,
        el::Integer, f::Real, g::Integer, h::Bool) = icxx"af::orb($a,$b,$c,$d,$el,$f,$g,$h);"
af_sift(a::AFAbstractArray, b::AFAbstractArray, c::AFAbstractArray, d::Integer,
        el::Real, f::Real, g::Real, h::Bool, i::Real,j::Real) = icxx"af::sift($a,$b,$c,$d,$el,$f,$g,$h,$i,$j);"
af_fast(a::AFAbstractArray, b::Real, c::Integer, d::Bool, el::Real, f::Integer) = icxx"af::fast($a,$b,$c,$d,$el,$f);"
af_harris(a::AFAbstractArray, b::Integer, c::Real, d::Real, el::Integer, f::Real) = icxx"af::harris($a,$b,$c,$d,$el,$f);"
af_susan(a::AFAbstractArray, b::Integer, c::Real, d::Real, el::Real, f::Integer) = icxx"af::susan($a,$b,$c,$d,$el,$f);"
af_hammingMatcher(a::AFAbstractArray, b::AFAbstractArray, c::AFAbstractArray,
                    d::AFAbstractArray, el::Integer, f::Integer) = icxx"af::hammingMatcher($a,$b,$c,$d,$el,$f);"
af_nearestNeighbour(a::AFAbstractArray, b::AFAbstractArray, c::AFAbstractArray,
                        d::AFAbstractArray, el::Integer, f::Integer,
                        g::Cxx.CppEnum{:af_match_type}) = icxx"af::nearestNeighbour($a,$b,$c,$d,$el,$f,$g);"

#Graphics
af_window() = icxx"af::Window();"
af_window(a::AbstractString) = icxx"af::Window($a);"
af_window(a::Integer, b::Integer, c::AbstractString) = icxx"af::Window($a,$b,$c);"
af_setTitle(a::AFWindow, b::AbstractString) = icxx"$a.setTitle($b);"
af_setImage(a::AFWindow, b::AFAbstractArray, c::AbstractString) = icxx"$a.image($b,$c);"
af_setPlot(a::AFWindow, b::AFAbstractArray, c::AFAbstractArray, d::AbstractString) = icxx"$a.plot($b,$c,$d);"
af_show(a::AFWindow) = icxx"$a.show();"
af_setSurface(a::AFWindow, b::AFAbstractArray) = icxx"$a.surface($b);"

#Unified Backend
const AF_BACKEND_CPU = icxx"AF_BACKEND_CPU;"
const AF_BACKEND_CUDA = icxx"AF_BACKEND_CUDA;"
const AF_BACKEND_OPENCL = icxx"AF_BACKEND_OPENCL;"

af_setBackend(backend::Cxx.CppEnum{:af_backend}) = icxx"af::setBackend($backend);"
function af_getActiveBackend()
    icxx"""
        af_backend b;
        af_get_active_backend(&b);
        return b;
        """
end
af_getAvailableBackends() = icxx"af::getAvailableBackends();"
af_getBackendCount() = icxx"af::getBackendCount();"
af_getBackendId(a::AFAbstractArray) = icxx"af::getBackendId($a);"
af_getDeviceId(a::AFAbstractArray)  = icxx"af::getDeviceId($a);"


text:
using DrWatson
@quickactivate "INN_Velocity-Migration"

using JLD, Statistics, PyPlot, ImageFiltering, InvertibleNetworks, LinearAlgebra, Dates
include(srcdir("FWI_64x200.jl"))
using .FWI_64x200

X_orig = load(datadir("vel_4k_samples_64x200_lin_vel.jld"), "m_all")
X_orig = Float32.(X_orig)
n1, n2, nc, nsamples = size(X_orig)
ntrain = Int(nsamples*.9)
X_train_orig=X_orig[:, :, :, 1:ntrain]

AN = ActNorm(ntrain)
X_train_orig = AN.forward(X_train_orig) # zero mean and unit std of the training data X

# define data and figure path
datapath = plotsdir("chint/test3_6000_16_16")
figfolder = "fwi/test3_6000_16_16/posterior_samples2"
mkpath(plotsdir(figfolder))

X_fixed = load(joinpath(datapath, "posterior_samples2.jld"), "X_fixed")
X_post  = load(joinpath(datapath, "posterior_samples2.jld"), "X_post")
X_fixed = AN.inverse(X_fixed) # unnormarlize based on the same parameter of the training data X
X_post = AN.inverse(X_post)   # unnormarlize based on the same parameter of the training data X
X_post_mean = mean(X_post; dims=4)

m = X_fixed[:, :, 1, 1] # true velocity in squared slowness
idx_wb = 9 # index of the water bottom
niterations = 15

####################################################################################################
## Scenario 1: initial velocity m0 = all posterior samples
println(string("Scenario 1: initial velocity m0 = posterior samples"))

test_size = size(X_post)[4] #100
m0_all = zeros(Float32, n1, n2, test_size)
fwi_all = zeros(Float32, n1, n2, test_size)
obj_all = zeros(Float32, niterations, test_size)

for i = 1:3 #test_size

    m0 = X_post[:, :, 1, i]
    m0[1:idx_wb, :] .= m[1,1]
    m0[idx_wb+1:end, :] = imfilter(m0[idx_wb+1:end, :], Kernel.gaussian(3)) # smoothed velocity
    fwi_all[:,:,i], obj_all[:,i] = fwi(m, m0, niterations)

    m0_all[:,:,i] = m0

    println(string("FWI using the ", i, "th sample of X posterior finished"))
end

save(plotsdir(figfolder, "fwi_post_samples.jld"), "m0_all", m0_all, "fwi_all", fwi_all, "obj_all", obj_all)

i=1
v = sqrt.(1f0 ./ m)
v0 = sqrt.(1f0 ./ m0_all[:,:,i])
fwi_v = sqrt.(1f0 ./ fwi_all[:,:,i])
obj = obj_all[:,i]

figure(figsize=[16,8])
ax1 = subplot(2, 2, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(2, 2, 2); imshow(v0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("Initial model, L2=", round(norm(v-v0,2), sigdigits=4)))
ax3 = subplot(2, 2, 3); imshow(fwi_v, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi_v,2), sigdigits=4)))
ax4 = subplot(2, 2, 4); plot(obj); title("Objective function value")
savefig(plotsdir(figfolder, "fwi_post_samples.png"))

####################################################################################################
## Scenario 2: initial velocity m0 = posterior mean
println(string("Scenario 2: initial velocity m0 = posterior mean"))

m0 = X_post_mean[:, :, 1, 1]
m0[1:idx_wb, :] .= m[1,1]
m0[idx_wb+1:end, :] = imfilter(m0[idx_wb+1:end, :], Kernel.gaussian(3)) # smoothed velocity

fwi2, obj2 = fwi(m, m0, niterations)
save(plotsdir(figfolder, "fwi_post_mean.jld"), "m", m, "m0", m0, "fwi", fwi2, "obj", obj2)

v = sqrt.(1f0 ./ m)
v0 = sqrt.(1f0 ./ m0)
fwi_v = sqrt.(1f0 ./ fwi2)
obj = obj2

figure(figsize=[16,8])
ax1 = subplot(2, 2, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(2, 2, 2); imshow(v0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("Initial model, L2=", round(norm(v-v0,2), sigdigits=4)))
ax3 = subplot(2, 2, 3); imshow(fwi_v, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi_v,2), sigdigits=4)))
ax4 = subplot(2, 2, 4); plot(obj); title("Objective function value")
savefig(plotsdir(figfolder, "fwi_post_mean.png"))

####################################################################################################
## Scenario 3: initial velocity m0 = migration velocity used in rtm
println(string("Scenario 3: initial velocity m0 = migration velocity used in rtm"))

# build a migration velocity m0 whose value increases linearly with depth
v0 = ones(n1, 1)
v0[idx_wb+1:end] = collect(range(3.5, 5.5, length = n1-idx_wb))
v0 = (1f0 ./ v0).^2 # convert to slowness squared
v0[1:idx_wb] .= m[1,1] # water layer velocity
m_lin = repeat(v0, 1, n2)

fwi3, obj3 = fwi(m, m_lin, niterations)
save(plotsdir(figfolder, "fwi_lin_vel.jld"), "m", m, "m0", m_lin, "fwi", fwi3, "obj", obj3)

v = sqrt.(1f0 ./ m)
v0 = sqrt.(1f0 ./ m_lin)
fwi_v = sqrt.(1f0 ./ fwi3)
obj = obj3

figure(figsize=[16,8])
ax1 = subplot(2, 2, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(2, 2, 2); imshow(v0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("Initial model, L2=", round(norm(v-v0,2), sigdigits=4)))
ax3 = subplot(2, 2, 3); imshow(fwi_v, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi_v,2), sigdigits=4)))
ax4 = subplot(2, 2, 4); plot(obj); title("Objective function value")
savefig(plotsdir(figfolder, "fwi_lin_vel.png"))

####################################################################################################
## Scenario 4: initial velocity m0 = smoothed true velocity
println(string("Scenario 4: initial velocity m0 = smoothed true velocity"))

m_smooth = m[1,1] .* ones(Float32, size(m)) # smoothed velocity
m_smooth[idx_wb+1:end, :] = imfilter(m[idx_wb+1:end, :], Kernel.gaussian(3)) # smoothed velocity

fwi4, obj4 = fwi(m, m_smooth, niterations)
save(plotsdir(figfolder, "fwi_smooth_vel.jld"), "m", m, "m0", m_smooth, "fwi", fwi4, "obj", obj4)

v = sqrt.(1f0 ./ m)
v0 = sqrt.(1f0 ./ m_smooth)
fwi_v = sqrt.(1f0 ./ fwi4)
obj = obj4

figure(figsize=[16,8])
ax1 = subplot(2, 2, 1); imshow(v, cmap="jet", aspect = "auto"); plt.colorbar(); title("True velocity (km/s)")
ax2 = subplot(2, 2, 2); imshow(v0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("Initial model, L2=", round(norm(v-v0,2), sigdigits=4)))
ax3 = subplot(2, 2, 3); imshow(fwi_v, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi_v,2), sigdigits=4)))
ax4 = subplot(2, 2, 4); plot(obj); title("Objective function value")
savefig(plotsdir(figfolder, "fwi_smooth_vel.png"))

####################################################################################################
# Load data for Scenario 1-4
m0_all, fwi_all, obj_all = load(string("./figs/", figfolder, "/fwi_post_samples.jld"), "m0_all", "fwi_all", "obj_all")
m, m0, fwi2, obj2 = load(string("./figs/", figfolder, "/fwi_post_mean.jld"), "m", "m0", "fwi", "obj")
m, m_lin, fwi3, obj3 = load(string("./figs/", figfolder, "/fwi_lin_vel.jld"), "m", "m0", "fwi", "obj")
m, m_smooth, fwi4, obj4 = load(string("./figs/", figfolder, "/fwi_smooth_vel.jld"), "m", "m0", "fwi", "obj")

# Figures in slowness squared (s^2/km^2)
figure(); imshow(m, cmap="jet", aspect = "auto"); plt.colorbar(); title(L"True velocity ($s^2/km^2$)")
savefig(plotsdir(figfolder, "vel_m.png"))

figure(figsize=[20,12])
ax1 = subplot(3, 4, 1); imshow(m0_all[idx_wb+1:end,:,1], cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed posterior sample1")
ax2 = subplot(3, 4, 2); imshow(m0_all[idx_wb+1:end,:,2], cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed posterior sample2")
ax3 = subplot(3, 4, 3); imshow(m0_all[idx_wb+1:end,:,3], cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed posterior sample3")
ax4 = subplot(3, 4, 4); imshow(m0[idx_wb+1:end, :], cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed posterior mean")
ax5 = subplot(3, 4, 5); imshow(fwi_all[idx_wb+1:end,:,1], cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(m-fwi_all[:,:,1],2), sigdigits=4)))
ax6 = subplot(3, 4, 6); imshow(fwi_all[idx_wb+1:end,:,2], cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(m-fwi_all[:,:,2],2), sigdigits=4)))
ax7 = subplot(3, 4, 7); imshow(fwi_all[idx_wb+1:end,:,3], cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(m-fwi_all[:,:,3],2), sigdigits=4)))
ax8 = subplot(3, 4, 8); imshow(fwi2[idx_wb+1:end, :], cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(m-fwi2,2), sigdigits=4)))
ax9 = subplot(3, 4, 9); plot(obj_all[:,1]); title("Objective function value")
ax10 = subplot(3, 4, 10); plot(obj_all[:,2]); title("Objective function value")
ax11 = subplot(3, 4, 11); plot(obj_all[:,3]); title("Objective function value")
ax12 = subplot(3, 4, 12); plot(obj2); title("Objective function value")
savefig(plotsdir(figfolder, "fwi_m_post.png"))

figure(figsize=[20,12])
ax1 = subplot(3, 4, 1); imshow(m0_all[idx_wb+1:end,:,1], cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed posterior sample")
ax2 = subplot(3, 4, 2); imshow(m0[idx_wb+1:end, :], cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed posterior mean")
ax3 = subplot(3, 4, 3); imshow(m_lin[idx_wb+1:end, :], cmap="jet", aspect = "auto"); plt.colorbar(); title("Migration velocity in RTM")
ax4 = subplot(3, 4, 4); imshow(m_smooth[idx_wb+1:end, :], cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed true velocity")
ax5 = subplot(3, 4, 5); imshow(fwi_all[idx_wb+1:end,:,1], cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(m-fwi_all[:,:,1],2), sigdigits=4)))
ax6 = subplot(3, 4, 6); imshow(fwi2[idx_wb+1:end, :], cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(m-fwi2,2), sigdigits=4)))
ax7 = subplot(3, 4, 7); imshow(fwi3[idx_wb+1:end, :], cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(m-fwi3,2), sigdigits=4)))
ax8 = subplot(3, 4, 8); imshow(fwi4[idx_wb+1:end, :], cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(m-fwi4,2), sigdigits=4)))
ax9 = subplot(3, 4, 9); plot(obj_all[:,1]); title("Objective function value")
ax10 = subplot(3, 4, 10); plot(obj2); title("Objective function value")
ax11 = subplot(3, 4, 11); plot(obj3); title("Objective function value")
ax12 = subplot(3, 4, 12); plot(obj4); title("Objective function value")
savefig(plotsdir(figfolder, "fwi_m.png"))

# Figures in velocity (km/s)
v = sqrt.(1f0 ./ m)

v0_all = sqrt.(1f0 ./ m0_all[:,:,1:3])
v0 = sqrt.(1f0 ./ m0)
v_lin = sqrt.(1f0 ./ m_lin)
v_smooth = sqrt.(1f0 ./ m_smooth)

fwi_v_all = sqrt.(1f0 ./ fwi_all[:,:,1:3])
fwi2_v = sqrt.(1f0 ./ fwi2)
fwi3_v = sqrt.(1f0 ./ fwi3)
fwi4_v = sqrt.(1f0 ./ fwi4)

figure(); imshow(v, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title(L"True velocity ($km/s$)")
savefig(plotsdir(figfolder, "vel.png"))

figure(figsize=[20,12])
ax1 = subplot(3, 4, 1); imshow(v0_all[:,:,1], vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed posterior sample1")
ax2 = subplot(3, 4, 2); imshow(v0_all[:,:,2], vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed posterior sample2")
ax3 = subplot(3, 4, 3); imshow(v0_all[:,:,3], vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed posterior sample3")
ax4 = subplot(3, 4, 4); imshow(v0, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed posterior mean")
ax5 = subplot(3, 4, 5); imshow(fwi_v_all[:,:,1], vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi_v_all[:,:,1],2), sigdigits=4)))
ax6 = subplot(3, 4, 6); imshow(fwi_v_all[:,:,2], vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi_v_all[:,:,2],2), sigdigits=4)))
ax7 = subplot(3, 4, 7); imshow(fwi_v_all[:,:,3], vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi_v_all[:,:,3],2), sigdigits=4)))
ax8 = subplot(3, 4, 8); imshow(fwi2_v, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi2_v,2), sigdigits=4)))
ax9 = subplot(3, 4, 9); plot(obj_all[:,1]); title("Objective function value")
ax10 = subplot(3, 4, 10); plot(obj_all[:,2]); title("Objective function value")
ax11 = subplot(3, 4, 11); plot(obj_all[:,3]); title("Objective function value")
ax12 = subplot(3, 4, 12); plot(obj2); title("Objective function value")
savefig(plotsdir(figfolder, "fwi_v_post.png"))

figure(figsize=[20,12])
ax1 = subplot(3, 4, 1); imshow(v0_all[:,:,1], vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed posterior sample")
ax2 = subplot(3, 4, 2); imshow(v0, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed posterior mean")
ax3 = subplot(3, 4, 3); imshow(v_lin, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title("Migration velocity in RTM")
ax4 = subplot(3, 4, 4); imshow(v_smooth, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title("Smoothed true velocity")
ax5 = subplot(3, 4, 5); imshow(fwi_v_all[:,:,1], vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi_v_all[:,:,1],2), sigdigits=4)))
ax6 = subplot(3, 4, 6); imshow(fwi2_v, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi2_v,2), sigdigits=4)))
ax7 = subplot(3, 4, 7); imshow(fwi3_v, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi3_v,2), sigdigits=4)))
ax8 = subplot(3, 4, 8); imshow(fwi4_v, vmin=1.8, vmax=5.0, cmap="jet", aspect = "auto"); plt.colorbar(); title(string("FWI result, L2=", round(norm(v-fwi4_v,2), sigdigits=4)))
ax9 = subplot(3, 4, 9); plot(obj_all[:,1]); title("Objective function value")
ax10 = subplot(3, 4, 10); plot(obj2); title("Objective function value")
ax11 = subplot(3, 4, 11); plot(obj3); title("Objective function value")
ax12 = subplot(3, 4, 12); plot(obj4); title("Objective function value")
savefig(plotsdir(figfolder, "fwi_v.png"))

ext: jl


text:
############################################################################################
"""
$(TYPEDEF)
Contains information for number of discretization steps in MCMC algorithm.

# Fields
$(TYPEDFIELDS)
"""
mutable struct StepNumberTune{T<:BaytesCore.UpdateBool}
    "If true, number of steps will be adapted."
    adaption::T
    "Number of discretization steps."
    steps::Int64
    "Maximum number of discretization steps."
    stepsᵐᵃˣ::Int64
    "Desired Integration time"
    ∫dt::Float64
    function StepNumberTune(
        adaption::A, steps::Int64, stepsᵐᵃˣ::Int64, ∫dt::Float64
    ) where {A<:UpdateBool}
        return new{typeof(adaption)}(adaption, steps, stepsᵐᵃˣ, ∫dt)
    end
end

############################################################################################
function update!(
    number::StepNumberTune, adaption::BaytesCore.UpdateTrue, ϵ::F
) where {F<:AbstractFloat}
    @unpack stepsᵐᵃˣ, ∫dt = number
    #!NOTE:  max(10^-4, ϵ) instead of ϵ if (cornercase) ϵ so small that after division number would no longer be Int64 ~ 10^-5 more than enough for max steps reach
    stepsᵖ = max(1, Int(floor(∫dt / max(10^-3, ϵ)))) #At least 1 Leapfrog step, target Integration time of 1
    number.steps = min(stepsᵐᵃˣ, stepsᵖ)
    return nothing
end
function update!(
    number::StepNumberTune, adaption::BaytesCore.UpdateFalse, ϵ::F
) where {F<:AbstractFloat}
    return nothing
end
function update!(number::StepNumberTune, ϵ::F) where {F<:AbstractFloat}
    return update!(number, number.adaption, ϵ)
end

############################################################################################
#export
export StepNumberTune, update!


text:
using DPC, DataStructures, Plots, DifferentialEquations, Distributions, DataFrames
# using DifferentialEquations, DataFrames, DataStructures, Plots, ProgressMeter
include("./integrator/utils.jl")

const mypalette = (orange="#fcaf3e",
                   green="#8ae234",
                   purple="#ad7fa8",
                   turquise="#729f7e",
                   yellow="#fce94f",
                   blue="#729fcf",
                   red="#ef2929",
                   gray="#cccccc")

## Set up neuron
const num_inputs = 1000
const num_segments = 1000
const synapses_per_segment = 20
const transmission_probability = 0.5
const θ_seg = 300
const plateau_duration = 100e-3
const spike_duration = 5e-3
const refractory_duration = 5.01e-3
const trange = (0.0,3.0)
const λ = 200.0
const background_rate = 10.0
const input_names = Symbol.("in",1:num_inputs)
const seg_names = Symbol.("seg",1:num_segments)
const neuron_id = :integrator
const volleysize=10
const θ_syn = 5

net = Network(
    id_type=Symbol,
    time_type=Float64,
    weight_type=BernoulliSynapseWeight{Float64},
    synaptic_input_type=Int,
    default_spike_duration=spike_duration,
    default_refractory_duration=refractory_duration,
    default_plateau_duration=plateau_duration
)

objects = Dict{Symbol,Any}()
n = Neuron(neuron_id, net;
            θ_syn=0,
            θ_seg=θ_seg
           )
objects[neuron_id] = n

for input_id in input_names
    inp = Input(input_id, net)
    objects[input_id] = inp
end

for (i, seg_id) in enumerate(seg_names)
    seg = Segment(seg_id, n; θ_syn=5, θ_seg=0)

    syn_inputs = input_names[mod.(i .+ (0:synapses_per_segment), Ref(1:num_inputs))]
    for (j, syn_inp) in enumerate(syn_inputs)
        syn_id = Symbol("syn$(i)$(j)_" * String(syn_inp))
        source = objects[syn_inp]
        target = seg

        syn = Synapse(syn_id, source, target;
                      weight=BernoulliSynapseWeight(1.0,transmission_probability))
        objects[syn_id] = syn
    end
end

# save_network(net,"integrator/network.yaml")

# Set up stimulus

"""
Generate contiguous volleys at a high rate in 3 distinct periods.
"""
stimulus(t) = λ * if 0.25 ≤ t < 0.5
    1.0
elseif 0.75 ≤ t < 1.0
    2.0
elseif 1.25 ≤ t < 1.5
    3.0
elseif 0 ≤ t < 3.0
    0.1
else
    0.0
end

"""
Integral of the stimulus function (for plotting purposes).
"""
int_stimulus(t) = λ * if t < 0
    0
elseif 0.0 ≤ t < 0.25
    0.1(t-0.0)
elseif 0.25 ≤ t < 0.5
    1.0(t-0.25)+0.025
elseif 0.5 ≤ t < 0.75
    0.1(t-0.5)+0.25+0.025
elseif 0.75 ≤ t < 1.00
    2.0(t-0.75)+0.025+0.25+0.025
elseif 1.00 ≤ t < 1.25
    0.1(t-1.00)+0.5+0.025+0.25+0.025
elseif 1.25 ≤ t < 1.5
    3.0(t-1.25)+0.025+0.5+0.025+0.25+0.025
elseif 1.5 ≤ t < 3.5
    0.1(t-1.5)+0.75+0.025+0.5+0.025+0.25+0.025
else
    0.05+1.5+0.05+1.0+0.05+0.5+0.05
end

# ## Simulate!
spike_trains,v = generateContiguousVolleys(trange, num_inputs, stimulus, volleysize, background_rate; dt_spike = spike_duration)

all_spikes = Event{Float64}[]
for (i, pop_spikes) in enumerate(spike_trains)
    for spike_t in pop_spikes
        push!(all_spikes, Event(:input_spikes, 0.0, spike_t, objects[input_names[i]]) )
    end
end

# don't ask
v = [v]
spike_trains = [spike_trains]

# # simuluate the neuron
function fltr(t, tp, id, x)
    if tp == :plateau_starts || tp == :plateau_ends || tp == :spikes
        return true
    end
    return false
end
logger=simulate!(net, all_spikes; logger! =Logger(net,filter=fltr))
logger = logger.data

# # calculate plateau counts
# plateau_counts = sum.(eachrow(logger[r"seg*"]))

# # extract plateau strating times
plateau_starts = groupby(select!(filter(:event => (ev -> (ev == :plateau_starts)), logger),
                                :t,
                                :object),
                         :object)
plateau_starts_vec = [[ (seg_name,) in keys(plateau_starts.keymap) ?
                        plateau_starts[(seg_name,)][!,:t] : Float64[] for seg_name in seg_names]]
# # extract the neuron's spike times
neuron_spikes = filter(:event => (ev -> (ev == :spikes)), logger)[:,:t]

# ## Plot!

# Plot the inputs, active segments, the ideal "rectangular integrator", spike threshold in regions of interest
p_input = plot( [Shape([xoff, xoff+0.125, xoff+0.125, xoff],[0,0,750,750]) for xoff ∈ [0.25, 0.75, 1.25]],
                 linecolor=nothing,
                 fillcolor="#EEEEEE",
                 title="dendritic integration of evidence",
                 xlabel="time [s]",
                 xlims=trange,
                 ylims=(0,750),
                 legend=:topright,
                 label="",
                 xticks=[0,0.25,0.5,0.75,1.0,1.25,1.5,2,3])

plot!(stimulus, trange..., linewidth=2, color=mypalette.orange, label="volleys/s")

hline!([θ_seg], linewidth=2, color=mypalette.red, label="spike threshold")

# plot!([trange[1];logger.t;trange[2]],
#       [0;plateau_counts;0],
#       linewidth=2,
#       color=mypalette.blue,
#       seriestype=:steppost,
#       xlims=trange,
#       label="active segments")

plot!(t->(int_stimulus(t)-int_stimulus(t-100e-3))/100e-3,
      trange...,
      linewidth=2,
      color=:black,
      linestyle=:dash,
      label="rect. integrator")

# Plot two panels of spike-trains
v_cropped = [[(time=vv.time,idxs=filter(x->x∈1:50, vv.idxs)) for vv ∈ v[1]]]
spike_trains_cropped = [spike_trains[1][1:50]]
p_spikes1 = plot_spike_raster(trange,
                              spike_trains_cropped,
                              spike_duration;
                              colors=[:gray],
                              highlight=v_cropped,
                              highligh_color=mypalette.orange,
                              title="input spike-trains",
                              xlims=trange,
                              xticks=false)
yticks!(collect(-9.5:-20.0:-49.5), string.(collect(10:20:50)))
plot!(xaxis=:off)

v_cropped = [[(time=vv.time,idxs=filter(x->x∈951:1000, vv.idxs).-950) for vv ∈ v[1]]]
spike_trains_cropped = [spike_trains[1][951:1000]]
p_spikes2 = plot_spike_raster(trange,
                              spike_trains_cropped,
                              spike_duration;
                              colors=[:gray],
                              highlight=v_cropped,
                              highligh_color=mypalette.orange,
                              xlims=trange,
                              xticks=false)
yticks!(collect(-9.5:-20.0:-49.5), string.(collect(960:20:1000)))

# Plot the "zoomed-in" regions of interest
subs = Plots.Plot[]
for xoff ∈ [0.25, 0.75, 1.25]
    sub_range = (xoff,xoff+0.125)
    pp = plot(#stimulus, sub_range..., color=mypalette.orange, label="volley rate",
              background_color_inside="#EEEEEE",
              framestyle=:box,
              ylims=(0,750),
              xlims=sub_range,
              legend=false,
              xticks=[xoff, xoff+0.0625, xoff+0.125])
    hline!([θ_seg],linewidth=3,color=mypalette.red)
    vline!([neuron_spikes],linewidth=1,color=mypalette.purple)

    # plot!([sub_range[1];logger.t;sub_range[2]],
    #       [0;plateau_counts;0],
    #       color=mypalette.blue,
    #       linewidth=3,
    #       seriestype=:steppost,
    #       label="num. of plateaus")
    plot!(t->(int_stimulus(t)-int_stimulus(t-100e-3))/100e-3,
          sub_range...,
          color=:black,
          linewidth=3,
          linestyle=:dash,
          label="rect. integrator")
    push!(subs,pp)
end

# Compose and save the figure
l = @layout [
    a1{0.15h}
    a2{0.15h}
    b{0.3h}
    [c{0.33w} d{0.33w} e{0.33w}]
]
plt=plot(p_spikes1, p_spikes2, p_input, subs..., layout=l, size=(500,600))
savefig("examples/integrator/figures/integrator.svg")


text:
using Revise
using SocialNet
using BenchmarkTools

# This file demonstrates the required features of the implementation
# The same tests are included in the test suite in ./test/

# Create a list of persons
person_list = ["cow", "horse", "pig", "chicken", "dog", "cat", "goldfish"]

# Create an emtpy social network
net = Network()

# Add all the persons to the network
# The dot notation broadcasts over a list
add_person!.(net, person_list)

# List all persons
@show persons(net)

# Add friend relation to pairs
friend_pair!(net, "cow", "horse")
friend_pair!(net, "cow", "pig")
friend_pair!(net, "cow", "chicken")
friend_pair!(net, "cow", "goldfish")

friend_pair!(net, "horse", "pig")
friend_pair!(net, "horse", "chicken")

friend_pair!(net, "dog", "cat")
friend_pair!(net, "dog", "goldfish")

# List all friends of a given person
@show friends(net, "cow")
@show friends(net, "horse")
@show friends(net, "chicken")
@show friends(net, "pig")

@show friends(net, "dog")
@show friends(net, "cat")
@show friends(net, "goldfish")

# Compute degree of separation between pairs of persons
@show degree_of_separation(net, "cow", "cow")
@show degree_of_separation(net, "cow", "horse")
@show degree_of_separation(net, "horse", "cat")
@show degree_of_separation(net, "horse", "dog")
@show degree_of_separation(net, "horse", "goldfish")
@show degree_of_separation(net, "pig", "cat")

# Remove friend relation from pairs
defriend_pair!(net, "cow", "goldfish")

# The distance between pig and cat is now infinite,
# represented by typemax(Int)
@show degree_of_separation(net, "pig", "cat")

nothing


text:
""" DSLs for defining categories, diagrams, and related structures.

Here "diagram" means diagram in the standard category-theoretic sense, not
string diagram or wiring diagram. DSLs for constructing wiring diagrams are
provided by other submodules.
"""
module DiagrammaticPrograms
export @graph, @fincat, @finfunctor, @diagram, @free_diagram,
  @migrate, @migration

using Base.Iterators: repeated
using MLStyle: @match

using ...GAT, ...Present, ...Graphs, ...CategoricalAlgebra
using ...Theories: munit
using ...CategoricalAlgebra.FinCats: mapvals, make_map
using ...CategoricalAlgebra.DataMigrations: ConjQuery, GlueQuery, GlucQuery
using ...Graphs.BasicGraphs: TheoryGraph
import ...CategoricalAlgebra.FinCats: vertex_name, vertex_named,
  edge_name, edge_named

# Graphs
########

@present TheoryNamedGraph <: TheoryGraph begin
  VName::AttrType
  EName::AttrType
  vname::Attr(V, VName)
  ename::Attr(E, EName)
end

""" Abstract type for graph with named vertices and edges.
"""
@abstract_acset_type AbstractNamedGraph <: AbstractGraph

""" Graph with named vertices and edges.

The default graph type used by [`@graph`](@ref), [`@fincat`](@ref),
[`@diagram`](@ref), and related macros.
"""
@acset_type NamedGraph(TheoryNamedGraph, index=[:src,:tgt,:ename],
                       unique_index=[:vname]) <: AbstractNamedGraph
# FIXME: The edge name should also be uniquely indexed, but this currently
# doesn't play nicely with nullable attributes.

vertex_name(g::AbstractNamedGraph, args...) = subpart(g, args..., :vname)
edge_name(g::AbstractNamedGraph, args...) = subpart(g, args..., :ename)

vertex_named(g::AbstractNamedGraph, name) = only(incident(g, name, :vname))
edge_named(g::AbstractNamedGraph, name)= only(incident(g, name, :ename))

""" Construct a graph in a simple, declarative style.

The syntax is reminiscent of Graphviz. Each line a declares a vertex or set of
vertices, or an edge. For example, the following defines a directed triangle:

```julia
@graph begin
  v0, v1, v2
  fst: v0 → v1
  snd: v1 → v2
  comp: v0 → v2
end
```

Vertices in the graph must be uniquely named, whereas edges names are optional.
"""
macro graph(graph_type, body)
  :(parse_graph($(esc(graph_type)), $(Meta.quot(body))))
end
macro graph(body)
  :(parse_graph(NamedGraph{Symbol,Union{Symbol,Nothing}}, $(Meta.quot(body))))
end

function parse_graph(G::Type, body::Expr; preprocess::Bool=true)
  g = G()
  if preprocess
    body = reparse_arrows(body)
  end
  for stmt in statements(body)
    @match stmt begin
      # v
      v::Symbol => add_vertex!(g, vname=v)
      # u, v, ...
      Expr(:tuple, vs...) => add_vertices!(g, length(vs), vname=vs)
      # u → v
      Expr(:call, :(→), u::Symbol, v::Symbol) => parse_edge!(g, u, v)
      # e : u → v
      Expr(:call, :(:), e::Symbol, Expr(:call, :(→), u::Symbol, v::Symbol)) =>
        parse_edge!(g, u, v, ename=e)
      # (e, f, ...) : u → v
      Expr(:call, (:), Expr(:tuple, es...), Expr(:call, :(→), u::Symbol, v::Symbol)) =>
        for e in es; parse_edge!(g, u, v, ename=e) end
      ::LineNumberNode => nothing
      _ => error("Cannot parse line in graph definition: $stmt")
    end
  end
  return g
end

function parse_edge!(g, s::Symbol, t::Symbol;
                     ename::Union{Symbol,Nothing}=nothing)
  e = add_edge!(g, vertex_named(g, s), vertex_named(g, t))
  if has_subpart(g, :ename)
    g[e,:ename] = ename
  end
  return e
end

function reparse_arrows(expr)
  # `f : x → y` is parsed by Julia as `(f : x) → y`, not `f : (x → y)`.
  @match expr begin
    Expr(:call, :(→), Expr(:call, :(:), f, x), y) =>
      Expr(:call, :(:), f, Expr(:call, :(→), x, y))
    Expr(head, args...) => Expr(head, (reparse_arrows(arg) for arg in args)...)
    _ => expr
  end
end

statements(expr) = (expr isa Expr && expr.head == :block) ? expr.args : [expr]

# Categories
############

""" Present a category by generators and relations.

The result is a finitely presented category (`FinCat`) represented by a graph,
possibly with path equations. For example, the simplex category truncated to one
dimension is:

```julia
@fincat begin
  V, E
  (δ₀, δ₁): V → E
  σ₀: E → V

  σ₀ ∘ δ₀ == id(V)
  σ₀ ∘ δ₁ == id(V)
end
```

The objects and morphisms must be uniquely named.
"""
macro fincat(body)
  parse_category(NamedGraph{Symbol,Symbol}, body)
end

function parse_category(G::Type, body::Expr; preprocess::Bool=true)
  g, eqs = G(), Pair[]
  if preprocess
    body = reparse_arrows(body)
  end
  for stmt in statements(body)
    @match stmt begin
      # v
      v::Symbol => add_vertex!(g, vname=v)
      # u, v, ...
      Expr(:tuple, vs...) => add_vertices!(g, length(vs), vname=vs)
      # e : u → v
      Expr(:call, :(:), e::Symbol, Expr(:call, :(→), u::Symbol, v::Symbol)) =>
        parse_edge!(g, u, v, ename=e)
      # (e, f, ...) : u → v
      Expr(:call, (:), Expr(:tuple, es...), Expr(:call, :(→), u::Symbol, v::Symbol)) =>
        for e in es; parse_edge!(g, u, v, ename=e) end
      # f == g
      Expr(:call, :(==), lhs, rhs) => push!(eqs, parse_path_equation(g, lhs, rhs))
      ::LineNumberNode => nothing
      _ => error("Cannot parse line in category definition: $stmt")
    end
  end
  isempty(eqs) ? FinCat(g) : FinCat(g, eqs)
end

""" Parse an equation between paths in a named graph.
"""
parse_path_equation(g, lhs, rhs) = parse_path(g, lhs) => parse_path(g, rhs)

""" Parse a path in a named graph.
"""
function parse_path(g, expr)
  function parse(expr)
    @match expr begin
      Expr(:call, :compose, args...) => mapreduce(parse, vcat, args)
      Expr(:call, :(⋅), f, g) ||
      Expr(:call, :(⨟), f, g) => vcat(parse(f), parse(g))
      Expr(:call, :(∘), f, g) => vcat(parse(g), parse(f))
      Expr(:call, :id, x::Symbol) => empty(Path, g, vertex_named(g, x))
      f::Symbol => Path(g, edge_named(g, f))
      _ => error("Invalid morphism expression $expr")
    end
  end
  parse(expr)
end

# Functors
##########

""" Define a functor between two finitely presented categories.

Such a functor is defined by sending the object and morphism generators of the
domain category to generic object and morphism expressions in the codomain
category. For example, the following functor embeds the schema for graphs into
the schema for circular port graphs by ignoring the ports:

```julia
@finfunctor TheoryGraph ThCPortGraph begin
  V => Box
  E => Wire
  src => src ⨟ box
  tgt => tgt ⨟ box
end
```
"""
macro finfunctor(dom_cat, codom_cat, body)
  :(parse_functor($(esc(dom_cat)), $(esc(codom_cat)), $(Meta.quot(body))))
end

function parse_functor(C::FinCat, D::FinCat, body::Expr)
  ob_rhs, hom_rhs = parse_ob_hom_maps(C, body)
  F = FinFunctor(mapvals(x -> parse_ob(D, x), ob_rhs),
                 mapvals(f -> parse_hom(D, f), hom_rhs), C, D)
  is_functorial(F, check_equations=false) ||
    error("Parsed functor is not functorial: $body")
  return F
end
function parse_functor(C::Presentation, D::Presentation, body::Expr; kw...)
  parse_functor(FinCat(C), FinCat(D), body; kw...)
end

function parse_ob_hom_maps(C::FinCat, body::Expr; allow_missing::Bool=false)
  assignments = Dict{Symbol,Union{Expr,Symbol}}()
  assign(lhs, rhs) = if haskey(assignments, lhs)
    error("Left-hand side $lhs assigned twice in $body")
  else
    assignments[lhs] = rhs
  end
  for stmt in statements(body)
    @match stmt begin
      # x => y
      Expr(:call, :(=>), lhs::Symbol, rhs) => assign(lhs, rhs)
      # (x, x′, ...) => y
      Expr(:call, :(=>), Expr(:tuple, args...), rhs) =>
        foreach(lhs -> assign(lhs, rhs), args)
      ::LineNumberNode => nothing
      _ => error("Invalid assignment statement $stmt")
    end
  end
  ob_rhs = make_map(ob_generators(C)) do x
    y = pop!(assignments, ob_generator_name(C, x), missing)
    (!ismissing(y) || allow_missing) ? y :
      error("Object $(ob_generator_name(C,x)) is not assigned")
  end
  hom_rhs = make_map(hom_generators(C)) do f
    g = pop!(assignments, hom_generator_name(C, f), missing)
    (!ismissing(g) || allow_missing) ? g :
      error("Morphism $(hom_generator_name(C,f)) is not assigned")
  end
  isempty(assignments) ||
    error(string("Unused assignment(s): ", join(keys(assignments), ", ")))
  (ob_rhs, hom_rhs)
end

""" Parse expression for object in a category.
"""
function parse_ob(C::FinCat{Ob,Hom}, expr) where {Ob,Hom}
  @match expr begin
    x::Symbol => ob_generator(C, x)
    Expr(:curly, _...) => parse_gat_expr(C, expr)::Ob
    _ => error("Invalid object expression $expr")
  end
end

""" Parse expression for morphism in a category.
"""
function parse_hom(C::FinCat{Ob,Hom}, expr) where {Ob,Hom}
  function parse(expr)
    @match expr begin
      Expr(:call, :compose, args...) =>
        mapreduce(parse, (fs...) -> compose(C, fs...), args)
      Expr(:call, :(⋅), f, g) ||
      Expr(:call, :(⨟), f, g) => compose(C, parse(f), parse(g))
      Expr(:call, :(∘), f, g) => compose(C, parse(g), parse(f))
      Expr(:call, :id, x::Symbol) => id(C, ob_generator(C, x))
      f::Symbol => hom_generator(C, f)
      Expr(:curly, _...) => parse_gat_expr(C, expr)::Hom
      _ => error("Invalid morphism expression $expr")
    end
  end
  parse(expr)
end

""" Parse GAT expression based on curly braces, rather than parentheses.
"""
function parse_gat_expr(C::FinCat, root_expr)
  pres = presentation(C)
  function parse(expr)
    @match expr begin
      Expr(:curly, head::Symbol, args...) =>
        invoke_term(pres.syntax, head, map(parse, args)...)
      x::Symbol => generator(pres, x)
      _ => error("Invalid GAT expression $root_expr")
    end
  end
  parse(root_expr)
end

# Diagrams
##########

""" Present a diagram in a given category.

Recall that a *diagram* in a category ``C`` is a functor ``F: J → C`` from a
small category ``J`` into ``C``. Given the category ``C``, this macro presents a
diagram in ``C``, i.e., constructs a finitely presented indexing category ``J``
together with a functor ``F: J → C``. This method of simultaneous definition is
often more convenient than defining ``J`` and ``F`` separately, as could be
accomplished by calling [`@fincat`](@ref) and then [`@finfunctor`](@ref).

As an example, the limit of the following diagram consists of the paths of
length two in a graph:

```julia
@diagram TheoryGraph begin
  v::V
  (e₁, e₂)::E
  (t: e₁ → v)::tgt
  (s: e₂ → v)::src
end
```

Morphisms in the indexing category can be left unnamed, which is convenient for
defining free diagrams (see also [`@free_diagram`](@ref)). For example, the
following diagram is isomorphic to the previous one:

```julia
@diagram TheoryGraph begin
  v::V
  (e₁, e₂)::E
  (e₁ → v)::tgt
  (e₂ → v)::src
end
```

Of course, unnamed morphisms cannot be referenced by name within the `@diagram`
call or in other settings, which can sometimes be problematic.
"""
macro diagram(cat, body)
  :(parse_diagram($(esc(cat)), $(Meta.quot(body))))
end

""" Present a free diagram in a given category.

Recall that a *free diagram* in a category ``C`` is a functor ``F: J → C`` where
``J`` is a free category on a graph, here assumed finite. This macro is
functionally a special case of [`@diagram`](@ref) but, for convenience, changes
the interpretation of equality expressions. Rather than interpreting them as
equations between morphisms in ``J``, equality expresions can be used to
introduce anonymous morphisms in a "pointful" style. For example, the limit of
the following diagram consists of the paths of length two in a graph:

```julia
@free_diagram TheoryGraph begin
  v::V
  (e₁, e₂)::E
  tgt(e₁) == v
  src(e₂) == v
end
```

Anonymous objects can also be introduced. For example, the previous diagram is
isomorphic to this one:

```julia
@free_diagram TheoryGraph begin
  (e₁, e₂)::E
  tgt(e₁) == src(e₂)
end
```

Some care must exercised when defining morphisms between diagrams with anonymous
objects, since they cannot be referred to by name.
"""
macro free_diagram(cat, body)
  :(parse_diagram($(esc(cat)), $(Meta.quot(body)), free=true))
end

function parse_diagram(C::FinCat, body::Expr; kw...)
  F_ob, F_hom, J = parse_diagram_data(
    C, x -> parse_ob(C,x), (f,x,y) -> parse_hom(C,f), body; kw...)
  F = FinFunctor(F_ob, F_hom, J, C)
  is_functorial(F, check_equations=false) ||
    error("Parsed diagram is not functorial: $body")
  return F
end
parse_diagram(pres::Presentation, body::Expr; kw...) =
  parse_diagram(FinCat(pres), body; kw...)

function parse_diagram_data(C::FinCat, parse_ob, parse_hom, body::Expr;
                            free::Bool=false, preprocess::Bool=true)
  g, eqs = NamedGraph{Symbol,Union{Symbol,Nothing}}(), Pair[]
  F_ob, F_hom = [], []
  function push_hom!(h, x, y; name=nothing)
    e = parse_edge!(g, x, y, ename=name)
    push!(F_hom, parse_hom(h, F_ob[src(g,e)], F_ob[tgt(g,e)]))
    return e
  end

  nanon = 0
  if preprocess
    body = reparse_arrows(body)
  end
  for stmt in statements(body)
    @match stmt begin
      # x => X
      # x::X
      Expr(:call, :(=>), x::Symbol, X) || Expr(:(::), x::Symbol, X) => begin
        add_vertex!(g, vname=x)
        push!(F_ob, parse_ob(X))
      end
      # (x, y, ...) => X
      # (x, y, ...)::X
      Expr(:call, :(=>), Expr(:tuple, xs...), X) ||
      Expr(:(::), Expr(:tuple, xs...), X) => begin
        add_vertices!(g, length(xs), vname=xs)
        append!(F_ob, repeated(parse_ob(X), length(xs)))
      end
      # (f: x → y) => h
      # (f: x → y)::h
      Expr(:call, :(=>), Expr(:call, :(:), f::Symbol,
                              Expr(:call, :(→), x::Symbol, y::Symbol)), h) ||
      Expr(:(::), Expr(:call, :(:), f::Symbol,
                       Expr(:call, :(→), x::Symbol, y::Symbol)), h) =>
        push_hom!(h, x, y, name=f)
      # (x → y) => h
      # (x → y)::h
      Expr(:call, :(=>), Expr(:call, :(→), x::Symbol, y::Symbol), h) ||
      Expr(:(::), Expr(:call, :(→), x::Symbol, y::Symbol), h) =>
        push_hom!(h, x, y)
      # h(x) == y
      # y == h(x)
      (Expr(:call, :(==), call::Expr, y::Symbol) ||
       Expr(:call, :(==), y::Symbol, call::Expr)) && if free end =>
         push_hom!(destructure_unary_call(call)..., y)
      # h(x) == k(y)
      Expr(:call, :(==), lhs::Expr, rhs::Expr) && if free end => begin
        (h, x), (k, y) = destructure_unary_call(lhs), destructure_unary_call(rhs)
        z = Symbol("##unnamed#$(nanon += 1)")
        add_vertex!(g, vname=z)
        push!(F_ob, nothing) # Assumes that codomain not needed to parse homs.
        e1, e2 = push_hom!(h, x, z), push_hom!(k, y, z)
        # Infer codomain from parsed homs.
        Z1, Z2 = codom(C, F_hom[e1]), codom(C, F_hom[e2])
        Z1 == Z2 || error("Codomain mismatch in $stmt: $Z1 != $Z2")
        F_ob[end] = Z1
      end
      # f == g
      Expr(:call, :(==), lhs, rhs) && if !free end =>
        push!(eqs, parse_path_equation(g, lhs, rhs))
      ::LineNumberNode => nothing
      _ => error("Cannot parse line in diagram definition: $stmt")
    end
  end
  J = isempty(eqs) ? FinCat(g) : FinCat(g, eqs)
  (F_ob, F_hom, J)
end

# Data migrations
#################

""" A diagram without a codomain category.

An intermediate data representation used internally by the parser for the
[`@migration`](@ref) macro.
"""
struct DiagramData{T,ObMap,HomMap,Shape<:FinCat}
  ob_map::ObMap
  hom_map::HomMap
  shape::Shape

  function DiagramData{T}(ob_map::ObMap, hom_map::HomMap, shape::Shape) where
      {T,ObMap,HomMap,Shape<:FinCat}
    new{T,ObMap,HomMap,Shape}(ob_map, hom_map, shape)
  end
end

Diagrams.ob_map(d::DiagramData, x) = d.ob_map[x]
Diagrams.hom_map(d::DiagramData, f) = d.hom_map[f]
Diagrams.shape(d::DiagramData) = d.shape

""" A diagram morphism without a domain or codomain.

Like [`DiagramData`](@ref), an intermediate data representation used internally
by the parser for the [`@migration`](@ref) macro.
"""
struct DiagramHomData{T,ObMap,HomMap}
  ob_map::ObMap
  hom_map::HomMap

  function DiagramHomData{T}(ob_map::ObMap, hom_map::HomMap) where {T,ObMap,HomMap}
    new{T,ObMap,HomMap}(ob_map, hom_map)
  end
end

""" Contravariantly migrate data from one acset to another.

This macro is shorthand for defining a data migration using the
[`@migration`](@ref) macro and then calling the `migrate` function. If the
migration will be used multiple times, it is more efficient to perform these
steps separately, reusing the functor defined by `@migration`.

For more about the syntax and supported features, see [`@migration`](@ref).
"""
macro migrate(tgt_type, src_acset, body)
  quote
    let T = $(esc(tgt_type)), X = $(esc(src_acset))
      migrate(T, X, parse_migration(Presentation(T), Presentation(X),
                                    $(Meta.quot(body))))
    end
  end
end

""" Define a contravariant data migration.

This macro provides a DSL to specify a contravariant data migration from
``C``-sets to ``D``-sets for given schemas ``C`` and ``D``. A data migration is
defined by a functor from ``D`` to a category of queries on ``C``. Thus, every
object of ``D`` is assigned a query on ``C`` and every morphism of ``D`` is
assigned a morphism of queries, in a compatible way. Example usages are in the
unit tests. What follows is a technical reference.

Several categories of queries are supported by this macro:

1. Trivial queries, specified by a single object of ``C``. In this case, the
   macro simply defines a functor ``D → C`` and is equivalent to
   [`@finfunctor`](@ref) or [`@diagram`](@ref).
2. *Conjunctive queries*, specified by a diagram in ``C`` and evaluated as a
   finite limit.
3. *Gluing queries*, specified by a diagram in ``C`` and evaluated as a finite
   colimit. An important special case is *linear queries*, evaluated as a
   finite coproduct.
4. *Gluc queries* (gluings of conjunctive queries), specified by a diagram of
   diagrams in ``C`` and evaluated as a colimit of limits. An important special
   case is *duc queries* (disjoint unions of conjunctive queries), evaluated as
   a coproduct of limits.

The query category of the data migration is not specified explicitly but is
inferred from the queries used in the macro call. Implicit conversion is
performed: trivial queries can be coerced to conjunctive queries or gluing
queries, and conjunctive queries and gluing queries can both be coerced to gluc
queries. Due to the implicit conversion, the resulting functor out of ``D`` has
a single query type and thus a well-defined codomain.

Syntax for the right-hand sides of object assignments is:

- a symbol, giving object of ``C`` (query type: trivial)
- `@product ...` (query type: conjunctive)
- `@join ...` (alias: `@limit`, query type: conjunctive)
- `@cases ...` (alias: `@coproduct`, query type: gluing or gluc)
- `@glue ...` (alias: `@colimit`, query type: gluing or gluc)

Thes query types supported by this macro generalize the kind of queries familiar
from relational databases. Less familiar is the concept of a morphism between
queries, derived from the concept of a morphism between diagrams in a category.
A query morphism is given by a functor between the diagrams' indexing categories
together with a natural transformation filling a triangle of the appropriate
shape. From a practical standpoint, the most important thing to remember is that
a morphism between conjunctive queries is contravariant with respect to the
diagram shapes, whereas a morphism between gluing queries is covariant. TODO:
Reference for more on this.
"""
macro migration(src_schema, body)
  :(parse_migration($(esc(src_schema)), $(Meta.quot(body))))
end
macro migration(tgt_schema, src_schema, body)
  :(parse_migration($(esc(tgt_schema)), $(esc(src_schema)), $(Meta.quot(body))))
end

""" Parse a contravariant data migration from a Julia expression.

The process kicked off by this internal function is somewhat complicated due to
the need to coerce queries and query morphisms to a common category. The
high-level steps of this process are:

1. Parse the queries and query morphisms into intermediate representations
   ([`DiagramData`](@ref) and [`DiagramHomData`](@ref)) whose final types are
   not yet determined.
2. Promote the query types to the tightest type encompassing all queries, an
   approach reminiscent of Julia's own type promotion system.
3. Convert all query and query morphisms to this common type, yielding `Diagram`
   and `DiagramHom` instances.
"""
function parse_migration(src_schema::Presentation, body::Expr;
                         preprocess::Bool=true)
  C = FinCat(src_schema)
  F_ob, F_hom, J = parse_query_diagram(C, body; free=false, preprocess=preprocess)
  make_migration_functor(F_ob, F_hom, J, C)
end
function parse_migration(tgt_schema::Presentation, src_schema::Presentation,
                         body::Expr; preprocess::Bool=true)
  D, C = FinCat(tgt_schema), FinCat(src_schema)
  if preprocess
    body = reparse_arrows(body)
  end
  ob_rhs, hom_rhs = parse_ob_hom_maps(D, body)
  F_ob = mapvals(expr -> parse_query(C, expr), ob_rhs)
  F_hom = mapvals(hom_rhs, keys=true) do f, expr
    parse_query_hom(C, expr, F_ob[dom(D,f)], F_ob[codom(D,f)])
  end
  make_migration_functor(F_ob, F_hom, D, C)
end

# Query parsing
#--------------

""" Parse expression defining a query.
"""
function parse_query(C::FinCat, expr)
  @match expr begin
    x::Symbol => ob_generator(C, x)
    Expr(:macrocall, form, args...) &&
        if form ∈ (Symbol("@limit"), Symbol("@join")) end => begin
      DiagramData{op}(parse_query_diagram(C, last(args))...)
    end
    Expr(:macrocall, form, args...) &&
        if form == Symbol("@product") end => begin
      d = DiagramData{op}(parse_query_diagram(C, last(args))...)
      is_discrete(shape(d)) ? d : error("Product query is not discrete: $expr")
    end
    Expr(:macrocall, form, args...) &&
        if form ∈ (Symbol("@colimit"), Symbol("@glue")) end => begin
      DiagramData{id}(parse_query_diagram(C, last(args))...)
    end
    Expr(:macrocall, form, args...) &&
        if form ∈ (Symbol("@coproduct"), Symbol("@cases")) end => begin
      d = DiagramData{id}(parse_query_diagram(C, last(args))...)
      is_discrete(shape(d)) ? d : error("Cases query is not discrete: $expr")
    end
    _ => error("Cannot parse query in migration definition: $expr")
  end
end
function parse_query_diagram(C::FinCat, expr::Expr;
                             free::Bool=true, preprocess::Bool=false)
  parse_diagram_data(C, X -> parse_query(C,X),
                     (f,x,y) -> parse_query_hom(C,f,x,y), expr;
                     free=free, preprocess=preprocess)
end

""" Parse expression defining a morphism of queries.
"""
parse_query_hom(C::FinCat{Ob}, expr, ::Ob, ::Union{Ob,Nothing}) where Ob =
  parse_hom(C, expr)

# Conjunctive fragment.

function parse_query_hom(C::FinCat, expr, d::DiagramData{op}, d′::DiagramData{op})
  ob_rhs, hom_rhs = parse_ob_hom_maps(shape(d′), expr)
  f_ob = mapvals(ob_rhs, keys=true) do j′, rhs
    parse_conj_query_ob_rhs(C, rhs, d, ob_map(d′, j′))
  end
  f_hom = mapvals(rhs -> parse_hom(shape(d), rhs), hom_rhs)
  DiagramHomData{op}(f_ob, f_hom)
end
function parse_query_hom(C::FinCat{Ob}, expr, c::Ob, d′::DiagramData{op}) where Ob
  ob_rhs, f_hom = parse_ob_hom_maps(shape(d′), expr, allow_missing=true)
  f_ob = mapvals(ob_rhs, keys=true) do j′, rhs
    ismissing(rhs) ? missing : (missing, parse_query_hom(C, rhs, c, ob_map(d′, j′)))
  end
  @assert all(ismissing, f_hom)
  DiagramHomData{op}(f_ob, f_hom)
end
function parse_query_hom(C::FinCat{Ob}, expr, d::DiagramData{op}, c′::Ob) where Ob
  DiagramHomData{op}([parse_conj_query_ob_rhs(C, expr, d, c′)], [])
end

# Gluing fragment.

function parse_query_hom(C::FinCat, expr, d::DiagramData{id}, d′::DiagramData{id})
  ob_rhs, hom_rhs = parse_ob_hom_maps(shape(d), expr)
  f_ob = mapvals(ob_rhs, keys=true) do j, rhs
    parse_glue_query_ob_rhs(C, rhs, ob_map(d, j), d′)
  end
  f_hom = mapvals(expr -> parse_hom(shape(d′), expr), hom_rhs)
  DiagramHomData{id}(f_ob, f_hom)
end
function parse_query_hom(C::FinCat{Ob}, expr, c::Union{Ob,DiagramData{op}},
                         d′::DiagramData{id}) where Ob
  DiagramHomData{id}([parse_glue_query_ob_rhs(C, expr, c, d′)], [])
end
function parse_query_hom(C::FinCat{Ob}, expr, d::DiagramData{id},
                         c′::Union{Ob,DiagramData{op}}) where Ob
  ob_rhs, f_hom = parse_ob_hom_maps(shape(d), expr, allow_missing=true)
  f_ob = mapvals(ob_rhs, keys=true) do j, rhs
    ismissing(rhs) ? missing : (missing, parse_query_hom(C, rhs, ob_map(d, j), c′))
  end
  @assert all(ismissing, f_hom)
  DiagramHomData{id}(f_ob, f_hom)
end

""" Parse RHS of object assignment in morphism out of conjunctive query.
"""
function parse_conj_query_ob_rhs(C::FinCat, expr, d::DiagramData{op}, c′)
  j_name, f_expr = @match expr begin
    x::Symbol => (x, nothing)
    Expr(:tuple, x::Symbol, f) => (x, f)
    Expr(:call, op, _...) && if op ∈ compose_ops end =>
      leftmost_arg(expr, (:(⋅), :(⨟)), all_ops=compose_ops)
    Expr(:call, name::Symbol, _) => reverse(destructure_unary_call(expr))
    _ => error("Cannot parse object assignment in migration: $expr")
  end
  j = ob_generator(shape(d), j_name)
  isnothing(f_expr) ? j :
    (j, parse_query_hom(C, f_expr, ob_map(d, j), c′))
end

""" Parse RHS of object assignment in morphism into gluing query.
"""
function parse_glue_query_ob_rhs(C::FinCat, expr, c, d′::DiagramData{id})
  j′_name, f_expr = @match expr begin
    x::Symbol => (x, nothing)
    Expr(:tuple, x::Symbol, f) => (x, f)
    Expr(:call, op, _...) && if op ∈ compose_ops end =>
      leftmost_arg(expr, (:(∘),), all_ops=compose_ops)
    _ => error("Cannot parse object assignment in migration: $expr")
  end
  j′ = ob_generator(shape(d′), j′_name)
  isnothing(f_expr) ? j′ :
    (j′, parse_query_hom(C, f_expr, c, ob_map(d′, j′)))
end

const compose_ops = (:(⋅), :(⨟), :(∘))

# Query construction
#-------------------

function make_migration_functor(F_ob, F_hom, D::FinCat, C::FinCat)
  diagram(make_query(C, DiagramData{id}(F_ob, F_hom, D)))
end

function make_query(C::FinCat{Ob}, d::DiagramData{T}) where {T, Ob}
  F_ob, F_hom, J = d.ob_map, d.hom_map, shape(d)
  F_ob = mapvals(x -> make_query(C, x), F_ob)
  query_type = mapreduce(typeof, promote_query_type, values(F_ob), init=Ob)
  @assert query_type != Any
  F_ob = mapvals(x -> convert_query(C, query_type, x), F_ob)
  F_hom = mapvals(F_hom, keys=true) do h, f
    make_query_hom(f, F_ob[dom(J,h)], F_ob[codom(J,h)])
  end
  Diagram{T}(if query_type <: Ob
    FinFunctor(F_ob, F_hom, shape(d), C)
  else
    # XXX: Why is the element type of `F_ob` sometimes too loose?
    D = TypeCat(typeintersect(query_type, eltype(values(F_ob))),
                eltype(values(F_hom)))
    FinDomFunctor(F_ob, F_hom, shape(d), D)

  end)
end

make_query(C::FinCat{Ob}, x::Ob) where Ob = x

function make_query_hom(f::DiagramHomData{op}, d::Diagram{op}, d′::Diagram{op})
  f_ob = mapvals(f.ob_map, keys=true) do j′, x
    x = @match x begin
      ::Missing => only_ob(shape(d))
      (::Missing, g) => (only_ob(shape(d)), g)
      _ => x
    end
    @match x begin
      (j, g) => Pair(j, make_query_hom(g, ob_map(d, j), ob_map(d′, j′)))
      j => j
    end
  end
  f_hom = mapvals(h -> ismissing(h) ? only_hom(shape(d)) : h, f.hom_map)
  DiagramHom{op}(f_ob, f_hom, d, d′)
end

function make_query_hom(f::DiagramHomData{id}, d::Diagram{id}, d′::Diagram{id})
  f_ob = mapvals(f.ob_map, keys=true) do j, x
    x = @match x begin
      ::Missing => only_ob(shape(d′))
      (::Missing, g) => (only_ob(shape(d′)), g)
      _ => x
    end
    @match x begin
      (j′, g) => Pair(j′, make_query_hom(g, ob_map(d, j), ob_map(d′, j′)))
      j′ => j′
    end
  end
  f_hom = mapvals(h -> ismissing(h) ? only_hom(shape(d′)) : h, f.hom_map)
  DiagramHom{id}(f_ob, f_hom, d, d′)
end

function make_query_hom(f::Hom, d::Diagram{T,C}, d′::Diagram{T,C}) where
    {T, Ob, Hom, C<:FinCat{Ob,Hom}}
  cat = codom(diagram(d))
  munit(DiagramHom{T}, cat, f, dom_shape=shape(d), codom_shape=shape(d′))
end
make_query_hom(f, x, y) = f

only_ob(C::FinCat) = only(ob_generators(C))
only_hom(C::FinCat) = (@assert is_discrete(C); id(C, only_ob(C)))

# Query promotion
#----------------

# Promotion of query types is modeled loosely on Julia's type promotion system:
# https://docs.julialang.org/en/v1/manual/conversion-and-promotion/

promote_query_rule(::Type, ::Type) = Union{}
promote_query_rule(::Type{<:ConjQuery{C}}, ::Type{<:Ob}) where {Ob,C<:FinCat{Ob}} =
  ConjQuery{C}
promote_query_rule(::Type{<:GlueQuery{C}}, ::Type{<:Ob}) where {Ob,C<:FinCat{Ob}} =
  GlueQuery{C}
promote_query_rule(::Type{<:GlucQuery{C}}, ::Type{<:Ob}) where {Ob,C<:FinCat{Ob}} =
  GlucQuery{C}
promote_query_rule(::Type{<:GlucQuery{C}}, ::Type{<:ConjQuery{C}}) where C =
  GlucQuery{C}
promote_query_rule(::Type{<:GlucQuery{C}}, ::Type{<:GlueQuery{C}}) where C =
  GlucQuery{C}

promote_query_type(T, S) = promote_query_result(
  T, S, Union{promote_query_rule(T,S), promote_query_rule(S,T)})
promote_query_result(T, S, ::Type{Union{}}) = typejoin(T, S)
promote_query_result(T, S, U) = U

convert_query(::FinCat, ::Type{T}, x::S) where {T, S<:T} = x

function convert_query(cat::C, ::Type{<:Diagram{T,C}}, x::Ob) where
  {T, Ob, C<:FinCat{Ob}}
  g = NamedGraph{Symbol,Symbol}(1, vname=nameof(x))
  munit(Diagram{T}, cat, x, shape=FinCat(g))
end
function convert_query(::C, ::Type{<:GlucQuery{C}}, d::ConjQuery{C}) where C
  munit(Diagram{id}, TypeCat(ConjQuery{C}, Any), d)
end
function convert_query(cat::C, ::Type{<:GlucQuery{C}}, d::GlueQuery{C}) where C
  J = shape(d)
  new_ob = make_map(ob_generators(J)) do j
    convert_query(cat, ConjQuery{C}, ob_map(d, j))
  end
  new_hom = make_map(hom_generators(J)) do h
    munit(Diagram{op}, cat, hom_map(d, h),
          dom_shape=new_ob[dom(J,h)], codom_shape=new_ob[codom(J,h)])
  end
  Diagram{id}(FinDomFunctor(new_ob, new_hom, J))
end
function convert_query(cat::C, ::Type{<:GlucQuery{C}}, x::Ob) where
    {Ob, C<:FinCat{Ob}}
  convert_query(cat, GlucQuery{C}, convert_query(cat, ConjQuery{C}, x))
end

# Utilities
###########

""" Left-most argument plus remainder of left-associated binary operations.
"""
function leftmost_arg(expr, ops; all_ops=nothing)
  isnothing(all_ops) && (all_ops = ops)
  function leftmost(expr)
    @match expr begin
      Expr(:call, op2, Expr(:call, op1, x, y), z) &&
          if op1 ∈ all_ops && op2 ∈ all_ops end => begin
        x, rest = leftmost(Expr(:call, op1, x, y))
        (x, Expr(:call, op2, rest, z))
      end
      Expr(:call, op, x, y) && if op ∈ ops end => (x, y)
      _ => (nothing, expr)
    end
  end
  leftmost(expr)
end

""" Destructure the expression `:(f(g(x)))` to `(:(f∘g), :x)`, for example.
"""
function destructure_unary_call(expr::Expr)
  @match expr begin
    Expr(:call, head, x::Symbol) => (head, x)
    Expr(:call, head, arg) => begin
      rest, x = destructure_unary_call(arg)
      (Expr(:call, :(∘), head, rest), x)
    end
  end
end

end

max_forks_repo_name: epatters/CompCat


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct TGTypeI <: TurbineGov
        R::Float64
        Ts::Float64
        Tc::Float64
        T3::Float64
        T4::Float64
        T5::Float64
        valve_position_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
        P_ref::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        internal::InfrastructureSystemsInternal
    end

Parameters of a Turbine Governor Type I.

# Arguments
- `R::Float64`: Droop parameter, validation range: `(0, nothing)`
- `Ts::Float64`: Governor time constant, validation range: `(0, nothing)`
- `Tc::Float64`: Servo time constant, validation range: `(0, nothing)`
- `T3::Float64`: Transient gain time constant, validation range: `(0, nothing)`
- `T4::Float64`: Power fraction time constant, validation range: `(0, nothing)`
- `T5::Float64`: Reheat time constant, validation range: `(0, nothing)`
- `valve_position_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}`: Valve position limits in MW
- `P_ref::Float64`: Reference Power Set-point, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`: The states of the TGTypeI model are:
	x_g1: Governor state,
	x_g2: Servo state,
	x_g3: Reheat state
- `n_states::Int`: TGTypeI has 3 states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct TGTypeI <: TurbineGov
    "Droop parameter"
    R::Float64
    "Governor time constant"
    Ts::Float64
    "Servo time constant"
    Tc::Float64
    "Transient gain time constant"
    T3::Float64
    "Power fraction time constant"
    T4::Float64
    "Reheat time constant"
    T5::Float64
    "Valve position limits in MW"
    valve_position_limits::NamedTuple{(:min, :max), Tuple{Float64, Float64}}
    "Reference Power Set-point"
    P_ref::Float64
    ext::Dict{String, Any}
    "The states of the TGTypeI model are:
	x_g1: Governor state,
	x_g2: Servo state,
	x_g3: Reheat state"
    states::Vector{Symbol}
    "TGTypeI has 3 states"
    n_states::Int
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function TGTypeI(R, Ts, Tc, T3, T4, T5, valve_position_limits, P_ref=1.0, ext=Dict{String, Any}(), )
    TGTypeI(R, Ts, Tc, T3, T4, T5, valve_position_limits, P_ref, ext, [:x_g1, :x_g2, :x_g3], 3, InfrastructureSystemsInternal(), )
end

function TGTypeI(; R, Ts, Tc, T3, T4, T5, valve_position_limits, P_ref=1.0, ext=Dict{String, Any}(), states=[:x_g1, :x_g2, :x_g3], n_states=3, internal=InfrastructureSystemsInternal(), )
    TGTypeI(R, Ts, Tc, T3, T4, T5, valve_position_limits, P_ref, ext, states, n_states, internal, )
end

# Constructor for demo purposes; non-functional.
function TGTypeI(::Nothing)
    TGTypeI(;
        R=0,
        Ts=0,
        Tc=0,
        T3=0,
        T4=0,
        T5=0,
        valve_position_limits=(min=0.0, max=0.0),
        P_ref=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`TGTypeI`](@ref) `R`."""
get_R(value::TGTypeI) = value.R
"""Get [`TGTypeI`](@ref) `Ts`."""
get_Ts(value::TGTypeI) = value.Ts
"""Get [`TGTypeI`](@ref) `Tc`."""
get_Tc(value::TGTypeI) = value.Tc
"""Get [`TGTypeI`](@ref) `T3`."""
get_T3(value::TGTypeI) = value.T3
"""Get [`TGTypeI`](@ref) `T4`."""
get_T4(value::TGTypeI) = value.T4
"""Get [`TGTypeI`](@ref) `T5`."""
get_T5(value::TGTypeI) = value.T5
"""Get [`TGTypeI`](@ref) `valve_position_limits`."""
get_valve_position_limits(value::TGTypeI) = value.valve_position_limits
"""Get [`TGTypeI`](@ref) `P_ref`."""
get_P_ref(value::TGTypeI) = value.P_ref
"""Get [`TGTypeI`](@ref) `ext`."""
get_ext(value::TGTypeI) = value.ext
"""Get [`TGTypeI`](@ref) `states`."""
get_states(value::TGTypeI) = value.states
"""Get [`TGTypeI`](@ref) `n_states`."""
get_n_states(value::TGTypeI) = value.n_states
"""Get [`TGTypeI`](@ref) `internal`."""
get_internal(value::TGTypeI) = value.internal

"""Set [`TGTypeI`](@ref) `R`."""
set_R!(value::TGTypeI, val) = value.R = val
"""Set [`TGTypeI`](@ref) `Ts`."""
set_Ts!(value::TGTypeI, val) = value.Ts = val
"""Set [`TGTypeI`](@ref) `Tc`."""
set_Tc!(value::TGTypeI, val) = value.Tc = val
"""Set [`TGTypeI`](@ref) `T3`."""
set_T3!(value::TGTypeI, val) = value.T3 = val
"""Set [`TGTypeI`](@ref) `T4`."""
set_T4!(value::TGTypeI, val) = value.T4 = val
"""Set [`TGTypeI`](@ref) `T5`."""
set_T5!(value::TGTypeI, val) = value.T5 = val
"""Set [`TGTypeI`](@ref) `valve_position_limits`."""
set_valve_position_limits!(value::TGTypeI, val) = value.valve_position_limits = val
"""Set [`TGTypeI`](@ref) `P_ref`."""
set_P_ref!(value::TGTypeI, val) = value.P_ref = val
"""Set [`TGTypeI`](@ref) `ext`."""
set_ext!(value::TGTypeI, val) = value.ext = val


text:
using SMC
using PyPlot
using Base.Test

@testset "logsumexp tests" begin

    # on vector
    vec = [1., 2.2, 3.1]
    actual = logsumexp(vec)
    expected = log(sum(exp.(vec)))
    @test isapprox(actual, expected)

    # on 2D array
    arr = [[1. 2. 3.]; [3.3 4.1 5.2]]
    # sum over rows
    actual = logsumexp(arr, 1)
    expected = [3.3955454645979626 4.215519523179755 5.305083319768697]
    @test isapprox(actual, expected)

    # sum over columns
    actual = logsumexp(arr, 2)
    expected = [3.40760596444438 5.593689178496464]'
    @test isapprox(actual, expected)
end

@testset "hmm tests" begin

    initial_state_prior = [0.4, 0.6]
    transition_model = [0.7 0.3;
                        0.1 0.9]
    observation_model = [0.6 0.3 0.1;
                         0.1 0.4 0.5]
    hmm = HiddenMarkovModel(initial_state_prior,
                            transition_model,
                            observation_model)
    observations = [1, 3, 2]

    @testset "forward probabilities" begin
        # p(x_t | y_{1:t-1}) for each t = 1,..,T
        actual = SMC.hmm_forward_pass(hmm, observations)
        # number of time steps x number of states
        @test size(actual) == (3, 2)

        # check against manually computed values
        expected = begin

            # p(x_2 = 1, y_1 = 1, y_2 = 3) = p(x_1 = 1, x_2 = 1, y_1 = 1, y_2 = 3) + p(x_1 = 2, x_2 = 1, y_1 = 1, y_2 = 3)
            shared_x2_eq_1 = ((0.4 * 0.7 * 0.6 * 0.1) + (0.6 * 0.1 * 0.1 * 0.1))

            # p(x_2 = 2, y_1 = 1, y_2 = 3) = p(x_1 = 1, x_2 = 2, y_1 = 1, y_2 = 3) + p(x_1 = 2, x_2 = 2, y_1 = 1, y_2 = 3)
            shared_x2_eq_2 = ((0.4 * 0.3 * 0.6 * 0.5) + (0.6 * 0.9 * 0.1 * 0.5)) #

            f11 = 0.4 # p(x_1 = 1)
            f12 = 0.6 # p(x_1 = 2)

            # p(x_2 = 1, y_1 = 1) = p(x_1 = 1, x_2 = 1, y_1 = 1) + p(x_1 = 2, x_2 = 1, y_1 = 1)
            f21 = (0.4 * 0.7 * 0.6) + (0.6 * 0.1 * 0.1)

            # p(x_2 = 2, y_1 = 1) = p(x_1 = 1, x_2 = 2, y_1 = 1) + p(x_1 = 2, x_2 = 2, y_1 = 1)
            f22 = (0.4 * 0.3 * 0.6) + (0.6 * 0.9 * 0.1)

            # p(x_3 = 1, y_1 = 1, y_2 = 3) = p(x_2 = 1, x_3 = 1, y_1 = 1, y_2 = 3) + p(x_2 = 2, x_3 = 1, y_1 = 1, y_2 = 3)
            f31 = shared_x2_eq_1 * 0.7 + shared_x2_eq_2 * 0.1

            # p(x_3 = 2, y_1 = 1, y_2 = 3) = p(x_2 = 1, x_3 = 2, y_1 = 1, y_2 = 3) + p(x_2 = 2, x_3 = 2, y_1 = 1, y_2 = 3)
            f32 = shared_x2_eq_1 * 0.3 + shared_x2_eq_2 * 0.9

            [f11 f12; f21 f22; f31 f32]
        end
        @test isapprox(actual, expected)

        # check log space agrees
        @test isapprox(exp.(SMC.hmm_log_forward_pass(hmm, observations)), expected)
    end

    @testset "marginal likelihoods" begin
        actual = SMC.hmm_marginal_likelihood(hmm, observations)

        expected = (
              (0.4 * 0.7 * 0.7 * 0.6 * 0.1 * 0.3) + # x = (1, 1, 1)
              (0.4 * 0.7 * 0.3 * 0.6 * 0.1 * 0.4) + # x = (1, 1, 2)
              (0.4 * 0.3 * 0.1 * 0.6 * 0.5 * 0.3) + # x = (1, 2, 1)
              (0.4 * 0.3 * 0.9 * 0.6 * 0.5 * 0.4) + # x = (1, 2, 2)
              (0.6 * 0.1 * 0.7 * 0.1 * 0.1 * 0.3) + # x = (2, 1, 1)
              (0.6 * 0.1 * 0.3 * 0.1 * 0.1 * 0.4) + # x = (2, 1, 2)
              (0.6 * 0.9 * 0.1 * 0.1 * 0.5 * 0.3) + # x = (2, 2, 1)
              (0.6 * 0.9 * 0.9 * 0.1 * 0.5 * 0.4)) # x = (2, 2, 2)

        @test isapprox(actual, expected)

        # check that log space agrees
        @test isapprox(actual, exp(hmm_log_marginal_likelihood(hmm, observations)))
    end

    @testset "posterior sampling" begin
        particle = hmm_posterior_sample(hmm, observations)

        # smoke test: check the type and dimensions and other basic properties
        @test typeof(particle) == Array{Int,1}
        @test length(particle) == 3
        @test all(particle .>= 1)
        @test all(particle .<= 2)

        # TODO: test the distribution? Maybe next year :)

    end

    @testset "sanity checks" begin

        # sanity check on extreme observation
        srand(1)
        hmm = HiddenMarkovModel([0.5, 0.5],                     # prior
                                [0.5 0.5; 0.5 0.5],             # transition
                                [0.9999 0.0001; 0.0001 0.9999]) # observation
        obs = [1]
        num_samples = 100
        num_1 = 0
        for i = 1:num_samples
            particle = hmm_posterior_sample(hmm, obs)
            num_1 += (particle[1] == 1)
        end
        @test num_1 >= 95

        # sanity check on extreme prior and transition
        srand(1)
        hmm = HiddenMarkovModel([0.0001, 0.9999],                     # prior
                                [0.9999 0.0001; 0.0001 0.9999],             # transition
                                [0.5 0.5; 0.5 0.5]) # observation
        obs = [1, 2]
        num_samples = 100
        num_1 = 0
        for i = 1:num_samples
            particle = hmm_posterior_sample(hmm, obs)
            num_1 += (particle[2] == 1)
        end
        @test num_1 <= 5

        # a symmetric extreme distribution
        hmm = HiddenMarkovModel([0.5, 0.5],         # prior
                                [0.9999 0.0001; 0.0001 0.9999],  # transition
                                [0.9999 0.0001; 0.0001 0.9999]) # observation

        # switch happens at step 9 -> 10
        obs_a = [1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2, 2, 2]
        obs_b = [2, 2, 2, 2, 2, 2, 2, 2, 2, 1, 1, 1, 1, 1, 1, 1, 1]

        # test for symmetry
        @test isapprox(hmm_log_marginal_likelihood(hmm, obs_a),
                       hmm_log_marginal_likelihood(hmm, obs_b))

        # check samples
        srand(1)
        num_samples = 100

        num_start_1 = 0
        num_end_2 = 0
        for i = 1:num_samples
            particle = hmm_posterior_sample(hmm, obs_a)
            num_start_1 += all(particle[1:9] .== 1)
            num_end_2 += all(particle[10:end] .== 2)
        end
        @test num_start_1 >= 95
        @test num_end_2 >= 95

        srand(1)
        num_start_2 = 0
        num_end_1 = 0
        for i = 1:num_samples
            particle = hmm_posterior_sample(hmm, obs_b)
            num_start_2 += all(particle[1:9] .== 2)
            num_end_1 += all(particle[10:end] .== 1)
        end
        @test num_start_2 >= 95
        @test num_end_1 >= 95
    end

end

@testset "smc tests" begin

    hmm = HiddenMarkovModel(
            [0.4, 0.6], # prior
            [0.7 0.3;
             0.1 0.9], # transition
            [0.6 0.3 0.1;
             0.1 0.4 0.5]) # observation
    observations = [1, 3, 2]

    # SMC with prior proposal
    srand(1)
    num_particles = 10000
    scheme = HMMPriorSMCScheme(hmm, observations, num_particles)
    output, log_ml_estimate = smc(scheme)
    @test length(output) == 3
    expected = hmm_log_marginal_likelihood(hmm, observations)
    @test isapprox(log_ml_estimate, expected, atol=0.1, rtol=0)

    # conditional SMC (take output from above as input)
    log_ml_estimate = conditional_smc(scheme, output)
    @test isapprox(log_ml_estimate, expected, atol=0.1, rtol=0)

    # SMC with conditional proposal
    srand(1)
    num_particles = 10000
    scheme = HMMConditionalSMCScheme(hmm, observations, num_particles)
    output, log_ml_estimate = smc(scheme)
    @test length(output) == 3
    expected = hmm_log_marginal_likelihood(hmm, observations)
    @test isapprox(log_ml_estimate, expected, atol=0.02, rtol=0)

    # conditional SMC (take output from above as input)
    log_ml_estimate = conditional_smc(scheme, output)
    @test isapprox(log_ml_estimate, expected, atol=0.1, rtol=0)
end

@testset "smc sandwiching the evidence with exact sampling" begin

    hmm = HiddenMarkovModel(
            [0.4, 0.6], # prior
            [0.7 0.3;
             0.1 0.9], # transition
            [0.6 0.3 0.1;
             0.1 0.4 0.5]) # observation
    observations = [1, 3, 2]

    # true marginal likelihood
    expected = hmm_log_marginal_likelihood(hmm, observations)

    function test_bounds(smc_scheme_constructor, all_num_particles::Array{Int,1},
                         num_reps::Int)
        lower_bound_raw_data = Array{Float64, 2}(length(all_num_particles), num_reps)
        upper_bound_raw_data = Array{Float64, 2}(length(all_num_particles), num_reps)
        for (i, num_particles) in enumerate(all_num_particles)
            scheme = HMMConditionalSMCScheme(hmm, observations, num_particles)
            for j = 1:num_reps
                output, lml_estimate = smc(scheme)
                lower_bound_raw_data[i, j] = lml_estimate
            end
            for j = 1:num_reps
                output = hmm_posterior_sample(hmm, observations)
                lml_estimate = conditional_smc(scheme, output)
                upper_bound_raw_data[i, j] = lml_estimate
            end
        end

        lower_bound_means = mean(lower_bound_raw_data, 2)
        lower_bound_stds = std(lower_bound_raw_data, 2)
        upper_bound_means = mean(upper_bound_raw_data, 2)
        upper_bound_stds = std(upper_bound_raw_data, 2)

        # they should be upper and lower bounds in all cases
        for i=1:length(all_num_particles)
            @test lower_bound_means[i] <= expected
            @test upper_bound_means[i] >= expected
        end

        # the estimates should get tighter and less noisy with more particles
        for i=1:length(all_num_particles)-1
            @test lower_bound_means[i] <= lower_bound_means[i+1]
            @test upper_bound_means[i] >= upper_bound_means[i+1]
            @test lower_bound_stds[i] >= lower_bound_stds[i+1]
            @test upper_bound_stds[i] >= upper_bound_stds[i+1]
        end

        # the final estimates should both be tight
        @test isapprox(lower_bound_means[end], expected, atol=0.01, rtol=0)
        @test isapprox(upper_bound_means[end], expected, atol=0.01, rtol=0)

    end

    # SMC with prior proposal
    srand(1)
    test_bounds(HMMConditionalSMCScheme, [1, 3, 10, 30, 100], 1000)

    # SMC with conditional proposal
    srand(1)
    test_bounds(HMMPriorSMCScheme, [1, 3, 10, 30, 100], 1000)

end

@testset "rendering hmm" begin
    hmm = HiddenMarkovModel(
        [0.4, 0.6], # prior
        [0.7 0.3;
         0.1 0.9], # transition
        [0.6 0.3 0.1;
         0.1 0.4 0.5]) # observation

    hmm_fname = "$(tempname()).png"
    latent_states_fname = "$(tempname()).png"
    observations_fname = "$(tempname()).png"

    # just check that rendering works in the current environment without crashign
    plt[:figure]()
    render_hmm!(hmm)
    plt[:savefig](hmm_fname)

    plt[:figure]()
    render_hmm_states!(hmm, [1, 1, 1, 2, 2, 1, 2, 1, 2, 2, 2], true, [1, 10], [1, 2])
    plt[:savefig](latent_states_fname)

    plt[:figure]()
    render_hmm_observations!(hmm, [1, 1, 1, 3, 3, 1, 2, 1, 2, 2, 2])
    plt[:savefig](observations_fname)

    Base.Filesystem.rm(hmm_fname)
    Base.Filesystem.rm(latent_states_fname)
    Base.Filesystem.rm(observations_fname)
end


text:
using Libdl

const gdip_version = "0.1.0"
download("https://github.com/comrob/gdip/archive/refs/tags/v$(gdip_version).tar.gz", "v$(gdip_version).tar.gz")
run(`tar xvf v$(gdip_version).tar.gz`)

pwd_dir = pwd()

cd("gdip-$(gdip_version)/gdip")
run(`bash install.sh`)
run(`cp lib/libGDIP.$(Libdl.dlext) $(pwd_dir)`)


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

using Markdown
using InteractiveUtils

# ╔═╡ 6fd616f6-1ea1-11eb-3814-8bfb4a096c49
using Pkg, DrWatson

# ╔═╡ 6ffe2628-1ea1-11eb-24ea-57f985146a72
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Plots, LaTeXStrings
end;

# ╔═╡ ed174bc4-1ea0-11eb-1e2f-a32874cec549
md"## Listing 3.5"

# ╔═╡ 1bec75ec-1f9d-11eb-1296-a5f2914d82de
f2(x) = (x<0 ? x+1 : 1-x)*(abs(x)<1 ? 1 : 0)

# ╔═╡ 7bf7e1a6-1f9d-11eb-02e3-fb81fb7fd472
begin
	a, b = -1.5, 1.5
	delta = 0.01
	xGrid = a:delta:b
	F(x) = sum([f2(u)*delta for u in a:delta:x])
	y = [F(u) for u in xGrid]
end;

# ╔═╡ 7c03a482-1f9d-11eb-3313-fd537ac60b9b
plot(xGrid, y, c=:blue, xlims=(a,b), ylims=(0,1),
	xlabel=L"x", ylabel=L"F(x)", legend=:none)

# ╔═╡ 70182cd8-1ea1-11eb-094d-8d0d49cf15f3
md"## End of listing 3.5"

# ╔═╡ Cell order:
# ╟─ed174bc4-1ea0-11eb-1e2f-a32874cec549
# ╠═6fd616f6-1ea1-11eb-3814-8bfb4a096c49
# ╠═6ffe2628-1ea1-11eb-24ea-57f985146a72
# ╠═1bec75ec-1f9d-11eb-1296-a5f2914d82de
# ╠═7bf7e1a6-1f9d-11eb-02e3-fb81fb7fd472
# ╠═7c03a482-1f9d-11eb-3313-fd537ac60b9b
# ╟─70182cd8-1ea1-11eb-094d-8d0d49cf15f3