text:
function ==(sent::ETree, received::ETree)

    if sent.name == received.name &&
       sent.attr == received.attr
        if isempty(sent.elements)
            return true
        end
        for elem_sent in sent.elements
            if isa(elem_sent, ETree)
                # there must be one of these in received
                elem_rec_index = 0

                for (i,node) in enumerate(received.elements)
                    if isa(node, ETree) && node.name == elem_sent.name
                        elem_rec_index = i
                        break
                    end
                end

                if elem_rec_index != 0 &&
                   elem_sent == received.elements[elem_rec_index]

                   return true
               end
            end
        end
    end

    false
end

message_payload_to_etree(message::SCPMessage) = message_payload_to_etree(string_from_buffer(message.payload))
function message_payload_to_etree(payload::String)
    str = "<MSG>" * payload * "</MSG>"
    root = xp_parse(str)
    index_of_first_etree = findfirst(e->isa(e, ETree), root.elements)
    @assert(index_of_first_etree != 0)
    root.elements[index_of_first_etree]::ETree
end

function message_contains_packet_with_element(message::ETree, name::String, elementname::String)
    if message.name == name
        payload = message.elements[findfirst(elem->isa(elem, ETree), message.elements)]
        if payload.name == elementname
            return true
        end
    end
    false
end
function message_contains_packet_with_element(message::SCPMessage, name::String, elementname::String)
    str = "<MSG>" * string_from_buffer(message.payload) * "</MSG>"
    for node in xp_parse(str).elements
        if isa(node, ETree)
            if message_contains_packet_with_element(node, name, elementname)
                return true
            end
        end
    end
    false
end

max_stars_repo_name: sisl/VirtualTestDrive

max_issues_repo_name: sisl/VirtualTestDrive

max_forks_repo_name: sisl/VirtualTestDrive


text:
abstract type AbstractSymmetry end
struct EvenSymmetry <: AbstractSymmetry end
struct OddSymmetry <: AbstractSymmetry end
struct NoSymmetry <: AbstractSymmetry end
const Even = EvenSymmetry()
const Odd = OddSymmetry()

ext: jl


text:
data = split(strip(read(stdin, String)), '\n')

x = y = 0
dx = 1
dy = 0
for instruction in data
    global x, y
    global dx, dy

    op, val = instruction[begin], parse(Int, instruction[2:end])

    if op == 'F'
        x += dx * val
        y += dy * val
    elseif op == 'N'
        y += val
    elseif op == 'S'
        y -= val
    elseif op == 'E'
        x += val
    elseif op == 'W'
        x -= val
    elseif op ∈ ('L', 'R')
        if op == 'R'
            val = 360 - val
        end

        if val == 90
            dx, dy = -dy, dx
        elseif val == 270
            dx, dy = dy, -dx
        else
            @assert val == 180 "val = $val"

            dx = -dx
            dy = -dy
        end
    end
end

println("Part 1 : $(abs(x) + abs(y))")

x = y = 0

# Waypoint
dx = 10
dy = 1
for instruction in data
    global x, y
    global dx, dy

    op, val = instruction[begin], parse(Int, instruction[2:end])

    if op == 'F'
        x += dx * val
        y += dy * val
    elseif op == 'N'
        dy += val
    elseif op == 'S'
        dy -= val
    elseif op == 'E'
        dx += val
    elseif op == 'W'
        dx -= val
    elseif op ∈ ('L', 'R')
        if op == 'R'
            val = 360 - val
        end

        if val == 90
            dx, dy = -dy, dx
        elseif val == 270
            dx, dy = dy, -dx
        else
            @assert val == 180 "val = $val"

            dx = -dx
            dy = -dy
        end
    end
end

println("Part 2 : $(abs(x) + abs(y))")


text:
using Chain
using ShiftedArrays

# PART 1
@chain "01/data" begin
    readlines # read into an array, each line is an element
    map(x -> parse(Int32, x), _) # convert strings to integers
    _ .- lag(_) # compute difference with previous value
    _[2:end] # remove first element, which will be missing
    filter(x -> x > 0, _) # filter to the increases in value
    length # now count
end

# PART 2
@chain "01/data" begin
    readlines # read into an array, each line is an element
    map(x -> parse(Int32, x), _) # convert strings to integers
    _ + lag(_) + lag(_, 2) # compute the sum of each value with its previous 2
    _ - lag(_) # now compute the difference with the previous element
    _[4:end] # remove first three elements, since we're summing triples
    filter(x -> x > 0, _) # filter to positive values, which are the increases in value
    length # now count
end


text:
struct DiffRuleVariable
    diffrule::DiffRule
    vars::OffsetVector{Variable,Vector{Variable}}
    funcold::Func
end
function DiffRuleVariable(∂m::∂Model, dr::DiffRule, func::Func)
    nparents = num_parents(dr)
    vars = OffsetVector{Variable}(undef, -nparents:last(last(dr.sections)))
    drv = DiffRuleVariable( dr, vars, func )
    vparents = parents(func)
    moldvars = ∂m.mold.vars
    for i ∈ -nparents:-1
        # drv.vars[i] = getvar!(∂m, moldvars[vparents[i + nparents + 1]])
        drv.vars[i] = getvar!(∂m, moldvars[vparents[i + nparents + 1]])
    end
    drv.vars[0] = var0 = getvar!(∂m.m, diffsym(moldvars[func.output[]]))
    drv
end

function add_section!(∂m::∂Model, drv::DiffRuleVariable, sectionid::Int)
    @unpack vars, funcold, diffrule = drv
    @unpack sections, returns, instructions, dependencies = diffrule
    section = sections[sectionid]
    if isone(sectionid)
        ret = returns[1]
    else
        str = ReverseDiffExpressionsBase.section_two_returns(diffrule)
        if sectionid == 2
            str || return
            ret = returns[2]
        else
            if (vars[sectionid - num_parents(diffrule) - 3]).tracked
                ret = returns[sectionid + str - 1]
            else
                return
            end
        end
    end
    for i ∈ section
        func = Func(instructions[i], false)
        # (instructions[i].instr === :constrain_pullback!) && uses!(func, ∂m.gradvar)
        for j ∈ dependencies[i]
            uses!(func, vars[j])
        end
        vars[i] = if i == ret && sectionid != 2
            if isone(sectionid)
                retv = getvar!(∂m, ∂m.mold.vars[funcold.output[]])
            else
                parent = vars[sectionid - num_parents(diffrule) - 3]
                retv = getvar!(∂m.m, diffsym(name(parent)))
            end
        else
            retv = getvar!(∂m.m, gensym(:temp))
        end
        retv.tracked = true
        returns!(func, retv)
        retold = ∂m.mold.vars[funcold.output[]]
        addfunc!(∂m.m, func)
    end
end

function add_fallback_diffrule!(∂m::∂Model, func::Func)
    @unpack m, mold = ∂m
    # add func as var, return DiffRuleVariable
    v = addvarref!(∂m.m, asexpr(func.instr))
    drv = DiffRuleVariable(∂m, ReverseDiffExpressionsBase.FALLBACK_RULES[length(parents(func))], func)

    retvsym = gensym(:rrule_LHS)
    rrulefunc = Func(Instruction(:ChainRules,:rrule), false)
    drv.vars[1] = rruleretv = addvar!(∂m.m, retvsym)
    rruleretv.tracked = true
    returns!(rrulefunc, rruleretv)
    uses!(rrulefunc, v)
    foreach(arg -> uses!(rrulefunc, arg), @view(drv.vars[begin:-1]))
    addfunc!(m, rrulefunc)
    funcfirst = Func(Instruction(:Base,:first), false)
    uses!(funcfirst, rruleretv)
    drv.vars[2] = retvar = getvar!(∂m, mold.vars[func.output[]])
    returns!(funcfirst, retvar)
    retvar.tracked = true
    addfunc!(m, funcfirst)
    drv
end
function add_diffrule!(∂m::∂Model, dr::DiffRule, func::Func)
    drv = DiffRuleVariable(∂m, dr, func)
    add_section!(∂m, drv, 1)
    drv
end

function add_view_diff!(∂m::∂Model, func::Func)
    gi_forward = Func(func.instr, false)
    gi_reverse = Func(func.instr, false)
    vparents = parents(func)
    moldvars = ∂m.mold.vars
    p = moldvars[first(vparents)]
    inputvf = getvar!(∂m, p); inputvf.tracked = true
    inputvr = getvar!(∂m.m, diffsym(name(p))); inputvr.tracked = true
    uses!(gi_forward, inputvf)
    uses!(gi_reverse, inputvr)
    for i ∈ @view(vparents[2:end])
        indold = moldvars[i]
        if isref(indold)
            ind = addvar!(∂m.m, Symbol(""))
            ind.ref = indold.ref
        else
            ind = addvar!(∂m.m, name(indold))
        end
        uses!(gi_forward, ind)
        uses!(gi_reverse, ind)
    end
    retv = moldvars[func.output[]]
    retvf = getvar!(∂m, retv); retvf.tracked = true
    returns!(gi_forward, retvf)
    retvr = getvar!(∂m.m, diffsym(name(retv))); retvr.tracked = true
    returns!(gi_reverse, retvr)
    push!(inputvr.paireddeps, identifier(retvr))
    addfunc!(∂m.m, gi_forward)
    addfunc!(∂m.m, gi_reverse)
    nothing
end

function add_probdist_diff!(∂m::∂Model, func::Func)
    ∂loglik = Func(func.instr, true)
    moldvars = ∂m.mold.vars
    for i ∈ parents(func)
        uses!(∂loglik, getvar!(∂m, moldvars[i]))
    end
    tup = addvar!(∂m.m, gensym(:tup))
    tup.tracked = true
    returns!(∂loglik, tup)
    addfunc!(∂m.m, ∂loglik)
    getconstindex!(∂m.m, ∂m.m.vars[0], tup, 1)
    ∂tup = addvar!(∂m.m, gensym(:∂tup))
    getconstindex!(∂m.m, ∂tup, tup, 2)
    for (ind,i) ∈ enumerate(parents(func))
        p = moldvars[i]
        if istracked(p)
            retv = getvar!(∂m.m, diffsym(name(p))); retv.tracked = true
            getconstindex!(∂m.m, retv, ∂tup, ind)
        end
    end
end
function add_constrain_diff!(∂m::∂Model, funcold::Func)
    drvc = DiffRuleVariable(∂m, DERIVATIVERULES[InstructionArgs(Instruction(:DistributionParameters,:constrain),3)], funcold)
    constrain_pullback = Func(Instruction(:constrain_pullback!), false)
    uses!(constrain_pullback, ∂m.gradvar)
    for i ∈ -3:-1
        uses!(constrain_pullback, drvc.vars[i])
    end
    tempv = getvar!(∂m.m, gensym(:constrainpullbacktup)); tempv.tracked = true
    returns!(constrain_pullback, tempv)
    addfunc!(∂m.m, constrain_pullback)
    retv = getvar!(∂m, ∂m.mold.vars[funcold.output[]]); retv.tracked = true
    getconstindex!(∂m.m, retv, tempv, 1)
    ∂retv = getvar!(∂m.m, diffsym(name(retv))); ∂retv.tracked = true
    getconstindex!(∂m.m, ∂retv, tempv, 2)

    constrain_reverse = Func(Instruction(:constrain_reverse!), false)
    uses!(constrain_reverse, ∂retv)
    uses!(constrain_reverse, drvc.vars[-1])
    nothingvar = addvar!(∂m.m, gensym(:nothing))
    returns!(constrain_reverse, nothingvar)
    push!(∂m.m.targets, identifier(nothingvar))
    addfunc!(∂m.m, constrain_reverse)
    nothing
end
function add_no_diff!(∂m::∂Model, funcold::Func)
    func = Func(funcold.instr, funcold.probdistapi)
    for i ∈ parents(funcold)
        uses!(func, getvar!(∂m, ∂m.mold.vars[i]))
    end
    returns!(func, getvar!(∂m, ∂m.mold.vars[funcold.output[]]))
    addfunc!(∂m.m, func)
    nothing
end
function has_tracked_parents(∂m::∂Model, func::Func)
    any(i -> istracked(∂m.mold.vars[i]), parents(func))
end

function add_func_diff!(∂m::∂Model, func::Func)
    has_tracked_parents(∂m, func) || return add_no_diff!(∂m, func)
    if func.probdistapi
        return add_probdist_diff!(∂m, func)
    elseif func.instr.instr === :constrain
        return add_constrain_diff!(∂m, func)
    elseif isindexfunc(func)
        if iszero(∂m.mold.vars[func.output[]].varid)
            return add_no_diff!(∂m, func)
        else
            return add_view_diff!(∂m, func)
        end
    end
    instrargs = InstructionArgs(func)
    dr = get(DERIVATIVERULES, instrargs, nothing)
    drv = isnothing(dr) ? add_fallback_diffrule!(∂m, func) : add_diffrule!(∂m, dr, func)
    for i ∈ 2:length(drv.diffrule.sections)
        add_section!(∂m, drv, i)
    end
