text:
"""
  Abstract type that defines the conversion between conservative variables
  and some other set of variables, typically the entropy variables.

  Required functions:

  ```
    convertToEntropy
    convertToConservative
    getA0
    getA0inv
  ```

  See those functions docstrings for details
"""
abstract type AbstractVariables end

"""
  Convert to entropy variables from conservative variables

  **Inputs**

   * obj: the [`AbstractVariables`](@ref) object specifying what variables
          to convert to
   * params: AbstractParamType object
   * qc: the conservative variables at the node, length `numDofPerNode`

  **Inputs/Outputs**

   * qe: vector to be overwritten with the entropy variables, same length
         as `qc`

  Aliasing: `qc` and `qe` may alias.  Implementations must ensure correct
            behavior in this case
"""
function convertToEntropy(obj::AbstractVariables, params::AbstractParamType,
                          qc::AbstractVector, qe::AbstractVector)

  error("abstract method for convertToEntropy() reached")
end

"""
  Convert to conservative variables from entropy variables

  **Inputs**

   * obj: the [`AbstractVariables`](@ref) object specifying what variables
          to convert to
   * params: AbstractParamType object
   * qe: the entropy variables at the node, length `numDofPerNode`

  **Inputs/Outputs**

   * qc: vector to be overwritten with the conservative variables, same length
         as `qe`

  Aliasing: `qc` and `qe` may alias.  Implementations must ensure correct
            behavior in this case
"""
function convertToConservative(obj::AbstractVariables, params::AbstractParamType,
                          qe::AbstractVector, qc::AbstractVector)

  error("abstract method for convertToConservative() reached")
end

"""
  Computes the Jacobian dq/dw, where `q` are the conservative variables and
  `w` are the entropy variables`.

  **Inputs**

   * obj: an [`AbstractVariables`](@ref) object
   * params: AbstractParamType object
   * qc: vector of conservative variables specifying the state at which to
         compute `A0`, length `numDofPerNode`

  **Inputs/Outputs**

   * A0: matrix to be overwritten with `dq/dw, `numDofPerNode` x `numDofPerNode`
"""
function getA0(obj::AbstractVariables, params::AbstractParamType,
               qc::AbstractVector, A0::AbstractMatrix)

  error("abstract method for getA0() reached")
end

"""
  Computes the Jacobian dw/dq = inv(dq/dw), where `q` are the conservative
  variables and `w` are the entropy variables`.

  **Inputs**

   * obj: an [`AbstractVariables`](@ref) object
   * params: AbstractParamType object
   * qc: vector of conservative variables specifying the state at which to
         compute `A0`, length `numDofPerNode`

  **Inputs/Outputs**

   * A0inv: matrix to be overwritten with `dw/dq, `numDofPerNode` x
           `numDofPerNode`
"""
function getA0inv(obj::AbstractVariables, params::AbstractParamType,
                  qc::AbstractVector, A0inv::AbstractMatrix)

  error("abstract method for getA0inv() reached")
end

"""
  Computes the derivative of [`getA0`](@ref) with respect to the conservative
  variables.

  **Inputs**

   * obj: an [`AbstractVariables`](@ref) object
   * params: `AbstractParamType` object
   * q: vector of conservative variables, length `numDofPerNode`
   * q_dot: derivative part of `q`, `numDofPerNode` x `nd`, where `nd` defines
            the number of derivative parts

  **Inputs/Outputs**

   * A0: same as [`getA0`](@ref)
   * A0_dot: `numDofPerNode` x `numDofPerNode` x `nd` array to be overwritten
             with the output.  The final dimension `nd` must be the same as
             the final dimension of `q_dot`.
"""
function getA0_diff(obj::AbstractVariables, params::AbstractParamType,
                    q::AbstractVector,
                    q_dot::AbstractMatrix,
                    A0::AbstractMatrix,
                    A0_dot::Abstract3DArray)

  error("abstract method for getA0_diff() reached")
end

"""
  Computes the reverse mode of getA0 wrt `q`

  **Inputs**

   * obj: an [`AbstractVariables`](@ref) object
   * params: `AbstractParamType` object
   * q: vector of conservative variables, length `numDofPerNode`
   * A0_bar: adjoint part of `A0`, same size as `A0`.

  **Inputs/Outputs**

   * q_bar: adjoint part of `q`, to be updated.  Same size as `q`
   * A0: same as [`getA0`](@ref)
"""
function getA0_revq(obj::AbstractVariables, params::AbstractParamType,
                     q::AbstractArray, q_bar::AbstractArray,
                     A0::AbstractMatrix, A0_bar::AbstractMatrix)

  error("abstract method for getA0_revq() reached")
end

#TODO: add: entropy function, entropy flux

#------------------------------------------------------------------------------
# Variables used by Ismail and Roe

"""
  Variables used by Ismail and Roe.  See "Affordable, entropy-consistent Euler
  Flux Functions II: Entropy Production at Shocks"
"""
struct IRVariables <: AbstractVariables
end

function convertToEntropy(obj::IRVariables, params::AbstractParamType,
                          qc::AbstractVector, qe::AbstractVector)

  convertToIR_(params, qc, qe)
end

function convertToConservative(obj::IRVariables, params::AbstractParamType,
                          qe::AbstractVector, qc::AbstractVector)

  convertToConservativeFromIR(params, qe, qc)
end

function getA0(obj::IRVariables, params::AbstractParamType, qc::AbstractVector,
               A0::AbstractMatrix)

  getIRA0(params, qc, A0)
end

function getA0inv(obj::IRVariables, params::AbstractParamType, qc::AbstractVector,
                  A0inv::AbstractMatrix)

  getIRA0inv(params, qc, A0inv)
end

function getA0_diff(obj::IRVariables, params::AbstractParamType,
                    q::AbstractVector,
                    q_dot::AbstractMatrix,
                    A0::AbstractMatrix,
                    A0_dot::Abstract3DArray)

  getIRA0_diff(params, q, q_dot, A0, A0_dot)
end

function getA0_revq(obj::IRVariables, params::AbstractParamType,
                     q::AbstractArray, q_bar::AbstractArray,
                     A0::AbstractMatrix, A0_bar::AbstractMatrix)

  getIRA0_revq(params, q, q_bar, A0, A0_bar)
end

#------------------------------------------------------------------------------
# Hughes entropy variables

"""
  Variables used by Hughes.  See "A New Finite Element Formulation For
  Computational Fluid Dynamics: Part I"
"""
struct HughesVariables <: AbstractVariables
end

function convertToEntropy(obj::HughesVariables, params::AbstractParamType,
                          qc::AbstractVector, qe::AbstractVector)

  convertToEntropy_(params, qc, qe)
end

function convertToConservative(obj::HughesVariables, params::AbstractParamType,
                          qe::AbstractVector, qc::AbstractVector)

  #TODO: cleanup the naming scheme
  convertToConservative(params, qe, qc)
end

function getA0(obj::HughesVariables, params::AbstractParamType, qc::AbstractVector,
               A0::AbstractMatrix)

  #TODO: cleanup the naming scheme
  calcA0(params, qc, A0)
end

function getA0inv(obj::HughesVariables, params::AbstractParamType, qc::AbstractVector,
                  A0inv::AbstractMatrix)

  calcA0inv(params, qc, A0inv)
end

#------------------------------------------------------------------------------
# Conservative variables
# Useful for testing to make conservative = entropy

"""
  Trick the solver into thinking the conservative and entropy variables are
  the same.
"""
struct ConservativeVariables <: AbstractVariables
end

function convertToEntropy(obj::ConservativeVariables, params::AbstractParamType,
                          qc::AbstractVector, qe::AbstractVector)

  for i=1:length(qc)
    qe[i] = qc[i]
  end
end

function convertToConservative(obj::ConservativeVariables, params::AbstractParamType,
                          qe::AbstractVector, qc::AbstractVector)

  for i=1:length(qc)
    qc[i] = qe[i]
  end
end

function getA0(obj::ConservativeVariables, params::AbstractParamType, qc::AbstractVector,
               A0::AbstractMatrix)

  fill!(A0, 0)
  for i=1:size(A0, 1)
    A0[i, i] = 1
  end
end

function getA0inv(obj::ConservativeVariables, params::AbstractParamType, qc::AbstractVector,
                  A0inv::AbstractMatrix)

  fill!(A0inv, 0)
  for i=1:size(A0inv, 1)
    A0inv[i, i] = 1
  end
end

function getA0_diff(obj::ConservativeVariables, params::AbstractParamType,
                    q::AbstractVector,
                    q_dot::AbstractMatrix,
                    A0::AbstractMatrix,
                    A0_dot::Abstract3DArray)

  getA0(obj, params, q, A0)
  fill!(A0_dot, 0)
end


text:
using Test#, ReferenceTests, BSON

plotyes = false
include("../scripts/car_travels.jl")

## Unit tests
@testset "update_position" begin
    @test update_position(0.0, 10, 1, 1, 200)[1] ≈ 10.0
    @test update_position(0.0, 10, 1, 1, 200)[2] == 1

    @test update_position(0.0, 10, -1, 1, 200)[1] ≈ -10.0
    @test update_position(0.0, 10, -1, 1, 200)[2] == 1

    @test update_position(0.0, 10, -1, 1, 200)[1] ≈ -10.0
    @test update_position(0.0, 10, -1, 1, 200)[2] == 1
end

## Reference Tests with ReferenceTests.jl
# We put both arrays X and T into a BSON.jl and then compare them

"Compare all dict entries"
comp(d1, d2) = keys(d1) == keys(d2) &&
    all([ v1≈v2 for (v1,v2) in zip(values(d1), values(d2))])

# run the model
T, X = car_travel_1D()

# Test just at some random indices. As for larger models,
# storing the full output array would create really large files!
inds = [18, 27, 45, 68, 71, 71, 102, 110, 123, 144]

d = Dict(:X=> X[inds], :T=>T[inds])
# @testset "Ref-tests" begin
#     @test_reference "reftest-files/X.bson" d by=comp
# end

ext: jl


text:
module Libtask

export CTask, consume, produce, TArray, get, tzeros, tfill, TRef

if !isfile(joinpath((@__DIR__) |> dirname, "deps/desp.jl"))
    import Pkg
    Pkg.build("Libtask")
end

include("../deps/deps.jl"); check_deps();
include("taskcopy.jl")
include("tarray.jl")

end


text:
module CameraTargetCornerTriggerModule

using ..Ahorn, Maple

@mapdef Trigger "HonlyHelper/CameraTargetCornerTrigger" CameraTargetCornerTrigger(x::Integer, y::Integer, width::Integer=Maple.defaultTriggerWidth, height::Integer=Maple.defaultTriggerHeight, lerpStrength::Number=0.0, positionMode::String="BottomLeft", xOnly::Bool=false, yOnly::Bool=false, deleteFlag::String="", nodes::Array{Tuple{Integer, Integer}, 1}=Tuple{Integer, Integer}[])

const placements = Ahorn.PlacementDict(
    "Camera Target Corner (HonlyHelper)" => Ahorn.EntityPlacement(
        CameraTargetCornerTrigger,
        "rectangle",
        Dict{String, Any}(),
        function(trigger)
            trigger.data["nodes"] = [(Int(trigger.data["x"]) + Int(trigger.data["width"]) + 8, Int(trigger.data["y"]))]
        end
    )
)

const positionModeOptions = String["BottomLeft", "BottomRight", "TopLeft", "TopRight"]

Ahorn.editingOptions(trigger::CameraTargetCornerTrigger) = Dict{String, Any}(
        "positionMode" => positionModeOptions
    )

function Ahorn.nodeLimits(trigger::CameraTargetCornerTrigger)
    return 1, 1
end

end

max_stars_repo_name: HollyMagala/HonlyHelper

max_issues_repo_name: HollyMagala/HonlyHelper

max_forks_repo_name: HollyMagala/HonlyHelper


text:
using Pkg
Pkg.activate(@__DIR__)
using RCall
using SigmaRidgeRegression
using StatsBase
using Random
using MLJ
using DrWatson
using JLD2
using Distributions

opt_no = parse(Int64, ARGS[1])
@show opt_no

# helper functions
function _merge(grp::GroupedFeatures; groups_goal=2)
    ngroups = length(grp.ps)
    new_ps = Vector{Int64}(undef, groups_goal)
    mod(ngroups, groups_goal) == 0 ||
        throw("Only implemented when #groups_goal divides ngroups")
    step_length = div(ngroups, groups_goal)
    cnt = 1
    for i in 1:groups_goal
        cnt_upper = cnt + (step_length - 1)
        new_ps[i] = sum(grp.ps[cnt:cnt_upper])
        cnt = cnt_upper
    end
    return GroupedFeatures(new_ps)
end

