text:
# This file holds the definition of the functions pertaining to the Attribute
#   type.

export  setAttributeName!,
        addPossibleAttributeValue!,
        removePossibleAttributeValue!,
        clearPossibleAttributeValues!,
        setPossibleAttributeValues!,
        addInitialAttributeValue!,
        removeInitialAttributeValue!,
        clearInitialAttributeValues!,
        setInitialAttributeValues!

"""
```
setAttributeName!(
    attribute::Attribute,
    name::String )
```
This function sets the name of the personnel attribute `attribute` to `name`.

This function returns `true`, indicating the name has been set successfully.
"""
function setAttributeName!( attribute::Attribute, name::String )::Bool

    attribute.name = name
    return true

end  # setAttributeName!( attribute, name )

"""
```
addPossibleAttributeValue!(
    attribute::Attribute,
    values::String... )
```
This function adds the values in `values` to the list of possible attribute values of the personnel attribute `attribute`, assuming they aren't in the list already.

This function returns `true` if any of the values has been added successfully, `false` if all of them were already in the list.
"""
function addPossibleAttributeValue!( attribute::Attribute,
    values::String... )::Bool

    tmpVals = unique( collect( values ) )
    inList = map( val -> val ∈ attribute.possibleValues, tmpVals )

    if all( inList )
        return false
    end  # if all( inList )

    append!( attribute.possibleValues, tmpVals[.!inList] )
    return true

end  # addPossibleAttributeValue!( attribute, values )

"""
```
removePossibleAttributeValue!(
    attribute::Attribute,
    values::String... )
```
This function removes the values in `values` from the list of possible attribute values of the personnel attribute `attribute`, assuming they are in the list. This function will also remove these values from the initialisation values if necessary.

This function returns `true` if any of the values has been removed successfully, `false` if none of them were in the list to begin with.
"""
function removePossibleAttributeValue!( attribute::Attribute,
    values::String... )::Bool

    toDelete = map( val -> val ∈ values, attribute.possibleValues )

    if !any( toDelete )
        return false
    end  # if !any( toDelete )

    deleteat!( attribute.possibleValues, toDelete )
    return true

end  # removePossibleAttributeValue!( attribute, values )

"""
```
clearPossibleAttributeValues!( attribute::Attribute )
```
This function clears the list of possible values of the personnel attribute `attribute`. This function also clears the initialisation values.

This function returns `true`, indicating the list of possible values has been cleared successfully.
"""
function clearPossibleAttributeValues!( attribute::Attribute )::Bool

    empty!( attribute.possibleValues )
    return true

end  # clearPossibleAttributeValues!( attribute::Attribute )

"""
```
setPossibleAttributeValues!(
    attribute::Attribute,
    values::Vector{String} )
```
This function sets the list of possible values of the personnel attribute `attribute` to the list given in `values`. This function will also remove those entries among the initialisation values that aren't found in the new list of possible values.

This function returns `true`, indicating the possible values have been set successfully.
"""
function setPossibleAttributeValues!( attribute::Attribute,
    values::Vector{String} )::Bool

    attribute.possibleValues = unique( values )
    return true

end  # setPossibleAttributeValues!( attribute, values )

"""
```
addInitialAttributeValue!(
    attribute::Attribute,
    valWeights::Tuple{String,T}... ) where T <: Real
```
This function adds each value/weight pair in `valWeights` to the initial values of the personnel attribute `attribute`. Values that aren't in the list of possible values get added as well.

Some remarks:
* If multiple pairs have the same value, this function gives a warning and makes no changes.
* Values with negative weights get ignored.
* If a value is already in the list of initial values, its weight gets updated.

This function returns `true` if it successfully adds initial values and/or updates their weights, and `false` if it makes no changes.
"""
function addInitialAttributeValue!( attribute::Attribute,
    valWeights::Tuple{String,T}... ) where T <: Real

    newVals = map( valWeight -> valWeight[1], collect( valWeights ) )
    newWeights = map( valWeight -> valWeight[2], collect( valWeights ) )
    addPossibleAttributeValue!( attribute, newVals... )

    if length( newVals ) != length( unique( newVals ) )
        @warn "Duplicate entries in the value/weight list, not making any changes."
        return false
    end  # if length( newVals ) != length( unique( newVals ) )

    # Filter out values with negative weights.
    isPosWeight = newWeights .> 0
    newVals = newVals[isPosWeight]
    newWeights = newWeights[isPosWeight]

    if isempty( newVals )
        return false
    end  # if isempty( newVals )

    # Add the new values to the list of initial values.
    isNewVal = map( newVal -> newVal ∉ attribute.initValues, newVals )
    append!( attribute.initValues, newVals[isNewVal] )
    append!( attribute.initValueWeights, zeros( Float64, count( isNewVal ) ) )

    # Add/update the weights where needed.
    valInds = findfirst.( map( newVal -> newVal .== attribute.initValues,
        newVals ) )
    attribute.initValueWeights[valInds] = newWeights

    createInitialValueDistribution( attribute )
    return true

end  # addInitialAttributeValue!( attribute, valWeights )

"""
```
addInitialAttributeValue!(
    attribute::Attribute,
    val::String,
    weight::Real )
```
This function adds the value `val` as an initial value with weight `weight` to the personnel attribute `attribute`.

All remarks for the version of this function with `Tuples` are valid here as well.

This function returns `true` if it successfully adds the initial value or updates its weight, and `false` if it makes no changes.
"""
addInitialAttributeValue!( attribute::Attribute, val::String,
    weight::Real )::Bool = addInitialAttributeValue!( attribute,
    (val, Float64( weight )) )

"""
```
removeInitialAttributeValue!(
    attribute::Attribute,
    values::String... )
```
This function removes the values in `values` from the list of initial values of the personnel attribute `attribute`, insofar as those values are on the list.

This function returns `true` if any values have been successfully removed from the list of initial values, and `false` if none were removed.
"""
function removeInitialAttributeValue!( attribute::Attribute,
    values::String... )::Bool

    removeVals = map( val -> val ∈ values, attribute.initValues )

    if !any( removeVals )
        return false
    end  # if !any( removeVals )

    deleteat!( attribute.initValues, removeVals )
    deleteat!( attribute.initValueWeights, removeVals )
    createInitialValueDistribution( attribute )
    return true

end  # removeInitialAttributeValue!( attribute, values )

"""
```
clearInitialAttributeValues!( attribute::Attribute )
```
This function clears the list of initial values of the personnel attribute `attribute`.

This function returns `true`, indicating the list of initial values has been cleared successfully.
"""
function clearInitialAttributeValues!( attribute::Attribute )::Bool

    empty!( attribute.initValues )
    empty!( attribute.initValueWeights )
    return true

end  # clearInitialAttributeValues!( attribute::Attribute )

"""
```
setInitialAttributeValues!(
    attribute::Attribute,
    valWeights::Dict{String,T} )
    where T <: Real
```
This function sets the list of interval values, with associated weights, of the personnel attribute `attribute` to the pairs in `valWeights`. Values with a negative weight are ignored. If the dictionary contains no valid entires, this function makes no changes. Hence, to clear the list of initial values, `clearInitialAttributeValues!` should be used.

This function returns `true` if the list of initial values is successfully set, and `false` if the dictionary contained no valid entries.
"""
function setInitialAttributeValues!( attribute::Attribute,
    valWeights::Dict{String,T} )::Bool where T <: Real

    newVals = filter( val -> valWeights[val] > 0,
        collect( keys( valWeights ) ) )

    if isempty( newVals )
        return false
    end  # if isempty( newVals )

    addPossibleAttributeValue!( attribute, newVals... )
    attribute.initValues = newVals
    attribute.initValueWeights = map( val -> Float64( valWeights[val] ),
        newVals )
    createInitialValueDistribution( attribute )
    return true

end  # setInitialAttributeValues!( attribute, valWeights )

"""
```
setInitialAttributeValues!(
    attribute::Attribute,
    valWeights::Tuple{String,Float64}... )
```
This function sets the list of interval values, with associated weights, of the personnel attribute `attribute` to the pairs in `valWeights`.

If the list of initial values contains duplicate entries, the function will issue a warning and not make any changes. All remarks associated with the version of this function with a `Dict` are valid for this version.

This function returns `true` if the list of initial values is successfully set, and `false` if the list of initial values contained no valid entries.
"""
function setInitialAttributeValues!( attribute::Attribute,
    valWeights::Tuple{String,Float64}... )::Bool

    # Check for duplicates.
    newVals = map( valWeight -> valWeight[1], collect( valWeights ) )

    if length( newVals ) != length( unique( newVals ) )
        @warn "Duplicate entries in the value/weight list, not making any changes."
        return false
    end  # if length( newVals ) != length( unique( newVals ) )

    # Convert the val/weight pairs to a dictionary.
    valWeightDict = Dict{String,Float64}()

    for valWeight in valWeights
        valWeightDict[valWeight[1]] = valWeight[2]
    end  # for valWeight in valWeights

    return setInitialAttributeValues!( attribute, valWeightDict )

end  # setInitialAttributeValues!( attribute, valWeights )

"""
```
setInitialAttributeValues!(
    attribute::Attribute,
    values::Vector{String},
    weights::Vector{T} )
    where T <: Real
```
This function sets the list of interval values, with associated weights, of the personnel attribute `attribute` to the list in `values`. The weights associated with each value are given through the vector `weights`.

If the vector of initial values and weights are of unequal length, the function will issue a warning and not make any changes.

If the list of initial values contains duplicate entries, the function will issue a warning and not make any changes.

All remarks associated with the version of this function with a `Dict` are valid for this version.

This function returns `true` if the list of initial values is successfully set, and `false` if the list of initial values contained no valid entries.
"""
function setInitialAttributeValues!( attribute::Attribute,
    values::Vector{String}, weights::Vector{T} )::Bool where T <: Real

    if length( values ) != length( weights )
        @warn "Mismatched lengths of vector of initial values and weights, not making any changes."
        return false
    end  # if length( values ) != length( weights )

    if length( values ) != length( unique( values ) )
        @warn "Duplicate entries in the initial value list, not making any changes."
        return false
    end  # if length( values ) != length( unique( values ) )

    # Convert the val/weight pairs to a dictionary.
    valWeightDict = Dict{String,Float64}()

    for ii in eachindex( values )
        valWeightDict[values[ii]] = weights[ii]
    end  # for ii in eachindex( values )

    return setInitialAttributeValues!( attribute, valWeightDict )

end  # setInitialAttributeValues!( attribute, values, weights )

function Base.show( io::IO, attribute::Attribute )::Nothing

    print( io, "Attribute: ", attribute.name )

    if isempty( attribute.possibleValues )
        print( io, "\n  Attribute can't take any values." )
        return
    end  # if isempty( attribute.possibleValues )

    print( io, "\n  Possible values: ", join( attribute.possibleValues, ", " ) )

    if isempty( attribute.initValues )
        print( io, "\n  No possible initial values" )
        return
    end  # if isempty( attribute.initValues )

    print( io, "\n  Initial values:" )
    initValueProbs = 100.0 * attribute.initValueWeights /
        sum( attribute.initValueWeights )

    for ii in eachindex( attribute.initValues )
        print( io, "\n    ", attribute.initValues[ii], " (weight ",
            attribute.initValueWeights[ii], " / ",
            round( initValueProbs[ii], sigdigits = 4 ), "%)" )
    end  # for ii in eachindex( attribute.initValues )

    return

end  # show( io, attribute )

include( joinpath( privPath, "attribute.jl" ) )


text:
const IMGUI_DEMO_INPUT_TEXT_STATES = Dict(
    "basic_input_text" => InputTexts.State(
        "input text",
        size=280,
        buffer="Hello, world!" * "\0"^267,
    ),
)

const IMGUI_DEMO_INPUT_TEXT_WITH_HINT_STATES = Dict(
    "basic_input_text_with_hint" => InputTextWithHints.State(
        "input text (w/ hint)",
        "enter text here",
    ),
)

const IMGUI_DEMO_INPUT_INT_STATES = Dict(
    "basic_input_int" => InputInts.State(
        "input int",
        123,
    ),
)

const IMGUI_DEMO_INPUT_FLOAT_STATES = Dict(
    "basic_input_float" => InputFloats.State(
        "input float",
        0.001f0,
        step=0.01f0,
        speed=1.0f0,
        format="%.3f",
    ),
    "basic_input_scientific" => InputFloats.State(
        "input scientific",
        Cfloat(1.e10),
        step=0.0f0,
        speed=0.0f0,
        format="%e",
    ),
)

const IMGUI_DEMO_INPUT_DOUBLE_STATES = Dict(
    "basic_input_double" => InputDoubles.State(
        "input double",
        999999.00000001,
        step=0.01f0,
        speed=1.0f0,
        format="%.8f",
    ),
)

ext: jl


text:
struct PhyloFields
    phylo::Phylo.Phylogeny
    nodespecies::Matrix{Bool}

    function PhyloFields(phylo, nodespecies)
        (Ntip(phylo) == size(nodespecies, 1) && Nnode(phylo) == size(nodespecies, 2)) || throw(DimensionMismatch("Dimension mismatch between nodespecies matrix and phylogeny"))
        new(phylo, nodespecies)
    end