end


text:
module ToyModelsMD

using Random

# model potentials
export MullerBrown,potential,force#,force!
export GaussianModel,GaussianModelTD,GaussianModel2D
export TrainableGaussianModel,TrainableGaussianModelTD
export NNModel
export MixedModel
# thermostats
export AbstractThermostat
export Langevin,LangevinND
export ActiveLangevin
# systems
export AbstractSystem,AbstractThermostattedSystem
export System,ThermostattedSystem,ActiveBrownianSystem
#export fderivative,gderivative
# Callbacks
export AbstractCallback,Callback
export MCVBCallback,initialize!,initializeall!
export MCVBTCallback,initialize!,initializeall!
export MixedCallback
# integrators
export StochasticEuler,runtraj!
# valuebaseline functions
export AbstractValueBaseline
export GaussianValueBaseline,GaussianValueBaselineTD,GaussianValueBaseline2D
export NNValueBaseline
## Monte Carlo samplers
#export MCSampler,BiasedMCSampler
#export sweep!,step!
# training algorithms
export mcvbtrainsgd!,mcvbtrainsgdpar!
export mcvbtrainadam!,mcvbtrainadampar!
export mcvbttrainsgdpar!

include("abstracttypes.jl")
include("potentials.jl")
include("thermostats.jl")
include("systems.jl")
include("activesystem.jl")
include("mullerbrown.jl")
include("gaussianmodel.jl")
include("trainablegaussianmodel.jl")
include("neuralnetwork.jl")
include("mixedmodel.jl")
include("callbacks.jl")
include("mixedcallback.jl")
include("integrators.jl")
include("mcvb.jl")
include("mcvbt.jl")
include("gaussianvaluebaseline.jl")
include("nnvaluebaseline.jl")
#include("montecarlo.jl")
include("trainforces.jl")
include("trainmcvbt.jl")

end


text:
using TimeIt
using AperturePhotometry

data = ones(10, 10)

error = 0.1 * ones(10, 10)
mask = falses(10, 10)

@timeit sum_circle(data, 5.0, 5.0, 4.0)
@timeit sum_circle(data, 5.0, 5.0, 4.0; mask=error, error=error, maskthresh=0.2)

@timeit sum_circann(data, 5.0, 5.0, 2.0, 4.0)
@timeit sum_circann(data, 5.0, 5.0, 2.0, 4.0; mask=error, error=error, maskthresh=0.2)

@timeit sum_ellipse(data, 5.0, 5.0, 3.0, 2.0, pi/4.)
@timeit sum_ellipse(data, 5.0, 5.0, 3.0, 2.0, pi/4.; mask=error, error=error, maskthresh=0.2)


text:
## Effects
# Effects are super useful to understand the actual modelfits. If you are an EEG-Researcher, you can think of effects as the "modelled ERPs", and the coefficients as the "difference-waves".
# In some way, we are fitting a model with coefficients and then try to get back the "original" ERPs - of course typically with some effect adjusted, overlap removed or similar - else why bother ;)

# Define some packages

using Unfold
using DataFrames
using Random
using CSV
using UnfoldMakie

# # Setup things
# Let's load some data and fit a model of a 2-level categorical and a continuous predictor with interaction.
include(joinpath(dirname(pathof(Unfold)), "../test/test_utilities.jl") ) # to load data

data, evts = loadtestdata("test_case_3a")

basisfunction = firbasis(τ = (-0.5, 1.5), sfreq = 20, name = "basisA")

evts.conditionA= ["off","on"][(evts.conditionA .+1)] # convert evts into categorical
f = @formula 0 ~ 1+conditionA*continuousA # 1

m = fit(UnfoldModel, Dict(Any=>(f,basisfunction)), evts, data,eventcolumn="type")

# Plot the results
plot_results(coeftable(m))

# As expected, we get four lines - the interaction is flat, the slope of the continuous is around 4, the categorical effect is at 3 and the intercept at 0 (everything is dummy coded by default)
# ### Effects
# A convenience function is [effects](@ref). It allows to specify effects on specific levels, while setting non-specified ones to a typical value (usually the mean)

eff = effects(Dict(:conditionA => ["off","on"]),m)
plot_results(eff;color=:conditionA)

# We can also generate continuous predictions
eff = effects(Dict(:continuousA => -1:0.1:1),m)
plot_results(eff;color=:continuousA)

# or split it up by condition
eff = effects(Dict(:conditionA=>["off","on"],:continuousA => -1:0.1:1),m)
plot_results(eff;color=:continuousA,col=:conditionA)

# ## What is typical anyway?
# The user can specify the typical function applied to the covariates/factors that are marginalized over. This offers even greater flexibility.
# Note that this is rarely necessary, in most applications the mean will be assumed.
eff_max = effects(Dict(:conditionA=>["off","on"]),m;typical=maximum) # typical continuous value fixed to 1
eff_max.typical .= :maximum
eff = effects(Dict(:conditionA=>["off","on"]),m)
eff.typical .= :mean # mean is the default

plot_results(vcat(eff,eff_max);color=:conditionA,col=:typical)

max_stars_repo_name: unfoldtoolbox/jlunfold

max_issues_repo_name: unfoldtoolbox/jlunfold

max_forks_repo_name: unfoldtoolbox/jlunfold


text:
using Stochastics
using Test

@testset "DescriptiveStatistics.jl" begin
    x = [1, 2]
    y = [3, 4]
    @test Stochastics.sample_mean(x) == 1.5
    @test Stochastics.sample_variance(x) == 0.5
    @test Stochastics.sum_of_squares(x) == 0.5
    @test Stochastics.sum_of_squares(x, y) == 0.5
    @test Stochastics.two_sample_covariance(x, y) == 0.5
    @test abs2(Stochastics.two_sample_correlation_coefficient(x, y) - 1) < 10e-8
    @test Stochastics.two_sample_variance(x, y) == 0.5
end

max_issues_repo_name: persson-io/Stochastics

max_forks_repo_name: persson-io/Stochastics


text:
import Pkg
Pkg.activate(normpath(@__DIR__))

using ProjectMoana
import Images: save
import Base: write
import ColorTypes: RGB

using ArgParse

function parse_commandline()
    s = ArgParseSettings(description = "This program generates an image of 10 spheres. Try me!",
                               usage = "usage: [--help] [--width W] [--height H] [--camera C] [--angle α] [--distance D] [--file_out FILENAME] [--render_alg ALG] [--a A] [--seq S]" ,
                              epilog = "Let's try again!")

    @add_arg_table s begin
        "--width"
            help = "width of the image"
            required = true
            arg_type = Int
        "--height"
            help = "height of the image"
            required = true
            arg_type = Int
        "--camera"
            help = "type of camera [O, for Orthogonal, or P, for Perspective]"
            required = false
            default = "P"
            arg_type = String
        "--angle"
            help = "angle of the z-axis rotation applied to camera (in degrees)"
            required = false
            default = 0.
            arg_type = Float64
        "--dist"
            help = "distance of a Perspective Camera"
            required = false
            default = 1.
            arg_type = Float64
        "--file_out"
            help = "name of the output file (without extension)"
            required = false
            default = "demo"
            arg_type = String
        "--render_alg"
            help = "type of rendering algortihm [O for On-Off, F for Flat, P for Path Tracer]"
            required = false
            default = "P"
            arg_type = String
        "--a"
            help = "a_factor for normalizing image luminosity during the convertion"
            required = false
            default = 1.
            arg_type = Float64
        "--seq"
            help = "sequence number for PCG generator"
            required = false
            default = 54
            arg_type = Int
        "--nrays"
            help = "Number of rays for antialasing"
            required = false
            default = 4
            arg_type = Int
   end

    return parse_args(s)
end

function main()
    params = parse_commandline()

    w = params["width"]
    h = params["height"]
    a = w/h
    d = params["dist"]
    camera_tr = rotation_z(params["angle"]*π/180.0) * translation(Vec(-2.0,0.,0.))
    image = HdrImage(w, h)
    file_out_pfm = "$(params["file_out"]).pfm"
    file_out_png = "$(params["file_out"]).png"
    algorithm = params["render_alg"]
    seq = convert(UInt64, params["seq"])
    samples_per_pixel = params["nrays"]

    samples_per_side = sqrt(samples_per_pixel)
    if samples_per_side^2 != samples_per_pixel
        println("Error, the number of rays per pixel ($samples_per_pixel) must be a perfect square")
        return
    end

# Creating WORLD with sky, sun, checkered gorud, Mars, Jupiter and mirror
    world = World()
    WHITE = RGB(1.,1.,1.)
    BLACK = RGB(0.,0.,0.)

    sky_color = Material(DiffuseBRDF(UniformPigment(RGB(.1,.5,.99))), UniformPigment(BLACK))
    sky = Plane(translation(Vec(0.,0.,500.)) * rotation_y( pi/3.), sky_color)
    add_shape(world, sky)

    sun_color = Material(DiffuseBRDF(UniformPigment(RGB(0.6,0.8,1.))), UniformPigment(WHITE))
    sun = Plane(translation(Vec(-100,0.,0.)) * rotation_y(pi/2.), sun_color )
    add_shape(world, sun)

    ground_color = Material(DiffuseBRDF(CheckeredPigment(RGB(0.3, 0.5, 0.1), RGB(1, 0.6, 0.8), 4)))
    ground = Plane(translation(Vec(0.,0.,-1.5)) * rotation_y(pi/12.) * rotation_x(pi/24.) * scaling(Vec(5,5,5)), ground_color )
    add_shape(world, ground)

    mirror_color = Material(SpecularBRDF())
    mirror = Sphere(translation(Vec(1.5,2,-0.5))*scaling(Vec(1.,1.5,1.)), mirror_color)
    add_shape(world, mirror)

    planet2_color = Material(DiffuseBRDF(ImagePigment(read_pfm_image("./examples/jupiter_texture.pfm"))))
    planet2 = Sphere(translation(Vec(6.,-4.,3.)) * scaling(Vec(3.,3.,3.)), planet2_color)
    add_shape(world, planet2)

    planet_color = Material(DiffuseBRDF(ImagePigment(read_pfm_image("./examples/mars_texture.pfm"))))
    planet = Sphere(translation(Vec(1.,-1.,-1.1)) * scaling(Vec(1.,1.,1.)), planet_color)
    add_shape(world, planet)

    println("World objects created.")

# Creating a Perspective of Orthogonal CAMERA
    if params["camera"] == "O"
        camera = OrthogonalCamera(a, camera_tr)
    elseif params["camera"] == "P"
        camera = PerspectiveCamera(a, camera_tr, d)
    end

# Creating an ImageTracer object
    tracer = ImageTracer(image, camera, samples_per_pixel)

    println("Observer initialized.")

    print("Computing ray intersection ")
    if params["render_alg"] == "F"
        println("using Flat renderer")
        renderer = Flat_Renderer(world, RGB(0.4,0.4,0.4))
        fire_all_rays(tracer, Flat, renderer)
    elseif params["render_alg"] == "P"
        println("using Path Tracer renderer")
        renderer = PathTracer_Renderer(world; background_color=RGB(0.,0.,0.), pcg=PCG(UInt64(42), seq), num_of_rays=2, max_depth=2, russian_roulette_limit=3)
        fire_all_rays(tracer, PathTracer, renderer)
    else
        println("using On/Off renderer")
        renderer = OnOff_Renderer(world)
        fire_all_rays(tracer, OnOff, renderer)
    end

    println("\nRay intersections evaluated.")

# Saving the PFM FILE
    write(file_out_pfm, tracer.image)
    println("$(file_out_pfm) has been written to disk.")

