text:
using Test
using PolynomialBasis
using ImplicitDomainQuadrature
# using Revise
using CutCellDG
include("../useful_routines.jl")

polyorder = 1
numqp = 2
elasticitybasis = LagrangeTensorProductBasis(2, polyorder)
levelsetbasis = HermiteTensorProductBasis(2)
levelset = InterpolatingPolynomial(1, levelsetbasis)

x0 = [0.0, 0.0]
meshwidths = [3.0, 1.0]
nelements = [3, 1]
dim,numinterp = size(interpolation_points(levelsetbasis))
points = interpolation_points(elasticitybasis)
dgmesh = CutCellDG.DGMesh(x0, meshwidths, nelements, points)
cgmesh = CutCellDG.CGMesh(x0, meshwidths, nelements, numinterp)

normal = [1.0, 0.0]
xI = [0.5, 0.0]
tol = 1e-3
perturbation = 1e-3

levelset = CutCellDG.LevelSet(
    x -> plane_distance_function(x, normal, xI),
    cgmesh,
    levelsetbasis,
)

cutmesh = CutCellDG.CutMesh(dgmesh, levelset)

cellquads = CutCellDG.CellQuadratures(cutmesh, levelset, numqp)
@test length(cellquads.quads) == 3
testcelltoquad = [
    2 1 1
    3 0 0
]
@test allequal(cellquads.celltoquad, testcelltoquad)


text:
# # Adaptive

# Here I show how one can use [Adaptive.jl](https://github.com/akels/Adaptive.jl) wrapper for [python adaptive package](https://github.com/python-adaptive/adaptive) to make adaptively sampled figures. Before we start, let's set up our environment by loading packages:
using Distributed
addprocs(2)
using TaskMaster
using Adaptive
using PyPlot
# where Adaptive.jl at the moment needs to be added directly from the GitHub repository.

# ## AdaptiveLearner1D

@everywhere f(x) = exp(-x^2)

fig = figure()

x = collect(range(-2,stop=2,length=200))
plot(x,f.(x),label=L"e^{-x^2}")

xx = collect(range(-2,stop=2,length=20))
plot(xx,f.(xx),".-",label="even sampling")

master = WorkMaster(f)
learner1d = AdaptiveLearner1D((-2,+2))
loop = Loop(master,learner1d)
evaluate!(loop,1:20)

plot(learner1d.x,learner1d.y,".-",label="AdaptiveLearner1D")

legend()
savefig("learner1d.svg")

# ![](learner1d.svg)

# ## AdaptiveLearner2D

@everywhere f(p) = exp(-p[1]^2 - p[2]^2)

master = WorkMaster(f)
learner2d = AdaptiveLearner2D([(-3,+3),(-3,+3)])
loop = Loop(master,learner2d)
evaluate!(loop,learner->learner.loss()<0.05)

fig = figure()

p,tri,v = learner2d.points, learner2d.vertices, learner2d.values

tricontourf(p[:,1],p[:,2],tri.-1,v)
triplot(p[:,1],p[:,2],tri.-1,"k.-")

savefig("learner2d.svg")

# ![](learner2d.svg)

ext: jl


text:
export uncertainty_exponent, basins_fractal_dimension

"""
    basins_fractal_dimension(basins; kwargs...) -> V_ε, N_ε ,d
Estimate the [Fractal Dimension](@ref) `d` of the boundary between basins of attraction using
the box-counting algorithm.

The output `N_ε` is a vector with the number of the balls of radius `ε` (in pixels)
that contain at least two initial conditions that lead to different attractors. `V_ε`
is a vector with the corresponding size of the balls. The ouput `d` is the estimation
of the box-counting dimension of the boundary by fitting a line in the `log.(N_ε)`
vs `log.(1/V_ε)` curve. However it is recommended to analyze the curve directly
for more accuracy.

## Keyword arguments
* `range_ε = 2:maximum(size(basins))÷20` is the range of sizes of the ball to
  test (in pixels).

## Description

It is the implementation of the popular algorithm of the estimation of the box-counting
dimension. The algorithm search for a covering the boundary with `N_ε` boxes of size
`ε` in pixels.
"""
function basins_fractal_dimension(basins::AbstractArray; range_ε = 3:maximum(size(basins))÷20)

    dims = size(basins)
    num_step = length(range_ε)
    N_u = zeros(Int, num_step) # number of uncertain box
    N = zeros(Int, num_step) # number of boxes
    V_ε = zeros(1, num_step) # resolution

    # Naive box counting estimator
    for (k,eps) in enumerate(range_ε)
        Nb, Nu = 0, 0
        # get indices of boxes
        bx_tuple = ntuple(i -> range(1, dims[i] - rem(dims[i],eps), step = eps), length(dims))
        box_indices = CartesianIndices(bx_tuple)
        for box in box_indices
            # compute the range of indices for the current box
            ind = CartesianIndices(ntuple(i -> range(box[i], box[i]+eps-1, step = 1), length(dims)))
            c = basins[ind]
            if length(unique(c))>1
                Nu = Nu + 1
            end
            Nb += 1
        end
        N_u[k] = Nu
        N[k] = Nb
        V_ε[k] = eps
    end
    N_ε = N_u
    # remove zeros in case there are any:
    ind = N_ε .> 0.0
    N_ε = N_ε[ind]
    V_ε = V_ε[ind]
    # get exponent via liner regression on `f_ε ~ ε^α`
    b, d = linreg(vec(-log10.(V_ε)), vec(log10.(N_ε)))
    return V_ε, N_ε, d
end

"""
    uncertainty_exponent(basins; kwargs...) -> ε, N_ε ,α
Estimate the uncertainty exponent[^Grebogi1983] of the basins of attraction. This exponent
is related to the final state sensitivity of the trajectories in the phase space.
An exponent close to `1` means basins with smooth boundaries whereas an exponent close
to `0` represent complety fractalized basins, also called riddled basins.

The output `N_ε` is a vector with the number of the balls of radius `ε` (in pixels)
that contain at least two initial conditions that lead to different attractors.
The ouput `α` is the estimation of the uncertainty exponent using the box-counting
dimension of the boundary by fitting a line in the `log.(N_ε)` vs `log.(1/ε)` curve.
However it is recommended to analyze the curve directly for more accuracy.

## Keyword arguments
* `range_ε = 2:maximum(size(basins))÷20` is the range of sizes of the ball to
  test (in pixels).

## Description

A phase space with a fractal boundary may cause a uncertainty on the final state of the
dynamical system for a given initial condition. A measure of this final state sensitivity
is the uncertainty exponent. The algorithm probes the basin of attraction with balls
of size `ε` at random. If there are a least two initial conditions that lead to different
attractors, a ball is tagged "uncertain". `f_ε` is the fraction of "uncertain balls" to the
total number of tries in the basin. In analogy to the fractal dimension, there is a scaling
law between, `f_ε ~ ε^α`. The number that characterizes this scaling is called the
uncertainty exponent `α`.

Notice that the uncertainty exponent and the box counting dimension of the boundary are
related. We have `Δ₀ = D - α` where `Δ₀` is the box counting dimension computed with
[`basins_fractal_dimension`](@ref) and `D` is the dimension of the phase space.
The algorithm first estimates the box counting dimension of the boundary and
returns the uncertainty exponent.

[^Grebogi1983]: C. Grebogi, S. W. McDonald, E. Ott and J. A. Yorke, Final state sensitivity: An obstruction to predictability, Physics Letters A, 99, 9, 1983
"""
function uncertainty_exponent(basins::AbstractArray; range_ε = 2:maximum(size(basins))÷20)
    V_ε, N_ε, d = basins_fractal_dimension(basins; range_ε)
    return V_ε, N_ε, length(size(basins)) - d
end


text:
# Matern 5/2 isotropic covariance function

"""
    Mat52Iso <: MaternIso

Isotropic Matern 5/2 kernel (covariance)
```math
k(x,x') = σ²(1+√5|x-x'|/ℓ + 5|x-x'|²/(3ℓ²))\\exp(- √5|x-x'|/ℓ)
```
with length scale ``ℓ`` and signal standard deviation ``σ``.
"""
mutable struct Mat52Iso <: MaternIso
    "Length scale"
    ℓ::Float64
    "Signal variance"
    σ2::Float64
    "Priors for kernel parameters"
    priors::Array

    """
        Mat52Iso(ll::Float64, lσ::Float64)

    Create `Mat52Iso` with length scale `exp(ll)` and signal standard deviation `exp(lσ)`.
    """
    Mat52Iso(ll::Float64, lσ::Float64) = new(exp(ll), exp(2 * lσ), [])
end

function set_params!(mat::Mat52Iso, hyp::VecF64)
    length(hyp) == 2 || throw(ArgumentError("Matern 5/2 only has two parameters"))
    mat.ℓ, mat.σ2 = exp(hyp[1]), exp(2 * hyp[2])
end
get_params(mat::Mat52Iso) = Float64[log(mat.ℓ), log(mat.σ2) / 2]
get_param_names(mat::Mat52Iso) = [:ll, :lσ]
num_params(mat::Mat52Iso) = 2

Statistics.cov(mat::Mat52Iso, r::Float64) =
    (s = √5 * r / mat.ℓ; mat.σ2 * (1 + s + s^2 / 3) * exp(-s))

@inline dk_dll(mat::Mat52Iso, r::Float64) =
    (s = √5 * r / mat.ℓ; mat.σ2 / 3 * s^2 * (1 + s) * exp(-s))


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

module ChabocheModule

using LinearAlgebra, ForwardDiff, Tensors, NLsolve, Parameters

import ..AbstractMaterial, ..AbstractMaterialState
import ..Utilities: Symm2, Symm4, isotropic_elasticity_tensor, lame, debang
import ..integrate_material!  # for method extension

# parametrically polymorphic for any type representing ℝ
export GenericChaboche, GenericChabocheDriverState, GenericChabocheParameterState, GenericChabocheVariableState

# specialization for Float64
export Chaboche, ChabocheDriverState, ChabocheParameterState, ChabocheVariableState

@with_kw mutable struct GenericChabocheDriverState{T <: Real} <: AbstractMaterialState
    time::T = zero(T)
    strain::Symm2{T} = zero(Symm2{T})
end

"""Parameter state for Chaboche material.

The classical viscoplastic material is a special case of this model with `C1 = C2 = 0`.

- `E`: Young's modulus
- `nu`: Poisson's ratio
- `R0`: initial yield strength
- `Kn`: plasticity multiplier divisor (drag stress)
- `nn`: plasticity multiplier exponent
- `C1`, `D1`: parameters governing behavior of backstress X1
- `C2`, `D2`: parameters governing behavior of backstress X2
- `Q`: hardening saturation state
- `b`: rate of convergence to hardening saturation
"""
@with_kw struct GenericChabocheParameterState{T <: Real} <: AbstractMaterialState
    E::T = 0
    nu::T = 0
    R0::T = 0
    Kn::T = 0
    nn::T = 0
    C1::T = 0
    D1::T = 0
    C2::T = 0
    D2::T = 0
    Q::T = 0
    b::T = 0
end

"""Problem state for Chaboche material.

- `stress`: stress tensor
- `X1`: backstress 1
- `X2`: backstress 2
- `plastic_strain`: plastic part of strain tensor
- `cumeq`: cumulative equivalent plastic strain (scalar, ≥ 0)
- `R`: yield strength
- `jacobian`: ∂σij/∂εkl
"""
@with_kw struct GenericChabocheVariableState{T <: Real} <: AbstractMaterialState
    stress::Symm2{T} = zero(Symm2{T})
    X1::Symm2{T} = zero(Symm2{T})
    X2::Symm2{T} = zero(Symm2{T})
    plastic_strain::Symm2{T} = zero(Symm2{T})
    cumeq::T = zero(T)
    R::T = zero(T)
    jacobian::Symm4{T} = zero(Symm4{T})
end

# TODO: Does this eventually need a {T}?
@with_kw struct ChabocheOptions <: AbstractMaterialState
    nlsolve_method::Symbol = :trust_region
end

@with_kw mutable struct GenericChaboche{T <: Real} <: AbstractMaterial
    drivers::GenericChabocheDriverState{T} = GenericChabocheDriverState{T}()
    ddrivers::GenericChabocheDriverState{T} = GenericChabocheDriverState{T}()
    variables::GenericChabocheVariableState{T} = GenericChabocheVariableState{T}()
    variables_new::GenericChabocheVariableState{T} = GenericChabocheVariableState{T}()
    parameters::GenericChabocheParameterState{T} = GenericChabocheParameterState{T}()
    dparameters::GenericChabocheParameterState{T} = GenericChabocheParameterState{T}()
    options::ChabocheOptions = ChabocheOptions()
end

ChabocheDriverState = GenericChabocheDriverState{Float64}
ChabocheParameterState = GenericChabocheParameterState{Float64}
ChabocheVariableState = GenericChabocheVariableState{Float64}
Chaboche = GenericChaboche{Float64}

"""
    state_to_vector(sigma::U, R::T, X1::U, X2::U) where U <: Symm2{T} where T <: Real

Adaptor for `nlsolve`. Marshal the problem state into a `Vector`.
"""
function state_to_vector(sigma::U, R::T, X1::U, X2::U) where U <: Symm2{T} where T <: Real
    return [tovoigt(sigma); R; tovoigt(X1); tovoigt(X2)]::Vector{T}
end

"""
    state_from_vector(x::AbstractVector{<:Real})

Adaptor for `nlsolve`. Unmarshal the problem state from a `Vector`.
"""
function state_from_vector(x::AbstractVector{T}) where T <: Real
    sigma::Symm2{T} = fromvoigt(Symm2{T}, @view x[1:6])
    R::T = x[7]
    X1::Symm2{T} = fromvoigt(Symm2{T}, @view x[8:13])
    X2::Symm2{T} = fromvoigt(Symm2{T}, @view x[14:19])
    return sigma, R, X1, X2
end

