text:
##########################
# FunctionCell Interface #
##########################
"""
	lr(r=0; trans=true, bias=true)

Constructs an untrained cell that when piped data inside (observations and
responses) calculates the regularized (ridge) Least-Squares regression
coefficients of the data based on the responses and regularization parameter
`r`.

The parameter `r` can be a real scalar, a vector or symetric matrix.

# Keyword arguments (same as in `MultivariateStats`)
  * `trans` whether to use the trasposed form, (default `true` e.g. rows are variables, columns are observations)
  * `bias` whether to include a bias term (default `true`)

Read the `MultivariateStats.jl` documentation for more information.
"""
lr(r=0; trans::Bool=true, bias::Bool=true) =
	FunctionCell(lr, (r,), ModelProperties(),
	      kwtitle("Linear regressor", [(:trans,trans),(:bias,bias)]);
	      trans=trans, bias=bias)

############################
# DataCell/Array Interface #
############################
"""
	lr(x,y, [;kwargs])

Trains a Least-Squares regression function cell. Returns the regression coefficients and bias (if available)
"""
# Training
lr(x::T where T<:CellDataU, r::Union{Real, AbstractArray}=0; trans::Bool=true, bias::Bool=true) =
	lr((getx!(x), zeros(nobs(x))), r; trans=trans, bias=bias)
lr(x::T where T<:CellData, r::Union{Real, AbstractArray}=0; trans::Bool=true, bias::Bool=true) =
	lr((getx!(x), gety(x)), r; trans=trans, bias=bias)
lr(x::Tuple{T,S} where T<:AbstractVector where S<:AbstractArray, r::Union{Real, AbstractArray}=0; trans::Bool=true, bias::Bool=true) =
	lr((mat(x[1], LearnBase.ObsDim.Constant{2}()), x[2]), r; trans=trans, bias=bias)
