text:
using CPUTime
using DifferentialEquations
using Plots

f = @ode_def begin
    dx = v
    dv = 1/m*(-c*v - k*x + u)
    du = -Kp*v + Ki*(xg-x) - Kd/m*(-c*v - k*x + u)
end m c k xg Kp Ki Kd

function do_1()
    """1つめのやつ"""
    m = 1.0
    c = 1.0
    k = 1.0
    xg = 1.0
    Kp = 3
    Ki = 1
    Kd = 0.05

    p = [m, c, k, xg, Kp, Ki, Kd]

    x0 = [0.0, 0.0, Kp*(xg - 1.0)]  # 初期値を設定
    tspan = (0.0, 20.0)  # 解く範囲を設定

    prob = ODEProblem(f ,x0 ,tspan, p)
    sol = solve(prob)  # ソルバで解く

    # plot
    plot(
        sol,
        linewidth=2,
        xaxis="time [sec]",
        yaxis="",
        label=["position" "velocity" "input"]
    )

end

function do_2()
    """2つめのやつ"""
    m = 1.0
    c = 1.0
    k = 1.0
    xg = 1.0

    n = 1000  # アニメの枚数
    Kp_min, Kp_max = 9.99, 10.0
    Ki_min, Ki_max = 4.99, 5.0
    Kd_min, Kd_max = 0.0, 5.0
    Kp = [i for i in Kp_min:(Kp_max-Kp_min)/n:Kp_max]
    Ki = [i for i in Ki_min:(Ki_max-Ki_min)/n:Ki_max]
    Kd = [i for i in Kd_min:(Kd_max-Kd_min)/n:Kd_max]

    sols = []
    tspan = (0.0, 5.0)  # 解く範囲を設定

    for j in 1:n
        x0 = [0.0, 0.0, Kp[j]*(xg - 1.0)]  # 初期値を設定
        p = [m, c, k, xg, Kp[j], Ki[j], Kd[j]]
        prob = ODEProblem(f ,x0 ,tspan, p)
        sol = solve(prob)  # ソルバで解く
        push!(sols, sol)
    end

    ano_xposi = (tspan[1] + tspan[2])/2
    anim = @animate for i in 1:1:n #length(sols)
        Kp_now = round(Kp[i], digits=1)
        Ki_now = round(Ki[i], digits=1)
        Kd_now = round(Kd[i], digits=1)
        s = "Kp = " * string(Kp_now) * ", Ki = " * string(Ki_now) * ", Kd = " * string(Kd_now)
        plot(
            sols[i],
            ylims = (-1, xg*1.5), annotate = (ano_xposi, -1, s),
        )
    end
    mp4(anim, "example_4_jl.mp4", fps = 30)

end

#@time do_1()

@time do_2()

ext: jl


text:
Hecke.add_verbose_scope(:Conjugacy)
Hecke.add_assert_scope(:Conjugacy)

include("Conjugacy/Husert.jl")
include("Conjugacy/Lifting.jl")
include("Conjugacy/Random.jl")
include("Conjugacy/Conjugacy.jl")

# Element which is failing:
# A = matrix(FlintQQ, 10, 10, [0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, -49, 0, -42, 14, -9, 6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, -1, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 5, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, -1])
#
#
# A = matrix(FlintQQ, 10, 10, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 6, 9, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 6, 9, -6, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, -1, -7, -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 6])
#
# This gives the right_order assertion:
#
# A = matrix(FlintQQ, 6, 6, [0, 1, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, -125, -75, -15, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, -5])
# B = matrix(FlintQQ, 6, 6, [-22910, 2266, 45832, 2658, 1285, -226465, -16600, 1660, 33621, 17622, 350, -168655, -11375, 1125, 22750, 1080, 625, -112400, 0, 0, 0, 7, 0, 0, 0, 0, 0, 0, 7, 0, -150, 15, 303, 126, 0, -1520])
#
# A = matrix(FlintQQ, 5, 5, [0, 1, 0, 0, 0, -81, -18, 0, 0, 0, 0, 0, -9, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 3, -9]])
#
# B = matrix(FlintQQ, 5, 5, 15054, 1, -105441, -149944, -5019, 21417, -15, -149919, -212937, -7151, 1254, 0, -8787, -12480, -418, 627, 0, -4389, -6249, -209, 108, 3, -756, -1179, -39])
#
# julia> Hecke._eltseq(A) |> print
# fmpz[0, 1, 0, 0, 0, 0, 1, 0, -343, -147, -21, 0, 0, 0, 0, -7]
# julia> Hecke._eltseq(B) |> print
# fmpz[-39277, -17864, -773376, 225560, -2841374, -1292556, -55957655, 16320401, -405524, -184499, -7988112, 2329751, -1622301, -738076, -31955520, 9319917]
#
# A = matrix(FlintQQ, 3, 3, [0, 1, 0, -1, 2, 0, 0, 0, 1])
#
# B = matrix(FlintQQ, 3, 3, [-2959757959, -15202918, 41804175660, -153800, -789, 2172300, -209552500, -1076375, 2959758751])
#
# B = matrix(FlintQQ, 3, 3, [-1, 2, 2, -1, 2, 1, -1, 1, 2])


text:
@testset "Pairwise Gap Percentage" begin

    @testset "Simple" begin
        file = joinpath(pwd(), "data", "simple.fasta")
        mat = [ 0. 0.
                0. 0. ]

        (gu, gi) = pairwisegapfraction(file, FASTA)
        @test gu == mat
        @test gi == mat
    end

    @testset "Gaps" begin

        file = joinpath(pwd(), "data", "gaps.txt")
        cl = hobohmI(read(file, Raw), 62)
        gu, gi = pairwisegapfraction(file, Raw)
        ncl = nclusters(cl)

        @test gu[1, 1] ≈ 0.0
        @test gi[1, 1] ≈ 0.0
        @test gu[1, 2] ≈ 100.0 * getweight(cl, 10)/ncl
        @test gi[1, 2] ≈ 0.0
        @test gu[10, 9] ≈ 100.0 * (ncl - getweight(cl, 1))/ncl
        @test gi[10, 9] ≈ 100.0 * (ncl - getweight(cl, 1) - getweight(cl, 2))/ncl
        @test gu[10, 10] ≈ 100.0 * (ncl - getweight(cl, 1))/ncl
        @test gu[10, 10] ≈ 100.0 * (ncl - getweight(cl, 1))/ncl
    end
end


text:
using Plots
p = plot(rand(5), rand(5))
display(p)

max_issues_repo_name: jondea/juliab

max_forks_repo_name: jondea/juliab


text:
function NbodyEnergy(u,Gm)
     N = length(Gm)
     zerouel = zero(eltype(u))
     T = zerouel
     U = zerouel
     for i in 1:N
        qi = u[2,:,i]
        vi = u[1,:,i]
        Gmi = Gm[i]
        T += Gmi*(vi[1]*vi[1]+vi[2]*vi[2]+vi[3]*vi[3])
        for j in (i+1):N
           qj = u[2,:,j]  # qj = u[2,:,j]
           Gmj = Gm[j]
           qij = qi - qj
           U -= Gmi*Gmj/norm(qij)
        end
     end
    1/2*T + U
end

function NbodyEnergy2(u,Gm)
#
#   special implementantion: considering:
#              u[1,:,:]: position
#              u[2,:,:]: velocity
#
     N = length(Gm)
     zerouel = zero(eltype(u))
     T = zerouel
     U = zerouel
     for i in 1:N
        qi = u[1,:,i]    # qi = u[2,:,i]
        vi = u[2,:,i]    # vi = u[1,:,i]
        Gmi = Gm[i]
        T += Gmi*(vi[1]*vi[1]+vi[2]*vi[2]+vi[3]*vi[3])
        for j in (i+1):N
           qj = u[1,:,j]  # qj = u[2,:,j]
           Gmj = Gm[j]
           qij = qi - qj
           U -= Gmi*Gmj/norm(qij)
        end
     end
    1/2*T + U
end

# OdeProblem

function NbodyODE!(du,u,Gm,t)
     N = length(Gm)
     du[1,:,:] .= zero(eltype(u))
     for i in 1:N
        qi = u[2,:,i]
        Gmi = Gm[i]
        du[2,:,i] = u[1,:,i]
        for j in (i+1):N
           qj = u[2,:,j]
           Gmj = Gm[j]
           qij = qi - qj
           auxij = (qij[1]*qij[1]+qij[2]*qij[2]+qij[3]*qij[3])^(-3/2)
           du[1,:,i] -= Gmj*auxij*qij
           du[1,:,j] += Gmi*auxij*qij
        end
     end
    return
end

# DynamicalODEProblem

function NbodyODEv!(dv,q,v,Gm,t)
#
#    dotv
#
     N = length(Gm)
     dv[:,:] .= zero(eltype(q))
     for i in 1:N
        qi = q[:,i]
        Gmi = Gm[i]
        for j in (i+1):N
           qj = q[:,j]
           Gmj = Gm[j]
           qij = qi - qj
           auxij = (qij[1]*qij[1]+qij[2]*qij[2]+qij[3]*qij[3])^(-3/2)
           dv[:,i] -= Gmj*auxij*qij
           dv[:,j] += Gmi*auxij*qij
        end
     end

    return

end

function NbodyODEq!(dq,q,v,Gm,t)
#
#    dotq
#
     N = length(Gm)
     for i in 1:N
        dq[:,i] = v[:,i]
     end

    return

end

## 2ndOrderProblem

function NbodyODE2nd!(ddu,du,u,Gm,t)
     N = length(Gm)
     ddu[:,:] .= zero(eltype(u))
     for i in 1:N
        qi = u[:,i]
        Gmi = Gm[i]
        for j in (i+1):N
           qj = u[:,j]
           Gmj = Gm[j]
           qij = qi - qj
           auxij = (qij[1]*qij[1]+qij[2]*qij[2]+qij[3]*qij[3])^(-3/2)
           ddu[:,i] -= Gmj*auxij*qij
           ddu[:,j] += Gmi*auxij*qij
        end
     end

    return

end


text:
include("../../../good_plots/iron_line_heatmap.jl")

vars = 5.0:85.0

hmp = iron_line_heatmap(
                        tolerance=1e-12,
                        dtmax=5,
                        size_multiplier=4,
                        nbins=200,
                        vars=vars
                        )
png(hmp, "iron_heatmap.png")


text:
using CIJ, Test

@testset "EC" begin
    @testset "Constructors" begin
        let
            # Check on size
            @test_throws ArgumentError EC([1 2; 3 4])
            @test_throws ArgumentError EC(rand(5,6))
            # No constructor
            @test_throws ArgumentError EC((1.0, 2.0))
            # Default precision
            @test EC(zeros(6,6)) isa EC{Float64}
            @test EC(zeros(Float32, 6, 6)) isa EC{Float32}
            @test EC(fill(1, 6, 6)) isa EC{float(Int)}
            @test EC(Array{Any,2}(rand(6,6))) isa EC{CIJ.DEFAULT_FLOAT}
            # Symmetrisation using upper half
            x = [10i + j for i in 1:6, j in 1:6]
            c = EC(x)
            @test c == [i<j ? 10i+j : 10j+i for i in 1:6, j in 1:6]
            # Warning for asymmetric matrices
            @test_logs (:warn, "input matrix not symmetrical: taking upper half"
                ) EC((x = rand(6,6); x[1,2] = -x[1,2]; x), warn=true)
            # Tuple and array construction
            c′ = EC(((10i + j for j in 1:6 for i in 1:6)...,))
            @test c == c′
            # Construction with another EC uses same underlying MMatrix
            c2 = EC(c)
            @test c == c2
            c[1,1] = rand()
            @test c2[1,1] == c[1,1]
            # Symmetrisation on setindex!
            c[1,2] = rand()
            @test c[2,1] == c[1,2]
            @test zero(EC) == zeros(CIJ.DEFAULT_FLOAT, 6, 6)
            @test zero(EC{CIJ.DEFAULT_FLOAT}) == zero(EC)
            @test zero(c) == zero(EC)
        end
    end

    @testset "Properties" begin
        let c = EC(rand(6,6))
            @test size(c) == (6, 6)
            @test length(c) == 36
        end
    end
end


text:
using DrWatson, Test

@testset "DrWatson" begin
    @testset "Naming" begin include("naming_tests.jl"); end
    @testset "Parse savename" begin include("parse_tests.jl"); end
    @testset "Project Setup" begin include("project_tests.jl"); end
    @testset "Saving tools" begin include("stools_tests.jl"); end
    @testset "Produce or Save" begin include("savefiles_tests.jl"); end
    @testset "Collect Results" begin include("update_results_tests.jl"); end
    @testset "Parameter Customization" begin include("customize_savename.jl"); end
end


text:
## Utils
function triangle_list_from_marker(flat_triangle_list)
    triangle_list = Array{Array{Int64,1},1}()

    for i in 1:3:length(flat_triangle_list)
        push!(triangle_list, flat_triangle_list[i:i+2])
    end

    return triangle_list
end

function flat_vertices(vertices::Array{Float64,2}, vertices_map::Array{Int64,1})
    vert_size = size(vertices)
    flat_vertices_vector = Vector{Cdouble}(undef, vert_size[1]*vert_size[2])

    # for vert_id in vertices_map
    for vert_id=1:vert_size[1]
        flat_vertices_vector[(vert_id*2)-1]=vertices[vert_id]
        flat_vertices_vector[(vert_id*2)]=vertices[vert_id+vert_size[1]]
    end

    return flat_vertices_vector
end

function flat_edges(edges::Array{Int64,2})
    edge_size = size(edges)
    flat_edges_vector = Vector{Cint}(undef, edge_size[1]*edge_size[2])

    for edge_id=1:edge_size[1]
        flat_edges_vector[(edge_id*2)-1]=edges[edge_id]
        flat_edges_vector[(edge_id*2)]=edges[edge_id+edge_size[1]]
    end

    return flat_edges_vector
end


text:
"""
Provides functionality for computing the Cuthill-McKee graph permutation on a sparse matrix in CSC format.
"""
module CuthillMcKee

import SparseArrays: SparseMatrixCSC

export rcmpermute, symrcm

"""
	 rcmpermute(A)
Computes and applies the reverse Cuthill-McKee permutation to the structurally-symmetric sparse matrix `A`. No checks are performed to ensure that the matrix is indeed structurally symmetric.
"""
function rcmpermute(A::SparseMatrixCSC)
	adj = coladjacency(A)
	sdeg = sortperm(coldegrees(A))
	p = symrcm(adj, sdeg)
	return A[p,p]
end

"""
	 symrcm(A,[ rev, inv, warnunconnected])
Computes the Cuthill-McKee permutation of the structurally-symmetric sparse matrix `A`. The default is to compute the reverse permutation, i.e. `rev` defaults to `true`. Optionally, when `inv` is set to `true`, the permutation is inverted. Per default the warning displaying the number of unconnected regions in the graph is disabled, to enable, set `warnunconnected` to `true`. No checks are performed to ensure that the matrix is indeed structurally symmetric.
"""
function symrcm(A::SparseMatrixCSC, args...)
	adj = coladjacency(A)
	sdeg = sortperm(coldegrees(A))
	return symrcm(adj, sdeg, args...)
end

"""
	 symrcm(adjac, sdegr,[ rev, inv, warnunconnected])
Computes the Cuthill-McKee permutation of the graph given as the adjacency list `adjac`. The adjacency list is required to be of type `Vector{Vector{Int}}` such that `adjac[i]` is a vector of the adjacent vertex numbers to vertex `i` in order of ascending vertex degree. The additional required argument `sdegr` is an array of vertex numbers in order of ascending vertex degree. The default is to compute the reverse permutation, i.e. `rev` defaults to `true`. Optionally, when `inv` is set to `true`, the permutation is inverted. Per default the warning displaying the number of unconnected regions in the graph is disabled, to enable, set `warnunconnected` to `true`.
"""
function symrcm(adjac::Array{Array{Int64,1}}, sdegr::Array{Int64,1}, rev=true, invert=false, warnunconnected=false)
	# Based on the description of the RCM algorithm by Ciprian Zavoianu.
	# See: http://ciprian-zavoianu.blogspot.com/2009/01/project-bandwidth-reduction.html
	nverts = length(adjac)
	checkbounds(sdegr,nverts)
	sd = copy(sdegr)
	R = Array{Int}(undef, 0)
	Q = Array{Int}(undef, 0)
	P = 0
	C = 0
	unconnected = 0
	inserted = falses(nverts)
	while !isempty(sd)
		found = false
		while !isempty(sd)
			P = popfirst!(sd)
			if !inserted[P]
				push!(R, P)
				inserted[P] = true
				found = true
				break
			end
		end
		if found
			empty!(Q)
			append!(Q, adjac[P])
			while !isempty(Q)
				while !isempty(Q)
					C = popfirst!(Q)
					if !inserted[C]
						push!(R, C)
						inserted[C] = true
						append!(Q, adjac[C])
						break
					end
				end
			end
		end
		if length(R) == nverts
			break
		else
			unconnected += 1
		end
	end
	if unconnected > 0 && warnunconnected
		warn("There are $(unconnected) unconnected regions.")
	end
	if rev
		reverse!(R)
	end
	if invert
		L = zeros(Int64,nverts)
		for i = 1:nverts
			L[R[i]] = i
		end
		return L
	else
		return R
	end
