text:
abstract type NBParam end

abstract type NBParamProb <: NBParam end
default_invlink(::NBParamProb) = logistic
abstract type NBParamMean <: NBParam end
default_invlink(::NBParamMean) = exp

"""
    NegativeBinomialLikelihood(param::NBParam, invlink::Union{Function,Link})

There are many possible parametrizations for the Negative Binomial likelihood.
We follow the convention laid out in p.137 of [^Hilbe'11] and provide some common parametrizations.
The `NegativeBinomialLikelihood` has a special structure; the first type parameter `NBParam`
defines in what parametrization the latent function is used, and contains the other (scalar) parameters.
`NBParam` itself has two subtypes:
- `NBParamProb` for parametrizations where `f -> p`, the probability of success of a Bernoulli event
- `NBParamMean` for parametrizations where `f -> μ`, the expected number of events

## `NBParam` predefined types

### `NBParamProb` types with `p = invlink(f)` the probability of success
- [`NBParamSuccess`](@ref): This is the definition used in [`Distributions.jl`](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial).
- [`NBParamFailure`](@ref): This is the definition used in [Wikipedia](https://en.wikipedia.org/wiki/Negative_binomial_distribution).

### `NBParamMean` types with `μ = invlink(f)` the mean/expected number of events
- [`NBParamI`](@ref): Mean is linked to `f` and variance is given by `μ(1 + α)`
- [`NBParamII`](@ref): Mean is linked to `f` and variance is given by `μ(1 + α * μ)`
- [`NBParamPower`](@ref): Mean is linked to `f` and variance is given by `μ(1 + α * μ^ρ)`

To create a new parametrization, you need to:
- create a new type `struct MyNBParam{T} <: NBParam; myparams::T; end`;
- dispatch `(l::NegativeBinomialLikelihood{<:MyNBParam})(f::Real)`, which must return a [`NegativeBinomial`](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial) from `Distributions.jl`.
`NegativeBinomial` follows the parametrization of [`NBParamSuccess`](@ref), i.e. the first argument is the number of successes
and the second argument is the probability of success.

## Examples

```julia-repl
julia> NegativeBinomialLikelihood(NBParamSuccess(10), logistic)(2.0)
NegativeBinomial{Float64}(r=10.0, p=0.8807970779778824)
julia> NegativeBinomialLikelihood(NBParamFailure(10), logistic)(2.0)
NegativeBinomial{Float64}(r=10.0, p=0.11920292202211757)

julia> d = NegativeBinomialLikelihood(NBParamI(3.0), exp)(2.0)
NegativeBinomial{Float64}(r=2.4630186996435506, p=0.25)
julia> mean(d) ≈ exp(2.0)
true
julia> var(d) ≈ exp(2.0) * (1 + 3.0)
true
```

[^Hilbe'11]: Hilbe, Joseph M. Negative binomial regression. Cambridge University Press, 2011.
"""
struct NegativeBinomialLikelihood{Tp<:NBParam,Tl<:AbstractLink} <: AbstractLikelihood
    params::Tp # Likelihood parametrization (and parameters)
    invlink::Tl
    function NegativeBinomialLikelihood(
        params::Tparam, invlink=default_invlink(params)
    ) where {Tparam<:NBParam}
        # we convert `invlink` into a `Link` object if it's not the case already
        invlink = link(invlink)
        return new{Tparam,typeof(invlink)}(params, invlink)
    end
end

@deprecate NegativeBinomialLikelihood(l=logistic; successes=1) NegativeBinomialLikelihood(
    NBParamI(successes), l
)

function (l::NegativeBinomialLikelihood)(::Real)
    return error(
        "not implemented for type $(typeof(l)). For your custom type to run you ",
        "need to implement `(l::NegativeBinomialLikelihood{<:MyNBParam})(f::Real)`. ",
        "For a full explanation, see `NegativeBinomialLikelihood` docs",
    )
end

@functor NegativeBinomialLikelihood

(l::NegativeBinomialLikelihood)(fs::AbstractVector{<:Real}) = Product(map(l, fs))

@doc raw"""
    NBParamSuccess(successes)

Negative Binomial parametrization with `successes` the number of successes and
`invlink(f)` the probability of success.
This corresponds to the definition used by [Distributions.jl](https://juliastats.org/Distributions.jl/latest/univariate/#Distributions.NegativeBinomial).

```math
  p(y|\text{successes}, p=\text{invlink}(f)) = \frac{\Gamma(y+\text{successes})}{y! \Gamma(\text{successes})} p^\text{successes} (1 - p)^y
```
"""
struct NBParamSuccess{T} <: NBParamProb
    successes::T
end

function (l::NegativeBinomialLikelihood{<:NBParamSuccess})(f::Real)
    return NegativeBinomial(l.params.successes, l.invlink(f))
end

@doc raw"""
    NBParamFailure(failures)

Negative Binomial parametrization with `failures` the number of failures and
`invlink(f)` the probability of success.
This corresponds to the definition used by [Wikipedia](https://en.wikipedia.org/wiki/Negative_binomial_distribution).

```math
  p(y|\text{failures}, p=\text{invlink}(f)) = \frac{\Gamma(y+\text{failures})}{y! \Gamma(\text{failures})} p^y (1 - p)^{\text{failures}}
```
"""
struct NBParamFailure{T} <: NBParamProb
    failures::T
end

function (l::NegativeBinomialLikelihood{<:NBParamFailure})(f::Real)
    return NegativeBinomial(l.params.failures, 1 - l.invlink(f))
end

# Helper function to convert mean and variance to p and r
_nb_mean_excessvar_to_r_p(μ::Real, ev::Real) = μ / ev, 1 / (1 + ev)

"""
    NBParamI(α)

Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α)` where `α > 0`.
"""
struct NBParamI{T} <: NBParamMean
    α::T
end

function (l::NegativeBinomialLikelihood{<:NBParamI})(f::Real)
    μ = l.invlink(f)
    ev = l.params.α
    return NegativeBinomial(_nb_mean_excessvar_to_r_p(μ, ev)...)
end

"""
    NBParamII(α)

Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α * μ)` where `α > 0`.
"""
struct NBParamII{T} <: NBParamMean
    α::T
end

function (l::NegativeBinomialLikelihood{<:NBParamII})(f::Real)
    μ = l.invlink(f)
    ev = l.params.α * μ
    return NegativeBinomial(_nb_mean_excessvar_to_r_p(μ, ev)...)
end

"""
    NBParamPower(α, ρ)

Negative Binomial parametrization with mean `μ = invlink(f)` and variance `v = μ(1 + α * μ^ρ)` where `α > 0` and `ρ > 0`.
"""
struct NBParamPower{Tα,Tρ} <: NBParamMean
    α::Tα
    ρ::Tρ
end

function (l::NegativeBinomialLikelihood{<:NBParamPower})(f::Real)
    μ = l.invlink(f)
    ev = l.params.α * μ^l.params.ρ
    return NegativeBinomial(_nb_mean_excessvar_to_r_p(μ, ev)...)
end


text:
using Random

include("utils.jl")

@syms 𝑥
@polyvar x

function test_deriv(x; n=10, min_deg=1, max_deg=6, sparcity=0.5)
    k = 0
    outcome = true
    for i = 1:n
        p = generate_rand_poly(x; min_deg, max_deg, sparcity)
        printstyled("P = ", unrationalize(p), '\n'; color=:green)
        try
            q = derivative(p)
            printstyled("∂p/∂x = ", unrationalize(q), '\n'; color=:red)
            k += 1
        catch e
            println(e)
        end
    end
    outcome
end

function test_factor(x; method=:schubert_kronecker)
    ps = Any[
        x^3 + 14x^2 + 56x + 64,
        18x^3 - 57x^2 + 53x - 12,
        6x^3 - 73x^2 - 86x + 273,
        6x^4 - x^3 + 4x^2 - x - 2,
        6x^4 - 11x^3 + 8x^2 - 33x - 30,
        5x^5 - 6x^4 - 24x^3 + 20x^2 + 7x - 2,
        15x^5 - 11x^4 + 47x^3 + 27x^2 - 38x + 8,
        x^4 - 4,
        x^4 - 8x^2 - 9,
        6x^4 - 7x^3 + 5x^2 - 20x + 17,
        x^6 - 1,
        x^6 + 1,
        x^5 + x + 1,
        6x^2 + 11x - 35,
        25x^4 - 16,
        6x^4 - 19x^3 + 24x^2 - 13x + 4,
        2x^6 + 2x^5 + 2x^4 + 4x^3 + 5x^2 - 3x - 2,
        8x^5 - 48x^4 + 90x^3 - 90x^2 + 117x - 27,
        # 30x^5 + 39x^4 + 35x^3 + 25x^2 + 9x + 2,
        x^5 - x^4 - 2x^3 + 2x^2 + x - 1,
        x^10 + 2x^9 + 2x^8 + 2x^7 + x^6 + x^5 + 2x^4 + x^3 + x^2 + 2x + 1,
        x^8 - 4x^6 + 16x^2 - 16,

        𝑥^3 + 14𝑥^2 + 56𝑥 + 64,
        18𝑥^3 - 57𝑥^2 + 53𝑥 - 12,
        6𝑥^3 - 73𝑥^2 - 86𝑥 + 273,
        6𝑥^4 - 𝑥^3 + 4𝑥^2 - 𝑥 - 2,
        6𝑥^4 - 11𝑥^3 + 8𝑥^2 - 33𝑥 - 30,
        5𝑥^5 - 6𝑥^4 - 24𝑥^3 + 20𝑥^2 + 7𝑥 - 2,
        15𝑥^5 - 11𝑥^4 + 47𝑥^3 + 27𝑥^2 - 38𝑥 + 8,
        𝑥^4 - 4,
        𝑥^4 - 8𝑥^2 - 9,
        6𝑥^4 - 7𝑥^3 + 5𝑥^2 - 20𝑥 + 17,
        𝑥^6 - 1,
        𝑥^6 + 1,
        𝑥^5 + 𝑥 + 1,
        6𝑥^2 + 11𝑥 - 35,
        25𝑥^4 - 16,
        6𝑥^4 - 19𝑥^3 + 24𝑥^2 - 13𝑥 + 4,
        2𝑥^6 + 2𝑥^5 + 2𝑥^4 + 4𝑥^3 + 5𝑥^2 - 3𝑥 - 2,
        8𝑥^5 - 48𝑥^4 + 90𝑥^3 - 90𝑥^2 + 117𝑥 - 27,
        # 30𝑥^5 + 39𝑥^4 + 35𝑥^3 + 25𝑥^2 + 9𝑥 + 2,
        𝑥^5 - 𝑥^4 - 2𝑥^3 + 2𝑥^2 + 𝑥 - 1,
        𝑥^10 + 2𝑥^9 + 2𝑥^8 + 2𝑥^7 + 𝑥^6 + 𝑥^5 + 2𝑥^4 + 𝑥^3 + 𝑥^2 + 2𝑥 + 1,
        𝑥^8 - 4𝑥^6 + 16𝑥^2 - 16,
    ]

    k = 0
    outcome = true
    for p = ps
        try
            printstyled(p, '\n'; color=:green)
            f = factor(p; method=method)
            outcome = outcome && (f != nothing)
            printstyled(poly(f), '\n'; color=:red)
            println(f)
            k += 1
        catch e
            println(e)
        end
    end
    outcome