# Automatic CONVERSION TO JPEG FILE
    normalize_image(tracer.image, params["a"])
    clamp_image(tracer.image)

    matrix_pixels = reshape(tracer.image.pixels, (tracer.image.width, tracer.image.height))

    save(file_out_png, matrix_pixels')

    println("$(file_out_png) has been automatically written to disk.")

end

main()

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


text:
@testset "VAMP" begin
    # test different constructors
    xdim, zdim, n, K = 4, 2, 6, 2
    v = VAMP(zeros(Float32, xdim, K))
    @test length(Flux.params(v)) == 1
    v = VAMP(K, xdim)
    @test length(Flux.params(v)) == 1

    # test gpu compatibility
    cutype = typeof(gpu(zeros(Float32,2,2)))
    v = VAMP(zeros(Float32, xdim, K)) |> gpu
    @test typeof(v.pseudoinputs) == cutype # works even when gpu not available

    # test constructors
    x = f32(hcat(ones(xdim, n), -ones(xdim, n)) + randn(xdim, n*2)/15) |> gpu

    # mean constructor - K slightly perturbed means of X
    v = init_vamp_mean(3, cpu(x), 0.1f0)
    @test size(v.pseudoinputs) == (xdim,3)
    @test typeof(v.pseudoinputs) == Array{Float32,2}
    v = init_vamp_mean(3, cpu(x)) |> gpu
    @test typeof(v.pseudoinputs) == cutype

    # rand constructor - randomly selects K samples from X
    v = init_vamp_sample(3, cpu(x))
    @test size(v.pseudoinputs) == (xdim,3)
    @test typeof(v.pseudoinputs) == Array{Float32,2}
    v = init_vamp_sample(3, cpu(x)) |> gpu
    @test typeof(v.pseudoinputs) == cutype

    # also test it as a part of a VAE model, test trainability
    hdim = 50
    enc = f32(Chain(Dense(xdim, hdim, relu), SplitLayer(hdim, [zdim,zdim])))
    dec = f32(Chain(Dense(zdim, hdim, relu), SplitLayer(hdim, [xdim,1])))
    q = ConditionalMvNormal(enc)
    p = ConditionalMvNormal(dec)
    v = init_vamp_mean(K, cpu(x), 0.1f0)
    m = VAE(v, q, p) |> gpu # gpu will make a copy of v
    @test typeof(m.prior.pseudoinputs) == typeof(gpu(zeros(Float32,2,2))) # works even when gpu not available

    # extract original param vals
    vamp_params_init = get_params(m.prior)
    model_params_init = get_params(m)
    @test length(model_params_init) == 13

    # test trainability
    k = IMQKernel(0.1f0)
    λ = .1f0
    lossf(x) = λ*mmd_mean(m, x, k) .- mean(logpdf(m.decoder, x, rand(m.encoder, x)))
    opt = ADAM(0.0001)
    data = [(x,) for _ in 1:2];
    Flux.train!(lossf, Flux.params(m), data, opt)
    @test all(param_change(vamp_params_init, m.prior)) # did the vamp params change?
    @test all(param_change(model_params_init, m)) # did the model params change?
end


text:
@service S3

bucket_name = "aws-jl-test-issues---" * _now_formatted()
S3.create_bucket(bucket_name)

@testset "issue 223" begin
    # https://github.com/JuliaCloud/AWS.jl/issues/223
    body = "Hello World!"
    file_name = "contains spaces"

    try
        S3.put_object(bucket_name, file_name, Dict("body" => body))
        resp = S3.get_object(bucket_name, file_name)

        @test String(resp) == body
    finally
        S3.delete_object(bucket_name, file_name)
        S3.delete_bucket(bucket_name)
    end
end

@testset "issue 227" begin
    @testset "s3 public bucket" begin
        # https://github.com/JuliaCloud/AWS.jl/issues/227
        config = AWSConfig(creds=nothing)
        resp = S3.get_object("www.invenia.ca", "index.html"; aws_config=config)

        @test !isempty(resp)
    end

    @testset "s3 private bucket" begin
        bucket_name = "aws-jl-test-issues---" * _now_formatted()
        file_name = "hello_world"

        try
            S3.create_bucket(bucket_name)
            S3.put_object(bucket_name, file_name)

            @test_throws AWSException S3.get_object(bucket_name, file_name; aws_config=AWSConfig(creds=nothing))
        finally
            S3.delete_object(bucket_name, file_name)
            S3.delete_bucket(bucket_name)
        end
    end

    @testset "lambda" begin
        @service Lambda

        @test_throws NoCredentials Lambda.list_functions(;aws_config=AWSConfig(creds=nothing))
    end
end


text:
# Copyright (c) 2021 Idiap Research Institute, http://www.idiap.ch/
#  Niccolò Antonello <nantonel@idiap.ch>

export Path

mutable struct Path{W,D <: Integer,I,O,
                    IS <: AbstractDict{I,D}, OS <: AbstractDict{O,D}}
  ilabels::Vector{D} # here we only use Int
  olabels::Vector{D}
  weight::W
  isym::IS     # input  symbol table idx::Int => sym::I
  osym::OS     # output symbol table idx::Int => sym::O
end

"""
  `Path(isym[,osym], iolabel::Pair, w=one(TropicalWeight{Float32}))`

Construct a path with input and output symbol table `isym` and (optional) `osym` (see [`WFST`](@ref)).

`iolabel` is a `Pair` of vectors.

```julia
julia> isym = ["a","b","c","d"];

julia> osym = ["α","β","γ","δ"];

julia> W = ProbabilityWeight{Float32};

julia> p = Path(isym,osym,["a","b","c"] => ["α","β","γ"], one(W))
["a", "b", "c"] → ["α", "β", "γ"], w:1.0f0

```

The weight of a path of a WFST results from the multiplication (``\\otimes``) of the weights of the different arcs that are transversed. See e.g. [`get_paths`](@ref) to extract paths from a WFST.

```julia
julia> p * Path(isym,osym,["d"] => ["δ"], W(0.5))
["a", "b", "c", "d"] → ["α", "β", "γ", "δ"], w:0.5f0

```

"""
function Path(isym::Dict, osym::Dict, iolabels::Pair, w = one(TropicalWeight{Float32}))
  Path([isym[s] for s in iolabels.first],
       [osym[s] for s in iolabels.second], w, isym, osym)
end
Path(isym, iolabels::Pair, w = one(TropicalWeight{Float32})) =
Path(isym, isym, iolabels, w)

Path(isym::AbstractVector, osym::AbstractVector, iolabels::Pair, w) =
Path(
     Dict( s => i for (i,s) in enumerate(isym)),
     Dict( s => i for (i,s) in enumerate(osym)),
     iolabels,w)

get_weight(p::Path) = p.weight
get_ilabel(p::Path) = p.ilabels
get_olabel(p::Path) = p.olabels

export get_isequence
"""
  `get_isequence(p::Path)`

Returns the input sequence of the path `p`.
"""
function get_isequence(p::Path)
  idx2sym = get_iisym(p)
  return [idx2sym[idx] for idx in get_ilabel(p)]
end

export get_osequence
"""
  `get_osequence(p::Path)`

Returns the output sequence of the path `p`.
"""
function get_osequence(p::Path)
  idx2sym = get_iosym(p)
  return [idx2sym[idx] for idx in get_olabel(p)]
end

function update_path!(path::Path{W,D}, ilabel::D, olabel::D, w) where {W,D}
  if !iseps(ilabel)
    push!(path.ilabels, ilabel)
  end
  if !iseps(olabel)
    push!(path.olabels, olabel)
  end
  update_weight!(path, W(w))
  return path
end
function update_path(path,ilabel,olabel,w)
  p = deepcopy(path)
  update_path!(p,ilabel,olabel,w)
  return p
end
update_weight!(path::Path{W}, w::W) where {W} = path.weight *= w
update_weight(path::Path{W}, w::W) where {W} =
Path(path.ilabels,path.olabels,w*path.weight,path.isym,path.osym)

import Base: *

function *(p1::Path, p2::Path)
  Path([p1.ilabels; p2.ilabels], [p1.olabels;p2.olabels], p1.weight*p2.weight, p1.isym, p1.osym)
end

import Base: ==
function ==(p1::Path,p2::Path)
  if (length(p1.ilabels) == length(p2.ilabels)) &&
    (length(p1.olabels) == length(p2.olabels))
    return all(p1.ilabels .== p2.ilabels) && all(p1.olabels .== p2.olabels) && (p1.weight == p2.weight)
  else
    return false
  end
end


text:
struct ANSICode{T<:AbstractString}
    open::T
    close::T
end

const CODES = Dict(
    :default => ANSICode("\e[22m", "\e[22m"),
    :bold => ANSICode("\e[1m", "\e[22m"),
    :b => ANSICode("\e[1m", "\e[22m"),
    :dim => ANSICode("\e[2m", "\e[22m"),
    :italic => ANSICode("\e[3m", "\e[23m"),
    :i => ANSICode("\e[3m", "\e[23m"),
    :underline => ANSICode("\e[4m", "\e[24m"),
    :u => ANSICode("\e[4m", "\e[24m"),
    :blink => ANSICode("\e[5m", "\e[25m"),
    :inverse => ANSICode("\e[7m", "\e[27m"),
    :hidden => ANSICode("\e[8m", "\e[28m"),
    :striked => ANSICode("\e[9m", "\e[29m"),
    "black" => 0,
    "red" => 1,
    "green" => 2,
    "yellow" => 3,
    "blue" => 4,
    "magenta" => 5,
    "cyan" => 6,
    "white" => 7,
    "default" => 9,
)

const CODES_16BIT_COLORS = Dict(
    "bright_black" => 8,
    "bright_red" => 9,
    "bright_green" => 10,
    "bright_yellow" => 11,
    "bright_blue" => 12,
    "bright_magenta" => 13,
    "bright_cyan" => 14,
    "bright_white" => 15,
    "grey0" => 16,
    "navy_blue" => 17,
    "dark_blue" => 18,
    "blue3" => 20,
    "blue1" => 21,
    "dark_green" => 22,
    "deep_sky_blue4" => 25,
    "dodger_blue3" => 26,
    "dodger_blue2" => 27,
    "green4" => 28,
    "spring_green4" => 29,
    "turquoise4" => 30,
    "deep_sky_blue3" => 32,
    "dodger_blue1" => 33,
    "green3" => 40,
    "spring_green3" => 41,
    "dark_cyan" => 36,
    "light_sea_green" => 37,
    "deep_sky_blue2" => 38,
    "deep_sky_blue1" => 39,
    "spring_green2" => 47,
    "cyan3" => 43,
    "dark_turquoise" => 44,
    "turquoise2" => 45,
    "green1" => 46,
    "spring_green1" => 48,
    "medium_spring_green" => 49,
    "cyan2" => 50,
    "cyan1" => 51,
    "dark_red" => 88,
    "deep_pink4" => 125,
    "purple4" => 55,
    "purple3" => 56,
    "blue_violet" => 57,
    "orange4" => 94,
    "grey37" => 59,
    "gray37" => 59,
    "medium_purple4" => 60,
    "slate_blue3" => 62,
    "royal_blue1" => 63,
    "chartreuse4" => 64,
    "dark_sea_green4" => 71,
    "pale_turquoise4" => 66,
    "steel_blue" => 67,
    "steel_blue3" => 68,
    "cornflower_blue" => 69,
    "chartreuse3" => 76,
    "cadet_blue" => 73,
    "sky_blue3" => 74,
    "steel_blue1" => 81,
    "pale_green3" => 114,
    "sea_green3" => 78,
    "aquamarine3" => 79,
    "medium_turquoise" => 80,
    "chartreuse2" => 112,
    "sea_green2" => 83,
    "sea_green1" => 85,
    "aquamarine1" => 122,
    "dark_slate_gray2" => 87,
    "dark_magenta" => 91,
    "dark_violet" => 128,
    "purple" => 129,
    "light_pink4" => 95,
    "plum4" => 96,
    "medium_purple3" => 98,
    "slate_blue1" => 99,
    "yellow4" => 106,
    "wheat4" => 101,
    "grey53" => 102,
    "gray53" => 102,
    "light_slate_grey" => 103,
    "light_slate_gray" => 103,
    "medium_purple" => 104,
    "light_slate_blue" => 105,
    "dark_olive_green3" => 149,
    "dark_sea_green" => 108,
    "light_sky_blue3" => 110,
    "sky_blue2" => 111,
    "dark_sea_green3" => 150,
    "dark_slate_gray3" => 116,
    "sky_blue1" => 117,
    "chartreuse1" => 118,
    "light_green" => 120,
    "pale_green1" => 156,
    "dark_slate_gray1" => 123,
    "red3" => 160,
    "medium_violet_red" => 126,
    "magenta3" => 164,
    "dark_orange3" => 166,
    "indian_red" => 167,
    "hot_pink3" => 168,
    "medium_orchid3" => 133,
    "medium_orchid" => 134,
    "medium_purple2" => 140,
    "dark_goldenrod" => 136,
    "light_salmon3" => 173,
    "rosy_brown" => 138,
    "grey63" => 139,
    "gray63" => 139,
    "medium_purple1" => 141,
    "gold3" => 178,
    "dark_khaki" => 143,
    "navajo_white3" => 144,
    "grey69" => 145,
    "gray69" => 145,
    "light_steel_blue3" => 146,
    "light_steel_blue" => 147,
    "yellow3" => 184,
    "dark_sea_green2" => 157,
    "light_cyan3" => 152,
    "light_sky_blue1" => 153,
    "green_yellow" => 154,
    "dark_olive_green2" => 155,
    "dark_sea_green1" => 193,
    "pale_turquoise1" => 159,
    "deep_pink3" => 162,
    "magenta2" => 200,
    "hot_pink2" => 169,
    "orchid" => 170,
    "medium_orchid1" => 207,
    "orange3" => 172,
    "light_pink3" => 174,
    "pink3" => 175,
    "plum3" => 176,
    "violet" => 177,
    "light_goldenrod3" => 179,
    "tan" => 180,
    "misty_rose3" => 181,
    "thistle3" => 182,
    "plum2" => 183,
    "khaki3" => 185,
    "light_goldenrod2" => 222,
    "light_yellow3" => 187,
    "grey84" => 188,
    "gray84" => 188,
    "light_steel_blue1" => 189,
    "yellow2" => 190,
    "dark_olive_green1" => 192,
    "honeydew2" => 194,
    "light_cyan1" => 195,
    "red1" => 196,
    "deep_pink2" => 197,
    "deep_pink1" => 199,
    "magenta1" => 201,
    "orange_red1" => 202,
    "indian_red1" => 204,
    "hot_pink" => 206,
    "dark_orange" => 208,
    "salmon1" => 209,
    "light_coral" => 210,
    "pale_violet_red1" => 211,
    "orchid2" => 212,
    "orchid1" => 213,
    "orange1" => 214,
    "sandy_brown" => 215,
    "light_salmon1" => 216,
    "light_pink1" => 217,
    "pink1" => 218,
    "plum1" => 219,
    "gold1" => 220,
    "navajo_white1" => 223,
    "misty_rose1" => 224,
    "thistle1" => 225,
    "yellow1" => 226,
    "light_goldenrod1" => 227,
    "khaki1" => 228,
    "wheat1" => 229,
    "cornsilk1" => 230,
    "grey100" => 231,
    "gray100" => 231,
    "grey3" => 232,
    "gray3" => 232,
    "grey7" => 233,
    "gray7" => 233,
    "grey11" => 234,
    "gray11" => 234,
    "grey15" => 235,
    "gray15" => 235,
    "grey19" => 236,
    "gray19" => 236,
    "grey23" => 237,
    "gray23" => 237,
    "grey27" => 238,
    "gray27" => 238,
    "grey30" => 239,
    "gray30" => 239,
    "grey35" => 240,
    "gray35" => 240,
    "grey39" => 241,
    "gray39" => 241,
    "grey42" => 242,
    "gray42" => 242,
    "grey46" => 243,
    "gray46" => 243,
    "grey50" => 244,
    "gray50" => 244,
    "grey54" => 245,
    "gray54" => 245,
    "grey58" => 246,
    "gray58" => 246,
    "grey62" => 247,
    "gray62" => 247,
    "grey66" => 248,
    "gray66" => 248,
    "grey70" => 249,
    "gray70" => 249,
    "grey74" => 250,
    "gray74" => 250,
    "grey78" => 251,
    "gray78" => 251,
    "grey82" => 252,
    "gray82" => 252,
    "grey85" => 253,
    "gray85" => 253,
    "grey89" => 254,
    "gray89" => 254,
    "grey93" => 255,
    "gray93" => 255,
)

# ---------------------------------------------------------------------------- #
#                                    COLORS                                    #
# ---------------------------------------------------------------------------- #
const NAMED_COLORS = [
    "default", "black", "red", "green", "yellow", "blue", "magenta", "cyan", "white"
]

const COLORS_16b = vcat(
    [
        "bright_black",
        "bright_red",
        "bright_green",
        "bright_yellow",
        "bright_blue",
        "bright_magenta",
        "bright_cyan",
        "bright_white",
        "grey0",
        "navy_blue",
        "dark_blue",
        "blue3",
        "blue1",
        "dark_green",
        "deep_sky_blue4",
        "dodger_blue3",
        "dodger_blue2",
        "green4",
        "spring_green4",
        "turquoise4",
        "deep_sky_blue3",
        "dodger_blue1",
        "green3",
        "spring_green3",
        "dark_cyan",
        "light_sea_green",
        "deep_sky_blue2",
        "deep_sky_blue1",
        "spring_green2",
        "cyan3",
        "dark_turquoise",
        "turquoise2",
        "green1",
        "spring_green1",
        "medium_spring_green",
        "cyan2",
        "cyan1",
        "dark_red",
        "deep_pink4",
        "purple4",
        "purple3",
        "blue_violet",
        "orange4",
        "grey37",
        "gray37",
        "medium_purple4",
        "slate_blue3",
        "royal_blue1",
        "chartreuse4",
        "dark_sea_green4",
        "pale_turquoise4",
        "steel_blue",
        "steel_blue3",
        "cornflower_blue",
        "chartreuse3",
        "cadet_blue",
        "sky_blue3",
        "steel_blue1",
        "pale_green3",
        "sea_green3",
        "aquamarine3",
        "medium_turquoise",
        "chartreuse2",
        "sea_green2",
        "sea_green1",
        "aquamarine1",
        "dark_slate_gray2",
        "dark_magenta",
        "dark_violet",
        "purple",
        "light_pink4",
        "plum4",
        "medium_purple3",
        "slate_blue1",
        "yellow4",
        "wheat4",
        "grey53",
        "gray53",
        "light_slate_grey",
        "light_slate_gray",
        "medium_purple",
        "light_slate_blue",
        "dark_olive_green3",
        "dark_sea_green",
        "light_sky_blue3",
        "sky_blue2",
        "dark_sea_green3",
        "dark_slate_gray3",
        "sky_blue1",
        "chartreuse1",
        "light_green",
        "pale_green1",
        "dark_slate_gray1",
        "red3",
        "medium_violet_red",
        "magenta3",
        "dark_orange3",
        "indian_red",
        "hot_pink3",
        "medium_orchid3",
        "medium_orchid",
        "medium_purple2",
        "dark_goldenrod",
        "light_salmon3",
        "rosy_brown",
        "grey63",
        "gray63",
        "medium_purple1",
        "gold3",
        "dark_khaki",
        "navajo_white3",
        "grey69",
        "gray69",
        "light_steel_blue3",
        "light_steel_blue",
        "yellow3",
        "dark_sea_green2",
        "light_cyan3",
        "light_sky_blue1",
        "green_yellow",
        "dark_olive_green2",
        "dark_sea_green1",
        "pale_turquoise1",
        "deep_pink3",
        "magenta2",
        "hot_pink2",
        "orchid",
        "medium_orchid1",
        "orange3",
        "light_pink3",
        "pink3",
        "plum3",
        "violet",
        "light_goldenrod3",
        "tan",
        "misty_rose3",
        "thistle3",
        "plum2",
        "khaki3",
        "light_goldenrod2",
        "light_yellow3",
        "grey84",
        "gray84",
        "light_steel_blue1",
        "yellow2",
        "dark_olive_green1",
        "honeydew2",
        "light_cyan1",
        "red1",
        "deep_pink2",
        "deep_pink1",
        "magenta1",
        "orange_red1",
        "indian_red1",
        "hot_pink",
        "dark_orange",
        "salmon1",
        "light_coral",
        "pale_violet_red1",
        "orchid2",
        "orchid1",
        "orange1",
        "sandy_brown",
        "light_salmon1",
        "light_pink1",
        "pink1",
        "plum1",
        "gold1",
        "navajo_white1",
        "misty_rose1",
        "thistle1",
        "yellow1",
        "light_goldenrod1",
        "khaki1",
        "wheat1",
        "cornsilk1",
        "grey100",
        "gray100",
        "grey3",
        "gray3",
        "grey7",
        "gray7",
        "grey11",
        "gray11",
        "grey15",
        "gray15",
        "grey19",
        "gray19",
        "grey23",
        "gray23",
        "grey27",
        "gray27",
        "grey30",
        "gray30",
        "grey35",
        "gray35",
        "grey39",
        "gray39",
        "grey42",
        "gray42",
        "grey46",
        "gray46",
        "grey50",
        "gray50",
        "grey54",
        "gray54",
        "grey58",
        "gray58",
        "grey62",
        "gray62",
        "grey66",
        "gray66",
        "grey70",
        "gray70",
        "grey74",
        "gray74",
        "grey78",
        "gray78",
        "grey82",
        "gray82",
        "grey85",
        "gray85",
        "grey89",
        "gray89",
        "grey93",
        "gray93",
    ],
    [string(i) for i in 1:255],
)

# ---------------------------------------------------------------------------- #
#                                     MODES                                    #
# ---------------------------------------------------------------------------- #
const NAMED_MODES = [
    "default",
    "bold",
    "b",
    "dim",
    "italic",
    "i",
    "underline",
    "u",
    "blink",
    "inverse",
    "hidden",
    "striked",
]


text:
# using Revise, Plots
using OrdinaryDiffEq, ForwardDiff, Test
	using PseudoArcLengthContinuation, LinearAlgebra, Parameters, Setfield
	const PALC = PseudoArcLengthContinuation

norminf = x -> norm(x, Inf)

function Fsl!(f, u, p, t)
	@unpack r, μ, ν, c3, c5 = p
	u1 = u[1]
	u2 = u[2]

	ua = u1^2 + u2^2

	f[1] = r * u1 - ν * u2 - ua * (c3 * u1 - μ * u2) - c5 * ua^2 * u1
	f[2] = r * u2 + ν * u1 - ua * (c3 * u2 + μ * u1) - c5 * ua^2 * u2
	return f
end

Fsl(x, p) = Fsl!(similar(x), x, p, 0.)

Fode(f, x, p, t) = Fsl!(f, x, p, t)
dFsl(x, dx, p) = ForwardDiff.derivative(t -> Fsl(x .+ t .* dx, p), 0.)
# JacOde(J, x, p, t) = copyto!(J, Jsl(x, p))

par_sl = (r = 0.5, μ = 0., ν = 1.0, c3 = 1.0, c5 = 0.0,)
u0 = [.001, .001]

function FslMono!(f, x, p, t)
	u = x[1:2]
	du = x[3:4]
	Fsl!(f[1:2], u, p, t)
	f[3:4] .= dFsl(u, du, p)
end

####################################################################################################
par_hopf = (@set par_sl.r = 0.1)
####################################################################################################
# continuation
optconteq = ContinuationPar(ds = -0.01, detectBifurcation = 2, pMin = -0.5)
br, _ = continuation(Fsl, u0, par_hopf, (@lens _.r), optconteq)
####################################################################################################
prob = ODEProblem(Fode, u0, (0., 100.), par_hopf)
probMono = ODEProblem(FslMono!, vcat(u0, u0), (0., 100.), par_hopf)
####################################################################################################
sol = @time solve(probMono, Tsit5(), abstol =1e-9, reltol=1e-6)
sol = @time solve(prob, Tsit5(), abstol =1e-9, reltol=1e-6)
# plot(sol[1,:], sol[2,:])
####################################################################################################
section(x) = x[1]
# standard simple shooting
M = 1
dM = 1
_pb = ShootingProblem(Fsl, par_hopf, prob, Rodas4(), 1, section; rtol = 1e-9)

initpo = [0.3, 0., 6.]
res = @time _pb(initpo, par_hopf)

# test of the differential of thew shooting method

_dx = rand(3)
resAD = ForwardDiff.derivative(z -> _pb(initpo .+ z .* _dx, par_hopf), 0.)
resFD = (_pb(initpo .+ 1e-8 .* _dx, par_hopf) - _pb(initpo, par_hopf)) * 1e8
resAN = _pb(initpo, par_hopf, _dx; δ = 1e-8)
@test norm(resAN - resFD, Inf) < 1e-1
@test norm(resAN - resAD, Inf) < 1e-1
####################################################################################################
# test shooting interface
_pb = ShootingProblem(Fsl, par_hopf, prob, Rodas4(), [initpo[1:end-1]]; rtol = 1e-9)
res = @time _pb(initpo, par_hopf)
res = @time _pb(initpo, par_hopf, initpo)

_pb2 = ShootingProblem(Fsl, par_hopf, prob, Rodas4(), probMono, Rodas4(autodiff=false), [initpo[1:end-1]]; rtol = 1e-9)
res = @time _pb2(initpo, par_hopf)
res = @time _pb2(initpo, par_hopf, initpo)
PALC.isSimple(_pb2)
####################################################################################################
# we test this using Newton - Continuation
ls = GMRESIterativeSolvers(tol = 1e-5, N = length(initpo))
optn = NewtonPar(verbose = false, tol = 1e-9,  maxIter = 20, linsolver = ls)
deflationOp = PALC.DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [zeros(3)])
outpo, _ = @time PALC.newton(_pb,
	initpo, par_hopf,
	optn,
	normN = norminf)

