text:
using SpecialFunctions

@testset "SpecialFunctions" for x in (1, -1, 0, 0.5, 10, -17.1, 1.5 + 0.7im)
    test_scalar(SpecialFunctions.erf, x)
    test_scalar(SpecialFunctions.erfc, x)
    test_scalar(SpecialFunctions.erfi, x)

    test_scalar(SpecialFunctions.airyai, x)
    test_scalar(SpecialFunctions.airyaiprime, x)
    test_scalar(SpecialFunctions.airybi, x)
    test_scalar(SpecialFunctions.airybiprime, x)

    test_scalar(SpecialFunctions.besselj0, x)
    test_scalar(SpecialFunctions.besselj1, x)

    test_scalar(SpecialFunctions.erfcx, x)
    test_scalar(SpecialFunctions.dawson, x)

    if x isa Real
        test_scalar(SpecialFunctions.invdigamma, x)
    end

    if x isa Real && 0 < x < 1
        test_scalar(SpecialFunctions.erfinv, x)
        test_scalar(SpecialFunctions.erfcinv, x)
    end

    if x isa Real && x > 0 || x isa Complex
        test_scalar(SpecialFunctions.bessely0, x)
        test_scalar(SpecialFunctions.bessely1, x)
        test_scalar(SpecialFunctions.gamma, x)
        test_scalar(SpecialFunctions.digamma, x)
        test_scalar(SpecialFunctions.trigamma, x)
        test_scalar(SpecialFunctions.lgamma, x)
    end
end


text:
#md # # Tutorial
#md # Notebook [download](tutorial.ipynb) [nbviewer](@__NBVIEWER_ROOT_URL__/tutorial.ipynb) [source](tutorial.jl)

using MGVI
#
using Distributions
using Random
using ValueShapes
using LinearAlgebra
using Optim

# We want to fit a 3-degree polynomial using two data sets (`a` and `b`). MGVI requires a model
# expressed as a function of the model parameters and returning an instance of the Distribution.
# In this example, since we have two sets of independent measurements, we express them as
# ValueShapes.NamedTupleDist.
#
#
# We assume errors are normally distributed with unknown covariance, which has to be learned as well.

const _x1_grid = [Float64(i)/10 for i in 1:25]
const _x2_grid = [Float64(i)/10 + 0.1 for i in 1:15]
const _common_grid = sort(vcat(_x1_grid, _x2_grid))

function _mean(x_grid, p)
    p[1]*10 .+ p[2]*40 .* x_grid .+ p[3]*600 .* x_grid.^2 .+ p[4]*80 .* x_grid.^3
end

function model(p)
    dist1 = Product(Normal.(_mean(_x1_grid, p), p[5]^2*60))
    dist2 = Product(Normal.(_mean(_x2_grid, p), p[5]^2*60))
    NamedTupleDist(a=dist1,
                   b=dist2)
end
# Here we define the ground truth of the parameters, as well as an initial guess.

const true_params =  [
 -0.3
 -1.5
 0.2
 -0.5
 0.3]

const starting_point = [
  0.2
  0.5
  -0.1
  0.3
 -0.6
];
#-
function pprintln(obj)
    show(stdout, "text/plain", obj)
    println()
end;
#-
using Plots
gr(size=(400, 300), dpi=700, fmt=:png)
#-
rng = MersenneTwister(157);

# We draw data directly from the model, using the true parameter values:

data = rand(rng, model(true_params), 1)[1];
#-
function _mean(x::Vector)
    _mean(_common_grid, x)
end

init_plots =() -> let
    truth = _mean(true_params)
    plot!(_common_grid, truth, markercolor=:blue, linecolor=:blue, label="truth")
    scatter!(_common_grid, _mean(starting_point), markercolor=:orange, markerstrokewidth=0, markersize=3, label="init")
    scatter!(vcat(_x1_grid, _x2_grid), reduce(vcat, data), markercolor=:black, markerstrokewidth=0, markersize=3, label="data")