end

struct PhyloAssemblage{T} <: SEAssemblage # A type to keep subtypes together, ensuring that they are all aligned at all times
    site::SELocations
    occ::SpeciesData{T}
    phy::PhyloFields

    function PhyloAssemblage(site, occ, phy)
        size(occ.commatrix.occurrences, 2) == size(coordinates(site), 1) || throw(DimensionMismatch("Length mismatch between occurrence matrix and coordinates"))
        Ntip(phy.Phylo) == size(occ.commatrix.occurrences, 1) || throw(DimensionMismatch("Occurrence matrix and phylogeny do not match in species numbers"))
        new(site, occ, phy)
    end
end

function PhyloAssemblage(site::SELocations, occ::SpeciesData, phylo::Phylogeny)
  nodespec = createNodeBySpeciesMatrix(phylo)
  PhyloAssemblage(site, occ, PhyloFields(phylo, nodespec))
end

PhyloAssemblage(assm::Assemblage, phylo::Phylogeny) = PhyloAssemblage(assm.site, assm.occ, phylo)

function createNodeBySpeciesMatrix(tree::Phylogenetics.Phylogeny)
   colname = tree.tipLabel
   rowname = ["$num" for num in nodeNumbers(tree)]
   nodespecies = zeros(Int, nNode(tree), nTip(tree))

   function loc(tree::Phylogenetics.Phylogeny, node::Int)
     if node <= nTip(tree)
       return node
     end
     ret = [loc(tree, descendants(node, tree)[1])..., loc(tree, descendants(node, tree)[2])...]
     nodespecies[node-nTip(tree), ret] = 1
     ret
   end
   loc(tree, basalNode(tree))

   nodespecies
end


text:
using StochasticSemiDiscretizationMethod
using Test

function createHayesProblem(a,β)
    AMx =  ProportionalMX(a*ones(1,1));
    τ1=1.
    BMx1 = DelayMX(τ1,zeros(1,1));
    cVec = Additive(1)
    noiseID = 1
    αMx1 = stCoeffMX(noiseID,ProportionalMX(zeros(1,1)))
    βMx11 = stCoeffMX(noiseID,DelayMX(τ1,β*ones(1,1)))
    σ = stAdditive(1,Additive(ones(1)))
    LDDEProblem(AMx,[BMx1],[αMx1],[βMx11],cVec,[σ])
end

function createSLDOProblem(A,B,ζ,α,β,σ)
    AMx =  ProportionalMX(@SMatrix [0. 1.;-A -2ζ]);
    τ1=2π
    BMx1 = DelayMX(τ1,@SMatrix [0. 0.; B 0.]);
    cVec = Additive(2)
    noiseID = 1
    αMx1 = stCoeffMX(noiseID,ProportionalMX(@SMatrix [0. 0.; α 0.]))
    βMx11 = stCoeffMX(noiseID,DelayMX(τ1,@SMatrix [0. 0.; β 0.]))
    σVec = stAdditive(1,Additive(@SVector [0., σ]))
    LDDEProblem(AMx,[BMx1],[αMx1],[βMx11],cVec,[σVec])
end

function tests()
    @testset "Testing the SemiDiscretizationMethod package with the examples" begin
		#Hayes
        @test begin
            hayes_lddep=createHayesProblem(-6.,2.); # LDDE problem for Hayes equation
            method=SemiDiscretization(0,0.1) # 0th order semi discretization with Δt=0.1
            τmax=1. # the largest τ of the system
            # Second Moment mapping
            mapping=DiscreteMapping_M2(hayes_lddep,method,τmax,n_steps=10,calculate_additive=true); #The discrete mapping of the system

            spectralRadiusOfMapping(mapping); # spectral radius ρ of the mapping matrix (ρ>1 unstable, ρ<1 stable)
            statM2=VecToCovMx(fixPointOfMapping(mapping), length(mapping.M1_Vs[1])); # stationary second moment matrix of the hayes equation (equilibrium position)
            true
        end
		#Stochastic Delay Oscillator
        @test begin
            SLDOP_lddep=createSLDOProblem(1.,0.1,0.1,0.1,0.1,0.5); # LDDE problem for Hayes equation
            method=SemiDiscretization(5,(2π+100eps())/10) # 5th order semi discretization with Δt=2π/10
            τmax=2π # the largest τ of the system
            # Second Moment mapping
            mapping=DiscreteMapping_M2(SLDOP_lddep,method,τmax,n_steps=10,calculate_additive=true); #The discrete mapping of the system

            spectralRadiusOfMapping(mapping); # spectral radius ρ of the mapping matrix (ρ>1 unstable, ρ<1 stable)
            statM2=VecToCovMx(fixPointOfMapping(mapping), length(mapping.M1_Vs[1])); # stationary second moment matrix of the hayes equation (equilibrium position)
            true
        end
    end
end
tests()


text:
#=
Filename: aiyagari_compute_policy.jl
Author: Victoria Gregory
Date: 8/29/2016

Computes and plots the optimal policy of a household
from the Aiyagari model, given prices.
=#

using LaTeXStrings
using Plots
pyplot()

# Example prices
r = 0.03
w = 0.956

# Create an instance of Household
am = Household(a_max=20.0, r=r, w=w)

# Use the instance to build a discrete dynamic program
am_ddp = DiscreteDP(am.R, am.Q, am.beta)

# Solve using policy function iteration
results = solve(am_ddp, PFI)

# Simplify names
z_size, a_size  = am.z_size, am.a_size
z_vals, a_vals = am.z_chain.state_values, am.a_vals
n = am.n

# Get all optimal actions across the set of
# a indices with z fixed in each column
a_star = reshape([a_vals[results.sigma[s_i]] for s_i in 1:n], a_size, z_size)

labels = [latexstring("\z = $(z_vals[1])") latexstring("\z = $(z_vals[2])")]
plot(a_vals, a_star, label=labels, lw=2, alpha=0.6)
plot!(a_vals, a_vals, label="", color=:black, linestyle=:dash)
plot!(xlabel="current assets", ylabel="next period assets", grid=false)


text:
module EponymTuples

export @eponymargs, @eponymtuple

var_and_type(ex::Symbol) = ex => Any

function var_and_type(ex::Expr)
    if ex.head == :(::) && (ex.args[1] isa Symbol)
        ex.args[1] => ex.args[2]
    else
        throw(ArgumentError("Can't handle expression $(ex), use a symbol or `symbol::type`."))
    end
end

function check_duplicate_vars(itr)
    seen = Set{Symbol}()
    for symbol in itr
        if symbol ∈ seen
            error(ArgumentError("duplicate variable name $(symbol)"))
        else
            push!(seen, symbol)
        end
    end
    nothing
end

"""
    @eponymargs(a, b::T, ...)

Expands to form like `((a, b)::NamedTuple{(:a, :b), <: Tuple{Any, T}})`, using
the variable names *both* for the `NamedTuple` and deconstruction in the
function arguments.

Valid arguments are variable names, or names followed by a type specifier (when
missing, `Any` is used).

# Example

```jldoctest
julia> using EponymTuples

julia> foo(@eponymargs(a, b)) = a + b
foo (generic function with 1 method)

julia> foo((a = 1, b = 2))
3
```
"""
macro eponymargs(args...)
    vars_and_types = map(var_and_type, args)
    vars = first.(vars_and_types)
    check_duplicate_vars(vars)
    types = map(esc ∘ last, (vars_and_types))
    :(($(map(esc, vars)...),)::$(esc(:NamedTuple)){$(vars),
                                                   <: $(esc(:Tuple)){$(types...)}})
end

var_and_value_form(ex::Symbol) = (e = esc(ex); :($(e) = $(e)))

function var_and_value_form(ex::Expr)
    if ex.head == :(=) && (ex.args[1] isa Symbol)
        :($(esc(ex.args[1])) = $(esc(ex.args[2])))
    else
        throw(ArgumentError("Can't handle expression $(ex), use a symbol or `symbol = value`."))
    end
end

"""
    @eponymtuple(a, b = bval, ...)

Expands to `(a = a, b = bval, ...)`, creating a named tuple.

Each argument is either a symbol, resulting in assigning its own value, or an
assignment, passed as is.

The recommended use is the first one, allowing the creation of named tuples
without repeating the name. Overriding the value is just for convenience, if you
are using it too much you are probably better off with the standard `NamedTuple`
constructor `(a = aval, b = bval, ...)`.
"""
macro eponymtuple(args...)
    forms = map(var_and_value_form, args)
    Expr(:tuple, forms...)
end

end # module


text:
export getAlpha, getBeta, getX,ParamLevelSetModelFunc
function ParamLevelSetModelFunc(Mesh::RegularMesh,m::Vector; computeJacobian = 1, bf::Int64 = 1,
		 sigma::Function = (m,n)->(n[:] .= 1.0),
		 Xc = convert(Array{Float32,2}, getCellCenteredGrid(Mesh)),u::Vector = zeros(prod(Mesh.n)),dsu::Vector = zeros(prod(Mesh.n)),
		 Jbuilder::SpMatBuilder{Int64,Float64} = getSpMatBuilder(Int64,Float64,prod(Mesh.n),length(m),computeJacobian*10*prod(Mesh.n)),
		 iifunc::Function=identity, numParamOfRBF::Int64 = 5)
	n = Mesh.n;
	if length(u) != prod(n)
		error("preallocated u is of wrong size");
	end
	u[:] .= 0.0;
	reset!(Jbuilder);
	h = Mesh.h;
	nRBFs = div(length(m),numParamOfRBF);
	xmin = Mesh.domain[[1;3;5]];
	beta_arr = Array{Float64}(undef,nRBFs);

	## in run 1 we calculate u. in run 2 we calculate J after we know the derivative of the heaviside.
	for k=1:nRBFs
		if numParamOfRBF==5
			computeFunc!(m,k,Xc,h,n,u,xmin);
		else
			# putting this inside the function takes a huge amount of time... go figure..
			(a1,a2,a3,a4,a5,a6) = getA(k,m);
			A = [a1 a2 a3 ; a2 a4 a5 ; a3 a5 a6];
			beta = sqrt(minimum(eigvals(A)));
			computeFunc_ellipsoids!(m,k,Xc,h,n,u,xmin,beta);
			beta_arr[k] = beta;
		end
	end
	sigma(u,dsu);
	if computeJacobian == 1
		if numParamOfRBF==5
			for k=1:nRBFs
				updateJacobian!(k,m,dsu,h,n,Jbuilder,u,iifunc,Xc,xmin);
			end
		else
			for k=1:nRBFs
				updateJacobian_ellipsoids!(k,m,dsu,h,n,Jbuilder,u,iifunc,Xc,xmin,beta_arr[k])
			end
		end
	end
	return u,Jbuilder;
end

function computeFunc!(m::Vector{Float64},k::Int64,Xc,h::Vector{Float64},n::Vector{Int64},
         u::Vector{Float64},xmin::Vector{Float64})
(x1,x2,x3) = getXt(k,m);
i = round.(Int64,(getX(k,m) - xmin)./h .+ 0.5);
# ## here we take a reasonable box around the bf.
# ## The constant depends on the type of RBF that we take.
beta = getBeta(k,m);
betaSq = beta^2;
thres = 0.81/betaSq;
boxL = ceil.(Int64,0.9./(h*abs(beta)));
imax = min.(i + boxL,n);
imin = max.(i - boxL,[1;1;1]);
alpha = getAlpha(k,m);

@inbounds for l = imin[3]:imax[3]
	@inbounds iishift3 = (l-1)*n[1]*n[2];
	@inbounds for j = imin[2]:imax[2]
		@inbounds iishift2 = iishift3 + (j-1)*n[1];
		@inbounds for q = imin[1]:imax[1]
			ii = iishift2 + q;
			@inbounds y = Xc[ii,1]-x1; y*=y;
			nX = y;
			@inbounds y = Xc[ii,2]-x2; y*=y;
			nX+=y;
			@inbounds y = Xc[ii,3]-x3; y*=y;
			nX+=y;
			if (nX <= thres) # 0.77~0.875^2
				nX*=betaSq;
				argii = radiust(nX);
				argii = psi1(argii);
				argii *= alpha;
				@inbounds u[ii] += argii;
			end
		end
	end
end
end

function updateJacobian!(k::Int64,m::Vector{Float64},dsu::Vector{Float64},h::Vector{Float64},n::Vector{Int64},Jbuilder::SpMatBuilder,u::Vector{Float64},iifunc::Function,Xc,xmin::Vector{Float64})
(x1,x2,x3) = getXt(k,m);
alpha = getAlpha(k,m);
beta = getBeta(k,m);
betaSQ = beta*beta;
invBeta = (1.0/beta);
thres = 0.81/betaSQ;

alphaBetaSq = alpha*betaSQ;
offset = convert(Int64,(k-1)*5 + 1);
md = 1e-3*maximum(abs.(dsu));
boxL = ceil.(Int64,0.9./(h*abs(beta)));
i = round.(Int64,([x1;x2;x3] - xmin)./h .+ 0.5);
imax = min.(i + boxL,n);
imin = max.(i - boxL,[1;1;1]);

