text:
module Quadmath
using Requires

export Float128, ComplexF128, Inf128

import Base: (*), +, -, /,  <, <=, ==, ^, convert,
          reinterpret, sign_mask, exponent_mask, exponent_one, exponent_half,
          significand_mask, exponent, significand,
          promote_rule, widen,
          string, print, show, parse,
          acos, acosh, asin, asinh, atan, atanh, cosh, cos, sincos,
          exp, expm1, log, log2, log10, log1p, sin, sinh, sqrt,
          tan, tanh,
          ceil, floor, trunc, round, fma, decompose,
          copysign, flipsign, max, min, hypot, abs,
          ldexp, frexp, modf, nextfloat, typemax, typemin, eps,
          isinf, isnan, isfinite, isinteger,
          floatmin, floatmax, precision, signbit, maxintfloat,
          Int32, Int64, Float64, BigFloat, BigInt

using Random

if Sys.isapple()
    const quadoplib = "libquadmath.0"
    const libquadmath = "libquadmath.0"
elseif Sys.isunix()
    const quadoplib = "libgcc_s.so.1"
    const libquadmath = "libquadmath.so.0"
elseif Sys.iswindows()
    if Sys.WORD_SIZE == 64
        const quadoplib = "libgcc_s_seh-1.dll"
    else
        const quadoplib = "libgcc_s_sjlj-1.dll"
    end
    const libquadmath = "libquadmath-0.dll"
end

macro ccall(expr)
    @assert expr isa Expr && expr.head == :(::)
    ret_type = expr.args[2]

    expr_call = expr.args[1]
    @assert expr_call isa Expr && expr_call.head == :call

    expr_fname = expr_call.args[1]

    if expr_fname isa Symbol
        fname = QuoteNode(expr_fname)
    elseif expr_fname isa Expr && expr_fname.head == :.
        fname = :(($(expr_fname.args[2]), $(esc(expr_fname.args[1]))))
    end

    expr_args = expr_call.args[2:end]

    @assert all(ex isa Expr && ex.head == :(::) for ex in expr_args)

    arg_names = [ex.args[1] for ex in expr_args]
    arg_types = [ex.args[2] for ex in expr_args]

    if Sys.isunix()
        :(ccall($fname, $(esc(ret_type)), ($(esc.(arg_types)...),), $(esc.(arg_names)...)))
    else
        if ret_type == :Cfloat128
            quote
                r = Ref{Cfloat128}()
                ccall($fname, Cvoid, (Ref{Cfloat128}, $(esc.(arg_types)...),), r, $(esc.(arg_names)...))
                r[]
            end
        else
            :(ccall($fname, $(esc(ret_type)), ($(esc.(arg_types)...),), $(esc.(arg_names)...)))
        end
    end
end

# we use this slightly cumbersome definition to ensure that the value is 128-bit aligned
# and passed on the xmm registers, matching the x86_64 ABI for __float128.
const Cfloat128 = NTuple{2,VecElement{Float64}}

struct Float128 <: AbstractFloat
    data::Cfloat128
end
convert(::Type{Float128}, x::Number) = Float128(x)

const ComplexF128 = Complex{Float128}

Base.cconvert(::Type{Cfloat128}, x::Float128) = x.data
Base.cconvert(::Type{Ref{Cfloat128}}, x::Float128) = Ref{Cfloat128}(x.data)

# reinterpret
function reinterpret(::Type{UInt128}, x::Float128)
    hi = reinterpret(UInt64, x.data[2].value)
    lo = reinterpret(UInt64, x.data[1].value)
    UInt128(hi) << 64 | lo
end
function reinterpret(::Type{Float128}, x::UInt128)
    fhi = reinterpret(Float64, (x >> 64) % UInt64)
    flo = reinterpret(Float64, x % UInt64)
    Float128((VecElement(flo), VecElement(fhi)))
end
reinterpret(::Type{Unsigned}, x::Float128) = reinterpret(UInt128, x)
reinterpret(::Type{Signed}, x::Float128) = reinterpret(Int128, x)

reinterpret(::Type{Int128}, x::Float128) =
    reinterpret(Int128, reinterpret(UInt128, x))
reinterpret(::Type{Float128}, x::Int128) =
    reinterpret(Float128, reinterpret(UInt128, x))

function __init__()
    @require SpecialFunctions="276daf66-3868-5448-9aa4-cd146d93841b" begin
        include("specfun.jl")
    end
end

sign_mask(::Type{Float128}) =        0x8000_0000_0000_0000_0000_0000_0000_0000
exponent_mask(::Type{Float128}) =    0x7fff_0000_0000_0000_0000_0000_0000_0000
exponent_one(::Type{Float128}) =     0x3fff_0000_0000_0000_0000_0000_0000_0000
exponent_half(::Type{Float128}) =    0x3ffe_0000_0000_0000_0000_0000_0000_0000
significand_mask(::Type{Float128}) = 0x0000_ffff_ffff_ffff_ffff_ffff_ffff_ffff

fpinttype(::Type{Float128}) = UInt128

# conversion
Float128(x::Float128) = x

# Float64
Float128(x::Float64) =
    Float128(@ccall(quadoplib.__extenddftf2(x::Cdouble)::Cfloat128))
Float64(x::Float128) =
    @ccall(quadoplib.__trunctfdf2(x::Cfloat128)::Cdouble)

# Float32
Float128(x::Float32) =
    Float128(@ccall(quadoplib.__extendsftf2(x::Cfloat)::Cfloat128))
Float32(x::Float128) =
    @ccall(quadoplib.__trunctfsf2(x::Cfloat128)::Cfloat)

# Float16
Float128(x::Float16) = Float128(Float32(x))
Float16(x::Float128) = Float16(Float32(x)) # TODO: avoid double rounding

# integer -> Float128
Float128(x::Int32) =
    Float128(@ccall(quadoplib.__floatsitf(x::Int32)::Cfloat128))

Float128(x::UInt32) =
    Float128(@ccall(quadoplib.__floatunsitf(x::UInt32)::Cfloat128))

Float128(x::Int64) =
    Float128(@ccall(quadoplib.__floatditf(x::Int64)::Cfloat128))

Float128(x::UInt64) =
    Float128(@ccall(quadoplib.__floatunditf(x::UInt64)::Cfloat128))

Float128(x::Int16) = Float128(Int32(x))
Float128(x::Int8) = Float128(Int32(x))
Float128(x::UInt16) = Float128(UInt32(x))
Float128(x::UInt8) = Float128(UInt32(x))

function Float128(x::UInt128)
    x == 0 && return Float128(0.0)
    n = 128-leading_zeros(x) # ndigits0z(x,2)
    if n <= 113
        y = ((x % UInt128) << (113-n)) & significand_mask(Float128)
    else
        y = ((x >> (n-114)) % UInt128) & 0x0001_ffff_ffff_ffff_ffff_ffff_ffff_ffff # keep 1 extra bit
        y = (y+1)>>1 # round, ties up (extra leading bit in case of next exponent)
        y &= ~UInt128(trailing_zeros(x) == (n-114)) # fix last bit to round to even
    end
    d = ((n+16382) % UInt128) << 112
    # reinterpret(Float128, d + y)
    d += y
    if Sys.iswindows()
        return reinterpret(Float128,d)
    else
        y1 = reinterpret(Float64,UInt64(d >> 64))
        y2 = reinterpret(Float64,(d % UInt64))
        return Float128((y2,y1))
    end
end

function Float128(x::Int128)
    x == 0 && return 0.0
    s = reinterpret(UInt128,x) & sign_mask(Float128) # sign bit
    x = abs(x) % UInt128
    n = 128-leading_zeros(x) # ndigits0z(x,2)
    if n <= 113
        y = ((x % UInt128) << (113-n)) & significand_mask(Float128)
    else
        y = ((x >> (n-114)) % UInt128) & 0x0001_ffff_ffff_ffff_ffff_ffff_ffff_ffff # keep 1 extra bit
        y = (y+1)>>1 # round, ties up (extra leading bit in case of next exponent)
        y &= ~UInt128(trailing_zeros(x) == (n-114)) # fix last bit to round to even
    end
    d = ((n+16382) % UInt128) << 112
    # reinterpret(Float128, s | d + y)
    d = s | d + y
    if Sys.iswindows()
        return reinterpret(Float128,d)
    else
        y1 = reinterpret(Float64,UInt64(d >> 64))
        y2 = reinterpret(Float64,(d % UInt64))
        Float128((y2,y1))
    end
end

# Float128 -> integer requires arithmetic, so is below

# Rational
Float128(x::Rational{T}) where T = Float128(numerator(x))/Float128(denominator(x))

Float128(x::Bool) = x ? Float128(1) : Float128(0)

# Comparison
(==)(x::Float128, y::Float128) =
    @ccall(quadoplib.__eqtf2(x::Cfloat128, y::Cfloat128)::Cint) == 0
(<)(x::Float128, y::Float128) =
    @ccall(quadoplib.__letf2(x::Cfloat128, y::Cfloat128)::Cint) == -1
(<=)(x::Float128, y::Float128) =
    @ccall(quadoplib.__letf2(x::Cfloat128, y::Cfloat128)::Cint) <= 0

# Arithmetic
(+)(x::Float128, y::Float128) =
    Float128(@ccall(quadoplib.__addtf3(x::Cfloat128, y::Cfloat128)::Cfloat128))
(-)(x::Float128, y::Float128) =
    Float128(@ccall(quadoplib.__subtf3(x::Cfloat128, y::Cfloat128)::Cfloat128))
(*)(x::Float128, y::Float128) =
    Float128(@ccall(quadoplib.__multf3(x::Cfloat128, y::Cfloat128)::Cfloat128))