getPeriod(_pb, outpo)
PALC.getAmplitude(_pb, outpo, par_hopf)

opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.01, ds= -0.01, pMax = 4.0, maxSteps = 30, detectBifurcation = 1, nev = 2, newtonOptions = @set optn.tol = 1e-7)#
	br_pok2, upo , _= @time PALC.continuation(
		_pb,
		outpo, par_hopf, (@lens _.r),
		opts_po_cont;
		tangentAlgo = BorderedPred(),
		verbosity = 0,
		plot = false,
		# plotSolution = (x, p; kwargs...) -> plot!(x[1:end-1]; kwargs...),
		printSolution = (u, p) -> norm(u[1:2]), normC = norminf)
# plot(br_pok2)
####################################################################################################
# test automatic branch switching
JFsl(x, p) = ForwardDiff.jacobian(t -> Fsl(t, p), x)
d1Fsl(x, p, dx) = ForwardDiff.derivative(t -> Fsl(x .+ t .* dx, p), 0.)
d2Fsl(x, p, dx1, dx2) = ForwardDiff.derivative(t -> d1Fsl(x .+ t .* dx2, p, dx1), 0.)
d3Fsl(x, p, dx1, dx2, dx3) = ForwardDiff.derivative(t -> d2Fsl(x .+ t .* dx3, p, dx1, dx2), 0.)
jet = (Fsl, JFsl, d2Fsl, d3Fsl)