@inbounds for l = imin[3]:imax[3]
	@inbounds iishift3 = (l-1)*n[1]*n[2];
	@inbounds for j = imin[2]:imax[2]
		@inbounds iishift2 = iishift3 + (j-1)*n[1];
		@inbounds for q = imin[1]:imax[1]
			ii = iishift2 + q;
			temp = dsu[ii];
			if temp >= md
				@inbounds y1 = x1 - Xc[ii,1]; nX =y1*y1;
				@inbounds y2 = x2 - Xc[ii,2]; nX+=y2*y2;
				@inbounds y3 = x3 - Xc[ii,3]; nX+=y3*y3;
				if (nX <= thres) # 0.77~0.875^2
					radii = radiust(nX*betaSQ);
					psi,dpsi = dpsi1_t(radii);
					psi*=temp;
					temp*= alphaBetaSq;
					temp/= radii;
					temp*=dpsi;
					nX *= temp;
					nX *= invBeta;
					y1*=temp; y2*=temp; y3*=temp;
					ii = iifunc(ii);
					setNext!(Jbuilder,ii,offset,psi);
					setNext!(Jbuilder,ii,offset+1,nX);
					setNext!(Jbuilder,ii,offset+2,y1);
					setNext!(Jbuilder,ii,offset+3,y2);
					setNext!(Jbuilder,ii,offset+4,y3);
				end
			end
		end
	end
end
end

function getAlpha(k::Int64,m::Vector{Float64},numParamOfRBF::Int64 = 5)
	return m[numParamOfRBF*(k-1) + 1];
end

function getBeta(k::Int64,m::Vector{Float64})
	return m[5*(k-1) + 2];
end

function getX(k::Int64,m::Vector{Float64},numParamOfRBF::Int64 = 5)
	return m[numParamOfRBF*(k-1).+((numParamOfRBF-2):numParamOfRBF)];
end

function getXt(k::Int64,m::Vector{Float64},numParamOfRBF::Int64 = 5)
	return (m[numParamOfRBF*(k-1)+numParamOfRBF-2],m[numParamOfRBF*(k-1)+numParamOfRBF-1],m[numParamOfRBF*(k-1)+numParamOfRBF]);
end

function setX(k::Int64,m::Vector{Float64},x::Array)
	m[5*(k-1)+(3:5)] = x;
end

export psi
function psi(r::Union{Array{Float64,1},Array{Float32,1}},bf = 1)
	rpp = max.(1.0-r,0.0);
	if bf==1
		psi = (rpp.^4).*(4.0*r + 1.0);
	elseif bf == 2
		psi = (rpp.^6).*((11.6666666)*r.^2 + 6.0*r + 1.0); ### this is the RBF divided by 3. No idea what went wrong there...
	elseif bf == 3
		psi = (rpp.^8).*(32.0*(r.^3) + 25.0.*(r.^2) + 8.0.*r + 1.0);
	end
	return psi;
end

# function dpsi(r::Union{Array{Float64,1},Array{Float32,1}},bf = 1)
	# rpp = max.(1.0-r,0.0);
	# if bf==1
		# dpsi = -4.0*(rpp.^3).*(4.0*r + 1.0) + (rpp.^4).*4.0;
	# elseif bf == 2
		# dpsi = -6.0*(rpp.^5).*((11.6666666)*r.^2 + 6.0*r + 1.0) + (rpp.^6).*((2.0*11.6666666)*r + 6.0); ### this is the RBF divided by 3. No idea what went wrong there...
	# elseif bf == 3
		# dpsi = -8.0*(rpp.^7).*(32.0*(r.^3) + 25.*(r.^2) + 8.0.*r + 1.0) + (rpp.^8).*(96.0*(r.^2) + 50.*r + 8.0);
	# end
	# return dpsi;
# end

function psi1(r::Float64)
	psi = 1.0;
	psi -= r;
	r *= 4.0;
	r += 1.0;
	psi*=psi;
	psi*=psi;
	psi*=r;
	return psi;
end

function dpsi1(r::Float64)
	dpsi =1.0;
	dpsi -= r;
	r *= -4.0;
	r -= 1.0;
	r += dpsi;
	dpsi *= (dpsi*dpsi);
	dpsi *= 4.0;
	dpsi *= r;
	return dpsi;
end

function dpsi1_t(r::Float64)
	psi = 1.0;
	psi -= r;
	r *= 4.0;
	r += 1.0;
	dpsi = psi;
	dpsi -= r;
	dpsi*=psi;
	dpsi *= 4.0;
	psi *= psi;
	dpsi *=psi;
	psi *= psi;
	psi*=r;
	return (psi,dpsi)
end

function radiust(t::Float64,eps::Float64 = 1e-3)
return sqrt(t + eps);#
end

export splitTheta,wrapTheta
function splitTheta(theta::Array)
	n = div(length(theta),5);
	if 5*n != length(theta)
		error("Inappropriate length of theta");
	end
	alpha = theta[1:5:end];
	beta = theta[2:5:end];
	Xs = [theta[3:5:end] theta[4:5:end] theta[5:5:end]];
	return (alpha,beta,Xs)
end

function wrapTheta( alpha:: Array{Float64},	beta:: Array{Float64}, Xs:: Array{Float64})
	theta = zeros(5*length(alpha));
	theta[1:5:end] = alpha;
	theta[2:5:end] = beta;
	theta[3:5:end] = Xs[:,1];
	theta[4:5:end] = Xs[:,2];
	theta[5:5:end] = Xs[:,3];
	return theta
end


text:
"""
Encapsulates the value of an evaluation calculated in one training epoch.
If the evaluation depends on a dataset,
the dataset's name can be specified also.
"""
struct MonitoringItem
   evaluation::String
   epoch::Int
   value::Float64
   datasetname::String
end

"""
A vector for collecting `MonitoringItem`s during training.
"""
const Monitor = Vector{MonitoringItem}

const monitoraislogr = "aislogr"
const monitoraisstandarddeviation = "aisstandarddeviation"
const monitorcordiff = "cordiff"
const monitorexactloglikelihood = "exactloglikelihood"
const monitorfreeenergy = "freeenergy"
const monitorloglikelihood = "loglikelihood"
const monitorlogproblowerbound = "logproblowerbound"
const monitormeandiff = "monitormeandiff"
const monitormeandiffpervariable = "monitormeandiffpervariable"
const monitorreconstructionerror = "reconstructionerror"
const monitorsd = "sd"
const monitorweightsnorm = "weightsnorm"

"""
    correlations(datadict)
Creates and returns a dictionary with the same keys as the given `datadict`.
The values of the returned dictionary are the correlations of the samples in
the datasets given as values in the `datadict`.
"""
function correlations(datadict::DataDict)
   DataDict(map(kv -> (kv[1] => cor(kv[2])), datadict))
end

"""
    means(datadict)
Creates and returns a dictionary with the same keys as the given `datadict`.
The values of the returned dictionary are the samples' means in the `datadict`.
"""
function means(datadict::DataDict)
   DataDict(map(kv -> (kv[1] => mean(kv[2], dims = 1)), datadict))
end

"""
    monitorcordiff!(monitor, rbm, epoch, cordict)
Generates samples and records the distance of their correlation matrix to the
correlation matrices for (original) datasets contained in the `cordict`.
"""
function monitorcordiff!(monitor::Monitor, bm::AbstractBM, epoch::Int,
      cordict::DataDict;
      nparticles::Int = 3000, burnin::Int = 10,
      xgenerated::Matrix{Float64} = sampleparticles(bm, nparticles, burnin)[1])

   samplecor = cor(xgenerated)
   for (datasetname, datacor) in cordict
      push!(monitor, MonitoringItem(monitorcordiff, epoch,
            norm(samplecor - datacor), datasetname))
   end

   monitor
end

"""
    monitorexactloglikelihood!(monitor, bm, epoch, datadict)
Computes the mean exact log-likelihood in the Boltzmann Machine model `bm`
for the data sets in the DataDict `datadict` and stores this information in
the Monitor `monitor`.
"""
function monitorexactloglikelihood!(monitor::Monitor, bm::AbstractBM,
      epoch::Int, datadict::DataDict)

   logz = BMs.exactlogpartitionfunction(bm)
   for (datasetname, x) in datadict
      push!(monitor, MonitoringItem(BMs.monitorexactloglikelihood, epoch,
            exactloglikelihood(bm, x, logz), datasetname))
   end

   monitor
end

"""
    monitorfreeenergy!(monitor, rbm, epoch, datadict)
Computes the free energy for the `datadict`'s data sets in the RBM model `rbm`
and stores the information in the `monitor`.
"""
function monitorfreeenergy!(monitor::Monitor, rbm::AbstractRBM,
      epoch::Int, datadict::DataDict)

   for (datasetname, x) in datadict
      push!(monitor, MonitoringItem(BMs.monitorfreeenergy, epoch,
            freeenergy(rbm, x), datasetname))
   end

   monitor
end

"""
    monitorloglikelihood!(monitor, rbm, epoch, datadict)
Estimates the log-likelihood of the `datadict`'s data sets in the RBM model
`rbm` with AIS and stores the values, together with information about the
variance of the estimator, in the `monitor`.

If there is more than one worker available, the computation is parallelized
by default. Parallelization can be turned on or off with the optional
boolean argument `parallelized`.

For the other optional keyword arguments, see `aislogimportanceweights`.

See also: `loglikelihood`.
"""
function monitorloglikelihood!(monitor::Monitor, rbm::AbstractRBM,
      epoch::Int, datadict::DataDict;
      parallelized::Bool = nworkers() > 1,
      # optional arguments for AIS:
      ntemperatures::Int = 100,
      temperatures::Array{Float64,1} = collect(0:(1/ntemperatures):1),
      nparticles::Int = 100,
      burnin::Int = 5)

   if parallelized
      logimpweights = BMs.batchparallelized(
            n -> BMs.aislogimpweights(rbm;
                  ntemperatures = ntemperatures, temperatures = temperatures,
                  nparticles = n, burnin = burnin),
            nparticles, vcat)
   else
      logimpweights = BMs.aislogimpweights(rbm;
            ntemperatures = ntemperatures, temperatures = temperatures,
            nparticles = nparticles, burnin = burnin)
   end

   logr = logmeanexp(logimpweights)
   sd = BMs.aisstandarddeviation(logimpweights)
   logz = BMs.logpartitionfunction(rbm, logr)
   for (datasetname, x) in datadict
      push!(monitor, MonitoringItem(BMs.monitorloglikelihood, epoch,
            BMs.loglikelihood(rbm, x, logz), datasetname))
   end
   push!(monitor,
         MonitoringItem(BMs.monitoraisstandarddeviation, epoch, sd, ""),
         MonitoringItem(BMs.monitoraislogr, epoch, logr, ""))

   monitor
end

"""
    monitorlogproblowerbound!(monitor, dbm, epoch, datadict)
Estimates the lower bound of the log probability of the `datadict`'s data sets
in the DBM `dbm` with AIS and stores the values, together with information about
the variance of the estimator, in the `monitor`.

If there is more than one worker available, the computation is parallelized
by default. Parallelization can be turned on or off with the optional
boolean argument `parallelized`.

For the other optional keyword arguments, see `aislogimpweights`.

See also: `logproblowerbound`.
"""
function monitorlogproblowerbound!(monitor::Monitor, dbm::MultimodalDBM,
      epoch::Int, datadict::DataDict;
      parallelized::Bool = nworkers() > 1,
      # optional arguments for AIS:
      ntemperatures::Int = 100,
      temperatures::Array{Float64,1} = collect(0:(1/ntemperatures):1),
      nparticles::Int = 100,
      burnin::Int = 5)

   if parallelized
      logimpweights = batchparallelized(
            n -> aislogimpweights(dbm;
                  ntemperatures = ntemperatures, temperatures = temperatures,
                  nparticles = n, burnin = burnin),
            nparticles, vcat)
   else
      logimpweights = aislogimpweights(dbm;
            ntemperatures = ntemperatures, temperatures = temperatures,
            nparticles = nparticles, burnin = burnin)
   end

   logr = logmeanexp(logimpweights)
   sd = aisstandarddeviation(logimpweights)
   logz = logpartitionfunction(dbm, logr)

   for (datasetname, x) in datadict
      push!(monitor, MonitoringItem(monitorlogproblowerbound, epoch,
            logproblowerbound(dbm, x, logz),
            datasetname))
   end
   push!(monitor,
         MonitoringItem(monitoraisstandarddeviation, epoch, sd, ""),
         MonitoringItem(monitoraislogr, epoch, logr, ""))

   monitor
end

function monitormeandiff!(monitor::Monitor, bm::AbstractBM, epoch::Int,
      meandict::DataDict;
      nparticles::Int = 3000, burnin::Int = 10,
      xgenerated::Matrix{Float64} = sampleparticles(bm, nparticles, burnin)[1])

   samplemean = mean(xgenerated, dims = 1)
   for (datasetname, datamean) in meandict
      push!(monitor, MonitoringItem(monitormeandiff, epoch,
               norm(samplemean - datamean), datasetname))
   end

   monitor
end

