text:
devices = Dict{Symbol, DeviceModel}(:Generators => DeviceModel(PSY.ThermalStandard, PSI.ThermalDispatch),
                                    :Loads =>  DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad))
branches = Dict{Symbol, DeviceModel}(:L => DeviceModel(PSY.Line, PSI.StaticLine),
                                     :T => DeviceModel(PSY.Transformer2W, PSI.StaticTransformer),
                                     :TT => DeviceModel(PSY.TapTransformer , PSI.StaticTransformer))
services = Dict{Symbol, PSI.ServiceModel}()
@testset "Operation Model kwargs with CopperPlatePowerModel base" begin
    template = OperationsTemplate(CopperPlatePowerModel, devices, branches, services);
    op_problem = OperationsProblem(TestOpProblem, template,
                                            c_sys5;
                                            optimizer = GLPK_optimizer,
                                           use_parameters = true)
    moi_tests(op_problem, true, 120, 120, 0, 0, 24, false)
#=
  j_model = op_problem.canonical.JuMPmodel
    @test (:params in keys(j_model.ext))
    @test JuMP.num_variables(j_model) == 120
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.Interval{Float64}) == 120
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.EqualTo{Float64}) == 24
    @test !((JuMP.VariableRef, MOI.ZeroOne) in JuMP.list_of_constraint_types(j_model))
    @test JuMP.objective_function_type(j_model) == JuMP.GenericAffExpr{Float64, VariableRef}
=#
    op_problem = OperationsProblem(TestOpProblem, template,
                                            c_sys14;
                                            optimizer = OSQP_optimizer)
    moi_tests(op_problem, false, 120, 120, 0, 0, 24, false)
    #=
    j_model = op_problem.canonical.JuMPmodel
    @test !(:params in keys(j_model.ext))
    @test JuMP.num_variables(j_model) == 120
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.Interval{Float64}) == 120
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.EqualTo{Float64}) == 24
    @test !((JuMP.VariableRef, MOI.ZeroOne) in JuMP.list_of_constraint_types(j_model))
    @test JuMP.objective_function_type(j_model) == JuMP.GenericQuadExpr{Float64, VariableRef}
    =#
    op_problem = OperationsProblem(TestOpProblem, template,
                                            c_sys5_re;
                                            use_forecast_data = false,
                                            optimizer = GLPK_optimizer)
    moi_tests(op_problem, false, 5, 5, 0, 0, 1, false)
    #=
    j_model = op_problem.canonical.JuMPmodel
    @test !(:params in keys(j_model.ext))
    @test JuMP.num_variables(j_model) == 5
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.Interval{Float64}) == 5
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.EqualTo{Float64}) == 1
    @test !((JuMP.VariableRef, MOI.ZeroOne) in JuMP.list_of_constraint_types(j_model))
    @test JuMP.objective_function_type(j_model) == JuMP.GenericAffExpr{Float64, VariableRef}
    =#

    op_problem = OperationsProblem(TestOpProblem, template,
                                            c_sys5_re;
                                            use_forecast_data = false,
                                            parameters = false,
                                            optimizer = GLPK_optimizer)
    moi_tests(op_problem, false, 5, 5, 0, 0, 1, false)
    #=
    j_model = op_problem.canonical.JuMPmodel
    @test !(:params in keys(j_model.ext))
    @test JuMP.num_variables(j_model) == 5
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.Interval{Float64}) == 5
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.LessThan{Float64}) == 0
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.GreaterThan{Float64}) == 0
    @test JuMP.num_constraints(j_model, JuMP.GenericAffExpr{Float64, VariableRef}, MOI.EqualTo{Float64}) == 1
    @test !((JuMP.VariableRef, MOI.ZeroOne) in JuMP.list_of_constraint_types(j_model))
    @test JuMP.objective_function_type(j_model) == JuMP.GenericAffExpr{Float64, VariableRef}
    =#
end

@testset "Operation Model Constructors with Parameters" begin
    networks = [PSI.CopperPlatePowerModel,
                PSI.StandardPTDFModel,
                DCPPowerModel,
                NFAPowerModel,
                ACPPowerModel,
                ACRPowerModel,
                ACTPowerModel,
                DCPLLPowerModel,
                LPACCPowerModel,
                SOCWRPowerModel,
                QCRMPowerModel,
                QCLSPowerModel];

    thermal_gens = [PSI.ThermalStandardUnitCommitment,
                    PSI.ThermalDispatch,
                    PSI.ThermalRampLimited,
                    PSI.ThermalDispatchNoMin];

        systems = [c_sys5,
                   c_sys5_re,
                   c_sys5_bat];

    for net in networks, thermal in thermal_gens, system in systems, p in [true, false]
        @testset "Operation Model $(net) - $(thermal) - $(system)" begin
            thermal_model = DeviceModel(PSY.ThermalStandard, thermal)
            devices = Dict{Symbol, DeviceModel}(:Generators => DeviceModel(PSY.ThermalStandard, thermal),
                                                :Loads =>       DeviceModel(PSY.PowerLoad, PSI.StaticPowerLoad))
            branches = Dict{Symbol, DeviceModel}(:L => DeviceModel(PSY.Line, PSI.StaticLine))
            template = OperationsTemplate(net, devices, branches, services);
            op_problem = OperationsProblem(TestOpProblem,
                                      template,
                                      system; PTDF = PTDF5, use_parameters = p);
        @test :nodal_balance_active in keys(op_problem.canonical.expressions)
        @test (:params in keys(op_problem.canonical.JuMPmodel.ext)) == p
        end

    end

end

#=
@testset "Build Operation Models" begin
    #SCED = PSI.SCEconomicDispatch(c_sys5; optimizer = GLPK_optimizer);
    #OPF = PSI.OptimalPowerFlow(c_sys5, ACPPowerModel, optimizer = ipopt_optimizer)
    #UC = PSI.UnitCommitment(c_sys5, DCPPowerModel; optimizer = GLPK_optimizer)

    #ED_rts_p = PSI.EconomicDispatch(c_rts, PSI.CopperPlatePowerModel; optimizer = GLPK_optimizer);
    #ED_rts = PSI.EconomicDispatch(c_rts, PSI.CopperPlatePowerModel; optimizer = GLPK_optimizer, use_parameters = false);
    # These other tests can be enabled when CDM parser get the correct HVDC type.
    #OPF_rts = PSI.OptimalPowerFlow(sys_rts, PSI.CopperPlatePowerModel, optimizer = ipopt_optimizer)
    #UC_rts = PSI.UnitCommitment(sys_rts, PSI.CopperPlatePowerModel; optimizer = GLPK_optimizer, use_parameters = false)
end

=#


text:
module TestShow
    using DataFrames
    using RDatasets
    movies = data("ggplot2", "movies")
    movies = movies[1:50, :]
    show(movies)
    show(movies, true)
    showall(movies)
    showall(movies, true)

    subdf = select(:(rating .> 6.4), movies)
    show(subdf)
    show(subdf, true)
    showall(subdf)
    showall(subdf, true)
end

ext: jl


text:
module TestMap

using Folds
using Test
include("utils.jl")

@testset begin
    @test Folds.map(identity, 1:10) == 1:10
    @test Folds.map(identity, reshape(1:6, 2, 3))::Matrix{Int} == reshape(1:6, 2, 3)
    @test Folds.map(tuple, 1:3, 4:6, 7:9) == collect(zip(1:3, 4:6, 7:9))
    @test Folds.map(identity, 1:10, SequentialEx()) == 1:10
    @test Folds.map(tuple, 1:3, 4:6, 7:9, SequentialEx()) == collect(zip(1:3, 4:6, 7:9))
    @test Folds.map(identity, reshape(1:6, 2, 3), SequentialEx())::Matrix{Int} ==
          reshape(1:6, 2, 3)
end

@testset "executor with keyword arguments" begin
    err = @test_error Folds.map(identity, 1:2, SequentialEx(); foo = :bar)
    @test "not accept any keyword arguments" ∈ᵉʳʳᵒʳ err
end

end  # module


text:
function rbc()
	∞ = Inf
	@variables α, β, ρ, δ, σ
	@variables t::Integer, k(..), z(..), c(..), q(..), i(..)
	@variables Ω_1
	x = [k, z]
	y = [c, q, i]
	p = [α, β, ρ, δ, σ, Ω_1]

	H = [1 / c(t) - (β / c(t + 1)) * (α * exp(z(t + 1)) * k(t + 1)^(α - 1) + (1 - δ)),
	c(t) + k(t + 1) - (1 - δ) * k(t) - q(t),
	q(t) - exp(z(t)) * k(t)^α,
	z(t + 1) - ρ * z(t),
	i(t) - (k(t + 1) - (1 - δ) * k(t))]

	steady_states = [k(∞) ~ (((1 / β) - 1 + δ) / α)^(1 / (α - 1)), z(∞) ~ 0,
				c(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)) -
						δ * (((1 / β) - 1 + δ) / α)^(1 / (α - 1)),
				q(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)),
				i(∞) ~ δ * (((1 / β) - 1 + δ) / α)^(1 / (α - 1))]

	steady_states_iv = [k(∞) ~ (((1 / β) - 1 + δ) / α)^(1 / (α - 1)), z(∞) ~ 0,
				c(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)) -
						δ * (((1 / β) - 1 + δ) / α)^(1 / (α - 1)),
				q(∞) ~ (((1 / β) - 1 + δ) / α)^(α / (α - 1)),
				i(∞) ~ δ * (((1 / β) - 1 + δ) / α)^(1 / (α - 1))]

	n_ϵ = 1
	n_x = length(x)
	n_y = length(y)
	n_p = length(p)
	Γ = reshape([σ], n_ϵ, n_ϵ)
	η = reshape([0; -1], n_x, n_ϵ) # η is n_x * n_ϵ matrix

	n_z = 2 # number of observables
	Q = zeros(n_z, n_x + n_y) # the order is [y, x]
	Q[1, 1] = 1.0 # c
	Q[2, 3] = 1.0 # i

	Ω = [Ω_1, Ω_1]

	return H, (; t, x, y, p, steady_states, steady_states_iv, Γ, η, Ω, Q), "rbc"
end

max_stars_repo_name: HighDimensionalEconLab/HMCExamples

max_issues_repo_name: HighDimensionalEconLab/HMCExamples

max_forks_repo_name: HighDimensionalEconLab/HMCExamples


text:
# Author: Ali Siahkoohi, alisk@gatech.edu
#         Mathias Louboutin
# Date: December 2020
# Copyright: Georgia Institute of Technology, 2020

module ExtendedConv

using Flux
using LinearAlgebra
using PyPlot: Figure
using JOLI
using SetIntersectionProjection
using SparseArrays
using Optim

import Base.+, Base.*
import DrWatson: _wsave

# Utilities
include("./utils/savefig.jl")

# Filter-wise normalization method
include("./visualization/visualization.jl")

include("./extension/extension.jl")
include("./extension/flux_optim.jl")

end


text:
using LinearAlgebra

mutable struct TaylorSeries{T,S<:Real}
    coeffs::Vector{T}
    offset::S
end
TaylorSeries(coeffs::Vector{T}) where T = TaylorSeries(coeffs, 0)

Base.length(t::TaylorSeries) = length(t.coeffs)
degree(t::TaylorSeries) = length(t) - 1
offset(t::TaylorSeries) = t.offset

function shift!(t::TaylorSeries{T}, b) where T
    a = offset(t)
    a == b && return t
    for k = 0:degree(t)
        for l = degree(t):-1:k+1
            if T <: AbstractArray
                axpy!((b-a)^(l-k)*binomial(l,k), t[l], t[k])
            else
                t[k] += t[l]*(b-a)^(l-k)*binomial(l,k)
            end
        end
    end
    t.offset = b
    return t
end
shift(t::TaylorSeries, b) = b == offset(t) ? t : shift!(copy(t), b)

Base.copy(t::TaylorSeries) = TaylorSeries([copy(t[k]) for k = 0:degree(t)], offset(t))

function Base.getindex(t::TaylorSeries, k)
    k < 0 && return BoundsError(t, k)
    @inbounds if k > degree(t)
        return zero(t.coeffs[1])
    else
        return t.coeffs[k+1]
    end
end
function Base.setindex!(t::TaylorSeries, v, k)
    k < 0 && return BoundsError(t, k)
    @inbounds if k > degree(t)
        t0 = t[0]
        while length(t) < k
            push!(t.coeffs, zero(t0))
        end
        push!(t.coeffs, v)
        return v
    else
        setindex!(t.coeffs, v, k+1)
        return v
    end
end

function (t::TaylorSeries)(x::Real)
    k = degree(t)
    xa = x - offset(t)
    (xa == zero(xa) || k == 0) && return copy(t[0])
    v = copy(t[k])
    if v isa AbstractArray
        while k > 0
            v = LinearAlgebra.axpby!(one(xa), t[k-1], xa, v)
            k -= 1
        end
    else
        while k > 0
            v = muladd(xa, v, t[k-1])
            k -= 1
        end
    end
    return v
end