function single_simulation(sim; Ks=Ks, save=true)
    res = []
    sim_name = randstring(16)

    groups = sim.groups

    _simulated_model = simulate(sim)
    X = MLJ.table(_simulated_model.X)
    Y = _simulated_model.Y
    resampling_idx = _simulated_model.resampling_idx

    bayes_λs =
        groups.ps .* var(sim.response_noise) ./ abs2.(sim.response_model.αs) ./ sim.ntrain
    bayes_ridge = MultiGroupRidgeRegressor(;
        groups=groups, λs=bayes_λs, center=false, scale=false
    )

    _mach = machine(bayes_ridge, X, Y)
    _eval = evaluate!(_mach; resampling=resampling_idx, measure=l2)
    mse_bayes = _eval.measurement[1]

    for K in Ks
        newgroups = _merge(groups; groups_goal=K)

        loo_single_ridge = LooRidgeRegressor(;
            ridge=SingleGroupRidgeRegressor(; groups=newgroups, center=false, scale=false)
        )

        loo_sigmaridge = LooSigmaRidgeRegressor(;
            groups=newgroups, center=false, scale=false
        )

        loo_multi_ridge = LooRidgeRegressor(;
            ridge=MultiGroupRidgeRegressor(; groups=newgroups, center=false, scale=false),
            rng=MersenneTwister(1),
        )

        holdout_glasso = TunedSeagull(; groups=newgroups, center=false, scale=false)

        models = [loo_sigmaridge, loo_single_ridge, loo_multi_ridge, holdout_glasso]

        tmp_mses = fill(Inf, length(models))

        for (model_idx, model) in enumerate(models)
            _mach = machine(model, X, Y)
            _eval = evaluate!(_mach; resampling=resampling_idx, measure=l2)
            tmp_mses[model_idx] = _eval.measurement[1]
        end

        push!(
            res,
            (
                mse_sigma=tmp_mses[1],
                mse_single=tmp_mses[2],
                mse_multi=tmp_mses[3],
                mse_glasso=tmp_mses[4],
                mse_bayes=mse_bayes,
                K=K,
                sim=sim,
                p=sim.groups.p,
                cov=sim.Σ,
                response=sim.response_model,
                sim_name=sim_name,
            ),
        )
    end
    if save
        @save "simulation_results/$(sim_name).jld2" res
    end
    return res
end

# Code that starts simulations

Ks = 2 .^ (1:5)

p = 32 * 25

ns = Int.([p / 2; p; 2p])

groups = GroupedFeatures(fill(25, 32))

ar1 = SigmaRidgeRegression.AR1Design(p, 0.8)
id = IdentityCovarianceDesign(p)

informative_response_model = RandomLinearResponseModel(; αs=(0:31) ./ 3.1, grp=groups)

all_opts = dict_list(Dict(:n => ns, :cov => [ar1, id]))
opt = all_opts[opt_no]

n = opt[:n]
@show n
Σ = opt[:cov]
@show Σ
nreps = 400
sim = GroupRidgeSimulationSettings(;
    groups=groups,
    Σ=Σ,
    response_noise=Normal(0, 5),
    response_model=informative_response_model,
    ntrain=n,
)

for i in Base.OneTo(nreps)
    @show i
    single_simulation(sim)
end


text:
module Qaintessent

using LinearAlgebra
using SparseArrays
using Memoize

include("util.jl")
export
    rdm

include("gates.jl")
export
    AbstractGate,
    X,
    Y,
    Z,
    XGate,
    YGate,
    ZGate,
    HadamardGate,
    SGate,
    TGate,
    SdagGate,
    TdagGate,
    RxGate,
    RyGate,
    RzGate,
    RotationGate,
    PhaseShiftGate,
    SwapGate,
    EntanglementXXGate,
    EntanglementYYGate,
    EntanglementZZGate,
    ControlledGate,
    controlled_not,
    MatrixGate,
    matrix,
    sparse_matrix,
    num_wires

include("circuitgate.jl")
export
    AbstractCircuitGate,
    CircuitGate,
    circuit_gate,
    req_wires

include("moment.jl")
export
    Moment

include("measurementoperator.jl")
export
    MeasurementOperator,
    mop

include("circuit.jl")
export
    Circuit,
    add_measurement!,
    distribution

include("density_matrix.jl")
export
    DensityMatrix,
    pauli_group_matrix,
    density_from_statevector,
    density_from_matrix

include("commute.jl")
    export
        iscommuting

include("apply.jl")
export
    apply

include("apply_density.jl")
export
    apply

include("gradients.jl")
include("gradients_density.jl")

include("view.jl")

include("openqasm/register.jl")
export
    qreg,
    creg,
    reg_check,
    set_creg!,
    add!,
    add_control!

include("openqasm/grammar.jl")
include("openqasm/gate_transformers.jl")

include("openqasm/transform_qasm.jl")
export
    qasm2cgc

include("openqasm/transform_cgc.jl")
export
    cgc2qasm

include("unitary2circuit.jl")
export
    unitary2circuit

include("measure.jl")
export
    measure

module MaxKColSubgraphQAOA
    include("qaoa/mixer_gates.jl")
    export
        RNearbyValuesMixerGate,
        ParityRingMixerGate,
        PartitionMixerGate,
        r_nearby_values_hamiltonian_onehot

    include("qaoa/phase_separator_gates.jl")
    export
        MaxKColSubgraphPhaseSeparationGate,
        Graph,
        max_k_col_subgraph_phase_separation_hamiltonian

    include("qaoa/qaoa_gradients.jl")
end

end


text:
@testset "Conversion" begin

    P,Q = Chebyshev, ChebyshevU
    x = variable()

    for _ in 1:10
        p = P(rand(1:5,5)); q=P(rand(1:5,5))
        @test convert(P, convert(Q, p))(x) ≈ p(x)
        @test convert(Q, convert(P,q))(x) ≈ q(x)
    end

end

@testset "Evaluation" begin

    p0 = ChebyshevU([1])  # 1
    p1 = ChebyshevU([0,1]) # 2x
    p2 = ChebyshevU([0,0,1]) # 4x^2 - 1
    p3 = ChebyshevU([0,0,0,1]) # 8x^3 - 4x
    p4 = ChebyshevU([0,0,0,0,1]) # 16x^4 - 12x^2 + 1

    x = 1/2
    @test p0(x) ≈ 1
    @test p1(x) ≈ 2x
    @test p2(x) ≈ 4x^2 - 1
    @test p3(x) ≈ 8x^3 - 4x
    @test p4(x) ≈ 16x^4 - 12x^2 + 1
    @test ChebyshevU([1,2,3,4])(x) ≈ 1p0(x) + 2p1(x) + 3p2(x) + 4p3(x)

end

@testset "roots" begin

    for n in 4:2:8
        p = ChebyshevU(vcat(zeros(n),1))
        @test cos.((n:-1:1)/(n+1) * pi) ≈ roots(p)
    end

end


text:
#Create typealiases like Mat4f0, Point2, Point2f0
for i=1:4
    for T=[:Point, :Vec]
        name 	= Symbol("$T$i")
        namef0 	= Symbol("$T$(i)f0")
        @eval begin
            typealias $name $T{$i}
            typealias $namef0 $T{$i, Float32}
            export $name
            export $namef0
        end
    end
    name   = Symbol("Mat$i")
    namef0 = Symbol("Mat$(i)f0")
    @eval begin
        typealias $name $Mat{$i,$i}
        typealias $namef0 $Mat{$i,$i, Float32}
        export $name
        export $namef0
    end
end

#Type aliases

"""
An alias for a one-simplex, corresponding to LineSegment{T} -> Simplex{2,T}.
"""
typealias LineSegment{T} Simplex{2,T}

typealias Triangle{T} Face{3, T, 0}
typealias GLFace{Dim} Face{Dim, Cuint, -1} #offset is relative to julia, so -1 is 0-indexed
typealias GLTriangle  Face{3, Cuint, -1}
typealias GLQuad      Face{4, Cuint, -1}

typealias Cube{T}   HyperCube{3, T}

"""
An alias for a HyperSphere of dimension 2. i.e. Circle{T} -> HyperSphere{2, T}
"""
typealias Circle{T} HyperSphere{2, T}

"""
An alias for a HyperSphere of dimension 3. i.e. Sphere{T} -> HyperSphere{3, T}
"""
typealias Sphere{T} HyperSphere{3, T}

typealias AbsoluteRectangle{T} HyperRectangle{2, T}

"""
AABB, or Axis Aligned Bounding Box, is an alias for a 3D HyperRectangle.
"""
typealias AABB{T} HyperRectangle{3, T}
@compat (::Type{AABB})(m...) = HyperRectangle(m...)

typealias HMesh HomogenousMesh

typealias UV{T} TextureCoordinate{2, T}
typealias UVW{T} TextureCoordinate{3, T}

"""
A `SimpleMesh` is an alias for a `HomogenousMesh` parameterized only by
vertex and face types.
"""
typealias SimpleMesh{VT, FT} HMesh{VT, FT, Void, Void, Void, Void, Void}
typealias PlainMesh{VT, FT} HMesh{Point{3, VT}, FT, Void, Void, Void, Void, Void}
typealias GLPlainMesh PlainMesh{Float32, GLTriangle}

typealias Mesh2D{VT, FT} HMesh{Point{2, VT}, FT, Void, Void, Void, Void, Void}
typealias GLMesh2D Mesh2D{Float32, GLTriangle}

typealias UVMesh{VT, FT, UVT} HMesh{Point{3, VT}, FT, Void, UV{UVT}, Void, Void, Void}
typealias GLUVMesh UVMesh{Float32, GLTriangle, Float32}

typealias UVWMesh{VT, FT, UVT} HMesh{Point{3, VT}, FT, Void, UVW{UVT}, Void, Void, Void}
typealias GLUVWMesh UVWMesh{Float32, GLTriangle, Float32}

typealias NormalMesh{VT, FT, NT} HMesh{Point{3, VT}, FT, Normal{3, NT}, Void, Void, Void, Void}
typealias GLNormalMesh NormalMesh{Float32, GLTriangle, Float32}

typealias UVMesh2D{VT, FT, UVT} HMesh{Point{2, VT}, FT, Void, UV{UVT}, Void, Void, Void}
typealias GLUVMesh2D UVMesh2D{Float32, GLTriangle, Float32}

typealias NormalColorMesh{VT, FT, NT, CT} HMesh{Point{3, VT}, FT, Normal{3, NT}, Void, CT, Void, Void}
typealias GLNormalColorMesh NormalColorMesh{Float32, GLTriangle, Float32, RGBA{Float32}}

typealias NormalVertexcolorMesh{VT, FT, NT, CT} HMesh{Point{3, VT}, FT, Normal{3, NT}, Void, Vector{CT}, Void, Void}
typealias GLNormalVertexcolorMesh NormalVertexcolorMesh{Float32, GLTriangle, Float32, RGBA{Float32}}

typealias NormalAttributeMesh{VT, FT, NT, AT, A_ID_T} HMesh{Point{3, VT}, FT, Normal{3, NT}, Void, Void, AT, A_ID_T}
typealias GLNormalAttributeMesh NormalAttributeMesh{Float32, GLTriangle, Float32, Vector{RGBA{U8}}, Float32}

typealias NormalUVWMesh{VT, FT, NT, UVT} HMesh{Point{3, VT}, FT, Normal{3, NT}, UVW{UVT}, Void, Void, Void}
typealias GLNormalUVWMesh NormalUVWMesh{Float32, GLTriangle, Float32, Float32}

typealias NormalUVMesh{VT, FT, NT, UVT} HMesh{Point{3, VT}, FT, Normal{3, NT}, UV{UVT}, Void, Void, Void}
typealias GLNormalUVMesh NormalUVMesh{Float32, GLTriangle, Float32, Float32}


text:
Wigner = PyCall.PyNULL()
Optics = PyCall.PyNULL()

function init_physics()
    copy!(Wigner, PyCall.pyimport_conda("sympy.physics.wigner", "sympy"))
    copy!(Optics, PyCall.pyimport_conda("sympy.physics.optics", "sympy"))
end

module Physics

using PyCall
using SymPy

physics = [(SymPy.Wigner, [:clebsch_gordan,
                           :dot_rot_grad_Ynm,
                           :gaunt,
                           :racah,
                           :wigner_3j,
                           :wigner_6j,
                           :wigner_9j]),
           (SymPy.Optics, [:RayTransferMatrix,
                           :FreeSpace,
                           :CurvedRefraction,
                           :FlatMirror])
]

for (m, meths) in physics
    for meth in meths
        meth_name = string(meth)
        @eval begin
            @doc """
            `$($meth_name)`: a SymPy function.
                The SymPy documentation can be found through: http://docs.sympy.org/latest/search.html?q=$($meth_name)
                """ ->
            ($meth)(args...;kwargs...) = pycall(getindex($m, $meth_name), PyAny, map(Sym, args)..., map(Sym, kwargs)...)
        end
        eval(Expr(:export, meth))
    end
end

end


text:
# -*- encoding: utf-8 -*-
#
# The MIT License (MIT)
#
# Copyright © 2021 Matteo Foglieni and Riccardo Gervasoni
#

@testset "test_OrthoNormalBasis" begin
    pcg = PCG()

    for i in 1:100000
        normal = Normal(random(pcg), random(pcg), random(pcg))
        e1, e2, e3 = create_onb_from_z(normal)

        @test e3 ≈ normal # z axis must be aligned with the normal

        # testing orthogonality
        @test are_close(e1 ⋅ e2, 0.)
        @test are_close(e1 ⋅ e3, 0.)
        @test are_close(e3 ⋅ e2, 0.)

        # testing normality
        @test are_close(squared_norm(e1), 1.)
        @test are_close(squared_norm(e2), 1.)
        @test are_close(squared_norm(e3), 1.)
    end
end


text:
#  Without using ApproxFun
# function integrate_B_full(n::Int, X, Y0; Y1 =1000000, θin = 0.0, num_coefs = 1000000)
#     Ys = LinRange(Y0,Y1,num_coefs);
#     σs = integration_scheme(Ys; scheme = :trapezoidal)
#     K(Y) = cos(Y*sin(θin) + n*atan(Y,X))*hankelh1(n,sqrt(X^2+Y^2))
#     return 2.0*(-1.0)^n*sum(K.(Ys).*σs)
# end

