text:
function insupport(d::TruncatedMvDistribution{D,R,S}, x::AbstractArray) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
    insupport(d.untruncated,x) && intruncationregion(d.region,x)
end

function rand(d::TruncatedMvDistribution{D,R,S}) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
    rand_naive(d)
end

function rand_naive(d::TruncatedMvDistribution{D,R,S}) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
    while true
        candidate = rand(d.untruncated)
        intruncationregion(d.region,candidate) && return candidate
    end
end

function length(d::TruncatedMvDistribution{D,R,S}) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
    length(d.untruncated)
end

function size(d::TruncatedMvDistribution{D,R,S}) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
    size(d.untruncated)
end

function pdf(d::TruncatedMvDistribution{D,R,S},x::AbstractArray; worst_tol = 1e-3) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
    d.state.tp_err < worst_tol || compute_tp(d)
    if intruncationregion(d.region,x)
        pdf(d.untruncated, x) / d.state.tp
    else
        zeros(length(d))
    end
end

function mean(d::TruncatedMvDistribution{D,R,S}; worst_tol = 1e-3 ) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
    d.state.μ_err < worst_tol || compute_mean(d)
    return d.state.μ
end

function cov(d::TruncatedMvDistribution{D,R,S}; worst_tol = 1e-3 ) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
    d.state.Σ_err < worst_tol || compute_cov(d)
    return d.state.Σ
end

function moment(d::TruncatedMvDistribution{D,R,S}, k::Vector{Int} ; worst_tol = 1e-3 ) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
    compute_moment(d,k)
end

function tp(d::TruncatedMvDistribution{D,R,S}; worst_tol = 1e-3) where {D <: MultivariateDistribution, R <: TruncationRegion, S <: TruncatedMvDistributionState}
    d.state.tp_err < worst_tol || compute_tp(d)
    d.state.tp
end


text:
module Move_AllVisualObjects

using Modia3D

filename = joinpath(Modia3D.path, "objects", "fish", "SiameseTiger0.3ds")
#filename = joinpath(Modia3D.path, "objects","engine", "crank", "crank.obj")
#file = FileShape.convexFile(filename; scaleFactor=MVector{3,Float64}(4.0,4.0,4.0))

# Material for Visualization
mat1 = Modia3D.Material(color="Blue", transparency=0.5)
mat2 = Modia3D.Material(transparency=0.2)

# Shapes
coordsys = Modia3D.CoordinateSystem(0.7)
r1 = [4.5, 0.0, 0.0]
dr = [-1.5,0.0,0.0]

@assembly VisualObjects begin
   world   = Modia3D.Object3D()
   shape1  = Modia3D.Object3D(world  , Modia3D.Box(0.9,0.5,0.3    , material=mat1); r=r1, fixed=false)
   shape2  = Modia3D.Object3D(shape1 , Modia3D.Sphere(0.7         , material=mat1); r=dr, fixed=false)
   shape3  = Modia3D.Object3D(shape2 , Modia3D.Cylinder(0.5,0.8   , material=mat1); r=dr, fixed=false)
   shape4  = Modia3D.Object3D(shape3 , Modia3D.Cone(0.3,0.7       , material=mat1); r=dr, fixed=false)
   shape5  = Modia3D.Object3D(shape4 , Modia3D.Capsule(0.4,0.45   , material=mat1); r=dr, fixed=false)
   shape6  = Modia3D.Object3D(shape5 , Modia3D.Spring(0.3,0.7     , material=mat1); r=dr, fixed=false)
   shape7  = Modia3D.Object3D(shape1 , Modia3D.GearWheel(0.5,0.8  , material=mat1); r=[0.0,0.0,-2.0], fixed=false)
   shape8  = Modia3D.Object3D(shape7 , Modia3D.Pipe(0.5,0.8       , material=mat1); r=dr, fixed=false)
   shape9  = Modia3D.Object3D(shape8 , Modia3D.Beam(0.4,0.5,0.3   , material=mat1); r=dr, fixed=false)
   shape10 = Modia3D.Object3D(shape9 , Modia3D.Grid(1.0,0.6)                      ; r=dr, fixed=false)
   shape11 = Modia3D.Object3D(shape10, coordsys                                   ; r=dr, fixed=false)
   shape12 = Modia3D.Object3D(shape11, Modia3D.FileMesh(filename,4, material=mat2); r=dr, fixed=false)
end
visualObjects = VisualObjects(sceneOptions=Modia3D.SceneOptions(visualizeFrames=true,
                                                                defaultFrameLength=0.7,
                                                                enableContactDetection=false))

# Kinematic simulation
Modia3D.initAnalysis!(visualObjects)

for time = range(0.0, stop=2.0, length=101)
  # update positional degrees of freedom
  r1[3] = 2*time
  Modia3D.set_r!(visualObjects.shape1,r1)
  Modia3D.updatePosition!(visualObjects)

  # Visualize shapes
  Modia3D.visualize!(visualObjects,time)
end

Modia3D.closeAnalysis!(visualObjects)

println("... success of Move_AllVisualObjects.jl!")

end

ext: jl


text:
export cuarrays, cubytes
using CUDA: CuArray

# Recursively search for CuArrays based on deepcopy_internal

cuarrays(x, c=CuArray[], d=IdDict{Any,Bool}()) = (_cuarrays(x,c,d); c)

_cuarrays(x::CuArray, c::Vector{CuArray}, d::IdDict{Any,Bool}) =
    if !haskey(d,x); d[x] = true; push!(c,x); if hasfield(typeof(x),:parent) && x.parent !== nothing; _cuarrays(x.parent,c,d); end; end

_cuarrays(x::Union{Module,String,Symbol,Core.MethodInstance,Method,GlobalRef,DataType,Union,UnionAll,Task,Regex},
          c::Vector{CuArray}, d::IdDict{Any,Bool}) = return

_cuarrays(x::Tuple, c::Vector{CuArray}, d::IdDict{Any,Bool}) =
    for xi in x; _cuarrays(xi, c, d); end

_cuarrays(x::Core.SimpleVector, c::Vector{CuArray}, d::IdDict{Any,Bool}) =
    if !haskey(d,x); d[x] = true; for xi in x; _cuarrays(xi, c, d); end; end

_cuarrays(x::Union{Dict,IdDict}, c::Vector{CuArray}, d::IdDict{Any,Bool}) =
    if !haskey(d,x); d[x] = true; for (k,v) in x; _cuarrays(k, c, d); _cuarrays(v, c, d); end; end

function _cuarrays(x::Array{T}, c::Vector{CuArray}, d::IdDict{Any,Bool}) where T
    if !isbitstype(T) && !haskey(d,x)
        d[x] = true
        for i = 1:(length(x)::Int)
            if ccall(:jl_array_isassigned, Cint, (Any, Csize_t), x, i-1) != 0
                xi = ccall(:jl_arrayref, Any, (Any, Csize_t), x, i-1)
                if !isbits(xi)
                    _cuarrays(xi, c, d)
                end
            end
        end
    end
end

function _cuarrays(@nospecialize(x), c::Vector{CuArray}, d::IdDict{Any,Bool})
    T = typeof(x)::DataType
    nf = nfields(x)
    (isbitstype(T) || nf == 0) && return
    if haskey(d, x)
        return
    end
    if T.mutable
        d[x] = true
    end
    for i in 1:nf
        if isdefined(x,i)
            _cuarrays(getfield(x,i), c, d)
        end
    end
end

cubytes(x)=sum(sizeof(a) for a in cuarrays(x))


text:
abstract type Action end
abstract type ByPermutations <: Action end
abstract type ByLinearTransformation <: Action end

abstract type InducedActionHomomorphism{A,T} end
#=
implements:
* features(action_hom)
* action(hom::InducedActionHomomorphism) → Action
* action(hom::InducedActionHomomorphism, g::GroupElement, feature) → the action of `g` on `feature`
=#

function action(
    hom::InducedActionHomomorphism{<:ByPermutations},
    g::GroupElement,
    v::AbstractVector,
)
    return v^induce(hom, g)
end

function action(
    hom::InducedActionHomomorphism{<:ByLinearTransformation},
    g::GroupElement,
    v::AbstractVector,
)
    return induce(hom, g) * v
end

function induce(ac::Action, hom::InducedActionHomomorphism, g::GroupElement)
    return throw(
        "No fallback is provided for $(typeof(ac)). You need to implement `induce(::$(typeof(ac)), ::$(typeof(hom)), ::$(typeof(g)))`.",
    )
end

induce(hom::InducedActionHomomorphism, g::GroupElement) =
    induce(action(hom), hom, g)

decompose(x::Any, hom::InducedActionHomomorphism) = throw(
    "No fallback is provided for $(typeof(x)). You need to implement `decompose(::$(typeof(x)), ::$(typeof(hom)))`.",
)
coeff_type(ac::ByLinearTransformation) = throw(
    "No fallback is provided for $(typeof(ac)). You need to implement `_coeff_type(::$(typeof(ac)))`.",
)

# Exceeding typemax(UInt16) here would mean e.g. that you're trying to block-diagonalize
# an SDP constraint of size 65535×65535, which is highly unlikely ;)
_int_type(::InducedActionHomomorphism) = UInt16

function induce(
    ac::ByLinearTransformation,
    hom::InducedActionHomomorphism,
    g::GroupElement,
)
    n = length(features(hom))

    I = Int[]
    J = Int[]
    V = coeff_type(ac)[]

    for (i, f) in enumerate(features(hom))
        k = action(action(hom), g, f)
        idcs, vals = decompose(k, hom)
        append!(I, fill(i, length(idcs)))
        append!(J, idcs)
        append!(V, vals)
    end
    return sparse(I, J, V)
end

# ala permutation action
function decompose(m::T, hom::InducedActionHomomorphism{A,T}) where {A,T}
    return [hom[m]], [one(coeff_type(action(hom)))]
end

struct ExtensionHomomorphism{A,T,I,V} <: InducedActionHomomorphism{A,T}
    ac::A
    features::V # supports linear indexing
    reversef::Dict{T,I}
end

# see the comment about UInt16 above
ExtensionHomomorphism(ac::Action, features) =
    ExtensionHomomorphism{UInt16}(ac, features)

function ExtensionHomomorphism{I}(ac::Action, features) where {I<:Integer}
    @assert typemax(I) >= length(features) "Use wider Integer type!"
    reversef =
        Dict{eltype(features),I}(f => idx for (idx, f) in pairs(features))
    return ExtensionHomomorphism(ac, features, reversef)
end

_int_type(::ExtensionHomomorphism{A,T,I}) where {A,T,I} = I

Base.getindex(ehom::ExtensionHomomorphism, i::Integer) = ehom.features[i]
Base.getindex(ehom::ExtensionHomomorphism{A,T}, f::T) where {A,T} =
    ehom.reversef[f]

# interface:
features(hom::ExtensionHomomorphism) = hom.features
action(hom::ExtensionHomomorphism) = hom.ac

function induce(::ByPermutations, hom::ExtensionHomomorphism, g::GroupElement)
    I = _int_type(hom)
    return Perm(vec(I[hom[action(action(hom), g, f)] for f in features(hom)]))
end

struct CachedExtensionHomomorphism{A,T,G,H,E<:InducedActionHomomorphism{A,T}} <:
       InducedActionHomomorphism{A,T}
    ehom::E
    cache::Dict{G,H}
end

for f in (:action, :features, :_int_type)
    @eval $f(h::CachedExtensionHomomorphism) = $f(h.ehom)
end

Base.getindex(h::CachedExtensionHomomorphism, x::Any) = h.ehom[x]

function CachedExtensionHomomorphism(
    G::Group,
    action::Action,
    basis;
    precompute = false,
)
    hom = ExtensionHomomorphism(action, basis)
    S = typeof(induce(hom, one(G)))
    chom = CachedExtensionHomomorphism{eltype(G),S}(hom)
    if precompute
        # to make sure that the future access to chom is read-only, i.e. thread-safe
        # one may choose to precompute the images
        for g in G
            induce(chom, g)
        end
    end
    return chom
end

CachedExtensionHomomorphism{G,H}(hom::InducedActionHomomorphism) where {G,H} =
    CachedExtensionHomomorphism(hom, Dict{G,H}())

function induce(ac::Action, hom::CachedExtensionHomomorphism, g::GroupElement)
    if !haskey(hom.cache, g)
        hom.cache[g] = induce(ac, hom.ehom, g)
    end
    return hom.cache[g]
end
#disambiguation:
function induce(
    ac::ByLinearTransformation,
    hom::CachedExtensionHomomorphism,
    g::GroupElement,
)
    if !haskey(hom.cache, g)
        hom.cache[g] = induce(ac, hom.ehom, g)
    end
    return hom.cache[g]
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 421056b1-5fd8-4a43-b9f9-29ab6655f01e
begin
    import Pkg
    Pkg.activate(mktempdir())

    # Registered packages
    pkgnames = [
        "ModelingToolkit",
        "OrdinaryDiffEq",
        "SpecialFunctions",
        "Plots",
        "Printf",
        "RecipesBase",
        "Romberg",
    ]
    pkgspecs = [Pkg.PackageSpec(name = pkgname) for pkgname in pkgnames]
    Pkg.add(pkgspecs)

    # Unregistered packages
    Pkg.add(Pkg.PackageSpec(url="https://github.com/SciML/MethodOfLines.jl"))

    using ModelingToolkit, SpecialFunctions, MethodOfLines, OrdinaryDiffEq
    using Plots, RecipesBase, Printf, Romberg
end

# ╔═╡ aef53a73-6a16-4259-bb35-1d9d735313a0
md"# Mathematical modeling of reactive melt infiltration (RMI)"

# ╔═╡ 8626400d-31e7-4b19-af63-47115eb893c7
md"## Setting up Julia language"

# ╔═╡ bf11b1e1-0459-4c9e-b2a2-cf325f6b7552
md"## Definition of variables and parameters"

# ╔═╡ 51457f69-ed29-4969-bb80-c531ccc7c695
md"#### Variables"

# ╔═╡ 43cd59d9-ad46-42d7-a7ba-443a2bd4f8de
@parameters t x

# ╔═╡ ff609196-76d6-4aba-9971-81835f800271
# c: vapor concentration
# s: surface coverage
# u: liquid volume fraction
# d: thickness of the product at liquid/substrate interface
@variables c(..) s(..) u(..) d(..) T(..)

# ╔═╡ 90323a52-6ec4-45bd-9ee0-11efbab19746
md"#### Model parameters"

# ╔═╡ 8c562f00-b437-4455-b3c2-a04dacb4e77f
# δ: thickness of VW-like crystals at vapor/substrate interface
# r: crystal growth rate
# v: evaporation rate
# D_gas: vapor diffusion coeffient
# w: wettability parameter
# k: permeability parameter
# D_fm: Reaction constant for FM-like reaction at liquid/substrate interface
@parameters δ, r, v, D_gas, w, k, D_fm

# ╔═╡ 6e282286-5366-469d-9f5f-96b85ab2f629
md"#### Physical constants"

# ╔═╡ 5b76a22f-2cd8-4402-b134-9349bd02b232
λ_s = 0.1 # solid thermal conductivity