lr(x::Tuple{T,S} where T<:AbstractMatrix where S<:AbstractArray, r::Union{Real, AbstractArray}=0; trans::Bool=true, bias::Bool=true ) = begin

	@assert nobs(x[1]) == nobs(x[2])

	# If y vector, leave as is, if matrix, traspose
	if x[2] isa Vector
		tmp = MultivariateStats.ridge(getobs(x[1]), x[2], r, trans=trans, bias=bias)
	else
		tmp = MultivariateStats.ridge(getobs(x[1]), x[2]', r, trans=trans, bias=bias)
	end

	if bias
		a = tmp[1:end-1,:]
		b = tmp[end,:]
	else
		a = tmp
		b = zeros(nvars(x[2]))
	end

	# Build model properties
	modelprops = ModelProperties(size(a,1), size(a,2))

	# Returned trained cell
	FunctionCell(lr, Model((a,b), modelprops), kwtitle("Linear regressor", [(:trans,trans),(:bias,bias)]))
end

# Execution
lr(x::T where T<:CellData, model::Model{<:Tuple{<:AbstractArray, <:AbstractArray}}) =
	datacell(lr(getx!(x), model), gety(x))
lr(x::T where T<:AbstractVector, model::Model{<:Tuple{<:AbstractArray, <:AbstractArray}}) =
	lr(mat(x, LearnBase.ObsDim.Constant{2}()), model)
lr(x::T where T<:AbstractMatrix, model::Model{<:Tuple{<:AbstractArray, <:AbstractArray}}) =
	model.data[1]'*x .+ model.data[2] # e.g. aᵀX+b


text:
module Utils
   export mycoolfunction

   function mycoolfunction()
      println("this is my cool function")
   end

   function mysecretfunction()
      println("this is my secret function")
   end

end

ext: jl


text:
# ahwtests.jl
# Tests the functions random_chromosome, execute_chromosome, and print_chromosome.
# The parameters of the run can be changed by editing them below.
# There is an optional integer command-line argument.  If present, it is used as
#     the random number seed.  When run with the same seed, the same random
#     chromosome will be constructed.
# With no command line argument, a different random chromosome will be constructed on each run.
# Example run from the "src" subdirectory:  "julia ../test/ahwtests.jl 43"

using CartesianGP
using Base.Test

# These are set in test/Chromosome.jl which is imported above
const numinputs = 3
const numoutputs = 3
const numperlevel = 1
const numlevels = 6
const numlevelsback = 4

if length(ARGS) > 0
    try
        # If the user gives an interger command line argument, use it for a random number seed
        seed = int(ARGS[1])
        srand(seed)
    catch
        # if the int() conversion fails, continue without setting the random number seed
    end
end

funcs = default_funcs()
p = Parameters(numinputs, numoutputs, numperlevel, numlevels, numlevelsback)
c0 = random_chromosome(p, funcs)
result = execute_chromosome(c0)
print_chromosome(c0)    # Print showing all nodes
print_chromosome(c0,true)  # Print showing only active nodes
print("Execution result: [")
for i in 1:numoutputs
   @printf("%#X ",result[i])
end
println("]")


text:
using DiffEqFlux, Flux, OrdinaryDiffEq, ReverseDiff

# Checks for Shapes and Non-Zero Gradients
u0 = rand(Float32, 6, 1)

hnn = HamiltonianNN(Chain(Dense(6, 12, relu), Dense(12, 1)))
p = hnn.p

@test size(hnn(u0)) == (6, 1)

@test ! iszero(ReverseDiff.gradient(p -> sum(hnn(u0, p)), p))

hnn = HamiltonianNN(Chain(Dense(6, 12, relu), Dense(12, 1)))
p = hnn.p

@test size(hnn(u0)) == (6, 1)

@test ! iszero(ReverseDiff.gradient(p -> sum(hnn(u0, p)), p))

# Test Convergence on a toy problem
t = range(0.0f0, 1.0f0, length = 64)
π_32 = Float32(π)
q_t = reshape(sin.(2π_32 * t), 1, :)
p_t = reshape(cos.(2π_32 * t), 1, :)
dqdt = 2π_32 .* p_t
dpdt = -2π_32 .* q_t

data = cat(q_t, p_t, dims = 1)
target = cat(dqdt, dpdt, dims = 1)

hnn = HamiltonianNN(Chain(Dense(2, 16, relu), Dense(16, 1)))
p = hnn.p

opt = ADAM(0.01)
loss(x, y, p) = sum((hnn(x, p) .- y) .^ 2)

initial_loss = loss(data, target, p)

epochs = 100
for epoch in 1:epochs
    gs = ReverseDiff.gradient(p -> loss(data, target, p), p)
    Flux.Optimise.update!(opt, p, gs)
end

final_loss = loss(data, target, p)

@test initial_loss > final_loss

# Test output and gradient of NeuralHamiltonianDE Layer
tspan = (0.0f0, 1.0f0)

model = NeuralHamiltonianDE(
    hnn, tspan, Tsit5(),
    save_everystep = false, save_start = true,
    saveat = range(tspan[1], tspan[2], length=10)
)
sol = Array(model(data[:, 1]))
@test size(sol) == (2, 10)

ps = Flux.params(model)
gs = Flux.gradient(() -> sum(Array(model(data[:, 1]))), ps)

@test ! iszero(gs[model.p])


text:
using RandomizedLinAlg

include("factorization.jl")
include("rlinalg.jl")
include("rsvd.jl")
include("rsvd_fnkz.jl")


text:
"""
Implementation of the streaming peicewise linear regression and hypothesis testing
    framework for partitioning and reconstructing simulation outputs.

Paper: Myers et al (2016). Partitioning a Large Simulation as It Runs.

Author: Wayne Wang
Email: wayneyw@umich.edu
Last modified: 08/02/2021
"""

function stream_lm(;
    p1::Int = 2, p2::Int = 4, buff_size::Int = 5, α::Real = 0.025, δ::Real = 0.001,
    sims_tot::Int = 200)
    """
    Streaming piecewise linear regression with modifed F test for change-point detection
    """
    # initialization
    curr_sims_tot = 0
    partitions, curr, buff, N_curr, lm_fit_buff, suff_stats_buff, lm_fit_curr, suff_stats_curr,
        lm_fit_curr_and_buff, suff_stats_curr_and_buff =
            init(buff_size)

    # compute linear partitions while simulation is running
    while curr_sims_tot <= sims_tot - buff_size - 1
        print(string(curr_sims_tot, "\n"))
        lm_fit_buff, suff_stats_buff, lm_fit_curr, suff_stats_curr,
        lm_fit_curr_and_buff, suff_stats_curr_and_buff =
            update_partitions!(partitions, curr, buff,
                suff_stats_curr, lm_fit_curr,
                suff_stats_buff, lm_fit_buff,
                suff_stats_curr_and_buff, lm_fit_curr_and_buff,
                N_curr, buff_size, p1, p2, α, δ)
    end

    return partitions
end

function init(buff_size::Int)
    """
    Compute initial statistics
    """
    partitions = [] #each element a tuple: (start/end time stamps, sufficient statistics: RSS, beta0, beta1)
    curr_and_buff = simulation_loader(2 * buff_size)
    curr = curr_and_buff[1:buff_size]
    buff = curr_and_buff[buff_size + 1:end]
    lm_fit_curr, suff_stats_curr = update_sufficient_stats(curr)
    curr = [curr[1][2], curr[end][2]] #discard simulation outputs, save only start/end time stamps
    N_curr = buff_size
    lm_fit_buff, suff_stats_buff = update_sufficient_stats(buff)
    lm_fit_curr_and_buff, suff_stats_curr_and_buff = update_sufficient_stats(curr_and_buff) #{curr \cup buff}

    return partitions, curr, buff, N_curr, lm_fit_buff, suff_stats_buff, lm_fit_curr, suff_stats_curr,
        lm_fit_curr_and_buff, suff_stats_curr_and_buff
end

function simulation_loader(N::Int;
    ys::AbstractVector = y, ts::AbstractVector = x)
    """
    Genrate N new simulation outputs
    return a time-series tuples {(y1,t1),…,(yN,tN)}
    """
    sims = []
    for i = curr_sims_tot + 1:curr_sims_tot + N
        push!(sims, (ys[i], ts[i]))
    end
    global curr_sims_tot += N
    return sims
end

function update_partitions!(partitions::AbstractVector,
    curr::AbstractVector,
    buff::AbstractVector,
    suff_stats_curr::AbstractVector,
    lm_fit_curr::AbstractVector,
    suff_stats_buff::AbstractVector,
    lm_fit_buff::AbstractVector,
    suff_stats_curr_and_buff::AbstractVector,
    lm_fit_curr_and_buff::AbstractVector,
    N_curr::Int, buff_size::Int, p1::Int, p2::Int, α::Real, δ::Real)
    """
    Update partitions, curr, and buff data/statistics/fit based on a modified F-test
    """
    # modified F-test
    N = N_curr + buff_size
    RSS_curr = lm_fit_curr[1]
    RSS_buff = lm_fit_buff[1]
    RSS1 = lm_fit_curr_and_buff[1]
    RSS2 = RSS_curr + RSS_buff
    F_stat = ((RSS1 - RSS2) / (p2 - p1)) / ((RSS2 + δ^2 * N) / (N - p2))
    F_crit = quantile(FDist(p2 - p1, N - p2), 1 - α)

    if F_stat > F_crit # reject null
        ## save partition info: start/end time stamps & regressioin fit
        push!(partitions, (copy(curr), copy(lm_fit_curr)))
        curr .= [buff[1][2], buff[end][2]]
        N_curr = buff_size
        buff .= simulation_loader(buff_size)
        ## compute suff. stats
        lm_fit_curr .= copy(lm_fit_buff)
        suff_stats_curr .= copy(suff_stats_buff)
        lm_fit_buff, suff_stats_buff = update_sufficient_stats(buff)
        ## update lm_fit_curr_and_buff for each new point in buff
        lm_fit_curr_and_buff, suff_stats_curr_and_buff =
                update_sufficient_stats(lm_fit_curr, suff_stats_curr, N_curr + buff_size, buff)
    else # fail to reject null
        curr[2] = buff[1][2]
        lm_fit_curr, suff_stats_curr = update_sufficient_stats(lm_fit_curr, suff_stats_curr, N_curr, [buff[1]])
        deleteat!(buff, 1)
        push!(buff, simulation_loader(1)[1])
        lm_fit_buff, suff_stats_buff = update_sufficient_stats(buff)
        lm_fit_curr_and_buff, suff_stats_curr_and_buff =
                update_sufficient_stats(lm_fit_curr_and_buff, suff_stats_curr_and_buff, N_curr + buff_size, [buff[end]])
    end

    return lm_fit_buff, suff_stats_buff, lm_fit_curr, suff_stats_curr,
        lm_fit_curr_and_buff, suff_stats_curr_and_buff
end

function update_sufficient_stats(sims::AbstractVector)
    """
    Compute sufficient statistics & fitted regression line for vectors
        of simulation outputs sims := {(y1,t1),…}
    """
    suff_stats = zeros(5) #sumt, sumt2, sumy, sumy2, summty
    lm_fit = zeros(3) #RSS, intercept, slope
    N = 0

    lm_fit, suff_stats = update_sufficient_stats(lm_fit, suff_stats, N, sims)

    return lm_fit, suff_stats
end

function update_sufficient_stats(lm_fit::AbstractVector, suff_stats::AbstractVector,
    N::Int, new_sims::AbstractVector)
    """
    Update sufficient statistics & fitted regression line for one or more newly
        added simulation outputs
    """
    # create copies to store the new stats
    lm_fit_new = similar(lm_fit) #RSS, intercept, slope
    suff_stats_new = copy(suff_stats) #sumt, sumt2, sumy, sumy2, summty

    for i = 1:length(new_sims)
        newy = new_sims[i][1]
        newt = new_sims[i][2]
        # update sufficient statistics
        suff_stats_new[1] += newt # sum of the time steps
        suff_stats_new[2] += newt * newt # sum of time steps squared
        suff_stats_new[3] += newy # sum of variable of interest (e.g. temp)
        suff_stats_new[4] += newy * newy # sum of variable of interest (e.g. temp)
        suff_stats_new[5] += newt * newy # sum of time cross variable
        N += 1
    end

    lm_fit_new .= compute_lm(N, suff_stats_new)

    return lm_fit_new, suff_stats_new
end

function compute_lm(N::Int, suff_stats::AbstractVector)
    """
    Compute linear regression fit (RSS, slope, intercept)
        using sufficient statistics
    """
    sumt = suff_stats[1]
    sumt2 = suff_stats[2]
    sumy = suff_stats[3]
    sumy2 = suff_stats[4]
    sumty = suff_stats[5]
    RSS = (sumty - sumt * sumy / N)^2
    RSS = RSS / (sumt2 - sumt * sumt / N)
    RSS = sumy2 - sumy * sumy / N - RSS
    #sigma2 = RSS / (N - 2)
    beta1 = sumty - sumt * sumy / N
    beta1 = beta1 / (sumt2 - sumt * sumt / N)
    beta0 = (sumy - beta1 * sumt) / N
    return [RSS, beta0, beta1]
end

max_stars_repo_name: ayanbiswas/PRISM

max_issues_repo_name: ayanbiswas/PRISM

max_forks_repo_name: ayanbiswas/PRISM


text:
struct BlockCOO{Tv,Ti} <: AbstractSparseM{Tv,Ti}
    m::Ti
    n::Ti
    val::Vector{<:AbstractMatrix{Tv}}
    rowind::Vector{Ti}
    colind::Vector{Ti}
end

function BlockCOO(m::Ti, n::Ti, elem::AbstractArray{Tuple{Ti,Ti,Tv},1}) where {Tv,Ti}
    rowind = Vector{Ti}()
    colind = Vector{Ti}()
    val = Vector{Tv}()
    prev_index::Tuple{Ti,Ti} = (0,0)
    for (i,j,u) in sort(elem)
        if m < i
            m = i
        end
        if n < j
            n = j
        end
        if prev_index != (i,j)
            push!(rowind, i)
            push!(colind, j)
            push!(val, u)
            prev_index = (i,j)
        else
            val[end] += u
        end
    end
    BlockCOO(m, n, val, rowind, colind)
end

block(A::SparseCOO{<:AbstractMatrix{Tv},Ti}) where {Tv,Ti} = BlockCOO{Tv,Ti}(A.m, A.n, copy(A.val), copy(A.rowind), copy(A.colind))
block(A::SparseCSR{<:AbstractMatrix{Tv},Ti}) where {Tv,Ti} = block(SparseCOO(A))
block(A::SparseCSC{<:AbstractMatrix{Tv},Ti}) where {Tv,Ti} = block(SparseCOO(A))
block(A::SparseArrays.SparseMatrixCSC{<:AbstractMatrix{Tv},Ti}) where {Tv,Ti} = block(SparseCOO(A))

SparseCSR(A::BlockCOO{Tv,Ti}) where {Tv,Ti} = _tocsr(_tocoo(A))
SparseCSC(A::BlockCOO{Tv,Ti}) where {Tv,Ti} = _tocsc(_tocoo(A))
SparseCOO(A::BlockCOO{Tv,Ti}) where {Tv,Ti} = _tocoo(A)
SparseArrays.sparse(A::BlockCOO{Tv,Ti}) where {Tv,Ti} = SparseArrays.sparse(_tocoo(A))

function Base.show(io::IO, A::BlockCOO{Tv,Ti}) where {Tv,Ti}
    m, n = size(A)
    println(io, string(A.m)*"x"*string(A.n)*" COO-BlockMatrix")
    for z = 1:SparseArrays.nnz(A)
        i = A.rowind[z]
        j = A.colind[z]
        println(io, "($(i), $(j)) $(A.val[z])")
    end
end

function _tocoo(A::BlockCOO{Tv,Ti}) where {Tv,Ti}
    m, n = size(A)
    bi = zeros(Ti, m)
    bj = zeros(Ti, n)
    for z = 1:SparseArrays.nnz(A)
        i = A.rowind[z]
        j = A.colind[z]
        m, n = size(A.val[z])
        if bi[i] == Ti(0)
            bi[i] = m
        else
            @assert bi[i] == m "assert error ($i, $j) ($m, $n)"
        end
        if bj[j] == Ti(0)
            bj[j] = n
        else
            @assert bj[j] == n "assert error ($i, $j) ($m, $n)"
        end
    end
    @assert all(bi .!= Ti(0))
    @assert all(bj .!= Ti(0))

    ## make
    m, n = sum(bi), sum(bj)
    ci = zero(bi)
    ci[1] = Ti(0)
    for i = 1:length(bi)-1
        ci[i+1] = ci[i] + bi[i]
    end
    cj = zero(bj)
    cj[1] = Ti(0)
    for j = 1:length(bj)-1
        cj[j+1] = cj[j] + bj[j]
    end
    val = Tv[]
    rowind = Ti[]
    colind = Ti[]
    for z = 1:SparseArrays.nnz(A)
        i = A.rowind[z]
        j = A.colind[z]
        coo = SparseCOO(A.val[z])
        append!(val, coo.val)
        append!(rowind, coo.rowind .+ ci[i])
        append!(colind, coo.colind .+ cj[j])
    end
    SparseCOO(m, n, val, rowind, colind)
end


text:
# ---
# title: 1397. Find All Good Strings
# id: problem1397
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/find-all-good-strings/description/>
# hidden: true
# ---
#
# Given the strings `s1` and `s2` of size `n`, and the string `evil`. _Return
# the number of **good** strings_.
#
# A **good** string has size `n`, it is alphabetically greater than or equal to
# `s1`, it is alphabetically smaller than or equal to `s2`, and it does not
# contain the string `evil` as a substring. Since the answer can be a huge
# number, return this modulo 10^9 + 7.
#
#
#
# **Example 1:**
#
#
#
#     Input: n = 2, s1 = "aa", s2 = "da", evil = "b"
#     Output: 51
#     Explanation: There are 25 good strings starting with 'a': "aa","ac","ad",...,"az". Then there are 25 good strings starting with 'c': "ca","cc","cd",...,"cz" and finally there is one good string starting with 'd': "da".
#
#
# **Example 2:**
#
#
#
#     Input: n = 8, s1 = "leetcode", s2 = "leetgoes", evil = "leet"
#     Output: 0
#     Explanation: All strings greater than or equal to s1 and smaller than or equal to s2 start with the prefix "leet", therefore, there is not any good string.
#
#
# **Example 3:**
#
#
#
#     Input: n = 2, s1 = "gx", s2 = "gz", evil = "x"
#     Output: 2
#
#
#
#
# **Constraints:**
#
#   * `s1.length == n`
#   * `s2.length == n`
#   * `s1 <= s2`
#   * `1 <= n <= 500`
#   * `1 <= evil.length <= 50`
#   * All strings consist of lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module MeasureBase

const logtwo = log(2.0)

using Random
import Random: rand!

using FillArrays
using ConcreteStructs
using MLStyle

export ≪
export sampletype

export AbstractMeasure

abstract type AbstractMeasure end

import PrettyPrinting

const Pretty = PrettyPrinting

sampletype(μ::AbstractMeasure) = typeof(testvalue(μ))

# sampletype(μ::AbstractMeasure) = sampletype(basemeasure(μ))

export logdensity
export basemeasure
export basekernel

using LogExpFunctions: logsumexp

"""
    logdensity(μ::AbstractMeasure{X}, x::X)

Compute the logdensity of the measure μ at the point x. This is the standard way
to define `logdensity` for a new measure. the base measure is implicit here, and
is understood to be `basemeasure(μ)`.

Methods for computing density relative to other measures will be
"""
function logdensity end

if VERSION < v"1.7.0-beta1.0"
    @eval begin
        struct Returns{T}
            value::T
        end

        (f::Returns)(x) = f.value
    end
end

include("kernel.jl")
include("parameterized.jl")
include("combinators/mapsto.jl")
include("combinators/half.jl")
include("exp.jl")
include("domains.jl")
include("utils.jl")
include("absolutecontinuity.jl")
include("macros.jl")

include("primitive.jl")
include("primitives/counting.jl")
include("primitives/lebesgue.jl")
include("primitives/dirac.jl")
include("primitives/trivial.jl")

include("combinators/factoredbase.jl")
include("combinators/weighted.jl")
include("combinators/superpose.jl")
include("combinators/product.jl")
include("combinators/for.jl")
include("combinators/power.jl")
include("combinators/affine.jl")
include("combinators/spikemixture.jl")
include("combinators/likelihood.jl")
include("combinators/pointwise.jl")
include("combinators/restricted.jl")
include("combinators/smart-constructors.jl")

include("rand.jl")

include("density.jl")

end


text:
using Hygese
using Test
using CVRPLIB, TSPLIB

@testset "Hygese.jl" begin
    # Write your tests here.

    @testset "CVRP" begin
        include("cvrp_test.jl")
    end

    @testset "TSP" begin
        include("tsp_test.jl")
    end

end


text:
for i = 1:100
    msg = "Fizz" ^ (i%3==0) * "Buzz" ^ (i%5==0)
    println(msg=="" ? i : msg)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Random
using Random: default_rng
using .Threads: nthreads

using Distributions, HypothesisTests, StatsBase, ProgressMeter, LoopVectorization,
    ThreadPools, GaitSymmetry, SafeBuffers
using GaitSymmetry: Ratio

export gensimdata, calcsymfun, estimatepower, estimatepower!, powersim

function gensimdata(rng, xdist, ydist, samples, T=promote_type(eltype(xdist), eltype(ydist)))
    stor = Matrix{T}(undef, (samples, 2))
    _gensimdata!(rng, stor, xdist, ydist)
end

function _gensimdata!(rng, stor, xdist, ydist)
    randn!(rng, stor)
    @views _randn_kernel!(stor[:,1], mean(xdist), std(xdist))
    @views _randn_kernel!(stor[:,2], mean(ydist), std(ydist))

    return stor
end

function _randn_kernel!(arr, μ, σ)
    @avx for i in eachindex(arr)
        arr[i] = abs(fma(arr[i], σ, μ))
    end
end

function calcsymfun(symfun, stor::AbstractArray{T}, len=size(stor, 1)) where T
    size(stor, 2) === 2 || throw(DimensionMismatch("expected 2 columns; got $(size(stor, 2))"))
    len ≤ size(stor, 1) || throw(DimensionMismatch("`len` must be ≤ `size(stor, 1)`"))
    res = Vector{T}(undef, len)
    _calcsymfun!(res, symfun, len, stor)
end

function _calcsymfun!(res, symfun, len, stor)
    @inbounds for i in 1:len-1
        x = stor[i,1]
        x′ = stor[i+1,1]
        y = stor[i,2]

        res[i] = symfun(x, x′) - symfun(x, y)
    end
    res[end] = symfun(stor[end,1], stor[1,1]) - symfun(stor[end,1], stor[end,2])

    return res
end

function estimatepower(symfun, xdist, ydist, n, samples)
    stor = gensimdata(default_rng(), xdist, ydist, n*samples)
    res = calcsymfun(symfun, stor)
    _estimatepower(res, n, samples)
end

function estimatepower(symfun, stor, n, samples)
    res = calcsymfun(symfun, stor, n*samples)
    _estimatepower(res, n, samples)
end

function estimatepower!(res, symfun, stor, n, samples)
    fullN = n*samples
    _calcsymfun!(res, symfun, fullN, stor)
    _estimatepower(res, n, samples)
end

function _estimatepower(res, n, samples)
    numsig = _numsignificant(res, n, samples)
    return numsig# / samples
end

function _numsignificant(res, n, samples, α=0.05)
    numsig = 0
    for i in 0:(samples-1)
        @inbounds slice = @view(res[(1:n).+i])
        avg, sd = mean_and_std(slice; corrected=true)
        if pvalue(OneSampleTTest(avg, sd, n)) ≤ α
            numsig += 1
        end
    end

    return numsig
end

function powersim(
    metrics,
    ratios,
    vars,
    Ns,
    samples,
    U=Float32;
    null=false,
    batch=1
)
    p = Progress(length(metrics)*length(ratios)*length(vars)*sum(Ns .* samples)+1;
        dt=5, desc="Batch $(lpad(batch, 2)) progress: ")

    ratios = convert(Array{U}, ratios)
    data = Array{U}(undef, (length(Ns), length(metrics), length(ratios), length(vars)))

    fullN = maximum(Ns)*samples
    respool = BufferPool(nthreads()+1, () -> Vector{U}(undef, fullN))
    storpool = BufferPool(nthreads()+1, () -> Matrix{U}(undef, (fullN, 2)))

    idxs = [ (σi, i) for σi in eachindex(vars) for i in eachindex(ratios) ]
    @qthreads for (σi, i) in idxs
        σ = vars[σi]
        xdist = Normal(1, σ)
        μ = ratios[i]
        μ_sd = μ*σ
        ydist = ifelse(null, xdist, Normal(μ, μ_sd))

        withbuffer(storpool) do stor
            _gensimdata!(default_rng(), stor, xdist, ydist)
            for mi in eachindex(metrics)
                if metrics[mi] === Alv20b
                    metric = Alv20b(U(σ))
                else
                    metric = metrics[mi]
                end

                withbuffer(respool) do res
                    _calcsymfun!(res, metric, fullN, stor)
                    for ni in eachindex(Ns)
                        data[ni, mi, i, σi] = _estimatepower(res, Ns[ni], samples)
                        next!(p; step=Ns[ni]*samples)
                    end
                end
            end
        end
    end
    finish!(p)

    return data
end


text:
module TestEGS
using Test
using PhaseSpaceIO
using PhaseSpaceIO.Testing
using PhaseSpaceIO: ptype, EGSHeader, write_header, write_particle, EGSParticle
using Accessors

@testset "test egs_iterator egs_writer" begin
    for _ in 1:100
        ZSLAB = rand(Bool) ? Nothing : Float32
        n = rand(1:1000)
        P = EGSParticle{ZSLAB}
        ps = P[arbitrary(P) for _ in 1:n]
        path = tempname() * ".egsphsp1"
        egs_writer(path,P) do w
            for p in ps
                write(w, p)
            end
        end
        @test ispath(path)
        ps_reload = egs_iterator(collect, path)
        ps_reload2 = phsp_iterator(collect, path)
        @test ps == ps_reload
        @test ps_reload2 == ps_reload
        @test eltype(ps) == eltype(ps_reload)
        rm(path)
    end
end

@testset "simple phase spaces" begin

    path = testdatapath("photon_electron_positron.egsphsp")
    h_truth = PhaseSpaceIO.EGSHeader{EGSParticle{Float32}}(3, 1, 1.0f0, 0.83312154f0, 1.0f0)
    ps_truth = [
                EGSParticle(
                    E=1.0, weight=5.2134085, x=2.0, y=-3.0,
                    u=0.99884456, v=-0.03811472, w=0.029271236,
                    new_history=false, zlast=0.23947382f0, latch=Latch(charge=0))
                EGSParticle(
                    E=0.8331216, weight=9.296892, x=-0.41874805, y=-0.04104328,
                    u=0.18903294, v=-0.29983893, w=-0.93507385,
                    new_history=true, zlast=1.3316472f0, latch=Latch(charge=-1))
                EGSParticle(
                    E=0.6077692, weight=0.8328271, x=0.922653, y=0.9263571,
                    u=0.22435355, v=-0.41052637, w=-0.8838176,
                    new_history=false, zlast=1.0724568f0, latch=Latch(charge=1))
    ]
    h_loaded, ps_loaded = phsp_iterator(path) do iter
        iter.header, collect(iter)
    end
    @test ps_truth[1] == ps_loaded[1]
    @test ps_truth[2] == ps_loaded[2]
    @test ps_truth[3] == ps_loaded[3]
    @test ps_truth == ps_loaded
    @test h_truth == h_loaded
end

@testset "finalizer" begin
    function write_sloppy(path, ps)
        P = eltype(ps)
        w = egs_writer(path, P)
        for p in ps
            write(w, p)
        end
    end
    path = tempname() * ".egsphsp"
    ps = [arbitrary(EGSParticle{Nothing})]
    write_sloppy(path, ps)
    GC.gc()
    ps2 = phsp_iterator(collect, path)
    @test length(ps) == length(ps2) == 1
    @test first(ps) ≈ first(ps2)
    rm(path)
end

@testset "Latch" begin
    l0 = Latch(zero(UInt32))

    # binary represenation taken from beamdp
    l = @set l0.creation = 5
    @test 0b00000101000000000000000000000000 === UInt32(l)

    l = @set l0.brems = true
    @test 0b00000000000000000000000000000001 === UInt32(l)

    l = @set l0.visited = @BitRegions(2,4,7,21)
    @test 0b00000000001000000000000010010100 === UInt32(l)

    l = @set l.creation = 3
    l = @set l.visited = @BitRegions(1,3,23)
    @test 0b00000011100000000000000000001010 === UInt32(l)
end

@testset "particle_type" begin
    p = EGSParticle(latch=Latch(multicross=false, charge=0, creation=0, visited=@BitRegions(), brems=false, ), new_history=true, E=0.74988735f0, x=-0.51910776f0, y=0.22886558f0, u=0.26358783f0, v=-0.9477993f0, w=-0.17943771f0, weight=4.7387686f0, zlast=nothing, )
    @test particle_type(p) === photon
    @test isphoton(p)
    p = @set p.latch.charge = -1
    @test iselectron(p)
    @test particle_type(p) === electron
    p = @set p.latch.charge = 1
    @test particle_type(p) === positron
    @test ispositron(p)
end

end #module


text:
export plot_simulation!

# convenience function for plotting simulation results
# use `output_map` to plot a linear combination of the state variables
function plot_simulation!(fig, sim::EnsembleSimulationSolution; vars=nothing, output_map=nothing, kwargs...)
    # The main problem is that plotting trajectories one by one changes the plot
    # limits. Hence we store the plot limits from an existing figure and restore
    # them after plotting all trajectories.

    # argument checking
    got_vars = !isnothing(vars)
    got_output_map = !isnothing(output_map)
    color = get(kwargs, :color, :red)
    label = get(kwargs, :lab, "")

    if !got_vars && !got_output_map
        throw(ArgumentError("either `vars` or `output_map` should be specified"))
    elseif got_vars
        if vars isa AbstractVector
            if length(vars) != 2
                error("unexpected length of `vars`: $(length(vars))")
            end
            vars = (vars[1], vars[2])
        end
        _plot_function = _plot_simulation_vars!
        opts = vars
    else
        _plot_function = _plot_simulation_output_map!
        opts = output_map
    end

    # obtain x and y limits
    xl = Plots.xlims(fig)
    yl = Plots.ylims(fig)

    _plot_function(fig, sim, opts, color=color, label=label)

    # restore x and y limits
    Plots.xlims!(fig, xl)
    Plots.ylims!(fig, yl)

    return fig
end

function _plot_simulation_vars!(fig, sim, vars; color, label)
    for simulation in trajectories(sim)
        for piece in simulation
            Plots.plot!(fig, piece, vars=vars, color=color, lab=label)
            label = ""  # overwrite to have exactly one label
        end
    end
end

function _plot_simulation_output_map!(fig, sim, output_map::Vector{<:Real}; color, label)

    # dimension check
    numvars = length(output_map)
    traj1 = sim.solutions[1].trajectory[1]
    n = length(traj1.u[1])
    if numvars == n
        c0 = 0.0
        coeffs = output_map
    elseif numvars == n + 1
        c0 = output_map[1]
        coeffs = output_map[2:end]
    else
        throw(ArgumentError("the length of the `output_map` should be $n or $(n+1), got $numvars"))
    end

    f(t, P) = t -> c0 + sum(coeffs[i] * P(t)[i] for i in 1:n)
    isfirst = true

    for simulation in trajectories(sim)
        for piece in simulation
            dt = piece.t
            trange = range(dt[1], dt[end], length=length(dt))
            if isfirst
                # plot first point only for the legend entry
                trange1 = trange[1]:trange[1]
                Plots.plot!(fig, trange1, f.(trange1, Ref(piece)), color=color, lab=label)
                label = ""  # overwrite to have exactly one label
                isfirst = false
            end
            Plots.plot!(fig, trange, f.(trange, Ref(piece)), color=color, lab="")
        end
    end
    return fig
end


text:
"""
    2D: global_frame(w::AbstractArray{<:Real,1}, cosa, sina)
    3D: global_frame(w::AbstractArray{<:Real,1}, dirccos::AbstractArray{<:Real,2})

Transform local flow variables to global frame
"""
function global_frame(w::T, cosa, sina) where {T<:AbstractArray{<:Real,1}}

    if eltype(w) <: Int
        G = similar(w, Float64)
    else
        G = similar(w)
    end

    if length(w) == 2
        G[1] = w[1] * cosa - w[2] * sina
        G[2] = w[1] * sina + w[2] * cosa
    elseif length(w) == 4
        G[1] = w[1]
        G[2] = w[2] * cosa - w[3] * sina
        G[3] = w[2] * sina + w[3] * cosa
        G[4] = w[4]
    else
        throw("local -> global: dimension dismatch")
    end

    return G

end

function global_frame(
    w::T,
    dirccos::X,
) where {T<:AbstractArray{<:Real,1},X<:AbstractArray{<:Real,2}}

    if eltype(w) <: Int
        G = similar(w, Float64)
    else
        G = similar(w)
    end

    if length(w) == 3
        G[1] = w[1] * dirccos[1, 1] + w[2] * dirccos[2, 1] + w[3] * dirccos[3, 1]
        G[2] = w[1] * dirccos[1, 2] + w[2] * dirccos[2, 2] + w[3] * dirccos[3, 2]
        G[3] = w[1] * dirccos[1, 3] + w[2] * dirccos[2, 3] + w[3] * dirccos[3, 3]
    elseif length(w) == 5
        G[1] = w[1]
        G[2] = w[2] * dirccos[1, 1] + w[3] * dirccos[2, 1] + w[4] * dirccos[3, 1]
        G[3] = w[2] * dirccos[1, 2] + w[3] * dirccos[2, 2] + w[4] * dirccos[3, 2]
        G[4] = w[2] * dirccos[1, 3] + w[3] * dirccos[2, 3] + w[4] * dirccos[3, 3]
        G[5] = w[5]
    else
        throw("local -> global: dimension dismatch")
    end

    return G

end

"""
    2D: local_frame(w::AbstractArray{<:Real,1}, cosa, sina)
    3D: local_frame(w::AbstractArray{<:Real,1}, dirccos::AbstractArray{<:Real,2})

Transform global flow variables to local frame
"""
function local_frame(w::T, cosa, sina) where {T<:AbstractArray{<:Real,1}}

    if eltype(w) <: Int
        L = similar(w, Float64)
    else
        L = similar(w)
    end

    if length(w) == 2
        L[1] = w[1] * cosa + w[2] * sina
        L[2] = w[2] * cosa - w[1] * sina
    elseif length(w) == 4
        L[1] = w[1]
        L[2] = w[2] * cosa + w[3] * sina
        L[3] = w[3] * cosa - w[2] * sina
        L[4] = w[4]
    else
        throw("global -> local: dimension dismatch")
    end

    return L

end

function local_frame(
    w::T,
    dirccos::X,
) where {T<:AbstractArray{<:Real,1},X<:AbstractArray{<:Real,2}}

    if eltype(w) <: Int
        L = similar(w, Float64)
    else
        L = similar(w)
    end

    if length(w) == 3
        L[1] = w[1] * dirccos[1, 1] + w[2] * dirccos[1, 2] + w[3] * dirccos[1, 3]
        L[2] = w[1] * dirccos[2, 1] + w[2] * dirccos[2, 2] + w[3] * dirccos[2, 3]
        L[3] = w[1] * dirccos[3, 1] + w[2] * dirccos[3, 2] + w[3] * dirccos[3, 3]
    elseif length(w) == 5
        L[1] = w[1]
        L[2] = w[2] * dirccos[1, 1] + w[3] * dirccos[1, 2] + w[4] * dirccos[1, 3]
        L[3] = w[2] * dirccos[2, 1] + w[3] * dirccos[2, 2] + w[4] * dirccos[2, 3]
        L[4] = w[2] * dirccos[3, 1] + w[3] * dirccos[3, 2] + w[4] * dirccos[3, 3]
        L[5] = w[5]
    else
        throw("global -> local: dimension dismatch")
    end

    return L

end

"""
    2D: unit_normal(p1::T, p2::T) where {T<:AbstractVector}
    3D: unit_normal(p1::T, p2::T, p3::T) where {T<:AbstractVector}

Calculate unit normal vector
"""
function unit_normal(p1::T, p2::T) where {T<:AbstractVector}
    Δ = p2 .- p1
    l = norm(Δ) + 1e-6

    return [-Δ[2], Δ[1]] ./ l
end

function unit_normal(p1::T, p2::T, p3::T) where {T<:AbstractVector}
    v1 = p2 .- p1
    v2 = p3 .- p1

    n = cross(v1, v2)
    l = norm(n) + 1e-6

    return n ./ l
end

"""
    point_distance(p1::T, p2::T) where {T<:AbstractVector}
    point_distance(p::T, p1::T, p2::T) where {T<:AbstractVector}

Calculate point-point/line/surface distance
"""
point_distance(p1::T, p2::T) where {T<:AbstractVector} = norm(p1 .- p2)

function point_distance(p::T, p1::T, p2::T) where {T<:AbstractVector}
    x0, y0 = p
    x1, y1 = p1
    x2, y2 = p2

    return abs((x2 - x1) * (y1 - y0) - (x1 - x0) * (y2 - y1)) /
           sqrt((x2 - x1)^2 + (y2 - y1)^2 + 1e-6)
end


text:
# (Ridge) regularization

############################################################################
### Types
############################################################################

struct SemRidge{P, W1, W2, FT, GT, HT} <: SemLossFunction
    α::P
    which::W1
    which_H::W2

    objective::FT
    gradient::GT
    hessian::HT
end

############################################################################
### Constructors
############################################################################

function SemRidge(;α_ridge, which_ridge, n_par, parameter_type = Float64, imply = nothing, kwargs...)
    if which_ridge isa Vector{Symbol}
        which_ridge = get_identifier_indices(which_ridge, imply)
    end
    which = [CartesianIndex(x) for x in which_ridge]
    which_H = [CartesianIndex(x, x) for x in which_ridge]
    return SemRidge(
        α_ridge,
        which,
        which_H,

        zeros(parameter_type, 1),
        zeros(parameter_type, n_par),
        zeros(parameter_type, n_par, n_par))
end

############################################################################
### functors
############################################################################

function (ridge::SemRidge)(par, F, G, H, model)

    if G
        ridge.gradient[ridge.which] .= 2*ridge.α*par[ridge.which]
    end

    if H
        @views @. ridge.hessian[ridge.which_H] += ridge.α*2.0
    end

    if F
        ridge.objective[1] = ridge.α*sum(par[ridge.which].^2)
    end

end

############################################################################
### Recommended methods
############################################################################

objective(lossfun::SemRidge) = lossfun.objective
gradient(lossfun::SemRidge) = lossfun.gradient
hessian(lossfun::SemRidge) = lossfun.hessian

update_observed(loss::SemRidge, observed::SemObs; kwargs...) = loss

############################################################################
### Pretty Printing
############################################################################

function Base.show(io::IO, struct_inst::SemRidge)
    print_type_name(io, struct_inst)
    print_field_types(io, struct_inst)
end


text:
# ---
# title: 1128. Number of Equivalent Domino Pairs
# id: problem1128
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/number-of-equivalent-domino-pairs/description/>
# hidden: true
# ---
#
# Given a list of `dominoes`, `dominoes[i] = [a, b]` is _equivalent_ to
# `dominoes[j] = [c, d]` if and only if either (`a==c` and `b==d`), or (`a==d`
# and `b==c`) - that is, one domino can be rotated to be equal to another
# domino.
#
# Return the number of pairs `(i, j)` for which `0 <= i < j < dominoes.length`,
# and `dominoes[i]` is equivalent to `dominoes[j]`.
#
#
#
# **Example 1:**
#
#
#
#     Input: dominoes = [[1,2],[2,1],[3,4],[5,6]]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= dominoes.length <= 40000`
#   * `1 <= dominoes[i][j] <= 9`
#
#
## @lc code=start
using LeetCode

function num_equiv_domino_pairs(dominoes::Vector{Tuple{Int,Int}})::Int
    counts = Dict{Tuple{Int,Int},Int}()
    n = 0
    for (k, v) in dominoes
        p = minmax(k, v)
        if haskey(counts, p)
            n += counts[p]
            counts[p] += 1
        else
            counts[p] = 1
        end
    end
    return n
end
## @lc code=end


text:
#   Copyright 2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
Example Julia source file in YASStyle format.
"""

# --- Exports

export factorial, add_one

# --- Imports

using Pkg

# --- Types

struct A
    field1::Any
end

# --- Methods

"""
    factorial(n)

Return ``n!``
"""
function factorial(n)
    result = 1
    for i in 1:n
        result *= i
    end

    return result
end

"""
    add_one(x)

Return ``x+1``.
"""
add_one(x) = x + 1

"""
    print_hello()

Print "Hello, World!".
"""
function print_hello()
    println("Hello, World!")

    return nothing
end


text:
# https://www.cs.cornell.edu/gries/Logic/Axioms.html
# The axioms of calculational propositional logic C are listed in the order in
# which they are usually presented and taught. Note that equivalence comes
# first. Note also that, after the first axiom, we take advantage of
# associativity of equivalence and write sequences of equivalences without
# parentheses. We use == for equivalence, | for disjunction, & for conjunction,
# ~ for negation (not), => for implication, and <= for consequence.
#
# Associativity of ==: ((p == q) == r) == (p == (q == r))
# Symmetry of ==: p == q == q == p
# Identity of ==: true == q == q
#
# Definition of false: false == ~true
# Distributivity of not: ~(p == q) == ~p == q
# Definition of =/=: (p =/= q) == ~(p == q)
#
# Associativity of |: (p | q) & r == p | (q | r)
# Symmetry of |: p | q == q | p
# Idempotency of |: p | p == p
# Distributivity of |: p |(q == r) == p | q == p | r
# Excluded Middle: p | ~p
#
# Golden rule: p & q == p == q == p | q
#
# Implication: p => q == p | q == q
# Consequence: p <= q == q => p

# Metatheory.options[:printiter] = true

using Metatheory

calc = @theory begin
    ((p == q) == r)     ==  (p == (q == r))
    (p == q)            ==  (q == p)
    (q == q)            =>  true

    ¬(p == q)           ==  (¬(p) == q)
    (p != q)            ==  ¬(p == q)

    ((p ∨ q) ∨ r)       ==  (p ∨ (q ∨ r))
    (p ∨ q)             ==  (q ∨ p)
    (p ∨ p)             =>  p
    (p ∨ (q == r))      ==  (p ∨ q == p ∨ r)
    (p ∨ ¬(p))          =>  true

    # DeMorgan
    ¬(p ∨ q)            ==  (¬p ∧ ¬q)
    ¬(p ∧ q)            ==  (¬p ∨ ¬q)

    (p ∧ q)             ==  ((p == q) == p ∨ q)

    (p => q)            ==  ((p ∨ q) == q)
    # (p => q)            ==  (¬p ∨ q)
    # (p <= q)            =>  (q => p)
end

fold = @theory begin
    (p::Bool == q::Bool)    |>     (p == q)
    (p::Bool ∨ q::Bool)     |>     (p || q)
    (p::Bool => q::Bool)    |>     ((p || q) == q)
    (p::Bool ∧ q::Bool)     |>     (p && q)
    ¬(p::Bool)              |>     (!p)
end

# t = or_alg ∪ and_alg ∪ neg_alg ∪ demorgan ∪ and_or_distrib ∪
#     absorption ∪ calc

t = calc ∪ fold

g = EGraph(:(((¬p == p) == false)))
saturate!(g, t)
extract!(g, astsize)

@test @areequal t true ((¬p == p) == false)
@test @areequal t true ((¬p == ¬p) == true)
@test @areequal t true ((¬p ∨ ¬p) == ¬p) (¬p ∨ p) ¬(¬p ∧ p)
@test @areequal t true ((p => (p ∨ p)) == true)
params = SaturationParams(timeout=12, eclasslimit=10000, schedulerparams=(1000, 5))

@test areequal(t, true, :(((p => (p ∨ p)) == ((¬(p) ∧ q) => q)) == true); params=params)

# Frege's theorem
# params = SaturationParams(timeout=12, eclasslimit=15000, scheduler=Schedulers.ScoredScheduler)
# params = SaturationParams(timeout=12, eclasslimit=15000, schedulerparams=(500, 2))
@test_skip areequal(t, true, :((p => (q => r)) => ((p => q) => (p => r))); params=params)

# Demorgan's
@test @areequal t true (¬(p ∨ q) == (¬p ∧ ¬q))

# Consensus theorem
# @test_skip
areequal(t, :((x ∧ y) ∨ (¬x ∧ z) ∨ (y ∧ z)), :((x ∧ y) ∨ (¬x ∧ z)); params=params)


text:
################################################################################
# Structure to keep track of performance information
################################################################################
"""
State tracking struct for Simulated Annealing placement.
"""
mutable struct SAState
    #= Information related to most recent move attempt =#
    "Current Temperature of the system."
    temperature::Float64
    "Current objective value."
    objective::Float64
    "Current distance limit."
    distance_limit::Float64
    "Current distance limit as an integer."
    distance_limit_int::Int64
    "Maximum distance limit for this architecture."
    max_distance_limit::Float64

    recent_move_attempts::Int64
    recent_successful_moves::Int64
    recent_accepted_moves::Int64
    recent_deviation::Float64

    "Current warming state"
    warming::Bool

    #= Global state information =#
    "Total number of move attempts made"
    total_moves::Int64
    "Number of successful move generations"
    successful_moves::Int64
    "Number of moves accepted"
    accepted_moves::Int64
    "Moves per second"
    moves_per_second::Float64
    "Average Deviaation"
    deviation::Float64
    "Creation time of the structure"
    start_time::Float64
    "Running Time - only an approximate measure"
    run_time::Float64
    #= Methods for dealing with the updates =#
    display_updates::Bool
    "Time of the last update"
    last_update_time::Float64
    "Update Interval"
    dt::Float64

    #= Auxiliary cost =#
    "Auxiliary cost"
    aux_cost::Float64

    function SAState(temperature, distance_limit, objective)
        # This is a hack to "cleanly" determine if we should display updates
        # or not.
        display_updates = true
        # if MicroLogging.Info >= MicroLogging._min_enabled_level[]
        #     display_updates = true
        # end

        SAStateUpdateInterval = 5
        return new(
            # Most recent run
            temperature,        # temperature
            Float64(objective), # objective
            distance_limit,     # distance_limit
            floor(Int, distance_limit), # distance_limit_int
            distance_limit,     # max_distance_limit (set to initial)
            0,                  # recent_move_attempts
            0,                  # recent_successful_moves
            0,                  # recent_accepted_moves
            0.0,                # recent_deviation
            true,               # warming

            # Global run
            0,      # total_moves
            0,      # suffessful_moves
            0,      # accepted_moves
            0.0,    # moves_per_second
            0.0,    # deviation
            time(), # creation time
            time(), # runtime

            # Update parameters
            display_updates,
            time(),                     # Time of creation
            SAStateUpdateInterval,      # Default update interval

            # auxiliary cost
            0.0,
        )
    end
end

function reset!(s::SAState)
    s.recent_move_attempts = 0
    s.recent_successful_moves = 0
    s.recent_accepted_moves = 0
    s.recent_deviation = 0.0

    return s.warming = true
end

function update!(state::SAState)
    # Update the global counters
    state.total_moves += state.recent_move_attempts
    state.successful_moves += state.recent_successful_moves
    state.accepted_moves += state.recent_accepted_moves
    # Compute number of moves per second.
    state.moves_per_second = state.successful_moves / (time() - state.start_time)
    # Update the exponential moving average of the objective function differnce.
    alpha = 0.2
    state.deviation = alpha * state.recent_deviation + (1.0 - alpha) * state.deviation

    # Determine whether or not to print out results
    current_time = time()
    state.run_time = current_time - state.start_time
    if current_time > state.last_update_time + state.dt
        show_stats(state)
        state.last_update_time = time()
    end

    # Compute the new distance limit integer and return "true" if it changed.
    new_distance_limit_int = max(1, round(Int, state.distance_limit))

    return_value = false
    if new_distance_limit_int != state.distance_limit_int
        return_value = true
    end
    state.distance_limit_int = new_distance_limit_int
    return return_value
end

function show_stats(state::SAState, first = false)
    state.display_updates || return nothing
    #=
    Fields of the SAState will be printed to the console in the order the
    order they show up here.

    The field name will be printed at the beginning of an placement iteration.
    Under-scores will be removed and the first letter of each word will
    be capitalized.

    The space occupied by the field name and the subsequent printed values
    will scale depending on the length of the string conversion of the
    field name. Thus, if the fields are printed out in a different order,
    everything should still look good.
    =#
    fields = (
        :temperature,
        :objective,
        #:total_moves,
        #:successful_moves,
        :accepted_moves,
        :moves_per_second,
        :deviation,
        :run_time,
        :distance_limit,
        #:aux_cost,
    )

    #=
    Key-word arguments for the "format" function. Add to this to fine tune
    how the result will be printed. Don't include the :width keyword argument
    as that will be automatically determined based on the width needed to
    print the field name.

    The "format" function belongs to the "Formatting" package.
    =#
    kwargs = Dict(
        :temperature => Dict(:precision => 5),
        :objective => Dict(:precision => 2),
        :moves_per_second => Dict(:precision => 3, :autoscale => :metric),
        :deviation => Dict(:precision => 3, :autoscale => :metric),
        :run_time => Dict(:precision => 2),
    )

    #=
    Spacing between words and fields. The white space allows easier
    differentiation for each field.
    =#
    spacing = 4

    #=
    If this is the first iteration, just print each field name and underline
    it with a bunch of "-".
    =#
    if first
        # Record the total number of characters written to the console.
        total_length = 0
        # Iterate through each field in order
        for f in fields
            # Replace the under scores in the field name with a space.
            string_f = replace(string(f), "_" => " ")
            # Get the length of the string and add in the spacing between
            # strings
            padded_length = length(string_f) + spacing
            # Pad the string to the desired length and capitalize the first
            # word of each sentence.
            string_to_print = titlecase(lpad(string_f, padded_length))
            # Print to console
            printstyled(string_to_print; bold = true, color = :light_green)
            # Record the number of characters written.
            total_length += padded_length
        end
        # Create a new line and print hyphens for the number of characters written.
        println("\n", "-"^total_length)

        #=
        Print out the numerical values for the requested fields.
        =#
    else
        # Iterate through fields in roder.
        for f in fields
            # Get the length of the field - append one space to the front.
            field_length = length(string(f)) + spacing
            # Get the format preferences dictionary
            kwarg = get(kwargs, f, Dict())
            # Check if we have to shorten the length - adjust because the "width'
            # parameter apparently doesn't take this into account
            haskey(kwarg, :autoscale) && (field_length -= 1)
            # Add the width field to the dictionary.
            kwarg[:width] = field_length
            # Print out the value of the field with the given options. Use
            # the function "format" from the "Formatting" package to make
            # things look nice.
            print(format(getfield(state, f); kwarg...))
        end
        # Terminate with a new line.
        print("\n")
    end
    return nothing
end


text:
export bezier, dbezier, bern_basis

"""
	bezier(npts, b[], dpts)

Calculate a Bezier curve over the `npts` control polygon points in `b`.

The method evaluates `dpts` points of the curve and gives back the tuple `(P, EV)`.
Here `P` is the vector of points memorized by columns while `EV` is the `1`-dimensional cellular complex.
This structure as it is formed could be directly imported in Plasm module by typing:
```julia
julia> P, EV = bezier(npts, b[], dpts);
julia> Plasm.view(P, EV)
```

---

# Arguments

- `npts::Int64`: the number of the control polygon vertices.
- `b::Array{Float64, 2}`: 2-dimensional Array containing the control polygon vertices in the x, y, z coordinates for columns.
- `dpts::Int64`: number of data points to be calculated on the curve.

---

# Examples

```jldoctest
julia> bezier(4,[1.0 2 4 3; 1 3 3 1; 0 0 0 0], 7)[1]
3×7 Array{Float64,2}:
 1.0  1.56481  2.18519  2.75  3.14815  3.26852  3.0
 1.0  1.83333  2.33333  2.5   2.33333  1.83333  1.0
 0.0  0.0      0.0      0.0   0.0      0.0      0.0
```

```jldoctest
julia> bezier(4,[1.0 2 4 3; 1 3 3 1; 0 0 0 0], 7)[2]
6-element Array{Array{Int64,1},1}:
 [1, 2]
 [2, 3]
 [3, 4]
 [4, 5]
 [5, 6]
 [6, 7]
```

---

_By Elia Onofri_
"""

function bezier(npts::Int64, b::Array{Float64,2}, dpts::Int64)::Tuple{Array{Float64,2}, Array{Array{Int64,1},1}}

	@assert npts == length(b[1,:]) ("ERROR: there are not npts = $(npts) points of the polygon in b[]")

	step::Float64 = 1/(dpts-1)
	n::Int64 = npts-1

	C::Array{Float64,2} = zeros(npts, npts)
	for i = 0:n
		for j = 0:n-i
			if (n-j-i)%2==0
				C[i+1, j+1] = binomial(n-j, n-j-i)
			else
				C[i+1, j+1] = 0.0-binomial(n-j, n-j-i)
			end
		end
	end

	D::Array{Float64,2} = zeros(npts, npts)
	for i = 1:npts
		D[i,i] = abs(C[npts+1-i,1])
	end

	G::Array{Float64,2} = transpose(b)

	T::Array{Float64,2} = ones(dpts, npts)

	t::Float64 = 0.0
	for i = 1:dpts
		for j = 0:n-1  #Per j=0 il valore è 1 ed è già impostato
			T[i, n-j]=T[i, npts-j]*t
		end
		t = t + step
	end

	EV::Array{Array{Int64,1},1} = Array{Int64,1}[] #1-dimensional cellular complex used for plotting the curve with Plasm package

	for i = 1 : dpts-1
		app = Array{Int64}(2)
		app[1] = i
		app[2] = i+1
		push!(EV,app)
	end

	P::Array{Float64,2} = transpose(T*C*D*G)

	return (P, EV)
end

#-----------------------------------------------------------------------

"""
	dbezier(npts, b[], dpts)

Calculate Bezier curve and its first and second derivatives over the `npts` control polygon points in `b`.

The method evaluates `dpts` points of the curve and gives back the tuple `(P, D1, D2, EV, DEV)`.

Here `P` is the vector of points memorized by columns while `EV` is the `1`-dimensional cellular complex associated to it.
Moreover `D1` and `D2` are the derivatives vectors and `DEV` is the `1`-dimensional cellular complex associated.
Please note that only `dpts-2` points are evaluated for the derivatives as there are not enought informations about the first and the last points.

This structure as it is formed could be directly imported in Plasm module by typing:
```julia
julia> P, _, _, EV, _ = bezier(npts, b[], dpts);
julia> Plasm.view(P, EV)
```

---

# Arguments

- `npts::Int64`: the number of the control polygon vertices.
- `b::Array{Float64, 2}`: 2-dimensional Array containing the control polygon vertices in the x, y, z coordinates for columns.
- `dpts::Int64`: number of data points to be calculated on the curve.

---

# Examples

```jldoctest
julia> dbezier(4,[1.0 2 4 3; 1 3 3 1; 0 0 0 0], 7)[1]
3×7 Array{Float64,2}:
 1.0  1.56481  2.18519  2.75  3.14815  3.26852  3.0
 1.0  1.83333  2.33333  2.5   2.33333  1.83333  1.0
 0.0  0.0      0.0      0.0   0.0      0.0      0.0
```

```jldoctest
julia> dbezier(4,[1.0 2 4 3; 1 3 3 1; 0 0 0 0], 7)[2]
3×5 Array{Float64,2}:
 3.66667  3.66667  3.0   1.66667  3.26852
 4.0      2.0      0.0  -2.0      1.83333
 0.0      0.0      0.0   0.0      0.0
```

```jldoctest
julia> dbezier(4,[1.0 2 4 3; 1 3 3 1; 0 0 0 0], 7)[3]
3×5 Array{Float64,2}:
   2.0   -2.0   -6.0  -10.0  3.26852
 -12.0  -12.0  -12.0  -12.0  1.83333
   0.0    0.0    0.0    0.0  0.0
```

```jldoctest
julia> dbezier(4,[1.0 2 4 3; 1 3 3 1; 0 0 0 0], 7)[4]
6-element Array{Array{Int64,1},1}:
 [1, 2]
 [2, 3]
 [3, 4]
 [4, 5]
 [5, 6]
 [6, 7]
```

```jldoctest
julia> dbezier(4,[1.0 2 4 3; 1 3 3 1; 0 0 0 0], 7)[5]
6-element Array{Array{Int64,1},1}:
 [1, 2]
 [2, 3]
 [3, 4]
 [4, 5]
```

---

_By Elia Onofri_
"""

function dbezier(npts::Int64, b::Array{Float64,2}, dpts::Int64)::Tuple{Array{Float64,2}, Array{Float64,2}, Array{Float64,2}, Array{Array{Int64,1},1}, Array{Array{Int64,1},1}}

	@assert npts == length(b[1,:]) ("ERROR: there are not npts = $(npts) points of the polygon in b[]")

	step::Float64 = 1/(dpts-1)
	n::Int64 = npts-1

	C::Array{Float64,2} = zeros(npts, npts)
	for i = 0:n
		for j = 0:n-i
			if (n-j-i)%2==0
				C[i+1, j+1] = binomial(n-j, n-j-i)
			else
				C[i+1, j+1] = 0.0-binomial(n-j, n-j-i)
			end
		end
	end

	D::Array{Float64,2} = zeros(npts, npts)
	for i = 1:npts
		D[i,i] = abs(C[npts+1-i,1])
	end

	G::Array{Float64,2} = transpose(b)

	T::Array{Float64,2} = ones(dpts, npts)

	t::Float64 = 0.0
	for i = 1:dpts
		for j = 0:n-1  #Per j=0 il valore è 1 ed è già impostato
			T[i, n-j] = T[i, npts-j]*t
		end
		t = t + step
	end

	Der1::Array{Float64,2} = ones(dpts-2, npts)
	Der2::Array{Float64,2} = ones(dpts-2, npts)

	t = step
	for i = 2:dpts-2
		den1::Float64 = t * (1 - t)
		den2::Float64 = den1^2
		for j = 0:n
			Der1[i-1, j+1] = (j - n * t) / den1
			Der2[i-1, j+1] = ((j - n * t)^2 - n * t^2 - j * (1 - 2 * t)) / den2
		end
		t = t + step
	end

	EV::Array{Array{Int64,1},1} = Array{Int64,1}[] #1-dimensional cellular complex used for plotting the curve with Plasm package

	for i = 1 : dpts-1
		app = Array{Int64}(2)
		app[1] = i
		app[2] = i+1
		push!(EV,app)
	end

	DEV::Array{Array{Int64,1},1} = EV[1:dpts-3] #1-dimensional cellular complex used for plotting the Derivatives with Plasm package

	F::Array{Float64,2} = T*C*D
	F1::Array{Float64,2} = Der1.*(F[range(2,dpts-2),:])
	F2::Array{Float64,2} = Der2.*(F[range(2,dpts-2),:])
	P::Array{Float64,2} = transpose(F*G)
	D1::Array{Float64,2} = transpose(F1*G)
	D2::Array{Float64,2} = transpose(F2*G)

	return (P, D1, D2, EV, DEV)
end

#-----------------------------------------------------------------------

"""
	bern_basis(n, i, t)

Calculate the Bernstein basis.

This method is meant to evaluate the `i`-th Bernstein Basis of order `n` in the parameter `t`.
Of course this method could be implemented in a more quick way, in order to evaluate all the basis of a specific order
but this is out of the scope of this library.

---

# Arguments

- `n::Int64`: the order of the basis (`n = npts-1`).
- `i::Int64`: the `i`-th basis of the given order. `i ∈ [0, …, n]`.
- `t::Float64`: the parameter of the basis `N_{n, i}(t)`.

---

# Examples

```jldoctest
julia> map(x -> bern_basis(3, 0, x), [0, 0.15, 0.35, 0.5, 0.65, 0.85, 1])
7-element Array{Float64,1}:
 1.0
 0.614125
 0.274625
 0.125
 0.042875
 0.003375
 0.0
```

```jldoctest
julia> map(x -> bern_basis(3, 3, x), [0, 0.15, 0.35, 0.5, 0.65, 0.85, 1])
7-element Array{Float64,1}:
 0.0
 0.003375
 0.042875
 0.125
 0.274625
 0.614125
 1.0
```

---

_By Elia Onofri_
"""

function bern_basis(n::Int64, i::Int64, t::Float64)::Float64
	@assert i >= 0 && i<= n ("ERROR: there are $(n+1) basis of order $(n). You cannot choose the $(i)-th one.")
	return binomial(n, i) * t^i * (1-t)^(n-i)
end


text:
export reaction_process_setup, langevin_process_setup, stoich_bounds, transform_state!, rescale_state!,
       value_function

stoichmat(rn::ReactionSystem) = prodstoichmat(rn) - substoichmat(rn)

reformat_jumps(S::Matrix, species_to_index::Dict, x::AbstractVector) = [x .+ S[:,i] for i in 1:size(S,2)]

∂(p,x) = differentiate.(p,x)
∂²(p,x,y) = differentiate(∂(p,x),y)
∂(X::BasicSemialgebraicSet) = [intersect(@set(X.p[i] == 0), [@set(X.p[k] >= 0) for k in 1:length(X.p) if k != i]..., X.V) for i in 1:length(X.p)]

function reformat_reactions(rxns::Vector{Reaction}, species_to_index::Dict, x::Vector{<:PV}, params = Dict())
    props = Polynomial{true,Float64}[]
    for r in rxns
        @unpack rate, substrates, substoich, only_use_rate = r
        rate = rate isa Sym ? params[rate] : rate
        a = rate*polynomial(MonomialVector(x,0))
        if !only_use_rate
            for (s, ν) in enumerate(substoich)
                idx = species_to_index[substrates[s]]
                a *= prod(x[idx] - i for i in 0:ν-1)/factorial(ν) # consistent with rxns
            end
        end
        push!(props, a)
    end
    return props
end

function stoich_bounds(rn::ReactionSystem, x0::Dict, solver)
    S = stoichmat(rn)
    B = nullspace(S)
    specs = species(rn)
    x_init = [x0[spec] for spec in specs]
    m = Model(solver)
    @variable(m, x[1:length(specs)] >= 0)
    @constraint(m, B'*x .== B'*x_init)
    scales = Dict()
    for i in 1:length(specs)
        @objective(m, Max, x[i])
        optimize!(m)
        if termination_status(m) == MOI.DUAL_INFEASIBLE
            @warn  "state space unbounded, autoscaling failed!"
            scales[specs[i]] = 1.0
        else
            scales[specs[i]] = min(1.0, objective_value(m))
        end
    end
    return scales
end

"""
    reaction_process_setup(rn::ReactionSystem, x0::Dict;
                    scales = Dict(s => 1.0 for s in species(rn)),
                    auto_scaling = false, solver = nothing)

    transforms a reaction network as defined by Catalyst.jl's ReactionSystem type
    into an equivalent ReactionProcess accounting for reaction invariants and
    scales.
"""
function reaction_process_setup(rn::ReactionSystem, x0::Dict;
                                scales = Dict(s => 1.0 for s in species(rn)),
                                auto_scaling = false, solver = nothing,
                                params::Dict = Dict())
    P = ReactionProcess(rn, params)
    project_into_subspace!(P, Float64[x0[s] for s in species(rn)])
    if auto_scaling
        if solver isa nothing
            error("Need specification of LP solver to perform automatic scaling")
        end
        scales = stoich_bounds(rn, x0, solver)
    end
    x_scale = [scales[P.state_to_species[x]] for x in P.JumpProcess.x]
    x0 = [x0[P.state_to_species[x]]/scales[P.state_to_species[x]] for x in P.JumpProcess.x]
    rescale_state!(P, x_scale)
    return P, x0
end

function langevin_process_setup(rn::ReactionSystem, x0::Dict;
                                scales = Dict(s => 1.0 for s in species(rn)),
                                auto_scaling = false, solver = nothing,
                                params::Dict = Dict())
    P = LangevinProcess(rn, params)
    project_into_subspace!(P, Float64[x0[s] for s in species(rn)])
    if auto_scaling
        if solver isa nothing
            error("Need specification of LP solver to perform automatic scaling")
        end
        scales = stoich_bounds(rn, x0, solver)
    end
    x_scale = [scales[P.state_to_species[x]] for x in P.DiffusionProcess.x]
    x0 = [x0[P.state_to_species[x]]/scales[P.state_to_species[x]] for x in P.DiffusionProcess.x]
    rescale_state!(P, x_scale)
    return P, x0
end

function reaction_process_setup(rn::ReactionSystem; scales = Dict(s => 1.0 for s in species(rn)), params::Dict = Dict())
    P = ReactionProcess(rn, params)
    x_scale = Float64[scales[P.state_to_species[x]] for x in P.JumpProcess.x]
    rescale_state!(P, x_scale)
    return P
end

function langevin_process_setup(rn::ReactionSystem;
                                scales = Dict(s => 1.0 for s in species(rn)), params::Dict = Dict())
    P = LangevinProcess(rn, params)
    x_scale = [scales[P.state_to_species[x]] for x in P.DiffusionProcess.x]
    rescale_state!(P, x_scale)
    return P
end

function transform_state!(P::JumpProcess, x::AbstractVector, z::AbstractVector; iv::AbstractVector = 1:length(P.x))
    Π = P.x => polynomial.(z) #[P.x[i] => z[i] for i in 1:length(P.x)]
    P.a = map(p -> p(Π), P.a)
    P.h = [map(p -> p(Π), h[iv]) for h in P.h]
    P.X = subs_X(P.X, Π) #P.X isa FullSpace ? P.X : intersect([@set(p(Π) >= 0) for p in P.X.p]...)
    P.x = x
end

function transform_state!(P::ReactionProcess, x::AbstractVector, z::AbstractVector; iv::AbstractVector = 1:length(P.JumpProcess.x))
    Π = P.JumpProcess.x => polynomial.(z) #[P.JumpProcess.x[i] => z[i] for i in 1:length(P.JumpProcess.x)]
    P.species_to_state = Dict(spec => P.species_to_state[spec](Π) for spec in keys(P.species_to_state))
    P.state_to_species = Dict(P.species_to_state[spec] => spec for spec in keys(P.species_to_state))
    transform_state!(P.JumpProcess, x, z; iv = iv)
end

function transform_state!(P::DiffusionProcess, x::AbstractVector, z::AbstractVector; iv::AbstractVector = 1:length(P.x))
    Π = P.x => polynomial.(z)
    P.f = map(p -> p(Π), P.f[iv])
    P.σ = map(p -> p(Π), P.σ[iv, iv])
    P.X = subs_X(P.X, Π) #P.X isa FullSpace ? P.X : intersect([@set(p(Π) >= 0) for p in P.X.p]...)
    P.x = x
end

function transform_state!(P::LangevinProcess, x::AbstractVector, z::AbstractVector; iv::AbstractVector = 1:length(P.JumpProcess.x))
    Π = P.DiffusionProcess.x => polynomial.(z)
    P.species_to_state = Dict(spec => P.species_to_state[spec](Π) for spec in keys(P.species_to_state))
    P.state_to_species = Dict(P.species_to_state[spec] => spec for spec in keys(P.species_to_state))
    transform_state!(P.DiffusionProcess, x, z; iv = iv)
end

function transform_state!(P::JumpDiffusionProcess, x::AbstractVector, z::AbstractVector; iv::AbstractVector = 1:length(P.x))
    Π = P.x => polynomial.(z) #[P.x[i] => z[i] for i in 1:length(P.x)]
    P.a = map(p -> p(Π), P.a)
    P.h = [map(p -> p(Π), h[iv]) for h in P.h]
    P.f = map(p -> p(Π), P.f[iv])
    P.σ = map(p -> p(Π), P.σ[iv, iv])
    P.X = subs_X(P.X, Π) #P.X isa FullSpace ? P.X : intersect([@set(p(Π) >= 0) for p in P.X.p]...)
    P.x = x
end

function rescale_state!(P::JumpProcess, x0::Vector{<:Real})
    transform_state!(P, P.x, P.x .* x0)
    for i in 1:length(P.h)
        P.h[i] ./= x0
    end
end

function rescale_state!(P::ReactionProcess, x0::Vector{<:Real})
    transform_state!(P, P.JumpProcess.x, polynomial.(P.JumpProcess.x .* x0))
    for i in 1:length(P.JumpProcess.h)
        P.JumpProcess.h[i] ./= x0
    end
end

function rescale_state!(P::DiffusionProcess, x0::Vector{<:Real})
    transform_state!(P, P.x, P.x .* x0)
    P.f ./= x0
    P.σ ./= x0*x0'
end

function rescale_state!(P::LangevinProcess, x0::Vector{<:Real})
    rescale_state!(P.DiffusionProcess, x0)
end

function rescale_state!(P::JumpDiffusionProcess, x0::Vector{<:Real})
    transform_state!(P, P.x, P.x .* x0)
    for i in 1:length(P.h)
        P.h[i] ./= x0
    end
    P.f ./= x0
    P.σ ./= x0*x0'
end

function partition_variables(B)
    m, n = size(B)
    @assert(m <= n)
    Q, R = qr(B)
    ds = []
    i, k = 1, 1
    while k <= n && i <= m
        if abs(R[i,k]) >= 1e-12
            push!(ds, k)
            i += 1
        end
        k += 1
    end
    is = setdiff(1:n, ds)
    return is, ds, Q, R
end

function project_into_subspace!(P::MarkovProcess, B, f::Vector{<:Real})
    iv, dv, Q, R = partition_variables(B)
    if !isempty(dv)
        z = Vector{Polynomial{true, Float64}}(undef, size(B,2))
        @polyvar(x[1:length(iv)])
        z[iv] .= polynomial.(x)
        z[dv] .= polynomial.(round.(R[:,dv]\(Q'*f), digits = 12) - round.((R[:,dv]\R[:,iv]), digits = 12)*x)
        transform_state!(P, x, z; iv=iv)
    end
end

function project_into_subspace!(P::Union{ReactionProcess,LangevinProcess}, x0::Vector{<:Real})
    B = transpose(nullspace(transpose(stoichmat(P.ReactionSystem))))
    if !isempty(B)
        f = B*x0
        project_into_subspace!(P, B, f)
    end
end

function init_moments(x, x0, d)
    n = length(x0)
    moms = Dict()
    for i in 1:(d+1)^n
        midx = invert_index(i, (d+1)*ones(Int64,n)) .- 1
        moms[prod(x.^midx)] = prod(x0.^midx)
    end
    return moms
end

init_moments(x, x0, d, P::Partition) =  Dict(P.get_vertex(x0) => init_moments(x, x0, d))

linearize_index(idx,rs) = idx[1] + (length(idx) > 1 ? sum((idx[i] - 1) * prod(rs[1:i-1]) for i in 2:length(idx)) : 0)

function invert_index(idx, rs)
    inv_idx = similar(rs)
    for i in length(rs):-1:2
        fac = prod(rs[1:i-1])
        n = div(idx - 1, fac)
        inv_idx[i] = n + 1
        idx -= n*fac
    end
    inv_idx[1] = idx
    return inv_idx
end

function expectation(w::Polynomial, μ::Dict)
    ex = AffExpr(0.0)
    for i in 1:length(w.a)
        if w.a[i] != 0 && μ[w.x[i]] != 0
            ex += w.a[i] * μ[w.x[i]]
        end
    end
    return ex
end

expectation(w::VariableRef, μ::Dict) = w*μ[1]
expectation(w::AbstractVector, μ::Dict) = sum(expectation(w[v], μ[v]) for v in keys(μ))
expectation(w::Dict, μ::Dict) = sum(expectation(w[v], μ[v]) for v in keys(μ))

function value_function(CP::ControlProcess, trange, bound::Bound)
    if trange[1] == 0
        trange = trange[2:end]
    end
    V_poly(t,x) = bound.w[get_piece(t, trange), bound.partition.get_vertex(x)]
    V_val(t,x) = V_poly(t,x)(CP.MP.iv => t, CP.MP.x => x)
    return V_poly, V_val
end

function get_piece(t, trange)
    if t > trange[end] || t < 0
        @warn string("t = ", t, " outside the domain of V")
        i = t < 0 ? 1 : length(trange)
    elseif t == trange[end]
        i = length(trange)
    else
        i = findfirst(ti -> ti > t, trange)
    end
    return i
end

function trivial_partition(X::AbstractSemialgebraicSet)
    G = MetaDiGraph()
    add_vertex!(G, :cell, X)
    return Partition(G, x -> 1)
end

function split_state_space(MP::MarkovProcess, X::BasicSemialgebraicSet)
    Xc = complement(X, MP.X)
    G = MetaDiGraph()
    add_vertex!(G, 1, :cell, X)
    add_vertex!(G, 2, :cell, Xc)
    add_edge!(G, 1, 2, ∂(X))
    return Partition(G, x -> check_membership(x, X) ? 1 : 2)
end

function dual_poly(w, t, trange)
    if trange[1] == 0
        trange = trange[2:end]
    end
    return Dict(key => subs(value(w[key]), t => key[1] > 1 ? (t - trange[key[1]-1])/(trange[key[1]] - trange[key[1]-1]) : t/trange[key[1]]) for key in keys(w))
end

function subs_X(X::BasicSemialgebraicSet, submap)
    eqs = Polynomial{true, Float64}[]
    for eq in equalities(X)
        push!(eqs, eq(submap))
    end
    ineqs = Polynomial{true, Float64}[]
    for eq in inequalities(X)
        push!(ineqs, eq(submap))
    end
    return BasicSemialgebraicSet(algebraicset(eqs), ineqs)
end

inequalities(::FullSpace) = []
polynomial(a::Real) = polynomial(Term{true}(a))


text:
using BandedMatrices, LinearAlgebra, ArrayLayouts, Random, Test, GenericLinearAlgebra
import BandedMatrices: MemoryLayout, SymmetricLayout, HermitianLayout, BandedColumns

@testset "Symmetric" begin
    A = Symmetric(brand(10,10,1,2))
    @test isbanded(A)
    @test BandedMatrix(A) == A
    @test bandwidth(A) == bandwidth(A,1) == bandwidth(A,2) ==  2
    @test bandwidths(A) == bandwidths(BandedMatrix(A)) == (2,2)
    @test MemoryLayout(typeof(A)) == SymmetricLayout{BandedColumns{DenseColumnMajor}}()

    @test A[1,2] == A[2,1]
    @test A[1,4] == 0
    x=rand(10)
    @test A*x ≈ Matrix(A)*x
    y = similar(x)
    @test all(A*x .=== muladd!(1.0,A,x,0.0,similar(x)) .===
                BLAS.sbmv!('U', 2, 1.0, parent(A).data, x, 0.0, similar(x)))

    @test A[2:10,1:9] isa BandedMatrix
    @test [A[k,j] for k=2:10, j=1:9] == A[2:10,1:9]
    A = Symmetric(brand(10,10,1,2),:L)
    @test isbanded(A)
    @test BandedMatrix(A) == A
    @test bandwidth(A) == bandwidth(A,1) == bandwidth(A,2) ==  1
    @test bandwidths(A) == bandwidths(BandedMatrix(A)) == (1,1)
    @test MemoryLayout(typeof(A)) == SymmetricLayout{BandedColumns{DenseColumnMajor}}()

    @test A[1,2] == A[2,1]
    @test A[1,3] == 0
    @test A[2:10,1:9] isa BandedMatrix
    @test [A[k,j] for k=2:10, j=1:9] == A[2:10,1:9]

    x=rand(10)
    @test A*x ≈ Matrix(A)*x
    y = similar(x)
    @test all(A*x .=== muladd!(1.0,A,x,0.0,y) .===
                BLAS.sbmv!('L', 1, 1.0, view(parent(A).data,3:4,:), x, 0.0, similar(x)))

    @test norm(A) ≈ norm(Matrix(A))
    @test cond(A) ≈ cond(Matrix(A))

    # (generalized) eigen & eigvals
    Random.seed!(0)

    A = brand(Float64, 100, 100, 2, 4)
    std = eigvals(Symmetric(Matrix(A)))
    @test eigvals(Symmetric(A)) ≈ std
    @test eigvals(Hermitian(A)) ≈ std
    @test eigvals(Hermitian(big.(A))) ≈ std

    A = brand(ComplexF64, 100, 100, 4, 0)
    @test Symmetric(A)[2:10,1:9] isa BandedMatrix
    @test Hermitian(A)[2:10,1:9] isa BandedMatrix
    @test [Symmetric(A)[k,j] for k=2:10, j=1:9] == Symmetric(A)[2:10,1:9]
    @test [Hermitian(A)[k,j] for k=2:10, j=1:9] == Hermitian(A)[2:10,1:9]

    std = eigvals(Hermitian(Matrix(A), :L))
    @test eigvals(Hermitian(A, :L)) ≈ std
    @test eigvals(Hermitian(big.(A), :L)) ≈ std

    A = Symmetric(brand(Float64, 100, 100, 2, 4))
    F = eigen(A)
    Λ, Q = F
    @test Q'Matrix(A)*Q ≈ Diagonal(Λ)
    FD = convert(Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}, F)
    @test FD.vectors'Matrix(A)*FD.vectors ≈ Diagonal(F.values)

    function An(::Type{T}, N::Int) where {T}
        A = Symmetric(BandedMatrix(Zeros{T}(N,N), (0, 2)))
        for n = 0:N-1
            parent(A).data[3,n+1] = T((n+1)*(n+2))
        end
        A
    end

    function Bn(::Type{T}, N::Int) where {T}
        B = Symmetric(BandedMatrix(Zeros{T}(N,N), (0, 2)))
        for n = 0:N-1
            parent(B).data[3,n+1] = T(2*(n+1)*(n+2))/T((2n+1)*(2n+5))
        end
        for n = 0:N-3
            parent(B).data[1,n+3] = -sqrt(T((n+1)*(n+2)*(n+3)*(n+4))/T((2n+3)*(2n+5)*(2n+5)*(2n+7)))
        end
        B
    end

    for T in (Float32, Float64)
        A = An(T, 100)
        B = Bn(T, 100)

        λ = eigvals(A, B)
        @test λ ≈ eigvals(Symmetric(Matrix(A)), Symmetric(Matrix(B)))

        err = λ*(T(2)/π)^2 ./ (1:length(λ)).^2 .- 1

        @test norm(err[1:40]) < 100eps(T)

        Λ, V = eigen(A, B)
        @test V'Matrix(A)*V ≈ Diagonal(Λ)
        @test V'Matrix(B)*V ≈ I
    end
end

@testset "Hermitian" begin
    T = ComplexF64
    A = Hermitian(brand(T,10,10,1,2))
    @test isbanded(A)
    @test BandedMatrix(A) == A
    @test bandwidth(A) == bandwidth(A,1) == bandwidth(A,2) ==  2
    @test bandwidths(A) == bandwidths(BandedMatrix(A)) == (2,2)
    @test MemoryLayout(typeof(A)) == HermitianLayout{BandedColumns{DenseColumnMajor}}()

    @test A[1,2] == conj(A[2,1])
    @test A[1,4] == 0
    x=rand(T,10)
    @test A*x ≈ Matrix(A)*x
    @test all(A*x .=== materialize!(MulAdd(one(T),A,x,zero(T),copy(x))) .===
                BLAS.hbmv!('U', 2, T(1.0), parent(A).data, x, T(0.0), similar(x)))

    A = Hermitian(brand(T,10,10,1,2),:L)
    @test isbanded(A)
    @test BandedMatrix(A) == A
    @test bandwidth(A) == bandwidth(A,1) == bandwidth(A,2) ==  1
    @test bandwidths(A) == bandwidths(BandedMatrix(A)) == (1,1)
    @test MemoryLayout(typeof(A)) == HermitianLayout{BandedColumns{DenseColumnMajor}}()

    @test A[1,2] == conj(A[2,1])
    @test A[1,3] == 0
    x=rand(T, 10)
    @test A*x ≈ Matrix(A)*x

    @test all(A*x .=== muladd!(one(T),A,x,zero(T),copy(x)) .===
                materialize!(MulAdd(one(T),A,x,zero(T),copy(x))) .===
                BLAS.hbmv!('L', 1, one(T), view(parent(A).data,3:4,:), x, zero(T), similar(x)))
end

@testset "LDLᵀ" begin
    for T in (Float16, Float32, Float64, BigFloat, Rational{BigInt})
        A = BandedMatrix{T}(undef,(10,10),(0,2))
        A[band(0)] .= 4
        A[band(1)] .= -one(T)/4
        A[band(2)] .= -one(T)/16
        SAU = Symmetric(A, :U)
        F = ldlt(SAU)
        b = collect(one(T):size(F, 1))
        x = Matrix(SAU)\b
        y = F\b
        @test x ≈ y
        A = BandedMatrix{T}(undef,(10,10),(2,0))
        A[band(0)] .= 4
        A[band(-1)] .= -one(T)/3
        A[band(-2)] .= -one(T)/9
        SAL = Symmetric(A, :L)
        x = Matrix(SAL)\b
        F = ldlt(SAL)
        y = F\b
        @test x ≈ y
        @test_throws DimensionMismatch F\[b;b]

        T ≠ Float16 && (@test det(F) ≈ det(SAL))
    end
    for T in (Int16, Int32, Int64, BigInt)
        A = BandedMatrix{T}(undef, (4,4), (1,1))
        A[band(0)] .= 3
        A[band(1)] .= 1
        F = ldlt(Symmetric(A))
        @test eltype(F) == float(T)
    end

    @testset "#175" begin
        A = BandedMatrix(0 => [1,2])
        F = ldlt(Symmetric(A))
        @test F.d == [1.,2.]
        @test F.D == Diagonal([1.,2.])
    end
end

@testset "Generalized eigenvalues $W{$T}($Ua,$Ub)($n,$wa-$wb)" for (T,W) in (
                                    (Float32, Symmetric),
                                    (Float64, Symmetric),
                                    #(Float32, Hermitian),
                                    #(Float64, Hermitian),
                                    #(ComplexF32, Hermitian),
                                    #(ComplexF64, Hermitian),
                                   ),
    (Ua, Ub) in  ((:L,:L), (:U,:U)),
    (wa, wb) in ((2, 3), (3, 2)), n in (4,)
    #
    function sbmatrix(W, T, U, w, n)
        r = U == :L ? (0:-1:-w+1) : (0:w-1)
        band(k) = k => ones(T, n - abs(k)) * 2.0^-abs(k)
        W(BandedMatrix(band.(r)...), U)
    end
    A = sbmatrix(W, T, Ua, wa, n)
    B = sbmatrix(W, T, Ub, wb, n)
    @test eigvals(A, B) ≈ eigvals(Matrix(A), Matrix(B))
end

@testset "Cholesky" begin
    for T in (Float64, BigFloat)
        A = Symmetric(BandedMatrix(0 => one(T) ./ [12, 6, 6, 6, 12],
                                1 => ones(T,4) ./ 24))
        Ac = cholesky(A)

        @test Ac isa Cholesky{T,<:BandedMatrix{T}}
        @test Ac.U ≈ cholesky(Matrix(A)).U

        b = rand(T,size(A,1))
        @test Ac\b ≈ Matrix(A)\b
        @test_broken Ac\b ≈ A\b
    end

    for T in (Float64, BigFloat)
        A = Symmetric(BandedMatrix(0 => one(T) ./ [12, 6, 6, 6, 12],
                                -1 => ones(T,4) ./ 24), :L)
        Ac = cholesky(A)

        @test Ac isa Cholesky{T,<:BandedMatrix{T}}
        @test Ac.L ≈ cholesky(Matrix(A)).L

        b = rand(T,size(A,1))
        @test Ac\b ≈ Matrix(A)\b
        @test_broken Ac\b ≈ A\b
    end

    let T = ComplexF64
        A = Hermitian(BandedMatrix(0 => one(T) ./ [12, 6, 6, 6, 12],
                                   1 => ones(T,4) ./ 24))
        Ac = cholesky(A)

        @test Ac isa Cholesky{T,<:BandedMatrix{T}}
        @test Ac.U ≈ cholesky(Matrix(A)).U

        b = rand(T,size(A,1))
        @test Ac\b ≈ Matrix(A)\b
        @test_broken Ac\b ≈ A\b
    end

    @testset "UnitRange" begin
        A = brand(100,100,3,2)
        S = Symmetric(A)
        @test S[3:10,4:11] == Symmetric(Matrix(A))[3:10,4:11]
    end

    @testset "Sym of degenerate bands" begin
        A = SymTridiagonal(zeros(8), fill(0.5,7))
        B = Symmetric(BandedMatrix(1 => fill(0.5,7)))
        @test A ≈ B
        @test BandedMatrices.inbands_getindex(B, 1, 1) == 0
        @test BandedMatrices.inbands_getindex(B, 1, 2) == BandedMatrices.inbands_getindex(B, 2, 1) == 0.5
    end
end