function integrate_B_full(n::Int, X, Y0; Y1 = 1000000, θin = 0.0 + 0.0im, num_coefs = 80000)
    K(Y) = cos(Y*sin(θin) + n*atan(Y,X))*hankelh1(n,sqrt(X^2+Y^2))

    return 2.0*(-1.0)^n * hquadrature(K,Y0,Y1; maxevals = num_coefs)[1]
    # approximate function with Chebyshev polynomial (to high precision) then integrate from Y0 to Y1
    # return 2.0*(-1.0)^n * sum(Fun(K,Y0..Y1, num_coefs))
end

# Y0 = sqrt(k^a12^2 - X^2)
function integrate_B(n::Int, X::T, Y0::T;
        θin::Union{T,Complex{T}} = 0.0, num_coefs::Int = 80000,
        Y1::T = max(T(2000)*X, T(4000)) # note Y1 is non-dimensional!
    ) where T<:AbstractFloat

    # assymptotically approximate the integral from Y1 to Inf (tested in integrate_hankels.nb)
    Binf = (1.0+1.0im)*exp(im*Y1*(1.0 - sin(θin)))*
        (1.0 + (-1.0)^n*exp(2.0im*Y1*sin(θin))*(1.0 - sin(θin)) + sin(θin))/(sqrt(pi*Y1)*cos(θin)^2)

    return Binf + integrate_B_full(n, X, Y0; Y1=Y1, θin=θin, num_coefs = num_coefs)
end

# for only whole-correction, this doesn't involve an integral
function integrate_S(n::Int, X::T; θin::Union{T,Complex{T}} = 0.0) where T <: AbstractFloat
    S = 2.0*(im^T(n))*exp(-im*n*θin)*exp(im*X*cos(θin))/cos(θin)
    if X<0 S = conj(S) end
    S
end

function BS_matrices(X::AbstractVector{T}, a12k::T; θin::Union{T,Complex{T}} = 0.0 + 0.0im,
        M::Int = 2, num_coefs::Int = 80000) where T<:AbstractFloat

    dX = X[2] - X[1]
    J = length(collect(X)) -1
    q = min(Int(floor(a12k/dX)),J)
    X = OffsetArray((-J:J)*dX, -J:J)

    B = OffsetArray{Complex{Float64}}(undef, -q:q, -2M:2M);
    for j = -q:q, m = -2M:2M
        if a12k^2 - X[j]^2 < - dX^2 error("evaluating B in the wrong domain") end
        B[j,m] = integrate_B(m, X[j], sqrt(abs(a12k^2 -X[j]^2)); θin = θin, num_coefs=num_coefs)
    end

    S = OffsetArray{Complex{Float64}}(undef, -J:J, -2M:2M);
    for j = -J:J, m = -2M:2M
        S[j,m] = integrate_S(m, X[j]; θin = θin)
    end

    return (B,S)
end

function intergrand_kernel(X::AbstractVector{T}, a12k::T; M::Int = 2,
        scheme::Symbol = :trapezoidal, kws...) where T<:AbstractFloat

    dX = X[2] - X[1]
    J = length(collect(X)) -1
    if abs(J*dX - X[end])/X[end] > 1e-10
        error("Expected X to be uniformly spaced from 0. Instead got X= $X.")
    end
    if ( abs(Int(round(a12k/dX)) - a12k/dX) > 1e-10 )
        @warn("There are no mesh points exactly on-top of the intergrands kinks. This could lead to poor accuracy.")
    end

    q = min(Int(floor(a12k/dX)),J)
    if q == 0 @warn("Mesh element larger than ka12. This is only accurate for low frequency.") end

    B, S = BS_matrices(X, a12k; M = M, kws...)

    # an numerical integration scheme for domain over X
    σ = OffsetArray(integration_scheme(X; scheme=scheme), 0:J)
    # an array of numerical integration schemes for the varying domain |j-l| <= q
    σs = OffsetArray([
        OffsetArray(
            integration_scheme(X[max(l-q,0)+1:min(l+q,J)+1]; scheme=scheme),
        max(l-q,0):min(l+q,J))
    for l=0:J], 0:J)

    intergrand_quad = [
        begin
            P = σ[j] * S[j-l,n-m]
            Q = (abs(j-l) <= q) ? σs[l][j] * (B[j-l,n-m] - S[j-l,n-m]) : zero(Complex{T})
            P + Q
        end
    for l=0:J, m=-M:M, j=0:J, n=-M:M]

    return intergrand_quad
end

max_stars_repo_name: arturgower/EffectiveWaves


text:
exprs = map(s -> replace(s, '*' => '⨣'), readlines("day18/input.txt"))
exprs = map(s -> replace(s, '+' => '⨱'), exprs)

# picking an operator that looks like a plus gives it the precedence of addition
⨣(a, b) = a * b
⨱(a, b) = a + b

println(reduce((v, l) -> v + eval(Meta.parse(l)), exprs; init=0))


text:
"""
    function test_2_eval()
---
This is to test the `evaluation` function on two traits.
Only `base` genotypes were used.
The first 400 are for production.
The rest 200 are for challenge.
"""
function test_2_eval()
    @load "dat/run/base.jld2"
    qtl = sim_QTL(base, 500, 500)
    snp1 = base[:hap][:, 1:600] # used for production
    snp2 = base[:hap][:, 601:end] # used for challenge

    # bv and phenotypes for the production pop
    b1, p1 = phenotype(snp1, qtl[1], .5)

    # and also simulate bv and phenotypes as if challenged
    b3, p3 = phenotype(snp1, qtl[2], .5, .5)

    # bv and phenotypes for the challenge pop.
    b2, p2 = phenotype(base[:hap][:, 601:end], qtl[2], .5, .5)

    # SNP effects for these two pops.
    g1 = alleles2gt(snp1)       # genotypes of 0, 1, and 2
    g2 = alleles2gt(snp2)
    m1, s1 = snp_blup(g1, p1)
    m2, s2 = snp_blup(g2, p2)

    ebv1 = g1's1 .+ m1
    ebv2 = g2's2 .+ m2
    ebv3 = g1's2 .+ m2
    println(cor(ebv1, b1), ' ', cor(ebv1, p1))
    println(cor(ebv2, b2), ' ', cor(b2, p2))
    println(cor(ebv3, b1), ' ', cor(ebv3, p1))
end

"""
    function test_eval()
---
Create F₁ to test genomic seleciton methods.
"""
function test_eval(; just_one = false)
    #base = deserialize("dat/run/ns.ser")
    @load "dat/run/base.jld2"

    nSire, nDam, nSib = 200, 400, 42
    h² = [.5, .5]
    par = begin
        nqtl = repeat([100, 500, 1000], inner=3)
        nclg = repeat([32, 22, 2], 3)
        [nqtl nclg]
    end
    ped = random_mate(nSire, nDam, nSib)
    nxt = gdrop(base[:hap], ped, base[:r])

    rst = open(joinpath(dat_dir, "run/test-evaluation.txt"), "w")
    lhs = nothing
    for (nqtl, nclg) in eachrow(par)
        qtl = sim_QTL(base, [nqtl, nqtl]...)
        ip, ic, hp, hc = test_sets(nDam, nSib, nclg)
        bv₁, p₁ = phenotype(nxt, qtl[1], h²[1])
        bv₂, p₂ = phenotype(nxt, qtl[2], h²[2], 0.5) # 50% death

        # genotypes and phenotypes
        gp = alleles2gt(nxt[:, hp]) # genotypes
        gc = alleles2gt(nxt[:, hc])
        pp, pc = p₁[ip], p₂[ic] # phenotypes

        results = zeros(6)
        results[1:2] = [cor(bv₁, p₁), cor(bv₂, p₂)]

        # evaluation on phenotypes
        _, sp = snp_blup(gp, pp)
        gebv = gp'sp
        bv = bv₁[ip]
        results[3:4] = [cor(gebv, bv), cor(gebv, pp)]

        # evaluation on TBV
        _, sp = snp_blup(gp, bv)
        gebv = gp'sp
        results[5:6] = [cor(gebv, bv), cor(gebv, pp)]
        nprd = nSib - nclg
        print("$nqtl, $nprd: ", round.(results; digits = 3), "\n")
        print(rst, "$nqtl, $nprd: ", round.(results; digits = 3), "\n")

        just_one && break
    end
    close(rst)
end


text:
using Pkg
Pkg.activate(".")
try
    Pkg.instantiate()
catch
    Pkg.add("CSV"); Pkg.add("DataFrames"); Pkg.add("Flux"); Pkg.add("HTTP"); Pkg.add("JSON"); Pkg.add("Statistics")
end
println("\e[32;1mCompleted setup!\e[0m")


text:
type INIDRandomSequence{S<:ValueSupport,T<:UnivariateDistribution} <: MultivariateDistribution{S}
    distributions::Vector{T}
end
function INIDRandomSequence{T<:UnivariateDistribution}(distributions::Vector{T})
    S = super(T).parameters[2]
    return INIDRandomSequence{S,T}(distributions)
end

length(X::INIDRandomSequence) = length(X.distributions)

function insupport(X::INIDRandomSequence,x::AbstractVector)
  length(x)==length(X) && all(dx->insupport(dx[1],dx[2]),zip(X.distributions,x))
end

function _rand!{T<:Real}(X::INIDRandomSequence, x::AbstractVector{T})
  for (i,d_i) in enumerate(X.distributions)
    x[i] = rand(d_i)
  end
  return x
end

function _logpdf{T<:Real}(X::INIDRandomSequence, x::AbstractVector{T})
  return sum(arg->logpdf(arg[1],arg[2]),zip(X.distributions,x))
end

mean(X::INIDRandomSequence) = Float64[mean(d_i) for d_i in X.distributions]
var(X::INIDRandomSequence) = Float64[var(d_i) for d_i in X.distributions]
cov(X::INIDRandomSequence) = Diagonal(var(X))
entropy(X::INIDRandomSequence) = sum(entropy,X.distributions)

function mgf(X::IIDRandomSequence, t::AbstractVector)
  prod([mgf(d_i,t_i) for (d_i,t_i) in zip(X.distributions,t)])
end
function cf(X::IIDRandomSequence, t::AbstractVector)
  prod([cf(d_i,t_i) for (d_i,t_i) in zip(X.distributions,t)])
end


text:
"""
    hd2ae(h, d, p)

Equatorial to horizon coordinates:  transform hour angle and
declination to azimuth and altitude.

### Given ###

- `ha`: hour angle (local)
- `dec`: declination
- `phi`: site latitude

### Returned ###

- `az`: azimuth
- `el`: altitude (informally, elevation)

### Notes ###

1. All the arguments are angles in radians.

2. Azimuth is returned in the range 0-2pi;  north is zero, and east
   is +pi/2.  Altitude is returned in the range +/- pi/2.

3. The latitude phi is pi/2 minus the angle between the Earth's
   rotation axis and the adopted zenith.  In many applications it
   will be sufficient to use the published geodetic latitude of the
   site.  In very precise (sub-arcsecond) applications, phi can be
   corrected for polar motion.

4. The returned azimuth az is with respect to the rotational north
   pole, as opposed to the ITRS pole, and for sub-arcsecond
   accuracy will need to be adjusted for polar motion if it is to
   be with respect to north on a map of the Earth's surface.

5. Should the user wish to work with respect to the astronomical
   zenith rather than the geodetic zenith, phi will need to be
   adjusted for deflection of the vertical (often tens of
   arcseconds), and the zero point of the hour angle ha will also
   be affected.

6. The transformation is the same as Vh = Rz(pi)*Ry(pi/2-phi)*Ve,
   where Vh and Ve are lefthanded unit vectors in the (az,el) and
   (ha,dec) systems respectively and Ry and Rz are rotations about
   first the y-axis and then the z-axis.  (n.b. Rz(pi) simply
   reverses the signs of the x and y components.)  For efficiency,
   the algorithm is written out rather than calling other utility
   functions.  For applications that require even greater
   efficiency, additional savings are possible if constant terms
   such as functions of latitude are computed once and for all.

7. Again for efficiency, no range checking of arguments is carried
   out.
"""
function hd2ae(h, d, p)
    a = Ref{Cdouble}()
    e = Ref{Cdouble}()
    ccall((:eraHd2ae, liberfa), Cvoid,
          (Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}),
          h, d, p, a, e)
    return a[], e[]
end

"""
    hd2pa(h, d, p)

Parallactic angle for a given hour angle and declination.

### Given ###

- `ha`: hour angle
- `dec`: declination
- `phi`: site latitude

### Returned ###

Returns the parallactic angle

### Notes ###

1. All the arguments are angles in radians.

2. The parallactic angle at a point in the sky is the position
   angle of the vertical, i.e. the angle between the directions to
   the north celestial pole and to the zenith respectively.

3. The result is returned in the range -pi to +pi.

4. At the pole itself a zero result is returned.

5. The latitude phi is pi/2 minus the angle between the Earth's
   rotation axis and the adopted zenith.  In many applications it
   will be sufficient to use the published geodetic latitude of the
   site.  In very precise (sub-arcsecond) applications, phi can be
   corrected for polar motion.

6. Should the user wish to work with respect to the astronomical
   zenith rather than the geodetic zenith, phi will need to be
   adjusted for deflection of the vertical (often tens of
   arcseconds), and the zero point of the hour angle ha will also
   be affected.

### Reference ###

- Smart, W.M., "Spherical Astronomy", Cambridge University Press,
  6th edition (Green, 1977), p49.
"""
function hd2pa(h, d, p)
    ccall((:eraHd2pa, liberfa), Cdouble, (Cdouble, Cdouble, Cdouble), h, d, p)
end