br_pok2, _ = continuation(jet...,br,1, opts_po_cont, ShootingProblem(1, par_hopf, prob, Rodas4());printSolution = (u, p) -> norm(u[1:2]), normC = norminf)

# test matrix-free computation of floquet coefficients
eil = EigKrylovKit(dim = 2, x₀=rand(2))
opts_po_contMF = @set opts_po_cont.newtonOptions.eigsolver = eil
br_pok2, _ = continuation(jet...,br,1, opts_po_contMF, ShootingProblem(1, par_hopf, prob, Rodas4());printSolution = (u, p) -> norm(u[1:2]), normC = norminf)
####################################################################################################
# Single Poincaré Shooting with hyperplane parametrization
normals = [[-1., 0.]]
centers = [zeros(2)]

probPsh = PoincareShootingProblem(2, par_hopf,prob, Rodas4(),probMono, Rodas4())
probPsh = PoincareShootingProblem(2, par_hopf,prob, Rodas4())

probPsh = PoincareShootingProblem(Fsl, par_hopf,
		prob, Rodas4(),
		probMono, Rodas4(),
		normals, centers; rtol = 1e-8)

hyper = probPsh.section

initpo_bar = PALC.R(hyper, [0,0.4], 1)

PALC.E(hyper, [1.0,], 1)
initpo_bar = [0.4]

probPsh(initpo_bar, par_hopf)

ls = GMRESIterativeSolvers(tol = 1e-7, N = length(initpo_bar), maxiter = 500, verbose = false)
	eil = EigKrylovKit(dim = 1, x₀=rand(1))
	optn = NewtonPar(verbose = false, tol = 1e-8,  maxIter = 140, linsolver = ls, eigsolver = eil)
	deflationOp = PALC.DeflationOperator(2.0, (x,y) -> dot(x, y), 1.0, [zero(initpo_bar)])
	outpo, _ = @time PALC.newton(probPsh,
			initpo_bar, par_hopf,
			optn; normN = norminf)
	println("--> Point on the orbit = ", PALC.E(hyper, outpo, 1))

getPeriod(probPsh, outpo, par_hopf)
PALC.getAmplitude(probPsh, outpo, par_hopf)

probPsh = PoincareShootingProblem(Fsl, par_hopf,
		prob, Rodas4(),
		# probMono, Rodas4(autodiff=false),
		normals, centers; rtol = 1e-8)

opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.015, ds= -0.01, pMax = 4.0, maxSteps = 50, newtonOptions = setproperties(optn;tol = 1e-9, eigsolver = eil), detectBifurcation = 1)
	br_pok2, upo , _= @time PALC.continuation(
		probPsh, outpo, par_hopf, (@lens _.r),
		opts_po_cont; verbosity = 0,
		tangentAlgo = BorderedPred(),
		plot = false,
		# plotSolution = (x, p;kwargs...) -> plot!(x; kwargs...),
		printSolution = (u, p) -> norm(u), normC = norminf)
# plot(br_pok2)
####################################################################################################
# normals = [[-1., 0.], [1, -1]]
# centers = [zeros(2), zeros(2)]
# initpo_bar = [1.04, -1.04/√2]

normals = [[-1., 0.], [1, 0]]
centers = [zeros(2), zeros(2)]
initpo_bar = [0.2, -0.2]

probPsh = PALC.PoincareShootingProblem(Fsl, par_hopf, prob, Tsit5(), normals, centers; rtol = 1e-6)
# version with analytical jacobian
probPsh2 = PALC.PoincareShootingProblem(Fsl, par_hopf, prob, Tsit5(), normals, centers; rtol = 1e-6, δ = 0)

hyper = probPsh.section

probPsh(initpo_bar, par_hopf; verbose = true)

ls = GMRESIterativeSolvers(tol = 1e-5, N = length(initpo_bar), maxiter = 500, verbose = false)
	eil = EigKrylovKit(dim = 1, x₀=rand(1))
	optn = NewtonPar(verbose = false, tol = 1e-9,  maxIter = 140, linsolver = ls, eigsolver = eil)
	deflationOp = PALC.DeflationOperator(2.0, (x,y) -> dot(x, y), 1.0, [zero(initpo_bar)])
	outpo, _ = @time PALC.newton(probPsh2, initpo_bar, par_hopf, optn; normN = norminf)
	outpo, _ = @time PALC.newton(probPsh, initpo_bar, par_hopf, optn; normN = norminf)
println("--> Point on the orbit = ", PALC.E(hyper, [outpo[1]], 1), PALC.E(hyper, [outpo[2]], 2))

getPeriod(probPsh, outpo, par_hopf)

opts_po_cont = ContinuationPar(dsmin = 0.0001, dsmax = 0.025, ds= -0.01, pMax = 4.0, maxSteps = 50, newtonOptions = (@set optn.tol = 1e-9), detectBifurcation = 1, nev = 2)
	br_pok2, upo , _= @time PALC.continuation(probPsh, outpo, par_hopf, (@lens _.r),
		opts_po_cont; verbosity = 0,
		tangentAlgo = BorderedPred(),
		plot = false,
		# plotSolution = (x, p;kwargs...) -> plot!(x; kwargs...),
		printSolution = (u, p) -> norm(u), normC = norminf)

####################################################################################################
normals = [[-1., 0.], [1, 0], [0, 1]]
centers = [zeros(2), zeros(2), zeros(2)]
initpo = [[0., 0.4], [0, -.3], [0.3, 0]]

probPsh = PoincareShootingProblem(Fsl, par_hopf, prob, Tsit5(), normals, centers; rtol = 1e-6)

hyper = probPsh.section
initpo_bar = reduce(vcat, [PALC.R(hyper, initpo[ii], ii) for ii in eachindex(centers)])

probPsh(initpo_bar, par_hopf; verbose = true)

ls = GMRESIterativeSolvers(tol = 1e-7, N = length(initpo_bar), maxiter = 10, verbose = false)
	optn = NewtonPar(verbose = false, tol = 1e-9,  maxIter = 50, linsolver = ls, eigsolver = eil)
	outpo, _ = @time PALC.newton(probPsh, initpo_bar, par_hopf, optn; normN = norminf)

for ii=1:length(normals)
	@show PALC.E(hyper, [outpo[ii]], ii)
end

getPeriod(probPsh, outpo, par_hopf)

opts_po_cont = ContinuationPar(dsmin = 0.0001, dsmax = 0.025, ds= -0.005, pMax = 4.0, maxSteps = 50, newtonOptions = setproperties(optn; tol = 1e-9), detectBifurcation = 1)
	br_hpsh, upo , _= @time PALC.continuation(probPsh, outpo, par_hopf, (@lens _.r),
		opts_po_cont; verbosity = 0, plot = false,
		# plotSolution = (x, p;kwargs...) -> plot!(x, subplot=3),
		printSolution = (u, p) -> norm(u), normC = norminf)
# plot(br_hpsh)
####################################################################################################
# test automatic branch switching
# calls with jacobian computed with finite differences
br_hpsh, _ = continuation(jet...,br,1, opts_po_cont, PoincareShootingProblem(1, par_hopf, prob, Rodas4());printSolution = (u, p) -> norm(u), normC = norminf)

br_hpsh, _ = continuation(jet...,br,1, opts_po_cont, PoincareShootingProblem(2, par_hopf, prob, Rodas4());printSolution = (u, p) -> norm(u), normC = norminf)

# calls with analytical jacobians
br_hpsh, _ = continuation(jet...,br,1, (@set opts_po_cont.detectBifurcation = 0), PoincareShootingProblem(1, par_hopf, prob, Rodas4(), probMono, Rodas4(),);printSolution = (u, p) -> norm(u), normC = norminf)

br_hpsh, _ = continuation(jet...,br,1, (@set opts_po_cont.detectBifurcation = 0), PoincareShootingProblem(2, par_hopf, prob, Rodas4(), probMono, Rodas4(),);printSolution = (u, p) -> norm(u), normC = norminf)

# plot(br_hpsh, vars = (:step, :p))
# plot(br_hpsh)


text:
module Impl

import Adapt
using Base: Dims
using Transducers
using Transducers: @return_if_reduced, next, complete
using Setfield: @set  # using Setfield instead of Accessors for older Julia

import Transducers: executor_type

using ..UnionArrays: Abstract, UnionArrays, buffereltypefor

include("utils.jl")
include("eltypespec.jl")
include("vectors.jl")
include("arrays.jl")

end  # module


text:
using PolynomialBasis
using CartesianMesh
using CutCellDG
using Test
include("../useful_routines.jl")

x0 = [-1.0]
meshwidth = [2.0]
nelements = [3]
polyorder = 3
solverbasis = LagrangeTensorProductBasis(1, polyorder)
refpoints = interpolation_points(solverbasis)

dgmesh = CutCellDG.DGMesh(x0,meshwidth,nelements,refpoints)
nc = CutCellDG.nodal_connectivity(dgmesh)
testnc = reshape(1:12,4,:)
@test allequal(nc,testnc)


text:
# This file is a part of KCenters.jl
# License is Apache 2.0: https://www.apache.org/licenses/LICENSE-2.0.txt

export partition, knr, sequence, invindex

"""
    partition(callback::Function, objects::AbstractVector{T}, refs::AbstractSearchContext; k::Int=1) where T

Groups items in `objects` using a nearest neighbor rule over `refs`.
The output is controlled using a callback function. The call is performed in `objects` order.

- `callback` is a function that is called for each `(objID, refItem)`
- `objects` is the input dataset
- `dist` a distance function ``(T, T) \\rightarrow \\mathbb{R}``
- `refs` the list of references
- `k` specifies the number of nearest neighbors to use
- `indexclass` specifies the kind of index to be used, a function receiving `(refs, dist)` as arguments,
    and returning the new metric index

Please note that each object can be related to more than one group ``k > 1`` (default ``k=1``).

"""
function partition(callback::Function, objects::AbstractVector{T}, refs::AbstractSearchContext; k::Int=1) where T
    res = KnnResult(k)
    for i in 1:length(objects)
        empty!(res)
        callback(i, search(refs, objects[i], res))
    end
end

"""
    invindex(objects::AbstractVector{T}, refs::AbstractSearchContext; k::Int=1) where T

Creates an inverted index from references to objects.
So, an object ``u`` is in ``r``'s posting list iff ``r``
is among the ``k`` nearest references of ``u``.

"""
function invindex(objects::AbstractVector{T}, refs::AbstractSearchContext; k::Int=1) where T
    π = [Vector{Int}() for i in 1:length(refs.db)]
    # partition((i, p) -> push!(π[p.id], i), dist, objects, refs, k=k)
    partition(objects, refs, k=k) do i, res
        for p in res
            push!(π[p.id], i)
        end
    end
    π
end

"""
    sequence(objects::AbstractVector{T}, refs::AbstractSearchContext) where T

Computes the nearest reference of each item in the dataset and return it as a sequence of identifiers
"""
function sequence(objects::AbstractVector{T}, refs::AbstractSearchContext) where T
    s = Vector{Int}(undef, length(objects))
    partition(objects, refs) do i, res
        s[i] = first(res).id
    end
    s
