text:
using Test
using Distributions
using DataStructures

import GlobalSensitivityAnalysis: scale_samples!, ishigami, sample

##
## 1. utils
##

# SobolData
N = 100
calc_second_order = false
parameters = OrderedDict(
    :param1 => Normal(0,1),
    :param2 => LogNormal(5, 20),
    :param3 => TriangularDist(0, 4, 1)
)

data1 = SobolData(params = parameters, calc_second_order = calc_second_order, N = N)
@test data1.params[:param1] == parameters[:param1]
@test data1.params[:param2] == parameters[:param2]
@test data1.params[:param2] == parameters[:param2]
@test data1.params[:param3] == parameters[:param3]
@test data1.N == N
@test data1.calc_second_order == calc_second_order

data2 = SobolData()
@test data2.params == nothing
@test data2.N == 1000
@test data2.calc_second_order == true

data3 = SobolData(params = parameters)
data4 = SobolData(calc_second_order = false)
data5 = SobolData(N = 100)

##
## 2. Sample Sobol
##

samples = sample(data1)
@test size(samples, 2) == length(data1.params)
if data1.calc_second_order
    @test size(samples,1) == data1.N * (2 * length(data1.params) + 2)
else
    @test size(samples,1) == data1.N * (length(data1.params) + 2)
end

samples3 = sample(data3)
@test size(samples3, 2) == length(data3.params)
if data3.calc_second_order
    @test size(samples3,1) == data3.N * (2 * length(data3.params) + 2)
else
    @test size(samples3,1) == data3.N * (length(data3.params) + 2)
end

@test_throws ErrorException sample(data2) # params are nothing
@test_throws ErrorException sample(data4) # params are nothing
@test_throws ErrorException sample(data5) # params are nothing

##
## 3a. Analyze Sobol
##

Y1 = ishigami(samples)
results = analyze(data1, Y1)
for Si in results[:firstorder]
    @test Si <= 1
end
for CI in results[:firstorder_conf]
    @test CI > 0
end
for St in results[:totalorder]
    @test St <= 1
end
for CI in results[:totalorder_conf]
    @test CI > 0
end
@test sum(results[:totalorder]) > sum(results[:firstorder])

Y3 = ishigami(samples3)
results = analyze(data3, Y3)
for Si in results[:firstorder]
    @test Si <= 1
end
for CI in results[:firstorder_conf]
    @test ismissing(CI) || CI > 0
end
for S2i in results[:secondorder]
    @test ismissing(S2i) || S2i <= 1
end
for CI in results[:secondorder_conf]
    @test ismissing(CI) || CI > 0
end
for St in results[:totalorder]
    @test St <= 1
end
for CI in results[:totalorder_conf]
    @test ismissing(CI) || CI > 0
end
@test sum(results[:totalorder]) > sum(results[:firstorder])

##
## 3b. Analyze Sobol Optional Keyword Args
##

data = SobolData(
    params = OrderedDict(:x1 => Normal(1, 0.2),
        :x2 => Uniform(0.75, 1.25),
        :x3 => LogNormal(0, 0.5)),
    N = 1000
)
samples = sample(data)
Y = ishigami(samples)
results = analyze(data, Y)

results = analyze(data, Y; progress_meter = false) # no progress bar should show

@test length(analyze(data, Y; N_override = 10)) == 6
results_override = analyze(data, Y, N_override = data.N)
results_original = analyze(data, Y)
@test results_override[:firstorder] == results_original[:firstorder]
@test results_override[:totalorder] == results_original[:totalorder]
@test_throws ErrorException analyze(data1, Y1; N_override = data.N + 1) # N_override > N

max_stars_repo_name: lrennels/SAJulia

max_issues_repo_name: lrennels/SAJulia


text:
module NLIDatasets

using DataDeps
using HTTP

function getzip(url, file)
    headers = ["User-Agent" => "NLIDatasets.jl",
               "Accept" => "*/*",
               "Accept-Encoding" => "gzip, deflate, br"]
    HTTP.download(url, file, headers)
end

function register_data(name, descr, url, sha; fetch=getzip, postfetch=unpack)
    kw = (fetch_method=fetch, post_fetch_method=postfetch)
    DataDeps.register(DataDep(name, descr, url, sha; kw...))
end

include("snli.jl")
include("multinli.jl")
include("xnli.jl")
include("scitail.jl")
include("hans.jl")
include("breaking_nli.jl")
include("anli.jl")

using .SNLI, .MultiNLI, .XNLI, .SciTail, .HANS, .BreakingNLI, .ANLI

end # module

ext: jl


text:
include("solvers.jl")

function main()
    solvers = [
        solve011 solve012;
        solve021 solve022;
        solve031 solve032;
        solve041 solve042;
        solve051 solve052;
        solve061 solve062;
        solve071 solve072;
        solve081 solve082;
        solve091 solve092;
    ]

    I_EXTRA_ARG = 3
    if length(ARGS) < I_EXTRA_ARG
        println("ARGS[1]    Day 1-25.")
        println("ARGS[2]    Part 1 or 2.")
        println("ARGS[$I_EXTRA_ARG:]   Args passed to solver, meaning is solver dependent.")
        println("")
        exit()
    end

    day = parse(Int, ARGS[1])
    part = parse(Int, ARGS[2])

    if !(firstindex(solvers) <= day <= lastindex(solvers))
        println("Day $day not implemented yet.")
        exit()
    elseif !(1 <= part <= 2)
        println("Part must be 1 or 2.")
        exit()
    end

    println("ARGS: ", ARGS)
    solve = solvers[day, part]

    solver_args = length(ARGS) >= I_EXTRA_ARG ? ARGS[I_EXTRA_ARG:end] : ()
    println("Solver Args: ", solver_args)

    result = solve(solver_args...)
    result === nothing ? println("Solve found nothing.") : println(result)
end

main()


text:
using Base.Test
using CUDArt
using CUDNN
@show CUDNN_VERSION

import CUDNN.cudnnConvolutionForward
import CUDNN.cudnnConvolutionBackwardFilter
import CUDNN.cudnnConvolutionBackwardData
import CUDNN.cudnnGetConvolutionNdForwardOutputDim
import CUDNN.cudnnPoolingForward
import CUDNN.cudnnPoolingBackward
include("conv_pool_cpu.jl")

padding=0
stride=1
# cd =  CD(; ndims=2, padding=0, stride=1, upscale=1, mode=CUDNN_CONVOLUTION, eltype=Float32)
x = reshape(Float32[1:16;], 4, 4, 1, 1); tx = CudaArray(x); @show x
w = reshape(Float32[1:4;], 2, 2, 1, 1); tw = CudaArray(w); @show w
y = zeros(Float32,3,3,1,1); ty = CudaArray(y); @show y
cudnnConvolutionForward(tx,tw,ty; padding=padding, stride=stride); y = to_host(ty); @show y
y2 = zeros(Float32,3,3,1,1); @show y2
cudnnConvolutionForward(x,w,y2; padding=padding, stride=stride); @show y2
@test_approx_eq y y2

dy = rand(Float32, size(y)); tdy = CudaArray(dy); @show dy
dw = zeros(Float32, size(w)); tdw = CudaArray(dw); @show dw
cudnnConvolutionBackwardFilter(tx,tdy,tdw); dw = to_host(tdw); @show dw
dw2 = zeros(Float32, size(w));  @show dw2
cudnnConvolutionBackwardFilter(x,dy,dw2); @show dw2
@test_approx_eq dw dw2

dx = zeros(Float32, size(x)); tdx = CudaArray(dx); @show dx
cudnnConvolutionBackwardData(tw, tdy, tdx); dx = to_host(tdx); @show dx
dx2 = zeros(Float32, size(x));  @show dx2
cudnnConvolutionBackwardData(w, dy, dx2); @show dx2
@test_approx_eq dx dx2

#=
x = rand(Float32,64,64,3,128); tx = CudaArray(x)
w = rand(Float32,8,8,3,10); tw = CudaArray(w)

dw = zeros(Float32, size(w)); tdw = CudaArray(dw);
dx = zeros(Float32, size(x)); tdx = CudaArray(dx);

dw2 = zeros(Float32, size(w))
dx2 = zeros(Float32, size(x));

y = zeros(Float32, cudnnGetConvolutionNdForwardOutputDim(tx,tw; padding=padding,stride=stride)); ty = CudaArray(y)
y2 = zeros(Float32, size(y))
dy = rand(Float32,size(y)); tdy = CudaArray(dy)

cudnnConvolutionForward(tx,tw,ty; padding=padding, stride=stride); y = to_host(ty)
cudnnConvolutionForward(x,w,y2; padding=padding, stride=stride);
@test_approx_eq y y2

cudnnConvolutionBackwardFilter(tx,tdy,tdw); dw = to_host(tdw)
cudnnConvolutionBackwardFilter(x,dy,dw2)
@test_approx_eq dw dw2

cudnnConvolutionBackwardData(tw, tdy, tdx); dx = to_host(tdx)
cudnnConvolutionBackwardData(w, dy, dx2);
@test_approx_eq dx dx2
=#

psize = 2
padding=0
stride=1
x = rand(Float32,4,4,1,1); tx = CudaArray(x); @show x
@show GetPoolingNdForwardOutputDim(x, window=psize, padding=padding, stride=stride, mode=0)
y = zeros(Float32, GetPoolingNdForwardOutputDim(x, window=psize, padding=padding, stride=stride, mode=0)); ty = CudaArray(y);
y2 = zeros(y)
cudnnPoolingForward(tx, ty; window=psize, padding=padding, stride=stride); y = to_host(ty); @show y
cudnnPoolingForward(x, y2; window=psize, padding=padding, stride=stride); @show y2
@test_approx_eq y y2

dx = zeros(Float32, size(x)); tdx = CudaArray(dx);
dx2 = zeros(Float32, size(x));
dy = rand(Float32, size(y)); tdy = CudaArray(dy); @show dy
cudnnPoolingBackward(ty, tdy, tx, tdx; window=psize, padding=padding, stride=stride, mode=0); dx = to_host(tdx); @show dx
cudnnPoolingBackward(y, dy, x, dx2; window=psize, padding=padding, stride=stride, mode=0); @show dx2
@test_approx_eq dx dx2


text:
module ShuffleProofs

#using Infiltrator

include("utils.jl") # Common functions

include("prover.jl")

include("parser.jl")
include("io.jl") # Some convinience methods

include("verifier.jl")

end # module


text:
push!(LOAD_PATH, "./")
push!(LOAD_PATH, "./policy")
push!(LOAD_PATH, "./pg")
push!(LOAD_PATH, "./opt")


text:
## Example for using TPA-LSTM

# Make sure all the required packages are available
cd(@__DIR__)
using Pkg; Pkg.activate(".")
Pkg.instantiate()

@info "Loading packages"
using FluxArchitectures
using Plots

# Load some sample data
@info "Loading data"
poollength = 10
horizon = 6
datalength = 2000
input, target = get_data(:solar, poollength, datalength, horizon) |> gpu

# Define the network architecture
@info "Creating model and loss"
inputsize = size(input, 1)
hiddensize = 10
layers = 2
filternum = 32
filtersize = 1

# Define the neural net
model = TPALSTM(inputsize, hiddensize, poollength, layers, filternum, filtersize) |> gpu

# MSE loss
function loss(x, y)
    Flux.reset!(model)
    return Flux.mse(model(x), permutedims(y))
end

# Callback for plotting the training
cb = function ()
    Flux.reset!(model)
    pred = model(input) |> permutedims |> cpu
    Flux.reset!(model)
    p1 = plot(pred, label="Predict")
    p1 = plot!(cpu(target), label="Data", title="Loss $(loss(input, target))")
    display(plot(p1))
end

# Training loop
@info "Start loss" loss = loss(input, target)
@info "Starting training"
Flux.train!(loss, Flux.params(model),Iterators.repeated((input, target), 50),
            ADAM(0.02), cb=cb)

@info "Finished"
@info "Final loss" loss = loss(input, target)

max_stars_repo_name: sdobber/FluxArchitectures

max_issues_repo_name: sdobber/FluxArchitectures

max_forks_repo_name: sdobber/FluxArchitectures


text:
# The SteinChampionLenardMillsKernel is based on the reference
# "An introduction to abstract splines" written in 1996 by the
# aforementioned authors. The kernel is a kernel for the Sobolev
# space H^2([0,1]) subject to the extra condition that f(0) = f(1) = 0.
# The kernel is of the form
#
# k(x, y) = (1/6) * [(x - y)_+^3 - x (1-y) (x^2 - 2y + y^2)]
#
# and has inner product
#
# <f, g> = \int_[0,1] f''(t) g''(t) dt
#
# In order to extend this arbitrary intervals, we dilate and translate
# the above kernel so that for the interval [0,a], the kernel
# becomes
#
# k_a(x, y) = k(x/a, y/a)
#
# and the inner product is normalized to become
#
# <f,g>_a = (1/a^3) * int_[0,a] f''(t) g''(t) dt
#
# Translating above is trivial.

# utility for grabbing the support bounds
mutable struct SupportBound
    lb::Float64
    ub::Float64
end

mutable struct SteinChampionLenardMillsKernel <: SteinTensorizedKernel
    support::Array{SupportBound,1}
end

# coerce arrays support bounds
SteinChampionLenardMillsKernel(supp::Array{Array{Float64,1},1}) =
    SteinChampionLenardMillsKernel([SupportBound(s[1], s[2]) for s in supp])

# Make it easy to make d dimensional supports
SteinChampionLenardMillsKernel(supp::Array{Float64,1}, d::Int) =
    SteinChampionLenardMillsKernel([supp for i in 1:d])

# Make it easy to make 1-d supports
SteinChampionLenardMillsKernel(supp::Array{Float64,1}) =
    SteinChampionLenardMillsKernel(supp, 1)

# Assumes two arrays are [lb] and [ub]
SteinChampionLenardMillsKernel(lb::Array{Float64,1}, ub::Array{Float64,1}) =
    SteinChampionLenardMillsKernel([SupportBound(lb[i], ub[i]) for i=1:length(lb)])