end;

# Before we start the optimization, let's have an initial look at the data.
# It is also interesting to see how our starting guess performs.

p = plot()
init_plots()

# Now we are ready to run one iteration of the MGVI.
# The output contains an updated parameter estimate (`first_iteration.result`),
# which we can compare to the true parameters.

first_iteration = mgvi_kl_optimize_step(rng,
                                        model, data,
                                        starting_point;
                                        jacobian_func=FwdRevADJacobianFunc,
                                        residual_sampler=ImplicitResidualSampler,
                                        optim_options=Optim.Options(iterations=10, show_trace=true),
                                        residual_sampler_options=(;cg_params=(;maxiter=10)))
pprintln(hcat(first_iteration.result, true_params))
p
#-
plot_iteration = (params, label) -> let
    #error_mat = mgvi_kl_errors(full_model, params)
    #display(error_mat)
    #errors = sqrt.(error_mat[diagind(error_mat)])
    #yerr = abs.(line(common_grid, params+errors) - line(common_grid, params-errors))
    #scatter!(common_grid, line(common_grid, params), markercolor=:green, label=label, yerr=yerr)
    for sample in eachcol(params.samples)
        scatter!(_common_grid, _mean(Vector(sample)), markercolor=:gray, markeralpha=0.3, markersize=2, label=nothing)
    end
    scatter!(_common_grid, _mean(params.result), markercolor=:green, label=label)
end;

# Now let's also plot the curve corresponding to the new parameters after the first iteration:

p = plot()
init_plots()
plot_iteration(first_iteration, "first")
p
#-
plot_iteration_light = (params, counter) -> let
    scatter!(_common_grid, _mean(params.result), markercolor=:green, markersize=3, markeralpha=2*atan(counter/18)/π, label=nothing)
end;

# From the plot above we see that one iteration is not enough. Let's do 5 more steps and plot the evolution of estimates.

init_plots()
plt = scatter()
next_iteration = first_iteration
for i in 1:5
    pprintln(minimum(next_iteration.optimized))
    pprintln(hcat(next_iteration.result, true_params))
    global next_iteration = mgvi_kl_optimize_step(rng,
                                                  model, data,
                                                  next_iteration.result;
                                                  jacobian_func=FwdRevADJacobianFunc,
                                                  residual_sampler=ImplicitResidualSampler,
                                                  optim_options=Optim.Options(iterations=10, show_trace=true),
                                                  residual_sampler_options=(;cg_params=(;maxiter=10)))
    plot_iteration_light(next_iteration, i)
end
pprintln(minimum(next_iteration.optimized))
pprintln(hcat(next_iteration.result, true_params))
plt

# Finally, let's plot the last estimate and compare it to the truth. Also, notice, that gray dots represent samples from
# the approximation.

p = plot()
init_plots()
plot_iteration(next_iteration, "last")
p

ext: jl


text:
module Example

if !isdefined(Core, :String)
    typealias String UTF8String
end

using FlatBuffers

@enum(Color, Red = 1, Green = 2, Blue = 8)
@default Color Red

@struct immutable Test
    a::Int16
    b::UInt8
    # _1::UInt8 # padding
end

type TestSimpleTableWithEnum
    color::Color
end

@default TestSimpleTableWithEnum color=Green

@struct immutable Vec3
    x::Float32
    y::Float32
    z::Float32
    # _::UInt32 # padding
    test1::Float64
    test2::Color
    # __::UInt8 # padding
    test3::Test
    # ___::UInt16 # padding
end

@align Vec3 16

type Stat
    id::String
    val::Int64
    count::UInt16
end

# Julia doesn't support forward referencing of types
# @union Any_ Union{Monster, TestSimpleTableWithEnum}

