text:
export Poly

function JuMP.value(t::AbstractTerm{<:JuMP.AbstractJuMPScalar})
    JuMP.value(coefficient(t)) * monomial(t)
end
function JuMP.value(p::AbstractPolynomialLike{<:JuMP.AbstractJuMPScalar})
    polynomial(JuMP.value.(terms(p)), MultivariatePolynomials.SortedUniqState())
end

abstract type AbstractPoly end

"""
    struct Poly{PB<:AbstractPolynomialBasis} <: AbstractPoly
        polynomial_basis::PB
    end

Polynomial variable ``v^\\top p`` where ``v`` is a vector of new decision variables and ``p`` is a vector of polynomials for the basis `polynomial_basis`.
"""
struct Poly{PB<:AbstractPolynomialBasis} <: AbstractPoly
    polynomial_basis::PB
end
Poly(x::AbstractVector{<:MultivariatePolynomials.AbstractPolynomialLike}) = Poly(MonomialBasis(x))

# Macro
function cvarchecks(_error::Function, info; extra_kwargs...)
    for (kwarg, _) in extra_kwargs
        _error("Unrecognized keyword argument $kwarg")
    end
    if !isnan(info.start)
        _error("Polynomial variable declaration does not support the start keyword argument.")
    end
    if info.has_fix
        _error("Polynomial variable declaration does not support the form ... == value.")
    end
    if info.has_lb && info.has_ub
        _error("Polynomial variable declaration does not support the form lb <= ... <= ub. Use ... >= 0 and separate constraints instead.")
    end
    if info.has_lb
        _error("Polynomial variable declaration does not support the form ... >= lb.")
    end
    if info.has_ub
        _error("Polynomial variable declaration does not support the form ... <= ub.")
    end
end
function _warnbounds(_error, p::AbstractPoly, info) end
function _warnbounds(_error, p::Poly, info)
    if info.has_lb
        _error("Free polynomial variable declaration does not support the form ... >= 0, use SOSPoly(x) instead of Poly(x) to create Sum of Squares polynomials. Note that SOSPoly(x) creates the polynomial x^T Q x with Q symmetric positive semidefinite while Poly(x) creates the polynomial a^T x so the meaning of the vector of monomial x changes from Poly to SOSPoly.")
    end
end

struct Variable{PT<:AbstractPoly} <: JuMP.AbstractVariable
    p::PT
    binary::Bool
    integer::Bool
end
function JuMP.build_variable(_error::Function, info::JuMP.VariableInfo, p::AbstractPoly; extra_kwargs...)
    cvarchecks(_error, info; extra_kwargs...)
    _warnbounds(_error, p, info)
    return Variable(p, info.binary, info.integer)
end

# Free polynomial
function JuMP.variable_type(model::JuMP.AbstractModel, p::Poly)
    return polytype(model, p, p.polynomial_basis)
end
function polytype(model::JuMP.AbstractModel, ::Poly, pb::AbstractPolynomialBasis)
    return MultivariatePolynomials.polynomialtype(pb, JuMP.VariableRef)
end

function JuMP.add_variable(model::JuMP.AbstractModel, v::Variable{<:Poly},
                           name::String="")
    function _newvar(i)
        vref = VariableRef(model)
        if v.binary
            JuMP.set_binary(vref)
        end
        if v.integer
            JuMP.set_integer(vref)
        end
        return vref
    end
    return polynomial(_newvar, v.p.polynomial_basis)
end


text:
# =============================
# Successor for discrete system
# =============================

"""
    successor(system::DiscreteIdentitySystem, x::AbstractVector)

Return the successor state of a `DiscreteIdentitySystem`.

### Input

- `system` -- `DiscreteIdentitySystem`
- `x`      -- state (it should be any vector type)

### Output

The same state `x`.
"""
function successor(system::DiscreteIdentitySystem, x::AbstractVector)
    !_is_conformable_state(system, x) && _argument_error(:x)
    return x
end

"""
    successor(system::ConstrainedDiscreteIdentitySystem, x::AbstractVector;
              [check_constraints]=true)

Return the successor state of a `ConstrainedDiscreteIdentitySystem`.

### Input

- `system`            -- `ConstrainedDiscreteIdentitySystem`
- `x`                 -- state (it should be any vector type)
- `check_constraints` -- (optional, default: `true`) check if the state belongs to
                         the state set

### Output

The same state `x`.
"""
function successor(system::ConstrainedDiscreteIdentitySystem, x::AbstractVector;
                   check_constraints::Bool=true)
    !_is_conformable_state(system, x) && _argument_error(:x)
    if check_constraints
        !_in_stateset(system, x) && _argument_error(:x,:X)
    end
    return x
end

"""
    successor(system::AbstractDiscreteSystem, x::AbstractVector;
              [check_constraints]=true)

Return the successor state of an `AbstractDiscreteSystem`.

### Input

- `system`            -- `AbstractDiscreteSystem`
- `x`                 -- state (it should be any vector type)
- `check_constraints` -- (optional, default: `true`) check if the state belongs to
                         the state set

### Output

The result of applying the system to state `x`.
"""
successor(system::AbstractDiscreteSystem, x::AbstractVector; kwargs...) =
    _instantiate(system, x; kwargs...)

"""
    successor(system::AbstractDiscreteSystem, x::AbstractVector, u::AbstractVector;
              [check_constraints]=true)

Return the successor state of an `AbstractDiscreteSystem`.

### Input

- `system`            -- `AbstractDiscreteSystem`
- `x`                 -- state (it should be any vector type)
- `u`                 -- input (it should be any vector type) or noise, if `system`
                         is not controlled
- `check_constraints` -- (optional, default: `true`) check if the state belongs to
                         the state set

### Output

The result of applying the system to state `x` and input `u`.

### Notes

If the system is not controlled but noisy, the input `u` is interpreted as noise.
    """
successor(system::AbstractDiscreteSystem, x::AbstractVector, u::AbstractVector; kwargs...) =
    _instantiate(system, x, u; kwargs...)

"""
    successor(system::AbstractDiscreteSystem,
              x::AbstractVector, u::AbstractVector, w::AbstractVector; [check_constraints]=true)

Return the successor state of an `AbstractDiscreteSystem`.

### Input

- `system`            -- `AbstractDiscreteSystem`
- `x`                 -- state (it should be any vector type)
- `u`                 -- input (it should be any vector type)
- `w`                 -- noise (it should be any vector type)
- `check_constraints` -- (optional, default: `true`) check if the state belongs to
                         the state set

### Output

The result of applying the system to state `x`, input `u` and noise `w`.
"""
successor(system::AbstractDiscreteSystem, x::AbstractVector, u::AbstractVector, w::AbstractVector; kwargs...) =
    _instantiate(system, x, u, w; kwargs...)

ext: jl


text:
# # Custom solvers
# In this example, we show how to define custom solvers. Our system
# will again be silicon, because we are not very imaginative
using DFTK, LinearAlgebra

a = 10.26
lattice = a / 2 * [[0 1 1.];
                   [1 0 1.];
                   [1 1 0.]]
Si = ElementPsp(:Si, psp=load_psp("hgh/lda/Si-q4"))
atoms = [Si => [ones(3)/8, -ones(3)/8]]

## We take very (very) crude parameters
model = model_LDA(lattice, atoms)
kgrid = [1, 1, 1]
Ecut = 5
basis = PlaneWaveBasis(model, Ecut; kgrid=kgrid);

# We define our custom fix-point solver: simply a damped fixed-point
function my_fp_solver(f, x0, max_iter; tol)
    mixing_factor = .7
    x = x0
    fx = f(x)
    for n = 1:max_iter
        inc = fx - x
        if norm(inc) < tol
            break
        end
        x = x + mixing_factor * inc
        fx = f(x)
    end
    (fixpoint=x, converged=norm(fx-x) < tol)
end;

# Our eigenvalue solver just forms the dense matrix and diagonalizes
# it explicitly (this only works for very small systems)
function my_eig_solver(A, X0; maxiter, tol, kwargs...)
    n = size(X0, 2)
    A = Array(A)
    E = eigen(A)
    λ = E.values[1:n]
    X = E.vectors[:, 1:n]
    (λ=λ, X=X, residual_norms=[], iterations=0, converged=true, n_matvec=0)
end;

# Finally we also define our custom mixing scheme. It will be a mixture
# of simple mixing (for the first 2 steps) and than default to Kerker mixing.
struct MyMixing
    α  # Damping parameter
end
MyMixing() = MyMixing(0.7)

function DFTK.mix(mixing::MyMixing, basis, ρin::RealFourierArray, ρout::RealFourierArray; n_iter, kwargs...)
    if n_iter <= 2
        ## Just do simple mixing
        ρin + mixing.α * (ρout - ρin)
    else
        ## Use the KerkerMixing from DFTK
        DFTK.mix(KerkerMixing(α=mixing.α), basis, ρin, ρout; kwargs...)
    end
end

# That's it! Now we just run the SCF with these solvers
scfres = self_consistent_field(basis;
                               tol=1e-8,
                               solver=my_fp_solver,
                               eigensolver=my_eig_solver,
                               mixing=MyMixing());
# Note that the default convergence criterion is on the difference of
# energy from one step to the other; when this gets below `tol`, the
# "driver" `self_consistent_field` artificially makes the fixpoint
# solver think it's converged by forcing `f(x) = x`. You can customize
# this with the `is_converged` keyword argument to
# `self_consistent_field`.