end

function test_fraction(x; n=10)
    k = 0
    outcome = true

    for i = 1:n
        p = prod(rand(1:10)*x^3 + rand(0:10)*x + rand(-15:15) for i=1:3)
        q = prod(rand(1:10)*x^2 + rand(0:10)*x + rand(-15:15) for i=1:3)

        printstyled("P/Q = ", p/q, '\n'; color=:green)

        try
            f = factor(p, q)
            print(f, '\n')
            A = poly(f)(x => π)
            B = p(x => π) / q(x => π)
            if ~isapprox(A, B)
                printstyled("incorrect: $A ≆ $B\n"; color=:red)
                outcome = false
            else
                k += 1
            end
            # f = factor(sym(p, x => 𝑥), sym(q, x => 𝑥))
            # printstyled(f, '\n'; color=:red)
        catch e
            println(e)
        end
    end
    outcome
end

function test_factor_extended(x, n, d=5; method=:roots_comb, seed=0)
    frags = [
        (x, i) -> (x - i),
        (x, i) -> (rand(2:11)*x - i),
        (x, i) -> (x^i + rand(1:100)),
        (x, i) -> (rand(1:20)*x^i + rand(1:100)),
        (x, i) -> (x^i - 1),
        (x, i) -> (x^(2^i) - i),
        (x, i) -> (rand(2:11)*x^(2^i) + i),
        (x, i) -> (x^2 + x + 1),
        (x, i) -> (x + 1),
        (x, i) -> (rand(1:20)*x + 1),
        (x, i) -> (x^3 + 7x^2 + 13x + 1),
        (x, i) -> (x^5 + x^4 + 11x^3 + 13x^2 + 5x + 1),
        (x, i) -> (rand(2:20)*x^3 + 7x^2 + 13x + 1),
        (x, i) -> (rand(2:20)*x^5 + x^4 + 11x^3 + 13x^2 + 5x + 1),
        (x, i) -> (rand(2:20)*x^7 + 1),
        (x, i) -> sum(rand(3:10)*x^j for j=1:rand(5:10)),
        (x, i) -> (x^rand(2:10)-1),
        (x, i) -> (x^rand(2:10)+1),
    ]

    if seed != 0
        Random.seed!(seed)
    end

    m = length(frags)

    for j = 1:n
        println("-------------------------------------")
        p = one(x)
        n₁ = 0
        for i = 1:d
            r = rand(1:m)
            q = frags[r](x, i)
            k = rand(1:2)
            print('(', q, ")^", k)
            if i < d
                print(" * ")
            else
                println()
            end
            p *= q^k
            n₁ += k
        end

        printstyled(p, "\n"; color=:blue)
        # printstyled("λ = ", landau_mignotte(p), '\n'; color=:blue)
        try
            f = factor(p; method=method)
            printstyled(f, "\n"; color=:magenta)

            if p isa AbstractPolynomialLike
                println("Δ = ", poly(f) - p)
            else
                println("Δ = ", expand(simplify(f) - p))
            end

            n₂ = sum(last(v) for v in f.factors)

            if n₂ < n₁
                printstyled("possible incomplete factorization\n"; color=:yellow)
            else
                printstyled("OK!\n"; color=:green)
            end
        catch e
            println(e)
        end
    end
    true
end

function swinnerton_dyer(x, n)
    a = zeros(n)
    r = 2
    for i = 1:n
        a[i] = sqrt(r)
        r = nextprime(r+1)
    end

    p = one(x)
    for i = 0:2^n-1
        k = sum(isodd(i >> (j-1)) ? a[j] : -a[j] for j = 1:n)
        p *= x + k
    end

    return integer_poly(p)
end

#############################################################################

function test_all()
    @test test_eq(x, (p,q)->(p+q)-p-q, "add")
    @test test_eq(x, (p,q)->p-(p÷q)*q-(p%q), "mul")
    @test test_eq(x, (p,q)->p % gcd(p,q)+q % gcd(p,q), "gcd"; max_deg=5)
    @test test_deriv(x)
    println("********* Schubert-Kronecker ************")
    @test test_factor(x; method=:schubert_kronecker)
    # println("********** Roundabout *******************")
    # @test test_factor(x; method=:roundabout)
    println("********* Roots Combinations*************")
    @test test_factor(x; method=:roots_comb)
    println("************* Roots SSP *****************")
    @test test_factor(x; method=:roots_SSP)
    @test test_factor_extended(x, 50; method=:roots_SSP)
    @test test_factor_extended(𝑥, 50; method=:roots_SSP)
    @test test_fraction(x)
end

@testset "arith" begin test_all() end

ext: jl


text:
include(joinpath(@__DIR__, "..", "..", "dynamics/quaternions.jl"))

function visualize!(vis, p, q; Δt = 0.1)
	setvisible!(vis["/Background"], true)
	setprop!(vis["/Background"], "top_color", RGBA(1.0, 1.0, 1.0, 1.0))
	setprop!(vis["/Background"], "bottom_color", RGBA(1.0, 1.0, 1.0, 1.0))
	setvisible!(vis["/Axes"], false)
	setvisible!(vis["/Grid"], false)

    setobject!(vis[:satellite],
    	Rect(Vec(-0.25, -0.25, -0.25),Vec(0.5, 0.5, 0.5)),
    	MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))

    arrow_x = ArrowVisualizer(vis[:satellite][:arrow_x])
    mat = MeshPhongMaterial(color=RGBA(1.0, 0.0, 0.0, 1.0))
    setobject!(arrow_x, mat)
    settransform!(arrow_x,
    	Point(0.0, 0.0, 0.0),
    	Vec(0.75, 0.0, 0.0),
    	shaft_radius=0.05,
    	max_head_radius=0.1)

    arrow_y = ArrowVisualizer(vis[:satellite][:arrow_y])
    mat = MeshPhongMaterial(color=RGBA(0.0, 1.0, 0.0, 1.0))
    setobject!(arrow_y, mat)
    settransform!(arrow_y,
    	Point(0.0, 0.0, 0.0),
    	Vec(0.0, 0.75, 0.0),
    	shaft_radius=0.05,
    	max_head_radius=0.1)

    arrow_z = ArrowVisualizer(vis[:satellite][:arrow_z])
    mat = MeshPhongMaterial(color=RGBA(0.0, 0.0, 1.0, 1.0))
    setobject!(arrow_z, mat)
    settransform!(arrow_z,
    	Point(0.0, 0.0, 0.0),
    	Vec(0.0, 0.0, 0.75),
    	shaft_radius=0.05,
    	max_head_radius=0.1)

    anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))

     # for t = 1:length(q)
    	#  MeshCat.atframe(anim, t) do
    	# 	 settransform!(vis["satellite"],
    	# 		   compose(Translation((q[t][1:3] + [-0.25; -0.25; -0.25])...),
    	# 				 LinearMap(quaternion_rotation_matrix(q[t][4:7]))))
    	#  end
     # end

	 for t = 1:length(q)
		MeshCat.atframe(anim, t) do
			settransform!(vis["satellite"],
				  compose(Translation(([-0.0; -0.0; -0.0])...),
						LinearMap(quaternion_rotation_matrix(q[t][1:4]))))
		end
	 end
	#
	#
    MeshCat.setanimation!(vis, anim)
end

vis = Visualizer()
# render(vis)
open(vis)

nq = 4
nv = 3
H = [zeros(1, 3); I]

J = Diagonal([1.0; 2.0; 3.0])

h = 0.01
T = 1000

function f(x, u)
	q = x[1:4]
	ω = x[5:7]
	τ = u[1:3]

	[0.5 * L_multiply(q) * H * ω;
	 J \ (τ - cross(ω, J * ω))]
 end

function f(x, u, h)
    x + h * f(x + 0.5 * h * f(x, u), u) # explicit midpoint
end

q1 = [1.0; 0.0; 0.0; 0.0]
ω1 = [0.0; 10.0; 0.1]

x1 = [q1; ω1]
u1 = zeros(3)

x = [copy(x1)]
for t = 1:T
	push!(x, f(x[end], u1, h))
end

visualize!(vis, nothing, x, Δt = h)

# implicit
function f(x⁺, x, u, h)
	x⁺ - (x + h * f(0.5 * (x + x⁺), u)) # implciit midpoint
end

function newton(x, u, h)
	x⁺ = copy(x)
	_r(z) = f(z, x, u, h)

	function cand(x, Δ, α)
		x⁺ = copy(x)
		x⁺[1:4] = L_multiply(x[1:4]) * φ(α * Δ[1:3])
		x⁺[5:7] += α * Δ[4:6]

		return x⁺
	end

	r = _r(x⁺)
	@show norm(r)

	for i = 1:100
		∇r = ForwardDiff.jacobian(_r, x⁺) * [attitude_jacobian(x⁺) zeros(4, 3); zeros(3,3) I]
		Δ = -(∇r' * ∇r) \ (∇r' * r)
		α = 1.0

		x_cand = cand(x⁺, Δ, α)
		r_cand = _r(x_cand)
		@show norm(r_cand)
		iter = 0

		while norm(r_cand) >= norm(r)
			α *= 0.5
			x_cand = cand(x⁺, Δ, α)
			r_cand = _r(x_cand)
			@show norm(r_cand)

			iter += 1
			if iter > 100
				@error "line search failure"
				break
			end
		end

		x⁺ = x_cand
		r = r_cand

		if norm(r) < 1.0e-8
			return x⁺
		end
	end
	@warn "newton failure"
	return x
end

newton(x1, u1, h)

x = [copy(x1)]
for t = 1:T
	push!(x, newton(x[end], u1, h))
end

visualize!(vis, nothing, x, Δt = h)

