text:
export FiniteDifferenceMethod, fdm, backward_fdm, forward_fdm, central_fdm, extrapolate_fdm

"""
    estimate_magitude(f, x::T) where T<:AbstractFloat

Estimate the magnitude of `f` in a neighbourhood of `x`, assuming that the outputs of `f`
have a "typical" order of magnitude. The result should be interpreted as a very rough
estimate. This function deals with the case that `f(x) = 0`.
"""
function estimate_magitude(f, x::T) where T<:AbstractFloat
    M = float(maximum(abs, f(x)))
    M > 0 && (return M)
    # Ouch, `f(x) = 0`. But it may not be zero around `x`. We conclude that `x` is likely a
    # pathological input for `f`. Perturb `x`. Assume that the pertubed value for `x` is
    # highly unlikely also a pathological value for `f`.
    Δ = convert(T, 0.1) * max(abs(x), one(x))
    return float(maximum(abs, f(x + Δ)))
end

"""
    estimate_roundoff_error(f, x::T) where T<:AbstractFloat

Estimate the round-off error of `f(x)`. This function deals with the case that `f(x) = 0`.
"""
function estimate_roundoff_error(f, x::T) where T<:AbstractFloat
    # Estimate the round-off error. It can happen that the function is zero around `x`, in
    # which case we cannot take `eps(f(x))`. Therefore, we assume a lower bound that is
    # equal to `eps(T) / 1000`, which gives `f` four orders of magnitude wiggle room.
    return max(eps(estimate_magitude(f, x)), eps(T) / 1000)
end

"""
    FiniteDifferences.DEFAULT_CONDITION

The default [condition number](https://en.wikipedia.org/wiki/Condition_number) used when
computing bounds. It provides amplification of the ∞-norm when passed to the function's
derivatives.
"""
const DEFAULT_CONDITION = 100

"""
    FiniteDifferenceMethod{G<:AbstractVector, C<:AbstractVector, E<:Function}

A finite difference method.

# Fields
- `grid::G`: Multiples of the step size that the function will be evaluated at.
- `q::Int`: Order of derivative to estimate.
- `coefs::C`: Coefficients corresponding to the grid functions that the function evaluations
    will be weighted by.
- `bound_estimator::Function`: A function that takes in the function and the evaluation
    point and returns a bound on the magnitude of the `length(grid)`th derivative.
"""
struct FiniteDifferenceMethod{G<:AbstractVector, C<:AbstractVector, E<:Function}
    grid::G
    q::Int
    coefs::C
    bound_estimator::E
end

"""
    FiniteDifferenceMethod(
        grid::AbstractVector,
        q::Int;
        condition::Real=DEFAULT_CONDITION
    )

Construct a finite difference method.

# Arguments
- `grid::Abstract`: The grid. See [`FiniteDifferenceMethod`](@ref).
- `q::Int`: Order of the derivative to estimate.
- `condition::Real`: Condition number. See [`DEFAULT_CONDITION`](@ref).

# Returns
- `FiniteDifferenceMethod`: Specified finite difference method.
"""
function FiniteDifferenceMethod(
    grid::AbstractVector,
    q::Int;
    condition::Real=DEFAULT_CONDITION
)
    p = length(grid)
    _check_p_q(p, q)
    return FiniteDifferenceMethod(
        grid,
        q,
        _coefs(grid, q),
        _make_default_bound_estimator(condition=condition)
    )
end

"""
    (m::FiniteDifferenceMethod)(
        f::Function,
        x::T;
        factor::Real=1,
        max_step::Real=0.1 * max(abs(x), one(x))
    ) where T<:AbstractFloat

Estimate the derivative of `f` at `x` using the finite differencing method `m` and an
automatically determined step size.

# Arguments
- `f::Function`: Function to estimate derivative of.
- `x::T`: Input to estimate derivative at.

# Keywords
- `factor::Real=1`: Factor to amplify the estimated round-off error by. This can be used
    to force a more conservative step size.
- `max_step::Real=0.1 * max(abs(x), one(x))`: Maximum step size.

# Returns
- Estimate of the derivative.

# Examples

```julia-repl
julia> fdm = central_fdm(5, 1)
FiniteDifferenceMethod:
  order of method:       5
  order of derivative:   1
  grid:                  [-2, -1, 0, 1, 2]
  coefficients:          [0.08333333333333333, -0.6666666666666666, 0.0, 0.6666666666666666, -0.08333333333333333]

julia> fdm(sin, 1)
0.5403023058681155

julia> fdm(sin, 1) - cos(1)  # Check the error.
-2.4313884239290928e-14

julia> FiniteDifferences.estimate_step(fdm, sin, 1.0)  # Computes step size and estimates the error.
(0.0010632902144695163, 1.9577610541734626e-13)
```
"""
@inline function (m::FiniteDifferenceMethod)(f::Function, x::Real; kw_args...)
    # Assume that converting to float is desired.
    return _call_method(m, f, float(x); kw_args...)
end
@inline function _call_method(
    m::FiniteDifferenceMethod,
    f::Function,
    x::T;
    factor::Real=1,
    max_step::Real=0.1 * max(abs(x), one(x))
) where T<:AbstractFloat
    # The automatic step size calculation fails if `m.q == 0`, so handle that edge case.
    iszero(m.q) && return f(x)
    h, _ = estimate_step(m, f, x, factor=factor, max_step=max_step)
    return _eval_method(m, f, x, h)
end

"""
    (m::FiniteDifferenceMethod)(f::Function, x::T, h::Real) where T<:AbstractFloat

Estimate the derivative of `f` at `x` using the finite differencing method `m` and a given
step size.

# Arguments
- `f::Function`: Function to estimate derivative of.
- `x::T`: Input to estimate derivative at.
- `h::Real`: Step size.

# Returns
- Estimate of the derivative.

# Examples

```julia-repl
julia> fdm = central_fdm(5, 1)
FiniteDifferenceMethod:
  order of method:       5
  order of derivative:   1
  grid:                  [-2, -1, 0, 1, 2]
  coefficients:          [0.08333333333333333, -0.6666666666666666, 0.0, 0.6666666666666666, -0.08333333333333333]

julia> fdm(sin, 1, 1e-3)
0.5403023058679624

julia> fdm(sin, 1, 1e-3) - cos(1)  # Check the error.
-1.7741363933510002e-13
```
"""
@inline function (m::FiniteDifferenceMethod)(f::Function, x::Real, h::Real)
    # Assume that converting to float is desired.
    return _eval_method(m, f, float(x), h)
end
@inline function _eval_method(
    m::FiniteDifferenceMethod,
    f::Function,
    x::T,
    h::Real
) where T<:AbstractFloat
    return sum(
        i -> convert(T, m.coefs[i]) * f(T(x + h * m.grid[i])),
        eachindex(m.grid)
    ) / h^m.q
end

# Check the method and derivative orders for consistency.
function _check_p_q(p::Integer, q::Integer)
    q >= 0 || throw(DomainError(q, "order of derivative (`q`) must be non-negative"))
    q < p || throw(DomainError(
        (q, p),
        "order of the method (q) must be strictly greater than that of the derivative (p)",
    ))
    # Check whether the method can be computed. We require the factorial of the method order
    # to be computable with regular `Int`s, but `factorial` will after 20, so 20 is the
    # largest we can allow.
    p > 20 && throw(DomainError(p, "order of the method (`p`) is too large to be computed"))
    return
end

const _COEFFS_CACHE = Dict{Tuple{AbstractVector{<:Real}, Integer}, Vector{Float64}}()

# Compute coefficients for the method and cache the result.
function _coefs(grid::AbstractVector{<:Real}, q::Integer)
    return get!(_COEFFS_CACHE, (grid, q)) do
        p = length(grid)
        # For high precision on the `\`, we use `Rational`, and to prevent overflows we use
        # `Int128`. At the end we go to `Float64` for fast floating point math, rather than
        # rational math.
        C = [Rational{Int128}(g^i) for i in 0:(p - 1), g in grid]
        x = zeros(Rational{Int128}, p)
        x[q + 1] = factorial(q)
        return Float64.(C \ x)
    end
end

# Estimate the bound on the derivative by amplifying the ∞-norm.
function _make_default_bound_estimator(; condition::Real=DEFAULT_CONDITION)
    default_bound_estimator(f, x) = condition * estimate_magitude(f, x)
    return default_bound_estimator
end

function Base.show(io::IO, m::MIME"text/plain", x::FiniteDifferenceMethod)
    @printf io "FiniteDifferenceMethod:\n"
    @printf io "  order of method:       %d\n" length(x.grid)
    @printf io "  order of derivative:   %d\n" x.q
    @printf io "  grid:                  %s\n" x.grid
    @printf io "  coefficients:          %s\n" x.coefs
end

"""
    function estimate_step(
        m::FiniteDifferenceMethod,
        f::Function,
        x::T;
        factor::Real=1,
        max_step::Real=0.1 * max(abs(x), one(x))
    ) where T<:AbstractFloat

Estimate the step size for a finite difference method `m`. Also estimates the error of the
estimate of the derivative.

# Arguments
- `m::FiniteDifferenceMethod`: Finite difference method to estimate the step size for.
- `f::Function`: Function to evaluate the derivative of.
- `x::T`: Point to estimate the derivative at.

# Keywords
- `factor::Real=1`. Factor to amplify the estimated round-off error by. This can be used
    to force a more conservative step size.
- `max_step::Real=0.1 * max(abs(x), one(x))`: Maximum step size.

# Returns
- `Tuple{T, <:AbstractFloat}`: Estimated step size and an estimate of the error of the
    finite difference estimate.
"""
function estimate_step(
    m::FiniteDifferenceMethod,
    f::Function,
    x::T;
    factor::Real=1,
    max_step::Real=0.1 * max(abs(x), one(x))
) where T<:AbstractFloat
    p = length(m.coefs)
    q = m.q

    # Estimate the round-off error.
    ε = estimate_roundoff_error(f, x) * factor

    # Estimate the bound on the derivatives.
    M = m.bound_estimator(f, x)

    # Set the step size by minimising an upper bound on the error of the estimate.
    C₁ = ε * sum(abs, m.coefs)
    C₂ = M * sum(n -> abs(m.coefs[n] * m.grid[n]^p), eachindex(m.coefs)) / factorial(p)
    # Type inference fails on this, so we annotate it, which gives big performance benefits.
    h::T = convert(T, min((q / (p - q) * (C₁ / C₂))^(1 / p), max_step))

    # Estimate the accuracy of the method.
    accuracy = h^(-q) * C₁ + h^(p - q) * C₂

    return h, accuracy
end

for direction in [:forward, :central, :backward]
    fdm_fun = Symbol(direction, "_fdm")
    grid_fun = Symbol("_", direction, "_grid")
    @eval begin function $fdm_fun(
            p::Int,
            q::Int;
            adapt::Int=1,
            condition::Real=DEFAULT_CONDITION,
            geom::Bool=false
        )
            _check_p_q(p, q)
            grid = $grid_fun(p)
            geom && (grid = _exponentiate_grid(grid))
            coefs = _coefs(grid, q)
            return FiniteDifferenceMethod(
                grid,
                q,
                coefs,
                _make_adaptive_bound_estimator($fdm_fun, p, q, adapt, condition, geom=geom),
            )
        end

        @doc """
    $($(Meta.quot(fdm_fun)))(
        p::Int,
        q::Int;
        adapt::Int=1,
        condition::Real=DEFAULT_CONDITION,
        geom::Bool=false
    )

Contruct a finite difference method at a $($(Meta.quot(direction))) grid of `p` linearly
spaced points.

# Arguments
- `p::Int`: Number of grid points.
- `q::Int`: Order of the derivative to estimate.

# Keywords
- `adapt::Int=1`: Use another finite difference method to estimate the magnitude of the
    `p`th order derivative, which is important for the step size computation. Recurse
    this procedure `adapt` times.
- `condition::Real`: Condition number. See [`DEFAULT_CONDITION`](@ref).
- `geom::Bool`: Use geometrically spaced points instead of linearly spaced points.

# Returns
- `FiniteDifferenceMethod`: The specified finite difference method.
        """ $fdm_fun
    end