(/)(x::Float128, y::Float128) =
    Float128(@ccall(quadoplib.__divtf3(x::Cfloat128, y::Cfloat128)::Cfloat128))

(-)(x::Float128) =
    Float128(@ccall(quadoplib.__negtf2(x::Cfloat128)::Cfloat128))

# Float128 -> Integer
unsafe_trunc(::Type{Int32}, x::Float128) =
    @ccall(quadoplib.__fixtfsi(x::Cfloat128)::Int32)

unsafe_trunc(::Type{Int64}, x::Float128) =
    @ccall(quadoplib.__fixtfdi(x::Cfloat128)::Int64)

unsafe_trunc(::Type{UInt32}, x::Float128) =
    @ccall(quadoplib.__fixunstfsi(x::Cfloat128)::UInt32)

unsafe_trunc(::Type{UInt64}, x::Float128) =
    @ccall(quadoplib.__fixunstfdi(x::Cfloat128)::UInt64)

function unsafe_trunc(::Type{UInt128}, x::Float128)
    xu = reinterpret(UInt128,x)
    k = (Int64(xu >> 112) & 0x07fff) - 16382 - 113
    xu = (xu & significand_mask(Float128)) | 0x0001_0000_0000_0000_0000_0000_0000_0000
    if k <= 0
        UInt128(xu >> -k)
    else
        UInt128(xu) << k
    end
end
function unsafe_trunc(::Type{Int128}, x::Float128)
    copysign(unsafe_trunc(UInt128,x) % Int128, x)
end
trunc(::Type{Signed}, x::Float128) = trunc(Int,x)
trunc(::Type{Unsigned}, x::Float128) = trunc(Int,x)
trunc(::Type{Integer}, x::Float128) = trunc(Int,x)

for Ti in (Int32, Int64, Int128, UInt32, UInt64, UInt128)
    if Ti <: Unsigned || sizeof(Ti) < sizeof(Float128)
        # Here `Float128(typemin(Ti))-1` is exact, so we can compare the lower-bound
        # directly. `Float128(typemax(Ti))+1` is either always exactly representable, or
        # rounded to `Inf` (e.g. when `Ti==UInt128 && Float128==Float32`).
        @eval begin
            function trunc(::Type{$Ti},x::Float128)
                if Float128(typemin($Ti)) - one(Float128) < x < Float128(typemax($Ti)) + one(Float128)
                    return unsafe_trunc($Ti,x)
                else
                    throw(InexactError(:trunc, $Ti, x))
                end
            end
            function (::Type{$Ti})(x::Float128)
                if (Float128(typemin($Ti)) <= x <= Float128(typemax($Ti))) && (round(x, RoundToZero) == x)
                    return unsafe_trunc($Ti,x)
                else
                    throw(InexactError($(Expr(:quote,Ti.name.name)), $Ti, x))
                end
            end
        end
    else
        # Here `eps(Float128(typemin(Ti))) > 1`, so the only value which can be truncated to
        # `Float128(typemin(Ti)` is itself. Similarly, `Float128(typemax(Ti))` is inexact and will
        # be rounded up. This assumes that `Float128(typemin(Ti)) > -Inf`, which is true for
        # these types, but not for `Float16` or larger integer types.
        @eval begin
            function trunc(::Type{$Ti},x::Float128)
                if Float128(typemin($Ti)) <= x < Float128(typemax($Ti))
                    return unsafe_trunc($Ti,x)
                else
                    throw(InexactError(:trunc, $Ti, x))
                end
            end
            function (::Type{$Ti})(x::Float128)
                if (Float128(typemin($Ti)) <= x < Float128(typemax($Ti))) && (round(x, RoundToZero) == x)
                    return unsafe_trunc($Ti,x)
                else
                    throw(InexactError($(Expr(:quote,Ti.name.name)), $Ti, x))
                end
            end
        end
    end
end

## math

## one argument
for f in (:acos, :acosh, :asin, :asinh, :atan, :atanh, :cosh, :cos,
          :exp, :expm1, :log, :log2, :log10, :log1p,
          :sin, :sinh, :sqrt, :tan, :tanh,
          :ceil, :floor, :trunc, )
    @eval function $f(x::Float128)
        Float128(@ccall(libquadmath.$(string(f,:q))(x::Cfloat128)::Cfloat128))
    end
end

abs(x::Float128) = Float128(@ccall(libquadmath.fabsq(x::Cfloat128)::Cfloat128))
round(x::Float128) = Float128(@ccall(libquadmath.rintq(x::Cfloat128)::Cfloat128))
round(x::Float128, r::RoundingMode{:Down}) = floor(x)
round(x::Float128, r::RoundingMode{:Up}) = ceil(x)
round(x::Float128, r::RoundingMode{:ToZero}) = round(x)

## two argument
(^)(x::Float128, y::Float128) =
    Float128(@ccall(libquadmath.powq(x::Cfloat128, y::Cfloat128)::Cfloat128))

# circumvent a failure in Base
function (^)(x::Float128, p::Integer)
    if p >= 0
        Base.power_by_squaring(x,p)
    else
        Base.power_by_squaring(inv(x),-p)
    end
end
copysign(x::Float128, y::Float128) =
    Float128(@ccall(libquadmath.copysignq(x::Cfloat128, y::Cfloat128)::Cfloat128))
hypot(x::Float128, y::Float128) =
    Float128(@ccall(libquadmath.hypotq(x::Cfloat128, y::Cfloat128)::Cfloat128))
atan(x::Float128, y::Float128) =
    Float128(@ccall(libquadmath.atan2q(x::Cfloat128, y::Cfloat128)::Cfloat128))
sincos(x::Float128) = (sin(x), cos(x))

## misc
@static if !Sys.iswindows()
    # disable fma on Windows until rounding mode issue fixed
    # https://github.com/JuliaMath/Quadmath.jl/issues/31
    fma(x::Float128, y::Float128, z::Float128) =
        Float128(@ccall(libquadmath.fmaq(x::Cfloat128, y::Cfloat128, z::Cfloat128)::Cfloat128))
end

isnan(x::Float128) = 0 != @ccall(libquadmath.isnanq(x::Cfloat128)::Cint)
isinf(x::Float128) = 0 != @ccall(libquadmath.isinfq(x::Cfloat128)::Cint)
isfinite(x::Float128) = 0 != @ccall(libquadmath.finiteq(x::Cfloat128)::Cint)

isinteger(x::Float128) = isfinite(x) && x === trunc(x)

signbit(x::Float128) = signbit(reinterpret(Int128, x))

flipsign(x::Float128, y::Float128) = signbit(y) ? -x : x

precision(::Type{Float128}) = 113

eps(::Type{Float128}) = reinterpret(Float128, 0x3f8f_0000_0000_0000_0000_0000_0000_0000)
floatmin(::Type{Float128}) = reinterpret(Float128, 0x0001_0000_0000_0000_0000_0000_0000_0000)
floatmax(::Type{Float128}) = reinterpret(Float128, 0x7ffe_ffff_ffff_ffff_ffff_ffff_ffff_ffff)
maxintfloat(::Type{Float128}) = Float128(0x0002_0000_0000_0000_0000_0000_0000_0000)

"""
    Inf128

Positive infinity of type [`Float128`](@ref).
"""
const Inf128 = reinterpret(Float128, 0x7fff_0000_0000_0000_0000_0000_0000_0000)
typemax(::Type{Float128}) = Inf128
typemin(::Type{Float128}) = -Inf128

ldexp(x::Float128, n::Cint) =
    Float128(@ccall(libquadmath.ldexpq(x::Cfloat128, n::Cint)::Cfloat128))
ldexp(x::Float128, n::Integer) =
    ldexp(x, clamp(n, typemin(Cint), typemax(Cint)) % Cint)

function frexp(x::Float128)
    ri = Ref{Cint}()
    f = Float128(@ccall(libquadmath.frexpq(x::Cfloat128, ri::Ptr{Cint})::Cfloat128))
    return f, Int(ri[])
end

function modf(x::Float128)
    isinf(x) && return (zero(Float128), x)
    ipart = trunc(x)
    fpart = x - ipart
    return fpart, ipart
end

significand(x::Float128) = frexp(x)[1] * 2
function exponent(x::Float128)
     !isfinite(x) && throw(DomainError("Cannot be NaN or Inf."))
     abs(x) > 0 && return frexp(x)[2] - 1
     throw(DomainError("Cannot be subnormal converted to 0."))
end

function nextfloat(f::Float128, d::Integer)
    F = typeof(f)
    fumax = reinterpret(Unsigned, F(Inf))
    U = typeof(fumax)

    isnan(f) && return f
    fi = reinterpret(Signed, f)
    fneg = fi < 0
    fu = unsigned(fi & typemax(fi))

    dneg = d < 0
    da = Base.uabs(d)
    if da > typemax(U)
        fneg = dneg
        fu = fumax
    else
        du = da % U
        if fneg ⊻ dneg
            if du > fu
                fu = min(fumax, du - fu)
                fneg = !fneg
            else
                fu = fu - du
            end
        else
            if fumax - fu < du
                fu = fumax
            else
                fu = fu + du
            end
        end
    end
    if fneg
        fu |= sign_mask(F)
    end
    reinterpret(F, fu)
end

Float128(::Irrational{:π}) =  reinterpret(Float128, 0x4000921fb54442d18469898cc51701b8)
Float128(::Irrational{:ℯ}) =  reinterpret(Float128, 0x40005bf0a8b1457695355fb8ac404e7a)
Float128(x::Irrational{T}) where {T} = Float128(BigFloat(x))

import Base.MPFR