"""
    integrate_material!(material::GenericChaboche{T}) where T <: Real

Chaboche material with two backstresses. Both kinematic and isotropic hardening.

See:

    J.-L. Chaboche. Constitutive equations for cyclic plasticity and cyclic
    viscoplasticity. International Journal of Plasticity 5(3) (1989), 247--302.
    https://doi.org/10.1016/0749-6419(89)90015-6

Further reading:

    J.-L. Chaboche. A review of some plasticity and viscoplasticity constitutive
    theories. International Journal of Plasticity 24 (2008), 1642--1693.
    https://dx.doi.org/10.1016/j.ijplas.2008.03.009

    J.-L. Chaboche, A. Gaubert, P. Kanouté, A. Longuet, F. Azzouz, M. Mazière.
    Viscoplastic constitutive equations of combustion chamber materials including
    cyclic hardening and dynamic strain aging. International Journal of Plasticity
    46 (2013), 1--22. https://dx.doi.org/10.1016/j.ijplas.2012.09.011
"""
function integrate_material!(material::GenericChaboche{T}) where T <: Real
    p = material.parameters
    v = material.variables
    dd = material.ddrivers
    d = material.drivers
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b = p
    lambda, mu = lame(E, nu)

    @unpack strain, time = d
    dstrain = dd.strain
    dtime = dd.time
    @unpack stress, X1, X2, plastic_strain, cumeq, R = v

    # elastic part
    jacobian = isotropic_elasticity_tensor(lambda, mu)  # dσ/dε, i.e. ∂σij/∂εkl
    stress += dcontract(jacobian, dstrain)  # add the elastic stress increment, get the elastic trial stress

    # resulting deviatoric plastic stress (accounting for backstresses Xm)
    seff_dev = dev(stress - X1 - X2)
    # von Mises yield function
    f = sqrt(1.5)*norm(seff_dev) - (R0 + R)  # using elastic trial problem state
    if f > 0.0
        g! = create_nonlinear_system_of_equations(material)
        x0 = state_to_vector(stress, R, X1, X2)
        res = nlsolve(g!, x0; method=material.options.nlsolve_method, autodiff=:forward)  # user manual: https://github.com/JuliaNLSolvers/NLsolve.jl
        converged(res) || error("Nonlinear system of equations did not converge!")
        x = res.zero
        stress, R, X1, X2 = state_from_vector(x)

        # using the new problem state
        seff_dev = dev(stress - X1 - X2)
        f = sqrt(1.5)*norm(seff_dev) - (R0 + R)

        dotp = ((f >= 0.0 ? f : 0.0)/Kn)^nn  # power law viscoplasticity (Norton-Bailey type)
        dp = dotp*dtime  # |dε_p|, using backward Euler (dotp is ∂ε_p/∂t at the end of the timestep)
        n = sqrt(1.5)*seff_dev/norm(seff_dev)  # Chaboche: a (tensorial) unit direction, s.t. 2/3 * (n : n) = 1; also n = ∂f/∂σ.

        plastic_strain += dp*n
        cumeq += dp   # cumulative equivalent plastic strain (note dp ≥ 0)

        # Compute the new Jacobian, accounting for the plastic contribution. Because
        #   x ≡ [σ R X1 X2]   (vector of length 19, with tensors encoded in Voigt format)
        # we have
        #   dσ/dε = (dx/dε)[1:6,1:6]
        # for which we can compute the LHS as follows:
        #   dx/dε = dx/dr dr/dε = inv(dr/dx) dr/dε ≡ (dr/dx) \ (dr/dε)
        # where r = r(x) is the residual, given by the function g!. AD can get us dr/dx automatically,
        # the other factor we will have to supply manually.
        drdx = ForwardDiff.jacobian(debang(g!), x)  # Array{19, 19}
        drde = zeros((length(x),6))                 # Array{19, 6}
        drde[1:6, 1:6] = tovoigt(jacobian)  # elastic Jacobian. Follows from the defn. of g!.
        jacobian = fromvoigt(Symm4, (drdx\drde)[1:6, 1:6])
    end
    variables_new = GenericChabocheVariableState{T}(stress = stress,
                                                    X1 = X1,
                                                    X2 = X2,
                                                    R = R,
                                                    plastic_strain = plastic_strain,
                                                    cumeq = cumeq,
                                                    jacobian = jacobian)
    material.variables_new = variables_new
    return nothing
end

"""
    create_nonlinear_system_of_equations(material::GenericChaboche{T}) where T <: Real

Create and return an instance of the equation system for the incremental form of
the evolution equations of the Chaboche material.

Used internally for computing the plastic contribution in `integrate_material!`.

The input `material` represents the problem state at the end of the previous
timestep. The created equation system will hold its own copy of that state.

The equation system is represented as a mutating function `g!` that computes the
residual:

```julia
    g!(F::V, x::V) where V <: AbstractVector{<:Real}
```

Both `F` (output) and `x` (input) are length-19 vectors containing
[sigma, R, X1, X2], in that order. The tensor quantities sigma, X1,
X2 are encoded in Voigt format.

The function `g!` is intended to be handed over to `nlsolve`.
"""
function create_nonlinear_system_of_equations(material::GenericChaboche{T}) where T <: Real
    p = material.parameters
    v = material.variables
    dd = material.ddrivers
    d = material.drivers
    @unpack E, nu, R0, Kn, nn, C1, D1, C2, D2, Q, b = p
    lambda, mu = lame(E, nu)

    # Old problem state (i.e. the problem state at the time when this equation
    # system instance was created).
    #
    # Note this does not include the elastic trial; this is the state at the
    # end of the previous timestep.
    @unpack strain, time = d
    dstrain = dd.strain
    dtime = dd.time
    @unpack stress, X1, X2, plastic_strain, cumeq, R = v

    jacobian = isotropic_elasticity_tensor(lambda, mu)

    # Compute the residual. F is output, x is filled by NLsolve.
    # The solution is x = x* such that g(x*) = 0.
    function g!(F::V, x::V) where V <: AbstractVector{<:Real}
        stress_new, R_new, X1_new, X2_new = state_from_vector(x)  # tentative new values from nlsolve

        seff_dev = dev(stress_new - X1_new - X2_new)
        f = sqrt(1.5)*norm(seff_dev) - (R0 + R_new)

        dotp = ((f >= 0.0 ? f : 0.0)/Kn)^nn
        dp = dotp*dtime
        n = sqrt(1.5)*seff_dev/norm(seff_dev)

        # The equations are written in an incremental form:
        #
        # Δσ = (∂σ/∂ε)_e : dε_e = (∂σ/∂ε)_e : (dε - dε_p)   (components 1:6)
        # ΔR = b (Q - R_new) |dε_p|                         (component 7)
        # ΔX1 = (2/3) C1 |dε_p| (n - (3/2) (D1/C1) X1_new)  (components 8:13)
        # ΔX2 = (2/3) C2 |dε_p| (n - (3/2) (D2/C2) X2_new)  (components 14:19)
        #
        # where
        #
        # Δ(...) = (...)_new - (...)_old
        #
        # Then move the terms on the RHS to the LHS to get the standard form, (stuff) = 0.
        # Also, below we avoid the multiplication and division that cancel each other
        # in the last terms of the equations for ΔX1 and ΔX2.
        #
        dstrain_plastic = dp*n
        dstrain_elastic = dstrain - dstrain_plastic
        tovoigt!(view(F, 1:6), stress_new - stress - dcontract(jacobian, dstrain_elastic))

        F[7] = R_new - R - b*(Q - R_new)*dp

        tovoigt!(view(F,  8:13), X1_new - X1 - dp*(2.0/3.0*C1*n - D1*X1_new))
        tovoigt!(view(F, 14:19), X2_new - X2 - dp*(2.0/3.0*C2*n - D2*X2_new))
        return nothing
    end
    return g!
end

end


text:
using MenuAdventures
using Documenter: deploydocs, makedocs

makedocs(
    sitename = "MenuAdventures.jl",
    modules = [MenuAdventures],
    doctest = false,
    pages = [
        "Public interface" => "index.md",
        "Internals" => "internals.md",
        "ExtraDirections extension" => "ExtraDirections.md",
        "ExtraActions extension" => "ExtraActions.md",
        "Onto extension" => "Onto.md",
        "Outfits extension" => "Outfits.md",
        "Parts extension" => "Parts.md",
        "Talking extension" => "Talking.md",
        "Testing extension" => "Testing.md"
    ]
)
deploydocs(repo = "github.com/bramtayl/MenuAdventures.jl.git")


text:
function testBulgeCPU()
    testBulge(false)
end

function testBulgeGPU()
    	testBulge(true)
end

function testBulge(useGPU)
    to = TimerOutput()

    myFloat = Float64
    input = getDefaultInput(myFloat, "../Athesis.toml")

    # calculate bulging level from input parameters
    bulge = input.Lx * input.Lx * input.constRecharge /
    (8.0 * input.K0 * input.boundaryPressure[1] * input.Δx * input.Δy)

    useGPU = false
    simulation = initSimulation(input, useGPU, myFloat, to)
    runSimulation!(simulation, to)

    ix = Int(ceil(simulation.grid.nx/2))
    iy = Int(ceil(simulation.grid.ny/2))
    iz = simulation.grid.nz
    head = simulation.state.hⁿ⁺¹[ix,iy,iz]

    diff = abs(head - (input.boundaryPressure[1] + bulge))

    return diff < 1.0e-4
end

function testBulgeGPUvsCPU()

    to = TimerOutput()
    myFloat = Float64
    input = getDefaultInput(myFloat, "../Athesis.toml")

    # CPU
    simulationCPU = initSimulation(input, false, myFloat, to)
    convergedCPU, nIterCPU = runSimulation!(simulationCPU, to)

    # GPU
    simulationGPU = initSimulation(input, true, myFloat, to)
    convergedGPU, nIterGPU = runSimulation!(simulationGPU, to)

    # equal convergence
    @test convergedCPU == convergedGPU
    @test nIterCPU == nIterGPU

    # same result
    absDiffHeads = abs.(simulationCPU.state.hⁿ⁺¹ - simulationGPU.state.hⁿ⁺¹)
    maxDiff, = findmax(absDiffHeads)

    @test maxDiff ≈ 0.0 atol=1.0e-09

    return true
end

@testset "Validation" begin
    @info "Running validation tests on CPU..."
    @test testBulgeCPU()

    if has_cuda()
        @info "Running validation tests on GPU..."
        @test testBulgeGPU()
        @test testBulgeGPUvsCPU()
    end
end


text:
ccall(:jl_exit_on_sigint, Void, (Cint,), 0)

push!(LOAD_PATH, "..")
using JukaiNLP: DepParser, Perceptron, Unlabeled, Labeled, FeedForward, UnlabeledFeedForward
using JukaiNLP: readconll, train!, decode, evaluate, toconll, initmodel!, MyAdaGrad
using JukaiNLP.DepParsing: Token, readconll_withprobs
using Merlin
# using JLD
using DocOpt

doc = """shift-reduce parser

Usage:
    demo.jl train (--labeled | --unlabeled) --nn --worddict <worddict> <train_path> <model_path> [<test_path>] [--iter=<iter>] [--embedfile=<embed_file>] [--batchsize=<batch>] [--evaliter=<eval>] [--nonlinearity=<nonlinear>] [--hiddenlayers=<hidden>] [--embedlayers=<embedlayers>] [--optimizer=<opt>] [--learnrate=<rate>] [--momentum=<momentum>] [--tagdict=<tagdict>] [--train-tag-probs=<trainprobs>] [--test-tag-probs=<testprobs>] [--use-topktags=<k>]
    demo.jl train (--labeled | --unlabeled) --struct-nn <train_path> <local_model> <model_path> [<test_path>] [--iter=<iter>] [--batchsize=<batch>] [--evaliter=<eval>] [--beamsize=<beam>] [--nonlinearity=<nonlinear>] [--optimizer=<opt>] [--learnrate=<rate>] [--momentum=<momentum>] [--tagdict=<tagdict>] [--train-tag-probs=<trainprobs>] [--test-tag-probs=<testprobs>] [--use-topktags=<k>]
    demo.jl train (--labeled | --unlabeled) --perceptron --worddict <worddict> <train_path> <model_path> [<test_path>] [--iter=<iter>]

    demo.jl test <test_path> <model_path>
    demo.jl (<sent> | -) <model_path>

Options:
    <train_path>    path to CoNLL format file to use for training
    <test_path>     path to CoNLL format file to use for evaluation
    <model_path>    path to output resulting parser model

"""
    # --iter          number of training iteration
    # --embedfile     path to pretrained embedding file (used in --nn)
    # --batchsize     batch size (used in --nn)
    # --evaliter      run evaluation on test data every after this number of iteration

args = docopt(doc)

worddict = args["<worddict>"]
parsertype = args["--labeled"] ? Labeled :
             args["--unlabeled"] ? Unlabeled :
             nothing
iter = args["--iter"] != nothing ? parse(Int, args["--iter"]) :
       args["--perceptron"] ? 20 :
       args["--nn"] ? 20000 : nothing
embedfile = args["--embedfile"] != nothing ? args["--embedfile"] : ""
trainfile = args["<train_path>"]
testfile = args["<test_path>"]
modelpath = args["<model_path>"]
localmodel = args["<local_model>"]
tagfile = args["--tagdict"]
trainprob = args["--train-tag-probs"]
testprob = args["--test-tag-probs"]
topk = args["--use-topktags"] != nothing ? parse(Int, args["--use-topktags"]) : 1
batchsize = args["--batchsize"] != nothing ? parse(Int, args["--batchsize"]) : 10000
beamsize = args["--beamsize"] != nothing ? parse(Int, args["--beamsize"]) : 32
evaliter = args["--evaliter"] != nothing ? parse(Int, args["--evaliter"]) : 100
nonlinear = args["--nonlinearity"] == nothing ? relu :
            args["--nonlinearity"] == "relu" ? relu :
            args["--nonlinearity"] == "tanh" ? tanh :
            args["--nonlinearity"] == "sigmoid" ? sigmoid :
            throw("no support for $(args["--nonlinearity"])")

hiddenlayers = args["--hiddenlayers"] == nothing ? [1024] :
            map(x -> parse(Int, x), split(args["--hiddenlayers"], "x"))
embedlayers = args["--embedlayers"] == nothing ? [50,50,50] :
            map(x -> parse(Int, x), split(args["--embedlayers"], "x"))

learnrate = args["--learnrate"] != nothing ? parse(Float64, args["--learnrate"]) : 0.001
momentum = args["--momentum"] != nothing ? parse(Float64, args["--momentum"]) : 0.9

opt = args["--optimizer"] == nothing ? SGD(learnrate, momentum=momentum) :
            lowercase(args["--optimizer"]) == "sgd" ? SGD(learnrate, momentum=momentum) :
            lowercase(args["--optimizer"]) == "adagrad" ? MyAdaGrad(learnrate) :
            throw("no support for $(args["--optimizer"])")