text:
@testset "pool" begin
    @testset "GlobalPool" begin
        n = 10
        X = rand(16, n)
        g = GNNGraph(random_regular_graph(n, 4))
        p = GlobalPool(+)
        @test p(g, X) ≈ NNlib.scatter(+, X, ones(Int, n))
    end

    @testset "TopKPool" begin
        N = 10
        k, in_channel = 4, 7
        X = rand(in_channel, N)
        for T = [Bool, Float64]
            adj = rand(T, N, N)
            p = TopKPool(adj, k, in_channel)
            @test eltype(p.p) === Float32
            @test size(p.p) == (in_channel,)
            @test eltype(p.Ã) === T
            @test size(p.Ã) == (k, k)
            y = p(X)
            @test size(y) == (in_channel, k)
        end
    end

    @testset "topk_index" begin
        X = [8,7,6,5,4,3,2,1]
        @test topk_index(X, 4) == [1,2,3,4]
        @test topk_index(X', 4) == [1,2,3,4]
    end
end


text:
using Ragel
using BufferedStreams

if VERSION >= v"0.5-"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

testfilepath(filename) = Pkg.dir("Ragel", "test", filename)

include("simple.jl")

@testset "simple parser" begin
    lines = map(ASCIIString, open(testfilepath("simple.txt"), Simple))
    @test lines == ["foo", "bar", "baz123"]
end


text:
function determine_qdim(equation::Union{SDE,PSDE,SPSDE})
    nd = equation.d
    ns = equation.ns
    ni = equation.n

    if nd==ns==ni==1
        NQ = 1
    elseif ns==ni==1
        NQ = 2
    elseif ns==1 || ni==1
        NQ = 3
    else
        NQ = 4
    end

    return NQ
end


text:
using LinQu, Test, Random
import LinQu.data # data() is extended

# define customized QGate
struct RandomSingleGate <: QGate
	rng::MersenneTwister
	seed::Int
	maxangle::Number

	RandomSingleGate(rng::MersenneTwister, seed::Int, maxangle::Number) = new(rng, seed, maxangle)
end
RandomSingleGate(seed::Int, maxangle::Number) = RandomSingleGate(MersenneTwister(seed), seed, maxangle)

reset!(gate::RandomSingleGate) = Random.seed!(gate.rng, gate.seed)

function data(gate::RandomSingleGate)
	params = rand(gate.rng, 4)
	params[2:end] ./= sum(params[2:end].^2)
	theta = params[1]* gate.maxangle
	return cos(theta/2)*I_DATA-sin(theta/2)im* (params[2]*X_DATA+params[3]*Y_DATA+params[4]*Z_DATA)
end

# Built-in interface still work for customized gate
state = MPSState(2)
state_copy = copy(state)
customizedgate = RandomSingleGate(5, π)
apply!(state, customizedgate, 1)
apply!(state, customizedgate, 2)
reset!(customizedgate)
apply!(state_copy, customizedgate, 1)
apply!(state_copy, customizedgate, 2)
@test toExactState(state) ≈ toExactState(state_copy)


text:
module RangeEnclosures

using Requires, Reexport
@reexport using IntervalArithmetic
const Interval_or_IntervalBox = Union{Interval, IntervalBox}

#=================
Available methods
==================#
const available_solvers = [:IntervalArithmetic,
                           :IntervalOptimisation,
                           :TaylorModels,
                           :SumOfSquares,
                           :BranchandBound]

const optional_solvers = [:AffineArithmetic]

include("intervals.jl")
include("taylormodels.jl")
include("sdp.jl")
include("branchandbound.jl")

#================
Optional methods
=================#
function __init__()
    @require AffineArithmetic = "2e89c364-fad6-56cb-99bd-ebadcd2cf8d2" include("affine.jl")
end

#================
API
=================#
include("enclose.jl")

export enclose

end # module


text:
"""
This functions are julia implementations of the functions in scikit-tda/tadasets.

Reference: https://github.com/scikit-tda/tadasets
"""

abstract type TopologicalDataType end
struct Infty <: TopologicalDataType end
struct SwissRoll <: TopologicalDataType end
struct Sphere <: TopologicalDataType
    d::Int64
end
struct Torus <: TopologicalDataType
    R::Float64
    r::Float64
end

function Sphere()
    return Sphere(2)
end
function Torus()
    return Torus(4.0, 1.0)
end

import Base.rand
"""
Infty

Sample random points on lemniscate.

Reference: https://math.stackexchange.com/a/7989/459895
"""
function rand(S::Infty, n::Int64; noise = 0.01)
    t = 2π * rand(Float64, n)
    x = cos.(t)
    y = sin.(t) .* cos.(t)
    return vcat(x',y') + noise * randn(2, n)
end

"""
Swiss Roll
Reference: http://people.cs.uchicago.edu/~dinoj/manifold/swissroll.html
"""
function rand(S::SwissRoll, n::Int64; noise = 0.01)
    t = 2π * rand(Float64, n)
    y = 2π * rand(Float64, n)
    x = t .* cos.(t)
    z = t .* sin.(t)
    return vcat(x',y',z') + noise * randn(3, n)
end

"""
d-dimensional Sphere

Reference: https://en.wikipedia.org/wiki/N-sphere#Generating_random_points
"""
function rand(S::Sphere, n::Int64; noise = 0.01)
    direction = randn(S.d, n)
    boundary = direction ./ sqrt.(sum(abs2, direction, dims = 1))
    return boundary + noise * randn(S.d, n)
end

"""
Torus
Reference: https://en.wikipedia.org/wiki/Torus#Geometry
"""
function rand(S::Torus, n::Int64; noise = 0.01)
    u = 2π * rand(Float64, n)
    v = 2π * rand(Float64, n)
    x = (S.R .+ S.r * cos.(v)) .* cos.(u)
    y = (S.R .+ S.r * cos.(v)) .* sin.(u)
    z = S.r * sin.(v)
    return vcat(x',y',z') + noise * randn(3, n)
end


text:
import ..Common.MARKERSIZE_SCALE
import ..Common.MARKERSHAPE
import ..Symmetric

# Cartesian Symmetric Simulation
@recipe function f(sim::Simulation{2,Symmetric,C,Tuple{LD}}) where {C,LD<:LatticeDomain{2,Symmetric,Cartesian}}
	lattice = sim.lattice_domain.lattice
	indices = sim.lattice_domain.indices
	imin, imax = extrema(map(ld->ld[1],indices))
	jmin, jmax = extrema(map(ld->ld[2],indices))
	x = [lattice[i,0].x for i ∈ imin:imax]
	y = [lattice[0,j].y for j ∈ jmin:jmax]

	n = permutedims(reshape(sqrt.(sim.ε[:]),length(x),length(y)))

    layout --> (1,3)
    legend --> false

	climr = maximum(real(n))
	climi = max(maximum(abs.(imag(n))),1e-5)
	climf = max(maximum(sim.F),1e-5)

    @series begin
		title --> L"{\rm Real}(n)"
        subplot := 1
        seriestype --> :heatmap
		seriescolor --> :sequential
		aspectratio --> 1
		clims --> (1,climr)
        x, y, real(n)
    end
	@series begin
		title --> L"{\rm Imag}(n)"
		subplot := 2
        seriestype --> :heatmap
		seriescolor --> :diverging
		aspectratio --> 1
		clims --> (-climi,climi)
        x, y, imag(n)
    end
	@series begin
		title --> L"{\rm Pump\ }F"
        subplot := 3
        seriestype --> :heatmap
		seriescolor --> :diverging
		aspectratio --> 1
		subplot := 3
		zcolor --> sim.F[:]
		clims --> (-climf,climf)
        y, x, reshape(sim.F[:],length(x),:)
    end
end

# Generic Simulation
@recipe function f(sim::Simulation{2})
    n = sqrt.(sim.ε[:])

    layout := (1,3)
    legend --> false

    @series begin
		title --> L"{\rm Real}(n)"
        subplot := 1
        seriestype --> :scatter
        markersize --> MARKERSIZE_SCALE/sqrt(length(sim.x))
        markerstrokealpha --> 0
        shape --> MARKERSHAPE
		seriescolor --> :sequential
		aspectratio --> 1
		zcolor --> real(n)
		clims --> (1,3)
        sim.x
    end
	@series begin
		title --> L"{\rm Imag}(n)"
        subplot := 2
        seriestype --> :scatter
        markersize --> MARKERSIZE_SCALE/sqrt(length(sim.x))
        markerstrokealpha --> 0
        shape --> MARKERSHAPE
		seriescolor --> :diverging
		aspectratio --> 1
		zcolor --> imag(n)
		clims --> (-1,1)
        sim.x
    end
	@series begin
		title --> L"{\rm Pump\ }F"
        subplot := 3
        seriestype --> :scatter
        markersize --> MARKERSIZE_SCALE/sqrt(length(sim.x))
        markerstrokealpha --> 0
        shape --> MARKERSHAPE
		seriescolor --> :diverging
		aspectratio --> 1
		subplot := 3
		zcolor --> sim.F[:]
		clims --> (-1,1)
        sim.x
    end
end

# Cartesian Symmetric ScalarField
@recipe function f(sim::Simulation{2,Symmetric,C,Tuple{LD}}, ψ::ScalarField{2}) where {C,LD<:LatticeDomain{2,Symmetric,Cartesian}}
	m = size(ψ,2)
	layout --> (1+m,3)
    legend --> false
	@series begin
		sim
	end
	for μ ∈ 1:m
		ψμ = ψ(μ)
		clim = max(maximum(abs.(real(ψμ))), maximum(abs.(real(ψμ))))
	    @series begin
			title --> L"{\rm Real}(\psi)"
	        subplot := 3μ+1
			sim, ψμ, real
	    end
		@series begin
			title --> L"{\rm Imag}(\psi)"
	        subplot := 3μ+2
	        sim, ψμ, imag
	    end
		@series begin
			title --> L"|\psi|^2"
	        subplot := 3μ+3
	        sim, ψμ, abs2
	    end
	end
end

# Cartesian Symmetric ScalarField by function
@recipe function f(sim::Simulation{2,Symmetric,C,Tuple{LD}}, ψ::ScalarField{2}, by::Function) where {C,LD<:LatticeDomain{2,Symmetric,Cartesian}}

	lattice = sim.lattice_domain.lattice
	indices = sim.lattice_domain.indices
	imin, imax = extrema(map(ld->ld[1],indices))
	jmin, jmax = extrema(map(ld->ld[2],indices))
	x = [lattice[i,0].x for i ∈ imin:imax]
	y = [lattice[0,j].y for j ∈ jmin:jmax]

	legend --> false

	for μ ∈ 1:size(ψ,2)
		ψμ = ψ(μ)
		clim = max(maximum(abs.(by.(ψμ))), maximum(abs.(by.(ψμ))))
	    @series begin
	        seriestype --> :heatmap
			aspectratio --> 1
			if by ∈ (abs, abs2)
				clims --> (0,clim)
				seriescolor --> :sequential
			else
				clims --> (-clim,clim)
				seriescolor --> :diverging
			end
	        x, y, by.(ψμ)
	    end
	end
end

# Generic ScalarField
@recipe function f(ψ::ScalarField{2})
	m = size(ψ,2)

	layout := (m,3)
    legend --> false
	for μ ∈ 1:m
		ψμ = ψ(μ)

		clim = max(maximum(abs.(real(ψμ))), maximum(abs.(real(ψμ))))
	    @series begin
			title --> L"{\rm Real}(\psi)"
	        subplot := 3*(μ-1)+1
	        seriestype --> :scatter
	        markersize --> MARKERSIZE_SCALE/sqrt(length(ψ.positions))
	        markerstrokealpha --> 0
	        shape --> MARKERSHAPE
			seriescolor --> :diverging
			aspectratio --> 1
			zcolor --> real(ψμ)
			clims --> (-clim,clim)
	        ψ.positions
	    end
		@series begin
			title --> L"{\rm Imag}(\psi)"
	        subplot := 3*(μ-1)+2
	        seriestype --> :scatter
	        markersize --> MARKERSIZE_SCALE/sqrt(length(ψ.positions))
	        markerstrokealpha --> 0
	        shape --> MARKERSHAPE
			seriescolor --> :diverging
			aspectratio --> 1
			zcolor --> imag(ψμ)
			clims --> (-clim,clim)
	        ψ.positions
	    end
		@series begin
			title --> L"|\psi|^2"
	        subplot := 3*(μ-1)+3
	        seriestype --> :scatter
	        markersize --> MARKERSIZE_SCALE/sqrt(length(ψ.positions))
	        markerstrokealpha --> 0
	        shape --> MARKERSHAPE
			seriescolor --> :sequential
			aspectratio --> 1
			zcolor --> abs2.(ψμ)
			clims --> (0, maximum(abs2.(ψμ)))
	        ψ.positions
	    end
	end
end

# Electric Field
# @recipe function f(e::ElectricField{1}, by::Function)
    # x = map(a->a.vec[1],e.pos)
    # perm = sortperm(x)
	# x = x[perm]

	# upper_lim = maximum(abs∘by,e.val)
	# if by ∈ (abs,abs2)
		# lower_lim = zero(upper_lim)
	# else
		# lower_lim = -upper_lim
	# end
#
# 	m = size(e,2)
#     layout --> (3,m)
#     legend --> false
# 	ylims --> (lower_lim,upper_lim)
# 	for μ ∈ 1:m
# 	    @series begin
# 			title --> latexstring("{\\rm Field}\\ $μ")
# 	        ylabel --> L"E_x"
# 	        subplot := 0m+μ
# 	        x, by.(e(μ).x[perm])
# 	    end
# 	    @series begin
# 	        ylabel --> L"E_y"
# 	        subplot := 1m+μ
# 	        x, by.(e(μ).y[perm])
# 	    end
# 	    @series begin
# 			xlabel --> L"x"
# 	        ylabel --> L"E_z"
# 	        subplot := 2m+μ
# 	        x, by.(e(μ).z[perm])
# 	    end
# 	end
# end
#
#
# @recipe function f(sim::Simulation{1},e::ElectricField{1}, by::Function)
#     x = map(a->a.vec[1],e.pos)
#     perm = sortperm(x)
#     x = x[perm]
#
# 	legend --> false
#
# 	m = 1+size(e,2)
# 	layout --> (3,m)
#
# 	n = sqrt.(sim.ε[perm])
# 	@series begin
# 		title --> L"{\rm Scattering\ Structure}"
# 		ylabel --> L"{\rm real}(n)"
#         subplot := 0m+1
#         x, real(n)
#     end
#     @series begin
# 		ylabel --> L"{\rm imag}(n)"
#         subplot := 1m+1
#         x, imag(n)
#     end
#     @series begin
# 		xlabel --> L"x"
# 		ylabel --> L"{\rm pump\ F}"
# 		subplot := 2m+1
#         x, sim.F[perm]
#     end
#
# 	upper_lim = maximum(abs∘by,e.val)
# 	if by ∈ (abs,abs2)
# 		lower_lim = zero(upper_lim)
# 	else
# 		lower_lim = -upper_lim
# 	end
# 	ylims --> (lower_lim,upper_lim)
#
# 	for μ ∈ 1:m-1
# 	    @series begin
# 			title --> latexstring("{\\rm Field}\\ $μ")
# 	        ylabel --> L"E_x"
# 	        subplot := 0m+μ+1
# 	        x, by.(e(μ).x[perm])
# 	    end
# 	    @series begin
# 	        ylabel --> L"E_y"
# 	        subplot := 1m+μ+1
# 	        x, by.(e(μ).y[perm])
# 	    end
# 	    @series begin
# 			xlabel --> L"x"
# 	        ylabel --> L"E_z"
# 	        subplot := 2m+μ+1
# 	        x, by.(e(μ).z[perm])
# 	    end
#
# 	end
# end


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

using FEMDynamics
using JuliaFEM
using FEMBase.Test

X = Dict(1 => [0.0, 0.0],
         2 => [2.0, 0.0],
         3 => [2.0, 2.0],
         4 => [0.0, 2.0])

element = Element(Quad4, [1, 2, 3, 4])
update!(element, "geometry", X)
update!(element, "youngs modulus", 288.0)
update!(element, "poissons ratio", 1/3)
update!(element, "density", 36.0)
update!(element, "displacement load 2", 36.0*9.81)
problem = Problem(Elasticity, "2x2 element", 2)
problem.properties.formulation = :plane_stress
add_elements!(problem, [element])
analysis = Analysis(LinearImplicit, "dropping element")
add_problems!(analysis, [problem])
run!(analysis)
midpnt_u = element("displacement", (0.0, 0.0), Inf)
midpnt_v = element("velocity", (0.0, 0.0), Inf)
@test isapprox(midpnt_u, [0.0, 1/2*9.81])
@test isapprox(midpnt_v, sqrt.(2.0*9.81*midpnt_u))


text:
module Jules

# global variables
Fs = 5000; # sampling frequency
T = 1/Fs;
# functions

function gentime(trial, sampling_freq=Fs)
    T = 1/sampling_freq
    L = size(trial, 1)
    return t = collect(0:L-1)*T
end

function genrange(start, finish)
    return convert(Integer, floor(start)):convert(Integer, floor(finish))
end

function genrangeseconds(start_sec, finish_sec, sampling_freq=Fs)
    return genrange(start_sec*sampling_freq, finish_sec*sampling_freq)
end

function fftplot(X; sampling_freq=Fs)
    T = 1/sampling_freq
    L = size(X, 1)
    t = collect(0:L-1)*T

    Y = fft(X)
    P2 = abs(Y/L)
    P1 = P2[1:div(L,2)+1]
    P1[2:end - 1] = 2*P1[2:end-1]

    freq = sampling_freq*collect(0:(L/2))/L

    subplot(2,1,1)
    plot(t,X)
    title("Time domain")
    xlabel("time (s)")
    ylabel("Amplitude")
    subplots_adjust(hspace=.5)
    subplot(2,1,2)
    plot(freq,P1)
    title("Single-Sided Amplitude Spectrum")
    xlabel("f (Hz)")
    ylabel("|P1(f)|")
    return (P1,freq)
end

function subdc(trial)
    trial - mean(trial)
end

function discreteder(trial, sampling_freq=Fs)
    T_s = 1/Fs
    b = zeros(size(trial))

    for i in 1:length(trial)-1
        b[i] = (trial[i+1] - trial[i])/T_s
    end

    return b
end

function getstats(trial_column)
    stddev = std(trial_column)
    mn = mean(trial_column)
    println("min: ", minimum(trial_column))
    println("max: ", maximum(trial_column))
    println("average: ", mn)
    println("standard dev: ", stddev)
    ran = -3.5*stddev+mn : 0.01 : 3.5*stddev+mn
    #plot(ran, gaussian(ran, mn, stddev))
    #hold
    #plot(stddev+mn,0,"o")
    #plot(-stddev+mn,0,"o")
end

function gaussian(x, avg, std_dev)
    return 1/(std_dev * sqrt(2*pi))*exp(-.5*((x-avg)/std_dev).^2)
end

function smooth(arr, num_point=5)
    if num_point % 2 == 0
        println("num_point must be odd")
        return
    end

    b = zeros(size(arr))
    span = div(num_point, 2)
    start, finish = 1, length(arr)

    for i in 1:length(arr)
        n = 0

        inspan = i-span < start ? i-start : span
        inspan = i+span > finish ? finish-i : inspan

        #println("start ", start, " finish ", finish, " i ", i, " inspan ", inspan)
        for j in max(start, i-inspan):min(finish, i+inspan)
            #println("j", j, " arr ", arr[j])
            b[i] += arr[j]
            n+=1
        end

        b[i]/=n
    end
    return b
end

#LOW-PASS WINDOWED-SINC FILTER

function windowsinc(input_arr, fc, kernel_length::Integer)
    if kernel_length % 2 == 0
        println("num_point must be odd")
        #only applies to 1 indexed languages
        return
    end
    #Set the cutoff frequency (between 0 and 0.5 relative to sampling freq)
    #kernel_length = 4/BW

    input = vcat(zeros(kernel_length), input_arr)
    kernel = zeros(kernel_length) #holds the filter kernel
    output = zeros(length(input))

    #Calculate the low-pass filter kernel

    kernel_length_half = div(kernel_length+1, 2)
    for i in 1:kernel_length
            if (i-kernel_length_half) == 0
                    kernel[i] = 2*pi*fc
            else
                    kernel[i] = sin(2*pi*fc * (i-kernel_length_half)) / (i-kernel_length_half)
                    kernel[i] = kernel[i] * (0.54 - 0.46*cos(2*pi*i/kernel_length) )
            end
    end

    #Normalize the low-pass filter kernel for unity gain at DC
    kernel_sum = sum(kernel)
    kernel = kernel ./ kernel_sum

    #Convolve the input signal & filter kernel
    for j in kernel_length+1:length(input)
            for i in 1:kernel_length
                    output[j] += input[j-i] * kernel[i]
            end
    end
    return output[kernel_length+1:end]
end

function acceltogs(x)
    return x*2.0/32768
end

function windowsinchp(input_arr, fc, kernel_length::Integer)
    if kernel_length % 2 == 0
        println("num_point must be odd")
        #only applies to 1 indexed languages
        return
    end
    #Set the cutoff frequency (between 0 and 0.5 relative to sampling freq)
    #kernel_length = 4/BW

    input = vcat(zeros(kernel_length), input_arr)
    kernel = zeros(kernel_length) #holds the filter kernel
    output = zeros(length(input))

    #Calculate the low-pass filter kernel

    kernel_length_half = div(kernel_length+1, 2)
    for i in 1:kernel_length
            if (i-kernel_length_half) == 0
                    kernel[i] = 2*pi*fc
            else
                    kernel[i] = sin(2*pi*fc * (i-kernel_length_half)) / (i-kernel_length_half)
                    kernel[i] = kernel[i] * (0.54 - 0.46*cos(2*pi*i/kernel_length) )
            end
    end

    #Normalize the low-pass filter kernel for unity gain at DC
    kernel_sum = sum(kernel)
    kernel = kernel./kernel_sum

    kernel = kernel .* -1
    kernel[kernel_length_half] += 1

    #Convolve the input signal & filter kernel
    for j in kernel_length+1:length(input)
            for i in 1:kernel_length
                    output[j] += input[j-i] * kernel[i]
            end
    end
    return output[kernel_length+1:end]
end

function findpeaks(x, tolerance=0.0001)
    peaks = []
    counter::Integer = 1
    for i in 2:length(x)-1
        if abs(x[i] - x[i+1]) > tolerance && abs(x[i] - x[i-1]) > tolerance
            append!(peaks, [(x[i], i)])
            counter += 1
        end
    end
    return peaks
end

function decimate(x, num::Integer)
    return x[1:num:length(x)]
end

function genfreqindex(finish, sampling_freq = Fs)
    return 1:finish/sampling_freq
end

function getenergy(signal)
    return sum(signal.^2)
end

function nernst(extra, intra, z, T=298)
    R = 8.314472
    F = 9.64853399 * 10.0^4
    E = (R*T)/(z*F)log(extra/intra)
end

function idft(X)
    a = conj(X)
    b = fft(a)
    return real(b)
end

function RK4(u, du, t::Float64, dt::Float64)
    k1 = du(t, u)
    k2 = du(t + dt/2, u + dt/2*k1)
    k3 = du(t + dt/2, u + dt/2*k2)
    k4 = du(t + dt, u + dt*k3)
    dU = k1 + 2*k2 + 2*k2 + k4
    return U = u + dU*dt/6.0
end

function Int4(u::Vector{Float64}, du, t, dt)
    A = 18.0/11.0
    B = -9.0/11.0
    C = 2.0/11.0
    D = 6.0/11.0

    # U = A*u[n-1] + B*[n-2] + C*[n-3] + D*du*dt
    return U = A*u[end-1] + B*u[end-2] + C*u[end-3] + D*du(t, u[end])*dt
end

function adams_bashforth_moulton(u::Float64, prev_du::Vector{Float64}, du, t::Float64, dt::Float64)
        f_k = du(t, u)
        predictor =  u + dt/24 * (-9*prev_du[end-2] + 37*prev_du[end-1] - 59*prev_du[end] + 55*f_k)
        corrector = u + dt/24 * (prev_du[end-1] - 4*prev_du[end] + 19*f_k + 9*du(t+dt, predictor))
        return corrector
end

function integ(method)
end

function constrain(a::Float64, lower_bound::Float64, upper_bound::Float64)
    if a > upper_bound
        return upper_bound
    end

    if a < lower_bound
        return lower_bound
    end

    return a
end

function interp(x_basis::Vector{Float64}, y_basis::Vector{Float64}, t::Float64)
    #print("time $t\n")
    next_index = findfirst(x -> x > t, x_basis)
    if t >= x_basis[end]
        next_index = length(x_basis)
    end
    prev_index = next_index -1
    nearest_time = x_basis[prev_index]
    #=
    print("nearest time $nearest_time\n")
    =#
    avg_constant = t - x_basis[prev_index]

    #=
    print("avg_constant ")
    print(avg_constant)
    print("\n")

    print("1 -avg_constant ")
    print(1- avg_constant)
    print("\n")
    print("1 -avg_constant ")
    print(1- avg_constant)
    print("\n")
    =#
    out = (1 - avg_constant)*y_basis[prev_index] + avg_constant*y_basis[next_index]

    #=
    print("out ")
    print(out)
    print("\n\n")
    =#

    return out
end

end


text:
using PGFPlotsX, Distributions, Cubature, LaTeXStrings, Interpolations

""" p(y|x) = N(f(x),σ)"""
p_y_x(y,x; σ=σ, slope=slope, context...) = pdf(Normal(f(x; slope=slope),σ),y)

""" p(x) = N(0,θ)"""
p_x(x; θ=θ, context...) = pdf(Normal(0,θ),x)

""" p(x,y) = p(y|x)p(x)"""
p_xy(x,y; context...) = p_y_x(y,x; context...)*p_x(x; context...)

""" p(y) = ∫ p(x,y) dx"""
p_y(y; θ=θ, context...) = hquadrature(x->p_xy(x,y; θ=θ, context...), -10*θ, 10*θ; reltol=ε)[1]

""" p(x|y) = p(x,y)/p(y)"""
p_x_y(x,y; context...) = p_xy(x,y; context...)/p_y(y; context...)

""" p(ŷ|y) = N(ŷ,σ)p(ŷ)/p(y)"""
p_ŷ_y(ŷ, y; context...) = p_x_y(f⁻(ŷ; context...),y; context...)*abs(f⁻′(ŷ; context...))

# How much does y reveal about x?
# Q: How much better than the prior do we know where x is after observing y?
# A: KLD from posterior to prior
# = ∫ p(y) ∫ p(x|y)log(p(x|y)/p(x)) dx dy
# = ∫∫ p(x,y) log(p(x|y)/p(x)) dxdy
# = ∫∫ p(x,y) log(p(x,y)/(p(x)p(y))) dxdy
# = MI(X,Y)

function observed_information(x,y; context...)
    P_x = p_x(x; context...)
    P_y = p_y(y; context...)
    P_y_x = p_y_x(y,x; context...)
    P_xy = P_y_x*P_x
    return P_xy*(log2(P_y_x)-log2(P_y))
end

MI(; context...) = hcubature(((x,y),)->observed_information(x,y; context...), [-10θ, -1-10σ], [10θ, 1+10σ]; reltol=ε)[1]

ε = 1e-9
σ = 0.1
θ = 1.0

f(x; slope=2.0, context...) = tanh(slope*x)
f⁻(ŷ; slope=2.0, context...) = atanh(ŷ)/slope
f⁻′(ŷ; slope=2.0, context...) = 1.0/(slope*(1.0-ŷ^2))

function test_mis()
    mi = zeros(100)
    slopes = 10 .^ LinRange(-2, 2, length(mi))
    for (i,slope) ∈ enumerate(slopes)
        mi[i] = MI(θ=θ, σ=σ, slope=slope)
        # println("Got MI = $(mi[i]) for slope=$(slope).")
    end
    return (slopes=slopes, mi=mi)
end

(slopes,mi) = test_mis()
slope_opt = slopes[argmax(mi)]

mi_fun = LinearInterpolation(slopes, mi)

push!(PGFPlotsX.CUSTOM_PREAMBLE, raw"\usetikzlibrary{backgrounds}")

xlims_plt = -3, 3
ylims_plt = -1.1, 1.1
x_plt = LinRange(xlims_plt...,200)
y_plt = LinRange(-0.995, 0.995, 200)
y_plt_prior = LinRange(-1.25, 1.25, 200)
y_true = [-0.90, 0.0, 0.90]
special_slopes = [10^-1, slope_opt, 10^1]

fig = begin
    fig = @pgf TikzPicture()

    ax = @pgf Axis(
        {
            axis_background_style={fill="none"},
            width=raw"\textwidth-1cm",
            height=raw"4cm",
            xlabel = raw"$\textrm{gain }f'(0)$",
            ylabel = raw"$\mathcal{I}(X;Y)$",
            xmode = "log",
            clip="false",
            xmin="10^-2",
            xmax="10^2",
            no_markers,
            name="curve",
        },
        # draw curve
        PlotInc({line_width="2pt", smooth}, Table(slopes, mi)),
        # draw asymptotes
        Plot({color="black", line_width="1pt", dashed}, Table([slopes[1],slopes[end]], [0,0])),
        Plot({color="black", line_width="1pt", dashed}, Table([slopes[1],slopes[end]], [1,1])),            # draw vertical lines
        # # draw vertical lines
        # (VLine({ color="gray", line_width="1pt"}, slope) for slope ∈ special_slopes)...,
        ("\\node (c-0-$(i)) at ($(slope),$(mi_fun(slope))) {\$\\pgfuseplotmark{*}\$}; " for (i,slope) ∈ enumerate(special_slopes))...
    )

    push!(fig,ax)

    for (i,slope) ∈ enumerate(special_slopes)
        context = Dict(:slope=>slope, :σ=>σ, :θ=>θ)

        # draw function plot
        ax_f = @pgf Axis(
            {
                axis_background_style={fill="black!5!white", fill_opacity=0.5},
                height="4cm",
                width=raw"0.4*\textwidth",
                xtick=[-1,0,1],
                ytick=[-1,0,1],
                xlabel=raw"$X$",
                ylabel=raw"$Y=f(X)+\eta$",
                xmin=xlims_plt[1],
                xmax=xlims_plt[2],
                ymin=ylims_plt[1],
                ymax=ylims_plt[2],
            },
            # draw function
            Plot(
                {
                    color="black",
                    line_width="2pt",
                },
                Table(x_plt, f.(x_plt; context...))
            ),
            # draw horizontal lines
            (HLine({ dashed, index_of_colormap=j, line_width="1pt", }, yy ) for (j,yy) ∈ enumerate(y_true))...,
            "\\coordinate (c-$(i)-1) at (axis description cs: 0,0);",
            "\\coordinate (c-$(i)-2) at (axis description cs: 1,0);",
            # add slope label
            "\\path (-1,-$(slope)) -- node[midway, sloped, anchor=south] {s=$(slope)} (1,$(slope));",
        )

        # draw x distributions
        ax_x = @pgf Axis(
            {
                at=[raw"($(curve.north west)+(0cm,3.5cm)$)",raw"($(curve.north)+(0cm,7cm)$)",raw"($(curve.north east)+(0cm,3.5cm)$)"][i],
                anchor=["south west", "south", "south east"][i],
                axis_background_style={fill="none"},
                height="2.5cm",
                width=raw"0.4*\textwidth",
                hide_axis,
                xlabel=raw"$p(X)^*, p(X|y_i)^*$",
                xmin=xlims_plt[1],
                xmax=xlims_plt[2],
                ymin=0,
            },
            Plot(
                {
                    draw = "none",
                    fill="gray",
                    fill_opacity=0.5,
                },
                Table(x_plt,p_x.(x_plt; context...)./maximum(p_x.(x_plt; context...)))
            ),
            # draw posteriors
            (Plot(
                {
                    draw = "none",
                    index_of_colormap=j,
                    fill=".",
                    fill_opacity=0.5
                },
                Table([xlims_plt[1];x_plt;xlims_plt[2]], [0;p_x_y.(x_plt,yⱼ; context...)./maximum(p_x_y.(x_plt,yⱼ; context...));0])
            ) for (j,yⱼ) ∈ enumerate(y_true))...,
        )

        # draw y distributions
        ax_y = @pgf Axis(
            {
                axis_background_style={fill="none"},
                height="4cm",
                width="2.5cm",
                hide_axis,
                xlabel=raw"$p(Y)^*,p(f(X)|y_i)^*$",
                xmin=0,
                ymin=ylims_plt[1],
                ymax=ylims_plt[2],
            },
            # draw prior
            Plot(
                {
                    draw = "none",
                    fill="gray",
                    fill_opacity=0.5,
                },
                Table([0.0; p_y.(y_plt; context...)./maximum(p_y.(y_plt; context...)); 0.0], [-1.0; y_plt; 1.0])
            ),
            # # draw posteriors
            # (Plot(
            #     {
            #         draw = "none",
            #         index_of_colormap=j,
            #         fill=".",
            #         fill_opacity=0.5
            #     },
            #     Table([0.0; p_ŷ_y.(y_plt,yⱼ; context...)./maximum(p_ŷ_y.(y_plt,yⱼ; context...)); 0], [-1.01; y_plt; 1.01])
            # ) for (j,yⱼ) ∈ enumerate(y_true))...,
        )

        gp = @pgf GroupPlot(
            {
                group_style = {
                    group_name="group$(i)",
                    group_size="2 by 2",
                    horizontal_sep="0pt",
                    vertical_sep="0pt",
                    xticklabels_at="edge bottom",
                    yticklabels_at="edge left"
                },
            },
            ax_x, {"group/empty plot",width="2.5cm",height="2.5cm", axis_background_style={fill="none"}}, ax_f, ax_y
        )
        #     plt[xidx][:yaxis][:flip]=true
        #     plt[fidx][:xaxis][:mirror]=true
        #     plt[yidx][:yaxis][:mirror]=true
        #     plt[yidx][:xaxis][:mirror]=true
        #     plt[eidx][:yaxis][:mirror]=true
        push!(fig, gp)
    end

    push!(fig,raw"\begin{scope}[on background layer]",("\\path[shade,bottom color=black!50!white,top color=black!10!white, fill opacity=0.2] (c-0-$(i).center) -- (c-$(i)-1) -- (c-$(i)-2) -- cycle;" for i ∈ 1:3)...,raw"\end{scope}")

    return fig
end

max_stars_repo_name: jleugeri/phd

max_issues_repo_name: jleugeri/phd

max_forks_repo_name: jleugeri/phd


text:
# TODO: use https
const IMAGE_SERVER = "http://image.eveonline.com"

img_get(endpoint, id, size, filetype = endpoint == "Character" ? "jpg" : "png") = eget("$IMAGE_SERVER/$endpoint/$(id)_$size.$filetype")


text:
module equations

using StatsModels, MixedModels, CategoricalArrays, CSV, StatsBase

function make_ran_matrix(x1::AbstractVector,x2::AbstractVector)
        isa(x1, CategoricalArray) ||
                       throw(ArgumentError("ran() only works with CategoricalArrays (got $(typeof(2)))"))
        isa(x2, CategoricalArray) ||
                       throw(ArgumentError("ran() only works with CategoricalArrays (got $(typeof(2)))"))

        u = unique(x2);
        filter!(x->x≠0,u)
        Z = Matrix{Bool}(undef, length(x1), length(u))
        for i in eachindex(u)
        	@. Z[:, i] = x1 .== u[i]
        end
           return unique(x1),Z
       end

ranMat(arg1,arg2,data1,data2) = make_ran_matrix(data1[!,Symbol(arg1)],data2[!,Symbol(arg2)])

function mme(f, userHints, userData, userPedData, blocks; paths2geno)
        terms4StatsModels = String.(split(repr(f.rhs), ('+')))
        terms4StatsModels = replace.(terms4StatsModels, ":" => "")
        terms4StatsModels = [filter(x -> !isspace(x), trm) for trm in terms4StatsModels]

        yVec = StatsModels.modelmatrix(f.lhs, userData)

        FE = Dict{Any,Any}()

        RE = Dict{Any,Any}()

	ME = Dict{Any,Any}()
	regionSizes = Array{Int64,1}(undef,0)

	#column id within pedigree
	idRE = []

        for i in 1:length(f.rhs)
		if (f.rhs[i] isa FunctionTerm) && (String(nameof(f.rhs[i].forig)) == "PR")
			println("$i has type BayesPR Type")
			println("terms4StatsModels[i]: $(terms4StatsModels[i])")
			arg1 = repr((f.rhs[i].args_parsed)[1])
			arg2 = parse(Int64,repr((f.rhs[i].args_parsed)[2]))
			path = paths2geno[Symbol(arg1)]
			thisM = CSV.read(path,CSV.Tables.matrix)
			#centering
			thisM .-= mean(thisM,dims=1)
			println("size of $arg1 data: $(size(thisM))")
			println("region size for $arg1: $arg2")
			ME[arg1] = thisM
                        thisM = 0 #I can directly merge to dict above
			push!(regionSizes, arg2)
                elseif (f.rhs[i] isa FunctionTerm) && (String(nameof(f.rhs[i].forig)) == "ran")
                        println("$i has type ran Type")
                        sym1 = repr((f.rhs[i].args_parsed)[1]) #now it is Symbol
                        sym2 = repr((f.rhs[i].args_parsed)[2]) #now it is from string
#                       arg2 = eval(Meta.parse(arg2)) #now it is from string to data. Later will be path
                        println("sym1: $sym1 sym2: $sym2")

			IDs,thisZ = ranMat(sym1, sym2, userData, userPedData)
			RE[sym1] = thisZ
			thisZ = 0
			push!(idRE,IDs)
                elseif (f.rhs[i] isa FunctionTerm) && (String(nameof(f.rhs[i].forig)) == "|")
                        println("$i has type | Type")
                        my_sch = schema(userData, userHints) #work on userData and userHints
                        my_ApplySch = apply_schema(terms(f.rhs[i]), my_sch, MixedModels.MixedModel)
			#####NO IDs for this effect!!! Will be added later!!!!#####################################################
                       	thisZ = modelcols(my_ApplySch, userData)
			RE[terms4StatsModels[i]] = thisZ
			thisZ = 0
			push!(RE,modelcols(my_ApplySch, userData))

                else
                println("$i has type $(typeof(f.rhs[i]))")
		thisX = StatsModels.modelmatrix(f.rhs[i], userData,hints= userHints)
		FE[terms4StatsModels[i]] = thisX
		thisX = 0
                end
        end

	#merging on the earliest, and then delete merged ones from the list
	#do not consider if there is reordering of variables
	#blocks must enter in the model in order!
	mergedOnes = []
	delThese = []
	newNamesFE = []
	for i in blocks
		blockThese = findall(x->x in i, namesFE)
		mergeTo = minimum(blockThese)
		FE[mergeTo] = hcat(FE[blockThese]...)
		delThese = vcat(delThese,blockThese[blockThese .!== minimum(blockThese)])
		push!(mergedOnes, blockThese)
	end

	delete!(FE, sort(delThese))

        return idRE, vec(yVec), FE, RE, ME, regionSizes
        end

end

max_stars_repo_name: animBreedTools/NextGP

max_issues_repo_name: animBreedTools/NextGP

max_forks_repo_name: animBreedTools/NextGP


text:
# module deals with querying eve-central's market data
#include("SDEUtils.jl")
module EveCentral
include("SDE.jl")
using .SDE
using DataFrames
using Requests
const baseURI = "http://api.eve-central.com/api"

function marketstat(itemNames :: Array{ASCIIString,1}, systemName :: AbstractString)
  const uri = "http://api.eve-central.com/api/marketstat/json"
  types = map(itemID, itemNames)
  types = map(string, types)

  q = Dict()
  q["typeid"] = join(types, ",")
  q["usesystem"] = string(solarSystemID(systemName))
  get(uri, query = q)
end
function marketstat(itemName :: AbstractString, systemName :: AbstractString)
  marketstat([itemName], systemName)
end
#historical data functions
function historicalData(item :: AbstractString, system :: AbstractString)
  item_id = itemID(item)
  system_id = solarSystemID(system)
  uri = "http://api.eve-central.com/api/history/for/type/$item_id/system/$system_id/bid/1"
  uri
end

function calcMargin(mktStat)
  1.0 - (mktStat["buy"]["max"] / mktStat["sell"]["min"])
end
function calcMargin(source, dest)
  1.0 - (source["buy"]["max"] / dest["sell"]["min"])
end
function calcVolumeMoved(mktStat)
  min(mktStat["buy"]["volume"], mktStat["sell"]["volume"])
end
function calcAvgProfit(mktStat, expectedVol :: AbstractFloat = 0.05)
  vol = calcVolumeMoved(mktStat)
  vol = expectedVol * vol
  buy = mktStat["buy"]["max"]
  sell = mktStat["sell"]["min"]
  profit = sell - buy
  profit * vol
end

export marketstat, calcMargin, calcAvgProfit, calcVolumeMoved

end


text:
function blkdiag(blocks::Array{Array{Pol{T}, 2}, 1}) where {T<:Number}
    m = sum(size.(blocks, 1))
    n = sum(size.(blocks, 2))
    J = zeros(Pol{T}, m, n)
    offset1 = 0
    offset2 = 0
    for blk in blocks
        J[offset1+1:offset1+size(blk, 1), offset2+1:offset2+size(blk, 2)] = blk
        offset1 += size(blk, 1)
        offset2 += size(blk, 2)
    end
    return J
end

"""
creates the matrix block of size n×(n+1) corresponding to column index n
"""
function colBlock(n::Int)
    var = Pol([0, 1])
    J = zeros(Pol{Int}, n, n+1)
    for ii in 1:n
        J[ii, ii] = var
        J[ii, ii+1] = -1
    end
    return J
end

colBlock(arr::AbstractArray{Int, 1}) = blkdiag(colBlock.(arr))

"""
creates the matrix block of size (n+1)×n corresponding to row index n
"""
function rowBlock(n::Int)
    var = Pol([0, 1])
    J = zeros(Pol{Int}, n+1, n)
    @inbounds for ii in 1:n
        J[ii, ii] = -1
        J[ii+1, ii] = var
    end
    return J
end

rowBlock(arr::AbstractArray{Int, 1}) = blkdiag(rowBlock.(arr))

"""
creates the block matrix corresponding to infinite eigenvalue of multiplicity n
"""
function infBlock(n::Int)
    var = Pol([0, 1])
    J = zeros(Pol{Int}, n, n)
    @inbounds for ii in 1:n-1
        J[ii, ii] = -1
        J[ii+1, ii] = var
    end
    J[n, n] = -1
    return J
end

infBlock(arr::AbstractArray{Int, 1}) = blkdiag(infBlock.(arr))

"""
creates the block matrix of size n×n corresponding to eigenvalue λ
"""
function eigBlock(λ::T, n::Int) where {T<:Number}
    var = Pol([0, 1])
    J = zeros(Pol{T}, n, n)
    @inbounds for ii in 1:n-1
        J[ii, ii] = var - λ
        J[ii+1, ii] = -1
    end
    J[n, n] = var - λ
    return J
end

eigBlock(λ::Number, arr::AbstractArray{Int, 1}) = blkdiag(eigBlock.(λ, arr))


text:
__precompile__()
module RaggedData

using ArgCheck: @argcheck
using AutoHashEquals: @auto_hash_equals
using DataStructures: OrderedDict
using DocStringExtensions: SIGNATURES
using Lazy: @forward
using Parameters: @unpack

import Base:
    count, sizehint!, indices, length, size, getindex, IndexStyle

"""
    RaggedCount(index, count)

An element in a `RaggedCounter`. `index` allows ordering of the elements as they
are encountered, `count` keeps track of the number of elements so far in the
data for a particular key.
"""
@auto_hash_equals mutable struct RaggedCount{S <: Integer}
    index::S
    count::S
end

"""
    RaggedCounter(T, S)

Count keys of type `T` with integers of type `S`. New keys can be added with
`(rc::RaggedCounter(key)`, which returns `key`. The order of new keys is
preserved.

`count(::RaggedCounter)` returns the total number of keys counter, while
`length(::RaggedCounter)` the number of unique keys. Use `sizehint!` for a
slight speed increase if you can estimate the number of elements in advance.

## Recommended usage

Use in the *first pass* of ingesting the data to count the number each key
occurs. Then generate a collator and indexer using `collate_index_keys` for the
*second pass* and then indexing.
"""
@auto_hash_equals struct RaggedCounter{T,S}
    dict::Dict{T,RaggedCount{S}}
end

RaggedCounter(T,S::Type{<:Integer}) =
    RaggedCounter{T,S}(Dict{T,RaggedCount{S}}())

function (rc::RaggedCounter{T,S})(key) where {T,S}
    @unpack dict = rc
    if haskey(dict, key)
        dict[key].count += one(S)
    else
        dict[key] = RaggedCount(S(length(dict)) + one(S), one(S))
    end
    key
end

count(rc::RaggedCounter) = sum(v.count for v in values(rc.dict))

@forward RaggedCounter.dict length, sizehint!

"""
    RaggedCollate

Use in the *second pass* of reading the data to return increasing indices
(`nextindex!`) that form an ordered partition of `1:count` after the whole data
is traversed.
"""
@auto_hash_equals struct RaggedCollate{T,S}
    dict::Dict{T,S}
    total::S
end

next_index!(coll::RaggedCollate{T,S}, key) where {T,S} = coll.dict[key] += one(S)

"""
    RaggedIndex(count_cumsum)

Create an index that maps key indices (with given counts) to indices of a flat
vector. See `collate_index_keys` for efficient creation from `RaggedCounter`s.

Example:
```jldoctest
julia> ri = RaggedIndex(cumsum(1:3))
RaggedData.RaggedIndex{Int64}([1, 3, 6])

julia> ri[1]
1:1

julia> ri[2]
2:3

julia> ri[3]
4:6
```
"""
@auto_hash_equals struct RaggedIndex{S} <: AbstractVector{UnitRange{S}}
    cumsum::Vector{S}
end

@forward RaggedIndex.cumsum length, size

count(ri::RaggedIndex) = ri.cumsum[end]

function getindex(ri::RaggedIndex{S}, i) where S
    UnitRange{S}(one(S) + (i == 1 ? zero(S) : ri.cumsum[i-1]), ri.cumsum[i])
end

"""
    sub_ix, sub_I = _subset(ix::RaggedIndex, I)

Calculate a `RangeIndex` and corresponding set of `UnitRange` indices such that
indexing the original is possible for the subset, ie

```julia
v[I][sub_ix[j]] == v[sub_I[j]]
```

∀ j ∈ indices(v).

Useful for subsets of `RaggedColumn` and `RaggedColumns`.
"""
function _subset(ix::RaggedIndex{S}, I) where {S}
    partial_sum = zero(S)
    sub_cumsum = similar(Array{S}, indices(I))
    sub_I = Vector{S}()
    for (j,i) in enumerate(I)
        k = ix[i]
        append!(sub_I, k)
        partial_sum += length(k)
        sub_cumsum[j] = partial_sum
    end
    RaggedIndex(sub_cumsum), sub_I
end

"""
    collate, index, keys = collate_index_keys(ragged_counter; collect_keys = false)

Return a collator `collate`, an index `index`, and possibly `keys` (may be an
empty vector if `!collect_keys`.
"""
function collate_index_keys(rc::RaggedCounter{T,S},
                            collect_keys = false) where {T,S}
    kv = collect(rc.dict)
    sort!(kv, by = x -> x.second.index)
    partial_sum = zero(S)
    dict = Dict{T, S}()
    sizehint!(dict) = length(kv)
    index = similar(kv, S)
    keys = ifelse(collect_keys, similar(kv, T), Vector{T}(0))
    @inbounds for i in indices(kv, 1) # dict, cumsum and keys in one sweep
        key = kv[i].first
        dict[key] = partial_sum
        partial_sum += kv[i].second.count
        index[i] = partial_sum
        if collect_keys
            keys[i] = key
        end
    end
    RaggedCollate(dict, partial_sum), RaggedIndex(index), keys
end

_vector_view_type(S,T) = SubArray{T,1,Array{T,1},Tuple{UnitRange{S}},true}

_vector_view_type(::RaggedIndex{S}, ::AbstractVector{T}) where {S, T} =
    _vector_view_type(S, T)

"""
    RaggedColumns(count_cumsum, columns)

A tuple of columns (vectors) indexes by a `RaggedIndex`.

`getindex(::RaggedColumns, ::Int)` returns a tuple of vectors for that index.

`getindex(::RaggedColumns, ::Any)` subsets the indexed vector.

Use `RaggedColumn(::RaggedColumns, i::Int)` or `RaggedColumns(::RaggedColumns, I)`
to obtain single or multiple columns.
"""
@auto_hash_equals struct RaggedColumns{Tix <: RaggedIndex,
                                       Tcolumns <: Tuple, S} <:
                                           AbstractVector{S}
    ix::Tix
    columns::Tcolumns
    function RaggedColumns(ix::Tix,
                           columns::Tcolumns) where {Tix <: RaggedIndex,
                                                     Tcolumns <:
                                                     Tuple{Vararg{AbstractVector}}}
        S = Tuple{map(c -> _vector_view_type(ix, c), columns)...}
        new{Tix, Tcolumns, S}(ix, columns)
    end
end

RaggedColumns(rc::RaggedColumns, I) = RaggedColumns(rc.ix, rc.columns[I])

RaggedColumns(rc::RaggedColumns, i::Int) =
    throw(ArgumentError("Use RaggedColumn(::RaggedColumns, ::Int) for a single column."))

@forward RaggedColumns.ix count, length

"""
    RaggedColumn(count_cumsum, column)

A single column indexed by a `RaggedIndex`.

`getindex(::RaggedColumn, ::Int)` returns a vector for that index.

`getindex(::RaggedColumn, ::Any)` subsets the indexed vector.
"""
@auto_hash_equals struct RaggedColumn{Tix <: RaggedIndex,
                                      Tcolumn <: AbstractVector, S} <:
                                          AbstractVector{S}
    ix::Tix
    column::Tcolumn
    function RaggedColumn(ix::Tix, column::Tcolumn) where {Tix <: RaggedIndex,
                                                           Tcolumn <: AbstractVector}
        S = _vector_view_type(ix, column)
        new{Tix, Tcolumn, S}(ix, column)
    end
end

RaggedColumn(rc::RaggedColumns, i::Int) = RaggedColumn(rc.ix, rc.columns[i])

@forward RaggedColumn.ix count, length

size(A::Union{RaggedColumn,RaggedColumns}) = (length(A.ix),)

IndexStyle(::Union{RaggedColumn,RaggedColumns}) = Base.IndexLinear()

function getindex(A::RaggedColumns, i::Int)
    j = A.ix[i]
    map(v -> @view(v[j]), A.columns)
end

function getindex(A::RaggedColumn, i::Int)
    j = A.ix[i]
    @view(A.column[j])
end

function getindex(A::RaggedColumns, I)
    sub_ix, sub_I = _subset(A.ix, to_indices(A.ix, (I,))...)
    RaggedColumns(sub_ix, map(v -> v[sub_I], A.columns))
end

function getindex(A::RaggedColumn, I)
    sub_ix, sub_I = _subset(A.ix, to_indices(A.ix, (I,))...)
    RaggedColumn(sub_ix, A.column[sub_I])
end

# permutation interface experiment

"""
    $SIGNATURES

Calculate counts distinct elements of `itr` as an `OrderedDict`. Counting is
performed using type `S`. See [`contiguous_ranges`](@ref) and
[`contiguous_invperm`](@ref).
"""
function ordered_counts(itr, S::Type{<: Integer} = Int)
    c = OrderedDict{eltype(itr), S}()
    for x in itr
        c[x] = get(c, x, zero(S)) + one(S)
    end
    c
end

"""
    $SIGNATURES

Return the contiguous ranges that [`contiguous_invperm`](@ref) would map data
with `counts` to.
"""
function contiguous_ranges(counts::OrderedDict{T, S}) where {T, S <: Integer}
    N = zero(S)
    [(z = N; N += c; (z+one(S)):N) for c in values(counts)]
end

"""
    $SIGNATURES

Similar to [`contiguous_invperm!`](@ref), but writes the permutation to a
preallocated vector, its first argument.
"""
function contiguous_invperm!(p::AbstractVector{S},
                             x::AbstractVector{T},
                             counts::OrderedDict{T,S}) where {T, S <: Integer}
    I = indices(x, 1)
    @argcheck I == indices(p, 1) "Incompatible indices."
    @argcheck sum(values(counts)) == length(I) "Length of x does not match counts."
    N = zero(S)
    c = OrderedDict((n = N; N += v; k => n) for (k, v) in counts)
    for i in I
        p[i] = (c[x[i]] += one(S))
    end
    p
end

"""
    $SIGNATURES

When `x` is a vector, and `counts[y]` is the number of elements `y .== x`,
return an inverse permutation vector `p` that groups the same values of `x`
together while keeping the order. That is to say, if

```julia
rs = contiguous_ranges(counts)
z = similar(x)
p = contiguous_invperm(x, counts)
z[p] .= x
```

then elements of `z[r]` are equal for `r in rs`.
"""
function contiguous_invperm(x::AbstractVector{T},
                            counts::OrderedDict{T,S}) where {T, S <: Integer}
    contiguous_invperm!(similar(x, S), x, counts)
end

end # module


text:
# Integrations with DFTK.jl
# This integration should ultimately live within the DFTK package itself

# -----------------------------------------------------------------------------
# Integration with AtomsBase
# -----------------------------------------------------------------------------

# TODO: should a version of this be integrated into AtomsBase?
# ? Is this logic even correct?
function fractional_position(pos::SVector{3,<:Unitful.Length}, box::SVector{3,<:SVector{3,<:Unitful.Length}})
    [pos[i] / box[i][i] for i ∈ 1:3]
end

# TODO: support multiple species using atom data
function dftk_atoms(system::AbstractSystem{3}, element::DFTK.Element)
    [element => Fix2(fractional_position, bounding_box(system)).(position(system))]
end

# -----------------------------------------------------------------------------
# Integration with InteratomicPotentials
# -----------------------------------------------------------------------------

@kwdef struct DFTKPotential <: ArbitraryPotential
    psp::ElementPsp
    lattice
    Ecut::Real
    kgrid::AbstractVector{<:Integer}
    n_bands::Union{Integer,Nothing} = nothing
    tol::Union{AbstractFloat,Nothing} = nothing
    damping::Union{AbstractFloat,Nothing} = nothing
    mixing::Union{Mixing,Nothing} = nothing
    previous_scfres::RefValue{Any} = Ref{Any}()
    potential_energy_cache::Dict{Float64,Float64} = Dict{Float64,Float64}()
end
function DFTKPotential(
    psp::ElementPsp,
    lattice,
    Ecut::Unitful.Energy,
    kgrid::AbstractVector{<:Integer},
    n_bands::Union{Integer,Nothing} = nothing,
    tol::Union{AbstractFloat,Nothing} = nothing,
    damping::Union{AbstractFloat,Nothing} = nothing,
    mixing::Union{Mixing,Nothing} = nothing,
    previous_scfres::RefValue{Any} = Ref{Any}(),
    potential_energy_cache::Dict{Float64,Float64} = Dict{Float64,Float64}()
)
    DFTKPotential(psp, austrip.(lattice), austrip(Ecut), kgrid, n_bands, tol, damping, mixing, previous_scfres, potential_energy_cache)
end

function InteratomicPotentials.potential_energy(system::AbstractSystem, potential::DFTKPotential)
    calculate_scf(system, potential).energies.total
end

function InteratomicPotentials.potential_energy(system::DynamicSystem, potential::DFTKPotential)
    get!(potential.potential_energy_cache, austrip(system.time)) do
        calculate_scf(system, potential).energies.total
    end
end

function InteratomicPotentials.force(system::AbstractSystem, potential::DFTKPotential)
    # TODO: support multiple species
    compute_forces_cart(calculate_scf(system, potential))[1]
end

function InteratomicPotentials.force(system::DynamicSystem, potential::DFTKPotential)
    scf = calculate_scf(system, potential)
    potential.potential_energy_cache[austrip(system.time)] = scf.energies.total
    # TODO: support multiple species
    compute_forces_cart(scf)[1]
end

function calculate_scf(system::AbstractSystem, potential::DFTKPotential)
    model = model_LDA(potential.lattice, dftk_atoms(system, potential.psp))
    basis = PlaneWaveBasis(model; Ecut = potential.Ecut, kgrid = potential.kgrid)

    extra_args = isassigned(potential.previous_scfres) ? (ψ = potential.previous_scfres[].ψ, ρ = potential.previous_scfres[].ρ) : (;)
    scfres = self_consistent_field(basis; extra_args..., (f => getfield(potential, f) for f ∈ (:n_bands, :tol, :damping, :mixing) if getfield(potential, f) !== nothing)...)
    potential.previous_scfres[] = scfres
end

# -----------------------------------------------------------------------------
# Miscelaneous functions
# -----------------------------------------------------------------------------

function analyze_convergence(system::AbstractSystem, potential::DFTKPotential, cutoffs::AbstractVector{<:Unitful.Energy})
    energies = Vector{Float64}()
    for Ecut ∈ cutoffs
        parameters = DFTKPotential(
            psp = potential.psp,
            lattice = potential.lattice,
            Ecut = Ecut,
            kgrid = potential.kgrid,
            n_bands = potential.n_bands,
            tol = potential.tol,
            damping = potential.damping,
            mixing = potential.mixing
        )
        @info "Ecut: $(Ecut)"
        scfres = calculate_scf(system, parameters)
        push!(energies, scfres.energies.total)
    end

    plot(
        title = "DFTK Analysis",
        xlab = "Ecut",
        ylab = "Total Energy",
        legend = false,
        cutoffs,
        energies * ENERGY_UNIT
    )
end


text:
# References:
# [1] ADS111x Ultra-Small, Low-Power, I2C-Compatible, 860-SPS, 16-Bit ADCs
#     SBAS444D –MAY 2009–REVISED JANUARY 2018

module PiADS111x

export ads_open

using PiGPIOC
import PiGPIOC.gpioInitialise
import PiGPIOC.i2cOpen
import PiGPIOC.i2cReadWordData
import PiGPIOC.i2cWriteWordData
# See http://abyz.me.uk/rpi/pigpio/cif.html

const ADS_CONVERSION_REGISTER = 0                       # [1, 9.6.1, p27]
const ADS_CONFIG_REGISTER     = 1                       # [1, 9.6.3, p28]
const ADS_CONFIG_DEFAULT      = 0x8583                  # [1, 9.6.3, p28]
const ADS_PGA_4V096           = 0x0001 <<  9            # [1, Table 8, p28]
const ADS_OS                  = 0x0000 << 15            # [1, Table 8, p28]
const ADS_COMP_QUE_DISABLE    = 0x0003 <<  0            # [1, Table 8, p28]

struct ADS1115 <: AbstractChannel{UInt16}
    i2c::Cint
end

"""
    ads_open(;bus=1, address=0x480)::ADS1115

Connect to ADS1115 at `address` on i2c `bus`.
"""
function ads_open(;bus=1, address=0x48)

    @assert bus >= 0

    res = gpioInitialise()
    @assert(res != PiGPIOC.PI_INIT_FAILED)

    i2c = i2cOpen(bus, address, 0)
    @assert i2c >= 0

    ads = ADS1115(i2c)
    ads_config(ads)
    ads
end

"""
    ads_config(::ADS1115 , [mux = 0])

Configure using `mux` bits. See MUX in [1, Table 8, p28].
"""
function ads_config(ads::ADS1115, mux=UInt16(0))
    conf = ADS_OS |
           ADS_PGA_4V096 |
           ADS_COMP_QUE_DISABLE |
           mux
    ads_write(ads, ADS_CONFIG_REGISTER, conf)
    x = ads_read(ads, ADS_CONFIG_REGISTER)
    @assert x == conf
    nothing
end

"""
    getindex(::ADS1115, i)

Run single Analog to Digital Conversion for pin `i` (0:3)
"""
function Base.getindex(ads::ADS1115, i)

    @assert i in 0:3 "ADS1115 input `i` must be in range 0:3"
    mux = UInt16(0x4000 | (i << 12)) # [1, Table 8, p28]
    ads_config(ads, mux)
    sleep(1/8) # 8 SPS [1, Table 8, p28]
    while ads_is_busy(ads)
        yield()
    end
    ads[]
end

"""
    getindex(::ADS1115)

Read latest value from ADS_CONVERSION_REGISTER.
"""
function Base.getindex(ads::ADS1115)
    n = signed(ads_read(ads, ADS_CONVERSION_REGISTER))
    V = (n * 125) / 1000000 # [1, Table 3, p17]
end

"""
    ads_write(::ADS1115, register, value)

Write 16-bit `value` to `register` on ADS1115 ADC.
"""
function ads_write(ads, register, v)
    @assert register in 0:3
    err = i2cWriteWordData(ads.i2c, register, bswap(UInt16(v)))
    @assert err == 0
    nothing
end

"""
    ads_read(::ADS1115, register))

Read 16-bit value from `register` on ADS1115 ADC.
"""
function ads_read(ads, register)
    @assert register in 0:3
    n = i2cReadWordData(ads.i2c, register)
    @assert n >= 0
    bswap(UInt16(n))
end

"""
    ads_is_busy(::ADS1115)

Is the ADS1115 ADC currently performing a conversion?
"""
ads_is_busy(ads) = (ads_read(ads, ADS_CONFIG_REGISTER) & ADS_OS) == 0

function ads_demo()
    ads = ads_open()
    while true
        @show [ads[n] for n in 0:3]
    end
end

end # module


text:
function get_extremes(x::Vector{T}, y::Vector{T}, z::Vector{T}) where {T<:AbstractFloat}
    minx, maxx = minimum(x), maximum(x)
    miny, maxy = minimum(y), maximum(y)
    minz, maxz = minimum(z), maximum(z)
    p::Vector{T} = [minx, minx, minx, minx, maxx, maxx, maxx, maxx]
    q::Vector{T} = [miny, miny, maxy, maxy, miny, miny, maxx, maxx]
    r::Vector{T} = [minz, maxz, minz, maxz, minz, maxz, minz, maxz]
    return p, q, r
end

function get_translation_parameters(hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer}
    x = hook.xoff
    y = hook.yoff
    z = hook.zoff
    return x, y, z
end

function get_translation_parameters(transform::Transform{T,S}) where {T<:AbstractFloat, S<:Integer}
    x = transform.matrix[1,4]
    y = transform.matrix[2,4]
    z = transform.matrix[3,4]
    return x, y, z
end

function get_extremes(assets::Array{Asset{T,S}, 1}, geometry_bvhs::Array{Bvh{T,S,A,E}, 1}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    n = size(assets, 1)
    boxes::Array{Float64, 2} = zeros(6, n)
    for i=1:n
        p = assets[i].oidx

        xmin = geometry_bvhs[p].xmin[1]
        ymin = geometry_bvhs[p].ymin[1]
        zmin = geometry_bvhs[p].zmin[1]

        xmax = geometry_bvhs[p].xmax[1]
        ymax = geometry_bvhs[p].ymax[1]
        zmax = geometry_bvhs[p].zmax[1]

        x1, y1, z1 = transform_forward(xmin, ymin, zmin, assets[i])
        x2, y2, z2 = transform_forward(xmin, ymin, zmax, assets[i])
        x3, y3, z3 = transform_forward(xmin, ymax, zmin, assets[i])
        x4, y4, z4 = transform_forward(xmin, ymax, zmax, assets[i])
        x5, y5, z5 = transform_forward(xmax, ymin, zmin, assets[i])
        x6, y6, z6 = transform_forward(xmax, ymin, zmax, assets[i])
        x7, y7, z7 = transform_forward(xmax, ymax, zmin, assets[i])
        x8, y8, z8 = transform_forward(xmax, ymax, zmax, assets[i])

        txmin = min(x1, x2, x3, x4, x5, x6, x7, x8)
        tymin = min(y1, y2, y3, y4, y5, y6, y7, y8)
        tzmin = min(z1, z2, z3, z4, z5, z6, z7, z8)

        txmax = max(x1, x2, x3, x4, x5, x6, x7, x8)
        tymax = max(y1, y2, y3, y4, y5, y6, y7, y8)
        tzmax = max(z1, z2, z3, z4, z5, z6, z7, z8)

        boxes[:,i] = [txmin, tymin, tzmin, txmax, tymax, tzmax]
    end
    xmin, xmax = minimum(boxes[1,:]), maximum(boxes[4,:])
    ymin, ymax = minimum(boxes[2,:]), maximum(boxes[5,:])
    zmin, zmax = minimum(boxes[3,:]), maximum(boxes[6,:])
    x = vcat(xmin, xmax)
    y = vcat(ymin, ymax)
    z = vcat(zmin, zmax)
    p, q, r = get_extremes(x, y, z)
    return p, q, r
end

# these functions below are not particularly accurate, so use with caution! They do not consider diameters of disks, cones, cylinders, or assemblies, and no orientations/scaling...

function get_extremes(disk::Disk{T,S,A,E}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    x = disk.p1x
    y = disk.p1y
    z = disk.p1z
    p, q, r = get_extremes(x, y, z)
    return p, q, r
end

function get_extremes(cone::Cone{T,S,A,E}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    x = cone.m1x
    y = cone.m1y
    z = cone.m1z
    p, q, r = get_extremes(x, y, z)
    return p, q, r
end

function get_extremes(cylinder::Cylinder{T,S,A,E}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    x = cylinder.m1x
    y = cylinder.m1y
    z = cylinder.m1z
    p, q, r = get_extremes(x, y, z)
    return p, q, r
end

function get_extremes(mesh::Mesh{T,S,A,E}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    x = vcat(mesh.v1x, mesh.v2x, mesh.v3x)
    y = vcat(mesh.v1y, mesh.v2y, mesh.v3y)
    z = vcat(mesh.v1z, mesh.v2z, mesh.v3z)
    p, q, r = get_extremes(x, y, z)
    return p, q, r
end

function get_extremes(geometry::Geometry{T,S,A,E}, hook::Hook{T,S}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
    x, y, z = get_extremes(geometry)
    x, y, z = transform_forward(x, y, z, hook)
    x, y, z = get_extremes(x, y, z)
    return x, y, z
end

#function get_extremes(assets::Array{Asset{T,S},1}, geometry_bvhs::AbstractArray{Bvh{T,S,A,E},1}) where {T<:AbstractFloat, S<:Integer, A<:AbstractArray{T,1}, E<:AbstractArray{S,1}}
#    n = size(assets,1)
#    x = Float64[]
#    y = Float64[]
#    z = Float64[]
#    for i=1:n
#        asset = assets[i]
#        bvh = geometry_bvhs[asset.oidx]
#        push!(x, bvh.xmin...)
#        push!(x, bvh.xmax...)
#        push!(y, bvh.ymin...)
#        push!(y, bvh.ymax...)
#        push!(z, bvh.zmin...)
#        push!(z, bvh.zmax...)
#    end
#    p, q, r = get_extremes(x, y, z)
#    return p, q, r
#end

#function get_extremes(assets::Array{Asset{T,S},1}) where {T<:AbstractFloat, S<:Integer}
#    n = size(assets,1)
#    x = Float64[]
#    y = Float64[]
#    z = Float64[]
#    for i=1:n
#        asset = assets[i]
#        a, b, c = get_translation_parameters(asset)
#        push!(x, a)
#        push!(y, b)
#        push!(z, c)
#    end
#    p, q, r = get_extremes(x, y, z)
#    return p, q, r
#end

max_stars_repo_name: martinvanleeuwen/RenderJay

max_issues_repo_name: martinvanleeuwen/RenderJay

max_forks_repo_name: martinvanleeuwen/RenderJay


text:
using LMMP
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
C = trunc.(Int,readdlm("C.txt"));
A = trunc.(Int,readdlm("A.txt"));
b = trunc.(Int,readdlm("b.txt"));
d = trunc.(Int,readdlm("d.txt"));

temp = LMMPSolver(A,b,C,d,XLB = zeros(400),XUB = ones(400))
@test temp[1] == 6256.816326530623


text:
## --- Custom display functions

    # Custom pretty-printing for colormaps
    function display(x::AllColormaps)
        println("AllColormaps:")
        for name in fieldnames(AllColormaps)
            println("  $name")
            display(getfield(x, name))
        end
    end

    # Custom pretty-printing for York fit results
    function display(x::YorkFit)
        print("York Fit y = a + bx:\n  intercept a: $(x.intercept) ± $(x.intercept_sigma) (1σ)")
        print("\n  slope b    : $(x.slope) ± $(x.slope_sigma) (1σ)\n  MSWD       : $(x.mswd)\n")
    end

## ---

max_stars_repo_name: brenhinkeller/StatGeochemBase

max_issues_repo_name: brenhinkeller/StatGeochemBase

max_forks_repo_name: brenhinkeller/StatGeochemBase


text:
#####
##### Momentum and tracer advective flux operators for upwind-biased advection schemes
#####
##### See topologically_conditional_interpolation.jl for an explanation of the underscore-prepended
##### functions _symmetric_interpolate_*, _left_biased_interpolate_*, and _right_biased_interpolate_*.
#####

const Upwind = AbstractUpwindBiasedAdvectionScheme

@inline upwind_biased_product(ũ, ψᴸ, ψᴿ) = ((ũ + abs(ũ)) * ψᴸ + (ũ - abs(ũ)) * ψᴿ) / 2

#####
##### Momentum advection operators
#####
##### Note the convention "momentum_flux_AB" corresponds to the advection _of_ B _by_ A.
#####

@inline function momentum_flux_uu(i, j, k, grid, scheme::Upwind, u)

    ũ  =    _symmetric_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme, u)
    uᴸ =  _left_biased_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme, u)
    uᴿ = _right_biased_interpolate_xᶜᵃᵃ(i, j, k, grid, scheme, u)

    return Axᵃᵃᶠ(i, j, k, grid) * upwind_biased_product(ũ, uᴸ, uᴿ)
end

@inline function momentum_flux_uv(i, j, k, grid, scheme::Upwind, u, v)

    ṽ  =    _symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme, v)
    uᴸ =  _left_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme, u)
    uᴿ = _right_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme, u)

    return Ayᵃᵃᶠ(i, j, k, grid) * upwind_biased_product(ṽ, uᴸ, uᴿ)
end

@inline function momentum_flux_uw(i, j, k, grid, scheme::Upwind, u, w)

    w̃  =    _symmetric_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme, w)
    uᴸ =  _left_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme, u)
    uᴿ = _right_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme, u)

    return Azᵃᵃᵃ(i, j, k, grid) * upwind_biased_product(w̃, uᴸ, uᴿ)
end

@inline function momentum_flux_vu(i, j, k, grid, scheme::Upwind, u, v)

    ũ  =    _symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme, u)
    vᴸ =  _left_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme, v)
    vᴿ = _right_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme, v)

    return Axᵃᵃᶠ(i, j, k, grid) * upwind_biased_product(ũ, vᴸ, vᴿ)