type Monster
    pos::Vec3
    mana::Int16
    hp::Int16
    name::String
    friendly::Bool # deprecated
    inventory::Vector{UInt8}
    color::Color
    # test_type::Any_
    # test::Vector{UInt8}
    test4::Vector{Test}
    testarrayofstring::Vector{String}
    testarrayoftables::Vector{Monster}
    # don't support nested circulr reference objects yet
    # enemy::Monster
    testnestedflatbuffer::Vector{UInt8}
    testempty::Stat
    testbool::Bool
    testhashs32_fnv1::Int32
    testhashu32_fnv1::UInt32
    testhashs64_fnv1::Int64
    testhashu64_fnv1::UInt64
    testhashs32_fnv1a::Int32
    testhashu32_fnv1a::UInt32
    testhashs64_fnv1a::Int64
    testhashu64_fnv1a::UInt64
    testarrayofbools::Vector{Bool}
    testf::Float32
    testf2::Float32
    testf3::Float32
end

@default Monster hp=100 mana=150 color=Blue friendly=false testf=Float32(3.14159) testf2=Float32(3)

end # module


text:
function DFSVisit(G::DictGraph, dfsDict::SortedDict, vertex::Int, count::Int)
    count  += 1;
    dfsDict[vertex][2] = t;
    dfsDict[vertex][1] = 0;

    for v in G[vertex]
        if dfsDict[v[1]][1] == -1
            dfsDict[v[1]][3] = vertex
            DFSVisit(G, dfsDict, v[1], count);
        end
    end

    dfsDict[vertex][1] = 1
    count += 1
    dfsDict[vertex][4] += count

    return;
end

function DFS(G::DictGraph; vertex::Int = -1)
    dfsDict = SortedDict(i => [-1 0 -1 0] for i in keys(G.relation))
    #                         [c  d  π f]
    count = -1

    if vertex > -1
        vertex in keys(G) || error("$vertex is not a vertex of this graph.")
        DFSVisit(G, dfsDict, vertex, count)
    end

    for v in keys(G.relation)
        dfsDict[v][1] == -1 && DFSVisit(G, dfsDict, v, count)
    end

    return dfsDict
end


text:
module NetworkModels

using Graphs
using Distributions
using Pipe
using LinearAlgebra

end


text:
module PilePenetration

using PoingrSimulator
using PoingrSimulator.Poingr
using PoingrSimulator.GeometricObjects
using PoingrSimulator.GeometricObjects: getline

using TOML

using Base: @_propagate_inbounds_meta, @_inline_meta

include("injection.jl")
include("postprocess.jl")

function main_simulation()::Cint
    if isempty(ARGS)
        inputtoml = "input.toml"
    else
        inputtoml = ARGS[1]
    end
    try
        main_simulation(inputtoml)
    catch
        Base.invokelatest(Base.display_error, Base.catch_stack())
        return 1
    end
    return 0
end