Base.eltype(t::TaylorSeries) = eltype(typeof(t))
Base.eltype(::Type{TaylorSeries{<:Real,T}}) where {T} = eltype(T)

function trunc!(t::TaylorSeries; kmax::Integer = length(F)-1, tol::Real = 0)
    ktol = findlast(x->norm(x)>=tol, t.coeffs) - 1
    kmax = min(kmax, ktol === nothing ? 0 : ktol)
    if degree(t) > kmax
        resize!(t.coeffs, kmax+1)
    end
    return t
end

function setdegree!(t::TaylorSeries, K::Int)
    while degree(t) < K
        push!(t.coeffs, zero(t[0]))
    end
    if degree(t) > K
        resize!(t.coeffs, K+1)
    end
    return t
end

function Base.copyto!(tdst::TaylorSeries, tsrc::TaylorSeries)
    @assert offset(tdst) == offset(tsrc)
    setdegree!(tdst, degree(tsrc))
    for j = 0:degree(tsrc)
        tdst[j] = copy(tsrc[j])
    end
    return tdst
end

Base.fill!(t::TaylorSeries{<:Number}, α) = (fill!(t.coeffs, α); return t)
Base.fill!(t::TaylorSeries{<:AbstractArray}, α) = (map!(x->fill!(x, α), t.coeffs, t.coeffs); return t)

LinearAlgebra.rmul!(t::TaylorSeries, α) = (rmul!(t.coeffs, α); return t)
LinearAlgebra.lmul!(α, t::TaylorSeries) = (lmul!(α, t.coeffs); return t)

function LinearAlgebra.mul!(tdst::TaylorSeries, α, tsrc::TaylorSeries)
    @assert offset(tdst) == offset(tsrc)
    setdegree!(tdst, degree(tsrc))
    mul!(tdst.coeffs, α, tsrc.coeffs)
    return tdst
end

function LinearAlgebra.mul!(tdst::TaylorSeries, tsrc::TaylorSeries, α)
    @assert offset(tdst) == offset(tsrc)
    setdegree!(tdst, degree(tsrc))
    mul!(tdst.coeffs, tsrc.coeffs, α)
    return tdst
end

function LinearAlgebra.axpy!(α, tx::TaylorSeries, ty::TaylorSeries)
    @assert offset(tx) == offset(ty)
    Kx = degree(tx)
    setdegree!(ty, max(tx, degree(ty)))
    Ky = degree(ty)
    if Ky > Kx
        LinearAlgebra.axpy!(α, tx.coeffs, view(ty.coeffs, 1:(Kx+1)))
    else
        LinearAlgebra.axpy!(α, tx.coeffs, ty.coeffs)
    end
    return ty
end
function LinearAlgebra.axpby!(α, tx::TaylorSeries, β, ty::TaylorSeries)
    @assert offset(tx) == offset(ty)
    Kx = degree(tx)
    setdegree!(ty, max(tx, degree(ty)))
    Ky = degree(ty)
    if Ky > Kx
        LinearAlgebra.axpby!(α, tx.coeffs, β, view(ty.coeffs, 1:(Kx+1)))
        lmul!(β, view(ty.coeffs, (Kx+2):(Ky+1)))
    else
        LinearAlgebra.axpby!(α, tx.coeffs, β, ty.coeffs)
    end
    return ty
end

function LinearAlgebra.axpy!(α, tx::TaylorSeries{<:AbstractArray},
                                ty::TaylorSeries{<:AbstractArray})
    @assert offset(tx) == offset(ty)
    Kx = degree(tx)
    setdegree!(ty, max(Kx, degree(ty)))
    for j = 0:Kx
        LinearAlgebra.axpy!(α, tx[j], ty[j])
    end
    return ty
end
function LinearAlgebra.axpby!(α, tx::TaylorSeries{<:AbstractArray},
                                β, ty::TaylorSeries{<:AbstractArray})
    @assert offset(tx) == offset(ty)
    Kx = degree(tx)
    setdegree!(ty, max(Kx, degree(ty)))
    Ky = degree(ty)
    for j = 0:Kx
        LinearAlgebra.axpby!(α, tx[j], β, ty[j])
    end
    for j = Kx+1:Ky
        LinearAlgebra.lmul!(β, ty[j])
    end
    return ty
end

Base.zero(t::TaylorSeries) = TaylorSeries([zero(t[0])], offset(t))
Base.one(t::TaylorSeries) = TaylorSeries([one(t[0])], offset(t))

function Base.similar(t::TaylorSeries{T}) where T
    if isbitstype(T)
        return TaylorSeries(similar(t.coeffs), offset(t))
    else
        return TaylorSeries(similar.(t.coeffs), offset(t))
    end
end

function Base.:+(t1::TaylorSeries, t2::TaylorSeries)
    @assert offset(t1) == offset(t2)
    K = min(degree(t1), degree(t2))
    coeffs = [t1[j]+t2[j] for j = 0:K]
    for j = K+1:degree(t1)
        push!(coeffs, copy(t1[j]))
    end
    for j = K+1:degree(t2)
        push!(coeffs, copy(t2[j]))
    end
    return TaylorSeries(coeffs, offset(t1))
end

function Base.:-(t1::TaylorSeries, t2::TaylorSeries)
    @assert offset(t1) == offset(t2)
    K = min(degree(t1), degree(t2))
    coeffs = [t1[j]-t2[j] for j = 0:K]
    for j = K+1:degree(t1)
        push!(coeffs, copy(t1[j]))
    end
    for j = K+1:degree(t2)
        push!(coeffs, -t2[j])
    end
    return TaylorSeries(coeffs, offset(t1))
end

Base.:-(t::TaylorSeries) = TaylorSeries(-t.coeffs, offset(t))

Base.:*(t::TaylorSeries, a::Number) = TaylorSeries([c*a for c in t.coeffs], offset(t))
Base.:*(a::Number, t::TaylorSeries) = TaylorSeries([a*c for c in t.coeffs], offset(t))
Base.:/(t::TaylorSeries, a::Number) = TaylorSeries([c/a for c in t.coeffs], offset(t))
Base.:\(a::Number, t::TaylorSeries) = TaylorSeries([a\c for c in t.coeffs], offset(t))

Base.conj(t::TaylorSeries) = TaylorSeries(map(conj, t.coeffs), offset(t))
Base.adjoint(t::TaylorSeries) = TaylorSeries(map(adjoint, t.coeffs), offset(t))
Base.transpose(t::TaylorSeries) = TaylorSeries(map(transpose, t.coeffs), offset(t))

function Base.:*(t1::TaylorSeries, t2::TaylorSeries; Kmax = degree(t1) + degree(t2))
    @assert offset(t1) == offset(t2)
    K = min(Kmax, degree(t1) + degree(t2))
    coeffs = [zero(t1[0])*zero(t2[0]) for k = 0:K]
    t = TaylorSeries(coeffs, offset(t1))
    return mul!(t, t1, t2, true, true)
end

function LinearAlgebra.mul!(t::TaylorSeries, t1::TaylorSeries, t2::TaylorSeries,
                                α = true, β = false)
    @assert offset(t) == offset(t1) == offset(t2)
    K1 = degree(t1)
    K2 = degree(t2)
    K = degree(t)
    for k = 0:K
        tk = t[k]
        if tk isa AbstractArray
            if β != 1
                rmul!(tk, β)
            end
            for k1 = max(0, k-K2):min(k, K1)
                tk1 = t1[k1]
                tk2 = t2[k-k1]
                if tk1 isa AbstractArray && tk2 isa AbstractArray
                    mul!(tk, tk1, tk2, α, true)
                elseif tk1 isa AbstractArray
                    axpy!(tk2*α, tk1, tk)
                elseif tk2 isa AbstractArray
                    axpy!(tk1*α, tk2, tk)
                else
                    t[k] += tk1*tk2*α
                end
            end
        else
            tk *= β
            for k1 = max(0, k-K2):min(k, K1)
                tk += t1[k1]*t2[k-k1]*α
            end
            t[k] = tk
        end
    end
    return t
end

function differentiate(t::TaylorSeries)
    if degree(t) >= 1
        return TaylorSeries([k*t[k] for k = 1:degree(t)], offset(t))
    else
        return zero(t)
    end
end
function integrate(t::TaylorSeries, (a,b)::Tuple{Real,Real})
    s = t[0]*(b-a)
    c = offset(t)
    for k = 1:degree(t)
        s += t[k]*((b-c)^(k+1) - (a-c)^(k+1))/(k+1)
    end
    return s
end

function localdot(t1::TaylorSeries, t2::TaylorSeries)
    @assert offset(t1) == offset(t2)
    K1 = degree(t1)
    K2 = degree(t2)
    coeffs = let K = K1 + K2
        [sum(dot(t1[k1],t2[k-k1]) for k1 = max(0,k-K2):min(k,K1)) for k=0:K]
    end
    return TaylorSeries(coeffs, offset(t1))
end

LinearAlgebra.tr(t::TaylorSeries) = TaylorSeries(map(tr, t.coeffs), offset(t))

Base.real(t::TaylorSeries) = TaylorSeries(map(real, t.coeffs), offset(t))
Base.imag(t::TaylorSeries) = TaylorSeries(map(imag, t.coeffs), offset(t))

function Base.inv(t::TaylorSeries; tol = 1e-10, Kmax = 10*degree(t))
    t0inv = inv(t[0])
    coeffs = [t0inv]
    K = degree(t)
    for k = 1:Kmax
        invtk = -t0inv*sum(coeffs[i+1]*t[k-i] for i = max(0,k-K):(k-1))
        if norm(invtk) < tol && norm(coeffs[end]) < tol
            break
        end
        push!(coeffs, invtk)
    end
    return TaylorSeries(coeffs, offset(t))
end


text:
# main.jl
#
# A projector QMC program for the TFIM

include("lattice.jl") #define the spatial lattice
include("updates.jl") #functions for the Monte Carlo updates
include("measurements.jl") #functions for the Monte Carlo updates

#Projector parameters
M = 300 #length of the projector operator_list is 2M
h_x = 0.5
J_ = 1.0
MCS = 10000 #the number of Monte Carlo steps

#*******  Globals
spin_left = fill(1,nSpin) #left and right trail spin state
spin_right = fill(1,nSpin)
LinkList =[]
LegType =[]
Associates =[]
operator_list = fill(0,(2*M,2))
#  (-2,i) is an off-diagonal site operator h(sigma^+_i + sigma^-_i)
#  (-1,i) is a diagonal site operator h
#  (0,0) is the identity operator I - NOT USED IN THE PROJECTOR CASE
#  (i,j) is a diagonal bond operator J(sigma^z_i sigma^z_j)
#*******  Globals

for i = 1:2000  #Equilibration
    DiagonalUpdate()
    LinkedList()
    ClusterUpdate()
end

M2 = 0
for i = 1:MCS #Monte Carlo Production Steps
    DiagonalUpdate()
    LinkedList()
	global M2 += Measure()
    ClusterUpdate()
end

println(M2/(nSpin*nSpin*MCS))


text:
#   *
function pointer_all(list::Array{Pointer,1},r::Int,c::Int)
    push!(list,Pointer(UInt64[0,0,0],[r,c],0,0,false,true))
    push!(list,Pointer(UInt64[0,0,0],[r,c],1,0,false,true))
    push!(list,Pointer(UInt64[0,0,0],[r,c],2,0,false,true))
    push!(list,Pointer(UInt64[0,0,0],[r,c],3,0,false,true))
    push!(list,Pointer(UInt64[0,0,0],[r,c],4,0,false,true))
    push!(list,Pointer(UInt64[0,0,0],[r,c],5,0,false,true))
end

#   \
function pointer_diag1(list::Array{Pointer,1},r::Int,c::Int)
    push!(list,Pointer(UInt64[0,0,0],[r,c],2,0,false,true))
    push!(list,Pointer(UInt64[0,0,0],[r,c],5,0,false,true))

end

#   /
function pointer_diag2(list::Array{Pointer,1},r::Int,c::Int)
    push!(list,Pointer(UInt64[0,0,0],[r,c],1,0,false,true))
    push!(list,Pointer(UInt64[0,0,0],[r,c],4,0,false,true))
end

#   _
function pointer_horiz(list::Array{Pointer,1},r::Int,c::Int)
    push!(list,Pointer(UInt64[0,0,0],[r,c],0,0,false,true))
    push!(list,Pointer(UInt64[0,0,0],[r,c],3,0,false,true))
end

#   >
function redir_r(list::Array{Pointer,1},ind::Int)
    list[ind].dir=0
end

#   d
function redir_ur(list::Array{Pointer,1},ind::Int)
    list[ind].dir=1
end

#   b
function redir_ul(list::Array{Pointer,1},ind::Int)
    list[ind].dir=2
end

#   <
function redir_l(list::Array{Pointer,1},ind::Int)
    list[ind].dir=3
end

#   p
function redir_dl(list::Array{Pointer,1},ind::Int)
    list[ind].dir=4
end

#   q
function redir_dr(list::Array{Pointer,1},ind::Int)
    list[ind].dir=5