if args["train"]

    if args["--nn"]
        parser = DepParser(parsertype, worddict)
        # trainsents = readconll_withprobs(parser, tagfile, trainfile, trainprob, topk=topk)
        # testsents = testfile == nothing ? Vector{Token}[] :
        # readconll_withprobs(parser, tagfile, testfile, testprob, topk=topk, train=false)
        trainsents = readconll(parser, trainfile)
        testsents = testfile == nothing ? Vector{Token}[] :
                        readconll(parser, testfile, train=false)
        train!(FeedForward, parser, trainsents, testsents, embed=embedfile,
                nonlinear=nonlinear, embedsizes=embedlayers, hiddensizes=hiddenlayers,
                topktags=topk!=1, opt=opt, iter=iter, batchsize=batchsize,
                evaliter=evaliter, outfile=modelpath)

    elseif args["--struct-nn"]
        parser = open(deserialize, localmodel)
        topk = typeof(parser.model.tag_f) == Linear ? -1 : 1
        trainsents = readconll_withprobs(parser, tagfile, trainfile, trainprob, topk=topk)
        testsents = testfile == nothing ? Vector{Token}[] :
        readconll_withprobs(parser, tagfile, testfile, testprob, topk=topk, train=false)
        train!(StructuredFeedForward, parser, trainsents, testsents, beamsize=beamsize,
                opt=opt, iter=iter, batchsize=batchsize, evaliter=evaliter, outfile=modelpath)

    elseif args["--perceptron"]
        parser = DepParser(parsertype, worddict)
        trainsents = readconll(parser, trainfile)
        testsents = testfile == nothing ? Vector{Token}[] :
                        readconll(parser, testfile)
        train!(Perceptron, parser, trainsents, testsents,
            iter=iter, outfile=modelpath)
    end

elseif args["test"]
    # parser = load(args["<model_path>"], "parser")
    parser = open(deserialize, modelpath)
    res = parser(testfile)
    for s in res
        toconll(s)
    end
    evaluate(parser, res)

elseif args["<sent>"] == "-"
    throw("yet to be supported")

elseif args["<sent>"] != nothing
    throw("yet to be supported")
end

max_stars_repo_name: hshindo/JukaiNLP

max_issues_repo_name: hshindo/JukaiNLP

max_forks_repo_name: hshindo/JukaiNLP


text:
module JRAF

    using Nemo
    using Cuba
    using WignerSymbols
    using Legendre
    using SphericalHarmonics
    using Base.Threads
    using DelimitedFiles
    # export JULIA_NUM_THREADS=8
    sprec = 750;
    ARBF = ArbField(sprec);
    ACBF = AcbField(sprec);
    CF = ComplexField(sprec);
    RF = RealField(sprec);
    export sprec, RF, CF, Plm, Nlm, computeYlm, SphericalHarmonics

    include("math.jl")
    include("angular_coefficients.jl")
    include("radial_coefficients.jl")
    include("special_functions.jl")
    include("gaux_p12_bsrep.jl")
    include("gaux_p123_bsrep.jl")
    include("gaux_p123_rec.jl")
    include("cgaux_p123_num.jl")
    include("sto_mol_integ_one_elect.jl")

end


text:
__precompile__()

module IScatterSpectrum

using SpecialFunctions,Polynomials

import PhysicalConstants.MKS: Boltzmann,ElectronCharge,MassElectron,SpeedOfLight,
                              UnifiedAtomicMass,VacuumPermittivity

friedconte(z) = im*sqrt(pi)*erfcx(-im*z)

struct Faddeeva
    n::Int
    l::Float64
    a::Poly{Float64}
    function Faddeeva(n=16)
        m  = 2*n
        l  = sqrt(n/sqrt(2))

        theta = (-m+1:m-1)*pi/m
        t = l*tan.(theta/2)
        f = [0; exp.(-t.^2).*(l*l + t.^2)]
        a = real.(fft(fftshift(f)))[2:n-1]/(2*m)   # no flip because of convention in polyval
        new(n, l, Poly(a))
    end
end

const deff = Faddeeva()

"""
    faddeeva(z, f=deff)

   The Faddeeva function is defined as:

   w(z) = exp(-z^2) * erfc(-im*z)

   where erfc(x) is the complex complementary error function. Also,

   friedconte(z) = im * sqrt(pi) * faddeeva(z)

   where friedconte(z) is the plasma dispersion function as defined in [2].

   w = faddeeva(z, f) reuses previously calculated coefficients in struct f
   where f is initialized like f=Faddeeva(n) with n specifying the number of terms
   to truncate the expansion (see (13) in [1]). n = 16 is used as default, i.e. for f=Faddeeva().

   Example:
       x = linspace(-10,10,1001); [X,Y] = meshgrid(x,x);
       W = faddeeva(complex(X,Y));

   Reference:
   [1] J.A.C. Weideman, "Computation of the Complex Error Function," SIAM
       J. Numerical Analysis, pp. 1497-1518, No. 5, Vol. 31, Oct., 1994
       Available Online: http://www.jstor.org/stable/2158232

   [2] B.D. Fried, S.D. Conte, The plasma dispersion function. New York Acadenic Press, 1961.
"""
function faddeeva(z, f::Faddeeva=deff)
    if real(z)==0      # purely imaginary input
        exp.(-z.^2).*erfc(imag(z))
    else               # complex valued input
        # make sure input in the upper half-plane (positive imag. values)
        zz = imag(z)<0? conj(z):z
        yy = (f.l + im*zz)/(f.l - im*zz)
        w  = 2*polyval(f.a, yy)/(f.l - im*zz)^2 + (1/sqrt(pi))/(f.l - im*zz)
        # convert the upper half-plane results to the lower half-plane if necessary
        imag(z)<0? conj(2*exp(-zz^2) - w):w
    end
end

const fnormf = SpeedOfLight*sqrt(MassElectron/(2*Boltzmann))/2
const debyef = (4*π/ElectronCharge/SpeedOfLight)^2*VacuumPermittivity*Boltzmann

"gfe(b in T) -> electron gyrofrequency in Hz"
gfe(b) = ElectronCharge*b/MassElectron

"""
    ScatterVolume(tfeff, system, b, alpha)

has parameters of the scattering volume that are constant in a fit.
"""
struct ScatterVolume{T<:Real}
    tfeff::T    # effective transmitter frequency in Hz
    fnormc::T   # factor for frequency normalization f -> theta
    debyec::T   # Debye length correction factor

    system::T   # spectrum scale factor (depends on tx power, volume size, ... )

    b::T        # magnetic field strength in Tesla
    gfe::T      # electron gyrofrequency
    phiel::T    # normalized electron gyrofrequency
    alpha::T    # angle between b and scattering vector in rad
    sinfac::T   # correction factor for magnetic field
    cosal::T    #
    cosal2::T

    function ScatterVolume{T}(tfeff::T, system::T, b::T, alpha::T) where T<:Real
        fnormc = fnormf/tfeff
        debyec = debyef*tfeff^2
        gf     = gfe(b)
        phiel  = fnormc*gf
        sinfac = 0.5*(sin(alpha)/phiel)^2
        cosal  = cos(alpha)
        new(tfeff, fnormc, debyec, system, b, gf, phiel, alpha, sinfac, cosal, cosal*cosal)
    end
end
ScatterVolume(tfeff::T, system::T, b::T, alpha::T) where {T<:Real} =
    ScatterVolume{T}(tfeff, system, b, alpha)

struct Ion{T<:Real}
    mass::T   # ion mass in amu
    srmr::T   # sqrt of mass ratio ion/electron
end
Ion(m=15.999) = Ion(m, sqrt(m*UnifiedAtomicMass/MassElectron))

"""
Plasma(ne, te, ti, ::ScatterVolume, [ions, fractn]

has parameters of the plasma that are constant over all frequencies of the spectrum.
"""
struct Plasma{T<:Real}
    ne::T
    te::T
    ti::T
    fractn1::T
    fractn::NTuple{2, T}    # up to two more ions
    ions::NTuple{3, Ion{T}}
    nion::Int

    spectrsf::T
    debyen::T   # normalized Debye length correction factor
    sqte::T
    sqti::T
    thetc::T
    tex::T
    tefac::T
    dpsidc::T

    # parameter for spectrum at zero frequency:
    cdtine::T
    ym::T

    function Plasma{T}(ne::T, te::T, ti::T, s::ScatterVolume{T}, fractn, ions, nion) where T<:Real
        sqti    = sqrt(ti)
        thetc   = s.fnormc/sqti
        tex     = exp(-te*s.sinfac)
        tefac   = 0.5/te + s.sinfac/tex
        fractn1 = one(ne) - fractn[1] - fractn[2]
        dpsidc  = thetc*ions[1].srmr
        cdtine  = ti*s.debyec + ne
        ym      = ti*ne + te*cdtine
        new(ne, te, ti, fractn1, fractn, ions, nion,
            s.system*ne, s.debyec/ne, sqrt(te), sqti, thetc,
            tex, tefac, dpsidc, cdtine, ym)
    end
end
Plasma(ne::T, te::T, ti::T,
       s::ScatterVolume{T}, nion=1, ions=(Ion(16*one(ne)), Ion(31*one(ne)), Ion(one(ne))),
       fractn=(zero(ne),zero(ne))) where {T<:Real} =
           Plasma{T}(ne, te, ti, s, fractn, ions, nion)

struct CollisionalPlasma{T<:Real}
    p::Plasma{T}
    cf::T
end

"""
    ionadm(theta, z=friedconte(theta))
-> unmagnetized ion admittance for real normalized frequency
"""
ionadm{T<:Real}(theta::T, z::Complex{T}=friedconte(theta)) = Complex(theta*imag(z), 1 + theta*real(z))

"""
    ionadm(freq, ::Plasma)
-> unmagnetized ion admittance for real unnormalized frequency and plasma parameters
"""
function ionadm{T<:Real}(freq::T, p::Plasma{T}, s::ScatterVolume{T})
    fnormt = s.fnormc/p.sqti
    thetai = fnormt*p.ions[1].srmr*freq
    yi     = ionadm(thetai)
    for k in 2:p.nion
        if k==2
            yi *= p.fractn1
        end
        thetai = fnormt*p.ions[k].srmr*freq
        yi    += p.fractn[k-1]*ionadm(thetai)
    end
    yi
end

"""
    ionadm(theta, z=friedconte(theta))
-> unmagnetized ion admittance for complex normalized frequency.
"""
function ionadm{T<:Complex}(theta::T, z::T=friedconte(theta))
    psizm = imag(theta)*(real(z)^2 + imag(z)^2)
    yd    = 1 - imag(theta)*(2*imag(z) - psizm)
    Complex(real(theta)*(imag(z) - psizm)/yd, 1 + real(theta)*real(z)/yd)
end

"""
    ionadm(freq, ::CollisionalPlasma)
-> unmagnetized ion admittance for unnormalized frequency and collisional plasma parameters
"""
function ionadm{T<:Real}(freq::T, coll::CollisionalPlasma{T}, s::ScatterVolume{T})
    fnormt = s.fnormc/coll.p.sqti
    thetai = fnormt*coll.p.ions[1].srmr*Complex(freq, coll.cf)
    yi     = ionadm(thetai)
    for k in 2:coll.p.nion
        if k==2
            yi *= coll.p.fractn1
        end
        thetai = fnormt*coll.p.ions[k].srmr*Complex(freq, coll.cf)
        yi    += coll.p.fractn[k-1]*ionadm(thetai)
    end
    yi
end

"""
    elecadm(theta, tex, z)
-> magnetized electron admittance for real normalized frequency
"""
function elecadm{T<:Real}(theta::T, tex::T, z::Complex{T}=friedconte(theta))
    tt = tex*theta
    Complex(tt*imag(z), 1 + tt*real(z)), imag(z)
end

# helper function, this part of the spectrum code is the same for non-collisional and collisional
"""
    espec(freq, ::Plasma, ::ScatterVolume, yi::Complex)
-> power spectral density for non-zero real frequency
"""
function espec{T<:Real}(freq::T, p::Plasma{T}, s::ScatterVolume{T}, yi::Complex{T})
    yiti   = yi/p.ti + im*p.debyen

    # electron admittance:
    thetae = s.fnormc/(p.sqte*s.cosal)*freq
    ye,    = elecadm(thetae, p.tex)
    yete   = ye/p.te

    the2   = thetae^2
    texthe = p.tex*thetae

    ym     = abs2(yete+yiti)
    yn     = abs2(yete)*real(yi) + abs2(yiti)*real(ye)

    (p.spectrsf/freq)*(yn/ym)
end

"""
    pwrspec(freq, ::Plasma, ::ScatterVolume)
-> power spectral density for non-zero real frequency
"""
function pwrspec{T<:Real}(freq::T, p::Plasma{T}, s::ScatterVolume{T})
    yi = ionadm(freq, p, s)
    espec(freq, p, s, yi)
end

"""
    pwrspec(freq, ::CollisionalPlasma, ::ScatterVolume)
-> power spectral density for non-zero frequency and collisional plasma
"""
function pwrspec{T<:Real}(freq::T, coll::CollisionalPlasma{T}, s::ScatterVolume{T})
    yi = ionadm(freq, coll, s)
    espec(freq, coll.p, s, yi)
end

# helper function
function yn0{T<:Real}(p::Plasma{T}, s::ScatterVolume{T}, yir::T)
    tineq  = p.sqti*p.ne^2*yir*p.ti
    sqtec  = p.cdtine*p.sqte*sqrt(pi)/s.cosal2 # zi(1) in fortran code = sqrt(π) for real zero frequency
    stecx  = sqtec*p.tex*p.te
    tineq + stecx*p.cdtine
end

"""
    pwrspec(::Plasma)
-> power spectral density for zero real frequency
"""
function pwrspec(p::Plasma, s::ScatterVolume)
    yir = sqrt(π)*p.ions[1].srmr
    for k in 2:p.nion
        if k==2
            yir *= p.fractn1
        end
        yir += p.fractn[k-1]*sqrt(π)*p.ions[k].srmr
    end
    yn = yn0(p, s, yir)
    s.fnormc*p.spectrsf*yn/p.ym^2
end

function yirc(coll::CollisionalPlasma, kion::Int)
    psiion = coll.p.dpsidc*coll.cf
    z      = friedconte(Complex(zero(coll.p.ne), psiion))
    psiz   = 1 - psiion*imag(z)
    psizm  = coll.p.ions[kion].srmr/psiz
    psizm*imag(z)
end