end

function _make_adaptive_bound_estimator(
    constructor::Function,
    p::Int,
    q::Int,
    adapt::Int,
    condition::Int;
    kw_args...
)
    if adapt >= 1
        estimate_derivative = constructor(
            p + 1, p, adapt=adapt - 1, condition=condition; kw_args...
        )
        return (f, x) -> estimate_magitude(x′ -> estimate_derivative(f, x′), x)
    else
        return _make_default_bound_estimator(condition=condition)
    end
end

_forward_grid(p::Int) = collect(0:(p - 1))

_backward_grid(p::Int) = collect((1 - p):0)

function _central_grid(p::Int)
    if isodd(p)
        return collect(div(1 - p, 2):div(p - 1, 2))
    else
        return vcat(div(-p, 2):-1, 1:div(p, 2))
    end
end

_exponentiate_grid(grid::Vector, base::Int=3) = sign.(grid) .* base .^ abs.(grid) ./ base

function _is_symmetric(vec::Vector; centre_zero::Bool=false, negate_half::Bool=false)
    half_sign = negate_half ? -1 : 1
    if isodd(length(vec))
        centre_zero && vec[end ÷ 2 + 1] != 0 && return false
        return vec[1:end ÷ 2] == half_sign .* reverse(vec[(end ÷ 2 + 2):end])
    else
        return vec[1:end ÷ 2] == half_sign .* reverse(vec[(end ÷ 2 + 1):end])
    end
end

function _is_symmetric(m::FiniteDifferenceMethod)
    grid_symmetric = _is_symmetric(m.grid, centre_zero=true, negate_half=true)
    coefs_symmetric =_is_symmetric(m.coefs, negate_half=true)
    return grid_symmetric && coefs_symmetric
end

"""
    extrapolate_fdm(
        m::FiniteDifferenceMethod,
        f::Function,
        x::T,
        h::Real=0.1 * max(abs(x), one(x));
        power=nothing,
        breaktol=Inf,
        kw_args...
    ) where T<:AbstractFloat

Use Richardson extrapolation to refine a finite difference method.

Takes further in keyword arguments for `Richardson.extrapolate`. This method
automatically sets `power = 2` if `m` is symmetric and `power = 1`. Moreover, it defaults
`breaktol = Inf`.

# Arguments
- `m::FiniteDifferenceMethod`: Finite difference method to estimate the step size for.
- `f::Function`: Function to evaluate the derivative of.
- `x::T`: Point to estimate the derivative at.
- `h::Real=0.1 * max(abs(x), one(x))`: Initial step size.

# Returns
- `Tuple{<:AbstractFloat, <:AbstractFloat}`: Estimate of the derivative and error.
"""
function extrapolate_fdm(
    m::FiniteDifferenceMethod,
    f::Function,
    x::T,
    h::Real=0.1 * max(abs(x), one(x));
    power::Int=1,
    breaktol::Real=Inf,
    kw_args...
) where T<:AbstractFloat
    (power == 1 && _is_symmetric(m)) && (power = 2)
    return extrapolate(h -> m(f, x, h), h; power=power, breaktol=breaktol, kw_args...)
end


text:
module TestDArray

using Dagger
using FLoops
using Folds
using FoldsDagger
using Test

function test_1d()
    @testset for n in 1:4
        A = distribute(1:4, Blocks(n))
        @test foldx_dagger(+, A) == sum(1:4)
        @test Folds.sum(A) == sum(1:4)
    end
end

function test_2d()
    @testset for n in 1:4, m in 1:3
        A = distribute(reshape(1:24, (4, 6)), Blocks(n, m))
        @test foldx_dagger(+, A) == sum(1:24)
        @test Folds.sum(A) == sum(1:24)
    end
end

end

ext: jl


text:
"""
Running the STG network
- For units to check out, synapses also get a conversion factor.
- The names of the variables in the system can be listed with `@show stg.states`.
- You can plot variables by their names with plot(sol;  vars = [STG.AB₊V]).
"""

using NeuronBuilder, ModelingToolkit, OrdinaryDiffEq, Plots

const Area = 0.0628 # Prinz/Liu 0.0628 mm2
const Cm = 10.0 # specific capacitance cₘ is a biological constant (around) 10 nF/mm^2

Prinz_conv = Cm / Area
synaptic_conv = 1e-3 / Area^2 #this gives μS/mm^2

τCa = 200.0
Ca∞ = 0.05

defaults = Dict(Voltage => BasicVoltageDynamics(),
    Calcium => Prinz.CalciumDynamics(τCa, Ca∞, 14.96 * 0.0628),
    Reversal{Calcium} => Prinz.CaReversalDynamics())

somatic_parameters = Dict(
    Reversal{Sodium} => 50.0,
    Reversal{Potassium} => -80.0,
    Reversal{Leak} => -50.0,
    Reversal{Proton} => -20.0,
    Voltage => -60.0,
    Calcium => 0.05,
    Reversal{Calcium} => 0.0)

#Using parameters from Prinz (2004) Similar network activity from disparate circuit parameters
#these are specifically in Table 2: AB2, LP4, PY1 for figure 3e
function AB_Neuron(num_connections)
    AB2_ch = [
        Prinz.Na(100.0 * Prinz_conv),
        Prinz.CaS(6.0 * Prinz_conv),
        Prinz.CaT(2.5 * Prinz_conv),
        Prinz.H(0.01 * Prinz_conv),
        Prinz.Ka(50.0 * Prinz_conv),
        Prinz.KCa(5.0 * Prinz_conv),
        Prinz.Kdr(100.0 * Prinz_conv),
        Prinz.leak(0.0 * Prinz_conv)
    ]
    b = BasicNeuron(NoGeometry(Cm), defaults, somatic_parameters, AB2_ch, :AB)
    return b(incoming_connections=num_connections)
end

function PY_Neuron(num_connections)
    PY1_ch = [Prinz.Na(100.0 * Prinz_conv),
        Prinz.CaS(2.0 * Prinz_conv),
        Prinz.CaT(2.4 * Prinz_conv),
        Prinz.H(0.05 * Prinz_conv),
        Prinz.Ka(50.0 * Prinz_conv),
        Prinz.KCa(0.0 * Prinz_conv),
        Prinz.Kdr(125.0 * Prinz_conv),
        Prinz.leak(0.01 * Prinz_conv)
    ]
    somatic_parameters[Voltage] = -55.0
    b = BasicNeuron(NoGeometry(Cm), defaults, somatic_parameters, PY1_ch, :PY)
    return b(incoming_connections=num_connections)
end

function LP_Neuron(num_connections)
    LP4_ch = [Prinz.Na(100.0 * Prinz_conv),
        Prinz.CaS(4.0 * Prinz_conv),
        Prinz.CaT(0.0 * Prinz_conv),
        Prinz.H(0.05 * Prinz_conv),
        Prinz.Ka(20.0 * Prinz_conv),
        Prinz.KCa(0.0 * Prinz_conv),
        Prinz.Kdr(25.0 * Prinz_conv),
        Prinz.leak(0.03 * Prinz_conv)
    ]
    somatic_parameters[Voltage] = -65.0
    b = BasicNeuron(NoGeometry(Cm), defaults, somatic_parameters, LP4_ch, :LP)
    return b(incoming_connections=num_connections)
end

function nodes(i)
    (i == 1) && return n -> AB_Neuron(n)
    (i == 2) && return n -> PY_Neuron(n)
    (i == 3) && return n -> LP_Neuron(n)
end

intoAB = [(EmptyConnection(),)
    (EmptyConnection(),) #hi
    (Glut(30.0 * synaptic_conv),)
]

intoPY = [
    (Chol(3.0 * synaptic_conv), Glut(10.0 * synaptic_conv))
    (EmptyConnection(),)
    (Glut(1.0 * synaptic_conv),)
]

intoLP = [
    (Chol(30.0 * synaptic_conv), Glut(30.0 * synaptic_conv))
    (Glut(30.0 * synaptic_conv),)
    (EmptyConnection(),)
]

Adjacency_matrix = cat(intoAB, intoPY, intoLP; dims=2)

function edges(i, j)
    return Adjacency_matrix[i, j]
end

tspan = (0.0, 10000.0)
stg = @time build_network(nodes, edges, 1:3, name=:stg)
prob = ODEProblem(stg, [], tspan, []; jac=true)

@time sol = solve(prob, AutoTsit5(Rosenbrock23()))

plot(sol, xlims=(5000, 10000), fg_legend=:transparent; vars=[stg.AB₊V, stg.LP₊V, stg.PY₊V, stg.PY₊CaS₊ICa], layout=(4, 1))

"""
Another way of adding connections without having to specify a connectivity matrix (no empty connections)
"""
AB = AB_Neuron(1)
LP = LP_Neuron(3)
PY = PY_Neuron(3)
neurons = [AB, PY, LP]

ABLP_chol = directed_synapse(AB, LP, Chol(30.0 * synaptic_conv))
ABPY_chol = directed_synapse(AB, PY, Chol(3.0 * synaptic_conv))
ABLP_glut = directed_synapse(AB, LP, Glut(30.0 * synaptic_conv))
ABPY_glut = directed_synapse(AB, PY, Glut(10.0 * synaptic_conv))

LPAB_glut = directed_synapse(LP, AB, Glut(30.0 * synaptic_conv))
LPPY_glut = directed_synapse(LP, PY, Glut(1.0 * synaptic_conv))

PYLP_glut = directed_synapse(PY, LP, Glut(30.0 * synaptic_conv))

connections = [LPAB_glut, ABPY_chol, ABPY_glut, LPPY_glut, ABLP_chol, ABLP_glut, PYLP_glut]

stg_2 = @time add_all_connections(neurons, connections)
prob_2 = ODEProblem(stg_2, [], tspan, []; jac=true)

@time sol_2 = solve(prob_2, AutoTsit5(Rosenbrock23()))

plot(sol_2, xlims=(5000, 10000), fg_legend = :transparent; vars=[stg_2.AB₊V, stg_2.LP₊V, stg_2.PY₊V, stg_2.PY₊CaS₊ICa], layout=(4, 1))


text:
using Test
using InvariantSets
using LinearAlgebra: I
using Polyhedra # needed for prerset_vrep, preset_hrep
using CDDLib # needed for preset_hrep

@testset "Preset computation" begin
    n = 2
    m = 1
    A = rand(n,n)
    B = rand(n,m)
    X = BallInf(zeros(n), 10.)
    U = BallInf(zeros(m), .1)
    K = [0.0890246  0.484408]
    F, g = tosimplehrep(X)

    P = preset(A, X)
    @test P ⊆ HPolytope(F*A,g) && HPolytope(F*A,g) ⊆ P
    # @test isequivalent(preset(A, X), HPolytope(F*A,g))

    D = state_constraints(X, U, K)
    D2 = InvariantSets.state_constraints_ops(X, U, K)
    @test D ⊆ D2 && D2 ⊆ D

    P_2 = InvariantSets._preset_vrep(A, B, U, D)
    P_3 = InvariantSets._preset_hrep(A, B, U, D)
    @test P_2 ⊆ P_3 && P_3 ⊆ P_2
end

@testset "Maximum invariant set" begin
    A2 = [1.0 1.0; 0.0 1.0]
    B2 = [1.0; 0.5]
    Q2 = I
    R2 = 90*I
    # K2 = dlqr(A2,B2,Q2,R2)
    K2 = [0.0890246  0.484408]
    X2 = BallInf(zeros(2), 2.)
    U2 = BallInf(zeros(1), 2.)
    S2 = A2-B2*K2

    max_PI_1 =  maximum_invariant_set(S2, X2)
    max_PI_2 =  maximum_invariant_set(S2, X2, 5)
    max_PI_3 =  maximum_invariant_set(S2, X2, 4)

    @test max_PI_1 ⊆ max_PI_2 && max_PI_2 ⊆ max_PI_1
    @test !(max_PI_1 ⊆ max_PI_3 && max_PI_3 ⊆ max_PI_1)

    max_PI_4 = maximum_invariant_set(A2, B2, X2, U2, K2)
    @test max_PI_1 ⊆ max_PI_4 && max_PI_4 ⊆ max_PI_1
