text:
module TopicModelHealthWebsites
using Distributions, StatsBase, Dates, SpecialFunctions, JLD
export HALT_LDA
# Latent Dirichlet Allocation with hierarchical asymmetric prior + local topics (1 per website)
function HALT_LDA(; K::Integer, datapath::String, savepath::String,
                n_iter = 100, n_burn = 0, n_lag = 5,
                c0alpha0_const = 10.0, beta = 0.05, gamma = 0.05, sigma0 = 5.0, sigmac = 0.05, a = 1.0, b = 1.0, c = 1.0, da = 1.0, db = 1.0)
    Ka = K + 1
    c0alpha0 = fill(c0alpha0_const, (Ka))
    data = load(datapath)
    words, V = data["words"], length(data["dictionary"])
    N = [[length(x) for x in y] for y in words]
    Mi = [length(x) for x in N]
    Mis = sum(Mi)
    M = length(Mi)
    topics, n, m0, m1, alpha = initialize(K, Ka, V, N, Mi, Mis, M, beta, gamma, c0alpha0, words)
    n_list, m0_list, m1_list, alpha_list, c0alpha0_list, c_list = sample_multiple(M, Mi, Mis, c, da, db, a, b, N, words, K, Ka, alpha, c0alpha0, gamma, beta, topics, V, n, m0, m1, sigma0, sigmac, n_iter, n_burn, n_lag)
    save(savepath, # MCMC samples
                    "n_list", n_list,
                    "m0_list", m0_list,
                    "m1_list", m1_list,
                    "alpha_list", alpha_list,
                    "c0alpha0_list", c0alpha0_list,
                    "c_list", c_list,
                    # Counts from data and number of topics
                    "K", K,
                    "Ka", Ka,
                    "V", V,
                    "M", M,
                    "Mi", Mi,
                    "N", N,
                    # Current state
                    "topics", topics,
                    "n", n,
                    "m0", m0,
                    "m1", m1,
                    "alpha", alpha,
                    "c0alpha0", c0alpha0,
                    # Fixed hyperparameters
                    "beta", beta,
                    "gamma", gamma,
                    "a", a,
                    "b", b,
                    "c", c,
                    "da", da,
                    "db", db,
                    # M-H
                    "sigma0", sigma0,
                    "sigmac", sigmac)