"""
    pwrspec(::CollisionalPlasma)
-> power spectral density for zero frequency and collisional plasma
"""
function pwrspec(coll::CollisionalPlasma, s::ScatterVolume)
    yir = yirc(coll, 1)
    for k in 2:p.nion
        if k==2
            yir *= p.fractn1
        end
        z    = friedconte(Complex(zero(coll.p.ne), p.thetc*coll.p.ions[k].srmr*coll.cf))
        yir += p.fractn[k-1]*yirc(coll, k)
    end

    yn = yn0(p, s, yir)
    s.fnormc*coll.p.spectrsf*yn/p.ym^2
end

end   # module IScatterSpectrum


text:
include("../src/prime_range.jl")

function test_prime_range_length()
	print("prime_range.length...")

  @test length(prime_range(0)) == 0
  @test length(prime_range(1)) == 0
  @test length(prime_range(2)) == 1
  @test length(prime_range(100)) == 25
	println("PASS")
end

function test_prime_range_elements()
	print("prime_range.elements...")

  @test prime_range(2) == [2]
  @test prime_range(10) == [2, 3, 5, 7]
  @test prime_range(100)[20] == 71
	println("PASS")
end

function test_prime_range()
	test_prime_range_length()
	test_prime_range_elements()

	println("")
end


text:
using Documenter
using Weave
using BoteSalvatICX
using Cairo
using Fontconfig

function addNISTHeaders(htmlfile::String)
    # read HTML
    html = transcode(String,read(htmlfile))
    # Find </head>
    i = findfirst(r"</[Hh][Ee][Aa][Dd]>", html)
    # Already added???
    j = findfirst("nist-header-footer", html)
    if isnothing(j) && (!isnothing(i))
        # Insert nist-pages links right before </head>
        res = html[1:i.start-1]*
            "<link rel=\"stylesheet\" href=\"https://pages.nist.gov/nist-header-footer/css/nist-combined.css\">\n"*
            "<script src=\"https://pages.nist.gov/nist-header-footer/js/jquery-1.9.0.min.js\" type=\"text/javascript\" defer=\"defer\"></script>\n"*
            "<script src=\"https://pages.nist.gov/nist-header-footer/js/nist-header-footer.js\" type=\"text/javascript\" defer=\"defer\"></script>\n"*
            html[i.start:end]
        write(htmlfile, res)
        println("Inserting NIST header/footer into $htmlfile")
    end
    return htmlfile
end

rm("build",force=true,recursive=true)

weaveit(name) = weave(joinpath("src", "$name"), out_path=joinpath("src", "$(splitext(name)[1]).md"), doctype="github")

names = ( "Example.ipynb", )

weaveit.(names)

makedocs(
    modules = [BoteSalvatICX],
    sitename = "BoteSalvatICX.jl",
    pages = [ "Home" => "index.md", "Using BoteSalvatICX" => "Example.md" ]
    )

map(name->rm(joinpath("src","$(splitext(name)[1]).md")), names)
addNISTHeaders(joinpath(@__DIR__, "build","index.html"))
addNISTHeaders.(map(name->joinpath(@__DIR__, "build", splitext(name)[1], "index.html"), names))
#deploydocs(
#    repo   = "github.com/USNISTGOV/BoteSalvatICX.jl.git",
#    branch = "nist-pages",
#    versions = "v^",
#)
# deploydocs(repo = "github.com/USNISTGOV/BoteSalvatICX.jl.git")


text:
"""
    s_metric(g;norm=true)

Return the normalised s-metric of `g`.

The s-metric is defined as the sum of the product of degrees between pair of nodes
for every edge in `g`. [Ref](https://arxiv.org/abs/cond-mat/0501169)
It is normalised by the maximum s_metric obtained from the family of graph
with similar degree distribution.
s_max is computed from an approximation formula as in https://journals.aps.org/pre/pdf/10.1103/PhysRevE.75.046102
If `norm=false`, no normalisation is performed.

# Examples
```jldoctest
julia> using LightGraphs

julia> s_metric(star_graph(4))
0.6
```
"""

function s_metric(g::AbstractGraph{T};norm=true) where T
    s = zero(T)
    for e in edges(g)
        s += degree(g,src(e)) * degree(g,dst(e))
    end
    if norm
        sm = sum(degree(g).^3)/2
        return s/sm
    else
        return s
    end
end


text:
using Sigma
using Base.Test

# Constrain an array of integers to be symmetric
function issymmetric(x)
  n = length(x)
  middle = div(n,2)+1
  x[1:middle-1] == x[length(x):-1:middle]
end

Xs = mvuniform(0,25,100)
sample = rand(Xs, issymmetric(Xs))


text:
module Utils

using Base.Meta
using GaloisFields
using GaloisFields: PrimeField
using Distributions
using Nemo
using AbstractAlgebra
using ..NTT
using ..NTT: modulus, degree
using Primes
using Mods
using StructArrays
using OffsetArrays
using BitIntegers
using Hecke

macro fields_as_locals(a)
    @assert isexpr(a, :(::))
    sym, T = a.args
    T = getfield(__module__, T)
    quote
        $([:(local $(esc(fname))) for fname in fieldnames(T)]...)
        let s = $(esc(sym))
            $([:($(esc(fname)) = s.$fname) for fname in fieldnames(T)]...)
        end
    end
end

# HACK - DiscreteUniform is the default for types, but it'd be nice to
# allow this.
Distributions.DiscreteUniform(T::Type) = T
Distributions.DiscreteUniform(R::Nemo.FmpzModRing) = R

# HACK
Base.BroadcastStyle(::Type{<:OffsetArray{<:Any, <:Any, T}}) where T<:StructArray = Base.BroadcastStyle(T)

# HACK
Base.convert(::Type{Integer}, x::PrimeField) = x.n
Base.convert(::Type{Integer}, x::Nemo.fmpz_mod) = x.data

# HACK
Base.convert(::Type{<:StructArray{T}}, A::AbstractArray{T}) where {T} =
    StructArray(A)

# Hack
Base.convert(::Type{Integer}, x::AbstractAlgebra.Generic.Res{fmpz}) = Nemo.lift(x)

# Hack
Base.div(a::Nemo.fmpz, b::BigInt, r::RoundingMode{:NearestTiesAway}) =
    div(BigInt(a), b, r)

# Hack
function Base.convert(T::Type{<:PrimeField}, x::PrimeField)
    x = convert(Integer, x)
    x > modulus(T) && throw(InexactError(:convert, T, x))
    T(x)
end

# Hack
Primes.isprime(x::Int256) = Primes.isprime(big(x))

# Hack
(::NmodPolyRing)(x::AbstractAlgebra.Generic.Res{fmpz_mod_poly}) = x

# Hack
Base.gcdx(x::nmod_poly, y::nmod_poly) = Hecke.rresx(x, y)

# Hack
Nemo.nmod_poly(a::Integer, b::Integer) = Nemo.nmod_poly(UInt64(a), UInt64(b))

# Hack (https://github.com/Nemocas/Nemo.jl/issues/668)
function (a::FqNmodFiniteField)(b::Nemo.nmod)
    a(data(b))
end

# Hack
function Base.oftype(R::AbstractAlgebra.Generic.Res{fmpz_mod_poly}, e::AbstractAlgebra.Generic.Res{nmod_poly})
    ℤx = Nemo.ZZ["x"][1]
    @assert Nemo.lift(ℤx, modulus(parent(R))) == Nemo.lift(ℤx, modulus(parent(e)))
    parent(R)(Nemo.lift(ℤx, data(e)))
end

struct UsageError
    msg::String
end

# For now
Base.widen(::Type{Int128}) = Int256

Base.promote_rule(::Type{Nemo.fmpz_mod}, ::Type{UInt64}) = Nemo.fmpz_mod

# HACK
Base.convert(::Type{OffsetArray{T, N, A1}}, A::OffsetArray{T, N, A1}) where {T,N,A1} = A
function Base.convert(::Type{OffsetArray{T, N, A1}}, A::OffsetArray{T, N, A2}) where {T,N,A1,A2}
    OffsetArray{T, N, A1}(convert(A1, A.parent), A.offsets)
end

#Hack
function Base.similar(::Type{<:StructArray{S, N}}, T::Type) where {S, N}
    C = Tuple{map(x->Array{x, N}, StructArrays.staticschema(T).parameters)...}
    I = StructArrays.index_type(C)
    StructArray{T, N, C, I}
end

end


text:
using ADCME
using PyPlot
using AdFem
using MAT

m = 20
n = 20
h = 1/n
bdnode = bcnode("all", m, n, h)
xy = gauss_nodes(m, n, h)
xy_fem = fem_nodes(m, n, h)

not_bdnode = ones(Bool, (m+1)*(n+1))
not_bdnode[bdnode] .= false
not_bdnode = findall(not_bdnode)

σ0 = 0.0
prior = 0.01
if length(ARGS)==2
    global σ0 = parse(Float64, ARGS[1])
    global prior = parse(Float64, ARGS[2])
end

DIR = "uq_nn2-$σ0-$prior"
rm(DIR, force=true, recursive=true)
mkdir(DIR)

function poisson1d(κ)
    κ = compute_space_varying_tangent_elasticity_matrix(κ, m, n, h)
    K = compute_fem_stiffness_matrix1(κ, m, n, h)
    K, _ = fem_impose_Dirichlet_boundary_condition1(K, bdnode, m, n, h)
    rhs = compute_fem_source_term1(ones(4m*n), m, n, h)
    rhs[bdnode] .= 0.0
    sol = K\rhs
end

function κfc(x, θ)
    xc = fem_to_fvm(x, m, n, h)
    xc = reshape(repeat(xc, 1, 4), (-1,))
    @assert length(xc)==4*m*n
    0.1 + abs(squeeze(fc(reshape(xc, (-1,1)), [20,20,20,1], θ)))
end

function param(x, θ)
    xc = fem_to_fvm(x, m, n, h)
    xc = reshape(repeat(xc, 1, 4), (-1,))
    return 0.1 + 1/(1+xc^2) + 100xc^2
end

function poisson1d_nn(fn)
    θ = missing
    function f(θ, u)
        κ = fn(u, θ)
        K = compute_space_varying_tangent_elasticity_matrix(κ, m, n, h)
        K = compute_fem_stiffness_matrix1(K, m, n, h)
        K, _ = fem_impose_Dirichlet_boundary_condition1(K, bdnode, m, n, h)
        rhs = K * u
        rhs = rhs[not_bdnode] - compute_fem_source_term1(ones(4m*n), m, n, h)[not_bdnode]
        rhs, gradients(rhs, u)
    end
    if fn==κfc
        θ = Variable(ae_init([1, 20, 20, 20, 1]))
    end
    u = newton_raphson_with_grad(f, constant(zeros((m+1)*(n+1))) ,θ)
    κnn = fn(u, θ)
    return u, κnn, θ
end

function jacobian(u, θ)
    κ = κfc(u, θ)
    K = compute_space_varying_tangent_elasticity_matrix(κ, m, n, h)
    K = compute_fem_stiffness_matrix1(K, m, n, h)
    K, _ = fem_impose_Dirichlet_boundary_condition1(K, bdnode, m, n, h)
    rhs = K * u
    rhs = rhs[not_bdnode] - compute_fem_source_term1(ones(4m*n), m, n, h)[not_bdnode]
    J = gradients(rhs, u)
    F = gradients(rhs, κ)
    -J\F
end

ADCME.options.newton_raphson.verbose=true
ADCME.options.newton_raphson.rtol = 1e-4
ADCME.options.newton_raphson.tol = 1e-4
y, κ, _ = poisson1d_nn(param)

sess = Session(); init(sess)
SOL = run(sess, y)

umax = maximum(SOL)
xc = LinRange(0, umax, 100)|>Array
u1 = @. 0.1 + 1/(1+xc^2) + 100xc^2

y, κnn, θ  = poisson1d_nn(κfc)

using Random; Random.seed!(233)
idx = fem_randidx(50, m, n, h)
obs = y[idx]
OBS = SOL[idx] .* ( 1. .+ σ0*randn(length(idx)))
loss = sum((obs-OBS)^2)

init(sess)
BFGS!(sess, loss, 200)

close("all")
figure(figsize=(10,4))
subplot(121)
visualize_scalar_on_fem_points(SOL, m, n, h)
subplot(122)
Yvalue = run(sess, y)
visualize_scalar_on_fem_points(Yvalue, m, n, h)
plot(xy_fem[idx,1], xy_fem[idx,2], "o", c="red", label="Observation")
legend()
savefig("$DIR/sol_compare.png")

close("all")
figure(figsize=(10,4))
subplot(121)
κvalue = run(sess, κ)
visualize_scalar_on_gauss_points(κvalue, m, n, h)
title("Exact \$K(x, y)\$")
subplot(122)
κnnvalue = run(sess, κnn)
visualize_scalar_on_gauss_points(κnnvalue, m, n, h)
title("Estimated \$K(x, y)\$")
savefig("$DIR/K_compare.png")

u0 = constant(LinRange(0,umax,100)|>Array )
z = 0.1 + abs(squeeze(fc(reshape(u0, (-1,1)), [20,20,20,1], θ)))
close("all")
plot(run(sess, u0), run(sess, z), "--")
plot(xc, u1, "-")
savefig("$DIR/ucurve.png")

u_est = run(sess, y)
H = jacobian(y, θ)
H = run(sess, H)[idx,:]
y_ = OBS
hs = run(sess, obs)
s = run(sess, κnn)

R = (prior)^2*diagm(0=>ones(length(obs)))
Q = (1e-2)^2*diagm(0=>ones(length(κnn)))
μ, Σ = uqnlin(y_, hs, H, R, s, Q)
# error()
σ = diag(Σ)

close("all")
u1 = fem_to_fvm(u_est, m, n, h)
z1 = μ[1:4:end]
σ1 = 2σ[1:4:end]
Is = sortperm(u1)
u1 = u1[Is]
z1 = z1[Is]
σ1 = σ1[Is]

xc = LinRange(0, umax, 100)|>Array
u2 = @. 0.1 + 1/(1+xc^2) + 100xc^2
plot(xc, u2, "-", linewidth=2, label="Referencee")

plot(run(sess, u0), run(sess, z), "--", linewidth=2, label="Estimated")
plot(u1, z1, "--", linewidth=2, label="Posterior Mean")

fill_between(u1, z1-σ1, z1+σ1, alpha=0.3, color="orange", label="Uncertainty Region")
xlim(minimum(u1), maximum(u1))
legend()
xlabel("\$u\$")
ylabel("\$K(u)\$")
savefig("$DIR/errorbar.png")