end

#   a
function rotate_cw(list::Array{Pointer,1},ind::Int)
    list[ind].dir=(list[ind].dir+1)%6
end

#   x
function rotate_ccw(list::Array{Pointer,1},ind::Int)
    list[ind].dir=(list[ind].dir+5)%6
end

#   s
function mirror_main1(list::Array{Pointer,1},ind::Int)
    list[ind].dir==0 ? list[ind].dir=4 :
    list[ind].dir==1 ? list[ind].dir=3 :
    list[ind].dir==3 ? list[ind].dir=1 :
    list[ind].dir==4 ? list[ind].dir=0 : nothing
end

#   t
function mirror_main2(list::Array{Pointer,1},ind::Int)
    list[ind].dir==0 ? list[ind].dir=2 :
    list[ind].dir==2 ? list[ind].dir=0 :
    list[ind].dir==3 ? list[ind].dir=5 :
    list[ind].dir==5 ? list[ind].dir=0 : nothing
end

#   u
function mirror_main3(list::Array{Pointer,1},ind::Int)
    list[ind].dir==1 ? list[ind].dir=5 :
    list[ind].dir==2 ? list[ind].dir=4 :
    list[ind].dir==4 ? list[ind].dir=2 :
    list[ind].dir==5 ? list[ind].dir=1 : nothing
end

#   j
function mirror_half1(list::Array{Pointer,1},ind::Int)
    list[ind].dir==0 ? list[ind].dir=3 :
    list[ind].dir==1 ? list[ind].dir=2 :
    list[ind].dir==2 ? list[ind].dir=1 :
    list[ind].dir==3 ? list[ind].dir=0 :
    list[ind].dir==4 ? list[ind].dir=5 :
    list[ind].dir==5 ? list[ind].dir=4 : nothing
end

#   k
function mirror_half2(list::Array{Pointer,1},ind::Int)
    list[ind].dir==0 ? list[ind].dir=1 :
    list[ind].dir==1 ? list[ind].dir=0 :
    list[ind].dir==2 ? list[ind].dir=5 :
    list[ind].dir==3 ? list[ind].dir=4 :
    list[ind].dir==4 ? list[ind].dir=3 :
    list[ind].dir==5 ? list[ind].dir=2 : nothing
end

#   l
function mirror_half3(list::Array{Pointer,1},ind::Int)
    list[ind].dir==0 ? list[ind].dir=5 :
    list[ind].dir==1 ? list[ind].dir=4 :
    list[ind].dir==2 ? list[ind].dir=3 :
    list[ind].dir==3 ? list[ind].dir=2 :
    list[ind].dir==4 ? list[ind].dir=1 :
    list[ind].dir==5 ? list[ind].dir=0 : nothing
end

#   O
function mirror_all(list::Array{Pointer,1},ind::Int)
    list[ind].dir==0 ? list[ind].dir=3 :
    list[ind].dir==1 ? list[ind].dir=4 :
    list[ind].dir==2 ? list[ind].dir=5 :
    list[ind].dir==3 ? list[ind].dir=0 :
    list[ind].dir==4 ? list[ind].dir=1 :
    list[ind].dir==5 ? list[ind].dir=2 : nothing
end

#   m
function catch_diag1(list::Array{Pointer,1},ind::Int)
    list[ind].dir==2 ? list[ind].alive=false :
    list[ind].dir==5 ? list[ind].alive=false : nothing
end

#   n
function catch_diag2(list::Array{Pointer,1},ind::Int)
    list[ind].dir==1 ? list[ind].alive=false :
    list[ind].dir==4 ? list[ind].alive=false : nothing
end

#   o
function catch_diag3(list::Array{Pointer,1},ind::Int)
    list[ind].dir==0 ? list[ind].alive=false :
    list[ind].dir==3 ? list[ind].alive=false : nothing
end

#   #
function catch_all(list::Array{Pointer,1},ind::Int)
    list[ind].alive=false
end

#   X (dir 012345 spread)
function pointer_spread(list::Array{Pointer,1},ind::Int,r::Int,c::Int)
    for d=0:5
        if d==list[ind].dir
            move(list,ind)
        elseif d==(list[ind].dir+3)%6
            nothing
        else
            push!(list,Pointer(deepcopy(list[ind].lstack),[r,c],d,0,false,true))
        end
    end
end

#   E (dir 03 spread)
function pointer_spread03(list::Array{Pointer,1},ind::Int,r::Int,c::Int)
    if list[ind].dir==0 || list[ind].dir==3
        nothing
    else
        list[ind].dir=0
        push!(list,Pointer(deepcopy(list[ind].lstack),[r,c],3,0,false,true))
    end
end

#   H (dir 14 spread)
function pointer_spread14(list::Array{Pointer,1},ind::Int,r::Int,c::Int)
    if list[ind].dir==1 || list[ind].dir==4
        nothing
    else
        list[ind].dir=1
        push!(list,Pointer(deepcopy(list[ind].lstack),[r,c],4,0,false,true))
    end
end

#   W (dir 25 spread)
function pointer_spread25(list::Array{Pointer,1},ind::Int,r::Int,c::Int)
    if list[ind].dir==2 || list[ind].dir==5
        nothing
    else
        list[ind].dir=2
        push!(list,Pointer(deepcopy(list[ind].lstack),[r,c],5,0,false,true))
    end
end

#   '
function iftopzero(list::Array{Pointer,1},ind::Int)
    if list[ind].lstack[end]==0
        move(list,ind)
    end
end

#   "
function iftopgrtzero(list::Array{Pointer,1},ind::Int)
    if list[ind].lstack[end]>0
        move(list,ind)
    end
end

#   K
function iftopeqsec(list::Array{Pointer,1},ind::Int)
    if list[ind].lstack[end]==list[ind].lstack[end-1]
        move(list,ind)
    end
end

#   L
function iftopgrtsec(list::Array{Pointer,1},ind::Int)
    if list[ind].lstack[end]>list[ind].lstack[end-1]
        move(list,ind)
    end
end

#   Q
function skipnext(list::Array{Pointer,1},ind::Int)
    move(list,ind)
end

#   v
function wait1(list::Array{Pointer,1},ind::Int)
    if list[ind].printstate==false
        w=list[ind].wait
        w=(w+1)%2
        w == 0 ? move(list,ind) : nothing
    end
end

#   ^
function wait2(list::Array{Pointer,1},ind::Int)
    if list[ind].printstate==false
        w=list[ind].wait
        w=(w+2)%3
        w==0 ? move(list,ind) : nothing
    end
end

#   J
function jumpto(list::Array{Pointer,1},ind::Int)
    list[ind].loc=[list[ind].lstack[end],list[ind].lstack[end-1]]
end

#   g
function globalgetfirst(list::Array{Pointer,1},ind::Int,gstack::Array{UInt64,1})
    length(gstack)>0 ? list[ind].lstack[end]=gstack[end] : nothing
end

#   f
function localfirst2global(list::Array{Pointer,1},ind::Int,gstack::Array{UInt64,1})
    push!(gstack,list[ind].lstack[end])
end

#   e
#   push lstack to gstack and initialize lstack to [0,0,0]
function localflush2global(list::Array{Pointer,1},ind::Int,gstack::Array{UInt64,1})
    append!(gstack,reverse(list[ind].lstack))
    fill!(list[ind].lstack,0)
end

#   U
#   push top 3 gstack on lstack in reverse order and pop top 3 gstack
function globalflush2local(list::Array{Pointer,1},ind::Int,gstack::Array{UInt64,1})
    list[ind].lstack[end-2]=gstack[end]
    list[ind].lstack[end-1]=gstack[end-1]
    list[ind].lstack[end]=gstack[end-2]
    for i=1:3
        pop!(gstack)
    end
end

#   ~
function localflip12!(list::Array{Pointer,1},ind::Int)
    push!(list[ind].lstack,splice!(list[ind].lstack,length(list[ind].lstack)-1))
end

#   @
function localflip13!(list::Array{Pointer,1},ind::Int)
    list[ind].lstack=reverse(list[ind].lstack)
end

#   F
function localallfirst!(list::Array{Pointer,1},ind::Int)
    fill!(list[ind].lstack,list[ind].lstack[end])
end

#   z
function localallzero!(list::Array{Pointer,1},ind::Int)
    fill!(list[ind].lstack,0)
end

#   y
function globalrotdown!(list::Array{Pointer,1},ind::Int,gstack::Array{UInt64,1})
    if length(gstack)>0
        depth=list[ind].lstack[end]
        steps=list[ind].lstack[end-1]
        if depth > length(gstack) || depth == 0
            depth = length(gstack)
        end
        append!(gstack,splice!(gstack,length(gstack)-depth+1:length(gstack)-depth+steps%depth))
    end
end

#   h
function globalrotup!(list::Array{Pointer,1},ind::Int,gstack::Array{UInt64,1})
    depth=list[ind].lstack[end]
    steps=list[ind].lstack[end-1]
    if depth > length(gstack) || depth == 0
        depth = length(gstack)
    end
    append!(gstack,splice!(gstack,(length(gstack)-depth+1):length(gstack)-steps%depth))
end

#   =
function globalduptop!(gstack::Array{UInt64,1})
    length(gstack)>0 ? gstack=push!(gstack,gstack[end]) : nothing
end

#   ?
function globalpop!(gstack::Array{UInt64,1})
    length(gstack)>0 ? pop!(gstack) : nothing
end

#   A
function globalstacklen!(gstack::Array{UInt64,1})
    gstack=push!(gstack,length(gstack))

end

#code manipulation

#   D (absolute addressing)
function adropto(list::Array{Pointer,1},ind::Int,arena::Honeycomb)
    r=list[ind].lstack[end-1]
    c=list[ind].lstack[end-2]
    dropstuff(list,ind,arena,Int128(r),Int128(c))
end

function dropstuff(list::Array{Pointer,1},ind::Int,arena::Honeycomb,r::Int128,c::Int128)
    rows=maximum(length(arena.a[:,1]))
    cols=maximum(length(arena.a[1,:]))
    # Padding
    r==0 ? arena.a=cat(1,transpose(fill(' ',cols)),arena.a) :
    r>rows ? (delta=Int(r-rows);arena.a=cat(1,arena.a,reshape(fill(' ',delta*cols),delta,cols))) : nothing
    rows=maximum(length(arena.a[:,1]))
    cols=maximum(length(arena.a[1,:]))
    c==0 ? arena.a=cat(2,(fill(' ',rows)),arena.a) :
    c>cols ? (delta=Int(c-cols);arena.a=cat(2,arena.a,reshape(fill(' ',rows*delta),rows,delta))) : nothing
    m=[0,0]
    (r==0 && c!=0) ? m=[1,0] :
    (r==0 && c==0) ? m=[1,1] :
    (c==0 && r!=0) ? m=[0,1] :
    (c!=0 && r!=0) ? m=[0,0] : nothing

    arena.a[list[ind].lstack[end-1]+m[1],list[ind].lstack[end-2]+m[2]]=Char(list[ind].lstack[end])

    if m!=[UInt64(0),UInt64(0)]
        @inbounds for k=1:length(list)
            slide(list,k,m)
        end
    end
end

#   G (absolute addressing)
function agetfrom(list::Array{Pointer,1},ind::Int,arena::Honeycomb)
    rows=maximum(length(arena.a[:,1]))
    cols=maximum(length(arena.a[1,:]))
    r=list[ind].lstack[end-1]
    c=list[ind].lstack[end-2]

    if r==0 || r>rows || c==0 || c>cols
        list[ind].lstack[end]=0
    else
        list[ind].lstack[end]=arena.a[r,c]
    end
end

#   Y (relative addressing)
function rdropto(list::Array{Pointer,1},ind::Int,arena::Honeycomb)
    rows=maximum(length(arena.a[:,1]))
    cols=maximum(length(arena.a[1,:]))
    r=list[ind].lstack[end-1]
    c=list[ind].lstack[end-2]
    x=list[ind].loc[1]
    y=list[ind].loc[2]
    ro=reinterpret(Int64,r)
    co=reinterpret(Int64,c)

    if x+ro<1 || y+co<1 || x+ro>rows || y+co>cols
        list[ind].lstack[end]=0
    else
        arena.a[x+ro,y+co]=list[ind].lstack[end]
    end
end

#   Z (relative addressing)
function rgetfrom(list::Array{Pointer,1},ind::Int,arena::Honeycomb)
    rows=maximum(length(arena.a[:,1]))
    cols=maximum(length(arena.a[1,:]))
    r=list[ind].lstack[end-1]
    c=list[ind].lstack[end-2]
    x=list[ind].loc[1]
    y=list[ind].loc[2]
    ro=reinterpret(Int64,r)
    co=reinterpret(Int64,c)

    if x+ro<1 || y+co<1 || x+ro>rows || y+co>cols
        list[ind].lstack[end]=0
    else
        list[ind].lstack[end]=arena.a[x+ro,y+co]
    end
end

# arithmetic

#   +
function add(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]=list[ind].lstack[end] + list[ind].lstack[end-1]
end