function BigFloat(x::Float128; precision=precision(BigFloat))
    if !isfinite(x) || iszero(x)
        @static if VERSION < v"1.1"
            return BigFloat(Float64(x), precision)
        else
            return BigFloat(Float64(x), precision=precision)
        end
    end

    @static if VERSION < v"1.1"
        b = setprecision(BigFloat, max(precision,113)) do
            BigFloat()
        end
    else
        b = BigFloat(precision=max(precision,113))
    end

    y, k = frexp(x)
    b.exp = Clong(k)
    b.sign = signbit(x) ? Cint(-1) : Cint(1)
    u = (reinterpret(UInt128, y) << 15) | 0x8000_0000_0000_0000_0000_0000_0000_0000
    i = cld(precision, sizeof(MPFR.Limb)*8)
    while u != 0
        w = (u >> (128-sizeof(MPFR.Limb)*8)) % MPFR.Limb
        unsafe_store!(b.d, w, i)
        i -= 1
        u <<= sizeof(MPFR.Limb)*8
    end
    # set remaining bits to zero
    while i > 0
        unsafe_store!(b.d, zero(MPFR.Limb), i)
        i -= 1
    end

    if precision < 113
        @static if VERSION < v"1.1"
            b2 = setprecision(BigFloat, precision) do
                BigFloat()
            end
            ccall((:mpfr_set, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32),
                  b2, b, MPFR.ROUNDING_MODE[])
            return b2
        else
            return BigFloat(b, precision=precision)
        end
    else
        return b
    end
end

function Float128(x::BigFloat)
    if !isfinite(x) || iszero(x)
        return Float128(Float64(x))
    end

    y,k = frexp(x)
    if k >= -16381
        prec = 113
    elseif k >= -16381-112
        prec = 113 + (k + 16381)
    elseif k == -16381-113 && abs(y) > 0.5
        z = reinterpret(Float128, UInt128(1))
    else
        z = reinterpret(Float128, UInt128(0))
    end

    @static if VERSION < v"1.1"
        y = setprecision(BigFloat, prec) do
            BigFloat()
        end
        ccall((:mpfr_set, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Int32),
                  y, x, MPFR.ROUNDING_MODE[])
    else
        y = BigFloat(x, precision=prec)
    end

    u = zero(UInt128)
    i = cld(prec, sizeof(MPFR.Limb)*8)
    j = 113
    while i > 0
        j -= sizeof(MPFR.Limb)*8
        u |= (unsafe_load(y.d, i) % UInt128) << j
        i -= 1
    end
    u &= significand_mask(Float128)
    u |= exponent_half(Float128)
    z = ldexp(reinterpret(Float128, u), y.exp)
    return copysign(z,x)
end

function BigInt(x::Float128)
    !isinteger(x) && throw(InexactError(BigInt, x))
    BigInt(BigFloat(x, precision=precision(Float128)))
end
function Float128(x::BigInt)
    @static if VERSION < v"1.1"
        y = setprecision(BigFloat, precision(Float128)) do
            BigFloat(x)
        end
    else
        y = BigFloat(x, precision=precision(Float128))
    end
    Float128(y)
end

promote_rule(::Type{Float128}, ::Type{Float16}) = Float128
promote_rule(::Type{Float128}, ::Type{Float32}) = Float128
promote_rule(::Type{Float128}, ::Type{Float64}) = Float128
promote_rule(::Type{Float128}, ::Type{<:Integer}) = Float128

#widen(::Type{Float64}) = Float128
widen(::Type{Float128}) = BigFloat

function decompose(x::Float128)::Tuple{Int128, Int, Int}
    isnan(x) && return 0, 0, 0
    isinf(x) && return ifelse(x < 0, -1, 1), 0, 0
    n = reinterpret(UInt128, x)
    s = (n & significand_mask(Float128)) % Int128
    e = ((n & exponent_mask(Float128)) >> 112) % Int
    s |= Int128(e != 0) << 112
    d = ifelse(signbit(x), -1, 1)
    s, e - 16495 + (e == 0), d
end

function Random.rand(rng::AbstractRNG, s::Random.SamplerTrivial{Random.CloseOpen01{Float128}})
    u = rand(rng, UInt128)
    x = (reinterpret(Float128, u & Base.significand_mask(Float128)
                     | Base.exponent_one(Float128))
         - one(Float128))
    return x
end

# TODO: need to do this better
function parse(::Type{Float128}, s::AbstractString)
    Float128(@ccall(libquadmath.strtoflt128(s::Cstring, C_NULL::Ptr{Ptr{Cchar}})::Cfloat128))
end

function string(x::Float128)
    lng = 64
    buf = Array{UInt8}(undef, lng + 1)
    lng = @ccall(libquadmath.quadmath_snprintf(buf::Ptr{UInt8}, (lng+1)::Csize_t, "%.35Qe"::Ptr{UInt8}, x::(Cfloat128...))::Cint)
    return String(resize!(buf, lng))
end

print(io::IO, b::Float128) = print(io, string(b))
show(io::IO, b::Float128) = print(io, string(b))

include("printf.jl")

end # module Quadmath


text:
using Test
using ProbabilisticCircuits
using DataFrames
using LogicCircuits