function newton_variational(x, u, h)
	q2 = copy(x)[1:4]
	ω2 = copy(x)[5:7]

	function f(ω1, ω2, q2, u, h)
		(J * ω2 * sqrt(4.0 / h^2.0 - ω2' * ω2)
			+ cross(ω2, J * ω2)
			- J * ω1 * sqrt(4.0 / h^2.0 - ω1' * ω1)
			+ cross(ω1, J * ω1)
			- 2.0 * u[1:3])
	end

	_r(z) = f(x[5:7], z, x[1:4], u, h)

	function cand(x, Δ, α)
		x⁺ = copy(x)
		# x⁺[1:4] = L_multiply(x[1:4]) * φ(α * Δ[1:3])
		x⁺[1:3] += α * Δ[1:3]

		return x⁺
	end

	r = _r(ω2)
	@show norm(r)

	for i = 1:100
		∇r = ForwardDiff.jacobian(_r, ω2)# * [attitude_jacobian(x⁺) zeros(4, 3); zeros(3,3) I]
		Δ = -(∇r' * ∇r) \ (∇r' * r)
		α = 1.0

		ω_cand = cand(ω2, Δ, α)
		r_cand = _r(ω_cand)
		@show norm(r_cand)
		iter = 0

		while norm(r_cand) >= norm(r)
			α *= 0.5
			x_cand = cand(ω2, Δ, α)
			r_cand = _r(ω_cand)
			@show norm(r_cand)

			iter += 1
			if iter > 100
				@error "line search failure"
				break
			end
		end

		ω2 = ω_cand
		r = r_cand

		if norm(r) < 1.0e-10
			q3 = 0.5 * h * L_multiply(q2) * [sqrt((2.0 / h)^2.0 - ω2' * ω2); ω2]
			return [q3; ω2]
		end
	end
	@warn "newton failure"
	return x
end

newton_variational(x1, u1, h)

x = [copy(x1)]
for t = 1:T
	push!(x, newton_variational(x[end], u1, h))
end

visualize!(vis, nothing, x, Δt = h)

q_variational = [xt[1:4] for xt in x]
w_variational = [xt[5:7] for xt in x]
plot(hcat(w_variational...)')
plot(hcat(q_variational...)')
plot([norm(q) for q in q_variational])

function newton_variational_2(q1, q2, u, h)
	q3 = copy(q2)

	# @show ω_finite_difference(q1, q2, h)
	# @show norm(q3)
	# @show q3
	function f(q1, q2, q3, u, h)
		ω1 = ω_finite_difference(q1, q2, h)
		ω2 = ω_finite_difference(q2, q3, h)

		(J * ω2 * sqrt(4.0 / h^2.0 - ω2' * ω2)
			+ cross(ω2, J * ω2)
			- J * ω1 * sqrt(4.0 / h^2.0 - ω1' * ω1)
			+ cross(ω1, J * ω1)
			- 2.0 * u[1:3])
	end

	_r(z) = f(q1, q2, z, u, h)

	function cand(q3, Δ, α)
		return L_multiply(q3) * φ(α * Δ[1:3])
	end

	r = _r(q3)
	@show norm(r)

	for i = 1:100
		∇r = ForwardDiff.jacobian(_r, q3) * attitude_jacobian(q3)
		Δ = -(∇r' * ∇r) \ (∇r' * r)
		α = 1.0

		q3_cand = cand(q3, Δ, α)
		r_cand = _r(q3_cand)
		@show norm(r_cand)
		iter = 0

		while norm(r_cand) > norm(r)
			α *= 0.5
			q3_cand = cand(q3, Δ, α)
			r_cand = _r(q3_cand)
			@show norm(r_cand)

			iter += 1
			if iter > 100
				@error "line search failure"
				break
			end
		end

		q3 = q3_cand
		r = r_cand

		# @show ω_finite_difference(q2, q3, h)
		# @show norm(q3)
		# @show q3
		if norm(r) < 1.0e-10
			return q3
		end
	end
	@warn "newton failure"
	return q2
end

x = [copy(q1), copy(Array(0.5 * h * L_multiply(q1) * [sqrt((2.0 / h)^2.0 - ω1' * ω1); ω1]))]

for t = 1:500
	push!(x, newton_variational_2(x[end-1], x[end], u1, h))
end

visualize!(vis, nothing, x, Δt = h)

function newton_variational_3(q1, q2, u, h)
	q3 = copy(q2)

	# @show ω_finite_difference(q1, q2, h)
	# @show norm(q3)
	# @show q3
	function f(q1, q2, q3, u, h)
		v1 = (q2[1:3] - q1[1:3]) ./ h
		v2 = (q3[1:3] - q2[1:3]) ./ h

		ω1 = ω_finite_difference(q1[4:7], q2[4:7], h)
		ω2 = ω_finite_difference(q2[4:7], q3[4:7], h)

		D1L1, D2L1 = [0.0; 0.0; -9.81], v1
		D1L2, D2L2 = [0.0; 0.0; -9.81], v2

		[0.5 * h * D1L1 + D2L1 + 0.5 * h * D1L2 - D2L2;
			(J * ω2 * sqrt(4.0 / h^2.0 - ω2' * ω2)
			+ cross(ω2, J * ω2)
			- J * ω1 * sqrt(4.0 / h^2.0 - ω1' * ω1)
			+ cross(ω1, J * ω1)
			- 2.0 * u[1:3])]
	end

	_r(z) = f(q1, q2, z, u, h)

	function cand(q3, Δ, α)
		return [q3[1:3] + α .* Δ[1:3];
			    L_multiply(q3[4:7]) * φ(α * Δ[4:6])]
	end

	r = _r(q3)
	@show norm(r)

	for i = 1:100
		∇r = ForwardDiff.jacobian(_r, q3) * [I zeros(3, 3);
		                                      zeros(4, 3) attitude_jacobian(q3[4:7])]
		# Δ = -(∇r' * ∇r) \ (∇r' * r)
		Δ = -∇r \ r

		α = 1.0

		q3_cand = cand(q3, Δ, α)
		r_cand = _r(q3_cand)
		@show norm(r_cand)
		iter = 0

		while norm(r_cand) > norm(r)
			α *= 0.5
			q3_cand = cand(q3, Δ, α)
			r_cand = _r(q3_cand)
			@show norm(r_cand)

			iter += 1
			if iter > 100
				@error "line search failure"
				break
			end
		end

		q3 = q3_cand
		r = r_cand

		# @show ω_finite_difference(q2, q3, h)
		# @show norm(q3)
		# @show q3
		if norm(r) < 1.0e-10
			return q3
		end
	end
	@warn "newton failure"
	return q2
end

p1 = [0.0; 0.0; 0.0]
v1 = [10.0; 0.0; 5.0]
x = [[copy(p1); copy(q1)],
	 [copy(p1) + h * copy(v1); copy(Array(0.5 * h * L_multiply(q1) * [sqrt((2.0 / h)^2.0 - ω1' * ω1); ω1]))]]

for t = 1:100
	push!(x, newton_variational_3(x[end-1], x[end], u1, h))
end
include(joinpath(module_dir(), "src/dynamics/rigid_body/visuals.jl"))

visualize!(vis, nothing, x, Δt = h)

struct RQuat <: Space
	n::Int
	r_idx
	Δr_idx
	quat_idx
	Δquat_idx
end

function rquat(dim, r_idx, Δr_idx, quat_idx, Δquat_idx)
	RQuat(dim, r_idx, Δr_idx, quat_idx, Δquat_idx)
end

rquat_space = rquat(6,
			collect((1:3)),
			collect((1:3)),
			collect((4:7)),
			collect((4:6)))

function candidate_point!(z̄::Vector{T}, s::RQuat, z::Vector{T}, Δ::Vector{T}, α::T) where T
    z̄[s.r_idx] .= z[s.r_idx] - α .* Δ[s.Δr_idx]
	z̄[s.quat_idx] .= L_multiply(z[s.quat_idx]) * φ(-1.0 * α .* Δ[s.Δquat_idx])
	return nothing
end

function f(q1, q2, q3, u, h)
	v1 = (q2[1:3] - q1[1:3]) ./ h
	v2 = (q3[1:3] - q2[1:3]) ./ h

	ω1 = ω_finite_difference(q1[4:7], q2[4:7], h)
	ω2 = ω_finite_difference(q2[4:7], q3[4:7], h)

	D1L1, D2L1 = [0.0; 0.0; -9.81], v1
	D1L2, D2L2 = [0.0; 0.0; -9.81], v2

	[0.5 * h * D1L1 + D2L1 + 0.5 * h * D1L2 - D2L2;
		-1.0 * (J * ω2 * sqrt(4.0 / h^2.0 - ω2' * ω2)
		+ cross(ω2, J * ω2)
		- J * ω1 * sqrt(4.0 / h^2.0 - ω1' * ω1)
		+ cross(ω1, J * ω1)
		- 2.0 * u[1:3])]
end

function G_func(x)
	q = x[4:7]
	[1.0 0.0 0.0 0.0 0.0 0.0;
     0.0 1.0 0.0 0.0 0.0 0.0;
	 0.0 0.0 1.0 0.0 0.0 0.0;
     zeros(4, 3) attitude_jacobian(q)]
end

function r!(r, z, θ, κ)
	q1 = θ[1:7]
	q2 = θ[7 .+ (1:7)]
	u1 = θ[2 * 7 .+ (1:3)]
	h = θ[2 * 7 + 3 + 1]
	r .= f(q1, q2, z, u1, h)
	nothing
end

function rz!(rz, z, θ)
	q1 = θ[1:7]
	q2 = θ[7 .+ (1:7)]
	u1 = θ[2 * 7 .+ (1:3)]
	h = θ[2 * 7 + 3 + 1]
	_f(y) = f(q1, q2, y, u1, h)
	rz .= ForwardDiff.jacobian(_f, z) * G_func(z)
end

rzz = zeros(6,6)
rz!(rzz, rand(7), rand(18))

# options
opts = ContactControl.InteriorPointOptions(
	r_tol = 1.0e-6, κ_tol = 1.0e-6,
	diff_sol = false,
	solver = :lu_solver)
#
# rr = zeros(6)
# r!(rr, z, θ, 0.0)
# rr
p1 = [0.0; 0.0; 0.0]
v1 = [10.0; 0.0; 5.0]
ω1 = [0.0; 10.0; 0.1]
x = [[copy(p1); copy(q1)],
	 [copy(p1) + h * copy(v1); copy(Array(0.5 * h * L_multiply(q1) * [sqrt((2.0 / h)^2.0 - ω1' * ω1); ω1]))]]

z = copy(x[end])
θ = [x[end-1]; x[end]; u1; h]

# solver
ip = ContactControl.interior_point(z, θ,
	s = rquat_space,
	r! = r!, rz! = rz!,
	rz = zeros(6, 6),
	opts = opts)

# solve
for t = 1:T
	ip.z .= copy(x[end])
	ip.θ .= [x[end-1]; x[end]; u1; h]

	status = ContactControl.interior_point_solve!(ip)

	push!(x, copy(ip.z))
end

include(joinpath(module_dir(), "src/dynamics/rigid_body/visuals.jl"))

visualize!(vis, nothing, x, Δt = h)


text:
"""
A wrapper for wasm_table_t specific to wasmtime, since the wasmtime c-api is currently missing
wasm_ref_as_func()
"""
mutable struct WasmTable <: AbstractVector{WasmFuncRef}
    store::WasmStore
    wasm_table_ptr::Ptr{wasm_table_t}
end
function WasmTable(store::WasmStore, limits::Pair{Int,Int})
    wasm_limits = wasm_limits_t(limits...)
    table_valtype = wasm_valtype_new(WASM_FUNCREF)
    @assert table_valtype != C_NULL "Failed to create table valtype"

    table_type = GC.@preserve wasm_limits wasm_tabletype_new(
        table_valtype,
        pointer_from_objref(wasm_limits),
    )
    @assert table_type != C_NULL "Failed to create table type"

    wasm_table_ptr = wasm_table_new(store, table_type, C_NULL)
    @assert wasm_table_ptr != C_NULL "Failed to create table"

    # finalizer(wasm_table_delete, WasmTable(store, wasm_table_ptr))
    WasmTable(store, wasm_table_ptr)
end

Base.unsafe_convert(::Type{Ptr{wasm_table_t}}, wasm_table::WasmTable) =
    wasm_table.wasm_table_ptr

Base.IndexStyle(::Type{WasmTable}) = IndexLinear()
Base.size(table::WasmTable) = (wasm_table_size(table) |> Int,)
function Base.getindex(table::WasmTable, i::Int)
    wasm_val_ref = Ref(wasm_val_t(tuple((zero(UInt8) for _ = 1:16)...)))

    index = i - 1

    wasmtime_context = @ccall libwasmtime.wasmtime_store_context(table.store::Ptr{wasm_store_t})::Ptr{Nothing}
    res = @ccall libwasmtime.wasmtime_table_get(
        wasmtime_context::Ptr{Nothing},
        table::Ptr{wasm_table_t},
        index::Cint,
        wasm_val_ref::Ptr{wasm_val_t},
    )::Bool
    @assert res "Failed to access index $i"
    @assert wasm_val_ref[].kind == WASM_FUNCREF

    WasmFuncRef(wasm_func_ptr[])
end
function Base.setindex!(table::WasmTable, v, i::Int)
    wasm_table_set(table, v, i - 1) # TODO: use result
end

map_to_extern(extern_table::WasmTable) = wasm_table_as_extern(extern_table)


text:
"""
	struct ExtractArray{T}
		item::T
	end

	convert array of values to one bag of values converted with `item`. Note that in order the function to work properly,
	calling `item` on a single item has to return Matrix.

```juliadoctest
julia> sc = ExtractArray(ExtractCategorical(Float32,2:4))
julia> sc([2,3,1,4]).data
3×4 Array{Float32,2}:
 1.0  0.0  0.0  0.0
 0.0  1.0  0.0  0.0
 0.0  0.0  0.0  1.0

```

```juliadoctest
julia> sc = ExtractArray(ExtractScalar())
julia> sc([2,3,4]).data
 2.0  3.0  4.0
```
"""
struct ExtractArray{T} <: AbstractExtractor
	item::T
end

extractsmatrix(s::ExtractArray) = false

# (s::ExtractArray)(v::V) where {V<:Nothing} = BagNode(reduce(catobs, s.item.([nothing])),[1:1])
(s::ExtractArray)(v::V) where {V<:Nothing} = BagNode(missing, [0:-1])
(s::ExtractArray)(v::V) where {V<:Missing} = BagNode(missing, [0:-1])
(s::ExtractArray)(v) = isempty(v) ? BagNode(missing, [0:-1]) : BagNode(reduce(catobs, s.item.(v)),[1:length(v)])
function Base.show(io::IO, m::ExtractArray; pad = [], key::String="")
	c = COLORS[(length(pad)%length(COLORS))+1]
	key *= isempty(key) ? "" : ": "
	paddedprint(io,"$(key)Array of\n", color = c)
	paddedprint(io, "  └── ", color=c, pad=pad)
	show(io,m.item, pad = [pad; (c, "      ")])
end

Base.hash(e::ExtractArray, h::UInt) = hash(e.item, h)
Base.:(==)(e1::ExtractArray, e2::ExtractArray) = e1.item == e2.item


text:
function mdl_avr_ode!(device_states,
                    output_ode,
                    device::DynGenerator{M, S, AVRFixed, TG, P})  where {M <: Machine,
                                                                         S <: Shaft,
                                                                         TG <: TurbineGov,
                                                                         P <: PSS}

    #Obtain parameters
    Vf = device.avr.Emf

    #Update inner vars
    device.inner_vars[Vf_var] = Vf

    return
end

function mdl_avr_ode!(device_states,
                      output_ode,
                      device::DynGenerator{M, S, AVRSimple, TG, P})  where {M <: Machine,
                                                                            S <: Shaft,
                                                                            TG <: TurbineGov,
                                                                            P <: PSS}

    #Obtain references
    V_ref = get_V_ref(device)

    local_ix = device.local_state_ix[device.avr]

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vf = internal_states[1]

    #Define external states for device
    V_th = sqrt(device.inner_vars[VR_gen_var]^2 + device.inner_vars[VI_gen_var]^2)

    #Get Parameters
    Kv = device.avr.Kv

    #Compute ODEs
    output_ode[local_ix[1]] = Kv*(V_ref - V_th)

    #Update inner_vars
    device.inner_vars[Vf_var] = Vf

    return
end

"""
Model of AVR Type I in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""

function mdl_avr_ode!(device_states,
                      output_ode,
                      device::DynGenerator{M, S, AVRTypeI, TG, P})  where {M <: Machine,
                                                                           S <: Shaft,
                                                                           TG <: TurbineGov,
                                                                           P <: PSS}

    #Obtain references
    V0_ref = get_V_ref(device)

    local_ix = device.local_state_ix[device.avr]

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vf = internal_states[1]
    Vr1 = internal_states[2]
    Vr2 = internal_states[3]
    Vm  = internal_states[4]

    #Define external states for device
    V_th = sqrt(device.inner_vars[VR_gen_var]^2 + device.inner_vars[VI_gen_var]^2)
    Vs = device.inner_vars[V_pss_var]

    #Get parameters
    Ka = device.avr.Ka
    Ke = device.avr.Ke
    Kf = device.avr.Kf
    Ta = device.avr.Ta
    Te = device.avr.Te
    Tf = device.avr.Tf
    Tr = device.avr.Tr
    Vr_max = device.avr.Vr_max
    Vr_min = device.avr.Vr_min
    Ae = device.avr.Ae
    Be = device.avr.Be

    #Compute auxiliary parameters
    Se_Vf = Ae*exp(Be*abs(Vf)) #16.13
    V_ref = V0_ref + Vs

    #Set anti-windup for Vr1. NOT WORKING FOR INITIALIZATION USING NLSOLVE
    #if Vr1 > Vr_max
    #    Vr1 = Vr_max
    #elseif Vr1 < Vr_min
    #    Vr1 = Vr_min
    #end

    #Compute 4 States AVR ODE:
    output_ode[local_ix[1]] = -(1.0/Te)*(Vf*(Ke + Se_Vf) - Vr1) #16.12c
    output_ode[local_ix[2]] = (1.0/Ta)*(Ka*(V_ref - Vm - Vr2 - (Kf/Tf)*Vf) - Vr1) #16.12a
    output_ode[local_ix[3]] = -(1.0/Tf)*((Kf/Tf)*Vf + Vr2) #16.12b
    output_ode[local_ix[4]] = (1.0/Tr)*(V_th - Vm) #16.11

    #Update inner_vars
    device.inner_vars[Vf_var] = Vf

    return
end

"""
Model of AVR Type II in Julia.
Refer to Power System Modelling and Scripting by F. Milano for the equations
"""

function mdl_avr_ode!(device_states,
                      output_ode,
                      device::DynGenerator{M, S, AVRTypeII, TG, P})  where {M <: Machine,
                                                                           S <: Shaft,
                                                                           TG <: TurbineGov,
                                                                           P <: PSS}

    #Obtain references
    V0_ref = get_V_ref(device)

    local_ix = device.local_state_ix[device.avr]

    #Define inner states for component
    internal_states = @view device_states[local_ix]
    Vf = internal_states[1]
    Vr1 = internal_states[2]
    Vr2 = internal_states[3]
    Vm  = internal_states[4]

    #Define external states for device
    V_th = sqrt(device.inner_vars[VR_gen_var]^2 + device.inner_vars[VI_gen_var]^2)
    Vs = device.inner_vars[V_pss_var]

    #Get parameters
    K0 = device.avr.K0
    T1 = device.avr.T1
    T2 = device.avr.T2
    T3 = device.avr.T3
    T4 = device.avr.T4
    Te = device.avr.Te
    Tr = device.avr.Tr
    Vr_max = device.avr.Vr_max
    Vr_min = device.avr.Vr_min
    Ae = device.avr.Ae
    Be = device.avr.Be

    #Compute auxiliary parameters
    Se_Vf = Ae*exp(Be*abs(Vf)) #16.13
    V_ref = V0_ref + Vs
    Vr = K0*Vr2 + (T4/T3)*(Vr1 + K0*(T2/T1)*(V_ref - Vm)) #16.21

    #Set anti-windup for Vr1. NOT WORKING FOR INITIALIZATION USING NLSOLVE
    #if Vr > Vr_max
    #    Vr = Vr_max
    #elseif Vr < Vr_min
    #    Vr = Vr_min
    #end

    #Compute 4 States AVR ODE:
    output_ode[local_ix[1]] = -(1.0/Te)*(Vf*(1.0 + Se_Vf) - Vr) #16.18
    output_ode[local_ix[2]] = (1.0/T1)*(K0*(1.0 - (T2/T1))*(V_ref - Vm) - Vr1) #16.14
    output_ode[local_ix[3]] = (1.0/(K0*T3))*((1.0 - (T4/T3))*(Vr1 + K0*(T2/T1)*(V_ref - Vm)) - K0*Vr2)  #16.20
    output_ode[local_ix[4]] = (1.0/Tr)*(V_th - Vm) #16.11

    #Update inner_vars
    device.inner_vars[Vf_var] = Vf

    return
end


text:
"""
    DynamicalSystemSetup

An abstract type for setup schemes for dynamical systems
"""
abstract type DynamicalSystemSetup end

export DynamicalSystemSetup


text:
using Documenter
using ConcurrentUtils

makedocs(
    sitename = "ConcurrentUtils",
    format = Documenter.HTML(),
    modules = [ConcurrentUtils],
    strict = [
        :autodocs_block,
        :cross_references,
        :docs_block,
        :doctest,
        :eval_block,
        :example_block,
        :footnote,
        :linkcheck,
        :meta_block,
        # :missing_docs,
        :parse_error,
        :setup_block,
    ],
    # Ref:
    # https://juliadocs.github.io/Documenter.jl/stable/lib/public/#Documenter.makedocs
)

deploydocs(
    repo = "github.com/JuliaConcurrent/ConcurrentUtils.jl",
    devbranch = "main",
    push_preview = true,
    # Ref:
    # https://juliadocs.github.io/Documenter.jl/stable/lib/public/#Documenter.deploydocs
)


text:
import Pkg

Pkg.add("Agents")
Pkg.add(name="MetaGraphs",version="0.6.7")
Pkg.add("AgentsPlots")
Pkg.add("Plots")
Pkg.add("LightGraphs")
Pkg.add("SimpleWeightedGraphs")
Pkg.add("GraphPlot")
Pkg.add("GraphRecipes")
Pkg.add("NetworkLayout")
Pkg.add("Tables")
Pkg.add("DataFrames")
Pkg.add("CSV")
Pkg.add("JSON")
Pkg.add("Serialization")
Pkg.add("DelimitedFiles")
Pkg.add("BritishNationalGrid")
Pkg.add("ZipFile")
Pkg.add("Shapefile")
Pkg.add("Random")
Pkg.add("Match")
Pkg.add("LinearAlgebra")
Pkg.add("StatsBase")
Pkg.add("Distributions")
Pkg.add("StatsPlots")
Pkg.add("SparseArrays")
Pkg.add("Laplacians")
Pkg.add("DataStructures")
Pkg.add("RollingFunctions")
Pkg.add("BenchmarkTools")
Pkg.add("Statistics")
Pkg.add("Logging")
Pkg.add("LoggingExtras")
Pkg.add("LoggingFacilities")
Pkg.add("Dates")
Pkg.add("TimeZones")
Pkg.add("PyCall")
Pkg.add(url="../NAssets.jl")
Pkg.precompile()


text:
mutable struct image_t <: LCMType
    utime::Int64
    width::Int32
    height::Int32
    row_stride::Int32
    pixelformat::Int32
    size::Int32
    data::Vector{UInt8}
    nmetadata::Int32
    metadata::Vector{image_metadata_t}
end

@lcmtypesetup(image_t,
    data => (size,),
    metadata => (nmetadata,)
)


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

mutable struct IoK8sApiCoreV1PodStatus <: SwaggerModel
    conditions::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PodCondition} } # spec name: conditions
    containerStatuses::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ContainerStatus} } # spec name: containerStatuses
    ephemeralContainerStatuses::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ContainerStatus} } # spec name: ephemeralContainerStatuses
    hostIP::Any # spec type: Union{ Nothing, String } # spec name: hostIP
    initContainerStatuses::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1ContainerStatus} } # spec name: initContainerStatuses
    message::Any # spec type: Union{ Nothing, String } # spec name: message
    nominatedNodeName::Any # spec type: Union{ Nothing, String } # spec name: nominatedNodeName
    phase::Any # spec type: Union{ Nothing, String } # spec name: phase
    podIP::Any # spec type: Union{ Nothing, String } # spec name: podIP
    podIPs::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1PodIP} } # spec name: podIPs
    qosClass::Any # spec type: Union{ Nothing, String } # spec name: qosClass
    reason::Any # spec type: Union{ Nothing, String } # spec name: reason
    startTime::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1Time } # spec name: startTime

    function IoK8sApiCoreV1PodStatus(;conditions=nothing, containerStatuses=nothing, ephemeralContainerStatuses=nothing, hostIP=nothing, initContainerStatuses=nothing, message=nothing, nominatedNodeName=nothing, phase=nothing, podIP=nothing, podIPs=nothing, qosClass=nothing, reason=nothing, startTime=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("conditions"), conditions)
        setfield!(o, Symbol("conditions"), conditions)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("containerStatuses"), containerStatuses)
        setfield!(o, Symbol("containerStatuses"), containerStatuses)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("ephemeralContainerStatuses"), ephemeralContainerStatuses)
        setfield!(o, Symbol("ephemeralContainerStatuses"), ephemeralContainerStatuses)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("hostIP"), hostIP)
        setfield!(o, Symbol("hostIP"), hostIP)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("initContainerStatuses"), initContainerStatuses)
        setfield!(o, Symbol("initContainerStatuses"), initContainerStatuses)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("message"), message)
        setfield!(o, Symbol("message"), message)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("nominatedNodeName"), nominatedNodeName)
        setfield!(o, Symbol("nominatedNodeName"), nominatedNodeName)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("phase"), phase)
        setfield!(o, Symbol("phase"), phase)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("podIP"), podIP)
        setfield!(o, Symbol("podIP"), podIP)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("podIPs"), podIPs)
        setfield!(o, Symbol("podIPs"), podIPs)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("qosClass"), qosClass)
        setfield!(o, Symbol("qosClass"), qosClass)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("reason"), reason)
        setfield!(o, Symbol("reason"), reason)
        validate_property(IoK8sApiCoreV1PodStatus, Symbol("startTime"), startTime)
        setfield!(o, Symbol("startTime"), startTime)
        o
    end
end # type IoK8sApiCoreV1PodStatus

const _property_map_IoK8sApiCoreV1PodStatus = Dict{Symbol,Symbol}(Symbol("conditions")=>Symbol("conditions"), Symbol("containerStatuses")=>Symbol("containerStatuses"), Symbol("ephemeralContainerStatuses")=>Symbol("ephemeralContainerStatuses"), Symbol("hostIP")=>Symbol("hostIP"), Symbol("initContainerStatuses")=>Symbol("initContainerStatuses"), Symbol("message")=>Symbol("message"), Symbol("nominatedNodeName")=>Symbol("nominatedNodeName"), Symbol("phase")=>Symbol("phase"), Symbol("podIP")=>Symbol("podIP"), Symbol("podIPs")=>Symbol("podIPs"), Symbol("qosClass")=>Symbol("qosClass"), Symbol("reason")=>Symbol("reason"), Symbol("startTime")=>Symbol("startTime"))
const _property_types_IoK8sApiCoreV1PodStatus = Dict{Symbol,String}(Symbol("conditions")=>"Vector{IoK8sApiCoreV1PodCondition}", Symbol("containerStatuses")=>"Vector{IoK8sApiCoreV1ContainerStatus}", Symbol("ephemeralContainerStatuses")=>"Vector{IoK8sApiCoreV1ContainerStatus}", Symbol("hostIP")=>"String", Symbol("initContainerStatuses")=>"Vector{IoK8sApiCoreV1ContainerStatus}", Symbol("message")=>"String", Symbol("nominatedNodeName")=>"String", Symbol("phase")=>"String", Symbol("podIP")=>"String", Symbol("podIPs")=>"Vector{IoK8sApiCoreV1PodIP}", Symbol("qosClass")=>"String", Symbol("reason")=>"String", Symbol("startTime")=>"IoK8sApimachineryPkgApisMetaV1Time")
Base.propertynames(::Type{ IoK8sApiCoreV1PodStatus }) = collect(keys(_property_map_IoK8sApiCoreV1PodStatus))
Swagger.property_type(::Type{ IoK8sApiCoreV1PodStatus }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1PodStatus[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1PodStatus }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1PodStatus[property_name]

function check_required(o::IoK8sApiCoreV1PodStatus)
    true
end

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


text:
# CUDA CuContext

immutable CuContext
	handle::Ptr{Void}
end

const CTX_SCHED_AUTO  = 0x00
const CTX_SCHED_SPIN  = 0x01
const CTX_SCHED_YIELD = 0x02
const CTX_SCHED_BLOCKING_SYNC = 0x04
const CTX_MAP_HOST = 0x08
const CTX_LMEM_RESIZE_TO_MAX = 0x10

function create_context(dev::CuDevice, flags::Integer)
	a = Array(Ptr{Void}, 1)
	@cucall(:cuCtxCreate, (Ptr{Ptr{Void}}, Cuint, Cint), a, flags, dev.handle)
	return CuContext(a[1])
end

create_context(dev::CuDevice) = create_context(dev, 0)

function destroy(ctx::CuContext)
	@cucall(:cuCtxDestroy, (Ptr{Void},), ctx.handle)
end

function push(ctx::CuContext)
	@cucall(:cuCtxPushCurrent, (Ptr{Void},), ctx.handle)
end

function pop(ctx::CuContext)
	a = Array(Ptr{Void}, 1)
	@cucall(:cuCtxPopCurrent, (Ptr{Ptr{Void}},), a)
	return CuContext(a[1])
end


text:
# Determine the mesh sensitivity for the default unexposed HTC and the recovered value of
# the exposed HTC. This is ran post-factum of the recovery calculation to justify it.
# In case the correct mesh was not picked the adjoint calculation muts be repeated with
# the converged value of the mesh size.

using NeuralThermalODE

const testfile = "Fig_67.csv_Result.csv"
const test_data = read_data(testfile)

# This is either the raw HTC or a spline fit.
const htc_file = "checkpoint_htc_360.csv"
const htc_recovered = read_htc(htc_file)

const sinkT = test_data[end, 2]
const ambT = test_data[1, 2]

tspan = (Float64(test_data[1, 1]), Float64(test_data[end, 1]))

const p = Vector(htc_recovered[2])

# Pick correct value.
const t = 0.00635

# Study the mesh sensitivity.
const N1 = 10
const N2 = 20
const N3 = 40
const N4 = 80

dx = t / N1
u = [ambT for i in 1:N1]
prob = ODEProblem(heat_transfer, u, tspan, p, saveat = 1.0)
sol1 = solve(prob, Tsit5(), saveat = 1.0)

dx = t / N2
u = [ambT for i in 1:N2]
prob = ODEProblem(heat_transfer, u, tspan, p, saveat = 1.0)
sol2 = solve(prob, Tsit5(), saveat = 1.0)

dx = t / N3
u = [ambT for i in 1:N3]
prob = ODEProblem(heat_transfer, u, tspan, p, saveat = 1.0)
sol3 = solve(prob, Tsit5(), saveat = 1.0)

dx = t / N4
u = [ambT for i in 1:N4]
prob = ODEProblem(heat_transfer, u, tspan, p, saveat = 1.0)
sol4 = solve(prob, Tsit5(), saveat = 1.0)

plot(sol1.t[:], sol1[end, :], label = "N = 10", color = :black, linestyle = :dash)
plot!(sol2.t[:], sol2[end, :], label = "N = 20", color = :black, linestyle = :dot)
plot!(sol3.t[:], sol3[end, :], label = "N = 40", color = :black, linestyle = :dashdot)
plot!(sol4.t[:], sol4[end, :], label = "N = 80", color = :black)
plot!(fmt = :svg,
      grid = false,
      xlabel = "Time (min)",
      ylabel = "Temperature (°C)",
      xlims = (0, 200),
      ylims = (-200, 50),
      xticks = 0:25:200,
      yticks = -200:25:50,
      minorticks = true)

savefig("mesh_sensitivity.svg")

max_stars_repo_name: NikStoyanov/NeuralThermalODE

max_issues_repo_name: NikStoyanov/NeuralThermalODE

max_forks_repo_name: NikStoyanov/NeuralThermalODE


text:
# Given some model parameters, synthesize and plot the channel maps and
# integrated spectrum

using ArgParse

s = ArgParseSettings()
@add_arg_table s begin
    # "--opt1"
    # help = "an option with an argument"
    # default = 0
    "--norad"
    help = "Use the image alread here."
    action = :store_true
    "config"
    help = "a YAML configuration file"
    required = true
end

parsed_args = parse_args(ARGS, s)

import YAML
config = YAML.load(open(parsed_args["config"]))

using constants
using image
using model
using HDF5

# import PyPlot.plt
# using LaTeXStrings

pp = config["parameters"]
params = ["M_star", "r_c", "T_10", "q", "gamma", "logM_CO", "ksi", "dpc", "incl", "PA", "vel", "mu_RA", "mu_DEC"]
nparam = length(params)
starting_param = Array(Float64, nparam)

for i=1:nparam
    starting_param[i] = pp[params[i]][1]
end

# Convert logM_CO to M_CO
starting_param[6] = 10^starting_param[6]

pars = Parameters(starting_param...)

grd = config["grid"]
grid = Grid(grd["nr"], grd["ntheta"], grd["r_in"], grd["r_out"], true)

# read the wavelengths for all 23 channels
fid = h5open(config["data_file"], "r")
lams = read(fid["lams"]) # [μm]
close(fid)

vel = pars.vel # [km/s]
# RADMC conventions for inclination and PA
incl = pars.incl # [deg]
PA = pars.PA # [deg] Position angle runs counter clockwise
npix = config["npix"] # number of pixels

# Doppler shift the dataset wavelength to rest-frame wavelength
beta = vel/c_kms # relativistic Doppler formula
shift_lams =  lams .* sqrt((1. - beta) / (1. + beta)) # [microns]

lam0 = cc/230.538e9 * 1e4 # [microns]

write_dust(pars, "", grid)
write_grid("", grid)
write_model(pars, "", grid)
write_lambda(shift_lams)

# Temporarily overwrite to dust mode
cp("radmc3d.inp.dust", "radmc3d.inp")

if !parsed_args["norad"]
    run(`radmc3d mctherm`) # next command will do this automatically
    run(`radmc3d image incl $incl posang $PA npix $npix lambda $lam0`)
end

# When done, revert to gas only mode
cp("radmc3d.inp.gas", "radmc3d.inp")

max_stars_repo_name: elnjensen/DiskJockey

max_issues_repo_name: elnjensen/DiskJockey

max_forks_repo_name: elnjensen/DiskJockey


text:
# ---
# title: 137. Single Number II
# id: problem137
# author: Indigo
# date: 2021-06-03
# difficulty: Medium
# categories: Bit Manipulation
# link: <https://leetcode.com/problems/single-number-ii/description/>
# hidden: true
# ---
#
# Given an integer array `nums` where every element appears **three times**
# except for one, which appears **exactly once**. _Find the single element and
# return it_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [2,2,3,2]
#     Output: 3
#
#
# **Example 2:**
#
#
#
#     Input: nums = [0,1,0,1,0,1,99]
#     Output: 99
#
#
#
#
# **Constraints:**
#
#   * `1 <= nums.length <= 3 * 104`
#   * `-231 <= nums[i] <= 231 - 1`
#   * Each element in `nums` appears exactly **three times** except for one element which appears **once**.
#
#
#
# **Follow up:**  Your algorithm should have a linear runtime complexity. Could
# you implement it without using extra memory?
#
#
## @lc code=start
using LeetCode

function single_number(nums::Vector{Int})::Int
    a, b = 0, 0
    for num in nums
        b = (b ⊻ num) & ~a
        a = (a ⊻ num) & ~b
    end
    return b
end
## @lc code=end


text:
# Gtk linked list #works

## Type hierarchy information
#compiles 100%
### an _LList is expected to have a data::Ptr{T} and next::Ptr{_LList{T}} element
### they are expected to be allocated and freed by GLib (e.g. with malloc/free)
import Base: popfirst!,Base.show,Base.size,Base.length,Base.pushfirst!,Base.push!
abstract type _LList{T} end

struct _GSList{T} <: _LList{T}
    data::Ptr{T}
    next::Ptr{_GSList{T}}
end
struct _GList{T} <: _LList{T}
    data::Ptr{T}
    next::Ptr{_GList{T}}
    prev::Ptr{_GList{T}}
end
eltype(::Type{_LList{T}}) where {T} = T
eltype(::Type{L}) where {L <: _LList} = eltype(supertype(L))

mutable struct GList{L <: _LList, T} <: AbstractVector{T}
    handle::Ptr{L}
    transfer_full::Bool
    function GList{L,T}(handle, transfer_full::Bool) where {L<:_LList,T}
        # if transfer_full == true, then also free the elements when finalizing the list
        # this function assumes the caller will take care of holding a pointer to the returned object
        # until it wants to be garbage collected
        @assert T == eltype(L)
        l = new{L,T}(handle, transfer_full)
        finalizer(empty!, l)
        return l
    end
end
GList(list::Type{T}) where {T} = GList(convert(Ptr{_GList{T}}, C_NULL), true)
GSList(list::Type{T}) where {T} = GList(convert(Ptr{_GSList{T}}, C_NULL), true)
GList(list::Ptr{L}, transfer_full::Bool = false) where {L <: _LList} = GList{L, eltype(L)}(list, transfer_full)

const  LList{L <: _LList} = Union{Ptr{L}, GList{L}}
eltype(::LList{L}) where {L <: _LList} = eltype(L)

_listdatatype(::Type{_LList{T}}) where {T} = T
_listdatatype(::Type{L}) where {L <: _LList} = _listdatatype(supertype(L))
deref(item::Ptr{L}) where {L <: _LList} = deref_to(L, unsafe_load(item).data) # extract something from the glist (automatically determine type)
deref_to(::Type{T}, x::Ptr) where {T} = unsafe_pointer_to_objref(x)::T # helper for extracting something from the glist (to type T)
deref_to(::Type{L}, x::Ptr) where {L <: _LList} = convert(eltype(L), deref_to(_listdatatype(L), x))
ref_to(::Type{T}, x) where {T} = gc_ref(x) # create a reference to something for putting in the glist
ref_to(::Type{L}, x) where {L <: _LList} = ref_to(_listdatatype(L), x)
empty!(li::Ptr{_LList}) = gc_unref(deref(li)) # delete an item in a glist
empty!(li::Ptr{L}) where {L <: _LList} = empty!(convert(Ptr{supertype(L)}, li))

## Standard Iteration protocol
start_(list::LList{L}) where {L} = unsafe_convert(Ptr{L}, list)
next_(::LList, s) = (deref(s), unsafe_load(s).next) # return (value, state)
done_(::LList, s) = (s == C_NULL)
iterate(list::LList, s=start_(list)) = done_(list, s) ? nothing : next_(list, s)

const  LListPair{L} = Tuple{LList, Ptr{L}}
function glist_iter(list::Ptr{L}, transfer_full::Bool = false) where L <: _LList
    # this function pairs every list element with the list head, to forestall garbage collection
    return (GList(list, transfer_full), list)
end
function next_(::LList, s::LListPair{L}) where L <: _LList
    return (deref(s[2]), (s[1], unsafe_load(s[2]).next))
end
done_(::LList, s::LListPair{L}) where {L <: _LList} = done_(s[1], s[2])

## Standard Array-like declarations
show(io::IO, ::MIME"text/plain", list::GList{L, T}) where {L, T} = show(io, list)
show(io::IO, list::GList{L, T}) where {L, T} = print(io, "GList{$L => $T}(length = $(length(list)), transfer_full = $(list.transfer_full))")

unsafe_convert(::Type{Ptr{L}}, list::GList) where {L <: _LList} = list.handle
size(list::LList) = (length(list),)
isempty(list::LList{L}) where {L} = (unsafe_convert(Ptr{L}, list) == C_NULL)

popfirst!(list::GList) = splice!(list, nth_first(list))
pop!(list::GList) = splice!(list, nth_last(list))
deleteat!(list::GList, i::Integer) = deleteat!(list, nth(list, i))

function splice!(list::GList, item::Ptr)
    x = deref(item)
    deleteat!(list, item)
    x
end

setindex!(list::GList, x, i::Real) = setindex!(list, x, nth(list, i))

## More Array-like declarations, this time involving ccall

### Non-modifying functions
length(list::LList{L}) where {L <: _GSList} = Int(ccall((:g_slist_length, libglib), Cuint, (Ptr{L},), list))
length(list::LList{L}) where {L <: _GList} = Int(ccall((:g_list_length, libglib), Cuint, (Ptr{L},), list))
copy(list::GList{L}) where {L <: _GSList} = typeof(list)(ccall((:g_slist_copy, libglib), Ptr{L}, (Ptr{L},), list), false)
copy(list::GList{L}) where {L <: _GList} = typeof(list)(ccall((:g_list_copy, libglib), Ptr{L}, (Ptr{L},), list), false)
check_undefref(p::Ptr) = (p == C_NULL ? error(UndefRefError()) : p)
nth_first(list::LList{L}) where {L <: _GSList} = unsafe_convert(Ptr{L}, list)
nth_first(list::LList{L}) where {L <: _GList} =
    check_undefref(ccall((:g_list_first, libglib), Ptr{L}, (Ptr{L},), list))
nth_last(list::LList{L}) where {L <: _GSList} =
    check_undefref(ccall((:g_slist_last, libglib), Ptr{L}, (Ptr{L},), list))
nth_last(list::LList{L}) where {L <: _GList} =
    check_undefref(ccall((:g_list_last, libglib), Ptr{L}, (Ptr{L},), list))
nth(list::LList{L}, i::Integer) where {L <: _GSList} =
    check_undefref(ccall((:g_slist_nth, libglib), Ptr{L}, (Ptr{L}, Cuint), list, i - 1))
nth(list::LList{L}, i::Integer) where {L <: _GList} =
    check_undefref(ccall((:g_list_nth, libglib), Ptr{L}, (Ptr{L}, Cuint), list, i - 1))
function getindex(list::LList{_GSList{T}}, i::Integer) where T
    p = check_undefref(ccall((:g_slist_nth_data, libglib), Ptr{T}, (Ptr{_GSList{T}}, Cuint), list, i - 1))
    return deref_to(_GSList{T}, p)
end
function getindex(list::LList{_GList{T}}, i::Integer) where T
    p = check_undefref(ccall((:g_list_nth_data, libglib), Ptr{T}, (Ptr{_GList{T}}, Cuint), list, i - 1))
    return deref_to(_GList{T}, p)
end
function get(list::LList{_GSList{T}}, i::Integer, default) where T
    p = ccall((:g_slist_nth_data, libglib), Ptr{T}, (Ptr{_GSList{T}}, Cuint), list, i - 1)
    p == C_NULL && return default
    return deref_to(_GSList{T}, p)
end
function get(list::LList{_GList{T}}, i::Integer, default) where T
    p = ccall((:g_list_nth_data, libglib), Ptr{T}, (Ptr{_GList{T}}, Cuint), list, i - 1)
    p == C_NULL && return default
    return deref_to(_GList{T}, p)
end

### Modifying functions (!) are only allowed on a GList
function empty!(list::GList{L}) where L <: _GSList
    if list.handle != C_NULL
        if list.transfer_full
            s = start_(list)
            while !done_(list, s)
                empty!(s)
                s = next_(list, s)[2]
            end
        end
        ccall((:g_slist_free, libglib), Nothing, (Ptr{L},), list)
        list.handle = C_NULL
    end
    return list
end
function empty!(list::GList{L}) where L <: _GList
    if list.handle != C_NULL
        if list.transfer_full
            s = start_(list)
            while !done_(list, s)
                empty!(s)
                s = next_(list, s)[2]
            end
        end
        ccall((:g_list_free, libglib), Nothing, (Ptr{L},), list)
        list.handle = C_NULL
    end
    return list
end
function append!(l1::GList{L}, l2::GList{L}) where L <: _GSList
    (l1.transfer_full | l2.transfer_full) && error("cannot combine lists with transfer_full = true")
    l1.handle = ccall((:g_slist_concat, libglib), Ptr{L}, (Ptr{L}, Ptr{L}), l1, l2)
    l2.handle = C_NULL
    return l1
end
function append!(l1::GList{L}, l2::GList{L}) where L <: _GList
    (l1.transfer_full | l2.transfer_full) && error("cannot combine lists with transfer_full = true")
    l1.handle = ccall((:g_list_concat, libglib), Ptr{L}, (Ptr{L}, Ptr{L}), l1, l2)
    l2.handle = C_NULL
    return l1
end
function reverse!(list::GList{L}) where L <: _GSList
    list.handle = ccall((:g_slist_reverse, libglib), Ptr{L}, (Ptr{L},), list)
    return list
end
function reverse!(list::GList{L}) where L <: _GList
    list.handle = ccall((:g_list_reverse, libglib), Ptr{L}, (Ptr{L},), list)
    return list
end
function insert!(list::GList{_GSList{T}}, i::Integer, item) where T
    list.handle = ccall((:g_slist_insert, libglib), Ptr{_GSList{T}},
        (Ptr{_GSList{T}}, Ptr{T}, Cint),
        list, ref_to(_GSList{T}, item), i - 1)
    return list
end
function insert!(list::GList{_GList{T}}, i::Integer, item) where T
    list.handle = ccall((:g_list_insert, libglib), Ptr{_GList{T}},
        (Ptr{_GList{T}}, Ptr{T}, Cint),
        list, ref_to(_GList{T}, item), i - 1)
    return list
end
function insert!(list::GList{_GSList{T}}, i::Ptr{_GSList{T}}, item) where T
    list.handle = ccall((:g_slist_insert_before, libglib), Ptr{_GSList{T}},
        (Ptr{_GSList{T}}, Ptr{_GSList{T}}, Ptr{T}),
        list, i, ref_to(_GSList{T}, item))
    return list
end
function insert!(list::GList{_GList{T}}, i::Ptr{_GList{T}}, item) where T
    list.handle = ccall((:g_list_insert_before, libglib), Ptr{_GList{T}},
        (Ptr{_GList{T}}, Ptr{_GList{T}}, Ptr{T}),
        list, i, ref_to(_GList{T}, item))
    return list
end
function pushfirst!(list::GList{_GSList{T}}, item) where T
    list.handle = ccall((:g_slist_prepend, libglib), Ptr{_GSList{T}}, (Ptr{_GSList{T}}, Ptr{T}), list, ref_to(_GSList{T}, item))
    return list
end
function pushfirst!(list::GList{_GList{T}}, item) where T
    list.handle = ccall((:g_list_prepend, libglib), Ptr{_GList{T}},
        (Ptr{_GList{T}}, Ptr{T}),
        list, ref_to(_GList{T}, item))
    return list
end
function push!(list::GList{_GSList{T}}, item) where T
    list.handle = ccall((:g_slist_append, libglib), Ptr{_GSList{T}},
        (Ptr{_GSList{T}}, Ptr{T}),
        list, ref_to(_GSList{T}, item))
    return list
end
function push!(list::GList{_GList{T}}, item) where T
    list.handle = ccall((:g_list_append, libglib), Ptr{_GList{T}},
        (Ptr{_GList{T}}, Ptr{T}),
        list, ref_to(_GList{T}, item))
    return list
end
function deleteat!(list::GList{L}, i::Ptr{L}) where L <: _GSList
    list.transfer_full && empty!(i)
    list.handle = ccall((:g_slist_delete_link, libglib), Ptr{L}, (Ptr{L}, Ptr{L}), list, i)
    return list
end
function deleteat!(list::GList{L}, i::Ptr{L}) where L <: _GList
    list.transfer_full && empty!(i)
    list.handle = ccall((:g_list_delete_link, libglib), Ptr{L}, (Ptr{L}, Ptr{L}), list, i)
    return list
end
function setindex!(list::GList{L}, item, i::Ptr{L}) where L <: _GSList
    list.transfer_full && empty!(i)
    idx = unsafe_load(i)
    idx = L(ref_to(L, item), idx.next)
    unsafe_store!(i, idx)
    return list
end
function setindex!(list::GList{L}, item, i::Ptr{L}) where L <: _GList
    list.transfer_full && empty!(i)
    idx = unsafe_load(i)
    idx = L(ref_to(L, item), idx.next, idx.prev)
    unsafe_store!(i, idx)
    return list
end

### Store most pointers without doing anything special
ref_to(::Type{P}, x) where {P <: Ptr} = x
deref_to(::Type{P}, x::Ptr) where {P <: Ptr} = x
empty!(li::Ptr{_LList{P}}) where {P <: Ptr} = nothing

### Store numbers directly inside the pointer bits (assuming convert(N, x) exists)
ref_to(::Type{N}, x) where {N <: Number} = Ptr{N}(x)
deref_to(::Type{N}, x::Ptr) where {N <: Number} = x
empty!(li::Ptr{_LList{N}}) where {N <: Number} = nothing

### Handle storing pointers to numbers
eltype(::Type{_LList{Ptr{N}}}) where {N <: Number} = N
deref_to(::Type{Ptr{N}}, p::Ptr) where {N <: Number} = unsafe_load(p)
ref_to(::Type{Ptr{N}}, x) where {N <: Number} = unsafe_store!(convert(Ptr{N}, g_malloc(N.size)), x)
empty!(li::Ptr{_GSList{Ptr{N}}}) where {N <: Number} = g_free(unsafe_load(li).data)
empty!(li::Ptr{_GList{Ptr{N}}}) where {N <: Number} = g_free(unsafe_load(li).data)

### Store (byte)strings as pointers
function deref_to(::Type{S}, p::Ptr) where {S <: String}
    p==C_NULL && return ""
    bytestring(convert(Ptr{UInt8}, p))
end
function ref_to(::Type{S}, x) where S <: String
    s = bytestring(x)
    p = ccall((:g_strdup, libglib), Ptr{UInt8}, (Cstring,), s)
    return p
end

function empty!(li::Ptr{_GSList{S}}) where {S <: String}
    gl=unsafe_load(li)
    g_free(gl.data)
    gl = _GSList{S}(C_NULL, gl.next) # set to NULL to prevent Julia from trying to access this while finalizing
    unsafe_store!(li, gl)
end

function empty!(li::Ptr{_GList{S}}) where {S <: String}
    gl=unsafe_load(li)
    g_free(gl.data)
    gl = _GList{S}(C_NULL, gl.next, gl.prev) # set to NULL to prevent Julia from trying to access this while finalizing
    unsafe_store!(li, gl)
end

max_stars_repo_name: adamwillisMastery/MiniGtk

max_issues_repo_name: adamwillisMastery/MiniGtk

max_forks_repo_name: adamwillisMastery/MiniGtk


text:
module Helmholtz
using LinearAlgebra
using SparseArrays
using jInv.Mesh;
using jInv.LinearSolvers
import jInv.Utils.clear!
import jInv.LinearSolvers.AbstractSolver
import jInv.LinearSolvers.solveLinearSystem

using Multigrid;

export HelmholtzParam,getShiftedHelmholtzParam
mutable struct HelmholtzParam
	Mesh   			:: RegularMesh;
	gamma  			:: Array{Float64};
	m               :: Array{Float64};
	omega			:: Union{Float64,ComplexF64}
	NeumannOnTop	:: Bool
	Sommerfeld 		:: Bool
end

import jInv.Utils.clear!

function clear!(HP::HelmholtzParam)
	clear!(HP.Mesh);
	m = zeros(0);
	gamma = zeros(0);
	omega = 0.0;
end

function getShiftedHelmholtzParam(p::HelmholtzParam,s::Float64)
	return HelmholtzParam(p.Mesh,p.gamma + s*real(p.omega), p.m,p.omega,p.NeumannOnTop,p.Sommerfeld);
end

include("GetHelmholtz.jl");

include("PlainNodalLaplacian.jl");
include("ShiftedLaplacianMultigridSolver.jl");
include("PointSourceADR/PointSourceADR.jl");
include("getPointSource.jl")
#include("Elastic/ElasticHelmholtz.jl");

end # module


text:
to_func_name(x::Symbol) = Symbol(lowercase(string(x)))

"""
     default_plot_signatures(funcname, funcname!, PlotType)
Creates all the different overloads for `funcname` that need to be supported for the plotting frontend!
Since we add all these signatures to different functions, we make it reusable with this function.
The `Core.@__doc__` macro transfers the docstring given to the Recipe into the functions.
"""
function default_plot_signatures(funcname, funcname!, PlotType)
    quote
        Core.@__doc__ function ($funcname)(args...; attributes...)
            plot($PlotType, args...; attributes...)
        end

        Core.@__doc__ function ($funcname!)(args...; attributes...)
            plot!($PlotType, args...; attributes...)
        end
    end
end

"""
# Plot Recipes in `AbstractPlotting`

There's two types of recipes. *Type recipes* define a simple mapping from a
user defined type to an existing plot type. *Full recipes* can customize the
theme and define a custom plotting function.

## Type recipes

Type recipe are really simple and just overload the argument conversion
pipeline. This can be done for all plot types or for a subset of plot types:

    # All plot types
    convert_arguments(P::Type{<:AbstractPlot}, x::MyType) = convert_arguments(P, rand(10, 10))
    # Only for scatter plots
    convert_arguments(P::Type{<:Scatter}, x::MyType) = convert_arguments(P, rand(10, 10))

Optionally you may define the default plot type so that `plot(x::MyType)` will
use this:

    plottype(::MyType) = Surface

## Full recipes with the `@recipe` macro

A full recipe for `MyPlot` comes in two parts. First is the plot type name,
arguments and theme definition which are defined using the `@recipe` macro.
Second is a custom `plot!` for `MyPlot`, implemented in terms of the atomic
plotting functions.

We use an example to show how this works:

    # arguments (x, y, z) && theme are optional
    @recipe(MyPlot, x, y, z) do scene
        Attributes(
            plot_color => :red
        )
    end

This macro expands to several things. Firstly a type definition:

    const MyPlot{ArgTypes} = Combined{myplot, ArgTypes}

The type parameter of `Combined` contains the function instead of e.g. a
symbol. This way the mapping from `MyPlot` to `myplot` is safer and simpler.
(The downside is we always need a function `myplot` - TODO: is this a problem?)

The following signatures are defined to make `MyPlot` nice to use:

    myplot(args...; kw_args...) = ...
    myplot!(scene, args...; kw_args...) = ...
    myplot(kw_args::Dict, args...) = ...
    myplot!(scene, kw_args::Dict, args...) = ...
    #etc (not 100% settled what signatures there will be)

A specialization of `argument_names` is emitted if you have an argument list
`(x,y,z)` provided to the recipe macro:

    argument_names(::Type{<: MyPlot}) = (:x, :y, :z)

This is optional but it will allow the use of `plot_object[:x]` to
fetch the first argument from the call
`plot_object = myplot(rand(10), rand(10), rand(10))`, for example.
Alternatively you can always fetch the `i`th argument using `plot_object[i]`,
and if you leave out the `(x,y,z)`, the default version of `argument_names`
will provide `plot_object[:arg1]` etc.

The theme given in the body of the `@recipe` invocation is inserted into a
specialization of `default_theme` which inserts the theme into any scene that
plots `MyPlot`:

    function default_theme(scene, ::MyPlot)
        Attributes(
            plot_color => :red
        )
    end

As the second part of defining `MyPlot`, you should implement the actual
plotting of the `MyPlot` object by specializing `plot!`:

    function plot!(plot::MyPlot)
        # normal plotting code, building on any previously defined recipes
        # or atomic plotting operations, and adding to the combined `plot`:
        lines!(plot, rand(10), color = plot[:plot_color])
        plot!(plot, plot[:x], plot[:y])
        plot
    end

It's possible to add specializations here, depending on the argument *types*
supplied to `myplot`. For example, to specialize the behavior of `myplot(a)`
when `a` is a 3D array of floating point numbers:

    const MyVolume = MyPlot{Tuple{<:AbstractArray{<: AbstractFloat, 3}}}
    argument_names(::Type{<: MyVolume}) = (:volume,) # again, optional
    function plot!(plot::MyVolume)
        # plot a volume with a colormap going from fully transparent to plot_color
        volume!(plot, plot[:volume], colormap = :transparent => plot[:plot_color])
        plot
    end

The docstring given to the recipe will be transferred to the functions it generates.

"""
macro recipe(theme_func, Tsym::Symbol, args::Symbol...)
    funcname_sym = to_func_name(Tsym)
    funcname! = esc(Symbol("$(funcname_sym)!"))
    PlotType = esc(Tsym)
    funcname = esc(funcname_sym)
    expr = quote
        $(funcname)() = not_implemented_for($funcname)
        const $(PlotType){$(esc(:ArgType))} = Combined{$funcname, $(esc(:ArgType))}
        Base.show(io::IO, ::Type{<: $PlotType}) = print(io, $(string(Tsym)), "{...}")
        $(default_plot_signatures(funcname, funcname!, PlotType))
        AbstractPlotting.default_theme(scene, ::Type{<: $PlotType}) = $(esc(theme_func))(scene)
        export $PlotType, $funcname, $funcname!
    end
    if !isempty(args)
        push!(expr.args, :($(esc(:(AbstractPlotting.argument_names)))(::Type{<: $PlotType}, len::Integer) = $args))
    end
    expr
end


text:
"""
A `Worker` object holds the necessary local states for a learner worker:

 - dset_tr: the local partition of the training dataset
 - dset_tt: the local partition of the test dataset
 - hp: hyper parameters specific to each algorithm (e.g. batch size)
 - lp: local parameters specific to each algorithm (e.g. algorithms in the
   dual might store dual variables locally)
"""
mutable struct Worker{HPType, LPType}
  dset_tr   :: Dataset # training set
  dset_tt   :: Dataset # test / validation set

  hp        :: HPType  # hyper parameters
  lp        :: LPType  # worker local parameters

  Worker{HPType, LPType}() where {HPType, LPType} = new()
end

"""
    spanw_workers(ctor, args...)

Create `Worker` objects on each worker process. Return a list
of (pid, worker future ref) tuples.
"""
function spawn_workers(ctor :: Function, args...)
  return map(1:nworkers()) do wid
    pid = wid + 1  # head node is 1, all other nodes are workers
    (pid, remotecall(ctor, pid, wid, nworkers(), args...))
  end
end

"""
    invoke_on_workers(f, workers, args...)

Invoke `f` on each workers in parallel, with `f(worker_ref, args...)`.
This function blocks until all workers finishes, and return a list
of results from each worker.
"""
function invoke_on_workers(f :: Function, workers, args...)
  rets = Array{Any}(length(workers))
  @sync begin
    for (i, (pid, worker_ref)) in enumerate(workers)
      @async begin
        rets[i] = remotecall_fetch(f, pid, worker_ref, args...)
      end
    end
  end

  # propagate remote exceptions
  for obj in rets
    if obj isa RemoteException
      throw(obj)
    end
  end
  return rets
end

"""
    accumulate_metric(workers, Type{MetricType}, w)

Collect and accumulate metrics from workers. This function assumes linear
predictors and assume only one metric object is needed. One can define this
in more general ways for general predictors and multiple metrics. Also it
might be more efficient to collect the metrics during the forward stage of
the training to avoid duplicated computation.
"""
function accumulate_metric{T <: Metric}(workers, :: Type{T}, w :: Vector{Float64})
  function eval_metric{T <: Metric}(worker_ref :: Future, :: Type{T}, w :: Vector{Float64})
    worker = fetch(worker_ref)
    return map((worker.dset_tr, worker.dset_tt)) do dset
      pred = dset.data' * w
      evaluate(T, pred, dset.labels)
    end
  end
  # rets is [(m_tr1, m_tt1), (m_tr2, m_tt2), ...]
  rets = invoke_on_workers(eval_metric, workers, T, w)
  return merge([x[1] for x in rets]), merge([x[2] for x in rets])
end


text:
module PandorasBoxTimefield

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/timefield" Timefield(x::Integer, y::Integer, start::Number=0.2, stop::Number=1.0, stopTime::Number=1.0, startTime::Number=3.0, animRate::Number=6.0, render::Bool=true, lingering::Bool=false, color::String="Teal", respectOtherTimeRates::Bool=false)

const placements = Ahorn.PlacementDict(
    "Timefield (Pandora's Box)" => Ahorn.EntityPlacement(
        Timefield,
        "rectangle"
    ),
)

const colors = sort(collect(keys(Ahorn.XNAColors.colors)))

Ahorn.editingOptions(entity::Timefield) = Dict{String, Any}(
    "color" => colors
)

Ahorn.minimumSize(entity::Timefield) = 8, 8
Ahorn.resizable(entity::Timefield) = true, true

Ahorn.selection(entity::Timefield) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::Timefield, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 8))
    height = Int(get(entity.data, "height", 8))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.5, 1.0, 1.0, 0.4), (0.5, 1.0, 1.0, 1.0))
end

end

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
module Radios

using Genie, Stipple, StippleUI, StippleUI.API
import Genie.Renderer.Html: HTMLString, normal_element, template

export radio

function __init__()
  Genie.Renderer.Html.register_normal_element("q__radio", context = Genie.Renderer.Html)
end

function radio( label::AbstractString = "",
                fieldname::Union{Symbol,Nothing} = nothing,
                args...;
                wrap::Function = StippleUI.DEFAULT_WRAPPER,
                kwargs...)
  wrap() do
    Genie.Renderer.Html.q__radio(args...; attributes([:label => label, :fieldname => fieldname, kwargs...], StippleUI.API.ATTRIBUTES_MAPPINGS)...)
  end
end

end


text:
using NeuralOperators
using Documenter
using DocumenterCitations

bib = CitationBibliography(joinpath(@__DIR__, "bibliography.bib"), sorting=:nyt)

DocMeta.setdocmeta!(NeuralOperators, :DocTestSetup, :(using NeuralOperators); recursive=true)

makedocs(
    bib,
    modules=[NeuralOperators],
    authors="JingYu Ning <foldfelis@gmail.com> and contributors",
    repo="https://github.com/SciML/NeuralOperators.jl/blob/{commit}{path}#{line}",
    sitename="NeuralOperators.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="http://neuraloperators.sciml.ai",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Introduction" => "introduction.md",
        "APIs" => "apis.md",
        "References" => "references.md",
    ],
)

deploydocs(;
    repo="github.com/SciML/NeuralOperators.jl",
)


text:
using DiscgolfRecord
using Test

@testset "DiscgolfRecord.jl" begin
    # Write your tests here.

    # Make sure we can preview courses
    @test_nowarn preview_course(COURSES["kit_carson"])
    @test_nowarn preview_course(COURSES["mast_park"])

    # Test the round reader
    sample_round_file = "20210719_mastpark.csv"
    @test_nowarn rnd = DiscgolfRecord.read_round(sample_round_file)

end