end

function coldegrees(A::SparseMatrixCSC)
	# assumes A is structurally symmetric
	cptr = A.colptr
	ncols = length(cptr)-1
	degr = zeros(Int,ncols)
	@inbounds for j = 1:ncols
		degr[j] = coldeg(A, j)
	end
	return degr
end

function coldeg(A::SparseMatrixCSC, j::Int)
	# assumes A is structurally symmetric
	cptr = A.colptr
	return cptr[j+1]-cptr[j]
end

function coladjacency(A::SparseMatrixCSC)
	# assumes A is structurally symmetric
	cptr = A.colptr
	rval = A.rowval
	ncols = length(cptr)-1
	adjac = Vector{Vector{Int}}(undef, ncols)
	sbyf = let A = A
		j->coldeg(A, j)
	end
	for j = 1:ncols
		strt = cptr[j]
		jdeg = coldeg(A, j)
		jadj = Vector{Int}(undef, jdeg)
		for i = 1:jdeg
			jadj[i] = rval[strt+i-1]
		end
		sort!(jadj, by=sbyf)
		adjac[j] = jadj
	end
	return adjac
end

end


text:
using Documenter, Indicators, Plots, GLMakie

makedocs(
    modules = [Indicators],
    sitename = "Indicators",
    authors="Jacob Amos",
    pages=["Home"=>"index.md",
           "Conventional" => ["Moving Averages" => "ma.md",
                              "Momentum Indicators" => "mom.md",
                              "Volatility Indicators" => "vol.md"],
           "Exotic" => ["Regressions"=>"reg.md",
                        "Trendlines" => "trendy.md",
                        "Chaos" => "chaos.md",
                        "Patterns" => "patterns.md"]],
    format = Documenter.HTML(),
    doctest=false,
    clean=true,
)

deploydocs(deps=Deps.pip("mkdocs", "python-markdown-math"),
           repo="github.com/dysonance/Indicators.jl.git",
           devbranch="master",
           devurl="dev",
           versions=["stable" => "v^", "v#.#", "dev" => "dev"])


text:
"""
    Bitbucket

Implements read-only access
"""
module Bitbucket

import ..GitForge: endpoint, into, postprocessor, ForgeAPIError
import ..GitForge.GitHub: ismemberorcollaborator
import Base.@kwdef, Base.Meta.quot
import StructTypes: Struct, UnorderedStruct, construct, constructfrom, StructType

using ..GitForge
using ..GitForge:
    @json,
    AStr,
    DoSomething,
    Endpoint,
    Forge,
    JSON,
    OnRateLimit,
    RateLimiter,
    HEADERS,
    ORL_THROW,
    @not_implemented
using ..GitForge.GitHub: NoToken, JWT, AbstractToken
using Dates, TimeZones, HTTP, JSON2, UUIDs, StructTypes
using Dates: Date, DateTime

export BitbucketAPI, Token, JWT

# using a const here in order to export it
const Token = GitForge.GitHub.Token
const DEFAULT_URL = "https://api.bitbucket.org/2.0"

auth_headers(::NoToken) = []
auth_headers(t::Token) = ["Authorization" => "Basic $(t.token)"]
auth_headers(t::JWT) = ["Authorization" => "Bearer $(t.token)"]

function JSON2.read(io::IO, ::Type{UUID}; kwargs...)
    str = JSON2.read(io, String)
    UUID(str[2:end-1])
end
JSON2.write(io::IO, uuid::UUID; kwargs...) = JSON2.write(io, "{$uuid}"; kwargs...)

"""
    BitbucketAPI(;
        token::AbstractToken=NoToken(),
        url::$AStr="$DEFAULT_URL",
        has_rate_limits::Bool=true,
        on_rate_limit::OnRateLimit=ORL_THROW,
    )

Create a Bitbucket API client.

## Keywords
- `token::AbstractToken=NoToken()`: Authorization token (or lack thereof).
- `url::$AStr="$DEFAULT_URL"`: Base URL of the target Bitbucket instance.
- `has_rate_limits::Bool=true`: Whether or not the Bitbucket server has rate limits.
- `on_rate_limit::OnRateLimit=ORL_THROW`: Behaviour on exceeded rate limits.
- `workspace::AbstractString=""`: slug for chosen workspace
"""
mutable struct BitbucketAPI <: Forge
    token::AbstractToken
    url::AbstractString
    hasrl::Bool
    orl::OnRateLimit
    rl_general::RateLimiter
    rl_search::RateLimiter
    workspace::AbstractString

    function BitbucketAPI(;
        token::AbstractToken=NoToken(),
        url::AStr=DEFAULT_URL,
        has_rate_limits::Bool=false,
        on_rate_limit::OnRateLimit=ORL_THROW,
        workspace::AbstractString="",
    )
        return new(token, url, has_rate_limits, on_rate_limit, RateLimiter(), RateLimiter(), workspace)
    end
end

GitForge.base_url(b::BitbucketAPI) = b.url
GitForge.request_headers(b::BitbucketAPI, ::Function) = [HEADERS; auth_headers(b.token)]
GitForge.postprocessor(::BitbucketAPI, ::Function) = JSON_Struct()
GitForge.has_rate_limits(b::BitbucketAPI, ::Function) = b.hasrl
GitForge.rate_limit_check(b::BitbucketAPI, ::Function) = GitForge.rate_limit_check(b.rl_general)
GitForge.on_rate_limit(b::BitbucketAPI, ::Function) = b.orl
GitForge.rate_limit_wait(b::BitbucketAPI, ::Function) = GitForge.rate_limit_wait(b.rl_general)
GitForge.rate_limit_period(b::BitbucketAPI, ::Function) = GitForge.rate_limit_period(b.rl_general)
GitForge.rate_limit_update!(b::BitbucketAPI, ::Function, r::HTTP.Response) =
    GitForge.rate_limit_update!(b.rl_general, r)

cvt_date(str::AbstractString) =
    ZonedDateTime(replace(str, r"(\.[0-9]{3})([0-9]*)\+" => s"\1+")).utc_datetime

function zdate(kw::NamedTuple, prop::Symbol)
    !haskey(kw, prop) && return (;)
    (; prop =>
        ZonedDateTime(replace(kw[prop], r"(\.[0-9]{3})([0-9]*)\+" => s"\1+")).utc_datetime)
end

struct BBStructClosure{T}
    obj::T
end

constructfield(::Type{FT}, v) where {FT} = constructfrom(FT, v)

# capture type Nothing fields here so the methods below won't apply to them
constructfield(::Type{Nothing}, v) = constructfrom(Nothing, v)

constructfield(::Type{FT}, str::AbstractString) where {FT <: Union{UUID, Nothing}} = UUID(str[2:end-1])

constructfield(::Type{FT}, v::AbstractString) where {FT <: Union{Date, Nothing}} =
    Date(ZonedDateTime(replace(v, r"(\.[0-9]{3})([0-9]*)\+" => s"\1+")))

constructfield(::Type{FT}, v::AbstractString) where {FT <: Union{DateTime, Nothing}} =
    DateTime(ZonedDateTime(replace(v, r"(\.[0-9]{3})([0-9]*)\+" => s"\1+")))

@inline function (f::BBStructClosure{T})(_i, nm, ::Type{FT}) where {T, FT}
    hasfield(T, nm) ? constructfield(FT, getfield(f.obj, nm)) : nothing
end

function collect_extras(type::Type, kw::NamedTuple)
    names = fieldnames(type)
    main = []
    extras = []
    for k in keys(kw)
        push!(k in names ? main : extras, k => kw[k])
    end
    (; main..., _extras = (; extras...))
end

sortkw(kw)=(;sort([pairs(kw)...], by=p->p[1])...)

struct JSON_Struct <: GitForge.PostProcessor
    f::Function
end
JSON_Struct() = JSON_Struct(identity)

function GitForge.postprocess(p::JSON_Struct, r::HTTP.Response, ::Type{T}) where T
    data = JSON2.read(IOBuffer(r.body))
    p.f(constructfrom(T, JSON2.read(IOBuffer(r.body))))
end

macro json_struct(type)
    quote
        StructTypes.StructType(::Type{$type}) = UnorderedStruct()

        StructTypes.constructfrom(::Struct, ::Type{$type}, ::Union{StructTypes.Mutable, StructTypes.Struct}, obj) =
            StructTypes.construct(BBStructClosure(obj), $type)

        function StructTypes.constructfrom(t::Struct, ::Type{$type}, kw::NamedTuple)
            kw = collect_extras($type, kw)
            constructfrom(UnorderedStruct(), $type, StructType(NamedTuple), kw)
        end
    end
end

include("pagination.jl")
include("users.jl")
include("repositories.jl")
include("pull_requests.jl")
include("workspaces.jl")
include("commits.jl")
include("branches.jl")
include("tags.jl")
#include("comments.jl")

end


text:
### Envs
## Basic environments
include("basics/basics.jl")

## Multicopters
include("multicopters/multicopters.jl")

## Control Allocators
include("allocators/allocators.jl")

## Controllers
include("controllers/controllers.jl")

## Integrated Envs
include("integrated_environments/integrated_environments.jl")


text:
# Using the BBP formula
#
# The Bailey–Borwein–Plouffe formula (BBP formula) is a formula for π

#=
Computing Bailey–Borwein–Plouffe formula for Pi
to ascertain the performance of FastRational{BigInt}
relative to Rational{BigInt} by BPP iterations

FastRational{BigInt}s perform better below n=328
Rational{BigInt}s     perform better above n=330

As BPP iterations increase, numerator and denominator grow.
After some digit count in numerator aor denominator, ~10_000,
The extra time spent in divgcd reductions becomes less than
the time spent performing arithmetic on very large integers;
so reducing the numerator and denominator begins to win.
=#

using FastRationals
using BenchmarkTools

const BT=BenchmarkTools.DEFAULT_PARAMETERS;
BT.overhead=BenchmarkTools.estimate_overhead();
BT.evals=1; ; BT.time_tolerance = 5.0e-7; BT.samples = 100;

const big1 = BigInt(1)
const big2 = BigInt(2)
const big4 = BigInt(4)
const big5 = BigInt(5)
const big6 = BigInt(6)
const big8 = BigInt(8)
const big16 = BigInt(16)

function bpp(::Type{T}, n) where {T}
    result = zero(T)
    for k = 0:n
       eightk = big8 * k
       cur = T(big4,eightk+1) -
             T(big2,eightk+4) -
             T(big1,eightk+5) -
             T(big1,eightk+6)
       cur = T(big1, big16^k) * cur
       result = result + cur
    end
    return result
end

systemqtime = @belapsed bpp(Rational{BigInt},   15);
fastqtime = @belapsed bpp(FastRational{BigInt}, 15);
bpp15 = round(systemqtime/fastqtime, digits=1)

systemqtime = @belapsed bpp(Rational{BigInt},   125);
fastqtime = @belapsed bpp(FastRational{BigInt}, 125);
bpp125 = round(systemqtime/fastqtime, digits=1)

systemqtime = @belapsed bpp(Rational{BigInt},   250);
fastqtime = @belapsed bpp(FastRational{BigInt}, 250);
bpp250 = round(systemqtime/fastqtime, digits=1)

# relspeeds meet at n=328

systemqtime = @belapsed bpp(Rational{BigInt},   500);
fastqtime = @belapsed bpp(FastRational{BigInt}, 500);
bpp500 = round(systemqtime/fastqtime, digits=1)

systemqtime = @belapsed bpp(Rational{BigInt},   1000);
fastqtime = @belapsed bpp(FastRational{BigInt}, 1000);
bpp1000 = round(systemqtime/fastqtime, digits=1)

systemqtime = @belapsed bpp(Rational{BigInt},   2000);
fastqtime = @belapsed bpp(FastRational{BigInt}, 2000);
bpp2000 = round(systemqtime/fastqtime, digits=1)

systemqtime = @belapsed bpp(Rational{BigInt},   3000);
fastqtime = @belapsed bpp(FastRational{BigInt}, 3000);
bpp3000 = round(systemqtime/fastqtime, digits=1)

relspeeds = (bpp15=bpp15, bpp125=bpp125, bpp250=bpp250,
             bpp500=bpp500, bpp1000=bpp1000, bpp2000=bpp2000,
             bpp3000=bpp3000)

relspeeds = (bpp15 = 0.3, bpp125 = 0.3, bpp250 = 0.7, bpp500 = 2.4, bpp1000 = 8.1, bpp2000 = 20.9, bpp3000 = 36.0)

xs = [15,125,250,500,1000,2000,3000];
ys = [relspeeds...,];

using GR
plot(xs, ys, size=(500,500))


text:
using Combinatorics

numbers = open("expense_report.txt") do f
    parse.(Int, readlines(f))
end

for pair in combinations(numbers, 2)
    if sum(pair) == 2020
        println(prod(pair))
    end
end

for pair in combinations(numbers, 3)
    if sum(pair) == 2020
        println(prod(pair))
    end
end


text:
using Compat

################################################################################
# Export definitions supported by our export script                            #
################################################################################
export Tox, Tox_Options

# typealias
export TOX_CHAT_CHANGE
export TOX_USER_STATUS
export TOX_MESSAGE_TYPE
export TOX_PROXY_TYPE
export TOX_SAVEDATA_TYPE
export TOX_ERR_OPTIONS_NEW
export TOX_ERR_NEW
export TOX_ERR_BOOTSTRAP
export TOX_CONNECTION
export TOX_ERR_SET_INFO
export TOX_ERR_FRIEND_ADD
export TOX_ERR_FRIEND_DELETE
export TOX_ERR_FRIEND_BY_PUBLIC_KEY
export TOX_ERR_FRIEND_GET_PUBLIC_KEY
export TOX_ERR_FRIEND_GET_LAST_ONLINE
export TOX_ERR_FRIEND_QUERY
export TOX_ERR_SET_TYPING
export TOX_ERR_FRIEND_SEND_MESSAGE
export TOX_FILE_KIND
export TOX_FILE_CONTROL
export TOX_ERR_FILE_CONTROL
export TOX_ERR_FILE_SEEK
export TOX_ERR_FILE_GET
export TOX_ERR_FILE_SEND
export TOX_ERR_FILE_SEND_CHUNK
export TOX_ERR_FRIEND_CUSTOM_PACKET
export TOX_ERR_GET_PORT

export tox_self_connection_status_cb
export tox_friend_name_cb
export tox_friend_status_message_cb
export tox_friend_status_cb
export tox_friend_connection_status_cb
export tox_friend_typing_cb
export tox_friend_read_receipt_cb
export tox_friend_request_cb
export tox_friend_message_cb
export tox_file_recv_control_cb
export tox_file_chunk_request_cb
export tox_file_recv_cb
export tox_file_recv_chunk_cb
export tox_friend_lossy_packet_cb
export tox_friend_lossless_packet_cb