"""
    hfk5z(rh, dh, date1, date2)

Transform a Hipparcos star position into FK5 J2000.0, assuming
zero Hipparcos proper motion.

### Given ###

- `rh`: Hipparcos RA (radians)
- `dh`: Hipparcos Dec (radians)
- `date1`, `date2`: TDB date (Note 1)

### Returned (all FK5, equinox J2000.0, date date1+date2) ###

- `r5`: RA (radians)
- `d5`: Dec (radians)
- `dr5`: FK5 RA proper motion (rad/year, Note 4)
- `dd5`: Dec proper motion (rad/year, Note 4)

### Notes ###

1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

   | `date1`   |     `date2` | Method      |
   |:----------|:------------|:------------|
   | 2450123.7 |         0.0 | JD          |
   | 2451545.0 |     -1421.3 | J2000       |
   | 2400000.5 |     50123.2 | MJD         |
   | 2450123.5 |         0.2 | date & time |

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

2. The proper motion in RA is dRA/dt rather than cos(Dec)*dRA/dt.

3. The FK5 to Hipparcos transformation is modeled as a pure rotation
   and spin;  zonal errors in the FK5 catalogue are not taken into
   account.

4. It was the intention that Hipparcos should be a close
   approximation to an inertial frame, so that distant objects have
   zero proper motion;  such objects have (in general) non-zero
   proper motion in FK5, and this function returns those fictitious
   proper motions.

5. The position returned by this function is in the FK5 J2000.0
   reference system but at date date1+date2.

6. See also [`fk52h`](@ref), [`h2fk5`](@ref), [`fk5hz`](@ref).

### Called ###

- [`s2c`](@ref): spherical coordinates to unit vector
- [`fk5hip`](@ref): FK5 to Hipparcos rotation and spin
- [`rxp`](@ref): product of r-matrix and p-vector
- [`sxp`](@ref): multiply p-vector by scalar
- [`rxr`](@ref): product of two r-matrices
- [`trxp`](@ref): product of transpose of r-matrix and p-vector
- [`pxp`](@ref): vector product of two p-vectors
- [`pv2s`](@ref): pv-vector to spherical
- [`anp`](@ref): normalize angle into range 0 to 2pi

### Reference ###

- F.Mignard & M.Froeschle, 2000, Astron.Astrophys. 354, 732-739.

"""
function hfk5z(rh, dh, date1, date2)
    r5 = Ref{Cdouble}()
    d5 = Ref{Cdouble}()
    dr5 = Ref{Cdouble}()
    dd5 = Ref{Cdouble}()
    ccall((:eraHfk5z, liberfa), Cvoid,
          (Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
          rh, dh, date1, date2, r5, d5, dr5, dd5)
    return r5[], d5[], dr5[], dd5[]
end

"""
    h2fk5(ra, dec, dra, ddec, px, rv)

Transform Hipparcos star data into the FK5 (J2000.0) system.

### Given (all Hipparcos, epoch J2000.0) ###

- `rh`: RA (radians)
- `dh`: Dec (radians)
- `drh`: Proper motion in RA (dRA/dt, rad/Jyear)
- `ddh`: Proper motion in Dec (dDec/dt, rad/Jyear)
- `pxh`: Parallax (arcsec)
- `rvh`: Radial velocity (km/s, positive = receding)

### Returned (all FK5, equinox J2000.0, epoch J2000.0) ###

- `r5`: RA (radians)
- `d5`: Dec (radians)
- `dr5`: Proper motion in RA (dRA/dt, rad/Jyear)
- `dd5`: Proper motion in Dec (dDec/dt, rad/Jyear)
- `px5`: Parallax (arcsec)
- `rv5`: Radial velocity (km/s, positive = receding)

### Notes ###

1. This function transforms Hipparcos star positions and proper
   motions into FK5 J2000.0.

2. The proper motions in RA are dRA/dt rather than
   cos(Dec)*dRA/dt, and are per year rather than per century.

3. The FK5 to Hipparcos transformation is modeled as a pure
   rotation and spin;  zonal errors in the FK5 catalog are not
   taken into account.

4. See also [`fk52h`](@ref), [`fk5hz`](@ref), [`hfk5z`](@ref).

### Called ###

- [`starpv`](@ref): star catalog data to space motion pv-vector
- [`fk5hip`](@ref): FK5 to Hipparcos rotation and spin
- [`rv2m`](@ref): r-vector to r-matrix
- [`rxp`](@ref): product of r-matrix and p-vector
- [`trxp`](@ref): product of transpose of r-matrix and p-vector
- [`pxp`](@ref): vector product of two p-vectors
- [`pmp`](@ref): p-vector minus p-vector
- [`pvstar`](@ref): space motion pv-vector to star catalog data

### Reference ###

- F.Mignard & M.Froeschle, Astron. Astrophys. 354, 732-739 (2000).

"""
function h2fk5(ra, dec, dra, ddec, px, rv)
    r = Ref{Cdouble}()
    d = Ref{Cdouble}()
    dr = Ref{Cdouble}()
    dd = Ref{Cdouble}()
    p = Ref{Cdouble}()
    v = Ref{Cdouble}()
    ccall((:eraH2fk5, liberfa), Cvoid,
            (Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Cdouble, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}, Ref{Cdouble}),
            ra, dec, dra, ddec, px, rv, r, d, dr, dd, p, v)
    return r[], d[], dr[], dd[], p[], v[]
end


text:
module SINDYPI_Julia

include("SINDYPI_main.jl")
include("Basis_generation.jl")
include("SparsityAlgorithm.jl")
include("Table_generation.jl")

export  SINDYPI
export Basis_SINDY_PI
export Sparse_matrix_SINDYPI
export Table_SINDY_PI

end # module


text:
import LinearAlgebra

mutable struct MatMul <: DiffableFunction
    grad_field::GradField
end

function forward(::MatMul, x, W)
    return x * W
end

function backward(f::MatMul, gy)
    x, W = f.grad_field.inputs
    gx = matmul(gy, transpose(W))
    gW = matmul(transpose(x), gy)
    return gx, gW
end

"""
    matmul(x, W)
Matrix multiplication.

# Example

```julia-repl
julia> x = Variable([1 1; 0 1])
name: nothing
values: [1 1; 0 1]
creator: User-Defined(nothing)

julia> W = Variable([1 0; 1 1])
name: nothing
values: [1 0; 1 1]
creator: User-Defined(nothing)

julia> matmul(x, W)
name: nothing
values: [2 1; 1 1]
creator: MatMul
```
"""
matmul(x, W) = MatMul(GradField())(x, W)

Base.:*(x::Variable{Matrix{T}}, W::Variable{Matrix{R}}) where {T <: Real, R <: Real} = matmul(x, W)


text:
table1 = [
  1 2
  0 3
]
table1 = table1[end:-1:1, :]
table2 = [
  15 12 11 10
  14 13  8  9
   1  2  7  6
   0  3  4  5
]
table2 = table2[end:-1:1, :]

## Test encode
zz = Set(Int64[])
for x in 0:15
  for y in 0:15
    z = encode_hilbert_zero(x, y)
    # @show (x, y, z)
    @test z ∉ zz
    @test 0 ≤ z < 256
    push!(zz, z)
  end
end

## Test decode
table2ind = CartesianIndices(table2)
for z in 0:15
  x, y = decode_hilbert_zero(z)
  found = table2ind[table2 .== z][1]
  # @show z
  # @show found
  # @show (x, y)
  @test found[2] - 1 == x
  @test found[1] - 1 == y
end

xy = Set(Tuple{Int64, Int64}[])
last = [-1, 0]
for z in 0:255
  x, y = decode_hilbert_zero(z)
  @test (x, y) ∉ xy
  push!(xy, (x, y))
  # The next point in the grid must be beside the previous.
  @test ((x - last[1])^2 == 1) || ((y - last[2])^2 == 1)
  last = [x, y]
end
for i in 0:15
  for j in 0:15
    @test (i, j) in xy
  end
end

## Test against each other.

for z in 1:400
  x, y = decode_hilbert_zero(z)
  zz = encode_hilbert_zero(x, y)
  @test zz == z
end


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

using ArraysOfArrays
using Test
using StatsBase
using Statistics

@testset "StatsBase support" begin
    r = rand(1,100)
    V = VectorOfSimilarVectors{Float64}(r)
    w = FrequencyWeights(rand(100))
    w0 = FrequencyWeights(vec(ones(100,1)))
    @test isapprox(@inferred(sum(V, w))[1], sum(r, w))
    @test isapprox(@inferred(mean(V,w))[1], mean(r, w))
    @test isapprox(@inferred(var(V,w, corrected=true))[1], var(r, w, corrected=true))
    @test isapprox(@inferred(std(V,w, corrected=true))[1], sqrt(var(V,w, corrected=true)[1]))
    @test isapprox(@inferred(std(V, corrected=true))[1], sqrt(var(V, w0, corrected=true)[1]))
    @test isapprox(@inferred(std(V, corrected=true))[1], std(V, w0, corrected=true)[1])
    @test isapprox(@inferred(cov(V, w))[1], var(r, w, corrected=true))
    @test isapprox(@inferred(cor(V,w))[1], 1.0)
end


text:
module Scheme

export lex, Obj, head, tail, displayobj, tokenize, cons, islist, dsplo, ptol, tokentoobj, ltop, tokenstolist, sread, repl, evaluate

#typestuff
struct Token
    type
    str
end

struct Sym
    name
end

eq(x::Sym, y::Sym) = x.name == y.name
eq(x, y) = x == y

struct Primitiveop
    f
end

struct Compositeop
    arglist
    body
end

const D = Union{Int, Float64, String, SubString{String}, Bool, Primitiveop, Compositeop, Sym, Nothing}

const Obj = Union{D, Tuple}

#list related functions
function islist(x::Obj)
    if typeof(x) <: Tuple
        if x[2] == nothing
            return true
        else
            return islist(x[2])
        end
    else
        return false
    end
end

function ptol(x::Obj)
    if typeof(x) == Nothing
        return []
    elseif typeof(x) <: D
        return x
    elseif islist(x)
        return vcat(x[1], ptol(x[2]))
    else
        print("cannot turn non-list nested pairs into list ")
    end
end

function ltop(x)
    y = reverse(x)
    r = nothing
    for i in y
        if typeof(i) <: Array
            r = cons(ltop(i), r)
        else
            r = cons(i, r)
        end
    end
    return r
end

function fevaluate end

ltop(x::D) = x

#builtins
cons(a, b) = (a, b)
car(x) = x[1]
cdr(x) = x[2]
caar(x) = car(car(x))
cadr(x) = car(cdr(x))
cddr(x) = cdr(cdr(x))
equal(x, y) = x == y
negate(x) = typeof(x) <: Bool ? !(x) : println("Object not boolean")
#ifthenelse(cond, tr, fal) = cond == true ? tr : fal
list(x...) = ltop(x)

function loadcont(x)
    io = open(x, "r")
    str = read(io, String)
    close(io)
    return str
end

function load(x)
    str = loadcont(x)
    objs = map(ltop, tokenstolist(lex(str)))
    map(fevaluate, objs)
    return nothing
end

symbtabel = Dict{Any, Any}(
 "cons" => Primitiveop(cons),
 "car" => Primitiveop(car),
 "cdr" => Primitiveop(cdr),
 "caar" => Primitiveop(caar),
 "cadr" => Primitiveop(cadr),
 "cddr" => Primitiveop(cddr),
 "=" => Primitiveop(equal),
 "!" => Primitiveop(negate),
 #"if" => Primitiveop(ifthenelse),
 "list" => Primitiveop(list),
 "+" => Primitiveop(+),
 "-" => Primitiveop(-),
 "*" => Primitiveop(*),
 "eval" => Primitiveop(fevaluate),
 "args" => Primitiveop(x->ltop(x.arglist)),
 "body" => Primitiveop(x->x.body),
 "load" => Primitiveop(load),
 "loadcont" => Primitiveop(loadcont),
 "type" => Primitiveop(x->print(typeof(x)))
 )

#evaluation logic
function fapply(func, args)
    if typeof(func) <: Primitiveop
        return func.f(args...)

    elseif typeof(func) <: Compositeop
        argbinds = Dict()
        i = 0
        for arg in func.arglist
            i += 1
            merge!(argbinds, Dict(arg.name=>args[i]))
        end

        fevaluate(func.body; exenv = argbinds)
    else
        print("error, first element is not callable")
    end
end

function fevaluate(x; exenv = Dict())
    if typeof(x) <: Sym
        return haskey(exenv, x.name) ? exenv[x.name] : symbtabel[x.name]
    elseif typeof(x) <: D
        return x
    elseif (islist(x))
        #binding symbols to values and creating functions
        if typeof(car(x)) <: Sym && car(x).name == "def"
            if typeof(cadr(x)) <: Sym
                merge!(symbtabel, Dict(car(cdr(x)).name => fevaluate(car(cddr(x)), exenv = exenv)))
                return fevaluate(car(cddr(x)), exenv = exenv)
            elseif islist(cadr(x))
                k = cadr(x)
                name = car(k).name
                args = cdr(k)
                body = car(cdr(cdr(x)))
                merge!(symbtabel, Dict(name => Compositeop(ptol(args), body)))
                return Compositeop(ptol(args), body)
            end

        #creating functions
        elseif typeof(car(x)) <: Sym && car(x).name == "lambda"
            return Compositeop(ptol(car(cdr(x))), car(cddr(x)))

        elseif typeof(car(x)) <: Sym && car(x).name == "quote"
            return cadr(x)

        elseif typeof(car(x)) <: Sym && car(x).name == "if"
            if fevaluate(cadr(x), exenv = exenv) == true
                return fevaluate(car(cddr(x)), exenv = exenv)
            else
                return fevaluate(cadr(cddr(x)), exenv = exenv)
            end
        else
            return fapply(fevaluate(car(x), exenv = exenv), map(i -> fevaluate(i, exenv = exenv), ptol(cdr(x))))
        end
    else
        println("Evaluation requires proper lists")
    end