end

"""
    knr(objects::AbstractVector{T}, refs::AbstractSearchContext) where T

Computes an array of k-nearest neighbors for `objects`
"""
function knr(objects::AbstractVector{T}, refs::AbstractSearchContext) where T
    s = Vector{Vector{Int}}(undef, length(objects))
    partition(objects, refs) do i, res
        s[i] = [p.id for p in res]
    end
    s
end


text:
function ab_premixing(::Type{T},mixing,Tc,pc,kij) where T <: ABCubicModel
    Ωa, Ωb = ab_consts(T)
    _Tc = Tc.values
    _pc = pc.values
    a = epsilon_LorentzBerthelot(SingleParam(pc, @. Ωa*R̄^2*_Tc^2/_pc),kij)
    b = sigma_LorentzBerthelot(SingleParam(pc, @. Ωb*R̄*_Tc/_pc))
    return a,b
end

function cubic_ab(model::ABCubicModel,V,T,z=SA[1.0],n=sum(z))
    invn2 = (one(n)/n)^2
    a = model.params.a.values
    b = model.params.b.values
    α = @f(α_function,model.alpha)
    c = @f(translation,model.translation)
    if length(z)>1
        ā,b̄,c̄ = @f(mixing_rule,model.mixing,α,a,b,c)
    else
        ā = a[1,1]*α[1]
        b̄ = b[1,1]
        c̄ = c[1,1]
    end
    return ā ,b̄, c̄
end

function second_virial_coefficient(model::ABCubicModel,T,z = SA[1.0])
    #@info "fast shortcut"
    a,b,c = cubic_ab(model,1e-4,T,z)
    return b-a/(R̄*T)
end

function lb_volume(model::CubicModel,z = SA[1.0])
    V = 1e-5
    T = 0.
    n = sum(z)
    invn = one(n)/n
    b = model.params.b.values
    c = @f(translation,model.translation)
    b̄ = dot(z,Symmetric(b),z)*invn*invn
    c̄ = dot(z,c)*invn
    return b̄-c̄
end
#dont use αa, just a, to avoid temperature dependence
function T_scale(model::CubicModel,z=SA[1.0])
    n = sum(z)
    invn2 = one(n)/(n*n)
    Ωa,Ωb = ab_consts(model)
    _a = model.params.a.values
    _b = model.params.b.values
    a = dot(z, Symmetric(_a), z)*invn2/Ωa
    b = dot(z, Symmetric(_b), z)*invn2/Ωb
    return a/b/R̄
end

function p_scale(model::CubicModel,z=SA[1.0])
    n = sum(z)
    invn2 = (1/n)^2
    Ωa,Ωb = ab_consts(model)
    _a = model.params.a.values
    _b = model.params.b.values
    a = invn2*dot(z, Symmetric(_a), z)/Ωa
    b = invn2*dot(z, Symmetric(_b), z)/Ωb
    return a/ (b^2) # Pc mean
end

function x0_crit_pure(model::CubicModel)
    lb_v = lb_volume(model)
    (1.0, log10(lb_v/0.3))
end

function volume(model::ABCubicModel,p,T,z=SA[1.0];phase=:unknown,threaded=false)
    lb_v   =lb_volume(model,z)
    RTp = R̄*T/p
    _poly,c̄ = cubic_poly(model,p,T,z)
    sols = Solvers.roots3(_poly)
    function imagfilter(x)
        absx = abs(imag(x))
        return absx < 8*eps(typeof(absx))
    end
    x1,x2,x3 = sols
    sols = (x1,x2,x3)
    xx = (x1,x2,x3)
    isreal = imagfilter.(xx)
    vvv = extrema(real.(xx))
    zl,zg = vvv
    vvl,vvg = RTp*zl,RTp*zg
    err() = @error("model $model Failed to converge to a volume root at pressure p = $p [Pa], T = $T [K] and compositions = $z")
    if sum(isreal) == 3 #3 roots
        vg = vvg
        _vl = vvl
        vl = ifelse(_vl>lb_v,_vl,vg) #catch case where solution is unphysical
    elseif  sum(isreal) == 1
        i = findfirst(imagfilter,sols)
        vl = real(sols[i])*RTp
        vg = real(sols[i])*RTp
    elseif  sum(isreal) == 0
        #try to use the default volume solver
        V0 = x0_volume(model,p,T,z;phase)
        v = _volume_compress(model,p,T,V0)
        isnan(v) && err()
        return v
    end

    function gibbs(v)
        _df,f =  ∂f(model,v,T,z)
        dv,dt = _df
        if abs((p+dv)/p) > 0.03
            return one(dv)/zero(dv)
        else
            return f + p*v
        end
    end
    #this catches the supercritical phase as well
    if vl ≈ vg
        return vl-c̄
    end

    if is_liquid(phase)
        return vl-c̄
    elseif is_vapour(phase)
        return vg-c̄
    else
        gg = gibbs(vg-c̄)
        gl = gibbs(vl-c̄)
        return ifelse(gg<gl,vg-c̄,vl-c̄)
    end
end

function ab_consts(model::CubicModel)
    return ab_consts(typeof(model))
end


text:
module DegreesOfFreedom

using Printf: @printf,@sprintf
using StaticArrays: SVector
using LaTeXStrings: latexstring
using ..Spatials: PID,AbstractBond
using ...Interfaces: rank,dimension
using ...Prerequisites: Float,decimaltostr,rawtype
using ...Prerequisites.TypeTraits: efficientoperations
using ...Prerequisites.CompositeStructures: CompositeDict
using ...Mathematics.VectorSpaces: VectorSpace
using ...Mathematics.AlgebraOverFields: SimpleID,ID,Element,Elements

import ..Spatials: pidtype,rcoord,icoord
import ...Interfaces: reset!,update!,sequence
import ...Mathematics.AlgebraOverFields: idpropertynames

export IID,Index,pid,iidtype,iid
export IndexToTuple,DirectIndexToTuple,directindextotuple,FilteredAttributes
export Internal,IDFConfig,Table
export LaTeX,OID,Operator,Operators,isHermitian,twist
export coordpresent,coordabsent
export latexformat,script,oidtype,otype
export Boundary

"""
    IID

The id of an internal degree of freedom.
"""
abstract type IID <: SimpleID end

"""
    Internal

The whole internal degrees of freedom at a single point.
"""
abstract type Internal{I<:IID} <: VectorSpace{I} end

"""
    show(io::IO,i::Internal)

Show an internal.
"""
Base.show(io::IO,i::Internal)=@printf io "%s(%s)" i|>typeof|>nameof join(("$name=$(getfield(i,name))" for name in i|>typeof|>fieldnames),",")

"""
    Index{P,I}

The complete index of a degree of freedom, which consist of the spatial part and the internal part.
"""
abstract type Index{P<:PID,I<:IID} <: SimpleID end

"""
    (INDEX::Type{<:Index})(pid::PID,iid::IID) -> INDEX

Get the corresponding index from a pid and an iid.
"""
@generated function (INDEX::Type{<:Index})(pid::PID,iid::IID)
    INDEX<:UnionAll ? :(INDEX(convert(Tuple,pid)...,convert(Tuple,iid)...)) : :(rawtype(INDEX)(convert(Tuple,pid)...,convert(Tuple,iid)...))
end

"""
    pidtype(index::Index)
    pidtype(::Type{<:Index{P}}) where {P<:PID}

Get the type of the spatial part of an index.
"""
pidtype(index::Index)=index|>typeof|>pidtype
pidtype(::Type{<:Index{P}}) where {P<:PID}=P

"""
    pid(index::Index) -> PID

Get the spatial part of an index.
"""
@generated function pid(index::Index)
    exprs=[:(getfield(index,$i)) for i=1:fieldcount(index|>pidtype)]
    return :(rawtype(pidtype(index))($(exprs...)))
end

"""
    iidtype(index::Index)
    iidtype(::Type{<:Index{<:PID,I}}) where {I<:IID}

Get the type of the internal part of an index.
"""
iidtype(index::Index)=index|>typeof|>iidtype
iidtype(::Type{<:Index{<:PID,I}}) where {I<:IID}=I

"""
    iid(index::Index) -> IID

Get the internal part of an index.
"""
@generated function iid(index::Index)
    exprs=[:(getfield(index,$i)) for i=fieldcount(index|>pidtype)+1:fieldcount(index)]
    return :(rawtype(iidtype(index))($(exprs...)))
end

"""
    adjoint(index::Index) -> typeof(index)

Get the adjoint of an index.
"""
Base.adjoint(index::Index)=typeof(index)(index|>pid,index|>iid|>adjoint)

"""
    union(::Type{P},::Type{I}) where {P<:PID,I<:IID}

Combine a concrete `PID` type and a concrete `IID` type to a concrete `Index` type.
"""
Base.union(::Type{P},::Type{I}) where {P<:PID,I<:IID}=Index{P,I}

"""
    IndexToTuple

The rules for converting an index to a tuple.

As a function, every instance should accept only one positional argument, i.e. the index to be converted to a tuple.
"""
abstract type IndexToTuple <:Function end
Base.:(==)(itt1::T,itt2::T) where T<:IndexToTuple = ==(efficientoperations,itt1,itt2)
Base.isequal(itt1::T,itt2::T) where T<:IndexToTuple = isequal(efficientoperations,itt1,itt2)

"""
    DirectIndexToTuple

Direct index to tuple.
"""
struct DirectIndexToTuple <: IndexToTuple end

"""
    (indextotuple::DirectIndexToTuple)(index::Index) -> Tuple

Convert an index to tuple directly.
"""
(indextotuple::DirectIndexToTuple)(index::Index)=convert(Tuple,index)

"""
    directindextotuple

Indicate that the conversion from an index to a tuple is direct.
"""
const directindextotuple=DirectIndexToTuple()

"""
    FilteredAttributes(::Type{I}) where I<:Index

A method that converts an arbitary index to a tuple, by iterating over the selected attributes in a specific order.
"""
struct FilteredAttributes{N} <: IndexToTuple
    attributes::NTuple{N,Symbol}
end
FilteredAttributes(attrs::Symbol...)=FilteredAttributes(attrs)
FilteredAttributes(::Type{I}) where I<:Index=FilteredAttributes(I|>fieldnames)

"""
    length(indextotuple::FilteredAttributes) -> Int
    length(::Type{<:FilteredAttributes{N}}) where N -> Int

Get the length of the filtered attributes.
"""
Base.length(indextotuple::FilteredAttributes)=indextotuple|>typeof|>length
Base.length(::Type{<:FilteredAttributes{N}}) where N=N

"""
    (indextotuple::FilteredAttributes)(index::Index) -> Tuple

Convert an index to tuple by the "filtered attributes" method.
"""
@generated function (indextotuple::FilteredAttributes)(index::Index)
    exprs=[:(getfield(index,indextotuple.attributes[$i])) for i=1:length(indextotuple)]
    return Expr(:tuple,exprs...)
end

"""
    filter(f::Function,indextotuple::FilteredAttributes)

Filter the attributes of a "filtered attributes" method.
"""
Base.filter(f::Function,indextotuple::FilteredAttributes)=FilteredAttributes(Tuple(attr for attr in indextotuple.attributes if f(attr)))

"""
    IDFConfig{I}(map::Function,pids::Union{AbstractVector{<:PID},Tuple{}}=()) where I<:Internal

Configuration of the internal degrees of freedom at a lattice.

`map` maps a `PID` to an `Internal`.
"""
struct IDFConfig{I<:Internal,M<:Function,P<:PID} <: CompositeDict{P,I}
    map::M
    contents::Dict{P,I}
end
function IDFConfig{I}(map::Function,pids::Union{AbstractVector{<:PID},Tuple{}}=()) where I<:Internal
    contents=Dict{pids|>eltype,I}()
    for pid in pids
        contents[pid]=map(pid)
    end
    IDFConfig(map,contents)
end

"""
    reset!(config::IDFConfig,pids) -> IDFConfig

Reset the idfconfig with new pids.
"""
function reset!(config::IDFConfig,pids)
    empty!(config)
    for pid in pids
        config[pid]=config.map(pid)
    end
    config
end

"""
    Table{I}(by::IndexToTuple) where I<:Index

Index-sequence table.
"""
struct Table{I<:Index,B<:IndexToTuple} <: CompositeDict{I,Int}
    by::B
    contents::Dict{I,Int}
end
Table{I}(by::IndexToTuple) where I<:Index=Table(by,Dict{I,Int}())