end

@inline function momentum_flux_vv(i, j, k, grid, scheme::Upwind, v)

    ṽ  =    _symmetric_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme, v)
    vᴸ =  _left_biased_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme, v)
    vᴿ = _right_biased_interpolate_yᵃᶜᵃ(i, j, k, grid, scheme, v)

    return Ayᵃᵃᶠ(i, j, k, grid) * upwind_biased_product(ṽ, vᴸ, vᴿ)
end

@inline function momentum_flux_vw(i, j, k, grid, scheme::Upwind, v, w)

    w̃  =    _symmetric_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme, w)
    vᴸ =  _left_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme, v)
    vᴿ = _right_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme, v)

    return Azᵃᵃᵃ(i, j, k, grid) * upwind_biased_product(w̃, vᴸ, vᴿ)
end

@inline function momentum_flux_wu(i, j, k, grid, scheme::Upwind, u, w)

    ũ  =    _symmetric_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme, u)
    wᴸ =  _left_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme, w)
    wᴿ = _right_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme, w)

    return Axᵃᵃᶠ(i, j, k, grid) * upwind_biased_product(ũ, wᴸ, wᴿ)
end

@inline function momentum_flux_wv(i, j, k, grid, scheme::Upwind, v, w)

    ṽ  =    _symmetric_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme, v)
    wᴸ =  _left_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme, w)
    wᴿ = _right_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme, w)

    return Ayᵃᵃᶠ(i, j, k, grid) * upwind_biased_product(ṽ, wᴸ, wᴿ)
