text:
#### Hyperdiffusion Model Functions
using DocStringExtensions
using LinearAlgebra
export HyperDiffusion, NoHyperDiffusion, StandardHyperDiffusion

abstract type HyperDiffusion end
vars_state(::HyperDiffusion, FT) = @vars()
vars_aux(::HyperDiffusion, FT) = @vars()
vars_gradient(::HyperDiffusion, FT) = @vars()
vars_gradient_laplacian(::HyperDiffusion, FT) = @vars()
vars_diffusive(::HyperDiffusion, FT) = @vars()
vars_hyperdiffusive(::HyperDiffusion, FT) = @vars()
function atmos_init_aux!(
    ::HyperDiffusion,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
) end
function atmos_nodal_update_aux!(
    ::HyperDiffusion,
    ::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function gradvariables!(
    ::HyperDiffusion,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function hyperdiffusive!(
    h::HyperDiffusion,
    hyperdiffusive::Vars,
    gradvars::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function flux_diffusive!(
    h::HyperDiffusion,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
) end
function diffusive!(
    h::HyperDiffusion,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end

"""
  NoHyperDiffusion <: HyperDiffusion
Defines a default hyperdiffusion model with zero diffusive fluxes.
"""
struct NoHyperDiffusion <: HyperDiffusion end

"""
  StandardHyperDiffusion{FT} <: HyperDiffusion
Horizontal hyperdiffusion methods for application in GCM and LES settings
Timescales are prescribed by the user while the diffusion coefficient is
computed as a function of the grid lengthscale.
"""
struct StandardHyperDiffusion{FT} <: HyperDiffusion
    τ_timescale::FT
end

vars_aux(::StandardHyperDiffusion, FT) = @vars(Δ::FT)
vars_gradient(::StandardHyperDiffusion, FT) =
    @vars(u::SVector{3, FT}, h_tot::FT)
vars_gradient_laplacian(::StandardHyperDiffusion, FT) =
    @vars(u::SVector{3, FT}, h_tot::FT)
vars_hyperdiffusive(::StandardHyperDiffusion, FT) =
    @vars(ν∇³u::SMatrix{3, 3, FT, 9}, ν∇³h_tot::SVector{3, FT})

function atmos_init_aux!(
    ::StandardHyperDiffusion,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.hyperdiffusion.Δ = lengthscale(geom)
end

function gradvariables!(
    h::StandardHyperDiffusion,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.hyperdiffusion.u = ρinv * state.ρu
    transform.hyperdiffusion.h_tot = transform.h_tot
end

function hyperdiffusive!(
    h::StandardHyperDiffusion,
    hyperdiffusive::Vars,
    hypertransform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ∇Δu = hypertransform.hyperdiffusion.u
    ∇Δh_tot = hypertransform.hyperdiffusion.h_tot
    # Unpack
    τ_timescale = h.τ_timescale
    # Compute hyperviscosity coefficient
    ν₄ = (aux.hyperdiffusion.Δ / 2)^4 / 2 / τ_timescale
    hyperdiffusive.hyperdiffusion.ν∇³u = ν₄ * ∇Δu
    hyperdiffusive.hyperdiffusion.ν∇³h_tot = ν₄ * ∇Δh_tot
end

function flux_diffusive!(
    h::StandardHyperDiffusion,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    hyperdiffusive::Vars,
    aux::Vars,
    t::Real,
)
    flux.ρu += state.ρ * hyperdiffusive.hyperdiffusion.ν∇³u
    flux.ρe += state.ρ * hyperdiffusive.hyperdiffusion.ν∇³h_tot
end

max_stars_repo_name: yairchn/CLIMA

max_issues_repo_name: yairchn/CLIMA

max_forks_repo_name: yairchn/CLIMA


text:
module Libm

const FloatTypes=Union{Float32,Float64}

include("utils.jl")
include("erf.jl")
include("log/tang.jl")

end

ext: jl


text:
#=

  Euler #39 in Julia.

  """
  If p is the perimeter of a right angle triangle with integral length sides,
  {a,b,c}, there are exactly three solutions for p = 120.

  {20,48,52}, {24,45,51}, {30,40,50}

  For which value of p <= 1000, is the number of solutions maximised?
  """

  This Julia model was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.11346191s
# This is surprinsingly slow!
function euler39a()
    n = 1000;
    # squares = Set((1:n).|>i->i*i) # slightly slower
    squares = Dict((1:n).|>i->i*i=>1)
    valid = []
    for x in keys(squares)
        for y in keys(squares)
            if x < y && haskey(squares,x+y) && (sqrt(x) + sqrt(x) + sqrt(x+y)) < 1000
                push!(valid,[x, y])
            end
        end
    end

    counts = Dict()
    for x in valid
        c = floor(Int,sqrt(x[1]) + sqrt(x[2]) + sqrt(x[1]+x[2]))
        get!(counts,c,0)
        counts[c] += 1
   end

    # find max count
    maxV = maximum(values(counts))
    t = [i for i in keys(counts) if counts[i] == maxV]
    return t[1]
end

run_euler(euler39a)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
module continuousEntropyEstimators

using Distances, NearestNeighbors, SpecialFunctions
# using Statistics
# TODO: implement version in 1D, possibly 2D that doesn't rely on trees
# TODO: implement version that runs in parallel
export entropyEstimator, kraskovEstimator, NNestimator, entropy, mutualInformation, conditionalMutualInformation

# because I like math symbols
Γ(x) = gamma(x)
abstract type entropyEstimator end
mutable struct kraskovEstimator<: entropyEstimator end
mutable struct NNestimator<: entropyEstimator end

"""
    entropy(data::Array{T,2},  estimator::kraskovEstimator=kraskovEstimator(); k::Int64=3, logbase::S=2, treeType=Chebyshev()) where {T<:AbstractFloat, S<:Real}
    entropy(data::Array{Int64}, logbase::Float64=2.0)

computes the entropy of the data, which should be a d×N array.

If data is an array of floats, it is treated as a differential entropy calculation, by default done using knn's with k=3. treeType should be a distance from Distances. Some options include

| Distance     | As a p norm |
|--------------|:-----------:|
| Chebyshev()  | ``∞``       |
| Euclidean()  | ``2``       |
| Cityblock()  | ``1``       |
| Minkowski(p) | ``p``       |

If data is an array of ints, it uses the plug-in estimator.

TODO: Calculate the entropy in the method of Gao et al for mixed distributions
"""
function entropy(data::Array{T,2},  estimator::kraskovEstimator=kraskovEstimator(); k::Int64=3, logbase::S=2, treeType=Chebyshev()) where {T<:AbstractFloat, S<:Real}
    d, N = size(data)
    # perturb all the points to avoid problems with deterministic functions
    X = copy(data) + 1e-10*randn(size(data))
    tree = KDTree(X, treeType)
    if treeType==Chebyshev()
        cb = d*log(logbase, 2)
    elseif treeType==Euclidean()
        p=2
    elseif treeType==Cityblock()
        p=1
    elseif typeof(treeType)==Minkowski{Int64} || typeof(treeType)==Minkowski{Float64}
        p = treeType.p
    end
    if treeType!=Chebyshev()
        print(treeType)
        cb = d*log(logbase,2) + d*log(logbase, Γ(1+1/p))- log(logbase, Γ(1+d/p))
    end
    entropyEstimate = digamma(N)-digamma(k) + cb + d*log(logbase, 2) # this last factor is left implicit in eq 20
    for i=1:N
        _, dist = knn(tree, X[:,i], k+1, true)
        # check to see if this point is discrete, in which case Gao et al suggest using a plug-in estimator
        # if dist[end]<1e-10
        entropyEstimate += d/N*log(logbase, dist[end])
    end
    return entropyEstimate
end
continuousEntropyEstimators.entropy
# function entropy(data::Array{Float64,2}, estimator::NNestimator=NNestimator(); logbase=2, treeType=Chebyshev())
#     d, N = size(data)
#     # perturb all the points to avoid problems with deterministic functions
#     X = copy(data) + 1e-10*randn(size(data))
#     tree = KDTree(X, treeType)
#     entropyEstimate = γ/log(logbase) + log(logbase, 2)
#     for i=1:N
#         _, dist = knn(tree, X[:,i], 2, true)
#         entropyEstimate += 1/N*log(logbase, N*dist[end])
#     end
#     return entropyEstimate
# end
# unique([1.5, 1.5, 3, 201])
# data = rand([1,3,2], 10)
# uniqueVal = unique(data)
# data[data.==uniqueVal[1]]
# maximum(data.==uniqueVal[1])
function entropy(data::Array{Int64}, logbase::Float64=2.0)
    pEstimate, _ = pluginEstimator(data)
    pEstimate[pEstimate.==0] = 1
    return -pEstimate'*log.(logbase, pEstimate)
end

"""
    pluginEstimator(data::Array{Int64,(1,2)})
estimate the probability distribution of a discrete variable by using the counts
"""
function pluginEstimator(data::Array{Int64,2})
    Y = [data[:,i] for i=1:size(data, 2)]
    uniqueVal = unique(Y)
    counts = zeros(size(uniqueVal))
    for i = 1:size(counts,1)
        for j = 1:length(Y)
            counts[i] += Y[j]==uniqueVal[i]
        end
    end
    counts = counts/length(data)
    return (counts, uniqueVal)
end

function pluginEstimator(data::Array{Int64,1})
    uniqueVal = unique(data)
    counts = zeros(size(uniqueVal))
    for i = 1:size(counts,1)
        counts[i] = sum(data.==uniqueVal[i])
    end
    counts = counts/length(data)
    return (counts, uniqueVal)
end
# abs(entropy(rand([1,3], 1000000))-1)<1e-5 # uniform variable on two values has unit entropy base 2
"""
    mutualInformation(Xc::Array{Float64,2}, Y::Array{Float64,2}, k=3; logbase=2, treeType=Chebyshev(), version::Int64=2)
    mutualInformation(Xc::Array{Float64,2}, Yc::Array{Int64}; k=3, logbase=2, treeType=Chebyshev())
    mutualInformation(Xc::Array{Int64,2}, Yc::Array{Float64}; k=3, logbase=2, treeType=Chebyshev())

Calculate the mutual information between Xc and Yc using knn nearest neighbor-type methods.

make sure that the data is `d`×`N`, where `d` is the dimension and `N` is the number of examples. `Xc` and `Yc` should match in the second dimension, of course. `version` corresponds to the versions in the Kraskov et. al paper, and is either 1 or 2. Default is chosen for high dimension. It uses neighborhoods restricted to `Xc` (resp `Yc`) to calculate the contribution from that direction. In the case that one of the variables is discrete and want it treated as such, make sure it is integer. It is calculated using ``I(X;Y) = H(X)-H(X|Y)`` If both are discrete, it just uses ``H(X) + H(Y) - H(X,Y)``
NOTE: only Chebyshev is implemented at this time
"""
function mutualInformation(Xc::Array{Float64,2}, Yc::Array{Float64,2}; k=3, logbase=2, treeType=Chebyshev(), version::Int64=2)
    @assert size(Xc,2)==size(Yc,2)
    @assert treeType==Chebyshev()
    dx, N = size(Xc)
    dy, N = size(Yc)
    X = copy(Xc) + 1e-10*randn(dx, N)
    Y = copy(Yc) + 1e-10*randn(dy, N)
    XY = [X; Y]
    tree = KDTree(XY, treeType)
    entropyEstimate = digamma(k) + digamma(N)
    if version==2
        entropyEstimate-=1/k
    end
    xTree = KDTree(X, treeType)
    yTree = KDTree(Y, treeType)
    for i=1:N
        neighborhood, dist = knn(tree, XY[:,i], k+1, true)
        neighbor = neighborhood[end]
        dist = dist[end]
        if version==1
            # count the points within x(y) distance of ε/2
            nx = length(inrange(xTree, X[:,i], dist, false)) # note that we both need to add one for the algorithm and we "accidentally" overcount by one, so there's no change
            ny = length(inrange(yTree, Y[:,i], dist, false)) # same w/y
        elseif version==2
            εx = chebyshev(X[:, neighbor],X[:,i]) # \varepsilon
            nx = length(inrange(xTree, X[:,i], εx, false)) # note that we both need to add one for the algorithm and we "accidentally" overcount by one, so there's no change
            εy = chebyshev(Y[:, neighbor],Y[:,i])
            ny = length(inrange(yTree, Y[:,i], dist, false)) # note that we both need to add one for the algorithm and we "accidentally" overcount by one, so there's no change
        else
            error("``(version)'s not an allowed value of version. Try 1 or 2")
        end
        entropyEstimate -= 1/N*digamma(nx)
    end
    return entropyEstimate
end

function mutualInformation(Xc::Array{Float64,2}, Yc::Array{Int64}; k=3, logbase=2, treeType=Chebyshev())
    @assert size(Xc,2)==size(Yc,2)
    @assert treeType==Chebyshev()
    dx, N = size(Xc)
    dy, N = size(Yc)
    X = copy(Xc) + 1e-10*randn(dx, N)
    Y = copy(Yc)
    # if Xis pure continuous, Y discrete, it's best to use I(X,Y) = H(Y)-H(Y|X)
    # First order of business is getting a plug-in estimate for Y
    pEstimate, yVals = pluginEstimator(Y)
    Nxgy =pEstimate*N # the number of Xs given the value of y in index
    eeXgY = digamma(k) + digamma.(Nxgy) # a vector of entropy estimates given the value of y
    entropyEstimate = entropy(X,k=k, logbase=logbase, treeType=treeType)
    #
    if dy>1
        Yp = [data[:,i] for i=1:size(data, 2)]
        for i = 1:length(pEstimate)
            for j = 1:length(Yp)
                tmp = pEstimate[i]*entropy(X[:, Yp[j]==yVals[i]])
                entropyEstimate-=pEstimate[i]*entropy(X[:, Yp[j]==yVals[i]])
            end
        end
    else
        # Y is one dimensional
        for i = 1:length(pEstimate)
            entropyEstimate-=pEstimate[i]*entropy(X[:, reshape(Y.==yVals[i],(:))], k=k, logbase=logbase, treeType=treeType)
        end
    end
    return abs(entropyEstimate)
end

function mutualInformation(Xc::Array{Int64}, Yc::Array{Int64}, logbase::Float64=2.0)
    entropy(Xc,logbase=logbase) + entropy(Yc,logbase=logbase) - entropy([Xc;Yc],logbase=logbase)
end

# Various kinds of conditional mutual information
"""
    conditionalMutualInformation(Xc::Array{X},Yc::Array{Y}, Zc::Array{Z}, k::Int64=3, logbase::Float64=2.0, treeType=Chebyshev(), version::Int64=2)) where {X<:Real, Y<:Real, Z<:Real}

Calculate the conditional mutual information between `Xc`, `Yc` given `Zc`, ``I(X;Y|Z) = I(X,Z; Y) - I(Z;Y)``.
"""
function conditionalMutualInformation(Xc::Array{Float64},Yc::Array{Float64}, Zc::Array{Float64}, k::Int64=3, logbase::Float64=2.0, treeType=Chebyshev(), version::Int64=2)
    mutualInformation([Xc; Zc], Yc,k=k,logbase=logbase, treeType=treeType, version=version) - mutualInformation(Zc, Yc, k=k, logbase=logbase, treeType=treeType, version=version)
end
function conditionalMutualInformation(Xc::Array{Float64},Yc::Array{Int64}, Zc::Array{Float64}, k::Int64=3, logbase::Float64=2.0, treeType=Chebyshev(), version::Int64=2)
    mutualInformation([Xc; Zc], Yc,k=k,logbase=logbase, treeType=treeType) - mutualInformation(Zc, Yc, k=k, logbase=logbase, treeType=treeType)
end
end # module

# entropy([rand([1,3,2,4,2],5000)'; rand([4,3,2,1,4],5000)'])
# data = [rand([1,3,2],500)'; rand([4,3,2],500)']
# X = [data[:,i] for i=1:size(data, 2)]
# uniqueVal = unique(X)
# uniqueVal[1]
# entropy([2.5*rand([1,3,2],500)'; 2.5*rand([4,3,2],500)'])
# d = 20
# X = randn(d,  1000)
# Y = X+.01*randn(d, 1000)
# @time mutualInformation(X, Y, version=2)
# @time entropy(X)
# @time entropy(Y)
# get_entropy(rand(1000))
# d = 3
# X = randn(d, 100000)
# @time entropy(X, kraskovEstimator(), k=3, treeType = Chebyshev())- 1/2*log2((e*2π)^d)
# @time entropy(X, NNestimator(), treeType = Chebyshev())
# 1/2*log2((e*2π)^d)
# using InformationMeasures
# 1.3.^range(-.001, -15, 5)
# 1.3.^range(3, 42, 5)
# 2.^range(4,14,5)
# 1/2*log(2,(2π)^2)
# Ns = [10, 50, 100, 500, 1000] #, 3000, 5000, 8000, 10000]
# esttimatesk8 = [zeros(x) for x in Ns]
# for (i,estimates) in enumerate(esttimatesk8)
#     print(i)
#     @time for j=1:length(estimates)
#         data = randn(10,Ns[i])
#         estimates[j] = entropy(data,8)
#     end
# end
# # simple binning
#
# using Cumulants
# c = cumulants(data')
# using InformationMeasures
# using Plots
# plot(Ns, [mean.(esttimates), mean.(esttimatesk8)])
# std.(esttimates)
# mean(estimates1)
# std(estimates1)
# entropyEstimate
# @assert N>=d
# data+=1e-10*randn(size(data))
#
#
# test = randn(10,10000)
# get_entropy(test)
# 1/2*log2((2π)^10)
# function sReLU(x)
#     if x>=0
#         return x*exp(-1/(1-x)^2)
#     else
#         return 0
#     end
# end
#
# using Plots
# t=-.5:.00001:10
# plot(t,[sReLU.(t), t])


text:
using MathOptInterface, ConstraintProgrammingExtensions

using Test

const CP = ConstraintProgrammingExtensions
const COIT = CP.Test
const COIDT = CP.DeprecatedTest
const COIB = CP.Bridges
const MOI = MathOptInterface
const MOIB = MOI.Bridges
const MOIBC = MOIB.Constraint
const MOIT = MOI.DeprecatedTest
const MOIU = MOI.Utilities

include("Bridges/models.jl")

@testset "ConstraintProgrammingExtensions" begin
    include("sets.jl")
    include("moi_traits.jl")
    include("Bridges/Bridges.jl")
    include("FlatZinc/FlatZinc.jl")
    include("DeprecatedTest/DeprecatedTest.jl")
    include("Test/Test.jl")
end


text:
#=
This is an example parameter file.

Parameters not set here will be set to their default values. For default values and the meaning of parameters see the definition of struct
type pars in module BrexDefs.

The script needs to define the following variables: par_comm, par_diff, PS, t_brex (see below for their example values)

!Warning: The types of constants uses in this script must agree with the abstract types given in the struct definition. This mainly means that parameters which are supposed to be floats need to be given as Float constants (for example k_min = 1.0 is fine, k_min = 1 is not, because the program is trying to save an Integer variable into an Float type - no automatic conversion is performed in the current version of the program)

There are 2 groups of parameters:
    (1) Parameters specific to a stationary equilibrium (aggregate state)
        - set using par_diff
    (2) Parameters non-specific to stationary equilibrium
        - set using par_comm
=#

#(1) Parameters specific to a stationary equilibrium (aggregate state):

par_comm = ""

#Example: Set the discount factor and preferences in all states [everything else will be equal to default values]
#par_comm = "
#β = 0.95,
#χ = 2.0
#"

#=
Example2: Leaving everything at default but changing the probability of different Brexit outcomes (first number is probability of no Brexit, probabilities have to sum to 1)
par_comm = "
PS = [0.0,0.5,0.5]
"
=#

#(Warning: putting a comma after the last entry results in an error!)

#(2) Parameters non-specific to stationary equilibrium:
#The number of elements of par_diff determines the number of states everywhere in the program.
#The first state corresponds to pre-brexit state (initial stationary equilibrium)

#Example: Tarrifs before Brexit are 0, in soft Brexit they are 0.05, in hard Brexit they are 0.1. TFP is also possibly affected, but no change in the baseline..
par_diff = [
    "τ = 0.0,A=0.5",
    "τ = 0.05,A=0.5",
    "τ = 0.1,A=0.5"
]

#Warning: Setting the values of some variables to be different between Brexit states will break the program - usually these are parameters related to the solution technique (see brexDefs.jl for details). In general, changing values of parameters which have economic interpretation should be fine.

#=If loadData = true, the program reads initial guesses for the Stationary equilibrium and transition paths from folder results/loadFolder

This can be very useful when the computation takes a long time, or when convergence is problematic (then we can use as initial guess stationary equilibria and transition paths obtained for not too different set of parameters)

These parameters are optional, if not defined here the program works fine.

Warning: If the folder does not contain the data or contains wrong data (for example data obtained for a different number of grid points or length of transition path) the program will crash or may result in unpredictable behaviour!
=#

loadData = false
loadFolder = "20200810_213249_example"


text:
module Grover

using QuantExQASM.Oracle
using QuantExQASM.Diffusion
using QuantExQASM.GateOps
using QuantExQASM.Circuit

export run_grover

"""
    apply_h!(cct::Circuit.Circ, tgt, reg::Union{String, Nothing}=nothing)

Internally used Hadamard implementation. Uses GateOps definition. Override for custom use.
"""
function apply_h!(cct::Circuit.Circ, tgt, reg::Union{String, Nothing}=nothing)
    Circuit.add_gatecall!(cct, GateOps.hadamard(tgt, reg))
end

"""
    apply_x!(cct::Circuit.Circ, tgt, reg::Union{String, Nothing}=nothing)

Internally used PauliX implementation. Uses GateOps definition. Override for custom use.
"""
function apply_x!(cct::Circuit.Circ, tgt, reg::Union{String, Nothing}=nothing)
    Circuit.add_gatecall!(c, GateOps.pauli_x(tgt, reg))
end

"""
    calc_iterations(num_states::Integer)

Calculate the required number of iterations to maximise the state's amplitude.
"""
function calc_iterations(num_states::Integer)
    return pi*sqrt(num_states)/4
end

"""
    state_init!(cct::Circuit.Circ, qubit_indices::Vector)

Initialises the state to the required target; defaults to ``H^{\\otimes n}\\vert \\psi \\rangle`` . Override for custom initialisation.
"""
function state_init!(cct::Circuit.Circ, qubit_indices::Vector)
    for i in qubit_indices
        apply_h!(cct, i)
    end
end

"""
    apply_grover_iteration!(cct::Circuit.Circ, qubit_indices::Vector)

Applies a single Grover iteration. To be used following `mark_state!`
"""
function apply_grover_iteration!(cct::Circuit.Circ, qubit_indices::Vector)
    return Diffusion.apply_diffusion(cct, qubit_indices[1:end-1], qubit_indices[end])
end

"""
    apply_grover_iteration!(cct::Circuit.Circ, qubit_indices::Vector, qubit_aux_indices::Vector)

Applies a single Grover iteration. To be used following `mark_state!`
"""
function apply_grover_iteration!(cct::Circuit.Circ, qubit_indices::Vector, qubit_aux_indices::Vector)
    return Diffusion.apply_diffusion(cct, qubit_indices[1:end-1], qubit_aux_indices, qubit_indices[end])
end

"""
    mark_state!(cct::Circuit.Circ, state::Integer, qubit_indices::Vector)

Applies the state marking procedure of the Grover iteration.
"""
function mark_state!(cct::Circuit.Circ, state::Integer, qubit_indices::Vector)
    return Oracle.bitstring_phase_oracle(cct::Circuit.Circ, state, qubit_indices[1:end-1], qubit_indices[end])
end

"""
    mark_state!(cct::Circuit.Circ, state::Integer, qubit_indices::Vector, qubit_aux_indices::Vector)

Applies the state marking procedure of the Grover iteration.
"""
function mark_state!(cct::Circuit.Circ, state::Integer, qubit_indices::Vector, qubit_aux_indices::Vector)
    return Oracle.bitstring_phase_oracle(cct::Circuit.Circ, state, qubit_indices[1:end-1], qubit_aux_indices, qubit_indices[end])
end

"""
    run_grover!(cct::Circuit.Circ, qubit_indices::Vector, state::Integer)

Generates a Grover search circuit sample, marking the state defined by `state`
and performing iterations to amplify the desired result upon measurement.
"""
function run_grover!(cct::Circuit.Circ, qubit_indices::Vector, state::Integer)
    num_states = 2^(length(qubit_indices))
    state_init!(cct, qubit_indices)
    for i in range(1,stop=ceil(Integer, calc_iterations(num_states)))
        mark_state!(cct, state, qubit_indices)
        apply_grover_iteration!(cct, qubit_indices)
    end
end

"""
    run_grover!(cct::Circuit.Circ, qubit_indices::Vector, state::Integer)

Generates a Grover search circuit sample, marking the state defined by `state`
and performing iterations to amplify the desired result upon measurement.
"""
function run_grover!(cct::Circuit.Circ, qubit_indices::Vector, qubit_aux_indices::Vector, state::Integer)
    num_states = 2^(length(qubit_indices))
    state_init!(cct, qubit_indices)
    for i in range(1,stop=ceil(Integer, calc_iterations(num_states)))
        mark_state!(cct, state, qubit_indices, qubit_aux_indices)
        apply_grover_iteration!(cct, qubit_indices, qubit_aux_indices)
    end
end

end


text:
module AdaGram

using SharedArrays
import Base.view
using LinearAlgebra
using ResumableFunctions
using Distributed
import Base.Iterators.Stateful
using Printf

using Libdl

_c_update_z = C_NULL
_c_inplace_update = C_NULL
_c_skip_gram = C_NULL

function __init__()
	_libpath = string(@__DIR__, "/../lib/superlib")
	_c_superlib = dlopen(_libpath)
	global _c_update_z = dlsym(_c_superlib, :update_z)
	global _c_inplace_update = dlsym(_c_superlib, :inplace_update)
	global _c_skip_gram = dlsym(_c_superlib, :skip_gram)
end

sigmoid(x) = 1. / (1. + exp(-x))
log_sigmoid(x) = -log(1. + exp(-x))

Tsf = Float32
Tw = Int32

include("softmax.jl")

struct Dictionary
	word2id::Dict{AbstractString, Tw}
	id2word::Array{AbstractString}

	function Dictionary(id2word::Array{AbstractString})
		word2id = Dict{AbstractString, Int}()
		for v in 1:length(id2word)
			setindex!(word2id, v, id2word[v])
		end
		return new(word2id, id2word)
	end
end

struct VectorModel
	frequencies::DenseArray{Int64}
	code::DenseArray{Int8, 2}
	path::DenseArray{Int32, 2}
	In::DenseArray{Tsf, 3}
	Out::DenseArray{Tsf, 2}
	alpha::Float64
	d::Float64
	counts::DenseArray{Float32, 2}
end

M(vm::VectorModel) = size(vm.In, 1) #dimensionality of word vectors
T(vm::VectorModel) = size(vm.In, 2) #number of meanings
V(vm::VectorModel) = size(vm.In, 3) #number of words

# view(x::SharedArray, i1::Subs, i2::Subs) = view(sdata(x), i1, i2)
# view(x::SharedArray, i1::Subs, i2::Subs, i3::Subs) = view(sdata(x), i1, i2, i3)

function shared_rand(dims::Tuple, norm::T) where {T <: AbstractFloat}
	S = SharedArray{T}(dims; init = S -> begin
			chunk = localindices(S)
			chunk_size = length(chunk)
			data = rand(chunk_size)
			data .-= 0.5
			data ./= norm
			S[chunk] = data
		end)
	return S
end

function shared_zeros(::Type{T}, dims::Tuple) where {T <: Number}
	S = SharedArray{T}(dims; init = S -> begin
			chunk = localindices(S)
			sdata(S[chunk]) .= T(0)
		end)
	return S
end

function VectorModel(max_length::Int64, V::Int64, M::Int64, T::Int64=1, alpha::Float64=1e-2,
		d::Float64=0.)
	path = shared_zeros(Int32, (max_length, V))
	code = shared_zeros(Int8, (max_length, V))

	code[:] .= -1

	In =  shared_zeros(Float32, (M, T, V))
	Out = shared_zeros(Float32, (M, V))

	counts = shared_zeros(Float32, (T, V))

	frequencies = shared_zeros(Int64, (V,))

	return VectorModel(frequencies, code, path, In, Out, alpha, d, counts)
end

function VectorModel(freqs::Array{Int64}, M::Int64, T::Int64=1, alpha::Float64=1e-2,
	d::Float64=0.)
	V = length(freqs)

	nodes = build_huffman_tree(freqs)
	outputs = convert_huffman_tree(nodes, V)

	max_length = maximum(map(x -> length(x.code), outputs))

	path = shared_zeros(Int32, (max_length, V))
	code = shared_zeros(Int8, (max_length, V))

	for v in 1:V
		code[:, v] .= -1
		for i in 1:length(outputs[v])
			code[i, v] = outputs[v].code[i]
			path[i, v] = outputs[v].path[i]
		end
	end

	In = shared_rand((M, T, V), Float32(M))
	Out = shared_rand((M, V), Float32(M))

	counts = shared_zeros(Float32, (T, V))

	frequencies = shared_zeros(Int64, (V,))
	frequencies[:] = freqs

	return VectorModel(frequencies, code, path, In, Out, alpha, d, counts)
end

view(vm::VectorModel, v::Integer, s::Integer) = view(vm.In, :, s, v)

function exp_normalize!(x)
	max_x = maximum(x)
	sum_x = 0.
	for i in 1:length(x)
		x[i] = exp(x[i] - max_x)
		sum_x += x[i]
	end
	for i in 1:length(x)
		x[i] /= sum_x
	end
end

include("kahan.jl")
include("skip_gram.jl")
include("stick_breaking.jl")
include("textutil.jl")
include("gradient.jl")
include("predict.jl")
include("util.jl")

export VectorModel, gradient!
export get_gradient, apply_gradient!
export V, T, M, L
export train_vectors!, inplace_train_vectors!
export vec, closest_words
export finalize!
export save_model, read_from_file, dict_from_file, build_from_file
export disambiguate, write_dictionary
export likelihood, parallel_likelihood
export expected_pi!, expected_pi
export load_model

import ArgParse
export ArgParse

end


text:
# BAM Overlap
# ===========

struct OverlapIterator{T}
    reader::Reader{T}
    refname::String
    interval::UnitRange{Int}
end

function Base.IteratorSize(::Type{OverlapIterator{T}}) where T
    return Base.SizeUnknown()
end

function Base.eltype(::Type{OverlapIterator{T}}) where T
    return Record
end

function GenomicFeatures.eachoverlap(reader::Reader, interval::Interval)
    return GenomicFeatures.eachoverlap(reader, interval.seqname, interval.first:interval.last)
end

function GenomicFeatures.eachoverlap(reader::Reader, interval)
    return GenomicFeatures.eachoverlap(reader, convert(Interval, interval))
end

function GenomicFeatures.eachoverlap(reader::Reader, refname::AbstractString, interval::UnitRange)
    return OverlapIterator(reader, String(refname), interval)
end

# Iterator
# --------

mutable struct OverlapIteratorState
    # reference index
    refindex::Int

    # possibly overlapping chunks
    chunks::Vector{GenomicFeatures.Indexes.Chunk}

    # current chunk index
    chunkid::Int

    # pre-allocated record
    record::Record
end

function Base.iterate(iter::OverlapIterator)
    refindex = findfirst(isequal(iter.refname), iter.reader.refseqnames)
    if refindex == 0
        throw(ArgumentError("sequence name $(iter.refname) is not found in the header"))
    end
    @assert iter.reader.index !== nothing
    chunks = GenomicFeatures.Indexes.overlapchunks(iter.reader.index.index, refindex, iter.interval)
    if !isempty(chunks)
        seek(iter.reader, first(chunks).start)
    end
    state = OverlapIteratorState(refindex, chunks, 1, Record())
    return iterate(iter, state)
end

function Base.iterate(iter::OverlapIterator, state)
    while state.chunkid ≤ lastindex(state.chunks)
        chunk = state.chunks[state.chunkid]
        while BGZFStreams.virtualoffset(iter.reader.stream) < chunk.stop
            read!(iter.reader, state.record)
            c = compare_intervals(state.record, (state.refindex, iter.interval))
            if c == 0
                return copy(state.record), state
            elseif c > 0
                # no more overlapping records in this chunk since records are sorted
                break
            end
        end
        state.chunkid += 1
        if state.chunkid ≤ lastindex(state.chunks)
            seek(iter.reader, state.chunks[state.chunkid].start)
        end
    end
    return nothing
end

function compare_intervals(record::Record, interval::Tuple{Int,UnitRange{Int}})
    rid = refid(record)
    if rid < interval[1] || (rid == interval[1] && rightposition(record) < first(interval[2]))
        # strictly left
        return -1
    elseif rid > interval[1] || (rid == interval[1] && position(record) > last(interval[2]))
        # strictly right
        return +1
    else
        # overlapping
        return 0
    end
end


text:
module moo_functions

function quickavg(xs::Vector{Float64}; first::Int64 = 1, last::Int64 = length(xs))::Float64
    out = 0.0
    for k in first:last
        out += xs[k]
    end
    return out/(last - first + 1)
end

function zdt1(x::Vector{Float64})::Tuple{Float64, Float64}
    g  = 1.0 + 9.0*quickavg(x, first = 2)
    f1 = x[1]
    f2 = g * (1.0 - sqrt(f1/g))
    return (f1, f2)
end

function zdt2(x::Vector{Float64})::Tuple{Float64, Float64}
    g  = 1.0 + 9.0*quickavg(x, first = 2)
    f1 = x[1]
    f2 = g * (1.0 - (f1/g)^2)
    return (f1, f2)
end

function zdt3(x::Vector{Float64})::Tuple{Float64, Float64}
    g  = 1.0 + 9.0*quickavg(x, first = 2)
    f1 = x[1]
    f2 = g * (1.0 - sqrt(f1/g) - f1/g * sin(10*pi*f1))
    return (f1, f2)
end

function zdt4(x::Vector{Float64})::Tuple{Float64, Float64}
    g  = 1.0 + 10*(length(x)-1)
    for k in 2:length(x)
        g += x[k]^2 - 10*cos(4*pi*x[k])
    end
    f1 = x[1]
    f2 = g * (1.0 - sqrt(f1/g))
    return (f1, f2)
end

function zdt6(x::Vector{Float64})::Tuple{Float64, Float64}
    g  = 1.0 + 9.0*quickavg(x, first = 1)^0.25
    f1 = 1.0 - exp(-4.0*x[1])*sin(6*pi*x[1])^6
    f2 = g*(1.0 - (f1/g)^2)
    return (f1, f2)
end

function zdt()
    return Dict{Int8, Function}(
        1 => zdt1,
        2 => zdt2,
        3 => zdt3,
        4 => zdt4,
        6 => zdt6
    )
end

function opt_hvs()
    return Dict{Int8, Float64}(
        1 => 120.0 + 2/3,
        2 => 120.0 + 1/3,
        3 => 128.77811613069076060,
        4 => 120.0 + 2/3,
        6 => 117.51857519692037009,
    )
end

end

max_stars_repo_name: OndrejKincl/evaTeaching-python

max_issues_repo_name: OndrejKincl/evaTeaching-python

max_forks_repo_name: OndrejKincl/evaTeaching-python


text:
using GR

function points_from_image(img, npts)
  w, h = size(img)
  xpts = Float64[]
  ypts = Float64[]
  cols = Int[]
  for i in 1:npts
    x = rand() * w
    y = rand() * h
    c = img[Int64(floor(x)) + 1, Int64(floor(h-y)) + 1]
    r = ( c        & 0xff) / 255.0
    g = ((c >> 8 ) & 0xff) / 255.0
    b = ((c >> 16) & 0xff) / 255.0
    if 0.2989 * r + 0.5870 * g + 0.1140 * b > 0.8
      if rand() < 0.1
        push!(xpts, x)
        push!(ypts, y)
        push!(cols, 1)
      end
      continue
    end
    push!(xpts, x)
    push!(ypts, y)
    push!(cols, inqcolorfromrgb(r, g, b))
  end
  xpts, ypts, cols
end

w, h, img = readimage("julia_logo.png")
x, y, cols = points_from_image(img, 30_000)

setwsviewport(0, 0.24, 0, 0.16)
setwswindow(0, 1, 0, 2/3)

setviewport(0, 1, 0, 2/3)
setwindow(0, w, 0, h)

setmarkersize(2/3)
setmarkertype(GR.MARKERTYPE_SOLID_CIRCLE)

settransparency(0.5)
polymarker(x, y)

n, tri = delaunay(x, y)
for i in 1:n
  if all(cols[tri[i,:]] .!= 1)
    settransparency(0.5)
    setfillcolorind(cols[tri[i,1]])
  else
    settransparency(0.1)
    setfillcolorind(1)
  end
  fillarea(x[tri[i,:]], y[tri[i,:]])
end
updatews()


text:
using Plots

function mainAlex(NN=20,TT=3650,pb=0.01,cb=1,cr=4)

    AA = dynamicCommunicators(NN,TT,pb,cb,cr)
    degOut = totalOut(AA,NN)
    Cbroad = dynamicCentrality(AA,NN,TT)
    scatterPlot(AA,size(AA)[1],degOut,Cbroad)

   # AA2 = ringNet()
   # degOut = totalOut(AA2,21)
   # Cbroad = dynamicCentrality(AA2,21,3)
   # scatterPlot(AA2,size(AA2)[1],degOut,Cbroad)
end

function scatterPlot(AA,NN,degOut,Cbroad)
    labels = [string(x) for x in 1:NN]
    scatter(degOut,Cbroad,markersize=1,series_annotations =labels,grid=false)
    xlabel!("total deg out")
    ylabel!("broadcast centrality")
    title!(string("node num=",NN))

end

function ringNet()
    AA2 = zeros(21,21,3)
    AA2[1,7,1] =1;AA2[7,1,1] =1;AA2[2,7,1] =1;AA2[7,2,1] =1;
    AA2[21,20,1] =1;AA2[20,21,1] =1;AA2[21,19,1] =1;AA2[19,21,1] =1;
    AA2[2,8,2]=1;AA2[8,2,2]=1;AA2[3,8,2]=1;AA2[8,3,2]=1;
    AA2[20,18,2]=1;AA2[18,20,2]=1;AA2[20,17,2]=1;AA2[17,20,2]=1;
    AA2[19,16,2]=1;AA2[16,19,2]=1;
    AA2[4,10,3]=1;AA2[10,4,3]=1;AA2[5,10,3]=1;AA2[10,5,3]=1;
    AA2[18,15,3]=1;AA2[15,18,3]=1;AA2[18,14,3]=1;AA2[14,18,3]=1;
    AA2[17,14,3]=1;AA2[14,17,3]=1;AA2[17,13,3]=1;AA2[13,17,3]=1;
    AA2[16,13,3]=1;AA2[13,16,3]=1;AA2[16,12,3]=1;AA2[12,16,3]=1;
    return AA2
end

function dynamicCommunicators(NN,TT,pb,cb,cr)
    imp = [x^2 for x in 1:NN]
    #imp = 2*ones(1,NN)
    #imp[NN] = 100
    println(imp)
    println("that was imp vector now normal")
    AA = zeros(NN,NN,TT)
    #INIT DATA>>
    for ii in 1:NN
        if(rand(1)[1] <= pb)
            tmp = deleteat!(collect(1:NN),ii)
            destinationN = tmp[rand(1:end)]
            AA[ii,destinationN,1] = 1
        end

    end
    #DATA LOOP: 1st put basal edges into <present> / 2nd put response from past to now
    for tt in 2:(TT)
        #println(string("tt=",tt))
        #BASAL LOOP
        for ii in 1:NN
            if(rand(1)[1] <= pb)
                tmp = deleteat!(collect(1:NN),ii)
                destinationN = tmp[rand(1:end)]
                AA[ii,destinationN,tt] = 1#FIRE NOW
            end

        end
        #RESPONSE LOOP
        for ii in 1:NN
            msgsTo_ii = AA[:,ii,tt-1]#RESPOND TO PREVIOUS STATE (tt-1)
            if(sum(AA[:,ii,tt-1]) >= 1)

                totalImportance = sum([msgsTo_ii[x]*imp[x] for x in 1:NN])
                #println(msgsTo_ii)
                #println(imp)
                #println(sum([msgsTo_ii[x]*imp[x] for x in 1:NN]))

                r_nNumerator = totalImportance
                r_nDenominator = 1 + (findmax(imp)[1] * sum(AA[:,ii,tt-1]))
                r_next = r_nNumerator / r_nDenominator
                println(r_next)
                if(r_next >=  rand(1)[1])#if so generate cr links

                    tmp = deleteat!(collect(1:NN),ii)
                    destinationNodes = tmp[randperm(length(tmp))[1:cr]]
                    AA[ii,destinationNodes,tt] = 1#FIRE NOW
                    println("fire!")

               end
            end
        end

    end

    return AA
end

function totalOut(AA,NN)
    coltmp = zeros(NN,1)
    outDegVec = sum(sum(AA,3),2)
    coltmp[:] = outDegVec[:]
    return coltmp
end

function dynamicCentrality(AA,NN,TT)
    println(":")
    daynorm = zeros(TT,1);
    for ii in 1:TT
        Atemp = AA[:,:,ii]
        daynorm[ii] = findmax(abs(eigvals(Atemp)))[1]
    end
    alpha = 0.9*(1/findmax(daynorm)[1]);#90percent of the max
    println(alpha)
    #alpha = 0.15;
    BB = eye(NN,NN);
    for kk = TT:-1:1
        Aday = AA[:,:,kk]
        #println(Aday)
        BBtmp = inv(eye(NN,NN) - alpha*Aday)#\BB
        #BBtmp = (eye(NN,NN) - alpha*Aday) \ BB
        BBtmp = abs(BBtmp)
        #BB = BBtmp / norm(BBtmp)
        BB = (BB*BBtmp) / norm(BB)
    end
    #sum across columns
    Cbroad = sum(BB',2)
    return Cbroad
end

max_stars_repo_name: mantzaris/competitiveCommunicators

max_issues_repo_name: mantzaris/competitiveCommunicators

max_forks_repo_name: mantzaris/competitiveCommunicators


text:
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
#  Project #5: The Coupled Cluster Singles and Doubles (CCSD) Energy
#~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# J. Chem. Phys. 94, 4334, doi:10.1063/1.460620

# We follow the convention that
# i, j, k,... represent occupied orbitals, with
# a, b, c,... unoccupied.
# p, q, r,... are generic indices which may represent
#             either occupied or unoccupied
# orbitals.

###################################################
#1: Preparing the Spin-Orbital Basis Integrals
###################################################
# Translate integrals from spatial mo basis to spin-orbital basis
function orb_to_spinorb(e::Array{T,1},     # eigE of spatial-MO orbitals
                        eri_mo::Array{T,1} # ERIs in spatial-MO basis
                       ) where T <: AbstractFloat
   # 1  2  3  4  5  6  7  8       (spin-orbital no.)
   # α, β, α, β, α, β, α, β, ...  (spin component
   # χ₁,χ₁,χ₂,χ₂,χ₃,χ₃,χ₄,χ₄, ... (spatial component)
   # In RHF closed shell both α and β have the same energy: map ϵ₁ ↦ ϵ₁,ϵ₂; ϵ₂ ↦ ϵ₃,ϵ₄; ...
   # cld(1,2) and cld(2,2) → 1, cld(3,2) and cld(4,2) → 2, ...
   nao  = length(e)
   nsmo = 2nao
   eri_smo = zeros(nsmo,nsmo,nsmo,nsmo)
   @inbounds (
   for p in 1:nsmo, q in 1:nsmo, r in 1:nsmo, s in 1:nsmo
      prqs = get_4idx( cld.((p,r,q,s), 2)...) # ceil division cld(p,2) == (p+1)÷2 with p ≥ 0
      psqr = get_4idx( cld.((p,s,q,r), 2)...) # broadcast and unpack arguments
      value1 = eri_mo[prqs] * (iseven(p) == iseven(r)) * (iseven(q) == iseven(s))
      value2 = eri_mo[psqr] * (iseven(p) == iseven(s)) * (iseven(q) == iseven(r))
      eri_smo[p,q,r,s] = value1 - value2
   end )
   # The fock matrix is diagonal in the spin-MO basis
   fs = Diagonal(zeros(nsmo))
   for i in 1:nsmo
      fs.diag[i] = e[cld(i,2)]
   end
   eri_smo, fs
end

###################################################
#2: Build the Initial-Guess Cluster Amplitudes
###################################################

function initial_amplitudes(occupied,virtual,
                            fs::Diagonal{T,Array{T,1}}, # Fock matrix in spin-MO basis
                            eri_smo::Array{T,4}         # ERIs        in spin-MO basis
                            ) where T <: AbstractFloat
   nsmo = length(occupied) + length(virtual)
   t1 = zeros(nsmo,nsmo)
   t2 = zeros(nsmo,nsmo,nsmo,nsmo)

   # MP₂ guess for T₂ amplitudes
   # t[a,b,i,j] = ⟨ij||ab⟩/ (ϵi + ϵj - ϵa - ϵb)
   # E_mp2      = 1/4 * ∑ijab ⟨ij||ab⟩ t[a,b,i,j]
   emp2 = 0.0
   @inbounds (
   for i in occupied, a in virtual,
                 j in occupied, b in virtual
      t2[a,b,i,j] += eri_smo[i,j,a,b] / (fs[i,i] + fs[j,j] - fs[a,a] - fs[b,b])
      emp2 += 0.5*0.5*t2[a,b,i,j]*eri_smo[i,j,a,b]
   end )
   t1, t2, emp2
end

###################################################
#3: Calculate the CC Intermediates
###################################################

# Effective doubles
function tau_(t1::Array{T,2},
              t2::Array{T,4},a::Int,b::Int,i::Int,j::Int) where T <: AbstractFloat
   # Stanton1991 (9)
   t2[a,b,i,j] + 0.5(t1[a,i]*t1[b,j] - t1[b,i]*t1[a,j])
end
function tau(t1::Array{T,2},
              t2::Array{T,4},a::Int,b::Int,i::Int,j::Int) where T <: AbstractFloat
   # Stanton1991 (10)
   t2[a,b,i,j] + t1[a,i]*t1[b,j] - t1[b,i]*t1[a,j]
end

function updateF_W!(Fae::Array{T,2},    # one-particle CC intermediates
                    Fmi::Array{T,2},    # one-particle CC intermediates
                    Fme::Array{T,2},    # one-particle CC intermediates
                    Wmnij::Array{T,4},  # two-particle CC intermediates
                    Wabef::Array{T,4},  # two-particle CC intermediates
                    Wmbej::Array{T,4},  # two-particle CC intermediates
                    occupied,virtual,   # occ and virt spin-MO spaces
                    t1::Array{T,2},     # CC amplitudes
                    t2::Array{T,4},     # CC amplitudes
                    fs::Diagonal{T,Array{T,1}}, # Fock matrix in spin-MO basis
                    eri_smo::Array{T,4}         # ERIs        in spin-MO basis
                   ) where T <: AbstractFloat
   nsmo = length(occupied) + length(virtual)
   # Stanton1991 (3)
   @inbounds (
   for a in virtual, e in virtual
      Fae[a,e] = (1 - (a == e))*fs[a,e]
      for m in occupied
         Fae[a,e] += -0.5*fs[m,e]*t1[a,m]
         for f in virtual
            Fae[a,e] += t1[f,m]*eri_smo[m,a,f,e]
            for n in occupied
               Fae[a,e] += -0.5*tau_(t1,t2,a,f,m,n)*eri_smo[m,n,e,f]
            end
         end
      end
   end )
   # Stanton1991 (4)
   @inbounds (
   for m in occupied, i in occupied
      Fmi[m,i] = (1 - (m == i))*fs[m,i]
      for e in virtual
         Fmi[m,i] += 0.5t1[e,i]*fs[m,e]
         for n in occupied
            Fmi[m,i] += t1[e,n]*eri_smo[m,n,i,e]
            for f in virtual
               Fmi[m,i] += 0.5tau_(t1,t2,e,f,i,n)*eri_smo[m,n,e,f]
            end
         end
      end
   end )
   # Stanton1991 (5)
   @inbounds (
   for e in virtual, m in occupied
      Fme[m,e] = fs[m,e]
      for f in virtual, n in occupied
         Fme[m,e] += t1[f,n]*eri_smo[m,n,e,f]
      end
   end )
   # Stanton1991 (6)
   @inbounds (
   for m in occupied, n in occupied, i in occupied, j in occupied
      Wmnij[m,n,i,j] = eri_smo[m,n,i,j]
      for e in virtual # P_(ij)
         Wmnij[m,n,i,j] += t1[e,j]*eri_smo[m,n,i,e] -
                           t1[e,i]*eri_smo[m,n,j,e]
         for f in virtual
            Wmnij[m,n,i,j] += 0.25tau(t1,t2,e,f,i,j)*eri_smo[m,n,e,f]
         end
      end
   end )
   # Stanton1991 (7)
   @inbounds (
   for a in virtual, b in virtual, e in virtual, f in virtual
      Wabef[a,b,e,f] = eri_smo[a,b,e,f]
      for m in occupied # P_(ab)
         Wabef[a,b,e,f] += -t1[b,m]*eri_smo[a,m,e,f] +
                            t1[a,m]*eri_smo[b,m,e,f]
         for n in occupied
            Wabef[a,b,e,f] += 0.25tau(t1,t2,a,b,m,n)*eri_smo[m,n,e,f]
         end
      end
   end )
   # Stanton1991 (8)
   @inbounds (
   for b in virtual, m in occupied, e in virtual, j in occupied
      Wmbej[m,b,e,j] = eri_smo[m,b,e,j]
      for f in virtual
         Wmbej[m,b,e,j] += t1[f,j]*eri_smo[m,b,e,f]
      end
      for n in occupied
         Wmbej[m,b,e,j] += -t1[b,n]*eri_smo[m,n,e,j]
         for f in virtual
            Wmbej[m,b,e,j] += -(0.5*t2[f,b,j,n] + t1[f,j]*t1[b,n])*eri_smo[m,n,e,f]
         end
      end
   end )
   Fae, Fmi, Fme, Wmnij, Wabef, Wmbej
end

###############################################
#4: Compute the Updated Cluster Amplitudes
###############################################

# Stanton1991 (1)
function updateT1(occupied,virtual,    # occ and virt spin-MO spaces
                  t1::Array{T,2},      # CC amplitudes
                  t2::Array{T,4},      # CC amplitudes
                  fs::Diagonal{T,Array{T,1}}, # Fock matrix in spin-MO basis
                  eri_smo::Array{T,4},        # ERIs        in spin-MO basis
                  Fae::Array{T,2},     # one-particle CC intermediates
                  Fmi::Array{T,2},     # one-particle CC intermediates
                  Fme::Array{T,2}      # one-particle CC intermediates
                 ) where T <: AbstractFloat
   nsmo = length(occupied) + length(virtual)
   T1 = zeros(nsmo,nsmo)
   @inbounds (
   for a in virtual, i in occupied
      T1[a,i] = fs[i,a] # 1st RHS term (leading term in the expansion)
      for e in virtual
         T1[a,i] += t1[e,i]*Fae[a,e] # 2nd RHS term
      end
      for m in occupied
         T1[a,i] += -t1[a,m]*Fmi[m,i] # 3rd RHS term
         for e in virtual
            T1[a,i] += t2[a,e,i,m]*Fme[m,e] # 4th RHS term
            for f in virtual
               T1[a,i] += -0.5t2[e,f,i,m]*eri_smo[m,a,e,f] # 6th RHS term
            end
            for n in occupied
               T1[a,i] += -0.5t2[a,e,m,n]*eri_smo[n,m,e,i] # 7th RHS term
            end
         end
      end
      for f in virtual,n in occupied
         T1[a,i] += -t1[f,n]*eri_smo[n,a,i,f] # 5th RHS term
      end
      # Denominator energies, Stanton1991 (12)
      T1[a,i] /= fs[i,i] - fs[a,a] # LHS
   end )
   T1
end

# Stanton1991 (2)
function updateT2(occupied,virtual,    # occ and virt spin-MO spaces
                  t1::Array{T,2},      # CC amplitudes
                  t2::Array{T,4},      # CC amplitudes
                  fs::Diagonal{T,Array{T,1}}, # Fock matrix in spin-MO basis
                  eri_smo::Array{T,4},        # ERIs        in spin-MO basis
                  Fae::Array{T,2},     # one-particle CC intermediates
                  Fmi::Array{T,2},     # one-particle CC intermediates
                  Fme::Array{T,2},     # one-particle CC intermediates
                  Wmnij::Array{T,4},   # two-particle CC intermediates
                  Wabef::Array{T,4},   # two-particle CC intermediates
                  Wmbej::Array{T,4}    # two-particle CC intermediates
                 ) where T <: AbstractFloat
   nsmo = length(occupied) + length(virtual)
   T2 = zeros(nsmo,nsmo,nsmo,nsmo)
   @inbounds (
   for a in virtual, i in occupied, b in virtual, j in occupied
      T2[a,b,i,j] = eri_smo[i,j,a,b] # 1st RHS (leading term in the expansion)
      for e in virtual # 2nd RHS term
         taeij, tbeij = t2[a,e,i,j], t2[b,e,i,j]
         T2[a,b,i,j] += taeij*Fae[b,e] -
                        tbeij*Fae[a,e]   # asym. permutes b <-> a
         for m in occupied # 3rd RHS term
            T2[a,b,i,j] += -0.5taeij*t1[b,m]*Fme[m,e] +
                            0.5tbeij*t1[a,m]*Fme[m,e]   # asym. permutes b <-> a
         end
      end
      for m in occupied # 4th RHS term
         tabim, tabjm = t2[a,b,i,m], t2[a,b,j,m]
         T2[a,b,i,j] += -tabim*Fmi[m,j] +
                         tabjm*Fmi[m,i]  # asym. permutes i <-> j
         for e in virtual # 5th RHS term
            T2[a,b,i,j] += -0.5tabim*t1[e,j]*Fme[m,e] +
                            0.5tabjm*t1[e,i]*Fme[m,e]  # asym. permutes i <-> -j
         end
      end
      for e in virtual # 10th RHS term
         T2[a,b,i,j] += t1[e,i]*eri_smo[a,b,e,j] -
                        t1[e,j]*eri_smo[a,b,e,i]  # asym. permutes i <-> j
         for f in virtual # 7th RHS term
            T2[a,b,i,j] += 0.5tau(t1,t2,e,f,i,j)*Wabef[a,b,e,f]
         end
      end
      for m in occupied # 11th RHS term
         T2[a,b,i,j] += -t1[a,m]*eri_smo[m,b,i,j] +
                         t1[b,m]*eri_smo[m,a,i,j]   # asym. permutes a <-> b
         for e in virtual # 8-9th RHS terms
            taeim, taejm = t2[a,e,i,m], t2[a,e,j,m]
            tbeim, tbejm = t2[b,e,i,m], t2[b,e,j,m]
            T2[a,b,i,j] +=  taeim*Wmbej[m,b,e,j] - t1[e,i]*t1[a,m]*eri_smo[m,b,e,j] +
                           -taejm*Wmbej[m,b,e,i] + t1[e,j]*t1[a,m]*eri_smo[m,b,e,i] + # i <-> j
                           -tbeim*Wmbej[m,a,e,j] - t1[e,i]*t1[b,m]*eri_smo[m,a,e,j] + #          a <-> b
                            tbejm*Wmbej[m,a,e,i] - t1[e,j]*t1[b,m]*eri_smo[m,a,e,i]   # i <-> j, a <-> b
         end
         for n in occupied # 6th RHS term
            T2[a,b,i,j] += 0.5tau(t1,t2,a,b,m,n)*Wmnij[m,n,i,j]
         end
      end
      # Denominator energies, Stanton1991 (13)
      T2[a,b,i,j] /=  fs[i,i] + fs[j,j] - fs[a,a] - fs[b,b] # LHS
   end )
   T2
end

###########################################
#5: Check for Convergence and Iterate
###########################################

function extant_E_CCSD(occupied, virtual, fs, t1, t2, eri_smo)
   # Calculate the extant CC correlation energy
   # E_CC = ∑ia fia t[a,i] + 1/4 ∑ijab ⟨ij||ab⟩t[a,b,i,j] + 1/2 ∑ijab ⟨ij||ab⟩ t[a,i] t[b,j]
   E_CCSD = 0.0
   @inbounds (
   for a in virtual, i in occupied
      E_CCSD += fs[i,a] * t1[a,i]
      for b in virtual, j in occupied
         E_CCSD += 0.25eri_smo[i,j,a,b] * t2[a,b,i,j] +
                         0.5eri_smo[i,j,a,b] * t1[a,i] * t1[b,j]
      end
   end )
   E_CCSD
end

function ccsd(e, C, eri, mol; verbose=true)

   nao      = size(C)[1]
   nsmo     = 2nao
   nelec    = mol.nelectron
   occupied =       1:nelec
   virtual  = nelec+1:nsmo
   eri_mo = ao2mo_smart(nao,C,eri,mol)

   eri_smo, fs = orb_to_spinorb(e,eri_mo)

   t1,t2, E_MP2 = initial_amplitudes(occupied,virtual,fs,eri_smo)

   E_CCSD  = 0.0
   ΔE_CCSD = 1.0
   iteri   = 0
   itermax = 60
   Fae = zeros(nsmo,nsmo)
   Fmi = zeros(nsmo,nsmo)
   Fme = zeros(nsmo,nsmo)
   Wmnij = zeros(nsmo,nsmo,nsmo,nsmo)
   Wabef = zeros(nsmo,nsmo,nsmo,nsmo)
   Wmbej = zeros(nsmo,nsmo,nsmo,nsmo)
   if verbose
      println(" Iter     Ecorr(CCSD)        ΔE(CCSD)")
      println("-----   --------------   --------------")
   end
   while ΔE_CCSD > 1e-9 && iteri < itermax
      iteri += 1
      Fae,Fmi,Fme,Wmnij,Wabef,Wmbej = updateF_W!(Fae,Fmi,Fme,Wmnij,Wabef,Wmbej,occupied,virtual,t1,t2,fs,eri_smo)
      t1 = updateT1(occupied,virtual,t1,t2,fs,eri_smo,Fae,Fmi,Fme)
      t2 = updateT2(occupied,virtual,t1,t2,fs,eri_smo,Fae,Fmi,Fme,Wmnij,Wabef,Wmbej)
      E_CCSD, oldE_CCSD = extant_E_CCSD(occupied,virtual,fs,t1,t2,eri_smo), E_CCSD
      ΔE_CCSD = abs(E_CCSD - oldE_CCSD)
      verbose && printfmt("{1:5d}   {2:13.10f}   {3:13.10f}\n",iteri,E_CCSD,ΔE_CCSD)
   end
   @assert iteri < itermax "NOT CONVERGED!!"
   E_CCSD, E_MP2
end


text:
#  ------------------------------------------------------------------------------------------
#  Copyright (c) Microsoft Corporation. All rights reserved.
#  Licensed under the MIT License (MIT). See LICENSE in the repo root for license information.
#  ------------------------------------------------------------------------------------------

module TestCore

    # Environment variables
    ENV["CREDENTIALS"] = ""
    ENV["PARAMETERS"] = "params.json"

    if haskey(ENV, "AZ_BATCH_TASK_WORKING_DIR")
        delete!(ENV, "AZ_BATCH_TASK_WORKING_DIR")
    end

    using AzureClusterlessHPC, PyCall, Test, SyntaxTree, Random, Serialization

    function __init__()

        global test_suite = get(ENV, "test_suite", "core") # "all", "core", "pyinterface"

        # AzureClusterlessHPC core
        global azureclusterlesshpc_core = ["core/test_batch_controller.jl",
                        "core/test_batch_macros.jl",
                        "core/test_batch_futures.jl"]

        # AzureClusterlessHPC python interface
        global azureclusterlesshpc_pyinterface = ["pyinterface/test_pyinterface.jl"]

    end

    function runtests()
        if test_suite == "all" || test_suite == "core"
            for t=azureclusterlesshpc_core
                @testset "Test $t" begin
                    @time include(t)
                end
            end
        end

        if test_suite == "all" || test_suite == "pyinterface"
            for t=azureclusterlesshpc_pyinterface
                @testset  "Test $t" begin
                    @time include(t)
                end
            end
        end
    end

end


text:
function plotoceans(ocns)
    plt = plot()
    for ocn in ocns
        OceanBasins.isocean(ocn) && plot!(plt, ocn.lon, ocn.lat, label=OceanBasins.name(ocn))
    end
    return plt
end

function filloceans(ocns)
    plt = plot()
    for ocn in ocns
        if OceanBasins.isocean(ocn)
            s = Shape(ocn.lon, ocn.lat)
            label = shorten_name(OceanBasins.name(ocn))
            plot!(plt, s, label=label, alpha=0.5)
        end
    end
    return plt
end

function shorten_name(name)
    name = replace(name, " Ocean" => "")
    name = replace(name, "South" => "S")
    name = replace(name, "North" => "N")
    name = replace(name, "eastern part" => "EP")
    name = replace(name, "western part" => "WP")
    name = replace(name, "Atlantic" => "ATL")
    name = replace(name, "Pacific" => "PAC")
    name = replace(name, "Indian" => "IND")
    name = replace(name, "Arctic" => "ARC")
    return name
end

#=
GMT plot for ReadMe
=#
using Libdl
push!(Libdl.DL_LOAD_PATH, "/usr/local/Cellar/gmt/6.0.0_5/lib")
using GMT
coast(region=:d, proj=:Robinson, frame=:g, res=:crude, area=10000, land=:lemonchiffon1, water=:lightsteelblue1, figsize=12)
using OceanBasins
const OCEANS = oceanpolygons()
N = 2000
lons = 360rand(N) .- 180
lats = 180rand(N) .- 90
isocns = [ispacific, isatlantic, isindian, isarctic, ismediterranean, isantarctic]
colors = [sum(iocn * isocns[iocn](lat,lon,OCEANS) for iocn in 1:length(isocns)) for (lat,lon) in zip(lats,lons)]
scatter!(lons, lats, title="Which ocean basin?", marker=:c, size=0.1, zcolor=colors, show=1, savefig="demo.png")


text:
export ft_getlatencytimer

function ft_getlatencytimer(ft_handle::Culong)
  timer = Ref{Cuchar}()
  ft_status = ccall((:FT_GetLatencyTimer, d2xx),
                     Cuint,
                     (Culong, Ref{Cuchar}),
                     ft_handle, timer)
  checkstatus(ft_status)
  return timer[]
end


text:
include("bnf_parse.jl")
include("bnf_transform.jl")
include("bnf_build.jl")

function bnf_rules(bnf::IO, startvariable::AbstractString, syntax::Symbol=:ebnf, addwhitespace::Bool=true, rulenameprefix="")
	ast = parse_bnf(bnf::IO, syntax)
	transform_bnf_ast(ast, startvariable)
	transform_ast(ast)
	build_bnf_rules(ast, addwhitespace, rulenameprefix)
end

function bnf_generator(io::IO, genname::Symbol, bnf::IO, startvariable::AbstractString, syntax::Symbol=:ebnf, addwhitespace::Bool=true)
	rules = bnf_rules(bnf, startvariable, syntax, addwhitespace)
	description = "string accepted by BNF starting with variable " * escape_string(startvariable)
	output_generator(io, genname, description, rules)
end

bnf_generator(genname::Symbol, bnf::IO, startvariable::AbstractString, syntax::Symbol=:ebnf, addwhitespace::Bool=true) = include_generator(genname, bnf_generator, bnf, startvariable, syntax, addwhitespace)


text:
using APES
using Base.Test
using Lint
using TypeCheck

# Run tests
include("APES.jl")
include("APES-speed.jl")

# Check code
lintpkg( "APES" )

# checkreturntypes(APES)
# println()
# checklooptypes(APES)
# println()
# checkmethodcalls(APES)


text:
module SampleJuliaProj

# Write your package code here.

using Combinatorics
export find_subset
greet() = "Hello World"
greet_random() = "Hello World Random"
# Write your package code here.
export func

"""
    func(x)

Returns double the number `x` plus `1`.
"""
func(x) = 2x + 1

"""
	find_subset(n::Int64,p,q,RV::AbstractString="X")
```@docs
Given `i` and `j` compute `Κ ⊆ \U1d4a9 \\{i,j}`; i.e., all non exclusive subsets. `i,j` can also be empty (i.e. []), in which case the non-empty superset gets listed. An optional prefix can be added (default is `X`).
#### Examples
```

Given `i` and `j` compute `Κ ⊆ \U1d4a9 \\{i,j}`; i.e., all non exclusive subsets. `i,j` can also be empty (i.e. []), in which case the non-empty superset gets listed. An optional prefix can be added (default is `X`).
#### Examples
```julia-repl
julia>find_subset(4,1,3,"X")
["" "X1" "X1,X2" "X1,X2,X5" "X1,X5" "X2" "X2,X5" "X5"]

julia>find_subset(4,1,3,"")
["" "1" "1,2" "1,2,5" "1,5" "2" "2,5" "5"]

julia> find_subset(5,4,3,"🍒")
["" "🍒3" "🍒3,🍒4" "🍒4"]

julia> find_subset(5,[],[],"🍓")
["🍓1" "🍓1,🍓2" "🍓1,🍓2,🍓3" "🍓1,🍓3" "🍓2" "🍓2,🍓3" "🍓3"]
```
"""
function find_subset(n::Int64,p,q,RV::AbstractString="X")
	s="$(RV)" .* string.(1:n)
	s1="$(RV)" .* string(p)
	s2="$(RV)" .* string(q)
	PowerSet = collect(powerset(s,1,n)) # Skip empty set
	setXk=[]
	for kk = 1:Int64(exp2(n)-1)
		subSet1 = vcat(PowerSet[kk,:][1])
		subSet2 = setdiff(subSet1,[s1,s2])
		setXk=[setXk; join(subSet2,",")]
	end

	return sort(unique(setXk))
end

end


text:
# SAM Reader
# ==========

type SAMReader <: Bio.IO.AbstractParser
    state::Ragel.State
    header::SAMHeader

    function SAMReader(input::BufferedInputStream)
        return new(Ragel.State(samparser_start, input), SAMHeader())
    end
end

function header(reader::SAMReader)
    return reader.header
end

function Base.eltype(::Type{SAMReader})
    return SAMRecord
end

function Base.eof(reader::SAMReader)
    return eof(reader.state.stream)
end

function Base.close(reader::SAMReader)
    close(reader.state.stream)
end

function Base.open(input::BufferedInputStream, ::Type{SAM})
    return SAMReader(input)
end


text:
using Random, Plots
include_ddp()

# Model
include_model("satellite")

# Time
T = 33
h = 0.1
tf = h * (T - 1)

# Reference trajectory
# https://en.wikipedia.org/wiki/Viviani%27s_curve
# # https://mathworld.wolfram.com/VivianisCurve.html
# t = range(-2.0 * π, stop = 0.0 * π, length = T)
# a = 0.05
# xf = a * (1.0 .+ cos.(t))
# yf = a * sin.(t)
# zf = 2.0 * a * sin.(0.5 .* t)
#
# plot(xf, zf, aspect_ratio = :equal)
# plot(xf, yf, aspect_ratio = :equal)
# plot(yf, zf, aspect_ratio = :equal)
# plot(xf, yf, zf, aspect_ratio = :equal)
# pf = [RotZ(0.0 * π) * [xf[t]; yf[t]; zf[t]] for t = 1:T]

θ = range(0.0, stop = 0.5 * π, length = 11)
p1 = 0.1 * cos.(θ)
p2 = 0.1 * sin.(θ)
plot(p1, p2)

ref1 = [[p1[t]; 0.0; p2[t]] for t = 1:11]
ref2 = [[0.0; p2[t]; p1[t]] for t = 1:11]
ref3 = [[p2[t]; p1[t]; 0.0] for t = 1:11]
pf = [ref1..., ref2..., ref3...]
xf = [pt[1] for pt in pf]
yf = [pt[2] for pt in pf]
zf = [pt[3] for pt in pf]

# Initial conditions, controls, disturbances
mrp1 = MRP(RotXYZ(0.0, 0.0, 0.0))
mrpT = MRP(RotXYZ(0.0, 0.0, 0.0))

x1 = [mrp1.x, mrp1.y, mrp1.z, 0.0, 0.0, 0.0]
xT = [mrpT.x, mrpT.y, mrpT.z, 0.0, 0.0, 0.0]

ū = [1.0e-3 * rand(model.m) for t = 1:T-1]
w = [zeros(model.d) for t = 1:T-1]

# Rollout
x̄ = rollout(model, x1, ū, w, h, T)

obj = StageCosts([NonlinearCost() for t = 1:T], T)

function g(obj::StageCosts, x, u, t)
	T = obj.T
	ω = view(x, 4:6)

	J = 0.0

	# tracking
	p = kinematics(model, x)
	J += transpose(p - pf[t]) * Diagonal(1000.0 * [1.0; 1.0; 1.0]) * (p - pf[t])

	# control limits
	if t < T
		# energy
		J += 0.5 * transpose(ω) * model.J * ω
		J += 1.0e-1 * transpose(u[1:3]) * u[1:3]
	end

	return J / T
end

# Constraints
p_con = [t == T ? model.n : (t == 11 || t == 22) ? 3 : 0 for t = 1:T]

info_t = Dict(:x11 => ref1[end], :x22 => ref2[end])
info_T = Dict(:xT => xT)

con_set = [StageConstraint(p_con[t], t < T ? info_t : info_T) for t = 1:T]

function c!(c, cons::StageConstraints, x, u, t)
	T = cons.T
	n = model.n

	if t == 11
		c[1:3] .= kinematics(model, x) - cons.con[t].info[:x11]
	end

	if t == 22
		c[1:3] .= kinematics(model, x) - cons.con[t].info[:x22]
	end

	if t == T
		c[1:n] .= x - cons.con[T].info[:xT]
	end
end

prob = problem_data(model, obj, con_set, copy(x̄), copy(ū), w, h, T)

# Solve
@time constrained_ddp_solve!(prob,
    max_iter = 1000, max_al_iter = 10,
	ρ_init = 1.0, ρ_scale = 10.0,
	con_tol = 1.0e-3,
	cache = false)

x, u = current_trajectory(prob)
x̄, ū = nominal_trajectory(prob)

# Visualize
include(joinpath(pwd(), "models/visualize.jl"))
vis = Visualizer()
render(vis)
visualize!(vis, model, x, Δt = h)

p = [kinematics(model, view(xt, 1:3)) for xt in x]
px = [pt[1] for pt in p]
py = [pt[2] for pt in p]
pz = [pt[3] for pt in p]

plot(yf[1:33], zf[1:33])
plot!([pt[2] for pt in p], [pt[3] for pt in p])

plot(xf[1:33], zf[1:33])
plot(xf[1:33], yf[1:33])
plot(xf)
plot(yf)
plot(zf)

plot(xf, yf, zf, aspect_ratio = :equal)
plot!(px, py, pz, aspect_ratio = :equal)

max_issues_repo_name: thowell/DirectMotionPlanning


text:
"""
    ForceIntensityModule

Module  to manage  distributed force intensity.
"""
module ForceIntensityModule

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
import ..VectorCacheModule: VectorCache, updateretrieve!, settime!

"""
    ForceIntensity{T<:Number, F<:Function}

Distributed force (force intensity) type.

The force intensity class. The physical units are force per unit volume,
where volume depends on to which manifold the force is applied:
- force/length^3 (when applied to a 3-D solid),
- force/length^2 (when applied to a surface),
- force/length^1 (when applied along a curve), or
- force/length^0 (when applied at a point).

Signature of the function to compute the value of the force  at any
given point `XYZ`, using the columns of the Jacobian matrix of the
element, `tangents`, and if necessary  also the finite element label,
`fe_label`:
```
getforce!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
```
The buffer `forceout` is filled with the value  of the force. The vector
`forceout` is also returned for convenience.
"""
struct ForceIntensity{T<:Number, F<:Function}
    cache::VectorCache{T, F}   # vector cache  where the current value of the force can be retrieved
end

"""
    ForceIntensity(::Type{T}, ndofn::FInt, computeforce!::F) where {T<:Number, F<:Function}

Construct force intensity when the function to compute the intensity
vector is given.

This constructor is intended for *time-independent* vector caches.

# Arguments
- `T` = the type of the elements of the force vector, typically floating-point or complex floating-point numbers,
- `ndofn` = number of elements of the force vector (the length of the force vector),
- `computeforce!` = callback function.
The function `computeforce!` needs to have a signature of
```
function computeforce!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    Calculate the force  and copy it into the buffer....
    return forceout
end
```
and it needs to  fill in the buffer `forceout` with the current force at the
location `XYZ`, using if appropriate the information supplied in the Jacobian
matrix `tangents`, and the label of the finite element, `fe_label`.
"""
function ForceIntensity(::Type{T}, ndofn::FInt, computeforce!::F) where {T<:Number, F<:Function}
    # Allocate the buffer to be ready for the first call
    return ForceIntensity(VectorCache(T, ndofn, computeforce!));
end

"""
    ForceIntensity(::Type{T}, ndofn::FInt, computeforce!::F, time::FFlt) where {T<:Number, F<:Function}

Construct force intensity when the function to compute the intensity
vector is given.

This constructor is intended for time-dependent force intensity caches.

# Arguments
- `T` = the type of the elements of the force vector, typically floating-point or complex floating-point numbers,
- `ndofn` = number of elements of the force vector (the length of the force vector),
- `computeforce!` = callback function,
- `time` = initial time.
The function `computeforce!` needs to have a signature of
```
function computeforce!(forceout::FFltVec, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt; time::FFlt = 0.0)
	Calculate the force  and copy it into the buffer....
	return forceout
end
```
and it needs to  fill in the buffer `forceout` with the current force at
the location `XYZ`, using if appropriate the information supplied in the
Jacobian matrix `tangents`, and the label of the finite element,
`fe_label`. The initial `time` is given.

The time needs to be set with `settime!` before calling `updateforce!`
as follows:
```
XYZ = reshape([0.0, 0.0], 2, 1)
tangents = reshape([0.0, 1.0], 2, 1)
fe_label = 0
setvector!(v, XYZ, tangents, fe_label; time::FFlt = 0.0) = begin
    return (time < 5.0 ?  v .= [10.0] : v .= [0.0])
end
vector = [10.0]
fi = ForceIntensity(FFlt, length(vector), setvector!, 0.0)
v = updateforce!(fi, XYZ, tangents, fe_label)
@test v == [10.0]
settime!(fi, 6.0)
v = updateforce!(fi, XYZ, tangents, fe_label)
@test v == [0.0]
```
"""
function ForceIntensity(::Type{T}, ndofn::FInt, computeforce!::F, time::FFlt) where {T<:Number, F<:Function}
    # Allocate the buffer to be ready for the first call
    return ForceIntensity(VectorCache(T, ndofn, computeforce!, time));
end

"""
    ForceIntensity(force::FVec{T}) where {T<:Number}

Construct force intensity when the constant `force` vector is given.
"""
function ForceIntensity(force::FVec{T}) where {T<:Number}
    return ForceIntensity(VectorCache(deepcopy(force)));
end

"""
    ForceIntensity(force::T) where {T<:Number}

Construct force intensity when the force is given as a scalar value.

The dimension of the force vector in this case is 1.
"""
function ForceIntensity(force::T) where {T<:Number}
    return ForceIntensity(T[force]);
end

"""
    updateforce!(self::ForceIntensity, ndofn::FInt, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)

Update the force intensity vector.

Returns a vector (stored in the cache `self.cache`).
"""
function updateforce!(self::ForceIntensity, XYZ::FFltMat, tangents::FFltMat, fe_label::FInt)
    return updateretrieve!(self.cache, XYZ, tangents, fe_label)
end

"""
    settime!(self::ForceIntensity, time::FFlt)

Set the current time for the force intensity.
"""
function settime!(self::ForceIntensity, time::FFlt)
	settime!(self.cache, time)
	return self
end

end


text:
# test file for low-level APIs

# constants
const width = 10
const height = 5
const pixelNum = width * height
const labelNum = 7

# grid graph
gridgraph = GCoptimizationGridGraph(width, height, labelNum)

for i = 0:pixelNum-1, l = 0:labelNum-1
    if i < 25
        l == 0 ? setDataCost(gridgraph, i, l, 0) : setDataCost(gridgraph, i, l, 10)
    else
        l == 5 ? setDataCost(gridgraph, i, l, 0) : setDataCost(gridgraph, i, l, 10)
    end
end

for l1 = 0:labelNum-1, l2 = 0:labelNum-1
    cost = (l1-l2)^2 <= 4 ? (l1-l2)^2 : 4
    setSmoothCost(gridgraph, l1, l2, cost)
end

@test compute_energy(gridgraph) == 250
@test expansion(gridgraph) == 44
@test expansion(gridgraph, 2) == 44

# data array & smooth array
data = Matrix{Int}(labelNum, pixelNum)
for i = 1:pixelNum, l = 1:labelNum
    if i <= 25
        if l == 1
            data[l,i] = 0
        else
            data[l,i] = 10
        end
    else
        if l == 6
            data[l,i] = 0
        else
            data[l,i] = 10
        end
    end
end
data = [data...;]

smooth = Matrix{Int}(labelNum, labelNum)
for l1 = 1:labelNum, l2 = 1:labelNum
    smooth[l1,l2] = (l1-l2)*(l1-l2) <= 4  ? (l1-l2)*(l1-l2) : 4;
end
smooth = [smooth...;]

gridgraphArray = GCoptimizationGridGraph(width, height, labelNum)
setDataCost(gridgraphArray, data)
setSmoothCost(gridgraphArray, smooth)

@test compute_energy(gridgraphArray) == 250
@test swap(gridgraphArray) == 44

# H & V
H = [i+(i+1)%3 for i in 0:pixelNum-1]
V = [i*(i+width)%7 for i in 0:pixelNum-1]

gridgraphHV = GCoptimizationGridGraph(width, height, labelNum)
setDataCost(gridgraphHV, data)
setSmoothCostVH(gridgraphHV, smooth, V, H)

@test compute_energy(gridgraphHV) == 250
@test swap(gridgraphHV, 2) == 170

# pass julia function to C++
function dataFn(p::GCOSiteID, l::GCOLabelID)::GCOEnergyTermType
    if p < 25
        return l == 0 ? 0 : 10
    else
        return l == 5 ? 0 : 10
    end
end

smoothFn(p1::GCOSiteID, p2::GCOSiteID, l1::GCOLabelID, l2::GCOLabelID)::GCOEnergyTermType = (l1-l2)*(l1-l2) <= 4 ? (l1-l2)*(l1-l2) : 4

gridgraphFn = GCoptimizationGridGraph(width, height, labelNum)

setDataCost(gridgraphFn, dataFn)
setSmoothCost(gridgraphFn, smoothFn)

@test compute_energy(gridgraphFn) == 250
@test swap(gridgraphFn, 2) == 44

# general graph
generalgraph = GCoptimizationGeneralGraph(pixelNum, labelNum)

setDataCost(generalgraph, data)
setSmoothCost(generalgraph, smooth)

for y = 0:height-1, x = 1:width-1
    setNeighbors(generalgraph, x+y*width, x-1+y*width)
end

for y = 1:height-1, x = 0:width-1
    setNeighbors(generalgraph, x+y*width, x+(y-1)*width)
end

@test compute_energy(generalgraph) == 250
@test expansion(generalgraph) == 44

# with non-one weight
generalgraphWeight = GCoptimizationGeneralGraph(pixelNum, labelNum)

setDataCost(generalgraphWeight, data)
setSmoothCost(generalgraphWeight, smooth)

for y = 0:height-1, x = 1:width-1
    p1 = x - 1 + y*width
    p2 = x + y*width
    setNeighbors(generalgraphWeight, p1, p2, p1+p2)
end

for y = 1:height-1, x = 0:width-1
    p1 = x + (y-1)*width
    p2 = x + y*width
    setNeighbors(generalgraphWeight, p1, p2, p1*p2)
end

@test compute_energy(generalgraphWeight) == 250
@test expansion(generalgraphWeight) == 244

# miscellaneous
@test numSites(gridgraph) == width * height
@test numLabels(gridgraph) == labelNum

@test whatLabel(gridgraph, 24) == 0
@test whatLabel(gridgraph, 25) == 5
labeling = Vector{Cint}(4)
whatLabel(gridgraph, 23, 4, labeling)
@test labeling == [0, 0, 5, 5]

@test giveDataEnergy(gridgraph) == 0
@test giveSmoothEnergy(gridgraph) == 44
@test giveLabelEnergy(gridgraph) == 0

alpha_expansion(gridgraph, 0)
alpha_beta_swap(gridgraph, 0, 5)

setVerbosity(gridgraph, 2)

# add test for setLabel
setLabel(gridgraph, 10, 2)
whatLabel(gridgraph, 10) == 2

# test setDataCost(gco::Cxx.CppPtr,  dataCostFuncExtra::Function, extraData::Ptr{Void})
#      setSmoothCost(gco::Cxx.CppPtr, smoothCostFuncExtra::Function, extraData::Ptr{Void})
immutable ForDataFn
    data::Vector{Cint}
    numLab::Cint
end

function dataFn(p::GCOSiteID, l::GCOLabelID, data::Ptr{Void})::GCOEnergyTermType
    myData = unsafe_pointer_to_objref(data)
    numLab = myData.numLab
    return myData.data[p*numLab+l+1]
end

immutable ForSmoothFn
    data::Vector{Cint}
    numLab::Cint
end

function smoothFn(p1::GCOSiteID, p2::GCOSiteID, l1::GCOLabelID, l2::GCOLabelID, data::Ptr{Void})::GCOEnergyTermType
    mySmooth = unsafe_pointer_to_objref(data)
    numLab = mySmooth.numLab
    return mySmooth.data[l1+l2*numLab+1]
end

gco = GCoptimizationGridGraph(width, height, labelNum)

toDataFn = ForDataFn(data, labelNum)
setDataCost(gco, dataFn, pointer_from_objref(toDataFn))

toSmoothFn = ForSmoothFn(smooth, labelNum)
setSmoothCost(gco, smoothFn, pointer_from_objref(toSmoothFn))

@test compute_energy(gco) == 250
@test expansion(gco) == 44


text:
import Distributions: pdf, rand, logpdf

struct LNRC{T1,T2,T3} <: ContinuousUnivariateDistribution
    μ::T1
    σ::T2
    ϕ::T3
    c::Int
end

Broadcast.broadcastable(x::LNRC) = Ref(x)

LNRC(;μ, σ, ϕ, c) = LNRC(μ, σ, ϕ, c)

function rand(dist::LNRC)
    @unpack μ,σ,ϕ,c = dist
    x = @. rand(LogNormal(μ, σ)) + ϕ
    rt,resp = findmin(x)
    return rt
end

function rand_sim(dists)
    total_rt = 0.0
    resps = fill(0, length(dists))
    for (i,d) in enumerate(dists)
        @unpack μ,σ,ϕ = d
        x = @. rand(LogNormal(μ, σ)) + ϕ
        rt,resp = findmin(x)
        total_rt += rt
        resps[i] = resp
    end
    return resps,total_rt
end

rand(dist::LNRC, N::Int) = [rand(dist) for i in 1:N]

function logpdf(d::T, t::Float64) where {T <: LNRC}
    @unpack μ,σ,ϕ,c = d
    LL = 0.0
    for (i,m) in enumerate(μ)
        if i == c
            LL += logpdf(LogNormal(m, σ), t - ϕ)
        else
            LL += logccdf(LogNormal(m, σ), t - ϕ)
        end
    end
    return LL
end

function logpdf(d::LNRC{T1,T2,Vector{T3}}, t::Float64) where {T1,T2,T3}
    @unpack μ,σ,ϕ,c = d
    LL = 0.0
    for (i,m) in enumerate(μ)
        if i == c
            LL += logpdf(LogNormal(m, σ), t - ϕ[i])
        else
            LL += logccdf(LogNormal(m, σ), t - ϕ[i])
        end
    end
    return LL
end

function pdf(d::T, t::Float64) where {T <: LNRC}
    @unpack μ,σ,ϕ,c = d
    density = 1.0
    for (i,m) in enumerate(μ)
        if i == c
            density *= pdf(LogNormal(m, σ), t - ϕ)
        else
            density *= (1 - cdf(LogNormal(m, σ), t - ϕ))
        end
    end
    return density
end


text:
struct WanMerweWeightingParameters <: AbstractWeightingParameters
    α::Float64
    β::Float64
    κ::Float64
end

function WanMerweWeightingParameters()
    WanMerweWeightingParameters(1e-3, 2, 0)
end

struct MeanSetWeightingParameters <: AbstractWeightingParameters
    ω₀::Float64
end

function MeanSetWeightingParameters()
    MeanSetWeightingParameters(1/3)
end

struct GaussSetWeightingParameters <: AbstractWeightingParameters
    κ::Float64
end

function GaussSetWeightingParameters()
    GaussSetWeightingParameters(3)
end

struct ScaledSetWeightingParameters <: AbstractWeightingParameters
    α::Float64
    β::Float64
    κ::Float64
end

function ScaledSetWeightingParameters()
    ScaledSetWeightingParameters(1e-3, 2.0, 1.0)
end

lambda(weight_params::WanMerweWeightingParameters, L) = weight_params.α^2 * (L + weight_params.κ) - L

function calc_mean_weights(weight_params::WanMerweWeightingParameters, num_states)
    λ = lambda(weight_params, num_states)
    weight_0 = λ / (num_states + λ)
    weight_i = 1 / (2 * (num_states + λ))
    weight_0, weight_i
end

function calc_cov_weights(weight_params::WanMerweWeightingParameters, num_states)
    weight_0, weight_i = calc_mean_weights(weight_params, num_states)
    weight_0 + 1 - weight_params.α^2 + weight_params.β, weight_i
end

function calc_cholesky_weight(weight_params::WanMerweWeightingParameters, num_states::Real)
    num_states + lambda(weight_params, num_states)
end

function calc_mean_weights(weight_params::ScaledSetWeightingParameters, num_states)
    weight_0 = (weight_params.α^2 * weight_params.κ - num_states) / (weight_params.α^2 * weight_params.κ)
    weight_i = 1 / (2 * weight_params.α^2 * weight_params.κ)
    weight_0, weight_i
end

function calc_cov_weights(weight_params::ScaledSetWeightingParameters, num_states)
    weight_0, weight_i = calc_mean_weights(weight_params, num_states)
    weight_0 + 1 - weight_params.α^2 + weight_params.β, weight_i
end

function calc_cholesky_weight(weight_params::ScaledSetWeightingParameters, num_states::Real)
    weight_params.α^2 * weight_params.κ
end

function calc_mean_weights(weight_params::MeanSetWeightingParameters, num_states)
    weight_params.ω₀, (1 - weight_params.ω₀) / (2 * num_states)
end

calc_cov_weights(weight_params::AbstractWeightingParameters, num_states) =
    calc_mean_weights(weight_params, num_states)

function calc_cholesky_weight(weight_params::MeanSetWeightingParameters, num_states::Real)
    num_states / (1 - weight_params.ω₀)
end

function calc_mean_weights(weight_params::GaussSetWeightingParameters, num_states)
    1 - num_states / weight_params.κ, 1 / (2 * weight_params.κ)
end

function calc_cholesky_weight(weight_params::GaussSetWeightingParameters, num_states::Real)
    weight_params.κ
end

struct SPTimeUpdate{X,P,O} <: AbstractTimeUpdate{X,P}
    state::X
    covariance::P
    χ::O
end

struct SPMeasurementUpdate{X,P,O,I,IC,K} <: AbstractMeasurementUpdate{X,P}
    state::X
    covariance::P
    𝓨::O
    innovation::I
    innovation_covariance::IC
    kalman_gain::K
end

abstract type AbstractSigmaPoints{T} <: AbstractMatrix{T} end

struct SigmaPoints{T, V <: AbstractVector{T}, L <: LowerTriangular{T}, W <: AbstractWeightingParameters} <: AbstractSigmaPoints{T}
    x0::V
    P_chol::L
    weight_params::W
    SigmaPoints{T, V, L, W}(x0, P_chol, weight_params) where {T<:Real, V<:AbstractVector{T}, L<:LowerTriangular{T}, W<:AbstractWeightingParameters} =
        size(x0, 1) == size(P_chol, 1) == size(P_chol, 2) ?
        new{T, V, L, W}(x0, P_chol, weight_params) :
        error("The length of the first dimension must be equal to the size of P_chol")
end

SigmaPoints(x0::V, P_chol::Cholesky{T}, weight_params::W) where {T<:Real, V<:AbstractVector{T}, W<:AbstractWeightingParameters} =
    SigmaPoints{T, V, typeof(P_chol.L), W}(x0, P_chol.L, weight_params)
SigmaPoints(x0::V, P_chol::L, weight_params::W) where {T<:Real, V<:AbstractVector{T}, L<:LowerTriangular{T}, W<:AbstractWeightingParameters} =
    SigmaPoints{T, V, L, W}(x0, P_chol, weight_params)

function calc_sigma_points(
    x::V,
    P::AbstractMatrix{T},
    weight_params::W
) where {T, V<:AbstractVector{T}, W<:AbstractWeightingParameters}
    weight = calc_cholesky_weight(weight_params, P)
    P_chol = cholesky(Hermitian(P .* weight, :L))
    SigmaPoints{T, V, typeof(P_chol.L), W}(x, P_chol.L, weight_params)
end

function calc_sigma_points(
    x::V,
    P::Cholesky{T},
    weight_params::W
) where {T, V<:AbstractVector{T}, W<:AbstractWeightingParameters}
    weight = calc_cholesky_weight(weight_params, P)
    weighted_P_chol = P.L * sqrt(weight)
    SigmaPoints{T, V, typeof(weighted_P_chol), W}(x, weighted_P_chol, weight_params)
end

function calc_sigma_points!(
    P_chol_temp::M,
    x::V,
    P::M,
    weight_params::W
) where {T, V<:AbstractVector{T}, M<:AbstractMatrix{T}, W<:AbstractWeightingParameters}
    weight = calc_cholesky_weight(weight_params, P)
    P_chol_temp .= P .* weight
    P_chol = cholesky!(Hermitian(P_chol_temp, :L))
    SigmaPoints{T, V, typeof(P_chol.L), W}(x, P_chol.L, weight_params)
end

function calc_sigma_points!(
    P_chol_temp::M,
    x::V,
    P::Cholesky{T},
    weight_params::W
) where {T, V<:AbstractVector{T}, M<:AbstractMatrix{T}, W<:AbstractWeightingParameters}
    weight = calc_cholesky_weight(weight_params, P)
    P_chol_temp .= (P.uplo === 'U' ? transpose(P.U) : P.L) .* sqrt(weight)
    SigmaPoints(x, LowerTriangular(P_chol_temp), weight_params)
end

struct TransformedSigmaPoints{T, V <: AbstractVector{T}, M <: AbstractMatrix{T}, W <: AbstractWeightingParameters} <: AbstractSigmaPoints{T}
    x0::V
    xi::M
    weight_params::W
    TransformedSigmaPoints{T, V, M, W}(x0, xi, weight_params) where {T<:Real, V<:AbstractVector{T}, M<:AbstractMatrix{T}, W<:AbstractWeightingParameters} =
        size(x0, 1) == size(xi, 1) ?
        new{T, V, M, W}(x0, xi, weight_params) :
        error("The length of the first dimension must be the same for all inputs")
end

TransformedSigmaPoints(x0::V, xi::M, weight_params::W) where {T<:Real, V<:AbstractVector{T}, M<:AbstractMatrix{T}, W<:AbstractWeightingParameters} =
    TransformedSigmaPoints{T, V, M, W}(x0, xi, weight_params)

function transform(F, χ::SigmaPoints{T}) where T
    𝓨_x0 = F(χ.x0)
    num_x = length(χ.x0)
    𝓨_xi = Matrix{T}(undef, length(𝓨_x0), 2 * length(χ.x0))
    xi_temp = copy(χ.x0)
    @inbounds for i = length(χ.x0):-1:1
        xi_temp[i:num_x] .= @view(χ.x0[i:num_x]) .+ @view(χ.P_chol.data[i:num_x, i])
        𝓨_xi[:, i] = F(xi_temp)
        xi_temp[i:num_x] .= @view(χ.x0[i:num_x]) .- @view(χ.P_chol.data[i:num_x, i])
        𝓨_xi[:, i + length(χ.x0)] = F(xi_temp)
    end
    TransformedSigmaPoints(𝓨_x0, 𝓨_xi, χ.weight_params)
end

function transform!(𝓨::TransformedSigmaPoints{T}, xi_temp, F!, χ::SigmaPoints{T}) where T
    F!(𝓨.x0, χ.x0)
    num_x = length(χ.x0)
    xi_temp .= χ.x0
    @inbounds for i = length(χ.x0):-1:1
        xi_temp[i:num_x] .= @view(χ.x0[i:num_x]) .+ @view(χ.P_chol.data[i:num_x, i])
        F!(@view(𝓨.xi[:, i]), xi_temp)
        xi_temp[i:num_x] .= @view(χ.x0[i:num_x]) .- @view(χ.P_chol.data[i:num_x, i])
        F!(@view(𝓨.xi[:, i + length(χ.x0)]), xi_temp)
    end
    TransformedSigmaPoints(𝓨.x0, 𝓨.xi, χ.weight_params)
end

function mean(𝓨::TransformedSigmaPoints)
    weight_0, weight_i = calc_mean_weights(𝓨)
    𝓨.x0 .* weight_0 .+ vec(sum(𝓨.xi, dims = 2)) .* weight_i
end

function mean!(y::Vector, 𝓨::TransformedSigmaPoints)
    weight_0, weight_i = calc_mean_weights(𝓨)
    sum!(y, 𝓨.xi)
    y .= y .* weight_i .+ 𝓨.x0 .* weight_0
end

function substract_mean(𝓨::TransformedSigmaPoints, y)
    𝓨 .- y
end

function substract_mean!(unbiased_𝓨::TransformedSigmaPoints, 𝓨::TransformedSigmaPoints, y)
    unbiased_𝓨 .= 𝓨 .- y
    TransformedSigmaPoints(unbiased_𝓨.x0, unbiased_𝓨.xi, 𝓨.weight_params)
end

function cov(unbiased_𝓨::TransformedSigmaPoints, Q::AbstractMatrix)
    weight_0, weight_i = calc_cov_weights(unbiased_𝓨)
    cov(unbiased_𝓨::TransformedSigmaPoints, Augment(Q)) .+ Q
end

function cov!(P, unbiased_𝓨::TransformedSigmaPoints, Q::AbstractMatrix)
    weight_0, weight_i = calc_cov_weights(unbiased_𝓨)
    P .= @~ unbiased_𝓨.x0 * unbiased_𝓨.x0' .* weight_0 .+ unbiased_𝓨.xi * unbiased_𝓨.xi' .* weight_i .+ Q
end

function cov(χ::SigmaPoints, unbiased_𝓨::TransformedSigmaPoints)
    weight_0, weight_i = calc_cov_weights(χ)
    num_states = length(χ.x0)
    (χ.P_chol * (@view(unbiased_𝓨.xi[:, 1:num_states]))' .-
        χ.P_chol * (@view(unbiased_𝓨.xi[:, num_states + 1:2 * num_states]))') .* weight_i
end

function cov!(P, χ::SigmaPoints, unbiased_𝓨::TransformedSigmaPoints)
    weight_0, weight_i = calc_cov_weights(χ)
    num_states = length(χ.x0)
    P .= @~ χ.P_chol * (@view(unbiased_𝓨.xi[:, 1:num_states]))'
    P .-= @~ χ.P_chol * (@view(unbiased_𝓨.xi[:, num_states + 1:2 * num_states]))'
    P .*= weight_i
end

function mean_and_cov(𝓨::TransformedSigmaPoints, Q)
    y = mean(𝓨)
    Ryy = cov(𝓨 .- y, Q)
    y, Ryy
end

function calc_mean_weights(χ::AbstractSigmaPoints)
    calc_mean_weights(χ.weight_params, (size(χ, 2) - 1) >> 1)
end

function calc_cov_weights(χ::AbstractSigmaPoints)
    calc_cov_weights(χ.weight_params, (size(χ, 2) - 1) >> 1)
end

function calc_cholesky_weight(weight_params::AbstractWeightingParameters, P)
    calc_cholesky_weight(weight_params, size(P, 2))
end

Base.size(S::SigmaPoints) = (length(S.x0), 2 * length(S.x0) + 1)

Base.getindex(S::SigmaPoints, inds::Vararg{Int,2}) =
    @inbounds if inds[2] == 1
        S.x0[inds[1]]
    elseif 1 < inds[2] <= length(S.x0) + 1
        S.x0[inds[1]] + S.P_chol[inds[1], inds[2] - 1]
    else
        S.x0[inds[1]] - S.P_chol[inds[1], inds[2] - length(S.x0) - 1]
    end

# make TransformedSigmaPoints broadcastable
Base.size(S::TransformedSigmaPoints) = (length(S.x0), size(S.xi, 2) + 1)

Base.getindex(S::TransformedSigmaPoints{T}, inds::Vararg{Int,2}) where {T} =
    @inbounds if inds[2] == 1
        S.x0[inds[1]]
    else
        S.xi[inds[1], inds[2] - 1]
    end

Base.setindex!(S::TransformedSigmaPoints{T}, val, inds::Vararg{Int,2}) where {T} =
    @inbounds if inds[2] == 1
        S.x0[inds[1]] = val
    else
        S.xi[inds[1], inds[2] - 1] = val
    end

Base.BroadcastStyle(::Type{<:TransformedSigmaPoints}) = Broadcast.ArrayStyle{TransformedSigmaPoints}()

Base.similar(bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{TransformedSigmaPoints}}, ::Type{T}) where T =
    TransformedSigmaPoints(similar(bc.args[1].x0), similar(bc.args[1].xi), bc.args[1].weight_params)

function Base.copyto!(dest::TransformedSigmaPoints, bc::Broadcast.Broadcasted{Broadcast.ArrayStyle{TransformedSigmaPoints}})
    for idx in eachindex(bc)
        @inbounds dest[idx] = bc[idx]
    end
    dest
end


text:
function all_close(x::Array, y::Array, rtol::Float64=1.e-5, atol::Float64=1.e-8)
    """
    Returns True if two arrays are element-wise equal within a tolerance.

    References
    ==========
    This function, and docstring, are taken directly from numpy.allclose

    See that documentation for more info

    """
    xinf = isinf(x)
    yinf = isinf(y)
    if any(xinf) || any(yinf)
        if !all(xinf .== yinf)
            return false
        end

        if !all(x[xinf] .== y[yinf])
            return false
        end

        x = x[~xinf]
        y = y[~xinf]
    end

    return all(.<=(abs(x - y), atol + rtol * abs(y)))
end

reload("../stable_reg.jl")
using MAT

# load in matlab results
ml = matread("num_stab_approxMATLAB.mat")

# Grab data
X = ml["X"]
Y = ml["Y"]

# compare to matlab answers

# OLS
@assert all_close(OLS(X, Y), ml["b_ols_yes"])
@assert all_close(OLS(X, Y, false), ml["b_ols_no"])

# SVD
@assert all_close(LS_SVD(X, Y), ml["b_ols_yes"])
@assert all_close(LS_SVD(X, Y, false), ml["b_ols_no"])

# LAD-PP
@assert all_close(LAD_PP(X, Y), ml["b_ladpp_yes"])
@assert all_close(LAD_PP(X, Y, false), ml["b_ladpp_no"])

# LAD-DP
@assert all_close(LAD_DP(X, Y), ml["b_laddp_yes"])
@assert all_close(LAD_DP(X, Y, false), ml["b_laddp_no"])

# Regularized least square Tikhonov
@assert all_close(RLS_T(X, Y), ml["b_ols_yes"])

# Regularized SVD Tikhonov
@assert all_close(RLS_TSVD(X, Y), ml["b_ols_yes"])

max_issues_repo_name: EconForge/Smolyak

max_forks_repo_name: EconForge/Smolyak


text:
using Birdsong
using Base.Test

import Birdsong: is_adjacent
include("test_syllablize.jl")
include("test_dftrace.jl")
include("test_summaryStats.jl")


text:
using DataFrames
using DuckDB
using Test
using Dates
using UUIDs

test_files = [
    "test_appender.jl",
    "test_basic_queries.jl",
    "test_config.jl",
    "test_connection.jl",
    "test_df_scan.jl",
    "test_prepare.jl",
    "test_transaction.jl",
    "test_sqlite.jl",
    "test_replacement_scan.jl",
    "test_table_function.jl",
    "test_old_interface.jl",
    "test_all_types.jl"
]

if size(ARGS)[1] > 0
    filtered_test_files = []
    for test_file in test_files
        if test_file == ARGS[1]
            push!(filtered_test_files, test_file)
        end
    end
    test_files = filtered_test_files
end

for fname in test_files
    println(fname)
    include(fname)
end


text:
using InverseLaplace

function test0()
    F = s -> 1/sqrt(s^2+1)
    t = [0.1, 1,10]
    N = 64
    sig =1.0
    b = 1.0
    f = weeks(F,t,N,sig,b)
    println(f)
    t1 = 0.1
    f = weeks(F,t1,N,sig,b)
    println(f)

    f = weekse(F,t1,N,sig,b)

end

function test1()
    F = s -> 1/sqrt(s^2+1)
    t = [0.1, 1,10]
    N = 64
    sig =1.0
    b = 1.0
    f = weeks(F,t,N,sig,b)
    fex = besselj0(t)
    println("abserr first: ",  abs(f-fex))
    (f, est) = weekse(F,t,N,sig,b)
    println("est error: ", est)
end

function test1a()
    F = s -> 1/sqrt(s^2+1)
    t = [0.1, 1,10]
    N = 64
    sig =1.0
    b = 1.0
    f = weeks(F,t,N,sig,b)
end

function test2()
    res = weeks( s -> s/(1+s^2) , [1.0,2.0,3.0],  18, 1.0,1.0)
    abs(res - [0.5403022953181169,-0.41614641727651525, -0.9899933071110818])
end

function test3()
    res = weekse( s -> s/(1+s^2) , [1.0,2.0,3.0],  18, 1.0,1.0)
    res
end

function test4(N=32)
    F(s) = exp(-4*sqrt(s))
    t = 1.0
    sig = 0.7
    b = 1.75
    res = weeks(F,[t],N,sig,b)[1]
    println(abs(res - 0.020675159969109663))
    fex = 2*exp(-4/t)/sqrt(pi*t^3)
    abserr = abs(res-fex)
    println("First abserr ", abserr)

    println("Difference from expected: ", abs(8.174615017609438e-6 - abserr))

    (so,bo) = wpar2(F,t,N,0,30,30)

    res = weeks(F,t,N,so,bo)

    abserr = abs(res-fex)
    println("Second abserr ", abserr)
    println("Difference from expected: ", abs(1.6406273631308643e-9  - abserr))
end

function test4a(N=32)
    F(s) = exp(-4*sqrt(s))
    t = 1.0
    sig = 0.7
    b = 1.75
    local res
    @time (so,bo) = wpar2(F,t,N,0,30,30)
    @time res = weeks(F,t,N,so,bo)
    res
end

function timeweeks(F, t::AbstractVector, N, sig, b)
    t1 = @elapsed  a0 = real(wcoef(F,N,sig,b))
    a = a0[N+1:2*N]
    t2 = @elapsed L = laguer(a,2*b*t)
    t3 = @elapsed f = L .* exp((sig-b)*t)
    (t1,t2,t3)
end

function runweeks()
    F = s -> 1/sqrt(s^2+1)
    t = [0.1, 1,10]
    N = 80
    sig =1.0
    b = 1.0
    ts1 = 0.0
    ts2 = 0.0
    ts3 = 0.0
    for i in 1:3*10^2
        (t1,t2,t3) = timeweeks(F,t,N,sig,b)
        ts1 += t1
        ts2 += t2
        ts3 += t3
    end
#    (ts1,ts2,ts3)
    (1,ts2/ts1,ts3/ts1)
end

function test5()
    F = s -> 1/sqrt(s^2+1)
    t = 100.0
    N = 100
    (so,bo) = wpar2(F,t,N,0,30,30)
    println("so $so, bo $bo")
    for t1 in 1//10:1//4:5
        println(t1)
        f = weeks(F,t1,N,so,bo)
        fex = besselj0(t1)
        ae  = abs(f-fex)
        ae2 = abs(ilt(F,BigFloat(t1), 64)- fex)
        @printf("t %.2f\t   ae %.4e    ae2 %.4e : fex %e ,  f %e \n", t1, ae, ae2, fex, f)
    end
end


text:
@testset "minimal_sum_test" begin

    @test sum([1.0, 2.0, 3.0]) === 6.0

    @test isnan(sum([1.0, 2.0, NaN, 3.0]))

    @test isnan(sum([1.0, -Inf, 2.0, Inf, 3.0]))

end

@testset "minimal_sum_abs_test" begin

    @test sum(abs.([1.0, -2.0, 3.0])) === 6.0

    @test isnan(sum(abs.([1.0, -2.0, NaN, 3.0])))

    @test sum(abs.([1.0, -Inf, 2.0, Inf, 3.0])) === Inf

end

@testset "minimal_sum_sqr_test" begin

    @test sum([1.0, 2.0, 3.0].^2) === 14.0

    @test isnan(sum([1.0, 2.0, NaN, 3.0].^2))

    @test sum([1.0, -Inf, 2.0, Inf, 3.0].^2) === Inf

end

@testset "minimal_dot_test" begin

    @test dot([1.0, 2.0, 3.0], [1.0, 2.0, 3.0]) === 14.0

    @test_broken dot([0x10000000000001p0, 0x1p104], [0x0fffffffffffffp0, -1.0]) === -1.0

    @test isnan(dot([1.0, 2.0, NaN, 3.0], [1.0, 2.0, 3.0, 4.0]))

    @test isnan(dot([1.0, 2.0, 3.0, 4.0], [1.0, 2.0, NaN, 3.0]))

    @test isnan(dot([1.0, 2.0, 0.0, 4.0], [1.0, 2.0, Inf, 3.0]))

    @test isnan(dot([1.0, 2.0, -Inf, 4.0], [1.0, 2.0, 0.0, 3.0]))

end


text:
# TODO: Make into proper tests and integrate in testset

import TensorKit: BraidingTensor

V1 = GradedSpace{FermionSpin}(0=>2, 1/2=>2, 1=>1, 3/2=>1)

V2 = GradedSpace{FibonacciAnyon}(:I => 2, :τ => 2)

V3 = GradedSpace{IsingAnyon}(:I => 2, :psi => 1, :sigma => 1)

for V in (V1, V2, V3)
    @show V

    t = TensorMap(randn, V * V' * V' * V, V * V')

    ττ = copy(BraidingTensor(V, V'))
    @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[-1 -2; 1 2]*t[1 2 -3 -4; -5 -6]
    @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[-1 -2; 1 2]*t[1 2 -3 -4; -5 -6]
    @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[2 1; -2 -1]*t[1 2 -3 -4; -5 -6]
    @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[-2 2; -1 1]*t[1 2 -3 -4; -5 -6]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V', V'))
    @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[-2 -3; 1 2]*t[-1 1 2 -4; -5 -6]
    @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[-2 -3; 1 2]*t[-1 1 2 -4; -5 -6]
    @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[2 1; -3 -2]*t[-1 1 2 -4; -5 -6]
    @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[-3 2; -2 1]*t[-1 1 2 -4; -5 -6]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V', V))
    @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[1 -2; 2 -3]*t[-1 1 2 -4; -5 -6]
    @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[1 -2; 2 -3]*t[-1 1 2 -4; -5 -6]
    @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[-3 2; -2 1]*t[-1 1 2 -4; -5 -6]
    @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[-2 -3; 1 2]*t[-1 1 2 -4; -5 -6]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V, V'))
    @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[-3 2; -2 1]*t[-1 1 2 -4; -5 -6]
    @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[-3 2; -2 1]*t[-1 1 2 -4; -5 -6]
    @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[1 -2; 2 -3]*t[-1 1 2 -4; -5 -6]
    @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[2 1; -3 -2]*t[-1 1 2 -4; -5 -6]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V, V))
    @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[2 1; -3 -2]*t[-1 1 2 -4; -5 -6]
    @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[2 1; -3 -2]*t[-1 1 2 -4; -5 -6]
    @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[-2 -3; 1 2]*t[-1 1 2 -4; -5 -6]
    @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[1 -2; 2 -3]*t[-1 1 2 -4; -5 -6]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V', V))
    @planar2 t1[-1 -2 -3 -4; -5 -6] := τ[-3 -4; 1 2]*t[-1 -2 1 2; -5 -6]
    @planar2 t2[-1 -2 -3 -4; -5 -6] := ττ[-3 -4; 1 2]*t[-1 -2 1 2; -5 -6]
    @planar2 t3[-1 -2 -3 -4; -5 -6] := τ[2 1; -4 -3]*t[-1 -2 1 2; -5 -6]
    @planar2 t4[-1 -2 -3 -4; -5 -6] := τ'[-4 2; -3 1]*t[-1 -2 1 2; -5 -6]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V', V))
    @planar2 t1[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2]*τ[1 2; -5 -6]
    @planar2 t2[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2]*ττ[1 2; -5 -6]
    @planar2 t3[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2]*τ[-6 -5; 2 1]
    @planar2 t4[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2]*τ'[2 -6; 1 -5]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V, V'))
    @planar2 t1[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2]*τ'[1 2; -5 -6]
    @planar2 t2[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2]*ττ'[1 2; -5 -6]
    @planar2 t3[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2]*τ'[-6 -5; 2 1]
    @planar2 t4[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 -4; 1 2]*τ[2 -6; 1 -5]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V, V))
    @planar2 t1[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2]*τ[-4 -6; 1 2]
    @planar2 t2[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2]*ττ[-4 -6; 1 2]
    @planar2 t3[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2]*τ[2 1; -6 -4]
    @planar2 t4[-1 -2 -3 -4; -5 -6] := t[-1 -2 -3 1; -5 2]*τ'[-6 2; -4 1]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V', V))
    @planar2 t1[();(-1,-2)] := τ[2 1; 3 4]*t[1 2 3 4; -1 -2]
    @planar2 t2[();(-1,-2)] := ττ[2 1; 3 4]*t[1 2 3 4; -1 -2]
    @planar2 t3[();(-1,-2)] := τ[4 3; 1 2]*t[1 2 3 4; -1 -2]
    @planar2 t4[();(-1,-2)] := τ'[1 4; 2 3]*t[1 2 3 4; -1 -2]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V, V))
    @planar2 t1[-1; -2] := τ[2 1; 3 4]*t[-1 1 2 3; -2 4]
    @planar2 t2[-1; -2] := ττ[2 1; 3 4]*t[-1 1 2 3; -2 4]
    @planar2 t3[-1; -2] := τ[4 3; 1 2]*t[-1 1 2 3; -2 4]
    @planar2 t4[-1; -2] := τ'[1 4; 2 3]*t[-1 1 2 3; -2 4]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V, V'))
    @planar2 t1[-1 -2] := τ[2 1; 3 4]*t[-1 -2 1 2; 4 3]
    @planar2 t2[-1 -2] := ττ[2 1; 3 4]*t[-1 -2 1 2; 4 3]
    @planar2 t3[-1 -2] := τ[4 3; 1 2]*t[-1 -2 1 2; 4 3]
    @planar2 t4[-1 -2] := τ'[1 4; 2 3]*t[-1 -2 1 2; 4 3]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V, V'))
    @planar2 t1[-1 -2; -3 -4] := τ[-1 3; 1 2]*t[1 2 3 -2; -3 -4]
    @planar2 t2[-1 -2; -3 -4] := ττ[-1 3; 1 2]*t[1 2 3 -2; -3 -4]
    @planar2 t3[-1 -2; -3 -4] := τ[2 1; 3 -1]*t[1 2 3 -2; -3 -4]
    @planar2 t4[-1 -2; -3 -4] := τ'[3 2; -1 1]*t[1 2 3 -2; -3 -4]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V', V'))
    @planar2 t1[-1 -2; -3 -4] := τ'[-2 3; 1 2]*t[-1 1 2 3; -3 -4]
    @planar2 t2[-1 -2; -3 -4] := ττ'[-2 3; 1 2]*t[-1 1 2 3; -3 -4]
    @planar2 t3[-1 -2; -3 -4] := τ'[2 1; 3 -2]*t[-1 1 2 3; -3 -4]
    @planar2 t4[-1 -2; -3 -4] := τ[3 2; -2 1]*t[-1 1 2 3; -3 -4]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V', V))
    @planar2 t1[-1 -2 -3; -4] := τ[-3 3; 1 2]*t[-1 -2 1 2; -4 3]
    @planar2 t2[-1 -2 -3; -4] := ττ[-3 3; 1 2]*t[-1 -2 1 2; -4 3]
    @planar2 t3[-1 -2 -3; -4] := τ[2 1; 3 -3]*t[-1 -2 1 2; -4 3]
    @planar2 t4[-1 -2 -3; -4] := τ'[3 2; -3 1]*t[-1 -2 1 2; -4 3]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V', V))
    @planar2 t1[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2]*τ[1 2; -4 3]
    @planar2 t2[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2]*ττ[1 2; -4 3]
    @planar2 t3[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2]*τ[3 -4; 2 1]
    @planar2 t4[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2]*τ'[2 3; 1 -4]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)

    ττ = copy(BraidingTensor(V, V'))
    @planar2 t1[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2]*τ'[1 2; -4 3]
    @planar2 t2[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2]*ττ'[1 2; -4 3]
    @planar2 t3[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2]*τ'[3 -4; 2 1]
    @planar2 t4[-1 -2 -3; -4] := t[-1 -2 -3 3; 1 2]*τ[2 3; 1 -4]
    @show norm(t1 - t2), norm(t1 - t3), norm(t1 - t4)
end