end

function parenmatch(str)
    counter = 0
    for l in str
        if l == '('
            counter += 1
        elseif l == ')'
            counter -= 1
        end

        if counter < 0
            return false
        end
    end
    counter == 0 ? true : false
end

function tokenize(str)
    tp = ""
    match(r"\".*\"", str) != nothing && match(r"\".*\"", str).match == str && (tp = "string")
    match(r"[0-9]+\.[0-9]+", str) != nothing && match(r"[0-9]+\.?[0-9]*", str).match == str && (tp = "float")
    match(r"[0-9]+", str) != nothing && match(r"[0-9]+", str).match == str && (tp = "int")
    (str == "true" || str == "false") && (tp = "bool")
    str == "(" && (tp = "lparen")
    str == ")" && (tp = "rparen")
    tp == "" ? tp = "id" : tp = tp
    return Token(tp, str)
end

function lex(inp)
    inp = replace(inp, "(" => " ( ")
    inp = replace(inp, ")" => " ) ")
    inp = split(inp, r"[\h, \v]+")
    inp = filter(i -> length(i) > 0, inp)
    tks = map(tokenize, inp)
end

function atom(x::Token)
    if x.type == "string"
        return String(x.str[2:end-1])
    elseif x.type == "float"
        return parse(Float64, x.str)
    elseif x.type == "int"
        return parse(Int, x.str)
    elseif x.type == "bool"
        return parse(Bool, x.str)
    elseif x.type == "id"
        return Sym(x.str)
    end
end

function recurreplacequote(lis)
    newlist = []
    templist = []
    i = 1
    while i != length(lis)+1
        if typeof(lis[i]) <: Sym && lis[i].name == "'"
            push!(templist, Sym("quote"))
            push!(templist, recurreplacequote(lis[i += 1]))
            push!(newlist, templist)
            i += 1
        elseif typeof(lis[i]) <: Array
            push!(newlist, recurreplacequote(lis[i]))
            i += 1
        else
            push!(newlist, lis[i])
            i += 1
        end
    end
    return newlist
end

function tokenstolist(tokens)
    rstack = []
    objstack = []
    for tok in tokens
        if tok.str == "("
            push!(objstack, [])
        elseif tok.str == ")"
            ob = pop!(objstack)
            objstack == [] ? push!(rstack, ob) : push!(objstack[end], ob)
        else
            objstack == [] ? push!(rstack, atom(tok)) : push!(objstack[end], atom(tok))
        end
    end
    return recurreplacequote(rstack)
end

function readexpression()
    buffer = ""

    while true
        buffer = buffer * readline();
        !(parenmatch(buffer)) || break
    end
    return buffer
end

sread() = map(ltop, tokenstolist(lex(readexpression())))

function displayobj(x::Obj)
    if typeof(x) <: Sym
        print(x.name)
    elseif typeof(x) <: Compositeop
        print("args: ")
        displayobj(ltop(x.arglist))
        println()
        print("body: ")
        displayobj(x.body)
    elseif typeof(x) <: D
        x == nothing ? print("()") :  print(x)
    elseif islist(x)
        print("(")
        while x != nothing
            #if islist(x[1])
            #    displayobj(x[1])
            #    x[2] != nothing && print(" ")
            #    x = x[2]
            #else
                displayobj(x[1])
                x[2] != nothing && print(" ")
                x = x[2]
            #end
        end
        print(")")
    elseif typeof(x) <: Tuple
        print("(")
        displayobj(x[1])
        print(" . ")
        displayobj(x[2])
        print(")")
    end
end

dsplo(x::Obj) = begin
    displayobj(x)
    println()
    end

function repl()
    while true
        print(">")
        #objs = map(evaluate, sread())
        #objs = sread()
        objs = map(fevaluate, sread())
        map(dsplo, objs)
        println()
    end
end
end

max_stars_repo_name: carl-gjerris/pseudo-scheme-interpreter

max_issues_repo_name: carl-gjerris/pseudo-scheme-interpreter

max_forks_repo_name: carl-gjerris/pseudo-scheme-interpreter


text:
using IntervalArithmetic, StaticArrays, IntervalLinearAlgebra, IntervalRootFinding

A = @SMatrix [4..6 -1..1 -1..1 -1..1;-1..1 -6.. -4 -1..1 -1..1;-1..1 -1..1 9..11 -1..1;-1..1 -1..1 -1..1 -11.. -9]

b = @SVector [-2..4, 1..8, -4..10, 2..12]

jac = Jacobi()
gs = GaussSeidel()
hbr = HansenBliekRohn()

@btime solve($A, $b, $gs)

@btime solve($A, $b, $jac)

@btime solve($A, $b, $hbr)

# comparison with IntervalRootFinding

@btime gauss_seidel_interval($A, $b)
@btime gauss_seidel_contractor($A, $b) #NOTE: THIS IS THE JACOBI method
@btime gauss_elimination_interval($A, $b)