# default is just supprt [0,1]
SteinChampionLenardMillsKernel() = SteinChampionLenardMillsKernel([[0.0,1.0]])

#### utility methods
# returns the parameters to translate and dilate the kernel to appropriate interval
function getaffinetransform(ker::SteinChampionLenardMillsKernel, i::Int)
    bt = ker.support[i]
    lb = bt.lb; ub = bt.ub
    lb, (ub - lb)
end

# this is kernel k with domain [0,1]
function k01(x::Float64, y::Float64)
    (max(x - y, 0)^3 - x * (1-y) * (x^2 - 2*y + y^2)) / 6
end

function gradxk01(x::Float64, y::Float64)
    (3 * (x-y)^2 * Int(x >= y) - (1-y) * (3*x^2 - 2*y + y^2)) / 6
end

function gradxyk01(x::Float64, y::Float64)
    0.5 * x^2 + 0.5 * y^2 - max(x,y) + (1/3)
end

### Mandatory methods to implement
function ki(ker::SteinChampionLenardMillsKernel, x::Float64, y::Float64, i::Int)
    c, a = getaffinetransform(ker, i)
    k01((x - c) / a, (y - c) / a)
end

function gradxki(ker::SteinChampionLenardMillsKernel, x::Float64, y::Float64, i::Int)
    c, a = getaffinetransform(ker, i)
    (1/a) * gradxk01((x - c) / a, (y - c) / a)
end

function gradxyki(ker::SteinChampionLenardMillsKernel, x::Float64, y::Float64, i::Int)
    c, a = getaffinetransform(ker, i)
    (1/a^2) * gradxyk01((x - c) / a, (y - c) / a)
end


text:
new_pardiso_handle() = zeros(Int, 64)

new_iparm(::AbstractSparseMatrix{T, Ti}) where {T, Ti} = zeros(Ti, 64)

new_iparm() = zeros(Int, 64)


text:
###############################################################################
#
#   fq_nmod.jl : Flint finite fields
#
###############################################################################

export fq_nmod, FqNmodFiniteField

###############################################################################
#
#   Type and parent object methods
#
###############################################################################

parent_type(::Type{fq_nmod}) = FqNmodFiniteField

elem_type(::Type{FqNmodFiniteField}) = fq_nmod

base_ring(a::FqNmodFiniteField) = Union{}

base_ring(a::fq_nmod) = Union{}

parent(a::fq_nmod) = a.parent

isdomain_type(::Type{fq_nmod}) = true

function check_parent(a::fq_nmod, b::fq_nmod)
   a.parent != b.parent && error("Operations on distinct finite fields not supported")
end

###############################################################################
#
#   Basic manipulation
#
###############################################################################

function Base.hash(a::fq_nmod, h::UInt)
   b = 0x78e5f766c8ace18d%UInt
   for i in 1:degree(parent(a)) + 1
      b = xor(b, xor(hash(coeff(a, i), h), h))
      b = (b << 1) | (b >> (sizeof(Int)*8 - 1))
   end
   return b
end

function coeff(x::fq_nmod, n::Int)
   n < 0 && throw(DomainError())
   return ccall((:nmod_poly_get_coeff_ui, :libflint), UInt,
                (Ref{fq_nmod}, Int), x, n)
end

function zero(a::FqNmodFiniteField)
   d = a()
   ccall((:fq_nmod_zero, :libflint), Void,
         (Ref{fq_nmod}, Ref{FqNmodFiniteField}), d, a)
   return d
end

function one(a::FqNmodFiniteField)
   d = a()
   ccall((:fq_nmod_one, :libflint), Void,
         (Ref{fq_nmod}, Ref{FqNmodFiniteField}), d, a)
   return d
end

function gen(a::FqNmodFiniteField)
   d = a()
   ccall((:fq_nmod_gen, :libflint), Void,
         (Ref{fq_nmod}, Ref{FqNmodFiniteField}), d, a)
   return d
end

iszero(a::fq_nmod) = ccall((:fq_nmod_is_zero, :libflint), Bool,
                     (Ref{fq_nmod}, Ref{FqNmodFiniteField}), a, a.parent)

isone(a::fq_nmod) = ccall((:fq_nmod_is_one, :libflint), Bool,
                    (Ref{fq_nmod}, Ref{FqNmodFiniteField}), a, a.parent)

isgen(a::fq_nmod) = a == gen(parent(a)) # there is no isgen in flint

isunit(a::fq_nmod) = ccall((:fq_nmod_is_invertible, :libflint), Bool,
                     (Ref{fq_nmod}, Ref{FqNmodFiniteField}), a, a.parent)

function characteristic(a::FqNmodFiniteField)
   d = fmpz()
   ccall((:__fq_nmod_ctx_prime, :libflint), Void,
         (Ref{fmpz}, Ref{FqNmodFiniteField}), d, a)
   return d
end

function order(a::FqNmodFiniteField)
   d = fmpz()
   ccall((:fq_nmod_ctx_order, :libflint), Void,
         (Ref{fmpz}, Ref{FqNmodFiniteField}), d, a)
   return d
end

function degree(a::FqNmodFiniteField)
   return ccall((:fq_nmod_ctx_degree, :libflint), Int,
                (Ref{FqNmodFiniteField},), a)
end

function deepcopy_internal(d::fq_nmod, dict::ObjectIdDict)
   z = fq_nmod(parent(d), d)
   z.parent = parent(d)
   return z
end

###############################################################################
#
#   Canonicalisation
#
###############################################################################

canonical_unit(x::fq_nmod) = x

###############################################################################
#
#   AbstractString I/O
#
###############################################################################

function show(io::IO, x::fq_nmod)
   cstr = ccall((:fq_nmod_get_str_pretty, :libflint), Ref{UInt8},
                (Ref{fq_nmod}, Ref{FqNmodFiniteField}), x, x.parent)

   print(io, unsafe_string(cstr))

   ccall((:flint_free, :libflint), Void, (Ref{UInt8},), cstr)
end

function show(io::IO, a::FqNmodFiniteField)
   print(io, "Finite field of degree ", degree(a))
   print(io, " over F_", characteristic(a))
end

needs_parentheses(x::fq_nmod) = x.length > 1

isnegative(x::fq_nmod) = false

show_minus_one(::Type{fq_nmod}) = true

###############################################################################
#
#   Unary operations
#
###############################################################################

function -(x::fq_nmod)
   z = parent(x)()
   ccall((:fq_nmod_neg, :libflint), Void,
       (Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}), z, x, x.parent)
   return z
end

###############################################################################
#
#   Binary operations
#
###############################################################################

function +(x::fq_nmod, y::fq_nmod)
   check_parent(x, y)
   z = parent(y)()
   ccall((:fq_nmod_add, :libflint), Void,
         (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}),
                                                       z, x, y, y.parent)
   return z
end

function -(x::fq_nmod, y::fq_nmod)
   check_parent(x, y)
   z = parent(y)()
   ccall((:fq_nmod_sub, :libflint), Void,
         (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}),
                                                       z, x, y, y.parent)
   return z
end

function *(x::fq_nmod, y::fq_nmod)
   check_parent(x, y)
   z = parent(y)()
   ccall((:fq_nmod_mul, :libflint), Void,
         (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}),
                                                       z, x, y, y.parent)
   return z
end

###############################################################################
#
#   Ad hoc binary operators
#
###############################################################################

function *(x::Int, y::fq_nmod)
   z = parent(y)()
   ccall((:fq_nmod_mul_si, :libflint), Void,
         (Ref{fq_nmod}, Ref{fq_nmod}, Int, Ref{FqNmodFiniteField}),
                                               z, y, x, y.parent)
   return z
end

*(x::fq_nmod, y::Int) = y*x

*(x::Integer, y::fq_nmod) = fmpz(x)*y

*(x::fq_nmod, y::Integer) = y*x

function *(x::fmpz, y::fq_nmod)
   z = parent(y)()
   ccall((:fq_nmod_mul_fmpz, :libflint), Void,
         (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fmpz}, Ref{FqNmodFiniteField}),
                                                    z, y, x, y.parent)
   return z
end

*(x::fq_nmod, y::fmpz) = y*x

+(x::fq_nmod, y::Integer) = x + parent(x)(y)

+(x::Integer, y::fq_nmod) = y + x

+(x::fq_nmod, y::fmpz) = x + parent(x)(y)

+(x::fmpz, y::fq_nmod) = y + x

-(x::fq_nmod, y::Integer) = x - parent(x)(y)

-(x::Integer, y::fq_nmod) = parent(y)(x) - y

-(x::fq_nmod, y::fmpz) = x - parent(x)(y)

-(x::fmpz, y::fq_nmod) = parent(y)(x) - y

###############################################################################
#
#   Powering
#
###############################################################################

function ^(x::fq_nmod, y::Int)
   if y < 0
      x = inv(x)
      y = -y
   end
   z = parent(x)()
   ccall((:fq_nmod_pow_ui, :libflint), Void,
         (Ref{fq_nmod}, Ref{fq_nmod}, Int, Ref{FqNmodFiniteField}),
                                               z, x, y, x.parent)
   return z
end

function ^(x::fq_nmod, y::fmpz)
   if y < 0
      x = inv(x)
      y = -y
   end
   z = parent(x)()
   ccall((:fq_nmod_pow, :libflint), Void,
         (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fmpz}, Ref{FqNmodFiniteField}),
                                                    z, x, y, x.parent)
   return z
end

###############################################################################
#
#   Comparison
#
###############################################################################

function ==(x::fq_nmod, y::fq_nmod)
   check_parent(x, y)
   ccall((:fq_nmod_equal, :libflint), Bool,
       (Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}), x, y, y.parent)
end

###############################################################################
#
#   Ad hoc comparison
#
###############################################################################

==(x::fq_nmod, y::Integer) = x == parent(x)(y)

==(x::fq_nmod, y::fmpz) = x == parent(x)(y)

==(x::Integer, y::fq_nmod) = parent(y)(x) == y

==(x::fmpz, y::fq_nmod) = parent(y)(x) == y

###############################################################################
#
#   Inversion
#
###############################################################################

function inv(x::fq_nmod)
   iszero(x) && throw(DivideError())
   z = parent(x)()
   ccall((:fq_nmod_inv, :libflint), Void,
       (Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}), z, x, x.parent)
   return z
end

###############################################################################
#
#   Exact division
#
###############################################################################

function divexact(x::fq_nmod, y::fq_nmod)
   check_parent(x, y)
   iszero(y) && throw(DivideError())
   z = parent(y)()
   ccall((:fq_nmod_div, :libflint), Void,
         (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}),
                                                       z, x, y, y.parent)
   return z
end

function divides(a::fq_nmod, b::fq_nmod)
   iszero(b) && error("Division by zero in divides")
   return true, divexact(a, b)
end

###############################################################################
#
#   Ad hoc exact division
#
###############################################################################

divexact(x::fq_nmod, y::Integer) = divexact(x, parent(x)(y))

divexact(x::fq_nmod, y::fmpz) = divexact(x, parent(x)(y))

divexact(x::Integer, y::fq_nmod) = divexact(parent(y)(x), y)

divexact(x::fmpz, y::fq_nmod) = divexact(parent(y)(x), y)

###############################################################################
#
#   Special functions
#
###############################################################################

function pth_root(x::fq_nmod)
   z = parent(x)()
   ccall((:fq_nmod_pth_root, :libflint), Void,
       (Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}), z, x, x.parent)
   return z
end

function trace(x::fq_nmod)
   z = fmpz()
   ccall((:fq_nmod_trace, :libflint), Void,
         (Ref{fmpz}, Ref{fq_nmod}, Ref{FqNmodFiniteField}), z, x, x.parent)
   return parent(x)(z)
end

function norm(x::fq_nmod)
   z = fmpz()
   ccall((:fq_nmod_norm, :libflint), Void,
         (Ref{fmpz}, Ref{fq_nmod}, Ref{FqNmodFiniteField}), z, x, x.parent)
   return parent(x)(z)
end

function frobenius(x::fq_nmod, n = 1)
   z = parent(x)()
   ccall((:fq_nmod_frobenius, :libflint), Void,
         (Ref{fq_nmod}, Ref{fq_nmod}, Int, Ref{FqNmodFiniteField}),
                                               z, x, n, x.parent)
   return z
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(z::fq_nmod)
   ccall((:fq_nmod_zero, :libflint), Void,
        (Ref{fq_nmod}, Ref{FqNmodFiniteField}), z, z.parent)
   return z
end

function mul!(z::fq_nmod, x::fq_nmod, y::fq_nmod)
   ccall((:fq_nmod_mul, :libflint), Void,
         (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}),
                                                       z, x, y, y.parent)
   return z
end

function addeq!(z::fq_nmod, x::fq_nmod)
   ccall((:fq_nmod_add, :libflint), Void,
         (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}),
                                                       z, z, x, x.parent)
   return z
end

function add!(z::fq_nmod, x::fq_nmod, y::fq_nmod)
   ccall((:fq_nmod_add, :libflint), Void,
         (Ref{fq_nmod}, Ref{fq_nmod}, Ref{fq_nmod}, Ref{FqNmodFiniteField}),
                                                       z, x, y, x.parent)
   return z
end

###############################################################################
#
#   Promotions
#
###############################################################################

promote_rule(::Type{fq_nmod}, ::Type{T}) where {T <: Integer} = fq_nmod

promote_rule(::Type{fq_nmod}, ::Type{fmpz}) = fq_nmod

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (a::FqNmodFiniteField)()
   z = fq_nmod(a)
   z.parent = a
   return z
end

(a::FqNmodFiniteField)(b::Integer) = a(fmpz(b))

function (a::FqNmodFiniteField)(b::Int)
   z = fq_nmod(a, b)
   z.parent = a
   return z
end

function (a::FqNmodFiniteField)(b::fmpz)
   z = fq_nmod(a, b)
   z.parent = a
   return z
end

function (a::FqNmodFiniteField)(b::fq_nmod)
   parent(b) != a && error("Coercion between finite fields not implemented")
   return b
end

###############################################################################
#
#   FlintFiniteField constructor
#
###############################################################################