# consts
export TOX_CHAT_CHANGE_PEER_ADD
export TOX_CHAT_CHANGE_PEER_DEL
export TOX_CHAT_CHANGE_PEER_NAME
export TOX_VERSION_MAJOR
export TOX_VERSION_MINOR
export TOX_VERSION_PATCH
export TOX_PUBLIC_KEY_SIZE
export TOX_SECRET_KEY_SIZE
export TOX_MAX_NAME_LENGTH
export TOX_MAX_STATUS_MESSAGE_LENGTH
export TOX_MAX_FRIEND_REQUEST_LENGTH
export TOX_MAX_MESSAGE_LENGTH
export TOX_MAX_CUSTOM_PACKET_SIZE
export TOX_HASH_LENGTH
export TOX_FILE_ID_LENGTH
export TOX_MAX_FILENAME_LENGTH
export TOX_USER_STATUS_NONE
export TOX_USER_STATUS_AWAY
export TOX_USER_STATUS_BUSY
export TOX_MESSAGE_TYPE_NORMAL
export TOX_MESSAGE_TYPE_ACTION
export TOX_PROXY_TYPE_NONE
export TOX_PROXY_TYPE_HTTP
export TOX_PROXY_TYPE_SOCKS
export TOX_SAVEDATA_TYPE_NONE
export TOX_SAVEDATA_TYPE_TOX_SAVE
export TOX_SAVEDATA_TYPE_SECRET_KEY
export TOX_ERR_OPTIONS_NEW_OK
export TOX_ERR_OPTIONS_NEW_MALLOC
export TOX_ERR_NEW_OK
export TOX_ERR_NEW_NULL
export TOX_ERR_NEW_MALLOC
export TOX_ERR_NEW_PORT_ALLOC
export TOX_ERR_NEW_PROXY_BAD_TYPE
export TOX_ERR_NEW_PROXY_BAD_HOST
export TOX_ERR_NEW_PROXY_BAD_PORT
export TOX_ERR_NEW_PROXY_NOT_FOUND
export TOX_ERR_NEW_LOAD_ENCRYPTED
export TOX_ERR_NEW_LOAD_BAD_FORMAT
export TOX_ERR_BOOTSTRAP_OK
export TOX_ERR_BOOTSTRAP_NULL
export TOX_ERR_BOOTSTRAP_BAD_HOST
export TOX_ERR_BOOTSTRAP_BAD_PORT
export TOX_CONNECTION_NONE
export TOX_CONNECTION_TCP
export TOX_CONNECTION_UDP
export TOX_ERR_SET_INFO_OK
export TOX_ERR_SET_INFO_NULL
export TOX_ERR_SET_INFO_TOO_LONG
export TOX_ERR_FRIEND_ADD_OK
export TOX_ERR_FRIEND_ADD_NULL
export TOX_ERR_FRIEND_ADD_TOO_LONG
export TOX_ERR_FRIEND_ADD_NO_MESSAGE
export TOX_ERR_FRIEND_ADD_OWN_KEY
export TOX_ERR_FRIEND_ADD_ALREADY_SENT
export TOX_ERR_FRIEND_ADD_BAD_CHECKSUM
export TOX_ERR_FRIEND_ADD_SET_NEW_NOSPAM
export TOX_ERR_FRIEND_ADD_MALLOC
export TOX_ERR_FRIEND_DELETE_OK
export TOX_ERR_FRIEND_DELETE_FRIEND_NOT_FOUND
export TOX_ERR_FRIEND_BY_PUBLIC_KEY_OK
export TOX_ERR_FRIEND_BY_PUBLIC_KEY_NULL
export TOX_ERR_FRIEND_BY_PUBLIC_KEY_NOT_FOUND
export TOX_ERR_FRIEND_GET_PUBLIC_KEY_OK
export TOX_ERR_FRIEND_GET_PUBLIC_KEY_FRIEND_NOT_FOUND
export TOX_ERR_FRIEND_GET_LAST_ONLINE_OK
export TOX_ERR_FRIEND_GET_LAST_ONLINE_FRIEND_NOT_FOUND
export TOX_ERR_FRIEND_QUERY_OK
export TOX_ERR_FRIEND_QUERY_NULL
export TOX_ERR_FRIEND_QUERY_FRIEND_NOT_FOUND
export TOX_ERR_SET_TYPING_OK
export TOX_ERR_SET_TYPING_FRIEND_NOT_FOUND
export TOX_ERR_FRIEND_SEND_MESSAGE_OK
export TOX_ERR_FRIEND_SEND_MESSAGE_NULL
export TOX_ERR_FRIEND_SEND_MESSAGE_FRIEND_NOT_FOUND
export TOX_ERR_FRIEND_SEND_MESSAGE_FRIEND_NOT_CONNECTED
export TOX_ERR_FRIEND_SEND_MESSAGE_SENDQ
export TOX_ERR_FRIEND_SEND_MESSAGE_TOO_LONG
export TOX_ERR_FRIEND_SEND_MESSAGE_EMPTY
export TOX_FILE_KIND_DATA
export TOX_FILE_KIND_AVATAR
export TOX_FILE_CONTROL_RESUME
export TOX_FILE_CONTROL_PAUSE
export TOX_FILE_CONTROL_CANCEL
export TOX_ERR_FILE_CONTROL_OK
export TOX_ERR_FILE_CONTROL_FRIEND_NOT_FOUND
export TOX_ERR_FILE_CONTROL_FRIEND_NOT_CONNECTED
export TOX_ERR_FILE_CONTROL_NOT_FOUND
export TOX_ERR_FILE_CONTROL_NOT_PAUSED
export TOX_ERR_FILE_CONTROL_DENIED
export TOX_ERR_FILE_CONTROL_ALREADY_PAUSED
export TOX_ERR_FILE_CONTROL_SENDQ
export TOX_ERR_FILE_SEEK_OK
export TOX_ERR_FILE_SEEK_FRIEND_NOT_FOUND
export TOX_ERR_FILE_SEEK_FRIEND_NOT_CONNECTED
export TOX_ERR_FILE_SEEK_NOT_FOUND
export TOX_ERR_FILE_SEEK_DENIED
export TOX_ERR_FILE_SEEK_INVALID_POSITION
export TOX_ERR_FILE_SEEK_SENDQ
export TOX_ERR_FILE_GET_OK
export TOX_ERR_FILE_GET_NULL
export TOX_ERR_FILE_GET_FRIEND_NOT_FOUND
export TOX_ERR_FILE_GET_NOT_FOUND
export TOX_ERR_FILE_SEND_OK
export TOX_ERR_FILE_SEND_NULL
export TOX_ERR_FILE_SEND_FRIEND_NOT_FOUND
export TOX_ERR_FILE_SEND_FRIEND_NOT_CONNECTED
export TOX_ERR_FILE_SEND_NAME_TOO_LONG
export TOX_ERR_FILE_SEND_TOO_MANY
export TOX_ERR_FILE_SEND_CHUNK_OK
export TOX_ERR_FILE_SEND_CHUNK_NULL
export TOX_ERR_FILE_SEND_CHUNK_FRIEND_NOT_FOUND
export TOX_ERR_FILE_SEND_CHUNK_FRIEND_NOT_CONNECTED
export TOX_ERR_FILE_SEND_CHUNK_NOT_FOUND
export TOX_ERR_FILE_SEND_CHUNK_NOT_TRANSFERRING
export TOX_ERR_FILE_SEND_CHUNK_INVALID_LENGTH
export TOX_ERR_FILE_SEND_CHUNK_SENDQ
export TOX_ERR_FILE_SEND_CHUNK_WRONG_POSITION
export TOX_ERR_FRIEND_CUSTOM_PACKET_OK
export TOX_ERR_FRIEND_CUSTOM_PACKET_NULL
export TOX_ERR_FRIEND_CUSTOM_PACKET_FRIEND_NOT_FOUND
export TOX_ERR_FRIEND_CUSTOM_PACKET_FRIEND_NOT_CONNECTED
export TOX_ERR_FRIEND_CUSTOM_PACKET_INVALID
export TOX_ERR_FRIEND_CUSTOM_PACKET_EMPTY
export TOX_ERR_FRIEND_CUSTOM_PACKET_TOO_LONG
export TOX_ERR_FRIEND_CUSTOM_PACKET_SENDQ
export TOX_ERR_GET_PORT_OK
export TOX_ERR_GET_PORT_NOT_BOUND
export TOX_ADDRESS_SIZE

################################################################################
# This is the code generated by Clang (Anonymous enums removed)                #
################################################################################

# begin enum TOX_CHAT_CHANGE
typealias TOX_CHAT_CHANGE Uint32
const TOX_CHAT_CHANGE_PEER_ADD = (UInt32)(0)
const TOX_CHAT_CHANGE_PEER_DEL = (UInt32)(1)
const TOX_CHAT_CHANGE_PEER_NAME = (UInt32)(2)
# end enum TOX_CHAT_CHANGE

const TOX_VERSION_MAJOR = 0
const TOX_VERSION_MINOR = 0
const TOX_VERSION_PATCH = 0

# Skipping MacroDefinition: TOX_VERSION_IS_API_COMPATIBLE ( MAJOR , MINOR , PATCH ) ( TOX_VERSION_MAJOR == MAJOR && ( TOX_VERSION_MINOR > MINOR || ( TOX_VERSION_MINOR == MINOR && TOX_VERSION_PATCH >= PATCH ) ) )
# Skipping MacroDefinition: TOX_VERSION_REQUIRE ( MAJOR , MINOR , PATCH ) typedef char tox_required_version [ TOX_IS_COMPATIBLE ( MAJOR , MINOR , PATCH ) ? 1 : - 1 ]
# Skipping MacroDefinition: TOX_VERSION_IS_ABI_COMPATIBLE ( ) tox_version_is_compatible ( TOX_VERSION_MAJOR , TOX_VERSION_MINOR , TOX_VERSION_PATCH )

const TOX_PUBLIC_KEY_SIZE = 32
const TOX_SECRET_KEY_SIZE = 32

# Skipping MacroDefinition: TOX_ADDRESS_SIZE ( TOX_PUBLIC_KEY_SIZE + sizeof ( uint32_t ) + sizeof ( uint16_t ) )

const TOX_MAX_NAME_LENGTH = 128
const TOX_MAX_STATUS_MESSAGE_LENGTH = 1007
const TOX_MAX_FRIEND_REQUEST_LENGTH = 1016
const TOX_MAX_MESSAGE_LENGTH = 1372
const TOX_MAX_CUSTOM_PACKET_SIZE = 1373
const TOX_HASH_LENGTH = 32
const TOX_FILE_ID_LENGTH = 32
const TOX_MAX_FILENAME_LENGTH = 255

# begin enum TOX_USER_STATUS
typealias TOX_USER_STATUS Uint32
const TOX_USER_STATUS_NONE = (UInt32)(0)
const TOX_USER_STATUS_AWAY = (UInt32)(1)
const TOX_USER_STATUS_BUSY = (UInt32)(2)
# end enum TOX_USER_STATUS

# begin enum TOX_MESSAGE_TYPE
typealias TOX_MESSAGE_TYPE Uint32
const TOX_MESSAGE_TYPE_NORMAL = (UInt32)(0)
const TOX_MESSAGE_TYPE_ACTION = (UInt32)(1)
# end enum TOX_MESSAGE_TYPE

# begin enum TOX_PROXY_TYPE
typealias TOX_PROXY_TYPE Uint32
const TOX_PROXY_TYPE_NONE = (UInt32)(0)
const TOX_PROXY_TYPE_HTTP = (UInt32)(1)
const TOX_PROXY_TYPE_SOCKS5 = (UInt32)(2)
# end enum TOX_PROXY_TYPE

# begin enum TOX_SAVEDATA_TYPE
typealias TOX_SAVEDATA_TYPE Uint32
const TOX_SAVEDATA_TYPE_NONE = (UInt32)(0)
const TOX_SAVEDATA_TYPE_TOX_SAVE = (UInt32)(1)
const TOX_SAVEDATA_TYPE_SECRET_KEY = (UInt32)(2)
# end enum TOX_SAVEDATA_TYPE

# begin enum TOX_ERR_OPTIONS_NEW
typealias TOX_ERR_OPTIONS_NEW Uint32
const TOX_ERR_OPTIONS_NEW_OK = (UInt32)(0)
const TOX_ERR_OPTIONS_NEW_MALLOC = (UInt32)(1)
# end enum TOX_ERR_OPTIONS_NEW

# begin enum TOX_ERR_NEW
typealias TOX_ERR_NEW Uint32
const TOX_ERR_NEW_OK = (UInt32)(0)
const TOX_ERR_NEW_NULL = (UInt32)(1)
const TOX_ERR_NEW_MALLOC = (UInt32)(2)
const TOX_ERR_NEW_PORT_ALLOC = (UInt32)(3)
const TOX_ERR_NEW_PROXY_BAD_TYPE = (UInt32)(4)
const TOX_ERR_NEW_PROXY_BAD_HOST = (UInt32)(5)
const TOX_ERR_NEW_PROXY_BAD_PORT = (UInt32)(6)
const TOX_ERR_NEW_PROXY_NOT_FOUND = (UInt32)(7)
const TOX_ERR_NEW_LOAD_ENCRYPTED = (UInt32)(8)
const TOX_ERR_NEW_LOAD_BAD_FORMAT = (UInt32)(9)
# end enum TOX_ERR_NEW

# begin enum TOX_ERR_BOOTSTRAP
typealias TOX_ERR_BOOTSTRAP Uint32
const TOX_ERR_BOOTSTRAP_OK = (UInt32)(0)
const TOX_ERR_BOOTSTRAP_NULL = (UInt32)(1)
const TOX_ERR_BOOTSTRAP_BAD_HOST = (UInt32)(2)
const TOX_ERR_BOOTSTRAP_BAD_PORT = (UInt32)(3)
# end enum TOX_ERR_BOOTSTRAP

# begin enum TOX_CONNECTION
typealias TOX_CONNECTION Uint32
const TOX_CONNECTION_NONE = (UInt32)(0)
const TOX_CONNECTION_TCP = (UInt32)(1)
const TOX_CONNECTION_UDP = (UInt32)(2)
# end enum TOX_CONNECTION

# begin enum TOX_ERR_SET_INFO
typealias TOX_ERR_SET_INFO Uint32
const TOX_ERR_SET_INFO_OK = (UInt32)(0)
const TOX_ERR_SET_INFO_NULL = (UInt32)(1)
const TOX_ERR_SET_INFO_TOO_LONG = (UInt32)(2)
# end enum TOX_ERR_SET_INFO

# begin enum TOX_ERR_FRIEND_ADD
typealias TOX_ERR_FRIEND_ADD Uint32
const TOX_ERR_FRIEND_ADD_OK = (UInt32)(0)
const TOX_ERR_FRIEND_ADD_NULL = (UInt32)(1)
const TOX_ERR_FRIEND_ADD_TOO_LONG = (UInt32)(2)
const TOX_ERR_FRIEND_ADD_NO_MESSAGE = (UInt32)(3)
const TOX_ERR_FRIEND_ADD_OWN_KEY = (UInt32)(4)
const TOX_ERR_FRIEND_ADD_ALREADY_SENT = (UInt32)(5)
const TOX_ERR_FRIEND_ADD_BAD_CHECKSUM = (UInt32)(6)
const TOX_ERR_FRIEND_ADD_SET_NEW_NOSPAM = (UInt32)(7)
const TOX_ERR_FRIEND_ADD_MALLOC = (UInt32)(8)
# end enum TOX_ERR_FRIEND_ADD

# begin enum TOX_ERR_FRIEND_DELETE
typealias TOX_ERR_FRIEND_DELETE Uint32
const TOX_ERR_FRIEND_DELETE_OK = (UInt32)(0)
const TOX_ERR_FRIEND_DELETE_FRIEND_NOT_FOUND = (UInt32)(1)
# end enum TOX_ERR_FRIEND_DELETE

# begin enum TOX_ERR_FRIEND_BY_PUBLIC_KEY
typealias TOX_ERR_FRIEND_BY_PUBLIC_KEY Uint32
const TOX_ERR_FRIEND_BY_PUBLIC_KEY_OK = (UInt32)(0)
const TOX_ERR_FRIEND_BY_PUBLIC_KEY_NULL = (UInt32)(1)
const TOX_ERR_FRIEND_BY_PUBLIC_KEY_NOT_FOUND = (UInt32)(2)
# end enum TOX_ERR_FRIEND_BY_PUBLIC_KEY

# begin enum TOX_ERR_FRIEND_GET_PUBLIC_KEY
typealias TOX_ERR_FRIEND_GET_PUBLIC_KEY Uint32
const TOX_ERR_FRIEND_GET_PUBLIC_KEY_OK = (UInt32)(0)
const TOX_ERR_FRIEND_GET_PUBLIC_KEY_FRIEND_NOT_FOUND = (UInt32)(1)
# end enum TOX_ERR_FRIEND_GET_PUBLIC_KEY

# begin enum TOX_ERR_FRIEND_GET_LAST_ONLINE
typealias TOX_ERR_FRIEND_GET_LAST_ONLINE Uint32
const TOX_ERR_FRIEND_GET_LAST_ONLINE_OK = (UInt32)(0)
const TOX_ERR_FRIEND_GET_LAST_ONLINE_FRIEND_NOT_FOUND = (UInt32)(1)
# end enum TOX_ERR_FRIEND_GET_LAST_ONLINE

# begin enum TOX_ERR_FRIEND_QUERY
typealias TOX_ERR_FRIEND_QUERY Uint32
const TOX_ERR_FRIEND_QUERY_OK = (UInt32)(0)
const TOX_ERR_FRIEND_QUERY_NULL = (UInt32)(1)
const TOX_ERR_FRIEND_QUERY_FRIEND_NOT_FOUND = (UInt32)(2)
# end enum TOX_ERR_FRIEND_QUERY