end

@inline function momentum_flux_ww(i, j, k, grid, scheme::Upwind, w)

    w̃  =    _symmetric_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme, w)
    wᴸ =  _left_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme, w)
    wᴿ = _right_biased_interpolate_zᵃᵃᶜ(i, j, k, grid, scheme, w)

    return Azᵃᵃᵃ(i, j, k, grid) * upwind_biased_product(w̃, wᴸ, wᴿ)
end

#####
##### Tracer advection operators
#####

@inline function advective_tracer_flux_x(i, j, k, grid, scheme::Upwind, u, c)

    @inbounds ũ = u[i, j, k]
    cᴸ =  _left_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme, c)
    cᴿ = _right_biased_interpolate_xᶠᵃᵃ(i, j, k, grid, scheme, c)

    Axᵃᵃᶠ(i, j, k, grid) * upwind_biased_product(ũ, cᴸ, cᴿ)
end

@inline function advective_tracer_flux_y(i, j, k, grid, scheme::Upwind, v, c)

    @inbounds ṽ = v[i, j, k]
    cᴸ =  _left_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme, c)
    cᴿ = _right_biased_interpolate_yᵃᶠᵃ(i, j, k, grid, scheme, c)

    return Ayᵃᵃᶠ(i, j, k, grid) * upwind_biased_product(ṽ, cᴸ, cᴿ)