#   -
function sub(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]=list[ind].lstack[end] - list[ind].lstack[end-1]
end

#   .
function mul(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]=list[ind].lstack[end] * list[ind].lstack[end-1]
end

#   :
function intdiv(list::Array{Pointer,1},ind::Int)
    try
        list[ind].lstack[end]=div(list[ind].lstack[end],list[ind].lstack[end-1])
    catch
        error("Instruction ':': Division by zero error!")
    end
end

#   %
function modulus(list::Array{Pointer,1},ind::Int)
    try
        list[ind].lstack[end]=list[ind].lstack[end] % list[ind].lstack[end-1]
    catch
        error("Instruction '%': Modulo by zero error!")
    end
end

#   0,1,2,3,4,5,6,7,8,9
function setnum(list::Array{Pointer,1},ind::Int,arena::Honeycomb)
    list[ind].lstack[end]=parse(UInt64,arena.a[list[ind].loc[1],list[ind].loc[2]])
end

#   P
function increment(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]+=1
end

#   M
function decrement(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]-=1
end

#   B
function pow(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]=list[ind].lstack[end]^list[ind].lstack[end-1]
end

# bitwise operations

#   &
function bitand(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]=list[ind].lstack[end]&list[ind].lstack[end-1]
end

#   |
function bitor(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]=list[ind].lstack[end]|list[ind].lstack[end-1]
end

#   $
function bitxor(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]=list[ind].lstack[end]$list[ind].lstack[end-1]
end

#   !
function bitnot(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]=~list[ind].lstack[end]
end

#   (
function bitshiftleft(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]=list[ind].lstack[end]<<list[ind].lstack[end-1]
end

#   )
function bitshiftright(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]=list[ind].lstack[end]>>>list[ind].lstack[end-1]
end

#   [
function bitrollleft(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]=list[ind].lstack[end]<<(list[ind].lstack[end-1]%64)+list[ind].lstack[end]>>>(64-list[ind].lstack[end-1]%64)
end

#   ]
function bitrollright(list::Array{Pointer,1},ind::Int)
    list[ind].lstack[end]=list[ind].lstack[end]>>>(list[ind].lstack[end-1]%64)+list[ind].lstack[end]<<(64-list[ind].lstack[end-1]%64)
end

#I/O

#   c
function ginputchar(gstack::Array{UInt64,1})
    input=readline(STDIN)
    inp=Char(input[1])
    gstack=push!(gstack,inp)
end

#   V
function ginputstring(gstack::Array{UInt64,1})
    input=readline(STDIN)
    inp=input[1:end]
    for c=1:strwidth(inp) #cut off \r\n
        gstack=push!(gstack,inp[chr2ind(inp,c)])
    end
    push!(gstack,0x0a) #append newline, \n as newline on all systems
end

#   i
function ginputint(gstack::Array{UInt64,1})
    input=readline(STDIN)
    gstack=push!(gstack,parse(UInt64,input))
end

#   C
function goutputchar(gstack::Array{UInt64,1})
    print(Char(gstack[end]))
end

#   I
function goutputint(gstack::Array{UInt64,1})
    print(Int128(gstack[end]))
end

#   ,
function linputchar(list::Array{Pointer,1},ind::Int)
    input=readline(STDIN)
    inp=(input[1])
    list[ind].lstack[end]=inp
end

#   T
function linputint(list::Array{Pointer,1},ind::Int)
    input=readline(STDIN)
    list[ind].lstack[end]=parse(UInt64,input)
end

#   }
function loutputchar(list::Array{Pointer,1},ind::Int)
    c=list[ind].lstack[end]
    print(Char(c))
end

#   {
function loutputint(list::Array{Pointer,1},ind::Int)
    outint=Int128(list[ind].lstack[end])
    print(outint)
end

#   `
function toggleoutput(list::Array{Pointer,1},ind::Int)
    list[ind].printstate= !list[ind].printstate
end

#   N
function newline()
    print("\n")
end

#   r
function readfile(list::Array{Pointer,1},ind::Int,gstack::Array{UInt64,1})
    code=()
    namebytes=list[ind].lstack[end]
    if namebytes>length(gstack)
        error("$(list[ind].loc)\nAmount of file name bytes longer than gstack!")
    else
        bytearray=reinterpret(UInt8,gstack)
        namechars=collect(Char,bytearray)[end:-1:end-BigInt(namebytes)+1]
        name=AbstractString(namechars)
        f=open(name,"r")
        try
            code=readbytes!(f)
        finally
            close(f)
        end
    end
    padbytesnum=8-length(code)%8
    padbytesarray=collect(zeros(UInt8,padbytesnum))
    padarray8=cat(1,padbytesarray,code)
    padarray64=reinterpret(UInt64,padarray8)
    gstack=cat(1,gstack,padarray64)
end

#   w
function writefile(list::Array{Pointer,1},ind::Int,gstack::Array{UInt64,1})
    namebytes=list[ind].lstack[end]
    filebytes=list[ind].lstack[end-1]
    bytearray=reinterpret(UInt8,gstack)
    if namebytes >= length(bytearray)
        error("$(list[ind].loc)\nFilename taking up all available bytes ($(namebytes)) in the stack ($(length(bytearray)))!")
    elseif namebytes==0
        error("$(list[ind].loc)\nZero length filename!")
    elseif namebytes+filebytes> length(bytearray)
        error("$(list[ind].loc)\nByte sum of file name and file content ($namebytes + $filebytes = $(namebytes+filebytes) bytes) bigger than the stack ($(length(bytearray)) bytes)!")
    else
        namechars=collect(Char,bytearray)[end:-1:end-BigInt(namebytes)+1]
        name=AbstractString(collect(namechars))
        filedata=collect(bytearray)[end-BigInt(namebytes):-1:end-BigInt(namebytes)-BigInt(filebytes)+1]
        f=open(name,"w")
        try
            write(f,filedata)
        finally
            close(f)
        end
    end
end

#   ;
function terminate()
    println("\nProgram finished!")
end


text:
# Test for the correct answer
@test abs(getconfidencelevel(m) - 0.53722) <= 1e-3

# Check that it doesn't work for other sets and then reset the set
setuncertaintyset(m, :PNorm, 1)
@test_throws ErrorException abs(getconfidencelevel(m) - 0.53722) <= 1e-3
setuncertaintyset(m, :Ellipsoid)


text:
using Test
using Cribbage

using Cribbage: score_fifteens, score_runs, score_pairs, score_jack, score_flush, _Hand

@testset "card" begin
    c = Card("A♠")
    @test c == Card("A♠")
    @test c == Card('A', '♠')
    @test c == A♠
    @test rank(c) == Rank('A')
    @test rank(c) == A
    @test suit(c) == Suit('♠')
    @test suit(c) == ♠
    @test mask(c) == 1
    @test value(c) == 1
    @test A♠ < 2♠
    @test A♠ < 4♡
    @test Q♠ < K♠
    @test K♠ > Q♠
    @test K♠ > Q♠
end

@testset "deck" begin
    @test length(deck()) == 52
    @test length(unique!(deck())) == 52
    @test length(DECK) == 52
end

@testset "hand" begin
    h = rand(_Hand)
    @test length(hand(h)) == 4

    s = rand(Hand)
    @test length(hand(s)) == 4
    @test cut(s) isa Card
end

@testset "scores" begin
    s = Hand([7♣, 8♣, 7♡, 8♡], 7♠)
    @test score_fifteens(s) == 12
    @test score_pairs(s) == 8
    @test score(s) == 20

    s = Hand([T♣, 4♣, A♡, 8♡], 7♠)
    @test score_fifteens(s) == 4
    @test score(s) == 4

    s = Hand([7♣, 8♣, 7♡, 8♡], 9♠)
    @test score_fifteens(s) == 8
    @test score_pairs(s) == 4
    @test score_runs(s) == 12
    @test score(s) == 24

    s = Hand([2♣, 3♣, 4♣, 5♣], 6♣)
    @test score_flush(s) == 5
    @test score_runs(s) == 5
    @test score_fifteens(s) == 4
    @test score(s) == 14

    s = Hand([2♣, 3♣, 4♣, 5♣], 6♡)
    @test score_flush(s) == 4
    @test score_runs(s) == 5
    @test score_fifteens(s) == 4
    @test score(s) == 13

    s = Hand([2♣, 3♣, 4♡, 5♣], 6♡)
    @test score_flush(s) == 0
    @test score_runs(s) == 5
    @test score_fifteens(s) == 4
    @test score(s) == 9

    s = Hand([2♣, 3♣, 4♣, 5♣], 4♡)
    @test score_flush(s) == 4
    @test score_runs(s) == 8
    @test score_fifteens(s) == 2
    @test score_pairs(s) == 2
    @test score(s) == 16

    s = Hand([T♣, J♣, Q♣, K♣], A♡)
    @test score_runs(s) == 5
    @test score_flush(s) == 4
    @test score(s) == 9

    s = Hand([Q♣, J♣, Q♣, K♣], A♡)
    @test score_runs(s) == 8
    @test score_flush(s) == 4
    @test score_pairs(s) == 2
    @test score(s) == 14

    s = Hand([Q♣, J♣, A♣, K♣], A♡)
    @test score_runs(s) == 8
    @test score_flush(s) == 4
    @test score_pairs(s) == 2
    @test score(s) == 14

    # Best hand
    s = Hand([5♣, 5♢, 5♠, J♡], 5♡)
    score(s) == 29
    @test score(s) == 29

    # Second-best hands
    s = Hand([5♣, 5♢, 5♠, 5♡], J♡)
    @test score(s) == 28

    s = Hand([5♣, 5♢, T♡, 5♡], 5♠)
    @test score(s) == 28
end


text:
module SuperTuxGen

using Base: write

# Use Julia's internal syntax trees as the representation for the SuperTux files
include("write_scheme.jl")

end # module


text:
# Command: Add Cookie
"""
    cookie!(session::Session, cookie::Cookie)::Nothing

Add Cookie.
"""
function cookie!(session::Session, cookie::Cookie)::Nothing
    @unpack addr, id = session
    response = HTTP.post(
        "$addr/session/$id/cookie",
        [("Content-Type" => "application/json")],
        JSON3.write("cookie" => cookie),
    )
    @assert response.status == 200
end


text:
################################################################################
# Drone Example
################################################################################
# using Algames
using StaticArrays
using LinearAlgebra
using MeshCat

vis  = Visualizer()
open(vis)
T = Float64

# Define the dynamics of the system
p = 4 # Number of players
model  = QuadrotorGame(p=p) # game with 3 players with double integrator dynamics in 3D
n = model.n
m = model.m

# Define the horizon of the problem
N = 50 # N time steps
dt = 0.40 # each step lasts 0.1 second
dt = 0.16 # each step lasts 0.1 second
probsize = ProblemSize(N,model) # Structure holding the relevant sizes of the problem

# Define the objective of each player
# We use a LQR cost
Q = [Diagonal(5*SVector{model.ni[i],T}([1*[0,0.3,0.3,1,1,1]; 5*[1,0.3,0.3,1,1,1]])) for i=1:p] # Quadratic state cost
R = [Diagonal(0.01*ones(SVector{model.mi[i],T})) for i=1:p] # Quadratic control cost
# Desrired state
xf = [SVector{model.ni[1],T}([1.0,-0.5,+0.3,0,0,0, 0.3,0,0, 0,0,0]),
	  SVector{model.ni[2],T}([1.5,-0.7,-0.3,0,0,0, 0.3,0,0, 0,0,0]),
	  SVector{model.ni[3],T}([0.8,-0.0,+0.5,0,0,0, 0.3,0,0, 0,0,0]),
      SVector{model.ni[3],T}([0.8,-0.0,-0.5,0,0,0, 0.3,0,0, 0,0,0]),
      ]
xf = xf[1:p]
# Desired control
uf = [- model.mass * model.gravity[end] / 4 / model.kf * ones(SVector{model.mi[i],T}) for i=1:p]
# Objectives of the game
game_obj = GameObjective(Q,R,xf,uf,N,model)
radius = 0.2*ones(p)
μ = 20.0*ones(p)
add_collision_cost!(game_obj, radius, μ)

# Define the constraints that each player must respect
game_con = GameConstraintValues(probsize)
# Add collision avoidance
radius = 0.08
add_spherical_collision_avoidance!(game_con, radius)
# Add wall constraint
room_walls = [
    Wall3D([-1.50, -1.00, -1.00], [1.50, -1.00, -1.00], [1.50,  1.00, -1.00], [0.00,  0.00, -1.00]),
    Wall3D([-1.50, -1.00,  1.00], [1.50, -1.00,  1.00], [1.50,  1.00,  1.00], [0.00,  0.00,  1.00]),
    Wall3D([-1.50, -1.00, -1.00], [1.50, -1.00, -1.00], [1.50, -1.00,  1.00], [0.00, -1.00,  0.00]),
    Wall3D([-1.50,  1.00, -1.00], [1.50,  1.00, -1.00], [1.50,  1.00,  1.00], [0.00,  1.00,  0.00]),
    ]