"""
    Table(indices::AbstractVector{<:Index},by::IndexToTuple=directindextotuple) -> Table

Convert an sequence of indices to the corresponding index-sequence table.

The input indices will be converted to tuples by the `by` function with the duplicates removed. The resulting unique tuples are sorted, which determines the sequence of the input indices. Note that two indices have the same sequence if their converted tupels are equal to each other.
"""
function Table(indices::AbstractVector{<:Index},by::IndexToTuple=directindextotuple)
    tuples=[by(index) for index in indices]
    permutation=sortperm(tuples,alg=Base.Sort.QuickSort)
    result=Table{indices|>eltype}(by)
    count=1
    for i=1:length(tuples)
        i>1 && tuples[permutation[i]]!=tuples[permutation[i-1]] && (count+=1)
        result[indices[permutation[i]]]=count
    end
    result
end

"""
    Table(config::IDFConfig,by::IndexToTuple=directindextotuple) -> Table

Get the index-sequence table of the whole internal Hilbert spaces at a lattice.
"""
function Table(config::IDFConfig,by::IndexToTuple=directindextotuple)
    result=union(config|>keytype,config|>valtype|>eltype)[]
    for (pid,internal) in config
        for iid in internal
            push!(result,(result|>eltype)(pid,iid))
        end
    end
    Table(result,by)
end

"""
    union(tables::Table...) -> Table

Unite several index-sequence tables.

See [`Table`](@ref) for more details.
"""
function Base.union(tables::Table...)
    @assert mapreduce(table->table.by,==,tables)
    indices=(tables|>eltype|>keytype)[]
    for table in tables
        for index in keys(table)
            push!(indices,index)
        end
    end
    Table(indices,tables[1].by)
end

"""
    reverse(table::Table) -> Dict{Int,Set{<:Index}}

Convert an index-sequence table to a sequence-indices table.

Since different indices may correspond to the same sequence, the reverse is a one-to-many map.
"""
function Base.reverse(table::Table)
    result=Dict{Int,Set{table|>keytype}}()
    for (index,seq) in table
        haskey(result,seq) || (result[seq]=Set{table|>keytype}())
        push!(result[seq],index)
    end
    result
end

"""
    reset!(table::Table,indices::AbstractVector{<:Index}) -> Table
    reset!(table::Table,config::IDFConfig) -> Table

Reset a table.
"""
function reset!(table::Table,indices::AbstractVector{<:Index})
    empty!(table)
    tuples=[table.by(index) for index in indices]
    permutation=sortperm(tuples,alg=Base.Sort.QuickSort)
    count=1
    for i=1:length(tuples)
        i>1 && tuples[permutation[i]]!=tuples[permutation[i-1]] && (count+=1)
        table[indices[permutation[i]]]=count
    end
    table
end
function reset!(table::Table,config::IDFConfig)
    indices=union(config|>keytype,config|>valtype|>eltype)[]
    for (pid,internal) in config
        for iid in internal
            push!(indices,(indices|>eltype)(pid,iid))
        end
    end
    reset!(table,indices)
end

"""
    LaTeX{SP,SB}(body) where {SP,SB}

LaTeX string representation.
"""
struct LaTeX{SP,SB,B}
    body::B
    function LaTeX{SP,SB}(body=nothing) where {SP,SB}
        @assert isa(SP,Tuple{Vararg{Symbol}}) && isa(SB,Tuple{Vararg{Symbol}}) "LaTeX error: SP and SB must be tuple of symbols."
        new{SP,SB,typeof(body)}(body)
    end
end
latexsuperscript(::Type{<:LaTeX{SP}}) where SP=SP
latexsubscript(::Type{<:LaTeX{SP,SB} where SP}) where SB=SB

@generated function oidcoord(vector::SVector{N,Float}) where N
    exprs=[:(vector[$i]===-0.0 ? 0.0 : vector[$i]) for i=1:N]
    return :(SVector($(exprs...)))
end

"""
    OID(index::Index,::Nothing,::Nothing,seq::Union{Nothing,Int})
    OID(index::Index,rcoord::SVector{N,Float},icoord::SVector{N,Float},seq::Union{Nothing,Int}) where N
    OID(index::Index,rcoord::Vector{Float},icoord::Vector{Float},seq::Union{Nothing,Int})
    OID(index::Index;rcoord::Union{Nothing,SVector,Vector{Float}}=nothing,icoord::Union{Nothing,SVector,Vector{Float}}=nothing,seq::Union{Nothing,Int}=nothing)

Operator id.
"""
struct OID{I<:Index,RC<:Union{Nothing,SVector},IC<:Union{Nothing,SVector},S<:Union{Nothing,Int}} <: SimpleID
    index::I
    rcoord::RC
    icoord::IC
    seq::S
    OID(index::Index,::Nothing,::Nothing,seq::Union{Nothing,Int})=new{typeof(index),Nothing,Nothing,typeof(seq)}(index,nothing,nothing,seq)
    function OID(index::Index,rcoord::SVector{N,Float},icoord::SVector{N,Float},seq::Union{Nothing,Int}) where N
        new{typeof(index),SVector{N,Float},SVector{N,Float},typeof(seq)}(index,oidcoord(rcoord),oidcoord(icoord),seq)
    end
end
OID(index::Index,rcoord::Vector{Float},icoord::Vector{Float},seq::Union{Nothing,Int})=OID(index,SVector{length(rcoord)}(rcoord),SVector{length(icoord)}(icoord),seq)
function OID(index::Index;rcoord::Union{Nothing,SVector,Vector{Float}}=nothing,icoord::Union{Nothing,SVector,Vector{Float}}=nothing,seq::Union{Nothing,Int}=nothing)
    OID(index,rcoord,icoord,seq)
end
totuple(::Nothing)=nothing
@generated totuple(v::SVector{N}) where N=Expr(:tuple,[:(v[$i]) for i=1:N]...)
Base.hash(oid::OID{<:Index},h::UInt)=hash((oid.index,totuple(oid.rcoord)),h)
Base.fieldnames(::Type{<:OID})=(:index,:rcoord,:icoord,:seq)
idpropertynames(::Type{<:ID{OID}})=(:indexes,:rcoords,:icoords,:seqs)

"""
    show(io::IO,oid::OID)

Show an operator id.
"""
function Base.show(io::IO,oid::OID)
    @printf io "OID(%s" oid.index
    oid.rcoord===nothing ? (@printf io ",:") : (@printf io ",[%s]" join(oid.rcoord,","))
    oid.icoord===nothing ? (@printf io ",:") : (@printf io ",[%s]" join(oid.icoord,","))
    @printf io ",%s)" oid.seq===nothing ? ":" : oid.seq
end

"""
    repr(oid::OID,l::LaTeX) -> String

LaTeX string representation of an oid.
"""
Base.repr(oid::OID,l::LaTeX)=@sprintf "%s^{%s}_{%s}" script(oid,l,Val(:B)) join(script(oid,l,Val(:SP)),',') join(script(oid,l,Val(:SB)),',')