end
function initialize(K, Ka, V, N, Mi, Mis, M, beta, gamma, c0alpha0, words) # Initialize topic and counts
    # theta
    theta_init = fill(1/Ka, Ka)
    # Generate topics
    topics = [[rand(Categorical(theta_init), N[i][j]) for j in 1:Mi[i]] for i in 1:M]
    topics_nump = topics
    w_nump = words
    alpha = [rand(Dirichlet(c0alpha0)) for i in 1:M]
    # Page-topic counts
    m1 = [[counts(topics_nump[i][j], Ka) for j in 1:Mi[i]] for i in 1:M]
    # Website-topic-word counts
    topics_nump = [collect(Iterators.flatten(x)) for x in topics_nump]
    w_nump = [collect(Iterators.flatten(x)) for x in w_nump]
    m0 = [counts(w_nump[i][findall(x -> x == Ka, topics_nump[i])], V) for i in 1:M]
    # Global-topic-word counts
    topics_nump = collect(Iterators.flatten(topics_nump))
    w_nump = collect(Iterators.flatten(w_nump))
    n = [counts(w_nump[findall(x -> x == k, topics_nump)], V) for k in 1:K]
    # Add hyperparameters
    n = hcat(n...)' .+ beta
    m0 = hcat(m0...)' .+ gamma
    m1 = [hcat(m1[x]...)' .+ 0.0 for x in 1:M] # needs to update in sampling at each iteration
    return(topics, n, m0, m1, alpha)
end
function sample_antoniak(;N, calpha)
    if N > 0
        lambda = sum([rand(Bernoulli(calpha/(calpha+ind))) for ind in 0:N-1])
    else
        lambda = 0
    end
    return(lambda)
end
# likelihood ratio for c0alpha0
function lpiratio_c0alpha0(c0alpha0_new, c0alpha0, alpha, da, db, M, k)
    lpiratio = M*(loggamma(sum(c0alpha0_new)) - loggamma(c0alpha0_new[k]) - loggamma(sum(c0alpha0)) + loggamma(c0alpha0[k])) +
        sum(log.([x[k] for x in alpha])) * (c0alpha0_new[k] - c0alpha0[k]) +
        (da-1)*(log(c0alpha0_new[k]) - log(c0alpha0[k])) -
        db*(c0alpha0_new[k] - c0alpha0[k])
    return(lpiratio)
end
# likelihood ratio for c
function lpiratio_c(c_new, c, a, b, N, M, Mi, Mis, m1, alpha)
    lpiratio = (a-1)*(log(c_new)-log(c)) -
        b*(c_new-c) +
        Mis*(loggamma(c_new)-loggamma(c)) -
        sum([sum([loggamma(c_new + N[i][j]) - loggamma(c + N[i][j]) for j in 1:Mi[i]]) for i in 1:M]) +
        sum([sum([sum((loggamma.(c_new*alpha[i] .+ m1[i][j,:]) - loggamma.(c_new*alpha[i])) -
            (loggamma.(c*alpha[i] .+ m1[i][j,:]) - loggamma.(c*alpha[i]))) for j in 1:Mi[i]]) for i in 1:M])
    return(lpiratio)
end
function sample_multiple(M, Mi, Mis, c, da, db, a, b, N, words, K, Ka, alpha, c0alpha0, gamma, beta, topics, V, n, m0, m1, sigma0, sigmac,
                        n_iter = 100, n_burn = 0, n_lag = 5)
    alpha_list = []
    c0alpha0_list = []
    c_list = []
    n_list = []
    m0_list = []
    m1_list = []
    for iter in 1:n_iter
        keep = (iter % n_lag == 0) * (iter >= n_burn)
        # Sample topics for each word
        for i in 1:M
            for j in 1:Mi[i]
                for h in 1:N[i][j]
                    # Data + unknown categories
                    w = words[i][j][h]
                    z = topics[i][j][h]
                    # Remove word s,t,u from counts
                    m1[i][j,z] -= 1
                    if z < Ka
                        n[z,w] -= 1
                    else
                        m0[i,w] -= 1
                    end
                    # Sample levels and topics
                    u_1 = rand(Uniform(0,1))
                    # Calculate the level 1 sampling probability
                    part_1 = n[:,w] ./ sum(n, dims = 2) .* (m1[i][j,1:K] .+ c*alpha[i][1:K])
                    # Calculate the level 0 sampling probability
                    part_0 = m0[i,w] / sum(m0[i,:]) * (m1[i][j,Ka] + c*alpha[i][Ka])
                    sum_p = sum(part_1) + part_0
                    # Draw topic
                    p_cum = part_0
                    if u_1 > p_cum/sum_p
                        for k = 1:K
                            p_cum += part_1[k]
                            if u_1 <= p_cum/sum_p
                                topics[i][j][h] = k
                                n[k,w] += 1
                                m1[i][j,k] += 1
                                break
                            end
                        end
                    else
                        topics[i][j][h] = Ka
                        m0[i,w] += 1
                        m1[i][j,Ka] += 1
                    end
                end
            end
        end
        # Sample site means
        for i in 1:M
            lambdas = sum([[sample_antoniak(N = m1[i][j,k], calpha = c*alpha[i][k]) for k in 1:Ka] for j in 1:Mi[i]])
            alpha[i,] = rand(Dirichlet(c0alpha0 .+ lambdas))
        end
        # Sample global means (univariate proposals)
        for k in 1:Ka
            c0alpha0_k = c0alpha0[k]
            c0alpha0_k_new = rand(truncated(Normal(c0alpha0_k, sigma0), 0, c0alpha0_k*2))
            c0alpha0_new = copy(c0alpha0)
            c0alpha0_new[k] = c0alpha0_k_new
            diff = lpiratio_c0alpha0(c0alpha0_new, c0alpha0, alpha, da, db, M, k) + logpdf(truncated(Normal(c0alpha0_k_new, sigma0), 0, c0alpha0_k_new*2), c0alpha0_k) - logpdf(truncated(Normal(c0alpha0_k, sigma0), 0, c0alpha0_k*2), c0alpha0_k_new)
            ar = exp(diff)
            if ar >= rand(Uniform(0,1))
                c0alpha0[k] = c0alpha0_k_new
            end
        end
        # Sample c
        c_new = rand(truncated(Normal(c, sigmac), 0, c*2))
        ar = exp(lpiratio_c(c_new, c, a, b, N, M, Mi, Mis, m1, alpha) + logpdf(truncated(Normal(c_new, sigmac), 0, c_new*2), c) - logpdf(truncated(Normal(c, sigmac), 0, c*2), c_new))
        if ar >= rand(Uniform(0,1))
            c = c_new
        end
        # Save
        if keep == true
            push!(alpha_list, copy(alpha))
            push!(c_list, copy(c))
            push!(c0alpha0_list, copy(c0alpha0))
            push!(n_list, copy(n))
            push!(m0_list, copy(m0))
            push!(m1_list, copy([copy(x) for x in m1]))
            print(Dates.format(now(), "HH:MM:SS") * ": Iteration: " * string(iter) * "\n")
        end
    end
    return([n_list, m0_list, m1_list, alpha_list, c0alpha0_list, c_list])
end
end # module

max_stars_repo_name: jwanghb/publichealth-websites

max_issues_repo_name: jwanghb/publichealth-websites

max_forks_repo_name: jwanghb/publichealth-websites


text:
include("constants.jl")

"""
Parse input string s as a float, otherwise return original string.
"""
function parsestring(s)
    x = tryparse(Float64, s)
    if isnothing(x)
        x = s
    end
    return x
end

"""
Return non-comment (#) lines from input file parsed as 2D array
"""
function readfile(filename)
    lines = readlines(filename)
    dat = [parsestring.(split(strip(line))) for line in lines if line[1] ≠ '#']
    return hcat(dat...)'
end

function loadmethyls(filename)
    println("load methyls: $filename")
    return readlines(filename)
end

function loadR2(filename, B0_MHz)
    println("load R2: $filename")
    println("   B0 (1H, in MHz) = $B0_MHz")

    tmp = readfile(filename)
    B0 = B0_MHz' * 1e6 * 2π / γH;
    R2_zq = tmp[:,1:8:end] .± tmp[:,2:8:end]
    R2_dq = tmp[:,3:8:end] .± tmp[:,4:8:end]
    R2_dqq = tmp[:,5:8:end] .± tmp[:,6:8:end]
    R2_qq = tmp[:,7:8:end] .± tmp[:,8:8:end]

    #return R2set(R2_zq, R2_dq, R2_dqq, R2_qq, B0)

    data = Vector{R2Data}()
    nmethyls, ~ = size(R2_zq)
    for i=1:nmethyls
        push!(data, R2Data(R2_zq[i,:], R2_dq[i,:], R2_dqq[i,:], R2_qq[i,:], B0))
    end
    return data
end

function loadS2tc(filename)
    println("load S2tc: $filename")

    S2tc = readfile(filename)
    return S2tc[:,1] .± S2tc[:,2]
end

function load1Hcsa(ηHfilename, S2tc, B0_MHz)
    println("load and calculate 1H CSA:")
    println("   ηH filename: $ηHfilename")
    println("   B0 (1H, in MHz) = $B0_MHz")

    ηH = readfile(ηHfilename)
    ηH = ηH[:,1] .± ηH[:,2]

    B0 = B0_MHz * 1e6 * 2π / γH;

    ΔσH = @. 1e6*ηH/(-(4/15)*(μ₀/(4π))*ħ*γH^2*γC*P2cosβ*rCH^-3*B0*S2tc*1e-9)
    return ΔσH
end

function load13Ccsa(filename)
    println("load 13C csa: $filename")

    ΔσC = readfile(filename)
    return ΔσC[:,1] .± ΔσC[:,2]
end

"""
calculate and subtract CSA contributions from Hahn-echo relaxation rates (for 1H coherence order pH)
"""
function subtractCSA(R2data::Vector{R2Data}, S2tc, ΔσH, ΔσC)
    result = similar(R2data)
    for i=1:length(R2data)
        R2s = R2data[i]
        ξHcsa = 2/(3*sqrt(5)) * γH * ΔσH[i] * 1e-6 * sqrt(S2tc[i] * 1e-9)
        ξCcsa = 2/(3*sqrt(5)) * γC * ΔσC[i] * 1e-6 * sqrt(S2tc[i] * 1e-9)

        B0 = R2s.B0'
        zq = @. R2s.ZQ - (ξCcsa - ξHcsa)^2 * B0^2
        dq = @. R2s.DQ - (ξCcsa + ξHcsa)^2 * B0^2
        dqq = @. R2s.QQ⁻ - (ξCcsa - 3ξHcsa)^2 * B0^2
        qq = @. R2s.QQ⁺ - (ξCcsa + 3ξHcsa)^2 * B0^2

        result[i] = R2Data(zq, dq, dqq, qq, B0)
    end
    return result

    # ξHcsa = @. 2/(3*sqrt(5)) * γH * ΔσH * 1e-6 * sqrt(S2tc * 1e-9)
    # ξCcsa = @. 2/(3*sqrt(5)) * γC * ΔσC * 1e-6 * sqrt(S2tc * 1e-9)

    # B0 = R2s.B0
    # zq = @. R2s.ZQ - (ξCcsa - ξHcsa)^2 * B0^2
    # dq = @. R2s.DQ - (ξCcsa + ξHcsa)^2 * B0^2
    # dqq = @. R2s.QQ⁻ - (ξCcsa - 3ξHcsa)^2 * B0^2
    # qq = @. R2s.QQ⁺ - (ξCcsa + 3ξHcsa)^2 * B0^2

    # return R2set(zq, dq, dqq, qq, B0)
end

"""
load CPMG data

# Parameters:
type = :MQ or :SQ1H
Trelax = relaxation time (in seconds)
B0_MHz = 1H Larmor frequency (in MHz)
"""
function loadCPMG(filename, methyls, type, B0_MHz, Trelax; min_err=nothing)
    nmethyls = length(methyls)
    println("load CPMG data: $filename ($nmethyls methyls)")
    println("   experiment type = $(String(type))")
    println("   B₀ = $B0_MHz MHz")
    println("   Trelax = $Trelax s")
    isnothing(min_err) || println("   minimum error = $min_err s⁻¹")

    B0 = B0_MHz * 1e6 * 2π / γH;
    if type == :MQ
        datatype = CPMGData_MQ
    elseif type == :SQ1H
        datatype = CPMGData_SQ1H
    end
    cpmgdata = Vector{datatype}(undef, nmethyls)

    lines = readlines(filename)
    linespermethyl = length(lines) ÷ nmethyls # integer division
    νCPMG = zeros(linespermethyl-4)
    #R2 = zeros(Measurement{Float64}, linespermethyl-4, nmethyls)
    for i=1:nmethyls
        headerline = lines[1 + linespermethyl*(i-1)]
        datalines = lines[3 + linespermethyl*(i-1):linespermethyl*i-2]
        methylname = replace(headerline[2:8], 'H'=>'C') # extract assuming format Rxxx(C/H)xx, e.g. L661HD1, I743CG2, etc.
        methylindex = findfirst(methyls .== methylname)
        data = hcat([parse.(Float64, split(line)) for line in datalines]...)'
        νCPMG = data[:,1]
        err = data[:,4]
        if ~isnothing(min_err)
            @.err[err<min_err] = min_err
        end
        #R2[:,methylindex] = data[:,2] .± err
        cpmgdata[methylindex] = datatype(νCPMG, data[:,2] .± err, B0, Trelax)
    end
    #return νCPMG, R2
    return cpmgdata
end

ext: jl

max_stars_repo_name: chriswaudby/qq-hahnecho

max_issues_repo_name: chriswaudby/qq-hahnecho

max_forks_repo_name: chriswaudby/qq-hahnecho


text:
using SwitchTimeOpt
using MathOptInterface
const MOI = MathOptInterface

using Plots
using Ipopt

# Define Solver options
maxiter = 300
maxtime = 1000.0
verbose = 5
tolerance = 1e-06

#Define solver
solver = Ipopt.Optimizer()
MOI.set(solver, MOI.RawOptimizerAttribute("tol"), tolerance)
MOI.set(solver, MOI.RawOptimizerAttribute("print_level"), verbose)
MOI.set(solver, MOI.RawOptimizerAttribute("max_cpu_time"), maxtime)
MOI.set(solver, MOI.RawOptimizerAttribute("max_iter"), maxiter)

### Define system parameters
#  Time vector
t0 = 0.0
tf = 10.0

# Integer inputs
uvec = repeat([1.0; 2.0], 8, 1)
uvec = copy(uvec')

# Number of switches
N = size(uvec,2) -1

# Cost function matrix
K = 1.0
C = sqrt(K)*[0.0 1.0 -1.0]
Q = C'*C

# Initial state
x0 = [2.0; 2.0; 3.0]

### Define system dynamics
function nldyn(x::Array{Float64,1}, u::Array{Float64,1})
  n = length(x)
  f = zeros(n)

  f[1] = -sqrt(x[1]) + u[1]
  f[2] = sqrt(x[1]) - sqrt(x[2])
  f[3] = -0.05
  return f

end

function nldyn_deriv(x::Array{Float64,1}, u::Array{Float64,1})
  df = [-1/(2*sqrt(x[1]))        0                      0;
        1/(2*sqrt(x[1]))        -1/(2*sqrt(x[2]))       0;
        0                        0                      0]
end

# Generate and solve problem
ngrid = 100

# Initialize the problem
m = stoproblem(
  x0,                 # Initial state
  nldyn,              # Nonlinear dynamics
  nldyn_deriv,        # Nonlinear dynamics derivative
  uvec,               # Vector of integer inputs
  ngrid = ngrid,      # Number of grid points
  t0 = t0,            # InitialtTime
  tf = tf,            # Final time
  Q = Q,              # Cost matrix
  solver = solver)

solve!(m)
objlin = getobjval(m)
tauopt = gettau(m)
status = getstat(m)
xsim, xpts, objode45, t = simulate(m, tauopt)

println("Status: ", status)
println("Objective value: ", objlin)

showPlots = true
savePlots = false
plotSolution(m, "Tank", "Tank", showPlots, savePlots)


text:
mutable struct Market{D<:Distribution, L<:Union{Bool, AbstractFloat}, F, Df, Db <: Distribution}
    stockout_cost::F
    demand_dist::Type{D}
    backorder::Float64
    lostsales::L
    horizon::Int
    forecasts::Vector{Float64}
    last_demand::Float64
    backorder_reset::Db
    forecast_reset::Df
    backorder_log::Vector{Float64}
    fillrate_log::Vector{Float64}
    cost_log::Vector{Float64}
    name::String
end

function Market(stockout_cost, demand_distribution::Type{<:Distribution}, horizon::Int, backorder_reset::NumDist, forecast_reset::State...; lostsales::Bool = false, name="market")
    @assert hasmethod(stockout_cost, Tuple{Market}) "stockout cost must have a method with `(::Market)` signature"
    @assert all(x -> x isa State, forecast_reset)
    @assert length(forecast_reset) == length(params(demand_distribution())) "Please provide a reset distribution for each parameter of $demand_distribution"
    bd = parametrify(backorder_reset)
    frd = Iterators.Stateful.(cycle.(parametrify.(forecast_reset)))
    forecasts = Float64[rand(param) for _ in 1:horizon for param in popfirst!.(frd)]

    Market{demand_distribution, typeof(lostsales), typeof(stockout_cost), typeof(frd), typeof(bd)}(
        stockout_cost, demand_distribution, Float64(rand(bd)), lostsales, horizon, forecasts, 0.0, bd, tuple(frd...), zeros(0), zeros(0), zeros(0), name)
end

function Market(stockout_cost::Number, demand_distribution::Type{<:Distribution}, horizon::Int, backorder_reset::NumDist, forecast_reset::State...; lostsales = false, name="market")
    Market(LinearStockoutCost(stockout_cost), demand_distribution, horizon, backorder_reset, forecast_reset..., lostsales = lostsales, name = name)
end

RLBase.state(ma::Market) = [ma.backorder; ma.forecasts]
state_size(ma::Market) = 1+ma.horizon*length(ma.forecast_reset)
function print_state(ma::Market)
    n_param = length(ma.forecasts) ÷ ma.horizon
    forecasts = ["$(ma.name) demand($j) t+$(i-1)" => p for (i,pars) in enumerate(partition(ma.forecasts, n_param)) for (j,p) in enumerate(pars)]
    return ma.lostsales ? forecasts : ["$(ma.name) backorder" => ma.backorder ; forecasts]
end

function demand!(ma::Market)
    param = ma.forecasts[1:length(ma.forecast_reset)]
    demand = rand(ma.demand_dist(param...))
    ma.last_demand = max(zero(demand), demand)
    return ma.backorder + ma.last_demand
end

function Base.push!(ma::Market, quantity, source)
    ma.backorder += ma.last_demand
    ma.backorder -= min(quantity, ma.backorder)
    return nothing
end

function reward!(ma::Market)
    deleteat!(ma.forecasts, 1:length(ma.forecast_reset))
    push!(ma.forecasts, rand.(popfirst!.(ma.forecast_reset))...)
    cost = ma.stockout_cost(ma)
    push!(ma.backorder_log, ma.backorder)
    push!(ma.fillrate_log, max(0, (1-ma.backorder/ma.last_demand)))
    push!(ma.cost_log, cost)
    ma.backorder *= (1 - ma.lostsales)
    return -cost
end

function RLBase.reset!(ma::Market)
    for fr in ma.forecast_reset
        Iterators.reset!(fr, fr.itr)
        reset!.(fr.itr.xs)
    end
    ma.forecasts = Float64[rand(param) for _ in 1:ma.horizon for param in popfirst!.(ma.forecast_reset)]
    ma.backorder = rand(ma.backorder_reset)
    empty!(ma.backorder_log)
    empty!(ma.cost_log)
    empty!(ma.fillrate_log)
    return nothing
end

inventory_position(ma::Market) = -ma.backorder

mutable struct LinearStockoutCost{T}
    b::T
end
(f::LinearStockoutCost)(ma::Market) = f.b*ma.backorder

Base.show(io::IO, market::Market{D,F,Df,Db}) where {D,F,Df,Db} = print(io, "Market($D, LS:$(market.lostsales))")

max_stars_repo_path: src/elementary blocks/market.jl

max_issues_repo_path: src/elementary blocks/market.jl

max_forks_repo_path: src/elementary blocks/market.jl


text:
using SimpleSolvers
using Test
using GeometricIntegrators
using GeometricIntegrators.Utils
using GeometricProblems.LorenzAttractor

SimpleSolvers.set_config(:nls_atol, 8eps())
SimpleSolvers.set_config(:nls_rtol, 2eps())

const Δt = 0.01
const nt = 1000

const ref = [-4.902687541134471, -3.743872921802973, 24.690858102790042]

@testset "$(rpad("Lorenz Attractor",80))" begin
    ode  = lorenz_attractor_ode()

    sol = integrate(ode, TableauGauss(1), Δt, nt)
    @test relative_maximum_error(sol.q, ref) < 4E-2

    sol = integrate(ode, TableauGauss(2), Δt, nt)
    @test relative_maximum_error(sol.q, ref) < 2E-5

end


text:
export Domain, SegmentDomain, tocanonical, fromcanonical, fromcanonicalD, ∂
export chebyshevpoints, fourierpoints, isambiguous, arclength
export components, component, ncomponents

# add indexing for all spaces, not just DirectSumSpace
# mimicking scalar vs vector

# prectype gives the precision, including for Vec
prectype(::Type{D}) where {D<:Domain} = float(eltype(eltype(D)))
prectype(d::Domain) = prectype(typeof(d))

#TODO: bivariate AnyDomain
struct AnyDomain <: Domain{UnsetNumber} end
struct EmptyDomain <: Domain{Nothing} end

isambiguous(::AnyDomain) = true
dimension(::AnyDomain) = 1

complexlength(::AnyDomain) = NaN
arclength(::AnyDomain) = NaN
arclength(::EmptyDomain) = false
arclength(::DomainSets.EmptySpace) = false

isempty(::AnyDomain) = false

reverseorientation(a::Union{AnyDomain,EmptyDomain}) = a

canonicaldomain(a::Union{AnyDomain,EmptyDomain}) = a

indomain(x::Domain,::EmptyDomain) = false

convert(::Type{Domain{T}}, ::AnyDomain) where T = Domain(T)

union(::AnyDomain, d::Domain) = d
union(d::Domain, ::AnyDomain) = d

union(::EmptyDomain, ::EmptyDomain) = EmptyDomain()
union(::EmptyDomain, a::Domain) = a
union(a::Domain, ::EmptyDomain) = a

##General routines
isempty(::EmptyDomain) = true

## Interval DomainSets

abstract type SegmentDomain{T} <: Domain{T} end
abstract type AbstractSegment{T} <: SegmentDomain{T} end
const IntervalOrSegment{T} = Union{AbstractInterval{T}, AbstractSegment{T}}
const IntervalOrSegmentDomain{T} = Union{AbstractInterval{T}, SegmentDomain{T}}

canonicaldomain(d::IntervalOrSegmentDomain) = ChebyshevInterval{real(prectype(d))}()

domainscompatible(a,b) = domainscompatible(domain(a),domain(b))
domainscompatible(a::Domain,b::Domain) = isambiguous(a) || isambiguous(b) ||
                    isapprox(a,b)

##TODO: Should fromcanonical be fromcanonical!?

#TODO consider moving these
leftendpoint(d::IntervalOrSegmentDomain{T}) where {T} = fromcanonical(d,-one(eltype(T)))
rightendpoint(d::IntervalOrSegmentDomain{T}) where {T} = fromcanonical(d,one(eltype(T)))

indomain(x,::AnyDomain) = true
function indomain(x,d::SegmentDomain)
    T=float(real(prectype(d)))
    y=tocanonical(d,x)
    ry=real(y)
    iy=imag(y)
    sc=norm(fromcanonicalD(d,ry<-1 ? -one(ry) : (ry>1 ? one(ry) : ry)))  # scale based on stretch of map on projection to interal
    dy=fromcanonical(d,y)
    # TODO: use isapprox once keywords are fast
    ((isinf(norm(dy)) && isinf(norm(x))) ||  norm(dy-x) ≤ 1000eps(T)*max(norm(x),1)) &&
        -one(T)-100eps(T)/sc ≤ ry ≤ one(T)+100eps(T)/sc &&
        -100eps(T)/sc ≤ iy ≤ 100eps(T)/sc
end

ncomponents(s::Domain) = 1
components(s::Domain) = [s]
function components(s::Domain,k)
    k ≠ 1 && throw(BoundsError())
    s
end

issubcomponent(a::Domain,b::Domain) = a in components(b)

##### canoncial
"""
    canonicaldomain(d)

returns a domain which we map to for operations. For example,
the canonical domain for an interval [a,b] is [-1,1]
"""
function canonicaldomain end

"""
    tocanonical(d, x)

maps the point `x` in `d` to a point in `canonical(d,x)`
"""
function tocanonical end

issubset(a::Domain,b::Domain) = a==b

## conveninece routines

ones(d::Domain) = ones(prectype(d),Space(d))
zeros(d::Domain) = zeros(prectype(d),Space(d))

function commondomain(P::AbstractVector)
    ret = AnyDomain()

    for op in P
        d = domain(op)
        @assert ret == AnyDomain() || d == AnyDomain() || ret == d

        if d != AnyDomain()
            ret = d
        end
    end

    ret
end

commondomain(P::AbstractVector,g::AbstractArray{T}) where {T<:Number} = commondomain(P)
commondomain(P::AbstractVector,g) = commondomain([P;g])

domain(::Number) = AnyDomain()

## rand

rand(d::IntervalOrSegmentDomain,k...) = fromcanonical.(Ref(d),2rand(k...)-1)

checkpoints(d::IntervalOrSegmentDomain) = fromcanonical.(Ref(d),[-0.823972,0.01,0.3273484])

## boundary

boundary(d::SegmentDomain) = [leftendpoint(d),rightendpoint(d)] #TODO: Points domain

## map domains
# we auto vectorize arguments
tocanonical(d::Domain,x,y,z...) = tocanonical(d,Vec(x,y,z...))
fromcanonical(d::Domain,x,y,z...) = fromcanonical(d,Vec(x,y,z...))

mappoint(d1::Domain,d2::Domain,x...) = fromcanonical(d2,tocanonical(d1,x...))
invfromcanonicalD(d::Domain,x...) = 1/fromcanonicalD(d,x...)

## domains in higher dimensions

## sorting
# we sort spaces lexigraphically by default

for OP in (:<,:(<=),:>,:(>=),:(isless))
    @eval $OP(a::Domain,b::Domain)=$OP(string(a),string(b))
end

## Other special domains

struct PositiveIntegers <: Domain{Int} end
struct Integers <: Domain{Int} end

const ℕ = PositiveIntegers()
const ℤ = Integers()


text:
struct AutoRegressiveModel{N, T} <: AbstractSampler
    nbits::Int
    masks::Vector{BitMatrix}
    W::NTuple{N, Matrix{T}}
    b::NTuple{N, Vector{T}}
end

function AutoRegressiveModel(nbits::Int, nhiddens)
    masks = createmasks(collect(1:nbits), nhiddens)

    nhiddens = [nbits, nhiddens..., nbits]
    W = ntuple(i -> glorot_uniform(nhiddens[i+1], nhiddens[i]), length(nhiddens)-1)
    b = ntuple(i -> zeros(nhiddens[i+1]), length(nhiddens)-1)

    AutoRegressiveModel(nbits, masks, W, b)
end

#follows https://github.com/karpathy/pytorch-made/blob/master/made.py#L68
function createmasks(order::Vector{Int}, hs::Vector{Int})
    D = length(order)
    L = length(hs)
    ## assign each unit in the hidden layer an integer
    hs_index = Dict()
    hs_index[0] = collect(1:D)
    for l in 1:L
        hs_index[l] = rand(minimum(hs_index[l-1]):D-1, hs[l])
    end
    ## construct the mask matrices
    masks = [hs_index[i] .>= hs_index[i-1]' for i in 1:L]
    push!(masks, hs_index[0] .> hs_index[L]')
    masks
end

function network(model::AutoRegressiveModel{N}, x::AbstractArray) where N
    for n in 1:N
        x = model.b[n] .+ (model.W[n] .* model.masks[n]) * x
        if n < N
            x = relu.(x)
        else
            x = sigmoid.(x)
        end
    end
    x
end

function get_logp(model::AutoRegressiveModel, x::AbstractArray)
    xhat = network(model, x)
    ϵ = 1E-7
    sum(x.*log.(xhat.+ϵ) + (1 .-x).*log.(1.0 .- xhat .+ ϵ), dims=1)
end

function gen_samples(model::AutoRegressiveModel, nbatch::Int)
    x = zeros(Int, model.nbits, nbatch)
    for i = 1:model.nbits
        xhat = network(model, x)
        x[i, :] .= xhat[i,:] .> rand(nbatch)
    end
    x
end

model_parameters(model::AutoRegressiveModel) = (model.W..., model.b...)

function model_dispatch!(model::AutoRegressiveModel{N}, θ) where N
    for n in 1:N
        model.W[n] .= θ[n]
        model.b[n] .= θ[N+n]
    end
end


text:
#
# Helper module for working with command line arguments.
# It understands arguments separated by space(s) symbol.
# Argument itself may consist of name and value. e.g.:
# "file:test.txt" or "arg=value" or "param/val". You may
# use different separators between name and the value.
# See ARG_VAL_SEPARATORS constant for details.
#
# Usage:
#     $julia test.jl 1 arg:val
#     #
#     # test.jl file...
#     #
#     args = CommandLine.create()
#     println(CommandLine.has(args, "1"))   # prints true
#     println(CommandLine.val(args, "1"))   # prints ""
#     println(CommandLine.val(args, "unk")  # prints false
#
# @author DeadbraiN
#
module CommandLine
  export create
  export val
  export has

  #
  # This separator set is used for splitting name and value
  # of one argument
  #
  const ARG_VAL_SEPARATORS = Set(":=/")
  #
  # Parses command line arguments and returns map of arg
  # names like keys and values like values. Example:
  # "file:test.txt" -> {"file" => "test.txt"}.
  # @return {Dict{String, String}}
  #
  function create()
    local args::Dict{String, String} = Dict{String, String}()
    local a::String

    for a in ARGS
      arr = split(a, ARG_VAL_SEPARATORS)
      args[arr[1]] = length(arr) > 1 ? arr[2] : ""
    end

    args
  end
  #
  # Returns specified argument's value. If no value or argument
  # doesn't exist, then returns empty string. In case if parse
  # argument is true, then it will try to parse it's value's type
  # @param args Arguments object returned by create()
  # @param arg Name of the argument
  # @param parse true if we have to try parse the value
  # @return {Any} Argument's value or empty string
  #
  function val(args::Dict{String, String}, arg::String, parse::Bool = false)
    if !haskey(args, arg) return "" end
    if parse return _parse(args, arg) end
    args[arg]
  end
  #
  # Checks if specified key (arg) is in command line
  # @param args Arguments object returned by create()
  # @param arg Name of the argument
  # @return {Bool}
  #
  function has(args::Dict{String, String}, arg::String)
    haskey(args, arg)
  end
  #
  # This function returns parsed command line argument value for specified
  # command. Supported types are in comments below...
  # @param args Command line arguments returned by create()
  # @param arg Name of the argument
  # @return {Any}
  #
  function _parse(args::Dict{String, String}, arg::String)
    local value::String = val(args, arg)

    if isempty(value) return true                                               # Bool
    elseif value == "true" return true                                          # Bool
    elseif value == "false" return false                                        # Bool
    elseif !isnull(tryparse(Int, value)) return parse(Int, value)               # Int
    elseif !isnull(tryparse(Float64, value)) return parse(Float64, value)       # Float64
    elseif value[1] === '[' return parse.(split(value[2:end-1], ","))           # Array{Type,1}
    elseif length(split(value, ".")) === 4 IPv4(value)                          # IPv4
    end

    value                                                                       # String
  end
end

max_stars_repo_name: tmptrash/jevo

max_issues_repo_name: tmptrash/jevo

max_forks_repo_name: tmptrash/jevo


text:
# Interactive seed selection in 2D.
#
# This script demonstrates how to select seed points interatively with PyPlot.

using LinearAlgebra
using FieldTracer
using ForwardDiff
using StaticArrays
using PyPlot

"Return x-type magnetic field from 2D potential using auto-differentiation."
function x_type_autodiff(X, Y, B₀, Bg, dx, dy)

   Az(x, y) = B₀/2*(x^2/dx^2 - y^2/dy^2)
   ∇Az(x, y, z) = ForwardDiff.gradient(p -> Az(p[1], p[2]), SVector(x,y,z))
   ẑ = SVector(0.0, 0.0, 1.0)
   B = Array{Float64,3}(undef, 3, length(X), length(Y))
   for j in axes(Y,1), i in axes(X,1)
      B[:,i,j] = ẑ × ∇Az(X[i], Y[j], 0.0)
      B[3,i,j] += Bg
   end
   B
end

function onclick(event, axs, Bx, By, X, Y)
   xs, ys = event.xdata, event.ydata
   u0 = [xs, ys]
   x1, y1 = trace(Bx, By, xs, ys, X, Y; maxstep=3000, ds=0.1, gridType="ndgrid")
   line = plot(x1, y1, "r--")
   add_arrow(line[1])
end

## Main

# Parameters
Bg = 1.0 # Guide field, [nT]
B₀ = 1.0

fig, axs = plt.subplots(figsize=(12,5))

Δ = 0.1           # grid interval
X = -5:Δ:5        # grid x range
Y = -3:Δ:3        # grid x range
dx, dy = 1.0, 0.1 # scale ratio of the analytic field

@time B = x_type_autodiff(X, Y, B₀, Bg, dx, dy)
axs.streamplot(collect(X), collect(Y), B[1,:,:]', B[2,:,:]')

xstart = [1.0, 3.0, 1.0] # seed starting point x
ystart = [0.0, 0.0, 1.0] # seed starting point y

for i in axes(xstart,1)
   x1, y1 = trace(B[1,:,:], B[2,:,:], xstart[i], ystart[i], X, Y;
      maxstep=3000, ds=0.1, gridType="ndgrid")
   axs.plot(x1, y1, "r--")
end
axs.plot(0, 0, color="k", marker="o", label="x-type null point")
axs.legend()
axs.axis("equal")

# Interactive method 1: draw a new streamline on mouse click
fig.canvas.mpl_connect("button_press_event",
   event -> onclick(event, axs, B[1,:,:], B[2,:,:], X, Y))

# Interactive mthod 2: select n seeds and then draw streamlines
#=
x = plt.ginput(3)
@show x
for i in 1:length(x)
   x1, y1 = trace2d(B[1,:,:], B[2,:,:], x[i][1], x[i][2], X, Y;
      maxstep=3000, ds=0.1, gridType="ndgrid")
   axs.plot(x1, y1, "r--")
end
=#


text:
using Base.Test
using HeavyGas

# Read a HGSYSTEM output file
function readhgfile(filename)
    x = Float64[]
    y = Float64[]
    f = open(filename)
    for ln in eachline(f)
        sl = split(ln)
        if !isempty(sl) && isa(parse(sl[1]), Number)
            push!(x,Float64(parse(sl[1])))
            push!(y,Float64(parse(sl[2])))
        end
    end
    close(f)
    return x, y
end

m = GasData(3, 10, 0.03, 15+273.15, 40.53, HeavyGas.D)
x, Meff, Beff, Sy, ypol, b, Heff, ueff = dispersion(0.1, 497.4, 2.5, m, E0 = 1.)

hgpath = joinpath(dirname(@__FILE__), "hgsystem") # Base path to HGSYSTEM reference files
x_hg, vfrac_hg = readhgfile(joinpath(hgpath, "HTAG01.PS1"))
x_hg, Sy_hg = readhgfile(joinpath(hgpath, "HTAG01.PS2"))
x_hg, Beff_hg = readhgfile(joinpath(hgpath, "HTAG01.PS3"))
x_hg, Sz_hg = readhgfile(joinpath(hgpath, "HTAG01.PS4"))
x_hg, b_hg = readhgfile(joinpath(hgpath, "HTAG01.PS6"))
x_hg, Heff_hg = readhgfile(joinpath(hgpath, "HTAG01.PS8"))

x_hg -= 2.5 # Offset of 2.5 m in HG ref data
const Vm_hg = 8314.3*(15+273.15)/101325
const α = 0.22664
const β = 1 + α
Meff_hg = 2*m.uref*gamma((1+α)/β)/(β*m.href^α*Vm_hg)*Beff_hg.*Sz_hg.^(1+α)

@test x_hg ≈ collect(x)
@test maximum(abs(vfrac_hg-(100*ypol))) < 2.
@test maximum(abs(Sy_hg-Sy)) < 0.2
@test maximum(abs(Beff_hg-Beff)) < 0.3
@test maximum(abs(b_hg-b)) < 0.3
@test maximum(abs(Heff_hg-Heff)) < 0.7
@test maximum(abs(Meff_hg-Meff)) < 12.

# using Plots
# pyplot()
# plot(x, Beff, label="Julia")
# plot!(x_hg, Beff_hg, label="HG")
# show()


text:
# Each mask is a tuple, (n, (t0, t1)) indicating that the neuron
# at index n is masked over the time interval (t0, t1)
const Mask = Tuple{Int64,Tuple{Float64,Float64}}

"""
Computes log-likelihood in masked regions.
"""
function log_like(
        model::SeqModel,
        spikes::Vector{Spike},
        masks::Vector{Mask}
    )

    globals = model.globals
    offsets = globals.neuron_response_offsets
    widths = globals.neuron_response_widths
    warps = model.priors.warp_values
    ll = 0.0

    # == FIRST TERM == #
    # -- Sum of Poisson Process intensity at all spikes -- #

    for x in spikes

        # Compute intensity.
        g = globals.bkgd_amplitude * exp(globals.bkgd_log_proportions[x.neuron])
        for event in model.sequence_events
            w = warps[event.sampled_warp]
            r = event.sampled_type
            g += event.sampled_amplitude * normpdf(
                event.sampled_timestamp + w * offsets[x.neuron, r],
                w * sqrt(widths[x.neuron, r]),
                x.timestamp
            )
        end

        # Add term to log-likelihood.
        ll += log(g)
    end

    # == SECOND TERM == #
    # -- Penalty on integrated intensity function -- #

    for (n, (start, stop)) in masks

        # Add contribution of background.
        ll -= (
            globals.bkgd_amplitude
            * exp(globals.bkgd_log_proportions[n])
            * (stop - start)
        )

        # Add contribution of each latent event.
        for event in model.sequence_events
            w = warps[event.sampled_warp]
            r = event.sampled_type
            g = Normal(
                w * offsets[n, r] + event.sampled_timestamp,
                w * sqrt(widths[n, r])
            )
            ll -= (
                event.sampled_amplitude
                * exp(globals.neuron_response_log_proportions[n, r])
                * (cdf(g, stop) - cdf(g, start))
            )
        end

    end

    return ll

end

"""
Computes log-likelihood of homogeneous poisson process
of spikes within a masked region. The `spikes` vector
should only contain spikes within some masked region.
"""
function homogeneous_baseline_log_like(
        spikes::Vector{Spike},
        masks::Vector{Mask}
    )

    # Count number of spikes and total time observed for each neuron.
    time_per_neuron = Dict{Int64,Float64}()
    spikes_per_neuron = Dict{Int64,Int64}()

    for (n, (t0, t1)) in masks

        # Add neuron n to dict keys
        if !(n in keys(time_per_neuron))
            time_per_neuron[n] = 0.0
            spikes_per_neuron[n] = 0
        end

        # Add masked time interval to neuron n
        time_per_neuron[n] += (t1 - t0)
    end

    # Count number of spikes for each neuron
    for x in spikes
        spikes_per_neuron[x.neuron] += 1
    end

    # Compute MLE estimate for a homogeneous PP.
    neurons = keys(time_per_neuron)
    mle_rates = Dict(n => max(eps(), spikes_per_neuron[n]) / time_per_neuron[n] for n in neurons)

    # Use MLE to compute total log-likelihood.
    ll = 0.0
    for n in neurons
        ll += spikes_per_neuron[n] * log(mle_rates[n])
        ll -= time_per_neuron[n] * mle_rates[n]
    end
    return ll

end

log_like(
    model::DistributedSeqModel,
    spikes::Vector{Spike},
    masks::Vector{Mask}
) = log_like(model.primary_model, spikes, masks)


text:
module AMDGPU

### Imports ###

using CEnum
using Setfield
using Libdl
using LLVM, LLVM.Interop
using GPUArrays
using GPUCompiler
using Adapt
using Requires
import LinearAlgebra
using Printf
import Core: LLVMPtr

### Exports ###

export HSAAgent, HSAQueue, HSAExecutable, HSAKernelInstance, HSASignal
export get_agents, profile, get_first_isa, launch!
export get_default_agent, get_default_queue

export ROCArray, ROCVector, ROCMatrix, ROCVecOrMat
export roc, roczeros, rocones, rocfill

### HSA Runtime ###

const libhsaruntime = "libhsa-runtime64.so.1"
include(joinpath(@__DIR__, "hsa", "HSA.jl"))
import .HSA: Agent, Queue, Executable, Status, Signal

struct Adaptor end

include("extras.jl")
include("error.jl")
include("agent.jl")
include("queue.jl")
include("signal.jl")
include("executable.jl")
include("kernel.jl")
include("memory.jl")

@enum DeviceRuntime HSA_rt OCL_rt
const RUNTIME = Ref{DeviceRuntime}(HSA_rt)
#=
if get(ENV, "AMDGPUNATIVE_OPENCL", "") != ""
    RUNTIME[] = OCL_rt
end
=#
include("runtime.jl")

# Device sources must load _before_ the compiler infrastructure
# because of generated functions.
include(joinpath("device", "tools.jl"))
include(joinpath("device", "pointer.jl"))
include(joinpath("device", "array.jl"))
include(joinpath("device", "gcn.jl"))
include(joinpath("device", "runtime.jl"))
include(joinpath("device", "llvm.jl"))
include(joinpath("device", "globals.jl"))

include("compiler.jl")
include("execution_utils.jl")
include("execution.jl")
include("exceptions.jl")
include("reflection.jl")

### ROCArray ###

include("array.jl")
#include("subarray.jl")
#include("utils.jl")
#include("indexing.jl")
include("broadcast.jl")
#include("matmul.jl")
#include("mapreduce.jl")
#include("gpuarray_interface.jl")

allowscalar(x::Bool) = nothing

### Initialization and Shutdown ###

const HSA_REFCOUNT = Threads.Atomic{UInt}(0)
hsaref!() = Threads.atomic_add!(HSA_REFCOUNT, UInt(1))
function hsaunref!()
    if Threads.atomic_cas!(HSA_REFCOUNT, UInt(1), UInt(0)) == UInt(1)
        HSA.shut_down()
    else
        Threads.atomic_sub!(HSA_REFCOUNT, UInt(1))
    end
end

function __init__()
    # Load binary dependencies
    include(joinpath(dirname(@__DIR__), "deps", "loaddeps.jl"))

    if hsa_configured
        # Make sure we load the library found by the last `] build`
        push!(Libdl.DL_LOAD_PATH, dirname(libhsaruntime_path))
        # TODO: Do the same (if possible) for the debug library

        # Initialize the HSA runtime
        status = HSA.init()
        if status == HSA.STATUS_SUCCESS
            hsaref!()
            # Register shutdown hook
            atexit() do
                hsaunref!()
            end

            # Populate the default agent
            agents = get_agents(:gpu)
            if length(agents) > 0
                DEFAULT_AGENT[] = first(agents)
            end

            # Load ROCm external libraries
            if ext_libs_configured
                libhip !== nothing         && include(joinpath(@__DIR__, "hip", "HIP.jl"))
                librocblas !== nothing     && include(joinpath(@__DIR__, "blas", "rocBLAS.jl"))
                librocfft !== nothing      && include(joinpath(@__DIR__, "fft", "rocFFT.jl"))
                #librocsparse !== nothing  && include("sparse/rocSPARSE.jl")
                #librocalution !== nothing && include("solver/rocALUTION.jl")
                librocrand !== nothing && include(joinpath(@__DIR__, "rand", "rocRAND.jl"))
                #libmiopen !== nothing     && include("dnn/MIOpen.jl")

                # Ensure external libraries are up to date
                function check_library(name, path)
                    path === nothing && return
                    if !ispath(path)
                        @warn "$name library has changed. Please run Pkg.build(\"AMDGPU\") and restart Julia."
                    end
                end
                check_library("rocBLAS", librocblas)
                check_library("rocSPARSE", librocsparse)
                check_library("rocALUTION", librocalution)
                check_library("rocFFT", librocfft)
                check_library("rocRAND", librocrand)
                check_library("MIOpen", libmiopen)
            else
                @warn """
                ROCm external libraries have not been built, runtime functionality will be unavailable.
                Please run Pkg.build("AMDGPU") and reload AMDGPU.
                """
            end
        else
            @warn "HSA initialization failed with code $status"
        end
        # we need to load it after rocRAND.jl
        include(joinpath(@__DIR__, "random.jl"))
    else
        @warn """
        HSA runtime has not been built, runtime functionality will be unavailable.
        Please run Pkg.build("AMDGPU") and reload AMDGPU.
        """
    end

    # Check whether device intrinsics are available
    if !ispath(device_libs_path)
        @warn """
        ROCm-Device-Libs have not been downloaded, device intrinsics will be unavailable.
        Please run Pkg.build("AMDGPU") and reload AMDGPU.
        """
    end

    # Load optional OpenCL integrations
    @require OpenCL="08131aa3-fb12-5dee-8b74-c09406e224a2" include("opencl.jl")

    # Load optional @requires packages
    @require ForwardDiff="f6369f11-7733-5829-9624-2563aa707210" include("forwarddiff.jl")
end

end # module


text:
@testset "Plots recipe" begin
    mtg = read_opf("files/simple_plant.opf", Dict)

    recipe = RecipesBase.apply_recipe(Dict{Symbol,Any}(), mtg)

    @test length(recipe) == length(mtg)
    #? NB: both should be equal because 6 lines (all except first point have an edge),
    #? and 1 scatter.

    @test recipe[1].plotattributes == Dict{Symbol,Any}(:label => "", :seriescolor => :black, :seriestype => :line)
    @test recipe[1].args == ([0.0, 0.0], [0.0, 0.2]) # coordinates of the two vertex for the edge
    @test recipe[6].args == ([0.0, 0.7071067811865475], [1.6, 2.3071067811865476]) # same, last one

    df_coordinates = PlantGeom.mtg_coordinates_df(mtg, force = true)

    @test recipe[7].args[1] == df_coordinates.XX
    @test recipe[7].plotattributes == Dict{Symbol,Any}(
        :color => [:black, :black, :black, :black, :black, :black, :black],
        :palette => colorschemes[:viridis],
        :label => "",
        :hover => [
            "name: `node_1`, link: `/`, symbol: `Scene`, index: `1`",
            "name: `node_2`, link: `/`, symbol: `Individual`, index: `2`",
            "name: `node_3`, link: `/`, symbol: `Axis`, index: `3`",
            "name: `node_4`, link: `/`, symbol: `Internode`, index: `4`",
            "name: `node_5`, link: `+`, symbol: `Leaf`, index: `5`",
            "name: `node_6`, link: `<`, symbol: `Internode`, index: `6`",
            "name: `node_7`, link: `+`, symbol: `Leaf`, index: `7`",
        ],
        :seriestype => :scatter,
        :colorbar_entry => false
    )
end


text:
using Revise
using Test
using NNFEM
using PyCall
using PyPlot
using JLD2
using ADCME
using LinearAlgebra
reset_default_graph()

include("nnutil.jl")

testtype = "NeuralNetwork1D"
include("NNTrussPull_Domain.jl")

#use only rho
prop = Dict("name"=> testtype, "rho"=> 8000.0, "E"=> 200e3, "nu"=> 0.45,
           "sigmaY"=>0.3e3, "K"=>1/9*200e3, "B"=> 0.0, "A0"=> 0.005, "nn"=>post_nn)

n_data = [1,2, 4,5]

function approximate_stress(tid::Int64, method::String)
    nodes, EBC, g, gt, FBC, fext, ft, npoints, node_to_point = BoundaryCondition(tid)
    #domain = Domain(nodes, elements, ndofs, EBC, g, FBC, fext)

    @load "Data/domain$(tid).jld2" domain
    setGeometryPoints!(domain, npoints, node_to_point)
    state = zeros(domain.neqs)
    ∂u = zeros(domain.neqs)
    globdat = GlobalData(state,zeros(domain.neqs), zeros(domain.neqs),∂u, domain.neqs, gt, ft)
    assembleMassMatrix!(globdat, domain)

    state_history, full_fext_history = read_data("$(@__DIR__)/Data/$(tid).dat")

    #update state history and fext_history on the homogenized domain
    fext_history = []
    setNeumannBoundary!(domain, FBC, fext)
    for i = 1:NT
        globdat.time = Δt*i
        updateDomainStateBoundary!(domain, globdat)
        push!(fext_history, domain.fext[:])
    end

    E_all, S_all = LSfittingStress(domain, globdat, hcat(state_history...), hcat(fext_history...), Δt, method)

end

nodes, _, _, _, _, _, _, _, _ = BoundaryCondition(tid)
elements = []
for i = 1:nx
    elnodes = [i, i+1]; coords = nodes[elnodes,:];
    push!(elements, FiniteStrainTruss(coords,elnodes, prop, ngp))
end

loss = constant(0.0)

method = "Constant"
for tid in n_data
    E_all, S_all = approximate_stress(tid, method)
    X = zeros(Float64, (NT-1)*ngp, 3)
    Y = zeros(Float64, (NT-1)*ngp)
    k = 1
    for it = 2:NT
        for igp = 1:ngp
            X[k,:] = [E_all[it,igp]; E_all[it-1,igp]; S_all[it-1,igp]]#ε, ε0, σ0
            Y[k] = S_all[it,igp]
            k = k + 1
        end
    end

    y = squeeze(nn(constant(X[:,1:1]), constant(X[:,2:2]), constant(X[:,3:3])))

    global loss
    loss += sum((y - Y)^2) #/stress_scale^2
end

sess = Session(); init(sess)
@show run(sess, loss)
# ADCME.load(sess, "Data/$(nntype)/NNPreLSfit.mat")
if !isdir("Data/$(nntype)")
    mkdir("Data/$(nntype)")
end

for i = 1:5
    BFGS!(sess, loss, 1000)
    ADCME.save(sess, "Data/$(nntype)/NNPreLSfit_nn$(idx)_ite$(i).mat")
end

tid = 3
strain, stress = read_strain_stress("Data/$(tid).dat")
X, Y = prepare_strain_stress_data1D(strain, stress )
y = squeeze(nn(constant(X[:,1:1]), constant(X[:,2:2]), constant(X[:,3:3])))
out = run(sess, y)
plot(X[:,1], out,"+", label="NN")
plot(X[:,1], Y, ".", label="Exact")
#legend()
savefig("nnpreLSfit_$(nntype)_truss1d_stress$tid.png")


text:
module MainPageController

using Genie.Renderer.Html

struct FileLocations
    script1::String
    script2::String
    css1::String
    css2::String
end

searchdir(path,key) = filter(x -> occursin(key, x), readdir(path))

# automatic discovery of file locations
function FileLocations()
    script1 = "/mainpage/js/" * searchdir("./public/mainpage/js", "chunk")[1]
    script2 = "/mainpage/js/" * searchdir("./public/mainpage/js", "chunk")[4]
    css1 = "/mainpage/css/" * searchdir("./public/mainpage/css", "chunk")[1]
    css2 = "/mainpage/css/" * searchdir("./public/mainpage/css", "chunk")[3]
    return FileLocations(script1, script2, css1, css2)
end

function mainpage()
    html(:mainpage, :mainpage, FL=FileLocations())
end

end


text:
# Bias-correction Pope (1990)
function bias_correction(V::VAR)
  (T,K) = size(V.Y')
  A = get_VAR1_rep(V)
  SIGMA = [V.Σ zeros(K,K*V.p-K); zeros(K*V.p-K,K*V.p)]
  vecSIGMAY = (eye((K*V.p)^2)-kron(A,A))\vec(SIGMA)    # Lutkepohl p.29 (2.1.39)
  SIGMAY = reshape(vecSIGMAY, K*V.p, K*V.p)
  I = eye(K*V.p, K*V.p)
  B = A'
  peigen = eigvals(A)
  sumeig = zeros(K*V.p, K*V.p)
  for h in 1:K*V.p
    sumeig += (peigen[h].\(I - peigen[h]*B)) #***
  end
  bias = SIGMA*(inv(I - B) + B/(I-B^2) + sumeig)/(SIGMAY) #***
  Abias = - bias/T
  bcstab = 9           # Arbitrary default value > 1
  delta = 1            # Adjustment factor
  bcA = cell() #### if not working remove and put bcA as global
  while bcstab >= 1
    # Adjust bias-correction proportionately
    bcA = A-delta*Abias
    bcmod = abs(eigvals(bcA))
    any(bcmod.>= 1) ?  bcstab = 1 : bcstab = 0
    delta += - 0.01
    (delta <= 0) && (bcstab = 0)
  end
  return bcA
end


text:
# reso.jl

#include("misc.jl")
#include("primitives.jl")
#include("types.jl")
#include("unify.jl")

# input form
# "[x,y].[+P(x,f(x),-Q(x)]"

# internal form
# CForm = ([:x,:y], [:(+(P(x,f(x)))), :(-(Q(x)))])
# CForm2 = (:C12, [:x,:y], [:(+(P(x,f(x)))), :(-(Q(x)))])

## result
# :NAP is Not APplicable

Clause = Array
Literal = Expr

EmptyClause = []

CForm = Tuple{Array,Array}
EmptyCForm = ([],EmptyClause)

## resolution of 2 clauses
function genvars(vars::Vlist)::Vlist
 map(v->gensym(v),vars)
end

function rename(vars, C::Clause, nvar)::Clause
 if isempty(vars); return C end
 apply(vars, C, nvar)
end

"""
resolution clash the i1'th of c1 with i2'th of c2
"""
function resolution(vars::Vlist, c1::Clause, c2::Clause, i1::Int, i2::Int)
# clause c1,c2 is an array.

 if length(c1) < i1; println("i1 over the c1's length");return :NAP end
 if length(c2) < i2; println("i2 over the c2's length");return :NAP end
 if i1 <= 0; println("i1 should be >=0");return :NAP end
 if i2 <= 0; println("i2 should be >=0");return :NAP end

 lit1 = c1[i1]
 lit2 = c2[i2]

 if lit1.args[1] == lit2.args[1]; return :NAP end
 lit1.args[1] = lit2.args[1]  # temporarily let signs same

 try
  sigma = unify(vars, lit1, lit2)
  c1=vcat(c1[1:(i1-1)],c1[(i1+1):end],c2[1:(i2-1)],c2[(i2+1):end])
  r1,s1 = apply(vars,c1,sigma), sigma
  return r1,s1
 catch e
  if isa(e, Loop)
   return :FAIL
  end
 end
 return :NAP
end

function resolution(var1::Vlist, c1::Clause, i1::Int, var2::Vlist, c2::Clause, i2::Int)
 vs1 = genvars(var1)
 vs2 = genvars(var2)
 var = vcat(vs1,vs2)
 rc1=rename(var1, c1, vs1)
 rc2=rename(var2, c2, vs2)
 rv = resolution(var, rc1, rc2, i1, i2)
 if typeof(rv) == Symbol;return rv end
 r1,sigma = rv
 return (var, r1), sigma
end

function resolution(c1::CForm, i1::Int, c2::CForm,i2::Int)
 cc1=deepcopy(c1)
 cc2=deepcopy(c2)
 rv =resolution(cc1[1],cc1[2],i1, cc2[1],cc2[2],i2)
 if typeof(rv) == Symbol; return rv end
 cf,sigma=rv
 return renamereadable(cf,sigma)
end

"""
reduction
"""
function reduction(vars::Vlist, c1::Clause, i1::Int)
 if length(c1) < i1; println("i1 over the c1's length");return :NAP end
 if i1 <= 0; println("i1 should be >=0");return :NAP end
 lit1 = c1[i1]
 atm1 = lit1.args[2]

 for i2 in 1:length(c1)
  if i2==i1; continue end
   lit2 = c1[i2]
   atm2 = lit2.args[2]
   if atm1.args[1] != atm2.args[1]; continue end
   if lit1.args[1] != lit2.args[1]; continue end

   try
    sigma = unify(vars, lit1, lit2)
    c1=vcat(c1[1:(i1-1)],c1[(i1+1):end])
    r1=apply(vars,c1,sigma)
    return renamereadable((vars, r1), sigma)
   catch e
    continue
   end
 end
 return ((vars,c1),[])
end

function reduction(c1::CForm, i1::Int)
 reduction(c1[1], c1[2], i1)
end

"""
satisfiable check
"""
function satisfiable(vars::Vlist, c1::Clause)
 for i1 in 1:length(c1)
  lit1 = c1[i1]
  atm1 = lit1.args[2]

  for i2 in i1+1:length(c1)
   lit2 = c1[i2]
   atm2 = lit2.args[2]
   if atm1.args[1] != atm2.args[1]; continue end
   if lit1.args[1] == lit2.args[1]; continue end
   lit2.args[1] = lit1.args[1] # temporaly let it same

   try
    sigma = unify(vars, lit1, lit2)
    return true
   catch e
    continue
   end
  end
 end
 return false

end

function satisfiable(c1::CForm)
 satisfiable(c1[1], c1[2])
end

### variable lists

numberingvar(vs::Array{Symbol},i::Int) = map(v->Symbol(v,i),vs)

function readablevars(n::Int)
 rv = [:u,:v,:w,:x,:y,:z]
 vars=[]
 for i in 1:ceil(Int, n/length(rv))
   append!(vars, numberingvar(rv, i))
 end
 vars[1:n]::Vlist
end

function containvar(vars::Vlist, form::Expr)
  cvars=[]::Vlist
  for arg in form.args
   if typeof(arg) == Symbol
    if arg in vars
      push!(cvars, arg)
    end
   else
    pvars =  containvar(vars, arg)
    append!(cvars, pvars)
   end
  end
  return union(cvars,cvars)
end

function containvar(vars::Vlist, cls::Clause)
  cvars=[]::Vlist
  for lit in cls
   pvars =  containvar(vars, lit)
   append!(cvars, pvars)
  end
  return union(cvars,cvars)
end

function fitvars(vars::Vlist, cls::Clause, sub::Tlist)
 fvars = containvar(vars, cls)
 fsub = reducesub(vars, fvars, sub)
 return((fvars, cls), fsub)
end

function renamereadable(clause::CForm, sub::Tlist)
 vars = clause[1]
 cls  = clause[2]
 rvars = readablevars(length(vars))
 rcls = apply(vars, cls, rvars)
 rsub = apply(vars, sub, rvars)
 fcls, fsub = fitvars(rvars, rcls, rsub)
 return fcls,fsub
end

function equalclause(clause1::CForm, clause2::CForm)
 rcls1,rsub1=renamereadable(clause1, [])
 rcls2,rsub2=renamereadable(clause2, [])

 var1=rcls1[1]
 cls1=rcls1[2]
 var2=rcls2[1]
 cls2=rcls2[2]

 if length(var1)!=length(var2);return false end
 if length(cls1)!=length(cls2);return false end

 scls2 = apply(var2, cls2, var1)
 if cls1 == scls2; return true end
 return false
end

function reducesub(vars::Vlist, fvars::Vlist, sub::Tlist)
 fsub  = []
 if isempty(sub); return [] end
 for i in 1:length(vars)
  if vars[i] in fvars
    push!(fsub, sub[i])
  end
 end
 return fsub
end

max_stars_repo_name: sazare/cheapViews

max_issues_repo_name: sazare/cheapViews

max_forks_repo_name: sazare/cheapViews


text:
module BasicStructs
using Parameters
export ResultMetrics,ConfigurtationStruct
"""
constants associated with image over which we will evaluate segmentations
"""
@with_kw struct ImageConstants
    mvspx ::Float64 # Voxelspacing x
    mvspy::Float64 # Voxelspacing y
    mvspz::Float64 #mean Voxelspacing z
    isZConst::Bool # true if slices thickness is the same in all image
    ZPositions::Vector{Float64} # array of true physical positions (in mm) of slices relative to the begining - used in case we have variable thickness of slices
    numberOfVox::Int64 # number of voxels in image
end #ImageConstants
"""
configuration struct that when passed will marks what kind of metrics we are intrested in

"""
@with_kw struct ConfigurtationStruct
    sliceWiseMetrics::Bool = false# if it will be marked as true metrics will be calculated not only globally but also
    dice::Bool = false #dice coefficient
    jaccard::Bool = false #jaccard coefficient
    gce::Bool = false #global consistency error
    vol::Bool = false# Volume metric
    randInd::Bool= false # Rand Index
    ic::Bool= false # interclass correlation
    kc::Bool= false # Kohen Cappa
    mi::Bool= false # mutual information
    vi::Bool= false # variation Of Information
    md::Bool= false # mahalanobis distance
    hd::Bool= false # hausdorff distance
end #ConfigurtationStruct

"""
Struct holding all resulting metrics - if some metric was not calculated its value is just -1
"""
@with_kw mutable struct ResultMetrics
    dice::Float64 = -1.0 #dice coefficient
    jaccard::Float64 =  -1.0 #jaccard coefficient
    gce::Float64 =  -1.0 #global consistency error
    vol::Float64 =  -1.0 # Volume metric
    randInd::Float64 = -1.0 # Rand Index
    ic::Float64 = -1.0 # interclass correlation
    kc::Float64 = -1.0 # Kohen Cappa
    mi::Float64 = -1.0 # mutual information
    vi::Float64 = -1.0 # variation Of Information
    md::Float64 = -1.0 # mahalanobis distance
    hd::Float64 = -1.0 # hausdorff distance

end #ResultMetrics

# """
# struct holding necess
# TN- true negative   TP - true positive  FN - false negative  FP - false positive
# """
# @with_kw struct TnTpFpFn
# # look into https://github.com/JuliaArrays/StructArrays.jl

# end #TnTpFpFn

end#BasicStructs


text:
module SLEEFInline

# export sin, cos, tan, asin, acos, atan, sincos, sinh, cosh, tanh,
#     asinh, acosh, atanh, log, log2, log10, log1p, ilogb, exp, exp2, exp10, expm1, ldexp, cbrt, pow

# fast variants (within 3 ulp)
# export sin_fast, cos_fast, tan_fast, sincos_fast, asin_fast, acos_fast, atan_fast, atan2_fast, log_fast, cbrt_fast

using Base.Math: uinttype, @horner, exponent_bias, exponent_mask, significand_bits, IEEEFloat, exponent_raw_max

const SLEEF = SLEEFInline

export SLEEF

## constants

const MLN2  = 6.931471805599453094172321214581765680755001343602552541206800094933936219696955e-01 # log(2)
const MLN2E = 1.442695040888963407359924681001892137426645954152985934135449406931 # log2(e)

const M_PI  = 3.141592653589793238462643383279502884 # pi
const PI_2 = 1.570796326794896619231321691639751442098584699687552910487472296153908203143099     # pi/2
const PI_4 = 7.853981633974483096156608458198757210492923498437764552437361480769541015715495e-01 # pi/4
const M_1_PI = 0.318309886183790671537767526745028724 # 1/pi
const M_2_PI = 0.636619772367581343075535053490057448 # 2/pi
const M_4_PI = 1.273239544735162686151070106980114896275677165923651589981338752471174381073817     # 4/pi

const MSQRT2 = 1.414213562373095048801688724209698078569671875376948073176679737990732478462102 # sqrt(2)
const M1SQRT2 = 7.071067811865475244008443621048490392848359376884740365883398689953662392310596e-01 # 1/sqrt(2)

const M2P13 = 1.259921049894873164767210607278228350570251464701507980081975112155299676513956 # 2^1/3
const M2P23 = 1.587401051968199474751705639272308260391493327899853009808285761825216505624206 # 2^2/3

const MLOG10_2 = 3.3219280948873623478703194294893901758648313930

const MDLN10E(::Type{Float64}) = Double(0.4342944819032518, 1.098319650216765e-17) # log10(e)
const MDLN10E(::Type{Float32}) = Double(0.4342945f0, -1.010305f-8)

const MDLN2E(::Type{Float64}) = Double(1.4426950408889634, 2.0355273740931033e-17) # log2(e)
const MDLN2E(::Type{Float32}) = Double(1.442695f0, 1.925963f-8)

const MDLN2(::Type{Float64}) = Double(0.693147180559945286226764, 2.319046813846299558417771e-17)  # log(2)
const MDLN2(::Type{Float32}) = Double(0.69314718246459960938f0, -1.904654323148236017f-9)

const MDPI(::Type{Float64})  = Double(3.141592653589793, 1.2246467991473532e-16) # pi
const MDPI(::Type{Float32})  = Double(3.1415927f0, -8.742278f-8)
const MDPI2(::Type{Float64}) = Double(1.5707963267948966, 6.123233995736766e-17) # pi/2
const MDPI2(::Type{Float32}) = Double(1.5707964f0, -4.371139f-8)

const MD2P13(::Type{Float64}) = Double(1.2599210498948732, -2.589933375300507e-17) # 2^1/3
const MD2P13(::Type{Float32}) = Double(1.2599211f0, -2.4018702f-8)

const MD2P23(::Type{Float64}) = Double(1.5874010519681996, -1.0869008194197823e-16) # 2^2/3
const MD2P23(::Type{Float32}) = Double(1.587401f0, 1.9520385f-8)

# Split pi into four parts (each is 26 bits)
const PI_A(::Type{Float64}) = 3.1415926218032836914
const PI_B(::Type{Float64}) = 3.1786509424591713469e-08
const PI_C(::Type{Float64}) = 1.2246467864107188502e-16
const PI_D(::Type{Float64}) = 1.2736634327021899816e-24

const PI_A(::Type{Float32}) = 3.140625f0
const PI_B(::Type{Float32}) = 0.0009670257568359375f0
const PI_C(::Type{Float32}) = 6.2771141529083251953f-7
const PI_D(::Type{Float32}) = 1.2154201256553420762f-10

const PI_XD(::Type{Float32}) = 1.2141754268668591976f-10
const PI_XE(::Type{Float32}) = 1.2446743939339977025f-13

# split 2/pi into upper and lower parts
const M_2_PI_H = 0.63661977236758138243
const M_2_PI_L = -3.9357353350364971764e-17

# Split log(10) into upper and lower parts
const L10U(::Type{Float64}) = 0.30102999566383914498
const L10L(::Type{Float64}) = 1.4205023227266099418e-13

const L10U(::Type{Float32}) = 0.3010253906f0
const L10L(::Type{Float32}) = 4.605038981f-6

# Split log(2) into upper and lower parts
const L2U(::Type{Float64}) = 0.69314718055966295651160180568695068359375
const L2L(::Type{Float64}) = 0.28235290563031577122588448175013436025525412068e-12

const L2U(::Type{Float32}) = 0.693145751953125f0
const L2L(::Type{Float32}) = 1.428606765330187045f-06

const TRIG_MAX(::Type{Float64}) = 1e14
const TRIG_MAX(::Type{Float32}) = 1f7

const SQRT_MAX(::Type{Float64}) = 1.3407807929942596355e154
const SQRT_MAX(::Type{Float32}) = 18446743523953729536f0

include("utils.jl")  # utility functions
include("double.jl") # Dekker style double double functions
include("priv.jl")   # private math functions
include("exp.jl")    # exponential functions
include("log.jl")    # logarithmic functions
include("trig.jl")   # trigonometric and inverse trigonometric functions
include("hyp.jl")    # hyperbolic and inverse hyperbolic functions
include("misc.jl")   # miscallenous math functions including pow and cbrt

# fallback definitions

for func in (:sin, :cos, :tan, :sincos, :asin, :acos, :atan, :sinh, :cosh, :tanh,
             :asinh, :acosh, :atanh, :log, :log2, :log10, :log1p, :exp, :exp2, :exp10, :expm1, :cbrt,
             :sin_fast, :cos_fast, :tan_fast, :sincos_fast, :asin_fast, :acos_fast, :atan_fast, :atan2_fast, :log_fast, :cbrt_fast)
    @eval begin
        $func(a::Float16) = Float16.($func(Float32(a)))
        $func(x::Real) = $func(float(x))
    end
end

for func in (:atan, :hypot)
    @eval begin
        $func(y::Real, x::Real) = $func(promote(float(y), float(x))...)
        $func(a::Float16, b::Float16) = Float16($func(Float32(a), Float32(b)))
    end
end
ldexp(x::Float16, q::Int) = Float16(ldexpk(Float32(x), q))

end


text:
module monty

using Random

"""
The player's strategy is to either switch after monty chooses a door,
or stay
"""
@enum Strategy switch stay

"""
Player either wins a car (a "win" allegedly)
or a goat.
"""
@enum Prize car goat

"""
There are three doors, behind two of which are goats while the third
hides a car.

The game is played in two rounds, first the player picks a door.
Monty (the game show host) then opens one of the other doors to reveal
a goat and offers the player the chance to switch to the third door,
or stay with their first pick.

This function runs a game for a player with a given strategy (stay or switch).
"""
function game(strategy::Strategy)::Prize
    # all possible combinations of goat and car behind doors 1, 2 and 3
    options = [(goat, goat, car),
               (goat, car, goat),
               (car, goat, goat)]

    # draw a random configuration
    doors = rand(options)
    # this is the door with the car
    door_with_car = findfirst(x -> x == car, doors)

    # player's first pick, they select a door at random
    first_pick = rand([1,2,3])
    other_doors = setdiff(Set([1,2,3]), first_pick)

    # from the other doors, monty picks the one with a goat to reveal
    monty_picks = setdiff(other_doors, door_with_car) |> first

    if strategy == switch
        # if the player switches, we choose the remaining door
        door_picked = setdiff(Set([1,2,3]), Set([first_pick, monty_picks])) |> first
    else
        # otherwise we stick with our first pick
        door_picked = first_pick
    end

    # and this is the prize!
    prize = doors[door_picked]
    return prize
end

"""
Runs `n_iterations` of the game, all with a given strategy and returns
the number of wins (player gets a car).
"""
function run(n_iterations::Int64, strategy::Strategy)
    mapreduce(
        x -> game(strategy),
        (acc, p) -> p == car ? acc + 1 : acc,
        1:n_iterations,
        init=0
    )
end

end # module

max_stars_repo_name: mattswoon/monty-hall

max_issues_repo_name: mattswoon/monty-hall

max_forks_repo_name: mattswoon/monty-hall


text:
function ftop(;dir=".",nf::Integer=20)
  @assert nf > 0
  efind =  `find $dir -type f -iname "*" -exec du -sh {} + `
  try
    run(pipeline(efind,`sort -rh`,`head -$nf`))
    println("Done.")
  catch
    println("No files found.")
  end
end


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

@doc doc"""
    AbstractMCMCCallback <: Function

*BAT-internal, not part of stable public API.*

Subtypes (e.g. `X`) must support

    (::X)(level::Integer, chain::MCMCIterator) => nothing
    (::X)(level::Integer, tuner::AbstractMCMCTuner) => nothing

to be compabtible with `mcmc_iterate!`, `mcmc_tune_burnin!`, etc.
"""
abstract type AbstractMCMCCallback <: Function end

@inline Base.convert(::Type{AbstractMCMCCallback}, x::AbstractMCMCCallback) = x

@inline Base.convert(::Type{Vector{<:AbstractMCMCCallback}}, V::Vector{<:AbstractMCMCCallback}) = V

Base.convert(::Type{Vector{<:AbstractMCMCCallback}}, V::Vector) =
            [convert(AbstractMCMCCallback, x) for x in V]

function mcmc_callback_vector(x, idxs::AbstractVector{<:Integer})
    V = convert(Vector{<:AbstractMCMCCallback}, x)
    if eachindex(V) != idxs
        throw(DimensionMismatch("Indices of callback vector incompatible with reference indices"))
    end
    V
end

mcmc_callback_vector(x::Tuple{}, idxs::AbstractVector{<:Integer}) =
    [MCMCNopCallback() for _ in idxs]

@doc doc"""
    MCMCCallbackWrapper{F} <: AbstractMCMCCallback

*BAT-internal, not part of stable public API.*

Wraps a callable object to turn it into an `AbstractMCMCCallback`.

Constructor:

    MCMCCallbackWrapper(f::Any)

`f` needs to support the call syntax of an `AbstractMCMCCallback`.
"""
struct MCMCCallbackWrapper{F} <: AbstractMCMCCallback
    f::F
end

@inline (wrapper::MCMCCallbackWrapper)(args...) = wrapper.f(args...)

Base.convert(::Type{AbstractMCMCCallback}, f::Function) = MCMCCallbackWrapper(f)

struct MCMCNopCallback <: AbstractMCMCCallback end

(cb::MCMCNopCallback)(level::Integer, obj::Any) = nothing

Base.convert(::Type{AbstractMCMCCallback}, ::Tuple{}) = MCMCNopCallback()

struct MCMCMultiCallback{N,TP<:NTuple{N,AbstractMCMCCallback}} <: AbstractMCMCCallback
    funcs::TP

    MCMCMultiCallback(fs::Vararg{AbstractMCMCCallback,N}) where {N} =
        new{N, typeof(fs)}(fs)

    function fMCMCMultiCallback(fs::Vararg{Any,N}) where {N}
        fs_conv = map(x -> convert(AbstractMCMCCallback, x), fs)
        new{N, typeof(fs_conv)}(fs_conv)
    end
end

function (cb::MCMCMultiCallback)(level::Integer, obj::Any)
    for f in cb.funcs
        f(level, obj)
    end
    nothing
end

Base.convert(::Type{AbstractMCMCCallback}, fs::Tuple) = MCMCMultiCallback(fs...)

struct MCMCAppendCallback{T,F<:Function} <: AbstractMCMCCallback
    appendable::T
    max_level::Int
    get_data_func!::F
    nonzero_weights::Bool
end

function (cb::MCMCAppendCallback)(level::Integer, subject::Any)
    if (level <= cb.max_level)
        cb.get_data_func!(cb.appendable, subject, cb.nonzero_weights)
    end
    nothing
end

Base.convert(::Type{AbstractMCMCCallback}, x::DensitySampleVector) = MCMCAppendCallback(x)

MCMCAppendCallback(x::DensitySampleVector, nonzero_weights::Bool = true) =
    MCMCAppendCallback(x, 1, get_samples!, nonzero_weights)


text:
using Rubin
using Tests
using Elliptic
using HypergeometricFunctions
using Polylogarithms
using SpecialFunctions

@test integrate(x^3*arccsch(a+b*x), x) == :(1//4*x^4*arccsch(a+b*x)+-1//4*a^4*b^-4*arccsch(a+b*x)+(1/2)*a*b^-4*(1+-2*a^2)*arctanh((1+(a+b*x)^-2)^(1/2))+-1//3*a*b^-4*(1+(a+b*x)^-2)^(1/2)*(a+b*x)^2+-1//12*b^-4*(1+(a+b*x)^-2)^(1/2)*(2+-17*a^2)*(a+b*x)+1//12*b^-2*x^2*(1+(a+b*x)^-2)^(1/2)*(a+b*x))
@test integrate(x^2*arccsch(a+b*x), x) == :(1//3*x^3*arccsch(a+b*x)+-1//6*b^-3*(1+-6*a^2)*arctanh((1+(a+b*x)^-2)^(1/2))+1//3*a^3*b^-3*arccsch(a+b*x)+-5//6*a*b^-3*(1+(a+b*x)^-2)^(1/2)*(a+b*x)+1//6*x*b^-2*(1+(a+b*x)^-2)^(1/2)*(a+b*x))
@test integrate(x*arccsch(a+b*x), x) == :((1/2)*x^2*arccsch(a+b*x)+(1/2)*b^-2*(1+(a+b*x)^-2)^(1/2)*(a+b*x)+-1*a*b^-2*arctanh((1+(a+b*x)^-2)^(1/2))+-1//2*a^2*b^-2*arccsch(a+b*x))
@test integrate(x^-1*arccsch(a+b*x), x) == :(-1//2*Polylogarithms.polylog(2,exp(2*arccsch(a+b*x)))+arccsch(a+b*x)*log(1+-1*a*(1+(1+a^2)^(1/2))^-1*exp(arccsch(a+b*x)))+arccsch(a+b*x)*log(1+-1*a*(1+-1*(1+a^2)^(1/2))^-1*exp(arccsch(a+b*x)))+-1*arccsch(a+b*x)*log(1+-1*exp(2*arccsch(a+b*x)))+Polylogarithms.polylog(2,a*(1+(1+a^2)^(1/2))^-1*exp(arccsch(a+b*x)))+Polylogarithms.polylog(2,a*(1+-1*(1+a^2)^(1/2))^-1*exp(arccsch(a+b*x))))
@test integrate(x^-2*arccsch(a+b*x), x) == :(-1*x^-1*arccsch(a+b*x)+-1*b*a^-1*arccsch(a+b*x)+2*b*a^-1*(1+a^2)^-1//2*arctanh((1+a^2)^-1//2*(a+tanh((1/2)*arccsch(a+b*x)))))
@test integrate(x^-3*arccsch(a+b*x), x) == :(-1//2*x^-2*arccsch(a+b*x)+(1/2)*a^-2*b^2*arccsch(a+b*x)+-1*a^-2*b^2*(1+a^2)^-3//2*(1+2*a^2)*arctanh((1+a^2)^-1//2*(a+tanh((1/2)*arccsch(a+b*x))))+(1/2)*b*a^-1*x^-1*(1+a^2)^-1*(1+(a+b*x)^-2)^(1/2)*(a+b*x))
@test integrate((a+b*arccsch(c+d*x))^2*(e+f*x)^3, x) == :(1//4*f^-1*(a+b*arccsch(c+d*x))^2*(e+f*x)^4+-2*b^2*d^-4*(d*e+-1*c*f)^3*Polylogarithms.polylog(2,exp(arccsch(c+d*x)))+2*b^2*d^-4*(d*e+-1*c*f)^3*Polylogarithms.polylog(2,-1*exp(arccsch(c+d*x)))+-1//3*b^2*d^-4*f^3*log(c+d*x)+-1//4*d^-4*f^-1*(a+b*arccsch(c+d*x))^2*(d*e+-1*c*f)^4+1//12*b^2*d^-4*f^3*(c+d*x)^2+x*b^2*d^-3*f^2*(d*e+-1*c*f)+b^2*d^-4*f^2*(d*e+-1*c*f)*Polylogarithms.polylog(2,exp(arccsch(c+d*x)))+-1*b^2*d^-4*f^2*(d*e+-1*c*f)*Polylogarithms.polylog(2,-1*exp(arccsch(c+d*x)))+3*f*b^2*d^-4*(d*e+-1*c*f)^2*log(c+d*x)+4*b*d^-4*(d*e+-1*c*f)^3*(a+b*arccsch(c+d*x))*arctanh(exp(arccsch(c+d*x)))+-2*b*d^-4*f^2*(a+b*arccsch(c+d*x))*(d*e+-1*c*f)*arctanh(exp(arccsch(c+d*x)))+-1//3*b*d^-4*f^3*(1+(c+d*x)^-2)^(1/2)*(a+b*arccsch(c+d*x))*(c+d*x)+1//6*b*d^-4*f^3*(1+(c+d*x)^-2)^(1/2)*(c+d*x)^3*(a+b*arccsch(c+d*x))+b*d^-4*f^2*(1+(c+d*x)^-2)^(1/2)*(c+d*x)^2*(a+b*arccsch(c+d*x))*(d*e+-1*c*f)+3*b*f*d^-4*(1+(c+d*x)^-2)^(1/2)*(d*e+-1*c*f)^2*(a+b*arccsch(c+d*x))*(c+d*x))
@test integrate((a+b*arccsch(c+d*x))^2*(e+f*x)^2, x) == :(1//3*f^-1*(a+b*arccsch(c+d*x))^2*(e+f*x)^3+-2*b^2*d^-3*(d*e+-1*c*f)^2*Polylogarithms.polylog(2,exp(arccsch(c+d*x)))+2*b^2*d^-3*(d*e+-1*c*f)^2*Polylogarithms.polylog(2,-1*exp(arccsch(c+d*x)))+-1//3*b^2*d^-3*f^2*Polylogarithms.polylog(2,-1*exp(arccsch(c+d*x)))+-1//3*d^-3*f^-1*(a+b*arccsch(c+d*x))^2*(d*e+-1*c*f)^3+1//3*x*b^2*d^-2*f^2+1//3*b^2*d^-3*f^2*Polylogarithms.polylog(2,exp(arccsch(c+d*x)))+2*f*b^2*d^-3*(d*e+-1*c*f)*log(c+d*x)+4*b*d^-3*(d*e+-1*c*f)^2*(a+b*arccsch(c+d*x))*arctanh(exp(arccsch(c+d*x)))+-2//3*b*d^-3*f^2*(a+b*arccsch(c+d*x))*arctanh(exp(arccsch(c+d*x)))+1//3*b*d^-3*f^2*(1+(c+d*x)^-2)^(1/2)*(c+d*x)^2*(a+b*arccsch(c+d*x))+2*b*f*d^-3*(1+(c+d*x)^-2)^(1/2)*(a+b*arccsch(c+d*x))*(c+d*x)*(d*e+-1*c*f))
@test integrate((a+b*arccsch(c+d*x))^2*(e+f*x), x) == :((1/2)*f^-1*(a+b*arccsch(c+d*x))^2*(e+f*x)^2+f*b^2*d^-2*log(c+d*x)+-2*b^2*d^-2*(d*e+-1*c*f)*Polylogarithms.polylog(2,exp(arccsch(c+d*x)))+2*b^2*d^-2*(d*e+-1*c*f)*Polylogarithms.polylog(2,-1*exp(arccsch(c+d*x)))+-1//2*d^-2*f^-1*(a+b*arccsch(c+d*x))^2*(d*e+-1*c*f)^2+4*b*d^-2*(a+b*arccsch(c+d*x))*(d*e+-1*c*f)*arctanh(exp(arccsch(c+d*x)))+b*f*d^-2*(1+(c+d*x)^-2)^(1/2)*(a+b*arccsch(c+d*x))*(c+d*x))
@test integrate((a+b*arccsch(c+d*x))^2, x) == :(d^-1*(a+b*arccsch(c+d*x))^2*(c+d*x)+-2*b^2*d^-1*Polylogarithms.polylog(2,exp(arccsch(c+d*x)))+2*b^2*d^-1*Polylogarithms.polylog(2,-1*exp(arccsch(c+d*x)))+4*b*d^-1*(a+b*arccsch(c+d*x))*arctanh(exp(arccsch(c+d*x))))
@test integrate((a+b*arccsch(c+d*x))^2*(e+f*x)^-1, x) == :(f^-1*(a+b*arccsch(c+d*x))^2*log(1+(f+(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+f^-1*(a+b*arccsch(c+d*x))^2*log(1+(f+-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+(1/2)*b^2*f^-1*Polylogarithms.polylog(3,exp(2*arccsch(c+d*x)))+-1*f^-1*(a+b*arccsch(c+d*x))^2*log(1+-1*exp(2*arccsch(c+d*x)))+-2*b^2*f^-1*Polylogarithms.polylog(3,-1*(f+(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+-2*b^2*f^-1*Polylogarithms.polylog(3,-1*(f+-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+-1*b*f^-1*(a+b*arccsch(c+d*x))*Polylogarithms.polylog(2,exp(2*arccsch(c+d*x)))+2*b*f^-1*(a+b*arccsch(c+d*x))*Polylogarithms.polylog(2,-1*(f+(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+2*b*f^-1*(a+b*arccsch(c+d*x))*Polylogarithms.polylog(2,-1*(f+-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x))))
@test integrate((a+b*arccsch(c+d*x))^2*(e+f*x)^-2, x) == :(-1*f^-1*(a+b*arccsch(c+d*x))^2*(e+f*x)^-1+d*f^-1*(a+b*arccsch(c+d*x))^2*(d*e+-1*c*f)^-1+-2*d*b^2*(d*e+-1*c*f)^-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-1//2*Polylogarithms.polylog(2,-1*(f+-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+2*d*b^2*(d*e+-1*c*f)^-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-1//2*Polylogarithms.polylog(2,-1*(f+(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+-2*b*d*(d*e+-1*c*f)^-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-1//2*(a+b*arccsch(c+d*x))*log(1+(f+-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+2*b*d*(d*e+-1*c*f)^-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-1//2*(a+b*arccsch(c+d*x))*log(1+(f+(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x))))
@test integrate((a+b*arccsch(c+d*x))^2*(e+f*x)^-3, x) == :(-1//2*f^-1*(a+b*arccsch(c+d*x))^2*(e+f*x)^-2+(1/2)*d^2*f^-1*(a+b*arccsch(c+d*x))^2*(d*e+-1*c*f)^-2+-2*b^2*d^2*(d*e+-1*c*f)^-2*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-1//2*Polylogarithms.polylog(2,-1*(f+-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+2*b^2*d^2*(d*e+-1*c*f)^-2*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-1//2*Polylogarithms.polylog(2,-1*(f+(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+f*b^2*d^2*(d*e+-1*c*f)^-2*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-1*log(f+(c+d*x)^-1*(d*e+-1*c*f))+b^2*d^2*f^2*(d*e+-1*c*f)^-2*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-3//2*Polylogarithms.polylog(2,-1*(f+-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+-1*b^2*d^2*f^2*(d*e+-1*c*f)^-2*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-3//2*Polylogarithms.polylog(2,-1*(f+(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+-2*b*d^2*(d*e+-1*c*f)^-2*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-1//2*(a+b*arccsch(c+d*x))*log(1+(f+-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+2*b*d^2*(d*e+-1*c*f)^-2*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-1//2*(a+b*arccsch(c+d*x))*log(1+(f+(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+b*d^2*f^2*(d*e+-1*c*f)^-2*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-3//2*(a+b*arccsch(c+d*x))*log(1+(f+-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+-1*b*d^2*f^2*(d*e+-1*c*f)^-2*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-3//2*(a+b*arccsch(c+d*x))*log(1+(f+(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^(1/2))^-1*(d*e+-1*c*f)*exp(arccsch(c+d*x)))+-1*b*f*d^2*(1+(c+d*x)^-2)^(1/2)*(f+(c+d*x)^-1*(d*e+-1*c*f))^-1*(d*e+-1*c*f)^-1*(d^2*e^2+f^2*(1+c^2)+-2*c*d*e*f)^-1*(a+b*arccsch(c+d*x)))
@test integrate(x^3*arccsch(x^(1/2)), x) == :(1//4*x^4*arccsch(x^(1/2))+-3//20*x^(1/2)*(-1x)^-1//2*(-1+-1x)^5//2+-1//4*x^(1/2)*(-1x)^-1//2*(-1+-1x)^(1/2)+-1//4*x^(1/2)*(-1x)^-1//2*(-1+-1x)^3//2+-1//28*x^(1/2)*(-1x)^-1//2*(-1+-1x)^7//2)
@test integrate(x^2*arccsch(x^(1/2)), x) == :(1//3*x^3*arccsch(x^(1/2))+1//3*x^(1/2)*(-1x)^-1//2*(-1+-1x)^(1/2)+1//15*x^(1/2)*(-1x)^-1//2*(-1+-1x)^5//2+2//9*x^(1/2)*(-1x)^-1//2*(-1+-1x)^3//2)
@test integrate(x*arccsch(x^(1/2)), x) == :((1/2)*x^2*arccsch(x^(1/2))+-1//2*x^(1/2)*(-1x)^-1//2*(-1+-1x)^(1/2)+-1//6*x^(1/2)*(-1x)^-1//2*(-1+-1x)^3//2)
@test integrate(arccsch(x^(1/2)), x) == :(x*arccsch(x^(1/2))+x^(1/2)*(-1x)^-1//2*(-1+-1x)^(1/2))
@test integrate(x^-1*arccsch(x^(1/2)), x) == :(arccsch(x^(1/2))^2+-1*Polylogarithms.polylog(2,exp(2*arccsch(x^(1/2))))+-2*arccsch(x^(1/2))*log(1+-1*exp(2*arccsch(x^(1/2)))))
@test integrate(x^-2*arccsch(x^(1/2)), x) == :(-1*x^-1*arccsch(x^(1/2))+(1/2)*x^-1//2*(-1x)^-1//2*(-1+-1x)^(1/2)+-1//2*x^(1/2)*(-1x)^-1//2*arctan((-1+-1x)^(1/2)))
@test integrate(x^-3*arccsch(x^(1/2)), x) == :(-1//2*x^-2*arccsch(x^(1/2))+-3//16*x^-1//2*(-1x)^-1//2*(-1+-1x)^(1/2)+1//8*x^-3//2*(-1x)^-1//2*(-1+-1x)^(1/2)+3//16*x^(1/2)*(-1x)^-1//2*arctan((-1+-1x)^(1/2)))
@test integrate(x^-4*arccsch(x^(1/2)), x) == :(-1//3*x^-3*arccsch(x^(1/2))+-5//48*x^(1/2)*(-1x)^-1//2*arctan((-1+-1x)^(1/2))+-5//72*x^-3//2*(-1x)^-1//2*(-1+-1x)^(1/2)+1//18*x^-5//2*(-1x)^-1//2*(-1+-1x)^(1/2)+5//48*x^-1//2*(-1x)^-1//2*(-1+-1x)^(1/2))
@test integrate(arccsch(x^-1), x) == :(-1*(1+x^2)^(1/2)+x*arcsinh(x))
@test integrate(x^-1*arccsch(a*x^n), x) == :((1/2)*n^-1*arccsch(a*x^n)^2+-1//2*n^-1*Polylogarithms.polylog(2,exp(2*arccsch(a*x^n)))+-1*n^-1*arccsch(a*x^n)*log(1+-1*exp(2*arccsch(a*x^n))))
@test integrate(x^-1*arccsch(a*x^5), x) == :(-1//10*Polylogarithms.polylog(2,exp(2*arccsch(a*x^5)))+1//10*arccsch(a*x^5)^2+-1//5*arccsch(a*x^5)*log(1+-1*exp(2*arccsch(a*x^5))))
@test integrate(arccsch(c*exp(a+b*x)), x) == :((1/2)*b^-1*arccsch(c*exp(a+b*x))^2+-1//2*b^-1*Polylogarithms.polylog(2,exp(2*arccsch(c*exp(a+b*x))))+-1*b^-1*arccsch(c*exp(a+b*x))*log(1+-1*exp(2*arccsch(c*exp(a+b*x)))))
@test integrate(x^4*exp(arccsch(a*x)), x) == :(1//4*a^-1*x^4+1//5*x^5*(1+a^-2*x^-2)^3//2+-2//15*a^-2*x^3*(1+a^-2*x^-2)^3//2)
@test integrate(x^3*exp(arccsch(a*x)), x) == :(-1//8*a^-4*arctanh((1+a^-2*x^-2)^(1/2))+1//3*a^-1*x^3+1//4*x^4*(1+a^-2*x^-2)^(1/2)+1//8*a^-2*x^2*(1+a^-2*x^-2)^(1/2))
@test integrate(x^2*exp(arccsch(a*x)), x) == :((1/2)*a^-1*x^2+1//3*x^3*(1+a^-2*x^-2)^3//2)
@test integrate(x*exp(arccsch(a*x)), x) == :(x*a^-1+(1/2)*a^-2*arctanh((1+a^-2*x^-2)^(1/2))+(1/2)*x^2*(1+a^-2*x^-2)^(1/2))
@test integrate(x^-1*exp(arccsch(a*x)), x) == :(-1*(1+a^-2*x^-2)^(1/2)+-1*a^-1*x^-1+arctanh((1+a^-2*x^-2)^(1/2)))
@test integrate(x^-2*exp(arccsch(a*x)), x) == :(-1//2*a*arccsch(a*x)+-1//2*a^-1*x^-2+-1//2*x^-1*(1+a^-2*x^-2)^(1/2))
@test integrate(x^-3*exp(arccsch(a*x)), x) == :(-1//3*a^-1*x^-3+-1//3*a^2*(1+a^-2*x^-2)^3//2)
@test integrate(x^-4*exp(arccsch(a*x)), x) == :(-1//4*a^-1*x^-4+-1//4*x^-3*(1+a^-2*x^-2)^(1/2)+1//8*a^3*arccsch(a*x)+-1//8*a^2*x^-1*(1+a^-2*x^-2)^(1/2))
@test integrate(x^-5*exp(arccsch(a*x)), x) == :(-1//5*a^-1*x^-5+-1//5*a^4*(1+a^-2*x^-2)^5//2+1//3*a^4*(1+a^-2*x^-2)^3//2)
@test integrate(x^4*exp(arccsch(a*x^2)), x) == :(1//3*a^-1*x^3+1//5*x^5*(1+a^-2*x^-4)^(1/2)+2//5*x*a^-2*(1+a^-2*x^-4)^(1/2)+-2//5*a^-2*x^-1*(1+a^-2*x^-4)^(1/2)*(a+x^-2)^-1+-1//5*a^-7//2*((a+x^-2)^-2*(a^2+x^-4))^(1/2)*(1+a^-2*x^-4)^-1//2*(a+x^-2)*Elliptic.F(2*arccot(x*a^(1/2)),1/2)+2//5*a^-7//2*((a+x^-2)^-2*(a^2+x^-4))^(1/2)*(1+a^-2*x^-4)^-1//2*(a+x^-2)*Elliptic.E(2*arccot(x*a^(1/2)),1/2))
@test integrate(x^3*exp(arccsch(a*x^2)), x) == :((1/2)*a^-1*x^2+1//4*a^-2*arctanh((1+a^-2*x^-4)^(1/2))+1//4*x^4*(1+a^-2*x^-4)^(1/2))
@test integrate(x^2*exp(arccsch(a*x^2)), x) == :(x*a^-1+1//3*x^3*(1+a^-2*x^-4)^(1/2)+-1//3*a^-5//2*((a+x^-2)^-2*(a^2+x^-4))^(1/2)*(1+a^-2*x^-4)^-1//2*(a+x^-2)*Elliptic.F(2*arccot(x*a^(1/2)),1/2))
@test integrate(x*exp(arccsch(a*x^2)), x) == :(a^-1*log(x)+(1/2)*x^2*(1+a^-2*x^-4)^(1/2)+-1//2*a^-1*arccsch(a*x^2))
@test integrate(exp(arccsch(a*x^2)), x) == :(x*(1+a^-2*x^-4)^(1/2)+-1*a^-1*x^-1+-2*x^-1*(1+a^-2*x^-4)^(1/2)*(a+x^-2)^-1+-1*a^-3//2*((a+x^-2)^-2*(a^2+x^-4))^(1/2)*(1+a^-2*x^-4)^-1//2*(a+x^-2)*Elliptic.F(2*arccot(x*a^(1/2)),1/2)+2*a^-3//2*((a+x^-2)^-2*(a^2+x^-4))^(1/2)*(1+a^-2*x^-4)^-1//2*(a+x^-2)*Elliptic.E(2*arccot(x*a^(1/2)),1/2))
@test integrate(x^-1*exp(arccsch(a*x^2)), x) == :((1/2)*arctanh((1+a^-2*x^-4)^(1/2))+-1//2*(1+a^-2*x^-4)^(1/2)+-1//2*a^-1*x^-2)
@test integrate(x^-2*exp(arccsch(a*x^2)), x) == :(-1//3*a^-1*x^-3+-1//3*x^-1*(1+a^-2*x^-4)^(1/2)+-1//3*a^-1//2*((a+x^-2)^-2*(a^2+x^-4))^(1/2)*(1+a^-2*x^-4)^-1//2*(a+x^-2)*Elliptic.F(2*arccot(x*a^(1/2)),1/2))
@test integrate(x^-3*exp(arccsch(a*x^2)), x) == :(-1//4*a*arccsch(a*x^2)+-1//4*a^-1*x^-4+-1//4*x^-2*(1+a^-2*x^-4)^(1/2))
@test integrate(x^-4*exp(arccsch(a*x^2)), x) == :(-1//5*a^-1*x^-5+-1//5*x^-3*(1+a^-2*x^-4)^(1/2)+-2*a^2*x^-1*(1+a^-2*x^-4)^(1/2)*(5a+5*x^-2)^-1+-1//5*a^(1/2)*((a+x^-2)^-2*(a^2+x^-4))^(1/2)*(1+a^-2*x^-4)^-1//2*(a+x^-2)*Elliptic.F(2*arccot(x*a^(1/2)),1/2)+2//5*a^(1/2)*((a+x^-2)^-2*(a^2+x^-4))^(1/2)*(1+a^-2*x^-4)^-1//2*(a+x^-2)*Elliptic.E(2*arccot(x*a^(1/2)),1/2))
@test integrate(x^-5*exp(arccsch(a*x^2)), x) == :(-1//6*a^-1*x^-6+-1//6*a^2*(1+a^-2*x^-4)^3//2)
@test integrate(x^4*exp(2*arccsch(a*x)), x) == :(1//5*x^5+-1//4*a^-5*arctanh((1+a^-2*x^-2)^(1/2))+2//3*a^-2*x^3+(1/2)*a^-1*x^4*(1+a^-2*x^-2)^(1/2)+1//4*a^-3*x^2*(1+a^-2*x^-2)^(1/2))
@test integrate(x^3*exp(2*arccsch(a*x)), x) == :(1//4*x^4+a^-2*x^2+2//3*a^-1*x^3*(1+a^-2*x^-2)^3//2)
@test integrate(x^2*exp(2*arccsch(a*x)), x) == :(1//3*x^3+a^-3*arctanh((1+a^-2*x^-2)^(1/2))+2*x*a^-2+a^-1*x^2*(1+a^-2*x^-2)^(1/2))
@test integrate(x*exp(2*arccsch(a*x)), x) == :((1/2)*x^2+-2*a^-2*arccsch(a*x)+2*a^-2*log(x)+2*x*a^-1*(1+a^-2*x^-2)^(1/2))
@test integrate(exp(2*arccsch(a*x)), x) == :(x+-2*a^-1*(1+a^-2*x^-2)^(1/2)+-2*a^-2*x^-1+2*a^-1*arctanh((1+a^-2*x^-2)^(1/2)))
@test integrate(x^-1*exp(2*arccsch(a*x)), x) == :(-1*arccsch(a*x)+-1*a^-2*x^-2+-1*a^-1*x^-1*(1+a^-2*x^-2)^(1/2)+log(x))
@test integrate(x^-3*exp(2*arccsch(a*x)), x) == :(-1//2*x^-2+-1//2*a^-2*x^-4+1//4*a^2*arccsch(a*x)+-1//2*a^-1*x^-3*(1+a^-2*x^-2)^(1/2)+-1//4*a*x^-1*(1+a^-2*x^-2)^(1/2))
@test integrate(x^-4*exp(2*arccsch(a*x)), x) == :(-1//3*x^-3+-2//5*a^-2*x^-5+-2//5*a^3*(1+a^-2*x^-2)^5//2+2//3*a^3*(1+a^-2*x^-2)^3//2)
@test integrate(x^-5*exp(2*arccsch(a*x)), x) == :(-1//4*x^-4+-1//3*a^-2*x^-6+-1//8*a^4*arccsch(a*x)+-1//3*a^-1*x^-5*(1+a^-2*x^-2)^(1/2)+-1//12*a*x^-3*(1+a^-2*x^-2)^(1/2)+1//8*a^3*x^-1*(1+a^-2*x^-2)^(1/2))
@test integrate(x^5*(1+c^2*x^2)^-1*exp(arccsch(c*x)), x) == :(c^-6*arctan(c*x)+-1*x*c^-5+1//3*c^-3*x^3+3//8*c^-6*arctanh((1+c^-2*x^-2)^(1/2))+-3//8*c^-4*x^2*(1+c^-2*x^-2)^(1/2)+1//4*c^-2*x^4*(1+c^-2*x^-2)^(1/2))
@test integrate(x^4*(1+c^2*x^2)^-1*exp(arccsch(c*x)), x) == :((1/2)*c^-3*x^2+-1//2*c^-5*log(1+c^2*x^2)+-2//3*x*c^-4*(1+c^-2*x^-2)^(1/2)+1//3*c^-2*x^3*(1+c^-2*x^-2)^(1/2))
@test integrate(x^3*(1+c^2*x^2)^-1*exp(arccsch(c*x)), x) == :(x*c^-3+-1*c^-4*arctan(c*x)+-1//2*c^-4*arctanh((1+c^-2*x^-2)^(1/2))+(1/2)*c^-2*x^2*(1+c^-2*x^-2)^(1/2))
@test integrate(x^2*(1+c^2*x^2)^-1*exp(arccsch(c*x)), x) == :((1/2)*c^-3*log(1+c^2*x^2)+x*c^-2*(1+c^-2*x^-2)^(1/2))
@test integrate(x*(1+c^2*x^2)^-1*exp(arccsch(c*x)), x) == :(c^-2*arctan(c*x)+c^-2*arctanh((1+c^-2*x^-2)^(1/2)))
@test integrate((1+c^2*x^2)^-1*exp(arccsch(c*x)), x) == :(c^-1*log(x)+-1*c^-1*arccsch(c*x)+-1//2*c^-1*log(1+c^2*x^2))
@test integrate(x^-1*(1+c^2*x^2)^-1*exp(arccsch(c*x)), x) == :(-1*(1+c^-2*x^-2)^(1/2)+-1*arctan(c*x)+-1*c^-1*x^-1)
@test integrate(x^-2*(1+c^2*x^2)^-1*exp(arccsch(c*x)), x) == :((1/2)*c*arccsch(c*x)+(1/2)*c*log(1+c^2*x^2)+-1*c*log(x)+-1//2*c^-1*x^-2+-1//2*x^-1*(1+c^-2*x^-2)^(1/2))
@test integrate(x^-3*(1+c^2*x^2)^-1*exp(arccsch(c*x)), x) == :(c*x^-1+c^2*(1+c^-2*x^-2)^(1/2)+c^2*arctan(c*x)+-1//3*c^-1*x^-3+-1//3*c^2*(1+c^-2*x^-2)^3//2)
@test integrate((d*x+a*d*b^-1)^-1*arccsch(a+b*x), x) == :((1/2)*d^-1*arccsch(a+b*x)^2+-1//2*d^-1*Polylogarithms.polylog(2,exp(2*arccsch(a+b*x)))+-1*d^-1*arccsch(a+b*x)*log(1+-1*exp(2*arccsch(a+b*x))))
@test integrate(x^3*arccsch(a+b*x^4), x) == :(1//4*b^-1*arctanh((1+(a+b*x^4)^-2)^(1/2))+1//4*b^-1*(a+b*x^4)*arccsch(a+b*x^4))
@test integrate(x^(-1+n)*arccsch(a+b*x^n), x) == :(b^-1*n^-1*arctanh((1+(a+b*x^n)^-2)^(1/2))+b^-1*n^-1*(a+b*x^n)*arccsch(a+b*x^n))

max_stars_repo_path: test/7 Inverse hyperbolic functions/7.6 Inverse hyperbolic cosecant/7.6.2 Inverse hyperbolic cosecant functions.jl


text:
using LearnByCheating
using Test

@testset "LearnByCheating.jl" begin
    # Write your tests here.
end


text:
# addprocs(3)

using Base.Test
using Compat
# using IncrementalInference

println("[TEST] out of module evalPotential...")
include("TestModuleFunctions.jl")
println("Success")

print("[TEST] Ensure memory return is working properly...")
include("typeReturnMemRef.jl")
println("Success")

print("[TEST] Ensure converter types can be run from extending namespaces...")
include("saveconvertertypes.jl")
println("Success")

println("[TEST] packing converters work...")
include("testgraphpackingconverters.jl")
println("Success")

println("[TEST] generic root finding by numeric solve of residual functions...")
include("testNumericRootGenericRandomized.jl")
println("Success")

println("[TEST] GenericWrapParam functors...")
include("testGenericWrapParam.jl")
println("Success")

println("[TEST] with simple local constraint examples Odo, Obsv2...")
include("testlocalconstraintexamples.jl")
println("Success")

println("[TEST] partial constraints...")
include("testpartialconstraint.jl")
println("Success")

println("[TEST] partial constraints...")
include("testnullhypothesis.jl")
println("Success")

println("[TEST] with local Graphs.jl dictionary and arrays only (multicore)...")
include("fourdoortest.jl")
println("Success")

println("[TEST] saving to and loading from .jld file")
savejld(fg, file="tempfg.jld" )
fgu = loadjld( file="tempfg.jld" )
Base.rm("tempfg.jld")
println("Success")

#


text:
# Autogenerated wrapper script for Qt5Svg_jll for aarch64-linux-gnu-cxx11
export libqt5svg

using Qt5Base_jll
JLLWrappers.@generate_wrapper_header("Qt5Svg")
JLLWrappers.@declare_library_product(libqt5svg, "libQt5Svg.so.5")
function __init__()
    JLLWrappers.@generate_init_header(Qt5Base_jll)
    JLLWrappers.@init_library_product(
        libqt5svg,
        "lib/libQt5Svg.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Literate

# ! Convert scripts to markdown using Literate
files = ["examples.jl", "api.jl"]

function lit_to_md(file)
    examples_path = joinpath("docs", "src", "examples")
    out_md_path = joinpath("docs", "src")
    Literate.markdown(
        joinpath(examples_path, file),
        out_md_path;
        documenter=true,
        execute=true
    )
end

map(lit_to_md, files)


text:
using Compat
using Base.Test
using PropertyGraph
using UUID

type TestGraphLoader <: GraphLoader
	# A graph loader that provides edges and vertices from a pair of dictionaries
	# used to test that registered GraphLoader is used as expected by PropertyGraph.jl
	graph::Graph

	created::Dict{UUID.Uuid, Container}
	sourceedges::Dict{UUID.Uuid, Edge}
	sourcevertices::Dict{UUID.Uuid, Vertex}

	function TestGraphLoader(g::Graph)
		# Constructs a GraphLoader for the specified Graph

		gl = new()
		gl.graph = g

		gl.created = Dict{UUID.Uuid,Container}()
		gl.sourceedges = Dict{UUID.Uuid,Edge}()
		gl.sourcevertices = Dict{UUID.Uuid,Vertex}()

		return gl
	end
end

function PropertyGraph.loadmissingdata(loader::TestGraphLoader, e::Edge)
	itemloaded::Bool = false

	if haskey(loader.sourceedges, e.id)
		edge = loader.sourceedges[e.id]
		clearproperties!(e)
		setpropertyvalues!(e, edge.attachedproperties)

		e.partiallyloaded = false

		loadmissingdata(loader, e.head)
		loadmissingdata(loader, e.tail)

		itemloaded = true
	end

	return itemloaded
end

function partialcopy(loader::TestGraphLoader, edge::Edge)
	if haskey(loader.created, edge.id)
		return loader.created[edge.id]
	else
		head = partialcopy(loader, edge.head)
		tail = partialcopy(loader, edge.tail)

		e = Edge(tail, head)
		e.partiallyloaded = true
		e.typelabel = edge.typelabel
		e.id = edge.id

		loader.created[edge.id] = e

		return e
	end
end

function partialcopy(loader::TestGraphLoader, vertex::Vertex)
	if haskey(loader.created, vertex.id)
		return loader.created[vertex.id]
	else
		v = Vertex()
		v.partiallyloaded = true
		v.typelabel = vertex.typelabel
		v.id = vertex.id

		loader.created[vertex.id] = v

		return v
	end
end

function PropertyGraph.loadmissingdata(loader::TestGraphLoader, v::Vertex)
	itemloaded::Bool = false

	if haskey(loader.sourcevertices, v.id)
		vertex = loader.sourcevertices[v.id]
		clearproperties!(v)
		setpropertyvalues!(v, vertex.attachedproperties)

		v.incomingedges = Set{Container}()
		v.outgoingedges = Set{Container}()

		for e in vertex.incomingedges
			push!(v.incomingedges, partialcopy(loader, e))
		end

		for e in vertex.outgoingedges
			push!(v.outgoingedges, partialcopy(loader, e))
		end

		v.partiallyloaded = false
		itemloaded = true
	end

	return itemloaded
end

function PropertyGraph.loadmissingdata(loader::TestGraphLoader, items::Set{Container})
	# given a partially loaded set of containers (edges or vertices), loading missing data
	allloaded::Bool = true

	for i in items
		if i.partiallyloaded
			itemloaded::Bool = loadmissingdata(loader, i)
			allloaded = allloaded && itemloaded
		end
	end

	return allloaded
end

function PropertyGraph.loadmissingdata(loader::GraphLoader, results::Set{CompositeQueryResultItem})
	# given a partially loaded set of containers (edges or vertices), load missing data
	allloaded::Bool = true

	for r in results
		if r.item.partiallyloaded
			allloaded = allloaded && loadmissingdata(loader, r.item)
		end
	end

	return allloaded
end

function PropertyGraph.loadvertices(loader::TestGraphLoader)
	vertices = Set{Vertex}()

	for v in values(loader.sourcevertices)
		push(vertices, partialcopy(loader, v))
	end

	return vertices
end

function PropertyGraph.loadvertices(loader::TestGraphLoader, where::Function)
	vertices = Set{Vertex}()

	for v in values(loader.sourcevertices)
		if where(CompositeQueryResultItem(v))
			push!(vertices, partialcopy(loader, v))
		end
	end

	return vertices
end

function PropertyGraph.loadvertices(loader::TestGraphLoader, whereExp::Expr)
	throw(NotImplementedException()) # must be implemented for specific graph loader
end

function PropertyGraph.loadedges(loader::TestGraphLoader)
	edges = Set{Edge}{}

	for e in values(loader.sourceedges)
		push!(edges, partialcopy(loader, e))
	end

	return edges
end

function PropertyGraph.loadedges(loader::TestGraphLoader, where::Function)
	edges = Set{Edge}{}

	for e in values(loader.sourceedges)
		if where(CompositeQueryResultItem(e))
			push!(edges, partialcopy(loader, e))
		end
	end

	return edges
end

function PropertyGraph.loadedges(loader::TestGraphLoader, whereExp::Expr)
	throw(NotImplementedException()) # must be implemented for specific graph loader
end