end

@testset "Maximum control invariant set" begin
    A2 = [1.0 1.0; 0.0 1.0]
    B2 = [1.0; 0.5]
    X2 = BallInf(zeros(2), 5.)
    U2 = BallInf(zeros(1), 1.)
    max_CPI_1 =  maximum_control_invariant_set(A2, B2, X2, U2; preset=InvariantSets._preset_vrep)
    max_CPI_2 =  maximum_control_invariant_set(A2, B2, X2, U2; preset=InvariantSets._preset_hrep)
    max_CPI_3 =  maximum_control_invariant_set(A2, B2, X2, U2, 5)
    @test max_CPI_1 ⊆ max_CPI_2 && max_CPI_2 ⊆ max_CPI_1
    @test !(max_CPI_1 ⊆ max_CPI_3 && max_CPI_3 ⊆ max_CPI_1)
end

max_stars_repo_name: mforets/InvariantSets

max_issues_repo_name: mforets/InvariantSets


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

using Markdown
using InteractiveUtils

# ╔═╡ e077ab16-22e5-11eb-15ed-5f0654bf7328
using Pkg, DrWatson

# ╔═╡ e75a9a06-22e5-11eb-2e36-8d4f62b830ed
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
	using Statistics, DataFrames, CSV
end;

# ╔═╡ ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
md"## Listing6.08"

# ╔═╡ 1aa41898-3afa-11eb-2f49-570853ce6435
begin
	df = DataFrame()
	for eps in [0.1, 0.05, 0.02, 0.01]
    	n = ceil(1/eps^2)
		append!(df, DataFrame(:epsilon => eps, :n => n))
	end
	df
end

# ╔═╡ 475ff888-22e6-11eb-2354-09f8f40a8e12
md"## End of listing6.08"

# ╔═╡ Cell order:
# ╟─ca5fbdc8-22e5-11eb-0a60-cb0a127a3ca5
# ╠═e077ab16-22e5-11eb-15ed-5f0654bf7328
# ╠═e75a9a06-22e5-11eb-2e36-8d4f62b830ed
# ╠═1aa41898-3afa-11eb-2f49-570853ce6435
# ╟─475ff888-22e6-11eb-2354-09f8f40a8e12


text:
using HostCPUFeatures
using Test

@testset "HostCPUFeatures.jl" begin

  println("vector_width.jl")
  @time @testset "vector_width.jl" begin
    for T ∈ (Float32,Float64)
      @test @inferred(HostCPUFeatures.pick_vector_width(T)) * @inferred(HostCPUFeatures.static_sizeof(T)) === @inferred(HostCPUFeatures.register_size(T)) === @inferred(HostCPUFeatures.register_size())
    end
    for T ∈ (Int8,Int16,Int32,Int64,UInt8,UInt16,UInt32,UInt64)
      @test @inferred(HostCPUFeatures.pick_vector_width(T)) * @inferred(HostCPUFeatures.static_sizeof(T)) === @inferred(HostCPUFeatures.register_size(T)) === @inferred(HostCPUFeatures.simd_integer_register_size())
    end
    @test HostCPUFeatures.static_sizeof(BigFloat) === HostCPUFeatures.static_sizeof(Int)
    @test HostCPUFeatures.static_sizeof(Float32) === HostCPUFeatures.static_sizeof(Int32) === HostCPUFeatures.StaticInt(4)

    @test @inferred(HostCPUFeatures.pick_vector_width(Float64, Int32, Float64, Float32, Float64)) * HostCPUFeatures.static_sizeof(Float64) === @inferred(HostCPUFeatures.register_size())
    @test @inferred(HostCPUFeatures.pick_vector_width(Float64, Int32)) * HostCPUFeatures.static_sizeof(Float64) === @inferred(HostCPUFeatures.register_size())

    @test @inferred(HostCPUFeatures.pick_vector_width(Float32, Float32)) * HostCPUFeatures.static_sizeof(Float32) === @inferred(HostCPUFeatures.register_size())
    @test @inferred(HostCPUFeatures.pick_vector_width(Float32, Int32)) * HostCPUFeatures.static_sizeof(Float32) === @inferred(HostCPUFeatures.simd_integer_register_size())

    FTypes = (Float32, Float64)
    Wv = ntuple(i -> @inferred(HostCPUFeatures.register_size()) >> (i+1), Val(2))
    for (T, N) in zip(FTypes, Wv)
      W = @inferred(HostCPUFeatures.pick_vector_width(T))
      # @test Vec{Int(W),T} == HostCPUFeatures.pick_vector(W, T) == HostCPUFeatures.pick_vector(T)
      @test W == @inferred(HostCPUFeatures.pick_vector_width(W, T))
      @test W === @inferred(HostCPUFeatures.pick_vector_width(W, T)) == @inferred(HostCPUFeatures.pick_vector_width(T))
      while true
        W >>= HostCPUFeatures.One()
        W == 0 && break
        W2, Wshift2 = @inferred(HostCPUFeatures.pick_vector_width_shift(W, T))
        @test W2 == HostCPUFeatures.One() << Wshift2 == @inferred(HostCPUFeatures.pick_vector_width(W, T)) == HostCPUFeatures.pick_vector_width(Val(Int(W)),T)  == W
        @test HostCPUFeatures.StaticInt(W) === HostCPUFeatures.pick_vector_width(Val(Int(W)), T) === HostCPUFeatures.pick_vector_width(W, T)
        for n in W+1:2W
          W3, Wshift3 = HostCPUFeatures.pick_vector_width_shift(HostCPUFeatures.StaticInt(n), T)
          @test W2 << 1 == W3 == 1 << (Wshift2+1) == 1 << Wshift3 == HostCPUFeatures.pick_vector_width(HostCPUFeatures.StaticInt(n), T) == HostCPUFeatures.pick_vector_width(Val(n),T) == W << 1
          # @test HostCPUFeatures.pick_vector(W, T) == HostCPUFeatures.pick_vector(W, T) == Vec{Int(W),T}
        end
      end
    end
    @test HostCPUFeatures.pick_vector_width(Float16) === HostCPUFeatures.pick_vector_width(Float32)
    # @test HostCPUFeatures.nextpow2(0) == 1

    @test HostCPUFeatures.unwrap(HostCPUFeatures.static(2)) === 2
    @test HostCPUFeatures.unwrap(HostCPUFeatures.static(1.2)) === 1.2
    @test HostCPUFeatures.unwrap(HostCPUFeatures.static(:a)) === :a
  end

end


text:
function finish_grid_bot(client, config, bot_path)
    trade_grid = JSON.parsefile(joinpath(bot_path, "trade_grid.json"));

    #
    ## Cancel all orders
    all_canceled = true
    for (price, v) in trade_grid
        if haskey(v, "order")
            order = v["order"]
            if !(order["status"] in [client.ORDER_STATUS_FILLED, client.ORDER_STATUS_CANCELED, client.ORDER_STATUS_EXPIRED])
                sleep(0.5)
                result = try client.get_order(symbol=order["symbol"], orderId=order["orderId"])
                    catch
                    nothing
                end
                if result !== nothing
                    sleep(0.5)
                    if !(result["status"] in [client.ORDER_STATUS_FILLED, client.ORDER_STATUS_CANCELED, client.ORDER_STATUS_EXPIRED])
                        result = try client.cancel_order(symbol=order["symbol"], orderId=order["orderId"])
                            catch
                            nothing
                        end
                    end
                else
                    all_canceled = false
                    continue
                end
                if result !== nothing
                    trade_grid[price]["order"] = result
                else
                    all_canceled = false
                    continue
                end
            end
        end
    end
    open(joinpath(bot_path, "trade_grid.json"), "w") do f
        JSON.print(f, trade_grid, 4)
    end
    if all_canceled == false
        return
    end
    trade_grid = JSON.parsefile(joinpath(bot_path, "trade_grid.json"));
    leftover_bought = 0
    for (k,v) in trade_grid
        if haskey(v, "order")
            leftover_bought += typeof(v["order"]["executedQty"]) == Float64 ? v["order"]["executedQty"] : parse(Float64, v["order"]["executedQty"])
        end
    end
    if leftover_bought == 0
        config["status"] = "finished"
    else
        config["status"] = "leftover"
    end
    open(joinpath(bot_path, "config.json"), "w") do f
        JSON.print(f, config, 4)
    end
end
export finish_grid_bot


text:
# A maximally fast upper branch incomplete gamma function when the first
# argument is a non-positive integer.
#
# (_gamma_upper_negative_integer)
@inline function _g_u_n_i(s::Int64, x, g0, expnx)
  n    = -s
  ser  = zero(x)
  _x   = one(x)
  _sgn = copy(_x)
  facn = convert(typeof(x), factorial(n))
  fac  = facn/n
  for k in 0:(n-1)
    term  = _sgn*fac*_x
    ser  += term
    _x   *= x
    _sgn *= -one(x)
    fac  /= (n-k-1)
  end
  ((expnx/_x)*ser + _sgn*g0)/facn
end

# The best speed I was able to accomplish with this thing was just to
# compartmentalize it into its own function. Definitely not ideal, but so it is.
#
# TODO (cg 2021/10/26 15:58): get rid of all the remaining factorial calls so
# that this won't literally break for order greater than, like, 19.
function exponential_improvement(v, x, l, m)
  onex = one(x)
  twox = onex+onex
  fv     = 4*v*v
  ser    = zero(x)
  _z     = x
  floatj = onex
  ak_num = fv - floatj
  factj  = onex
  twofloatj = onex
  eightj  = 8
  expx    = exp(2*x)
  expnx   = exp(-2*x)
  expintx = -expinti(-2*x)
  ser     = expx*factorial(l-1)*_g_u_n_i(1-l, 2*x, expintx, expnx)/(2*pi)
  for j in 1:(m-1)
    # add to the series:
    s   = Int(l-j)
    _g  = expx*factorial(Int(s-1))*_g_u_n_i(1-s, 2*x, expintx, expnx)/(2*pi)
    term       = ak_num/(factj*_z*eightj)*_g
    ser       += term
    # update ak and _z:
    floatj    += onex
    twofloatj += twox
    factj     *= floatj
    ak_num    *= (fv - twofloatj^2)
    _z        *= x
    eightj    *= 8
  end
  _sgn = isodd(l) ? -one(x) : one(x)
  ser*_sgn*twox*cospi(v)
end

function _besselk_as(v, x, order, use_remainder=true, modify=false)
  onex = one(x)
  twox = onex+onex
  fv     = 4*v*v
  ser    = zero(x)
  _z     = x
  ser    = onex #zero(x) #onex/_z
  floatj = onex
  ak_num = fv - floatj
  factj  = onex
  twofloatj = onex
  eightj = 8
  for j in 1:order
    # add to the series:
    term       = ak_num/(factj*_z*eightj)
    ser       += term
    # update ak and _z:
    floatj    += onex
    twofloatj += twox
    factj     *= floatj
    ak_num    *= (fv - twofloatj^2)
    _z        *= x
    eightj    *= 8
  end
  if use_remainder
    _rem = exponential_improvement(v, x, Int(order+1), Int(order))
    ser += _rem
  end
  # if you're modifying as (x^v)*besselk(v,x), since the series part is pretty
  # stable numerically, what we want to deal with is the (x^v)*exp(-x). That's
  # the problem of potentially huge*tiny.
  if modify
    mulval = exp(v*log(x)-x)
  else
    mulval = exp(-x)
  end
  return sqrt(pi/(x*twox))*mulval*ser
end


text:
using Random
using ComoniconTypes
using ComoniconTestUtils
using Test

Random.seed!(42)

@testset "ComoniconTestUtils.jl" begin
    ComoniconTestUtils.test_function("a", "b", "c"; option_a=1, option_b=2)
    @test_args ["a", "b", "c"]
    @test_kwargs [:option_a=>1, :option_b=>2]