# begin enum TOX_ERR_SET_TYPING
typealias TOX_ERR_SET_TYPING Uint32
const TOX_ERR_SET_TYPING_OK = (UInt32)(0)
const TOX_ERR_SET_TYPING_FRIEND_NOT_FOUND = (UInt32)(1)
# end enum TOX_ERR_SET_TYPING

# begin enum TOX_ERR_FRIEND_SEND_MESSAGE
typealias TOX_ERR_FRIEND_SEND_MESSAGE Uint32
const TOX_ERR_FRIEND_SEND_MESSAGE_OK = (UInt32)(0)
const TOX_ERR_FRIEND_SEND_MESSAGE_NULL = (UInt32)(1)
const TOX_ERR_FRIEND_SEND_MESSAGE_FRIEND_NOT_FOUND = (UInt32)(2)
const TOX_ERR_FRIEND_SEND_MESSAGE_FRIEND_NOT_CONNECTED = (UInt32)(3)
const TOX_ERR_FRIEND_SEND_MESSAGE_SENDQ = (UInt32)(4)
const TOX_ERR_FRIEND_SEND_MESSAGE_TOO_LONG = (UInt32)(5)
const TOX_ERR_FRIEND_SEND_MESSAGE_EMPTY = (UInt32)(6)
# end enum TOX_ERR_FRIEND_SEND_MESSAGE

# begin enum TOX_FILE_KIND
typealias TOX_FILE_KIND Uint32
const TOX_FILE_KIND_DATA = (UInt32)(0)
const TOX_FILE_KIND_AVATAR = (UInt32)(1)
# end enum TOX_FILE_KIND

# begin enum TOX_FILE_CONTROL
typealias TOX_FILE_CONTROL Uint32
const TOX_FILE_CONTROL_RESUME = (UInt32)(0)
const TOX_FILE_CONTROL_PAUSE = (UInt32)(1)
const TOX_FILE_CONTROL_CANCEL = (UInt32)(2)
# end enum TOX_FILE_CONTROL

# begin enum TOX_ERR_FILE_CONTROL
typealias TOX_ERR_FILE_CONTROL Uint32
const TOX_ERR_FILE_CONTROL_OK = (UInt32)(0)
const TOX_ERR_FILE_CONTROL_FRIEND_NOT_FOUND = (UInt32)(1)
const TOX_ERR_FILE_CONTROL_FRIEND_NOT_CONNECTED = (UInt32)(2)
const TOX_ERR_FILE_CONTROL_NOT_FOUND = (UInt32)(3)
const TOX_ERR_FILE_CONTROL_NOT_PAUSED = (UInt32)(4)
const TOX_ERR_FILE_CONTROL_DENIED = (UInt32)(5)
const TOX_ERR_FILE_CONTROL_ALREADY_PAUSED = (UInt32)(6)
const TOX_ERR_FILE_CONTROL_SENDQ = (UInt32)(7)
# end enum TOX_ERR_FILE_CONTROL

# begin enum TOX_ERR_FILE_SEEK
typealias TOX_ERR_FILE_SEEK Uint32
const TOX_ERR_FILE_SEEK_OK = (UInt32)(0)
const TOX_ERR_FILE_SEEK_FRIEND_NOT_FOUND = (UInt32)(1)
const TOX_ERR_FILE_SEEK_FRIEND_NOT_CONNECTED = (UInt32)(2)
const TOX_ERR_FILE_SEEK_NOT_FOUND = (UInt32)(3)
const TOX_ERR_FILE_SEEK_DENIED = (UInt32)(4)
const TOX_ERR_FILE_SEEK_INVALID_POSITION = (UInt32)(5)
const TOX_ERR_FILE_SEEK_SENDQ = (UInt32)(6)
# end enum TOX_ERR_FILE_SEEK

# begin enum TOX_ERR_FILE_GET
typealias TOX_ERR_FILE_GET Uint32
const TOX_ERR_FILE_GET_OK = (UInt32)(0)
const TOX_ERR_FILE_GET_NULL = (UInt32)(1)
const TOX_ERR_FILE_GET_FRIEND_NOT_FOUND = (UInt32)(2)
const TOX_ERR_FILE_GET_NOT_FOUND = (UInt32)(3)
# end enum TOX_ERR_FILE_GET

# begin enum TOX_ERR_FILE_SEND
typealias TOX_ERR_FILE_SEND Uint32
const TOX_ERR_FILE_SEND_OK = (UInt32)(0)
const TOX_ERR_FILE_SEND_NULL = (UInt32)(1)
const TOX_ERR_FILE_SEND_FRIEND_NOT_FOUND = (UInt32)(2)
const TOX_ERR_FILE_SEND_FRIEND_NOT_CONNECTED = (UInt32)(3)
const TOX_ERR_FILE_SEND_NAME_TOO_LONG = (UInt32)(4)
const TOX_ERR_FILE_SEND_TOO_MANY = (UInt32)(5)
# end enum TOX_ERR_FILE_SEND

# begin enum TOX_ERR_FILE_SEND_CHUNK
typealias TOX_ERR_FILE_SEND_CHUNK Uint32
const TOX_ERR_FILE_SEND_CHUNK_OK = (UInt32)(0)
const TOX_ERR_FILE_SEND_CHUNK_NULL = (UInt32)(1)
const TOX_ERR_FILE_SEND_CHUNK_FRIEND_NOT_FOUND = (UInt32)(2)
const TOX_ERR_FILE_SEND_CHUNK_FRIEND_NOT_CONNECTED = (UInt32)(3)
const TOX_ERR_FILE_SEND_CHUNK_NOT_FOUND = (UInt32)(4)
const TOX_ERR_FILE_SEND_CHUNK_NOT_TRANSFERRING = (UInt32)(5)
const TOX_ERR_FILE_SEND_CHUNK_INVALID_LENGTH = (UInt32)(6)
const TOX_ERR_FILE_SEND_CHUNK_SENDQ = (UInt32)(7)
const TOX_ERR_FILE_SEND_CHUNK_WRONG_POSITION = (UInt32)(8)
# end enum TOX_ERR_FILE_SEND_CHUNK

# begin enum TOX_ERR_FRIEND_CUSTOM_PACKET
typealias TOX_ERR_FRIEND_CUSTOM_PACKET Uint32
const TOX_ERR_FRIEND_CUSTOM_PACKET_OK = (UInt32)(0)
const TOX_ERR_FRIEND_CUSTOM_PACKET_NULL = (UInt32)(1)
const TOX_ERR_FRIEND_CUSTOM_PACKET_FRIEND_NOT_FOUND = (UInt32)(2)
const TOX_ERR_FRIEND_CUSTOM_PACKET_FRIEND_NOT_CONNECTED = (UInt32)(3)
const TOX_ERR_FRIEND_CUSTOM_PACKET_INVALID = (UInt32)(4)
const TOX_ERR_FRIEND_CUSTOM_PACKET_EMPTY = (UInt32)(5)
const TOX_ERR_FRIEND_CUSTOM_PACKET_TOO_LONG = (UInt32)(6)
const TOX_ERR_FRIEND_CUSTOM_PACKET_SENDQ = (UInt32)(7)
# end enum TOX_ERR_FRIEND_CUSTOM_PACKET

# begin enum TOX_ERR_GET_PORT
typealias TOX_ERR_GET_PORT Uint32
const TOX_ERR_GET_PORT_OK = (UInt32)(0)
const TOX_ERR_GET_PORT_NOT_BOUND = (UInt32)(1)
# end enum TOX_ERR_GET_PORT

immutable Array_128_Uint8
    d1::Uint8
    d2::Uint8
    d3::Uint8
    d4::Uint8
    d5::Uint8
    d6::Uint8
    d7::Uint8
    d8::Uint8
    d9::Uint8
    d10::Uint8
    d11::Uint8
    d12::Uint8
    d13::Uint8
    d14::Uint8
    d15::Uint8
    d16::Uint8
    d17::Uint8
    d18::Uint8
    d19::Uint8
    d20::Uint8
    d21::Uint8
    d22::Uint8
    d23::Uint8
    d24::Uint8
    d25::Uint8
    d26::Uint8
    d27::Uint8
    d28::Uint8
    d29::Uint8
    d30::Uint8
    d31::Uint8
    d32::Uint8
    d33::Uint8
    d34::Uint8
    d35::Uint8
    d36::Uint8
    d37::Uint8
    d38::Uint8
    d39::Uint8
    d40::Uint8
    d41::Uint8
    d42::Uint8
    d43::Uint8
    d44::Uint8
    d45::Uint8
    d46::Uint8
    d47::Uint8
    d48::Uint8
    d49::Uint8
    d50::Uint8
    d51::Uint8
    d52::Uint8
    d53::Uint8
    d54::Uint8
    d55::Uint8
    d56::Uint8
    d57::Uint8
    d58::Uint8
    d59::Uint8
    d60::Uint8
    d61::Uint8
    d62::Uint8
    d63::Uint8
    d64::Uint8
    d65::Uint8
    d66::Uint8
    d67::Uint8
    d68::Uint8
    d69::Uint8
    d70::Uint8
    d71::Uint8
    d72::Uint8
    d73::Uint8
    d74::Uint8
    d75::Uint8
    d76::Uint8
    d77::Uint8
    d78::Uint8
    d79::Uint8
    d80::Uint8
    d81::Uint8
    d82::Uint8
    d83::Uint8
    d84::Uint8
    d85::Uint8
    d86::Uint8
    d87::Uint8
    d88::Uint8
    d89::Uint8
    d90::Uint8
    d91::Uint8
    d92::Uint8
    d93::Uint8
    d94::Uint8
    d95::Uint8
    d96::Uint8
    d97::Uint8
    d98::Uint8
    d99::Uint8
    d100::Uint8
    d101::Uint8
    d102::Uint8
    d103::Uint8
    d104::Uint8
    d105::Uint8
    d106::Uint8
    d107::Uint8
    d108::Uint8
    d109::Uint8
    d110::Uint8
    d111::Uint8
    d112::Uint8
    d113::Uint8
    d114::Uint8
    d115::Uint8
    d116::Uint8
    d117::Uint8
    d118::Uint8
    d119::Uint8
    d120::Uint8
    d121::Uint8
    d122::Uint8
    d123::Uint8
    d124::Uint8
    d125::Uint8
    d126::Uint8
    d127::Uint8
    d128::Uint8
end

zero(::Type{Array_128_Uint8}) = begin  # /home/sbordt/.julia/v0.4/Clang/src/wrap_c.jl, line 264:
        Array_128_Uint8(fill(zero(Uint8),128)...)
    end

################################################################################
# Added by hand                                                                #
################################################################################
immutable Tox end

const TOX_ADDRESS_SIZE = TOX_PUBLIC_KEY_SIZE + sizeof(Uint32) + sizeof(Uint16)

typealias tox_self_connection_status_cb Void
typealias tox_friend_name_cb Void
typealias tox_friend_status_message_cb Void
typealias tox_friend_status_cb Void
typealias tox_friend_connection_status_cb Void
typealias tox_friend_typing_cb Void
typealias tox_friend_read_receipt_cb Void
typealias tox_friend_request_cb Void
typealias tox_friend_message_cb Void
typealias tox_file_recv_control_cb Void
typealias tox_file_chunk_request_cb Void
typealias tox_file_recv_cb Void
typealias tox_file_recv_chunk_cb Void
typealias tox_friend_lossy_packet_cb Void
typealias tox_friend_lossless_packet_cb Void

immutable Tox_Options
    ipv6_enabled::Bool
    udp_enabled::Bool
    proxy_type::TOX_PROXY_TYPE
    proxy_host::Ptr{Cchar}
    proxy_port::Uint16
    start_port::Uint16
    end_port::Uint16
    tcp_port::Uint16
    savedata_type::TOX_SAVEDATA_TYPE
    savedata_data::Ptr{Uint8}
    savedata_length::Csize_t
end

function Tox_Options()
    # some default values that might make sense
    Tox_Options(false, false, 0, C_NULL, 0, 0, 0, 0, TOX_SAVEDATA_TYPE_NONE, C_NULL, 0 )
end


text:
#!/usr/bin/julia

using Flux

W = rand(2, 5)
b = rand(2)

predict(x) = (W * x) .+ b
loss(x, y) = sum((predict(x) .- y).^2)

x, y = rand(5), rand(2) # Dummy data
println("Data: ")
println(x, y)

l = loss(x, y) # ~ 3

θ = Flux.Params([W, b])
println("\nInitial params: ", θ)

grads = gradient(() -> loss(x, y), θ)

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

# η = 0.1 # Learning Rate
opt = Descent(0.1) # Gradient descent with learning rate 0.1

for i in 1:10
  for p in (W, b)
    # Flux.Optimise.update!(p, -η * grads[p]) # Can specify your own optimizer
    Flux.Optimise.update!(opt, p, grads[p])
  end
  println("Loss: ", loss(x, y))
end

println("\nOptimized params: ", θ)
println("Loss: ", loss(x, y))


text:
using GeneralizedDynamicsFromData
using OrderedCollections
using HDF5
using FileIO

experiment_name = "fritzhugh_nagumo_iv_study"
repetitions = 50

weight_decay = 1e-4
η = 1e-1
η_decay_rate = 0.5
η_decay_step = 100
η_limit = 1e-4
noise = 1e-4

loss_weights = [1/2, 1/2]

construct_optimiser() = Flux.Optimiser(WeightDecay(weight_decay),
                                       ExpDecay(η,η_decay_rate,η_decay_step,η_limit),
                                       ADAM())
construct_loss(θ, y, predict) = polar_loss(θ, y, loss_weights, predict)

net_config = OrderedDict([
    :inputs => 2,
    :outputs => 1,
    :neurons => 16,
    :layers => 1,
    :non_lin => tanh,
    :initialization => Flux.glorot_normal
])

initial_values = [ Float64[-2.0, -0.5 ],
                   Float64[-2.0,  0.5 ],
                   Float64[-1.0,  1.5 ],
                   Float64[ 1.0,  0.5 ],
                   Float64[-1.0, -0.5 ],
                   Float64[ 0.5,  1.55] ]

params = [ [0.9, 0.5, 1.2, 1.25],    # monostable
           [0.9, 0.5, 1.0, 1.25] ]   # bistable

for (id, u0) in enumerate(initial_values)

  problem= Dict([:equation => fritzhugh_nagumo,
                  :parameters => params[1],
                  :u0 => u0,
                  :tspan => (0.0f0, 5.0f0),
                  :ts => 0.1,
                  :solver => Tsit5,
                  :optimizer => construct_optimiser,
                  :max_iter => 1000,
                  :loss => construct_loss
                ])

  summary, _ = repeat_experiment(problem,
                                  net_config,
                                  repetitions;
                                  longterm = 20.0,
                                  ε = noise,
                                  progress=false)

  min_losses = convert(Array{Float64,1}, summary["losses"])
  longterm_predictions = cat([sol for sol in summary["longterm_predictions"]]...; dims=3)
  longterm_solution = cat([sol for sol in summary["longterm_solution"]]...; dims=3)

  fn = joinpath("./data/", experiment_name*"_monostable_"*"$(id).h5")

  h5open(fn, "w") do file
      file["parameters"] = params[1]
      file["u0"] = u0
      file["losses"] = min_losses
      file["longterm_predictions"] = longterm_predictions
      file["longterm_solution"] = longterm_solution
  end

end

for (id, u0) in enumerate(initial_values)

  problem= Dict([:equation => fritzhugh_nagumo,
                  :parameters => params[2],
                  :u0 => u0,
                  :tspan => (0.0f0, 5.0f0),
                  :ts => 0.1,
                  :solver => Tsit5,
                  :optimizer => construct_optimiser,
                  :max_iter => 1000,
                  :loss => construct_loss
                ])

  summary, _ = repeat_experiment(problem,
                                  net_config,
                                  repetitions;
                                  longterm = 20.0,
                                  ε = noise,
                                  progress=false)

  min_losses = convert(Array{Float64,1}, summary["losses"])
  longterm_predictions = cat([sol for sol in summary["longterm_predictions"]]...; dims=3)
  longterm_solution = cat([sol for sol in summary["longterm_solution"]]...; dims=3)

  fn = joinpath("./data/", experiment_name*"_bistable_"*"$(id).h5")

  h5open(fn, "w") do file
      file["parameters"] = params[2]
      file["u0"] = u0
      file["longterm_predictions"] = longterm_predictions
      file["longterm_solution"] = longterm_solution
  end

end

max_stars_repo_name: pnieters/GeneralizedDynamicsFromData