end

@inline function advective_tracer_flux_z(i, j, k, grid, scheme::Upwind, w, c)

    @inbounds w̃ = w[i, j, k]
    cᴸ =  _left_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme, c)
    cᴿ = _right_biased_interpolate_zᵃᵃᶠ(i, j, k, grid, scheme, c)

    return Azᵃᵃᵃ(i, j, k, grid) * upwind_biased_product(w̃, cᴸ, cᴿ)
end


text:
using Random, LinearAlgebra, BlockArrays
using BlockArrays: _BlockArray
using Stheno: BlockLowerTriangular, BlockUpperTriangular
using LinearAlgebra: ldiv!, \

function solve_BlockTriangular_Vector_tests(rng, T, x, TriangleType)
    Ps, Tmat, xvec = blocksizes(T, 1), TriangleType(Matrix(T)), Vector(x)

    @test ldiv!(T, copy(x)) isa BlockVector
    @test Vector(ldiv!(T, copy(x))) ≈ ldiv!(Tmat, copy(xvec))
    @test T \ x isa BlockVector
    @test ldiv!(T, copy(x)) == T \ x

    z̄ = _BlockArray([randn(rng, P) for P in Ps], Ps)

    z, back = Zygote.pullback(\, T, x)
    L̄, x̄ = back(z̄)
    @test L̄ isa BlockMatrix
    @test x̄ isa BlockVector

    z_vec, back_mat = Zygote.pullback(\, Tmat, xvec)
    L̄_mat, x̄_vec = back_mat(Vector(z̄))
    @test z_vec ≈ Vector(z)
    @test L̄_mat ≈ Matrix(L̄)
    @test x̄_vec ≈ Vector(x̄)