matwrite("$DIR/data.mat",
    Dict(
        "SOL"=>SOL,
        "Yvalue"=>Yvalue,
        "kappavalue"=>κvalue,
        "u_est"=>u_est,
        "H"=>H,
        "yobs"=>y_,
        "hs"=>hs,
        "s"=>s,
        "sigma"=>σ,
        "mu"=>μ
    ))

max_stars_repo_name: ADCMEMarket/ADCMEImages

max_issues_repo_name: ADCMEMarket/ADCMEImages

max_forks_repo_name: ADCMEMarket/ADCMEImages


text:
module runsim

export Simulation

import ..exclusionprocess: ExclusionProcess
import ..epstate: EPState
import ..boundarycondition: BoundaryCondition

"""
Basic simulation type.
"""
mutable struct Simulation
    ep::Type{S} where {S<:ExclusionProcess}
    bc::Type{T} where {T<:BoundaryCondition}
    state::EPState
    steps::Int
end

end

max_stars_repo_name: AustenLamacraft/ExclusionProcesses

max_issues_repo_name: AustenLamacraft/ExclusionProcesses

max_forks_repo_name: AustenLamacraft/ExclusionProcesses


text:
module VLMinorityGameKit

# include -
include("Include.jl")

# method and type exports -
# types -
export VLAbstractGameWorld
export VLAbstractGameAgent
export VLAbstractGameStrategy

export VLBasicMinorityGameWorld
export VLBasicMinorityGameAgent
export VLBasicMinorityGameStrategy

export VLThermalMinorityGameWorld
export VLThermalMinorityGameAgent
export VLThermalMinorityGameStrategy

export VLGCMinorityGameStrategy
export VLGCMinorityGameAgent
export VLGCMinorityGameWorld

export VLMinorityGameSimulationResult

# methods -
export execute
export build_minority_game_world

# analysis functions -
export compute_price_return
export compute_autocorrelation_array

# for debug -
export design

end # module


text:
using WAV, DataFrames, CSV, PlotlyJS, DSP, PyPlot
include("spectrograms.jl")

function plot_classes(classes, frequencies)
    plot(bar(;x=classes, y=frequencies))
end

function create_audio_metadata(labels)
    audiodata = DataFrame(sampleRate=Float32[], channels=Int8[], samples=Int32[], length=Float32[], class=String[])
    looked = 0
    n = length(labels)
    for (id, lbl) in labels
        println("Analyzing file $id.wav")
        data, srate = wavread("train/Train/$id.wav")
        channels = size(data,2)
        samples = size(data,1)
        len = samples / srate
        row = [srate channels samples len lbl]
        push!(audiodata, row)
        looked += 1
        println("$(convert(Int8, round(looked / n * 100)))% Analyzed")
    end
    return audiodata
end

# Read class data
df = CSV.read("train/train.csv", header=true)
labels = gen_label_dict("train/train.csv")

# Plot class distribution
classdata = by(df, :Class, nrow)
# plot_classes(classdata[:,1], classdata[:,2])

# audio_metadata = create_audio_metadata(labels)
audio_metadata = CSV.read("metadata.csv", header=true)
# describe(audio_metadata, stats=:all)
# CSV.write("metadata.csv", audio_metadata)

# Plot a single wav file
# data, fs = wavread("train/0.wav")
# trace1 = scatter(;x=1:size(y,1), y=y[:,1], mode="lines")
# trace2 = scatter(;x=1:size(y,1), y=y[:,2], mode="lines")
# plot([trace1, trace2])
spectr = create_spectrogram("train/Train/0.wav")
imshow(spectr)
# imshow(log10.(power(spec)), extent=[first(t), last(t),
#              fs*first(f), fs*last(f)], aspect="auto")


text:
"""
    some_elm(::Type{T})

Some arbitrary element of type `T`
"""
some_elm(::Type{T}) where T = begin
    try
        return zero(T)
    catch
        try
            return T(0)
        catch
            try
                rand(T)
            catch
                try
                    Vector{T}(undef, 1)[1]
                catch
                    throw("the type $T is not supported by JDF.jl yet. Try to update JDF.jl. If it still doesn't work after update, please submit an issue at https://github.com/xiaodaigh/JDF.jl/issues")
                end
            end
        end
    end
end

"""
    savejdf(outdir, dataframe)

    savejdf(dataframe, outdir)

Save a `DataFrame` to the `outdir`. On Julia > v1.3, a multithreaded version is
used.

The columns of the dataframe can be of the following vector types columns are
supported

    * `isbits` types e.g. `Int*`, `UInt*`, `Float*`
    * `Bool`
    * `Strings`
    * `WeakRefStrings.StringVector`
    * `CategoricalArrays`
    * `Union{Missing, T}`` for `T` support above

"""
savejdf(df::AbstractDataFrame, outdir::AbstractString) = savejdf(outdir, df)

savejdf(outdir, df::AbstractDataFrame; verbose = false) = begin
    if VERSION < v"1.3.0-rc1"
        return ssavejdf(outdir, df)
    end

    pmetadatas = Any[missing for i = 1:length(DataFrames.names(df))]

    if !isdir(outdir)
        mkpath(outdir)
    end

    # use a bounded channel to limit
    c1 = Channel{Bool}(Threads.nthreads())
    atexit(()->close(c1))

    for (i, n) in enumerate(DataFrames.names(df))
        if verbose
            println(n)
        end
        put!(c1, true)
        pmetadatas[i] = @spawn begin
            io = BufferedOutputStream(open(joinpath(outdir, string(n)) ,"w"))
            res = compress_then_write(df[!, i], io)
            close(io)
            res
        end
        take!(c1)
    end

    metadatas = fetch.(pmetadatas)

    fnl_metadata = (
        names = DataFrames.names(df),
        rows = DataFrames.size(df, 1),
        metadatas = metadatas,
        version = v"0.2"
    )

    open(joinpath(outdir, "metadata.jls"), "w") do io
        serialize(io, fnl_metadata)
    end
    fnl_metadata
end

"""
    serially save a DataFrames to the outdir
"""
ssavejdf(outdir, df::AbstractDataFrame) = begin
    pmetadatas = Any[missing for i = 1:length(DataFrames.names(df))]

    if !isdir(outdir)
        mkpath(outdir)
    end

    for i = 1:length(DataFrames.names(df))
        io = BufferedOutputStream(open(joinpath(outdir, string(DataFrames.names(df)[i])), "w"))
        pmetadatas[i] = compress_then_write(df[!, i], io)
        close(io)
    end

    fnl_metadata = (
        names = DataFrames.names(df),
        rows = DataFrames.size(df, 1),
        metadatas = pmetadatas,
        version = v"0.2"
    )

    open(joinpath(outdir, "metadata.jls"), "w") do io
        serialize(io, fnl_metadata)
    end
    fnl_metadata
end

# figure out from metadata how much space is allocated
get_bytes(metadata) = begin
    if metadata.type == String
        return max(
            metadata.string_compressed_bytes,
            metadata.string_len_bytes,
        )
    elseif metadata.type == Missing
        return 0
    elseif metadata.type >: Missing
        return max(
            get_bytes(metadata.Tmeta),
            get_bytes(metadata.missingmeta),
        )
    else
        return metadata.len
    end
end

hasfieldnames(::Type{T}) where {T} = fieldnames(T) >= 1


text:
using Distributed, Test, oneAPI

# GPUArrays has a testsuite that isn't part of the main package.
# Include it directly.
import GPUArrays
gpuarrays = pathof(GPUArrays)
gpuarrays_root = dirname(dirname(gpuarrays))
include(joinpath(gpuarrays_root, "test", "testsuite.jl"))
testf(f, xs...; kwargs...) = TestSuite.compare(f, oneArray, xs...; kwargs...)

# JuliaGPU/oneAPI.jl#108: fill! broken for Int16/Float16
const eltypes = [Int32, Int64,
                 Complex{Int16}, Complex{Int32}, Complex{Int64},
                 Float32,
                 ComplexF32]