"""
    adjoint(oid::OID) -> typeof(oid)
    adjoint(oid::ID{OID,N}) where N -> typeof(oid)

Get the adjoint of an operator id.
"""
Base.adjoint(oid::OID)=OID(oid.index',oid.rcoord,oid.icoord,oid.seq)
@generated Base.adjoint(oid::ID{OID,N}) where N=Expr(:call,:ID,[:(oid[$i]') for i=N:-1:1]...)

"""
    isHermitian(oid::ID{OID,N}) where N -> Bool

Judge whether an operator id is Hermitian.
"""
function isHermitian(oid::ID{OID,N}) where N
    for i=1:((N+1)÷2)
        oid[i]'==oid[N+1-i] || return false
    end
    return true
end

"""
    script(oid::OID,::Val{attr}) where attr -> String

Get the `:rcoord/:icoord` script of an oid.
"""
@generated function script(oid::OID,::Val{attr}) where attr
    @assert attr in (:rcoord,:icoord) "script error: not supported attr($attr)."
    fieldtype(oid,attr)<:Nothing ? 'N' : (attr=QuoteNode(attr);:("[$(join(getfield(oid,$attr),','))]"))
end

"""
    script(oid::OID,l::LaTeX,::Val{:B}) -> Any
    script(oid::OID,l::LaTeX,::Val{:SP}) -> Tuple
    script(oid::OID,l::LaTeX,::Val{:SB}) -> Tuple

Get the body/superscript/subscript of the latex string representation of an oid.
"""
@generated script(oid::OID,l::LaTeX,::Val{:B})=fieldtype(l,:body)<:Nothing ? :(oid.index.scope) : :(l.body)
@generated script(oid::OID,l::LaTeX,::Val{:SP})=Expr(:tuple,[:(script(oid,Val($super))) for super in QuoteNode.(l|>latexsuperscript)]...)
@generated script(oid::OID,l::LaTeX,::Val{:SB})=Expr(:tuple,[:(script(oid,Val($sub))) for sub in QuoteNode.(l|>latexsubscript)]...)

"""
    coordpresent

Indicate that the `:icoord` and `:rcoord` attributes in an oid should not be nothing.
"""
const coordpresent=Val(true)

"""
    coordabsent

Indicate that the `:icoord` and `:rcoord` attributes in an oid should be nothing
"""
const coordabsent=Val(false)

"""
    oidtype(I::Type{<:IID},B::Type{<:AbstractBond},::Type{Nothing},::Val{true})
    oidtype(I::Type{<:IID},B::Type{<:AbstractBond},::Type{<:Table},::Val{true})
    oidtype(I::Type{<:IID},B::Type{<:AbstractBond},::Type{Nothing},::Val{false})
    oidtype(I::Type{<:IID},B::Type{<:AbstractBond},::Type{<:Table},::Val{false})

Get the compatible oid type.
"""
oidtype(I::Type{<:IID},B::Type{<:AbstractBond},::Type{Nothing},::Val{true})=OID{union(B|>pidtype,I),SVector{B|>dimension,Float},SVector{B|>dimension,Float},Nothing}
oidtype(I::Type{<:IID},B::Type{<:AbstractBond},::Type{<:Table},::Val{true})=OID{union(B|>pidtype,I),SVector{B|>dimension,Float},SVector{B|>dimension,Float},Int}
oidtype(I::Type{<:IID},B::Type{<:AbstractBond},::Type{Nothing},::Val{false})=OID{union(B|>pidtype,I),Nothing,Nothing,Nothing}
oidtype(I::Type{<:IID},B::Type{<:AbstractBond},::Type{<:Table},::Val{false})=OID{union(B|>pidtype,I),Nothing,Nothing,Int}

"""
    angle(id::ID{OID},vectors::AbstractVector{<:AbstractVector{Float}},values::AbstractVector{Float}) -> Float

Get the total twist phase of an id.
"""
@generated function Base.angle(id::ID{OID},vectors::AbstractVector{<:AbstractVector{Float}},values::AbstractVector{Float})
    Expr(:call,:+,[:(angle(id[$i],vectors,values)) for i=1:rank(id)]...)
end

"""
    Operator{V<:Number,I<:ID{OID}} <: Element{V,I}

Abstract type for an operator.
"""
abstract type Operator{V<:Number,I<:ID{OID}} <: Element{V,I} end
function (O::Type{<:Operator})( value,
                                indexes::NTuple{N,Index};
                                rcoords::Union{Nothing,NTuple{N,SVector{M,Float}}}=nothing,
                                icoords::Union{Nothing,NTuple{N,SVector{M,Float}}}=nothing,
                                seqs::Union{Nothing,NTuple{N,Int}}=nothing
                                ) where {N,M}
    rcoords===nothing && (rcoords=ntuple(i->nothing,N))
    icoords===nothing && (icoords=ntuple(i->nothing,N))
    seqs===nothing && (seqs=ntuple(i->nothing,N))
    return O(value,ID(OID,indexes,rcoords,icoords,seqs))
end

const latexformats=Dict{Symbol,LaTeX}()
"""
    latexformat(T::Type{<:Operator}) -> LaTeX
    latexformat(T::Type{<:Operator},l::LaTeX) -> LaTeX

Get/Set the latex format for a subtype of `Operator`.
"""
latexformat(T::Type{<:Operator})=latexformats[nameof(T)]
latexformat(T::Type{<:Operator},l::LaTeX)=latexformats[nameof(T)]=l

"""
    show(io::IO,opt::Operator)
    show(io::IO,::MIME"text/latex",opt::Operator)

Show an operator.
"""
Base.show(io::IO,opt::Operator)=@printf io "%s(%s,%s)" nameof(typeof(opt)) decimaltostr(opt.value) opt.id
Base.show(io::IO,::MIME"text/latex",opt::Operator)=show(io,MIME"text/latex"(),latexstring(repr(opt)))

"""
    adjoint(opt::Operator{N}) where N -> Operator

Get the adjoint of an operator.
"""
Base.adjoint(opt::Operator)=rawtype(typeof(opt))(opt.value',opt.id')

"""
    repr(opt::Operator,l::Union{LaTeX,Nothing}=nothing) -> String

Get the latex string representation of an operator.
"""
Base.repr(opt::Operator,::Nothing=nothing)=repr(opt,latexformat(typeof(opt)))
function Base.repr(opt::Operator,l::LaTeX)
    rank(opt)==0 && return replace(valuetolatextext(opt.value)," "=>"")
    return @sprintf "%s%s" valuetostr(opt.value) join(NTuple{rank(opt),String}(repr(opt.id[i],l) for i=1:rank(opt)),"")
end
function valuetostr(v)
    v==1 && return ""
    v==-1 && return "-"
    result=valuetolatextext(v)
    if occursin(" ",result) || (isa(v,Complex) && real(v)≠0 && imag(v)≠0)
        bra,ket=occursin("(",result) ? ("[","]") : ("(",")")
        result=@sprintf "%s%s%s" bra replace(result," "=>"") ket
    end
    return result
end
valuetolatextext(value::Union{Real,Complex})=decimaltostr(value)
function valuetolatextext(value)
    io=IOBuffer()
    showable(MIME"text/latex"(),value) ? show(IOContext(io,:limit=>false),MIME"text/latex"(),value) : show(IOContext(io,:limit=>false),MIME"text/plain"(),value)
    return replace(replace(String(take!(io)),"\\begin{equation*}"=>""),"\\end{equation*}"=>"")
end

"""
    isHermitian(opt::Operator) -> Bool

Judge whether an operator is Hermitian.
"""
isHermitian(opt::Operator)=isa(opt.value,Real) && isHermitian(opt.id)

"""
    rcoord(opt::Operator) -> SVector

Get the whole rcoord of an operator.
"""
@generated function rcoord(opt::Operator)
    rank(opt)==1 && return :(opt.id[1].rcoord)
    rank(opt)==2 && return :(opt.id[1].rcoord-opt.id[2].rcoord)
    error("rcoord error: not supported rank($(rank(opt))) of $(nameof(opt)).")
end

"""
    icoord(opt::Operator) -> SVector

Get the whole icoord of an operator.
"""
@generated function icoord(opt::Operator)
    rank(opt)==1 && return :(opt.id[1].icoord)
    rank(opt)==2 && return :(opt.id[1].icoord-opt.id[2].icoord)
    error("icoord error: not supported rank($(rank(opt))) of $(nameof(opt)).")
end

"""
    sequence(opt::Operator,table=nothing) -> NTuple{rank(opt),Int}

Get the sequence of the oids of an operator according to a table.
"""
@generated function sequence(opt::Operator,table=nothing)
    table<:Nothing && return :(opt.id.seqs::NTuple{rank(opt),Int})
    return Expr(:tuple,[:(get(table,opt.id[$i].index,nothing)) for i=1:rank(opt)]...)
end

"""
    otype

Get the compatible operator type from a term type, a bond type and a table type.
"""
function otype end

"""
    Operators(opts::Operator...)

A set of operators.

Type alias for `Elements{<:ID{OID},<:Operator}`.
"""
const Operators{I<:ID{OID},O<:Operator}=Elements{I,O}
Operators(opts::Operator...)=Elements(opts...)

"""
    adjoint(opts::Operators) -> Operators

Get the adjoint of a set of operators.
"""
function Base.adjoint(opts::Operators)
    result=Operators{opts|>keytype,opts|>valtype}()
    for opt in values(opts)
        nopt=opt|>adjoint
        result[nopt.id]=nopt
    end
    return result
end

"""
    repr(opts::Operators,l::Union{LaTeX,Nothing}=nothing) -> String

Get the latex string representation of a set of operators.
"""
function Base.repr(opts::Operators,l::Union{LaTeX,Nothing}=nothing)
    result=String[]
    for (i,opt) in enumerate(values(opts))
        rep=repr(opt,l)
        i>1 && rep[1]!='-' && push!(result,"+")
        push!(result,rep)
    end
    return join(result,"")
end

"""
    show(io::IO,::MIME"text/latex",opts::Operators)

Show latex formed operators.
"""
Base.show(io::IO,::MIME"text/latex",opts::Operators)=show(io,MIME"text/latex"(),latexstring(repr(opts)))

"""
    summary(io::IO,opts::Operators)

Print a brief description of a set of operators to an io.
"""
Base.summary(io::IO,opts::Operators)=@printf io "Operators{%s}" valtype(opts)

"""
    isHermitian(opts::Operators) -> Bool

Judge whether a set of operators as a whole is Hermitian.
"""
isHermitian(opts::Operators)=opts==opts'

"""
    twist(operator::Operator,vectors::AbstractVector{<:AbstractVector{Float}},values::AbstractVector{Float}) -> Operator

Twist an operator.
"""
function twist(operator::Operator,vectors::AbstractVector{<:AbstractVector{Float}},values::AbstractVector{Float})
    replace(operator,value=operator.value*exp(1im*angle(operator.id,vectors,values)))
end

"""
    Boundary{Names}(values::AbstractVector{Float},vectors::AbstractVector{<:AbstractVector{Float}}) where Names
    Boundary()

Boundary twist of operators.
"""
struct Boundary{Names,V<:AbstractVector{Float}}
    values::Vector{Float}
    vectors::Vector{V}
    function Boundary{Names}(values::AbstractVector{Float},vectors::AbstractVector{<:AbstractVector{Float}}) where Names
        @assert length(Names)==length(values)==length(vectors) "Boundary error: dismatched names, values and vectors."
        new{Names,eltype(vectors)}(convert(Vector{Float},values),vectors)
    end
end
const boundaryemptyvalues=Float[]
const boundaryemptyvectors=SVector{0,Float}[]
Boundary()=Boundary{()}(boundaryemptyvalues,boundaryemptyvectors)

"""
    ==(bound1::Boundary,bound2::Boundary) -> Bool
    isequal(bound1::Boundary,bound2::Boundary) -> Bool
"""
Base.:(==)(bound1::Boundary,bound2::Boundary) = ==(efficientoperations,bound1,bound2)
Base.isequal(bound1::Boundary,bound2::Boundary)=isequal(efficientoperations,bound1,bound2)

"""
    (bound::Boundary)(operator::Operator) -> Operator
    (bound::Boundary{()})(operator::Operator) -> Operator

Get the boundary twisted operator.
"""
(bound::Boundary)(operator::Operator)=twist(operator,bound.vectors,bound.values)
(bound::Boundary{()})(operator::Operator)=operator

"""
    angle(bound::Boundary,operator::Operator) -> Float
    angle(bound::Boundary{()},operator::Operator) -> Int

Get the boundary twist phase of an operator.
"""
Base.angle(bound::Boundary,operator::Operator)=angle(operator.id,bound.vectors,bound.values)
Base.angle(bound::Boundary{()},operator::Operator)=0

"""
    update!(bound::Boundary{Names},args...;kwargs...) where Names -> Boundary

Update the values of the boundary twisted phase.
"""
@generated function update!(bound::Boundary{Names},args...;kwargs...) where Names
    return Expr(:block,[:(bound.values[$i]=get(kwargs,$name,bound.values[$i])) for (i,name) in enumerate(QuoteNode.(Names))]...,:(bound))
end

end #module


text:
# This file is a part of BAT.jl, licensed under the MIT License (MIT).

struct HistogramAsUvDistribution{T <: AbstractFloat} <: ContinuousUnivariateDistribution
    h::Histogram{<:Real, 1}
    inv_weights::Vector{T}
    edges::Vector{T}
    volumes::Vector{T}

    probabilty_edges::Vector{T}
    probabilty_volumes::Vector{T}
    probabilty_inv_volumes::Vector{T}

    acc_prob::Vector{T}

    μ::T
    var::T
    cov::Matrix{T}
    σ::T
end

function HistogramAsUvDistribution(h::Histogram{<:Real, 1}, T::DataType = Float64)
    nh = normalize(h)
    probabilty_widths::Vector{T} = h.weights * inv(sum(h.weights))
    probabilty_edges::Vector{T} = Vector{Float64}(undef, length(probabilty_widths) + 1)
    probabilty_edges[1] = 0
    @inbounds for (i, w) in enumerate(probabilty_widths)
        probabilty_edges[i+1] = probabilty_edges[i] + probabilty_widths[i]
    end
    probabilty_edges[end] = 1
    volumes = diff(h.edges[1])
    mean = Statistics.mean(StatsBase.midpoints(nh.edges[1]), ProbabilityWeights(nh.weights))
    var = Statistics.var(StatsBase.midpoints(nh.edges[1]), ProbabilityWeights(nh.weights), mean = mean)

    acc_prob::Vector{T} = zeros(T, length(nh.weights))
    for i in 2:length(acc_prob)
        acc_prob[i] += acc_prob[i-1] + nh.weights[i-1] * volumes[i-1]
    end

    d::HistogramAsUvDistribution{T} = HistogramAsUvDistribution{T}(
        nh,
        inv.(nh.weights),
        nh.edges[1],
        volumes,
        probabilty_edges,
        probabilty_widths,
        inv.(probabilty_widths),
        acc_prob,
        mean,
        var,
        fill(var, 1, 1),
        sqrt(var)
    )
end

function Base.rand(rng::AbstractRNG, d::HistogramAsUvDistribution{T})::T where {T <: AbstractFloat}
    r::T = rand()
    next_inds::UnitRange{Int} = searchsorted(d.probabilty_edges, r)
    next_ind_l::Int = next_inds.start
    next_ind_r::Int = next_inds.stop
    if next_ind_l > next_ind_r
        next_ind_l = next_inds.stop
        next_ind_r = next_inds.start
    end
    ret::T = d.edges[next_ind_l]
    if next_ind_l < next_ind_r
        ret += d.volumes[next_ind_l] * (d.probabilty_edges[next_ind_r] - r) * d.probabilty_inv_volumes[next_ind_l]
    end
    return ret
end

function Distributions.pdf(d::HistogramAsUvDistribution{T}, x::Real)::T where {T <: AbstractFloat}
    i::Int = StatsBase.binindex(d.h, x)
    return @inbounds d.h.weights[i]
end

function Distributions.logpdf(d::HistogramAsUvDistribution{T}, x::Real)::T where {T <: AbstractFloat}
    return log(pdf(d, x))
end

function Distributions.cdf(d::HistogramAsUvDistribution{T}, x::Real)::T where {T <: AbstractFloat}
    i::Int = StatsBase.binindex(d.h, x)
    p::T = @inbounds sum(d.h.weights[1:i-1] .* d.volumes[1:i-1])
    p += (x - d.edges[i]) * d.h.weights[i]
    return p
end

function Distributions.minimum(d::HistogramAsUvDistribution{T})::T where {T <: AbstractFloat}
    d.edges[1]
end

function Distributions.maximum(d::HistogramAsUvDistribution{T})::T where {T <: AbstractFloat}
    d.edges[end]
end

function Distributions.insupport(d::HistogramAsUvDistribution{T}, x::Real)::Bool where {T <: AbstractFloat}
    d.edges[1] <= x <= d.edges[end]
end

function Distributions.quantile(d::HistogramAsUvDistribution{T}, x::Real)::T where {T <: AbstractFloat}
    r::UnitRange{Int} = searchsorted(d.acc_prob, x)
    idx::Int = min(r.start, r.stop)
    p::T = d.acc_prob[ idx ]
    q::T = d.edges[idx]
    missing_p::T = x - p
    inv_weight::T = d.inv_weights[idx]
    if !isinf(inv_weight)
        q += missing_p * inv_weight
    end
    return min(q, maximum(d))
end

Base.eltype(d::HistogramAsUvDistribution{T}) where {T <: AbstractFloat}= T

_np_bounds(d::HistogramAsUvDistribution) =
    HyperRectBounds(Vector{eltype(d)}([quantile(d, 0)]), Vector{eltype(d)}([quantile(d, 1)]), hard_bounds)

Statistics.mean(d::HistogramAsUvDistribution) = d.μ
Statistics.var(d::HistogramAsUvDistribution) = d.var
Statistics.cov(d::HistogramAsUvDistribution) = d.cov


text:
using Documenter
using TensorKit

makedocs(modules=[TensorKit],
            sitename="TensorKit.jl",
            authors = "Jutho Haegeman",
            format = Documenter.HTML(; prettyurls = get(ENV, "CI", nothing) == "true"),
            pages = [
                "Home" => "index.md",
                "Manual" => ["man/intro.md", "man/spaces.md", "man/sectors.md", "man/tensors.md"],
                "Library" => ["lib/spaces.md"]
            ])

deploydocs(repo = "github.com/Jutho/TensorKit.jl.git")