# ╔═╡ 2ac5c5f2-9882-420a-94f3-17cc9614bf53
λ_l = λ_s # liquid thermal conductivity

# ╔═╡ e8aca2ec-14e4-4dcd-83e1-de898a312eb4
H = 0.1 # enthalpy of reaction

# ╔═╡ 8a2af572-4a87-4f30-a031-56efbda43a70
md"## Governing equations"

# ╔═╡ 3b75ae00-df1a-497b-8553-cec55597fb46
Dx = Differential(x)

# ╔═╡ 547313be-fc25-4868-9331-428976ae14ad
Dxx = Dx^2

# ╔═╡ ede6012f-9cbd-4dbc-b10e-c55b14336c38
Dt = Differential(t)

# ╔═╡ 84ff092f-bde9-4696-91c9-5a5b97a20cde
ev(c, u) = v * (1-c) * u # evaporation model

# ╔═╡ f9d3a857-5fca-4329-ad6e-3a0ab4a53331
vw(c, s) = r * c * (1-s) # VW-like crystal growth at vapor/substrate interface

# ╔═╡ 471e3c14-33ae-4103-8586-41385073f230
cp(w, s) = 1 + erf((s-1)/w) # capillary pressure model

# ╔═╡ 9486d7f2-aa1a-4bdf-828d-6c4dbd1744c3
fm(u, d) = u * D_fm / d # FM-like crystal growth at liquid/substrate interface

# ╔═╡ 8a706e21-9a86-45f6-b329-92f3439d44b1
λ(u) = (1-u) * λ_s + u * λ_l # thermal conductivity

# ╔═╡ ca6af1e9-2bf5-4751-aa95-8b685f4bb987
eqs = [
    # Gas concentration
    Dt(c(t,x)) ~ D_gas * Dxx(c(t,x)) - δ * vw(c(t,x), s(t,x)) + ev(c(t,x), u(t,x)),

    # Surface coverage
    Dt(s(t,x)) ~ vw(c(t,x), s(t,x)),

    # Liquid volume fraction
    Dt(u(t,x)) ~ Dx(k*cp(w, s(t,x))*Dx(u(t,x))),

    # Thickness of the product at liquid/substrate interface
    Dt(d(t,x)) ~ fm(u(t,x), d(t,x)),

    # Temperature
    Dt(T(t,x)) ~ Dx(λ(u(t,x))*Dx(T(t,x))) + H * fm(u(t,x), d(t,x)),
]

# ╔═╡ 59d5fcdd-28cf-44ce-8bd3-bd790c095852
plot([s -> cp(w, s) for w in [10, 1, 0.1]], 0, 1,
     labels = map(w -> "w = " * w, ["10" "1" "0.1"]),
     xlabel = "Surface coverage",
     ylabel = "capillary pressure")

# ╔═╡ a7599967-4633-4b2e-a355-8fcaa0446580
md"## Boundary conditions"

# ╔═╡ b084b8bd-2fe6-4d38-9090-34613ef41842
bcs = [
    c(0,x) ~ 0.0,
    Dx(c(t,0)) ~ 0.0,
    Dx(c(t,1)) ~ 0.0,

    s(0,x) ~ 0.0,
    Dx(s(t,0)) ~ 0.0,
    Dx(s(t,1)) ~ 0.0,

    u(0,x) ~ 0.0,
    u(t,0) ~ 1.0,
    Dx(u(t,1)) ~ 0.0,

    d(0,x) ~ δ,
    Dx(d(t,0)) ~ 0.0,
    Dx(d(t,1)) ~ 0.0,

    T(0,x) ~ 1,
    Dx(T(t,0)) ~ 0.0,
    Dx(T(t,1)) ~ 0.0,
]

# ╔═╡ c8a0e235-9c6d-4d3d-a6c0-75a96db1139b
md"## Solving PDEs"

# ╔═╡ 25ad6ce3-86fc-478f-b957-66857422e2b9
domains = [t ∈ (0.0, 1.0),
           x ∈ (0.0, 1.0)]

# ╔═╡ 1b5c9ffb-70d5-45d8-9042-4892df89008c
@named pde = PDESystem(expand_derivatives.(eqs),
                       bcs,
                       domains,
                       [t,x],
                       [c(t,x), s(t,x), u(t,x), d(t,x), T(t,x)],
                       [δ => 0.1,       # thickness of VW-like reaction
                        r => 10,        # crystal growth rate
                        v => 40,        # evaporation rate
                        D_gas => 0.01,  # vapor diffusion coeffient
                        w => 0.3,       # wettability parameter
                        k => 10,        # permeability parameter
                        D_fm => 10]);   # Reaction constant for FM-like reaction

# ╔═╡ 4860b5f5-dcb8-40ca-b378-702bcfb40738
dx = 0.01

# ╔═╡ 3f2f9dba-ce74-460a-ae59-78a331de9424
disc = MOLFiniteDifference([x => dx], t);

# ╔═╡ 967e96b7-b29b-480d-9383-2faf8709b13c
prob = discretize(pde, disc)

# ╔═╡ c62ef778-aa47-4079-a580-e9cb5234fe73
sol = solve(prob, Tsit5());

# ╔═╡ b75aac9a-ede8-44db-980e-df35a48b18a3
md"## Plotting results"

# ╔═╡ 93588680-d150-47d1-8cde-8f2c39441881
@recipe function plot(range::StepRangeLen, sol::ODESolution)
    ts = sol.t
    N = length(ts)
    xs = collect(range[2:end-1])
    M = length(xs)

    layout := @layout [c s
                       u l
                       d T]

    t_out = collect(0:0.2:1)
    ls = zeros(length(t_out))
    k = 1
    for i in 1:N
        if ts[i] ≥ t_out[k]
            @series begin
                subplot := 1
                ylabel --> "c"
                xlabel --> "x"
                ylims --> (0, 1)
                label --> @sprintf "t = %1.1f" ts[i]
                xs, sol.u[i][1:M]
            end
            @series begin
                subplot := 2
                ylabel --> "s"
                xlabel --> "x"
                ylims --> (0, 1)
                label --> @sprintf "t = %1.1f" ts[i]
                xs, sol.u[i][M+1:2M]
            end
            @series begin
                subplot := 3
                ylabel --> "u"
                xlabel --> "x"
                ylims --> (0, 1)
                label --> @sprintf "t = %1.1f" ts[i]
                xs, sol.u[i][2M+1:3M]
            end
            @series begin
                subplot := 5
                ylabel --> "d"
                xlabel --> "x"
                label --> @sprintf "t = %1.1f" ts[i]
                xs, sol.u[i][3M+1:4M]
            end
            @series begin
                subplot := 6
                ylabel --> "T"
                xlabel --> "x"
                label --> @sprintf "t = %1.1f" ts[i]
                xs, sol.u[i][4M+1:5M]
            end

            # Calculate infiltration length
            ls[k], _ = romberg(range[2:end-1], sol.u[i][2M+1:3M])

            k += 1
        end
    end
    @series begin
        subplot := 4
        ylabel --> "Infiltration length"
        xlabel --> "t"
        legend --> false
        markershape --> :circle
        t_out, ls
    end
    xlims --> (0, 1)
end

# ╔═╡ eac0b5d5-866d-4a59-b33a-e81eb5c11a4e
plot(0:dx:1, sol, size=(600,800))

# ╔═╡ Cell order:
# ╟─aef53a73-6a16-4259-bb35-1d9d735313a0
# ╟─8626400d-31e7-4b19-af63-47115eb893c7
# ╠═421056b1-5fd8-4a43-b9f9-29ab6655f01e
# ╟─bf11b1e1-0459-4c9e-b2a2-cf325f6b7552
# ╟─51457f69-ed29-4969-bb80-c531ccc7c695
# ╠═43cd59d9-ad46-42d7-a7ba-443a2bd4f8de
# ╠═ff609196-76d6-4aba-9971-81835f800271
# ╟─90323a52-6ec4-45bd-9ee0-11efbab19746
# ╠═8c562f00-b437-4455-b3c2-a04dacb4e77f
# ╟─6e282286-5366-469d-9f5f-96b85ab2f629
# ╠═5b76a22f-2cd8-4402-b134-9349bd02b232
# ╠═2ac5c5f2-9882-420a-94f3-17cc9614bf53
# ╠═e8aca2ec-14e4-4dcd-83e1-de898a312eb4
# ╟─8a2af572-4a87-4f30-a031-56efbda43a70
# ╠═3b75ae00-df1a-497b-8553-cec55597fb46
# ╠═547313be-fc25-4868-9331-428976ae14ad
# ╠═ede6012f-9cbd-4dbc-b10e-c55b14336c38
# ╠═ca6af1e9-2bf5-4751-aa95-8b685f4bb987
# ╠═84ff092f-bde9-4696-91c9-5a5b97a20cde
# ╠═f9d3a857-5fca-4329-ad6e-3a0ab4a53331
# ╠═471e3c14-33ae-4103-8586-41385073f230
# ╠═9486d7f2-aa1a-4bdf-828d-6c4dbd1744c3
# ╠═8a706e21-9a86-45f6-b329-92f3439d44b1
# ╠═59d5fcdd-28cf-44ce-8bd3-bd790c095852
# ╟─a7599967-4633-4b2e-a355-8fcaa0446580
# ╠═b084b8bd-2fe6-4d38-9090-34613ef41842
# ╟─c8a0e235-9c6d-4d3d-a6c0-75a96db1139b
# ╠═25ad6ce3-86fc-478f-b957-66857422e2b9
# ╠═1b5c9ffb-70d5-45d8-9042-4892df89008c
# ╠═4860b5f5-dcb8-40ca-b378-702bcfb40738
# ╠═3f2f9dba-ce74-460a-ae59-78a331de9424
# ╠═967e96b7-b29b-480d-9383-2faf8709b13c
# ╠═c62ef778-aa47-4079-a580-e9cb5234fe73
# ╟─b75aac9a-ede8-44db-980e-df35a48b18a3
# ╠═93588680-d150-47d1-8cde-8f2c39441881
# ╠═eac0b5d5-866d-4a59-b33a-e81eb5c11a4e

max_stars_repo_name: hasundue/pluto

max_issues_repo_name: hasundue/pluto

max_forks_repo_name: hasundue/pluto


text:
using ..Atmos
using ..TurbulenceClosures
using ..TurbulenceConvection

# Method definitions for diagnostics collection for all the components
# of `AtmosModel`.

const AtmosComponentTypes = Union{
    AbstractMoistureModel,
    PrecipitationModel,
    RadiationModel,
    TracerModel,
    TurbulenceClosureModel,
    TurbulenceConvectionModel,
}

dv_PointwiseDiagnostic(
    ::AtmosConfigType,
    ::PointwiseDiagnostic,
    ::AtmosComponentTypes,
    ::AtmosModel,
    ::States,
    ::AbstractFloat,
    ::Dict{Symbol, Any},
) = 0
dv_HorizontalAverage(
    ::AtmosConfigType,
    ::HorizontalAverage,
    ::AtmosComponentTypes,
    ::AtmosModel,
    ::States,
    ::AbstractFloat,
    ::Dict{Symbol, Any},
) = 0


text:
using DataAssim

include("resample_test.jl")


text:
module GLFWBackend

using InteractiveUtils
using ModernGL
using GLFW
using ..LibNuklear

# constants
const NK_GLFW_TEXT_MAX = 256
const NK_GLFW_DOUBLE_CLICK_LO = 0.02
const NK_GLFW_DOUBLE_CLICK_HI = 0.2
const NK_GLFW3_DEFAULT = 0
const NK_GLFW3_INSTALL_CALLBACKS = 1

export NK_GLFW3_DEFAULT, NK_GLFW3_INSTALL_CALLBACKS

include("types.jl")
export nk_glfw_vertex, nk_glfw_device, nk_glfw

include("callbacks.jl")
export nk_glfw3_char_callback, nk_gflw3_scroll_callback, nk_glfw3_mouse_button_callback

include("impl.jl")
export nk_glfw3_device_create, nk_glfw3_device_destroy, nk_glfw3_init, nk_glfw3_shutdown
export nk_glfw3_font_stash_begin, nk_glfw3_font_stash_end
export nk_glfw3_device_upload_atlas, nk_glfw3_render, nk_glfw3_new_frame
export nk_glfw3_create_texture, nk_glfw3_update_texture, nk_glfw3_delete_texture, create_nk_image

function __init__()
    global glfw = nk_glfw()
end
export glfw

end # module


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE
using Test, Tensors
dtime = 0.25
# mat = Material(IdealPlastic, tuple())
# mat.properties.youngs_modulus = 200.0e3
# mat.properties.poissons_ratio = 0.3
# mat.properties.yield_stress = 100.0
parameters = IdealPlasticParameterState(youngs_modulus = 200.0e3,
                                        poissons_ratio = 0.3,
                                        yield_stress = 100.0)