@testset "BMC tests with SamplePSDD" begin
    # Set up a logic constraint ϕ as a BDD and scope size n.
    function case(ϕ::Bdd, n::Integer; atol::Real = 0)
        # All possible valuations (including impossible ones).
        M = all_valuations(collect(1:n))
        # Get only possible worlds.
        W = M[findall(ϕ.(eachrow(M))),:]
        # Assign random probabilities for each world in W.
        R = rand(1:20, size(W, 1))
        # Construct a dataset that maps the distribution of R (world W[i] repeats R[i] times).
        D = DataFrame(vcat([repeat(W[i,:], 1, R[i])' for i ∈ 1:size(W, 1)]...), :auto)
        # Learn ensemble of PSDDs from ϕ and D.
        B = bmc_sample_psdd(5, ϕ, 16, D, 3, 4; pseudocount = 0.0, maxiter = 100, verbose = false)
        T = DataFrame(M, :auto)
        # Test consistency.
        @test (EVI(B, T) .> -Inf) == ϕ.(eachrow(M))
        # Test probabilities.
        evi = exp.(EVI(B, T))
        @test isapprox(evi[findall(>(0), evi)], (R/sum(R)); atol)
    end

    # some tests commented out to speed up tests
    case((1 ∧ 2) ∨ (3 ∧ ¬4) ∨ (¬1 ∧ 5), 5)
    # case((1 → 3) ∧ (5 → ¬2), 5)
    # case((1 ∧ 2 ∧ 3) ∨ (4 ∧ 5), 5)
    case(exactly(3, collect(1:5)), 5)
    # case(atleast(3, collect(1:5)), 5)
    # case(atmost(3, collect(1:5)), 5)
end

ext: jl


text:
using NearestNeighborDescent: sample_neighbors, make_knn_heaps, deheap_knns

@testset "deheap_knns tests" begin
    true_ids = [1 2 3;
                2 3 4]
    true_dists = [0.1 0.2 0.3;
                  0.2 0.3 0.4]
    heap1 = BinaryMaxHeap{NNTuple{Int, Float64}}()
    push!(heap1, NNTuple(1, 0.1))
    push!(heap1, NNTuple(2, 0.2))
    push!(heap1, NNTuple(3, 0.3))
    heap2 = BinaryMaxHeap{NNTuple{Int, Float64}}()
    push!(heap2, NNTuple(2, 0.2))
    push!(heap2, NNTuple(3, 0.3))
    push!(heap2, NNTuple(4, 0.4))
    heap3 = BinaryMaxHeap{NNTuple{Int, Float64}}()
    push!(heap3, NNTuple(3, 0.3))
    push!(heap3, NNTuple(4, 0.4))
    push!(heap3, NNTuple(5, 0.5))
    @test top(heap1) == NNTuple(3, 0.3)
    @test top(heap2) == NNTuple(4, 0.4)
    @test top(heap3) == NNTuple(5, 0.5)
    heaps = [heap1, heap2, heap3]
    ids, dists = @inferred deheap_knns(heaps, 2)
    @test ids == true_ids
    @test dists == true_dists
end

@testset "sample neighbors tests" begin
    @testset "sample_rate = 1. tests" begin
        points = collect(1:10)
        n_neighbors = 4

        # zero neighbors
        idxs = sample_neighbors(length(points), 0)
        @test length(idxs) == 0

        # k < n
        idxs = sample_neighbors(length(points), n_neighbors)
        @test length(idxs) == n_neighbors
        @test issubset(idxs, points)

        # k > n
        idxs = sample_neighbors(length(points), length(points)+5)
        @test length(idxs) == length(points)
    end
    @testset "sample_rate = .5 tests" begin
        points = collect(1:10)
        n_neighbors = 4
        ρ = .5

        # zero neighbors
        idxs = sample_neighbors(length(points), 0, ρ)
        @test length(idxs) == 0

        # k < n
        idxs = sample_neighbors(length(points), n_neighbors, ρ)
        @test ρ*n_neighbors ≥ length(idxs)
        @test issubset(idxs, points)

        # k > n
        idxs = sample_neighbors(length(points), 2*length(points), ρ)
        @test length(idxs) == length(points)
    end
    @testset "exclude set tests" begin
        points = collect(1:10)

        # exclude 1
        idxs = sample_neighbors(length(points),
                                          length(points),
                                          exclude=[1])
        @test idxs ⊊ points
        @test !(1 ∈ idxs)

        # exclude all
        idxs = sample_neighbors(length(points),
                                          length(points),
                                          exclude=points)
        @test length(idxs) == 0

    end
end

@testset "make_knn_heaps tests" begin
    @testset "Float64 tests" begin
        data = [rand(3) for _ in 1:10]
        n_neighbors = 3

        knn_heaps = make_knn_heaps(data, n_neighbors, Euclidean())

        @test length(knn_heaps) == length(data)
        for p in 1:length(knn_heaps)
            @test length(knn_heaps[p]) == n_neighbors
            for t = 1:length(knn_heaps[p])
                @test knn_heaps[p].valtree[t].idx ≠ p
                @test knn_heaps[p].valtree[t].idx ≤ length(data)
                d = evaluate(Euclidean(), data[p], data[knn_heaps[p].valtree[t].idx])
                @test knn_heaps[p].valtree[t].dist == d
            end
        end
    end
    @testset "Int tests" begin
        data = [rand([0, 1], 3) for _ in 1:10]
        n_neighbors = 2
        knn_heaps = make_knn_heaps(data, n_neighbors, Hamming())

        @test length(knn_heaps) == length(data)
        for p in 1:length(knn_heaps)
            @test length(knn_heaps[p]) == n_neighbors
            for t = 1:length(knn_heaps[p])
                @test knn_heaps[p].valtree[t].idx ≠ p
                @test knn_heaps[p].valtree[t].idx ≤ length(data)
                d = evaluate(Hamming(), data[p], data[knn_heaps[p].valtree[t].idx])
                @test knn_heaps[p].valtree[t].dist == d
            end
        end

    end
end

@testset "neighbors tests" begin
    # create heaps
    knn_heaps = [BinaryMaxHeap{NNTuple{Int, Float64}}() for _ in 1:5]
    push!(knn_heaps[1], NNTuple(2, 100.))
    push!(knn_heaps[1], NNTuple(5, 100.))
    push!(knn_heaps[2], NNTuple(4, 100.))
    push!(knn_heaps[3], NNTuple(2, 100.))
    push!(knn_heaps[3], NNTuple(4, 100.))
    push!(knn_heaps[4], NNTuple(1, 100.))
    push!(knn_heaps[5], NNTuple(3, 100.))
    push!(knn_heaps[5], NNTuple(4, 100.))
    old_fw, fw, old_bw, bw = _neighbors(knn_heaps)

    @testset "new fw neighbors tests" begin
        @test fw[1] == [2, 5]
        @test fw[2] == [4]
        @test fw[3] == [2, 4]
        @test fw[4] == [1]
        @test fw[5] == [3, 4]
        @test all(isempty(v) for v in old_fw)
    end
    @testset "new bw neighbors tests" begin
        @test bw[1] == [4]
        @test bw[2] == [1, 3]
        @test bw[3] == [5]
        @test bw[4] == [2, 3, 5]
        @test bw[5] == [1]
        @test all(isempty(v) for v in old_bw)
    end
    old_fw, fw, old_bw, bw = _neighbors(knn_heaps)
    @testset "old fw neighbors tests" begin
        @test old_fw[1] == [2, 5]
        @test old_fw[2] == [4]
        @test old_fw[3] == [2, 4]
        @test old_fw[4] == [1]
        @test old_fw[5] == [3, 4]
        @test all(isempty(v) for v in fw)
    end
    @testset "old bw neighbors tests" begin
        @test old_bw[1] == [4]
        @test old_bw[2] == [1, 3]
        @test old_bw[3] == [5]
        @test old_bw[4] == [2, 3, 5]
        @test old_bw[5] == [1]
        @test all(isempty(v) for v in bw)
    end
end

@testset "_heappush! tests" begin
    @testset "max_cand tests" begin
        h = BinaryMaxHeap{NNTuple{Int, Float64}}()
        t = NNTuple(1, 1.)
        _heappush!(h, t, 0)
        @test length(h) == 0

        _heappush!(h, t, 1)
        @test length(h) == 1
        @test top(h) == t

        d = NNTuple(2, .5)
        _heappush!(h, d, 1)
        @test length(h) == 1
        @test top(h) == d
    end

    @testset "return val tests" begin
        h = BinaryMaxHeap{NNTuple{Int, Float64}}()
        # max_cand
        @test _heappush!(h, NNTuple(1, rand()), 0) == 0
        # empty heap push
        @test _heappush!(h, NNTuple(1, 1.), 1) == 1
        # length == max AND further away, no push
        @test _heappush!(h, NNTuple(2, 2.), 1) == 0
        # length == max BUT closer, push
        @test _heappush!(h, NNTuple(3, .5), 1) == 1
        @test top(h).idx == 3
        @test top(h).dist == .5
        @test length(h) == 1
        # length < max AND further, push
        @test _heappush!(h, NNTuple(4, 4.), 2) == 1
        @test top(h).idx == 4
        @test top(h).dist == 4.
        # tuple already in heap, no push
        @test _heappush!(h, NNTuple(3, .5), 3) == 0
        @test length(h) == 2
        @test top(h).idx == 4
        @test top(h).dist == 4.
    end
end


text:
mutable struct Normal_Record <: GMC_NS_Model_Record
    trajectory::Int64
    i::Int64
    pos::Vector{Float64}
    path::String
    log_Li::Float64
end

struct Normal_Model <: GMC_NS_Model
    trajectory::Integer #id integer assigned at construction
    i::Integer #id of parent model, if any

    θ::Vector{Float64} #model's parameter Vector
    log_Li::Float64 #model's likelihood, calculated by constructor

    pos::Vector{Float64}
    v::Vector{Float64} #model's velocity in parameter space
end

function construct_normal_model(trajectory::Integer, i::Integer, θ::Vector{Float64}, pos::Vector{Float64}, v::Vector{Float64}, obs::Vector{Float64}; v_init=false)
    μ,λ=θ
    mod_normal=Normal(μ,sqrt(1/λ))
    log_lh=lps(logpdf.(mod_normal, obs))
    v_init && (v=rand(MvNormal(length(θ),1.)))

    Normal_Model(trajectory, i, θ, log_lh, pos, v)
end

function Base.show(io::IO, m::Normal_Model; xsteps=100)
    μ,λ=m.θ
    σ=sqrt(1/λ)
    n=Normal(μ,σ)
    X=[μ-4σ:((μ+4σ)-(μ-4σ))/xsteps:μ+4σ...]
    y=pdf.(n,X)

    show(io, lineplot(X,y,xlabel="X",ylabel="p",title="Normal Model $(m.trajectory).$(m.i), log_Li $(m.log_Li)"))
    println()
    println("θ: $(m.θ)")
    println("v: $(m.v)")
end


text:
module test_julia_global

using Test

n = 1
@test n == 1

function f()
    global n += 2
end

f()
@test n == 3

@test test_julia_global.n == 3

end # module test_julia_global

module test_julia_global_local

using Test

local n = 1
@test_throws UndefVarError n == 1

function f()
    global n
    n += 2
end

@test_throws UndefVarError f()
@test_throws UndefVarError n == 3

end # module test_julia_global_local


text:
"""
    CentralReconstruction{K}

A simple, central, polynomial reconstruction (linear procedure) without adaptive
techniques using polynomials of degee `K-1`.
"""
struct CentralReconstruction{K} <: AbstractReconstruction
    function CentralReconstruction{K}() where {K}
        @argcheck (K>0 && isodd(K)) ArgumentError("`K` must be odd (K = $K).")
        new{K}()
    end
end

function CentralReconstruction(::Val{K}=Val{3}()) where {K}
    CentralReconstruction{K}()
end

function CentralReconstruction(K::Integer)
    CentralReconstruction{K}()
end

@inline stencil_width(::CentralReconstruction{K}) where {K} = K
@inline stencil_width_val(::CentralReconstruction{K}) where {K} = Val{K}()

@inline order(::CentralReconstruction{K}) where {K} = K

function (recons::CentralReconstruction{3})(edge_u, cell, u_m1, u_0, u_p1, balance_law, meshx)
    edge_u[1,cell] = ( 5*u_0 - u_p1 + 2*u_m1 ) / 6
    edge_u[2,cell] = ( 5*u_0 + 2*u_p1 - u_m1 ) / 6

    nothing
end

function (recons::CentralReconstruction{5})(edge_u, cell, u_m2, u_m1, u_0, u_p1, u_p2, balance_law, meshx)
    edge_u[1,cell] = ( 47*u_0 - 13*u_p1 + 2*u_p2 + 27*u_m1 - 3*u_m2 ) / 60
    edge_u[2,cell] = ( 47*u_0 + 27*u_p1 - 3*u_p2 - 13*u_m1 + 2*u_m2 ) / 60

    nothing
end

function (recons::CentralReconstruction{7})(edge_u, cell, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, balance_law, meshx)
    edge_u[1,cell] = ( 319*u_0 - 101*u_p1 + 25*u_p2 - 3*u_p3 + 214*u_m1 - 38*u_m2 + 4*u_m3 ) / 420
    edge_u[2,cell] = ( 319*u_0 + 214*u_p1 - 38*u_p2 + 4*u_p3 - 101*u_m1 + 25*u_m2 - 3*u_m3 ) / 420

    nothing
end

function interpolate!(uval, ξ, u_m1, u_0, u_p1, balance_law, mesh, recons::CentralReconstruction{3})
    for i in eachindex(uval)
        uval[i] = @evalpoly(ξ[i], 5*u_0/6 - u_p1/6 + u_m1/3, u_0 - u_m1, -u_0 + u_p1/2 + u_m1/2)
    end

    nothing
end

function interpolate!(uval, ξ, u_m2, u_m1, u_0, u_p1, u_p2, balance_law, mesh, recons::CentralReconstruction{5})
    for i in eachindex(uval)
        uval[i] = @evalpoly(ξ[i], 47*u_0/60 - 13*u_p1/60 + u_p2/30 + 9*u_m1/20 - u_m2/20, 5*u_0/4 - u_p1/12 - 5*u_m1/4 + u_m2/12, -u_0 + 3*u_p1/4 - u_p2/8 + u_m1/4 + u_m2/8, -u_0/2 + u_p1/6 + u_m1/2 - u_m2/6, u_0/4 - u_p1/6 + u_p2/24 - u_m1/6 + u_m2/24)
    end

    nothing
end

function interpolate!(uval, ξ, u_m3, u_m2, u_m1, u_0, u_p1, u_p2, u_p3, balance_law, mesh, recons::CentralReconstruction{7})
    for i in eachindex(uval)
        uval[i] = @evalpoly(ξ[i], 319*u_0/420 - 101*u_p1/420 + 5*u_p2/84 - u_p3/140 + 107*u_m1/210 - 19*u_m2/210 + u_m3/105, 49*u_0/36 - 5*u_p1/36 + u_p2/90 - 49*u_m1/36 + 5*u_m2/36 - u_m3/90, -23*u_0/24 + 7*u_p1/8 - 19*u_p2/80 + 7*u_p3/240 + u_m1/16 + 21*u_m2/80 - u_m3/30, -7*u_0/9 + 11*u_p1/36 - u_p2/36 + 7*u_m1/9 - 11*u_m2/36 + u_m3/36, 23*u_0/72 - 13*u_p1/48 + 5*u_p2/48 - u_p3/72 - u_m1/6 + u_m2/48 + u_m3/144, u_0/12 - u_p1/24 + u_p2/120 - u_m1/12 + u_m2/24 - u_m3/120, -u_0/36 + u_p1/48 - u_p2/120 + u_p3/720 + u_m1/48 - u_m2/120 + u_m3/720)
    end

    nothing
end


text:
module _DispEn2D
export DispEn2D
using Clustering: kmeans, assignments
using Statistics: std, mean
using StatsFuns: normcdf
    """
        Disp2D, RDE = DispEn2D(Mat)

    Returns the bidimensional dispersion entropy estimate (`Disp2D`) and reverse
    bidimensional dispersion entropy (`RDE`) estimated for the data matrix (`Mat`)
    using the default parameters:  time delay = 1, symbols = 3, logarithm = natural,
    data transform = normalised cumulative density function (`'ncdf'`), logarithm = natural,
    template matrix size = [floor(H/10) floor(W/10)], (where H and W represent
    the height (rows) and width (columns) of the data matrix `Mat`) \n
    ** The minimum number of rows and columns of Mat must be > 10.**

        Disp2D, RDE = DispEn2D(Mat::AbstractArray{T,2} where T<:Real;
                            m::Union{Int,Tuple{Int,Int}}=floor.(Int, size(Mat)./10), tau::Int=1,
                                c::Int=3, Typex::String="ncdf", Logx::Real=exp(1), Norm::Bool=false, Lock::Bool=true)

    Returns the bidimensional dispersion entropy (`Disp2D`) and reverse
    bidimensional distribution entropy (`RDE`) estimate for the data matrix (`Mat`)
    using the specified 'keyword' arguments:

    # Arguments:
    `m`     - Template submatrix dimensions, an integer scalar (i.e. the same
              height and width) or a two-element tuple of integers
              [height, width] with a minimum value > 1.
              [default: [floor(H/10) floor(W/10)]] \n
    `tau`   - Time Delay, a positive integer       [default: 1]  \n
    `c`     - Number of symbols, an integer > 1
    `Typex` - Type of symbolic mapping transform, one of the following:
              {`linear`, `kmeans`, `ncdf`, `equal`}
              See the `EntropyHub Guide` for more info on these transforms.
    `Logx`  - Logarithm base, a positive scalar        [default: natural]\n
              ** enter 0 for natural logarithm.**\n
    `Norm`  - Normalisation of `Disp2D` value, a boolean:
                - [false]   no normalisation - default
                - [true]    normalises w.r.t number of possible dispersion patterns.
    `Lock`  - By default, DispEn2D only permits matrices with a maximum
              size of 128 x 128 to prevent memory errors when storing data on RAM.
              e.g. For Mat = [200 x 200], m = 3, and tau = 1, DispEn2D
              creates a vector of 753049836 elements. To enable matrices
              greater than [128 x 128] elements, set `Lock` to false.
              [default: 'true']
              `WARNING: unlocking the permitted matrix size may cause your Julia
               IDE to crash.`

    # See also `DispEn`, `DistEn2D`, `SampEn2D`, `FuzzEn2D`, `MSEn`

    # References:
        [1] Hamed Azami, et al.,
            "Two-dimensional dispersion entropy: An information-theoretic
            method for irregularity analysis of images."
            Signal Processing: Image Communication,
            75 (2019): 178-187.

    """
    function DispEn2D(Mat::AbstractArray{T,2} where T<:Real; m::Union{Int,Tuple{Int,Int}}=floor.(Int, size(Mat)./10), tau::Int=1,
                        c::Int=3, Typex::String="ncdf", Logx::Real=exp(1), Norm::Bool=false, Lock::Bool=true)

    Logx == 0  ? Logx = exp(1) : nothing

    NL, NW = size(Mat)
    ((NL > 128 || NW > 128) && Lock) ?
        error("To prevent memory errors, matrix width & length must have <= 128 elements.
            To estimatDispEn2D  for the current matrix ($NL,$NW) change Lock to 'false'.
            Caution: unlocking the safe matrix size may cause the Julia IDE to crash.") :
            nothing

    length(m)==1 ?  (mL = m; mW = m) : (mL = m[1];  mW = m[2])
    (NL > 10 && NW > 10) ? nothing :
    error("Number of rows and columns in Mat must be > 10")
    (minimum(m) > 1) ? nothing : error("m:   must be an integer > 1, or 2-element integer tuple w/ values > 1")
    (tau > 0) ? nothing : error("tau:   must be an integer > 0")
    (c > 1) ? nothing :   error("c:     must be an integer > 1")
    (lowercase(Typex) in ["linear", "kmeans", "ncdf", "equal"]) ? nothing :
        error("Typex:    must be one of the following strings - 'linear','kmeans','ncdf','equal'")
    (Logx>0) ? nothing :  error("Logx:  must be a positive number > 0")

    if lowercase(Typex) == "linear"
        Edges = range(minimum(Mat),maximum(Mat),length=c+1)
        Zi = map(x -> sum(Edges[1:c].<=x), Mat)

    elseif lowercase(Typex) == "kmeans"
        Temp = kmeans(transpose(Mat[:]), c; maxiter=200)
        Zx = assignments(Temp)
        Clux = Temp.centers
        xx = sortperm(Clux[:]);   Zi = zeros(Int,length(Mat))
        for k = 1:c
            Zi[Zx.==xx[k]] .= k;
        end
        Zi = reshape(Zi,size(Mat))

    elseif lowercase(Typex) == "ncdf"
        Zx = normcdf.(mean(Mat),std(Mat,corrected=false),Mat);
        Zi = map(x -> sum(range(0,1,length=c+1)[1:c].<=x), Zx);

    elseif lowercase(Typex) == "equal"
        ix = sortperm(Mat[:]);
        xx = Int.(round.(range(0,length(Mat),length=c+1)))
        Zi = zeros(Int, length(Mat))
        for k = 1:c
            Zi[ix[xx[k]+1:xx[k+1]]] .= k
        end
        Zi = reshape(Zi,size(Mat))
    end

    NL = NL - (mL-1)*tau
    NW = NW - (mW-1)*tau
    X = zeros(Int,NL*NW,mL*mW)
    p = 0
    for k = 1:NL
        for n = 1:NW
            p += 1
            X[p,:] = Zi[k:tau:(mL-1)*tau+k,n:tau:(mW-1)*tau+n][:]
        end
    end

    p != NL*NW ? @warn("Potential error with submatrix division.") : nothing

    T = unique(X,dims=1)
    Nx = size(T)[1]
    Counter = zeros(Nx)
    for n = 1:Nx
        Counter[n] = sum(all(X .- transpose(T[n,:]) .==0, dims=2))
    end

    Ppi = Counter[Counter.!= 0]/size(X)[1]

    big(c)^(mL*mW) > 10^16 ? error("RDE cannot be estimated with c = $c and
    a submatrix of size $mL x $mW. Required floating point precision exceeds 10^16.
    Consider reducing the template submatrix size (m) or the number of symbols (c).") : nothing

    #RDE = sum(Ppi.^2) - (1/Nx)
    RDE = sum((Ppi .- (1/(c^(mL*mW)))).^2);

    if round(sum(Ppi),digits=4) != 1
        @warn("Potential Error calculating probabilities")
    end

    Disp2D = -sum(Ppi.*log.(Logx, Ppi))
    if Norm
        Disp2D = Disp2D/log(Logx, c^(mL*mW))
        RDE = RDE./(1 - (1/(c^(mL*mW))))
    end

    return Disp2D, RDE
    end

end

"""
Copyright 2021 Matthew W. Flood, EntropyHub

Licensed under the Apache License, Version 2.0 (the "License");
you may not use this file except in compliance with the License.
You may obtain a copy of the License at

    http://www.apache.org/licenses/LICENSE-2.0

Unless required by applicable law or agreed to in writing, software
distributed under the License is distributed on an "AS IS" BASIS,
WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
See the License for the specific language governing permissions and
limitations under the License.

For Terms of Use see https://github.com/MattWillFlood/EntropyHub
"""


text:
using ITensors, Test

@testset "sweepnext function" begin
  @testset "one site" begin
    N = 6
    count = 1
    output = [
      (1, 1),
      (2, 1),
      (3, 1),
      (4, 1),
      (5, 1),
      (6, 1),
      (6, 2),
      (5, 2),
      (4, 2),
      (3, 2),
      (2, 2),
      (1, 2),
    ]
    for (b, ha) in sweepnext(N; ncenter=1)
      @test (b, ha) == output[count]
      count += 1
    end
    @test count == 2 * N + 1
  end

  @testset "two site" begin
    N = 6
    count = 1
    output = [
      (1, 1), (2, 1), (3, 1), (4, 1), (5, 1), (5, 2), (4, 2), (3, 2), (2, 2), (1, 2)
    ]
    for (b, ha) in sweepnext(N)
      @test (b, ha) == output[count]
      count += 1
    end
    @test count == 2 * (N - 1) + 1
  end

  @testset "three site" begin
    N = 6
    count = 1
    output = [(1, 1), (2, 1), (3, 1), (4, 1), (4, 2), (3, 2), (2, 2), (1, 2)]
    for (b, ha) in sweepnext(N; ncenter=3)
      @test (b, ha) == output[count]
      count += 1
    end
    @test count == 2 * (N - 2) + 1
  end
end

nothing


text:
using SparseArrays, JavaCall

# DBtable contains the table binding information, as well as the
# d4mQuery Java object for the table.
struct DBtable
    DB::DBserver
    name::String
    security::String
    numLimit::Int
    numRow::Int
    columnfamily::String
    putBytes::Number
    d4mQuery
    tableOps
end

# DBtablePair contains the table pair binding information, as well as the
# d4mQuery Java object for the table.
# The DBtablePair contains bindings to both a table and its transpose.
struct DBtablePair
    DB::DBserver
    name1::String
    name2::String
    security::String
    numLimit::Int
    numRow::Int
    columnfamily::String
    putBytes::Number
    d4mQuery
    tableOps
end

# Deletes specified DB table
# TODO: add check whether want to delete
function delete(table::DBtable)
    name = table.name
    DB = table.DB
    if  ispresent(DB, name)
        println("Deleting " * name * " in " * DB.instanceName)
        jcall(table.tableOps, "deleteTable", Nothing, (JString,), name)
    else
        println("Table (name " * name * ") already does not exist in " * DB.instanceName)
    end
    return nothing
end

function delete(table::DBtablePair)
    name1 = table.name1
    name2 = table.name2
    DB = table.DB

    if  ispresent(DB, name1)
        println("Deleting " * name1 * " in " * DB.instanceName)
        jcall(table.tableOps, "deleteTable", Nothing, (JString,), name1)
    else
        println("Table 1 (name " * name1 * ") already does not exist in " * DB.instanceName)
    end
    if  ispresent(DB, name2)
        println("Deleting " * name2 * " in " * DB.instanceName)
        jcall(table.tableOps, "deleteTable", Nothing, (JString,), name2)
    else
        println("Table 2 (name " * name2 * ") already does not exist in " * DB.instanceName)
    end

    return nothing
end

 # Helper function to convert Assoc fields to an Accumulo-friendly string
function toDBstring(input)
    StringArray = Array{T} where T <: AbstractString
    NumericArray = Array{T} where T <: Number
    UnionArray = Array{T} where T <: Union{AbstractString,Number}

    if isa(input, StringArray)
        output = join(input, "\n") * "\n"
    elseif isa(input, NumericArray)
        output = join(string.(input), "\n") * "\n"
    elseif isa(input, UnionArray)
        output = join(convert(Array{AbstractString}, input), "\n") * "\n"
    elseif isa(input, Colon)
        output = ":"
    elseif isa(input, StartsWith)
        # StartsWith supports one element, with a delimiter at the end,
        # as well as a delimited list with a delimiter at the end

        str = input.inputString

        output = ""
        del = str[end:end]
        idx1 = 1
        idx2 = findnext(del, str, idx1 + 1)

        while ~isa(idx2, Nothing)
            output = output * str[idx1:idx2[1]] * ":" * del * str[idx1:idx2[1] - 1] * Char(127) * del
            idx1 = idx2[1] + 1
            idx2 = findnext(del, str, idx1)
        end
    else
        output = input
    end

    return output
end

# Base getindex function- i and j are d4m formatted strings (delimitered)
DBtableType = Union{DBtable,DBtablePair}
function getindex(table::DBtableType, i::AbstractString, j::AbstractString)

    jcall(table.d4mQuery, "setCloudType", Nothing, (JString,), table.DB.dbType)
    jcall(table.d4mQuery, "setLimit", Nothing, (jint,), table.numLimit)
    jcall(table.d4mQuery, "reset", Nothing, (), )

    dbResultSet = @jimport "edu.mit.ll.d4m.db.cloud.D4mDbResultSet"

    if i != ":" || j == ":" || isa(table, DBtable)

        if isa(table, DBtablePair)
            jcall(table.d4mQuery, "setTableName", Nothing, (JString,), table.name1)
        end
        jcall(table.d4mQuery, "doMatlabQuery", dbResultSet, (JString, JString, JString, JString,), i, j, table.columnfamily, table.security)

        r = jcall(table.d4mQuery, "getRowReturnString", JString, (), )
        c = jcall(table.d4mQuery, "getColumnReturnString", JString, (), )

    else # Search transpose table if column query
        jcall(table.d4mQuery, "setTableName", Nothing, (JString,), table.name2)
        jcall(table.d4mQuery, "doMatlabQuery", dbResultSet, (JString, JString, JString, JString,), j, i, table.columnfamily, table.security)

        c = jcall(table.d4mQuery, "getRowReturnString", JString, (), )
        r = jcall(table.d4mQuery, "getColumnReturnString", JString, (), )
    end

    v = jcall(table.d4mQuery, "getValueReturnString", JString, (), )

    return deepCondense(Assoc(r, c, v))

end

UnionArray = Array{Union{AbstractString,Number}}
ValidQueryTypes = Union{Colon,AbstractString,Array,UnionArray,StartsWith}

getindex(table::DBtableType, i::ValidQueryTypes, j::ValidQueryTypes) = getindex(table, toDBstring(i), toDBstring(j))

# Query iterator functionality
function getiterator(table::DBtableType, nelements::Number)

    ops = @jimport "edu.mit.ll.d4m.db.cloud.D4mDbTableOperations"
    opsObj = ops((JString, JString, JString, JString,), table.DB.instanceName, table.DB.host, table.DB.user, table.DB.pass)

    if isa(table, DBtablePair)
        d4mQuery = @jimport "edu.mit.ll.d4m.db.cloud.D4mDataSearch"
        queryObj = d4mQuery((JString, JString, JString, JString, JString,), table.DB.instanceName, table.DB.host, table.name1, table.DB.user, table.DB.pass)

        Ti = DBtablePair(table.DB, table.name1, table.name2, table.security,
        nelements, table.numRow, table.columnfamily, table.putBytes,
        queryObj, opsObj)
    else
        d4mQuery = @jimport "edu.mit.ll.d4m.db.cloud.D4mDataSearch"
        queryObj = d4mQuery((JString, JString, JString, JString, JString,), table.DB.instanceName, table.DB.host, table.name, table.DB.user, table.DB.pass)

        Ti = DBtable(table.DB, table.name, table.security,
        nelements, table.numRow, table.columnfamily, table.putBytes,
        queryObj, opsObj)
    end

    return Ti
end

# The getindex function to get the next few elements for an iterator table object
function getindex(table::DBtableType)

    #T.d4mQuery.next();
    #dbResultSet = @jimport "edu.mit.ll.d4m.db.cloud.D4mDbResultSet"
    jcall(table.d4mQuery, "next", Nothing, (), )

    # check if transpose table query- if so flip the r and c results.
    tablename = jcall(table.d4mQuery, "getTableName", JString, (), )

    if isa(table, DBtablePair) && tablename == table.name2
        c = jcall(table.d4mQuery, "getRowReturnString", JString, (), )
        r = jcall(table.d4mQuery, "getColumnReturnString", JString, (), )
    else
        r = jcall(table.d4mQuery, "getRowReturnString", JString, (), )
        c = jcall(table.d4mQuery, "getColumnReturnString", JString, (), )
    end

    v = jcall(table.d4mQuery, "getValueReturnString", JString, (), )

    return deepCondense(Assoc(r, c, v))

end

# Helper function for ingest
function searchall(str::String, c::Char)
    idx = search(str, c)
    allIdx = idx

    while idx < endof(str)
        idx = search(str, c, idx + 1)
        allIdx = [allIdx idx]
    end

    return allIdx
end

# putTriple is the main db ingest function
function putTriple(table::DBtableType, r::UnionArray, c::UnionArray, v::UnionArray)

    # Find chunk size for ingest
    chunkBytes = table.putBytes
    numTriples = length(r);
    avgBytePerTriple = (sum(length.(r)) + sum(length.(c)) + sum(length.(v))) / numTriples;
    chunkSize = min(max(1, round(Integer, chunkBytes / avgBytePerTriple)), numTriples);

    # In Matlab D4M this is inside the loop- do we need to make a new object each time?
    dbInsert = @jimport "edu.mit.ll.d4m.db.cloud.D4mDbInsert"

    if isa(table, DBtablePair)
        insertObj = dbInsert((JString, JString, JString, JString, JString,), table.DB.instanceName, table.DB.host, table.name1, table.DB.user, table.DB.pass)
        insertObjT = dbInsert((JString, JString, JString, JString, JString,), table.DB.instanceName, table.DB.host, table.name2, table.DB.user, table.DB.pass)
    else
        insertObj = dbInsert((JString, JString, JString, JString, JString,), table.DB.instanceName, table.DB.host, table.name, table.DB.user, table.DB.pass)
    end

    # Insert each chunk
    for i in 1:chunkSize:numTriples
        iNext = min(i + chunkSize, numTriples);
        rr = toDBstring(r[i:iNext]);
        cc = toDBstring(c[i:iNext]);
        vv = toDBstring(v[i:iNext]);

        jcall(insertObj, "doProcessing", Nothing, (JString, JString, JString, JString, JString,), rr, cc, vv, table.columnfamily, table.security)

        # Insert transpose into transpose table if DBtablePair
        if isa(table, DBtablePair)
            jcall(insertObjT, "doProcessing", Nothing, (JString, JString, JString, JString, JString,), cc, rr, vv, table.columnfamily, table.security)
        end
    end

end

function ingestprep(input)
    StringArray = Array{T} where T <: AbstractString
    NumericArray = Array{T} where T <: Number
    UnionArray = Array{Union{AbstractString,Number}}

    if isa(input, StringArray)
        output = convert(UnionArray, input)
    elseif isa(input, NumericArray)
        output = string.(input)
    elseif isa(input, AbstractString)
        output = split(input, input[end], limit = NumStr(input))
    else
        output = input
    end

    return output
end

UnionArray = Array{Union{AbstractString,Number}}
#StringOrNumArray = Union{AbstractString,Array,Number}
ValidType = Union{UnionArray,AbstractString,Array,Number}

putTriple(table::DBtableType, r::ValidType, c::ValidType, v::ValidType) = putTriple(table, ingestprep(r), ingestprep(c), ingestprep(map(x-> isinteger2(x) ? convert(Integer, x) : x, v)))
# ^ automatically convert any integral floats (of value) to actual integers
# Graphulo has some problems with floats

function isinteger2(x)
    if isa(x, AbstractString)
        return false
    elseif isinteger(x)
        return true
    end
    return false
end

# The put function deconstructs A and calls putTriple
function put(table::DBtableType, A::Assoc)
    r, c, v = find(A)
    putTriple(table, r, c, v)
end

#addColCombiner: Adds combiners to specific column names.
# Inputs: T = database table object
    # colNames = string list of column names. If ':,' (default) then sets all columns to combining.
    # combineType = combiner name "min", "max", or "sum";
    #      IF D4M_API_JAVA.jar is installed on the Accumulo instance, you can
    #      specify "min_decimal", "max_decimal", or "sum_decimal"
    #      to obtain the ability to combine on decimals
function addColCombiner(table::DBtable, colNames = ":,", combineType = "sum")
    jcall(table.tableOps, "designateCombiningColumns", Nothing, (JString, JString, JString, JString,), table.name, colNames, combineType, table.columnfamily)
end

# nnz returns the number of entries in specified table
function nnz(table::DBtableType)
    if isa(table, DBtablePair)
        tname = table.name1
    else
        tname = table.name
    end

    arrayList = @jimport "java.util.ArrayList"
    list = @jimport "java.util.List"

    tableList = arrayList((), )
    jcall(tableList, "add", jboolean, (JObject,), tname)

    jcall(table.tableOps, "getNumberOfEntries", jlong, (list,), tableList)
end

function addsplits(table::DBtable, splitstring)
    jcall(table.tableOps, "addSplits", Nothing, (JString, JString,), table.name, splitstring)
end

function addsplits(table::DBtablePair, splitstring, splitstringT)

    jcall(table.tableOps, "addSplits", Nothing, (JString, JString,), table.name1, splitstring)
    jcall(table.tableOps, "addSplits", Nothing, (JString, JString,), table.name2, splitstringT)

end

function getsplits(table::DBtableType)

    if isa(table, DBtablePair)
        splits = jcall(table.tableOps, "getSplitsString", JString, (JString,), table.name1)
        splitsT = jcall(table.tableOps, "getSplitsString", JString, (JString,), table.name2)

        return splits, splitsT
    else
        splits = jcall(table.tableOps, "getSplitsString", JString, (JString,), table.name)

        return splits
    end
end

DBtableTypeorString = Union{DBtableType, AbstractString}
function makedegreetable(A::DBtableTypeorString, Rname::AbstractString; countColumns = true::Bool, colq = ""::AbstractString)
    # /**
    # * Create a degree table from an existing table.
    # * @param table Name of original table.
    # * @param Degtable Name of degree table. Created if it does not exist.
    # *                 Use a combiner if you want to sum in the new degree entries into an existing table.
    # * @param countColumns True means degrees are the number of entries in each row.
    # *                     False means degrees are the sum or weights of entries in each row.
    # * @param colq The name of the degree column in Degtable. Default is "".

    # Returns a database table struct.

    # public long generateDegreeTable(String table, String Degtable, boolean countColumns, String colq)

    if isa(A, DBtablePair)
        Aname = A.name1
    elseif isa(A, DBtable)
        Aname = A.name
    else
        Aname = A
    end

    jcall(A.DB.Graphulo, "generateDegreeTable", jlong,
        (JString, JString, jboolean, JString,),
        Aname, Rname, countColumns, colq)

end


text:
"""
Helping functions
"""

# Converts vector of words to vector of indices
function indices(wordVect::Vector, vocab::Vector, unk::String="_unk_")
    function index(x, unk)
        idx = something(findfirst(isequal(x), vocab), 0)
        idx > 0 || return findfirst(isequal(unk), vocab)
        return idx
    end
    return broadcast(x -> index(x, unk), wordVect)
end

# Padding multiple sequences w r t the max size sequence
function pre_pad_sequences(sequences::Vector, pad::String="_pad_")
    max_len = maximum([length(x) for x in sequences])
    return [[fill(pad, max_len-length(sequence)); sequence] for sequence in sequences]
end

function post_pad_sequences(sequences::Vector, pad::String="_pad_")
    max_len = maximum([length(x) for x in sequences])
    return [[sequence; fill(pad, max_len-length(sequence))] for sequence in sequences]
end

# To initialize funciton for model LSTM weights
init_weights(extreme::AbstractFloat, dims...) = randn(Float32, dims...) .* sqrt(Float32(extreme))

# Generator, whenever it should be called two times since it gives X in first and y in second call
function generator(c::Channel, corpus::AbstractDocument; batchsize::Integer=64, bptt::Integer=70)
    X_total = post_pad_sequences(chunk(tokens(corpus), batchsize))
    n_batches = Int(floor(length(X_total[1])/bptt))
    put!(c, n_batches)
    for i=1:n_batches
        start = bptt*(i-1) + 1
        batch = [Flux.batch(X_total[k][j] for k=1:batchsize) for j=start:start+bptt]
        put!(c, batch[1:end-1])
        put!(c, batch[2:end])
    end
end

"""
    get_buckets(c::Corpus, bucketsize::Integer; order::Bool=true)

Simple Sequence-Bucketing

This function will return the groups of `Document`s with close sequence lengths from the given `Corpus`.
Use this function with `data_loader` function to get a `Channel`. Also, the order of sequences lengths
can be set using keyword argument `order`, which is, `true` for ascending order (default) and `false` for
descending order of lengths of sequences.

# Example:

julia> buckets = get_buckets(corpus, 32);

"""
function get_buckets(c::Corpus, labels::Vector, bucketsize::Integer; order::Bool=true)
    lengths = length.(tokens.(documents(c)))
    sorted_lens = order ? sortperm(lengths) : reverse(sortperm(lengths))
    c, labels = c[sorted_lens], labels[sorted_lens]
    buckets = []
    for i=1:bucketsize:length(c)
        (length(c)-i) < (bucketsize-1) && (push!(buckets, zip(c[i:end], labels[i:end]));continue)
        push!(buckets, zip(c[i:i+bucketsize-1], labels[i:i+bucketsize-1]))
    end
    return buckets
end

"""
    data_loader(buckets::AbstractArray, classes::Vector; pad_func::Function=pre_pad_sequences)
    data_loader(dataset::Corpus, labels::Vector, batchsize::Integer; pad_func::Function=pre_pad_sequences)

This funciton can be use to make `Channel` to load data for the training of `TextClassifier`.
Do preprocessing of the text before using this function. And get all the examples into a `Corpus` type
and all the corresponding labels into a `Vector`.

NOTE: Remember that the first call to the `Channel` will output an integer which specifies the number of batches the `Channel` can give.

It can be used in two ways:
 - With sequence bucketing: use `get_buckets` function first on the prepared dataset to get buckets
   and pass those buckets to this function. Now, make a Vector of all possible classes and pass this to the `data_loader`
   and keep this classes variable safe since the output units of the classifier will be corresponding to this `Vector`.
   If already a list of all the classes is available that can also be passed to this function.

   Remember, the output buckets from the `get_buckets` the length of the last bucket may not be equal to the batchsize specified,
   since the number of examples couldn't be division by batchsize specified. That can be passed to data_loader
   but it might cause problems later while training, so that batch can be removed if it is not a big loss of training data.
   Or pass number of examples that can be divided into batches of specified batchsize.

julia> classes = unique(labels)     # labels contains the labels for the examples
julia> buckets = get_buckets(data, labels, 32);

# To check whether the buckets are all of same size:
julia> print(length(buckets[end]))

julia> loader = data_loader(buckets, classes)

 - Without sequence bucketing: pass the `Corpus` and labels direclty to this funciton and specify batchsize.
   This will return a Channel, which will give a batch at every call with labels.

julia> loader = data_loader(data, labels, 32)

The padding can be controlled by setting the `pad_func` keyword argument to
either `pre_pad_sequences` or `post_pad_sequence`, former is for pre padding the sequence
and that latter is for post padding the sequences.
"""
function data_loader(buckets::AbstractArray, classes::Vector; pad_func::Function=pre_pad_sequences)
    shuffle!(buckets)
    Channel(csize=1) do docs
        n_batches = length(buckets)
        put!(docs, n_batches)
        for b in buckets
            shuffle!(b)
            X, Y = [], []
            for (cur_text, label) in b
                toks = tokenize(cur_text)
                push!(X, toks)
                y = Flux.onehotbatch([label], classes)
                push!(Y, y)
            end#for
            X = pad_func(X, "_pad_")
            batchsize = length(X)
            put!(docs, [Flux.batch(X[k][j] for k=1:batchsize) for j=1:length(X[1])])
            put!(docs, cat(Y..., dims=2))
        end #for
    end #channel
end

function data_loader(dataset::Corpus, labels::Vector, batchsize::Integer; pad_func::Function=pre_pad_sequences)
    classes = unique(labels)
    n_batches = Int(floor(length(dataset)/batchsize))
    Channel(csize=1) do loader
        put!(loader, n_batches)
        for i=1:n_batches
            X = tokens.(dataset[(i-1)*batchsize+1:i*batchsize])
            Y = Flux.onehotbatch(labels[(i-1)*batchsize+1:i*batchsize], classes)
            X = pad_func(X, "_pad_")
            put!(loader, [Flux.batch(X[k][j] for k=1:batchsize) for j=1:length(X[1])])
            put!(loader, Y)
        end
    end
end


text:
function main!(r)
    counters = move_to_corner!(r)
    println(counters)
    mark_each_corner!(r)
    nachalo!(r, counters)
end

function move_to_corner!(r)
    counters = []
    while !(isborder(r,Nord) && isborder(r,West))
        counter_N = 0
        while !isborder(r,Nord)
            move!(r, Nord)
            counter_N+=1
        end
        pushfirst!(counters, counter_N)
        counter_W = 0
        while !isborder(r,West)
            move!(r, West)
            counter_W+=1
        end
        pushfirst!(counters, counter_W)
    end
    return counters
end

function mark_each_corner!(r)
    for i in 1:4
        while !isborder(r, HorizonSide(4-i))
            move!(r, HorizonSide(4-i))
        end
        putmarker!(r)
    end
end

function nachalo!(r, counters)
    println(size(counters))
    for i in 1:size(counters, 1)
        if i%2==0
            for j in 0:counters[i]-1
                move!(r, Sud)
            end

        end
        if i%2==1
            for j in 0:counters[i]-1
                move!(r, Ost)
            end
        end
    end
end


text:
# adapted from https://github.com/znxlwm/tensorflow-MNIST-GAN-DCGAN

using ADCME
using PyCall
using PyPlot
using DelimitedFiles

tt = pyimport("tensorflow.examples.tutorials.mnist")

function lrelu(x, th=0.2)
    return tf.maximum(th * x, x)
end

function generator(x, isTrain=true)
    local o
    variable_scope("generator") do

        # 1st hidden layer
        conv1 = tf.layers.conv2d_transpose(x, 1024, [4, 4], strides=(1, 1), padding="valid")
        lrelu1 = lrelu(tf.layers.batch_normalization(conv1, training=isTrain), 0.2)

        # 2nd hidden layer
        conv2 = tf.layers.conv2d_transpose(lrelu1, 512, [4, 4], strides=(2, 2), padding="same")
        lrelu2 = lrelu(tf.layers.batch_normalization(conv2, training=isTrain), 0.2)

        # 3rd hidden layer
        conv3 = tf.layers.conv2d_transpose(lrelu2, 256, [4, 4], strides=(2, 2), padding="same")
        lrelu3 = lrelu(tf.layers.batch_normalization(conv3, training=isTrain), 0.2)

        # 4th hidden layer
        conv4 = tf.layers.conv2d_transpose(lrelu3, 128, [4, 4], strides=(2, 2), padding="same")
        lrelu4 = lrelu(tf.layers.batch_normalization(conv4, training=isTrain), 0.2)

        # output layer
        conv5 = tf.layers.conv2d_transpose(lrelu4, 1, [4, 4], strides=(2, 2), padding="same")
        o = tf.nn.tanh(conv5)
    end
    return o
end

function discriminator(x, isTrain=true)
    local o, conv5
    variable_scope("discriminator") do
        # 1st hidden layer
        conv1 = tf.layers.conv2d(x, 128, [4, 4], strides=(2, 2), padding="same")
        lrelu1 = lrelu(conv1, 0.2)

        # 2nd hidden layer
        conv2 = tf.layers.conv2d(lrelu1, 256, [4, 4], strides=(2, 2), padding="same")
        lrelu2 = lrelu(tf.layers.batch_normalization(conv2, training=isTrain), 0.2)

        # 3rd hidden layer
        conv3 = tf.layers.conv2d(lrelu2, 512, [4, 4], strides=(2, 2), padding="same")
        lrelu3 = lrelu(tf.layers.batch_normalization(conv3, training=isTrain), 0.2)

        # 4th hidden layer
        conv4 = tf.layers.conv2d(lrelu3, 1024, [4, 4], strides=(2, 2), padding="same")
        lrelu4 = lrelu(tf.layers.batch_normalization(conv4, training=isTrain), 0.2)

        # output layer
        conv5 = tf.layers.conv2d(lrelu4, 1, [4, 4], strides=(1, 1), padding="valid")
        o = tf.nn.sigmoid(conv5)

    end
    return o, conv5
end

let fixed_z_ = randn(25, 1, 1, 100)
    global show_result
    function show_result(filename)
        test_images = run(sess, G_z, feed_dict=Dict(z=> fixed_z_, isTrain=>false))

        size_figure_grid = 5
        fig, ax = subplots(size_figure_grid, size_figure_grid, figsize=(5, 5))
        for i = 1:5
            for j = 1:5
                ax[i, j].get_xaxis().set_visible(false)
                ax[i, j].get_yaxis().set_visible(false)
            end
        end

        for k = 1:size_figure_grid*size_figure_grid
            i = div(k-1, size_figure_grid)+1
            j = (k%size_figure_grid) + 1
            ax[i, j].cla()
            ax[i, j].imshow(reshape(test_images[k,:,:,:], 64, 64), cmap="gray")
        end

        label = "Epoch $filename"
        # fig.text(0.5, 0.04, label, ha="center")

        savefig(filename)
        close("all")
    end
end

# training parameters
batch_size = 100
lr = 0.0002
train_epoch = 20

if !isdir("MNIST_data")
    mkdir("MNIST_data")
    download("http://yann.lecun.com/exdb/mnist/train-images-idx3-ubyte.gz", "MNIST_data/train-images-idx3-ubyte.gz")
    download("http://yann.lecun.com/exdb/mnist/t10k-labels-idx1-ubyte.gz", "MNIST_data/t10k-labels-idx1-ubyte.gz")
    download("http://yann.lecun.com/exdb/mnist/t10k-images-idx3-ubyte.gz", "MNIST_data/t10k-images-idx3-ubyte.gz")
    download("http://yann.lecun.com/exdb/mnist/train-labels-idx1-ubyte.gz", "MNIST_data/train-labels-idx1-ubyte.gz")
end
if !isdir("figures")
    mkdir("figures")
end

mnist = tt.input_data.read_data_sets("MNIST_data/", one_hot=true, reshape=[])

x = placeholder(Float32, shape=(nothing, 64, 64, 1))
z = placeholder(Float32, shape=(nothing, 1, 1, 100))
isTrain = placeholder(Bool)

G_z = generator(z, isTrain)
# networks : discriminator
D_real, D_real_logits = discriminator(x, isTrain)
D_fake, D_fake_logits = discriminator(G_z, isTrain)

# loss for each network
D_loss_real = mean(tf.nn.sigmoid_cross_entropy_with_logits(logits=D_real_logits, labels=tf.ones([batch_size, 1, 1, 1])))
D_loss_fake = mean(tf.nn.sigmoid_cross_entropy_with_logits(logits=D_fake_logits, labels=tf.zeros([batch_size, 1, 1, 1])))
D_loss = D_loss_real + D_loss_fake
G_loss = mean(tf.nn.sigmoid_cross_entropy_with_logits(logits=D_fake_logits, labels=tf.ones([batch_size, 1, 1, 1])))

D_vars = [x for x in get_collection() if startswith(x.name, "discriminator")]
G_vars = [x for x in get_collection() if startswith(x.name, "generator")]

control_dependencies(get_collection(tf.GraphKeys.UPDATE_OPS)) do
    global D_optim = AdamOptimizer(lr, beta1=0.5).minimize(D_loss, var_list=D_vars)
    global G_optim = AdamOptimizer(lr, beta1=0.5).minimize(G_loss, var_list=G_vars)
end

sess = Session(); init(sess)
train_set = run(sess, tf.image.resize_images(mnist[1].images, [64, 64]))
train_set = (train_set .- 0.5) / 0.5
num_examples = size(train_set, 1)

loss_d = []
loss_g = []
for epoch = 1:train_epoch
    for iter = 1:div(num_examples, batch_size)
        x_ = train_set[(iter-1)*batch_size+1:iter*batch_size,:,:,:]
        z_ = randn(batch_size, 1, 1, 100)
        loss_d_, _ = run(sess, [D_loss, D_optim], feed_dict=Dict(x=>x_, z=>z_, isTrain=>true))

        z_ = randn(batch_size, 1, 1, 100)
        loss_g_, _ = run(sess, [G_loss, G_optim], feed_dict=Dict(x=>x_, z=>z_, isTrain=>true))

        println("[$epoch, #$iter] $loss_d_, $loss_g_")
        push!(loss_d, loss_d_)
        push!(loss_g, loss_g_)
        show_result("figures/$(epoch)_$(iter)")
    end

end

close("all")
plot(loss_d, label="Discriminator")
plot(loss_g, label="Generator")
xlabel("Iterations")
ylabel("Loss")
savefig("figures/loss.png")

writedlm("figures/lossd.txt", loss_d)
writedlm("figures/lossg.txt", loss_g)


text:
# Unit test for Campbell Scientific logger file reading
@testset "Campbell logger file reading" begin
	data,units = readcampbell(joinpath(dirname(@__DIR__),"test","input","campbell_data.dat"));
	@test size(data) == (5,5)
	@test size(units) == (4,)
	@test data[:datetime][1] == DateTime(2008,07,15,16,15,00)
	@test sum(data[:Val1]) == 0.5
	@test data[:Val3][end] == 0.3
	@test units[1] == "RN"
	@test units[2] == "unit1"
	@test units[end] == "unit3"
end