function main_simulation(inputtoml_file::AbstractString)
    proj_dir = dirname(inputtoml_file)
    input = TOML.parsefile(inputtoml_file)

    General  = input["General"]
    Advanced = input["Advanced"]
    Pile     = input["Pile"]

    inputtoml = Dict{String, Any}(
        "General" => Dict{String, Any}(
            "type"              => "PenetrateIntoGround",
            "coordinate_system" => "axisymmetric",
            "domain"            => General["domain"],        # [[xmin, xmax], [ymin, ymax]] (m)
            "grid_space"        => General["grid_space"],    # (m)
            "gravity"           => General["gravity"],       # (m/s2)
            "total_time"        => General["total_time"],    # (sec)
            "show_progress"     => General["show_progress"],
        ),
        "Output" => Dict{String, Any}(
            "interval"       => General["output_interval"],   # (sec)
            "directory"      => General["output_folder_name"],
            "history"        => false,
            "snapshots"      => true,
            "paraview"       => true,
            "paraview_grid"  => false,
            "copy_inputfile" => false,
        ),
        "Advanced" => Dict{String, Any}(
            "CFL"                       => Advanced["CFL"],
            "contact_threshold_scale"   => Advanced["contact_threshold_scale"],
            "contact_penalty_parameter" => Advanced["contact_penalty_parameter"],
            "npoints_in_cell"           => 2, # in each dimension
            "reorder_pointstate"        => false, # don't reorder to prevent error
        ),
    )
    if haskey(General, "restart")
        inputtoml["General"]["restart"] = General["restart"]
    end

    inputtoml["SoilLayer"] = map(input["SoilLayer"]) do SoilLayer
        friction_with_pile_inside  = SoilLayer["friction_with_pile_inside"]
        friction_with_pile_outside = SoilLayer["friction_with_pile_outside"]
        Dict{String, Any}(
            "mohr_coulomb_type"       => "circumscribed",
            "thickness"               => SoilLayer["thickness"],                        # (m)
            "density"                 => SoilLayer["unit_weight"] / General["gravity"], # (kg/m3)
            "poissons_ratio"          => SoilLayer["poissons_ratio"],
            "youngs_modulus"          => SoilLayer["youngs_modulus"],                   # (N/m2)
            "cohesion"                => SoilLayer["cohesion"],                         # (N/m2)
            "friction_angle"          => SoilLayer["friction_angle"],                   # (degree)
            "dilatancy_angle"         => SoilLayer["dilatancy_angle"],                  # (degree)
            "tension_cutoff"          => SoilLayer["tension_cutoff"],
            "friction_with_rigidbodies" => [[friction_with_pile_inside, friction_with_pile_inside,
                                             friction_with_pile_outside, friction_with_pile_outside,
                                             friction_with_pile_outside, friction_with_pile_outside]],
        )
    end

    inputtoml["RigidBody"] = begin
        radius_head    = Pile["diameter_head"] / 2 # (m)
        radius_tip     = Pile["diameter_tip"] / 2  # (m)
        pile_length    = Pile["pile_length"]       # (m)
        tapered_length = Pile["tapered_length"]    # (m)
        thickness      = Pile["thickness"]         # (m) NOTE: This must be at least `2 * grid_space`
        Any[
            Dict{String, Any}(
                "type"        => "Polygon",
                "velocity"    => [0.0, -Pile["velocity"]],  # (m/s)
                "coordinates" => [[radius_head, pile_length], [radius_head, tapered_length], [radius_tip, 0.0],
                                  [radius_tip+thickness, 0.0], [radius_head+thickness, tapered_length], [radius_head+thickness, pile_length]]
            )
        ]
    end

    inputtoml["Pile"] = input["Pile"]
    if haskey(input, "OutputHistory")
        inputtoml["OutputHistory"] = input["OutputHistory"]
    end

    if haskey(inputtoml["General"], "restart")
        dir_ext = splitext(inputtoml["Output"]["directory"])
        output_dir = joinpath(proj_dir, string(dir_ext[1], "_restarted_from_", inputtoml["General"]["restart"], dir_ext[2]))
    else
        output_dir = joinpath(proj_dir, inputtoml["Output"]["directory"])
    end
    mkpath(output_dir)
    cp(inputtoml_file, joinpath(output_dir, "input.toml"); force = true)

    io = IOBuffer()
    TOML.print(io, inputtoml)
    PoingrSimulator.main(proj_dir, String(take!(io)), Injection)
end

end # module


text:
## Run MMALA using a Bayesian probit model with a Normal prior N(0, priorVar*I)
## on the Swiss banknote data

using Distributions, GeometricMCMC

include("vaso.jl")
include("probitNormalPrior.jl")

# Create Model instance
model = Model(nPars, data, logPrior, logLikelihood, gradLogPosterior,
  tensor, derivTensor, randPrior);