const float16_supported = oneL0.module_properties(device()).fp64flags & oneL0.ZE_DEVICE_MODULE_FLAG_FP16 == oneL0.ZE_DEVICE_MODULE_FLAG_FP16
if float16_supported
    append!(eltypes, [#=Float16,=# ComplexF16])
end
const float64_supported = oneL0.module_properties(device()).fp64flags & oneL0.ZE_DEVICE_MODULE_FLAG_FP64 == oneL0.ZE_DEVICE_MODULE_FLAG_FP64
if float64_supported
    append!(eltypes, [Float64, ComplexF64])
end
TestSuite.supported_eltypes(::Type{<:oneArray}) = eltypes

using Random

## entry point

function runtests(f, name)
    old_print_setting = Test.TESTSET_PRINT_ENABLE[]
    Test.TESTSET_PRINT_ENABLE[] = false

    try
        # generate a temporary module to execute the tests in
        mod_name = Symbol("Test", rand(1:100), "Main_", replace(name, '/' => '_'))
        mod = @eval(Main, module $mod_name end)
        @eval(mod, using Test, Random, oneAPI)

        let id = myid()
            wait(@spawnat 1 print_testworker_started(name, id))
        end

        ex = quote
            GC.gc(true)
            Random.seed!(1)
            oneAPI.allowscalar(false)

            @timed @testset $"$name" begin
                $f()
            end
        end
        data = Core.eval(mod, ex)
        #data[1] is the testset

        # process results
        cpu_rss = Sys.maxrss()
        passes,fails,error,broken,c_passes,c_fails,c_errors,c_broken =
            Test.get_test_counts(data[1])
        if data[1].anynonpass == false
            data = ((passes+c_passes,broken+c_broken),
                    data[2],
                    data[3],
                    data[4],
                    data[5])
        end
        res = vcat(collect(data), cpu_rss)

        GC.gc(true)
        res
    finally
        Test.TESTSET_PRINT_ENABLE[] = old_print_setting
    end
end

## auxiliary stuff

# NOTE: based on test/pkg.jl::capture_stdout, but doesn't discard exceptions
macro grab_output(ex)
    quote
        mktemp() do fname, fout
            ret = nothing
            open(fname, "w") do fout
                redirect_stdout(fout) do
                    ret = $(esc(ex))
                end
            end
            ret, read(fname, String)
        end
    end
end

# @test_throw, peeking into the load error for testing macro errors
macro test_throws_macro(ty, ex)
    return quote
        Test.@test_throws $(esc(ty)) try
            $(esc(ex))
        catch err
            if VERSION < v"1.7-"
                @test err isa LoadError
                @test err.file === $(string(__source__.file))
                @test err.line === $(__source__.line + 1)
                rethrow(err.error)
            else
                rethrow(err)
            end
        end
    end
end

# Run some code on-device
macro on_device(ex...)
    code = ex[end]
    kwargs = ex[1:end-1]

    @gensym kernel
    esc(quote
        let
            function $kernel()
                $code
                return
            end

            oneAPI.@sync @oneapi $(kwargs...) $kernel()
        end
    end)
end

# helper function for sinking a value to prevent the callee from getting optimized away
@inline sink(i::Int32) =
    Base.llvmcall("""%slot = alloca i32
                     store volatile i32 %0, i32* %slot
                     %value = load volatile i32, i32* %slot
                     ret i32 %value""", Int32, Tuple{Int32}, i)
@inline sink(i::Int64) =
    Base.llvmcall("""%slot = alloca i64
                     store volatile i64 %0, i64* %slot
                     %value = load volatile i64, i64* %slot
                     ret i64 %value""", Int64, Tuple{Int64}, i)

nothing # File is loaded via a remotecall to "include". Ensure it returns "nothing".


text:
#MDDatasets: Additional math tools
#-------------------------------------------------------------------------------
#=NOTE:
These tools should eventually be moved to a separate unit.
=#

struct Point2D{TX<:Number, TY<:Number}
	x::TX
	y::TY
end

#Allows one to specify limits of a 1D range
#TODO: do we want to enforce min<=max???
#TODO: Add parameter to indicate if limits can go negative, overlap, ...??
struct Limits1D{T<:Number}
	min::T
	max::T
end
#Auto-detect type (Limits1D(min=4)):
Limits1D(min::TL, max::TH) where {TL<:Number,TH<:Number} = Limits1D{promote_type(TL,TH)}(min, max)
Limits1D(min::T, ::Nothing) where T<:Number = Limits1D(min, typemax(T))
Limits1D(::Nothing, max::T) where T<:Number = Limits1D(typemin(T), max)
Limits1D(;min=nothing, max=nothing) = Limits1D(min, max)
Limits1D(r::AbstractRange) = Limits1D(minimum(r), maximum(r)) #TODO: is it preferable to use rng[1/end]?

#Constructor with forced type (Limits1D{Float32}(min=4)):
(::Type{Limits1D{T}})(min::Number, ::Nothing) where T<:Number = Limits1D{T}(convert(T, min), typemax(T))
(::Type{Limits1D{T}})(::Nothing, max::Number) where T<:Number = Limits1D{T}(typemin(T), convert(T, max))
(::Type{Limits1D{T}})(;min=nothing, max=nothing) where T<:Number = Limits1D{T}(min, max)
(::Type{Limits1D{T}})(r::AbstractRange) where T<:Number = Limits1D{T}(convert(T,minimum(r)), convert(T,maximum(r)))

#==Enhance base functions
===============================================================================#
#Maximum/minimum finite values:
minfinite(::Type{T}) where T<:Integer = typemin(T)
maxfinite(::Type{T}) where T<:Integer = typemax(T)
minfinite(::Type{T}) where T<:AbstractFloat = nextfloat(convert(T, -Inf))
maxfinite(::Type{T}) where T<:AbstractFloat = prevfloat(convert(T, Inf))

#Use Limits1D with clamp:
Base.clamp(v::Union{Number,Vector}, r::Limits1D) = clamp(v, r.min, r.max)
Base.clamp!(v::Union{Number,Vector}, r::Limits1D) = clamp!(v, r.min, r.max)

#sanitize: clamp down infinite/large values & replace NaNs
#-------------------------------------------------------------------------------
_sanitize_type(::T, ::TL, ::TH) where {T, TL, TH} = promote_type_nonvoid(T, TL, TH)
_sanitize_dfltnan(::Type{T}) where T<:Integer = 0 #Does not matter (input never NaN)
_sanitize_dfltnan(::Type{T}) where T<:AbstractFloat = convert(T, NaN) #Don't modify NaN
function _sanitize_defaults(::Type{T}, min, max, nan) where T
	return tuple(convert(T, substvoid(min, minfinite(T))),
		convert(T, substvoid(max, maxfinite(T))),
		convert(T, substvoid(nan, _sanitize_dfltnan(T)))
	)
end

function _sanitize(x, min, max, nan)
	return ifelse(isnan(x), nan,
		ifelse(x>max, max,
		ifelse(x<min, min, x)))
end
function _sanitize(::Type{T}, x::Real, min, max, nan) where T
	(min, max, nan) = _sanitize_defaults(T, min, max, nan)
	return _sanitize(convert(T, x), min, max, nan)
end
sanitize(x::Real; min=nothing, max=nothing, nan=nothing) =
	_sanitize(_sanitize_type(x,min,max), x, min, max, nan)

function _sanitize(::Type{T}, x::Vector, min, max, nan) where T
	(min, max, nan) = _sanitize_defaults(T, min, max, nan)
	return T[_sanitize(xi, min, max, nan) for xi in x]
end
sanitize(x::Vector; min=nothing, max=nothing, nan=nothing) =
	_sanitize(promote_type_nonvoid(eltype(x),typeof(min),typeof(max)), x, min, max, nan)

#==Interpolation
===============================================================================#
#Interpolate between two points.
function _interpolate(p1::Point2D, p2::Point2D, x::Number)
	m = (p2.y-p1.y) / (p2.x-p1.x)
	return m*(x-p1.x)+p1.y
end
_interpolate(p1::Point2D, p2::Point2D, x::Nothing) =
	throw(ArugmentError("Must provide an x-value: interpolate(p1, p2, x=VAL)"))
interpolate(p1::Point2D, p2::Point2D; x=nothing) = _interpolate(p1, p2, x)

#Last line


text:
using SimpleSolvers
using Test
using GeometricIntegrators
using GeometricIntegrators.Integrators.VPRK
using GeometricIntegrators.Utils
using GeometricProblems.HarmonicOscillator

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

const Δt = 0.1
const nt = 1000

const ref = [-0.012420428712136283, -0.35344429048339926]

@testset "$(rpad("Harmonic Oscillator",80))" begin

    ode  = harmonic_oscillator_ode()
    pode = harmonic_oscillator_pode()
    iode = harmonic_oscillator_iode()
    dae  = harmonic_oscillator_dae()
    idae = harmonic_oscillator_idae()
    pdae = harmonic_oscillator_pdae()

    ode_equs = get_functions(ode)
    @test_nowarn ode_equs[:v](ode.t₀, ode.q₀[begin], zero(ode.q₀[begin]))

    iode_equs = get_functions(iode)
    @test_nowarn iode_equs[:ϑ](iode.t₀, iode.q₀[begin], iode.p₀[begin], zero(iode.q₀[begin]))
    @test_nowarn iode_equs[:f](iode.t₀, iode.q₀[begin], iode.p₀[begin], zero(iode.p₀[begin]))
    @test_nowarn iode_equs[:v̄](iode.t₀, iode.q₀[begin], zero(iode.q₀[begin]))
    @test_nowarn iode_equs[:f̄](iode.t₀, iode.q₀[begin], iode.p₀[begin], zero(iode.q₀[begin]))

    pode_equs = get_functions(pode)
    @test_nowarn pode_equs[:v](pode.t₀, pode.q₀[begin], pode.p₀[begin], zero(pode.q₀[begin]))
    @test_nowarn pode_equs[:f](pode.t₀, pode.q₀[begin], pode.p₀[begin], zero(pode.p₀[begin]))

    dae_equs = get_functions(dae)
    @test_nowarn dae_equs[:v](dae.t₀, dae.q₀[begin], zero(dae.q₀[begin]))
    @test_nowarn dae_equs[:u](dae.t₀, dae.q₀[begin], dae.λ₀[begin], zero(dae.q₀[begin]))
    @test_nowarn dae_equs[:ϕ](dae.t₀, dae.q₀[begin], zero(dae.λ₀[begin]))
    @test_nowarn dae_equs[:v̄](dae.t₀, dae.q₀[begin], zero(dae.q₀[begin]))

    idae_equs = get_functions(idae)
    @test_nowarn idae_equs[:ϑ](idae.t₀, idae.q₀[begin], idae.λ₀[begin], zero(idae.q₀[begin]))
    @test_nowarn idae_equs[:f](idae.t₀, idae.q₀[begin], idae.λ₀[begin], zero(idae.p₀[begin]))
    @test_nowarn idae_equs[:u](idae.t₀, idae.q₀[begin], pdae.p₀[begin], idae.λ₀[begin], zero(idae.q₀[begin]))
    @test_nowarn idae_equs[:g](idae.t₀, idae.q₀[begin], pdae.p₀[begin], idae.λ₀[begin], zero(idae.p₀[begin]))
    @test_nowarn idae_equs[:ϕ](idae.t₀, idae.q₀[begin], idae.p₀[begin], zero(idae.λ₀[begin]))
    @test_nowarn idae_equs[:v̄](idae.t₀, idae.q₀[begin], zero(idae.q₀[begin]))
    @test_nowarn idae_equs[:f̄](idae.t₀, idae.q₀[begin], idae.p₀[begin], zero(idae.q₀[begin]))

    pdae_equs = get_functions(pdae)
    @test_nowarn pdae_equs[:v](pdae.t₀, pdae.q₀[begin], pdae.p₀[begin], zero(pdae.q₀[begin]))
    @test_nowarn pdae_equs[:f](pdae.t₀, pdae.q₀[begin], pdae.p₀[begin], zero(pdae.p₀[begin]))
    @test_nowarn pdae_equs[:u](pdae.t₀, pdae.q₀[begin], pdae.p₀[begin], pdae.λ₀[begin], zero(pdae.q₀[begin]))
    @test_nowarn pdae_equs[:g](pdae.t₀, pdae.q₀[begin], pdae.p₀[begin], pdae.λ₀[begin], zero(pdae.p₀[begin]))
    @test_nowarn pdae_equs[:ϕ](pdae.t₀, pdae.q₀[begin], pdae.p₀[begin], zero(pdae.λ₀[begin]))
    @test_nowarn pdae_equs[:v̄](pdae.t₀, pdae.q₀[begin], pdae.p₀[begin], zero(pdae.q₀[begin]))
    @test_nowarn pdae_equs[:f̄](pdae.t₀, pdae.q₀[begin], pdae.p₀[begin], zero(pdae.p₀[begin]))

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

    int = IntegratorVPRKpMidpoint(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, ref) < 1E-4

    int = IntegratorVPRKpSymmetric(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, ref) < 1E-4

end


text:
function imview(colortype::Type{<:ColorTypes.Colorant}, imgvalues...)
    return PermutedDimsArray(
        ImageCore.colorview(
            colortype,
            (
                ImageCore.normedview(
                    ImageCore.Normed{eltype(img),8 * sizeof(eltype(img))},
                    img,
                ) for img in imgvalues
            )...,
        ),
        (2, 1),
    )
end

function imview(
    gci::GDALColorInterp,
    imgvalues::AbstractMatrix;
    colortable::ColorTable = ColorTable(C_NULL),
)
    return if gci == GCI_GrayIndex
        imview(ColorTypes.Gray, imgvalues)
    elseif gci == GCI_Undefined
        imview(ColorTypes.Gray, imgvalues)
    elseif gci == GCI_RedBand
        zerovalues = zeros(eltype(imgvalues), size(imgvalues))
        imview(GPI_RGB, imgvalues, zerovalues, zerovalues)
    elseif gci == GCI_GreenBand
        zerovalues = zeros(eltype(imgvalues), size(imgvalues))
        imview(GPI_RGB, zerovalues, imgvalues, zerovalues)
    elseif gci == GCI_BlueBand
        zerovalues = zeros(eltype(imgvalues), size(imgvalues))
        imview(GPI_RGB, zerovalues, zerovalues, imgvalues)
    elseif gci == GCI_PaletteIndex
        if colortable.ptr == C_NULL
            error(
                """
                `imview` is only supported for `GCI_PaletteIndex` with non-null
                colortables.
                """,
            )
        end
        gpi = paletteinterp(colortable)
        if gpi == GPI_Gray
            imview(GPI_Gray, imgvalues)
        elseif gpi == GPI_RGB
            colorentries = GDAL.GDALColorEntry[
                getcolorentryasrgb(colortable, i - 1) for
                i in 1:GDAL.gdalgetcolorentrycount(colortable.ptr)
            ]
            c1 = Matrix{UInt8}(undef, size(imgvalues)...)
            c2 = Matrix{UInt8}(undef, size(imgvalues)...)
            c3 = Matrix{UInt8}(undef, size(imgvalues)...)
            c4 = Matrix{UInt8}(undef, size(imgvalues)...)
            for i in eachindex(imgvalues)
                c1[i] = UInt8(colorentries[imgvalues[i]+1].c1)
                c2[i] = UInt8(colorentries[imgvalues[i]+1].c2)
                c3[i] = UInt8(colorentries[imgvalues[i]+1].c3)
                c4[i] = UInt8(colorentries[imgvalues[i]+1].c4)
            end
            imview(GPI_RGB, c1, c2, c3, c4)
        else
            error("""
                  Unsupported GPI: $(paletteinterp(colortable)). Please file an
                  issue at https://github.com/yeesian/ArchGDAL.jl/issues if it
                  should be supported.
                  """)
        end
    else
        error("""
              Unknown GCI: $gci. Please file an issue at
              https://github.com/yeesian/ArchGDAL.jl/issues if it should be
              supported.
              """)
    end
end

function imview(gpi::GDALPaletteInterp, imgvalues::AbstractMatrix)
    return if gpi == GPI_Gray
        imview(ColorTypes.Gray, imgvalues)
    else
        error("""
              Unsupported GPI: $gpi. Please file an issue at
              https://github.com/yeesian/ArchGDAL.jl/issues if it should be
              supported.
              """)
    end
end

function imview(
    gpi::GDALPaletteInterp,
    c1::AbstractMatrix,
    c2::AbstractMatrix,
    c3::AbstractMatrix,
)
    return if gpi == GPI_Gray
        imview(ColorTypes.Gray, c1, c2, c3)
    elseif gpi == GPI_RGB
        imview(ColorTypes.RGB, c1, c2, c3)
    else
        error("""
              Unsupported GPI: $gpi. If it should be supported, please file an
              issue at https://github.com/yeesian/ArchGDAL.jl/issues with the
              desired output.
              """)
    end
end

function imview(
    gpi::GDALPaletteInterp,
    c1::AbstractMatrix,
    c2::AbstractMatrix,
    c3::AbstractMatrix,
    c4::AbstractMatrix,
)
    return if gpi == GPI_Gray
        imview(ColorTypes.Gray, c1, c2, c3, c4)
    elseif gpi == GPI_RGB
        imview(ColorTypes.RGBA, c1, c2, c3, c4)
    else
        error("""
              Unsupported GPI: $gpi. If it should be supported, please file an
              issue at https://github.com/yeesian/ArchGDAL.jl/issues with the
              desired output.
              """)
    end
end

function imread(
    colortype::Union{GDALPaletteInterp,GDALColorInterp},
    dataset::AbstractDataset,
    i::Integer,
    args...,
)
    return imread(colortype, getband(dataset, i), args...)
end

function imread(gpi::GDALPaletteInterp, rb::AbstractRasterBand, args...)
    return imview(gpi, read(rb, args...))
end

function imread(gci::GDALColorInterp, rb::AbstractRasterBand, args...)
    return getcolortable(rb) do colortable
        return imview(gci, read(rb, args...), colortable = colortable)
    end
end

function imread(rb::AbstractRasterBand, args...)
    return getcolortable(rb) do colortable
        return imview(
            getcolorinterp(rb),
            read(rb, args...),
            colortable = colortable,
        )
    end
end

function imread(dataset::AbstractDataset, i::Integer, args...)
    return imread(getband(dataset, i), args...)
end

function _colorindices(dataset::AbstractDataset, indices)
    gci = unique(
        GDALColorInterp[getcolorinterp(getband(dataset, i)) for i in indices],
    )
    gciorder = sort(gci)
    colortype = if gciorder == [GCI_GrayIndex]
        GCI_GrayIndex
    elseif gciorder == [GCI_PaletteIndex]
        GCI_PaletteIndex
    elseif gciorder == [GCI_RedBand, GCI_GreenBand, GCI_BlueBand]
        GPI_RGB
    elseif gciorder == [GCI_RedBand, GCI_GreenBand, GCI_BlueBand, GCI_AlphaBand]
        GPI_RGB
    else
        error("""
              Unknown GCI: $gciorder. Please file an issue at
              https://github.com/yeesian/ArchGDAL.jl/issues if it should be
              supported.
              """)
    end
    return colortype, Tuple(indices[sortperm(gci)])
end

function imread(dataset::AbstractDataset, indices, args...)
    colortype, idxs = _colorindices(dataset, indices)
    return if colortype == GCI_PaletteIndex
        getcolortable(getband(dataset, 1)) do colortable
            return imview(
                colortype,
                (read(getband(dataset, i), args...) for i in idxs)...,
                colortable = colortable,
            )
        end
    else
        imview(colortype, (read(getband(dataset, i), args...) for i in idxs)...)
    end
end

imread(dataset::AbstractDataset) = imread(dataset, 1:nraster(dataset))

function imread(filename::AbstractString)
    return read(filename) do dataset
        return imread(dataset)
    end
end


text:
#
# cdefs.jl --
#
# Constants and structures defined in "xpa.h".  These "private" definitions are
# put in a sub-module to make them not directly accessible.
#
module CDefs

export XPACommRec, XPARec

#
# CONSTANTS
#

# Sizes.
const SZ_LINE = 4096
const XPA_NAMELEN = 1024

# This is the number of actual commands we have above
const XPA_CMDS = 4

#
# CALLBACK FUNCTIONS
#

# typedef int (*SendCb)(void *client_data, void *call_data, char *paramlist,
#                       char **buf, size_t *len);
struct SendCb end

# typedef int (*ReceiveCb)(void *client_data, void *call_data,
#                          char *paramlist, char *buf, size_t len );
struct ReceiveCb end

# typedef int (*InfoCb)(void *client_data, void *call_data,
#                       char *paramlist);
struct InfoCb end

# typedef void *(*SelAdd)(void *client_data, int fd);
abstract type SelAdd end

# typedef void (*SelDel)(void *client_data);
abstract type SelDel end

# typedef void (*SelOn)(void *client_data);
abstract type SelOn end

# typedef void  (*SelOff)(void *client_data);
abstract type SelOff end

# typedef void (*MyFree)(void *buf);
struct MyFree end

#
# OPAQUE STRUCTURES
#

# Opaque structure for `struct nsrec`.
abstract type NSRec end

# Opaque structure for `struct cliprec`.
abstract type ClipRec end

# Opaque structure for `struct xpainputrec`.
abstract type XPAInputRec end

# Opaque structure for `struct xpaclientrec`.
abstract type XPAClientRec end

# Opaque structure for `struct xpacmdrec`.
abstract type XPACmdRec end

#
# EXPORTED STRUCTURES
#

# XPA communication structure for each connection.
struct XPACommRec
    next::Ptr{XPACommRec}
    status::Cint
    message::Cint
    n::Cint
    cmd::Cint
    mode::Cint
    telnet::Cint
    usebuf::Cint
    useacl::Cint
    id::Ptr{UInt8} # (C string)
    info::Ptr{UInt8} # (C string)
    target::Ptr{UInt8} # (C string)
    paramlist::Ptr{UInt8} # (C string)
    cmdfd::Cint
    datafd::Cint
    cendian::Ptr{UInt8} # (C string)
    ack::Cint
    # buf and len passed to callbacks
    buf::Ptr{UInt8} # (C string)
    len::Csize_t
    # for AF_INET
    cmdip::Cuint
    cmdport::Cint
    dataport::Cint
    # for AF_UNIX
    cmdname::Ptr{UInt8} # (C string)
    dataname::Ptr{UInt8} # (C string)
    acl::NTuple{XPA_CMDS+1,Cint}
    # for handling fd's in non-select event loops
    selcptr::Ptr{Cvoid} # cmdfd struct for seldel
    seldptr::Ptr{Cvoid} # datafd struct for seldel
    # pointer to associated name server
    ns::Ptr{NSRec}
    # myfree routine
    myfree::Ptr{MyFree}
    myfree_ptr::Ptr{Cvoid}
end

struct XPARec
    version::Ptr{UInt8} # xpa version string
    status::Cint        # status of this xpa
    type::Ptr{UInt8}    # (C string) "g", "s", "i" are server types; "c" for client
    #
    # THE SERVER SIDE
    #
    next::Ptr{XPARec}
    xclass::Ptr{UInt8} # (C string)
    name::Ptr{UInt8}   # (C string)
    help::Ptr{UInt8}   # (C string)
    # send callback info
    send_callback::Ptr{SendCb}
    send_data::Ptr{Cvoid}
    send_mode::Cint
    # receive callback info
    receive_callback::Ptr{ReceiveCb}
    receive_data::Ptr{Cvoid}
    receive_mode::Cint
    # info callback info
    info_callback::Ptr{InfoCb}
    info_data::Ptr{Cvoid}
    info_mode::Cint
    # list of sub-commands for this access point
    commands::Ptr{XPACmdRec}
    # communication info
    fd::Cint                      # listening socket file descriptor
    method::Ptr{UInt8}            # (C string) method string: host:ip or unix_filename
    nshead::Ptr{NSRec}            # name servers associated with this access point
    commhead::Ptr{XPACommRec}     # linked list of communcation records
    cliphead::Ptr{ClipRec}        # linked list of cliboard records
    filename::Ptr{UInt8}          # (C string) file name (unix sockets) for listening
    sendian::Ptr{UInt8}           # (C string) endian-ness of server
    # request-specific info
    comm::Ptr{XPACommRec}         # current comm if we are processing a request
    # select loop info
    seldel::Ptr{SelDel}           # routine to remove xpa socket from select loop
    seladd::Ptr{SelAdd}           # routine to add xpa command sockets to select loop
    selon ::Ptr{SelOn}            # routine to enable xpa command sockets
    seloff::Ptr{SelOff}           # routine to disable xpa command sockets
    selptr::Ptr{Cvoid}            # additional info for seldelete()
    #
    # THE CLIENT SIDE
    #
    persist::Cint                 # flag whether this is a persistent client
    nclient::Cint                 # number of clients -- used in processing headers
    client_mode::Cint             # global client mode
    clienthead::Ptr{XPAClientRec} # linked list of active clients
    ifd::Cint                     # input fd for XPASetFd()
    inpbytes::Csize_t             # total number of bytes in input lists
    inphead::Ptr{XPAInputRec}     # linked list of input structs
end

end # module


text:
using Test
using Distributed

addprocs(2)
using IncrementalInference
@everywhere using IncrementalInference

@testset "test multiprocess solveTree!" begin

fg = generateGraph_Kaess()

solveTree!(fg)

end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

function GitReference(repo::GitRepo, refname::AbstractString)
    ref_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_reference_lookup, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Cstring),
                   ref_ptr_ptr, repo.ptr, refname)
    return GitReference(ref_ptr_ptr[])