function monitormeandiffpervariable!(monitor::Monitor, bm::AbstractBM, epoch::Int,
      meandict::DataDict;
      nparticles::Int = 3000, burnin::Int = 10,
      xgenerated::Matrix{Float64} = sampleparticles(bm, nparticles, burnin)[1],
      variables::Vector{Int} = collect(1:size(xgenerated, 2)))

   samplemean = mean(xgenerated, dims = 1)
   for (datasetname, datamean) in meandict
      for i in variables
         push!(monitor, MonitoringItem(monitormeandiffpervariable, epoch,
               samplemean[i] - datamean[i], datasetname * "/Var" * string(i)))
      end
   end

   monitor
end

function monitormeanandcordiff!(monitor::Monitor, bm::AbstractBM,
      epoch::Int;
      meandict::DataDict = DataDict(),
      cordict::DataDict = DataDict(),
      nparticles::Int = 3000, burnin::Int = 10,
      xgenerated::Matrix{Float64} = sampleparticles(bm, nparticles, burnin)[1])

   monitormeandiff!(monitor, bm, epoch, meandict; xgenerated = xgenerated)
   monitorcordiff!(monitor, bm, epoch, cordict; xgenerated = xgenerated)
end

"""
    monitorreconstructionerror!(monitor, rbm, epoch, datadict)
Computes the reconstruction error for the data sets in the `datadict`
and the `rbm` and stores the values in the `monitor`.
"""
function monitorreconstructionerror!(monitor::Monitor, rbm::AbstractRBM,
      epoch::Int, datadict::DataDict)

   for (datasetname, x) in datadict
      push!(monitor, MonitoringItem(monitorreconstructionerror, epoch,
               reconstructionerror(rbm, x), datasetname))
   end

   monitor
end

function monitorsd!(monitor::Monitor,
      gbrbm::Union{GaussianBernoulliRBM, GaussianBernoulliRBM2}, epoch::Int)

   nvisible = length(gbrbm.sd)
   for i = 1:nvisible
      push!(monitor,
            MonitoringItem(monitorsd, epoch, gbrbm.sd[i], string(i)))
   end

   monitor
end

"""
    monitorweightsnorm!(monitor, rbm, epoch)
Computes the L2-norm of the weights matrix and the bias vectors of the `rbm`
and stores the values in the `monitor`.
These values can give a hint how much the updates are changing the parameters
during learning.
"""
function monitorweightsnorm!(monitor::Monitor, rbm::AbstractRBM, epoch::Int)
   push!(monitor,
         MonitoringItem(monitorweightsnorm, epoch,
               norm(rbm.weights), "Weights"),
         MonitoringItem(monitorweightsnorm, epoch,
               norm(rbm.visbias), "Visible bias"),
         MonitoringItem(monitorweightsnorm, epoch,
               norm(rbm.hidbias), "Hidden bias"))
end


text:
function p=pFromTheta(Th,Param)
p=(Param.Rd*Th/Param.p0^Param.kappa).^(1/(1-Param.kappa));
end


text:
using Documenter, MathPhysicalConstants

makedocs(
    modules = [MathPhysicalConstants],
    format = Documenter.HTML(),
    sitename = "MathPhysicalConstants",
    pages    = Any[
        "Introduction" => "index.md",
        "Installation" => "installation.md",
        "Usage"        => "usage.md",
        "api"          => "api.md",
    ]
)