door_cylinders = [
	CylinderWall([0.00, -1.00, -1.00], :z, 2.00, 0.95),
	CylinderWall([0.00,  1.00, -1.00], :z, 2.00, 0.95),
	CylinderWall([0.00, -1.00, -1.00], :y, 2.00, 0.95),
	CylinderWall([0.00, -1.00,  1.00], :y, 2.00, 0.95),
	]
add_wall_constraint!(game_con, room_walls)
add_wall_constraint!(game_con, door_cylinders)

build_wall!(vis, room_walls, α=0.05, name=:Room)
build_cylinder!(vis, door_cylinders, radius=radius, α=0.05, name=:Door1)
build_door_frame!(vis, size=10, α=1.0)

# Define the initial state of the system
x0 = [
	[-1.0, -0.4, -0.6, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
	[-0.8,  0.4,  0.4, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
	[-1.4,  0.2,  0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
	[-1.6, -0.3, -0.3, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
	]
x0 = SVector{model.n,T}(reshape(vcat([x0[i]' for i=1:p]...), model.n))

# Define the Options of the solver
opts = Options()
# Define the game problem
prob = GameProblem(N,dt,x0,model,opts,game_obj,game_con)

# Solve the problem
@time newton_solve!(prob)

# Visualize the Results
using Plots
Algames.plot!(prob.model, prob.pdtraj.pr)
Algames.plot!(prob.stats)

build_traj!(vis, model, prob.pdtraj.pr, α=1.0, name=:Traj_orange)
build_xf!(vis, model, xf, α=1.0, name=:Xf)

anim = visualize_drone!(vis, model, prob.pdtraj.pr)

plot(hcat([Algames.state(prob.pdtraj.pr[i]) for i=1:N]...)')
plot(hcat([Algames.control(prob.pdtraj.pr[i]) for i=1:N]...)')

# filename = "4_drones_door"
# MeshCat.convert_frames_to_video(
#     "/home/simon/Downloads/$filename.tar",
#     "/home/simon/Documents/$filename.mp4", overwrite=true)
#
# convert_video_to_gif(
#     "/home/simon/Documents/$filename.mp4",
#     "/home/simon/Documents/$filename.gif", overwrite=true)
#
# using FFMPEG
#


text:
using ValidatedNumerics, ValidatedNumerics.RootFinding
using ForwardDiff

using FactCheck

const D = RootFinding.derivative

include("wilkinson.jl")

function generate_wilkinson(n)#, T=BigFloat)   # SLOW

    p = poly(collect(1:n))
    p_prime = polyder(p)

    coeffs = p.a #[convert(T, coeff) for coeff in p.a]
    # coeffs_prime = [convert(T, coeff) for coeff in p_prime.a]

    function f(x)
        total = coeffs[1]
        for i in 2:length(coeffs)
            total += coeffs[i]*x^(i-1)
        end
        total
    end

    f
end

setprecision(Interval, Float64)
float_pi = @interval(pi)

setprecision(Interval, 10000)
big_pi = @interval(pi)
# Using precision "only" 256 leads to overestimation of the true roots for `cos`
# i.e the Newton method gives more accurate results!

half_pi = big_pi / 2
three_halves_pi = 3*big_pi/2

# Format:  (function, derivative, lower_bound, upper_bound, [true_roots])
function_list = [
                    (sin, cos,    -5,  5,    [ -big_pi, @interval(0), big_pi ] ) ,
                    (cos, x->D(cos, x), -7.5, 7.5, [ -three_halves_pi, -half_pi, half_pi, three_halves_pi ] ),
                    (W₃, x->D(W₃, x),   -10, 10,   [ @interval(1), @interval(2), @interval(3) ] ),
                    (W₇, x->D(W₇, x),   -10, 10,   [ @interval(i) for i in 1:7 ] ),
                    (x->exp(x)-2, y->D(x->exp(x),y),  -20, 20,  [log(@biginterval(2))] ),
                    (x->asin(sin(x)) - 0.1, y->D(x->asin(sin(x)),y), 0, 1, [@biginterval(0.1)])
                ]

facts("Testing root finding") do

    for rounding_type in (:wide, :narrow)
        context("Interval rounding: $rounding_type") do
            setrounding(Interval, rounding_type)

            for prec in ( (BigFloat,53), (BigFloat,256), (Float64,64) )
                context("Precision: $prec") do
                    setprecision(Interval, prec)

                    for method in (newton, krawczyk)
                        context("Method $method") do

                            for func in function_list

                                f, f_prime, a_lower, a_upper, true_roots = func
                                a = @interval(a_lower, a_upper)

                                context("Function $f; interval $a") do

                                    for autodiff in (false, true)

                                        context("With autodiff=$autodiff") do

                                            if autodiff
                                                roots = method(f, a)
                                            else
                                                roots = method(f, f_prime, a)
                                            end

                                            @fact length(roots) --> length(true_roots)

                                            for i in 1:length(roots)
                                                root = roots[i]

                                                @fact isa(root, Root{prec[1]}) --> true
                                                @fact is_unique(root) --> true
                                                @fact true_roots[i] ⊆ root.interval --> true
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

f(x) = x^2 - 2

roots = newton(f, @interval(10, 11))

facts() do
    @fact length(roots) --> 0
end

setprecision(Interval, Float64)

facts("find_roots tests") do
    f(x) = x^2 - 2

    roots = find_roots_midpoint(f, -5, 5)
    @fact length(roots) --> 3
    @fact length(roots[1]) --> 2

    roots = find_roots(f, -5, 5)
    @fact length(roots) --> 2

    setprecision(Interval, 256)

    for method in (newton, krawczyk)
        new_roots = method(f, roots)
        @fact length(new_roots) == length(roots) --> true

        for (root1, root2) in zip(new_roots, roots)
            @fact root1 ⊆ root2 --> true
        end
    end

end

facts("Multiple roots") do
    setprecision(Interval, Float64)
    let
        f(x) = (x-1) * (x^2 - 2)^3 * (x^3 - 2)^4

        roots = newton(f, -5..5.1, maxlevel=1000)

        @fact length(roots) --> 4
        @fact roots[1].status --> :unknown
        @fact roots[1].interval --> Interval(-1.4142135623730954, -1.414213562373095)
        @fact roots[3].interval --> Interval(1.259921049894873, 1.2599210498948734)

    end
end

# Example of a function with a double root at 0 from Burden & Faires, 9th ed, p.84
# exp(x) - x - 1


text:
using Metaheuristics

if VERSION < v"0.7.0"
    using Base.Test
    srand(31415926534)
else
    using Test
end

# write your own tests here
@testset "Multi objective" begin

    function run_methods(problem)

        ff, bounds, pf = Metaheuristics.TestProblems.get_problem(problem)
        D = size(bounds, 2)

        # number of function evaluations
        f_calls = 0

        f(x) = begin
            f_calls += 1
            ff(x)
        end

        options = Options( seed = 1, iterations=10000, f_calls_limit = 25000)
        nobjectives = length(pf[1].f)
        npartitions = nobjectives == 2 ? 100 : 12

        methods = [
                SMS_EMOA(N = 50, n_samples=500, options=options),
                NSGA2(options=options),
                MOEAD_DE(gen_ref_dirs(nobjectives, npartitions), options=Options( seed = 1, iterations = 500)),
                NSGA3(options=options),
              ]

        for method in methods
            f_calls = 0
            result = ( optimize(f, bounds, method) )
            show(IOBuffer(), "text/html", result)
            show(IOBuffer(), "text/plain", result.population)
            show(IOBuffer(), "text/html", result.population)
            show(IOBuffer(), result.population[1])
            @test Metaheuristics.PerformanceIndicators.igd(result.population, pf) <= 0.2
            @test Metaheuristics.PerformanceIndicators.spacing(result) < 0.2
            @test Metaheuristics.PerformanceIndicators.covering(pf, result.population) <= 1.0
            @test Metaheuristics.PerformanceIndicators.covering(result, result) ≈ 0.0

            # number of function evaluations should be reported correctly
            @test f_calls == result.f_calls

            # test obtaining non-dominated solutions
            pf1 = pareto_front(result)
            pf2 = pareto_front(result.population)

            @test size(pf1, 1) == size(pf2,1) &&
                  Metaheuristics.PerformanceIndicators.igd(pf1, pf2) ≈ 0.0
        end
    end

    for problem in [:ZDT3, :DTLZ2, :MTP]
        run_methods(problem)
    end

end


text:
parse_float(x) = parse(Float64, x)

function run_lblrtm(work_dir::String, ν_s::Float64, ν_e::Float64)::Tuple{Array{Float64,1},Array{Float64,1}}
    ν_s -= 25
    ν_e += 25

    @info "Entering $(pwd()) ..."
    cd(work_dir)
    @info "Entered."

    @info "Checking system links..."
    requirements = ["lblrtm", "TAPE3"]
    for link in requirements
        @assert(islink(link) || isfile(link), "$link does not exist.")
    end
    @info "System links OK."

    @info "Removing last run's files..."
    rm("TAPE5", force=true)
    rm("TAPE6", force=true)
    rm("TAPE7", force=true)
    rm("TAPE8", force=true)
    rm("TAPE10", force=true)
    @info "Removed."

    IHIRAC = 1
    ILBLF4 = 2
    ICNTNM = 1
    IAERSL = 1
    IEMIT = 1
    ISCAN = 0
    IFILTR = 0
    IPLOT = 0
    ITEST = 0
    IATM = 1
    IMRG = 0
    ILAS = 0
    IOD = 1
    IXSECT = 0
    NPTS = 80
    ISOTPL = 0
    IBRD = 0

    SAMPLE = 4
    DVSET = 0.0 # Spectral resolution for calculation
    ALFAL0 = 0.04
    AVMASS = 36.0
    DPTMIN = 0.0002
    DPTFAC = 0.001
    IFNFLG = 0
    DVOUT = 0.1 # Spectral resolution for output

    TBOUND = 300.0 # Boundary temperature
    SREMIS = 1.0 # Emissivity coefficient - Rand 0

    MODEL = 6 # U.S.Standard
    ITYPE = 2
    IBMAX = 0
    ZERO = 0
    NOPRNT = 0

    H1 = 100.0
    H2 = 0.0
    ANGLE = 180.0
    @info "Writing TAPE5..."
    open("TAPE5", "w") do tape5
        write(tape5, "\$ LBLRTM INPUT\n") # Record 1.1
        write(tape5, " HI=$IHIRAC F4=$ILBLF4 CN=$ICNTNM AE=$IAERSL EM=$IEMIT SC=$ISCAN FI=$IFILTR PL=$IPLOT TS=$ITEST AM=$IATM MG=$IMRG LA=$ILAS OD=$IOD XS=$IXSECT   00   00\n") # Record 1.2
        write(tape5, @sprintf("%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%10.3f%5d%10.3f\n", ν_s, ν_e, SAMPLE, DVSET, ALFAL0, AVMASS, DPTMIN, DPTFAC, IFNFLG, DVOUT)) # Record 1.3
        write(tape5, @sprintf("%10.3f%10.3f\n", TBOUND, SREMIS)) # Record 1.4
        write(tape5, @sprintf("%5d%5d%5d%5d%5d\n", MODEL, ITYPE, IBMAX, ZERO, NOPRNT)) # Record 3.1
        write(tape5, @sprintf("%10.3f%10.3f%10.3f\n", H1, H2, ANGLE)) # Record 3.2
        write(tape5, "\n") # Record 3.3A
        write(tape5, "\n") # Record 4.1
        write(tape5, "-1.\n")
        write(tape5, @sprintf(
"""
\$ Transfer to ASCII plotting data
 HI=0 F4=0 CN=0 AE=0 EM=0 SC=0 FI=0 PL=1 TS=0 AM=0 MG=0 LA=0 MS=0 XS=0    0    0
# Plot title not used
%10.3f%10.3f   10.2000  100.0000    5    0   12    0     1.000 0  0    0
    0.0000    1.2000    7.0200    0.2000    4    0    1    1    0    0 0    3 27
%10.3f%10.3f   10.2000  100.0000    5    0   12    0     1.000 0  0    0
    0.0000    1.2000    7.0200    0.2000    4    0    1    1    1    0 0    3 28
-1.
""", ν_s, ν_e, ν_s, ν_e))
    end
    @info "TAPE5 OK."

    # Run LBLRTM
    @info "Running LBLRTM..."
    @time run(`./lblrtm`)
    @info "LBLRTM completed."

    # Collect result
    @info "Collecting result..."
    wavenumber = []
    brightness_temperature = []
    open("TAPE28") do io
        see_header = false
        see_data = false
        for line in eachline(io)
            if !see_header
                if findfirst("WAVENUMBER", line) ≠ nothing
                    see_header = true
                end
            elseif !see_data
                see_data = true
            else
                wb, bt = map(parse_float, split(strip(line), r"\s+"))
                append!(wavenumber, wb)
                append!(brightness_temperature, bt)
            end
        end
    end
    @info "Results collected."

    return (wavenumber, brightness_temperature)
end


text:
using PoincareInvariants.SecondFinDiffPlans

@safetestset "getpoints, getpointspec and getpointnum" begin
    using PoincareInvariants

    @test getpointspec(50, SecondFinDiffPlan) == (9, 9)
    @test getpointspec(87, SecondFinDiffPlan) == (11, 11)
    @test getpointspec((53, 42), SecondFinDiffPlan) == (53, 43)

    @test getpointnum((3, 5)) == 15
    @test getpointnum((423789, 326)) == 423789 * 326
    @test getpointnum((123, 908343)) == 123 * 908343

    @inferred Matrix{Float64} getpoints((x, y) -> (x, x+y), Float64, (7, 5), SecondFinDiffPlan)
    @inferred Vector{Float32} getpoints((x, y) -> sin(x), Float32, (9, 9), SecondFinDiffPlan)

    @test getpoints((x, y) -> (x, y), Float64, (5, 3), SecondFinDiffPlan) ≈ [
        0    0  ;
        0    0.5;
        0    1  ;
        0.25 0  ;
        0.25 0.5;
        0.25 1  ;
        0.5  0  ;
        0.5  0.5;
        0.5  1  ;
        0.75 0  ;
        0.75 0.5;
        0.75 1  ;
        1.0  0  ;
        1.0  0.5;
        1.0  1  ]

    @test getpoints((x, y) -> (4x, 6y), Float64, (5, 7), SecondFinDiffPlan) ≈ [
        0 0;
        0 1;
        0 2;
        0 3;
        0 4;
        0 5;
        0 6;
        1 0;
        1 1;
        1 2;
        1 3;
        1 4;
        1 5;
        1 6;
        2 0;
        2 1;
        2 2;
        2 3;
        2 4;
        2 5;
        2 6;
        3 0;
        3 1;
        3 2;
        3 3;
        3 4;
        3 5;
        3 6;
        4 0;
        4 1;
        4 2;
        4 3;
        4 4;
        4 5;
        4 6]

    @test getpoints((x, y) -> y, Float64, (11, 11), SecondFinDiffPlan) ≈ Float64[
        0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
        0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
        0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
        0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
        0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
        0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
        0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
        0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
        0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
        0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
        0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0,
    ]

    ps = getpointspec(6342, SecondFinDiffPlan)
    testpnts = getpoints((x, y) -> (x, y), Float64, ps, SecondFinDiffPlan)
    pnts4 = getpoints(Float64, ps, SecondFinDiffPlan) do x, y
        x, y, x * y, x + y
    end
    @test pnts4 isa Matrix{Float64}
    @test size(pnts4) == (getpointnum(ps), 4)

    @test pnts4[:, 1] ≈ testpnts[:, 1]
    @test pnts4[:, 2] ≈ testpnts[:, 2]
    @test pnts4[:, 3] ≈ testpnts[:, 1] .* testpnts[:, 2]
    @test pnts4[:, 4] ≈ testpnts[:, 1] .+ testpnts[:, 2]
end

@safetestset "Differentiation" begin
    using ..SecondFinDiffPlans: getpoints, SecondFinDiffPlan, differentiate

    # Finite difference derivatives should be exact for a quadratic
    let f(x, y)  = 1 + 0.75*x + y - 0.5*x*x - y*y - 1.5*x*y
        fx(x, y) = 0.75 - x - 1.5*y
        fy(x, y) = 1 - 2*y - 1.5*x

        for (nx, ny, maxerr) in [(3, 3, 50), (21, 33, 500), (331, 123, 5000)]
            vals = getpoints(f, Float64, (nx, ny), SecondFinDiffPlan)

            @test all((ix, iy) for ix in 1:nx, iy in 1:ny) do (ix, iy)
                ∂x, ∂y = differentiate(vals, ix, iy, (nx, ny))
                x ,  y = (ix - 1) / (nx - 1), (iy - 1) / (ny - 1)
                tx, ty = fx(x, y), fy(x, y)
                return abs(tx - ∂x) / eps() < maxerr && abs(ty - ∂y) / eps() < maxerr
            end
        end
    end

    let f(x, y) = sinpi(x) * exp(y)
        fx(x, y) = π * cospi(x) * exp(y)
        fy(x, y) = sinpi(x) * exp(y)

        nx, ny = 1321, 1235
        vals = getpoints(f, Float64, (nx, ny), SecondFinDiffPlan)
        maxerr = 5e-5

        @test all((ix, iy) for ix in 1:nx, iy in 1:ny) do (ix, iy)
            ∂x, ∂y = differentiate(vals, ix, iy, (nx, ny))
            x ,  y = (ix - 1) / (nx - 1), (iy - 1) / (ny - 1)
            tx, ty = fx(x, y), fy(x, y)
            return abs(tx - ∂x) < maxerr && abs(ty - ∂y) < maxerr
        end
    end
end

@safetestset "Simpson weights" begin
    using ..SecondFinDiffPlans: getpoints, SecondFinDiffPlan, getsimpweight

    @test [getsimpweight(Float64, ix, iy, (3, 5)) * 9 * 4 * 2 for iy in 1:5, ix in 1:3] ≈ [
        1  4 1;
        4 16 4;
        2  8 2;
        4 16 4;
        1  4 1]

    @test [getsimpweight(Float64, ix, iy, (5, 7)) for iy in 1:7, ix in 1:5] ≈ [
        1  4  2  4  1;
        4 16  8 16  4;
        2  8  4  8  2;
        4 16  8 16  4;
        2  8  4  8  2;
        4 16  8 16  4;
        1  4  2  4  1
    ] ./ (9 * 4 * 6)

    let f(x, y) = 1 + 2x + 3y + 4x^2 + 5x*y + 6y^2
        nx, ny = (3, 3)
        vals = getpoints(f, Float64, (nx, ny), SecondFinDiffPlan)
        weights = [getsimpweight(Float64, ix, iy, (nx, ny)) for iy in 1:ny, ix in 1:nx] |> vec
        @test sum(weights .* vals) ≈ 97 / 12 atol=10eps()
    end
end

@safetestset "compute!" begin
    using ..SecondFinDiffPlans: getpoints, SecondFinDiffPlan, getsimpweight
    using PoincareInvariants
    using LinearAlgebra: dot

    f(x, y) = [
         1 +  2*x +  3*y +  4*x^2 +  5*x*y +  6*y^2,
         7 +  8*x +  9*y + 10*x^2 + 11*x*y + 12*y^2,
        13 + 14*x + 15*y + 16*x^2 + 17*x*y + 18*y^2,
        19 + 20*x + 21*y + 22*x^2 + 23*x*y + 24*y^2]

    fx(x, y) = [
         2 +  8*x +  5*y,
         8 + 20*x + 11*y,
        14 + 32*x + 17*y,
        20 + 44*x + 23*y]

    fy(x, y) = [
         3 +  5*x + 12*y,
         9 + 11*x + 24*y,
        15 + 17*x + 36*y,
        21 + 23*x + 48*y]

    ω(z, ::Any, ::Any) = [0 0 z[1] z[2]; 0 0 z[3] z[4]; -z[1] -z[3] 0 0; -z[2] -z[4] 0 0]
    nx, ny = 11, 17
    integrand = getpoints(Float64, (nx, ny), SecondFinDiffPlan) do x, y
        dot(fy(x, y), ω(f(x, y), 0, nothing), fx(x, y))
    end

    ws = [getsimpweight(Float64, ix, iy, (nx, ny)) for iy in 1:ny, ix in 1:nx] |> vec
    testI = dot(ws, integrand)

    pinv = SecondPoincareInvariant{Float64, 4}(ω, (nx, ny), SecondFinDiffPlan)

    @test pinv.plan isa SecondFinDiffPlan

    points = getpoints(f, Float64, (nx, ny), SecondFinDiffPlan)
    I = compute!(pinv, points, 0, nothing)

    @test I ≈ testI atol=5e-11
end


text:
# LibTOAST.jl: interface to the TOAST++ library
# Copyright (C) 2019 Samuel Powell

@testset "CoeffTypes" begin

  vtx,  = data(mesh2D)
  nprm = sin.(5.0.*vtx[:,1]./15) .+ cos.(4.0.*vtx[:,2]./15)

  nx, ny = 250, 250
  pixelmap = PixelMap(mesh2D, (nx,ny))

  @test all(zero(NodalCoeff, mesh2D) .== 0)
  @test all(zero(SolutionCoeff, pixelmap) .== 0)
  @test all(zero(RasterCoeff, pixelmap) .== 0)
  @test all(zero(IntermediateCoeff, pixelmap) .== 0)

  @test all(one(NodalCoeff, mesh2D) .== 1)
  @test all(one(SolutionCoeff, pixelmap) .== 1)
  @test all(one(RasterCoeff, pixelmap) .== 1)
  @test all(one(IntermediateCoeff, pixelmap) .== 1)

  nc0 = NodalCoeff(mesh2D, nprm)

  sc1 = SolutionCoeff(pixelmap, nc0)
  gc1 = IntermediateCoeff(pixelmap, nc0)
  bc1 = RasterCoeff(pixelmap, nc0)

  #nc1 = NodalCoeff(bc1)
  #nc2 = NodalCoeff(sc1)
  #nc3 = NodalCoeff(gc1)

  sc2 = SolutionCoeff(gc1)
  sc3 = SolutionCoeff(bc1)

  gc2 = IntermediateCoeff(sc1)
  gc3 = IntermediateCoeff(bc1)

  bc2 = RasterCoeff(sc1)
  bc3 = RasterCoeff(gc1)

  maperr(c1, c2) = 100*norm(c2.-c1)./norm(c1)
  maptol = .5

  # Tset functions on the mesh
  #@test maperr(nc0,nc1) < maptol
  #@test maperr(nc0,nc2) < maptol
  #@test maperr(nc0,nc3) < maptol

  # Test functions in the raster
  @test maperr(bc1,bc2) < maptol
  @test maperr(bc1,bc3) < maptol

  # Test functions in the solution
  @test maperr(sc1,sc2) < maptol
  @test maperr(sc1,sc3) < maptol

  # Test functions in the intermediate
  @test maperr(gc1,gc2) < maptol
  @test maperr(gc1,gc3) < maptol

end


text:
using Hecke
bla = normpath(joinpath(dirname(pathof(Hecke)), "..", "docs", "make.jl"))
local_build = true
include(bla)

function _open(filename; wait = false)
    @static if Sys.isapple()
        run(`open $(filename)`; wait = wait)
    elseif Sys.islinux() || Sys.isbsd()
        run(`xdg-open $(filename)`; wait = wait)
    elseif Sys.iswindows()
        cmd = get(ENV, "COMSPEC", "cmd.exe")
        run(`$(cmd) /c start $(filename)`; wait = wait)
    else
        @warn("Opening files the default application is not supported on this OS.",
              KERNEL = Sys.KERNEL)
    end
end

if !Hecke.html_build[]
  cd(normpath(joinpath(dirname(pathof(Hecke)), "..", "docs"))) do
    run(`mkdocs build --no-directory-urls`, wait = true)
  end

  bla = normpath(joinpath(dirname(pathof(Hecke)), "..", "docs", "site", "index.html"))
  _open(bla)
else
  bla = normpath(joinpath(dirname(pathof(Hecke)), "..", "docs", "build", "index.html"))
  _open(bla)
end


text:
export rectangulararray

#
# switch to generator expression:
# f(i) for i in 1:n

"""
    rectangulararray(x[,y[,z]])
    rectangulararray(<keyword arguments>)

Used with `transform` to create rectangular arrays of elements.

`rectangulararray` returns a iterable object which returns a sequence of
transformation matrices, which when passed to `transform` will produce
a rectangular array of elements.  `x`, `y`, and `z` are iterable objects
which specify the offsets along the x, y, and z axis.  Unspecified
arguments are zero.
"""
function rectangulararray(x,y=0.0,z=0.0)
  tm = txyz(0,0,0)
  ch = Channel{typeof(tm)}(1)
  @async begin
    for a in x
      tm[1,4] = a
      for b in y
        tm[2,4] = b
        for c in z
          tm[3,4] = c
          put!(ch,copy(tm))
        end
      end
    end
    close(ch)
  end
  return ch
end
rectangulararray(;x=0.0, y=0.0, z=0.0) = rectangulararray(x,y,z)


text:
#Define main function for the ground state optimization
#and import the necessary routines

#Main function for ground state optimization
function GS_optimization(GsHp::GS_hp, dim::Int64)

    println("# Starting ground state optimization for the Heisenberg Hamiltonian, with $nspins = $L_x x $L_y spins and α=$alpha")
    println("# Number of workers = $(nprocs())")
    println("# The following hyper-parameters are used:")
    println("# Number of sweeps = $(GsHp.samples),  iteration step = $(GsHp.n_iter),  learning rate = $(GsHp.step_size)")

    #If the number of workers exceeds the number of variational parameters,
    #only a number of workers equals to the number of variational parameters will be used.
    #The exceeding workers will be removed.
    if nprocs() > dim
        println("WARNING: number of workers exceeds number of variational parameters. Only $dim workers will be used.")
        diff = nprocs() - dim
        rmprocs(nprocs()-diff+1:nprocs())
    end

    #Store number of workers for later use
    NW = Int64[]; push!(NW, NumberWorkers())

    #Define matrix storing energy and wavefunction gradients for each sampled state
    #@everywhere ensure a copy of ENWfDer is created for each worker
    @everywhere ENWfDer = Array{Complex{Float64}}(undef,@getfrom 1 Nsample,dim + 1)

    #Initialize network parameters to random numbers
    RBM_Par = init_RBM_par(alpha,nspins)
    #Sort network parameters to W_RBM independent parameters and vice versa
    W_RBM = Array{Complex{Float64}}(undef,dim)
    ParToW(RBM_Par,W_RBM); WToPar(RBM_Par,W_RBM)

    #Set value of J_x, J_y
    @everywhere J = [J_x; J_y]
################################################################################
#############    Start main iteration  #########################################
    for NIter = 1:GsHp.n_iter

        #Sampling of states
        #At the end of the loop, the energy sampled are returned in EnergyStates
        #The derivatives of the wavefunction are stored in ENWFDer
        EnergyStates = @distributed (vcat) for i=1:NW[1]
              sampling(RBM_Par,J,i)
         end

         #Evaluation entries of S-matrix stored in the vector SVector
        SVector = @distributed (+) for i=1:NW[1]
            S_gen(dim)./NW[1]
         end

         #Sort entries SVector into the S matrix
         SMatrix = Gen_S(SVector)

         #Regularization of the S-matrix with a small shift of the diagonal element of S
         SMatrix = SMatrix .+ Matrix{Float64}(I,dim,dim).*max(100*(0.9)^NIter,0.0001)

         #Inversion of the S-matrix with normal inversion routine (not pseudo-inverse)
         InvS = inv(SMatrix)

         #Store the total energy for each state
         energy = EnergyStates[:]
         #Calculate the mean value of the energy and energy variance and store them
         mean_E = mean(energy); var_E = var(energy)
         push!(Energy_,real(mean_E)); push!(Variance_,norm(var_E))

         ##calculate the derivative of the energy respect to the variational parameters and assign them to the array dF
         EGrad = @distributed (+) for i=1:NW[1]
            Egrad(dim)./NW[1]
         end

         #Update variational parameters according to SR optimization rule
         for i=1:dim
             W_RBM[i] = W_RBM[i] - GsHp.step_size*(transpose(InvS[i,:])*EGrad)
         end

         #Update neural network parameters after optimization
         for i in procs()
             fetch(@spawnat i WToPar(RBM_Par,W_RBM))
         end

         #print value of energy every 10 steps
         if NIter % 10 == false
             println("Iteration step #$NIter")
             println("Energy: $(real(mean_E)) +- $(norm(var_E))")
             println(" ")
         end

     end

###########################   End main iteration    ############################
################################################################################
    println("# ground state optimization completed")
    println(" ")

    #Store the final independent variational parameters (real and imaginary parts)
    writedlm("W_RBM_$(nspins)_$(alpha)_real.jl", real(W_RBM)); writedlm("W_RBM_$(nspins)_$(alpha)_imag.jl", imag(W_RBM))

end

max_stars_repo_name: GiamFbn/MyULTRAFAST

max_issues_repo_name: GiamFbn/MyULTRAFAST

max_forks_repo_name: GiamFbn/MyULTRAFAST


text:
module IndexCore

using StaticRanges

using Base: OneTo,
            to_index,
            tail,
            axes,
            broadcasted,
            AbstractCartesianIndex,
            @propagate_inbounds,
            @_propagate_inbounds_meta

using StaticRanges:
    can_set_length,
    can_set_first,
    can_set_last,
    ForwardOrdering,
    ReverseOrdering,
    ContinuousTrait,
    Continuous,
    DiscreteTrait,
    Discrete,
    similar_type,
    OneToRange,
    AbstractStepRange,
    AbstractStepRangeLen,
    AbstractLinRange,
    StaticUnitRange,
    Length,
    Size

export
    # Types
    AbstractIndex,
    AbstractLengthCheck,
    AllUniqueTrait,
    AllUnique,
    Index,
    LengthCheckedTrait,
    LengthChecked,
    LengthNotCheckedTrait,
    LengthNotChecked,
    NotUniqueTrait,
    SimpleIndex,
    Uniqueness,
    UnkownUnique,
    UnkownUniqueTrait,
    # methods
    append_axes,
    append_axes!,
    append_axis,
    append_axis!,
    append_keys,
    append_values,
    axes2length,
    axes2size,
    cat_axes,
    cat_axis,
    cat_names,
    cat_keys,
    cat_values,
    check_dimensions,
    check_index_length,
    check_index_uniqueness,
    check_index_params,
    combine_index,
    combine_indices,
    combine_keys,
    combine_names,
    combine_values,
    covcor_axes,
    dimnames,
    drop_axes,
    filter_axes,
    hcat_axes,
    index_by,
    indices,
    inverse_axes,
    keys_type,
    matmul_axes,
    permute_axes,
    pop_index,
    pop_index!,
    popfirst_index,
    popfirst_index!,
    push_index!,
    pushfirst_index!,
    reduce_axes,
    reduce_axis,
    reindex,
    reshape_axes,
    reshape_axes!,
    resize_axis!,
    to_dims,
    unsafe_reindex,
    values_type,
    vcat_axes,
    dimn

include("abstractindex.jl")
include("param_checks.jl")
include("index.jl")
#include("subindex.jl")    # rewrite, doc, ex, tests
include("simpleindex.jl")
include("iterate.jl")      # tests
include("names.jl")
include("size.jl")         # doc, ex, tests
include("reindex.jl")
include("to_index.jl")
include("to_indices.jl")
#include("to_linear.jl")   # rewrite, doc, ex, tests
include("checkindex.jl")   # tests
include("getindex.jl")     # tests
include("combine.jl")      # ex
include("drop_axes.jl")    # tests
include("reduce_axes.jl")  # tests
include("reshape_axes.jl") # ex, tests
include("inverse_axes.jl")
include("covcor_axes.jl")
include("permute_axes.jl") # tests
include("matmul_axes.jl")  # tests
include("filter_axes.jl")  # tests
include("cat_axes.jl")     # tests
include("append_axes.jl")  # tests
include("make_unique.jl")  # doc, ex, tests, TODOs
include("promotion.jl")
include("push_pop.jl")     # tests
#include("prmote_shape.jl")

end


text:
# Subplot geometries

abstract SubplotGeometry <: Gadfly.GeometryElement

# Adding elements to subplots in a generic way.

function add_subplot_element(subplot::SubplotGeometry, arg::Base.Callable)
    add_subplot_element(subplot, arg())
end

function add_subplot_element(subplot::SubplotGeometry, arg::Gadfly.Layer)
    push!(subplot.layers, arg)
end

function add_subplot_element(subplot::SubplotGeometry, arg::Vector{Gadfly.Layer})
    for layer in arg
        push!(subplot.layers, layer)
    end
end

function add_subplot_element(p::SubplotGeometry, arg::Gadfly.GeometryElement)
    if !isempty(p.layers) && isa(p.layers[end].geom, Geom.Nil)
        p.layers[end].geom = arg
    else
        layer = Layer()
        layer.geom = arg
        push!(p.layers, layer)
    end
end

function add_subplot_element(subplot::SubplotGeometry, arg::Gadfly.StatisticElement)
    push!(subplot.statistics, arg)
end

function add_subplot_element(subplot::SubplotGeometry, arg::Gadfly.ScaleElement)
    push!(subplot.scales, arg)
    if (isa(arg, Scale.ContinuousColorScale) ||
        isa(arg, Scale.DiscreteColorScale)) && !haskey(subplot.guides, Guide.ColorKey)
        subplot.guides[Guide.ColorKey] = Guide.colorkey()
    end
end

function add_subplot_element(subplot::SubplotGeometry, arg::Gadfly.GuideElement)
    subplot.guides[typeof(arg)] = arg
end

function add_subplot_element{T <: Gadfly.Element}(subplot::SubplotGeometry,
                                                  arg::Type{T})
    add_subplot_element(subplot, arg())
end

function add_subplot_element(subplot::SubplotGeometry, coord::Gadfly.CoordinateElement)
    subplot.coord = coord
end

function add_subplot_element(subplot::SubplotGeometry, arg)
    error("Subplots do not support elements of type $(typeof(arg))")
end

type SubplotGrid <: SubplotGeometry
    layers::Vector{Gadfly.Layer}
    statistics::Vector{Gadfly.StatisticElement}
    scales::Vector{Gadfly.ScaleElement}
    guides::Dict{Type, Gadfly.GuideElement}
    coord::Gadfly.CoordinateElement
    free_x_axis::Bool
    free_y_axis::Bool

    # Current plot has no way of passing existing aesthetics. It always produces
    # these using scales.
    function SubplotGrid(elements::Gadfly.ElementOrFunctionOrLayers...;
                         free_x_axis=false, free_y_axis=false)
        subplot = new(Gadfly.Layer[], Gadfly.ScaleElement[], Gadfly.StatisticElement[],
                      Dict{Type, Gadfly.GuideElement}(), Coord.cartesian(),
                      free_x_axis, free_y_axis)

        for element in elements
            add_subplot_element(subplot, element)
        end

        # TODO: Handle default guides and statistics
        subplot
    end
end

function layers(geom::SubplotGrid)
    return geom.layers
end

const subplot_grid = SubplotGrid

function element_aesthetics(geom::SubplotGrid)
    vars = [:xgroup, :ygroup]
    for layer in geom.layers
        append!(vars, element_aesthetics(layer.geom))
    end
    vars
end

function default_scales(geom::SubplotGrid)
    scales = Gadfly.ScaleElement[]
    for stat in geom.statistics
        append!(scales, default_scales(stat))
    end

    for layer in geom.layers
        append!(scales, default_scales(layer.geom))
        append!(scales, default_scales(default_statistic(layer.geom)))
    end

    return scales
end

function element_coordinate_type(::SubplotGrid)
    return Gadfly.Coord.subplot_grid
end

# Render a subplot grid geometry, which consists of rendering and arranging
# many smaller plots.
function render(geom::SubplotGrid, theme::Gadfly.Theme,
                superplot_aes::Gadfly.Aesthetics,
                subplot_layer_aess::Vector{Gadfly.Aesthetics},
                subplot_layer_datas::Vector{Gadfly.Data},
                scales::Dict{Symbol, Gadfly.ScaleElement})

    # inherit aesthetics from the plot when needed but not provided
    for (layer, layer_aes) in zip(geom.layers, subplot_layer_aess)
        inherited_aes = element_aesthetics(layer.geom)
        push!(inherited_aes, :xgroup, :ygroup)
        for var in inherited_aes
            if getfield(layer_aes, var) === nothing
                setfield!(layer_aes, var, getfield(superplot_aes, var))
            end
        end
    end

    # work out the grid size
    m = 1
    n = 1
    for layer_aes in subplot_layer_aess
        if layer_aes.xgroup != nothing
            m = max(m, maximum(layer_aes.xgroup))
        end
        if layer_aes.ygroup != nothing
            n = max(n, maximum(layer_aes.ygroup))
        end
    end

    layer_aes_grid = Array(Matrix{Gadfly.Aesthetics}, length(geom.layers))
    for (i, (layer, aes)) in enumerate(zip(geom.layers, subplot_layer_aess))
        layer_aes_grid[i] = Gadfly.by_xy_group(aes, aes.xgroup, aes.ygroup, m, n)
    end

    layer_data_grid = Array(Matrix{Gadfly.Data}, length(geom.layers))
    for (i, (layer, data, aes)) in enumerate(zip(geom.layers, subplot_layer_datas,
                                                 subplot_layer_aess))
        layer_data_grid[i] = Gadfly.by_xy_group(data, aes.xgroup, aes.ygroup, m, n)
    end

    coord = geom.coord
    plot_stats = Gadfly.StatisticElement[stat for stat in geom.statistics]
    layer_stats = [isempty(layer.statistics) ?
                        Gadfly.StatisticElement[Geom.default_statistic(layer.geom)] : layer.statistics
                   for layer in geom.layers]

    for i in 1:n, j in 1:m
        Scale.apply_scales(geom.scales,
                           Gadfly.Aesthetics[layer_aes_grid[k][i, j]
                                             for k in 1:length(geom.layers)],
                           Gadfly.Data[layer_data_grid[k][i, j]
                                       for k in 1:length(geom.layers)]...)

        for (k, stats) in enumerate(layer_stats)
            Stat.apply_statistics(stats, scales, coord, layer_aes_grid[k][i, j])
        end
    end

    # apply geom-wide statistics
    geom_aes = Gadfly.concat([layer_aes_grid[k][i,j]
                              for i in 1:n, j in 1:m, k in 1:length(geom.layers)]...)
    geom_stats = Gadfly.StatisticElement[]

    has_stat_xticks = false
    has_stat_yticks = false
    for guide in values(geom.guides)
        stat = default_statistic(guide)
        if !isa(stat, Gadfly.Stat.identity)
            if isa(stat, Gadfly.Stat.TickStatistic)
                if stat.out_var == "x"
                    has_stat_xticks = true
                elseif stat.out_var == "y"
                    has_stat_yticks = true
                end
            end
            push!(geom_stats, stat)
        end
    end

    if !geom.free_x_axis && !has_stat_xticks
        push!(geom_stats, Stat.xticks())
    end

    if !geom.free_y_axis && !has_stat_yticks
        push!(geom_stats, Stat.yticks())
    end

    Stat.apply_statistics(geom_stats, scales, coord, geom_aes)
    aes_grid = [geom_aes for i in 1:n, j in 1:m]

    # if either axis is on a free scale, we need to apply row/column-wise
    # tick statistics.
    if geom.free_x_axis
        for j in 1:m
            col_aes = Gadfly.concat([layer_aes_grid[k][i, j]
                                     for i in 1:n, k in 1:length(geom.layers)]...)
            Gadfly.inherit!(col_aes, geom_aes)
            Stat.apply_statistic(Stat.xticks(), scales, coord, col_aes)

            aes_grid[:, j] = col_aes
        end
    end

    if geom.free_y_axis
        for i in 1:n
            row_aes = Gadfly.concat([layer_aes_grid[k][i, j]
                                     for j in 1:m, k in 1:length(geom.layers)]...)
            Gadfly.inherit!(row_aes, geom_aes)
            Stat.apply_statistic(Stat.yticks(), scales, coord, row_aes)

            for j in 1:m
                aes_grid[i, j] = Gadfly.inherit(row_aes, aes_grid[i, j])
            end
        end
    end

    for k in length(geom.layers), i in 1:n, j in 1:m
        Gadfly.inherit!(layer_aes_grid[k][i, j], aes_grid[i, j])
    end

    # TODO: this assumes a rather ridged layout

    hascolorkey = haskey(geom.guides, Guide.ColorKey)

    if hascolorkey
        xprop = [(isodd(j) ? 1.0 : NaN) for j in 3:2*m+2]
        tbl = table(n + 2, 2*m + 2, 1:n, 3:2*m+2,
                    x_prop=xprop, y_prop=ones(n),
                    fixed_configs=Any[
                        [(i, 1) for i in 1:n],
                        [(i, 2) for i in 1:n],
                        [(n+1, j) for j in 3:2:m+2],
                        [(n+2, j) for j in 3:2:m+2]])
    else
        tbl = table(n + 2, m + 2, 1:n, 3:m+2,
                    x_prop=ones(m), y_prop=ones(n),
                    fixed_configs=Any[
                        [(i, 1) for i in 1:n],
                        [(i, 2) for i in 1:n],
                        [(n+1, j) for j in 3:m+2],
                        [(n+2, j) for j in 3:m+2]])
    end

    xtitle = "x"
    for v in [:x, :xmin, :xmax]
        if haskey(superplot_aes.titles, v)
            xtitle = superplot_aes.titles[v]
            break
        end
    end

    ytitle = "y"
    for v in [:y, :ymin, :ymax]
        if haskey(superplot_aes.titles, v)
            ytitle = superplot_aes.titles[v]
            break
        end
    end

    xlabels = superplot_aes.xgroup_label(1:m)
    ylabels = superplot_aes.ygroup_label(1:n)
    subplot_padding = 2mm

    # This assumes non of the layers themselves are subplot geometries
    layer_subplot_aess = Vector{Gadfly.Aesthetics}[Array(Gadfly.Aesthetics, 0)
                                                   for _ in 1:length(geom.layers)]
    layer_subplot_datas = Vector{Gadfly.Data}[Array(Gadfly.Data, 0)
                                                   for _ in 1:length(geom.layers)]

    for i in 1:n, j in 1:m
        p = Plot()
        p.theme = theme
        p.layers = geom.layers
        guides = Gadfly.GuideElement[]

        for guide in values(geom.guides)
            if typeof(guide) in [Guide.XTicks, Guide.YTicks, Guide.XLabel, Guide.YLabel]
                continue
            end
            push!(guides, guide)
        end

        p.scales = collect(ScaleElement, values(scales))

        # default guides
        push!(guides, Guide.background())

        if i == n
            push!(guides, get(geom.guides, Guide.XTicks, Guide.xticks()))

            if !is(superplot_aes.xgroup, nothing)
                push!(guides, get(geom.guides, Guide.XLabel, Guide.xlabel(xlabels[j])))
            end
        else
            push!(guides, Guide.xticks(label=false))
        end

        joff = 0
        if j == 1
            joff += 1
            push!(guides, get(geom.guides, Guide.YTicks, Guide.yticks()))
            if !is(superplot_aes.ygroup, nothing)
                joff += 1
                push!(guides, get(geom.guides, Guide.YLabel, Guide.ylabel(ylabels[i])))
            end
        else
            push!(guides, Guide.yticks(label=false))
        end

        subtbl = Gadfly.render_prepared(
                            p, Gadfly.Coord.cartesian(),
                            aes_grid[i, j],
                            Gadfly.Aesthetics[layer_aes_grid[k][i, j]
                                              for k in 1:length(geom.layers)],
                            layer_stats,
                            layer_subplot_aess,
                            layer_subplot_datas,
                            scales,
                            guides,
                            table_only=true)

        # copy over the correct units, since we are reparenting the children
        for u in 1:size(subtbl, 1), v in 1:size(subtbl, 2)
            for child in subtbl[u, v]
                if isnull(child.units)
                    child.units = subtbl.units
                end
            end
        end

        # All of the below needs to be rewritten to take into account the
        # possibility that subplots will have their own external guides.
        # We could special case it.
        if hascolorkey
            tbl[i, 2 + 2*j - 1] = pad(subtbl[1, 1 + joff],
                                      j > 1 ? subplot_padding : 0mm,
                                      subplot_padding,
                                      subplot_padding,
                                      i < n ? subplot_padding : 0mm)
            tbl[i, 2 + 2*j] = subtbl[1, 2 + joff]
        else
            tbl[i, 2 + j] = pad(subtbl[1, 1 + joff],
                                j > 1 ? subplot_padding : 0mm,
                                subplot_padding,
                                subplot_padding,
                                i < n ? subplot_padding : 0mm)
        end

        # bottom guides
        if i == n
            for k in 2:size(subtbl, 1)
                tbl[i + k - 1, hascolorkey ? 2 + 2*j - 1: 2 + j] =
                    pad(subtbl[k, 1 + joff],
                        j > 1 ? subplot_padding : 0mm,
                        subplot_padding,
                        0mm, 0mm)
            end
        end

        # left guides
        if j == 1
            for k in 1:(size(subtbl, 2)-1)
                tbl[i, k] =
                    pad(subtbl[1, k],
                        0mm, 0mm,
                        subplot_padding,
                        i < n ? subplot_padding : 0mm)
            end
        end
    end

    return compose!(context(), tbl)
end


text:
using PPInterpolation, Test
using ForwardDiff

@testset "Akima" begin
    x = [3.0, 5.0, 6.0, 8.0, 9.0, 11.0, 12.0, 14.0, 15.0]
    f = [10.0, 10.0, 10.0, 10.0, 10.5, 15.0, 50.0, 60.0, 85.0]
    varieties = [C2Hyman89(), C2MP(), HuynRational(), VanLeer(), FritschButland(), Brodlie(), VanAlbada()]
    for v in varieties
        spline = makeCubicPP(x, f, PPInterpolation.SECOND_DERIVATIVE, 0.0, PPInterpolation.SECOND_DERIVATIVE, 0.0, v)
        nSample = 1000
        previousVal = spline(x[1])
        for i = 1:nSample
            ti = x[1] + (x[end] - x[1]) * (i) / (nSample)
            zval = spline(ti)
            if i > 2
                @test zval >= previousVal
            end
            previousVal = zval
        end
    end
    #     varieties = [ C2Hyman89(), C2MP(), HuynRational(), VanLeer(), FritschButland(), Brodlie(), VanAlbada() ]^C
    #  p = plot()
    #  for v in varieties
    #        spline = makeCubicPP(x,f,PPInterpolation.SECOND_DERIVATIVE,0.0,PPInterpolation.SECOND_DERIVATIVE,0.0,v)
    #        plot!(p, t, spline.(t), label=string(v))
    #        end
    #  plot!(legend=:top)

end

@testset "FrtschCarlson" begin
    x = [7.99, 8.09, 8.19, 8.7, 9.2, 10.0, 12.0, 15.0, 20.0]
    f = [0.0, 2.76429e-5, 4.37498e-2, 0.169183, 0.469428, 0.943740, 0.998636, 0.999919, 0.999994]

    varieties = [C2Hyman89(), C2MP(), HuynRational(), VanLeer(), FritschButland(), Brodlie(), VanAlbada()]
    for v in varieties
        spline = makeCubicPP(x, f, PPInterpolation.SECOND_DERIVATIVE, 0.0, PPInterpolation.SECOND_DERIVATIVE, 0.0, v)
        nSample = 1000
        previousVal = spline(x[1])
        for i = 1:nSample
            ti = x[1] + (x[end] - x[1]) * (i) / (nSample)
            zval = spline(ti)
            if i > 2
                @test zval >= previousVal
            end
            previousVal = zval
        end
    end
end
@testset "HaganCurveAuto" begin
    #for quintic: first pass compute cubic first derivative, second pass, use cubic b, and optimize second derivative.
    t = [0.0, 0.1, 1.0, 2.0, 3.0, 4.0, 9.0, 20.0, 30.0]
    y = [0.06, 0.06, 0.06, 0.05, 0.05, 0.05, 0.065, 0.06, 0.06]
    t = [0.0, 0.1, 1.0, 4.0, 9.0, 20.0, 30.0]
    y = [0.081, 0.081, 0.07, 0.044, 0.07, 0.04, 0.03]
    t = [0, 0.1, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 2.0, 2.5, 3.0, 3.5]
    y = [0.012, 0.008, 0.0104, 0.00995, 0.010633333333333333, 0.010849999999999999, 0.010179999999999998, 0.009816666666666665, 0.008612499999999999, 0.007789999999999998, 0.007741666666666665, 0.007635714285714285]

    # t = [0.0, 0.1, 1, 2, 3, 4, 5, 6, 7, 30.0]
    # f = [0.01, 0.01, 0.01, 0.01, 0.01, 0.03, 0.03, 0.03, 0.03, 0.02]
    # y = zeros(length(t))
    # y[1] = f[1]
    # y[2:end] = (y[1:end-1] .* t[1:end-1] + (t[2:end]-t[1:end-1]) .* f[2:end]) ./ t[2:end]
    # y[3:end-2] = f[3:end-2] ./ t[3:end-2]
    # println(y)
    yt = zeros(length(t))
    yt[1:end] = y[1:end] .* t[1:end]
    spline = makeCubicPP(t, yt, PPInterpolation.SECOND_DERIVATIVE, 0.0, PPInterpolation.SECOND_DERIVATIVE, 0.0, C2())
    nSample = 1000
    for i = 1:nSample
        ti = (t[end]) * (i) / (nSample)
        val = evaluateDerivative(spline, ti)
        zval = spline(ti) / ti
        println(ti, " CubicPP ", val, " ", zval)
        autoval = ForwardDiff.derivative(spline, ti)
        @test isapprox(val, autoval, atol = 1e-15)
    end
    #AUTODIFF // where is autodiff useful (hagan wave read),  BigFloat
end


text:
using Test
using Distances: Euclidean, CosineDist
using Random
using SparseArrays
using LinearAlgebra
using SimilaritySearch
using UMAP
using UMAP: initialize_embedding, fuzzy_simplicial_set, compute_membership_strengths, smooth_knn_dist_opt_binsearch, spectral_layout, optimize_embedding, combine_fuzzy_sets, fit_ab, SMOOTH_K_TOLERANCE

include("utils_tests.jl")
include("umap_tests.jl")


text:
"""
    function IF(C ;order=4, fs=1)
Calculate instantaneous frequencies
"""
function IF(C; order=4, fs=1)

	n,N = size(C)
    t   = 1:n
	Phi = zeros(n-1,N)

	@inbounds for i =1:N
        H1 = Spline1D(t, C[:,i], k = order)
        Phi[:,i] = (fs/2π) .* diff(unwrap(angle.(hilbert(H1(t)))))
    end

	return Phi
end

function unwrap!(v::AbstractVector, lim=π)
  for i in 2:length(v)
    while unwrapped[i] - unwrapped[i-1] >= lim
      unwrapped[i] -= 2*lim
    end
    while unwrapped[i] - unwrapped[i-1] <= -lim
      unwrapped[i] += 2*lim
    end
  end
  return unwrapped
end

unwrap(v) = unwrap!(copy(v))

max_stars_repo_name: alexmorley/EmpiricalModeDecomposition

max_issues_repo_name: alexmorley/EmpiricalModeDecomposition

max_forks_repo_name: alexmorley/EmpiricalModeDecomposition