end

function GitReference(repo::GitRepo, obj_oid::Oid, refname::AbstractString = Consts.HEAD_FILE;
                      force::Bool=false, msg::AbstractString="")
    ref_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_reference_create, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{UInt8}, Ptr{Oid}, Cint, Cstring),
                   ref_ptr_ptr, repo.ptr, refname, Ref(obj_oid), Cint(force),
                   isempty(msg) ? Cstring_NULL : msg)
    return GitReference(ref_ptr_ptr[])
end

function head(repo::GitRepo)
    head_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_repository_head, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}), head_ptr_ptr, repo.ptr)
    return GitReference(head_ptr_ptr[])
end

function shortname(ref::GitReference)
    isempty(ref) && return ""
    name_ptr = ccall((:git_reference_shorthand, :libgit2), Cstring, (Ptr{Void},), ref.ptr)
    name_ptr == C_NULL && return ""
    return unsafe_string(name_ptr)
end

function reftype(ref::GitReference)
    return ccall((:git_reference_type, :libgit2), Cint, (Ptr{Void},), ref.ptr)
end

function fullname(ref::GitReference)
    isempty(ref) && return ""
    reftype(ref) == Consts.REF_OID && return ""
    rname = ccall((:git_reference_symbolic_target, :libgit2), Cstring, (Ptr{Void},), ref.ptr)
    rname == C_NULL && return ""
    return unsafe_string(rname)
end

function name(ref::GitReference)
    isempty(ref) && return ""
    name_ptr = ccall((:git_reference_name, :libgit2), Cstring, (Ptr{Void},), ref.ptr)
    name_ptr == C_NULL && return ""
    return unsafe_string(name_ptr)
end

function branch(ref::GitReference)
    isempty(ref) && return ""
    str_ptr_ptr = Ref(LibGit2.Cstring_NULL)
    @check ccall((:git_branch_name, :libgit2), Cint,
                  (Ptr{Cstring}, Ptr{Void},), str_ptr_ptr, ref.ptr)
    return unsafe_string(str_ptr_ptr[])
end

function ishead(ref::GitReference)
    isempty(ref) && return false
    err = ccall((:git_branch_is_head, :libgit2), Cint,
                  (Ptr{Void},), ref.ptr)
    return err == 1
end

function isbranch(ref::GitReference)
    isempty(ref) && return false
    err = ccall((:git_reference_is_branch, :libgit2), Cint,
                  (Ptr{Void},), ref.ptr)
    return err == 1
end

function istag(ref::GitReference)
    isempty(ref) && return false
    err = ccall((:git_reference_is_tag, :libgit2), Cint,
                  (Ptr{Void},), ref.ptr)
    return err == 1
end

function isremote(ref::GitReference)
    isempty(ref) && return false
    err = ccall((:git_reference_is_remote, :libgit2), Cint,
                  (Ptr{Void},), ref.ptr)
    return err == 1
end

function peel{T <: GitObject}(::Type{T}, ref::GitReference)
    git_otype = getobjecttype(T)
    obj_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    err = ccall((:git_reference_peel, :libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Cint), obj_ptr_ptr, ref.ptr, git_otype)
    if err == Int(Error.ENOTFOUND)
        return Oid()
    elseif err != Int(Error.GIT_OK)
        if obj_ptr_ptr[] != C_NULL
            finalize(GitAnyObject(obj_ptr_ptr[]))
        end
        throw(Error.GitError(err))
    end
    return T(obj_ptr_ptr[])
end

function ref_list(repo::GitRepo)
    with(StrArrayStruct()) do sa
        sa_ref = Ref(sa)
        @check ccall((:git_reference_list, :libgit2), Cint,
                      (Ptr{StrArrayStruct}, Ptr{Void}), sa_ref, repo.ptr)
        convert(Vector{AbstractString}, sa_ref[])
    end
end

function create_branch(repo::GitRepo,
                       bname::AbstractString,
                       commit_obj::GitCommit;
                       force::Bool=false)
    ref_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_branch_create, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Cstring, Ptr{Void}, Cint),
                   ref_ptr_ptr, repo.ptr, bname, commit_obj.ptr, Cint(force))
    return GitReference(ref_ptr_ptr[])
end

function delete_branch(branch::GitReference)
    @check ccall((:git_branch_delete, :libgit2), Cint, (Ptr{Void},), branch.ptr)
end

function head!(repo::GitRepo, ref::GitReference)
    ref_name = name(ref)
    @check ccall((:git_repository_set_head, :libgit2), Cint,
                  (Ptr{Void}, Cstring), repo.ptr, ref_name)
    return ref
end

function lookup_branch(repo::GitRepo,
                       branch_name::AbstractString,
                       remote::Bool=false)
    ref_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    branch_type = remote ? Consts.BRANCH_REMOTE : Consts.BRANCH_LOCAL
    err = ccall((:git_branch_lookup, :libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{UInt8}, Cint),
                  ref_ptr_ptr, repo.ptr, branch_name, branch_type)
    if err == Int(Error.ENOTFOUND)
        return nothing
    elseif err != Int(Error.GIT_OK)
        if ref_ptr_ptr[] != C_NULL
            finalize(GitReference(ref_ptr_ptr[]))
        end
        throw(Error.GitError(err))
    end
    return GitReference(ref_ptr_ptr[])
end

function upstream(ref::GitReference)
    isempty(ref) && return nothing
    ref_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    err = ccall((:git_branch_upstream, :libgit2), Cint,
                  (Ref{Ptr{Void}}, Ptr{Void},), ref_ptr_ptr, ref.ptr)
    if err == Int(Error.ENOTFOUND)
        return nothing
    elseif err != Int(Error.GIT_OK)
        if ref_ptr_ptr[] != C_NULL
            finalize(GitReference(ref_ptr_ptr[]))
        end
        throw(Error.GitError(err))
    end
    return GitReference(ref_ptr_ptr[])
end

function owner(ref::GitReference)
    repo_ptr = ccall((:git_reference_owner, :libgit2), Ptr{Void},
                      (Ptr{Void},), ref.ptr)
    return GitRepo(repo_ptr)
end

function target!(ref::GitReference, new_oid::Oid; msg::AbstractString="")
    ref_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_reference_set_target, :libgit2), Cint,
             (Ptr{Ptr{Void}}, Ptr{Void}, Ptr{Oid}, Cstring),
             ref_ptr_ptr, ref.ptr, Ref(new_oid), isempty(msg) ? Cstring_NULL : msg)
    return GitReference(ref_ptr_ptr[])
end

function GitBranchIter(r::GitRepo, flags::Cint=Cint(Consts.BRANCH_LOCAL))
    bi_ptr = Ref{Ptr{Void}}(C_NULL)
    @check ccall((:git_branch_iterator_new, :libgit2), Cint,
                  (Ptr{Ptr{Void}}, Ptr{Void}, Cint), bi_ptr, r.ptr, flags)
    return GitBranchIter(bi_ptr[])
end

function Base.start(bi::GitBranchIter)
    ref_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    btype = Ref{Cint}()
    err = ccall((:git_branch_next, :libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Cint}, Ptr{Void}),
                  ref_ptr_ptr, btype, bi.ptr)
    err != Int(Error.GIT_OK) && return (nothing, -1, true)
    return (GitReference(ref_ptr_ptr[]), btype[], false)
end

Base.done(bi::GitBranchIter, state) = Bool(state[3])

function Base.next(bi::GitBranchIter, state)
    ref_ptr_ptr = Ref{Ptr{Void}}(C_NULL)
    btype = Ref{Cint}()
    err = ccall((:git_branch_next, :libgit2), Cint,
                 (Ptr{Ptr{Void}}, Ptr{Cint}, Ptr{Void}),
                  ref_ptr_ptr, btype, bi.ptr)
    err != Int(Error.GIT_OK) && return (state[1:2], (nothing, -1, true))
    return (state[1:2], (GitReference(ref_ptr_ptr[]), btype[], false))
end

Base.iteratorsize(::Type{GitBranchIter}) = Base.SizeUnknown()

function Base.map(f::Function, bi::GitBranchIter)
    res = nothing
    s = start(bi)
    while !done(bi, s)
        val = f(s[1:2])
        if res === nothing
            res = Array{typeof(val)}(0)
        end
        push!(res, val)
        val, s = next(bi, s)
    end
    return res
end

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
@testset "test contigency" begin

@testset "pg response - $(i)" for (i,network) in enumerate(networks)

    network = deepcopy(network)
    PowerModels.update_data!(network, solutions[i])

    for (i,gen) in network["gen"]
        gen["pg_base"] = gen["pg"]
    end

    pg_target = 0.0
    for (i,gen) in network["gen"]
        if gen["gen_status"] != 0
            pg_target += gen["pg"]
        end
    end

    for (i,gen) in network["gen"]
        if gen["gen_status"] == 0
            continue
        end

        gen_bus = network["bus"]["$(gen["gen_bus"])"]
        area_gens = network["area_gens"][gen_bus["area"]]

        # NOTE this is a slow operation
        network_cont = deepcopy(network)

        cont_gen = network_cont["gen"][i]
        cont_gen["gen_status"] = 0
        pg_missing = cont_gen["pg"]

        #println("")
        #println("generation removed: $(pg_missing)")

        apply_pg_response!(network_cont, area_gens, pg_missing)

        pg_total = 0.0
        for (i,gen) in network_cont["gen"]
            if gen["gen_status"] != 0
                pg_total += gen["pg"]
            end
        end

        pg_comp = comp_pg_response_total(network_cont, area_gens)

        #println("delta $delta")
        #println("gen $(pg_target) / $(pg_total) / $(pg_comp)")

        @test isapprox(pg_target, pg_total)
        @test isapprox(pg_total, pg_comp)
    end