end

find_leaf(cmd::Entry, inputs, current=1) = find_leaf(cmd.root, inputs, current)
find_leaf(cmd::LeafCommand, inputs, current=1) = cmd

function find_leaf(cmd::NodeCommand, inputs, current=1)
    @test haskey(cmd.subcmds, inputs[current])
    return find_leaf(cmd.subcmds[inputs[current]], inputs, current+1)
end

@testset "random_command + random_input" begin
    for _ in 1:10
        cmd = rand_command()
        inputs = rand_input(cmd)
        leaf = find_leaf(cmd, inputs)

        for each in inputs
            if startswith(each, "--")
                name = first(split(each, '='))
                name = lstrip(name, '-')
                @test haskey(leaf.options, name) || haskey(leaf.flags, name)
            end
        end
    end
end


text:
#=
const unit_excludes = String[
    "number_threads",                 # EAGO won't support number of threads in near future
    "raw_status_string",              # TODO: ADD internal status states to EAGO
    "solve_unbounded_model",          # CBC returns infeasible or unbounded for linear...
    "solve_zero_one_with_bounds_3",   # GLPK has a non-standard return code
    "solve_result_index",             # TODO: Should throw error when querying for multiple results... (expected behavior?)
    "solve_qcp_edge_cases",           # Not box constrained NLP type problems...
    "solve_qp_zero_offdiag",
    "solve_qp_edge_cases",
    "solve_affine_deletion_edge_cases"   # TODO: Fix this
]

const contlinear_excludes = String[
    "partial_start",  # EAGO doesn't support VariablePrimalStart
    "linear1",   # TODO: Fix this

    #=
    "linear13",
    "linear8a",
    "linear14",
    "linear6",
    "linear4",
    "linear3",
    "linear9",
    "linear8c",
    "linear2",
    "linear12",
    "linear7",
    "linear8b",
    "linear10b",
    "linear10",
    "linear15",
    "linear5",
    "linear11"
    =#
]

const intlinear_excludes = String[
    "indicator1",  # doesn't currently support indicator sets
    "indicator2",  # can't check using Cbc until https://github.com/jump-dev/Cbc.jl/issues/151 is resolved
    "indicator3",
    "indicator4",

    "int2",        # currently doesn't support sos1 or sos2 constraints
    "int3",   # TODO: Fix this

    # Passing
    #"semiinttest",
    #"semiconttest",
    #"int1",
    #"knapsack"
]

const contconic_excludes = String[
    "dualexp",  # Not directly bridged to common cones
    "dualpow",
    "logdet",
    "rootdet",
    "sdp",
    "normnuc",
    "exp",
    "soc",
    "normspec",
    "relentr",
    "rsoc",
    "pow",
    "geomean"
]

const contquadratic_excludes = String[
    "ncqcp",
    "qp",
    "socp",
    "qcp",
]

function test_moi(T::Type{<:Real}; solver_options...)

    optimizer = MOIU.CachingOptimizer(MOIU.UniversalFallback(MOIU.Model{T}()), EAGO.Optimizer())
    MOI.set(optimizer, MOI.RawParameter("verbosity"), 0)

    tol = 2sqrt(sqrt(eps(T)))
    config = MOIT.Config(T;
    atol = tol,
    rtol = tol,
    #solve = true,
    #query = true,
    #modify_lhs = false,
    #duals = false,
    #infeas_certificates = false,
    )

    @testset "unit tests" begin
        MOIT.unittest(MOIB.full_bridge_optimizer(optimizer, T), config, unit_excludes)
    end

    @testset "continuous linear tests" begin
        MOIT.contlineartest(MOIB.full_bridge_optimizer(optimizer, T), config, contlinear_excludes)
    end

    @testset "mixed-integer linear tests" begin
        MOIT.intlineartest(MOIB.full_bridge_optimizer(optimizer, T), config, intlinear_excludes)
    end

    @testset "continuous conic tests" begin
        MOIT.contconictest(MOIB.full_bridge_optimizer(optimizer, T), config, contconic_excludes)
    end

    @testset "continuous quadratic tests" begin
        MOIT.contquadratictest(MOIB.full_bridge_optimizer(optimizer, T), config, contquadratic_excludes)
    end
end

# Test with mip_solver = Cbc as it supports SOS1 & SOS2 constraints
# TODO: Use bridges for SOS1 & SOS2 constraint if unsupported
# Need to test with GLPK as well to ensure subsolver supports constraint
# coefficient modification.
test_moi(Float64)
=#

module TestEAGO

import EAGO
using MathOptInterface
using Test

const MOI = MathOptInterface
const OPTIMIZER = MOI.instantiate(MOI.OptimizerWithAttributes(EAGO.Optimizer, MOI.Silent() => true))
const BRIDGED = MOI.instantiate(MOI.OptimizerWithAttributes(EAGO.Optimizer, MOI.Silent() => true), with_bridge_type = Float64)
const CONFIG = MOI.Test.Config(atol = 1e-3, rtol = 1e-3, optimal_status = MOI.OPTIMAL,
                               exclude = Any[MOI.DualObjectiveValue, MOI.ConstraintBasisStatus, MOI.VariableName, MOI.ConstraintName, MOI.delete,
                                             MOI.ConstraintDual, MOI.ListOfModelAttributesSet, MOI.add_constrained_variables])

"""
    runtests()

This function runs all functions in the this Module starting with `test_`.
"""
function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
end

"""
    test_runtests()

This function runs all the tests in MathOptInterface.Test.

Pass arguments to `exclude` to skip tests for functionality that is not
implemented or that your solver doesn't support.
"""
function test_runtests()
    MOI.Test.runtests(BRIDGED, CONFIG,
                      exclude = [# IPOPT Inherited test exclusions
                                "test_model_ScalarFunctionConstantNotZero",
                                "test_solve_TerminationStatus_DUAL_INFEASIBLE",
                                "test_linear_VectorAffineFunction_empty_row",
                                "test_solve_DualStatus_INFEASIBILITY_CERTIFICATE_",
                                "test_model_LowerBoundAlreadySet",
                                "test_model_UpperBoundAlreadySet",
                                "test_model_copy_to_UnsupportedAttribute",
                                "test_model_copy_to_UnsupportedConstraint",
                                "test_objective_set_via_modify",
                                "test_model_ModelFilter_ListOfConstraintIndices",
                                "test_model_ModelFilter_ListOfConstraintTypesPresent",
                                # Cbc default test exlucisons
                                "_Indicator_",
                                "test_linear_SOS1_integration",
                                "test_linear_SOS2_integration",
                                "test_solve_SOS2_add_and_delete",
                                "test_conic_NormInfinityCone_INFEASIBLE",
                                "test_conic_NormOneCone_INFEASIBLE",
                                "test_solve_TerminationStatus_DUAL_INFEASIBLE",

                                # EAGO test exclusions
                                "test_attribute_NumberOfThreads",
                                "test_modification_",
                                "test_linear_integration_delete_variables",

                                "test_quadratic_SecondOrderCone_basic",

                                "test_quadratic_constraint_GreaterThan",
                                "test_quadratic_constraint_LessThan",
                                "test_quadratic_constraint_basic",
                                "test_quadratic_constraint_minimize",
                                "test_quadratic_duplicate_terms",
                                "test_quadratic_integration",
                                "test_quadratic_nonconvex_constraint_integration",
                                "test_quadratic_nonhomogeneous",

                                "test_constraint_qcp_duplicate_diagonal",
                                "test_constraint_qcp_duplicate_off_diagonal",
                                "test_objective_get_ObjectiveFunction_ScalarAffineFunction",
                                "test_objective_qp_ObjectiveFunction_edge_cases",
                                "test_objective_qp_ObjectiveFunction_zero_ofdiag",
                                ],
                      exclude_tests_after = v"0.10.5")
end

end

@testset "MOI" begin
    TestEAGO.runtests()
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
"""
`KeyPair{T}(𝑑::BigInt, 𝑄::Point)` represents a private-public key pair in which
- `𝑑` is the secret key
- `𝑄` the public key
= `T` represents key pair scheme

`KeyPair{:ECDSA}(𝑑::Integer)` instantiate a `KeyPair` such as 𝑄 = 𝑑𝐺 and where
- {𝑑 ∈ ℤ | 𝑑 < 𝑛}
- 𝐺 is the secp256k1 base point.
- 𝑛 is the order of 𝐺
"""
struct KeyPair{T}
    𝑑::BigInt
    𝑄::Point
end

"""
Signature{T}(𝑟::BigInt, 𝑠::BigInt) represents a following scheme `T`
- `T` can be set to `:ECDSA`
"""
struct Signature{T}
    𝑟::BigInt
    𝑠::BigInt
end

"Formats Signature as (r, s) in hexadecimal format"
function show(io::IO, z::Signature)
    print(io, "scep256k1 signature(𝑟, 𝑠):\n", string(z.𝑟, base = 16), ",\n", string(z.𝑠, base = 16))
end

==(x::Signature, y::Signature) = x.𝑟 == y.𝑟 && x.𝑠 == y.𝑠

"""
    serialize(x::Signature) -> Vector{UInt8}

Serialize a `Signature` to DER format
"""
function serialize(x::Signature)
    rbin = bytes(x.𝑟)
    # if rbin has a high bit, add a 00
    if rbin[1] >= 128
        rbin = pushfirst!(rbin, 0x00)
    end
    prepend!(rbin, bytes(length(rbin)))
    pushfirst!(rbin, 0x02)

    sbin = bytes(x.𝑠)
    # if sbin has a high bit, add a 00
    if sbin[1] >= 128
        sbin = pushfirst!(sbin, 0x00)
    end
    prepend!(sbin, bytes(length(sbin)))
    pushfirst!(sbin, 0x02)

    result = sbin
    prepend!(result, rbin)
    prepend!(result, bytes(length(result)))
    return pushfirst!(result, 0x30)
end