mat = IdealPlastic(parameters=parameters)
times = [mat.drivers.time]
stresses = [copy(tovoigt(mat.variables.stress))]
stresses_expected = [[50.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                     [100.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                     [100.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                     [50.0, 0.0, 0.0, 0.0, 0.0, 0.0],
                     [-100.0, 0.0, 0.0, 0.0, 0.0, 0.0]]
dstrain11 = 1e-3*dtime
strains_expected = [[dstrain11, -0.3*dstrain11, -0.3*dstrain11, 0.0, 0.0, 0.0],
                    [2*dstrain11, -0.3*dstrain11*2, -0.3*dstrain11*2, 0.0, 0.0, 0.0],
                    [3*dstrain11, -0.3*dstrain11*2 - 0.5*dstrain11, -0.3*dstrain11*2 - 0.5*dstrain11, 0.0, 0.0, 0.0],
                    [2*dstrain11, -0.3*dstrain11 - 0.5*dstrain11, -0.3*dstrain11 - 0.5*dstrain11, 0.0, 0.0, 0.0],
                    [-2*dstrain11, 0.3*dstrain11*2, 0.3*dstrain11*2, 0.0, 0.0, 0.0]]
dtimes = [dtime, dtime, dtime, dtime, 1.0]
dstrains11 = [dstrain11, dstrain11, dstrain11, -dstrain11, -4*dstrain11]
for i in 1:length(dtimes)
    dstrain11 = dstrains11[i]
    dtime = dtimes[i]
    uniaxial_increment!(mat, dstrain11, dtime)
    # uniaxial_increment!(mat, dstrain11, dtime; dstrain = copy(tovoigt(mat.ddrivers.strain))*dstrain11/mat.ddrivers.strain[1,1]*dtime/mat.ddrivers.time)
    # mat.time += mat.dtime
    # mat.strain .+= mat.dstrain
    # mat.stress .+= mat.dstress
    update_material!(mat)
    # push!(times, mat.drivers.time)
    # push!(stresses, copy(tovoigt(mat.variables.stress)))
    #@info(tovoigt(mat.variables.stress), stresses_expected[i])
    @test isapprox(tovoigt(mat.variables.stress), stresses_expected[i])
    @test isapprox(tovoigt(mat.drivers.strain; offdiagscale=2.0), strains_expected[i])
end


text:
using Coincheck
using Base.Test

using HTTP

@testset "HttpUtil" begin
    importall Coincheck.HttpUtil
    @test convert_to_url("http://localhost", "foo") == "http://localhost/foo"
    @test convert_to_url("http://localhost", "foo", Dict("one" => 1, "two" => 2)) == "http://localhost/foo?two=2&one=1"
    @test convert_to_url("http://localhost/", "foo") == "http://localhost/foo"
    @test convert_to_url("http://localhost", "/foo") == "http://localhost/foo"
    @test convert_to_url("http://localhost", "foo/") == "http://localhost/foo"
end

@testset "call_http_api Tests" begin
    server = HTTP.Server((req, rep) -> begin
        req.uri == HTTP.URI("/success") && return HTTP.Response("""{"success": true}""")
        req.uri == HTTP.URI("/failure") && return HTTP.Response("""failure""")
        HTTP.Response("")
    end)
    client = Coincheck.HttpClient("http://127.0.0.1:8081")

    try
        @spawn HTTP.serve(server, HTTP.IPv4(127,0,0,1), 8081)

        @test Coincheck.call_http_api(client, "success") == Dict("success" => true)
        @test_throws ErrorException Coincheck.call_http_api(client, "failure")
    finally
        put!(server.in, 9)
    end

end


text:
#=
Buttons on right side determine what is displayed on canvas1
=#

function rb1_cb(widgetptr::Ptr,user_data::Tuple{Gui_Handles,Int64})

    han, event_id = user_data

    if han.c_right_top != event_id
        han.c_right_top=event_id

        if event_id==3
            if han.c_right_bottom != 5
                setproperty!(han.rb2[5],:active,true)
            end
        elseif event_id==4
            if han.c_right_bottom != 6
                setproperty!(han.rb2[6],:active,true)
            end
        else
            if (han.c_right_bottom == 5)|(han.c_right_bottom == 6)
                setproperty!(han.rb2[8],:active,true)
            end
        end
    end

    if han.c_right_top == 1
	set_slider(han,16)
    elseif han.c_right_top == 2
	set_slider(han,32)
    elseif (han.c_right_top == 3)|(han.c_right_top == 4)
	set_slider(han,64)
    end

    clear_c(han)
    nothing
end

function rb2_cb(widgetptr::Ptr,user_data::Tuple{Gui_Handles,Int64})

    han, event_id = user_data

    if han.c_right_bottom != event_id
        han.c_right_bottom=event_id

        if event_id==5
            if han.c_right_top!=3
                setproperty!(han.rb1[3],:active,true)
            end
        elseif event_id==6
            if han.c_right_top != 4
                setproperty!(han.rb1[4],:active,true)
            end
        elseif event_id == 7
            ctx=Gtk.getgc(han.c)
            han.c.back.ptr = CairoRGBSurface(width(ctx),height(ctx)).ptr
            han.c.backcc = CairoContext(han.c.back)
        else
            if (han.c_right_top == 3)|(han.c_right_top == 4)
                setproperty!(han.rb1[1],:active,true)
            end
        end
    end

    clear_c(han)
    nothing
end

function set_slider(han::Gui_Handles,chan_num::Int64)
    han.num16=div(han.sc.spike-1,chan_num)+1
    han.num=rem(han.sc.spike-1,chan_num)+1
    han.chan_per_display=chan_num
    if han.num16*chan_num>length(han.enabled)
        setproperty!(han.adj2,:upper,han.num16*chan_num-length(han.enabled))
        setproperty!(han.adj,:upper,div(length(han.enabled),chan_num)+1)
    else
        setproperty!(han.adj2,:upper,chan_num)
        setproperty!(han.adj,:upper,div(length(han.enabled),chan_num))
    end
    setproperty!(han.adj2, :value, han.num)
    setproperty!(han.adj, :value, han.num16)
    nothing
end


text:
"""
    holdout([rng], N, T, λvals, ρvals, [weights,] xDistsTest, xDistsTrain,
            noiseDist, signalModels, kernelgenerator; showprogress)
    holdout([rng], N, T, λvals, ρvals, [weights,] xDistsTest, νDistsTest,
            xDistsTrain, νDistsTrain, noiseDist, signalModels, kernelgenerator;
            showprogress)

Select λ and ρ via a holdout process.

# Arguments
- `rng::AbstractRNG = Random.GLOBAL_RNG`: Random number generator to use
- `N::Integer`: Number of test points
- `T::Integer`: Number of training points
- `λvals::AbstractVector{<:Real}`: Values of λ to search over \\[nλ\\]
- `ρvals::AbstractVector{<:Real}`: Values of ρ to search over \\[nρ\\]
- `weights::AbstractVector{<:Real}`: Weights for calculating holdout cost
  \\[L\\]; omit if L = 1
- `xDistsTest`: Distributions of latent parameters \\[L\\] or scalar (if L = 1);
  `xDists` can be any object such that `rand(xDists, ::Integer)` is defined (or
  a collection of such objects)
- `νDistsTest`: Distributions of known parameters \\[K\\] or scalar (if K = 1);
  `νDists` can be any object such that `rand(νDists, ::Integer)` is defined (or
  a collection of such objects); omit this parameter if K = 0
- `xDistsTrain`: Distributions of latent parameters \\[L\\] or scalar
  (if L = 1); `xDists` can be any object such that `rand(xDists, ::Integer)` is
  defined (or a collection of such objects)
- `νDistsTrain`: Distributions of known parameters \\[K\\] or scalar (if K = 1);
  `νDists` can be any object such that `rand(νDists, ::Integer)` is defined (or
  a collection of such objects); omit this parameter if K = 0
- `noiseDist`: Distribution of noise (assumes same noise distribution for both
  real and imaginary channels in complex case); `noiseDist` can be any object
  such that `rand(noiseDist, ::Integer)` is defined
- `signalModels::Union{<:Function,<:AbstractVector{<:Function}}`: Signal models
  used to generate noiseless data \\[numSignalModels\\]; each signal model
  accepts as inputs L latent parameters (scalars) first, then K known parameters
  (scalars); user-defined parameters (e.g., scan parameters in MRI) should be
  built into the signal model
- `kernelgenerator::Function`: Function that creates a `Kernel` object given a
  vector `Λ` of lengthscales
- `showprogress::Bool = false`: Whether to show progress

## Note
- L is the number of unknown or latent parameters to be estimated
- K is the number of known parameters
- nλ is the number of λ values to try
- nρ is the number of ρ values to try

# Return
- `λ::Real`: Bandwidth scaling parameter
- `ρ::Real`: Tikhonov regularization parameter
- `Ψ::AbstractMatrix{<:Real}`: Holdout costs for λvals and ρvals [nλ,nρ]
"""
function holdout(
    rng::AbstractRNG,
    N::Integer,
    T::Integer,
    λvals::AbstractVector{<:Real},
    ρvals::AbstractVector{<:Real},
    weights::AbstractVector{<:Real},
    xDistsTest::AbstractVector,
    xDistsTrain::AbstractVector,
    noiseDist,
    signalModels::Union{<:Function,<:AbstractArray{<:Function,1}},
    kernelgenerator::Function;
    showprogress::Bool = false
)

    # Make sure length of weights matches xDistsTest and xDistsTrain
    length(weights) == length(xDistsTest) == length(xDistsTrain) ||
        throw(DimensionMismatch("lengths of weights, xDistsTest, and " *
                                "xDistsTrain should be the same"))

    # Generate synthetic test data
    (y, x) = generatenoisydata(rng, N, xDistsTest, noiseDist, signalModels)
    (ytrain, xtrain) = generatenoisydata(rng, T, xDistsTrain, noiseDist,
                                         signalModels)

    # Loop through each value of λ and ρ
    nλ = length(λvals)
    nρ = length(ρvals)
    Ψ  = zeros(nλ, nρ)
    for idxλ = 1:nλ

        showprogress && println("idxλ = $idxλ/$nλ")

        λ = λvals[idxλ]

        # Generate length scales
        if ndims(y) == 1 || size(y, 1) == 1
            Λ = λ * max(mean(abs.(y)), eps()) # scalar (D = 1)
        else
            Λ = λ * max.(dropdims(mean(abs.(y), dims = 2), dims = 2), eps()) # [D]
        end

        # Create the kernel
        kernel = kernelgenerator(Λ)

        for idxρ = 1:nρ

            showprogress && println("    idxρ = $idxρ/$nρ")

            ρ = ρvals[idxρ]

            # Train PERK
            trainData = PERK.krr_train(rng, xtrain, ytrain, kernel, ρ)

            # Run PERK
            xhat = perk(y, trainData, kernel) # [L,N]

            # Calculate Ψ(λ,ρ), the holdout cost
            werr = ((xhat - x) ./ x) .* sqrt.(weights) # [L,N]
            Ψ[idxλ,idxρ] = sqrt(norm(werr) / N)

        end

    end

    # Return values of λ and ρ that minimize Ψ
    (idxλ, idxρ) = Tuple(argmin(Ψ))
    λ = λvals[idxλ]
    ρ = ρvals[idxρ]

    return (λ, ρ, Ψ)

end

function holdout(
    rng::AbstractRNG,
    N::Integer,
    T::Integer,
    λvals::AbstractVector{<:Real},
    ρvals::AbstractVector{<:Real},
    weights::AbstractVector{<:Real},
    xDistsTest::AbstractVector,
    νDistsTest,
    xDistsTrain::AbstractVector,
    νDistsTrain,
    noiseDist,
    signalModels::Union{<:Function,<:AbstractArray{<:Function,1}},
    kernelgenerator::Function;
    showprogress::Bool = false
)

    # Make sure length of weights matches xDistsTest and xDistsTrain
    length(weights) == length(xDistsTest) == length(xDistsTrain) ||
        throw(DimensionMismatch("lengths of weights, xDistsTest, and " *
                                "xDistsTrain should be the same"))

    # Make sure lengths of νDistsTest and νDistsTrain match
    length(νDistsTest) == length(νDistsTrain) ||
        throw(DimensionMismatch("lengths of νDistsTest and νDistsTrain " *
                                "should be the same"))

    # Generate synthetic test data and training data
    (y, x, ν) = generatenoisydata(rng, N, xDistsTest, νDistsTest, noiseDist,
                                  signalModels)
    (ytrain, xtrain, νtrain) = generatenoisydata(rng, T, xDistsTrain,
                                           νDistsTrain, noiseDist, signalModels)

    # Combine y and ν
    q = combine(y, ν) # [D+K,N]
    qtrain = combine(ytrain, νtrain) # [D+K,T]

    # Loop through each value of λ and ρ
    nλ = length(λvals)
    nρ = length(ρvals)
    Ψ  = zeros(nλ, nρ)
    for idxλ = 1:nλ

        showprogress && println("idxλ = $idxλ/$nλ")

        λ = λvals[idxλ]

        # Generate length scales
        Λ = λ * max.(dropdims(mean(abs.(q), dims = 2), dims = 2), eps()) # [D+K]

        # Create the kernel
        kernel = kernelgenerator(Λ)

        for idxρ = 1:nρ

            showprogress && println("    idxρ = $idxρ/$nρ")

            ρ = ρvals[idxρ]

            # Train PERK
            trainData = PERK.krr_train(rng, xtrain, qtrain, kernel, ρ)

            # Run PERK
            xhat = perk(y, ν, trainData, kernel) # [L,N]

            # Calculate Ψ(λ,ρ), the holdout cost
            werr = ((xhat - x) ./ x) .* sqrt.(weights) # [L,N]
            Ψ[idxλ,idxρ] = sqrt(norm(werr) / N)

        end

    end

    # Return values of λ and ρ that minimize Ψ
    (idxλ, idxρ) = Tuple(argmin(Ψ))
    λ = λvals[idxλ]
    ρ = ρvals[idxρ]

    return (λ, ρ, Ψ)

end

function holdout(
    rng::AbstractRNG,
    N::Integer,
    T::Integer,
    λvals::AbstractVector{<:Real},
    ρvals::AbstractVector{<:Real},
    xDistsTest,
    xDistsTrain,
    noiseDist,
    signalModels::Union{<:Function,<:AbstractArray{<:Function,1}},
    kernelgenerator::Function;
    showprogress::Bool = false
)

    weights = [1]
    holdout(rng, N, T, λvals, ρvals, weights, [xDistsTest], [xDistsTrain],
            noiseDist, signalModels, kernelgenerator,
            showprogress = showprogress)

end

function holdout(
    rng::AbstractRNG,
    N::Integer,
    T::Integer,
    λvals::AbstractVector{<:Real},
    ρvals::AbstractVector{<:Real},
    xDistsTest,
    νDistsTest,
    xDistsTrain,
    νDistsTrain,
    noiseDist,
    signalModels::Union{<:Function,<:AbstractArray{<:Function,1}},
    kernelgenerator::Function;
    showprogress::Bool = false
)

    weights = [1]
    holdout(rng, N, T, λvals, ρvals, weights, [xDistsTest], νDistsTest,
            [xDistsTrain], νDistsTrain, noiseDist, signalModels,
            kernelgenerator, showprogress = showprogress)

end

holdout(N::Integer, args...) = holdout(Random.GLOBAL_RNG, N, args...)


text:
function myppmwrite(img, file::String)
    s = open(file, "w")
    write(s, "P6\n")
    n, m = size(img)
    write(s, "$m $n 255\n")
    for i=1:n, j=1:m
        p = img[i,j]
        write(s, UInt64(p))
        write(s, UInt64(p))
        write(s, UInt64(p))
    end
    close(s)
end

max_stars_repo_name: rabarar/juliaFractals

max_issues_repo_name: rabarar/juliaFractals

max_forks_repo_name: rabarar/juliaFractals


text:
module Operators

import LinearAlgebra

using StaticArrays

import ..slab, ..slab_args, ..column, ..column_args
import ..DataLayouts: DataLayouts, Data2D, DataSlab2D
import ..Geometry: Geometry, Covariant12Vector, Contravariant12Vector, ⊗
import ..Spaces: Spaces, Quadratures, AbstractSpace
import ..Topologies
import ..Meshes
import ..Fields: Fields, Field

using ..RecursiveApply

include("spectralelement.jl")
include("numericalflux.jl")
include("finitedifference.jl")
include("remapping.jl")

end # module


text:
###
### **arb.h** -- real numbers
###

### Types, macros and constants

### Memory management
arbcall"void arb_init(arb_t x)"
arbcall"void arb_clear(arb_t x)"
#mo arbcall"arb_ptr _arb_vec_init(slong n)" # clashes with similar method for acb
arbcall"void _arb_vec_clear(arb_ptr v, slong n)"
arbcall"void arb_swap(arb_t x, arb_t y)"
arbcall"slong arb_allocated_bytes(const arb_t x)"
arbcall"slong _arb_vec_allocated_bytes(arb_srcptr vec, slong len)"
#mo arbcall"double _arb_vec_estimate_allocated_bytes(slong len, slong prec)" # clashes with similar method for acb

### Assignment and rounding
arbcall"void arb_set(arb_t y, const arb_t x)"
arbcall"void arb_set_arf(arb_t y, const arf_t x)"
arbcall"void arb_set_si(arb_t y, slong x)"
arbcall"void arb_set_ui(arb_t y, ulong x)"
arbcall"void arb_set_d(arb_t y, double x)"
#ni arbcall"void arb_set_fmpz(arb_t y, const fmpz_t x)"
#ni arbcall"void arb_set_fmpz_2exp(arb_t y, const fmpz_t x, const fmpz_t e)"
arbcall"void arb_set_round(arb_t y, const arb_t x, slong prec)"
#ni arbcall"void arb_set_round_fmpz(arb_t y, const fmpz_t x, slong prec)"
#ni arbcall"void arb_set_round_fmpz_2exp(arb_t y, const fmpz_t x, const fmpz_t e, slong prec)"
#ni arbcall"void arb_set_fmpq(arb_t y, const fmpq_t x, slong prec)"
arbcall"int arb_set_str(arb_t res, const char * inp, slong prec)"
arbcall"char * arb_get_str(const arb_t x, slong n, ulong flags)"

### Assignment of special values
arbcall"void arb_zero(arb_t x)"
arbcall"void arb_one(arb_t f)"
arbcall"void arb_pos_inf(arb_t x)"
arbcall"void arb_neg_inf(arb_t x)"
arbcall"void arb_zero_pm_inf(arb_t x)"
arbcall"void arb_indeterminate(arb_t x)"
arbcall"void arb_zero_pm_one(arb_t x)"
arbcall"void arb_unit_interval(arb_t x)"

### Input and output
#mo arbcall"void arb_print(const arb_t x)" # clashes with Base.print
#ns arbcall"void arb_fprint(FILE * file, const arb_t x)"
arbcall"void arb_printd(const arb_t x, slong digits)"
#ns arbcall"void arb_fprintd(FILE * file, const arb_t x, slong digits)"
arbcall"void arb_printn(const arb_t x, slong digits, ulong flags)"
#ns arbcall"void arb_fprintn(FILE * file, const arb_t x, slong digits, ulong flags)"
arbcall"char * arb_dump_str(const arb_t x)"
arbcall"int arb_load_str(arb_t x, const char * str)"
#ns arbcall"int arb_dump_file(FILE * stream, const arb_t x)"
#ns arbcall"int arb_load_file(arb_t x, FILE * stream)"

### Random number generation
#ns arbcall"void arb_randtest(arb_t x, flint_rand_t state, slong prec, slong mag_bits)"
#ns arbcall"void arb_randtest_exact(arb_t x, flint_rand_t state, slong prec, slong mag_bits)"
#ns arbcall"void arb_randtest_precise(arb_t x, flint_rand_t state, slong prec, slong mag_bits)"
#ns arbcall"void arb_randtest_wide(arb_t x, flint_rand_t state, slong prec, slong mag_bits)"
#ns arbcall"void arb_randtest_special(arb_t x, flint_rand_t state, slong prec, slong mag_bits)"
#ni arbcall"void arb_get_rand_fmpq(fmpq_t q, flint_rand_t state, const arb_t x, slong bits)"

### Radius and interval operations
arbcall"void arb_get_mid_arb(arb_t m, const arb_t x)"
arbcall"void arb_get_rad_arb(arb_t r, const arb_t x)"
arbcall"void arb_add_error_arf(arb_t x, const arf_t err)"
arbcall"void arb_add_error_mag(arb_t x, const mag_t err)"
arbcall"void arb_add_error(arb_t x, const arb_t err)"
arbcall"void arb_add_error_2exp_si(arb_t x, slong e)"
#ni arbcall"void arb_add_error_2exp_fmpz(arb_t x, const fmpz_t e)"
arbcall"void arb_union(arb_t z, const arb_t x, const arb_t y, slong prec)"
arbcall"int arb_intersection(arb_t z, const arb_t x, const arb_t y, slong prec)"
arbcall"void arb_nonnegative_part(arb_t res, const arb_t x)"
arbcall"void arb_get_abs_ubound_arf(arf_t u, const arb_t x, slong prec)"
arbcall"void arb_get_abs_lbound_arf(arf_t u, const arb_t x, slong prec)"
arbcall"void arb_get_ubound_arf(arf_t u, const arb_t x, slong prec)"
arbcall"void arb_get_lbound_arf(arf_t u, const arb_t x, slong prec)"
arbcall"void arb_get_mag(mag_t z, const arb_t x)"
arbcall"void arb_get_mag_lower(mag_t z, const arb_t x)"
arbcall"void arb_get_mag_lower_nonnegative(mag_t z, const arb_t x)"
#ni arbcall"void arb_get_interval_fmpz_2exp(fmpz_t a, fmpz_t b, fmpz_t exp, const arb_t x)"
arbcall"void arb_set_interval_mag(arb_t x, const mag_t a, const mag_t b, slong prec)"
arbcall"void arb_set_interval_arf(arb_t x, const arf_t a, const arf_t b, slong prec)"
arbcall"void arb_set_interval_mpfr(arb_t x, const mpfr_t a, const mpfr_t b, slong prec)"
arbcall"void arb_set_interval_neg_pos_mag(arb_t x, const mag_t a, const mag_t b, slong prec)"
arbcall"void arb_get_interval_arf(arf_t a, arf_t b, const arb_t x, slong prec)"
arbcall"void arb_get_interval_mpfr(mpfr_t a, mpfr_t b, const arb_t x)"
arbcall"slong arb_rel_error_bits(const arb_t x)"
arbcall"slong arb_rel_accuracy_bits(const arb_t x)"
arbcall"slong arb_rel_one_accuracy_bits(const arb_t x)"
arbcall"slong arb_bits(const arb_t x)"
arbcall"void arb_trim(arb_t y, const arb_t x)"
#ni arbcall"int arb_get_unique_fmpz(fmpz_t z, const arb_t x)"
arbcall"void arb_floor(arb_t y, const arb_t x, slong prec)"
arbcall"void arb_ceil(arb_t y, const arb_t x, slong prec)"
#ni arbcall"void arb_get_fmpz_mid_rad_10exp(fmpz_t mid, fmpz_t rad, fmpz_t exp, const arb_t x, slong n)"
#mo arbcall"int arb_can_round_arf(const arb_t x, slong prec, arf_rnd_t rnd)" # clashes with arb_can_round_mpfr
#mo arbcall"int arb_can_round_mpfr(const arb_t x, slong prec, mpfr_rnd_t rnd)" # clashes with arb_can_round_arf

### Comparisons
arbcall"int arb_is_zero(const arb_t x)"
arbcall"int arb_is_nonzero(const arb_t x)"
arbcall"int arb_is_one(const arb_t f)"
arbcall"int arb_is_finite(const arb_t x)"
arbcall"int arb_is_exact(const arb_t x)"
arbcall"int arb_is_int(const arb_t x)"
arbcall"int arb_is_int_2exp_si(const arb_t x, slong e)"
arbcall"int arb_equal(const arb_t x, const arb_t y)"
arbcall"int arb_equal_si(const arb_t x, slong y)"
arbcall"int arb_is_positive(const arb_t x)"
arbcall"int arb_is_nonnegative(const arb_t x)"
arbcall"int arb_is_negative(const arb_t x)"
arbcall"int arb_is_nonpositive(const arb_t x)"
arbcall"int arb_overlaps(const arb_t x, const arb_t y)"
arbcall"int arb_contains_arf(const arb_t x, const arf_t y)"
#ni arbcall"int arb_contains_fmpq(const arb_t x, const fmpq_t y)"
#ni arbcall"int arb_contains_fmpz(const arb_t x, const fmpz_t y)"
arbcall"int arb_contains_si(const arb_t x, slong y)"
arbcall"int arb_contains_mpfr(const arb_t x, const mpfr_t y)"
arbcall"int arb_contains(const arb_t x, const arb_t y)"
arbcall"int arb_contains_int(const arb_t x)"
arbcall"int arb_contains_zero(const arb_t x)"
arbcall"int arb_contains_negative(const arb_t x)"
arbcall"int arb_contains_nonpositive(const arb_t x)"
arbcall"int arb_contains_positive(const arb_t x)"
arbcall"int arb_contains_nonnegative(const arb_t x)"
arbcall"int arb_contains_interior(const arb_t x, const arb_t y)"
arbcall"int arb_eq(const arb_t x, const arb_t y)"
arbcall"int arb_ne(const arb_t x, const arb_t y)"
arbcall"int arb_lt(const arb_t x, const arb_t y)"
arbcall"int arb_le(const arb_t x, const arb_t y)"
arbcall"int arb_gt(const arb_t x, const arb_t y)"
arbcall"int arb_ge(const arb_t x, const arb_t y)"

### Arithmetic
arbcall"void arb_neg(arb_t y, const arb_t x)"
arbcall"void arb_neg_round(arb_t y, const arb_t x, slong prec)"
arbcall"void arb_abs(arb_t y, const arb_t x)"
arbcall"void arb_sgn(arb_t y, const arb_t x)"
arbcall"int arb_sgn_nonzero(const arb_t x)"
arbcall"void arb_min(arb_t z, const arb_t x, const arb_t y, slong prec)"
arbcall"void arb_max(arb_t z, const arb_t x, const arb_t y, slong prec)"
arbcall"void arb_add(arb_t z, const arb_t x, const arb_t y, slong prec)"
arbcall"void arb_add_arf(arb_t z, const arb_t x, const arf_t y, slong prec)"
arbcall"void arb_add_ui(arb_t z, const arb_t x, ulong y, slong prec)"
arbcall"void arb_add_si(arb_t z, const arb_t x, slong y, slong prec)"
#ni arbcall"void arb_add_fmpz(arb_t z, const arb_t x, const fmpz_t y, slong prec)"
#ni arbcall"void arb_add_fmpz_2exp(arb_t z, const arb_t x, const fmpz_t m, const fmpz_t e, slong prec)"
arbcall"void arb_sub(arb_t z, const arb_t x, const arb_t y, slong prec)"
arbcall"void arb_sub_arf(arb_t z, const arb_t x, const arf_t y, slong prec)"
arbcall"void arb_sub_ui(arb_t z, const arb_t x, ulong y, slong prec)"
arbcall"void arb_sub_si(arb_t z, const arb_t x, slong y, slong prec)"
#ni arbcall"void arb_sub_fmpz(arb_t z, const arb_t x, const fmpz_t y, slong prec)"
arbcall"void arb_mul(arb_t z, const arb_t x, const arb_t y, slong prec)"
arbcall"void arb_mul_arf(arb_t z, const arb_t x, const arf_t y, slong prec)"
arbcall"void arb_mul_si(arb_t z, const arb_t x, slong y, slong prec)"
arbcall"void arb_mul_ui(arb_t z, const arb_t x, ulong y, slong prec)"
#ni arbcall"void arb_mul_fmpz(arb_t z, const arb_t x, const fmpz_t y, slong prec)"
arbcall"void arb_mul_2exp_si(arb_t y, const arb_t x, slong e)"
#ni arbcall"void arb_mul_2exp_fmpz(arb_t y, const arb_t x, const fmpz_t e)"
arbcall"void arb_addmul(arb_t z, const arb_t x, const arb_t y, slong prec)"
arbcall"void arb_addmul_arf(arb_t z, const arb_t x, const arf_t y, slong prec)"
arbcall"void arb_addmul_si(arb_t z, const arb_t x, slong y, slong prec)"
arbcall"void arb_addmul_ui(arb_t z, const arb_t x, ulong y, slong prec)"
#ni arbcall"void arb_addmul_fmpz(arb_t z, const arb_t x, const fmpz_t y, slong prec)"
arbcall"void arb_submul(arb_t z, const arb_t x, const arb_t y, slong prec)"
arbcall"void arb_submul_arf(arb_t z, const arb_t x, const arf_t y, slong prec)"
arbcall"void arb_submul_si(arb_t z, const arb_t x, slong y, slong prec)"
arbcall"void arb_submul_ui(arb_t z, const arb_t x, ulong y, slong prec)"
#ni arbcall"void arb_submul_fmpz(arb_t z, const arb_t x, const fmpz_t y, slong prec)"
arbcall"void arb_inv(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_div(arb_t z, const arb_t x, const arb_t y, slong prec)"
arbcall"void arb_div_arf(arb_t z, const arb_t x, const arf_t y, slong prec)"
arbcall"void arb_div_si(arb_t z, const arb_t x, slong y, slong prec)"
arbcall"void arb_div_ui(arb_t z, const arb_t x, ulong y, slong prec)"
#ni arbcall"void arb_div_fmpz(arb_t z, const arb_t x, const fmpz_t y, slong prec)"
#ni arbcall"void arb_fmpz_div_fmpz(arb_t z, const fmpz_t x, const fmpz_t y, slong prec)"
arbcall"void arb_ui_div(arb_t z, ulong x, const arb_t y, slong prec)"
arbcall"void arb_div_2expm1_ui(arb_t z, const arb_t x, ulong n, slong prec)"

### Dot product
arbcall"void arb_dot_precise(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec)"
arbcall"void arb_dot_simple(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec)"
arbcall"void arb_dot(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec)"
arbcall"void arb_approx_dot(arb_t res, const arb_t s, int subtract, arb_srcptr x, slong xstep, arb_srcptr y, slong ystep, slong len, slong prec)"

### Powers and roots
arbcall"void arb_sqrt(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_sqrt_arf(arb_t z, const arf_t x, slong prec)"
#ni arbcall"void arb_sqrt_fmpz(arb_t z, const fmpz_t x, slong prec)"
arbcall"void arb_sqrt_ui(arb_t z, ulong x, slong prec)"
arbcall"void arb_sqrtpos(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_hypot(arb_t z, const arb_t x, const arb_t y, slong prec)"
arbcall"void arb_rsqrt(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_rsqrt_ui(arb_t z, ulong x, slong prec)"
arbcall"void arb_sqrt1pm1(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_root_ui(arb_t z, const arb_t x, ulong k, slong prec)"
#mo arbcall"void arb_root(arb_t z, const arb_t x, ulong k, slong prec)" # alias to arb_root_ui
arbcall"void arb_sqr(arb_t y, const arb_t x, slong prec)"
#ni arbcall"void arb_pow_fmpz_binexp(arb_t y, const arb_t b, const fmpz_t e, slong prec)"
#ni arbcall"void arb_pow_fmpz(arb_t y, const arb_t b, const fmpz_t e, slong prec)"
arbcall"void arb_pow_ui(arb_t y, const arb_t b, ulong e, slong prec)"
arbcall"void arb_ui_pow_ui(arb_t y, ulong b, ulong e, slong prec)"
arbcall"void arb_si_pow_ui(arb_t y, slong b, ulong e, slong prec)"
#ni arbcall"void arb_pow_fmpq(arb_t y, const arb_t x, const fmpq_t a, slong prec)"
arbcall"void arb_pow(arb_t z, const arb_t x, const arb_t y, slong prec)"

### Exponentials and logarithms
arbcall"void arb_log_ui(arb_t z, ulong x, slong prec)"
#ni arbcall"void arb_log_fmpz(arb_t z, const fmpz_t x, slong prec)"
arbcall"void arb_log_arf(arb_t z, const arf_t x, slong prec)"
arbcall"void arb_log(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_log_ui_from_prev(arb_t log_k1, ulong k1, arb_t log_k0, ulong k0, slong prec)"
arbcall"void arb_log1p(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_log_base_ui(arb_t res, const arb_t x, ulong b, slong prec)"
arbcall"void arb_log_hypot(arb_t res, const arb_t x, const arb_t y, slong prec)"
arbcall"void arb_exp(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_expm1(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_exp_invexp(arb_t z, arb_t w, const arb_t x, slong prec)"

### Trigonometric functions
arbcall"void arb_sin(arb_t s, const arb_t x, slong prec)"
arbcall"void arb_cos(arb_t c, const arb_t x, slong prec)"
arbcall"void arb_sin_cos(arb_t s, arb_t c, const arb_t x, slong prec)"
arbcall"void arb_sin_pi(arb_t s, const arb_t x, slong prec)"
arbcall"void arb_cos_pi(arb_t c, const arb_t x, slong prec)"
arbcall"void arb_sin_cos_pi(arb_t s, arb_t c, const arb_t x, slong prec)"
arbcall"void arb_tan(arb_t y, const arb_t x, slong prec)"
arbcall"void arb_cot(arb_t y, const arb_t x, slong prec)"
#ni arbcall"void arb_sin_cos_pi_fmpq(arb_t s, arb_t c, const fmpq_t x, slong prec)"
#ni arbcall"void arb_sin_pi_fmpq(arb_t s, const fmpq_t x, slong prec)"
#ni arbcall"void arb_cos_pi_fmpq(arb_t c, const fmpq_t x, slong prec)"
arbcall"void arb_tan_pi(arb_t y, const arb_t x, slong prec)"
arbcall"void arb_cot_pi(arb_t y, const arb_t x, slong prec)"
arbcall"void arb_sec(arb_t res, const arb_t x, slong prec)"
arbcall"void arb_csc(arb_t res, const arb_t x, slong prec)"
arbcall"void arb_csc_pi(arb_t res, const arb_t x, slong prec)"
arbcall"void arb_sinc(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_sinc_pi(arb_t z, const arb_t x, slong prec)"

### Inverse trigonometric functions
arbcall"void arb_atan_arf(arb_t z, const arf_t x, slong prec)"
arbcall"void arb_atan(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_atan2(arb_t z, const arb_t b, const arb_t a, slong prec)"
arbcall"void arb_asin(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_acos(arb_t z, const arb_t x, slong prec)"

### Hyperbolic functions
arbcall"void arb_sinh(arb_t s, const arb_t x, slong prec)"
arbcall"void arb_cosh(arb_t c, const arb_t x, slong prec)"
arbcall"void arb_sinh_cosh(arb_t s, arb_t c, const arb_t x, slong prec)"
arbcall"void arb_tanh(arb_t y, const arb_t x, slong prec)"
arbcall"void arb_coth(arb_t y, const arb_t x, slong prec)"
arbcall"void arb_sech(arb_t res, const arb_t x, slong prec)"
arbcall"void arb_csch(arb_t res, const arb_t x, slong prec)"

### Inverse hyperbolic functions
arbcall"void arb_atanh(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_asinh(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_acosh(arb_t z, const arb_t x, slong prec)"

### Constants
arbcall"void arb_const_pi(arb_t z, slong prec)"
arbcall"void arb_const_sqrt_pi(arb_t z, slong prec)"
arbcall"void arb_const_log_sqrt2pi(arb_t z, slong prec)"
arbcall"void arb_const_log2(arb_t z, slong prec)"
arbcall"void arb_const_log10(arb_t z, slong prec)"
arbcall"void arb_const_euler(arb_t z, slong prec)"
arbcall"void arb_const_catalan(arb_t z, slong prec)"
arbcall"void arb_const_e(arb_t z, slong prec)"
arbcall"void arb_const_khinchin(arb_t z, slong prec)"
arbcall"void arb_const_glaisher(arb_t z, slong prec)"
arbcall"void arb_const_apery(arb_t z, slong prec)"

### Lambert W function
arbcall"void arb_lambertw(arb_t res, const arb_t x, int flags, slong prec)"

### Gamma function and factorials
arbcall"void arb_rising_ui_bs(arb_t z, const arb_t x, ulong n, slong prec)"
arbcall"void arb_rising_ui_rs(arb_t z, const arb_t x, ulong n, ulong step, slong prec)"
arbcall"void arb_rising_ui_rec(arb_t z, const arb_t x, ulong n, slong prec)"
arbcall"void arb_rising_ui(arb_t z, const arb_t x, ulong n, slong prec)"
arbcall"void arb_rising(arb_t z, const arb_t x, const arb_t n, slong prec)"
#ni arbcall"void arb_rising_fmpq_ui(arb_t z, const fmpq_t x, ulong n, slong prec)"
arbcall"void arb_fac_ui(arb_t z, ulong n, slong prec)"
arbcall"void arb_doublefac_ui(arb_t z, ulong n, slong prec)"
arbcall"void arb_bin_ui(arb_t z, const arb_t n, ulong k, slong prec)"
arbcall"void arb_bin_uiui(arb_t z, ulong n, ulong k, slong prec)"
arbcall"void arb_gamma(arb_t z, const arb_t x, slong prec)"
#ni arbcall"void arb_gamma_fmpq(arb_t z, const fmpq_t x, slong prec)"
#ni arbcall"void arb_gamma_fmpz(arb_t z, const fmpz_t x, slong prec)"
arbcall"void arb_lgamma(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_rgamma(arb_t z, const arb_t x, slong prec)"
arbcall"void arb_digamma(arb_t y, const arb_t x, slong prec)"

### Zeta function
arbcall"void arb_zeta_ui_vec_borwein(arb_ptr z, ulong start, slong num, ulong step, slong prec)"
arbcall"void arb_zeta_ui_asymp(arb_t x, ulong s, slong prec)"
arbcall"void arb_zeta_ui_euler_product(arb_t z, ulong s, slong prec)"
arbcall"void arb_zeta_ui_bernoulli(arb_t x, ulong s, slong prec)"
arbcall"void arb_zeta_ui_borwein_bsplit(arb_t x, ulong s, slong prec)"
arbcall"void arb_zeta_ui_vec(arb_ptr x, ulong start, slong num, slong prec)"
arbcall"void arb_zeta_ui_vec_even(arb_ptr x, ulong start, slong num, slong prec)"
arbcall"void arb_zeta_ui_vec_odd(arb_ptr x, ulong start, slong num, slong prec)"
arbcall"void arb_zeta_ui(arb_t x, ulong s, slong prec)"
arbcall"void arb_zeta(arb_t z, const arb_t s, slong prec)"
arbcall"void arb_hurwitz_zeta(arb_t z, const arb_t s, const arb_t a, slong prec)"

### Bernoulli numbers and polynomials
arbcall"void arb_bernoulli_ui(arb_t b, ulong n, slong prec)"
#ni arbcall"void arb_bernoulli_fmpz(arb_t b, const fmpz_t n, slong prec)"
arbcall"void arb_bernoulli_ui_zeta(arb_t b, ulong n, slong prec)"
arbcall"void arb_bernoulli_poly_ui(arb_t res, ulong n, const arb_t x, slong prec)"
arbcall"void arb_power_sum_vec(arb_ptr res, const arb_t a, const arb_t b, slong len, slong prec)"

### Polylogarithms
arbcall"void arb_polylog(arb_t w, const arb_t s, const arb_t z, slong prec)"
arbcall"void arb_polylog_si(arb_t w, slong s, const arb_t z, slong prec)"

### Other special functions
#ni arbcall"void arb_fib_fmpz(arb_t z, const fmpz_t n, slong prec)"
arbcall"void arb_fib_ui(arb_t z, ulong n, slong prec)"
arbcall"void arb_agm(arb_t z, const arb_t x, const arb_t y, slong prec)"
arbcall"void arb_chebyshev_t_ui(arb_t a, ulong n, const arb_t x, slong prec)"
arbcall"void arb_chebyshev_u_ui(arb_t a, ulong n, const arb_t x, slong prec)"
arbcall"void arb_chebyshev_t2_ui(arb_t a, arb_t b, ulong n, const arb_t x, slong prec)"
arbcall"void arb_chebyshev_u2_ui(arb_t a, arb_t b, ulong n, const arb_t x, slong prec)"
#ni arbcall"void arb_bell_sum_bsplit(arb_t res, const fmpz_t n, const fmpz_t a, const fmpz_t b, const fmpz_t mmag, slong prec)"
#ni arbcall"void arb_bell_sum_taylor(arb_t res, const fmpz_t n, const fmpz_t a, const fmpz_t b, const fmpz_t mmag, slong prec)"
#ni arbcall"void arb_bell_fmpz(arb_t res, const fmpz_t n, slong prec)"
arbcall"void arb_bell_ui(arb_t res, ulong n, slong prec)"
#ni arbcall"void arb_euler_number_fmpz(arb_t res, const fmpz_t n, slong prec)"
arbcall"void arb_euler_number_ui(arb_t res, ulong n, slong prec)"
#ni arbcall"void arb_partitions_fmpz(arb_t res, const fmpz_t n, slong prec)"
arbcall"void arb_partitions_ui(arb_t res, ulong n, slong prec)"

### Internals for computing elementary functions
#ni arbcall"void _arb_atan_taylor_naive(mp_ptr y, mp_limb_t * error, mp_srcptr x, mp_size_t xn, ulong N, int alternating)"
#ni arbcall"void _arb_atan_taylor_rs(mp_ptr y, mp_limb_t * error, mp_srcptr x, mp_size_t xn, ulong N, int alternating)"
#ni arbcall"void _arb_exp_taylor_naive(mp_ptr y, mp_limb_t * error, mp_srcptr x, mp_size_t xn, ulong N)"
#ni arbcall"void _arb_exp_taylor_rs(mp_ptr y, mp_limb_t * error, mp_srcptr x, mp_size_t xn, ulong N)"
#ni arbcall"void _arb_sin_cos_taylor_naive(mp_ptr ysin, mp_ptr ycos, mp_limb_t * error, mp_srcptr x, mp_size_t xn, ulong N)"
#ni arbcall"void _arb_sin_cos_taylor_rs(mp_ptr ysin, mp_ptr ycos, mp_limb_t * error, mp_srcptr x, mp_size_t xn, ulong N, int sinonly, int alternating)"
#ni arbcall"int _arb_get_mpn_fixed_mod_log2(mp_ptr w, fmpz_t q, mp_limb_t * error, const arf_t x, mp_size_t wn)"
#ni arbcall"int _arb_get_mpn_fixed_mod_pi4(mp_ptr w, fmpz_t q, int * octant, mp_limb_t * error, const arf_t x, mp_size_t wn)"
arbcall"slong _arb_exp_taylor_bound(slong mag, slong prec)"
arbcall"void arb_exp_arf_bb(arb_t z, const arf_t x, slong prec, int m1)"
#ni arbcall"void _arb_exp_sum_bs_simple(fmpz_t T, fmpz_t Q, flint_bitcnt_t * Qexp, const fmpz_t x, flint_bitcnt_t r, slong N)"
#ni arbcall"void _arb_exp_sum_bs_powtab(fmpz_t T, fmpz_t Q, flint_bitcnt_t * Qexp, const fmpz_t x, flint_bitcnt_t r, slong N)"
arbcall"void arb_exp_arf_rs_generic(arb_t res, const arf_t x, slong prec, int minus_one)"
#ni arbcall"void _arb_atan_sum_bs_simple(fmpz_t T, fmpz_t Q, flint_bitcnt_t * Qexp, const fmpz_t x, flint_bitcnt_t r, slong N)"
#ni arbcall"void _arb_atan_sum_bs_powtab(fmpz_t T, fmpz_t Q, flint_bitcnt_t * Qexp, const fmpz_t x, flint_bitcnt_t r, slong N)"
arbcall"void arb_atan_arf_bb(arb_t z, const arf_t x, slong prec)"
arbcall"void arb_sin_cos_arf_generic(arb_t s, arb_t c, const arf_t x, slong prec)"
arbcall"void arb_sin_cos_arf_bb(arb_t s, arb_t c, const arf_t x, slong prec)"
arbcall"void arb_sin_cos_wide(arb_t s, arb_t c, const arb_t x, slong prec)"
arbcall"void arb_sin_cos_generic(arb_t s, arb_t c, const arb_t x, slong prec)"

### Vector functions
arbcall"void _arb_vec_zero(arb_ptr vec, slong n)"
arbcall"int _arb_vec_is_zero(arb_srcptr vec, slong len)"
arbcall"int _arb_vec_is_finite(arb_srcptr x, slong len)"
arbcall"void _arb_vec_set(arb_ptr res, arb_srcptr vec, slong len)"
arbcall"void _arb_vec_set_round(arb_ptr res, arb_srcptr vec, slong len, slong prec)"
arbcall"void _arb_vec_swap(arb_ptr vec1, arb_ptr vec2, slong len)"
arbcall"void _arb_vec_neg(arb_ptr B, arb_srcptr A, slong n)"
arbcall"void _arb_vec_sub(arb_ptr C, arb_srcptr A, arb_srcptr B, slong n, slong prec)"
arbcall"void _arb_vec_add(arb_ptr C, arb_srcptr A, arb_srcptr B, slong n, slong prec)"
arbcall"void _arb_vec_scalar_mul(arb_ptr res, arb_srcptr vec, slong len, const arb_t c, slong prec)"
arbcall"void _arb_vec_scalar_div(arb_ptr res, arb_srcptr vec, slong len, const arb_t c, slong prec)"
#ni arbcall"void _arb_vec_scalar_mul_fmpz(arb_ptr res, arb_srcptr vec, slong len, const fmpz_t c, slong prec)"
arbcall"void _arb_vec_scalar_mul_2exp_si(arb_ptr res, arb_srcptr src, slong len, slong c)"
arbcall"void _arb_vec_scalar_addmul(arb_ptr res, arb_srcptr vec, slong len, const arb_t c, slong prec)"
arbcall"void _arb_vec_get_mag(mag_t bound, arb_srcptr vec, slong len, slong prec)"
arbcall"slong _arb_vec_bits(arb_srcptr x, slong len)"
arbcall"void _arb_vec_set_powers(arb_ptr xs, const arb_t x, slong len, slong prec)"
#ni arbcall"void _arb_vec_add_error_arf_vec(arb_ptr res, arf_srcptr err, slong len)"
#ni arbcall"void _arb_vec_add_error_mag_vec(arb_ptr res, mag_srcptr err, slong len)"
arbcall"void _arb_vec_indeterminate(arb_ptr vec, slong len)"
arbcall"void _arb_vec_trim(arb_ptr res, arb_srcptr vec, slong len)"
#ni arbcall"int _arb_vec_get_unique_fmpz_vec(fmpz * res,  arb_srcptr vec, slong len)"


text:
function update_parameters!(net::Net{GPUBackend}, method::Nesterov, learning_rate,
    last_momentum, momentum, param_blob, hist_blob, gradient, data_type)

  # param_blob += -last_momentum* hist_blob (update with vt-1)
  CuBLAS.axpy(net.backend.cublas_ctx, length(hist_blob), convert(data_type, -last_momentum), hist_blob.ptr, 1, param_blob.ptr, 1)

  # hist_blob = last_momentum * hist_blob
  CuBLAS.scal(net.backend.cublas_ctx, length(hist_blob), convert(data_type, last_momentum), hist_blob.ptr, 1)

  # hist_blob = -learning_rate * gradient + hist_blob (calc vt)
  CuBLAS.axpy(net.backend.cublas_ctx, length(hist_blob), convert(data_type, -learning_rate), gradient.ptr, 1, hist_blob.ptr, 1)

  # param_blob += (1+momentum) * hist_blob (update with vt)
  CuBLAS.axpy(net.backend.cublas_ctx, length(hist_blob), convert(data_type, 1 + momentum), hist_blob.ptr, 1, param_blob.ptr, 1)
end


text:
module SPARK

    using Documenter: @doc
    using LinearAlgebra: Diagonal
    using RungeKutta
    using SimpleSolvers

    import GeometricEquations: LODE, LDAE, HDAE, IDAE, PDAE, get_functions

    using ..GeometricBase
    using ..Config

    import ..Solutions: AtomicSolutionPDAE, SolutionVector, update!

    import ..Integrators

    import ..Integrators: PDAEIntegrator, InitialGuessIODE, InitialGuessPODE, Parameters
    import ..Integrators: IDAEIntegratorCache, IntegratorCache, CacheDict, CacheType
    import ..Integrators: AbstractTableau, AbstractCoefficients,
                          @CoefficientsRK, @HeaderCoefficientsRK
    import ..Integrators: create_internal_stage_vector, create_nonlinear_solver,
                          update_vector_fields!, update_solution!, update_multiplier!,
                          initialize!
    import ..Integrators: equation, equations, tableau, timestep,
                          eachdim, eachstage, nstages

    import ..Utils: @define

    export CoefficientsARK, CoefficientsPRK, CoefficientsMRK, CoefficientsIRK,
           CoefficientsSPARK

    export AbstractIntegratorSPARK
    export AbstractTableauSPARK, TableauSPARK

    export IntegratorSPARK, TableauSPARK

    export IntegratorHPARK, TableauHPARK
    export IntegratorVPARK, TableauVPARK

    export IntegratorHSPARK, TableauHSPARK
    export IntegratorHSPARKprimary, TableauHSPARKprimary
    export IntegratorHSPARKsecondary, TableauHSPARKsecondary

    export IntegratorVSPARK, TableauVSPARK
    export IntegratorVSPARKprimary, TableauVSPARKprimary
    export IntegratorVSPARKsecondary, TableauVSPARKsecondary

    export IntegratorSLRK, TableauSLRK

    include("spark/abstract_integrator_spark.jl")
    include("spark/coefficients.jl")

    include("spark/integrators_spark_cache.jl")
    include("spark/integrators_spark_common.jl")
    include("spark/integrators_spark_tableau.jl")
    include("spark/integrators_spark_parameters.jl")

    include("spark/integrators_spark.jl")

    include("spark/integrators_vpark.jl")
    include("spark/integrators_hpark.jl")

    include("spark/integrators_vspark_common.jl")
    include("spark/integrators_vspark.jl")
    include("spark/integrators_vspark_primary.jl")
    include("spark/integrators_vspark_secondary.jl")

    include("spark/integrators_hspark_common.jl")
    include("spark/integrators_hspark.jl")
    include("spark/integrators_hspark_primary.jl")
    include("spark/integrators_hspark_secondary.jl")

    include("spark/integrators_slrk.jl")

end


text:
function _precompile_()
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
    precompile(URIParser.parse_url, (Compat.ASCIIString,))
    precompile(URIParser.parse_authority, (Compat.ASCIIString, Bool,))
    precompile(URIParser.escape_with, (Compat.ASCIIString, Compat.String,))
    precompile(URIParser.is_host_char, (Char,))
    precompile(URIParser.is_url_char, (Char,))
    precompile(URIParser.isnum, (Char,))
    precompile(URIParser.is_mark, (Char,))
    precompile(URIParser.is_userinfo_char, (Char,))
    precompile(URIParser.escape, (Compat.ASCIIString,))
    precompile(URIParser.ishex, (Char,))
end


text:
# ----------------------- #
# Display-esque functions #
# ----------------------- #

const _WINDOW_PLOTLY_CONFIG = """
"""

const _MATHJAX_CONFIG = """
"""

const PLOTLYJS_VERSION = "2.3.0"

const CDN_URL = "https://cdn.plot.ly/plotly-$(PLOTLYJS_VERSION).min.js"

function _requirejs_config()
    """
    $(_WINDOW_PLOTLY_CONFIG)
    $(_MATHJAX_CONFIG)

    """
end

"""
    to_html(
        io::IO,
        p::Plot;
        autoplay::Bool=true,
        include_plotlyjs::Union{String,Missing}="cdn",
        include_mathjax::Union{String,Missing}="cdn",
        post_script::Union{String,Missing}=missing,
        full_html::Bool=true,
        animation_opts::Union{Dict,Missing}=missing,
        default_width::String="100%",
        default_height::String="100%"
    )

- io: IO stream to write to
- p: Plot to save
- autoplay: Should animations start automatically
- include_plotlyjs: How to include plotly.js. Options are
    - cdn: include a
            """
        elseif endswith(including, ".js")
            # check if this is a file
            if embed_plotlyjs
                if !isfile(including)
                    error("cannot embed $(including) file not found")
                end
                # if so, read it and insert inline
                load_plotlyjs = """

                """
            else
                # assume it is a url that can be found
                load_plotlyjs = """
                $(_WINDOW_PLOTLY_CONFIG)

                """
            end
        else  # assume cdn
            load_plotlyjs = """
            $(_WINDOW_PLOTLY_CONFIG)

            """
        end
    end

    # handle mathjax
    load_mathjax = ""
    if isa(include_mathjax, String)
        inmath = lowercase(include_mathjax)
        math_url = ""
        if endswith(inmath, ".js")
            math_url = inmath
        else
            # assume cdn
            math_url = "https://cdnjs.cloudflare.com/ajax/libs/mathjax/2.7.5/MathJax.js"
        end
        load_mathjax = ""
    end

    html_div = """

        $(load_mathjax)
        $(load_plotlyjs)

    """

    if full_html
        return print(io, """
        <meta chartset="utf-8" />

        $(html_div)

        """)
    end
    return print(io, html_div)
end  # function

function savejson(p::Plot, fn::AbstractString)
    ext = split(fn, ".")[end]
    if ext == "json"
        open(f -> print(f, json(p)), fn, "w")
        return p
    else
        msg = "PlotlyBase can only save figures as JSON. For all other"
        msg *= " file types, please use PlotlyJS.jl"
        throw(ArgumentError(msg))
    end
end

# jupyterlab/nteract integration
Base.Multimedia.istextmime(::MIME"application/vnd.plotly.v1+json") = true
function Base.show(io::IO, ::MIME"application/vnd.plotly.v1+json", p::Plot)
    JSON.print(io, p)
end

function Base.show(io::IO, ::MIME"text/plain", p::Plot)
    println(io, """
    data: $(json(map(_describe, p.data), 2))
    layout: "$(_describe(p.layout))"
    """)
end

Base.show(io::IO, ::MIME"text/html", p::Plot; kwargs...) = to_html(io, p; kwargs...)

# integration with vscode and Juno
function Base.show(io::IO, ::MIME"application/prs.juno.plotpane+html", p::Plot)
    show(io, MIME("text/html"), p; include_mathjax="cdn", include_plotlyjs="cdn")
end

Base.show(io::IO, p::Plot) = show(io, MIME("text/plain"), p)

import REPL

"""
opens a browser tab with the given html file
"""
function launch_browser(tmppath::String)
    if Sys.isapple()
        run(`open $tmppath`)
    elseif Sys.iswindows()
        run(`cmd /c start $tmppath`)
    elseif Sys.islinux()
        run(`xdg-open $tmppath`)
    end
end

function Base.display(::REPL.REPLDisplay, p::Plot)
    tmppath = string(tempname(), ".plotlyjs-jl.html")
    open(tmppath, "w") do io
        show(io, MIME("text/html"), p; include_mathjax="cdn", include_plotlyjs="cdn")
    end
    launch_browser(tmppath) # Open the browser
end


text:
using DataFrames, Statistics, GLM, Test

df = CSV.read(joinpath(dirname(@__FILE__), "..", "dataset/Cigar.csv"))
df[:pState] = categorical(df[:State])
df[:pYear] = categorical(df[:Year])

function glm_helper(formula::Formula, df::DataFrame)
    model_response(ModelFrame(formula, df)) - predict(glm(formula, df, Normal(), IdentityLink()))
end
function glm_helper(formula::Formula, df::DataFrame, wts::Symbol)
    model_response(ModelFrame(formula, df)) - predict(glm(formula, df, Normal(), IdentityLink(), wts = convert(Array{Float64}, df[wts])))
end

test = (
    convert(Array{Float64}, partial_out(df, @model(Sales + Price ~ NDI))[1]),
    convert(Array{Float64}, partial_out(df, @model(Sales + Price ~ NDI, fe = pState))[1]),
    convert(Array{Float64}, partial_out(df, @model(Sales + Price ~ 1, fe = pState))[1]),
    convert(Array{Float64}, partial_out(df, @model(Sales + Price ~ 1))[1]),
    mean(convert(Array{Float64}, partial_out(df, @model(Sales + Price ~ NDI, add_mean = true))[1]), dims = 1[1]),
    mean(convert(Array{Float64}, partial_out(df, @model(Sales + Price ~ NDI, fe = pState, add_mean = true))[1]), dims = 1[1]),
    mean(convert(Array{Float64}, partial_out(df, @model(Sales + Price ~ 1, fe = pState, add_mean = true))[1]), dims = 1[1]),
    mean(convert(Array{Float64}, partial_out(df, @model(Sales + Price ~ 1, add_mean = true))[1]), dims = 1[1]),
    convert(Array{Float64}, partial_out(df, @model(Sales + Price ~ NDI, weights = Pop))[1]),
    convert(Array{Float64}, partial_out(df, @model(Sales + Price ~ NDI, fe = pState, weights = Pop))[1]),
    convert(Array{Float64}, partial_out(df, @model(Sales + Price ~ 1, fe = pState, weights = Pop))[1]),
    convert(Array{Float64}, partial_out(df, @model(Sales + Price ~ 1, weights = Pop))[1]),
    )

answer = (
    hcat(glm_helper(@formula(Sales ~ NDI), df), glm_helper(@formula(Price ~ NDI), df)),
    hcat(glm_helper(@formula(Sales ~ NDI + pState), df), glm_helper(@formula(Price ~ NDI + pState), df)),
    hcat(glm_helper(@formula(Sales ~ pState), df), glm_helper(@formula(Price ~ pState), df)),
    hcat(glm_helper(@formula(Sales ~ 1), df), glm_helper(@formula(Price ~ 1), df)),
    hcat(mean(df[:Sales]), mean(df[:Price])),
    hcat(mean(df[:Sales]), mean(df[:Price])),
    hcat(mean(df[:Sales]), mean(df[:Price])),
    hcat(mean(df[:Sales]), mean(df[:Price])),
    hcat(glm_helper(@formula(Sales ~ NDI), df, :Pop), glm_helper(@formula(Price ~ NDI), df, :Pop)),
    hcat(glm_helper(@formula(Sales ~ NDI + pState), df, :Pop), glm_helper(@formula(Price ~ NDI + pState), df, :Pop)),
    hcat(glm_helper(@formula(Sales ~ pState), df, :Pop), glm_helper(@formula(Price ~ pState), df, :Pop)),
    hcat(glm_helper(@formula(Sales ~ 1), df, :Pop), glm_helper(@formula(Price ~ 1), df, :Pop))
    )

for i in 1:12
    @test test[i] ≈ answer[i] atol = 1e-5
end


text:
"""
`assemble_system!` - assemble the matrix for the left hand side of the problem
"""
function assemble_system!(problem::RandomODEProblem)
    # construct the first row of the discrete Laplacian
    problem.A[1,1] = -2
    problem.A[1,2] = 1

    # construct middle rows
    for j ∈ 2:problem.nₓ - 1
        problem.A[j, j- 1] = 1
        problem.A[j, j] = -2
        problem.A[j, j + 1] = 1
    end

    # construct last row
    problem.A[problem.nₓ, problem.nₓ - 1] = 1
    problem.A[problem.nₓ, problem.nₓ] = -2

    # adjust for random potential
    problem.A = spdiagm(problem.V) - problem.A/problem.Δx^2

    return problem
end


text:
using CollegeStratBase, CollegeStratMakie
using Test

csp = CollegeStratMakie;

function fig_test_setup(figName)
    fPath = joinpath(csp.pkg_test_dir(), figName);
    notesPath = csp.fig_notes_path(fPath);
    isfile(fPath)  &&  rm(fPath);
    isfile(notesPath)  &&  rm(notesPath);

    plot_defaults();
    return fPath, notesPath
end

function subplot_layout_test()
    @testset "Subplot layout" begin
        for nPlots = 2 : 10
            nr, nc = subplot_layout(nPlots);
            isValid = true;
            for j = 1 : nPlots
                ir, ic = subplot_pos(j, nr, nc);
                isValid = isValid  && (1 <= ic <= nc)  && (1 <= ir <= nr);
            end
            @test isValid
        end
    end
end

function fig_save_test()
    @testset "Save figure" begin
        fPath, notesPath = fig_test_setup("fig_save_test.pdf");
        p, _ = csp.test_bar_graph();

        fNotes = ["This is a test figure.", "With figure notes"];
        figsave(p, fPath; figNotes = fNotes);
        @test isfile(fPath)
        @test isfile(notesPath)
    end
end

function bar_graph_test()
    @testset "Bar graph" begin
        fPath, notesPath = fig_test_setup("bar_graph_test.pdf");
        n = 5;
        groupLabelV = ["G$j"  for j = 1 : n];
        p, _ = bar_graph(groupLabelV, LinRange(2, 3, n);
            xlabel = "x", ylabel = "y", yerror = LinRange(0.2, 0.3, n));
        figsave(p, fPath);
        @test isfile(fPath);
    end
end

function grouped_bar_test()
    @testset "Grouped bar graph" begin
        fPath, notesPath = fig_test_setup("grouped_bar_test.pdf");
        # 4 groups of 6 bars each
        nr = 4; nc = 6;
        # Bar labels (one for each column)
        grpLabelV = ["Grp $j" for j = 1 : nr];
        labelV = ["Lbl $j" for j = 1 : nc];
        dataM = (10 .*(1:nr)) .+ (1 : nc)';
        yerrorM = 0.1 .* dataM;
        # fig = blank_plot();
        # ax = make_axis(fig; xlabel = "x", ylabel = "y");
            # xticks = (1:nr, ));
        fig, _ = grouped_bar_graph(grpLabelV, dataM;
            yerror = yerrorM,
            legendPos = (1,2), xlabel = "x", ylabel = "y");
        # Cannot create fig[1,2] first
        # not having a title may not work either
        grouped_bar_legend(fig[1,2], labelV; title = "Bars");
        figsave(fig, fPath; figNotes = ["Notes"]);
        @test isfile(fPath)
        @test isfile(notesPath)
    end
end

function histogram_test()
    @testset "Histogram" begin
        fPath, notesPath = fig_test_setup("histogram_test.pdf");
        p, _ = histogram_plot(log.(LinRange(1.0, 2.0, 100)); xlabel = "x")
        figsave(p, fPath);
        @test isfile(fPath);
    end
end

function contour_test()
    @testset "Contour" begin
        fPath, notesPath = fig_test_setup("contour_test.pdf");
        xV = LinRange(1.0, 4.0, 10);
        yV = LinRange(2.0, 3.0, 8);
        zM = sin.(xV) .+ sin.(yV)';
        p, _ = contour_plot(xV, yV, zM; xlabel = "x", ylabel = "y");
        figsave(p, fPath);
        @test isfile(fPath);
    end
end

function line_plot_test()
    @testset "Line plot" begin
        fPath, notesPath = fig_test_setup("line_plot_test.pdf");
        nx = 7;
        xV = LinRange(-2.0, 1.5, nx);
        yM = LinRange(0.3, 0.5, nx) .+ LinRange(1.0, 0.5, 4)';
        p, ax = line_plot(xV, yM; xlabel = "x",
            labelV = ["Lbl $j" for j = 1 : 4], legPos = :below);
        figsave(p, fPath);
        @test isfile(fPath);

        fPath2, _ = fig_test_setup("add_line_test.pdf");
        yV = xV .+ 0.5;
        p, ax = line_plot(xV, yV; ylabel = "y", color = get_colors(1, 2),
            ylims = (0, nothing));
        add_line!(ax, xV, yV .+ 1.0; color = get_colors(2, 2));
        # ylims!(ax, 0, nothing);
        figsave(p, fPath2);
        @test isfile(fPath2);
    end
end

function line_error_bands_test()
    @testset "Line error bands" begin
        nx = 7;
        ny = 4;
        xV = LinRange(-2.0, 1.5, nx);
        yM = LinRange(0.3, 0.5, nx) .+ LinRange(1.0, 0.5, ny)';
        p, ax = line_plot(xV, yM; xlabel = "x",
            labelV = ["Lbl $j" for j = 1 : ny], legPos = :below);
        for iy = 1 : ny
            add_error_band!(ax, xV, yM[:,iy], LinRange(0.05, 0.1, nx);
                color = (get_colors(iy, ny), 0.2));
        end
        fPath3, _ = fig_test_setup("line_error_bands.pdf");
        figsave(p, fPath3);
        @test isfile(fPath3);
    end
end

function scatter_plot_test()
    @testset "Line plot" begin
        fPath, notesPath = fig_test_setup("scatter_plot_test.pdf");
        nx = 7;
        xV = LinRange(-2.0, 1.5, nx);
        yM = LinRange(0.3, 0.5, nx) .+ LinRange(1.0, 0.5, 4)';
        p, _ = scatter_plot(xV, yM; xlabel = "x");
        figsave(p, fPath);
        @test isfile(fPath);

        fPath2, _ = fig_test_setup("add_scatter_test.pdf");
        yV = xV .+ 0.5;
        p, ax = scatter_plot(xV, yV; ylabel = "y");
        add_scatter!(ax, xV, yV .+ 1.0; color = get_colors(2, 2));
        figsave(p, fPath2);
        @test isfile(fPath2);
    end
end

@testset "All" begin
    subplot_layout_test();
    fig_save_test();
    grouped_bar_test();
    bar_graph_test();
    line_plot_test();
    line_error_bands_test();
    scatter_plot_test();
    contour_test();
    histogram_test();
    include("subplot_test.jl");
    include("regression_test.jl");
    include("xyz_test.jl");
end

# ------------------


text:
using StaticArrays: SVector
using LinearAlgebra: I
using Combinatorics: combinations

function Point(coords::AbstractVector{<:Real})
    N = length(coords)
    return Point{N}(SVector{N, Float64}(coords))
end

Point(coords::Vararg{<:Real}) = Point(collect(coords))

Point(b::Barycentric) = Point(barycentric_matrix(b.simplex) * b.coords)

Point(b::SimpleBarycentric) = Point(Barycentric(b))

function Simplex(points::AbstractVector{Point{N}}) where N
    K = length(points)
    return Simplex{N, K}(SVector{K, Point{N}}(points))
end

Simplex(points::Vararg{Point{N}}) where N = Simplex(collect(points))

function Simplex(c::Cell)
    @assert length(c.points) == c.K
    return Simplex(c.points)
end

Simplex(s::SimpleSimplex) = Simplex(s.points)

SimpleSimplex(points::AbstractVector{Point{N}}) where N = SimpleSimplex{N}(collect(points))

SimpleSimplex(s::Simplex) = SimpleSimplex(s.points)

SimpleSimplex(s::Vector{SimpleBarycentric{N}}) where N = SimpleSimplex(map(Point, s))

function Barycentric(s::Simplex{N, K}, coords::AbstractVector{<:Real}) where {N, K}
    @assert length(coords) == K
    return Barycentric{N, K}(s, SVector{K, Float64}(coords))
end

Barycentric(s::Simplex, coords::Vararg{<:Real}) = Barycentric(s, collect(coords))

"""
    Barycentric(m::Metric{N}, s::Simplex{N}, p::Point{N}) where N

Find the Barycentric for the projection of a point onto the affine subspace spanned by a
simplex.
"""
function Barycentric(m::Metric{N}, s::Simplex{N}, p::Point{N}) where N
    A, B = barycentric_projection_matrices(m, s)
    coords = A * p.coords + B
    return Barycentric(s, coords)
end

Barycentric(b::SimpleBarycentric) = Barycentric(Simplex(b.simplex), b.coords)

SimpleBarycentric(s::SimpleSimplex{N}, coords::AbstractVector{<:Real}) where N =
    SimpleBarycentric{N}(s, collect(coords))

SimpleBarycentric(b::Barycentric) = SimpleBarycentric(SimpleSimplex(b.simplex), b.coords)

export subsimplices
"""
    subsimplices(s::Simplex{N, K}, k::Int) where {N, K}

Find all subsimplices of `s` with dimension `k-1`. If `k == 0` or `k > K`, return an
empty vector.
"""
function subsimplices(s::Simplex{N, K}, k::Int) where {N, K}
    @assert k >= 0
    if (k == 0) || (k > K)
        return Simplex{N, k}[]
    else
        return [Simplex(ps...) for ps in combinations(s.points, k)]
    end
end

export barycentric_matrix
"""
    barycentric_matrix(s::Simplex)

Find the matrix `M` so that `Mx` is the point corresponding to the barycentric
coordinates `x`.
"""
barycentric_matrix(s::Simplex) = hcat([p.coords for p in s.points]...)

export barycentric_projection_matrices
"""
    barycentric_projection_matrices(m::Metric{N}, s::Simplex{N, K}) where {N, K}

For a vector `p`, extremize `norm_square(m, p - barycentric_matrix(s) * x)` over `x`
with the constraint `sum(x) == 1`. The result is the barycentric coordinates of the
projection of `p` onto the affine subspace spanned by `s` and it can be written in the form
`x = A * p + B`. Return `A`, `B`.
"""
function barycentric_projection_matrices(m::Metric{N}, s::Simplex{N, K}) where {N, K}
    b = barycentric_matrix(s)
    # solve sum(x) == 1 by saying x = c * y + d
    c = transpose(hcat(Matrix{Float64}(I, K-1, K-1), -ones(K-1)))
    d = zeros(K); d[end] = 1
    bc = b * c; bcT = transpose(bc)
    A = c * ((bcT * m.mat * bc) \ (bcT * m.mat))
    B = d - A * b * d
    return A, B
end

export distance_quadratic
"""
    distance_quadratic(m::Metric{N}, s::Simplex{N})

The squared distance from a point with barycentric coordinates `x` to a point `p`
can be written in the form `xᵀAx + Bᵀx` where `B` depends on `p` but `A` does not. Return A.
"""
distance_quadratic(m::Metric{N}, s::Simplex{N}) where N = -.5 * circumsphere_matrix(m, s)

export distance_linear
"""
    distance_linear(m::Metric{N}, s::Simplex{N}, p::Point{N})

The squared distance from a point with barycentric coordinates `x` to a point `p`
can be written in the form `xᵀAx + Bᵀx` where `B` depends on `p` but `A` does not. Return B.
"""
distance_linear(m::Metric{N}, s::Simplex{N}, p::Point{N}) where N = [norm_square(
    m, p.coords - q.coords) for q in s.points]

export distance_square
"""
    distance_square(m::Metric{N}, b::Barycentric{N}, p::Point{N}) where N

Find the squared distance between the points `Point(b)` and `p`. Note that this
value can be negative if the metric is not positive semi-definite.
"""
function distance_square(m::Metric{N}, b::Barycentric{N}, p::Point{N}) where N
    A = distance_quadratic(m, b.simplex)
    B = distance_linear(m, b.simplex, p)
    return transpose(b.coords) * A * b.coords + transpose(B) * b.coords
end

export circumsphere_matrix
"""
    circumsphere_matrix(m::Metric{N}, s::Simplex{N}) where N

If `x` are the barycentric coordinates of the circumcenter of `s` then there is a matrix M
such that `M * x = 2 * R² * ones(K)` where `R²` is the square circumradius. Return M.
"""
function circumsphere_matrix(m::Metric{N}, s::Simplex{N}) where N
    norm_squares = [[norm_square(m, p.coords - q.coords) for p in s.points] for q in s.points]
    return transpose(hcat(norm_squares...))
end

export circumsphere_barycentric
"""
    circumsphere_barycentric(m::Metric{N}, s::Simplex{N, K}) where {N, K}

Find the Barycentric of the circumcenter and the squared circumradius.
"""
function circumsphere_barycentric(m::Metric{N}, s::Simplex{N, K}) where {N, K}
    if K == 1
        return Barycentric(s, [1.0]), 0
    else
        coords = circumsphere_matrix(m, s) \ ones(K)
        sum_coords = sum(coords)
        return Barycentric(s, coords / sum_coords), .5/sum_coords
    end
end

export circumsphere
"""
    circumsphere(m::Metric{N}, s::Simplex{N}) where N

Find the circumcenter and the squared circumradius.
"""
function circumsphere(m::Metric{N}, s::Simplex{N}) where N
    b, R² = circumsphere_barycentric(m, s)
    return Point(b), R²
end

export circumcenter
"""
    circumcenter(m::Metric{N}) where N

Return a function taking a `Simplex{N}` to the Barycentric of its circumcenter.
"""
circumcenter(m::Metric{N}) where N = s::Simplex{N} -> circumsphere_barycentric(m, s)[1]

export centroid
"""
    centroid(s::Simplex{N, K}) where {N, K}

Compute the Barycentric of the centroid of a simplex.
"""
centroid(s::Simplex{N, K}) where {N, K} = Barycentric(s, SVector{K, Float64}(ones(K)/K))

"""
    barycentric_subspace(m::Metric{N}, s::Simplex{N, K}, f::Simplex{N, J},
        p::Point{N}) where {N, K, J}

Let `s` and `f` be simplices where `f` lies in the affine subspace spanned by `s`. For a
point `p`, the subset of points with the same projection onto `f` as that of `p` is a line.
This line can be represented in the form `Mx = v` where `x` are barycentric coordinates
with respect to `s`. Return M and v.
"""
function barycentric_subspace(m::Metric{N}, s::Simplex{N, K}, f::Simplex{N, J},
    p::Point{N}) where {N, K, J}
    @assert J <= K
    # require that the projection of x onto f is the same as for p
    A, _ = barycentric_projection_matrices(m, f)
    M = A * barycentric_matrix(s)
    v = A * p.coords
    # since the sum of barycentric coordinates is 1, one equation is redundant
    # and one variable is redundant. Replace the last equation with sum(x) == 1
    M[end, :] = ones(K)
    v[end] = 1
    return M, v
end

"""
    rotate_about_face(m::Metric{N}, s::Simplex{N, K}, p::Point{N}, i::Int) where {N, K}

Let `s` be a simplex, `q = s.points[i]`, and `f` the face of `s` opposite `q`. Rotate the
point `p` about `f` into the affine subspace spanned by `s` and on the halfspace bounded
by `f` that does not contain `q`.
"""
function rotate_about_face(m::Metric{N}, s::Simplex{N, K}, p::Point{N}, i::Int) where {N, K}
    f = Simplex(s.points[filter(j -> j != i, 1:K)])
    a, b = parametrize_subspace(barycentric_subspace(m, s, f, p)...)
    q = f.points[1]
    A, B = distance_quadratic(m, s), distance_linear(m, s, q)
    C = -norm_square(m, p.coords - q.coords)
    x1, x2 = intersect_linear_quadratic(A, B, C, a, b[:]) # b[:] reshapes it to a vector
    x = x1[i] > 0 ? x2 : x1
    return Barycentric(s, x)
end

"""
    pairwise_delaunay(m::Metric{N}, s1::Simplex{N, K}, s2::Simplex{N, K}) where {N, K}

Given two simplices `s1` and `s2` that share a face, rotate `s2` into the affine subspace
spanned by `s1`. Let `p` be the image of the the point of `s2` that is not in `s1` under
this rotation and let `c` and `R²` be the circumcenter and squared circumradius of `s1`.
Return the squared distance from `p` to `c` minus `R²`.
"""
function pairwise_delaunay(m::Metric{N}, s1::Simplex{N, K}, s2::Simplex{N, K}) where {N, K}
    face_points = Point{N}[]
    i = 0
    for j in 1:K
        if s1.points[j] in s2.points
            push!(face_points, s1.points[j])
        else
            i = j
        end
    end
    @assert length(face_points) == K-1
    @assert i > 0
    p = first_setdiff(s2.points, face_points)
    b = rotate_about_face(m, s1, p, i)
    return transpose(b.coords) * distance_quadratic(m, s1) * b.coords
end


text:
@everywhere using Pkg
@everywhere Pkg.activate(".")
@everywhere using CmdStan, Distributions, Distributed, MCMCChains

ProjDir = @__DIR__

@everywhere mutable struct Job
  name::AbstractString
  model::AbstractString
  data::Dict
  output_format::Symbol
  tmpdir::AbstractString
  nchains::Int
  num_samples::Int
  num_warmup::Int
  save_warmup::Bool
  delta::Float64
  sm::Union{Nothing, Stanmodel}
end
function Job(
    name::AbstractString,
    model::AbstractString,
    data::Dict;
    output_format=:particles,
    tmpdir=pwd()
  )

  nchains=4
  num_samples=1000
  num_warmup=1000
  save_warmup=false
  delta=0.85

  sm = Stanmodel(
    CmdStan.Sample(
      num_samples=num_samples,
      num_warmup=num_warmup,
      save_warmup=save_warmup,
      adapt=CmdStan.Adapt(delta=delta)
    ),
    name=name, model=model, nchains=nchains,
    tmpdir=tmpdir, output_format=:mcmcchains
  )

  Job(name, model, data, output_format,
    tmpdir, nchains, num_samples, num_warmup,
    save_warmup, delta, sm)
end

# Assume we have 2 models, m1 & m2
m1 = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

m2 = "
data {
  int<lower=1> N;
  int<lower=0,upper=1> y[N];
}
parameters {
  real<lower=0,upper=1> theta;
}
model {
  theta ~ beta(1,1);
  y ~ bernoulli(theta);
}
";

# Assume we have 4 sets of data
d = Vector{Dict}(undef, 4)
d[1] = Dict("N" => 10, "y" => [0, 1, 0, 0, 0, 0, 0, 0, 0, 1])
d[2] = Dict("N" => 10, "y" => [0, 1, 0, 0, 0, 1, 1, 1, 0, 1])
d[3] = Dict("N" => 10, "y" => [0, 1, 0, 1, 0, 0, 0, 0, 0, 1])
d[4] = Dict("N" => 10, "y" => [0, 1, 1, 1, 0, 1, 1, 1, 0, 1])

tmpdir = "$(ProjDir)/tmpcmd"
@everywhere jobs = Vector{Job}(undef, 4)
jobs[1] = Job("m1.1", m1, d[1]; tmpdir=tmpdir*"1")
jobs[2] = Job("m1.2", m1, d[2]; tmpdir=tmpdir*"2")
jobs[3] = Job("m2.1", m2, d[3]; tmpdir=tmpdir*"3")
jobs[4] = Job("m2.2", m2, d[4]; tmpdir=tmpdir*"4")

@everywhere function runjob(i, jobs)
  p = []
  rc, chn, _ = stan(jobs[i].sm, jobs[i].data)
  if rc == 0
    push!(p, chn)
    println("Job $i completed")
  else
    println("Job[i] failed, adjust p indeces accordingly.")
  end
  p
end

@time res = pmap(i -> runjob(i, jobs), 1:length(jobs));
for i in 1:length(jobs)
  display(res[i])
end


text:
export PauliID, pauli_I, pauli_X, pauli_Y, pauli_Z
export QuantumCircuit, QulacsQuantumCircuit
export ParametricQuantumCircuit, QulacsParametricQuantumCircuit
export add_parametric_multi_Pauli_rotation_gate!, set_parameter!
export add_CNOT_gate!, add_H_gate!, add_X_gate!, add_S_gate!,add_Sdag_gate!, add_Z_gate!
export get_parameter_count, get_parameter
export add_parametric_RX_gate!
export add_parametric_RY_gate!
export add_parametric_RZ_gate!
export add_RY_gate!
export add_RZ_gate!
export add_SWAP_gate!
export add_gate!
export get_thetas

export VariationalQuantumCircuit
export QulacsVariationalQuantumCircuit

export Gate, QulacsGate, RX, RY, RZ, get_matrix, add_U1_gate!

export QuantumState, QulacsQuantumState
export set_computational_basis!, create_hf_state
export FermionOperator, jordan_wigner, get_number_preserving_sparse_operator
export get_vector, get_n_qubit

export QubitOperator, OFQubitOperator
export get_term_count, get_n_qubit, terms_dict, is_hermitian
export get_expectation_value, get_transition_amplitude
export inner_product

export get_expectation_value, create_operator_from_openfermion
export hermitian_conjugated
export create_observable

export imag_time_evolve
export compute_gtau
export compute_A
export compute_C
export compute_thetadot
#export _create_quantum_state
export taus_list
export read_taus_list
export compute_next_thetas_direct
export compute_next_thetas_vqs
export compute_fubini
export compute_F
export overlap, compute_B, compute_B2, compute_B2x, compute_B3 #vqs.jl

max_stars_repo_name: Sakurairihito/QCMaterialNew

max_issues_repo_name: sakurairihito/QCMaterialNew

max_forks_repo_name: sakurairihito/QCMaterialNew


text:
using DelimitedFiles

filename = "./log-conv-euler-1D-nodewise_CFL.5"
file = readdlm(filename)

num_scheme = 1#2
Narr = [2;5]
Karr = [50;100;200;400;800]
block_size = 4
num_lines = size(file,1)
num_block = div(num_lines,block_size)
L1_err_dict = Dict()
for b in 1:num_block
    l = 1+block_size*(b-1)
    N = parse(Int,file[l,3][1])
    K = file[l,6]
    L1_err = file[l+1,4]
    @show N,K
    @show L1_err
    if !haskey(L1_err_dict,(N,K))
        L1_err_dict[(N,K)] = [L1_err]
    else
        push!(L1_err_dict[(N,K)],L1_err)
    end
end

num_row = length(Karr)-1
num_col = num_scheme*length(Narr)
conv_table = zeros(num_row,num_col)
for i = 1:num_row
    for j = 1:num_col
        K = Karr[i]
        N = Narr[div(j-1,num_scheme)+1]
        s = mod1(j,num_scheme)
        conv_table[i,j] = log2(L1_err_dict[(N,K)][s]/L1_err_dict[(N,2*K)][s])
    end
end

max_stars_repo_name: yiminllin/ESDG-PosLimit

max_issues_repo_name: yiminllin/ESDG-PosLimit

max_forks_repo_name: yiminllin/ESDG-PosLimit


text:
#=
analyze benchmark results
~/julia/julia analyze.jl
=#

using Printf
using CSV
using DataFrames

include(joinpath(@__DIR__, "../examples/Examples.jl"))
import Main.Examples

bench_file = joinpath(@__DIR__, "raw", "bench.csv")

output_dir = mkpath(joinpath(@__DIR__, "analysis"))
tex_dir = mkpath(joinpath(output_dir, "tex"))
stats_dir = mkpath(joinpath(output_dir, "stats"))
csv_dir = mkpath(joinpath(output_dir, "csvs"))

# uncomment examples to process
examples_params = [
    # PSD:
    "completablepsd" => ([:d], [1], [:nat, :ext], [:iters]),
    # WSOS:
    "polyfacilitylocation" => ([:n], [1], [:nat, :ext], [:iters]),
    "polyregression" => ([:m], [3], [:nat, :ext], [:iters]),
    "twostagestochastic" => ([:halfdeg], [2], [:nat, :ext], [:iters]),
    # spectral norm:
    "matrixcompletion" => ([:nrow], [3], [:nat, :ext], [:iters]),
    "matrixdecomposition" => ([:n], [2], [:nat, :ext], [:iters]),
    "matrixregression" => ([:p], [2], [:nat, :ext], [:iters]),
    # spectral function:
    "inversecovariance" => ([:d], [1], [:nat, :ext], [:iters]),
    "vectorregression" => ([:n], [1], [:nat, :noext], [:iters]),
    "knapsack" => ([:n], [1], [:cont_geo, :cont_noext_geo], [:iters], "cont_geo"),
    "knapsack" => ([:n], [1], [:cont_log, :cont_noext_log], [:iters], "cont_log"),
    "knapsack" => ([:n], [1], [:cont_inv, :cont_noext_inv], [:iters], "cont_inv"),
    "knapsack" => ([:n], [1], [:nat_geo, :noext_geo, :ext_geo], [:iters], "geo"),
    "experimentdesign" => ([:d], [1], [:nat_rtdet, :ext_rtdet], [:iters], "rtdet"),
    "experimentdesign" => ([:d], [1], [:nat_entr, :ext_entr], [:iters], "entr"),
    # MIP formulations:
    "ballpacking" => ([:num_pts], [4], [:sos2, :logib, :cc], [:nodes]),
    "modulardesign" => ([:n], [3], [:nat, :noext], [:iters], "convex"),
    "modulardesign" => ([:n], [3], [:sos2, :logib, :cc], [:iters], "nonconvex"),
]

function analyze()
    col_types = Dict(:example => String, :inst_set => String, :status => String)
    all_df = CSV.read(bench_file, DataFrame, types = col_types)
    @show all_df

    for (ex_name, ex_params) in examples_params
        println()
        @info("starting $ex_name with params: $ex_params")

        sets = string.(ex_params[3])
        ex_df = filter(
            t -> (t.inst_num > 1) && (t.example == ex_name) && (t.inst_set in sets),
            all_df,
        )
        if isempty(ex_df)
            @info("no data for $ex_name")
            continue
        end

        if length(ex_params) == 5
            ex_name *= "_" * ex_params[5]
        end

        ex_df_wide = make_wide_csv(ex_df, ex_name, ex_params)
        make_table_tex(ex_name, ex_params, ex_df_wide)
        make_plot_csv(ex_name, ex_params, ex_df_wide)

        @info("finished")
    end

    println()
    @info("finished all")
end

status_map = Dict(
    "OPTIMAL" => "co",
    "TIME_LIMIT" => "tl",
    "ITERATION_LIMIT" => "il",
    "OTHER_ERROR" => "er",
)

function relative_tol_satisfied(a::T, b::T, tol::T = 1e-3) where {T <: Real}
    return (abs(a - b) / (1 + max(abs(a), abs(b))) < tol)
end

function make_wide_csv(ex_df, ex_name, ex_params)
    @info("making wide csv for $ex_name")
    inst_keys = ex_params[1]
    str(x) = split(x[2:(end - 1)], ", ")

    replace!(ex_df.status, status_map...)

    for (name, pos) in zip(inst_keys, ex_params[2])
        transform!(
            ex_df,
            :inst_data => ByRow(x -> string(eval(Meta.parse(x))[pos])) => name,
        )
    end

    # check objectives if solver claims optimality
    show_cols = [:inst_set, :inst_num, :status, :primal_obj]
    for group_df in groupby(ex_df, inst_keys)
        # check all pairs of verified converged results
        co_idxs = findall(group_df[:, :status] .== "co")
        (length(co_idxs) < 2) && continue

        first_optval = group_df[co_idxs[1], :primal_obj]
        other_optvals = group_df[co_idxs[2:end], :primal_obj]
        if !all(relative_tol_satisfied.(other_optvals, first_optval))
            dat = group_df[!, :inst_data][1]
            println("$ex_name $dat primal optimal objective values disagree:")
            show(group_df[!, show_cols], eltypes = false)
            println("\n")
        end
    end

    unstacked = [
        unstack(
            ex_df,
            inst_keys,
            :inst_set,
            v,
            renamecols = x -> Symbol(v, :_, x),
            allowduplicates = true,
        ) for v in [:status, :solve_time, :nodes, :iters, :num_subp, :num_cuts]
    ]
    ex_df_wide = outerjoin(unstacked..., on = inst_keys)

    CSV.write(joinpath(stats_dir, ex_name * "_wide.csv"), ex_df_wide)

    return ex_df_wide
end

process_entry(::Missing) = "\$\\ast\$"

function process_entry(x::Float64)
    isnan(x) && return "\$\\ast\$"
    @assert x > 0
    if x < 10
        return @sprintf("%.1f", x)
    else
        return @sprintf("%.0f.", x)
    end
end

process_entry(x) = string(x)

function make_table_tex(ex_name, ex_params, ex_df_wide)
    @info("making table tex for $ex_name")
    inst_sets = string.(ex_params[3])
    inst_keys = ex_params[1]
    col_names = ex_params[4]

    sep = " & "
    ex_tex = open(joinpath(tex_dir, ex_name * "_table.tex"), "w")

    header_str = "% " * prod(string(s) * sep for s in inst_keys)
    for s in inst_sets
        header_str *= s * prod(sep * string(c) for c in col_names) * sep * "time" * sep
    end
    println(ex_tex, header_str)

    for row in eachrow(ex_df_wide)
        row_str = process_entry(row[1])
        for i in 2:length(inst_keys)
            row_str *= sep * process_entry(row[i])
        end
        for inst_set in inst_sets
            row_str *= sep * process_entry(row[Symbol(:status_, inst_set)])
            for col_name in col_names
                row_str *= sep * process_entry(row[Symbol(col_name, :_, inst_set)])
            end
            row_str *= sep * process_entry(row[Symbol(:solve_time_, inst_set)])
        end
        row_str *= " \\\\"
        println(ex_tex, row_str)
    end

    close(ex_tex)

    return
end

function transform_plot_cols(ex_df_wide, inst_set)
    old_cols = Symbol.([:status_, :solve_time_], inst_set)
    # old_cols = Symbol.([:status_, :iters_], inst_set)
    return transform!(
        ex_df_wide,
        old_cols =>
            ByRow((x, y) -> ((!ismissing(x) && x == "co") ? y : missing)) => inst_set,
    )
end

function make_plot_csv(ex_name, ex_params, ex_df_wide)
    @info("making plot csv for $ex_name")
    inst_sets = string.(ex_params[3])
    inst_keys = ex_params[1]

    for inst_set in inst_sets
        transform_plot_cols(ex_df_wide, inst_set)
    end

    plot_file_start = joinpath(csv_dir, ex_name * "_plot")

    axis_name = last(inst_keys)
    if length(inst_keys) == 1
        success_df = select(ex_df_wide, axis_name, inst_sets...)
        CSV.write(plot_file_start * ".csv", success_df)
    else
        group_name = first(inst_keys)
        success_df = select(ex_df_wide, axis_name, group_name, inst_sets...)
        for (group_id, group_df) in pairs(groupby(success_df, group_name))
            CSV.write(
                plot_file_start * "_$(group_id[1]).csv",
                select(group_df, Not(group_name)),
            )
        end
    end

    return
end

analyze();


text:
#===============================================================================
# PID controller written in pure Julia.
#
# Author: Carl Cortright
# Date: 1/1/2017
#
===============================================================================#

module PIDLib

export PID, Update

type PID
  P::Float64
  I::Float64
  D::Float64
  setpoint::Float64
  errors::AbstractArray{Float64, 1}
  epochs::Int64

  PID(P, I, D, setpoint) = new(P, I, D, setpoint, zeros(Float64, 10000), 1)
end

#
# Updates the PID mechanism and
#
function Update(pid::PID, value::Float64)::Float64
  error = value - pid.setpoint
  pid.errors[pid.epochs + 1] = error
  pid.epochs += 1
  return (pid.P * error) + (pid.I * RiemannSum(pid)) + (pid.D * Derivative(pid))
end

#
# Takes the "integral" over the errors from 0 to t
#
function RiemannSum(pid::PID)
  sum = 0.0
  for i in 1:length(pid.errors)
    sum += pid.errors[i]
  end
  return sum
end

#
# Takes the current Derivative of the errors
#
function Derivative(pid::PID)
  return pid.errors[end] - pid.errors[end - 1]
end

end  # module pid

max_issues_repo_name: carlcortright/Playground


text:
using FractionalDiffEq, Plots
rightfun(x)=sin(x^2)
prob = MultiTermsFODEProblem([1 8 26 73 90], [3.5 3.1 2.3 1.2 0.5], rightfun, [30 90], [1 0.3], [0, 0], (0, 10))
sol = solve(prob, 0.02, ClosedFormHankelM())
plot(sol)