max_issues_repo_name: pnieters/GeneralizedDynamicsFromData

max_forks_repo_name: pnieters/GeneralizedDynamicsFromData


text:
ℓ = multivariate_normal(ones(1))
κ = GaussianKineticEnergy(1)
Q = evaluate_ℓ(ℓ, [1.0])

reporters_1 = [
    NoProgressReport(),
    ProgressMeterReport(),
]

reporters_2 = [
    LogProgressReport(),
]

for reporter in deepcopy(vcat(reporters_1, reporters_2))
    results = mcmc_with_warmup(RNG, ℓ, 10_000; reporter = reporter)
end

for reporter in deepcopy(vcat(reporters_1, reporters_2))
    DynamicHMC.report(reporter, "")

    mcmc_reporter_1 = DynamicHMC.make_mcmc_reporter(reporter, 1_000; currently_warmup = true)
    DynamicHMC.report(mcmc_reporter_1, "")
    DynamicHMC.report(mcmc_reporter_1, 1)

    mcmc_reporter_2 = DynamicHMC.make_mcmc_reporter(reporter, 1_000; currently_warmup = false)
    DynamicHMC.report(mcmc_reporter_2, "")
    DynamicHMC.report(mcmc_reporter_2, 1)
end

for reporter in deepcopy(reporters_1)
    DynamicHMC.report(reporter, "")
    DynamicHMC.report(reporter, 1)

    mcmc_reporter_1 = DynamicHMC.make_mcmc_reporter(reporter, 1_000; currently_warmup = true)
    DynamicHMC.report(mcmc_reporter_1, "")
    DynamicHMC.report(mcmc_reporter_1, 1)

    mcmc_reporter_2 = DynamicHMC.make_mcmc_reporter(reporter, 1_000; currently_warmup = false)
    DynamicHMC.report(mcmc_reporter_2, "")
    DynamicHMC.report(mcmc_reporter_2, 1)
end


text:
using DataFrames, GLM

function construct_variables(cause_cnt, data_cnt)
    exact = zeros(Float64, cause_cnt * data_cnt)
    weight = zeros(Float64, cause_cnt * data_cnt)
    prevalence = zeros(Float64, cause_cnt * data_cnt)
    total = zeros(Float64, cause_cnt * data_cnt)
    bc = zeros(Float64, cause_cnt * data_cnt)
    bc2 = zeros(Float64, cause_cnt * data_cnt)
    w = rand(cause_cnt)
    for i in 1:data_cnt
        p = rand(cause_cnt)
        e = exact_burden(w, p)
        assign = ((i - 1) * cause_cnt + 1):(i * cause_cnt)
        exact[assign] = e
        weight[assign] = w
        prevalence[assign] = p
        bc[assign] = (w .* p) / sum(w .* p)
        bc2[assign] = (w .* p) / sum(e)
        total[assign] .= sum(e)
    end
    weight, prevalence, exact, total, bc, bc2
end

w, p, e, t, bc, bc2 = construct_variables(20, 20)
data = DataFrame(E = e, W = w, P = p, BC = bc, BC2 = bc2, T = t)
ols = lm(@formula(E ~ BC + W) , data)
scatter(e, predict(ols))

w1 = 0.3 * rand(10)
p1 = 0.6 * rand(10)
b1 = exact_burden(w1, p1)

w2 = 0.5 * rand(10)
p2 = 0.2 * rand(10)
b2 = exact_burden(w2, p2)

b12 = vcat(b1, b2)
b3 = exact_burden(vcat(w1, w2), vcat(p1, p2))
sum(b3)
1 - (1 - sum(b1)) * (1 - sum(b2))

# This is an oddly successful way to combine the two smaller exact
# calculations into a larger calculation. It's off by a constant.
b1p = b1 * (1 - sum(b2))
b2p = b2 * (1 - sum(b1))
scatter(b3[1:10], b1p, markercolor = :green)
scatter!(b3[11:20], b2p, markercolor = :blue)
xlabel!("Exact Burden")
ylabel!("Estimated Burden")

guess = vcat(b1p, b2p)
guess = guess * sum(b3) / sum(guess)
rel = (guess - b3) ./ b3
# Max relative error is 0.0363.
maximum(abs.(rel))


text:
#===================================================================================================
NSGA-II

[1] Deb, Kalyanmoy, et al. "A fast and elitist multiobjective genetic algorithm: NSGA-II."
IEEE transactions on evolutionary computation 6.2 (2002): 182-197.
===================================================================================================#

"""
    nondominated_sort(fitness::AbstractMatrix{<:Real}, senses::AbstractVector{<:Integer})::Vector{Vector{Int}}

Fast nondominated sorting on page 3 of [1]. The fronts are returned, where the first one is the Pareto
optimal front, and the i-th front contains the indices of individuals in that front.
"""
function nondominated_sort(fitness::AbstractMatrix{<:Real}, senses::AbstractVector{<:Integer})::Vector{Vector{Int}}
    n = size(fitness, 2) # number of individuals
    d_counter = zeros(Int, n)
    dominated = [Int[] for _ = 1:n]
    fronts = Vector{Int}[]
    @inbounds for ip = 1:n, iq = ip+1:n
        p = view(fitness, :, ip)
        q = view(fitness, :, iq)
        if dominate(p, q, senses)
            push!(dominated[ip], iq)
            d_counter[iq] += 1
        elseif dominate(q, p, senses)
            push!(dominated[iq], ip)
            d_counter[ip] += 1
        end
    end
    # the first front
    push!(fronts, Int[])
    f1 = fronts[1]
    @inbounds for i = 1:n
        if d_counter[i] == 0
            push!(f1, i)
        end
    end
    n_in_fronts = length(f1)
    # the other fronts
    current = 1
    while true
        new_front = Int[]
        sizehint!(new_front, n - n_in_fronts)
        for ip in fronts[current]
            for iq in dominated[ip]
                d_counter[iq] -= 1
                if d_counter[iq] == 0
                    push!(new_front, iq)
                end
            end
        end
        if isempty(new_front)
            break
        else
            push!(fronts, new_front)
            n_in_fronts += length(new_front)
        end
        current += 1
    end
    return fronts
end

"""
    get_nondomination_rank(fronts::AbstractVector{<:AbstractVector{<:Integer}})::Vector{Int}

Given the sorted front, determine which front each individual resides (starting from 1).
The value at the i-th position of the returned array is the rank of the i-th individual in the
population.
"""
function get_nondomination_rank(fronts::AbstractVector{<:AbstractVector{<:Integer}})::Vector{Int}
    n = sum(length, fronts)
    ranks = zeros(Int, n)
    for (r, front) in enumerate(fronts)
        ranks[front] .= r
    end
    return ranks
end

"""
    assign_crowding_distance(fitness::AbstractMatrix{<:Real},
        fronts::AbstractVector{<:AbstractVector{<:Integer}})::Vector{Real}

Compute the crowding distance for each individual in the whole population characterized by its
`fitness` and Pareto front in  `fronts`. Refer to the algorithm on the bottom of page 4 of Ref. [1].

A vector of crowding distances `cds` is returned, in which `cds[i]` denotes the crowding distance of
the individual `i`.
"""
function assign_crowding_distance(fitness::AbstractMatrix{<:Real},
    fronts::AbstractVector{<:AbstractVector{<:Integer}})::Vector{Float64}
    m, n = size(fitness)
    cds = zeros(n)  # crowding distance
    fronts = deepcopy(fronts)
    # handle each front separately
    for front in fronts # front contains the index of an individual in the population
        if length(front) <= 2
            cds[front] .= Inf
            continue
        end
        for io = 1:m  # each objective
            sort!(front; by = idx -> fitness[io, idx])
            scale = fitness[io, front[end]] - fitness[io, front[1]]
            # all individuals in this front have the identical `io`-th objective value, ignore it
            if scale > 0
                cds[front[1]] = Inf
                cds[front[end]] = Inf
                for i = 2:length(front)-1
                    cds[front[i]] += (fitness[io, front[i+1]] - fitness[io, front[i-1]]) / scale
                end
            end
        end
    end
    return cds
end

"""
   crowded_compare(i_rank, i_dist, j_rank, j_dist)::Bool

Crowed-comparison operator. Return -1 if we prefer the individual `i`; otherwise, 1 means
the individual `j` is preferable. If both individuals `i` and `j` are equally good, 0 is returned.
"""
function crowded_compare(i_rank, i_dist, j_rank, j_dist)::Int
    if i_rank < j_rank
        return -1
    elseif i_rank > j_rank
        return 1
    end
    # now i and j belong to the same Pareto front: prefer the less crowded one
    if i_dist > j_dist
        return -1
    elseif i_dist < j_dist
        return 1
    end
    return 0
end


text:
struct SO3 <: Manifold end

dim(::Type{SO3}) = 3
embdim(::Type{SO3}) = 9

abstract type ExpMap <: ChartID end
isglobal(::Chart{<:ExpMap,SO3}) = false
struct ExpMap0 <: ExpMap end
struct ExpMapX <: ExpMap end
struct ExpMapY <: ExpMap end
struct ExpMapZ <: ExpMap end

skew(w) = SA[0.   -w[3]  w[2]
             w[3]  0.   -w[1]
            -w[2]  w[1]  0.]
unskew(ŵ) = SA[ŵ[3,2], ŵ[1,3], ŵ[2,1]]

function logSO3(R)
    trR = tr(R)
    if trR < 3
        θ = acos((trR-1)/2)
        # ŵ = (R - R')/(2*sin(θ))
        # w = unskew(ŵ)
        w = SA[R[3,2] - R[2,3], R[1,3] - R[3,1], R[2,1] - R[1,2]]/(2*sin(θ))
        return θ*w
    else
        return @SVector zeros(3)
    end
end

function expSO3(w)
    w == (@SVector zeros(3)) && (@show "expSO3 singularity"; return SMatrix{3,3,eltype(w)}(I))
    ŵ = skew(w)
    wnorm = norm(w)
    SMatrix{3,3,eltype(w)}(I) + ŵ*sin(wnorm)/wnorm + ŵ^2*(1-cos(wnorm))/wnorm^2
end

tr_vec(pe::SVector{9,S}) where S = pe[1] + pe[5] + pe[9]
chart_boundary_test(pe, ::Chart{ExpMap0,SO3}) = tr_vec(pe)
chart_boundary_test(pe, ::Chart{ExpMapX,SO3}) = tr(RotX(π/2)*reshape(pe, Size(3,3)))
chart_boundary_test(pe, ::Chart{ExpMapY,SO3}) = tr(RotY(π/2)*reshape(pe, Size(3,3)))
chart_boundary_test(pe, ::Chart{ExpMapZ,SO3}) = tr(RotZ(π/2)*reshape(pe, Size(3,3)))

chart_to_emb(p0, ::Chart{ExpMap0,SO3}) = reshape(expSO3(p0), Size(9))
chart_to_emb(px, ::Chart{ExpMapX,SO3}) = reshape(RotX(-π/2)*expSO3(px), Size(9))
chart_to_emb(py, ::Chart{ExpMapY,SO3}) = reshape(RotY(-π/2)*expSO3(py), Size(9))
chart_to_emb(pz, ::Chart{ExpMapZ,SO3}) = reshape(RotZ(-π/2)*expSO3(pz), Size(9))

emb_to_chart(pe, ::Chart{ExpMap0,SO3}) = logSO3(reshape(pe, Size(3,3)))
emb_to_chart(pe, ::Chart{ExpMapX,SO3}) = logSO3(RotX(π/2)*reshape(pe, Size(3,3)))
emb_to_chart(pe, ::Chart{ExpMapY,SO3}) = logSO3(RotY(π/2)*reshape(pe, Size(3,3)))
emb_to_chart(pe, ::Chart{ExpMapZ,SO3}) = logSO3(RotZ(π/2)*reshape(pe, Size(3,3)))

function zero_jacobian_component(::Chart{ExpMap0,SO3})
    SA[0 0  0  0 0 1 0 -1 0
       0 0 -1  0 0 0 1  0 0
       0 1  0 -1 0 0 0  0 0.]
end

function chart_to_emb_jacobian(p, C::Chart{ExpMap0,SO3})
    if norm(p - (@SVector zeros(3))) < 1e-12
        return zero_jacobian_component(C)'
    else
        return ForwardDiff.jacobian(p -> chart_to_emb(p, C), p)
    end
end

function emb_to_chart_jacobian(pe, C::Chart{ExpMap0,SO3})
    if (chart_boundary_test(pe, C) - 3) < 1e-12
        return zero_jacobian_component(C)/2
    else
        return ForwardDiff.jacobian(pe -> emb_to_chart(pe, C), pe)
    end
end

chart_choice_coord_rep(::Chart{<:ExpMap,SO3}) = EmbRep()
function choose_chart_emb(::EmbRep, pe, Cn::Chart{<:ExpMap,SO3})
    C0, Cx = Chart{ExpMap0,SO3}(), Chart{ExpMapX,SO3}()
    Cy, Cz = Chart{ExpMapY,SO3}(), Chart{ExpMapZ,SO3}()

    lo_bound, hi_bound = -0.5, 2.5
    val = chart_boundary_test(pe, Cn)
    if val < lo_bound || val > hi_bound
        val0 = chart_boundary_test(pe, C0)
        valx = chart_boundary_test(pe, Cx)
        valy = chart_boundary_test(pe, Cy)
        valz = chart_boundary_test(pe, Cz)

        val0 = val0 > hi_bound ? -1. : val0
        valx = valx > hi_bound ? -1. : valx
        valy = valy > hi_bound ? -1. : valy
        valz = valz > hi_bound ? -1. : valz

        val_max, ind_max = findmax(SA[val0, valx, valy, valz])
        return (C0, Cx, Cy, Cz)[ind_max]
    else
        return Cn
    end
end

function apply_rotation_emb(xme, R)
    Rtarget = reshape(xme, Size(3,3))
    reshape(R*Rtarget, Size(9))
end

function apply_inverse_rotation_jacobian_emb(xme)
    R = reshape(xme, Size(3,3))
    Rinv = R'
    ForwardDiff.jacobian(xme -> apply_rotation_emb(xme, Rinv), xme)
end

function apply_rotation_chart(xm, R, Ci::Chart{<:ExpMap,SO3})
    xme = chart_to_emb(xm, Ci)
    xme_out = apply_rotation_emb(xme, R)
    xm_out = emb_to_chart(xme_out, Ci)
end

function apply_rotation_chart_emb(xm, R, Ci::Chart{<:ExpMap,SO3})
    xme = chart_to_emb(xm, Ci)
    xme_out = apply_rotation_emb(xme, R)
end

function apply_inverse_rotation_jacobian_chart(xm, Ci::Chart{<:ExpMap,SO3})
    xme = chart_to_emb(xm, Ci)
    R = reshape(xme, Size(3,3))
    Rinv = R'
    xme_out = reshape(SMatrix{3,3,Float64}(I), Size(9))
    J1 = ForwardDiff.jacobian(xm -> apply_rotation_chart_emb(xm, Rinv, Ci), xm)
    J2 = emb_to_chart_jacobian(xme_out, Ci)
    J2*J1
end

function default_metric(pn, task::Task{<:TaskMap{M,SO3}}, CN::Chart{<:ExpMap,SO3}) where M
    C0 = Chart{ExpMap0,SO3}()
    JL = apply_inverse_rotation_jacobian_chart(pn, CN)
    Jϕ = chart_transition_jacobian(pn, CN, C0)
    J = Jϕ*JL
    J'*SMatrix{3,3,Float64}(I)*J
end

default_weight_metric(pn, wn, task::Task{<:TaskMap{M,SO3}}, CN::Chart{<:ExpMap,SO3}) where M =
    default_metric(pn, task, CN)