end

function solve_BlockTriangular_Matrix_tests(rng, T, X, TriangleType)
    Ps, Qs = blocksizes(X, 1), blocksizes(X, 2)
    Tmat, Xmat = TriangleType(Matrix(T)), Matrix(X)

    @test ldiv!(T, copy(X)) isa BlockMatrix
    @test Matrix(ldiv!(T, copy(X))) ≈ ldiv!(Tmat, copy(Xmat))
    @test T \ X isa BlockMatrix
    @test ldiv!(T, copy(X)) == T \ X

    Z̄ = _BlockArray([randn(rng, P, Q) for P in Ps, Q in Qs], Ps, Qs)

    Z, back = Zygote.pullback(\, T, X)
    T̄, X̄ = back(Z̄)
    @test T̄ isa BlockMatrix
    @test X̄ isa BlockMatrix

    Z_vec, back_mat = Zygote.pullback(\, Tmat, Xmat)
    T̄_mat, X̄_mat = back_mat(Matrix(Z̄))
    @test Z_vec ≈ Matrix(Z)
    @test T̄_mat ≈ Matrix(T̄)
    @test X̄_mat ≈ Matrix(X̄)
end

@testset "triangular" begin
    @testset "adjoint / transpose" begin
        rng, Ps = MersenneTwister(123456), [5, 4]
        X = _BlockArray([randn(rng, P, P′) for P in Ps, P′ in Ps], Ps, Ps)
        Xmat = Matrix(X)
        for foo in [adjoint, transpose]
            @test foo(UpperTriangular(X)) isa BlockLowerTriangular
            @test Matrix(foo(UpperTriangular(X))) == collect(foo(UpperTriangular(Xmat)))
            @test foo(LowerTriangular(X)) isa BlockUpperTriangular
            @test Matrix(foo(LowerTriangular(X))) == collect(foo(LowerTriangular(Xmat)))
        end
    end
    @testset "construction" begin
        rng, Ps = MersenneTwister(123456), [5, 6, 7]
        @testset "BlockLowerTriangular" begin
            L = LowerTriangular(randn(rng, sum(Ps), sum(Ps)))
            @test BlockArray(L, Ps, Ps) isa BlockLowerTriangular
            @test BlockArray(L, Ps, Ps) == L
        end
        @testset "BlockUpperTriangular" begin
            U = UpperTriangular(randn(rng, sum(Ps), sum(Ps)))
            @test BlockArray(U, Ps, Ps) isa BlockUpperTriangular
            @test BlockArray(U, Ps, Ps) == U
        end
    end
    @testset "copy" begin
        rng, Ps = MersenneTwister(123456), [5, 4]
        X = BlockArray(randn(rng, sum(Ps), sum(Ps)), Ps, Ps)

        @testset "BlockLowerTriangular" begin
            L = LowerTriangular(X)
            @test copy(L) isa BlockLowerTriangular
            @test copy(L) == L
        end
        @testset "BlockUpperTriangular" begin
            U = UpperTriangular(X)
            @test copy(U) isa BlockUpperTriangular
            @test copy(U) == U
        end
    end
    @testset "mul! / *" begin
        rng, Ps, Qs = MersenneTwister(123456), [5, 4, 3], [7, 6, 5]
        A = _BlockArray([randn(rng, P, P′) for P in Ps, P′ in Ps], Ps, Ps)
        L, U = LowerTriangular(A), UpperTriangular(A)
        x = _BlockArray([randn(rng, P) for P in Ps], Ps)
        X = _BlockArray([randn(rng, P, Q) for P in Ps, Q in Qs], Ps, Qs)

        @testset "LowerTriangular" begin
            @testset "Matrix-Vector" begin
                dense_BlockMatrix_BlockVector_mul_tests(rng, L, x)
            end
            @testset "Matrix-Matrix" begin
                dense_BlockMatrix_BlockMatrix_mul_tests(rng, L, X)
            end
        end
        @testset "UpperTriangular" begin
            @testset "Matrix-Vector" begin
                dense_BlockMatrix_BlockVector_mul_tests(rng, U, x)
            end
            @testset "Matrix-Matrix" begin
                dense_BlockMatrix_BlockMatrix_mul_tests(rng, U, X)
            end
        end
    end
    @testset "ldiv! / \\" begin
        rng, Ps, Qs = MersenneTwister(123456), [5, 4, 3], [7, 6, 5]
        A = _BlockArray([randn(rng, P, P′) for P in Ps, P′ in Ps], Ps, Ps)
        L, U = LowerTriangular(A), UpperTriangular(A)
        x = _BlockArray([randn(rng, P) for P in Ps], Ps)
        X = _BlockArray([randn(rng, P, Q) for P in Ps, Q in Qs], Ps, Qs)

        @testset "LowerTriangular" begin
            @testset "Matrix-Vector" begin
                solve_BlockTriangular_Vector_tests(rng, L, x, LowerTriangular)
            end
            @testset "Matrix-Matrix" begin
                solve_BlockTriangular_Matrix_tests(rng, L, X, LowerTriangular)
            end
        end
        @testset "UpperTriangular" begin
            @testset "Matrix-Vector" begin
                solve_BlockTriangular_Vector_tests(rng, U, x, UpperTriangular)
            end
            @testset "Matrix-Matrix" begin
                solve_BlockTriangular_Matrix_tests(rng, U, X, UpperTriangular)
            end
        end
    end