end

cuts_ratec_nd_first_lazy_gen = [0, 1]
cuts_ratec_nd_first_lazy_branch = [0, 1]
@testset "cuts ratec_nd_first_lazy - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_opf_cheap(network, DCPPowerModel, lp_solver)
    @test isapprox(result["termination_status"], OPTIMAL)
    update_active_power_data!(network, result["solution"])

    cuts = check_contingencies_branch_power_bpv(network, total_cut_limit=1000, gen_flow_cuts=[], branch_flow_cuts=[])

    @test isapprox(length(cuts.gen_cuts), cuts_ratec_nd_first_lazy_gen[i])
    @test isapprox(length(cuts.branch_cuts), cuts_ratec_nd_first_lazy_branch[i])
end

cuts_ratec_gen = [0, 10]
cuts_ratec_branch = [0, 10]
@testset "cuts ratec - $(i)" for (i,network) in enumerate(networks)
    network = deepcopy(network)
    network["gen_flow_cuts"] = []
    network["branch_flow_cuts"] = []

    result = run_opf_cheap(network, DCPPowerModel, lp_solver)
    @test isapprox(result["termination_status"], OPTIMAL)
    update_active_power_data!(network, result["solution"])

    cuts = check_contingencies_branch_power(network, total_cut_limit=1000, gen_flow_cuts=[], branch_flow_cuts=[])

    @test isapprox(length(cuts.gen_cuts), cuts_ratec_gen[i])
    @test isapprox(length(cuts.branch_cuts), cuts_ratec_branch[i])
end

end # close test group


text:
using Test
using LatticeTools

@testset "sitemap" begin
    unitcell = makeunitcell([1.0 0.0; 0.0 1.0]; SiteType=String)
    addsite!(unitcell, "Ox", FractCoord([0,0], [0.5, 0.0]))
    addsite!(unitcell, "Oy", FractCoord([0,0], [0.0, 0.5]))
    lattice = make_lattice(unitcell, [4 0; 0 4])

    tsym = FiniteTranslationSymmetry(lattice)
    psym = project(PointSymmetryDatabase.get(13), [1 0 0; 0 1 0])

    @test findsitemap(unitcell, TranslationOperation([1,0])) == [(1,[0,0]), (2,[0,0])]
    @test findsitemap(unitcell, TranslationOperation([2,0])) == [(1,[0,0]), (2,[0,0])]

    @test findsitemap(unitcell, tsym) == collect(
        Iterators.repeated([(1,[0,0]), (2,[0,0])], 16)
    )

    rotC4 = element(psym, 3)
    @test findsitemap(unitcell, rotC4) == [ (2, [0,0]), (1, [-1,0]) ]
end


text:
"""
    sic(v::VAREstimate)

Return the Schwarz (Bayesian) Information Criterion for VAR model `v`.
"""
function sic(v::VAREstimate)
    obs, K = size(residuals(v))
    nparam = size(coef(v), 1)
    ΣU = v.ΣU
    ΣUml = ((obs - nparam) / obs) * ΣU
    return logdet(ΣUml) + (log(obs) / obs) * nparam * K
end

"""
    hqc(v::VAREstimate)

Return the Hannan-Quinn Criterion for VAR model `v`.
"""
function hqc(v::VAREstimate)
    obs, K = size(residuals(v))
    nparam = size(coef(v), 1)
    ΣU = v.ΣU
    ΣUml = ((obs - nparam) / obs) * ΣU
    return logdet(ΣUml) + (2 * log(log(obs)) / obs) * nparam * K
end

"""
    StabilityCheck(v::VAREstimate) -> StabilityCheck

Check the stability of the estimated VAR model `v`.

# Fields
- `isstable::Bool`
- `eigenvals::Array{Number}`
- `eigenmod::Array{Float64}`
"""
struct StabilityCheck
    isstable::Bool
    eigenvals::Array{Number}
    eigenmod::Array{Float64}
end

function StabilityCheck(v::VAREstimate)
    # See Lutkepohl (2006) p.15ff
    p = v.lags
    B = coef(v)
    nparam, K = size(B)
    startidx = nparam - K * p + 1
    A = B[startidx:end, :]
    Acomp = companion(A, p)
    eigenvals = eigvals(Acomp)
    eigenmod = abs.(eigenvals)
    isstable = all(eigenmod .< 1)
    StabilityCheck(isstable, eigenvals, eigenmod)
end

function show(io::IO, obj::StabilityCheck)
    # K = length(stab.ynames)
    E = obj.eigenvals
    Emod = round.(unique(obj.eigenmod), digits = 3)
    println(io, typeof(obj))
    if obj.isstable
        println(io, "VAR is stable")
    else
        println(io, "VAR is not stable")
    end
    println(io, "Modulus of eigenvalues (dropping repeated values):")
    println(io, sort(Emod))
end

"""
    PortmanteauTest(v::VAREstimate, h::Int) -> PortmanteauTest

Conduct a multivariate portmanteau test for residual autocorrelation up to lag
    `h` for VAR model `v`.
Implements the adjusted test described on p.171 of Lütkepohl (2006).

# Fields
- `U::Matrix`
- `h::Int`
- `Q::Float64`
- `df::Int`
- `pval::Float64`
"""
struct PortmanteauTest
    U::Matrix
    h::Int
    Q::Float64
    df::Int
    pval::Float64
end

function PortmanteauTest(v::VAREstimate, h::Int)
    U = residuals(v)
    T, K = size(U)
    p = v.lags
    C = permutedims(crosscov(U, U, collect(0:h), demean = false), [3, 2, 1])
    Q = 0.0
    for i = 1:h
        Q += tr(C[:, :, i+1]' * (C[:, :, 1] \ C[:, :, i+1] / C[:, :, 1])) * (T^2 / (T - i))
    end
    df = K^2 * (h - p)
    pval = ccdf(Chisq(df), Q)
    PortmanteauTest(U, h, Q, df, pval)
end

function show(io::IO, obj::PortmanteauTest)
    println(io, typeof(obj))
    println(io, "Multivariate test for residual autocorrelation.")
    println(io, "- Null is no autocorrelation")
    println(io, "- Lags: $(obj.h)")
    println(io, "- Q stat: $(round(obj.Q; digits=3))")
    println(io, "- p-value: $(round(obj.pval; digits=3))")
    outcome = obj.pval > 0.05 ? "Fail to reject" : "Reject"
    println(io, "$outcome the null at the 5% level")
end

"""
    LMCorrTest(v::VAREstimate, h::Int; smallsample = false) -> LMCorrTest

Conduct a multivariate LM test (Breusch-Godfrey test) for residual autocorrelation
    up to lag `h` for VAR model `v`.
If `smallsample = true`, the alternative test statistic detailed in Lütkepohl (2006),
    p.173 is used.

# Fields
- `h::Int`
- `Q::Float64`
- `df1::Any`
- `df2::Any`
- `pval::Float64`
- `smallsample::Bool`
"""
struct LMCorrTest
    h::Int
    Q::Float64
    df1::Any
    df2::Any
    pval::Float64
    smallsample::Bool
end

function LMCorrTest(v::VAREstimate, h::Int; smallsample = false)
    U = residuals(v)
    T, K = size(U)
    p = v.lags
    nparam = size(v.B, 1)
    Z = v.Z

    U2 = vcat(zeros(p + (h - p), K), U)
    Ulag = VectorAutoregressions.lag_matrix(U2, h)
    ZUlag = [Z Ulag]
    AD = (ZUlag' * ZUlag) \ (ZUlag' * U)
    E = U - ZUlag * AD

    ΣUml = ((T - nparam) / T) * v.ΣU
    ΣEml = (E' * E) / T

    if smallsample
        # Test stat from Lütkepohl 2006, p.173
        s = sqrt((K^4 * h^2 - 4) / (K^2 + (K^2 * h^2) - 5))
        N = T - (K*p) - 1 - (K*h) - (K - K*h +1)/2
        df1 = h * K^2
        df2 = (N * s) - ((K^2 * h) / 2) + 1
        Q = ((det(ΣUml) / det(ΣEml))^(1 / s) - 1) * (df2 / df1)
        pval = ccdf(FDist(df1, df2), Q)
    else
        # Test stat from Killian and Lütkepohl 2017, p.54
        Q = T * (K - tr(ΣUml \ ΣEml))
        df1 = h * K^2
        df2 = nothing
        pval = ccdf(Chisq(df1), Q)
    end
    return LMCorrTest(h, Q, df1, df2, pval, smallsample)
end

function show(io::IO, obj::LMCorrTest)
    println(io, typeof(obj))
    println(io, "Multivariate LM test for residual autocorrelation.")
    println(io, "- Null is no autocorrelation")
    println(io, "- Lags: $(obj.h)")
    println(io, "- Q stat: $(round(obj.Q; digits=3))")
    println(io, "- p-value: $(round(obj.pval; digits=3))")
    if obj.smallsample
        println(io, "- Small sample test")
    end
    outcome = obj.pval > 0.05 ? "Fail to reject" : "Reject"
    println(io, "$outcome the null at the 5% level")
end

struct ResidualNormalityTest
    λ_s::Float64
    λ_k::Float64
    λ_sk::Float64
    pval1::Float64
    pval2::Float64
    pval3::Float64
end

function ResidualNormalityTest(v::VAREstimate)
    U = residuals(v)
    ΣU = v.ΣU
    λs, λk, λsk, pval1, pval2, pval3 = multivariate_normality_test(U, ΣU)
    ResidualNormalityTest(λs, λk, λsk, pval1, pval2, pval3)
end

function multivariate_normality_test(U, Σ)
    # see Lutkepohl pp. 177-181
    uobs, K = size(U)
    P   = cholesky(Σ).L
    W   = transpose(P\U')
    W3  = W.^3
    W4  = W.^4
    b1  = dropdims(sum(W3, dims=1)/uobs, dims=1)
    b2  = dropdims(sum(W4, dims=1)/uobs, dims=1)
    λs = dot(b1, b1)*(uobs/6)
    λk = dot(b2 .- 3, b2 .- 3)*(uobs/24)
    λsk = λs + λk

    pval1 = 1 - cdf(Chisq(K), λs)
    pval2 = 1 - cdf(Chisq(K), λk)
    pval3 = 1 - cdf(Chisq(2K), λsk)

    return λs, λk, λsk, pval1, pval2, pval3
end


text:
#=--- CpuId / test / mock.jl ----------------------------------------------=#

#
# This file contains a number of cpuid dumps from known CPUs.
# These dumps are used to test the analysis function of CpuId.jl.
#

import CpuId
import CpuId.CpuInstructions

"""
Entry for mocking `cpuid`, maps input and output of the `cpuid` instruction
on a known CPU; viz. `(eax, ecx) -> (eax, ebx, ecx, edx)`.
"""
const _mockdb_entry = Pair{ Dict{NTuple{2,UInt32}, NTuple{4,UInt32}}, Dict{Symbol,Any} }

"""
List of cpuid dumps.
Each entry is represented by a dictionary which maps input and output of the
`cpuid` instruction on a known CPU; viz. `(eax, ecx) -> (eax, ebx, ecx, edx)`.
"""
const _mockdb = _mockdb_entry[]

"""
Temporarily overwrite the low-level cpuid function, to return pre-recorded results.
"""
function mock_cpuid(idx::Integer)

    global _mockdb

    !(1 ≤ idx ≤ length(_mockdb)) &&
        error("The cpuid mocking database does only have entries " *
              "with indices 1:$(length(_mockdb))")

    _fake_cpuid(eax=0, ecx=0)::NTuple{4,UInt32} =
        get( first(_mockdb[idx]), (UInt32(eax), UInt32(ecx))
           , (zero(UInt32), zero(UInt32), zero(UInt32), zero(UInt32),) )

    CpuInstructions.eval( :(cpuid(eax=0, ecx=0) = $_fake_cpuid(eax, ecx) ) )
end

"""
Dump a table of cpuid answers for mocking.  Ouput in a copy & paste ready format.
"""
function dump_cpuid_table()

    println("Dump of all cpuid leafs, used for mocking 'cpuid'...\n")

    println("# ", strip(CpuId.cpubrand()), " with",
            CpuId.hypervised() ? " "*string(CpuId.hvvendor()) : "out", " hypervisor" )
    println("push!( _mockdb, (Dict(")

    for minleaf in [0x0000_0000,0x2000_0000,0x4000_0000,0x8000_0000]
        # get the maximum leaf
        maxleaf = first(CpuId.cpuid(minleaf))
        maxleaf < minleaf && continue
        for leaf in minleaf:maxleaf
            println("    (",(leaf, 0x0000_0000), " => ", CpuId.cpuid(leaf), "),")
        end
    end

    leaf = 0x0000_0004
    if CpuId.hasleaf(leaf)
        for subleaf in 0x0000_0001:0x0000_000f
            eax, ebx, ecx, edx = CpuId.cpuid(leaf, subleaf)
            println("    (",(leaf, subleaf), " => ", (eax, ebx, ecx, edx), "),")
            eax & 0x1f == 0 && break
        end
    end

    leaf = 0x0000_000b
    if CpuId.hasleaf(leaf)
        for subleaf in 0x0000_0001:0x0000_000f
            eax, ebx, ecx, edx = CpuId.cpuid(leaf, subleaf)
            println("    (",(leaf, subleaf), " => ", (eax, ebx, ecx, edx), "),")
            ebx & 0xffff == 0x0000 && break
        end
    end

    leaf = 0x8000_001d
    if CpuId.hasleaf(leaf)
        for subleaf in 0x0000_0001:0x0000_000f
            eax, ebx, ecx, edx = CpuId.cpuid(leaf, subleaf)
            println("    (",(leaf, subleaf), " => ", (eax, ebx, ecx, edx), "),")
            eax & 0x1f == 0 && break
        end
    end

    println("  ) => Dict{Symbol,Any}(")
    # print the results of certain identification functions
    println("    :cpuvendor       => :", cpuvendor(),       ",")
    println("    :cpuarchitecture => :", cpuarchitecture(), ",")
    println("    :cpucores        => ", cpucores(),         ",")
    println("    :cputhreads      => ", cputhreads(),       ",")
    println("    :cachesize       => ", cachesize(),        ",")
    println("    :cachelinesize   => ", cachelinesize(),    ",")
    println("    :simdbits        => ", simdbits(),         ",")

    println("  )))\n\nDone.\n")

end

#=--- end of file ---------------------------------------------------------=#