"""
    Signature(x::Vector{UInt8}; scheme::Symbol) -> Signature

Parse a DER binary to a `Signature{scheme}`
- `scheme` is optional and set to `:ECDSA` by default.
"""
function Signature(x::Vector{UInt8}; scheme::Symbol=:ECDSA)
    io = IOBuffer(x)
    prefix = read(io, 1)[1]
    if prefix != 0x30
        throw(PrefixError())
    end
    len = read(io, 1)[1]
    if len + 2 != length(x)
        throw(LengthError())
    end
    prefix = read(io, 1)[1]
    if prefix != 0x02
        throw(PrefixError())
    end
    rlength = read(io, 1)[1]
    r = to_int(read(io, rlength))
    prefix = read(io, 1)[1]
    if prefix != 0x02
        throw(PrefixError())
    end
    slength = read(io, 1)[1]
    s = to_int(read(io, slength))
    if length(x) != 6 + rlength + slength
        throw(LengthError())
    end
    return Signature{scheme}(r, s)
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -0.27492067550218746 <= q <= 2.084201684402027)

                     add_NL_constraint(m, :(log(1 + exp(-0.0416569430659357 + 0.13675627334278007*log(1 + exp(0.1752604319892206 + -0.8557289632032381*log(1 + exp(-0.028885055816546323 + 0.15091441453709065*log(1 + exp(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3]))) + -0.8760277243880852*log(1 + exp(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))) + 0.21149394336428795*log(1 + exp(-0.29077418744555983 + -0.11119024401630329*log(1 + exp(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3]))) + 0.4817613487326513*log(1 + exp(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))))) + 0.703085236566185*log(1 + exp(-0.4914467550189556 + 0.5078124686195609*log(1 + exp(-0.028885055816546323 + 0.15091441453709065*log(1 + exp(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3]))) + -0.8760277243880852*log(1 + exp(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))) + 0.7523099505151869*log(1 + exp(-0.29077418744555983 + -0.11119024401630329*log(1 + exp(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3]))) + 0.4817613487326513*log(1 + exp(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))))))) + log(1 + exp(0.6991595157542867 + 0.5132010992296534*log(1 + exp(0.1752604319892206 + -0.8557289632032381*log(1 + exp(-0.028885055816546323 + 0.15091441453709065*log(1 + exp(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3]))) + -0.8760277243880852*log(1 + exp(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))) + 0.21149394336428795*log(1 + exp(-0.29077418744555983 + -0.11119024401630329*log(1 + exp(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3]))) + 0.4817613487326513*log(1 + exp(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))))) + 0.6682334830395109*log(1 + exp(-0.4914467550189556 + 0.5078124686195609*log(1 + exp(-0.028885055816546323 + 0.15091441453709065*log(1 + exp(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3]))) + -0.8760277243880852*log(1 + exp(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))) + 0.7523099505151869*log(1 + exp(-0.29077418744555983 + -0.11119024401630329*log(1 + exp(-0.5352161553133312 + 0.05336126885466008*$(x[1]) + 0.7318437425770079*$(x[2]) + -0.7787417630340121*$(x[3]))) + 0.4817613487326513*log(1 + exp(0.0006254708492061667 + 0.44360632638752495*$(x[1]) + -0.9819405794216607*$(x[2]) + 0.2295552892872319*$(x[3]))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m


text:
# Autogenerated wrapper script for czmq_jll for x86_64-apple-darwin
export libczmq

using ZeroMQ_jll
JLLWrappers.@generate_wrapper_header("czmq")
JLLWrappers.@declare_library_product(libczmq, "@rpath/libczmq.4.dylib")
function __init__()
    JLLWrappers.@generate_init_header(ZeroMQ_jll)
    JLLWrappers.@init_library_product(
        libczmq,
        "lib/libczmq.4.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
# Autogenerated wrapper script for FFplay_jll for i686-linux-gnu
export ffplay

using FFMPEG_jll
using SDL2_jll
## Global variables
PATH = ""
LIBPATH = ""
LIBPATH_env = "LD_LIBRARY_PATH"
LIBPATH_default = ""

# Relative path to `ffplay`
const ffplay_splitpath = ["bin", "ffplay"]

# This will be filled out by __init__() for all products, as it must be done at runtime
ffplay_path = ""

# ffplay-specific global declaration
function ffplay(f::Function; adjust_PATH::Bool = true, adjust_LIBPATH::Bool = true)
    global PATH, LIBPATH
    env_mapping = Dict{String,String}()
    if adjust_PATH
        if !isempty(get(ENV, "PATH", ""))
            env_mapping["PATH"] = string(PATH, ':', ENV["PATH"])
        else
            env_mapping["PATH"] = PATH
        end
    end
    if adjust_LIBPATH
        LIBPATH_base = get(ENV, LIBPATH_env, expanduser(LIBPATH_default))
        if !isempty(LIBPATH_base)
            env_mapping[LIBPATH_env] = string(LIBPATH, ':', LIBPATH_base)
        else
            env_mapping[LIBPATH_env] = LIBPATH
        end
    end
    withenv(env_mapping...) do
        f(ffplay_path)
    end
end

"""
Open all libraries
"""
function __init__()
    global artifact_dir = abspath(artifact"FFplay")

    # Initialize PATH and LIBPATH environment variable listings
    global PATH_list, LIBPATH_list
    # From the list of our dependencies, generate a tuple of all the PATH and LIBPATH lists,
    # then append them to our own.
    foreach(p -> append!(PATH_list, p), (FFMPEG_jll.PATH_list, SDL2_jll.PATH_list,))
    foreach(p -> append!(LIBPATH_list, p), (FFMPEG_jll.LIBPATH_list, SDL2_jll.LIBPATH_list,))

    global ffplay_path = normpath(joinpath(artifact_dir, ffplay_splitpath...))

    push!(PATH_list, dirname(ffplay_path))
    # Filter out duplicate and empty entries in our PATH and LIBPATH entries
    filter!(!isempty, unique!(PATH_list))
    filter!(!isempty, unique!(LIBPATH_list))
    global PATH = join(PATH_list, ':')
    global LIBPATH = join(vcat(LIBPATH_list, [joinpath(Sys.BINDIR, Base.LIBDIR, "julia"), joinpath(Sys.BINDIR, Base.LIBDIR)]), ':')
end  # __init__()


text:
# Load / save test

using ValueHistoriesLH, Test

function load_uni_test()
    @testset "Univalue history" begin
        h = make_test_history(3, [2,3]; keySpacing = 3);
        fPath = joinpath(test_dir(),  "save_history_test");

        sName = :test1;
        sp = ValueHistoriesLH.series_path(fPath, sName);
        @test endswith(sp, "_test1.json3")

        if isfile(sp)
            rm(sp);
        end
        save_history(h, sp; showMsg = true);
        @test isfile(sp)

        # Remember that loaded history makes matrices into vectors
        h2 = load_history(sp);
        @test length(h) == length(h2)
        for j = 1 : length(h)
            v1 = values(h)[j];
            v2 = values(h2)[j];
            @test isapprox(v1,  reshape(v2, size(v1)))
        end
    end
end

function load_mv_test()
    @testset "MV History" begin
        h = make_test_mvhistory(3);
        fPath = joinpath(test_dir(),  "save_mvhistory_test");
        if isfile(fPath)
            rm(fPath);
        end
        save_mvhistory(h, fPath);
        @test isfile(ValueHistoriesLH.meta_path(fPath))
        v = ValueHistoriesLH.load_meta(fPath);
        @test isa(v, Vector{Symbol})

        h2 = load_mvhistory(fPath);
        @test isa(h2, MVHistory);

        for key in keys(h)
            v1 = values(h, key);
            v2 = values(h2, key);
            @test isapprox(v1, reshape(v2, size(v1)))
        end
    end
end

@testset "Load / save" begin
    load_uni_test();
    load_mv_test();
end

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


text:
#=
Example of importing a module. This example uses the module contained in the
folder named MyModule.
=#

# Load the module file
include("MyModule/MyModule.jl")

# Import the module
import .MyModule

# Use the module functions
MyModule.myfuncA()
MyModule.myfuncB()

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
# # Honeycomb Lattice in k-space

using LinearAlgebra
using Plots
using LatticeTools

simplifyname(s::AbstractString) = (
    s
    |> (x-> replace(x, r"(.+?)" => s"\1"))
    |> (x-> replace(x, r"(.+?)" => s"[\1]"))
)

function makewithin(extent::AbstractVector{<:Real})
    a, b, c, d = extent
    (x::Real, y::Real) -> ((a <= x <= b) && (c <= y <= d))
end
mkpath("example_little_symmetry_kspace_honeycomb")

# ## Honeycome Lattice

# ### Define UnitCell
latticevectors = [1 -0.5; 0 sqrt(3)*0.5];
unitcell = makeunitcell(latticevectors; SiteType=String);
addsite!(unitcell, "A", carte2fract(unitcell, [0.5, 0.5 / sqrt(3)]));
addsite!(unitcell, "B", carte2fract(unitcell, [0.5, -0.5 / sqrt(3)]));

# ### Plot Lattice
let bravais_lattice = [], a_sites = [], b_sites = []
    extent = [-2, 2, -2, 2]
    within = makewithin(extent)
    for i1 in -5:5, i2 in -5:5
        R = latticevectors * [i1, i2]
        within(R...) && push!(bravais_lattice, R)
        for (orb_name, orb_fc) in unitcell.sites
            orb_cc = fract2carte(unitcell, orb_fc)
            r = R + orb_cc
            within(r...) && push!(orb_name == "A" ? a_sites : b_sites, r)
        end
    end
    bravais_lattice = hcat(bravais_lattice...)
    a_sites = hcat(a_sites...)
    b_sites = hcat(b_sites...)

    img = plot(aspect_ratio=1)
    scatter!(bravais_lattice[1,:], bravais_lattice[2,:], color="black", markershape=:circle, markersize=3, label="Bravais")
    scatter!(a_sites[1,:], a_sites[2,:], color="red", markerstrokewidth=0, markersize=5, label="A")
    scatter!(b_sites[1,:], b_sites[2,:], color="blue", markerstrokewidth=0, markersize=5, label="B")
    xlims!(extent[1], extent[2])
    ylims!(extent[3], extent[4])
    savefig(img, "example_little_symmetry_kspace_honeycomb/realspace.svg")
end

# ![](example_little_symmetry_kspace_honeycomb/realspace.svg)

# ## √3 × √3 Supercell

# ### Make Superlattice
shape = [2 -1; 1 1]
lattice = makelattice(unitcell, shape)

# ### Plot Lattice
let bravais_lattice = [], a_sites = [], b_sites = []
    extent = [-4, 4, -4, 4]
    within = makewithin(extent)

    for i1 in -5:5, i2 in -5:5
        R = lattice.supercell.latticevectors * [i1, i2]
        within(R...) && push!(bravais_lattice, R)
        for (orb_name, orb_fc) in lattice.supercell.sites
            orb_cc = fract2carte(lattice.supercell, orb_fc)
            r = R + orb_cc
            within(r...) && push!(orb_name[1] == "A" ? a_sites : b_sites, r)
        end
    end
    bravais_lattice = hcat(bravais_lattice...)
    a_sites = hcat(a_sites...)
    b_sites = hcat(b_sites...)

    img = plot(aspect_ratio=1)
    scatter!(bravais_lattice[1,:], bravais_lattice[2,:], color="black", markershape=:circle, markersize=3, label="Bravais (supercell)")
    scatter!(a_sites[1,:], a_sites[2,:], color="red", markerstrokewidth=0, markersize=5, label="A")
    scatter!(b_sites[1,:], b_sites[2,:], color="blue", markerstrokewidth=0, markersize=5, label="B")
    xlims!(extent[1], extent[2])
    ylims!(extent[3], extent[4])
    savefig(img, "example_little_symmetry_kspace_honeycomb/realspace-root3xroot3.svg")
end

# ![](example_little_symmetry_kspace_honeycomb/realspace-root3xroot3.svg)

# ### Symmetries
tsym = FiniteTranslationSymmetry(lattice)
psym = project(PointSymmetryDatabase.get(25), [1 0 0; 0 1 0])
println("point symmetry (Hermann Mauguin): $(psym.hermann_mauguin)")
println("translation symmetry compatible with point symmetry? ", iscompatible(tsym, psym))

println("irreducible representations:")
for tsym_irrep_index in 1:num_irreps(tsym)
    kf = tsym.fractional_momenta[tsym_irrep_index]
    kc = lattice.unitcell.reducedreciprocallatticevectors * kf
    psym_little = little_symmetry(tsym, tsym_irrep_index, psym)
    print("- irrep index: $tsym_irrep_index\n")
    print("  fractional momentum: [")
    for (i, x) in enumerate(kf)
        i != 1 && print(", ")
        print(x)
    end
    print("]\n")
    print("  momentum: [")
    for (i, x) in enumerate(kc)
        i != 1 && print(", ")
        print(x)
    end
    print("]\n")
    print("  is point symmetry compatible with momentum: $(iscompatible(tsym, tsym_irrep_index, psym))\n")
    print("  little point symmetry: $(psym_little.hermann_mauguin)\n")
end

# ### Plot Momentum Space
let
    extent = [-1.1, 1.1, -1.4, 1.4]
    within = makewithin(extent)

    reciprocallatticepoints = []
    kpoints = []
    littlegroupnames = String[]
    generatornames = String[]
    kpointnames = String[]

    for tsym_irrep_index in 1:num_irreps(tsym)
        psym_little = little_symmetry(tsym, tsym_irrep_index, psym)
        kf = tsym.fractional_momenta[tsym_irrep_index]
        kc = lattice.unitcell.reducedreciprocallatticevectors * kf
        for i1 in -2:2, i2 in -2:2
            G = lattice.unitcell.reducedreciprocallatticevectors * [i1, i2]
            within(G...) && push!(reciprocallatticepoints, G)
            k = G + kc
            if within(k...)
                push!(kpoints, k)
                push!(littlegroupnames, psym_little.hermann_mauguin)
                push!(generatornames, join(simplifyname.(psym_little.element_names[psym_little.generators]), "\n"))
                push!(kpointnames, "$tsym_irrep_index")
            end
        end
    end
    reciprocallatticepoints = hcat(reciprocallatticepoints...)
    kpoints = hcat(kpoints...)

    img = plot(size=(400, 500), xlims=(extent[1]-0.1, extent[2]+0.1), ylims=(extent[3]-0.1, extent[4]+0.1), aspect_ratio=1)
    scatter!(reciprocallatticepoints[1,:], reciprocallatticepoints[2,:],
             markershape=:circle, markersize=6, markercolor=RGBA(1,1,1,0), markerstrokecolor=RGBA(1,0,0,1), label="reciprocal lattice")
    scatter!(kpoints[1,:], kpoints[2,:],
             markershape=:circle, markersize=4, markercolor=RGBA(0,0,1,0.5), markerstrokecolor=RGBA(0,0,1,0.5),
             series_annotations=[Plots.text("k[$x]", 8, :left, :bottom) for x in kpointnames], label="momentum points")
    scatter!(kpoints[1,:], kpoints[2,:],
             markershape=:circle, markersize=0, color=RGBA(0,0,0,0),
             series_annotations=[Plots.text("$x ", 8, :right, :top) for x in littlegroupnames], label="")
    scatter!(kpoints[1,:], kpoints[2,:],
             markershape=:circle, markersize=0, color=RGBA(0,0,0,0),
             series_annotations=[Plots.text("$x", 6, :left, :top) for x in generatornames], label="")
    savefig(img, "example_little_symmetry_kspace_honeycomb/momentumspace-root3xroot3.svg")
end

# ![](example_little_symmetry_kspace_honeycomb/momentumspace-root3xroot3.svg)

# ## 2√3 × 2√3 Supercell

# ### Make Superlattice
shape = [2 2; -2 4]
lattice = makelattice(unitcell, shape)

# ### Plot Lattice
let bravais_lattice = [], a_sites = [], b_sites = []
    extent = [-8, 8, -8, 8]
    within = makewithin(extent)

    for i1 in -5:5, i2 in -5:5
        R = lattice.supercell.latticevectors * [i1, i2]
        within(R...) && push!(bravais_lattice, R)
        for (orb_name, orb_fc) in lattice.supercell.sites
            orb_cc = fract2carte(lattice.supercell, orb_fc)
            r = R + orb_cc
            within(r...) && push!(orb_name[1] == "A" ? a_sites : b_sites, r)
        end
    end
    bravais_lattice = hcat(bravais_lattice...)
    a_sites = hcat(a_sites...)
    b_sites = hcat(b_sites...)

    img = plot(aspect_ratio=1)
    scatter!(bravais_lattice[1,:], bravais_lattice[2,:], color="black", markershape=:circle, markersize=5, label="Bravais (supercell)")
    scatter!(a_sites[1,:], a_sites[2,:], color="red", markerstrokewidth=0, markersize=3, label="A")
    scatter!(b_sites[1,:], b_sites[2,:], color="blue", markerstrokewidth=0, markersize=3, label="B")
    xlims!(extent[1], extent[2])
    ylims!(extent[3], extent[4])
    savefig(img, "example_little_symmetry_kspace_honeycomb/realspace-2root3x2root3.svg")
end

# ![](example_little_symmetry_kspace_honeycomb/realspace-2root3x2root3.svg)

# ### Symmetries
tsym = FiniteTranslationSymmetry(lattice)
psym = project(PointSymmetryDatabase.get(25), [1 0 0; 0 1 0])
println("point symmetry (Hermann Mauguin): $(psym.hermann_mauguin)")
println("translation symmetry compatible with point symmetry? ", iscompatible(tsym, psym))

println("irreducible representations:")
for tsym_irrep_index in 1:num_irreps(tsym)
    kf = tsym.fractional_momenta[tsym_irrep_index]
    kc = lattice.unitcell.reducedreciprocallatticevectors * kf
    psym_little = little_symmetry(tsym, tsym_irrep_index, psym)
    print("- irrep index: $tsym_irrep_index\n")
    print("  fractional momentum: [")
    for (i, x) in enumerate(kf)
        i != 1 && print(", ")
        print(x)
    end
    print("]\n")
    print("  momentum: [")
    for (i, x) in enumerate(kc)
        i != 1 && print(", ")
        print(x)
    end
    print("]\n")
    print("  is point symmetry compatible with momentum: $(iscompatible(tsym, tsym_irrep_index, psym))\n")
    print("  little point symmetry: $(psym_little.hermann_mauguin)\n")
end

# ### Plot Momentum Space
let
    extent = [-1.2, 1.2, -1.6, 1.6]
    extent = [-1.1, 1.1, -1.4, 1.4]
    within = makewithin(extent)

    reciprocallatticepoints = []
    kpoints = []
    littlegroupnames = String[]
    generatornames = String[]
    kpointnames = String[]

    for tsym_irrep_index in 1:num_irreps(tsym)
        psym_little = little_symmetry(tsym, tsym_irrep_index, psym)
        kf = tsym.fractional_momenta[tsym_irrep_index]
        kc = lattice.unitcell.reducedreciprocallatticevectors * kf
        for i1 in -2:2, i2 in -2:2
            G = lattice.unitcell.reducedreciprocallatticevectors * [i1, i2]
            within(G...) && push!(reciprocallatticepoints, G)
            k = G + kc
            if within(k...)
                push!(kpoints, k)
                push!(littlegroupnames, psym_little.hermann_mauguin)
                push!(generatornames, join(simplifyname.(psym_little.element_names[psym_little.generators]), "\n"))
                push!(kpointnames, "$tsym_irrep_index")
            end
        end
    end
    reciprocallatticepoints = hcat(reciprocallatticepoints...)
    kpoints = hcat(kpoints...)

    img = plot(size=(600, 750), xlims=(extent[1]-0.1, extent[2]+0.1), ylims=(extent[3]-0.1, extent[4]+0.1), aspect_ratio=1)
    scatter!(reciprocallatticepoints[1,:], reciprocallatticepoints[2,:],
             markershape=:circle, markersize=6, markercolor=RGBA(1,1,1,0), markerstrokecolor=RGBA(1,0,0,1), label="reciprocal lattice")
    scatter!(kpoints[1,:], kpoints[2,:],
             markershape=:circle, markersize=4, markercolor=RGBA(0,0,1,0.5), markerstrokecolor=RGBA(0,0,1,0.5),
             series_annotations=[Plots.text("k[$x]", 8, :left, :bottom) for x in kpointnames], label="momentum points")
    scatter!(kpoints[1,:], kpoints[2,:],
             markershape=:circle, markersize=0, color=RGBA(0,0,0,0),
             series_annotations=[Plots.text("$x ", 8, :right, :top) for x in littlegroupnames], label="")
    scatter!(kpoints[1,:], kpoints[2,:],
             markershape=:circle, markersize=0, color=RGBA(0,0,0,0),
             series_annotations=[Plots.text("$x", 6, :left, :top) for x in generatornames], label="")
    savefig(img, "example_little_symmetry_kspace_honeycomb/momentumspace-2root3x2root3.svg")
end

# ![](example_little_symmetry_kspace_honeycomb/momentumspace-2root3x2root3.svg)

# ## 6 × 6 Supercell

# ### Make Superlattice
shape = [6 0; 0 6]
lattice = makelattice(unitcell, shape)

# ### Plot Lattice
let bravais_lattice = [], a_sites = [], b_sites = []
    extent = [-8, 8, -8, 8]
    within = makewithin(extent)

    for i1 in -5:5, i2 in -5:5
        R = lattice.supercell.latticevectors * [i1, i2]
        within(R...) && push!(bravais_lattice, R)
        for (orb_name, orb_fc) in lattice.supercell.sites
            orb_cc = fract2carte(lattice.supercell, orb_fc)
            r = R + orb_cc
            within(r...) && push!(orb_name[1] == "A" ? a_sites : b_sites, r)
        end
    end
    bravais_lattice = hcat(bravais_lattice...)
    a_sites = hcat(a_sites...)
    b_sites = hcat(b_sites...)

    img = plot(aspect_ratio=1)
    scatter!(bravais_lattice[1,:], bravais_lattice[2,:], color="black", markershape=:circle, markersize=5, label="Bravais (supercell)")
    scatter!(a_sites[1,:], a_sites[2,:], color="red", markerstrokewidth=0, markersize=3, label="A")
    scatter!(b_sites[1,:], b_sites[2,:], color="blue", markerstrokewidth=0, markersize=3, label="B")
    xlims!(extent[1], extent[2])
    ylims!(extent[3], extent[4])
    savefig(img, "example_little_symmetry_kspace_honeycomb/realspace-6x6.svg")
end

# ![](example_little_symmetry_kspace_honeycomb/realspace-6x6.svg)

# ### Symmetries
tsym = FiniteTranslationSymmetry(lattice)
psym = project(PointSymmetryDatabase.get(25), [1 0 0; 0 1 0])
println("point symmetry (Hermann Mauguin): $(psym.hermann_mauguin)")
println("translation symmetry compatible with point symmetry? ", iscompatible(tsym, psym))

println("irreducible representations:")
for tsym_irrep_index in 1:num_irreps(tsym)
    kf = tsym.fractional_momenta[tsym_irrep_index]
    kc = lattice.unitcell.reducedreciprocallatticevectors * kf
    psym_little = little_symmetry(tsym, tsym_irrep_index, psym)
    print("- irrep index: $tsym_irrep_index\n")
    print("  fractional momentum: [")
    for (i, x) in enumerate(kf)
        i != 1 && print(", ")
        print(x)
    end
    print("]\n")
    print("  momentum: [")
    for (i, x) in enumerate(kc)
        i != 1 && print(", ")
        print(x)
    end
    print("]\n")
    print("  is point symmetry compatible with momentum: $(iscompatible(tsym, tsym_irrep_index, psym))\n")
    print("  little point symmetry: $(psym_little.hermann_mauguin)\n")
end

# ### Plot Momentum Space
let
    extent = [-1.2, 1.2, -1.6, 1.6]
    extent = [-1.1, 1.1, -1.4, 1.4]
    within = makewithin(extent)

    reciprocallatticepoints = []
    kpoints = []
    littlegroupnames = String[]
    generatornames = String[]
    kpointnames = String[]

    for tsym_irrep_index in 1:num_irreps(tsym)
        psym_little = little_symmetry(tsym, tsym_irrep_index, psym)
        kf = tsym.fractional_momenta[tsym_irrep_index]
        kc = lattice.unitcell.reducedreciprocallatticevectors * kf
        for i1 in -2:2, i2 in -2:2
            G = lattice.unitcell.reducedreciprocallatticevectors * [i1, i2]
            within(G...) && push!(reciprocallatticepoints, G)
            k = G + kc
            if within(k...)
                push!(kpoints, k)
                push!(littlegroupnames, psym_little.hermann_mauguin)
                push!(generatornames, join(simplifyname.(psym_little.element_names[psym_little.generators]), "\n"))
                push!(kpointnames, "$tsym_irrep_index")
            end
        end
    end
    reciprocallatticepoints = hcat(reciprocallatticepoints...)
    kpoints = hcat(kpoints...)

    img = plot(size=(800, 1000), xlims=(extent[1]-0.1, extent[2]+0.1), ylims=(extent[3]-0.1, extent[4]+0.1), aspect_ratio=1)
    scatter!(reciprocallatticepoints[1,:], reciprocallatticepoints[2,:],
             markershape=:circle, markersize=6, markercolor=RGBA(1,1,1,0), markerstrokecolor=RGBA(1,0,0,1), label="reciprocal lattice")
    scatter!(kpoints[1,:], kpoints[2,:],
             markershape=:circle, markersize=4, markercolor=RGBA(0,0,1,0.5), markerstrokecolor=RGBA(0,0,1,0.5),
             series_annotations=[Plots.text("k[$x]", 8, :left, :bottom) for x in kpointnames], label="momentum points")
    scatter!(kpoints[1,:], kpoints[2,:],
             markershape=:circle, markersize=0, color=RGBA(0,0,0,0),
             series_annotations=[Plots.text("$x ", 8, :right, :top) for x in littlegroupnames], label="")
    scatter!(kpoints[1,:], kpoints[2,:],
             markershape=:circle, markersize=0, color=RGBA(0,0,0,0),
             series_annotations=[Plots.text("$x", 6, :left, :top) for x in generatornames], label="")
    savefig(img, "example_little_symmetry_kspace_honeycomb/momentumspace-6x6.svg")
end

# ![](example_little_symmetry_kspace_honeycomb/momentumspace-6x6.svg)


text:
function update_paramsfromM(M, ϵ, ζ)
    κ = (M+2ϵ) / ϵ
    ζ̂ = ζ / (3κ)
    τ = √(κ) / (√(κ) + 1)
    T = 4 * κ^4 / (1-√(τ))^2
    return κ, ζ̂, τ, T
end

"""
dₖ = solve_tCG(M, x, gradf_k, hessf_x_h, ν=1e-2, ϵ_residual = 1e-5)

Conjugate gradient method to solve the linear problem:
find d s.t. hessf_x_h(d) + gradf_k = 0

An implementation of ``capped CG'' following
    *A Newton-CG algorithm with complexity guarantees for smooth
    unconstrained optimization*, Royer, Neil, Wright.
"""
function solve_tCG_capped(Man, x, gradfₖ, hessf_x_h; ϵ=1e-8, ζ = 1e-8, maxiter=1e5, check_tvectors=false, printlev=0)

    # Short notations
    norm_x(η) = norm(Man, x, η)
    inner_x(η, ξ) = inner(Man, x, η, ξ)
    H(η) = hessf_x_h(η)
    H̄(η) = hessf_x_h(η) + 2ϵ*η

    # Operator norm of H
    M::Float64 = 0.0

    κ, ζ̂, τ, T = update_paramsfromM(M, ϵ, ζ)

    r₀ = gradfₖ
    yⱼ = zero_tangent_vector(Man, x)
    rⱼ, rⱼ_prev = deepcopy(gradfₖ), deepcopy(gradfₖ)
    pⱼ, pⱼ_prev = deepcopy(-gradfₖ), deepcopy(-gradfₖ)
    j = 0

    # ##
    # @show ϵ
    # n = size(x, 1)
    # @show Man.nnz_coords

    # Hess = zeros(10, 10)
    # Hess_bar = zeros(10, 10)

    # for i in 1:n
    #     e_i = zeros(n)
    #     e_i[i] = 1.0

    #     Hess[:, i] = H(e_i)
    #     Hess_bar[:, i] = H̄(e_i)
    # end

    # @show (eigvals(Hess))
    # @show (eigvals(Hess_bar))

    # d = rand(10)
    # @assert norm(Hess*d - H(d)) < 1e-14
    # ##

    if inner_x(pⱼ, H̄(pⱼ)) < ϵ * inner_x(pⱼ, pⱼ)
        return pⱼ, :NegativeCurvature, j
    elseif norm_x(H(pⱼ)) > M * norm_x(pⱼ)
        M = norm_x(H(pⱼ)) / norm_x(pⱼ)
        κ, ζ̂, τ, T = update_paramsfromM(M, ϵ, ζ)
    end

    (printlev>0) && @printf "\nj     norm(rⱼ)             norm(vⱼ)             ⟨vⱼ, hessf(x)[vⱼ]⟩   ν * norm(vⱼ)^2\n"
    while true

        rⱼ_prev = deepcopy(rⱼ)
        # yⱼ_prev = deepcopy(yⱼ)
        pⱼ_prev = deepcopy(pⱼ)

        # Standard CG operations
        αⱼ = inner_x(rⱼ, rⱼ) / inner_x(pⱼ, H̄(pⱼ))
        yⱼ += αⱼ * pⱼ
        rⱼ += αⱼ * H̄(pⱼ)

        #TODO: storing the inner product is enough here
        βⱼ = inner_x(rⱼ, rⱼ) / inner_x(rⱼ_prev, rⱼ_prev)
        pⱼ = -rⱼ + βⱼ * pⱼ_prev
        j += 1

        # Update estimate of H's operator norm
        for v in [pⱼ, yⱼ, rⱼ]
            norm_Hv = norm_x(H(v))
            norm_v = norm_x(v)
            if norm_Hv > M * norm_v
                M = norm_Hv / norm_v
                κ, ζ̂, τ, T = update_paramsfromM(M, ϵ, ζ)
            end
        end

        # Plotting
        if printlev>0
            # @printf "%5i %.10e    %.10e    % .10e     %.10e\n" j norm_x(rⱼ) norm_x(vⱼ) inner(M, x, vⱼ, hessf_x_vⱼ) ν * norm_x(vⱼ)^2
            @printf "%5i %.10e\n" j norm_x(rⱼ)
        end

        # Termination criteria
        if inner_x(yⱼ, H̄(yⱼ)) < ϵ * inner_x(yⱼ, yⱼ)
            return yⱼ, :NegativeCurvature, j
        elseif norm_x(rⱼ) ≤ ζ̂ * norm_x(r₀)
            return yⱼ, :Solution, j
        elseif inner_x(pⱼ, H̄(pⱼ)) < ϵ * inner_x(pⱼ, pⱼ)
            return pⱼ, :NegativeCurvature, j
        elseif norm_x(rⱼ) > √(T) * τ^(j/2) * norm_x(r₀)
            # αⱼ = inner_x(rⱼ, rⱼ) / inner_x(pⱼ, H(pⱼ))
            # yⱼ += αⱼ * pⱼ
            # TODO : implement last item, by executing again iterations of CG (instead of storing)
            @warn "Weird unfound negative curvature direction. Not implemented."
            @assert false
        end

        if j >= 400
            return -yⱼ * sign(inner_x(yⱼ, gradfₖ)), :MaxIter, j
        end
        # (printlev>0) && @printf "%5i %.10e    %.10e    % .10e     %.10e\n" j norm_x(rⱼ) norm_x(vⱼ) inner(M, x, vⱼ, hessf_x_vⱼ) ν * norm_x(vⱼ)^2
        # if norm_x(rⱼ) < ϵ_residual || inner(M, x, vⱼ, hessf_x_vⱼ) < ν * norm_x(vⱼ)^2 || j > maxiter

        # j += 1
    end

    return
end


text:
using DiffEqCallbacks
using OrdinaryDiffEq
using SimulationLogs
using Test

function lotka!(du, u, p, t)
    @log t
    @log x, y = u
    @log α, β, δ, γ = p

    @log total_population = x + y

    du[1] = α*x - β*x*y
    du[2] = δ*x*y - γ*y
end
function lotka(u, p, t)
    du = similar(u)
    lotka!(du, u, p , t)
    return du
end

lotka_sols = let
    p = [2/3, 4/3, 1, 1]
    u0 = [0.9, 0.9]
    tspan = (0.0, 10.0)

    callback = PresetTimeCallback(1:10, integrator->(integrator.p[end] += 0.01))
    prob_oop = ODEProblem(lotka, u0, tspan, p)
    prob_iip = ODEProblem(lotka!, u0, tspan, p)

    sol_oop = logged_solve(prob_oop, Tsit5())
    sol_iip = logged_solve(prob_iip, Tsit5())
    sol_oop_cb = logged_solve(prob_oop, Tsit5(); callback=callback)
    sol_iip_cb = logged_solve(prob_iip, Tsit5(); callback=callback)
    sol_oop_cbset = logged_solve(prob_oop, Tsit5(); callback=CallbackSet(callback))
    sol_iip_cbset = logged_solve(prob_iip, Tsit5(); callback=CallbackSet(callback))

    sols_iip = [sol_iip, sol_iip_cb, sol_iip_cb]
    sols_oop = [sol_oop, sol_oop_cb, sol_oop_cb]

    (iip=sols_iip, oop=sols_oop)
end

@testset "Lotka example" begin
    test_all(f, i=:) = all(f, lotka_sols.iip[i]) && all(f, lotka_sols.oop[i])
    test_cb(f) = test_all(f, 2:3)

    @testset "Logged values" begin
        @test test_all(sol->sol.log.γ[1]==1)
        @test test_all(sol->sol.log.x==sol[1,:])
        @test test_all(sol->sol.log.y==sol[2,:])
        @test test_all(sol->sol.log.t==sol.t)

        @test test_all(sol->all(sol.log.α .== sol.prob.p[1]))
        @test test_cb(sol->sol.log.γ[end] == 1.1)
    end

    this_sol = lotka_sols.iip[1]

    @testset "Indexing" begin
        @test this_sol[1] == this_sol.u[1]
        @test this_sol[:] == this_sol.u
        @test this_sol[:, :] == reduce(hcat, this_sol.u)
    end

    @testset "Interpolation" begin
        @test this_sol(10) ≈ this_sol.u[end]
    end

    @testset "Properties" begin
        @test this_sol.log isa SimulationLog
        @test this_sol.alg == this_sol.sol.alg
        @test Set(propertynames(this_sol)) == Set((:log, propertynames(this_sol.sol)...))
    end

    @testset "Printing" begin
        sol_string = @capture_out(print(this_sol))
        sol_sol_string = @capture_out(print(this_sol.sol))
        @test contains(sol_string, sol_sol_string)
    end
end

max_stars_repo_name: jonniedie/SimulationLogs

max_issues_repo_name: jonniedie/SimulationLogs

max_forks_repo_name: jonniedie/SimulationLogs


text:
# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ
function ITensors.MPO(::Model"ising", s; J, h)
  N = length(s)
  a = OpSum()
  for n in 1:(N - 1)
    a .+= -J, "X", n, "X", n + 1
  end
  for n in 1:N
    a .+= -h, "Z", n
  end
  return splitblocks(linkinds, MPO(a, s))
end

# H = -J Σⱼ XⱼXⱼ₊₁ - h Σⱼ Zⱼ
function ITensors.OpSum(::Model"ising", n1, n2; J, h)
  opsum = OpSum()
  if J != 0
    opsum += -J, "X", n1, "X", n2
  end
  if h != 0
    opsum += -h / 2, "Z", n1
    opsum += -h / 2, "Z", n2
  end
  return opsum
end

# H = -J X₁X₂ - h Z₁
# XXX: use `op` instead of `ITensor`
function ITensors.ITensor(::Model"ising", s1::Index, s2::Index; J, h)
  # -J * op("X", s1) * op("X", s2) - h * op("Z", s1) * op("Id", s2)
  opsum = OpSum()
  n = 1
  opsum += -J, "X", n, "X", n + 1
  opsum += -h, "Z", n
  return prod(MPO(opsum, [s1, s2]))
end

# P. Pfeuty, The one-dimensional Ising model with a transverse field, Annals of Physics 57, p. 79 (1970)
function reference(::Model"ising", ::Observable"energy"; h, J=1.0)
  (J == h) && return -4 / π
  f(k; λ) = sqrt(1 + λ^2 + 2λ * cos(k))
  return -h / (2π * J) * ∫(k -> f(k; λ=(J / h)), -π, π)
end

# https://journals.aps.org/pr/abstract/10.1103/PhysRev.65.117
# https://en.wikipedia.org/wiki/Ising_model#Onsager's_exact_solution
function reference(::Model"ising_classical", ::Observable"free_energy"; β, J=1.0)
  βJ = β * J
  function f(θ₁, θ₂)
    return log(cosh(2βJ)^2 - sinh(2βJ) * cos(θ₁) - sinh(2βJ) * cos(θ₂))
  end
  I = ∫(θ₂ -> ∫(θ₁ -> f(θ₁, θ₂), 0, 2π), 0, 2π)
  return -(log(2) + I / (8π^2)) / β
end

function reference(::Model"ising_classical", ::Observable"critical_inverse_temperature")
  # log1p(x) = log(1 + x)
  return log1p(sqrt(2)) / 2
end


text:
@assert typeof(88) == Int64
for t in [Int64, Integer, Signed, Real, Number, Any]
  @assert isa(88, t)
end

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
using CairoMakie, Distributions, Plots, PSIS
using Documenter

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

makedocs(;
    modules=[PSIS],
    authors="Seth Axen <seth.axen@gmail.com> and contributors",
    repo="https://github.com/arviz-devs/PSIS.jl/blob/{commit}{path}#{line}",
    sitename="PSIS.jl",
    format=Documenter.HTML(; prettyurls=get(ENV, "CI", "false") == "true", assets=String[]),
    pages=[
        "Home" => "index.md",
        "Plotting" => "plotting.md",
        "API" => "api.md",
        "Internal" => "internal.md",
    ],
)

deploydocs(; repo="github.com/arviz-devs/PSIS.jl.git", devbranch="main")


text:
include("../roblib.jl")

function mark_ladder(r)
    num_steps = through_rectangles_into_angle(r, (Down, Left))

    cells_to_mark = mark_and_enumerate!(r, Right) - 1
    while isborder(r, Up) == false
        movements!(r, Left)
        move!(r, Up)
        putmarkers!(r, Right, cells_to_mark)
        cells_to_mark -= 1
    end

    through_rectangles_into_angle(r, (Down, Left))
    movements!(r, (Up, Right), reverse(num_steps))
end

"""
    mark_and_enumerate(r, side::HorizonSide)

-- Ставит маркеры в указанном направлении, пока не наткнется на ограду, а также считает количетво шагов сделанных в эту сторону и возвращает их
"""
function mark_and_enumerate!(r, side::HorizonSide)
    steps = 0

    putmarker!(r)
    while isborder(r, side) == false
        move!(r, side)
        putmarker!(r)
        steps += 1
    end

    return steps
end


text:
using MultivariateSeries

r  = 4
w0 = rand(r)
A0 = rand(5,r)
B0 = rand(3,r)
C0 = rand(4,r)

T0 = tensor(w0, A0, B0, C0)
w, A, B, C = decompose(T0, eps_rkf(1.e-10), mode=2)

T = tensor(w, A, B, C)
@assert norm(T-T0)< 1.e-6


text:
__precompile__(true)

module CompressedSensing

import JuMP
using Clp
using LossFunctions

export l1reconstruct
export evaluateloss

"""
    l1reconstruct(x, A)

Return estimated x by ``\ell_1`` reconstruction.
"""
function l1reconstruct(unknown_x::Array{Float64, 1}, A::Array{Float64, 2})
    N = length(unknown_x)
    if size(A, 2) ≠ N
        error("doesn't match the size of x and A")
    end
    M = size(A, 1)
    y = A * unknown_x

    l1 = JuMP.Model(solver=ClpSolver())
    JuMP.@variable(l1, t[1:N])
    JuMP.@variable(l1, x[1:N] >= 0.0)

    JuMP.@objective(l1, Min, sum(t[n] for n=1:N))

    JuMP.@constraint(l1, [n=1:N], -t[n] <= x[n])
    JuMP.@constraint(l1, [n=1:N], x[n] <= t[n])
    JuMP.@constraint(l1, [m=1:M], y[m] == sum(A[m,n] * x[n] for n=1:N))
    status = JuMP.solve(l1)
    JuMP.getvalue(x)
end

"""
    evaluateloss(x, A; alg, loss)

Return evaluated loss by an algorithm with a loss.
Default alg=l1reconstruct, loss=L2DisctLoss().
"""
function evaluateloss(x::Array{Float64, 1}, A::Array{Float64, 2}; alg=l1reconstruct, loss=L2DistLoss())
    estimatedx = alg(x, A)
    loss(estimated_x, A)
end

end # module


text:
using Documenter
using SimJulia

makedocs(
  modules = [SimJulia],
  clean   = true,
  format   = Documenter.Formats.HTML,
  sitename = "SimJulia.jl",
  pages    = [
    "Home" => "index.md",
    "Intro" => "intro.md",
    "Manual" => "topics.md",
    "Library" => "api.md"
    ]

)

deploydocs(
  repo = "github.com/BenLauwens/SimJulia.jl.git",
  julia  = "0.5",
  osname = "linux",
  target = "build",
  deps = nothing,
  make = nothing,
)


text:
import Base: isempty, ∈, ∪

export UnionSet,
       swap

# ========================================
# Binary set union
# ========================================

"""
    UnionSet{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}}

Type that represents the set union of two convex sets.

### Fields

- `X` -- convex set
- `Y` -- convex set
"""
struct UnionSet{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}}
    X::S1
    Y::S2

    # default constructor with dimension check
    function UnionSet(X::S1, Y::S2) where{N<:Real, S1<:LazySet{N}, S2<:LazySet{N}}
        @assert dim(X) == dim(Y) "sets in a union must have the same dimension"
        return new{N, S1, S2}(X, Y)
    end
end

isoperationtype(::Type{<:UnionSet}) = true
isconvextype(::Type{<:UnionSet}) = false

# EmptySet is the neutral element for UnionSet
@neutral(UnionSet, EmptySet)

# Universe is the absorbing element for UnionSet
@absorbing(UnionSet, Universe)

"""
    ∪

Alias for `UnionSet`.
"""
∪(X::LazySet, Y::LazySet) = UnionSet(X, Y)

"""
    swap(cup::UnionSet)

Return a new `UnionSet` object with the arguments swapped.

### Input

- `cup` -- union of two convex sets

### Output

A new `UnionSet` object with the arguments swapped.
"""
function swap(cup::UnionSet)
    return UnionSet(cup.Y, cup.X)
end

"""
    dim(cup::UnionSet)

Return the dimension of the set union of two convex sets.

### Input

- `cup` -- union of two convex sets

### Output

The ambient dimension of the union of two convex sets.
"""
function dim(cup::UnionSet)
    return dim(cup.X)
end

"""
    σ(d::AbstractVector{N}, cup::UnionSet{N}; [algorithm]="support_vector") where {N<:Real}

Return the support vector of the union of two convex sets in a given direction.

### Input

- `d`         -- direction
- `cup`       -- union of two convex sets
- `algorithm` -- (optional, default: "support_vector"): the algorithm to compute
                 the support vector; if "support_vector", use the support
                 vector of each argument; if "support_function" use the support
                 function of each argument and evaluate the support vector of only
                 one of them

### Output

The support vector in the given direction.

### Algorithm

The support vector of the union of two convex sets ``X`` and ``Y`` can be obtained
as the vector that maximizes the support function of either ``X`` or ``Y``, i.e.
it is sufficient to find the ``\\argmax(ρ(d, X), ρ(d, Y)])`` and evaluate its support
vector.

The default implementation, with option `algorithm="support_vector"`, computes
the support vector of ``X`` and ``Y`` and then compares the support function using
a dot product. If it happens that the support function can be more efficiently
computed (without passing through the support vector), consider using the alternative
`algorithm="support_function"` implementation, which evaluates the support function
of each set directly and then calls only the support vector of either ``X`` *or*
``Y``.
"""
function σ(d::AbstractVector{N}, cup::UnionSet{N}; algorithm="support_vector") where {N<:Real}
    X, Y = cup.X, cup.Y
    if algorithm == "support_vector"
        σX, σY = σ(d, X), σ(d, Y)
        return dot(d, σX) > dot(d, σY) ? σX : σY
    elseif algorithm == "support_function"
        m = argmax([ρ(d, X), ρ(d, Y)])
        return m == 1 ? σ(d, X) : σ(d, Y)
    else
        error("algorithm $algorithm for the support vector of a `UnionSet` is unknown")
    end
end

"""
    ρ(d::AbstractVector{N}, cup::UnionSet{N}) where {N<:Real}

Return the support function of the union of two convex sets in a given direction.

### Input

- `d`   -- direction
- `cup` -- union of two convex sets

### Output

The support function in the given direction.

### Algorithm

The support function of the union of two convex sets ``X`` and ``Y`` is the
maximum of the support functions of ``X`` and ``Y``.
"""
function ρ(d::AbstractVector{N}, cup::UnionSet{N}) where {N<:Real}
    X, Y = cup.X, cup.Y
    return max(ρ(d, X), ρ(d, Y))
end

"""
    an_element(cup::UnionSet{N}) where {N<:Real}

Return some element of a union of two convex sets.

### Input

- `cup` -- union of two convex sets

### Output

An element in the union of two convex sets.

### Algorithm

We use `an_element` on the first wrapped set.
"""
function an_element(cup::UnionSet{N}) where {N<:Real}
    return an_element(cup.X)
end

"""
    ∈(x::AbstractVector{N}, cup::UnionSet{N}) where {N<:Real}

Check whether a given point is contained in a union of two convex sets.

### Input

- `x`   -- point/vector
- `cup` -- union of two convex sets

### Output

`true` iff ``x ∈ cup``.
"""
function ∈(x::AbstractVector{N}, cup::UnionSet{N}) where {N<:Real}
    return x ∈ cup.X || x ∈ cup.Y
end

"""
    isempty(cup::UnionSet)

Check whether a union of two convex sets is empty.

### Input

- `cup` -- union of two convex sets

### Output

`true` iff the union is empty.
"""
function isempty(cup::UnionSet)
    return isempty(cup.X) && isempty(cup.Y)
end

"""
    isbounded(cup::UnionSet)

Determine whether a union of two convex sets is bounded.

### Input

- `cup` -- union of two convex sets

### Output

`true` iff the union is bounded.
"""
function isbounded(cup::UnionSet)
    return isbounded(cup.X) && isbounded(cup.Y)
end


text:
using TopOpt, LinearAlgebra, StatsFuns
using TopOpt.TrussTopOptProblems
using TopOptMakie
using VTKDataIO
using PyCall
using VTKDataTypes
using DashVtk

using Dash
using DashHtmlComponents

#dash = pyimport("dash")
#dash_html = pyimport("dash_html_components")
dash_vtk = pyimport("dash_vtk")

# 2D
ndim = 2
node_points, elements, mats, crosssecs, fixities, load_cases =
    parse_truss_json(joinpath(@__DIR__, "tim_$(ndim)d.json"));
ndim, nnodes, ncells = length(node_points[1]), length(node_points), length(elements)
loads = load_cases["0"]
problem = TrussProblem(
    Val{:Linear}, node_points, elements,
    loads, fixities, mats, crosssecs,
);

xmin = 0.0001 # minimum density
x0 = fill(1.0, ncells) # initial design
p = 4.0 # penalty
V = 0.5 # maximum volume fraction

solver = FEASolver(Displacement, Direct, problem, xmin = xmin)
comp = TopOpt.Compliance(problem, solver)

function obj(x)
    # minimize compliance
    return comp(x)
end
function constr(x)
    # volume fraction constraint
    return sum(x) / length(x) - V
end

m = Model(obj)
addvar!(m, zeros(length(x0)), ones(length(x0)))
Nonconvex.add_ineq_constraint!(m, constr)

options = MMAOptions(
    maxiter=1000, tol = Tolerance(kkt = 1e-4, f = 1e-4),
)
TopOpt.setpenalty!(solver, p)
@time r = Nonconvex.optimize(
    m, MMA87(dualoptimizer = ConjugateGradient()),
    x0, options = options,
);

@show obj(r.minimizer)
@show constr(r.minimizer)

topology = r.minimizer;

function fg(topology, problem)
    mesh = VTKUnstructuredData(problem)
    topology = round.(topology)
    inds = findall(isequal(0), topology)
    deleteat!(mesh.cell_connectivity, inds)
    deleteat!(mesh.cell_types, inds)
    return mesh
end

pymesh = PyVTK(fg(topology, problem))
dash_vtk_utils = pyimport("dash_vtk.utils")
dash_vtk_utils.to_mesh_state
state = dash_vtk_utils.to_mesh_state(pymesh)

content = vtk_view([
        vtk_geometryrepresentation([
        vtk_mesh(state=state)
    ]),
])

app = dash()

app.layout = html_div() do
    html_div(
    style = Dict(
        "width" => "100%",
        "height" => "400px"
    ),
    children=[content],
    )
end

run_server(app, "0.0.0.0", debug = true)