end


text:
using Statistics
using Logging
using Printf
using DataDeps
using JLD2

using Oceananigans
using Oceananigans.Units

using Oceananigans.Diagnostics: accurate_cell_advection_timescale

Logging.global_logger(OceananigansLogger())

#####
##### Progress monitor
#####

mutable struct Progress
    interval_start_time :: Float64
end

function (p::Progress)(sim)
    wall_time = (time_ns() - p.interval_start_time) * 1e-9

    @info @sprintf("Time: %s, iteration: %d, max(|u⃗|): (%.2e, %.2e) m/s, extrema(η): (min=%.2e, max=%.2e), CFL: %.2e, wall time: %s",
                   prettytime(sim.model.clock.time),
                   sim.model.clock.iteration,
                   maximum(abs, sim.model.velocities.u),
                   maximum(abs, sim.model.velocities.v),
                   minimum(sim.model.free_surface.η),
                   maximum(sim.model.free_surface.η),
                   sim.parameters.cfl(sim.model),
                   prettytime(wall_time))

    p.interval_start_time = time_ns()

    return nothing
end

#####
##### Script starts here
#####

ENV["DATADEPS_ALWAYS_ACCEPT"] = "true"

dd = DataDep("cubed_sphere_32_grid",
    "Conformal cubed sphere grid with 32×32 grid points on each face",
    "https://github.com/CliMA/OceananigansArtifacts.jl/raw/main/cubed_sphere_grids/cubed_sphere_32_grid.jld2",
    "b1dafe4f9142c59a2166458a2def743cd45b20a4ed3a1ae84ad3a530e1eff538" # sha256sum
)

DataDeps.register(dd)

cs32_filepath = datadep"cubed_sphere_32_grid/cubed_sphere_32_grid.jld2"

# Central (λ, φ) for each face of the cubed sphere
central_longitude = (0, 90,  0, 180, -90,   0)
central_latitude  = (0,  0, 90,   0,   0, -90)

function cubed_sphere_surface_gravity_waves(; face_number)

    H = 4kilometers
    grid = ConformalCubedSphereGrid(cs32_filepath, Nz=1, z=(-H, 0))

    ## Model setup

    model = HydrostaticFreeSurfaceModel(
                      grid = grid,
        momentum_advection = nothing,
              free_surface = ExplicitFreeSurface(gravitational_acceleration=0.1),
                  coriolis = nothing,
                   closure = nothing,
                   tracers = nothing,
                  buoyancy = nothing
    )

    ## Initial condition:
    ## Very small sea surface height perturbation so the resulting dynamics are well-described
    ## by a linear free surface.

    A  = 1e-5 * H  # Amplitude of the perturbation
    λ₀ = central_longitude[face_number]
    φ₀ = central_latitude[face_number]
    Δλ = 15  # Longitudinal width
    Δφ = 15  # Latitudinal width

    η′(λ, φ) = A * exp(- (λ - λ₀)^2 / Δλ^2) * exp(- (φ - φ₀)^2 / Δφ^2)

    Oceananigans.set!(model, η=η′)

    ## Simulation setup

    Δt = 10minutes

    cfl = CFL(Δt, accurate_cell_advection_timescale)

    simulation = Simulation(model,
                        Δt = Δt,
                 stop_time = 25days,
        iteration_interval = 1,
                  progress = Progress(time_ns()),
                parameters = (; cfl)
    )

    fields_to_check = (
        u = model.velocities.u,
        v = model.velocities.v,
        η = model.free_surface.η,
        Gu = model.timestepper.Gⁿ.u,
        Gv = model.timestepper.Gⁿ.v,
        Gη = model.timestepper.Gⁿ.η
    )

    simulation.diagnostics[:state_checker] =
        StateChecker(model, fields=fields_to_check, schedule=IterationInterval(1))

    output_fields = merge(model.velocities, (η=model.free_surface.η,))

    simulation.output_writers[:fields] =
        JLD2OutputWriter(model, output_fields,
            schedule = TimeInterval(1hour),
              prefix = "cubed_sphere_surface_gravity_waves_face$face_number",
               force = true)

    run!(simulation)

    return simulation
end

include("animate_on_map_projection.jl")

function run_cubed_sphere_surface_gravity_waves_validation()

    for f in 1:6
        cubed_sphere_surface_gravity_waves(face_number=f)
    end

    projections = [
        ccrs.NearsidePerspective(central_longitude=0,   central_latitude=30),
        ccrs.NearsidePerspective(central_longitude=180, central_latitude=-30)
    ]

    for f in 1:6
        animate_surface_gravity_waves(face_number=f, projections=projections)
    end
end


text:
# A Modeling Language Package for Solving Data Envelopment Analysis Problems
module FrontierEfficiencyAnalysis_v2

# An algebraic modeling language for Julia
import JuMP
# Solver-independent functions and low-level interface for Mathematical Programming
import MathProgBase
# An example solver
import Gurobi

# function
export solveDEA

# solving dea problems by a smart algorithm
# it is computationally efficient and can solve large scale DEA problems via any solvers (even with LP size limitation)
include("solveDEA.jl")

end


text:
solution_ml = sem_fit(model_ml)
bs = se_bootstrap(solution_ml; n_boot = 20)
se = se_hessian(solution_ml)

update_partable!(partable, solution_ml, se, :se)
update_partable!(partable, solution_ml, bs, :se_boot)


text:
module CassetteZipper

using ..Ahorn, Maple

@pardef CassetteZipMover(x1::Integer, y1::Integer, x2::Integer=x1 + 16, y2::Integer=y1, width::Integer=8, height::Integer=8, theme="Normal", index::Integer=0, tempo::Number=1) = Entity("brokemiahelper/cassetteZipMover", x=x1, y=y1, nodes=Tuple{Int, Int}[(x2, y2)], width=width, height=height, theme=theme, index=index, tempo=tempo)

colorNames = Dict{String, Int}(
    "Blue" => 0,
    "Rose" => 1,
    "Bright Sun" => 2,
    "Malachite" => 3
)

const placements = Ahorn.PlacementDict(
    "Cassette Zip Mover ($index - $color)" => Ahorn.EntityPlacement(
        CassetteZipMover,
        "rectangle",
        Dict{String, Any}(
            "index" => index,
            "theme" => "Normal"
        ),
        function(entity)
            entity.data["nodes"] = [(Int(entity.data["x"]) + Int(entity.data["width"]) + 8, Int(entity.data["y"]))]
        end
    ) for (color, index) in colorNames
)

Ahorn.nodeLimits(entity::CassetteZipMover) = 0, 1

Ahorn.minimumSize(entity::CassetteZipMover) = 16, 16
Ahorn.resizable(entity::CassetteZipMover) = true, true

Ahorn.editingOptions(entity::CassetteZipMover) = Dict{String, Any}(
    "index" => colorNames,
    "theme" => Maple.zip_mover_themes
)