text:
@testset "TopListArchive" begin
    @testset "TopListIndividual" begin
        # test equality
        @test isequal(BlackBoxOptim.TopListIndividual([3.0, 2.0], 2.0, 0), BlackBoxOptim.TopListIndividual([3.0, 2.0], 2.0, 0))
        @test (BlackBoxOptim.TopListIndividual([3.0, 2.0], 2.0, 0) == BlackBoxOptim.TopListIndividual([3.0, 2.0], 2.0, 0))
        @test (BlackBoxOptim.TopListIndividual([3.0, 2.0], 2.0, 0) != BlackBoxOptim.TopListIndividual([3.0, 2.0], 2.0, 0)) == false
        @test (BlackBoxOptim.TopListIndividual([3.0, 2.0], 1.0, 0) != BlackBoxOptim.TopListIndividual([3.0, 2.0], 2.0, 0))
        @test (BlackBoxOptim.TopListIndividual([1.0, 2.0], 2.0, 0) != BlackBoxOptim.TopListIndividual([3.0, 2.0], 2.0, 0))
    end

    @testset "Constructing a small archive and adding to it" begin

        a = TopListArchive(MinimizingFitnessScheme, 1, 3)

        @test capacity(a) == 3
        @test length(a)   == 0
        @test isnan(best_fitness(a))

        BlackBoxOptim.add_candidate!(a, 1.0, [0.0])
        @test capacity(a)         == 3
        @test length(a)           == 1
        @test best_fitness(a)     == 1.0
        @test best_candidate(a)   == [0.0]
        @test last_top_fitness(a) == 1.0
        @test delta_fitness(a)    == Inf

        BlackBoxOptim.add_candidate!(a, 2.0, [1.0])
        @test capacity(a)       == 3
        @test length(a)         == 2
        @test best_fitness(a)   == 1.0
        @test best_candidate(a) == [0.0]
        @test last_top_fitness(a) == 2.0
        @test delta_fitness(a)    == Inf

        BlackBoxOptim.add_candidate!(a, 0.5, [2.0])
        @test capacity(a)         == 3
        @test length(a)           == 3
        @test best_fitness(a)   == 0.5
        @test best_candidate(a) == [2.0]
        @test last_top_fitness(a) == 2.0
        @test delta_fitness(a)    == 0.5

        BlackBoxOptim.add_candidate!(a, 0.8, [4.0])
        @test capacity(a)         == 3
        @test length(a)           == 3
        @test best_fitness(a)   == 0.5
        @test best_candidate(a) == [2.0]
        @test last_top_fitness(a) == 1.0
        @test delta_fitness(a)    == 0.5

        expected = ((0.8 - 0.5) / ((1 - 0.05)^(-2/1) - 1))
        @test width_of_confidence_interval(a, 0.05) ≈ expected
        @test fitness_improvement_potential(a, 0.05) ≈ (expected / 0.50)

        BlackBoxOptim.add_candidate!(a, 0.4, [1.9])
        @test capacity(a)       == 3
        @test length(a)         == 3
        @test best_fitness(a)   == 0.4
        @test best_candidate(a) == [1.9]
        @test last_top_fitness(a) == 0.8
        @test isapprox(delta_fitness(a), 0.1)

        expected = ((0.5 - 0.4) / ((1 - 0.05)^(-2/1) - 1))
        @test width_of_confidence_interval(a, 0.05) ≈ expected
        @test fitness_improvement_potential(a, 0.05) ≈ (expected / 0.40)

        # identical candidate is not inserted
        BlackBoxOptim.add_candidate!(a, 0.5, [2.0])
        @test capacity(a)       == 3
        @test length(a)         == 3
        @test last_top_fitness(a)   == 0.8

        # different candidates with the same fitness are inserted
        BlackBoxOptim.add_candidate!(a, 0.5, [3.0])
        @test length(a)             == 3
        @test last_top_fitness(a)   == 0.5
    end

    @testset "magnitude_class for positive fitness values" begin

        @test BlackBoxOptim.magnitude_class(1.0)      == (1.0, 0.0)
        @test BlackBoxOptim.magnitude_class(10.0)     == (1.0, 1.0)
        @test BlackBoxOptim.magnitude_class(1e2)      == (1.0, 2.0)
        @test BlackBoxOptim.magnitude_class(1e102)    == (1.0, 102.0)
        @test BlackBoxOptim.magnitude_class(1e-1)     == (1.0, -1.0)
        @test BlackBoxOptim.magnitude_class(1e-2)     == (1.0, -2.0)
        @test BlackBoxOptim.magnitude_class(1e-9)     == (1.0, -9.0)
        @test BlackBoxOptim.magnitude_class(1e-12)    == (1.0, -12.0)

        @test BlackBoxOptim.magnitude_class(2.0)      == (1.0, 0.3)
        @test BlackBoxOptim.magnitude_class(20.0)     == (1.0, 1.3)
        @test BlackBoxOptim.magnitude_class(200.0)    == (1.0, 2.3)

        @test BlackBoxOptim.magnitude_class(5.0)      == (1.0, 0.6)
        @test BlackBoxOptim.magnitude_class(50.0)     == (1.0, 1.6)
        @test BlackBoxOptim.magnitude_class(500.0)    == (1.0, 2.6)

    end

    @testset "magnitude_class for negative fitness values" begin

        @test BlackBoxOptim.magnitude_class(-1.0)      == (-1.0, 0.0)
        @test BlackBoxOptim.magnitude_class(-10.0)     == (-1.0, 1.0)
        @test BlackBoxOptim.magnitude_class(-1e2)      == (-1.0, 2.0)
        @test BlackBoxOptim.magnitude_class(-1e102)    == (-1.0, 102.0)
        @test BlackBoxOptim.magnitude_class(-1e-1)     == (-1.0, -1.0)
        @test BlackBoxOptim.magnitude_class(-1e-2)     == (-1.0, -2.0)
        @test BlackBoxOptim.magnitude_class(-1e-9)     == (-1.0, -9.0)
        @test BlackBoxOptim.magnitude_class(-1e-12)    == (-1.0, -12.0)

        @test BlackBoxOptim.magnitude_class(-2.0)      == (-1.0, 0.3)
        @test BlackBoxOptim.magnitude_class(-20.0)     == (-1.0, 1.3)
        @test BlackBoxOptim.magnitude_class(-200.0)    == (-1.0, 2.3)

        @test BlackBoxOptim.magnitude_class(-5.0)      == (-1.0, 0.6)
        @test BlackBoxOptim.magnitude_class(-50.0)     == (-1.0, 1.6)
        @test BlackBoxOptim.magnitude_class(-500.0)    == (-1.0, 2.6)

    end

    @testset "archive copies the individuals" begin
        a = TopListArchive(MinimizingFitnessScheme, 2, 3)

        indiv = [0.0, 2.0]
        BlackBoxOptim.add_candidate!(a, 1.0, indiv)
        # equal but not identical
        @test (best_candidate(a) == indiv)
        @test (best_candidate(a) === indiv) == false

        # modify the vector
        indiv[1] = 5.0
        # test that archived version is still the same
        @test best_candidate(a) == [0.0, 2.0]
    end

    @testset "for maximizing fitness" begin
        a = TopListArchive(MaximizingFitnessScheme, 1, 3)

        BlackBoxOptim.add_candidate!(a, 1.0, [0.0])
        @test best_fitness(a)     == 1.0
        @test best_candidate(a)   == [0.0]
        @test last_top_fitness(a) == 1.0
        @test delta_fitness(a)    == Inf

        BlackBoxOptim.add_candidate!(a, 2.0, [1.0])
        @test best_fitness(a)   == 2.0
        @test best_candidate(a) == [1.0]
        @test last_top_fitness(a) == 1.0
        @test delta_fitness(a)    == 1.0

        BlackBoxOptim.add_candidate!(a, 0.5, [2.0])
        @test best_fitness(a)   == 2.0
        @test best_candidate(a) == [1.0]
        @test last_top_fitness(a) == 0.5
        @test delta_fitness(a)    == 1.0

        BlackBoxOptim.add_candidate!(a, 1.5, [4.0])
        @test best_fitness(a)   == 2.0
        @test best_candidate(a) == [1.0]
        @test last_top_fitness(a) == 1.0
        @test delta_fitness(a)    == 1.0

        BlackBoxOptim.add_candidate!(a, 2.5, [5.0])
        @test best_fitness(a)   == 2.5
        @test best_candidate(a) == [5.0]
        @test last_top_fitness(a) == 1.5
        @test delta_fitness(a)    == 0.5
    end
end


text:
import Distances
import Base.Test
import StatsBase.entropy

function entropy{T<:Real}(D::Dict{Int64,Array{T,2}})
    E = Dict{Int64,Float64}()
    Et = 0
    S = 0
    Ec = 0
    N = length(D)
    for (k,v) in D
        n = size(v,2) #get the number of trials
        ee = entropy(v)
        E[k] = ee
        Et += n*ee
        S += n
        Ec -= log2(n)
    end
    Ec += N*log2(S)
    Ec /= S
    Et /= S
    return Et,Ec,E
end

function entropy{T<:Real}(x::Array{T,1})
    n = length(x)
    d = zeros(T,n)
    dd = 0.0
    dm = Inf
    for i=1:n
        dm = Inf
        for j=1:n
            if i==j
                continue
            end
            dd = x[i]-x[j]
            dd = sqrt(dd*dd)
            dm = dd < dm ? dd : dm
        end
        d[i] = dm
    end
    H = (1/n)*sum(log2(d)) + log2(2*(n-1)) + eulergamma/log(2)
    return H
end

function entropy{T<:Real}(X::Array{T,2})
    #find the minimum distance between pairs of columns
    ndims = size(X,1)
    N = size(X,2)
    d = Distances.pairwise(Distances.Euclidean(),X)
    #find the minimum distance for each pair
    d[d.==0] = Inf
    md = minimum(d,2)
    S = spherical_volume(ndims)
    H = (ndims/N)*sum(log2(md)) + log2(S*(N-1)/ndims) + eulergamma/log(2)
    return H
end

function spherical_volume(r::Real)
    return r*pi^(r/2)/gamma(r/2+1)
end


text:
#Use XR' = G' = [G1 G2 G3...] to reduce columns of A in
# MXA' + *X* =F
# here G is a vector of Funs

cont_reduce_dofs!{T<:Fun,NT<:Number}( A::AbstractArray{NT},M::Operator,G::Vector{T},F::Fun )=cont_reduce_dofs!( A,M,G,ProductFun(F))
function cont_reduce_dofs!{T<:Fun,NT<:Number}( A::AbstractArray{NT},M::Operator,G::Vector{T},F::ProductFun )
        # first multiply to get MXR' = M*G' = [M*G1 M*G2 ...]
        # then kill the row by subtracting
        # MXR'[:,k]*A'[k,:]  from MXA'
        # i.e., subtacting A[:,k]*R[k,:] from A
        # and M*G'[:,k]*A'[k,:] from F
        # i.e. M*G[k]*A[:,k]' from

    pad!(F,:,max(size(F,2),size(A,1)))

    for k = 1:length(G)
        MG = M*G[k]         # coefficients in the range space of M
        for j=1:min(length(F.coefficients),size(A,1))
            axpy!(-A[j,k],MG,F.coefficients[j]) # equivalent to X+=a*Y
        end
    end

    F
end

# G is  ∞ x K array
# A is ∞ x K list of opcols
# M is ∞ x ∞ operator
# used by kron
cont_reduce_dofs!{NT<:Number,T<:Number}( A::AbstractArray{NT},M::AbstractArray{T},G::Array,F::Fun )=cont_reduce_dofs!( A,M,G,ProductFun(F) )
function cont_reduce_dofs!{NT<:Number,T<:Number}( A::AbstractArray{NT},M::AbstractArray{T},G::Array,F::ProductFun )
    MGA=M*pad(G,size(M,1),size(G,2))*full(A).'
    pad!(F,:,max(size(F,2),size(MGA,2)))
    for j=1:size(MGA,2)
        axpy!(-1,MGA[:,j],F.coefficients[j])
    end
    F
end

cont_reduce_dofs!{T<:Fun,NT<:Number,MT<:Number}( A::AbstractArray{NT},M::AbstractArray{MT},G::Vector{T},F::Fun )=cont_reduce_dofs!(A,M,G,ProductFun(F))
function cont_reduce_dofs!{T<:Fun,NT<:Number,MT<:Number}( A::AbstractArray{NT},M::AbstractArray{MT},G::Vector{T},F::ProductFun )
        # first multiply to get MXR' = M*G' = [M*G1 M*G2 ...]
        # then kill the row by subtracting
        # MXR'[:,k]*A'[k,:]  from MXA'
        # i.e., subtacting A[:,k]*R[k,:] from A
        # and M*G'[:,k]*A'[k,:] from F
        # i.e. M*G[k]*A[:,k]' from

    for k = 1:length(G)
        MG = M*pad(G[k].coefficients,size(M,2))         # coefficients in the range space of M
        for j=1:length(F.coefficients)
            axpy!(-A[j,k],MG,F.coefficients[j]) # equivalent to X+=a*Y
        end
    end

    F
end

function cont_reduce_dofs!{T<:Fun}(S::OperatorSchur,L::Operator,M::Operator,G::Vector{T},F::ProductFun)
    F=cont_reduce_dofs!(S.Lcols,L,G,F)
    cont_reduce_dofs!(S.Mcols,M,G,F)
end

function cont_reduce_dofs!{T<:Fun}(S::OperatorSchur,L::Operator,M::Operator,G::Vector{T},F::Array)
    for Fk in F
        cont_reduce_dofs!(S,L,M,G,Fk)
    end
    F
end

function cont_reduce_dofs!{T<:Fun}(S::OperatorSchur,L::Operator,M::Operator,G::Matrix{T},F::Matrix)
    @assert size(G,2)==size(F,2)
    @assert size(F,1)==1

    for k=1:size(G,2)
        cont_reduce_dofs!(S,L,M,G[:,k],F[1,k])
    end
    F
end

cont_reduce_dofs!{M<:AbstractArray}(Ax::Vector{M},Ay::Vector,G,F::Fun)=cont_reduce_dofs!(Ax,Ay,G,ProductFun(F))
function cont_reduce_dofs!{M<:AbstractArray}(Ax::Vector{M},Ay::Vector,G,F::ProductFun)
    @assert length(Ax)==length(Ay)
    for k=1:length(Ax)
        F=cont_reduce_dofs!(Ax[k],Ay[k],G,F)
    end

    F
end

regularize_bcs(S::OperatorSchur,Gy) = length(Gy)==0?Gy:S.bcQ*Gy

# Solve Bx*Y=Gx and P*Y*R' + S*Y*T' = F
# where R and T are upper triangular

##TODO: Support complex in boundary conditions
cont_constrained_lyapuptriang{OSS,T}(OS::PDEOperatorSchur{OSS,T},Gx,F::ProductFun;kwds...)=cont_constrained_lyapuptriang(promote_type(T,eltype(F)),OS,Gx,F;kwds...)
#cont_constrained_lyapuptriang{N}(::Type{N},OS::PDEOperatorSchur,Gx,F::Array)=cont_constrained_lyapuptriang(N,OS,Gx,F,100000)