function FlintFiniteField(char::Int, deg::Int, s::AbstractString; cached = true)
   S = Symbol(s)
   parent_obj = FqNmodFiniteField(fmpz(char), deg, S, cached)

   return parent_obj, gen(parent_obj)
end

function FlintFiniteField(pol::nmod_poly, s::AbstractString; cached = true)
   S = Symbol(s)
   parent_obj = FqNmodFiniteField(pol, S, cached)

   return parent_obj, gen(parent_obj)
end


text:
@testset "$TEST $G" begin

g = CompleteGraph(5, G)
klist = cores(g)
@test klist == fill(4, 5)

g = PathGraph(3)
klist = cores(g)
@test klist == fill(1, 3)

g = CycleGraph(3)
klist = cores(g)
@test klist == fill(2, 3)

add_vertex!(g)
add_edge!(g, 3, 4)
klist = cores(g)
@test klist == [2, 2, 2, 1]

h, vmap = kcore(g, 2)
@test h == CycleGraph(3)
@test length(vmap) == 0
@test vmap == VertexMap(g, vertextype(g))
end # testset


text:
isrotated(plotattributes) = get(plotattributes, :rotated, false)

function make_label(var_name, selection)
    label = String(ArviZ.arviz.plots.plot_utils.make_label(var_name, selection))
    return replace(label, '\n' => ' ')
end

function update_plot_extra_kwargs!(plotattributes, group, group_extra_kwargs)
    extra_kwargs = Dict{Any,Any}(get(plotattributes, :extra_kwargs, Dict{Any,Any}()))
    group_kwargs = Dict{Any,Any}(get(extra_kwargs, group, Dict{Any,Any}()))
    extra_kwargs[group] = merge(group_extra_kwargs, group_kwargs)
    plotattributes[:extra_kwargs] = extra_kwargs
    return plotattributes
end

get_bins(x) = ArviZ.arviz.stats.density_utils.get_bins(x)

function calculate_point_estimate(point_estimate::Union{String,Symbol}, values; kwargs...)
    return ArviZ.arviz.plots.plot_utils.calculate_point_estimate(point_estimate, values; kwargs...)
end
function calculate_point_estimate(point_estimate, values; kwargs...)
    return calculate_point_estimate("$point_estimate", values; kwargs...)
end

# adapted from Plots
function widen(lmin, lmax, scale=:identity)
    f, invf = RecipesPipeline.scale_func(scale), RecipesPipeline.inverse_scale_func(scale)
    span = f(lmax) - f(lmin)
    eps = max(1e-16, 0.03span)
    return invf(f(lmin) - eps), invf(f(lmax) + eps)
end


text:
using Gurobi
function sudoku_gurobi(X)
# SUDOKU_GUROBI Solve a sedoku puzzle using Gurobi for an Binary LP
#
# S=sudoku_gurobi(X)
#   X has a set of initial conditions where anything ~= 1 ... 9 is an
#   unknown
# S is the solved sudoku puzzle

    assert(all(size(X) == (9, 9)))
    # create variables
    nvars = 9*9*9

    # index the variables for easier access
    vars = zeros(Int64, 9, 9, 9)
    ind = 1
    for i = 1:9
        for j = 1:9
            for v = 1:9
                vars[i, j, v] = ind
                ind = ind + 1
            end
        end
    end

    #obj = ones(Int64, 1, nvars) # no objective, just constraints

    # add known cells as constraints

    A = spzeros(nvars, 0)
    b = ones(1, 0)

    for i = 1:9
        for j = 1:9
            if X[i, j] < 1 || X[i, j] > 9
                # this is a free variable, ignore
            else
                # this is a constraint, add it
                cons = zeros(nvars)
                cons[vars[i, j, X[i, j]]] = 1
                A = [A sparse(cons)] # append a row
                b = [b 1]
            end
        end
    end

    #
    # add sudoku constraints
    #

    # each cell only has one value
    for i = 1:9
        for j = 1:9
            cons = zeros(Int64, nvars)
            for v = 1:9
                cons[vars[i, j, v]] = 1
            end
            A = [A sparse(cons)]
            b = [b 1]
        end
    end

    # each row has 1-9
    for i = 1:9
        for v = 1:9
            cons = zeros(nvars)
            # the value v only appears in each row once
            for j = 1:9
                cons[vars[i, j, v]] = 1
            end
            A = [A sparse(cons)]
            b = [b 1]
        end
    end

    # each column has 1-9
    for j = 1:9
        for v = 1:9
            cons = zeros(nvars)
            # the value v only appears in each row once
            for i = 1:9
                cons[vars[i, j, v]] = 1
            end
            A = [A sparse(cons)]
            b = [b 1]
        end
    end

    # subgrid has 1-9
    for v = 1:9
        for i0 = 0:2
            for j0 = 0:2
                cons = zeros(nvars)
                for i = 1:3
                    for j = 1:3
                        cons[vars[i+3*i0, j+3*j0, v]] = 1
                    end
                end
                A = [A sparse(cons)]
                b = [b 1]
            end
        end
    end

    # setup Gurobi call
    env = Gurobi.Env()

    model = Gurobi.Model(env, "sudoku_blp", :minimize)

    # add binary variables
    add_bvars!(model, zeros(nvars))

    # have the variables incorporated into the model
    update_model!(model)

    add_constrs_t!(model, A, '=', vec(b))

    # run optimization
    optimize(model)

    # get results
    x = get_solution(model)

    # assign solution
    S = zeros(Int64, 9 ,9)
    for i = 1:9
        for j = 1:9
            for v = 1:9
                if x[vars[i, j, v]] > 0
                    S[i, j] = v
                end
            end
        end
    end
    return S
end

##
X = [0 0 8 0 0 0 6 0 0;
     0 2 0 0 0 0 0 1 0;
     5 0 0 6 0 0 0 0 9;
     0 0 7 0 4 0 0 0 3;
     0 4 0 2 0 9 0 6 0;
     6 0 0 0 1 0 8 0 0;
     9 0 0 0 0 7 0 0 4;
     0 5 0 0 0 0 0 7 0;
     0 0 1 0 0 0 5 0 0];
S = sudoku_gurobi(X)

##

X = [0 0 4 9 3 0 0 0 5
     2 8 0 0 0 0 0 6 0
     0 1 0 0 2 0 0 4 0
     0 0 5 0 4 0 0 0 7
     0 2 0 0 0 5 0 0 0
     0 0 0 0 0 6 0 0 0
     0 0 0 1 0 0 0 8 2
     8 0 0 0 9 0 0 0 1
     0 7 0 0 0 0 0 0 0];
S = sudoku_gurobi(X)

##
X = [0 6 0 0 0 0 9 0 0;
     1 7 0 0 0 0 0 0 0;
     0 0 5 0 1 0 0 8 6;
     0 0 0 0 0 0 7 0 0;
     0 0 0 0 8 6 0 9 0;
     0 4 0 0 9 0 0 0 2;
     0 3 8 0 0 0 0 1 0;
     6 0 0 0 0 5 0 2 0;
     0 0 0 0 0 0 3 0 0];
S = sudoku_gurobi(X)

##
# http://www.nytimes.com/crosswords/game/sudoku/hard


text:
# SPDX-License-Identifier: MIT

module Semirings

export Semiring
export Semifield
export ProbabilitySemifield
export LogSemifield
export TropicalSemiring

include("abstractsemiring.jl")
include("logsemifield.jl")
include("probsemifield.jl")
include("tropsemiring.jl")

end


text:
"""
Constructs `Table` out of `FeatureLayer`, where every row is a `Feature` consisting of Geometry and attributes.
```
ArchGDAL.Table(T::Union{IFeatureLayer, FeatureLayer})
```
"""
struct Table{T<:Union{IFeatureLayer, FeatureLayer}}
    layer::T
end

getlayer(t::Table) = Base.getfield(t, :layer)

function Tables.schema(layer::AbstractFeatureLayer)
    field_names, geom_names, featuredefn, fielddefns = schema_names(layer)
    ngeom = ArchGDAL.ngeom(featuredefn)
    geomdefns = (ArchGDAL.getgeomdefn(featuredefn, i) for i in 0:ngeom-1)
    field_types = (_FIELDTYPE[gettype(fielddefn)] for fielddefn in fielddefns)
    geom_types = (IGeometry for i in 1:ngeom)
    Tables.Schema((field_names..., geom_names...), (field_types..., geom_types...))
end

Tables.istable(::Type{<:Table}) = true
Tables.rowaccess(::Type{<:Table}) = true
Tables.rows(t::Table) = t

function Base.iterate(t::Table, st = 0)
    layer = getlayer(t)
    st >= nfeature(layer) && return nothing
    if iszero(st)
        resetreading!(layer)
    end
    return nextnamedtuple(layer), st + 1
end

function Base.getindex(t::Table, idx::Integer)
    layer = getlayer(t)
    setnextbyindex!(layer, idx-1)
    return nextnamedtuple(layer)
end

Base.IteratorSize(::Type{<:Table}) = Base.HasLength()
Base.size(t::Table) = nfeature(getlayer(t))
Base.length(t::Table) = size(t)
Base.IteratorEltype(::Type{<:Table}) = Base.HasEltype()
Base.propertynames(t::Table) = Tables.schema(getlayer(t)).names
Base.getproperty(t::Table, s::Symbol) = [getproperty(row, s) for row in t]

"""
Returns the feature row of a layer as a `NamedTuple`

Calling it iteratively will work similar to `nextfeature` i.e. give the consecutive feature as `NamedTuple`
"""
function nextnamedtuple(layer::IFeatureLayer)
    field_names, geom_names = schema_names(layer)
    return nextfeature(layer) do feature
        prop = (getfield(feature, name) for name in field_names)
        geom = (getgeom(feature, idx-1) for idx in 1:length(geom_names))
        NamedTuple{(field_names..., geom_names...)}((prop..., geom...))
    end
end

function schema_names(layer::AbstractFeatureLayer)
    featuredefn = layerdefn(layer)
    fielddefns = (getfielddefn(featuredefn, i) for i in 0:nfield(layer)-1)
    field_names = (Symbol(getname(fielddefn)) for fielddefn in fielddefns)
    geom_names = (Symbol(getname(getgeomdefn(featuredefn, i-1))) for i in 1:ngeom(layer))
    return (field_names, geom_names, featuredefn, fielddefns)
end

function Base.show(io::IO, t::Table)
    println(io, "Table with $(nfeature(getlayer(t))) features")
end
Base.show(io::IO, ::MIME"text/plain", t::Table) = show(io, t)


text:
using Optim