# Create setDriftStep function for adjusting the drift step of MMALA
driftSteps = [1, 1e-8, 1e-7, 1e-6, 1e-5, 1e-4, 1e-3, 1e-2, 1e-1, 0.25]

setDriftStep(i::Int, acceptanceRatio::Float64, nMcmc::Int,
  nBurnin::Int, currentStep::Float64) = setMmalaDriftStep(i, acceptanceRatio,
  nMcmc, nBurnin, currentStep, driftSteps)

# Create instance malaOpts of MMALA options
mmalaOpts = MmalaOpts(55000, 5000, setDriftStep);

# Run MMALA
mcmc = Array(Float64, nPars, nPars)
z = Array(Float64, nPars, nPars)
linearZvMcmc = Array(Float64, nPars, nPars)
linearCoef = Array(Float64, nPars, nPars)
quadraticZvMcmc = Array(Float64, nPars, nPars)
quadraticCoef = Array(Float64, convert(Int, nPars*(nPars+3)/2), nPars)

try
  mcmc, z = mmala(model, mmalaOpts);

  # Compute ZV-MMALA mean estimators based on linear polynomial
  linearZvMcmc, linearCoef = linearZv(mcmc, z);

  # Compute ZV-MMALA mean estimators based on quadratic polynomial
  quadraticZvMcmc, quadraticCoef = quadraticZv(mcmc, z);
catch msg
  println(msg)
end


text:
using CUDAapi

using Test

import Libdl

@testset "library types" begin
    @test CUDAapi.PATCH_LEVEL == CUDAapi.libraryPropertyType(2)
    @test CUDAapi.C_32U == CUDAapi.cudaDataType(Complex{UInt32})
end

@testset "properties" begin
    CUDAapi.devices_for_cuda(v"8.0")
    CUDAapi.devices_for_llvm(v"5.0")
    CUDAapi.isas_for_cuda(v"8.0")
    CUDAapi.isas_for_llvm(v"5.0")
end

# helper macro to test for non-nothingness
macro test_something(ex)
    quote
        local rv = $(esc(ex))
        @test rv !== nothing
        rv
    end
end

@testset "discovery" begin
    CUDAapi.find_binary([Sys.iswindows() ? "CHKDSK" : "true"])
    CUDAapi.find_library([Sys.iswindows() ? "NTDLL" : "c"])

    dirs = find_toolkit()
    @test !isempty(dirs)

    ver = find_toolkit_version(dirs)

    @testset "CUDA tools and libraries" begin
        @test_something find_cuda_binary("nvcc", dirs)
        @test_something find_cuda_library("cudart", dirs)
        @test_something find_cuda_library("nvtx", dirs)
        @test_something find_libdevice([v"3.0"], dirs)
        @test_something find_libcudadevrt(dirs)
    end

    if haskey(ENV, "CI")
        # CI deals with plenty of CUDA versions, which makes discovery tricky.
        # dump a relevant tree of files to help debugging
        function traverse(dir, level=0)
            for entry in readdir(dir)
                print("  "^level)
                path = joinpath(dir, entry)
                if isdir(path)
                    println("└ $entry")
                    traverse(path, level+1)
                else
                    println("├ $entry")
                end
            end
        end
        for dir in dirs
            println("File tree of toolkit directory $dir:")
            traverse(dir)
        end
    end
end

@testset "availability" begin
    @test isa(has_cuda(), Bool)
    @test isa(has_cuda_gpu(), Bool)
end

@testset "call" begin
    # ccall throws if the lib doesn't exist, even if not called
    foo(x) = (x && ccall((:whatever, "nonexisting"), Cvoid, ()); 42)
    @test_throws ErrorException foo(false)

    # @runtime_ccall prevents that
    bar(x) = (x && @runtime_ccall((:whatever, "nonexisting"), Cvoid, ()); 42)
    @test bar(false) == 42
    # but should still error nicely if actually calling the library
    @test_throws ErrorException bar(true)
end