cont_constrained_lyap{OSS<:DiagonalOperatorSchur}(OS::PDEOperatorSchur{OSS},Gxin,Gyin,F::Fun;kwds...)=cont_constrained_lyap(OS,Gxin,Gyin,ProductFun(F);kwds...)
function cont_constrained_lyap{OSS<:DiagonalOperatorSchur}(OS::PDEOperatorSchur{OSS},Gxin,Gyin,F::ProductFun;kwds...)
    if isempty(OS.Bx)
        return cont_constrained_lyap_nobcs(OS,F;kwds...)
    end
    n = size(OS.S,1)
    F=pad(F,size(F,1),n)
    Gx=pad(coefficients(Gxin).',:,n)

    TYP=promote_type(eltype(OS),eltype(F))
    Y=Array(Fun{typeof(domainspace(OS,1)),TYP},n)

    for k=1:n
        op=OS.Rdiags[k]
        rhs=Any[Gx[:,k]...;F.coefficients[k]]
        Y[k]=chop!(linsolve(op,rhs;kwds...),eps())
    end

    Y
end

function cont_constrained_lyap_nobcs{OSS<:DiagonalOperatorSchur}(OS::PDEOperatorSchur{OSS},F::ProductFun;kwds...)
    n = size(OS.S,1)
    F=pad(F,size(F,1),n)

    TYP=promote_type(eltype(OS),eltype(F))
    Y=Array(Fun{typeof(domainspace(OS,1)),TYP},n)

    for k=1:n
        op=OS.Rdiags[k]
        rhs=F.coefficients[k]
        Y[k]=chop!(linsolve(op,rhs;kwds...),eps())
    end

    Y
end

function cont_constrained_lyapuptriang{N,OSS<:OperatorSchur}(::Type{N},OS::PDEOperatorSchur{OSS},Gx,F::ProductFun;kwds...)
    n = min(size(OS.S.T,2),max(size(F,2),size(Gx,2)))

    if !isempty(Gx)
        Gx=pad(Gx,size(Gx,1),n)
    end

    rs=rangespace(OS)
    @assert space(F)==rs

    Y=Array(Fun{typeof(domainspace(OS,1)),N},n) # the solution
    PY=Array(Fun{typeof(rs[1]),N},n) # the first x-op times the solution
    SY=Array(Fun{typeof(rs[1]),N},n) # the second x-op times the solution

    k=n
    m=n  # max length

    rhs=Array(Any,size(Gx,1)+1)
    TT=isempty(OS.Bx)?eltype(OS):promote_type(eltype(OS),mapreduce(eltype,promote_type,OS.Bx))

    blkops=Array(Operator{TT},2length(OS.Bx)+2,2)
    if !isempty(OS.Bx)
        blkops[1:2length(OS.Bx),:]=blkdiag(OS.Bx,OS.Bx)
    end
    blkrhs=Array(Any,2size(Gx,1)+2)

    while k≥1
        if k==1 || (OS.S.T[k,k-1] == 0 && OS.S.R[k,k-1] == 0        )  # triangular setting
            rhs[end] = k≤length(F.coefficients)?F.coefficients[k]:zeros(rs[1])

            if k < n
                for j=k+1:n
                    axpy!(-OS.S.R[k,j],PY[j],rhs[end]) # equivalent to X+=a*Y
                    axpy!(-OS.S.T[k,j],SY[j],rhs[end])
                end
            end

            for j=1:size(Gx,1)
                rhs[j]=Gx[j,k]
            end

            ops=OS.Rdiags[k]
            chop!(rhs[end],eps())
            Y[k]=chop!(linsolve(ops,rhs;kwds...),eps())

            if k > 1
                PY[k]=OS.Lx*Y[k];SY[k]=OS.Mx*Y[k]
            end

            k-=1
        else # quasitriangular
            rhs1 = k-1≤length(F.coefficients)?F.coefficients[k-1]:zeros(rs[1])
            rhs2 = k≤length(F.coefficients)?F.coefficients[k]:zeros(rs[1])

            if k < n
                for j=k+1:n
                    axpy!(-OS.S.R[k-1,j],PY[j],rhs1)  # +=
                    axpy!(-OS.S.T[k-1,j],SY[j],rhs1)  # +=

                    axpy!(-OS.S.R[k,j],PY[j],rhs2)  # +=
                    axpy!(-OS.S.T[k,j],SY[j],rhs2)  # +=
                end
            end

            blkops[end-1:end,:]=OS.S.R[k-1:k,k-1:k].*OS.Lx+OS.S.T[k-1:k,k-1:k].*OS.Mx

            for j=1:size(Gx,1)
                blkrhs[j]=Gx[j,k-1]
                blkrhs[j+size(Gx,1)]=Gx[j,k]
            end
            blkrhs[end-1]=rhs1;blkrhs[end]=rhs2

            y=vec(linsolve(blkops,blkrhs;kwds...))
            Y[k-1]=chop!(y[1],eps());Y[k]=chop!(y[2],eps())

            PY[k-1]=OS.Lx*Y[k-1]; PY[k]=OS.Lx*Y[k]
            SY[k-1]=OS.Mx*Y[k-1]; SY[k]=OS.Mx*Y[k]

            k-=2
        end
    end

    Y
end

# Solve for multiple RHS
# TODO: Generalize def
function cont_constrained_lyapuptriang{N,OSS<:OperatorSchur,PF<:ProductFun}(::Type{N},
                                                                            OS::PDEOperatorSchur{OSS},
                                                                            Gx,
                                                                            F::Array{PF,2};kwds...)
    n = min(size(OS.S.T,2),max(mapreduce(Fk->size(Fk,2),max,F),mapreduce(Fk->size(Fk,2),max,Gx)))
    rs=rangespace(OS)

    Fm=size(F,2)

    Y=Array(Fun{typeof(domainspace(OS,1)),N},n,Fm) # the solution
    PY=Array(Fun{typeof(rs[1]),N},n,Fm) # the first x-op times the solution
    SY=Array(Fun{typeof(rs[1]),N},n,Fm) # the second x-op times the solution

    k=n
    m=n  # max length

    nbcs=size(first(Gx),1)
    rhs=Array(Any,nbcs+1,Fm)

    TT=isempty(OS.Bx)?eltype(OS):promote_type(eltype(OS),mapreduce(eltype,promote_type,OS.Bx))

     while k≥1
        @assert k==1 || (OS.S.T[k,k-1] == 0 && OS.S.R[k,k-1] == 0        )

        for j=1:Fm
            rhs[1:nbcs,j]=Gx[j][:,k]
            rhs[end,j]=F[j].coefficients[k]
            if k < n
                for kk=k+1:n
                    axpy!(-OS.S.R[k,kk],PY[kk,j],rhs[end,j]) # equivalent to X+=a*Y
                    axpy!(-OS.S.T[k,kk],SY[kk,j],rhs[end,j])
                end
            end
        end
        ops=OS.Rdiags[k]
        Y[k,:]=vec(linsolve(ops,rhs;kwds...))
        for j=1:Fm
            chop!(Y[k,j],eps())
        end

        if k > 1
            PY[k,:]=vec(OS.Lx*Y[k:k,:])
            SY[k,:]=vec(OS.Mx*Y[k:k,:])  # L*Array returns an array space op
        end
        k-=1
    end
    Y
end

## Solve ∞-dimensional lyap equation
#
#   Bx*X=[Gx[1](y);Gx[2](y);...]=Gx'
#   X*By' = [Gy[1](x) Gy[2](x) ...]
#   Lx*X*Ly' + Mx*X*My' = F
#
# by discretizing in y

#TODO: F should be adaptive rather than array
#TODO: Describe precisely the permutation structure of OS so this algorithm actually makes sense
#      to others (and me too!!)

cont_constrained_lyap{OSS<:OperatorSchur}(OS::PDEOperatorSchur{OSS},
                                          Gx::Vector,
                                          Gyin::Vector,
                                          F::Fun;kwds...)=cont_constrained_lyap(OS,Gx,Gyin,ProductFun(F);kwds...)
function cont_constrained_lyap{OSS<:OperatorSchur}(OS::PDEOperatorSchur{OSS},
                                                   Gx::Vector,
                                                   Gyin::Vector,
                                                   F::ProductFun;kwds...)
    Gy=regularize_bcs(OS.S,Gyin)
    F=pad!(F,:,min(size(OS.S.Q,1),size(F,2)))

    F=cont_reduce_dofs!(OS.S,OS.Lx,OS.Mx,Gy,F)

     # Q2 says how to rearrange the columns of F so that the operator is upper triangular
    Q2 = OS.S.Q
    F=ProductFun(Q2[1:length(F.coefficients),:]'*F.coefficients,space(F))

    ny=size(OS.S,2)
    Ky=numbcs(OS.S)

    ## we've discretized, in y, and rhs for Bx is a function of y
    # so we need to discetize it as well
    # and permute columns by P

    if !isempty(Gx)
        # bcP recombines boundary conditions
        Gx=pad(coefficients(Gx).',:,ny)*OS.S.bcP
        # remove unused DOFs and rearrange columns
        Gx=Gx[:,Ky+1:end]*OS.S.Z
    else
        Gx=[]
    end

    Y=cont_constrained_lyapuptriang(OS,Gx,F;kwds...)
    # Y is a Vector{Fun}, so that Y[k][j] corresponds to matrix element M[k,j]
    # This means acting on Y is acting on *columns* of M

    X22=OS.S.Z*Y  #think of it as transpose

    X11=Gy-OS.S.bcs[:,Ky+1:end]*X22
    X=[X11;X22]
    X=OS.S.bcP*X        # this is equivalent to acting on columns by P'
end

function cont_constrained_lyap{OSS<:OperatorSchur,PF<:ProductFun}(OS::PDEOperatorSchur{OSS},
                                                                  Gx::Array,
                                                                  Gyin::Array,
                                                                  F::Array{PF};kwds...)
    @assert size(F,1)==1

    Gy=regularize_bcs(OS.S,Gyin)
    cont_reduce_dofs!(OS.S,OS.Lx,OS.Mx,Gy,F)
    Q2 = OS.S.Q
    for k=1:size(F,2)
        F[1,k]=ProductFun(Q2[1:length(F[1,k].coefficients),:].'*F[1,k].coefficients,space(F[1,k]))
    end
    ny=size(OS.S,2)
    Ky=numbcs(OS.S)

    if !isempty(Gx)
        GxM=Array(Matrix{Float64},size(Gx,2))
        for k=1:size(Gx,2)
            # bcP recombines boundary conditions
            GxM[k]=pad(coefficients(Gx[:,k]).',:,ny)*OS.S.bcP
            # remove unused DOFs and rearrange columns
            GxM[k]=GxM[k][:,Ky+1:end]*OS.S.Z
        end
    else
        GxM=[]
    end
    nx=10000
    Y=cont_constrained_lyapuptriang(Float64,OS,GxM,F;kwds...)
    #TODO: Ky or Kx?
    X=Array(eltype(Y),size(Y,1)+Ky,size(Y,2))
    for k=1:size(Y,2)
        X[Ky+1:end,k]=OS.S.Z*Y[:,k]
        X[1:Ky,k]=Gy[:,k]-OS.S.bcs[:,Ky+1:end]*X[Ky+1:end,k]
        X[:,k]=OS.S.bcP*X[:,k]
    end

    X
end

cont_constrained_lyap{OSS<:OperatorSchur,PF<:Fun}(OS::PDEOperatorSchur{OSS},
                                      Gx::Array,
                                      Gyin::Array,
                                      F::Array{PF};kwds...)=cont_constrained_lyap(OS,Gx,Gyin,map(ProductFun,F);kwds...)


text:
using ClosureCoefficients
using SparseArrays
using StatsBase
using Test

function undir_test1()
    I = [1, 1, 2, 3, 3]
    J = [2, 3, 3, 4, 5]
    A = sparse(I, J, 1, 5, 5)

    global_clcf = (2 + 2 + 2 + 0 + 0) / (4 + 4 + 2 + 3 + 3)
    avg_clcf    = mean([1.0 / 2, 1.0 / 2, 1.0, 0.0, 0.0])
    local_clcfs =      [1.0 / 2, 1.0 / 2, 1.0, 0.0, 0.0]
    wedges      = [4, 4, 2, 3, 3]
    triangles   = [1, 1, 1, 0, 0]

    C = copy(A)
    C[1, 1] = 1
    for B in [max.(A, A'), A, C, max.(C, C')]
        ret = undir_clcfs(B)
        @test ret.global_clcf ≈ global_clcf
        @test ret.avg_clcf    ≈ avg_clcf
        @test ret.local_clcfs ≈ local_clcfs
        @test ret.wedges     == wedges
        @test ret.triangles  == triangles
    end
end

function dir_test1()
    # Simple all bidirectional edges
    I = [1, 1, 2, 2, 3, 3]
    J = [2, 3, 1, 3, 1, 2]
    A = sparse(I, J, 1, 3, 3)

    ret = dir_clcfs(A)
    for (k, v) in ret
        @test v.global_clcf ≈ 1.0
        @test v.avg_clcf    ≈ 1.0
        @test v.local_clcfs ≈ [1.0, 1.0, 1.0]
        @test v.wedges    == [2, 2, 2]
        @test v.triangles == [2, 2, 2]
    end
end

function dir_test2()
    # Feed forward loop
    I = [1, 1, 2]
    J = [2, 3, 3]
    A = sparse(I, J, 1, 3, 3)

    ret = dir_clcfs(A)
    @test ret["oo_o"].global_clcf ≈ 1.0
    @test ret["oo_o"].avg_clcf    ≈ 1.0
    @test ret["oo_o"].local_clcfs ≈ [1.0, 0.0, 0.0]
    @test ret["oo_o"].wedges    == [1, 0, 0]
    @test ret["oo_o"].triangles == [1, 0, 0]

    @test ret["oo_i"].global_clcf ≈ 0.0
    @test ret["oo_i"].avg_clcf    ≈ 0.0
    @test ret["oo_i"].local_clcfs ≈ [0.0, 0.0, 0.0]
    @test ret["oo_i"].wedges    == [1, 0, 0]
    @test ret["oo_i"].triangles == [0, 0, 0]

    @test ret["oi_o"].global_clcf ≈ 0.5
    @test ret["oi_o"].avg_clcf    ≈ 0.5
    @test ret["oi_o"].local_clcfs ≈ [1.0, 0.0, 0.0]
    @test ret["oi_o"].wedges    == [1, 1, 0]
    @test ret["oi_o"].triangles == [1, 0, 0]

    @test ret["oi_i"].global_clcf ≈ 0.5
    @test ret["oi_i"].avg_clcf    ≈ 0.5
    @test ret["oi_i"].local_clcfs ≈ [0.0, 1.0, 0.0]
    @test ret["oi_i"].wedges    == [1, 1, 0]
    @test ret["oi_i"].triangles == [0, 1, 0]

    @test ret["io_o"].global_clcf ≈ 0.5
    @test ret["io_o"].avg_clcf    ≈ 0.5
    @test ret["io_o"].local_clcfs ≈ [0.0, 1.0, 0.0]
    @test ret["io_o"].wedges    == [0, 1, 1]
    @test ret["io_o"].triangles == [0, 1, 0]

    @test ret["io_i"].global_clcf ≈ 0.5
    @test ret["io_i"].avg_clcf    ≈ 0.5
    @test ret["io_i"].local_clcfs ≈ [0.0, 0.0, 1.0]
    @test ret["io_i"].wedges    == [0, 1, 1]
    @test ret["io_i"].triangles == [0, 0, 1]

    @test ret["ii_o"].global_clcf ≈ 0.0
    @test ret["ii_o"].avg_clcf    ≈ 0.0
    @test ret["ii_o"].local_clcfs ≈ [0.0, 0.0, 0.0]
    @test ret["ii_o"].wedges    == [0, 0, 1]
    @test ret["ii_o"].triangles == [0, 0, 0]

    @test ret["ii_i"].global_clcf ≈ 1.0
    @test ret["ii_i"].avg_clcf    ≈ 1.0
    @test ret["ii_i"].local_clcfs ≈ [0.0, 0.0, 1.0]
    @test ret["ii_i"].wedges    == [0, 0, 1]
    @test ret["ii_i"].triangles == [0, 0, 1]
end

function dir_test3()
    # Cycle
    I = [1, 2, 3]
    J = [2, 3, 1]
    A = sparse(I, J, 1, 3, 3)

    ret = dir_clcfs(A)
    @test ret["oo_o"].global_clcf ≈ 0.0
    @test ret["oo_o"].avg_clcf    ≈ 0.0
    @test ret["oo_o"].local_clcfs ≈ [0.0, 0.0, 0.0]
    @test ret["oo_o"].wedges    == [1, 1, 1]
    @test ret["oo_o"].triangles == [0, 0, 0]

    @test ret["oo_i"].global_clcf ≈ 1.0
    @test ret["oo_i"].avg_clcf    ≈ 1.0
    @test ret["oo_i"].local_clcfs ≈ [1.0, 1.0, 1.0]
    @test ret["oo_i"].wedges    == [1, 1, 1]
    @test ret["oo_i"].triangles == [1, 1, 1]

    @test ret["oi_o"].global_clcf ≈ 0.0
    @test ret["oi_o"].avg_clcf    ≈ 0.0
    @test ret["oi_o"].local_clcfs ≈ [0.0, 0.0, 0.0]
    @test ret["oi_o"].wedges    == [0, 0, 0]
    @test ret["oi_o"].triangles == [0, 0, 0]

    @test ret["oi_i"].global_clcf ≈ 0.0
    @test ret["oi_i"].avg_clcf    ≈ 0.0
    @test ret["oi_i"].local_clcfs ≈ [0.0, 0.0, 0.0]
    @test ret["oi_i"].wedges    == [0, 0, 0]
    @test ret["oi_i"].triangles == [0, 0, 0]

    @test ret["io_o"].global_clcf ≈ 0.0
    @test ret["io_o"].avg_clcf    ≈ 0.0
    @test ret["io_o"].local_clcfs ≈ [0.0, 0.0, 0.0]
    @test ret["io_o"].wedges    == [0, 0, 0]
    @test ret["io_o"].triangles == [0, 0, 0]

    @test ret["io_i"].global_clcf ≈ 0.0
    @test ret["io_i"].avg_clcf    ≈ 0.0
    @test ret["io_i"].local_clcfs ≈ [0.0, 0.0, 0.0]
    @test ret["io_i"].wedges    == [0, 0, 0]
    @test ret["io_i"].triangles == [0, 0, 0]

    @test ret["ii_o"].global_clcf ≈ 1.0
    @test ret["ii_o"].avg_clcf    ≈ 1.0
    @test ret["ii_o"].local_clcfs ≈ [1.0, 1.0, 1.0]
    @test ret["ii_o"].wedges    == [1, 1, 1]
    @test ret["ii_o"].triangles == [1, 1, 1]

    @test ret["ii_i"].global_clcf ≈ 0.0
    @test ret["ii_i"].avg_clcf    ≈ 0.0
    @test ret["ii_i"].local_clcfs ≈ [0.0, 0.0, 0.0]
    @test ret["ii_i"].wedges    == [1, 1, 1]
    @test ret["ii_i"].triangles == [0, 0, 0]
end

function data_test()
    A1 = load_example_data("FW-Florida.txt", symm=false)
    dc = dir_clcfs(A1)
    # TODO(arb): check these numbers more precisely with data from Hao
    @test ≈(dc["ii_i"].avg_clcf, 0.21, atol=0.03)
    @test ≈(dc["oo_i"].avg_clcf, 0.05, atol=0.03)
    A2 = load_example_data("arxiv-AstroPh.txt", symm=true)
    uc = undir_clcfs(A2)
    @test ≈(uc.avg_clcf, 0.250, atol=0.001)
    @test ≈(uc.global_clcf, 0.318, atol=0.001)
end

function run_all()
    undir_test1()
    dir_test1()
    dir_test2()
    dir_test3()
    data_test()
end
run_all()


text:
module day22

p1, p2 = (x -> (parse.(Int, x[2:26]), parse.(Int, x[29:53])))(readlines("day22_input.txt"))

function playround(a, b)
    ta, tb = popfirst!.([a, b])
    ta > tb ? append!(a, [ta, tb]) : append!(b, [tb, ta])
end

function winningscore(a, b)
    while minimum(length.([a, b])) > 0
        playround(a, b)
    end
    sum(prod.(enumerate(reverse(length(a) == 0 ? b : a))))
end

println("Part 1: score = $(winningscore(copy(p1), copy(p2)))")

# Part 2
function recround(a, b)
    ta, tb = popfirst!.([a, b])
    if length(a) < ta || length(b) < tb
        ta > tb ? append!(a, [ta, tb]) : append!(b, [tb, ta])
        return a, b
    end

    awins = recgame(copy(a[1:ta]), copy(b[1:tb]))
    awins ? append!(a, [ta, tb]) : append!(b, [tb, ta])
    a, b
end

# Return true if a wins, false if b wins
function recgame(a, b)
    roundmem = Set()
    while minimum(length.([a, b])) > 0
        if [a, b] ∈ roundmem return true end
        push!(roundmem, [copy(a), copy(b)])
        a, b = recround(a, b)
    end
    length(a) != 0
end

awins = recgame(p1, p2)
println("Part 2: score = $(sum(prod.(enumerate(reverse(awins ? p1 : p2)))))")

end

max_stars_repo_name: ericgreveson/adventofcode

max_issues_repo_name: ericgreveson/adventofcode

max_forks_repo_name: ericgreveson/adventofcode


text:
"""
    CVI.jl

Description:
    Includes all of the CVI/ICVI definitions.
"""

# include("CONN.jl")
include("XB.jl")
include("DB.jl")
include("PS.jl")
include("CH.jl")
include("cSIL.jl")
include("GD53.jl")
include("GD43.jl")
include("WB.jl")
include("rCIP.jl")


text:
"""
$(TYPEDEF)
Kinetic energy specifications. Implements the methods
- `Base.size`
- [`kinetic_energy`](@ref)
- [`calculate_p♯`](@ref)
- [`∇kinetic_energy`](@ref)
- [`rand_p`](@ref)
For all subtypes, it is implicitly assumed that kinetic energy is symmetric in
the momentum `p`,
```julia
kinetic_energy(κ, p, q) == kinetic_energy(κ, .-p, q)
```
When the above is violated, the consequences are undefined.
"""
abstract type KineticEnergy end

"""
$(TYPEDEF)
Euclidean kinetic energies (position independent).
"""
abstract type EuclideanKineticEnergy <: KineticEnergy end

"""
$(TYPEDEF)
Gaussian kinetic energy, with ``K(q,p) = p ∣ q ∼ 1/2 pᵀ⋅M⁻¹⋅p + log|M|`` (without constant),
which is independent of ``q``.
The inverse covariance ``M⁻¹`` is stored.
!!! note
    Making ``M⁻¹`` approximate the posterior variance is a reasonable starting point.
"""
struct GaussianKineticEnergy{P,T,L} <: EuclideanKineticEnergy
    "M⁻¹"
    M⁻¹::Diagonal{T,PtrVector{P,T,L,false}}
    "W such that W*W'=M. Used for generating random draws."
    W::Diagonal{T,PtrVector{P,T,L,false}}
end

# """
# $(SIGNATURES)
# Gaussian kinetic energy with the given inverse covariance matrix `M⁻¹`.
# """
# GaussianKineticEnergy(M⁻¹::AbstractMatrix) = GaussianKineticEnergy(M⁻¹, cholesky(inv(M⁻¹)).L)

"""
$(SIGNATURES)
Gaussian kinetic energy with the given inverse covariance matrix `M⁻¹`.
"""
function GaussianKineticEnergy(sptr::StackPointer, M⁻¹::Diagonal{T,PtrVector{P,T,L,false}}) where {P,T,L}
    sptr, W = PtrVector{P,T,L}(sptr)
    M⁻¹d = M⁻¹.diag
    @avx for l ∈ eachindex(W)
        W[l] = one(T) / sqrt( M⁻¹d[l] )
    end
    sptr, GaussianKineticEnergy(M⁻¹, Diagonal( W ))
end

"""
$(SIGNATURES)
Gaussian kinetic energy with a diagonal inverse covariance matrix `M⁻¹=m⁻¹*I`.
"""
function GaussianKineticEnergy(sptr::StackPointer, ::Static{D}, m⁻¹::T = 1.0) where {D,T}
    sptr, M⁻¹ = PtrVector{D,T}(sptr)
    sptr, W   = PtrVector{D,T}(sptr)
    @fastmath mroot = 1 / sqrt(m⁻¹)
    @avx for d in eachindex(W)
        M⁻¹[d] = m⁻¹
        W[d] = mroot
    end
    # @show M⁻¹, W
    # @show pointer(M⁻¹), pointer(W)
    sptr, GaussianKineticEnergy(Diagonal(M⁻¹), Diagonal(W))
end

function regularized_cov_block_quote(W::Int, T, reps_per_block::Int, stride::Int, mask_last::Bool = false, mask::Unsigned = 0xff)
    # loads from ptr_sample
    # stores in ptr_s² and ptr_invs
    # needs vNinv, mulreg, and addreg to be defined
    reps_per_block -= 1
    size_T = sizeof(T)
    WT = size_T*W
    V = Vec{W,T}
    quote
        $([Expr(:(=), Symbol(:μ_,i), :(SIMDPirates.vload($V, ptr_sample + $(WT*i), $([mask for _ ∈ 1:((i==reps_per_block) & mask_last)]...)))) for i ∈ 0:reps_per_block]...)
        $([Expr(:(=), Symbol(:Σδ_,i), :(SIMDPirates.vbroadcast($V,zero($T)))) for i ∈ 0:reps_per_block]...)
        $([Expr(:(=), Symbol(:Σδ²_,i), :(SIMDPirates.vbroadcast($V,zero($T)))) for i ∈ 0:reps_per_block]...)
        for n ∈ 1:N-1
            $([Expr(:(=), Symbol(:δ_,i), :(SIMDPirates.vsub(SIMDPirates.vload($V, ptr_sample + $(WT*i) + n*$stride*$size_T),$(Symbol(:μ_,i))))) for i ∈ 0:reps_per_block]...)
            $([Expr(:(=), Symbol(:Σδ_,i), :(SIMDPirates.vadd($(Symbol(:δ_,i)),$(Symbol(:Σδ_,i))))) for i ∈ 0:reps_per_block]...)
            $([Expr(:(=), Symbol(:Σδ²_,i), :(SIMDPirates.vmuladd($(Symbol(:δ_,i)),$(Symbol(:δ_,i)),$(Symbol(:Σδ²_,i))))) for i ∈ 0:reps_per_block]...)
        end
        $([Expr(:(=), Symbol(:ΣδΣδ_,i), :(SIMDPirates.vmul($(Symbol(:Σδ_,i)),$(Symbol(:Σδ_,i))))) for i ∈ 0:reps_per_block]...)
        $([Expr(:(=), Symbol(:s²nm1_,i), :(SIMDPirates.vfnmadd($(Symbol(:ΣδΣδ_,i)),vNinv,$(Symbol(:Σδ²_,i))))) for i ∈ 0:reps_per_block]...)
        $([Expr(:(=), Symbol(:regs²_,i), :(SIMDPirates.vmuladd($(Symbol(:s²nm1_,i)), vmulreg, vaddreg))) for i ∈ 0:reps_per_block]...)
        $([Expr(:(=), Symbol(:reginvs_,i), :(SIMDPirates.rsqrt($(Symbol(:regs²_,i))))) for i ∈ 0:reps_per_block]...)
        $([:(vstore!(ptr_s² + $(WT*i), $(Symbol(:regs²_,i)), $([mask for _ ∈ 1:((i==reps_per_block) & mask_last)]...))) for i ∈ 0:reps_per_block]...)
        $([:(vstore!(ptr_invs + $(WT*i), $(Symbol(:reginvs_,i)), $([mask for _ ∈ 1:((i==reps_per_block) & mask_last)]...))) for i ∈ 0:reps_per_block]...)
    end
end

@generated function GaussianKineticEnergy(sp::StackPointer, sample::AbstractMatrix{T}, λ::T, ::Val{D}) where {D,T}
    W, Wshift = VectorizationBase.pick_vector_width_shift(M, T)
    Wm1 = W-1
    M = (D + Wm1) & ~Wm1
    AL = VectorizationBase.align(M*size_T)
    quote
        ptr_s² = pointer(sp, $T)
        regs² = PtrVector{$M,$T}(ptr_s²)
        ptr_invs = ptr_s² + $AL
        invs = PtrVector{$M,$T}(ptr_invs)
        sp + $(2AL), GaussianKineticEnergy!(regs², invs, sample, λ)
    end
end

GaussianKineticEnergy!(K::GaussianKineticEnergy, sample::AbstractMatrix, λ) = GaussianKineticEnergy!(K.M⁻¹.diag, K.W.diag, sample, λ)

@generated function GaussianKineticEnergy!(
    regs²::PtrVector{D,T,L}, invs::PtrVector{D,T,L}, sample::AbstractMatrix{T}, λ::T
# ) where {D,L,T}
) where {D,T,L}
    W, Wshift = VectorizationBase.pick_vector_width_shift(L, T)
    Wm1 = W-1
    # M = (D + Wm1) & ~Wm1
    M = L
    V = Vec{W,T}
    MdW = (M + Wm1) >> Wshift
    Wrem = M & Wm1
    size_T = sizeof(T)
    WT = size_T*W
    need_to_mask = Wrem > 0
    reps_per_block = 4
    blocked_reps, blocked_rem = divrem(MdW, reps_per_block)
    if (MdW % (blocked_reps + 1)) == 0
        blocked_reps, blocked_rem = blocked_reps + 1, 0
        reps_per_block = MdW ÷ blocked_reps
    end
    if (need_to_mask && blocked_rem == 0) || blocked_reps == 1
        blocked_rem += reps_per_block
        blocked_reps -= 1
    end
    # AL = VectorizationBase.align(M*size_T)
    q = quote
        # ptr_s² = pointer(sp, $T)
        # regs² = PtrVector{$M,$T}(ptr_s²)
        # ptr_invs = ptr_s² + $AL
        # invs = PtrVector{$M,$T}(ptr_invs)
        ptr_s² = pointer(regs²)
        ptr_invs = pointer(invs)
        # @show size(sample)
        # display(sample)
        N = size(sample, 2)
        ptr_sample = pointer(sample)
        @fastmath begin
            Ninv = one($T) / N
            mulreg = N / ((N + λ)*(N - 1))
            addreg = $(T(1e-3)) * λ / (N + λ)
        end
        vNinv = SIMDPirates.vbroadcast($V, Ninv)
        vmulreg = SIMDPirates.vbroadcast($V, mulreg)
        vaddreg = SIMDPirates.vbroadcast($V, addreg)
    end
    if blocked_reps > 0
        loop_block = regularized_cov_block_quote(W, T, reps_per_block, M, false)
        body_quote = quote
            $loop_block
            ptr_sample += $WT*$reps_per_block
            ptr_s² += $WT*$reps_per_block
            ptr_invs += $WT*$reps_per_block
        end
        if blocked_reps > 1
            body_quote = quote
                for _ ∈ 1:$blocked_reps
                    $body_quote
                end
            end
        end
        push!(q.args, body_quote)
    end
    if blocked_rem > 0
        push!(q.args, regularized_cov_block_quote(W, T, blocked_rem, M, need_to_mask, VectorizationBase.mask(T,Wrem)))
    end
    # push!(q.args, :(@show regs²))
    # push!(q.args, :(@show invs))
    push!(q.args, :(GaussianKineticEnergy(Diagonal(regs²), Diagonal(invs))))
    # display(q)
    q
end

function Base.show(io::IO, κ::GaussianKineticEnergy{D,T}) where {D,T}
    print(io::IO, "Gaussian kinetic energy ($D,$T), √diag(M⁻¹): $(.√(diag(κ.M⁻¹)))")
end

## NOTE about implementation: the 3 methods are callable without a third argument (`q`)
## because they are defined for Gaussian (Euclidean) kinetic energies.

Base.size(κ::GaussianKineticEnergy{P}) where {P} = (P,P)
# Base.size(κ::GaussianKineticEnergy, args...) = size(κ.M⁻¹, args...)

####
#### Hamiltonian
####

struct Hamiltonian{P,T,N,L}
    "The kinetic energy specification."
    κ::GaussianKineticEnergy{P,T,N}
    """
    The (log) density we are sampling from. Supports the `LogDensityProblem` API.
    Technically, it is the negative of the potential energy.
    """
    ℓ::L
    """
    $(SIGNATURES)
    Construct a Hamiltonian from the log density `ℓ`, and the kinetic energy specification
    `κ`. `ℓ` with a vector are expected to support the `LogDensityProblems` API, with
    gradients.
    """
    function Hamiltonian(κ::GaussianKineticEnergy{P,T,N}, ℓ::L) where {P,T,N,L<:AbstractProbabilityModel{P}}
        # @argcheck capabilities(ℓ) ≥ LogDensityOrder(1)
        # @argcheck dimension(ℓ) == size(κ, 1)
        new{P,T,N,L}(κ, ℓ)
    end
end

Base.show(io::IO, H::Hamiltonian) = print(io, "Hamiltonian with $(H.κ)")

"""
$(TYPEDEF)
A log density evaluated at position `q`. The log densities and gradient are saved, so that
they are not calculated twice for every leapfrog step (both as start- and endpoints).
Because of caching, a `EvaluatedLogDensity` should only be used with a specific Hamiltonian,
preferably constructed with the `evaluate_ℓ` constructor.
In composite types and arguments, `Q` is usually used for this type.
"""
struct EvaluatedLogDensity{P,T,L}
    "Position."
    q::PtrVector{P,T,L,false}
    "ℓ(q). Saved for reuse in sampling."
    ℓq::T
    "∇ℓ(q). Cached for reuse in sampling."
    ∇ℓq::PtrVector{P,T,L,false}
    function EvaluatedLogDensity(q::PtrVector{P,T,L,false}, ℓq::T, ∇ℓq::PtrVector{P,T,L,false}) where {P,T<:Real,L}
#        @argcheck length(q) == length(∇ℓq)
        new{P,T,L}(q, ℓq, ∇ℓq)
    end
end

# general constructors below are necessary to sanitize input from eg Diagnostics, or an
# initial position given as integers, etc

function EvaluatedLogDensity(q::AbstractVector, ℓq::Real, ∇ℓq::AbstractVector)
    q, ∇ℓq = promote(q, ∇ℓq)
    EvaluatedLogDensity(q, ℓq, ∇ℓq)
end

EvaluatedLogDensity(q, ℓq::Real, ∇ℓq) = EvaluatedLogDensity(collect(q), ℓq, collect(∇ℓq))

"""
$(TYPEDEF)
A point in phase space, consists of a position (in the form of an evaluated log density `ℓ`
at `q`) and a momentum.
"""
struct PhasePoint{D,T,L}
    "Evaluated log density."
    Q::EvaluatedLogDensity{D,T,L}
    "Momentum."
    p::PtrVector{D,T,L,false}
    flag::UInt32
    # function PhasePoint(Q::EvaluatedLogDensity, p::S) where {T,S}
        # @argcheck length(p) == length(Q.q)
        # new{T,S}(Q, p)
    # end
end
PhasePoint(Q, p) = PhasePoint(Q, p, 0x00000000)