colors = Dict{Integer, Ahorn.colorTupleType}(
    1 => (240, 73, 190, 255) ./ 255,
	2 => (252, 220, 58, 255) ./ 255,
	3 => (56, 224, 78, 255) ./ 255,
)

defaultColor = (73, 170, 240, 255) ./ 255

borderMultiplier = (0.9, 0.9, 0.9, 1)

function Ahorn.selection(entity::CassetteZipMover)
    x, y = Ahorn.position(entity)
    nx, ny = Int.(entity.data["nodes"][1])

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

    return [Ahorn.Rectangle(x, y, width, height), Ahorn.Rectangle(nx + floor(Int, width / 2) - 5, ny + floor(Int, height / 2) - 5, 10, 10)]
end

function getTextures(entity::CassetteZipMover)
    theme = lowercase(get(entity, "theme", "normal"))

    if theme == "moon"
        return "objects/zipmover/moon/block", "objects/zipmover/moon/light01", "objects/zipmover/moon/cog"
    end

    return "objects/zipmover/block", "objects/zipmover/light01", "objects/zipmover/cog"
end

ropeColor = (102, 57, 49) ./ 255

function renderCassetteZipMover(ctx::Ahorn.Cairo.CairoContext, entity::CassetteZipMover)
    x, y = Ahorn.position(entity)
    nx, ny = Int.(entity.data["nodes"][1])

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

    block, light, cog = getTextures(entity)
    lightSprite = Ahorn.getSprite(light, "Gameplay")

    tilesWidth = div(width, 8)
    tilesHeight = div(height, 8)

    cx, cy = x + width / 2, y + height / 2
    cnx, cny = nx + width / 2, ny + height / 2

    length = sqrt((x - nx)^2 + (y - ny)^2)
    theta = atan(cny - cy, cnx - cx)

    Ahorn.Cairo.save(ctx)

    Ahorn.translate(ctx, cx, cy)
    Ahorn.rotate(ctx, theta)

    Ahorn.setSourceColor(ctx, ropeColor)
    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1);

    # Offset for rounding errors
    Ahorn.move_to(ctx, 0, 4 + (theta <= 0))
    Ahorn.line_to(ctx, length, 4 + (theta <= 0))

    Ahorn.move_to(ctx, 0, -4 - (theta > 0))
    Ahorn.line_to(ctx, length, -4 - (theta > 0))

    Ahorn.stroke(ctx)

    Ahorn.Cairo.restore(ctx)

    Ahorn.drawRectangle(ctx, x + 2, y + 2, width - 4, height - 4, (0.0, 0.0, 0.0, 1.0))
    Ahorn.drawSprite(ctx, cog, cnx, cny, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))

    for i in 2:tilesWidth - 1
        Ahorn.drawImage(ctx, block, x + (i - 1) * 8, y, 8, 0, 8, 8, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))
        Ahorn.drawImage(ctx, block, x + (i - 1) * 8, y + height - 8, 8, 16, 8, 8, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))
    end

    for i in 2:tilesHeight - 1
        Ahorn.drawImage(ctx, block, x, y + (i - 1) * 8, 0, 8, 8, 8, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))
        Ahorn.drawImage(ctx, block, x + width - 8, y + (i - 1) * 8, 16, 8, 8, 8, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))
    end

    Ahorn.drawImage(ctx, block, x, y, 0, 0, 8, 8, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))
    Ahorn.drawImage(ctx, block, x + width - 8, y, 16, 0, 8, 8, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))
    Ahorn.drawImage(ctx, block, x, y + height - 8, 0, 16, 8, 8, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))
    Ahorn.drawImage(ctx, block, x + width - 8, y + height - 8, 16, 16, 8, 8, tint=get(colors, Int(get(entity.data, "index", 0)), defaultColor))

    Ahorn.drawImage(ctx, lightSprite, x + floor(Int, (width - lightSprite.width) / 2), y)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::CassetteZipMover, room::Maple.Room)
    renderCassetteZipMover(ctx, entity)
end

end

max_stars_repo_name: Brokemia/BrokemiaHelper

max_issues_repo_name: Brokemia/BrokemiaHelper

max_forks_repo_name: Brokemia/BrokemiaHelper


text:
using Raylib

const GLSL_VERSION = 330

function build_font(filedata, baseSize, glyphCount, fonttype, padding, packmeth)
    null = convert(Ptr{Cint}, C_NULL)
    glyphs_ptr = Raylib.LoadFontData(pointer(filedata), length(filedata), 16, null, glyphCount, Int(fonttype))
    recs = Ref{Ptr{Raylib.RayRectangle}}()
    atlas = Raylib.GenImageFontAtlas(glyphs_ptr, recs, 95, 16, padding, packmeth)
    recs_ptr = recs[]
    texture = Raylib.LoadTextureFromImage(atlas)
    Raylib.UnloadImage(atlas)

    font = Raylib.RayFont(baseSize, 95, 0, texture, recs_ptr, glyphs_ptr)
    return font
end

function main()
    MAX_FONTS = 8
    screenWidth = 800
    screenHeight = 450

    Raylib.InitWindow(screenWidth, screenHeight, "raylib [text] example - SDF fonts")

    msg = "Signed Distance Fields"

    filedata = Raylib.RayFileData(joinpath(@__DIR__, "resources/anonymous_pro_bold.ttf"))

    basesize = 16
    glyphcount = 95

    default_font = build_font(filedata, basesize, 95, Raylib.FONT_DEFAULT, 4, 0)
    sdf_font = build_font(filedata, basesize, 0, Raylib.FONT_SDF, 0, 0)

    filedata = nothing

    shader = Raylib.LoadShader("", joinpath(@__DIR__, "resources/shaders/glsl$(GLSL_VERSION)/sdf.fs"))
    Raylib.SetTextureFilter(sdf_font.texture, Int(Raylib.TEXTURE_FILTER_BILINEAR))

    font_pos = Raylib.rayvector(40, screenHeight/2 - 50)
    textSize = Raylib.rayvector(0, 0)
    fontSize = 16.0
    current_font = false

    GC.gc()
    Raylib.SetTargetFPS(60)

    while !Raylib.WindowShouldClose()

        fontSize += 8Raylib.GetMouseWheelMove()
        fontSize = max(fontSize, 6)

        current_font = Raylib.IsKeyDown(Raylib.KEY_SPACE)

        textSize = if current_font
            Raylib.MeasureTextEx(sdf_font, msg, fontSize, 0)
        else
            Raylib.MeasureTextEx(default_font, msg, fontSize, 0)
        end

        font_pos = Raylib.rayvector(
            Raylib.GetScreenWidth(),
            Raylib.GetScreenHeight()
        ) ./ 2 .- textSize ./ 2 .+ (0, 80)

        Raylib.BeginDrawing()
        Raylib.ClearBackground(Raylib.RAYWHITE)

        if current_font
            Raylib.BeginShaderMode(shader)
            Raylib.DrawTextEx(sdf_font, msg, font_pos, fontSize, 0, Raylib.BLACK)
            Raylib.EndShaderMode()
            Raylib.DrawTexture(sdf_font.texture, 10, 10, Raylib.BLUE)

            Raylib.DrawText("SDF!", 320, 20, 80, Raylib.RED)
        else
            Raylib.DrawTextEx(default_font, msg, font_pos, fontSize, 0, Raylib.BLACK)
            Raylib.DrawTexture(default_font.texture, 10, 10, Raylib.BLUE)

            Raylib.DrawText("default font", 315, 40, 30, Raylib.GRAY)
        end

        Raylib.DrawText("FONT SIZE: 16.0", Raylib.GetScreenWidth() - 240, 20, 20, Raylib.DARKGRAY)
        Raylib.DrawText("RENDER SIZE: $(fontSize)", Raylib.GetScreenWidth() - 240, 50, 20, Raylib.DARKGRAY)
        Raylib.DrawText("Use MOUSE WHEEL to SCALE TEXT!", Raylib.GetScreenWidth() - 240, 90, 10, Raylib.DARKGRAY)
        Raylib.DrawText("hold space to use sdf font version!", 340, Raylib.GetScreenHeight() - 30, 20, Raylib.MAROON)
        Raylib.EndDrawing()
    end

    Raylib.UnloadFont(default_font)
    Raylib.UnloadFont(sdf_font)
    Raylib.UnloadShader(shader)
    Raylib.CloseWindow()
end


text:
using EzXML, OrdinaryDiffEq, DataFrames, CSV

function clone(repo)
    run(`git clone $url`)
end

function run_all_repos(csv_file, root, filelist=readdir(root); skip=0)
    df = DataFrame(file=String[], len = Int[], res = Int[], deps = Int[], msg = String[])
    n = 0

    for d in filelist
        path = joinpath(root, d)
        if isdir(path)
            n += run_repo(path, df; dry_run=n<skip)
            CSV.write(csv_file, df)
        end
    end
end

function run_repo(repo, df; file_limit=500000, dry_run=false)
    printstyled("processing repo $repo\n"; color=:blue)
    files = list_top_cellml_files(repo)

    for f in files
        l = filesize(f)
        msg = "OK!"
        k = 0
        m = 0
        if dry_run
            printstyled("\tskipping file $f ($l bytes)\n"; color=:blue)
        else
            printstyled("\tprocessing file $f ($l bytes)\n"; color=:blue)
            try
                if l > file_limit
                    msg = "Too big!"
                    k = 9
                    printstyled("\tfile $f ($l bytes) is too big and is skipped\n"; color=:red)
                else
                    ml = CellModel(f)
                    k = 1
                    m = length(ml.doc.xmls)
                    prob  = ODEProblem(ml, (0,1000.0))
                    k = 2
                    sol = solve(prob, TRBDF2(), dtmax=0.5)
                    k = 3
                end
            catch e
                println(e)
                msg = string(e)
            finally
                push!(df, (f, l, k, m, msg))
                printstyled("$f done with a code $k\n"; color=:green)
            end
        end
    end

    return length(files)
end

function list_top_cellml_files(repo)
    files = [f for f in readdir(repo) if endswith(f,".cellml")]

    if length(files) == 1
        return joinpath.(repo, files)
    end

    imported = Set{String}()
    for f in files
        xml = readxml(joinpath(repo, f))
        for n in CellMLToolkit.list_imports(xml)
            push!(imported, n["xlink:href"])
        end
    end

    if !isempty(imported)
        printstyled("imported files are $imported\n"; color=:yellow)
    end

    joinpath.(repo, [f for f in files if f ∉ imported])
end

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

const cellml_files_with_import = [
    "25d",
    "bugbuster",
    "butera_rinzel_smith_II_1999",
    "cardiovascular_circulation_dvad",
    "cardiovascular_circulation_ivad",
    "cardiovascular_circulation_systemic",
    "cardiovascular_circulation_westkessel",
    "cooling_hunter_crampin_CellML_2008",
    "electrecoblu",
    "faville_pullan_sanders_koh_lloyd_smith_2009",
    "garmendia-torres_goldbeter_jacquet_2007",
    "goldbeter_1991",
    "goldbeter_gonze_pourquie_2007",
    "grandi_pasqualini_bers_2010",
    "guyton_aldosterone_2008",
    "guyton_angiotensin_2008",
    "guyton_antidiuretic_hormone_2008",
    "guyton_atrial_natriuretic_peptide_2008",
    "guyton_autonomics_2008",
    "pasek_simurda_christe_2006",
    "saucerman_bers_2008",
    "shi_hose_2009",
    "smith_abdala_koizumi_rybak_paton_2007",
    "terkildsen_niederer_crampin_hunter_smith_2008",
    "virtualbacterialphotography",
    "westermark_lansner_2003"
]


text:
function calc_block_coords(cur_block::AbstractBlock, cur_rotation::Int=0)

  @assert cur_rotation in [ 0 , -1 , +1 ]

  cur_piece = cur_block.piece

  if iszero(cur_rotation)
    cur_row = cur_block.y
    cur_col = cur_block.x
  else
    cur_row = cur_block.x
    cur_col = cur_block.y

    cur_row *= -sign(cur_rotation)
    cur_col *= +sign(cur_rotation)
  end

  if iseven(cur_piece.width)
    ( cur_row < 0 ) && ( cur_row += 1 )
    ( cur_col < 0 ) && ( cur_col += 1 )
  end

  cur_row += cur_piece.i
  cur_col += cur_piece.j

  return (cur_row, cur_col)

end


text:
function FieldGeodeticPoint(arg0::RealFieldElement, arg1::RealFieldElement, arg2::RealFieldElement)
    return FieldGeodeticPoint((RealFieldElement, RealFieldElement, RealFieldElement), arg0, arg1, arg2)
end

function equals(obj::FieldGeodeticPoint, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_altitude(obj::FieldGeodeticPoint)
    return jcall(obj, "getAltitude", RealFieldElement, ())
end

function get_east(obj::FieldGeodeticPoint)
    return jcall(obj, "getEast", FieldVector3D, ())
end

function get_latitude(obj::FieldGeodeticPoint)
    return jcall(obj, "getLatitude", RealFieldElement, ())
end

function get_longitude(obj::FieldGeodeticPoint)
    return jcall(obj, "getLongitude", RealFieldElement, ())
end

function get_nadir(obj::FieldGeodeticPoint)
    return jcall(obj, "getNadir", FieldVector3D, ())
end

function get_north(obj::FieldGeodeticPoint)
    return jcall(obj, "getNorth", FieldVector3D, ())
end

function get_south(obj::FieldGeodeticPoint)
    return jcall(obj, "getSouth", FieldVector3D, ())
end

function get_west(obj::FieldGeodeticPoint)
    return jcall(obj, "getWest", FieldVector3D, ())
end

function get_zenith(obj::FieldGeodeticPoint)
    return jcall(obj, "getZenith", FieldVector3D, ())
end

function hash_code(obj::FieldGeodeticPoint)
    return jcall(obj, "hashCode", jint, ())
end

function to_string(obj::FieldGeodeticPoint)
    return jcall(obj, "toString", JString, ())
end