text:
"""
	grdmask(cmd0::String="", arg1=nothing, kwargs...)

1. It reads one or more pathfiles that each define a closed polygon.
2. The pathfiles simply represent data point locations and the mask is set to the inside or outside
   value depending on whether a node is within a maximum distance from the nearest data point.

Full option list at [`grdmask`]($(GMTdoc)grdmask.html)

Parameters
----------

- $(GMT.opt_R)
- **I** | **inc** :: [Type => Str | Number]

    *x_inc* [and optionally *y_inc*] is the grid spacing.
    ($(GMTdoc)grdmask.html#i)
- **A** | **steps** | **straight_lines** :: [Type => Str | Number]		``Arg = m|p|x|y``

    If the input data are geographic then the sides in the polygons will be approximated by great circle arcs.
    When using this option sides will be regarded as straight lines.
    ($(GMTdoc)grdmask.html#a)
- **G** | **outgrid** :: [Type => Str]

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdmask(....) form.
    ($(GMTdoc)grdmask.html#g)
- **N** | **out_edge_in** :: [Type => Str | List]    ``Arg = [z|Z|p|P]values``

    Sets the out/edge/in that will be assigned to nodes that are outside the polygons, on the edge, or inside.
    Values can be any number, including the textstring NaN [Default is 0/0/1].
    ($(GMTdoc)grdmask.html#n)
- **S** | **search_radius** :: [Type => Str | List]    ``Arg = search_radius[unit] |xlim/ylim``

    Set nodes to inside, on edge, or outside depending on their distance to the nearest data point.
- $(GMT.opt_V)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_g)
- $(GMT.opt_h)
- $(GMT.opt_i)
- $(GMT.opt_j)
- $(GMT.opt_n)
- $(GMT.opt_r)
- $(GMT.opt_x)
- $(GMT.opt_swap_xy)
"""
function grdmask(cmd0::String="", arg1=nothing; kwargs...)

	length(kwargs) == 0 && occursin(" -", cmd0) && return monolitic("grdmask", cmd0, arg1)

	d = KW(kwargs)
	cmd = parse_common_opts(d, "", [:R :V_params :a :e :f :g :j :n :yx :r :x :params])
	cmd = parse_these_opts(cmd, d, [[:A :steps :straight_lines], [:I :inc], [:G :outgrid],
	                                [:N :out_edge_in], [:S :search_radius]])

	return common_grd(d, "grdmask " * cmd, arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdmask(arg1, cmd0::String=""; kw...) = grdmask(cmd0, arg1; kw...)


text:
"""
Abstract Factor type
"""
abstract type Factor end

"""
Factor type defined by energy function

`` \\psi_{a}(x_{\\partial a}) = exp[ -E(x_{\\partial a})]``

Parameters:
* `E::Function`: ``E(x_{\\partial a}) ``
* `idx::Vector{Int64}`: ``\\partial a`` indeces

"""
struct EnergyFun <: Factor
    E::Function
    idx::Vector{Int64}
end

"""
Factor type for Ising model

`` \\psi_{ij}(x_{i},x_{j}) = exp[ Jx_{i}x_{j} + h_{i}x_{i} +h_{j}x_{j}]``

Parameters:
* `J::Float64`: coupling
* `hi::Float64`: field on ``x_{i}``
* `hi::Float64`: field on ``x_{j}``
* `idx::Vector{Int64}`: ``(i,j)``

moments are computed analytically

"""
struct IsingPair <: Factor
    J::Float64
    hi::Float64
    hj::Float64
    idx::Vector{Int64}
	IsingPair(J,hi,hj,idx) = length(idx)!=2 ? error("idx must have length 2") : new(J,hi,hj,idx)
end

"""
Basic constructor for Ising Factor Graph

Input:
* `N::Int64`: number of spins
* `h::Vector{Float64}`: Array of fields ``\\in \\mathbb{R}^{N}``
* `J::SparseMatrixCSC{Float64,Int64}`: symmetric coupling matrix ``\\in \\mathbb{R}^{N\\times N}``
* `t::Symbol`: choose between `EnergyFun` and `IsingPair`

Output: `Array{Factor}`: array of factor nodes `` \\{ \\psi_{ij} \\}``

"""
function IsingFG(N::Int64,h::Vector{Float64},J::SparseMatrixCSC{Float64,Int64}, t::Symbol; Adj::SparseMatrixCSC{Int64,Int64} = (J.!=0)*1)
	t ∉ (:EnergyFun,:IsingPair) && error("invalid option")
	degree = sum(Adj,dims = 2)[:]
	if t == :EnergyFun
		return collect([EnergyFun((x::Vector{Int64}-> -( J[i,j]*x[1]*x[2]+x[1]*(h[i]/degree[i])+ x[2]*(h[j]/degree[j]))),[i,j] ) for i=1:N, j=1:N if Adj[i,j]!=0 && i<j])
	elseif t == :IsingPair
		return collect([IsingPair(J[i,j],h[i]/degree[i],h[j]/degree[j],[i;j] ) for i=1:N, j=1:N if Adj[i,j]!=0 && i<j])
	end
end

"""
DCState

Istantaneous state of Density Consistency run

Fields:
* `μtl :: Vector{Vector}`		 array of first moments of tilted distributions (one for each factor)
* `Σtl :: Vector{Matrix}`		 array of second moments of tilted distributions (one for each factor)
* `μt  :: Vector`				 mean vector w.r.t tilted distribution (size N)
* `Σt  :: Matrix`				 correlation matrix w.r.t tilted distribution (size NxN)
* `µ   :: Vector`				 mean vector of gaussian distribution (size N)
* `Σ   :: Matrix`				 correlation matrix of the gaussian distribution  (size NxN)
* `h   :: Vector{Vector}`		 array of gaussian field vector (one for each factor)
* `S   :: Vector{Matrix}`		 array of inverse correlation matrix (one for each factor)
* `yc  :: Vector{Vector}`		 array of cavity fields vector (one for each factor)
* `Sc  :: Vector{Matrix}`		 array of cavity inverse correlation matrix (one for each factor)
"""
struct DCState
    μtl :: Vector{Vector}
    Σtl :: Vector{Matrix}
    μt  :: Vector
    Σt  :: Matrix
    µ   :: Vector
    Σ   :: Matrix
    h   :: Vector{Vector}
    S   :: Vector{Matrix}
    yc  :: Vector{Vector}
    Sc  :: Vector{Matrix}
end

function DCState(N::Int,m::Vector{Int};h0::Float64 = 1e-3,f::Function=n->h0*(rand(n) .- 0.5))
    M=length(m)
    return DCState([zeros(m[a]) for a=1:M],
                   [zeros(m[a],m[a]) for a=1:M],
                   zeros(N),
                   zeros(N,N),
				   zeros(N),
                   zeros(N,N),
                   [f(m[a]) for a=1:M],
                   [eye(m[a]) for a=1:M],
                   [zeros(m[a]) for a=1:M],
                   [zeros(m[a],m[a]) for a=1:M])
end

max_stars_repo_name: giovact/DensityConsistency

max_issues_repo_name: giovact/DensityConsistency

max_forks_repo_name: giovact/DensityConsistency


text:
module Covid_Forecast

# Write your package code here.
using Pkg

Pkg.add("Conda")
Pkg.add("webio_jupyter_extension")
Pkg.add("")
Pkg.add("")
Pkg.add("")

#using Conda

using WebIO

 #https://juliagizmos.github.io/WebIO.jl/latest/providers/ijulia/

#using Interact
#using CSV
#using DataFrames
#Conda.pip_interop(true)
#Conda.pip("install", "webio_jupyter_extension")

#using IJulia
#notebook()

end


text:
dataframes_installed = haskey(Pkg.installed(), "DataFrames")

if dataframes_installed
  using DataFrames
end

using PairwiseListMatrices
using Base.Test

print("""

Test PairwiseListMatrices
=======================
""")

let list = [1,-2,3],
    list_diag_triu =  [ 1 -2
                        0 3 ],
    list_diag_tril =  [ 1 0
                       -2 3 ],
    list_triu = [ 0 1 -2
                  0 0 3
                  0 0 0 ],
    list_tril = [ 0 0 0
                  1 0 0
                 -2 3 0 ],
    list_diag_sym = [ 1 -2
                     -2 3 ],
    list_sym = [ 0 1 -2
                 1 0 3
                -2 3 0 ]

  print("""

  Test creations from list
  ------------------------
  """)

  pld_triu = UpperTriangular(PairwiseListMatrix(list, true))
  pld_tril = LowerTriangular(PairwiseListMatrix(list, true))
  pl_triu  = UpperTriangular(PairwiseListMatrix(list))
  pl_tril  = LowerTriangular(PairwiseListMatrix(list))
  pld_sym  = PairwiseListMatrix(list, true)
  pl_sym   = PairwiseListMatrix(list)

  @test pld_triu == list_diag_triu
  @test pld_tril == list_diag_tril
  @test pl_triu  == list_triu
  @test pl_tril  == list_tril
  @test pld_sym  == list_diag_sym
  @test pl_sym   == list_sym

  print("""

  Getters
  -------
  """)

  @test getlist(pld_sym) == pld_sym.list
  @test getlist(pl_sym)  ==  pl_sym.list

  @test getdiag(pld_sym) == pld_sym.diag
  @test getdiag(pl_sym)  ==  pl_sym.diag

  print("""

  Length of the list
  """)

  @test lengthlist(2, true)  == lengthlist(pld_sym)
  @test lengthlist(3, false) == lengthlist(pl_sym)

  print("""

  Convert
  -------
  """)

  @test list_diag_sym == convert(PairwiseListMatrix{Int, true},  list_diag_sym)
  @test list_diag_sym == convert(PairwiseListMatrix{Int, false}, list_diag_sym)

  @test convert(Matrix{Float64}, list_sym)  == convert(PairwiseListMatrix{Float64, true},  list_sym)
  @test convert(Matrix{Float64}, list_sym)  == convert(PairwiseListMatrix{Float64, false}, list_sym)

  @test convert(Matrix{Int}, pl_sym)  == list_sym
  @test convert(Matrix{Int}, pld_sym) == list_diag_sym

  @test convert(PairwiseListMatrix{Int, false}, pld_sym) == pld_sym
  @test convert(PairwiseListMatrix{Int, true}, pl_sym)   == pl_sym

  print("""

  Test unary operations
  ---------------------
  """)

  @test -pld_triu == -list_diag_triu
  @test -pld_tril == -list_diag_tril
  @test -pl_triu  == -list_triu
  @test -pl_tril  == -list_tril
  @test -pld_sym  == -list_diag_sym
  @test -pl_sym   == -list_sym

  @test abs(pld_triu) == abs(list_diag_triu)
  @test abs(pld_tril) == abs(list_diag_tril)
  @test abs(pl_triu ) == abs(list_triu)
  @test abs(pl_tril ) == abs(list_tril)
  @test abs(pld_sym ) == abs(list_diag_sym)
  @test abs(pl_sym  ) == abs(list_sym)

  print("""

  Test binary operations
  ----------------------
  """)

  @test pld_triu - pld_triu == list_diag_triu - list_diag_triu
  @test pld_tril - pld_tril == list_diag_tril - list_diag_tril
  @test pl_triu  - pl_triu  == list_triu - list_triu
  @test pl_tril  - pl_tril  == list_tril - list_tril
  @test pld_sym  - pld_sym  == list_diag_sym - list_diag_sym
  @test pl_sym   - pl_sym   == list_sym - list_sym

  @test pld_triu + pld_triu == list_diag_triu + list_diag_triu
  @test pld_tril + pld_tril == list_diag_tril + list_diag_tril
  @test pl_triu  + pl_triu  == list_triu + list_triu
  @test pl_tril  + pl_tril  == list_tril + list_tril
  @test pld_sym  + pld_sym  == list_diag_sym + list_diag_sym
  @test pl_sym   + pl_sym   == list_sym + list_sym

  @test pld_triu .* pld_triu == list_diag_triu .* list_diag_triu
  @test pld_tril .* pld_tril == list_diag_tril .* list_diag_tril
  @test pl_triu  .* pl_triu  == list_triu .* list_triu
  @test pl_tril  .* pl_tril  == list_tril .* list_tril
  @test pld_sym  .* pld_sym  == list_diag_sym .* list_diag_sym
  @test pl_sym   .* pl_sym   == list_sym .* list_sym

  @test pld_sym  .* 2 == list_diag_sym .* 2
  @test pl_sym   .* 2 == list_sym .* 2
  @test pld_sym  .+ 2 == list_diag_sym .+ 2
  @test pl_sym   .+ 2 == list_sym .+ 2
  @test pld_sym  .- 2 == list_diag_sym .- 2
  @test pl_sym   .- 2 == list_sym .- 2
  @test pld_sym  ./ 2 == list_diag_sym ./ 2
  @test pl_sym   ./ 2 == list_sym ./ 2
  @test pld_sym  .* 2 == list_diag_sym .* 2
  @test pl_sym   + 2 == list_sym + 2
  @test pld_sym  + 2 == list_diag_sym + 2
  @test pl_sym   - 2 == list_sym - 2
  @test pld_sym  - 2 == list_diag_sym - 2
  @test pl_sym   / 2 == list_sym / 2
  @test pld_sym  / 2 == list_diag_sym / 2

  print("""

  Transpose
  ---------
  """)

  @test transpose(pld_triu) == pld_tril == list_diag_triu.'
  @test transpose(pld_tril) == pld_triu == list_diag_tril.'
  @test transpose(pl_triu) == pl_tril == list_triu.'
  @test transpose(pl_tril) == pl_triu == list_tril.'

  @test ctranspose(pld_triu) == pld_tril == list_diag_triu'
  @test ctranspose(pld_tril) == pld_triu == list_diag_tril'
  @test ctranspose(pl_triu) == pl_tril == list_triu'
  @test ctranspose(pl_tril) == pl_triu == list_tril'

  @test transpose(pld_sym) == pld_sym == list_diag_sym.'
  @test transpose(pl_sym) == pl_sym == list_sym.'

  @test ctranspose(pld_sym) == pld_sym == list_diag_sym'
  @test ctranspose(pl_sym) == pl_sym == list_sym'

  @test transpose!(pl_sym) == pl_sym
  @test ctranspose!(pl_sym) == pl_sym

  print("""

  Linear algebra
  --------------
  """)

  @test pld_triu * pld_triu == list_diag_triu * list_diag_triu
  @test pld_tril * pld_tril == list_diag_tril * list_diag_tril
  @test pl_triu  * pl_triu  == list_triu      * list_triu
  @test pl_tril  * pl_tril  == list_tril      * list_tril
  @test pld_sym  * pld_sym  == list_diag_sym  * list_diag_sym
  @test pl_sym   * pl_sym   == list_sym       * list_sym
  @test Symmetric(pl_sym) * Symmetric(pl_sym) == Symmetric(list_sym) * Symmetric(list_sym)

  @test pld_triu / pld_triu == list_diag_triu / list_diag_triu
  @test pld_tril / pld_tril == list_diag_tril / list_diag_tril
  @test pld_sym  / pld_sym  == list_diag_sym  / list_diag_sym
  @test pl_sym   / pl_sym   == list_sym       / list_sym

  @test svd(pld_triu) == svd(list_diag_triu)
  @test svd(pld_tril) == svd(list_diag_tril)
  @test svd(pl_triu ) == svd(list_triu)
  @test svd(pl_tril ) == svd(list_tril)
  @test svd(pld_sym ) == svd(list_diag_sym)
  @test svd(pl_sym  ) == svd(list_sym)

  print("""

  Stats
  -----
  """)

  @test mean(pld_triu) == mean(list_diag_triu)
  @test mean(pld_tril) == mean(list_diag_tril)
  @test mean(pl_triu ) == mean(list_triu)
  @test mean(pl_tril ) == mean(list_tril)
  @test mean(pld_sym ) == mean(list_diag_sym)
  @test mean(pl_sym  ) == mean(list_sym)

  @test mean(pld_sym,1) == mean(list_diag_sym, 1)
  @test mean(pl_sym ,1) == mean(list_sym, 1)

  @test mean(pld_sym,2) == mean(list_diag_sym, 2)
  @test mean(pl_sym ,2) == mean(list_sym, 2)

  @test std(pld_triu) == std(list_diag_triu)
  @test std(pld_tril) == std(list_diag_tril)
  @test std(pl_triu ) == std(list_triu)
  @test std(pl_tril ) == std(list_tril)
  @test std(pld_sym ) == std(list_diag_sym)
  @test std(pl_sym  ) == std(list_sym)

  @test std(pld_sym,1) == std(list_diag_sym, 1)
  @test std(pl_sym, 1) == std(list_sym, 1)

  @test std(pld_sym,2) == std(list_diag_sym, 2)
  @test std(pl_sym, 2) == std(list_sym, 2)

  @test cor(pld_sym ) == cor(list_diag_sym)
  @test cor(pl_sym  ) == cor(list_sym)

  print("""

  triu
  ----
  """)

  @test triu(pld_sym ) == triu(list_diag_sym)
  @test triu(pl_sym  ) == triu(list_sym)

  @test triu(pld_sym,1) == triu(list_diag_sym, 1)
  @test triu(pl_sym, 1) == triu(list_sym, 1)

  @test_throws ErrorException triu!(pl_sym)
  @test_throws ErrorException triu!(pl_sym, 1)

end

let list = PairwiseListMatrix([true, true, true])

  @test list - true == [ -1  0  0
                         0 -1  0
                         0  0 -1 ]
  @test list + true == [ 1  2  2
                         2  1  2
                         2  2  1 ]
end

print("""

List with Labels
----------------
""")

let list = [1,2,3], labels=['a', 'b', 'c'], labels_diag = ['A', 'B'],
    list_diag_triu =  UpperTriangular(PairwiseListMatrix(list, labels_diag, true)),
    list_diag_tril =  LowerTriangular(PairwiseListMatrix(list, labels_diag, true)),
    list_triu = UpperTriangular(PairwiseListMatrix(list, labels)),
    list_tril = LowerTriangular(PairwiseListMatrix(list, labels)),
    list_diag_sym = PairwiseListMatrix(list, labels_diag, true),
    list_sym = PairwiseListMatrix(list, labels)

  print("""
  labels(!) & copy
  """)

  @test PairwiseListMatrices.labels(list_diag_sym) == labels_diag

  copy_lm = copy(list_sym)
  labels!(copy_lm, ['A', 'B', 'C'])
  @test PairwiseListMatrices.labels(copy_lm) == ['A', 'B', 'C']

  @test PairwiseListMatrices.labels(list_sym) == labels

  @test diag(list_diag_sym) == [1, 3]
  @test diag(list_sym) == [0, 0, 0]

  print("""
  getlabel & getindex
  """)

  @test getlabel(list_diag_sym, 'A', 'B') == list_diag_sym[1,2] == 2
  @test getlabel(list_diag_sym, 'B', 'A') == list_diag_sym[2,1] == 2

  @test getlabel(list_sym, 'a', 'b') == list_sym[1,2] == 1
  @test getlabel(list_sym, 'b', 'a') == list_sym[2,1] == 1

  print("""
  setlabel! & setindex!
  """)

  list_diag_sym[1,2] = 30
  @test list_diag_sym[1,2] == 30
  list_diag_sym[2] = 40
  @test list_diag_sym[2] == 40

  setlabel!(list_diag_sym, 20, 'A', 'B')
  @test list_diag_sym[1,2] == 20
  setlabel!(list_diag_sym, 10, 'B', 'A')
  @test list_diag_sym[2,1] == 10

  setlabel!(list_sym, 20, 'a', 'b')
  @test list_sym[1,2] == 20
  setlabel!(list_sym, 10, 'b', 'a')
  @test list_sym[2,1] == 10

  list_diag_sym[1,2] = 10
  @test list_diag_sym[2,1] == 10

  list_sym[1,2] = 15
  @test list_sym[2,1] == 15
  list_sym[2] = 30
  @test list_sym[2] == 30

  list_diag_sym[1,2] = 15
  @test list_diag_sym[2,1] == 15
  list_diag_sym[2] = 30
  @test list_diag_sym[2] == 30

end

print("""

Vector{PairwiseListMatrix}
--------------------------
""")

let list_samples = [ PairwiseListMatrix(Int, 4, false) for i in 1:100 ],
  list_diag_samples = [ PairwiseListMatrix(Int, 4, true) for i in 1:100 ],
  full_samples = Matrix{Int}[ full(mat) for mat in list_samples ],
  full_diag_samples = Matrix{Int}[ full(mat) for mat in list_diag_samples ]

  @test sum(list_samples) == sum(full_samples)
  @test sum(list_diag_samples) == sum(full_diag_samples)

  @test mean(list_samples) == mean(full_samples)
  @test mean(list_diag_samples) == mean(full_diag_samples)
end

@test sum(PairwiseListMatrix{Int, false}[ PairwiseListMatrix([1, 1, 1]) ]) == PairwiseListMatrix([1, 1, 1])
@test sum(PairwiseListMatrix{Int, true}[ PairwiseListMatrix([1, 1, 1], true) ]) == PairwiseListMatrix([1, 1, 1], true)

@test_throws ErrorException sum(PairwiseListMatrix{Int, true}[])
@test_throws ErrorException sum(PairwiseListMatrix{Int, true}[ PairwiseListMatrix([1, 1, 1], true), PairwiseListMatrix([1, 1, 1, 1, 1, 1], true) ])
@test_throws ErrorException std(PairwiseListMatrix{Int, true}[ PairwiseListMatrix([1, 1, 1], true) ])

let list_samples = [ ones(PairwiseListMatrix(Int, 4, false)) for i in 1:100 ],
  list_diag_samples = [ zeros(PairwiseListMatrix(Int, 4, true)) for i in 1:100 ]

  @test std(list_samples) == std(list_diag_samples)
end

print("""

Mean without diagonal
---------------------
""")

let list = PairwiseListMatrix([10,20,30], ["A", "B", "C"])

  @test vec( mean_nodiag(list, 1) ) == [15., 20., 25.]
  @test vec( mean_nodiag(list, 2) ) == [15., 20., 25.]
  @test mean_nodiag(list) == 20.
end

let list = PairwiseListMatrix([0,10,20,0,30,0], ["A", "B", "C"], true)

  @test vec( mean_nodiag(list, 1) ) == [15., 20., 25.]
  @test vec( mean_nodiag(list, 2) ) == [15., 20., 25.]
  @test mean_nodiag(list) == 20.
end

print("""

To and from table
-----------------
""")

let table = [ "A" "B" 10
              "A" "C" 20
              "B" "C" 30 ],
  list = PairwiseListMatrix([10,20,30], ["A", "B", "C"]),
  list_diag = PairwiseListMatrix([0,10,20,0,30,0], ["A", "B", "C"], true)

  @test to_table(list, false) == table
  @test from_table(table, Int64, false) == list

  @test to_table(list) == [ "A" "A" 0
                            "A" "B" 10
                            "A" "C" 20
                            "B" "B" 0
                            "B" "C" 30
                            "C" "C" 0]
  @test to_table(list) == to_table(list_diag)
  @test to_table(list, false) == to_table(list_diag, false)
end

print("""

./ Int & Float64
----------------
""")

let  list = PairwiseListMatrix([.5, .4, .3])

  result =  list ./ 4
  @test isa(result, PairwiseListMatrix{Float64,false})
  @test sum(result) == 0.6

  result =  list ./ 4.
  @test isa(result, PairwiseListMatrix{Float64,false})
  @test sum(result) == 0.6
end

let  list = PairwiseListMatrix([.5, .4, .3], true)

  result =  list ./ 4
  @test isa(result, PairwiseListMatrix{Float64,true})
  @test_approx_eq sum(result) (0.5 + 0.4*2.0 + 0.3)/4

  result =  list ./ 4.
  @test isa(result, PairwiseListMatrix{Float64,true})
  @test_approx_eq sum(result) (0.5 + 0.4*2.0 + 0.3)/4.0

end

print("""

Z score
-------
""")

let list = PairwiseListMatrix{Int, false}[ PairwiseListMatrix([1, 1, 1]), PairwiseListMatrix([3, 3, 3]) ],
  mat = PairwiseListMatrix([2, 2, 2])

  @test std(list) == PairwiseListMatrix([1., 1., 1.])
  @test mean(list) == mat
  @test PairwiseListMatrices.zscore(list, mat) == PairwiseListMatrix([0., 0., 0.])

  @test_throws ErrorException PairwiseListMatrices.zscore(list, PairwiseListMatrix([2, 2, 2], true))
  @test_throws ErrorException PairwiseListMatrices.zscore(list, PairwiseListMatrix([2, 2, 2, 2, 2, 2]))
end

let list = PairwiseListMatrix{Int, true}[ PairwiseListMatrix([1, 1, 1], true), PairwiseListMatrix([3, 3, 3], true) ],
  mat = PairwiseListMatrix([2, 2, 2], true)

  @test std(list) == PairwiseListMatrix([1., 1., 1.], true)
  @test mean(list) == mat
  @test PairwiseListMatrices.zscore(list, mat) == PairwiseListMatrix([0., 0., 0.], true)

  @test_throws ErrorException PairwiseListMatrices.zscore(list, PairwiseListMatrix([2, 2, 2]))
  @test_throws ErrorException PairwiseListMatrices.zscore(list, PairwiseListMatrix([2, 2, 2, 2, 2, 2], true))
end

print("""

macros
------
""")

let PLM = PairwiseListMatrix([1,2,3], false)
  @iteratelist PLM Base.Test.@test list[k] == k
end

values   = [1,2,3,4,5,6]
PLMtrue  = PairwiseListMatrix(values, true)
PLMfalse = PairwiseListMatrix(values, false)
full_t   = full(PLMtrue)
full_f   = full(PLMfalse)

@iteratelist PLMtrue  Base.Test.@test list[k] == Main.values[k]
@iteratelist PLMfalse Base.Test.@test list[k] == Main.values[k]

@iteratediag PLMtrue  Base.Test.@test false
@iteratediag PLMfalse Base.Test.@test diag[k] == 0

@iterateupper PLMtrue  true  list[k] = Main.values[k]
@iterateupper PLMfalse false list[k] = Main.values[k]

@iterateupper PLMtrue  true  list[k] = Main.full_t[i,j]
@iterateupper PLMtrue  false list[k] = Main.full_t[i,j]

@iterateupper PLMtrue  true  list[k] = Main.full_f[i,j]
@iterateupper PLMfalse false list[k] = Main.full_f[i,j]

print("""

write...
--------
""")

let values   = [1,2,3,4,5,6],
    PLMtrue  = PairwiseListMatrix(values, true),
    PLMfalse = PairwiseListMatrix(values, false),
    filename = "test.temp"

  try
    writecsv(filename, PLMtrue, true)
    @test all( readcsv(filename) .== to_table(PLMtrue, true) )
  finally
    rm(filename)
  end

  try
    writecsv(filename, PLMtrue, false)
    @test all( readcsv(filename) .== to_table(PLMtrue, false) )
  finally
    rm(filename)
  end

  try
    writecsv(filename, PLMfalse, true)
    @test all( readcsv(filename) .== to_table(PLMfalse, true) )
  finally
    rm(filename)
  end

  try
    writecsv(filename, PLMfalse, false)
    @test all( readcsv(filename) .== to_table(PLMfalse, false) )
  finally
    rm(filename)
  end

end

print("""

join
----
""")

let left = PairwiseListMatrix(collect(1.:15.), Char['a','b','c','d','e'], true),
    right = PairwiseListMatrix(collect(1.:10.), Char['a','e','i','o','u'], false)

  a, b = join(left, right)
  @test (a, b) == join(left, right, kind=:inner)
  @test labels(a) == labels(b)
  @test labels(a) == ['a', 'e']
  @test size(a) == (2,2)
  @test size(b) == (2,2)

  a, b = join(left, right, kind=:left)
  @test labels(a) == labels(b)
  @test labels(a) == labels(left)
  @test a == left
  @test b != right

  a, b = join(left, right, kind=:right)
  @test labels(a) == labels(b)
  @test labels(b) == labels(right)
  @test a != left
  @test b == right

  a, b = join(left, right, kind=:outer)
  @test labels(a) == labels(b)
  @test labels(a) == ['a', 'b', 'c', 'd', 'e', 'i', 'o', 'u']
  @test size(a) == (8,8)
  @test size(b) == (8,8)
  @test a != left
  @test b != right

end

# #### DataFrames #### #

if dataframes_installed
  print("""

  # DataFrames
  # ==========
  """)

  let values   = [1,2,3,4,5,6],
    PLMtrue  = PairwiseListMatrix(values, true),
    PLMfalse = PairwiseListMatrix(values, false)

    df = to_dataframe(PLMtrue)
    @test PLMtrue == from_dataframe(df, true)

    df = to_dataframe(PLMtrue, false)
    @test triu(PLMtrue,1) == triu(from_dataframe(df, false),1)

    df = to_dataframe(PLMfalse)
    @test PLMfalse == from_dataframe(df, true)

    df = to_dataframe(PLMfalse, false)
    @test PLMfalse == from_dataframe(df, false)

  end

end


text:
module SynapticFilter

using CUDA: AnyCuMatrix
using Distances: msd
using Distributions: Poisson
using LinearAlgebra: dot, diagind
using Statistics: mean

include("models.jl")
export
    update!,
    State,
    NeuronObs,
    PoissonExpModel, BlockPoissonExpModel,
    ExpGainModel, AdaptiveExpGainModel,
    OUModel

include("filters.jl")
export
    FilterState,
    BlockSF, DiagSF, FullSF,
    GradientRule

include("simulation.jl")
export
    Simulation, SimulationState

include("performance.jl")
export
    ErrorMeasure, ComputeMSE

end#module


text:
#####
##### Functions to deal with a cross-section interpolator
#####

using ..Architectures: GPU

"""
    $(FUNCTIONNAME)(hitran::HitranTable,
                    broadening::AbstractBroadeningFunction,
                    wave_grid::AbstractRange{<:Real},
                    p_grid::AbstractRange{<:Real},
                    t_grid::AbstractRange{<:Real};
                    wavelength_flag::Bool=false,
                    wing_cutoff::Real=40,
                    vmr::Real=0,
                    CEF::AbstractComplexErrorFunction=HumlicekWeidemann32SDErrorFunction(),
                    architecture::AbstractArchitecture=default_architecture )

Using a HitranModel, create an InterpolationModel by interpolating the lineshape function
at the given pressure and temperature grids

"""
function make_interpolation_model(
                                  # Required
                                  hitran::HitranTable,
                                  broadening::AbstractBroadeningFunction,
                                  wave_grid::AbstractRange{<:Real},
                                  p_grid::AbstractRange{<:Real},
                                  t_grid::AbstractRange{<:Real};
                                  # Optionals
                                  wavelength_flag::Bool=false,
                                  wing_cutoff::Integer=40,
                                  vmr::Real=0,
                                  CEF::AbstractComplexErrorFunction=HumlicekWeidemann32SDErrorFunction(),
                                  architecture::AbstractArchitecture=default_architecture     # Computer `Architecture` on which `Model` is run
                                )

    # Warn user if using incompatible/untested CEF
    if architecture isa GPU && !(CEF isa HumlicekWeidemann32SDErrorFunction)
        @warn "Cross-section calculations on GPU may or may not work with this CEF (use HumlicekWeidemann32SDErrorFunction if you encounter issues)"
    end

    # Convert from wavelength to wavenumber if necessary
    ν_grid = wavelength_flag ? reverse(nm_per_m ./ wave_grid) : wave_grid

    # Empty matrix to store the calculated cross-sections
    cs_matrix = zeros(length(p_grid), length(t_grid), length(ν_grid));

    # Calculate all the cross-sections at the pressure and temperature grids
    @showprogress 1 "Computing Cross Sections for Interpolation..." for i in 1:length(p_grid)
        for j in 1:length(t_grid)
            # make_hitran_model(hitran_data, Voigt(), wing_cutoff = 40, CEF=HumlicekWeidemann32SDErrorFunction(), architecture=CPU())
            model = make_hitran_model(hitran, broadening, wing_cutoff=wing_cutoff, CEF=CEF, architecture=architecture)
            cs_matrix[i,j,:] = Array(compute_absorption_cross_section(model, collect(ν_grid), p_grid[i], t_grid[j], wavelength_flag=wavelength_flag))
        end
    end

    # Perform the interpolation
    itp = interpolate(cs_matrix, BSpline(Cubic(Line(OnGrid()))))

    # Get the molecule and isotope numbers from the HitranTable
    mol = hitran.mol[1]
    iso = all(x->x==hitran.iso[1], hitran.iso) ? hitran.iso[1] : -1

    # Return the interpolation model with all the proper parameters
    return InterpolationModel(itp, mol, iso, broadening, ν_grid, p_grid, t_grid, wing_cutoff, vmr, CEF, architecture)
end

""" Convenience function to save an InterpolationModel at a specified filepath (Using JLD2) """
function save_interpolation_model(itp_model::InterpolationModel, filepath::String)
    @save filepath itp_model
end

""" Convenience function to load an InterpolationModel from a specified filepath (Using JLD2) """
function load_interpolation_model(filepath::String)
    itp_model::InterpolationModel = @load filepath itp_model
    return itp_model
end


text:
const FHIR_TYPE_MAP = Dict{String, String}(
    "ResourceList" => "ResourceList",
    "base64Binary" => "String",
    "boolean"      => "Bool",
    "canonical"    => "String",
    "code"         => "String",
    "date"         => "Date",
    "dateTime"     => "FHIRDateTime",
    "decimal"      => "FHIRNumber",
    "id"           => "String",
    "instant"      => "ZonedDateTime",
    "integer"      => "SafeInt32",
    "markdown"     => "String",
    "oid"          => "String",
    "positiveInt"  => "SafeInt32",
    "string"       => "String",
    "time"         => "Time",
    "unsignedInt"  => "SafeInt32",
    "uri"          => "String",
    "url"          => "String",
    "uuid"         => "String",
    "xhtml"        => "String",
)

@inline function has_julia_type(fhirtype::AbstractString)
    return (fhirtype in collect(keys(FHIR_TYPE_MAP)))::Bool
end

@inline function get_juliatype(fhirtype::AbstractString)
    if fhirtype in collect(keys(FHIR_TYPE_MAP))
        return FHIR_TYPE_MAP[fhirtype]::String
    end
    throw(KeyError("\"$(fhirtype)\" not found"))
end

@inline function get_juliatype_maybe(fhirtype::AbstractString)
    if fhirtype in collect(keys(FHIR_TYPE_MAP))
        return FHIR_TYPE_MAP[fhirtype]::String
    end
    return fhirtype
end


text:
export Vandermonde

using LinearAlgebra: Adjoint, Transpose

"""
    Vandermonde(c::AbstractVector)

Create a "lazy" `n × n`
[`Vandermonde` matrix](http://en.wikipedia.org/wiki/Vandermonde_matrix)
where ``A_{ij} = c_i^{j-1}``,
requiring only `O(n)` storage for the vector `c`.

The `transpose` and `adjoint` operations are also lazy.

```jldoctest van
julia> a = 1:5; A = Vandermonde(a)
5×5 Vandermonde{Int64}:
 1  1   1    1    1
 1  2   4    8   16
 1  3   9   27   81
 1  4  16   64  256
 1  5  25  125  625

julia> A'
5×5 adjoint(::Vandermonde{Int64}) with eltype Int64:
 1   1   1    1    1
 1   2   3    4    5
 1   4   9   16   25
 1   8  27   64  125
 1  16  81  256  625
```

The backslash operator `\\` is overloaded to solve Vandermonde and adjoint
Vandermonde systems in `O(n^2)` time using the algorithm of
Björck & Pereyra (1970), https://doi.org/10.2307/2004623.

```jldoctest van
julia> A \\ a
5-element Vector{Float64}:
 0.0
 1.0
 0.0
 0.0
 0.0

julia> A' \\ A[2,:]
5-element Vector{Float64}:
 0.0
 1.0
 0.0
 0.0
 0.0
```
"""
struct Vandermonde{T} <: AbstractMatrix{T}
    c :: Vector{T}
end

function Vandermonde(c::AbstractVector{T}) where T
    return Vandermonde{T}(convert(Vector{T}, c))
end

@inline Base.@propagate_inbounds function getindex(
    V::Vandermonde,
    i::Int,
    j::Int,
)
    @boundscheck checkbounds(V, i, j)
    return (@inbounds V.c[i])^(j-1)
end

size(V::Vandermonde) = (length(V.c), length(V.c))

# solvers

function \(V::Adjoint{T1,<:Vandermonde{T1}}, y::AbstractVecOrMat{T2}) where {T1, T2}
    T = vandtype(T1,T2)
    x = Array{T}(undef, size(y))
    copyto!(x, y)
    pvand!(adjoint(V.parent.c), x)
    return x
end

function \(V::Transpose{T1,<:Vandermonde{T1}}, y::AbstractVecOrMat{T2}) where {T1, T2}
    T = vandtype(T1,T2)
    x = Array{T}(undef, size(y))
    copyto!(x, y)
    pvand!(V.parent.c, x)
    return x
end

function \(V::Vandermonde{T1}, y::AbstractVecOrMat{T2}) where {T1, T2}
    T = vandtype(T1,T2)
    x = Array{T}(undef, size(y))
    copyto!(x, y)
    dvand!(V.c, x)
    return x
end

"""
    vandtype(T1::Type, T2::Type)
Determine the return type of `Vandermonde{T1} \\ Vector{T2}`.
"""
function vandtype(T1::Type, T2::Type)
    T = promote_type(T1, T2)
    S = typeof(oneunit(T) / oneunit(T1))
    return S
end

"""
    pvand!(a, b) -> b

Solve system ``A^T*x = b`` in-place.

``A^T`` is transpose of Vandermonde matrix ``A_{ij} = a_i^{j-1}``.

Algorithm by Bjorck & Pereyra,
Mathematics of Computation, Vol. 24, No. 112 (1970), pp. 893-903,
https://doi.org/10.2307/2004623
"""
function pvand!(alpha, B)
    n = length(alpha)
    n == size(B,1) ||
        throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
    nrhs = size(B,2)
    @inbounds begin
        for j=1:nrhs
            for k=1:n-1
                for i=n:-1:k+1
                    B[i,j] -= alpha[k] * B[i-1,j]
                end
            end
            for k=n-1:-1:1
                for i=k+1:n
                    B[i,j] = B[i,j] / (alpha[i] - alpha[i-k])
                end
                for i=k:n-1
                    B[i,j] -= B[i+1,j]
                end
            end
        end
    end
    return B
end

"""
    dvand!(a, b) -> b

Solve system ``A*x = b`` in-place.

`A` is Vandermonde matrix ``A_{ij} = a_i^{j-1}``.

Algorithm by Bjorck & Pereyra,
Mathematics of Computation, Vol. 24, No. 112 (1970), pp. 893-903,
https://doi.org/10.2307/2004623
"""
function dvand!(alpha, B)
    n = length(alpha)
    n == size(B,1) ||
        throw(DimensionMismatch("matrix has dimensions ($n,$n) but right hand side has $(size(B,1)) rows"))
    nrhs = size(B,2)
    @inbounds begin
        for j=1:nrhs
            for k=1:n-1
                for i=n:-1:k+1
                    B[i,j] = (B[i,j] - B[i-1,j]) / (alpha[i] - alpha[i-k])
                end
            end
            for k=n-1:-1:1
                for i=k:n-1
                    B[i,j] -= alpha[k] * B[i+1,j]
                end
            end
        end
    end
    return B
end


text:
#======
A. HAC
=======#

function lrvar(
    k::HAC,
    m::AbstractMatrix;
    prewhite::Bool = false,
    demean::Bool=true,
    scale::Real=1,
)
    mm = demean ? m .- mean(m, dims = 1) : m
    scale *= inv(size(m,1))
    return _lrvar(k, mm, prewhite, scale)
end

function _lrvar(k::HAC, mm, prewhite::Bool, scale)
    n, p = size(mm)
    m, D = prewhiter(mm, Val{prewhite})
    bw = optimalbandwidth(k, m; prewhite=prewhite)
    V = triu!(m'*m)
    @inbounds for j in covindices(k, n)
        κⱼ = kernel(k, j/bw)
        LinearAlgebra.axpy!(κⱼ, Γ(m, j), V)
    end
    LinearAlgebra.copytri!(V, 'U')
    dewhiter!(V, m, D, Val{prewhite})
    return Symmetric(V.*convert(eltype(V), scale))
end

function lrvarmatrix(
    k::HAC,
    m::AbstractMatrix,
    factorization=Cholesky;
    prewhite::Bool=false,
    demean::Bool=true,
    scale::Real=1,
)
    mm = demean ? m .- mean(m, dims = 1) : m
    scale *= inv(size(m,1))
    return _lrvarmatrix(k, mm, prewhite, scale, factorization)
end

function _lrvarmatrix(
    k::HAC,
    mm::AbstractMatrix,
    prewhite::Bool,
    scale::Real,
    ::Type{Cholesky},
)
    V = _lrvar(k, mm, prewhite, scale)
    CovarianceMatrix(cholesky(V, check=true), k, V)
end

function _lrvarmatrix(
    k::HAC,
    mm::AbstractMatrix,
    prewhite::Bool,
    scale::Real,
    ::Type{SVD},
)
    V = _lrvar(k, mm, prewhite, scale)
    CovarianceMatrix(svd(V.data), k, V)
end

#=======
B. VARHC
=======#
"""

"""
function lrvar(
    k::VARHAC,
    m::AbstractMatrix,
    demean::Bool=true,
    scale::Real=1
)
    mm = demean ? m .- mean(m, dims = 1) : m
    return _lrvar(k, mm, scale)
end

function _lrvar(k::VARHAC, m, scale)
    maxlag, lagstrategy, selectionstrategy = k.maxlag, k.lagstrategy, k.selectionstrategy
    strategy = selectionstrategy == :aic ? 1 : (selectionstrategy == :bic ? 2 : 3)
    return Symmetric(varhac(m,maxlag,lagstrategy,strategy).*scale)
end

#=======
C. HC
=======#
function lrvar(
    k::HC,
    m::AbstractMatrix;
    demean::Bool=true,
    scale::Real=1
)
    mm = demean ? m .- mean(m, dims=1) : m
    return _lrvar(k, mm, scale)
end

function _lrvar(k::HC, m::AbstractMatrix{T}, scale) where T
    V = m'*m
    F = promote_type(T, eltype(V))
    return Symmetric(V.*convert(F, inv(size(m,1))*scale))
end

function lrvarmatrix(
    k::HC,
    m::AbstractMatrix,
    factorization=Cholesky;
    demean::Bool=true,
    scale::Real=1
)
    mm = demean ? m .- mean(m, dims=1) : m
    return _lrvarmatrix(k, mm, scale, factorization)
end

function _lrvarmatrix(k::HC, m::AbstractMatrix, scale, ::Type{Cholesky})
    V = _lrvar(k, m, scale)
    return CovarianceMatrix(cholesky(V, check=true), k, V)
end

function _lrvarmatrix(k::HC, m::AbstractMatrix, scale, ::Type{SVD})
    V = _lrvar(k, m, scale)
    return CovarianceMatrix(svd(V.data), k, V)
end
#======
D. CRHC
=======#
function lrvar(
    k::CRHC,
    m::AbstractMatrix;
    demean::Bool=true,
    scale::Real=1,
)
    mm = demean ? m .- mean(m, dims=1) : m
    scale *= inv(size(m,1))*scale
    return _lrvar(k, mm, scale)
end

function _lrvar(k::CRHC, m::AbstractMatrix{T}, scale) where T
    cache = installcache(k, m)
    Shat = clusterize!(cache)
    F = promote_type(T, eltype(Shat))
    return Symmetric(Shat.*convert(F, scale))
end

function lrvarmatrix(
    k::CRHC,
    m::AbstractMatrix,
    factorization=Cholesky;
    demean::Bool=true,
    scale::Real=1,
)
    mm = demean ? m .- mean(m, dims=1) : m
    return _lrvarmatrix(k, mm, scale, factorization)
end

function _lrvarmatric(k::CRHC, m::AbstractMatrix, scale, ::Type{Cholesky})
    V = _lrvar(k, m, scale)
    CovarianceMatric(cholesky(V, check=true), k, V)
end

function _lrvarmatric(k::CRHC, m::AbstractMatrix, scale, ::Type{SVD})
    V = _lrvar(k, m, scale)
    CovarianceMatric(SVD(V.data), k, V)
end
#=========
Finalizers
==========#
# factorizer(::Type{SVD}, V) = svd(V.data)
# factorizer(::Type{Cholesky}, V) = cholesky(V, check = false)

# finalize(k, V, M, F) = finalize(k, V, M, F, 1)

# function finalize(k, V, ::Type{M}, F, scale) where M<:Matrix
#     return Symmetric(V.*scale)
# end

# function finalize(k, V, ::Type{M}, F, scale) where M<:CovarianceMatrix
#     V .= V.*scale
#     CovarianceMatrix(factorizer(F, V), k, Symmetric(V))
# end


text:
module CEnum

abstract Cenum{T}
Base.(:|){T<:Cenum}(a::T, b::T) = T(Int(a) | Int(b))
Base.(:&){T<:Cenum}(a::T, b::T) = T(Int(a) & Int(b))
# typemin and typemax won't change for an enum, so we might as well inline them per type
function Base.typemax{T<:Cenum}(::Type{T})
    last(enum_values(T))
end
function Base.typemin{T<:Cenum}(::Type{T})
    first(enum_values(T))
end
Base.convert{T<:Integer}(::Type{T}, x::Cenum) = convert(T, Base.box(Int32, x))
Base.write(io::IO, x::Cenum) = write(io, Int32(x))
Base.read{T<:Cenum}(io::IO, ::Type{T}) = T(read(io, Int32))

enum_values{T<:Cenum}(::T) = enum_values(T)
enum_names{T<:Cenum}(::T) = enum_names(T)

function is_member{T<:Cenum}(::Type{T}, x::Integer)
    is_member(T, enum_values(T), x)
end
@inline is_member{T<:Cenum}(::Type{T}, r::UnitRange, x::Integer) = x in r
@inline function is_member{T<:Cenum}(::Type{T}, values::Tuple, x::Integer)
    lo, hi = typemin(T), typemax(T)
    x<lo || x>hi && return false
    for val in values
        val == x && return true
        val > x && return false # is sorted
    end
    return false
end

function enum_name{T<:Cenum}(x::T)
    index = findfirst(enum_values(T), Int(x))
    if index != 0
        return enum_names(T)[index]
    end
    error("Invalid enum: $x, name not found")
end
function Base.show(io::IO, x::Cenum)
    print(io, enum_name(x), "($(Int(x)))")
end

function islinear(array)
    isempty(array) && return false # false, really? it's kinda undefined?
    lastval = first(array)
    for val in rest(array, 2)
        val-lastval == 1 || return false
    end
    return true
end

macro cenum(name, args...)
    if !isa(name, Symbol)
        error("Name must be symbol or Name{Type}. Found: $name")
    end
    lastval = -1
    name_values = map([args...]) do arg
        if isa(arg, Symbol)
            lastval += 1
            val = lastval
            sym = arg
        elseif arg.head == :(=) || arg.head == :kw
            sym,val = arg.args
        else
            error("Expression of type $arg not supported. Try only symbol or name = value")
        end
        (sym, val)
    end
    sort!(name_values, by=last) # sort for values
    values = map(last, name_values)

    if islinear(values) # optimize for linear values
        values = :($(first(values)):$(last(values)))
    else
        values = :(tuple($(values...)))
    end
    value_block = Expr(:block)
    typename = esc(name)
    for (ename, value) in name_values
        push!(value_block.args, :(const $(esc(ename)) = $typename($value)))
    end

    expr = quote
        bitstype 32 $typename <: CEnum.Cenum{UInt32}
        function Base.convert(::Type{$typename}, x::Integer)
            is_member($typename, x) || Base.enum_argument_error($(Expr(:quote, name)), x)
            Base.box($typename, convert(Int32, x))
        end
        CEnum.enum_names(::Type{$typename}) = tuple($(map(x-> Expr(:quote, first(x)), name_values)...))
        CEnum.enum_values(::Type{$typename}) = $values
        $value_block
    end
    expr
end
export @cenum

end # module


text:
"""
# Greedy algorithms

- Minimum spanning trees
    - A greedy approach
    - Trees
    - The cut property
    - Kruskal's algorithm
    - A data structure for disjoint sets
    - Union by rank
    - Path compression
    - A randomized algorithm for minimum cut
    - Prim's algorithm
- Huffman encoding
    - Entropy
- Horn formulas
- Set cover
"""
module Chapter5

end # module