deploydocs(
    repo = "github.com/LaGuer/MathPhysicalConstants.jl.git",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
##############################################################################
##
## Sorting
##
##############################################################################

#########################################
## Permutation & Ordering types/functions
#########################################
# Sorting in Julia works through Orderings, where each ordering is a type
# which defines a comparison function lt(o::Ord, a, b).

# UserColOrdering: user ordering of a column; this is just a convenience container
#                  which allows a user to specify column specific orderings
#                  with "order(column, rev=true,...)"

type UserColOrdering{T<:ColumnIndex}
    col::T
    kwargs
end

# This is exported, and lets a user define orderings for a particular column
order{T<:ColumnIndex}(col::T; kwargs...) = UserColOrdering{T}(col, kwargs)

# Allow getting the column even if it is not wrapped in a UserColOrdering
_getcol(o::UserColOrdering) = o.col
_getcol(x) = x

###
# Get an Ordering for a single column
###
function ordering(col_ord::UserColOrdering, lt::Function, by::Function, rev::Bool, order::Ordering)
    for (k,v) in col_ord.kwargs
        if     k == :lt;    lt    = v
        elseif k == :by;    by    = v
        elseif k == :rev;   rev   = v
        elseif k == :order; order = v
        else
            error("Unknown keyword argument: ", string(k))
        end
    end

    Order.ord(lt,by,rev,order)
end

ordering(col::ColumnIndex, lt::Function, by::Function, rev::Bool, order::Ordering) =
             Order.ord(lt,by,rev,order)

# DTPerm: defines a permutation on a particular DataTable, using
#         a single ordering (O<:Ordering) or a list of column orderings
#         (O<:AbstractVector{Ordering}), one per DataTable column
#
#         If a user only specifies a few columns, the DataTable
#         contained in the DTPerm only contains those columns, and
#         the permutation induced by this ordering is used to
#         sort the original (presumably larger) DataTable

immutable DTPerm{O<:Union{Ordering, AbstractVector}, DT<:AbstractDataTable} <: Ordering
    ord::O
    dt::DT
end

function DTPerm{O<:Ordering, DT<:AbstractDataTable}(ords::AbstractVector{O}, dt::DT)
    if length(ords) != ncol(dt)
        error("DTPerm: number of column orderings does not equal the number of DataTable columns")
    end
    DTPerm{typeof(ords), DT}(ords, dt)
end

DTPerm{O<:Ordering, DT<:AbstractDataTable}(o::O, dt::DT) = DTPerm{O,DT}(o,dt)

# get ordering function for the i-th column used for ordering
col_ordering{O<:Ordering}(o::DTPerm{O}, i::Int) = o.ord
col_ordering{V<:AbstractVector}(o::DTPerm{V}, i::Int) = o.ord[i]

Base.@propagate_inbounds Base.getindex(o::DTPerm, i::Int, j::Int) = o.dt[i, j]
Base.@propagate_inbounds Base.getindex(o::DTPerm, a::DataTableRow, j::Int) = a[j]

function Sort.lt(o::DTPerm, a, b)
    @inbounds for i = 1:ncol(o.dt)
        ord = col_ordering(o, i)
        va = o[a, i]
        vb = o[b, i]
        lt(ord, va, vb) && return true
        lt(ord, vb, va) && return false
    end
    false # a and b are equal
end

###
# Get an Ordering for a DataTable
###

################
## Case 1: no columns requested, so sort by all columns using requested order
################
## Case 1a: single order
######
ordering(dt::AbstractDataTable, lt::Function, by::Function, rev::Bool, order::Ordering) =
    DTPerm(Order.ord(lt, by, rev, order), dt)

######
## Case 1b: lt, by, rev, and order are Arrays
######
function ordering{S<:Function, T<:Function}(dt::AbstractDataTable,
                                            lt::AbstractVector{S}, by::AbstractVector{T},
                                            rev::AbstractVector{Bool}, order::AbstractVector)
    if !(length(lt) == length(by) == length(rev) == length(order) == size(dt,2))
        throw(ArgumentError("Orderings must be specified for all DataTable columns"))
    end
    DTPerm([Order.ord(_lt, _by, _rev, _order) for (_lt, _by, _rev, _order) in zip(lt, by, rev, order)], dt)
end

################
## Case 2:  Return a regular permutation when there's only one column
################
## Case 2a: The column is given directly
######
ordering(dt::AbstractDataTable, col::ColumnIndex, lt::Function, by::Function, rev::Bool, order::Ordering) =
    Perm(Order.ord(lt, by, rev, order), dt[col])

######
## Case 2b: The column is given as a UserColOrdering
######
ordering(dt::AbstractDataTable, col_ord::UserColOrdering, lt::Function, by::Function, rev::Bool, order::Ordering) =
    Perm(ordering(col_ord, lt, by, rev, order), dt[col_ord.col])

################
## Case 3:  General case: cols is an iterable of a combination of ColumnIndexes and UserColOrderings
################
## Case 3a: None of lt, by, rev, or order is an Array
######
function ordering(dt::AbstractDataTable, cols::AbstractVector, lt::Function, by::Function, rev::Bool, order::Ordering)

    if length(cols) == 0
        return ordering(dt, lt, by, rev, order)
    end

    if length(cols) == 1
        return ordering(dt, cols[1], lt, by, rev, order)
    end

    # Collect per-column ordering info

    ords = Ordering[]
    newcols = Int[]

    for col in cols
        push!(ords, ordering(col, lt, by, rev, order))
        push!(newcols, index(dt)[(_getcol(col))])
    end

    # Simplify ordering when all orderings are the same
    if all([ords[i] == ords[1] for i = 2:length(ords)])
        return DTPerm(ords[1], dt[newcols])
    end

    return DTPerm(ords, dt[newcols])
end

######
# Case 3b: cols, lt, by, rev, and order are all arrays
######
function ordering{S<:Function, T<:Function}(dt::AbstractDataTable, cols::AbstractVector,
                                            lt::AbstractVector{S}, by::AbstractVector{T},
                                            rev::AbstractVector{Bool}, order::AbstractVector)

    if !(length(lt) == length(by) == length(rev) == length(order))
        throw(ArgumentError("All ordering arguments must be 1 or the same length."))
    end

    if length(cols) == 0
        return ordering(dt, lt, by, rev, order)
    end

    if length(lt) != length(cols)
        throw(ArgumentError("All ordering arguments must be 1 or the same length as the number of columns requested."))
    end

    if length(cols) == 1
        return ordering(dt, cols[1], lt[1], by[1], rev[1], order[1])
    end

    # Collect per-column ordering info

    ords = Ordering[]
    newcols = Int[]

    for i in 1:length(cols)
        push!(ords, ordering(cols[i], lt[i], by[i], rev[i], order[i]))
        push!(newcols, index(dt)[(_getcol(cols[i]))])
    end

    # Simplify ordering when all orderings are the same
    if all([ords[i] == ords[1] for i = 2:length(ords)])
        return DTPerm(ords[1], dt[newcols])
    end

    return DTPerm(ords, dt[newcols])
end

######
## At least one of lt, by, rev, or order is an array or tuple, so expand all to arrays
######
function ordering(dt::AbstractDataTable, cols::AbstractVector, lt, by, rev, order)
    to_array(src::AbstractVector, dims) = src
    to_array(src::Tuple, dims) = [src...]
    to_array(src, dims) = fill(src, dims)

    dims = length(cols) > 0 ? length(cols) : size(dt,2)
    ordering(dt, cols,
             to_array(lt, dims),
             to_array(by, dims),
             to_array(rev, dims),
             to_array(order, dims))
end

#### Convert cols from tuple to Array, if necessary
ordering(dt::AbstractDataTable, cols::Tuple, args...) = ordering(dt, [cols...], args...)

###########################
# Default sorting algorithm
###########################

# TimSort is fast for data with structure, but only if the DataTable is large enough
# TODO: 8192 is informed but somewhat arbitrary

Sort.defalg(dt::AbstractDataTable) = size(dt, 1) < 8192 ? Sort.MergeSort : SortingAlgorithms.TimSort

# For DataTables, we can choose the algorithm based on the column type and requested ordering
function Sort.defalg{T<:Real}(dt::AbstractDataTable, ::Type{T}, o::Ordering)
    # If we're sorting a single numerical column in forward or reverse,
    # RadixSort will generally be the fastest stable sort
    if isbits(T) && sizeof(T) <= 8 && (o==Order.Forward || o==Order.Reverse)
        SortingAlgorithms.RadixSort
    else
        Sort.defalg(dt)
    end
end
Sort.defalg(dt::AbstractDataTable,        ::Type,            o::Ordering) = Sort.defalg(dt)
Sort.defalg(dt::AbstractDataTable, col    ::ColumnIndex,     o::Ordering) = Sort.defalg(dt, eltype(dt[col]), o)
Sort.defalg(dt::AbstractDataTable, col_ord::UserColOrdering, o::Ordering) = Sort.defalg(dt, col_ord.col, o)
Sort.defalg(dt::AbstractDataTable, cols,                     o::Ordering) = Sort.defalg(dt)

function Sort.defalg(dt::AbstractDataTable, o::Ordering; alg=nothing, cols=[])
    alg != nothing && return alg
    Sort.defalg(dt, cols, o)
end

########################
## Actual sort functions
########################

Base.issorted(dt::AbstractDataTable; cols=Any[], lt=isless, by=identity, rev=false, order=Forward) =
    issorted(eachrow(dt), ordering(dt, cols, lt, by, rev, order))

# sort and sortperm functions

for s in [:(Base.sort), :(Base.sortperm)]
    @eval begin
        function $s(dt::AbstractDataTable; cols=Any[], alg=nothing,
                    lt=isless, by=identity, rev=false, order=Forward)
            if !(isa(by, Function) || eltype(by) <: Function)
                msg = "'by' must be a Function or a vector of Functions. Perhaps you wanted 'cols'."
                throw(ArgumentError(msg))
            end
            ord = ordering(dt, cols, lt, by, rev, order)
            _alg = Sort.defalg(dt, ord; alg=alg, cols=cols)
            $s(dt, _alg, ord)
        end
    end
end

Base.sort(dt::AbstractDataTable, a::Algorithm, o::Ordering) = dt[sortperm(dt, a, o),:]
Base.sortperm(dt::AbstractDataTable, a::Algorithm, o::Union{Perm,DTPerm}) = sort!([1:size(dt, 1);], a, o)
Base.sortperm(dt::AbstractDataTable, a::Algorithm, o::Ordering) = sortperm(dt, a, DTPerm(o,dt))


text:
module TestPackage

f(x) = 2x

export f

end


text:
using StackOverflow

"""
    checkmyprivileges(username::String)

Takes in a username and returns the privileges said user has.
"""
function checkmyprivileges(username::String = "")

    r = HTTP.request("GET", "https://api.stackexchange.com/2.2/users?order=desc&sort=reputation&inname=$(username)&site=stackoverflow")
    #TODO: May return multiple names.

    json_obj = StackOverflow.convert_HTTP_Response_To_JSON(r)

    for (k, v) in json_obj
        if occursin("items", k)
            for item in v
                for (key, value) in item
                    if key == "reputation"

                        println("Your current rep is $(value) and your privileges are as follows: \n")
                        if value >= 1

                            println("\ndescription: Ask a question or contribute an answer")
                            println("short_description: create posts")
                        end

                        if value >= 5

                            println("\ndescription: Discuss the site itself,  bugs, feedback, and governance")
                            println("short_description: participate in meta")
                        end

                        if value >= 10

                            println("\ndescription: Post more links, answer protected questions")
                            println("short_description: answer protected questions")

                            println("\ndescription: Create answers that can be easily edited by most users")
                            println("short_description: create wiki posts")
                        end

                        if value >= 15

                            println("\ndescription: Bring content to the attention of the community via flags")
                            println("short_description: flag posts")

                            println("\ndescription: Indicate when questions and answers are useful")
                            println("short_description: vote up")
                        end

                        if value >= 20

                            println("\ndescription: Participate in this site's chat rooms")
                            println("short_description: talk in chat")
                        end

                        if value >= 50

                            println("\ndescription: Leave comments on other people's posts")
                            println("short_description: comment everywhere")
                        end

                        if value >= 75

                            println("\ndescription: Offer some of your reputation as bounty on a question")
                            println("short_description: set bounties")
                        end

                        if value >= 100

                            println("\ndescription: Collaborate on the editing and improvement of wiki posts")
                            println("short_description: edit community wikis")

                            println("\ndescription: Create new chat rooms")
                            println("short_description: create chat rooms")
                        end

                        if value >= 125

                            println("\ndescription: Indicate when questions and answers are not useful")
                            println("short_description: vote down")
                        end

                        if value >= 200

                            println("\ndescription: Some ads are now automatically disabled")
                            println("short_description: see reduced ads")
                        end

                        if value >= 250

                            println("\ndescription: View and cast close/reopen votes on your own questions")
                            println("short_description: view close votes")
                        end

                        if value >= 500

                            println("\ndescription: Access first posts and late answers review queues")
                            println("short_description: access review queues")
                        end

                        if value >= 1000

                            println("\ndescription: See vote breakout, expandable usercard")
                            println("short_description: see votes, expandable usercard")

                            println("\ndescription: Create chat rooms where only specific users may talk")
                            println("short_description: create gallery chat rooms")
                        end

                        if value >= 1500

                            println("\ndescription: Add new tags to the site")
                            println("short_description: create new tags")
                        end

                        if value >= 2000

                            println("\ndescription: Edits to any question or answer are applied immediately")
                            println("short_description: edit questions and answers")
                        end

                        if value >= 2500

                            println("\ndescription: Decide which tags have the same meaning as others")
                            println("short_description: create tag synonyms")
                        end

                        if value >= 3000

                            println("\ndescription: Help decide whether posts are off-topic or duplicates")
                            println("short_description: cast close & reopen votes")
                        end

                        if value >= 5000

                            println("\ndescription: Approve edits to tag wikis made by regular users")
                            println("short_description: approve tag wiki edits")
                        end

                        if value >= 10000

                            println("\ndescription: Access reports, delete questions, review reviews")
                            println("short_description: access moderator tools")
                        end

                        if value >= 15000

                            println("\ndescription: Mark questions as protected")
                            println("short_description: protect questions")
                        end

                        if value >= 20000

                            println("\ndescription: Expanded editing, deletion and undeletion privileges")
                            println("short_description: access 'trusted user' tools")
                        end

                        if value >= 25000

                            println("\ndescription: Access to internal and Google site analytics")
                            println("short_description: access to site analytics")
                        end

                        if value == 0
                            println("You have no privileges currently.")
                        end

                    end
                end
            end
        end
    end

end

"""
    getprivileges()

Prints out the various privilege levels in Stack Overflow and the rep
required to have said privileges.
"""
function getprivileges()
    r = HTTP.request("GET", "https://api.stackexchange.com/2.2/privileges?site=stackoverflow")
    json_obj = Stackoverflow.convert_HTTP_Response_To_JSON(r)

    for (k, v) in json_obj
        if occursin("items", k)
            for item in v
                for (key, value) in item
                    println("$(key): $(value)")
                end
                println("\n")
            end
        end
    end

end


text:
include("types.jl")

function parse_rep(name)
    pps = 0.0
    ppm = 0.0
    gamerate = 0.0
    time = 0
    totalPieceActiveTime = 0
    totalPieceLocked = 0
    lps = 0.0
    lpm = 0.0
    totalLines = 0
    gmt = now()

    rplay_ok = false

    open(name) do f
        for line in readlines(f)
            if contains(line, "0.ruleopt.strRuleName=STANDARD-FAST-B")
                rplay_ok = true
            elseif contains(line, "timestamp.gmt")
                nss = replace(replace(line, r"timestamp.gmt=", ""), r"\n", "")
                df = "yyyy-mm-dd-HH-MM-SS"
                #=println(df)=#
                #=println(nss)=#
                #=println(DateTime(nss, df))=#
                df = DateTime(nss, df)
                #=println(Dates.format(df, "yyyy-mm-dd HH:MM:SS"))=#
                gmt = df
            elseif contains(line, "timestamp.time")
                #=nss = replace(line, r"timestamp.time=", "")=#
                #=ns  = chop(replace(nss, r"\\:", "."))=#

                #=ns_split = split(ns, ".")=#

                #=xh_, xm_, xs_ = ns_split=#

                #=xh = parse(Int32, xh_)=#
                #=xm = parse(Int32, xm_)=#
                #=xs = parse(Int32, xs_)=#

                #=println("time ", xh, xm, xs)=#
            elseif contains(line, "statistics.time=")
                nss = replace(replace(line, r"0.statistics.time=", ""), r"\n", "")
                time = parse(Int32, nss)
                #=println("time ", time)=#
            elseif contains(line, "statistics.totalPieceActiveTime=")
                nss = replace(replace(line, r"0.statistics.totalPieceActiveTime=", ""), r"\n", "")
                totalPieceActiveTime = parse(Int32, nss)
                #=println("totalPieceActiveTime ", totalPieceActiveTime)=#
            elseif contains(line, "statistics.gamerate=")
                nss = replace(replace(line, r"0.statistics.gamerate=", ""), r"\n", "")
                gamerate = parse(Float64, nss)
                #=println("gamerate ", gamerate)=#
            elseif contains(line, "statistics.totalPieceLocked=")
                nss = replace(replace(line, r"0.statistics.totalPieceLocked=", ""), r"\n", "")
                totalPieceLocked = parse(Int32, nss)
                #=println("pieces ", totalPieceLocked)=#
            elseif contains(line, "statistics.lpm=")
                nss = replace(replace(line, r"0.statistics.lpm=", ""), r"\n", "")
                lpm = parse(Float64, nss)
                #=println("lpm ", lpm)=#
            elseif contains(line, "statistics.lps=")
                nss = replace(replace(line, r"0.statistics.lps=", ""), r"\n", "")
                lps = parse(Float64, nss)
                #=println("lps ", lps)=#
            elseif contains(line, "statistics.pps=")
                nss = replace(replace(line, r"0.statistics.pps=", ""), r"\n", "")
                pps = parse(Float64, nss)
                #=println("pps ", pps)=#
            elseif contains(line, "statistics.ppm=")
                nss = replace(replace(line, r"0.statistics.ppm=", ""), r"\n", "")
                ppm = parse(Float64, nss)
                #=println("ppm ", ppm)=#
            elseif contains(line, "result.totallines=")
                nss = replace(line, r"result.totallines=", "")
                totalLines = parse(Int32, nss)
                #=println("pps ", totalLines)=#
            end
        end
    end

    return game_data(rplay_ok, pps, ppm, gamerate, lps, lpm, time, totalPieceActiveTime, totalPieceLocked, totalLines, gmt)
end

function main()
    name = ARGS[1]

    base_pwd = pwd()
    cd(name)

    first = true
    best = 0

    for file in readdir(name)
        if contains(file, ".rep")
            data = parse_rep(file)

            if data.ok && data.totalLines == 40
                if first
                    first = false
                    best = data.totalPieceActiveTime / 60.0
                end

                best = min(best, data.totalPieceActiveTime / 60.0)

                #=println(data)=#
                @printf("%s %f %f %f %f\n", Dates.format(data.gmt, "yyyy-mm-dd-HH:MM:SS"), data.totalPieceActiveTime / 60.0, data.pps, data.lps, best)
            end
        end
    end

    cd(base_pwd)
end

main()


text:
using Test
using SymArrays
using SymArrays: symarrlength, sub2ind
using TensorOperations
using Random

@testset "SymArrays.jl" begin
    @testset "SymArray" begin
        # Write your own tests here.
        @test symarrlength((3,6,4,3),(3,2,1,3)) == 8400

        S = SymArray{(2,),Float64}(5,5);
        @test sub2ind(S,2,5) == sub2ind(S,5,2)
        @test sub2ind(S,3,4) != sub2ind(S,5,3)
        # sub2ind has to have same number of arguments as size of N
        @test_throws MethodError sub2ind(S,3,5,1)
        # indexing the array allows having additional 1s at the end
        @test S[3,5,1] == S[3,5]
        @test_throws BoundsError S[3,5,3]

        S = SymArray{(3,1,2,2),Float64}(3,3,3,2,4,4,4,4);
        @test size(S) == (3,3,3,2,4,4,4,4)
        @test length(S) == 2000
        # iterating over all indices should give only the distinct indices,
        # i.e., give the same number of terms as the array length
        @test sum(1 for s in S) == length(S)
        @test sum(1 for I in eachindex(S)) == length(S)

        # calculating the linear index when iterating over Cartesian indices should give sequential access to the array
        @test 1:length(S) == [(sub2ind(S,Tuple(I)...) for I in eachindex(S))...]

        # test that permuting exchangeable indices accesses the same array element
        i1  = sub2ind(S,1,2,3,1,4,3,2,1)
        @test sub2ind(S,2,3,1,1,4,3,2,1) == i1
        @test sub2ind(S,2,3,1,1,3,4,2,1) == i1
        @test sub2ind(S,2,3,1,1,3,4,1,2) == i1
        @test sub2ind(S,2,1,3,1,4,3,1,2) == i1
        # make sure that swapping independent indices gives a different array element
        @test sub2ind(S,2,1,3,1,4,1,3,2) != i1
        @test sub2ind(S,1,2,3,2,4,3,2,1) != i1
        @test sub2ind(S,1,2,3,1,1,3,2,4) != i1

        SI = eachindex(S)
        @test first(SI) == CartesianIndex(1,1,1,1,1,1,1,1)

        A = rand(5,5)
        @assert A' != A
        @test SymArray{(1,1)}(A) == A
        @test SymArray{(2,)}(A) != A

        A = A + A'
        # the standard equality test uses iteration over the arrays, (a,b) in zip(A,B),
        # which only accesses the actually different indices in SymArrays
        @test SymArray{(2,)}(A) != A
        # broadcasting goes over all indices
        @test all(SymArray{(2,)}(A) .== A)

        # views on existing vector-like types
        x = rand(3)
        S = SymArray{(2,)}(x,2,2)
        @test S.data === x

        x = 5:10
        S = SymArray{(2,)}(x,3,3)
        @test S.data === x

        A = rand(10)
        S = SymArray{(1,)}(A,10)
        @test S.data === A

        # but if called without sizes, the copy constructor should be used
        S = SymArray{(1,)}(A)
        @test S.data !== A
    end

    @testset "Contractions" begin
        N,M,O = 10, 15, 18
        for T in (Float64,ComplexF64)
            A = rand(T,N)

            for U in (Float64,ComplexF64)
                for (n,dims) in enumerate([(N,M,O),(O,N,M),(M,O,N)])
                    B = rand(U,dims...)
                    n==1 && (@tensor D1[j,k] := A[i]*B[i,j,k])
                    n==2 && (@tensor D1[j,k] := A[i]*B[j,i,k])
                    n==3 && (@tensor D1[j,k] := A[i]*B[j,k,i])
                    D2 = contract(A,B,Val(n))
                    @test D1 ≈ D2
                    if T==U
                        contract!(D2,A,B,Val(n))
                        @test D1 ≈ D2
                    else
                        @test_throws MethodError contract!(D2,A,B,Val(n))
                    end
                end

                S = SymArray{(3,),U}(N,N,N)
                rand!(S.data)
                B = collect(S)
                @test contract(A,S,Val(1)) == contract(A,S,Val(2))
                @test contract(A,S,Val(1)) == contract(A,S,Val(3))
                for n = 1:3
                    @test collect(contract(A,S,Val(n))) ≈ contract(A,B,Val(n))
                end

                S = SymArray{(2,1),U}(N,N,N)
                rand!(S.data)
                B = collect(S)
                @test contract(A,S,Val(1)) == contract(A,S,Val(2))
                for n = 1:3
                    @test collect(contract(A,S,Val(n))) ≈ contract(A,B,Val(n))
                end
                @test !(collect(contract(A,S,Val(1))) ≈ collect(contract(A,S,Val(3))))

                S = SymArray{(2,),U}(N,N)
                rand!(S.data)
                B = collect(S)
                @test contract(A,S,Val(1)) ≈ contract(A,S,Val(2))
                for n = 1:2
                    @test collect(contract(A,S,Val(n))) ≈ contract(A,B,Val(n))
                end
            end
        end
    end
end


text:
function estimator_statistics(
    estimate;
    dataset,
    demo_gt,
    position_indices,
    trajectory_distance = Distances.meanad,
    parameter_distance = Distances.cosine_dist,
)

    function trajectory_component_errors(trajectory)
        if haskey(trajectory, :x)
            (;
                x_error = trajectory_distance(demo_gt.x, trajectory.x),
                position_error = trajectory_distance(
                    demo_gt.x[position_indices, :],
                    trajectory.x[position_indices, :],
                ),
            )
        else
            (; x_error = Inf, position_error = Inf)
        end
    end

    observation = dataset[estimate.observation_idx]
    x_observation_error, position_observation_error = trajectory_component_errors(observation)
    x_estimation_error, position_estimation_error = trajectory_component_errors(estimate)

    parameter_estimation_error =
        map(demo_gt.player_cost_models_gt, estimate.player_weights) do cost_model_gt, weights_est
            @assert sum(weights_est) ≈ 1
            parameter_distance(CostUtils.normalize(cost_model_gt.weights), weights_est)
        end |> Statistics.mean

    (;
        estimate.estimator_name,
        estimate.observation_idx,
        estimate.converged,
        observation.σ,
        x_observation_error,
        position_observation_error,
        x_estimation_error,
        position_estimation_error,
        parameter_estimation_error,
    )
end


text:
type VcovSimple <: AbstractVcovMethod end
type VcovSimpleData <: AbstractVcovMethodData end
VcovMethodData(::VcovSimple, ::AbstractDataFrame) = VcovSimpleData()
function vcov!(::VcovSimpleData, x::VcovData)
    invcrossmatrix = inv(x.crossmatrix)
    scale!(invcrossmatrix, sumabs2(x.residuals) /  x.df_residual)
    return invcrossmatrix
end
shat!(::VcovSimpleData, x::VcovData) = scale(x.crossmatrix, sumabs2(x.residuals))


text:
using QXContexts
using Test
using TestSetExtensions
using Logging

@testset ExtendedTestSet "QXContexts.jl" begin
    @includetests ARGS
end


text:
export Marginal, getdata, is_clamped, is_initial, as_marginal
export SkipClamped, SkipInitial, SkipClampedAndInitial, IncludeAll

using Distributions
using Rocket

import Rocket: getrecent
import Base: ndims, precision, length, size

struct Marginal{D}
    data       :: D
    is_clamped :: Bool
    is_initial :: Bool
end

Base.show(io::IO, marginal::Marginal) = print(io, string("Marginal(", getdata(marginal), ")"))

getdata(marginal::Marginal)    = marginal.data
is_clamped(marginal::Marginal) = marginal.is_clamped
is_initial(marginal::Marginal) = marginal.is_initial

# TupleTools.prod is a more efficient version of Base.all for NTuple here
is_clamped(marginals::Tuple)  = TupleTools.prod(map(is_clamped, marginals))
is_initial(marginals::Tuple)  = TupleTools.prod(map(is_initial, marginals))

## Statistics

Distributions.pdf(marginal::Marginal, x)    = Distributions.pdf(getdata(marginal), x)
Distributions.logpdf(marginal::Marginal, x) = Distributions.logpdf(getdata(marginal), x)

MacroHelpers.@proxy_methods Marginal getdata [
    Distributions.mean,
    Distributions.median,
    Distributions.mode,
    Distributions.shape,
    Distributions.scale,
    Distributions.rate,
    Distributions.var,
    Distributions.std,
    Distributions.cov,
    Distributions.invcov,
    Distributions.logdetcov,
    Distributions.entropy,
    Distributions.params,
    Base.precision,
    Base.length,
    Base.ndims,
    Base.size,
    Base.eltype,
    mean_cov,
    mean_var,
    mean_invcov,
    mean_precision,
    weightedmean_cov,
    weightedmean_var,
    weightedmean_invcov,
    weightedmean_precision,
    probvec,
    weightedmean,
    inversemean,
    logmean,
    meanlogmean,
    mirroredlogmean,
    loggammamean
]

## Utility functions

getdata(marginals::Tuple) = map(getdata, marginals)

as_marginal(marginal::Marginal) = marginal

## Marginal observable

abstract type MarginalSkipStrategy end

struct SkipClamped           <: MarginalSkipStrategy end
struct SkipInitial           <: MarginalSkipStrategy end
struct SkipClampedAndInitial <: MarginalSkipStrategy end
struct IncludeAll            <: MarginalSkipStrategy end

apply_skip_filter(observable, ::SkipClamped)           = observable |> filter(v -> !is_clamped(v))
apply_skip_filter(observable, ::SkipInitial)           = observable |> filter(v -> !is_initial(v))
apply_skip_filter(observable, ::SkipClampedAndInitial) = observable |> filter(v -> !is_initial(v) && !is_clamped(v))
apply_skip_filter(observable, ::IncludeAll)            = observable

struct MarginalObservable <: Subscribable{ Marginal }
    subject :: Rocket.RecentSubjectInstance{ Marginal, Subject{ Marginal, AsapScheduler, AsapScheduler } }
    stream  :: LazyObservable{ Marginal }
end

MarginalObservable() = MarginalObservable(RecentSubject(Marginal), lazy(Marginal))

as_marginal_observable(observable::MarginalObservable, skip_strategy::MarginalSkipStrategy) = apply_skip_filter(observable, skip_strategy)
as_marginal_observable(observable)                                                          = as_marginal_observable(observable, IncludeAll())

function as_marginal_observable(observable, skip_strategy::MarginalSkipStrategy)
    output = MarginalObservable()
    connect!(output, observable)
    return as_marginal_observable(output, skip_strategy)
end

Rocket.getrecent(observable::MarginalObservable) = Rocket.getrecent(observable.subject)
Rocket.getrecent(observables::Tuple)             = Rocket.getrecent.(observables)
Rocket.getrecent(::Nothing)                      = nothing

@inline Rocket.on_subscribe!(observable::MarginalObservable, actor) = subscribe!(observable.stream, actor)

@inline Rocket.subscribe!(observable::MarginalObservable, actor::Rocket.Actor{ <: Marginal })           = Rocket.on_subscribe!(observable.stream, actor)
@inline Rocket.subscribe!(observable::MarginalObservable, actor::Rocket.NextActor{ <: Marginal })       = Rocket.on_subscribe!(observable.stream, actor)
@inline Rocket.subscribe!(observable::MarginalObservable, actor::Rocket.ErrorActor{ <: Marginal })      = Rocket.on_subscribe!(observable.stream, actor)
@inline Rocket.subscribe!(observable::MarginalObservable, actor::Rocket.CompletionActor{ <: Marginal }) = Rocket.on_subscribe!(observable.stream, actor)

@inline Rocket.subscribe!(observable::MarginalObservable, actor::Rocket.Subject{ <: Marginal })                 = Rocket.on_subscribe!(observable.stream, actor)
@inline Rocket.subscribe!(observable::MarginalObservable, actor::Rocket.BehaviorSubjectInstance{ <: Marginal }) = Rocket.on_subscribe!(observable.stream, actor)
@inline Rocket.subscribe!(observable::MarginalObservable, actor::Rocket.PendingSubjectInstance{ <: Marginal })  = Rocket.on_subscribe!(observable.stream, actor)
@inline Rocket.subscribe!(observable::MarginalObservable, actor::Rocket.RecentSubjectInstance{ <: Marginal })   = Rocket.on_subscribe!(observable.stream, actor)
@inline Rocket.subscribe!(observable::MarginalObservable, actor::Rocket.ReplaySubjectInstance{ <: Marginal })   = Rocket.on_subscribe!(observable.stream, actor)

function connect!(marginal::MarginalObservable, source)
    set!(marginal.stream, source |> multicast(marginal.subject) |> ref_count())
    return nothing
end

function setmarginal!(marginal::MarginalObservable, value)
    next!(marginal.subject, Marginal(value, false, true))
    return nothing
end

## Marginal Mapping structure
## https://github.com/JuliaLang/julia/issues/42559
## Explanation: Julia cannot fully infer type of the lambda callback function in activate! method in node.jl file
## We create a lambda-like callable structure to improve type inference and make it more stable
## However it is not fully inferrable due to dynamic tags and variable constraints, but still better than just a raw lambda callback

struct MarginalMapping{F, T, N, M, A, R}
    vtag            :: T
    msgs_names      :: N
    marginals_names :: M
    meta            :: A
    factornode      :: R
end

marginal_mapping_fform(::MarginalMapping{F}) where F = F
marginal_mapping_fform(::MarginalMapping{F}) where F <: Function = F.instance

function MarginalMapping(::Type{F}, vtag::T, msgs_names::N, marginals_names::M, meta::A, factornode::R) where { F, T, N, M, A, R }
    return MarginalMapping{F, T, N, M, A, R}(vtag, msgs_names, marginals_names, meta, factornode)
end

function MarginalMapping(::F, vtag::T, msgs_names::N, marginals_names::M, meta::A, factornode::R) where { F <: Function, T, N, M, A, R}
    return MarginalMapping{F, T, N, M, A, R}(vtag, msgs_names, marginals_names, meta, factornode)
end

function (mapping::MarginalMapping)(dependencies)

    messages  = dependencies[1]
    marginals = getrecent(dependencies[2])

    # Marginal is clamped if all of the inputs are clamped
    is_marginal_clamped = __check_all(is_clamped, messages) && __check_all(is_clamped, marginals)

    # Marginal is initial if it is not clamped and all of the inputs are either clamped or initial
    is_marginal_initial = !is_marginal_clamped && (__check_all(is_clamped_or_initial, messages) && __check_all(is_clamped_or_initial, marginals))

    marginal = marginalrule(
        marginal_mapping_fform(mapping),
        mapping.vtag,
        mapping.msgs_names,
        messages,
        mapping.marginals_names,
        marginals,
        mapping.meta,
        mapping.factornode
    )

    return Marginal(marginal, is_marginal_clamped, is_marginal_initial)
end

Base.map(::Type{T}, mapping::M) where { T, M <: MarginalMapping } = Rocket.MapOperator{T, M}(mapping)


text:
using Documenter, SDFReader

makedocs(
    sitename = "SDFReader",
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    pages = [
        "index.md",
    ],
    modules = [SDFReader]
)

deploydocs(
    repo = "github.com/SebastianM-C/SDFReader.jl",
    push_preview = true
)


text:
using ITensors
using ITensorTDVP
using KrylovKit
using Observers
using Random
using Test

@testset "Basic TDVP" begin
  N = 10
  cutoff = 1e-12

  s = siteinds("S=1/2", N)

  os = OpSum()
  for j in 1:(N - 1)
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
    os += "Sz", j, "Sz", j + 1
  end

  H = MPO(os, s)

  ψ0 = randomMPS(s; linkdims=10)

  # Time evolve forward:
  ψ1 = tdvp(H, -0.1im, ψ0; nsweeps=1, cutoff, nsite=1)

  @test ψ1 ≈ tdvp(-0.1im, H, ψ0; nsweeps=1, cutoff, nsite=1)
  @test ψ1 ≈ tdvp(H, ψ0, -0.1im; nsweeps=1, cutoff, nsite=1)

  @test norm(ψ1) ≈ 1.0

  ## Should lose fidelity:
  #@test abs(inner(ψ0,ψ1)) < 0.9

  # Average energy should be conserved:
  @test real(inner(ψ1', H, ψ1)) ≈ inner(ψ0', H, ψ0)

  # Time evolve backwards:
  ψ2 = tdvp(H, +0.1im, ψ1; nsweeps=1, cutoff)

  @test norm(ψ2) ≈ 1.0

  # Should rotate back to original state:
  @test abs(inner(ψ0, ψ2)) > 0.99
end

@testset "TDVP: Sum of Hamiltonians" begin
  N = 10
  cutoff = 1e-10

  s = siteinds("S=1/2", N)

  os1 = OpSum()
  for j in 1:(N - 1)
    os1 += 0.5, "S+", j, "S-", j + 1
    os1 += 0.5, "S-", j, "S+", j + 1
  end
  os2 = OpSum()
  for j in 1:(N - 1)
    os2 += "Sz", j, "Sz", j + 1
  end

  H1 = MPO(os1, s)
  H2 = MPO(os2, s)
  Hs = [H1, H2]

  ψ0 = randomMPS(s; linkdims=10)

  ψ1 = tdvp(Hs, -0.1im, ψ0; nsweeps=1, cutoff, nsite=1)

  @test ψ1 ≈ tdvp(-0.1im, Hs, ψ0; nsweeps=1, cutoff, nsite=1)
  @test ψ1 ≈ tdvp(Hs, ψ0, -0.1im; nsweeps=1, cutoff, nsite=1)

  @test norm(ψ1) ≈ 1.0

  ## Should lose fidelity:
  #@test abs(inner(ψ0,ψ1)) < 0.9

  # Average energy should be conserved:
  @test real(sum(H -> inner(ψ1', H, ψ1), Hs)) ≈ sum(H -> inner(ψ0', H, ψ0), Hs)

  # Time evolve backwards:
  ψ2 = tdvp(Hs, +0.1im, ψ1; nsweeps=1, cutoff)

  @test norm(ψ2) ≈ 1.0

  # Should rotate back to original state:
  @test abs(inner(ψ0, ψ2)) > 0.99
end

@testset "Custom solver in TDVP" begin
  N = 10
  cutoff = 1e-12

  s = siteinds("S=1/2", N)

  os = OpSum()
  for j in 1:(N - 1)
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
    os += "Sz", j, "Sz", j + 1
  end

  H = MPO(os, s)

  ψ0 = randomMPS(s; linkdims=10)

  function solver(PH, t, psi0; kwargs...)
    solver_kwargs = (;
      ishermitian=true, tol=1e-12, krylovdim=30, maxiter=100, verbosity=0, eager=true
    )
    psi, info = exponentiate(PH, t, psi0; solver_kwargs...)
    return psi, info
  end

  ψ1 = tdvp(solver, H, -0.1im, ψ0; cutoff, nsite=1)

  @test ψ1 ≈ tdvp(solver, -0.1im, H, ψ0; cutoff, nsite=1)
  @test ψ1 ≈ tdvp(solver, H, ψ0, -0.1im; cutoff, nsite=1)

  @test norm(ψ1) ≈ 1.0

  ## Should lose fidelity:
  #@test abs(inner(ψ0,ψ1)) < 0.9

  # Average energy should be conserved:
  @test real(inner(ψ1', H, ψ1)) ≈ inner(ψ0', H, ψ0)

  # Time evolve backwards:
  ψ2 = tdvp(H, +0.1im, ψ1; cutoff)

  @test norm(ψ2) ≈ 1.0

  # Should rotate back to original state:
  @test abs(inner(ψ0, ψ2)) > 0.99
end

@testset "Accuracy Test" begin
  N = 4
  tau = 0.1
  ttotal = 1.0
  cutoff = 1e-12

  s = siteinds("S=1/2", N; conserve_qns=false)

  os = OpSum()
  for j in 1:(N - 1)
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
    os += "Sz", j, "Sz", j + 1
  end
  H = MPO(os, s)
  HM = prod(H)

  Ut = exp(-im * tau * HM)

  psi = productMPS(s, n -> isodd(n) ? "Up" : "Dn")
  psix = prod(psi)

  Sz_tdvp = Float64[]
  Sz_exact = Float64[]

  c = div(N, 2)
  Szc = op("Sz", s[c])

  Nsteps = Int(ttotal / tau)
  for step in 1:Nsteps
    psix = noprime(Ut * psix)
    psix /= norm(psix)

    psi = tdvp(
      H,
      -im * tau,
      psi;
      cutoff,
      normalize=false,
      solver_tol=1e-12,
      solver_maxiter=500,
      solver_krylovdim=25,
    )
    push!(Sz_tdvp, real(expect(psi, "Sz"; sites=c:c)[1]))
    push!(Sz_exact, real(scalar(dag(prime(psix, s[c])) * Szc * psix)))
    F = abs(scalar(dag(psix) * prod(psi)))
  end

  @test norm(Sz_tdvp - Sz_exact) < 1e-5
end

@testset "TEBD Comparison" begin
  N = 10
  cutoff = 1e-12
  tau = 0.1
  ttotal = 1.0

  s = siteinds("S=1/2", N; conserve_qns=true)

  os = OpSum()
  for j in 1:(N - 1)
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
    os += "Sz", j, "Sz", j + 1
  end

  H = MPO(os, s)

  gates = ITensor[]
  for j in 1:(N - 1)
    s1 = s[j]
    s2 = s[j + 1]
    hj =
      op("Sz", s1) * op("Sz", s2) +
      1 / 2 * op("S+", s1) * op("S-", s2) +
      1 / 2 * op("S-", s1) * op("S+", s2)
    Gj = exp(-1.0im * tau / 2 * hj)
    push!(gates, Gj)
  end
  append!(gates, reverse(gates))

  psi = productMPS(s, n -> isodd(n) ? "Up" : "Dn")
  phi = copy(psi)
  c = div(N, 2)

  #
  # Evolve using TEBD
  #

  Nsteps = convert(Int, ceil(abs(ttotal / tau)))
  Sz1 = zeros(Nsteps)
  En1 = zeros(Nsteps)
  Sz2 = zeros(Nsteps)
  En2 = zeros(Nsteps)

  for step in 1:Nsteps
    psi = apply(gates, psi; cutoff)
    #normalize!(psi)

    nsite = (step <= 3 ? 2 : 1)
    phi = tdvp(
      H, -tau * im, phi; nsweeps=1, cutoff, nsite, normalize=true, exponentiate_krylovdim=15
    )

    Sz1[step] = expect(psi, "Sz"; sites=c:c)[1]
    Sz2[step] = expect(phi, "Sz"; sites=c:c)[1]
    En1[step] = real(inner(psi', H, psi))
    En2[step] = real(inner(phi', H, phi))
  end

  #
  # Evolve using TDVP
  #
  struct TDVPObserver <: AbstractObserver end

  Sz2 = zeros(Nsteps)
  En2 = zeros(Nsteps)
  function ITensors.measure!(obs::TDVPObserver; sweep, bond, half_sweep, psi, kwargs...)
    if bond == 1 && half_sweep == 2
      Sz2[sweep] = expect(psi, "Sz"; sites=c)
      En2[sweep] = real(inner(psi', H, psi))
    end
  end

  phi = productMPS(s, n -> isodd(n) ? "Up" : "Dn")

  phi = tdvp(
    H,
    -im * ttotal,
    phi;
    time_step=-im * tau,
    cutoff,
    normalize=false,
    (observer!)=TDVPObserver(),
  )

  #display(En1)
  #display(En2)
  #display(Sz1)
  #display(Sz2)
  #@show norm(Sz1 - Sz2)
  #@show norm(En1 - En2)

  @test norm(Sz1 - Sz2) < 1e-3
  @test norm(En1 - En2) < 1e-3
end

@testset "Imaginary Time Evolution" begin
  N = 10
  cutoff = 1e-12
  tau = 1.0
  ttotal = 50.0

  s = siteinds("S=1/2", N)

  os = OpSum()
  for j in 1:(N - 1)
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
    os += "Sz", j, "Sz", j + 1
  end

  H = MPO(os, s)

  psi = randomMPS(s; linkdims=2)

  trange = 0.0:tau:ttotal
  for (step, t) in enumerate(trange)
    nsite = (step <= 10 ? 2 : 1)
    psi = tdvp(H, -tau, psi; cutoff, nsite, normalize=true, exponentiate_krylovdim=15)
  end
  #@show maxlinkdim(psi)

  @test inner(psi', H, psi) < -4.25
end

@testset "Observers" begin
  N = 10
  cutoff = 1e-12
  tau = 0.1
  ttotal = 1.0

  s = siteinds("S=1/2", N; conserve_qns=true)

  os = OpSum()
  for j in 1:(N - 1)
    os += 0.5, "S+", j, "S-", j + 1
    os += 0.5, "S-", j, "S+", j + 1
    os += "Sz", j, "Sz", j + 1
  end
  H = MPO(os, s)

  c = div(N, 2)

  #
  # Using the ITensors observer system
  #
  struct TDVPObserver <: AbstractObserver end

  Nsteps = convert(Int, ceil(abs(ttotal / tau)))
  Sz1 = zeros(Nsteps)
  En1 = zeros(Nsteps)
  function ITensors.measure!(obs::TDVPObserver; sweep, bond, half_sweep, psi, kwargs...)
    if bond == 1 && half_sweep == 2
      Sz1[sweep] = expect(psi, "Sz"; sites=c)
      En1[sweep] = real(inner(psi', H, psi))
    end
  end

  psi1 = productMPS(s, n -> isodd(n) ? "Up" : "Dn")
  tdvp(
    H,
    -im * ttotal,
    psi1;
    time_step=-im * tau,
    cutoff,
    normalize=false,
    (observer!)=TDVPObserver(),
  )

  #
  # Using Observers.jl
  #

  function measure_sz(; psi, bond, half_sweep)
    if bond == 1 && half_sweep == 2
      return expect(psi, "Sz"; sites=c)
    end
    return nothing
  end

  function measure_en(; psi, bond, half_sweep)
    if bond == 1 && half_sweep == 2
      return real(inner(psi', H, psi))
    end
    return nothing
  end

  obs = Observer("Sz" => measure_sz, "En" => measure_en)

  step_measure_sz(; psi) = expect(psi, "Sz"; sites=c)

  step_measure_en(; psi) = real(inner(psi', H, psi))

  step_obs = Observer("Sz" => step_measure_sz, "En" => step_measure_en)

  psi2 = MPS(s, n -> isodd(n) ? "Up" : "Dn")
  tdvp(
    H,
    -im * ttotal,
    psi2;
    time_step=-im * tau,
    cutoff,
    normalize=false,
    (observer!)=obs,
    (step_observer!)=step_obs,
  )

  Sz2 = results(obs)["Sz"]
  En2 = results(obs)["En"]

  Sz2_step = results(step_obs)["Sz"]
  En2_step = results(step_obs)["En"]

  @test Sz1 ≈ Sz2
  @test En1 ≈ En2
  @test Sz1 ≈ Sz2_step
  @test En1 ≈ En2_step
end

nothing


text:
function getPieceMoves!(moves::Moves, bitboard::UInt64, type::Symbol,
                        friends::UInt64, enemy::ChessSet, white::Bool, b::Board,
                        k_in_check::Bool=false)
    if bitboard == EMPTY; return end
    if type == :pawn
        getPawnMoves!(moves, bitboard, b.taken, friends, enemy, white, b.enpassant)
    else
        b0, d = bitshift(bitboard)
        for i = 0:d
            if b0<<i & bitboard != EMPTY
                src = b0<<i
                piece_moves = EMPTY

                if type == :knight
                    piece_moves = KNIGHT[src]
                elseif type == :king
                    piece_moves = KING[src]
                    if ~k_in_check; getCastlingMoves!(moves, src, b, white) end
                elseif type == :rook
                    piece_moves = getMagicAttack(src, b.taken, true)
                elseif type == :bishop
                    piece_moves = getMagicAttack(src, b.taken, false)
                elseif type == :queen
                    piece_moves = getMagicAttack(src, b.taken, true)
                    piece_moves |= getMagicAttack(src, b.taken, false)
                end
                piece_moves &= ~friends

                if piece_moves == EMPTY; continue end

                b0_m, d_m = bitshift(piece_moves)
                for j = 0:d_m
                    if b0_m<<j & piece_moves != EMPTY
                        target = b0_m<<j
                        enemy_type = getTypeAt(enemy, target)
                        if enemy_type == :king; continue end
                        target & enemy.friends != EMPTY ? take = Piece(enemy_type, target) : take = NONE
                        push!(moves, Move(type, src, target, take, EMPTY, :none, NOCASTLING))
                    end
                end
            end
        end
    end
end

function getAttack(attack::UInt64, b::UInt64, type::Symbol, taken::UInt64)
    if b == EMPTY; return attack end
    b0, d = bitshift(b)
    for i = 0:d
        if b0<<i & b != EMPTY
            src = b0<<i
            if type == :knight
                return attack |= KNIGHT[src]
            elseif type == :rook
                return attack |= getMagicAttack(src, taken, true)
            elseif type == :bishop
                return attack |= getMagicAttack(src, taken, false)
            elseif type == :queen
                attack |= getMagicAttack(src, taken, false)
                return attack |= getMagicAttack(src, taken, true)
            end
        end
    end
end

function getAttack(b::Board, white::Bool)
    white ? cs = b.white : cs = b.black
    attack = KING[cs.K]
    attack = getAttack(attack, cs.R, :rook, b.taken)
    attack = getAttack(attack, cs.B, :bishop, b.taken)
    attack = getAttack(attack, cs.Q, :queen, b.taken)
    return getPawnAttack(attack, cs.P, white)
end

function getMoves(b::Board, white::Bool)
    white ? friends = b.white.friends : friends = b.black.friends
    white ? enemy = b.black : enemy = b.white
    white ? cs = b.white : cs = b.black

    moves = Moves()
    king_in_check = inCheck(b, ~b.active)
    for (bitboard, s) in zip([cs.P, cs.N, cs.B, cs.R, cs.Q, cs.K],
                             [:pawn, :knight, :bishop, :rook, :queen, :king])
        getPieceMoves!(moves, bitboard, s, friends, enemy, white, b, king_in_check)
    end
    return filterMoves(b, moves)
end

function filterMoves(b::Board, moves::Moves)
    filtered = Moves()
    for m in moves.moves
        if m.type == :none || m.take.type == :king; continue end
        b1 = makeMove(b, m)
        king_in_check = inCheck(b1, b.active)
        if ~king_in_check
            push!(filtered, m)
        end
    end
    return filtered
end

function moveFromTo(bitboard, move)
    return bitboard ⊻= move.from | move.to
end

function updateSet(cs::ChessSet, move::Move)
    newp = cs.P
    newn = cs.N
    newb = cs.B
    newr = cs.R
    newq = cs.Q
    newk = cs.K
    if move.type == :pawn
        if move.promotion != :none
            newp = removeFrom(newp, move.from)
            if move.promotion == :queen
                newq |= move.to
            elseif move.promotion == :rook
                newr |= move.to
            elseif move.promotion == :knight
                newn |= move.to
            elseif move.promotion == :bishop
                newb |= move.to
            end
        else
            newp = moveFromTo(newp, move)
        end
    elseif move.type == :knight
        newn = moveFromTo(newn, move)
    elseif move.type == :bishop
        newb = moveFromTo(newb, move)
    elseif move.type == :rook
        newr = moveFromTo(newr, move)
    elseif move.type == :queen
        newq = moveFromTo(newq, move)
    elseif move.type == :king
        newk = moveFromTo(newk, move)
    end
    newf = newp|newn|newb|newr|newq|newk
    return ChessSet(newp, newn, newb, newr, newq, newk, newf)
end

function removeFrom(bitboard::UInt64, square::UInt64)
    return bitboard ⊻= square
end

function updateSet(cs::ChessSet, piece::Piece)
    newp = cs.P
    newn = cs.N
    newb = cs.B
    newr = cs.R
    newq = cs.Q
    newk = cs.K
    if piece.type == :pawn
        newp = removeFrom(newp, piece.square)
    elseif piece.type == :knight
        newn = removeFrom(newn, piece.square)
    elseif piece.type == :bishop
        newb = removeFrom(newb, piece.square)
    elseif piece.type == :rook
        newr = removeFrom(newr, piece.square)
    elseif piece.type == :queen
        newq = removeFrom(newq, piece.square)
    elseif piece.type == :king
        newk = removeFrom(newk, piece.square)
    end
    newf = newp|newn|newb|newr|newq|newk
    return ChessSet(newp, newn, newb, newr, newq, newk, newf)
end

function makeMove(board::Board, move::Move)
    castling = board.castling
    if board.active #true is white, false is black
        new_white = updateSet(board.white, move)

        if move.take != NONE
            new_black = updateSet(board.black, move.take)
            if move.take.type == :rook
                if move.take.square == A8
                    castling ⊻= Cq
                elseif move.take.square == H8
                    castling ⊻= Ck
                end
            end
        else
            new_black = board.black
        end

        if move.castling != NOCASTLING
            if move.castling == CQ
                new_white = updateSet(new_white, Move(:rook, A1, D1, NONE, EMPTY, :none, NOCASTLING))
                castling ⊻= CQ
            elseif move.castling == CK
                new_white = updateSet(new_white, Move(:rook, H1, F1, NONE, EMPTY, :none, NOCASTLING))
                castling ⊻= CK
            end
        end
        if castling&CK != NOCASTLING && (move.type == :king || (move.type == :rook && move.from == H1))
            castling ⊻= CK
        end
        if castling&CQ != NOCASTLING && (move.type == :king || (move.type == :rook && move.from == A1))
            castling ⊻= CQ
        end
    else
        new_black = updateSet(board.black, move)

        if move.take != NONE
            new_white = updateSet(board.white, move.take)
            if move.take.type == :rook
                if  move.take.square == A1
                    castling ⊻= CQ
                elseif move.take.square == H1
                    castling ⊻= CK
                end
            end
        else
            new_white = board.white
        end

        if move.castling != NOCASTLING
            if move.castling == Cq
                new_black = updateSet(new_black, Move(:rook, A8, D8, NONE, EMPTY, :none, NOCASTLING))
                castling ⊻= Cq
            elseif move.castling == Ck
                new_black = updateSet(new_black, Move(:rook, H8, F8, NONE, EMPTY, :none, NOCASTLING))
                castling ⊻= Ck
            end
        end
        if castling&Ck != NOCASTLING && (move.type == :king || (move.type == :rook && move.from == H8))
            castling ⊻= Ck
        end
        if castling&Cq != NOCASTLING && (move.type == :king || (move.type == :rook && move.from == A8))
            castling ⊻= Cq
        end
    end

    taken = new_white.friends | new_black.friends
    return Board(new_white, new_black, taken, ~board.active, castling,
                 move.enpassant, board.halfmove, board.fullmove)
end

function moveNow(b::Board, m::AbstractString)
    src = m[1:2]
    dst = m[3:4]
    if length(m) == 5
        prm = m[end]
    else
        prm = ""
    end

    avail_moves = Cassandra.getMoves(b, b.active)
    for am in avail_moves
        if Cassandra.UINT2PGN[am.from] == src && Cassandra.UINT2PGN[am.to] == dst
             if (prm == 'q' && am.promotion == :queen) ||
                (prm == 'r' && am.promotion == :rook) ||
                (prm == 'b' && am.promotion == :bishop) ||
                (prm == 'n' && am.promotion == :knight) ||
                (prm == "" && am.promotion == :none)
                b1 = Cassandra.makeMove(b, am)
                return b1
            end
        end
    end
end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

# Porting note: no need to store `hash` as we can just use the standard hash function
struct Label
    scope::Optional{String}
    name::String
end

Label(name::AbstractString) = Label(nothing, name)
Label(scope::AbstractString, name::AbstractString) = Label(scope, name)

function scoped_name(label::Label)
    return label.scope === nothing ?
        label.name :
        string(label.scope, ":", label.name)
end

Base.show(io::IO, label::Label) = print(io, scoped_name(label))


text:
@testset "Boltzmann" begin
    @test RTFOOL.boltzmann(1.0, 1.0) == e^-1.00
    @test RTFOOL.boltzmann(1.0, 2.0) == e^-2.00
    @test RTFOOL.boltzmann(0.5, 1.5) == e^-0.75

    @test RTFOOL.boltzmann(1.0, [1.0]) == [1.0]
    @test RTFOOL.boltzmann(1.0, [2.0]) == [1.0]
    @test RTFOOL.boltzmann(1.0, [1.0, 1.0]) == [0.5, 0.5]
    @test RTFOOL.boltzmann(1.0, [1.0, 2.0]) == [e/(1+e), 1/(1+e)]
    @test RTFOOL.boltzmann(0.5, [0.5, 1.5]) ≈ [e/(e + √e), √e/(e + √e)]
    @test RTFOOL.boltzmann(1/3, [1, 2, 3]) ==
        [e^-(1/3), e^-(2/3), e^-1] / (e^-(1/3) + e^-(2/3) + e^-1)

    @test RTFOOL.boltzmann(1.0, [1.0], [2]) == [1.0]
    @test RTFOOL.boltzmann(1.0, [2.0], [3]) == [1.0]
    @test RTFOOL.boltzmann(2.0, [1.0, 1.0], [1, 2]) == [1.0/3, 2.0/3]
    @test RTFOOL.boltzmann(1.0, [1.0, 2.0], [2, 1]) == [2e/(1+2e), 1/(1+2e)]
    @test RTFOOL.boltzmann(0.5, [0.5, 1.5], [2, 3]) ≈ [2e/(2e + 3*√e), 3*√e/(2e + 3*√e)]
    @test RTFOOL.boltzmann(1/3, [1, 2, 3], [1, 2, 3]) ==
        [e^-(1/3), 2e^-(2/3), 3e^-1] / (e^-(1/3) + 2e^-(2/3) + 3e^-1)

    let space = StateSpace{1}([(1,)], [1.0], BigInt[2])
        @test RTFOOL.boltzmann(1.0, space) == [1.0]
    end
    let space = StateSpace{1}([(1,)], [2.0], BigInt[3])
        @test RTFOOL.boltzmann(1.0, [2.0], [3]) == [1.0]
    end
    let space = StateSpace{2}([(1,1), (2,2)], [1.0, 1.0], BigInt[1,2])
        @test RTFOOL.boltzmann(2.0, space) == [1.0/3, 2.0/3]
    end
    let space = StateSpace{2}([(1,1),(2,2)], [1.0, 2.0], BigInt[2,1])
        @test RTFOOL.boltzmann(1.0, space) == [2e/(1+2e), 1/(1+2e)]
    end
    let space = StateSpace{2}([(1,1),(2,2)], [0.5, 1.5], BigInt[2,3])
        @test RTFOOL.boltzmann(0.5, space) ≈ [2e/(2e + 3*√e), 3*√e/(2e + 3*√e)]
    end
    let space = StateSpace{3}([(1,1,1), (1,1,2), (1,2,2)], [1,2,3], BigInt[1,2,3])
        @test RTFOOL.boltzmann(1/3, space) ≈
            [e^-(1/3), 2e^-(2/3), 3e^-1] / (e^-(1/3) + 2e^-(2/3) + 3e^-1)
    end
end


text:
function eftAddGTE{T<:StdFloat}(a::T,b::T,c::T,d::T)

   w=x=y=z=zero(T)

   w,d = eftAddGTE(c,d)
   w,c = eftAddGTE(b,w)
   a,b = eftAddGTE(a,w)

   w,x = a,b
   if x != zero(T)
       x,y = eftAddGTE(x,c)
       if y != zero(T)
           y,z = eftAddGTE(y,d)
       else
           x,y = eftAddGTE(x,d)
       end
   else
       w,x = eftAddGTE(w,c)
       if x != zero(T)
           x,y = eftAddGTE(x,d)
       else
           w,x = eftAddGTE(w,d)
       end
   end

   w,x,y,z
end

function eftAddGTEas3{T<:StdFloat}(a::T,b::T,c::T,d::T)

   w=x=y=zero(T)

   w,d = eftAddGTE(c,d)
   w,c = eftAddGTE(b,w)
   a,b = eftAddGTE(a,w)

   w,x = a,b
   if x != zero(T)
       x,y = eftAddGTE(x,c)
       if y != zero(T)
           y += d
       else
           x,y = eftAddGTE(x,d)
       end
   else
       w,x = eftAddGTE(w,c)
       if x != zero(T)
           x,y = eftAddGTE(x,d)
       else
           w,x = eftAddGTE(w,d)
       end
   end

   w,x,y
end

                                    # hi1 lo1  hi2  lo2
function eftAddGTEx2as2{T<:StdFloat}(a::T,b::T,c::T,d::T)

   w=zero(T)

   w,c = eftAdd(b,c)
   a,b = eftAdd(a,w)

   if b != zero(T)
       b += c + d  # return a, b+c+d
   else
       a,b = eftAddGTE(a,c)
       if b != zero(T)
           b +=  d
       else
           a,b = eftAddGTE(a,d)
       end
   end

   a,b
end

function eftAdd{T<:StdFloat}(a::T,b::T,c::T,d::T)

   w=x=y=z=zero(T)

   w,d = eftAdd(c,d)
   w,c = eftAdd(b,w)
   a,b = eftAdd(a,w)

   w,x = a,b
   if x != zero(T)
       x,y = eftAddGTE(x,c)
       if y != zero(T)
           y,z = eftAddGTE(y,d)
       else
           x,y = eftAddGTE(x,d)
       end
   else
       w,x = eftAddGTE(w,c)
       if x != zero(T)
           x,y = eftAddGTE(x,d)
       else
           w,x = eftAddGTE(w,d)
       end
   end

   w,x,y,z
end

function eftAddAs3{T<:StdFloat}(a::T,b::T,c::T,d::T)

   w=x=y=zero(T)

   w,d = eftAdd(c,d)
   w,c = eftAdd(b,w)
   a,b = eftAdd(a,w)

   w,x = a,b
   if x != zero(T)
       x,y = eftAddGTE(x,c)
       if y != zero(T)
           y += d
       else
           x,y = eftAddGTE(x,d)
       end
   else
       w,x = eftAddGTE(w,c)
       if x != zero(T)
           x,y = eftAddGTE(x,d)
       else
           w,x = eftAddGTE(w,d)
       end
   end

   w,x,y
end

function eftAddAs2{T<:StdFloat}(a::T,b::T,c::T,d::T)

   w=x=zero(T)

   w,d = eftAdd(c,d)
   w,c = eftAdd(b,w)
   a,b = eftAdd(a,w)

   w,x = a,b
   if x != zero(T)
       x += c + d
   else
       w,x = eftAddGTE(w,c)
       if x != zero(T)
           x +=  d
       else
           w,x = eftAddGTE(w,d)
       end
   end

   w,x
end