"""
Linear Regression Cost function with Regularisation

Parameters:
	params      = vector containing rolled up X and Θ
		X       = array of features for each movie
		Θ       = array of weights for each user
	Y           = array of ratings for each movie
	R           = array of the existance of ratings for each movie
	numUsers    = number of users
	numMovies   = number of movies
	numFeatures = number of features
	λ           = regularisation constant

Returns:
	J = linear regression cost of the given weights with the given regularisation constant
"""
function cost(params, Y, R, numUsers, numMovies, numFeatures, λ)
	# Unfold params -> X, Θ
	X, Θ = unroll(params, numUsers, numMovies, numFeatures)

	# Calculate cost
	J = 1/2 * sum(R .* (X * Θ' - Y).^2)
	# add regularisation
	J += λ/2 * (sum(Θ.^2) + sum(X.^2))

	if rand() > 0.95
		@show J
	end

	return J
end

"""
Linear Regression Gradient with Regularisation

Parameters:
	params      = vector containing rolled up X and Θ
		X       = array of features for each movie
		Θ       = array of weights for each user
	gradient    = vector containing rolled up Xgrad and ΘGrad from previous iteration
	Y           = array of ratings for each movie
	R           = array of the existance of ratings for each movie
	numUsers    = number of users
	numMovies   = number of movies
	numFeatures = number of features
	λ           = regularisation constant

Modifies:
	gradient = XGrad and ΘGrad are updated based on the given weights and regularisation costant
"""
function grad!(params, gradient, Y, R, numUsers, numMovies, numFeatures, λ)
	# Unfold params -> X, Θ
	X, Θ = unroll(params, numUsers, numMovies, numFeatures)

	XGrad = zeros(size(X))
	ΘGrad = zeros(size(Θ))

	#batchSize = min(1000, numMovies - 1)
	#start = batchSize >= numMovies ? 1 : rand(1:numMovies - batchSize)
	# Calculate XGrad
	#for i = start:start + batchSize
	for i = 1:numMovies
		idx = find(R[i, :])
		ΘTemp = Θ[idx, :]
		YTemp = Y[i, idx];
		XGrad[i, :] = (X[i, :] * ΘTemp' - YTemp) * ΘTemp + λ * X[i,:];
	end

	#batchSize = min(200, numUsers - 1)
	#start = batchSize >= numUsers ? 1 : rand(1:numUsers - batchSize)
	# Calculate ΘGrad
	#for i = start:start + batchSize
	for i = 1:numUsers
		idx = find(R[:, i])
		XTemp = X[idx, :];
		YTemp = Y[idx, i];
		ΘGrad[i, :] = XTemp' * (XTemp * Θ[i, :]' - YTemp) + λ * Θ[i, :]';
	end

	gradient[:] = rollup(XGrad, ΘGrad)
end

function costGrad!(params, gradient, Y, R, numUsers, numMovies, numFeatures, λ)
	# Unfold params -> X, Θ
	X, Θ = unroll(params, numUsers, numMovies, numFeatures)

	XGrad = zeros(size(X))
	ΘGrad = zeros(size(Θ))

	start = rand(1:size(X)[1]-200)

	#tic()
	# Calculate XGrad
	for i = start:start+200
		idx = find(R[i, :])
		ΘTemp = Θ[idx, :]
		YTemp = Y[i, idx];
		XGrad[i, :] = (X[i, :] * ΘTemp' - YTemp) * ΘTemp + λ * X[i,:];
	end

	#toc()
	start = rand(1:size(Θ)[1]-200)
	#tic()
	# Calculate ΘGrad
	for i = start:start+200
		idx = find(R[:, i])
		XTemp = X[idx, :];
		YTemp = Y[idx, i];
		ΘGrad[i, :] = XTemp' * (XTemp * Θ[i, :]' - YTemp) + λ * Θ[i, :]';
	end

	#toc()
	#tic()
	gradient[:] = rollup(XGrad, ΘGrad)

	# Calculate cost
	J = 1/2 * sum(R .* (X * Θ' - Y).^2)
	# add regularisation
	J += λ/2 * (sum(Θ.^2) + sum(X.^2))

	if rand() > 0.99
		@show J
	end

	#toc()
	return J
end

"""
Unroll parameter vector into X and Θ

Parameters:
	params = vector containing rolled up X and Θ
	numUsers    = number of users
	numMovies   = number of movies
	numFeatures = number of features

Returns:
	X = array of features for each movie
	Θ = array of weights for each user
"""
function unroll(params, numUsers, numMovies, numFeatures)
	X = reshape(params[1 : numMovies * numFeatures], numMovies, numFeatures)
	Θ = reshape(params[numMovies * numFeatures + 1 : end], numUsers, numFeatures)
	return X, Θ
end

"""
Roll up X and Θ into single vector

Parameters:
	X = array of features for each movie
	Θ = array of weights for each user

Returns:
	params = vector containing rolled up X and Θ
"""
function rollup(X, Θ)
	return params = vcat(X..., Θ...)
end

"""
Entry point

Params:
	dataPath = path to() CSVs
"""
function main(dataPath)
	# Load data from CSV - using wc to speed up loading by preallocating
	# N.B. requires cat and wc to be installed - a shell should not be nessecary though
	lines = parse(readall(pipeline(`cat $[dataPath]ratings.csv`, `wc -l`)))
	@time ratings = readdlm(dataPath * "ratings.csv", ',', header=true, dims=(lines, 4))
	ratingsHeader = ratings[2]
	ratings = ratings[1]

	lines = parse(readall(pipeline(`cat $[dataPath]movies.csv`, `wc -l`)))
	@time movies = readdlm(dataPath * "movies.csv", ',', header=true, dims=(lines, 3))
	moviesHeader = movies[2]
	movies = movies[1]

	lines = parse(readall(pipeline(`cat $[dataPath]links.csv`, `wc -l`)))
	@time links = readdlm(dataPath * "links.csv", ',', header=true, dims=(lines, 3))
	linksHeader = links[2]
	links = links[1]

	lines = parse(readall(pipeline(`cat $[dataPath]tags.csv`, `wc -l`)))
	@time tags = readdlm(dataPath * "tags.csv", ',', header=true, dims=(lines, 4))
	tagsHeader = tags[2]
	tags = tags[1]

	# Count users, movies and features
	numUsers = length(union(ratings[:,1]))
	numMovies = size(movies)[1]
	numFeatures = 18 # Leave this for now, 18 initiall chosen b/c 18 listed genres

	# Random initialisation
	X = rand(numMovies, numFeatures)
	Θ = rand(numUsers, numFeatures)

	Y = zeros(Int8, numMovies, numUsers)
	for i = 1:size(ratings)[1]
		Y[Int(ratings[i,2]), Int(ratings[i,1])] = Int8(ratings[i,3] * 2)
	end

	R = Y .!= 0

	@time @show cost(rollup(X, Θ), Y, R, numUsers, numMovies, numFeatures, 0.1)
	@time grad!(rollup(X, Θ), rollup(zeros(X), zeros(Θ)),Y, R, numUsers, numMovies, numFeatures, 0.1)

	@time @show results = optimize(x -> cost(x, Y, R, numUsers, numMovies, numFeatures, 0.1), (x,y) -> grad!(x,y, Y, R, numUsers, numMovies, numFeatures, 0.1), rollup(X,Θ), iterations=10_000)
	#@time @show optimize(DifferentiableFunction(x -> cost(x, Y, R, numUsers, numMovies, numFeatures, 0.1), (x,y) -> grad!(x,y, Y, R, numUsers, numMovies, numFeatures, 0.1), (x,y) -> costGrad!(x,y, Y, R, numUsers, numMovies, numFeatures, 0.1)), rollup(X,Θ), iterations=100)
#=
	gradient = rollup(ones(size(X)), ones(size(Θ)))
	params = rollup(X, Θ)
	α = 0.0001
	tic()
	for i = 1:10_000
		grad!(params, gradient, Y, R, numUsers, numMovies, numFeatures, 0.1)
		params -= α * gradient
		params[1:10]
		if i % 100 == 0
			toc()
			gradient[1:10]
			@show cost(params, Y, R, numUsers, numMovies, numFeatures, 0.1)
			tic()
		end
	end=#

	X, Θ = unroll(results.minimum, numUsers, numMovies, numFeatures)
	@show X[1,:] * Θ[1,:]'
	@show Y[1,1]
end

main("../ml-latest-small/fixed/")


text:
using SparseADRules:imul!, idot
const approx_rtol = 100*eps()

@testset "matrix multiplication" begin
    for T in (Float64, ComplexF64)
        for i = 1:5
            A = sprand(T, 10, 5, 0.5)
            b = rand(T, 5)
            outv = A*b
            c = zero(outv)
            ioutv = imul!(copy(c), A, b, 1 ,1)[1]
            @test check_inv(imul!, (copy(c), A, b, 1, 1))
            @test ≈(ioutv, outv, rtol=approx_rtol)

            B = rand(T, 5, 3)
            outm = A*B
            C = zero(outm)
            ioutm = imul!(copy(C), A, B, 1, 1)[1]
            @test check_inv(imul!, (copy(C), A, B, 1, 1))
            @test ≈(ioutm, outm, rtol=approx_rtol)
        end
    end
end

@testset "adjoint matrix multiplication" begin
    for T in (Float64, ComplexF64)
        for i = 1:5
            A = sprand(T, 5, 5, 0.2)
            b = rand(T, 5)
            outv = A'*b
            c= zero(outv)
            ioutv = imul!(copy(c), A', b, 1, 1)[1]
            @test check_inv(imul!, (copy(c), A', b, 1, 1))
            @test ≈(ioutv, outv, rtol=approx_rtol)

            B = rand(T, 5, 3)
            outm = A'*B
            C = zero(outm)
            ioutm = imul!(copy(C), A', B, 1, 1)[1]
            @test check_inv(imul!, (copy(C), A', B, 1, 1))
            @test ≈(ioutm, outm, rtol=approx_rtol)
        end
    end
end

@testset "dense matrix-sparse matrix multiplication" begin
    for T in (Float64, ComplexF64)
        for i = 1:5
            B = rand(T, 10, 10)
            A = sprand(T, 10, 5, 0.5)
            outm = B*A
            C = zero(outm)
            ioutm = imul!(copy(C), B, A, 1 ,1)[1]
            @test check_inv(imul!, (copy(C), B, A, 1, 1))
            @test ≈(ioutm, outm, rtol=approx_rtol)
        end
    end
end

@testset "adjoint dense matrix-sparse matrix multiplication" begin
    for T in (Float64, ComplexF64)
        for i = 1:5
            B = rand(T, 10, 5)
            A = sprand(T, 10, 5, 0.5)
            outm = B'*A
            C = zero(outm)
            ioutm = imul!(copy(C), B', A, 1 ,1)[1]
            @test check_inv(imul!, (copy(C), B', A, 1, 1))
            @test ≈(ioutm, outm, rtol=approx_rtol)
        end
    end
end

@testset "dense matrix - adjoint sparse matrix multiplication" begin
    for T in (Float64, ComplexF64)
        for i = 1:5
            B = rand(T, 10, 10)
            A = sprand(T, 5, 10, 0.5)
            outm = B*A'
            C = zero(outm)
            ioutm = imul!(copy(C), B, A', 1 ,1)[1]
            @test check_inv(imul!, (copy(C), B, A', 1, 1))
            @test ≈(ioutm, outm, rtol=approx_rtol)
        end
    end
end

@testset "dense vector - sparse matrix dot" begin
    for T in (Float64, ComplexF64)
        for i = 1:50
            x = rand(T, 10)
            A = sprand(T, 10, 5, 0.2)
            y = rand(T, 5)
            outd = dot(x, A, y)
            r = zero(T)
            ioutd = idot(r, x, A, y)[1]
            @test check_inv(idot, (r, x, A, y))
            @test ≈(ioutd, outd, rtol=approx_rtol)
        end
    end
end

@testset "sparse vector - sparse matrix dot" begin
    for T in (Float64, ComplexF64)
        for i = 1:50
            x = sprand(T, 10, 0.2)
            A = sprand(T, 10, 5, 0.2)
            y = sprand(T, 5, 0.3)
            outd = dot(x, A, y)
            r = zero(T)
            ioutd = idot(r, x, A, y)[1]
            @test check_inv(idot, (r, x, A, y))
            @test ≈(ioutd, outd, rtol=approx_rtol)
        end
    end
end

@testset "idot - sparse sparse" begin
    for T in (Float64, ComplexF64)
        for i = 1:5
            A = sprand(T, 10, 5, 0.2)
            B = sprand(T, 10, 5, 0.2)
            outd = dot(A, B)
            r = zero(T)
            ioutd = idot(r, A, B)[1]
            @test check_inv(idot, (r, A, B))
            @test ≈(ioutd, outd, rtol=approx_rtol)
        end
    end
end


text:
#include("../src/ttv.jl")
#include("/Users/ericagol/Computer/Julia/regress.jl")

#@testset "ttv_cartesian" begin

# This routine takes derivative of transit times with respect
# to the initial cartesian coordinates of bodies. [x]
#n = 8
n = 3
t0 = 7257.93115525-7300.0
#h  = 0.12
h  = 0.04
#h  = 0.02
#tmax = 600.0
tmax = 40000.0
#tmax = 100.0
#tmax = 10.0

# Read in initial conditions:
elements = readdlm("elements.txt",',')

# Make an array, tt,  to hold transit times:
# First, though, make sure it is large enough:
ntt = zeros(Int64,n)
for i=2:n
  ntt[i] = ceil(Int64,tmax/elements[i,2])+3
end
println("ntt: ",ntt)
tt  = zeros(n,maximum(ntt))
tt2 = zeros(n,maximum(ntt))
# Save a counter for the actual number of transit times of each planet:
count = zeros(Int64,n)
# Call the ttv function:
rstar = 1e12

# Now, compute derivatives (with respect to initial cartesian positions/masses):
dtdq0 = zeros(n,maximum(ntt),7,n)
@time dtdelements = ttv_elements!(n,t0,h,tmax,elements,tt,count,dtdq0,rstar)
# Now with half the timestep:
#dtdq2 = zeros(n,maximum(ntt),7,n)
#@time dtdelements = ttv_elements!(n,t0,h/2,tmax,elements,tt2,count,dtdq2,rstar)

tt_big = big.(tt); elementsbig = big.(elements); rstarbig = big(rstar)
tmaxbig = big(tmax); t0big = big(t0); hbig = big(h)

# Compute the derivatives in BigFloat precision to see whether finite difference
# derivatives or Float64 derivatives are imprecise at the start:
dtdq0_big = zeros(BigFloat,n,maximum(ntt),7,n)
hbig = big(h); tt_big = big.(tt); elementsbig = big.(elements); rstarbig = big(rstar)
@time dtdelements_big = ttv_elements!(n,t0big,hbig,tmaxbig,elementsbig,tt_big,count,dtdq0_big,rstarbig)
# Now with half the timestep:
dtdq2_big = zeros(BigFloat,n,maximum(ntt),7,n)
tt2_big = big.(tt2)
#@time dtdelements_big = ttv_elements!(n,t0big,hbig/2,tmaxbig,elementsbig,tt2_big,count,dtdq2_big,rstarbig)

using PyPlot

clf()
# Plot the difference in the TTVs:
for i=2:3
  diff1 = abs.(tt[i,2:count[i]].-tt_big[i,2:count[i]])/h;
  loglog(tt[i,2:count[i]]-tt[i,1],diff1);
#  diff2 = abs.(tt2[i,2:count[i]].-tt2_big[i,2:count[i]])/elements[i,2];
#  loglog(tt[i,2:count[i]]-tt[i,1],diff2);
end
xlabel("Time since start")
ylabel(L"(time_{dbl}-time_{big})/timestep")
#loglog([1.0,1024.0],2e-15*[1,2^15],":")
loglog([1.0,40000.0],0.5e-16*([1.0,40000.0]/h).^1.5,":")

clf()
nmed = 10
for i=2:3, k=1:7, l=1:3
  if maximum(abs.(dtdq0[i,2:count[i],k,l])) > 0
#    diff1 = abs.(asinh.(dtdq0_big[i,2:count[i],k,l])-asinh.(dtdq0[i,2:count[i],k,l]));
#    diff1 = abs.(dtdq0_big[i,2:count[i],k,l]-dtdq0[i,2:count[i],k,l])/maximum(dtdq0[i,2:count[i],k,l]);
    diff1 = abs.(dtdq0_big[i,2:count[i],k,l]-dtdq0[i,2:count[i],k,l])
    diff2 = copy(diff1);
    diff3 = abs.(dtdelements_big[i,2:count[i],k,l]-dtdelements[i,2:count[i],k,l])
#    diff3 = abs.(asinh.(dtdelements_big[i,2:count[i],k,l])-asinh.(dtdelements[i,2:count[i],k,l]))
    diff4 = copy(diff3);
    ntt1 = size(diff1)[1];
    for it=nmed:ntt1-nmed
      # Divide by a median smoothed dtdq0:
      diff2[it] = diff1[it]/maximum(abs.(dtdq0[i,it+1-nmed+1:it+1+nmed,k,l]))
      maxit = maximum(abs.(dtdelements[i,it+1-nmed+1:it+1+nmed,k,l]))
      if maxit > 0
        diff4[it] = diff3[it]/maximum(abs.(dtdelements[i,it+1-nmed+1:it+1+nmed,k,l]))
      end
    end;
    diff2[1:nmed-1] .= diff2[nmed]; diff2[ntt1-nmed+1:ntt1] .= diff2[ntt1-nmed];
    diff4[1:nmed-1] .= diff4[nmed]; diff4[ntt1-nmed+1:ntt1] .= diff4[ntt1-nmed];
#    diff3 = abs.(asinh.(dtdq2_big[i,2:count[i],k,l])-asinh.(dtdq2[i,2:count[i],k,l]));
    loglog(tt[i,2:count[i]]-tt[i,1],diff2,linestyle=":");
    loglog(tt[i,2:count[i]]-tt[i,1],diff4,linestyle="--");
#    loglog(tt[i,2:count[i]]-tt[i,1],diff3);
    println("xvs: ",i," ",k," ",l," asinh error h  : ",convert(Float64,maximum(diff2))); #read(STDIN,Char);
    println("els: ",i," ",k," ",l," asinh error h  : ",convert(Float64,maximum(diff4))); #read(STDIN,Char);
#    println(i," ",k," ",l," asinh error h/2: ",convert(Float64,maximum(diff3))); #read(STDIN,Char);
  end
end
loglog([1.0,40000.0],0.5e-16*([1.0,40000.0]/h).^1.5)

#axis([1,1024,1e-19,1e-9])
axis([1,31600,1e-19,2.5e-7])
xlabel("Time since start")
#ylabel("median(asinh(dt/dq)_dbl - asinh(dt/dq)_big,20)")
ylabel("(dt/dq_{dbl}-dt/dq_{big})/maximum(dt/dq_{dbl}[i-9:i+10])")

# Plot a line that scales as time^{3/2}:

# loglog([1.0,1024.0],1e-12*[1,2^15],":",linewidth=3)
# loglog([1.0,1024.0],1e-12*[1,2^15],":",linewidth=3)

println("Max diff asinh(dtdq0): ",maximum(abs.(asinh.(dtdq0_big)-asinh.(dtdq0))))
@test isapprox(asinh.(dtdq0),asinh.(convert(Array{Float64,4},dtdq0_big));norm=maxabs)
#end


text:
using Documenter
using DocumenterLaTeX
using ProbabilisticCircuits
using Literate

#######################################
# 1/ generate the top-level README.md
#######################################

source_dir = "$(@__DIR__)/src"

"replace script includes with file content in Literate code"
function replace_includes(str)
    pat = r"include\(\"(.*)\"\)"
    m = match(pat, str)
    while !isnothing(m)
        str = replace(str, "$(m.match)" =>
                read("$source_dir/$(m[1])", String))
        m = match(pat, str)
    end
    str
end

"hide `#plot` lines in Literate code"
function hide_plots(str)
    str = replace(str, r"#plot (.)*[\n\r]" => "")
    replace(str, r"#!plot (.*)[\n\r]" => s"\g<1>\n")
end

"show `#plot` lines in Literate code"
function show_plots(str)
    str = replace(str, r"#!plot (.)*[\n\r]" => "")
    replace(str, r"#plot (.*)[\n\r]" => s"\g<1>\n")
end

Literate.markdown("$source_dir/README.jl", "$(@__DIR__)/../"; documenter=false, credit=false, execute=true,
    preprocess = hide_plots ∘ replace_includes)

# The DOCSARGS environment variable can be used to pass additional arguments to make.jl.
# This is useful on CI, if you need to change the behavior of the build slightly but you
# can not change the .travis.yml or make.jl scripts any more (e.g. for a tag build).
if haskey(ENV, "DOCSARGS")
    for arg in split(ENV["DOCSARGS"])
        (arg in ARGS) || push!(ARGS, arg)
    end
end

const  pages = [
    "Home" => "index.md",
    "Installation" => "installation.md",
    "Manual" => [
        "manual/demo.md",
        "manual/queries.md",
        "manual/learning.md"
    ],
    "API" => [
        "api/public.md",
        "api/types.md",
        "Internals" => map(
            s -> "api/internals/$(s)",
            sort(readdir(joinpath(@__DIR__, "src/api/internals")))
        ),
    ],
];

const format = if ("pdf" in ARGS)
    LaTeX(platform  = "native")
else
   Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = !("local" in ARGS),
        canonical = "https://juice-jl.github.io/ProbabilisticCircuits.jl/stable/",
        assets = ["assets/favicon.ico"],
        analytics = "UA-136089579-2",
        highlights = ["yaml"],
        collapselevel = 1,
    )
end

makedocs(
    sitename = "ProbabilisticCircuits.jl",
    format = format,
    doctest = true,
    modules = [ProbabilisticCircuits],
    pages = pages,
    linkcheck_ignore = [
        # We'll ignore links that point to GitHub's edit pages, as they redirect to the
        # login screen and cause a warning:
        r"https://github.com/([A-Za-z0-9_.-]+)/([A-Za-z0-9_.-]+)/edit(.*)"
    ]
)

# Documenter can also automatically deploy documentation to gh-pages.
# See "Hosting Documentation" and deploydocs() in the Documenter manual
# for more information.
deploydocs(
    target = "build",
    repo = "github.com/Juice-jl/ProbabilisticCircuits.jl.git",
    branch = "gh-pages",
    devbranch = "master",
    devurl = "dev",
    versions = ["stable" => "v^", "v#.#"],
)


text:
using Revise, BeamPropagation, Distributions, StatsBase, StaticArrays, PhysicalConstants.CODATA2018, Plots, BenchmarkTools

const h = PlanckConstant.val
const ħ = h / 2π
const λ = @with_unit 626 "nm"
const k = 2π / λ
const m = @with_unit 57 "u"
const Δv = ħ * k / m

vz_μ = @with_unit 125 "m/s"
vz_σ = @with_unit 25 "m/s"
vxy_μ = @with_unit 0 "m/s"
vxy_σ = @with_unit 25 "m/s"
exit_radius = @with_unit 4 "mm"

const r = (Normal(0, exit_radius/2), Normal(0, exit_radius/2), Normal(0, 0))
const v = (Normal(vxy_μ, vxy_σ), Normal(vxy_μ, vxy_σ), Normal(vz_μ, vz_σ))
const a = (Normal(0, 0), Normal(0, 0), Normal(0, 0))

VBRs = Weights([
    0.9457,     # to 000
    0.0447,     # 100
    3.9e-3,     # 0200
    2.7e-3,     # 200
    9.9e-4,     # 0220
    11.3e-4,    # 0110
    3.9e-4,     # 1200
    1.5e-4,     # 300
    1.3e-4,     # 1220
    0.7e-4,     # 110 (only N=1, assuming roughly 2/3 to 1/3 rotational branching)
    0.4e-4,     # 110, N=2
    5.7e-5,     # 220
    4.3e-5      # other states
])

addressed = [
    true,   # to 000
    true,   # 100
    true,   # 0200
    true,   # 200
    true,   # 0220
    true,   # 0110
    true,   # 1200
    true,   # 300
    true,   # 1220
    true,   # 110
    false,   # 110, N=2
    false,   # 220
    false   # other states
]

transverse = [
    false,  # to 000
    false,  # 100
    false,  # 0200
    false,  # 200
    false,  # 0220
    false,  # 0110
    false,  # 1200
    false,  # 300
    true,   # 1220
    true,   # 110
    false,   # 110, N=2
    false,   # 220
    false   # other states
]

@inline function random_unit3Dvector()
    θ = rand(Uniform(0, 2π))
    z = rand(Uniform(-1, 1))
    return @SVector [sqrt(1-z^2)*cos(θ), sqrt(1-z^2)*sin(θ), z]
end

@inline function transverse_on(z)
    if 0.175 < z < 0.20     # 1st transverse region 17.5 - 20 cm after cell
        return true
    elseif 0.45 < z < 0.50  # 2nd transverse region 45 - 50 cm after cell
        return true
    elseif 0.56 < z < 0.61  # 2nd transverse region 56 - 61 cm after cell
        return true
    else
        return false
    end
end

const detect_rad  = @with_unit 0.5 "cm"
const detect_zloc = @with_unit 70 "cm"
const detect_zlen = @with_unit 0.5 "cm"
const dead_rad = @with_unit 1 "cm"
const dead_len = detect_zloc + detect_zlen

n = Int64(40e6)
scattering_rate = @with_unit 1.2 "MHz"
save_every      = 100
delete_every    = 15
dt              = 1 / scattering_rate
max_steps       = Int64(3.5e4)

@inline function simple_prop(r, v)
    dist_detect = detect_zloc - r[3]
    x_final = r[1] + v[1] * dist_detect / v[3]
    y_final = r[2] + v[2] * dist_detect / v[3]
    return sqrt(x_final^2 + y_final^2)
end
@inline discard(r, v) = (simple_prop(r, v) > dead_rad) || (r[3] > dead_len)
@inline is_detectable(r) = v_perp(r) < detect_rad && (detect_zloc + detect_zlen > r[3] > detect_zloc)

@inline function save(i, r, v, a, state, s)
    s.vzs[i] = v[3]
    s.detectable[i] = simple_prop(r, v) < detect_rad
    return nothing
end

@inline function f(i, r, v, a, state, dt, p, s)
    if is_detectable(r) && ~s.detectable[i]
        s.detectable[i] = true
    end
    if p.addressed[state] & (~p.transverse[state] | transverse_on(r[3]))
        state′ = sample(1:13, p.VBRs)
        v′ = @SVector [v[1], v[2], v[3] - Δv]
        v′ += Δv .* random_unit3Dvector()
        s.photons[i] += 1
    else
        state′ = state
        v′ = v
    end
    s.vzs[i] = v′[3]
    s.states[i] = state′
    return (state′, v′, a)
end

vzs           = zeros(Float64, n)
photons       = zeros(Int64, n)
detectable    = zeros(Bool, n)
states        = ones(Int64, n)

p = @params (VBRs, transverse, addressed)
s = @params (vzs, photons, detectable, states)

@time s₀, sₜ = propagate!(n, r, v, a, f, save, discard, delete_every, dt, max_steps, p, s)

bright = [addressed[final_state] for final_state in sₜ.states]

barhist(s₀.vzs[s₀.detectable], bins=0:25/4:225, alpha=0.5, xlim=[0,225])
barhist!(sₜ.vzs[sₜ.detectable], bins=0:25/4:225, alpha=0.5)
barhist!(sₜ.vzs[sₜ.detectable .& bright], bins=0:25/4:225, alpha=0.5)

sum(abs.(s₀.vzs[s₀.detectable]) .< 10) / n
sum(abs.(sₜ.vzs[sₜ.detectable]) .< 10) / n
sum(abs.(sₜ.vzs[sₜ.detectable .& bright]) .< 10) / n

v_class = 0 .< sₜ.vzs .< 100
histogram2d(s₀.vzs[sₜ.detectable .& bright .& v_class], sₜ.vzs[sₜ.detectable .& bright.& v_class], bins=20)


text:
# Active patterns
module Active
using MLStyle
using MLStyle.Infras
using MLStyle.MatchCore
using MLStyle.Qualification
using MLStyle.TypeVarExtraction

export @active, def_active_pattern

function def_active_pattern(qualifier_ast, case, active_body, mod)
    (case_name, IDENTS, param) = @match case begin
        :($(case_name :: Symbol)($param)) => (case_name, nothing, param)
        :($(case_name :: Symbol){$(idents...)}($param)) => (case_name, idents, param)
    end
    TARGET = mangle(mod)

    if !isdefined(mod, case_name)
        mod.eval(quote struct $case_name end end)
    end
    qualifier = get_qualifier(qualifier_ast, mod)
    case_obj = getfield(mod, case_name)
    if IDENTS === nothing
        def_app_pattern(mod,
            predicate = (hd_obj, args) -> hd_obj === case_obj,
            rewrite = (tag, hd_obj, args, mod) -> begin
                    (test_var, pat) = @match length(args) begin
                        0 => (false, true)
                        1 => (nothing, args[1])
                        _ => (nothing, Expr(:tuple, args...))
                    end

                    function (body)
                        @format [tag, test_var, param, TARGET, active_body, body] quote
                            let  TARGET =
                                let param = tag
                                    active_body
                                end
                                TARGET === test_var ?  failed : body
                            end
                        end
                    end ∘ mk_pattern(TARGET, pat, mod)
                end,
            qualifiers = Set([qualifier]))
    else
        n_idents = length(IDENTS)
        def_gapp_pattern(mod,
            predicate = (spec_vars, hd_obj, args) -> hd_obj === case_obj && length(spec_vars) === n_idents,
            rewrite   = (tag, forall, spec_vars, hd_obj, args, mod) -> begin
                    (test_var, pat) = @match length(args) begin
                        0 => (false, true)
                        1 => (nothing, args[1])
                        _ => (nothing, Expr(:tuple, args...))
                    end
                    assign_elts_and_active_body =
                        let arr = [:($IDENT = $(spec_vars[i]))
                                    for (i, IDENT)
                                    in enumerate(IDENTS)]
                            Expr(:let, Expr(:block, arr...), Expr(:block, active_body))
                        end
                    function (body)
                        @format [tag, test_var, param, TARGET, assign_elts_and_active_body, body] quote
                            let TARGET =
                                let param = tag
                                    assign_elts_and_active_body
                                end
                                TARGET === test_var ?  failed : body
                            end
                        end
                    end ∘ mk_pattern(TARGET, pat, mod)
            end,
        qualifiers = Set([qualifier]))
    end
    nothing
end

"""
Simple active pattern implementation.
You can give a qualifier in the first argument of `@active` to customize its visibility in other modules.
```julia
    @active F(x) begin
        if x > 0
            nothing
        else
            :ok
        end
    end

    @match -1 begin
        F(:ok) => false
        _ => true
    end # true

    @active public IsEven(x) begin
        x % 2 === 0
    end

    @match 4 begin
        IsEven() => :ok
        _ => :err
    end # :ok
```
"""
macro active(qualifier, case, active_body)
   def_active_pattern(qualifier, case, active_body, __module__)
end

macro active(case, active_body)
    def_active_pattern(:public, case, active_body, __module__)
end

end


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

mutable struct IoK8sApiBatchV1JobSpec <: SwaggerModel
    activeDeadlineSeconds::Any # spec type: Union{ Nothing, Int64 } # spec name: activeDeadlineSeconds
    backoffLimit::Any # spec type: Union{ Nothing, Int32 } # spec name: backoffLimit
    completions::Any # spec type: Union{ Nothing, Int32 } # spec name: completions
    manualSelector::Any # spec type: Union{ Nothing, Bool } # spec name: manualSelector
    parallelism::Any # spec type: Union{ Nothing, Int32 } # spec name: parallelism
    selector::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1LabelSelector } # spec name: selector
    template::Any # spec type: Union{ Nothing, IoK8sApiCoreV1PodTemplateSpec } # spec name: template
    ttlSecondsAfterFinished::Any # spec type: Union{ Nothing, Int32 } # spec name: ttlSecondsAfterFinished

    function IoK8sApiBatchV1JobSpec(;activeDeadlineSeconds=nothing, backoffLimit=nothing, completions=nothing, manualSelector=nothing, parallelism=nothing, selector=nothing, template=nothing, ttlSecondsAfterFinished=nothing)
        o = new()
        validate_property(IoK8sApiBatchV1JobSpec, Symbol("activeDeadlineSeconds"), activeDeadlineSeconds)
        setfield!(o, Symbol("activeDeadlineSeconds"), activeDeadlineSeconds)
        validate_property(IoK8sApiBatchV1JobSpec, Symbol("backoffLimit"), backoffLimit)
        setfield!(o, Symbol("backoffLimit"), backoffLimit)
        validate_property(IoK8sApiBatchV1JobSpec, Symbol("completions"), completions)
        setfield!(o, Symbol("completions"), completions)
        validate_property(IoK8sApiBatchV1JobSpec, Symbol("manualSelector"), manualSelector)
        setfield!(o, Symbol("manualSelector"), manualSelector)
        validate_property(IoK8sApiBatchV1JobSpec, Symbol("parallelism"), parallelism)
        setfield!(o, Symbol("parallelism"), parallelism)
        validate_property(IoK8sApiBatchV1JobSpec, Symbol("selector"), selector)
        setfield!(o, Symbol("selector"), selector)
        validate_property(IoK8sApiBatchV1JobSpec, Symbol("template"), template)
        setfield!(o, Symbol("template"), template)
        validate_property(IoK8sApiBatchV1JobSpec, Symbol("ttlSecondsAfterFinished"), ttlSecondsAfterFinished)
        setfield!(o, Symbol("ttlSecondsAfterFinished"), ttlSecondsAfterFinished)
        o
    end
end # type IoK8sApiBatchV1JobSpec

const _property_map_IoK8sApiBatchV1JobSpec = Dict{Symbol,Symbol}(Symbol("activeDeadlineSeconds")=>Symbol("activeDeadlineSeconds"), Symbol("backoffLimit")=>Symbol("backoffLimit"), Symbol("completions")=>Symbol("completions"), Symbol("manualSelector")=>Symbol("manualSelector"), Symbol("parallelism")=>Symbol("parallelism"), Symbol("selector")=>Symbol("selector"), Symbol("template")=>Symbol("template"), Symbol("ttlSecondsAfterFinished")=>Symbol("ttlSecondsAfterFinished"))
const _property_types_IoK8sApiBatchV1JobSpec = Dict{Symbol,String}(Symbol("activeDeadlineSeconds")=>"Int64", Symbol("backoffLimit")=>"Int32", Symbol("completions")=>"Int32", Symbol("manualSelector")=>"Bool", Symbol("parallelism")=>"Int32", Symbol("selector")=>"IoK8sApimachineryPkgApisMetaV1LabelSelector", Symbol("template")=>"IoK8sApiCoreV1PodTemplateSpec", Symbol("ttlSecondsAfterFinished")=>"Int32")
Base.propertynames(::Type{ IoK8sApiBatchV1JobSpec }) = collect(keys(_property_map_IoK8sApiBatchV1JobSpec))
Swagger.property_type(::Type{ IoK8sApiBatchV1JobSpec }, name::Symbol) = Union{Nothing,eval(Meta.parse(_property_types_IoK8sApiBatchV1JobSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiBatchV1JobSpec }, property_name::Symbol) =  _property_map_IoK8sApiBatchV1JobSpec[property_name]

function check_required(o::IoK8sApiBatchV1JobSpec)
    (getproperty(o, Symbol("template")) === nothing) && (return false)
    true
end

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


text:
export scedctransfer

function df_subset(df::DataFrame,col::String,colsymbol::Symbol)
        col = regex_helper(col)
        ind = findall(occursin.(col,df[!,colsymbol]))
        df = df[ind,:]
	return df
end

function df_subset(df::DataFrame,col::Nothing,colsymbol::Symbol)
	 return df
end

function df_remove(df::DataFrame,col::String,colsymbol::Symbol)
        ind = findall(.!occursin.(col,df[!,colsymbol]))
        df = df[ind,:]
	return df
end

function regex_helper(reg::String)
    if reg == '*'
        # pass for all
    elseif occursin('?',reg)
        reg = replace(reg, '?' => '.')
        reg = Regex(reg)
    elseif occursin('*',reg)
        if reg[end] == '*'
                reg = '^' * strip(reg,'*')
        elseif reg[1] == '*'
                reg = strip(reg,'*') * '$'
        end
            reg = Regex(reg)
    end
    return reg
end

function s3_file_map(filein::String,fileout::String,aws::Dict)
    s3_get_file(aws, "scedc-pds", filein, fileout)
    println("Downloading file: $filein       \r")
end

function filenamecleaner(filename::String)
	filename = strip(filename,'b')
	return replace(filename,"'"=>"")
end

"""

  yearday2Date(year,day)

Convert year and day of year date format to `Date` object.
"""
function yearday2Date(year::String,day::String)
    day = parse(Int,day)
    return Date(year) + Day(day -1)
end

"""

  Date2yearday(d)

Convert `Date` object to yearday string, e.g. 2017354.
"""
function Date2yearday(d::Date)
    days = (d - Date(Year(d))).value + 1
    n = ndigits(days)
    return string(Year(d).value) * ('0' ^ (3 - n)) * string(days)
end

function indexpath(d::Date)
    days = (d - Date(Year(d))).value + 1
    n = ndigits(days)
    outstring = "continuous_waveforms/index/csv/year="
    outstring *= string(Year(d).value) * "/year_doy="
    outstring *= string(Year(d).value) * '_' * ('0' ^ (3 - n)) * string(days)
    outstring *= "/index.csv"
    return outstring
end

function scedcpath(d::Date)
    days = (d - Date(Year(d))).value + 1
    n = ndigits(days)
    outstring = "continuous_waveforms/" * string(Year(d).value) *'/'
    outstring *= string(Year(d).value) * '_' * ('0' ^ (3 - n)) * string(days)
    return outstring
end

"""

  scedctransfer(OUTDIR,startdate)

Transfer data from a single day in the SCEDC open data set from S3 to EC2.

# Arguments
- `OUTDIR::String`: The output directory.
- `startdate::Date`: The start day of the download.
- `enddate::Date`: The (optional) end day of the download.
- `network::String`: Network to download from. If network = "*" or is unspecified,
                       data is downloaded from all available networks.
- `station::String`: Station to download, e.g. "RFO". If station = "*" or is unspecified,
                       data is downloaded from all available stations.
- `channel::String`: Channels to download, e.g. "HH*". If channel = "*" or is unspecified,
                       data is downloaded from all available channels.
- `location::String`: Locations to download, e.g. "00". If channel = "*" or is unspecified,
                       data is downloaded from all available locations. NOTE: most files do
                       not have a location.
- `minlatitude::Float64`: Minimum latitude in data search.
- `maxlatitude::Float64`: Maximum latitude in data search.
- `minlongitude::Float64`: Minimum longitude in data search.
- `maxlongitude::Float64`: Maximum longitude in data search.
"""
function scedctransfer(OUTDIR::String,
                  startdate::Date;
			  	  enddate::Union{Date,Nothing}=nothing,
                  network::Union{String,Nothing}=nothing,
                  station::Union{String,Nothing}=nothing,
                  location::Union{String,Nothing}=nothing,
                  channel::Union{String,Nothing}=nothing,
                  minlatitude::Union{Float64,Nothing}=nothing,
                  maxlatitude::Union{Float64,Nothing}=nothing,
                  minlongitude::Union{Float64,Nothing}=nothing,
                  maxlongitude::Union{Float64,Nothing}=nothing,
                  latitude::Union{Float64,Nothing}=nothing,
                  longitude::Union{Float64,Nothing}=nothing,
                  stationXML::Bool=false)

    tstart = now()
    # connect to S3
    println("Connecting to AWS...      $(now())")
    aws = aws_config(region = "us-west-2")
	@eval @everywhere aws=$aws
    firstdate = Date(2000,1,1)

	if isnothing(enddate)
		enddate = startdate
	end

    # check dates
    if startdate > now()
        @warn("Date must be earlier than today. Aborting download.")
        return nothing
    end

    if enddate < firstdate
        @warn("End date must be later than $firstdate. Aborting download.")
		return nothing
    end

    # get filedf
	date_range = startdate:Day(1):enddate
	paths = indexpath.(date_range)

	println("Starting Download...      $(now())")
	println("Using $(nworkers()) cores...")

	for ii in eachindex(paths)
	    filedf = CSV.read(IOBuffer(s3_get(aws,"scedc-pds",paths[ii])))

	    # subset dataframe
	    # filter by lat/lon
	    if !isnothing(minlatitude)
	        filedf = filedf[filedf[:lat] .> minlatitude,:]
	    end

	    if !isnothing(maxlatitude)
	        filedf = filedf[filedf[:lat] .< maxlatitude,:]
	    end

	    if !isnothing(minlongitude)
	        filedf = filedf[filedf[:lon] .> minlongitude,:]
	    end

	    if !isnothing(maxlongitude)
	        filedf = filedf[filedf[:lon] .< maxlongitude,:]
	    end

	    # filter stations
	    filedf = df_subset(filedf,network,:net)
	    filedf = df_subset(filedf,station,:sta)
	    filedf = df_subset(filedf,channel,:seedchan)
	    filedf = df_subset(filedf,location,:loc)

	    # return if nothing in dataframe
	    if size(filedf,1) == 0
	        @warn("No data available for request $(date_range[ii])!")
			continue
	    end

	    # create directory for instrument responses
	    if stationXML
	        println("Downloading stationXML... $(now())")
	        XMLDIR = joinpath(OUTDIR,"FDSNstationXML")
	        mkpath(XMLDIR)
	        networks = filedf[!,:net]
	        stations = filedf[!,:sta]
	        xmlfiles = networks .* "_" .* stations .* ".xml"
	        ind = indexin(unique(xmlfiles), xmlfiles)
	        networks, stations, xmlfiles = networks[ind], stations[ind], xmlfiles[ind]
	        xml_in = [joinpath("FDSNstationXML",networks[ii],xmlfiles[ii]) for ii = 1:length(xmlfiles)]
	        xml_out = [joinpath(XMLDIR,networks[ii],xmlfiles[ii]) for ii = 1:length(xmlfiles)]
	        xml_dir = unique([dirname(f) for f in xml_out])
	        for ii = 1:length(xml_dir)
	            if !isdir(xml_dir[ii])
	                mkpath(xml_dir[ii])
	            end
	        end

	        # check if requested channels have an instrument response
	        stations2remove = []
	        for ii = 1:length(xml_in)
	            if s3_exists(aws,"scedc-pds",xml_in[ii])
	                s3_get_file(aws,"scedc-pds",xml_in[ii],xml_out[ii])
	            else
	                push!(stations2remove,xmlfiles[ii])
	            end
	        end

	        stations2remove = [split(s[1:end-4],'_')[2] for s in stations2remove]
	        for ii = 1:length(stations2remove)
	            df = df_remove!(filedf,stations2remove[ii],:sta)
	        end
	    end

	    # return if nothing in dataframe
	    if size(filedf,1) == 0
			@warn("No data available for request $(date_range[ii])!")
			continue
	    end

	    # query files
	    dpath = scedcpath(date_range[ii])
	    filenames = filenamecleaner.(filedf[!,:ms_filename])
	    files2download = [joinpath(dpath,f) for f in filenames]

	    # create directory structure
	    OUTDIR = expanduser(OUTDIR)
	    out_files = [joinpath(OUTDIR,dpath,f) for f in filenames]
	    file_dir = unique([dirname(f) for f in out_files])
	    for ii = 1:length(file_dir)
	        if !isdir(file_dir[ii])
	            mkpath(file_dir[ii])
	        end
	    end

	    # download files
	    @eval @everywhere files2download=$files2download
	    @eval @everywhere out_files=$out_files
	    pmap(s3_file_map,files2download,out_files,fill(aws,length(out_files)))
	end
    println("Download Complete!        $(now())          ")
    tend = now()
    println("Download took $(Dates.canonicalize(Dates.CompoundPeriod(tend - tstart)))")
    return nothing
end


text:
module TestCAPI
using Test
using ONNXRunTime.CAPI
import ONNXRunTime as OX

@testset "Session" begin
    api = GetApi()
    env = CreateEnv(api, name="myenv")
    @testset "increment2x3" begin
        path = OX.testdatapath("increment2x3.onnx")
        session_options = CreateSessionOptions(api)
        @test_throws Exception CreateSession(api, env, "does_not_exits.onnx", session_options)
        session = CreateSession(api, env, path, session_options)
        @test SessionGetInputCount(api, session) == 1
        @test SessionGetOutputCount(api, session) == 1
        mem = CreateCpuMemoryInfo(api)
        allocator = CreateAllocator(api, session, mem)
        @test SessionGetInputName(api, session, 0, allocator) == "input"
        @test_throws OX.OrtException SessionGetInputName(api, session, 1, allocator)
        @test SessionGetOutputName(api, session, 0, allocator) == "output"
        @test_throws OX.OrtException SessionGetOutputName(api, session, 1, allocator)
        input_vec = randn(Float32, 6)
        input_array = [
            input_vec[1] input_vec[2] input_vec[3];
            input_vec[4] input_vec[5] input_vec[6];
        ]
        input_tensor = CreateTensorWithDataAsOrtValue(api, mem, input_vec, (2,3))
        run_options = CreateRunOptions(api)
        input_names = ["input"]
        output_names = ["output"]
        inputs = [input_tensor]
        outputs = Run(api, session, run_options, input_names, inputs, output_names)
        @test length(outputs) == 1
        output_tensor = first(outputs)
        @test output_tensor isa OrtValue
        output_array = GetTensorMutableData(api, output_tensor)
        @test typeof(output_array) <: AbstractMatrix{Float32}
        @test size(output_array) == (2,3)
        @test output_array ≈ 1 .+ input_array
    end
end

@testset "tensor roundtrip" begin
    api = GetApi()
    mem = CreateCpuMemoryInfo(api)
    data = randn(2,3)
    tensor = CreateTensorWithDataAsOrtValue(api, mem, vec(OX.reversedims(data)), size(data))
    @test IsTensor(api, tensor)
    info = GetTensorTypeAndShape(api, tensor)
    onnxelty = GetTensorElementType(api, info)
    @test onnxelty isa ONNXTensorElementDataType
    @test OX.juliatype(onnxelty) == eltype(data)
    @test GetDimensionsCount(api, info) == 2
    @test GetDimensions(api, info) == [2,3]
    data2 = GetTensorMutableData(api, tensor)
    @test length(data2) == length(data)
    @test size(data2) == size(data)
    @test eltype(data2) == eltype(data)
    @test data2 == data
end

end


text:
using Test
using LogicCircuits
using ProbabilisticCircuits
using DiscriminativeCircuits
using DataFrames
using CUDA
using Random
using Tables

function test_expectation_brute_force(pc::ProbCircuit, lc::LogisticCircuit, data, CLASSES::Int)
    Random.seed!(124)

    EPS = 1e-4;
    COUNT = size(data)[1]
    # Compute True expectation brute force
    true_exp = zeros(COUNT, CLASSES)
    for i in 1:COUNT
        row = Array(data[i, :])
        cur_data_all = generate_all(row)

        calc_p = log_likelihood_per_instance(pc, cur_data_all)
        calc_p = exp.(calc_p)

        calc_f = class_weights_per_instance(lc, cur_data_all)
        true_exp[i, :] = sum(calc_p .* calc_f, dims=1)
        true_exp[i, :] ./= sum(calc_p) #p_observed
    end
    # Compute Circuit Expect
    calc_exp, cache = Expectation(pc, lc, data);
    for i = 1:COUNT
        for j = 1:CLASSES
            @test true_exp[i,j] ≈ calc_exp[i,j] atol= EPS;
        end
    end
    # Compute Bottom Up Expectation
    calc_exp_2, exp_flow = ExpectationUpward(pc, lc, data);
    for i = 1:COUNT
        for j = 1:CLASSES
            @test true_exp[i,j] ≈ calc_exp_2[i,j] atol= EPS;
        end
    end
    # Compute BitCircuit Expectations (CPU)
    bit_exps, fvalues, gvalues, pbc_cpu = ExpectationBit(pc, lc, data; return_aux=true)
    bit_exps_batch2 = ExpectationBit(pbc_cpu, pc, lc, data, fvalues, gvalues; return_aux=false)
    for i = 1:COUNT
        for j = 1:CLASSES
            @test true_exp[i,j] ≈ bit_exps[i,j] atol= EPS;
            @test true_exp[i,j] ≈ bit_exps_batch2[i,j] atol= EPS;
        end
    end

    # Compute BitCircuit Expectations (GPU)
    if CUDA.functional()
        data_gpu = to_gpu(data);
        bit_exps_gpu, fvalues, gvalues, pbc_gpu = ExpectationBit(pc, lc, data_gpu; return_aux=true)
        bit_exps_gpu_2 = ExpectationBit(pbc_gpu, pc, lc, data_gpu, fvalues, gvalues; return_aux=false)

        # Move to cpu
        bit_exps_cpu = to_cpu(bit_exps_gpu)
        bit_exps_cpu_2 = to_cpu(bit_exps_gpu_2)

        for i = 1:COUNT
            for j = 1:CLASSES
                @test true_exp[i,j] ≈ bit_exps_cpu[i,j] atol= EPS;
                @test true_exp[i,j] ≈ bit_exps_cpu_2[i,j] atol= EPS;
            end
        end
    end

end

function test_moment_brute_force(pc::ProbCircuit, lc::LogisticCircuit, data, CLASSES::Int, moment::Int)
    EPS = 1e-4;
    COUNT = size(data)[1]
    # Compute True moment brute force
    true_mom = zeros(COUNT, CLASSES)
    for i in 1:COUNT
        row = Vector(data[i, :])
        cur_data_all = generate_all(row)

        calc_p = log_likelihood_per_instance(pc, cur_data_all)
        calc_p = exp.(calc_p)

        calc_f = class_weights_per_instance(lc, cur_data_all)
        true_mom[i, :] = sum(calc_p .* (calc_f .^ moment), dims=1)
        true_mom[i, :] ./= sum(calc_p) #p_observed
    end

    # Compute Circuit Moment
    calc_mom, cache = Moment(pc, lc, data, moment);
    for i = 1:COUNT
        for j = 1:CLASSES
            @test (true_mom[i,j] / (calc_mom[i,j] )) ≈ 1.0 atol= EPS;
        end
    end
end

@testset "Expectation Brute Force Test" begin
    #  Small (4 Var)
    psdd_file       = "little_4var.psdd"
    logistic_file   = "little_4var.circuit"
    CLASSES = 2
    N = 4

    pc = zoo_psdd(psdd_file);
    lc = zoo_dc(logistic_file)
    data = DataFrame([0 0 0 0;
            0 1 1 0;
            0 0 1 1;
            missing missing missing missing;
            missing 0 1 missing;
            0 1 missing 1;
            1 missing 0 missing;
            missing 0 1 missing;
            missing missing 0 1;
            missing missing missing 1;
            missing missing missing 0;
            ], :auto)

    data = DataFrame(map(x -> (ismissing(x) ? missing : Bool(x)) , Tables.matrix(data)), :auto)
    test_expectation_brute_force(pc, lc, data, CLASSES)

    # Big circuit (15 Var)
    psdd_file       = "exp-D15-N1000-C4.psdd"
    logistic_file   = "exp-D15-N1000-C4.circuit"
    CLASSES = 4
    N = 15
    COUNT = 10

    pc = zoo_psdd(psdd_file);
    lc = zoo_dc(logistic_file)
    data = DataFrame(rand( (missing,true,false), (COUNT, N) ), :auto)

    test_expectation_brute_force(pc, lc, data, CLASSES)
end

@testset "Moment Brute Force Test" begin
    # Small (4 Var)
    psdd_file       = "little_4var.psdd"
    logistic_file   = "little_4var.circuit";
    CLASSES = 2
    N = 4
    COUNT = 100

    pc = zoo_psdd(psdd_file);
    lc = zoo_dc(logistic_file)
    data = DataFrame(rand( (missing,true,false), (COUNT, N) ), :auto)

    test_moment_brute_force(pc, lc, data, CLASSES, 1)
    test_moment_brute_force(pc, lc, data, CLASSES, 2)
    test_moment_brute_force(pc, lc, data, CLASSES, 3)
    test_moment_brute_force(pc, lc, data, CLASSES, 4)
    test_moment_brute_force(pc, lc, data, CLASSES, 10)
    test_moment_brute_force(pc, lc, data, CLASSES, 15)

    # Big Var
    psdd_file       = "exp-D15-N1000-C4.psdd"
    logistic_file   = "exp-D15-N1000-C4.circuit";
    CLASSES = 4
    N = 15
    COUNT = 10

    pc = zoo_psdd(psdd_file);
    lc = zoo_dc(logistic_file)
    data = DataFrame(rand( (missing,true,false), (COUNT, N) ), :auto)

    test_moment_brute_force(pc, lc, data, CLASSES, 1)
    test_moment_brute_force(pc, lc, data, CLASSES, 2)
    test_moment_brute_force(pc, lc, data, CLASSES, 3)
    test_moment_brute_force(pc, lc, data, CLASSES, 4)
    test_moment_brute_force(pc, lc, data, CLASSES, 10)
    test_moment_brute_force(pc, lc, data, CLASSES, 15)
end

max_stars_repo_name: Juice-jl/DiscriminativeCircuits

max_issues_repo_name: Juice-jl/DiscriminativeCircuits

max_forks_repo_name: Juice-jl/DiscriminativeCircuits


text:
# function ChainRulesCore.rrule(::Interpolations.BSplineInterpolation, T, A, it, axs)
#     y = Interpolations.BSplineInterpolation(T, A, it, axs)
#     function bspline_pb(Δy)
#         return NoTangent(), NoTangent(), Δy, NoTangent(), NoTangent()
#     end
#     return y, bspline_pb
# end TODO: Remove the rrule later.

function ChainRulesCore.rrule(::Type{Interpolations.BSplineInterpolation{T, N, TA, IT, TAX}}, A, axs, it) where {T, N, TA, IT, TAX}
    y = Interpolations.BSplineInterpolation{T, N, TA, IT, TAX}(A, axs, it)
    function bspline_const_pb(Δy)
        return NoTangent(), Δy, NoTangent(), NoTangent()
    end
    return y, bspline_const_pb
end

function ChainRulesCore.rrule(::Type{Interpolations.Extrapolation{T, N, ITPT, IT, ET}}, itp, et::Union{Interpolations.Flat, Interpolations.Reflect, Interpolations.Line}) where {T,N,ITPT,IT,ET}
    y = Interpolations.Extrapolation{T,N,ITPT,IT,ET}(itp, et)
    function extrap_const_pb(Δy)
        return NoTangent(), Δy, NoTangent()
    end
    return y, extrap_const_pb
end

Zygote.@nograd Interpolations.tweight
Zygote.@nograd Interpolations.tcoef

function ChainRulesCore.rrule(::typeof(Interpolations.copy_with_padding), T, A, it)
    y = Interpolations.copy_with_padding(T, A, it)
    function copy_with_padding_pullback(Δy)
        return NoTangent(), NoTangent(), Δy, NoTangent()
    end
    return y, copy_with_padding_pullback
end

function ChainRulesCore.rrule(::Type{Interpolations.FilledExtrapolation}, itp, fv)
    y = Interpolations.FilledExtrapolation(itp, fv)
    function filledextra_const_pb(Δy)
        return NoTangent(), Δy, NoTangent()
    end
    return y, filledextra_const_pb
end

function ChainRulesCore.rrule(::Type{SArray{D, T, ND, L}}, x...) where {D, T, ND, L}
    y = SArray{D, T, ND, L}(x...)
    function sarray_pb(Δy)
        Δy = map(t->eltype(x...)(t...), Δy)
        return NoTangent(), Δy
    end
    return y, sarray_pb
end

function ChainRulesCore.rrule(p::Type{RotMatrix{N, T, L}}, x) where {N, T, L}
    y = p(x)
    function rotmat_pb(Δy)
        Δy = p(Δy)
        return NoTangent(), Δy
    end
    return y, rotmat_pb
end

function ChainRulesCore.rrule(::typeof(ImageTransformations._getindex), A::AbstractExtrapolation, x)
    y = A(Tuple(x)...)
    function _getindex_pb(Δy)
        # Δy :: NamedTuple{(:r, :g, :b)} or something similar
        Δy = eltype(A)(Δy...)
        gr = Interpolations.gradient(A, x...)
        n = nfields(Δy) > 0 ? nfields(Δy) : 1
        return NoTangent(), NoTangent(), (Δy .⋅ gr) * n
    end
    return y, _getindex_pb
end

function ChainRulesCore.rrule_via_ad(::Zygote.ZygoteRuleConfig, h::DiffImages.Homography, x::SVector{2, K}) where K
    Zygote.pullback((x, y)->x(y), h, x)
end

function ChainRulesCore.rrule(::typeof(ImageTransformations.warp!), out, img::AbstractExtrapolation, tform)
    out = ImageTransformations.warp!(out, img, tform)
    function warp!_pb(Δy)
        ∇out = NoTangent()
        ∇img = @not_implemented("To be implemented.")
        flds = fieldnames(typeof(tform))
        ∇tform = canonicalize(Tangent{typeof(tform)}())
        Δy = collect(Δy)
        for p in CartesianIndices(out)
            _, ∇τ = rrule(ImageTransformations._getindex, img, p.I)
            _, _, ∇τ = ∇τ(Δy[p])
            _, ∇ϕ = rrule_via_ad(Zygote.ZygoteRuleConfig(), tform, SVector(p.I))
            ∇h, _ = ∇ϕ(∇τ)
            ∇tform += Tangent{typeof(tform)}(;NamedTuple{flds}(∇h)...)
        end
        return NoTangent(), ∇out, ∇img, ∇tform
    end
    return out, warp!_pb
end


text:
using Test
using PilePenetration

using TOML
using CSV

using Serialization
using NaturalSort

@testset "Short pile penetration" begin
    PilePenetration.main_simulation("input.toml"); println()
    output_dir = TOML.parsefile("input.toml")["General"]["output_folder_name"]

    # check results
    output = CSV.File("output.csv") # expected output
    history = CSV.File(joinpath(output_dir, "history.csv"))
    for name in propertynames(output)
        output_col = output[name]
        history_col = history[name]
        @test history_col ≈ output_col atol=1e-3
    end

    # check input.toml
    @test TOML.parsefile("input.toml") == TOML.parsefile(joinpath(output_dir, "input.toml"))

    # check paraview files
    @test all(i -> isfile(joinpath(output_dir, string("paraview/output", i, ".vtm"))), 0:10)
    # check snapshots file
    root, _, files = only(walkdir(joinpath(output_dir, "snapshots")))
    count = 0
    for file in sort(files; lt = natural)
        @test file == "snapshot$count"
        data = deserialize(joinpath(root, file))
        @test data isa NamedTuple
        @test keys(data) == (:grid, :pointstate, :rigidbody, :rigidbody0, :t)
        count += 1
    end
end

@testset "Restart simulation" begin
    PilePenetration.main_simulation("input_restart.toml"); println()
    dict = TOML.parsefile("input_restart.toml")
    dir_ext = splitext(dict["General"]["output_folder_name"])
    output_dir = string(dir_ext[1], "_restarted_from_", dict["General"]["restart"], dir_ext[2])

    # check results
    output = CSV.File("output.csv") # expected output
    history = CSV.File(joinpath(output_dir, "history.csv"))
    for name in propertynames(output)
        output_last = output[name][end]
        history_last = history[name][end]
        @test history_last ≈ output_last atol=1e-3
    end
end

@testset "Post-processing" begin
    @testset "normal" begin
        output_dir = "output.tmp"
        inputfile = joinpath(output_dir, "postprocess.toml")
        cp("postprocess.toml", inputfile; force = true)
        PilePenetration.main_postprocess(inputfile)

        # check results
        output = CSV.File("output.csv") # expected output
        history = CSV.File(joinpath(output_dir, "postprocess.tmp", "history.csv"))
        for name in propertynames(output)
            output_col = output[name]
            history_col = history[name]
            @test history_col ≈ output_col atol=1e-3
        end
    end
    @testset "restart" begin
        output_dir = "output_restarted_from_5.tmp"
        inputfile = joinpath(output_dir, "postprocess.toml")
        cp("postprocess.toml", inputfile; force = true)
        PilePenetration.main_postprocess(inputfile)

        # check results
        output = CSV.File("output.csv") # expected output
        history = CSV.File(joinpath(output_dir, "postprocess.tmp", "history.csv"))
        for name in propertynames(output)
            output_last = output[name][end]
            history_last = history[name][end]
            @test history_last ≈ output_last atol=1e-3
        end
    end
end


text:
# # Fractal Growth

# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../fractal.mp4" type="video/mp4">
# </video>
# ```

# This model follows the process known as diffusion-limited aggregation to simulate the growth of fractals.
# It is a kinetic process that consists of randomly diffusing particles giving rise to fractal-like structures
# resembling those observed naturally. This examplet is based off of
# ["Particularly Stuck" example](https://www.complexity-explorables.org/explorables/particularly-stuck/)
# in Complexity Explorables.

# The environment is a two dimensional, continuous space world. Agents are particles that diffuse and aggregate to form
# fractals. Initially, there are particles of random size distributed across the space, and one static particle in the center
# that forms the seed for the fractal growth. As moving particles collide with the seed or any particle that previously
# collided with the seed, it gets stuck and contributes to the fractal. As a particle gets stuck, another one is created at a
# circular border around the center to feed the growth.

# It is also available from the `Models` module as [`Models.fractal_growth`](@ref).
using Agents, LinearAlgebra
using Random # hide

# We use the [`@agent`](@ref) macro to conveniently define a `Particle` agent. Each agent
# has a radius, representing the particle size, a boolean to define whether it is stuck and part of the fractal,
# and an axis around which it spins (elaborated on later). In addition, since we use the [`ContinuousAgent`](@ref)
# type, the [`@agent`](@ref) macro also provides each agent with fields for `id`, `pos` (its position in space) and
# `vel` (its velocity).
@agent Particle ContinuousAgent{2} begin
    radius::Float64
    is_stuck::Bool
    spin_axis::Array{Float64,1}
end

# A custom constructor allows convenient creation of agents.
Particle(
    id::Int,
    radius::Float64,
    spin_clockwise::Bool;
    pos = (0.0, 0.0),
    is_stuck = false,
) = Particle(id, pos, (0.0, 0.0), radius, is_stuck, [0.0, 0.0, spin_clockwise ? -1.0 : 1.0])

# We also define a few utility functions for ease of implementation.
# `rand_circle` returns a random point on the unit circle. `particle_radius`
# generates a random radius for a particle, within given range defined by `min_radius`
# and `max_radius`. If `max_radius < min_radius`, it returns `min_radius`: allowing
# a fixed particle size to be specified.
rand_circle(rng) = (θ = rand(rng, 0.0:0.1:359.9); (cos(θ), sin(θ)))
particle_radius(min_radius::Float64, max_radius::Float64, rng) =
    min_radius <= max_radius ? rand(rng, min_radius:0.01:max_radius) : min_radius

# The `initialize_model` function returns a new model containing particles placed
# randomly in the given space and one seed particle at the center.
function initialize_model(;
    initial_particles::Int = 100, # initial particles in the model, not including the seed
    ## size of the space in which particles exist
    space_extents::NTuple{2,Float64} = (150.0, 150.0),
    speed = 0.5, # speed of particle movement
    vibration = 0.55, # amplitude of particle vibration
    attraction = 0.45, # velocity of particles towards the center
    spin = 0.55, # tangential velocity with which particles orbit the center
    ## fraction of particles orbiting clockwise. The rest are anticlockwise
    clockwise_fraction = 0.0,
    min_radius = 1.0, # minimum radius of any particle
    max_radius = 2.0, # maximum radius of any particle
)
    properties = Dict(
        :speed => speed,
        :vibration => vibration,
        :attraction => attraction,
        :spin => spin,
        :clockwise_fraction => clockwise_fraction,
        :min_radius => min_radius,
        :max_radius => max_radius,
        :spawn_count => 0,
    )
    ## space is periodic to allow particles going off one edge to wrap around to the opposite
    space = ContinuousSpace(space_extents, 1.0; periodic = true)
    model = ABM(Particle, space; properties)
    center = space_extents ./ 2.0
    for i in 1:initial_particles
        particle = Particle(
            i,
            particle_radius(min_radius, max_radius, model.rng),
            rand(model.rng) < clockwise_fraction,
        )
        ## `add_agent!` automatically gives the particle a random position in the space
        add_agent!(particle, model)
    end
    ## create the seed particle
    particle = Particle(
        initial_particles + 1,
        particle_radius(min_radius, max_radius, model.rng),
        true;
        pos = center,
        is_stuck = true,
    )
    ## `add_agent_pos!` will use the position of the agent passed in, instead of assigning it
    ## to a random value
    add_agent_pos!(particle, model)
    return model
end

# The `agent_step!` function simulates particle motion for those who are not yet `stuck`.
# For each particle, we first perform a crude distance check to all other particles.
# If the current particle intersects any particle in the fractal, it also becomes
# part of the fractal and is not simulated further. Agent velocity has a radial component that
# attracts it towards the center, a tangential component that makes it orbit around the center,
# and a random component that simulates vibration of the particle. The velocity is scaled
# to be inversely proportional to the square of the particle's radius, so that larger particles
# move slower. The `speed` parameter is implemented as the time difference between successive
# steps of the simulation. A larger value causes particles to move more per step, but leads to
# inaccuracies as particles do not move through the intervening space.
function agent_step!(agent::Particle, model)
    agent.is_stuck && return

    for id in nearby_ids(agent.pos, model, agent.radius)
        if model[id].is_stuck
            agent.is_stuck = true
            ## increment count to make sure another particle is spawned as this one gets stuck
            model.spawn_count += 1
            return
        end
    end
    ## radial vector towards the center of the space
    radial = model.space.extent ./ 2.0 .- agent.pos
    radial = radial ./ norm(radial)
    ## tangential vector in the direction of orbit of the particle
    tangent = Tuple(cross([radial..., 0.0], agent.spin_axis)[1:2])
    agent.vel =
        (
            radial .* model.attraction .+ tangent .* model.spin .+
            rand_circle(model.rng) .* model.vibration
        ) ./ (agent.radius^2.0)
    move_agent!(agent, model, model.speed)
end

# The `model_step!` function serves the sole purpose of spawning additional particles
# as they get stuck to the growing fractal.
function model_step!(model)
    while model.spawn_count > 0
        particle = Particle(
            nextid(model),
            particle_radius(model.min_radius, model.max_radius, model.rng),
            rand(model.rng) < model.clockwise_fraction;
            pos = (rand_circle(model.rng) .+ 1.0) .* model.space.extent .* 0.49,
        )
        add_agent_pos!(particle, model)
        model.spawn_count -= 1
    end
end

# # Running the model

# We run the model using the `InteractiveDynamics` package with `GLMakie` backend so
# the fractal growth can be visualised as it happens. `InteractiveDynamics` provides the `abmvideo` function to easily record a video of the simulation running.
Random.seed!(42) # hide
model = initialize_model()

using InteractiveDynamics
import CairoMakie
CairoMakie.activate!() # hide

# Particles that are stuck and part of the fractal are shown in red, for visual distinction
particle_color(a::Particle) = a.is_stuck ? :red : :blue
# The visual size of particles corresponds to their radius, and has been calculated
# for the default value of `space_extents` of the `initialize_model` function. It will
# not look accurate on other values.
particle_size(a::Particle) = 7.5 * a.radius

abmvideo(
    "fractal.mp4",
    model,
    agent_step!,
    model_step!;
    ac = particle_color,
    as = particle_size,
    am = '●',
    spf = 10,
    frames = 600,
    framerate = 25,
    title = "Fractal Growth",
    scatterkwargs = (strokewidth = 0.5, strokecolor = :white),
)
nothing # hide
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="../fractal.mp4" type="video/mp4">
# </video>
# ```

# Using `InteractiveDynamics` simulation parameters can also be tweaked
# dynamically. This makes use of the [`InteractiveDynamics.abm_data_exploration`](@ref) function.

# ```julia
# using InteractiveDynamics
# using GLMakie
# model = initialize_model()
# ```

# `params` defines the range in which different parameter values can be adjusted through
# sliders.
# ```julia
# params = (
#     :attraction => 0.0:0.01:2.0,
#     :speed => 0.0:0.01:2.0,
#     :vibration => 0.0:0.01:2.0,
#     :spin => 0.0:0.01:2.0,
#     :clockwise_fraction => 0.0:0.01:1.0,
#     :min_radius => 0.5:0.01:3.0,
#     :max_radius => 0.5:0.01:3.0,
# )
#
# particle_size(a::Particle) = 4 * a.radius
# abm_data_exploration(
#     model,
#     agent_step!,
#     model_step!,
#     params;
#     ac = particle_color,
#     as = particle_size,
#     am = '⚪',
# )
# ```
#
# ```@raw html
# <video width="auto" controls autoplay loop>
# <source src="https://raw.githubusercontent.com/JuliaDynamics/JuliaDynamics/master/videos/agents/fractal_interact.mp4?raw=true" type="video/mp4">
# </video>
# ```


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

using Markdown
using InteractiveUtils

# ╔═╡ cb402178-8d6c-11eb-0162-0929b93a7cdb
sqrt(2)

# ╔═╡ 1b81c7a4-8d6d-11eb-1fff-050c4510b871
f = sqrt∘sum

# ╔═╡ 3652cce0-8d6d-11eb-204f-6d4f9f3e16ff
f([1,2])

# ╔═╡ 3ce8e698-8d6d-11eb-1666-5d62cbf21d1b
sum([1,2,3]) |> abs |> sqrt

# ╔═╡ 4cae5ade-8d6d-11eb-0169-7154cca2b874
([1,2,3])<| sum

# ╔═╡ 8b124a94-8d6d-11eb-0bb1-9f777e1d98f6
map(x -> x +1 ,[1,2,3]) == [2,3,4]

# ╔═╡ 24c947fa-8d6e-11eb-1b1b-dfdea765d7e6
foreach(println,[1,2,3])

# ╔═╡ 571228d8-8d6e-11eb-381f-4bb319df187f
foldl(=>,1:4) == (1 => 2 => 3 => 4)

# ╔═╡ 9943f8f0-8d6e-11eb-1bd0-4f90bdfb4543
foldl(=>,1:4)

# ╔═╡ 8e198494-8d6f-11eb-012d-034e88313e4e
begin
	a = 1
	b = 2
end

# ╔═╡ faa90d6c-8d6f-11eb-06b1-530a08ed4846
⊗(a,b) = kron(a,b)

# ╔═╡ 168e6272-8d70-11eb-2f8b-4dabdadaa454
⊗([1,2],[2,3])

# ╔═╡ 28272208-8d70-11eb-12ad-8fc5c7f16c98
foldl(⊗,[[1,2],[1,2],[1,2]])

# ╔═╡ 7c054e56-8d70-11eb-043d-cfebcc74b872
g(x) = (x==1||x==2) ? "x=1" : "x=2"

# ╔═╡ 62a9c638-8d72-11eb-253a-cb0f833dc233
g(2)

# ╔═╡ 9b18de58-8d77-11eb-335b-4b1775144c05

# ╔═╡ ae4013f2-8d77-11eb-13f6-3fd35a98cf09

# ╔═╡ Cell order:
# ╠═cb402178-8d6c-11eb-0162-0929b93a7cdb
# ╠═1b81c7a4-8d6d-11eb-1fff-050c4510b871
# ╠═3652cce0-8d6d-11eb-204f-6d4f9f3e16ff
# ╠═3ce8e698-8d6d-11eb-1666-5d62cbf21d1b
# ╠═4cae5ade-8d6d-11eb-0169-7154cca2b874
# ╠═8b124a94-8d6d-11eb-0bb1-9f777e1d98f6
# ╠═24c947fa-8d6e-11eb-1b1b-dfdea765d7e6
# ╠═571228d8-8d6e-11eb-381f-4bb319df187f
# ╠═9943f8f0-8d6e-11eb-1bd0-4f90bdfb4543
# ╠═8e198494-8d6f-11eb-012d-034e88313e4e
# ╠═faa90d6c-8d6f-11eb-06b1-530a08ed4846
# ╠═168e6272-8d70-11eb-2f8b-4dabdadaa454
# ╠═28272208-8d70-11eb-12ad-8fc5c7f16c98
# ╠═7c054e56-8d70-11eb-043d-cfebcc74b872
# ╠═62a9c638-8d72-11eb-253a-cb0f833dc233
# ╠═9b18de58-8d77-11eb-335b-4b1775144c05
# ╠═ae4013f2-8d77-11eb-13f6-3fd35a98cf09

max_stars_repo_name: Fromeworld/JuliaSeminar

max_issues_repo_name: Fromeworld/JuliaSeminar

max_forks_repo_name: Fromeworld/JuliaSeminar