text:
using PlanktonIndividuals, Serialization

grid = RectilinearGrid(size = (16, 16, 1), x = (0,32), y = (0,32), z = (0,-32))

model = PlanktonModel(CPU(), grid)

function tot_mass(nut, g)
    mass = zeros(g.Nx, g.Ny, g.Nz)
    for i in 1:g.Nx
        for j in 1:g.Ny
            for k in 1:g.Nz
                mass[i,j,k] = nut[i+g.Hx, j+g.Hy, k+g.Hz] * PlanktonIndividuals.Grids.volume(i+g.Hx, j+g.Hy, k+g.Hz, g)
            end
        end
    end
    return sum(mass)
end

TP = tot_mass(model.nutrients.PO4.data, grid) +
     tot_mass(model.nutrients.DOP.data, grid) +
     tot_mass(model.nutrients.POP.data, grid)
TP = TP + sum(model.individuals.phytos.sp1.data.Pq .+
              model.individuals.phytos.sp1.data.Bm .* model.individuals.phytos.sp1.p.R_PC)

uvel = zeros(16,16,1,11)
vvel = zeros(16,16,1,11)
wvel = zeros(16,16,2,11)

for i in 1:11
    uvel[:,:,:,i] .= randn(16,16,1) .* 1e-4
    vvel[:,:,:,i] .= randn(16,16,1) .* 1e-4
    wvel[:,:,:,i] .= randn(16,16,2) .* 1e-4
end

sim = PlanktonSimulation(model, ΔT = 60, iterations = 10, vels=(u=uvel, v=vvel, w=wvel))

update!(sim)

TPt = tot_mass(model.nutrients.PO4.data, grid) +
      tot_mass(model.nutrients.DOP.data, grid) +
      tot_mass(model.nutrients.POP.data, grid)
TPt = TPt + sum(model.individuals.phytos.sp1.data.Pq .+
                model.individuals.phytos.sp1.data.Bm .* model.individuals.phytos.sp1.p.R_PC)

@testset "PlanktonIndividuals 2D tests:" begin
    @test isapprox(TP,TPt; atol=1e1)
end


text:
using StatsBase
using StatsModels
using GLM
using SparseArrays
using DataFrames
#using RDatasets: dataset
using Test

using RobustModels

# Import data
#data = dataset("robustbase", "Animals2")

Animal = ["Lesser short-tailed shrew", "Little brown bat", "Big brown bat", "Mouse", "Musk shrew", "Star-nosed mole", "E. American mole", "Ground squirrel", "Tree shrew", "Golden hamster", "Mole", "Galago", "Rat", "Chinchilla", "Owl monkey", "Desert hedgehog", "Rock hyrax-a", "European hedgehog", "Tenrec", "Artic ground squirrel", "African giant pouched rat", "Guinea pig", "Mountain beaver", "Slow loris", "Genet", "Phalanger", "N.A. opossum", "Tree hyrax", "Rabbit", "Echidna", "Cat", "Artic fox", "Water opossum", "Nine-banded armadillo", "Rock hyrax-b", "Yellow-bellied marmot", "Verbet", "Red fox", "Raccoon", "Rhesus monkey", "Potar monkey", "Baboon", "Roe deer", "Goat", "Kangaroo", "Grey wolf", "Chimpanzee", "Sheep", "Giant armadillo", "Human", "Grey seal", "Jaguar", "Brazilian tapir", "Donkey", "Pig", "Gorilla", "Okapi", "Cow", "Horse", "Giraffe", "Asian elephant", "African elephant", "Triceratops", "Dipliodocus", "Brachiosaurus"]
Brain = [0.14, 0.25, 0.3, 0.4, 0.33, 1.0, 1.2, 4.0, 2.5, 1.0, 3.0, 5.0, 1.9, 6.4, 15.5, 2.4, 12.3, 3.5, 2.6, 5.7, 6.6, 5.5, 8.1, 12.5, 17.5, 11.4, 6.3, 12.3, 12.1, 25.0, 25.6, 44.5, 3.9, 10.8, 21.0, 17.0, 58.0, 50.4, 39.2, 179.0, 115.0, 179.5, 98.2, 115.0, 56.0, 119.5, 440.0, 175.0, 81.0, 1320.0, 325.0, 157.0, 169.0, 419.0, 180.0, 406.0, 490.0, 423.0, 655.0, 680.0, 4603.0, 5712.0, 70.0, 50.0, 154.5]
Body = [0.005, 0.01, 0.023, 0.023, 0.048, 0.06, 0.075, 0.101, 0.104, 0.12, 0.122, 0.2, 0.28, 0.425, 0.48, 0.55, 0.75, 0.785, 0.9, 0.92, 1.0, 1.04, 1.35, 1.4, 1.41, 1.62, 1.7, 2.0, 2.5, 3.0, 3.3, 3.385, 3.5, 3.5, 3.6, 4.05, 4.19, 4.235, 4.288, 6.8, 10.0, 10.55, 14.83, 27.66, 35.0, 36.33, 52.16, 55.5, 60.0, 62.0, 85.0, 100.0, 160.0, 187.1, 192.0, 207.0, 250.0, 465.0, 521.0, 529.0, 2547.0, 6654.0, 9400.0, 11700.0, 87000.0]
data = DataFrame(; Animal=Animal, Brain=Brain, Body=Body)

data.logBrain = log.(data.Brain)
data.logBody = log.(data.Body)
form = @formula(logBrain ~ 1 + logBody)
X = hcat(ones(size(data, 1)), data.logBody)
sX = SparseMatrixCSC(X)
y = data.logBrain

# Include tests
include("estimators.jl")
include("interface.jl")
include("linearfit.jl")
include("mquantile.jl")
include("robustridge.jl")
include("qreg.jl")
include("univariate.jl")

ext: jl


text:
push!(LOAD_PATH, "./gridmodule/")
push!(LOAD_PATH, "./femmodule/")

using BaseFEM
using BaseGrid
using LinearAlgebra

# 1D ring
#m=grid1D("./mesh1d/mesh2")
#tau = 1e-3
#eps = 0.1

# 3D ball
#m=grid3D("./mesh3d/ball")
#tau = 1e-3
#eps = 0.1

# greenland example
#m=grid2D("./mesh2d/greenland")
#tau = 1e-1
#eps = 3.00

# sphere example
m=grid2D("./mesh2d/sphere_hq05")
tau = 1e-6
eps  = 0.01

top  = m.topology
ndof = ndof_P1(top)
geo  = restrain_geometry(m.geometry,ndof)
A,M=assembly_P1(geo,top)

# initial conditions
u      = 2.0*rand(ndof).-1.0
time = 0.0
S = factorize([M tau*A;-A M])
mprint = 0

for i=1:1000
	global time,u,mprint
	println("iteration $(i) at time $(time)")

	dWdu = (4*u.^3-4*u)/eps

	rhs1 = M*u
	rhs2 = M*dWdu

	sol = S\[rhs1;rhs2]
	u   = sol[1:ndof]
	time = time + tau

	# output to ParaView vtk file
	if mod(i,10) == 0
		mprint = mprint + 1
		outvtk_scalar("output/example_cahnhilliard"*string(mprint)*".vtk",geo,top,"sol",u)
	end
end
# alternative meshes
# m=grid2D("./mesh2d/test")
# m=grid2D("./mesh2d/sphere")

max_stars_repo_name: dpeschka/jPDE

max_issues_repo_name: dpeschka/jPDE

max_forks_repo_name: dpeschka/jPDE


text:
using ClusterTrees
using StaticArrays
using Test

const P = SVector{3,Float64}

using DelimitedFiles
Q = readdlm(joinpath(@__DIR__,"points.dlm"), Float64)
points = [SVector{3,Float64}(Q[i,:]) for i in axes(Q,1)]

root_center = P(0.5,0.5,0.5)
root_size = 0.5
root_node = ClusterTrees.LevelledTrees.HNode(ClusterTrees.PointerBasedTrees.Node(ClusterTrees.LevelledTrees.Data(0,Int[]), 0, 0, 0, 0), 0)
tree = ClusterTrees.LevelledTrees.LevelledTree([root_node], 1, root_center, root_size, Int[1])

smallest_box_size = 0.1
root_sector = 2
root_sfc_state = 2
root_depth = 1
for i in 1:length(points)
    router = ClusterTrees.LevelledTrees.Router(smallest_box_size, points[i])
    root_state = root(tree), root_sector, root_center, root_size, root_sfc_state, root_depth
    ClusterTrees.update!(tree, root_state, i, router)  do tree, node, data
        push!(ClusterTrees.data(tree,node).values, data)
    end
end

# ClusterTrees.print_tree(tree)

num_bins = 3
bins = [P[] for i in 1:num_bins]

num_printed = 0
num_points = 0
num_nodes = length(tree.nodes)
for (i,node) in enumerate(ClusterTrees.DepthFirstIterator(tree, root(tree)))
    # println(node, ": ", tree.nodes[node])
    b = div((i-1)*num_bins, num_nodes) + 1
    append!(bins[b], points[ClusterTrees.data(tree,node).values])
    global num_printed += 1
    global num_points += length(data(tree,node).values)
end

ordered_points = reduce(append!, bins)
ordered_points = [p[i] for p in ordered_points, i = 1:3]
@test num_printed == length(tree.nodes)
@test num_points == length(points)

num_points = 0
ordered_points = P[]
for b in ClusterTrees.LevelledTrees.LevelIterator(tree, length(tree.levels))
    global num_points += length(ClusterTrees.data(tree,b).values)
    append!(ordered_points, points[ClusterTrees.data(tree,b).values])
end
@test num_points == length(points)

@test length(tree.levels) == 4
@test ClusterTrees.LevelledTrees.numlevels(tree) == 4
lvl1 = collect(ClusterTrees.LevelledTrees.LevelIterator(tree,1))
@test lvl1 == [1]
lvl2 = collect(ClusterTrees.LevelledTrees.LevelIterator(tree,2))
@test lvl2 == [10,13,2,5,23]
lvl3 = collect(ClusterTrees.LevelledTrees.LevelIterator(tree,3))
@test lvl3 == [19,11,21,16,14,3,8,6,24]
lvl4 = collect(ClusterTrees.LevelledTrees.LevelIterator(tree,4))
@test lvl4 == [20,18,12,22,17,15,4,9,7,25]


text:
using .MCMCChains
export pointwise_log_likelihoods

"""
    pointwise_log_likelihoods(ll_fun::Function, chains::Chains, data)

Compute the pointwise log likelihoods.

# Arguments
  - $LIKELIHOOD_FUNCTION_ARG
  - `chain::Chains`: A chain object from MCMCChains.
  - `data`: An array of data points used to estimate the parameters of the model.

# Returns
  - `Array`: a three dimensional array of pointwise log-likelihoods. Dimensions are ordered
    as `[data, step, chain]`.
"""
function pointwise_log_likelihoods(
    ll_fun::Function, chain::Chains, data::AbstractVector; kwargs...
)
    samples = Chains(chain, :parameters).value
    return pointwise_log_likelihoods(ll_fun, samples, data; kwargs...)
end

"""
    function psis_loo(
        ll_fun::Function,
        chain::Chains,
        data,
        args...;
        kwargs...
    ) -> PsisLoo

Use Pareto-Smoothed Importance Sampling to calculate the leave-one-out cross validation
score from an MCMCChains object.

# Arguments

  - $LIKELIHOOD_FUNCTION_ARG
  - `chain::Chain`: A chain object from MCMCChains.
  - `data`: A vector of data points used to estimate the parameters of the model.
  - $ARGS [`psis_loo`](@ref).
  - $KWARGS [`psis_loo`](@ref).

See also: [`psis_loo`](@ref), [`loo`](@ref), [`PsisLoo`](@ref).
"""
function psis_loo(ll_fun::Function, chain::Chains, data::AbstractVector, args...; kwargs...)
    pointwise_log_likes = pointwise_log_likelihoods(ll_fun, chain, data)
    return psis_loo(pointwise_log_likes, args...; kwargs...)
end

"""
    psis(ll_fun::Function, chain::Chains, data[, kwargs...]; args...) -> Psis

Implements Pareto-smoothed importance sampling (PSIS) based on MCMCChain object.

# Arguments

  - $LIKELIHOOD_FUNCTION_ARG
  - `chain::Chain`: A chain object from MCMCChains.
  - `data`: A vector of data points used to estimate the parameters of the model.
  - $ARGS [`psis`](@ref).
  - $KWARGS [`psis`](@ref).

See also: [`psis`](@ref), [`psis_loo`](@ref), [`PsisLoo`](@ref).
"""
function psis(ll_fun::Function, chain::Chains, data::AbstractVector, args...; kwargs...)
    pointwise_log_likes = pointwise_log_likelihoods(ll_fun, chain, data)
    return psis(-pointwise_log_likes, args...; kwargs...)
end


text:
#mutable struct DisjointSetNode
#    rank::Int
#    p::DisjointSetNode
#    DisjointSetNode() = (obj = new(0); obj.p = obj;)
#end

function findset(x::DisjointSetNode)
    if (x != x.p)
        x.p = findset(x.p)
    end
    return x.p
end

function my_union!(x::DisjointSetNode, y::DisjointSetNode)
    link(findset(x), findset(y))
end

function link(x::DisjointSetNode, y::DisjointSetNode)
    if (x.rank > y.rank)
        y.p = x
    else
        x.p = y
        if (y.rank == x.rank)
            y.rank = y.rank + 1
        end
    end
end

function hammingdistance(s1::String, s2::String)
    dist = 0 # accumulated "distance" between strings
    for i = 1:length(s1)
        if (Int(s1[i]) - Int(s2[i]) != 0) # if letters at element i dont match
            dist = dist + 1 # increase distance strings
        end
    end
    return dist
end

function findclusters(E::Vector{Tuple{Int, Int, Int}}, n::Int, k::Int)
    A = []
    sort!(E)
    nodes = []
    for i = 1:n
        push!(nodes, DisjointSetNode())
    end
    counter = 1
    while (counter <= n-k)
        edge = E[counter] # edge is tuple (w,u,v)
        if findset(nodes[edge[2]]) != findset(nodes[edge[3]]) # compare rootnode of node u and v
            push!(A,edge)
            my_union!(nodes[edge[2]], nodes[edge[3]])
        end
        counter = counter + 1
    end
    clusters = []
    c = 1 # cluster iterator
    push!(clusters, [])
    push!(clusters[1],[E[1][2]])
    clusterParent = findset(nodes[E[1][2]])
    println(nodes)
    for i = 2:n
        if findset(nodes[E[i][2]])==clusterParent
            push!(clusters[1][c], E[i][2])
        else
            c = c + 1
        end
    end
    if length(clusters[1]) < 2
        return clusters[1]
    else
    	return clusters
    end
end

max_stars_repo_name: evenlwanvik-student/algdat

max_issues_repo_name: evenlwanvik-student/algdat

max_forks_repo_name: evenlwanvik-student/algdat


text:
# This file was generated, do not modify it. # hide
using StatsBase #hide

@show describe(D.localrank)


text:
using Distributions
using Base.Test

d = NormalInverseGaussian(1.7, 1.8, 1.2, 2.3)

@test_approx_eq_eps params(d)[1] 1.7 0.000000001
@test_approx_eq_eps params(d)[2] 1.8 0.000000001
@test_approx_eq_eps params(d)[3] 1.2 0.000000001
@test_approx_eq_eps params(d)[4] 2.3 0.000000001

# The solution was computed using this R code:
# dnig(4.2, mu=1.7, alpha=1.8, beta=1.2, delta=2.3)
@test_approx_eq_eps pdf(d, 4.2) 0.2021958 0.0000001

# The solution was computed using this R code:
# dnig(4.8, mu=1.7, alpha=1.8, beta=1.2, delta=2.3, log=TRUE)
@test_approx_eq_eps logpdf(d, 4.8) -1.909973 0.000001

# The solution was computed using this R code:
# mean(rnig(100000000, mu=1.7, alpha=1.8, beta=1.2, delta=2.3))
@test_approx_eq_eps mean(d) 3.757509 0.001

# The solution was computed using this R code:
# var(rnig(100000000, mu=1.7, alpha=1.8, beta=1.2, delta=2.3))
@test_approx_eq_eps var(d) 3.085488 0.001

# The solution was computed using this R code:
# skewness(rnig(100000000, mu=1.7, alpha=1.8, beta=1.2, delta=2.3))
@test_approx_eq_eps skewness(d) 1.138959 0.001


text:
using Printf, Test, Plots

include("../../hybmc/mathutils/Helpers.jl")

include("../../hybmc/termstructures/YieldCurve.jl")
include("../../hybmc/models/HullWhiteModel.jl")
include("../../hybmc/models/AssetModel.jl")
include("../../hybmc/models/DeterministicModel.jl")
include("../../hybmc/models/HybridModel.jl")

include("../../hybmc/simulations/McSimulation.jl")
include("../../hybmc/simulations/Payoffs.jl")

using LinearAlgebra

function HwModel(rate=0.01,vol=0.0050,mean=0.03)
    curve = YieldCurve(rate)
    times = [ 10.0 ]
    vols  = [ vol  ]
    return HullWhiteModel(curve, mean, times, vols)
end

function setUp()
    domAlias = "EUR"
    domModel = HwModel(0.01, 0.0050, 0.03)
    forAliases = [ "USD", "GBP", "JPY" ]
    spotS0     = [   1.0,   2.0,   3.0 ]
    spotVol    = [   0.3,   0.3,   0.3 ]
    rates      = [  0.02,  0.03,  0.04 ]
    ratesVols  = [ 0.006, 0.007, 0.008 ]
    mean       = [  0.01,  0.01,  0.01 ]
    #
    forAssetModels = [
        AssetModel(S0, vol) for (S0, vol) in zip(spotS0,spotVol) ]
    forRatesModels = [
        HwModel(r,v,m) for (r,v,m) in zip(rates,ratesVols,mean) ]
    #
    nCorrs = 2 * size(forAliases)[1] + 1
    corr = Matrix{Float64}(I, nCorrs, nCorrs)
    corr[2,3] = 0.6  # USD quanto
    corr[3,2] = 0.6
    corr[4,5] = 0.8  # GBP quanto
    corr[5,4] = 0.8
    corr[6,7] = -0.7  # JPY quanto
    corr[7,6] = -0.7
    #
    return HybridModel(domAlias,domModel,forAliases,forAssetModels,forRatesModels,corr)
end

@testset "test_ModelSetup" begin
    model = setUp()
    @test stateSize(model) == 11
    @test factors(model)  == 7
    @test size(initialValues(model))[1] == stateSize(model)
    @test model.index == Dict([ ("USD",1), ("GBP",2), ("JPY",3) ])
    @test model.modelsStartIdx == [2, 3, 5, 6, 8, 9]
    C = model.L * transpose(model.L)
    @test norm(C-model.correlations) == 0.0
end

@testset "test_HybridModel" begin
    model = setUp()
    x = initialValues(model)
    @test asset(model,0.0,x,"USD") == model.forAssetModels[1].X0
    @test asset(model,0.0,x,"GBP") == model.forAssetModels[2].X0
    @test asset(model,0.0,x,"JPY") == model.forAssetModels[3].X0
    #
    @test zeroBond(model,0.0,5.0,x,"EUR") == discount(model.domRatesModel.yieldCurve,5.0)
    @test zeroBond(model,0.0,5.0,x,"USD") == discount(model.forRatesModels[1].yieldCurve,5.0)
    @test zeroBond(model,0.0,5.0,x,"GBP") == discount(model.forRatesModels[2].yieldCurve,5.0)
    @test zeroBond(model,0.0,5.0,x,"JPY") == discount(model.forRatesModels[3].yieldCurve,5.0)
end

@testset "test_HybridModelEvolution" begin
    model = setUp()
    times = [0.0]
    nPaths = 1
    seed = 314159265359
    # risk-neutral simulation
    mcSim = McSimulation(model,times,nPaths,seed,false)
    p = Path(mcSim,1)
    #
    @test asset(p,0.0,"USD") == model.forAssetModels[1].X0
    @test asset(p,0.0,"GBP") == model.forAssetModels[2].X0
    @test asset(p,0.0,"JPY") == model.forAssetModels[3].X0
    #
    @test zeroBond(p,0.0,5.0,"EUR") == discount(model.domRatesModel.yieldCurve,5.0)
    @test zeroBond(p,0.0,5.0,"USD") == discount(model.forRatesModels[1].yieldCurve,5.0)
    @test zeroBond(p,0.0,5.0,"GBP") == discount(model.forRatesModels[2].yieldCurve,5.0)
    @test zeroBond(p,0.0,5.0,"JPY") == discount(model.forRatesModels[3].yieldCurve,5.0)
end

@testset "test_HybridModelWithDeterministicRates" begin
    curve0 = YieldCurve(0.03)
    dcfModel0 = DcfModel(curve0,"EUR")
    #
    times = [0.0, 1.0]
    nPaths = 1
    seed = 314159265359
    # hybrid adjuster
    hybAdjTimes = [0.0, 1.0, 2.0]
    # simulate deterministic model only
    mcSim = McSimulation(dcfModel0,times,nPaths,seed,false)
    @test zeroBond(Path(mcSim,1),1.0,10.0,nothing) == discount(curve0,10.0)/discount(curve0,1.0)
    # simulate deterministic domestic model
    hwModel = HwModel(0.01,0.0050,0.03)
    asModel = AssetModel(1.0,0.30)
    corr =  Matrix{Float64}(I, 2, 2)
    model = HybridModel("EUR",dcfModel0,["USD"],[asModel],[hwModel],corr)
    hybVolAdj = hybridVolAdjuster(model,hybAdjTimes)
    model = HybridModel(model,hybAdjTimes,hybVolAdj)
    mcSim = McSimulation(model,times,nPaths,seed,false)
    p = Path(mcSim,1)
    @test asset(p,0.0,"USD") == model.forAssetModels[1].X0
    @test zeroBond(p,0.0,5.0,"EUR") == discount(curve0,5.0)
    @test zeroBond(p,0.0,5.0,"USD") == discount(model.forRatesModels[1].yieldCurve,5.0)
    # simulate deterministic foreign model
    dcfModel0 = DcfModel(curve0,"USD")
    model = HybridModel("EUR",hwModel,["USD"],[asModel],[dcfModel0],corr)
    hybVolAdj = hybridVolAdjuster(model,hybAdjTimes)
    model = HybridModel(model,hybAdjTimes,hybVolAdj)
    mcSim = McSimulation(model,times,nPaths,seed,false)
    p = Path(mcSim,1)
    @test asset(p,0.0,"USD") == model.forAssetModels[1].X0
    @test zeroBond(p,0.0,5.0,"EUR") == discount(model.domRatesModel.yieldCurve,5.0)
    @test zeroBond(p,0.0,5.0,"USD") == discount(curve0,5.0)
        # simulate deterministic domestic and foreign curve
    curve1 = YieldCurve(0.05)
    dcfModel0 = DcfModel(curve0,"EUR")
    dcfModel1 = DcfModel(curve1,"USD")
    corr =  Matrix{Float64}(I, 1, 1)
    model = HybridModel("EUR",dcfModel0,["USD"],[asModel],[dcfModel1],corr)
    hybVolAdj = hybridVolAdjuster(model,hybAdjTimes)
    model = HybridModel(model,hybAdjTimes,hybVolAdj)
    mcSim = McSimulation(model,times,nPaths,seed,false)
    p = Path(mcSim,1)
    @test asset(p,0.0,"USD") == model.forAssetModels[1].X0
    @test zeroBond(p,0.0,5.0,"EUR") == discount(curve0,5.0)
    @test zeroBond(p,0.0,5.0,"USD") == discount(curve1,5.0)
end

function fwd(mcSim::McSimulation,p::Payoff)
    samples = discountedAt(p,mcSim)
    fwd = mean(samples) /
        discount(mcSim.model.domRatesModel.yieldCurve,obsTime(p))
    err = std(samples) / sqrt(size(samples)[1]) /
        discount(mcSim.model.domRatesModel.yieldCurve,obsTime(p))
    return (fwd, err)
end

@testset "test_hybridSimulation" begin
    model = setUp()
    times = append!(Array(range(0,stop=10.0,length=11)), [10.5])
    nPaths = 2^13
    seed = 141592653593
    # risk-neutral simulation
    mcSim = McSimulation(model,times,nPaths,seed,false)
    #
    T = 10.0
    P = Pay(Fixed(1.0),T)
    (fw, err) = fwd(mcSim,P)
    # domestic numeraire
    @printf("1.0   @ %4.1lfy %8.6lf - mc_err = %8.6lf\n", T,fw,err)
    @test abs(fw - 1.0) < 1.5*err
    # foreign assets
    for (k, alias) in enumerate(model.forAliases)
        p = Asset(T,alias)
        xT = model.forAssetModels[k].X0 *
            discount(model.forRatesModels[k].yieldCurve,T) /
            discount(model.domRatesModel.yieldCurve,T)
        (fw, err) = fwd(mcSim,p)
        @printf("%s   @ %4.1lfy %8.6lf vs %8.6lf (curve) - mc_err = %8.6lf\n", alias,T,fw,xT,err)
        @test abs(fw - xT) < 2.0*err
    end
    # domestic Libor rate
    Tstart = 10.0
    Tend = 10.5
    L = Pay(LiborRate(T,Tstart,Tend,alias="EUR"),Tend)
    (fw, err) = fwd(mcSim,L)
    Lref = (discount(mcSim.model.domRatesModel.yieldCurve,Tstart) /
            discount(mcSim.model.domRatesModel.yieldCurve, Tend) - 1) /
            (Tend - Tstart)
    @printf("L_EUR @ %4.1lfy %8.6lf vs %8.6lf (curve) - mc_err = %8.6lf\n", T,fw,Lref,err)
    @test abs(fw - Lref) < 2.0*err
    # foreign Lbor rates
    for (k, alias) in enumerate(model.forAliases)
        L = Pay(LiborRate(T,Tstart,Tend,alias=alias)*Asset(Tend,alias),Tend)
        (fw, err) = fwd(mcSim,L)
        fw *= discount(mcSim.model.domRatesModel.yieldCurve, Tend) /
              discount(mcSim.model.forRatesModels[k].yieldCurve, Tend) /
              mcSim.model.forAssetModels[k].X0
        err *= discount(mcSim.model.domRatesModel.yieldCurve, Tend) /
               discount(mcSim.model.forRatesModels[k].yieldCurve, Tend) /
               mcSim.model.forAssetModels[k].X0
        Lref = (discount(mcSim.model.forRatesModels[k].yieldCurve,Tstart) /
                discount(mcSim.model.forRatesModels[k].yieldCurve, Tend) - 1) /
                (Tend - Tstart)
        @printf("L_%s @ %4.1lfy %8.6lf vs %8.6lf (curve) - mc_err = %8.6lf\n", alias,T,fw,Lref,err)
        @test abs(fw - Lref) < 2.0*err
    end
end

@testset "test_HybridVolAdjusterCalculation" begin
    # we set up a hybrid model consistent to QuantLib
    domAlias = "EUR"
    domModel = HwModel(0.01, 0.0050, 0.01)
    forAliases = [ "USD", "GBP" ]
    forAssetModels = [
        AssetModel(1.0, 0.30),
        AssetModel(2.0, 0.15) ]
    forRatesModels = [
        HwModel(0.02, 0.0060, 0.02),
        HwModel(0.02, 0.0070, 0.03) ]
    nCorrs = 2 * size(forAliases)[1] + 1
    corr = Matrix{Float64}(I, nCorrs, nCorrs)
    # [ EUR, USD-EUR, USD, GBP-EUR, GBP ] - +1 for Julia indeces
    #   0    1        2    3        4
    # USD-EUR - EUR
    corr[0+1,1+1] = 0.5
    corr[1+1,0+1] = 0.5
    # USD-EUR - USD
    corr[1+1,2+1] = -0.5
    corr[2+1,1+1] = -0.5
    # EUR - USD
    corr[0+1,2+1] = -0.5
    corr[2+1,0+1] = -0.5
    # GBP-EUR - EUR
    corr[0+1,3+1] = -0.5
    corr[3+1,0+1] = -0.5
    # GBP-EUR - GBP
    corr[3+1,4+1] = 0.5
    corr[4+1,3+1] = 0.5
    # EUR - GBP
    corr[0+1,4+1] = -0.8
    corr[4+1,0+1] = -0.8
    # USD - GBP
    corr[2+1,4+1] = 0.0
    corr[4+1,2+1] = 0.0
    model = HybridModel(domAlias,domModel,forAliases,forAssetModels,forRatesModels,corr)
    hybAdjTimes = Array(range(0.0, stop=20.0, length=21))
    hybVolAdj = hybridVolAdjuster(model,hybAdjTimes)
    model = HybridModel(model,hybAdjTimes,hybVolAdj)
    # plot(hybAdjTimes, transpose(hybVolAdj))
    # savefig("hybVolAdj.png")
    times = Array(range(0.0, stop=10.0, length=11))
    nPaths = 2^13
    seed = 314159265359
    # risk-neutral simulation
    mcSim = McSimulation(model,times,nPaths,seed,false)
    #
    T = 10.0
    for (k, alias) in enumerate(model.forAliases)
        # ATM forward
        xT = model.forAssetModels[k].X0 *
             discount(model.forRatesModels[k].yieldCurve, T) /
             discount(model.domRatesModel.yieldCurve, T)
        K = Fixed(xT)
        Z = Fixed(0.0)
        C = Pay(Max(Asset(T,alias)-K,Z),T)
        (fw, err) = fwd(mcSim,C)
        vol = BlackImpliedVol(fw,xT,xT,T,1.0)
        vega = BlackVega(xT,xT,vol,T)
        err /= vega
        volRef = model.forAssetModels[k].sigma
        @printf("C_%s @ %4.1lfy %8.6lf vs %8.6lf (curve) - mc_err = %8.6lf\n", alias,T,vol,volRef,err)
        @test abs(vol - volRef) < 3.0*err
        P = Pay(Max(K-Asset(T,alias),Z),T)
        (fw, err) = fwd(mcSim,P)
        vol = BlackImpliedVol(fw,xT,xT,T,-1.0)
        vega = BlackVega(xT,xT,vol,T)
        err /= vega
        volRef = model.forAssetModels[k].sigma
        @printf("P_%s @ %4.1lfy %8.6lf vs %8.6lf (curve) - mc_err = %8.6lf\n", alias,T,vol,volRef,err)
        @test abs(vol - volRef) < 3.0*err
    end
end

max_stars_repo_name: sschlenkrich/HybridMonteCarlo

max_issues_repo_name: sschlenkrich/HybridMonteCarlo

max_forks_repo_name: sschlenkrich/HybridMonteCarlo


text:
"""
Author: Jared J. Thomas
Date:   May 29, 2020

Test based on:
[1] An Aero-acoustic Noise Distribution Prediction Methodology for Offshore Wind Farms
by Jiufa Cao, Weijun Zhu, Xinbo Wu, Tongguang Wang, and Haoran Xu
"""

using FLOWFarm; const ff = FLOWFarm
using DelimitedFiles
using PyPlot

# load data
data = readdlm("inputfiles/velocity_def_row_of_10_turbs.txt",  ',', skipstart=4)

# load problem set up
include("./model_sets/model_set_2.jl")

# calculate wind turbine velocities and corresponding aerodynamic operational states
turbine_inflow_velcities, turbine_ct, turbine_ai, turbine_local_ti = ff.turbine_velocities_one_direction(turbine_x, turbine_y, turbine_z, rotor_diameter, hub_height, turbine_yaw,
sorted_turbine_index, ct_model, rotor_sample_points_y, rotor_sample_points_z, windresource,
model_set, wind_farm_state_id=1)

# calculate the power production of each wind turbine
ff.turbine_powers_one_direction(generator_efficiency, cut_in_speed, cut_out_speed, rated_speed, rated_power, rotor_diameter, turbine_inflow_velcities, air_density, power_model)

# set up the point locations to be used in generating the plot
stepsize = 5
xrange = 1:stepsize:10*rotor_diameter[1]*nturbines
yrange = -1*rotor_diameter[1]*nturbines:stepsize:1*rotor_diameter[1]*nturbines
zrange = hub_height[1]:stepsize:hub_height[1]

# calculate the wind speed at each point in the flow field
velh = ff.calculate_flow_field(xrange, yrange, zrange,
model_set, turbine_x, turbine_y, turbine_z, turbine_yaw, turbine_ct, turbine_ai,
rotor_diameter, hub_height, turbine_local_ti, sorted_turbine_index, turbine_inflow_velcities,
windresource; wind_farm_state_id=1)

# visualize the resulting data
flowfieldplot = contourf(xrange, yrange, velh, cmap="Blues_r")


text:
# ---
# title: 561. Array Partition I
# id: problem561
# author: Indigo
# date: 2021-02-16
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/array-partition-i/description/>
# hidden: true
# ---
#
# Given an integer array `nums` of `2n` integers, group these integers into `n`
# pairs `(a1, b1), (a2, b2), ..., (an, bn)` such that the sum of `min(ai, bi)`
# for all `i` is **maximized**. Return _the maximized sum_.
#
#
#
# **Example 1:**
#
#
#
#     Input: nums = [1,4,3,2]
#     Output: 4
#     Explanation: All possible pairings (ignoring the ordering of elements) are:
#     1. (1, 4), (2, 3) -> min(1, 4) + min(2, 3) = 1 + 2 = 3
#     2. (1, 3), (2, 4) -> min(1, 3) + min(2, 4) = 1 + 2 = 3
#     3. (1, 2), (3, 4) -> min(1, 2) + min(3, 4) = 1 + 3 = 4
#     So the maximum possible sum is 4.
#
# **Example 2:**
#
#
#
#     Input: nums = [6,2,6,5,1,2]
#     Output: 9
#     Explanation: The optimal pairing is (2, 1), (2, 5), (6, 6). min(2, 1) + min(2, 5) + min(6, 6) = 1 + 2 + 6 = 9.
#
#
#
#
# **Constraints:**
#
#   * `1 <= n <= 104`
#   * `nums.length == 2 * n`
#   * `-104 <= nums[i] <= 104`
#
#
## @lc code=start
using LeetCode

function array_pair_sum(nums::Vector{Int})
    sort!(nums)
    sum(nums[1:2:end])
end
## @lc code=end


text:
# =========================================================================================
"""
    function φΩ{T}(
        ξ    ::Vector{T},
        model::NonlocalXieModel1{T}
    )

Computes the interior nonlocal electrostatic potential ``φ_Ω`` for the given observation
point ``ξ``.

# Unit
``V = \\frac{C}{F}``

# Return type
`T`

!!! warning
    This function does not verify whether ξ is located inside of the sphere!
"""
function φΩ(ξ::Vector{T}, model::NonlocalXieModel1{T}) where T
    a  = model.radius
    λ  = model.params.λ
    εΩ = model.params.εΩ
    εΣ = model.params.εΣ
    ε∞ = model.params.ε∞
    A₃ = model.A₃
    r  = norm(ξ)

    φ = zero(T)
    for (qi, q) in enumerate(model.charges)   # Eq. (17a)

        # if ξ is close to the origin, all terms for n > 0 become negligible
        if r < 1e-10
            φ += A₃[1, qi]
            continue
        end

        # if q is close to the origin, compute nonlocal Born potential
        # (see born/potentials.jl)
        if norm(q.pos) < 1e-10
            ν = √(εΣ/ε∞) * a/λ
            φ += q.val * (1/r + 1/a/εΣ * (1 - εΣ + (εΣ - ε∞)/ε∞ * sinh(ν)/ν * exp(-ν)))
            continue
        end

        # otherwise, use Eq. (18)
        P = legendre(model.len, cos(ξ, q.pos, r))
        φj = zero(T)
        for n in 0:model.len-1
            φj += A₃[n+1, qi] * r^n * P(n)
        end
        φ += φj * q.val
    end

    (φ + φmol(ξ, model.charges) / 4π / εΩ) * T(ec/ε0)
end

# =========================================================================================
"""
    function φΣ{T}(
        ξ    ::Vector{T},
        model::NonlocalXieModel1{T}
    )

Computes the exterior nonlocal electrostatic potential ``φ_Σ`` for the given observation
point ``ξ``.

# Unit
``V = \\frac{C}{F}``

# Return type
`T`

!!! warning
    This function does not verify whether ξ is located outside of the sphere!
"""
function φΣ(ξ::Vector{T}, model::NonlocalXieModel1{T}) where T
    a  = model.radius
    λ  = model.params.λ
    εΣ = model.params.εΣ
    ε∞ = model.params.ε∞
    A₁ = model.A₁
    A₂ = model.A₂
    κ  = λ \ √(εΣ/ε∞)
    r  = norm(ξ)

    kᵣ = spherical_besselk(model.len, κ * r)

    φ = zero(T)
    for (qi, q) in enumerate(model.charges)   # Eq. (17a)

        # if q is close to origin, compute nonlocal Born potential
        # (see born/potentials.jl)
        if norm(q.pos) < 1e-10
            ν = √(εΣ/ε∞) * a/λ
            φ += q.val / εΣ / r * (1 + (εΣ - ε∞)/ε∞ * sinh(ν)/ν * exp(-ν * r/a))
            continue
        end

        # otherwise, use Eq. (18)
        P = legendre(model.len, cos(ξ, q.pos, r))
        φj = zero(T)
        for n in 0:model.len-1
            φj += (ε∞ - εΣ)/ε∞ * A₂[n+1, qi] * kᵣ(n) * P(n) +
                  A₁[n+1, qi] / r^(n+1) * P(n)
        end
        φ += φj * q.val
    end

    φ * T(ec/ε0)
end


text:
# module Ecto.MirrorTypes

abstract MirrorModel

type MirrorField
    name::Symbol
    typ::Type
    opts::Dict
end

type MirrorPredicate
    body::Vector
end

type MirrorStatement
end

type MirrorNaming
   plain::Symbol
   mirror::Symbol
   schema::Symbol
   function MirrorNaming(::Type{Val{:plain}}, T::Type)
       plain = T.name.name
       MirrorNaming(plain)
   end
   function MirrorNaming(::Type{Val{:mirror}}, T::Type)
       mirror = T.name.name
       plain = Symbol(string(mirror)[2:end])
       MirrorNaming(plain)
   end
   function MirrorNaming(::Type{Val{:schema}}, M::Module)
       modulename = Base.module_name(M)
       plain = Symbol(replace(string(modulename), r"Schema", "") )
       MirrorNaming(plain)
   end
   function MirrorNaming(plain::Symbol)
       new(plain, Symbol(:_, plain), Symbol(plain, :Schema))
   end
end

type MirrorStruct
    naming::MirrorNaming
    fields::Vector{MirrorField}
end

type MirrorStructError
    message
end

function schema_to_mirrorstruct(modul::Module)::MirrorStruct
    if isdefined(modul, :schema) && isdefined(Ecto.Schema, :__attributes__)
        naming = MirrorNaming(Val{:schema}, modul)
        ecotofields = Ecto.Schema.get_attribute(modul, :ecto_fields)
        fields = Vector{MirrorField}()
        for (k,v) in ecotofields
            name = k
            typ = Ecto.Typ.julia_type(Val{v})
            opts = Dict()
            field = MirrorField(name, typ, opts)
            push!(fields, field)
        end
        MirrorStruct(naming, fields)
    else
        throw(MirrorStructError(""))
    end
end

function plain_to_mirrorstruct(T::Type)::MirrorStruct
    naming = MirrorNaming(Val{:plain}, T)
    names = fieldnames(T)
    fields = Vector{MirrorField}()
    if !in(:id, names)
        field = MirrorField(:id, Int, Dict())
        push!(fields, field)
    end
    for name in names
        typ = fieldtype(T, name)
        opts = Dict()
        field = MirrorField(name, typ, opts)
        push!(fields, field)
    end
    MirrorStruct(naming, fields)
end

function mirror_to_mirrorstruct{MM<:MirrorModel}(mirror::MM)::MirrorStruct
    T = typeof(mirror)
    naming = MirrorNaming(Val{:mirror}, T)
    names = fieldnames(T)
    fields = Vector{MirrorField}()
    if !in(:id, names)
        field = MirrorField(:id, Int, Dict())
        push!(fields, field)
    end
    for name in names
        field = getfield(mirror, name)
        push!(fields, field)
    end
    MirrorStruct(naming, fields)
end

function mirrorstruct_to_schema(struct::MirrorStruct)::Module
    modul = MirrorTypes
    typename = struct.naming.schema
    isdefined(modul, typename) && return getfield(modul, typename)
    plural = Ecto.English.pluralize(string(struct.naming.plain))
    lines = String[]
    push!(lines, "module $typename")
    push!(lines, string("using ", Ecto.Schema))
    push!(lines, """schema("$plural") do""")
    for field in struct.fields
        push!(lines, string("    field(", repr(field.name), ", ", repr(Ecto.Typ.elixir_type(field.typ)), ")"))
    end
    push!(lines, "end")
    push!(lines, "end")
    schema_declaration = join(lines, '\n')
    # println(schema_declaration)
    eval(modul, parse(schema_declaration))
end

function mirrorstruct_to_plain(struct::MirrorStruct)::Type
    modul = MirrorTypes
    typename = struct.naming.plain
    isdefined(modul, typename) && return getfield(modul, typename)
    lines = String[]
    push!(lines, "type $typename")
    for field in struct.fields
        push!(lines, string("    ", field.name, "::", field.typ))
    end
    push!(lines, "end")
    plain_declaration = join(lines, '\n')
    # println(plain_declaration)
    eval(modul, parse(plain_declaration))
    getfield(modul, typename)
end

function mirrorstruct_to_mirror(struct::MirrorStruct)
    modul = MirrorTypes
    typename = struct.naming.mirror
    isdefined(modul, typename) && return getfield(modul, typename)(struct.fields...)
    lines = String[]
    push!(lines, "type $typename <: Ecto.MirrorTypes.MirrorModel")
    for field in struct.fields
        push!(lines, string("    ", field.name, "::", MirrorField))
    end
    push!(lines, "end")
    mirror_declaration = join(lines, '\n')
    # println(mirror_declaration)
    eval(modul, parse(mirror_declaration))
    T = getfield(modul, typename)
    T(struct.fields...)
end


text:
module ESNICE_examples
using Statistics
using FinEtools
using FinEtools.MeshImportModule
using FinEtools.MeshExportModule
using LinearAlgebra: dot
using Arpack
using LinearAlgebra
using SparseArrays
using PGFPlotsX
using Test
using StatsBase

function ESNICE_energies()
    E = 1e6*phun("PA");
    nu = 0.0;
    L = 2*phun("M");
    hs = L * collect(10 .^ range(-4.0, stop=0.0, length = 10))
    mag = 0.001

    rs = Float64[]
    PEs = Float64[]
    APEs = Float64[]
    for h = hs
        xs = collect(linearspace(0.0, L, 2))
        ys = collect(linearspace(0.0, h, 2))
        zs = collect(linearspace(0.0, h, 2))
        global fens; global fes
        fens, fes = T4blockx(xs, ys, zs, :a)
        fens.xyz[:, 3] .-= h/2
        MR = DeforModelRed3D
        global geom = NodalField(fens.xyz)
        global u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field
        numberdofs!(u)

        l1 = selectnode(fens; plane = [1.0 0.0 0.0 0.0], thickness = h/1000)
        for i = l1
            x,y,z = geom.values[i, :]
            u.values[i, 1] = z * mag
        end
        l1 = selectnode(fens; plane = [1.0 0.0 0.0 L], thickness = h/1000)
        for i = l1
            x,y,z = geom.values[i, :]
            u.values[i, 1] = -z * mag
        end
        @show h, u.values
        material = MatDeforElastIso(MR, E, nu)
        global femm = FEMMDeforLinear(MR, IntegDomain(fes, TetRule(1)), material)
        associategeometry!(femm, geom)

        K = stiffness(femm, geom, u)
        U = gathersysvec(u)
        PE = dot(U, K * U) / 2.0
        push!(PEs, PE)
        I = h * h^3 / 12
        APE = 2 * E * I * mag^2 / L
        push!(APEs, APE)

        ars = []
        for i = 1:count(fes)
            res  = FinEtools.FEMMDeforLinearESNICEModule.aspectratio(geom.values[collect(fes.conn[i]), :])
            ar = sort([res[1], res[2], res[3], res[4]])
            push!(ars, mean([ar[2:3]...]))
        end
        # @show h/L, ars
        # push!(rs, minimum(ars))
        push!(rs, h/L)
    end

    @show rPE = PEs ./ APEs

    # Least-squares fit
    A = hcat([-log10(r) for r in rs], [-1 for r in rPE])
    b = [log10(r-1) for r in rPE]
    p = A \ b
    @show p
    a = p[1]
    b = 10^p[2]
    @show a, b

    @pgf a = Axis({
            xlabel = "Aspect ratio",
            ylabel = "Relative Potential Energy",
            grid="major",
            legend_pos  = "north east"
        },
        Plot({color="red"}, Table([:x => log10.(vec(rs)), :y => log10.(vec(rPE))])),
        Plot({"only marks", mark="x"}, Table([:x => log10.(vec(rs)), :y => log10.(vec([1/(b*r^a)+1 for r in rs]))]))
        )
    display(a)

    # fld = fieldfromintegpoints(femm, geom, u, :Cauchy, 1)
    # File =  "mt4energy2.vtk"
    # vtkexportmesh(File, fens, fes; scalars=[("sigmax", fld.values)], vectors=[("u", u.values)])
    # @async run(`"paraview.exe" $File`)
    true
end

function ESNICE_vibration()
    E = 70000*phun("MPa");
    nu = 0.33;
    rho = 2700*phun("KG/M^3");
    radius = 0.5*phun("ft");
    neigvs = 20                   # how many eigenvalues
    OmegaShift = (10.0*2*pi)^2;

    MR = DeforModelRed3D
    output = import_ABAQUS("alum_cyl.inp")
    fens, fes = output["fens"], output["fesets"][1]
    fens.xyz .*= phun("mm") # The input is provided in SI(mm) units
    fens, fes = T10toT4(fens, fes)

    geom = NodalField(fens.xyz)
    u = NodalField(zeros(size(fens.xyz,1),3)) # displacement field

    numberdofs!(u)

    material = MatDeforElastIso(MR, rho, E, nu, 0.0)

    femm = FEMMDeforLinearESNICET4(MR, IntegDomain(fes, NodalSimplexRule(3)), material)
    associategeometry!(femm,  geom)
    @show minimum(vec(femm.nphis)), maximum(vec(femm.nphis))
    @pgf a = Axis({
            xlabel = "Entity",
            ylabel = "Stabilization factor",
            grid="major",
            legend_pos  = "north east"
        },
        Plot({mark="circle"}, Table([:x => vec(1:count(fes)), :y => vec(femm.ephis)]))
        )
    display(a)
    K  = stiffness(femm, geom, u)
    M = mass(femm, geom, u)
    d,v,nev,nconv = eigs(K+OmegaShift*M, M; nev=neigvs, which=:SM)
    d = d .- OmegaShift;
    fs = real(sqrt.(complex(d)))/(2*pi)
    println("Eigenvalues: $fs [Hz]")

    vectors = []
    for i = 7:length(fs)
        scattersysvec!(u, v[:, i])
        push!(vectors, ("Mode_$i", deepcopy(u.values)))
    end
    File  =   "alum_cyl_mode_shapes.vtk"
    vtkexportmesh(File, connasarray(fes), fens.xyz, FinEtools.MeshExportModule.T4; vectors = vectors)
    @async run(`"paraview.exe" $File`)

    true
end # function

function allrun()
    println("#####################################################")
    println("# ESNICE_energies ")
    ESNICE_energies()
    println("#####################################################")
    println("# ESNICE_vibration ")
    ESNICE_vibration()
    return true
end # function allrun

end # module LE11NAFEMS_examples


text:
function get_unitcell_hyperhoneycomb() :: Unitcell

    # define list of basis sites
    basis    = Vector{SVector{3, Float64}}(undef, 4)
    basis[1] = SVector{3, Float64}(0.0,  0.0, 0.0)
    basis[2] = SVector{3, Float64}(1.0,  1.0, 0.0)
    basis[3] = SVector{3, Float64}(1.0,  2.0, 1.0)
    basis[4] = SVector{3, Float64}(0.0, -1.0, 1.0)

    # define list of Bravais vectors
    vectors    = Vector{SVector{3, Float64}}(undef, 3)
    vectors[1] = SVector{3, Float64}(-1.0, 1.0, -2.0)
    vectors[2] = SVector{3, Float64}(-1.0, 1.0,  2.0)
    vectors[3] = SVector{3, Float64}( 2.0, 4.0,  0.0)

    # define list of bonds for each basis site
    bonds     = Vector{Vector{SVector{4, Int64}}}(undef, length(basis))
    bonds1    = Vector{SVector{4, Int64}}(undef, 3)
    bonds1[1] = SVector{4, Int64}( 0,  0,  0,  1)
    bonds1[2] = SVector{4, Int64}( 0,  0,  0,  3)
    bonds1[3] = SVector{4, Int64}( 1,  0,  0,  3)
    bonds[1]  = bonds1
    bonds2    = Vector{SVector{4, Int64}}(undef, 3)
    bonds2[1] = SVector{4, Int64}( 0,  0,  0, -1)
    bonds2[2] = SVector{4, Int64}( 0,  0,  0,  1)
    bonds2[3] = SVector{4, Int64}( 0, -1,  0,  1)
    bonds[2]  = bonds2
    bonds3    = Vector{SVector{4, Int64}}(undef, 3)
    bonds3[1] = SVector{4, Int64}( 0,  0,  0, -1)
    bonds3[2] = SVector{4, Int64}( 0,  1,  0, -1)
    bonds3[3] = SVector{4, Int64}( 0,  0,  1,  1)
    bonds[3]  = bonds3
    bonds4    = Vector{SVector{4, Int64}}(undef, 3)
    bonds4[1] = SVector{4, Int64}( 0,  0,  0, -3)
    bonds4[2] = SVector{4, Int64}(-1,  0,  0, -3)
    bonds4[3] = SVector{4, Int64}( 0,  0, -1, -1)
    bonds[4]  = bonds4

    # build unitcell
    uc = Unitcell(basis, vectors, bonds)

    return uc
end


text:
@block SimonDanisch ["dataset_examples"] begin
    @cell "WorldClim visualization" [visualization, dataset, bigdata] begin
        using FileIO, GeometryTypes, Colors, GDAL

        #=
        This example requires the GDAL package, from https://github.com/JuliaGeo/GDAL.jl
        For more information about GDAL, see the official documentation at: https://gdal.org/
        =#

        # register GDAL drivers
        GDAL.allregister()

        # function to check if a file is a .tif file
        istiff(x) = endswith(x, ".tif")

        # set up 7zip
        exe7z = is_windows() ? joinpath(JULIA_HOME, "7z.exe") : joinpath(JULIA_HOME, "7z")
        unzip(in, out) = run(`$exe7z x -y $in -o$out`)

        # function to read the raster data from the GeoTIFF
        function loadf0(x)
            img = GDAL.open(x, GDAL.GA_ReadOnly)
            band = GDAL.getrasterband(img, 1)
            xsize = GDAL.getrasterbandxsize(band)
            ysize = GDAL.getrasterbandysize(band)
            data = Array{Float32}(xsize, ysize)
            GDAL.rasterio(band, GDAL.GF_Read, 0, 0, xsize, ysize, data, xsize, ysize, GDAL.GDT_Float32, 0, 0)
            data'
        end

        # since we cannot re-distribute the dataset, this function grabs the dataset from the host server
        function load_dataset(name)
            # get the dataset from:
            # http://worldclim.org/version2
            """
            This is WorldClim 2.0 Beta version 1 (June 2016) downloaded from http://worldclim.org
            They represent average monthly climate data for 1970-2000.
            The data were created by Steve Fick and Robert Hijmans
            You are not allowed to redistribute these data.
            """
            if !isfile("$name.zip")
                # This might fail on windows - just try again a couple of times -.-
                download("http://biogeo.ucdavis.edu/data/worldclim/v2.0/tif/base/wc2.0_10m_$name.zip", "$name.zip")
            end
            if !isdir(name)
                unzip("$name.zip", name)
            end
            loadf0.(filter(istiff, joinpath.(name, readdir(name))))
        end

        # load the actual datasets!
        water = load_dataset("prec")
        temperature = load_dataset("tmax")

        # calculate geometries
        m = GLNormalUVMesh(Sphere(Point3f0(0), 1f0), 200)
        p = decompose(Point3f0, m)
        uv = decompose(UV{Float32}, m)
        norms = decompose(Normal{3, Float32}, m)

        # plot the temperature as color map on the globe
        cmap = [:darkblue, :deepskyblue2, :deepskyblue, :gold, :tomato3, :red, :darkred]
        scene = Makie.mesh(m, color = temperature[10], colorrange = (-50, 50), colormap = cmap, shading = true, show_axis = false)
        temp_plot = scene[end];

        # function to scale precipitation to a suitable marker size
        function to_msize(uv, water, xysize = 0.008, normalization = 908f0 * 4f0)
            markersize = map(uv) do uv
                wh = size(water) .- 1
                x, y = round.(Int, Tuple(uv) .* wh) .+ 1
                val = water[x, y] / normalization
                val = val < 0.0 ? 0f0 : val
                Vec3f0(xysize, xysize, val)
            end
        end

        # plot precipitation as vertical bars
        meshscatter!(scene,
            p, rotations = norms,
            marker = Rect3D(Vec3f0(0), Vec3f0(1)), markersize = to_msize(uv, water[1]), raw = true,
            color = (:skyblue2, 0.9)
        )
        wplot = scene[end]

        # update eye position
        eye_position, lookat, upvector = Vec3f0(0.5, 0.8, 2.5), Vec3f0(0), Vec3f0(0, 1, 0)
        update_cam!(scene, eye_position, lookat, upvector)
        scene

        # save animation
        r = record(scene, @outputfile(mp4), 0:(11*4)) do i
            # Make simulation slower. TODO figure out how do this nicely with ffmpeg
            if i % 4 == 0
                i2 = (i ÷ 4) + 1
                temp_plot[:color] = temperature[i2]
                wplot[:markersize] = to_msize(uv, water[i2])
            end
        end
        return r
    end

end


text:
struct JointDistribution <: RandomUQInput
    marginals::Array{RandomVariable}
    copula::Copula

    function JointDistribution(marginals::Array{RandomVariable}, copula::Copula)
        length(marginals) == dimensions(copula) ||
            error("Dimension mismatch between copula and marginals")

        return new(marginals, copula)
    end
end

function sample(jd::JointDistribution, n::Integer=1)
    u = sample(jd.copula, n)

    samples = DataFrame()

    for (i, rv) in enumerate(jd.marginals)
        samples[!, rv.name] = quantile.(rv.dist, u[:, i])
    end

    return samples
end

function to_physical_space!(jd::JointDistribution, x::DataFrame)
    correlated_cdf = to_copula_space(jd.copula, Matrix{Float64}(x[:, names(jd)]))
    for (i, rv) in enumerate(jd.marginals)
        x[!, rv.name] = quantile.(rv.dist, correlated_cdf[:, i])
    end
    return nothing
end

function to_standard_normal_space!(jd::JointDistribution, x::DataFrame)
    for rv in jd.marginals
        x[!, rv.name] = cdf.(rv.dist, x[:, rv.name])
    end
    uncorrelated_stdnorm = to_standard_normal_space(
        jd.copula, Matrix{Float64}(x[:, names(jd)])
    )
    for (i, rv) in enumerate(jd.marginals)
        x[!, rv.name] = uncorrelated_stdnorm[:, i]
    end
    return nothing
end

names(jd::JointDistribution) = vec(map(x -> x.name, jd.marginals))

mean(jd::JointDistribution) = mean(jd.marginals)

dimensions(jd::JointDistribution) = dimensions(jd.copula)

max_stars_repo_name: FriesischScott/julia-uq

max_issues_repo_name: FriesischScott/julia-uq

max_forks_repo_name: FriesischScott/julia-uq


text:
@testset "7.6 Inverse hyperbolic cosecant" begin
    include("7.6.1 u (a+b arccsch(c x))^n.jl")
    include("7.6.2 Inverse hyperbolic cosecant functions.jl")
end

max_stars_repo_path: test/7 Inverse hyperbolic functions/7.6 Inverse hyperbolic cosecant/.jl


text:
using Documenter
import GraphIdx

DocMeta.setdocmeta!(GraphIdx, :DocTestSetup, :(import GraphIdx);
                    recursive=false, warn=true)
DocMeta.setdocmeta!(GraphIdx.Tree, :DocTestSetup, :(import GraphIdx.Tree);
                    recursive=false, warn=true)
doctest(GraphIdx; manual=false)

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
# ---
# title: 1640. Check Array Formation Through Concatenation
# id: problem1640
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Array, Sort
# link: <https://leetcode.com/problems/check-array-formation-through-concatenation/description/>
# hidden: true
# ---
#
# You are given an array of **distinct** integers `arr` and an array of integer
# arrays `pieces`, where the integers in `pieces` are **distinct**. Your goal is
# to form `arr` by concatenating the arrays in `pieces` **in any order**.
# However, you are **not** allowed to reorder the integers in each array
# `pieces[i]`.
#
# Return `true` _if it is possible_ _to form the array_`arr` _from_`pieces`.
# Otherwise, return `false`.
#
#
#
# **Example 1:**
#
#
#
#     Input: arr = [85], pieces = [[85]]
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: arr = [15,88], pieces = [[88],[15]]
#     Output: true
#     Explanation: Concatenate [15] then [88]
#
#
# **Example 3:**
#
#
#
#     Input: arr = [49,18,16], pieces = [[16,18,49]]
#     Output: false
#     Explanation: Even though the numbers match, we cannot reorder pieces[0].
#
#
# **Example 4:**
#
#
#
#     Input: arr = [91,4,64,78], pieces = [[78],[4,64],[91]]
#     Output: true
#     Explanation: Concatenate [91] then [4,64] then [78]
#
# **Example 5:**
#
#
#
#     Input: arr = [1,3,5,7], pieces = [[2,4,6,8]]
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `1 <= pieces.length <= arr.length <= 100`
#   * `sum(pieces[i].length) == arr.length`
#   * `1 <= pieces[i].length <= arr.length`
#   * `1 <= arr[i], pieces[i][j] <= 100`
#   * The integers in `arr` are **distinct**.
#   * The integers in `pieces` are **distinct**  (i.e., If we flatten pieces in a 1D array, all the integers in this array are distinct).
#
#
## @lc code=start
using LeetCode

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


text:
using EnsembleInference
using Test

@testset "EnsembleInference.jl" begin
    include("utility.jl")
    include("dirac.jl")
    include("haar.jl")
    include("diffusion_normal.jl")
    include("representation.jl")
end


text:
module test_flux_basic

using Test
using Flux.Tracker

W = # rand(2, 5)
[0.857747   0.291713  0.179873  0.938979  0.51022
 0.0852085  0.977716  0.246164  0.460672  0.772312]

b = # rand(2)
[0.663369
 0.132996]

predict(x) = W*x .+ b
loss(x, y) = sum((predict(x) .- y).^2)

x, y = # rand(5), rand(2)
([0.496864, 0.947507, 0.874288, 0.251528, 0.192234], [0.901991, 0.0802404])

@test loss(x, y) ≈ 3.166070569359008

W = param(W)
b = param(b)
l = loss(x, y)

@test W isa TrackedArray
@test b isa TrackedArray
@test l isa Tracker.TrackedReal
@test l.tracker.grad == 0

@test W.grad == [0.0  0.0  0.0  0.0  0.0
                 0.0  0.0  0.0  0.0  0.0]

back!(l)
@test W.grad ≈ [0.949491  1.81065  1.67074  0.480662  0.367353
                1.49163   2.84449  2.62468  0.755107  0.577102] atol=0.00001

W.data .-= 0.1(W.grad)
@test W.data ≈ [0.762798   0.110648   0.0127994  0.890913  0.473485
               -0.0639541  0.693267  -0.0163043  0.385161  0.714602] atol=0.00001

end # module test_flux_basic


text:
using BSONqs
using Test

const BSON = BSONqs

roundtrip_equal(x, flags...) = BSON.roundtrip(x, flags...) == x

abstract type Bar end

mutable struct Foo
  x
end

struct T{A}
  x::A
end

struct S end

struct Baz <: Bar end

@testset "BSON" begin

@testset "Primitive Types" begin
  @test roundtrip_equal(nothing)
  @test roundtrip_equal(1)
  @test roundtrip_equal(Dict(:a => 1,:b => 2))
  @test roundtrip_equal(UInt8[1,2,3])
  @test roundtrip_equal("b")
  @test roundtrip_equal([1,"b"])
  @test roundtrip_equal([[1, 2, 3], "foo"])
  @test roundtrip_equal(Tuple, :compat)
end

@testset "Complex Types" begin
  @test roundtrip_equal(:foo)
  @test roundtrip_equal(Int64)
  @test roundtrip_equal(Complex{Float32})
  @test roundtrip_equal(Complex, :compat)
  @test roundtrip_equal(Array, :compat)
  @test roundtrip_equal([1,2,3])
  @test roundtrip_equal(rand(2,3))
  @test roundtrip_equal(Array{Real}(rand(2,3)), :compat)
  @test roundtrip_equal(1+2im)
  @test roundtrip_equal(Nothing[])
  @test roundtrip_equal(S[])
  @test roundtrip_equal(fill(nothing, (3,2)))
  @test roundtrip_equal(fill(S(), (1,3)))
  @test roundtrip_equal(Set([1,2,3]))
  @test roundtrip_equal(Dict("a"=>1))
  @test roundtrip_equal(T(()))
  @test roundtrip_equal(T{Bar}(Baz()))
  @test roundtrip_equal(T{Symbol}(:wibble))
end

@testset "Circular References" begin
  x = [1,2,3]
  (x1, x2) = BSON.roundtrip((x,x))
  @test x1 == x
  @test x1 === x2

  d = Dict{Symbol,Any}(:a=>1)
  d[:d] = d
  d = BSON.roundtrip(d)
  @test d[:a] == 1
  @test d[:d] === d

  x = Foo(1)
  x.x = x
  x = BSON.roundtrip(x)
  @test x.x === x
end

@testset "Specify root document types" begin
  @test roundtrip_equal(T{String}("foo"), :set_root_type)
  @test roundtrip_equal(rand(7), :set_root_type)
end

end


text:
# EXAMPLE: Robertson's reaction from pages 3-4 of Solving Ordinary Differential Equations II
#     (Stiff and Differential-Algebraic Problems)

struct RobertsonSystem
    unused_var::Float64  # to make generalization more obvious to readers
    de::Function
    function RobertsonSystem(de::Function)
        return new(1.0, de)
    end
end

function stiff_de!(xx::Vector{T}, x::Vector{T}, s::RobertsonSystem, t::Float64=0.0) where {T}
    c1 = 0.04
    c2 = 1.0e4
    c3 = 3.0e7
    y1 = x[1]
    y2 = x[2]
    y3 = x[3]
    xx[1] = -c1 * y1 + c2 * y2 * y3
    xx[2] =  c1 * y1 - c2 * y2 * y3 - c3 * y2^2
    xx[3] =                           c3 * y2^2
    return nothing
end

my_stiff_struct = RobertsonSystem(stiff_de!)
x0 = [1.0, 0.0, 0.0]  # see Eq. 1.4
initial_h = 1.0e-4

rr = makeRadauIntegrator(my_stiff_struct, x0, 1.0e-16, 2, 1)
update_h!(rr, initial_h)  # force it to start with h = 1.0e-4
rr.rule.s = 2  # force it to start with Radau5
rr.step.h_max = Inf  # don't want timesteps truncated
n_steps = 10

xx = [x0]  # to avoid putting the loop in a function
t = [0.0]  # to avoid putting the loop in a function
for k = 1:n_steps
    x0 = xx[1]
    h, x_final = solveRadau(rr, x0)
    xx[1] = x_final
    t[1] += h
end

println("time: ", t[1])
println("second term: ", xx[1][2])
println("compare to Fig 1.3")

@testset "Robertson's reaction" begin
    @test 3.45e-5 < xx[1][2] < 3.7e-5
    @test t[1] == initial_h * (2^(n_steps) - 1)  # it can at most double h every step
end


text:
__precompile__()
module UnregisteredWithoutProject

if !isfile(joinpath(@__DIR__, "..", "deps", "deps.jl"))
    error("UnregisteredWithoutProject is not installed correctly")
end

export f
f(x) = x

end


text:
# abstract type for scalar-valued kernels
abstract type ScalarKernel <: AbstractKernel end

# scalar-valued exponential kernel
struct ExponentialKernel{T<:Real,M<:Metric} <: ScalarKernel
    """Inverse length-scale."""
    γ::T
    """Distance function."""
    distance::M

    function ExponentialKernel{T,M}(γ::T, distance::M) where {T,M}
        γ > zero(γ) ||
            throw(ArgumentError("inverse length-scale has to be positive"))

        new{T,M}(γ, distance)
    end
end

ExponentialKernel(γ::Real, distance::Metric = Euclidean()) =
    ExponentialKernel{typeof(γ),typeof(distance)}(γ, distance)
ExponentialKernel(distance::Metric = Euclidean()) = ExponentialKernel(1, distance)

(kernel::ExponentialKernel)(x, y) = exp(- kernel.γ * evaluate(kernel.distance, x, y))

const LaplacianKernel = ExponentialKernel

# scalar-valued squared exponential kernel
struct SquaredExponentialKernel{T<:Real,M<:Metric} <: ScalarKernel
    """Squared inverse length-scale."""
    γ::T
    """Distance function."""
    distance::M

    function SquaredExponentialKernel{T,M}(γ::T, distance::M) where {T,M}
        γ > zero(γ) ||
            throw(ArgumentError("squared inverse length-scale has to be positive"))

        new{T,M}(γ, distance)
    end
end

SquaredExponentialKernel(γ::Real, distance::Metric = Euclidean()) =
    SquaredExponentialKernel{typeof(γ),typeof(distance)}(γ, distance)
SquaredExponentialKernel(distance::Metric = Euclidean()) =
    SquaredExponentialKernel(1, distance)

(kernel::SquaredExponentialKernel)(x, y) =
    exp(- kernel.γ * evaluate(kernel.distance, x, y)^2)

const GaussianKernel = SquaredExponentialKernel


text:
"""
# LocalRegistry

Create and maintain local registries for Julia packages. This package
is intended to provide a simple but manual workflow for maintaining
small local registries (private or public) without making any
assumptions about how the registry or the packages are hosted.

Registry creation is done by the function `create_registry`.

Registration of new and updated packages is done by the function `register`.
"""
module LocalRegistry

using RegistryTools: RegistryTools, gitcmd, Compress,
                     check_and_update_registry_files, ReturnStatus, haserror,
                     find_registered_version
using UUIDs: uuid4
using Pkg: Pkg, TOML

export create_registry, register

# Note: The `uuid` keyword argument is intentionally omitted from the
# documentation string since it's not intended for users.
"""
    create_registry(name, repo)
    create_registry(path, repo)

Create a registry with the given `name` or at the local directory
`path`, and with repository URL `repo`. The first argument is
interpreted as a path if it has more than one path component and
otherwise as a name. If a path is given, the last path component is
used as the name of the registry. If a name is given, it is created in
the standard registry location. In both cases the registry path must
not previously exist.

*Keyword arguments*

    create_registry(...; description = nothing, push = false, gitconfig = Dict())

* `description`: Optional description of the purpose of the registry.
* `push`: If `false`, the registry will only be prepared locally. Review the result and `git push` it manually.
* `gitconfig`: Optional configuration parameters for the `git` command.
"""
function create_registry(name_or_path, repo; description = nothing,
                         gitconfig::Dict = Dict(), uuid = nothing,
                         push = false)
    if length(splitpath(name_or_path)) > 1
        path = abspath(expanduser(name_or_path))
    else
        path = joinpath(first(DEPOT_PATH), "registries", name_or_path)
    end
    name = basename(path)
    if isempty(name)
        # If `path` ends with a slash, `basename` becomes empty and
        # the slash needs to be peeled off with `dirname`.
        name = basename(dirname(path))
    end
    if isdir(path) || isfile(path)
        error("$path already exists.")
    end
    mkpath(path)

    # The only reason for the `uuid` function argument is to allow
    # deterministic testing of the package.
    if isnothing(uuid)
        uuid = string(uuid4())
    end

    registry_data = RegistryTools.RegistryData(name, uuid, repo = repo,
                                               description = description)
    RegistryTools.write_registry(joinpath(path, "Registry.toml"), registry_data)

    git = gitcmd(path, gitconfig)
    run(`$git init -q`)
    run(`$git add Registry.toml`)
    run(`$git commit -qm 'Create registry.'`)
    run(`$git remote add origin $repo`)
    if push
        run(`$git push -u origin master`)
    end
    @info "Created registry in directory $(path)"

    return path
end

"""
    register(package, registry = registry)

Register the new `package` in the registry `registry`. The version
number and other information is obtained from the package's
`Project.toml`.

    register(package)

Register a new version of `package`.  The version number is obtained
from the package's `Project.toml`.

    register()

Register a new version of the package in the currently active project.

Notes:
 * By default this will, in all cases, `git push` the updated registry
   to its remote repository. If you prefer to do the push manually,
   use the keyword argument `push = false`.
 * The package must live in a git working copy, e.g. having been
   cloned by `Pkg.develop`.

`package` can be specified in the following ways:
* By package name. The package must be available in the active `Pkg`
  environment.
* By path. This is distinguished from package name by having more than
  one path component. A path in the current working directory can be
  specified by starting with `"./"`.
* By module. It needs to first be loaded by `using` or `import`.

`registry` can be specified by name or by path in the same way as
`package`. If omitted or `nothing`, it will be automatically located
by `package`.

*Keyword arguments*

    register(package; registry = nothing, commit = true, push = true,
             repo = nothing, gitconfig = Dict())

* `registry`: Name or path of registry.
* `commit`: If `false`, only make the changes to the registry but do not commit.
* `push`: If `true`, push the changes to the registry repository automatically. Ignored if `commit` is false.
* `repo`: Specify the package repository explicitly. Otherwise looked up as the `git remote` of the package the first time it is registered.
* `gitconfig`: Optional configuration parameters for the `git` command.
"""
function register(package::Union{Nothing, Module, AbstractString} = nothing;
                  registry::Union{Nothing, AbstractString} = nothing,
                  kwargs...)
    do_register(package, registry; kwargs...)
    return
end

# Differs from the above by looser type restrictions on `package` and
# `registry`. Also returns false if there was nothing new to register
# and true if something new was registered.
function do_register(package, registry;
                     commit = true, push = true, repo = nothing,
                     gitconfig::Dict = Dict())
    # Find and read the `Project.toml` for the package. First look for
    # the alternative `JuliaProject.toml`.
    package_path = find_package_path(package)
    local pkg
    for project_file in Base.project_names
        pkg = Pkg.Types.read_project(joinpath(package_path, project_file))
        if !isnothing(pkg.name)
            break
        end
    end
    if isnothing(pkg.name)
        error("$(package) does not have a Project.toml or JuliaProject.toml file")
    end

    # If the package directory is dirty, a different version could be
    # present in Project.toml.
    if is_dirty(package_path, gitconfig)
        error("Package directory is dirty. Stash or commit files.")
    end

    registry_path = find_registry_path(registry, pkg)
    if is_dirty(registry_path, gitconfig)
        error("Registry directory is dirty. Stash or commit files.")
    end

    # Compute the tree hash for the package and the subdirectory
    # location within the git repository. For normal packages living
    # at the top level of the repository, `subdir` will be the empty
    # string.
    tree_hash, subdir = get_tree_hash(package_path, gitconfig)

    # Check if this exact package has already been registered. Do
    # nothing and don't error in that case.
    if find_registered_version(pkg, registry_path) == tree_hash
        @info "This version has already been registered and is unchanged."
        return false
    end

    # Use the `repo` argument or, if this is a new package
    # registration, check for the git remote. If `repo` is `nothing`
    # and this is an existing package, the repository information will
    # not be updated.
    if isnothing(repo)
        if !has_package(registry_path, pkg)
            package_repo = get_remote_repo(package_path, gitconfig)
        else
            package_repo = ""
        end
    else
        package_repo = repo
    end

    @info "Registering package" package_path registry_path package_repo uuid=pkg.uuid version=pkg.version tree_hash subdir
    clean_registry = true

    git = gitcmd(registry_path, gitconfig)
    HEAD = readchomp(`$git rev-parse --verify HEAD`)
    status = ReturnStatus()
    try
        check_and_update_registry_files(pkg, package_repo, tree_hash,
                                        registry_path, String[], status,
                                        subdir = subdir)
        if !haserror(status)
            if commit
                commit_registry(pkg, package_path, package_repo, tree_hash, git)
                if push
                    run(`$git push`)
                end
            end
            clean_registry = false
        end
    finally
        if clean_registry
            run(`$git reset --hard $(HEAD)`)
            run(`$git clean -f -d`)
        end
    end

    # Registration failed. Explain to the user what happened.
    if haserror(status)
        error(explain_registration_error(status))
    end

    return true
end

function explain_registration_error(status)
    for triggered_check in status.triggered_checks
        check = triggered_check.id
        data = triggered_check
        if check == :version_exists
            return "Version $(data.version) has already been registered and the content has changed."
        elseif check == :change_package_name
            return "Changing package name is not supported."
        elseif check == :change_package_uuid
            return "Changing package UUID is not allowed."
        elseif check == :package_self_dep
            return "The package depends on itself."
        elseif check == :name_mismatch
            return "Error in (Julia)Project.toml: UUID $(data.uuid) refers to package '$(data.reg_name)' in registry but Project.toml has '$(data.project_name)'."
        elseif check == :wrong_stdlib_uuid
            return "Error in (Julia)Project.toml: UUID $(data.project_uuid) for package $(data.name) should be $(data.stdlib_uuid)"
        elseif check == :package_url_missing
            return "No repo URL provided for a new package."
        elseif check == :unexpected_registration_error
            return "Unexpected registration error."
        end
    end
end

# This does the same thing as `LibGit2.isdirty(LibGit2.GitRepo(path))`
# but also works when `path` is a subdirectory of a git
# repository. Only dirt within the subdirectory is considered.
function is_dirty(path, gitconfig)
    git = gitcmd(path, gitconfig)
    # TODO: There should be no need for the `-u` option but without it
    # a bogus diff is reported in the tests.
    return !isempty(read(`$git diff-index -u HEAD -- .`))
end

# If the package is omitted, the active project must correspond to a
# package.
function find_package_path(::Nothing)
    path = ""
    if VERSION < v"1.4"
        env = Pkg.Types.EnvCache()
        if !isnothing(env.pkg)
            path = dirname(env.project_file)
        end
    else
        # Pkg.project() was introduced in Julia 1.4 as an experimental
        # feature. Effectively this does the same thing as the code
        # above but is hopefully more future safe.
        project = Pkg.project()
        if project.ispackage
            path = dirname(project.path)
        end
    end

    if path == ""
        error("The active project is not a package.")
    end

    return path
end

# If the package is provided as a module, directly find the package
# path from the loaded code. This works both if the module is loaded
# from the current package environment or found in LOAD_PATH.
function find_package_path(package::Module)
    return dirname(dirname(pathof(package)))
end

# A string argument is either interpreted as a path or as a package
# name, decided by the number of components returned by `splitpath`.
#
# If the package is given by name, it must be available in the current
# package environment as a developed package.
function find_package_path(package_name::AbstractString)
    if length(splitpath(package_name)) > 1
        if !isdir(package_name)
            error("Package path $(package_name) does not exist.")
        end
        return abspath(expanduser(package_name))
    end

    ctx = Pkg.Types.Context()
    if !haskey(ctx.env.project.deps, package_name)
        error("Unknown package $package_name.")
    end
    pkg_uuid = ctx.env.project.deps[package_name]
    pkg_path = ctx.env.manifest[pkg_uuid].path
    if isnothing(pkg_path)
        error("Package must be developed to be registered.")
    elseif !isabspath(pkg_path)
        # If the package is developed with --local, pkg_path is
        # relative to the project path.
        pkg_path = joinpath(dirname(ctx.env.manifest_file), pkg_path)
    end

    # `pkg_path` might be a relative path, in which case it is
    # relative to the directory of Manifest.toml. If `pkg_path`
    # already is an absolute path, this call does not affect it.
    pkg_path = abspath(dirname(ctx.env.manifest_file), pkg_path)

    return pkg_path
end

function find_registry_path(registry::AbstractString, ::Pkg.Types.Project)
    return find_registry_path(registry)
end

function find_registry_path(registry::AbstractString)
    if length(splitpath(registry)) > 1
        return abspath(expanduser(registry))
    end

    all_registries = collect_registries()
    matching_registries = filter(r -> r.name == registry, all_registries)
    if isempty(matching_registries)
        error("Registry $(registry) not found.")
    end

    return first(matching_registries).path
end

function find_registry_path(::Nothing, pkg::Pkg.Types.Project)
    all_registries = collect_registries()
    all_registries_but_general = filter(r -> r.name != "General",
                                        all_registries)

    matching_registries = filter(all_registries) do reg_spec
        reg_data = Pkg.TOML.parsefile(joinpath(reg_spec.path, "Registry.toml"))
        haskey(reg_data["packages"], string(pkg.uuid))
    end

    if isempty(matching_registries)
        if length(all_registries_but_general) == 1
            return first(all_registries_but_general).path
        else
            error("Package $(pkg.name) not found in any registry. Please specify in which registry you want to register it.")
        end
    elseif length(matching_registries) > 1
        error("Package $(pkg.name) is registered in more than one registry, please specify in which you want to register the new version.")
    end

    return first(matching_registries).path
end

# This replaces the use of `Pkg.Types.collect_registries` which was
# removed in Julia 1.7.
#
# TODO: Once Julia versions before 1.7 are no longer supported,
# consider switching over to use `Pkg.Registry.reachable_registries`
# where this is called.
function collect_registries()
    registries = []
    for depot in Pkg.depots()
        isdir(depot) || continue
        reg_dir = joinpath(depot, "registries")
        isdir(reg_dir) || continue
        for name in readdir(reg_dir)
            file = joinpath(reg_dir, name, "Registry.toml")
            isfile(file) || continue
            push!(registries, (name = name, path = joinpath(reg_dir, name)))
        end
    end
    return registries
end

function has_package(registry_path, pkg::Pkg.Types.Project)
    registry = Pkg.TOML.parsefile(joinpath(registry_path, "Registry.toml"))
    return haskey(registry["packages"], string(pkg.uuid))
end

function get_tree_hash(package_path, gitconfig)
    git = gitcmd(package_path, gitconfig)
    subdir = readchomp(`$git rev-parse --show-prefix`)
    tree_hash = readchomp(`$git rev-parse HEAD:$subdir`)
    # Get rid of trailing slash.
    if isempty(basename(subdir))
        subdir = dirname(subdir)
    end
    return tree_hash, subdir
end

function get_remote_repo(package_path, gitconfig)
    git = gitcmd(package_path, gitconfig)
    remote_names = split(readchomp(`$git remote`), '\n', keepempty=false)
    repos = String[]
    foreach(remote_names) do remote_name
        push!(repos,readchomp(`$git remote get-url $(remote_name)`))
    end
    length(repos) === 0 && error("No repo URL found. Try calling `register` with the keyword `repo` to provide a URL.")
    length(repos) > 1 && error("Multiple repo URLs found. Try calling `register` with the keyword `repo` to provide a URL.\n$(repos)")
    return repos[1]
end

function commit_registry(pkg::Pkg.Types.Project, package_path, package_repo, tree_hash, git)
    @debug("commit changes")
    message = """
    New version: $(pkg.name) v$(pkg.version)

    UUID: $(pkg.uuid)
    Repo: $(package_repo)
    Tree: $(string(tree_hash))
    """
    run(`$git add --all`)
    run(`$git commit -qm $message`)
end

end


text:
#operators constructors

struct Op
    dim::Int
    cmtot::SparseMatrixCSC{Float64,Int64}
    cdtot::SparseMatrixCSC{Float64,Int64}
    le::Int
    basis::SparseMatrixCSC{Float64,Int64}
    Op(dim) = new(dim, operators(dim)...)
end

# after being defined, one can
# retrieve the dimension by typping
# dim(o)
dim(o::Op) = o.dim
cmtot(o::Op) = o.cmtot
cdtot(o::Op) = o.cdtot
le(o::Op) = o.le
basis(o::Op) = o.basis

#These are the fermionic operators we wanted
#By calling cm(op, 1) we get the sparse matrix
#corresponding to the destruction of the
#first fermionic mode
a(o::Op, i::Int) = cmtot(o)[1:le(o),((i-1)*le(o)+1):i*le(o)]
ad(o::Op, i::Int) = cdtot(o)[((i-1)*le(o)+1):i*le(o),1:le(o)]
ada(o::Op, i::Int, j::Int) = ad(o,i)*a(o,j)
aad(o::Op, i::Int, j::Int) = a(o,i)*ad(o,j)
aa(o::Op, i::Int, j::Int) = a(o,i)*a(o,j)
adad(o::Op, i::Int, j::Int) = ad(o,i)*ad(o,j)

function vacuum(o::Op)
    l = 2^dim(o)
    vac = spzeros(l)
    vac[1] = 1.0
    return vac
end


text:
using DataFrames
using TimeSeries
using Base.Dates
using Blpapi
using Gadfly
using GLM
using JuliaInXL

type connObject

  IP
  Port
  session

end

function toXLDate(d::Date)
    return Int(d - Date(1899,12,30))
end

function createConn(IP = "localhost", Port = 8194)
  try
    session = createSession(IP, Port)
    conn = connObject(IP, Port, session)
  catch e
    e_msg = string(sprint(io->Base.showerror(io,e)), sprint(io->Base.show_backtrace(io,catch_backtrace())))
    println(e_msg)
    return e_msg
  end
end

type instrument

  ticker
  fields
  startDate
  endDate

end

type basicInstrument

    ticker
    fieldnames

end

function arrayAusgabe(n)

  A = rand(Float64, Int(n), Int(n))
  return A

end

function TS(con = connection, inst = ibm)

  response = bdh(con.session, [inst.ticker], [inst.fields], inst.startDate, inst.endDate; periodicitySelection = "DAILY")
  data = response[inst.ticker, inst.fields]

  dformat = Dates.DateFormat("y-m-d")
  dates = similar(data[:,1], Date)
  prices = similar(data[:,2], Float64)

  for i in (1:length(data[:,1]))
    dates[i] = Date(data[i,1], dformat)
    prices[i] = data[i,2]
  end

  k = sortperm(dates)
  dates = dates[k]
  prices = prices[k]

  data = TimeArray(dates, prices, [inst.ticker])

end

function BarFrame(con, inst, barsize)

  Response = bar(con.session, inst.ticker, inst.fields, string(inst.startDate), string(inst.endDate); interval=barsize)

  df = DataFrame()
  df[:Dates] = Response.datetimes
  df[:Open] = zeros(length(Response.datetimes))
  df[:High] = zeros(length(Response.datetimes))
  df[:Low] = zeros(length(Response.datetimes))
  df[:Close] = zeros(length(Response.datetimes))
  df[:Volume] = zeros(length(Response.datetimes))
  df[:RotFact] = zeros(length(Response.datetimes))

  for i in 1:length(Response.datetimes)
    df[i,:Open] = Response[Response.datetimes[i]].openVar
    df[i,:High] = Response[Response.datetimes[i]].highVar
    df[i,:Low] = Response[Response.datetimes[i]].lowVar
    df[i,:Close] = Response[Response.datetimes[i]].closeVar
    df[i,:Volume] = Response[Response.datetimes[i]].volumeVar
  end

  for i in collect(nrow(df):-1:2)
      if df[i,:High] > df[i-1,:High] && df[i,:Low] > df[i-1,:Low]
          df[i,:RotFact] = mean([df[i,:High], df[i,:Low]])*df[i,:Volume]
      elseif df[i,:High] < df[i-1,:High] && df[i,:Low] < df[i-1,:Low]
          df[i,:RotFact] = mean([df[i,:High], df[i,:Low]])*df[i,:Volume]*-1
      else
      df[i,:RotFact] = 0
      end
  end

  df

end

function intrRotFact(commodity, days_back)

    df = bbgTS.VolAnalytics[commodity]
    df = Array{Any}(df[Date(df[:,:Dates]) .== Date(now())-Day(days_back), [:Dates, :Volume, :RotFact]])
    df[:,1] = toXLDate.(df[:,1])
    return df

end

function intrRotFactXL(commodity, days_back)

    inst = bbgTS.contracts[commodity]
    inst.startDate = now() - Dates.Month(2)
    inst.endDate = now()
    df = BarFrame(bbgTS.conInternal, inst, bbgTS.n)
    df = Array{Any}(df[Date(df[:,:Dates]) .== Date(now())-Dates.Day(days_back), [:Dates, :Volume, :RotFact]])
    df[:,1] = toXLDate.(df[:,1])
    return df

end

# t = now()
# startDate = t - Dates.Month(2)
# conInternal = createConn()
# n = 30

# VolAnalytics = Dict([("Gas", BarFrame(conInternal, instrument("FNA Comdty", "TRADE", startDate, t), n)),
#                 ("Brent", BarFrame(conInternal, instrument("COA Comdty", "TRADE", startDate, t), n)),
#                 ("Bund", BarFrame(conInternal, instrument("RXA Comdty", "TRADE", startDate, t), n))])

# contracts = Dict([("Gas", instrument("FNA Comdty", "TRADE", startDate, t)),
#                     ("Brent", instrument("COA Comdty", "TRADE", startDate, t)),
#                     ("Bund", instrument("RXA Comdty", "TRADE", startDate, t))])
#
#
#
#
#
# instContracts = Dict([("NBP", basicInstrument("FNA Comdty", "TRADE")),
#                         ("Brent", basicInstrument("COA Comdty", "TRADE")),
#                         ("Bund", basicInstrument("RXA Comdty", "TRADE")),
#                         ("Sugar", basicInstrument("QWA Comdty", "TRADE")),
#                         ("Cocoa", basicInstrument("QCA Comdty", "TRADE"))])


text:
function xyz_from_zmat(zmat_string::String)
    # Match an atom id (String for atomic symbol or integer for atomic number) with one or more blank spaces on the left of it
    re_id = r"\s*(\d+|\w{1,2})"

    # Match an integer (mapping to an atom) followed by a symbol or value (bond length, angle or dihedral). e.g. "1 R" or "2 106.5"
    re_map = r"\s*(\d+)\s+?([^\s]+)"

    # Match a definition e.g. "R = 1.024"
    re_def = r"(\w+)\s*[=]\s*(\d+)"

    line_num = 0
    entries = []
    definitions = []
    for line in split(strip(zmat_string), "\n")

        line_num += 1

        if isempty(strip(line, [' ']))
            continue
        end

        # Required format for first line:  "ATOM_ID"
        if line_num == 1
            m = occursin(re_id, line) ? match(re_id, line) : throw(ArgumentError("Failed to process first ZMAT entry:\n $(line)"))
            push!(entries, m.captures...)

        # Required format for second line:  "ATOM_ID INT VALUE"
        # e.g. "H 1 R" "1 1 1.025"
        elseif line_num == 2
            re_2line = re_id*re_map
            m = occursin(re_2line, line) ? match(re_2line, line) : throw(ArgumentError("Failed to process second ZMAT entry:\n $(line)"))
            println(m)
            push!(entries, m.captures...)

        # Required format for third line:  "ATOM_ID INT VALUE INT VALUE"
        # e.g. "H 1 R 2 A" "1 1 1.025 2 120"
        elseif line_num == 3
            re_3line = re_id*re_map*re_map
            m = occursin(re_3line, line) ? match(re_3line, line) : throw(ArgumentError("Failed to process third ZMAT entry:\n $(line)"))
            push!(entries, m.captures...)

        # Remaining lines are expected to be "ATOM_ID INT VALUE INT VALUE INT VALUE" or "VAR = VALUE"
        else
            re_line = re_id*re_map*re_map*re_map
            if occursin(re_def, line)
                m = match(re_def, line)
                push!(definitions, m.captures...)
            elseif occursin(re_line, line)
                m = match(re_line, line)
                push!(entries, m.captures...)
            else
                throw(ArgumentError("Failed to process line $line_num on ZMAT:\n $(line)"))
            end
        end
    end

    line_format = FormatExpr("{:<3}   {}")
    out = ""
    return entries, definitions
end


text:
__precompile__()

module Porta

# The default tolerance for comparing the equality of two values.
const TOLERANCE = 1e-7

logicdir = "logic"
linearalgebradir = "linearalgebra"
geometrydir = "geometry"
bundlesdir = "bundles"
uidir = "ui"

include(joinpath(logicdir, "propositionallogic.jl"))
include(joinpath(logicdir, "predicatelogic.jl"))

include(joinpath(linearalgebradir, "abstractvectorspace.jl"))
include(joinpath(linearalgebradir, "determinant.jl"))
include(joinpath(linearalgebradir, "real3.jl"))
include(joinpath(linearalgebradir, "real4.jl"))

include(joinpath(geometrydir, "s1.jl"))
include(joinpath(geometrydir, "s2.jl"))
include(joinpath(geometrydir, "s3.jl"))
include(joinpath(geometrydir, "biquaternions.jl"))
include(joinpath(geometrydir, "stereographicprojection.jl"))
include(joinpath(geometrydir, "rotations.jl"))
include(joinpath(geometrydir, "body.jl"))

include(joinpath(bundlesdir, "clifford.jl"))

include(joinpath(uidir, "ui.jl"))
include(joinpath(uidir, "arrow.jl"))
include(joinpath(uidir, "cylinder.jl"))
include(joinpath(uidir, "sphere.jl"))
include(joinpath(uidir, "torus.jl"))
include(joinpath(uidir, "triad.jl"))
include(joinpath(uidir, "whirl.jl"))
include(joinpath(uidir, "frame.jl"))

end # module


text:
### Test error exits
f = Yield ~ 1 + (1 | Batch)
mf = ModelFrame(f, ds)
Λ = [LowerTriangular(eye(1)) for i = 1]
y = convert(Vector{Float64}, ds[:Yield])
trms = push!(Any[], remat(:(1|Batch), ds), ones(30, 1), reshape(y, (length(y), 1)))

@test_throws DomainError LinearMixedModel(f, mf, trms, Λ, fill(-1., 30))

@test_throws ArgumentError lmm(Yield ~ 1, ds)

modl = lmm(Yield ~ 1 + (1|Batch), ds);

@test_throws DimensionMismatch MixedModels.getθ!(Array(Float64, (2,)), modl)
@test_throws DimensionMismatch setθ!(modl, [0., 1.])

@test_throws DimensionMismatch LinAlg.Ac_ldiv_B!(UpperTriangular(fm3.R[1, 1]), ones(30, 1))

const tri = LowerTriangular(eye(3))
const hblk = MixedModels.HBlkDiag(ones(2,3,2))
const dd3 = Diagonal(ones(3))

@test_throws DimensionMismatch setθ!(tri, ones(2))
@test_throws DimensionMismatch MixedModels.tscale!(tri, ones(2))
@test_throws DimensionMismatch MixedModels.tscale!(ones(2,2), tri)
@test_throws DimensionMismatch MixedModels.tscale!(tri, hblk)
@test_throws ArgumentError Base.cholfact!(hblk)
@test_throws DimensionMismatch MixedModels.downdate!(dd3, Diagonal(ones(2)))
@test_throws DimensionMismatch MixedModels.downdate!(dd3, sparse(ones(4,1)))
@test_throws DimensionMismatch MixedModels.downdate!(ones(2,2), dd3, ones(3,3))
@test_throws DimensionMismatch MixedModels.downdate!(ones(2,2), dd3, ones(2,2))
@test_throws DimensionMismatch MixedModels.inject!(ones(2,2), dd3)
@test_throws DimensionMismatch MixedModels.inject!(ones(3,2), dd3)
@test_throws DimensionMismatch LinAlg.A_ldiv_B!(dd3, Diagonal(ones(2)))

const speye4 = speye(4)
@test_throws ArgumentError MixedModels.inject!(speye(4), sparse(ones(4,4)))
@test_throws DimensionMismatch LinAlg.A_ldiv_B!(dd3, speye4)
@test_throws DimensionMismatch LinAlg.Ac_ldiv_B!(UpperTriangular(hblk), ones(1,1))
@test_throws DimensionMismatch MixedModels.unscaledre!(ones(1), modl.wttrms[1], ones(1, 6))

@test_throws DimensionMismatch MixedModels.downdate!(ones(1,1), speye4, ones(4,1))
@test_throws DimensionMismatch MixedModels.downdate!(ones(1,1), speye4, speye4)
@test_throws DimensionMismatch MixedModels.downdate!(ones(1,1), speye(3), speye4)
@test_throws DimensionMismatch MixedModels.downdate!(ones(1,1), speye4)

@test_throws DimensionMismatch MixedModels.inject!(speye4, speye(3))

@test_throws DimensionMismatch MixedModels.unscaledre!(zeros(size(slp, 1)),
  remat(:(1 + Days | Subject), slp), ones(1,30))

@test_throws DimensionMismatch MixedModels.tscale!(LowerTriangular(eye(2)), Diagonal(ones(30)))
@test_throws DimensionMismatch MixedModels.tscale!(Diagonal(ones(30)), LowerTriangular(eye(2)))
@test_throws ArgumentError MixedModels.lrt(modl)


text:
valid_moves(piece::Bishop) = _valid_moves_diagonal(piece)

valid_moves(piece::Queen) = vcat(
  _valid_moves_straight(piece), _valid_moves_diagonal(piece)
)

function valid_moves(piece::Pawn)
  is_white = piece.player.is_white
  other_player = filter(tmp_player -> tmp_player != piece.player, piece.player.game.players)[1]

  if is_white
    @assert piece.row > 1
    direction = +1
  else
    @assert piece.row < 8
    direction = -1
  end

  @assert 1 <= piece.row+direction <= 8

  board = piece.player.game.board
  moves = []

  work_i, work_j = piece.row+direction, piece.col
  if isnothing(piece.player.game.board[work_i,work_j])
    push!(moves, [work_i, work_j])

    is_init = (
      ( is_white && piece.row == 2 ) ||
      ( !is_white && piece.row == 7 )
    )

    if is_init
      work_i, work_j = piece.row+2*direction, piece.col
      isnothing(piece.player.game.board[work_i,work_j]) && push!(moves, [work_i, work_j])
    end
  end

  work_i = piece.row+direction
  for work_j in piece.col .+ [-1,+1]
    ( 1 <= work_j <= 8 ) || continue

    anti_pawn = other_player.anti_pawn
    work_piece = board[work_i,work_j]

    can_move = !isnothing(anti_pawn) && anti_pawn.row == work_i && anti_pawn.col == work_j
    can_move |= !isnothing(work_piece) && work_piece.player != piece.player

    can_move && push!(moves, [work_i, work_j])
  end

  moves
end

function valid_moves(piece::Knight)
  board = piece.player.game.board
  moves = []

  for work_i in piece.row .+ [+1,-1]
    ( 1 <= work_i <= 8 ) || continue
    for work_j in piece.col .+ [+2,-2]
      ( 1 <= work_j <= 8 ) || continue
      is_valid = isnothing(piece.player.game.board[work_i,work_j])
      is_valid || ( is_valid = ( board[work_i,work_j].player != piece.player ) )
      is_valid && push!(moves, [work_i, work_j])
    end
  end

  for work_i in piece.row .+ [+2,-2]
    ( 1 <= work_i <= 8 ) || continue
    for work_j in piece.col .+ [+1,-1]
      ( 1 <= work_j <= 8 ) || continue
      is_valid = isnothing(piece.player.game.board[work_i,work_j])
      is_valid || ( is_valid = ( board[work_i,work_j].player != piece.player ) )
      is_valid && push!(moves, [work_i, work_j])
    end
  end

  moves
end

function valid_moves(piece::Rook)
  moves = _valid_moves_straight(piece)
  ( piece.player.can_castle_left || piece.player.can_castle_right ) || return moves
  piece.player.has_castled && return moves

  if piece.player.is_white
    work_i = 1
  else
    work_i = 8
  end

  king = piece.player.pieces[findfirst(work_piece -> isa(work_piece,King), piece.player.pieces)]

  ( piece.row == work_i ) || return moves
  ( king.row == work_i ) || return moves

  @assert king.col == 5

  if piece.col < king.col
    j_list = piece.col+1:king.col-1
  else
    j_list = king.col+1:piece.col-1
  end

  for work_j in j_list
    isnothing(piece.player.game.board[work_i,work_j]) || return moves
  end

  push!(moves, [king.row,king.col])

  moves
end

function valid_moves(piece::King)
  board = piece.player.game.board
  moves = []

  for work_i in piece.row .+ (-1:+1)
    ( 1 <= work_i <= 8 ) || continue
    for work_j in piece.col .+ (-1:+1)
      ( 1 <= work_j <= 8 ) || continue

      if !isnothing(piece.player.game.board[work_i,work_j]) && board[work_i,work_j].player == piece.player
        continue
      end

      push!(moves, [work_i, work_j])
    end
  end

  moves
end

function _valid_moves_diagonal(piece::Piece)
  board = piece.player.game.board
  moves = []

  for i_dir = [-1,+1]
    for j_dir = [-1,+1]
      work_i = piece.row + i_dir
      work_j = piece.col + j_dir

      while 1 <= work_i <= 8 && 1 <= work_j <= 8
        if isnothing(piece.player.game.board[work_i,work_j])
          push!(moves, [work_i, work_j])
        elseif board[work_i,work_j].player == piece.player
          break
        else
          push!(moves, [work_i, work_j])
          break
        end

        work_i += i_dir
        work_j += j_dir
      end
    end
  end

  moves
end

function _valid_moves_straight(piece::Piece)
  board = piece.player.game.board
  moves = []

  work_j = piece.col

  for work_i in piece.row-1:-1:1
    push!(moves, [work_i, work_j])
    isnothing(board[work_i,work_j]) && continue
    ( board[work_i,work_j].player == piece.player ) || break
    pop!(moves) ; break
  end

  for work_i in piece.row+1:+1:8
    push!(moves, [work_i, work_j])
    isnothing(board[work_i,work_j]) && continue
    ( board[work_i,work_j].player == piece.player ) || break
    pop!(moves) ; break
  end

  work_i = piece.row

  for work_j in piece.col-1:-1:1
    push!(moves, [work_i, work_j])
    isnothing(board[work_i,work_j]) && continue
    ( board[work_i,work_j].player == piece.player ) || break
    pop!(moves) ; break
  end

  for work_j in piece.col+1:+1:8
    push!(moves, [work_i, work_j])
    isnothing(board[work_i,work_j]) && continue
    ( board[work_i,work_j].player == piece.player ) || break
    pop!(moves) ; break
  end

  moves
end


text:
ENV["PKG_CONFIG_PATH"] = "/Users/dc/anaconda/envs/python35/lib/pkgconfig"

using OpenCV
using Images, ImageView
using Cxx

import ImageView

# C++ OpenCV code
cxx"""

#include "opencv2/objdetect.hpp"
#include "opencv2/highgui.hpp"
#include "opencv2/imgproc.hpp"

#include <iostream>
#include <stdio.h>

using namespace std;
using namespace cv;

/// Load cascade file
CascadeClassifier load_face_cascade(String face_cascade_name) {

    CascadeClassifier face_cascade;
    face_cascade.load(face_cascade_name);

    return face_cascade;
}

/// detect and draw a bounding box on the image
void detect_face(cv::Mat frame, CascadeClassifier face_cascade) {

    std::vector<Rect> faces;
    cv::Mat frame_gray;
    std::vector<int>vec(4);

    cvtColor( frame, frame_gray, COLOR_BGR2GRAY );
    equalizeHist( frame_gray, frame_gray );

    face_cascade.detectMultiScale( frame_gray, faces, 1.1, 2, 0|CASCADE_SCALE_IMAGE, Size(30, 30) );

    for( size_t i = 0; i < faces.size(); i++ ) {
        Point center( faces[i].x + faces[i].width/2, faces[i].y + faces[i].height/2 );
        ellipse( frame, center, Size( faces[i].width/2, faces[i].height/2), 0, 0, 360, Scalar( 255, 0, 255 ), 4, 8, 0 );
    }
}

/// detect and return bounding box coordinates
std::vector<int> detect_face_coords(cv::Mat frame, CascadeClassifier face_cascade) {

    std::vector<Rect> faces;
    cv::Mat frame_gray;
    std::vector<int>vec(4);

    cvtColor( frame, frame_gray, COLOR_BGR2GRAY );
    equalizeHist( frame_gray, frame_gray );

    face_cascade.detectMultiScale( frame_gray, faces, 1.1, 2, 0|CASCADE_SCALE_IMAGE, Size(30, 30) );

    if (faces.size() > 0) {
        vec[0] = faces[0].x;
        vec[1] = faces[0].width;
        vec[2] = faces[0].y;
        vec[3] = faces[0].height;
    }
    else {
        vec[0] = 0; vec[1] = 0;
        vec[2] = 0; vec[3] = 0;
    }

    return vec;
}

"""

cascade_path = download("https://raw.githubusercontent.com/opencv/opencv/master/data/haarcascades/haarcascade_frontalface_alt.xml")
face_cascade = @cxx load_face_cascade(pointer(cascade_path));

# create a bounding box around face and save results in the root folder
filename = joinpath(pwd(), "sample-images", "beautiful-1274051_640_100_1.jpg");
img_opencv = imread(filename);
@cxx detect_face(img_opencv, face_cascade);
imwrite(joinpath(pwd(), "result.jpg"), img_opencv)

# retrieve coordinates and crop the image
filename = joinpath(pwd(), "sample-images", "beautiful-1274051_640_100_1.jpg");
img_opencv = imread(filename);
img_images = load(filename);

coords_cxx = @cxx detect_face_coords(img_opencv, face_cascade);
coords = map(x -> Int(at(coords_cxx, x)), 0:3);

face = img_images[coords[3]:coords[3] + coords[4], coords[1]:coords[1] + coords[2]]
ImageView.imshow(face)


text:
"""
    Index
"""
struct Index{name,K,V,Ks,Vs} <: AbstractIndex{name,K,V,Ks,Vs}
    _keys::Ks
    _values::Vs

    function Index{name,K,V,Ks,Vs}(ks::Ks, vs::Vs, uc::Uniqueness, lc::AbstractLengthCheck) where {name,K,V,Ks<:AbstractVector{K},Vs<:AbstractUnitRange{V}}
        check_index_length(ks, vs, lc)
        check_index_uniqueness(ks, uc)
        return new{name,K,V,Ks,Vs}(ks, vs)
    end
end

###
### Constructors
###
function Index{name,K,V,Ks,Vs}(
    ks::Ks2,
    vs::Vs2,
    uc::Uniqueness=UnkownUnique,
    lc::AbstractLengthCheck=LengthNotChecked
   ) where {name,K,V,Ks<:AbstractVector{K},Vs<:AbstractUnitRange{V},Ks2,Vs2}
    return Index{name}(Ks(ks), Vs(vs), uc, lc)
end

function Index(ks::AbstractVector, uc::Uniqueness=UnkownUnique)
    if is_static(ks)
        return Index(ks, OneToSRange(length(ks)), uc, LengthChecked)
    elseif is_fixed(ks)
        return Index(ks, OneTo(length(ks)), uc, LengthChecked)
    else
        return Index(ks, OneToMRange(length(ks)), uc, LengthChecked)
    end
end

function Index{name}(
    ks::AbstractVector,
    uc::Uniqueness=UnkownUnique
   ) where {name,K}
    if is_static(ks)
        return Index{name}(ks, OneToSRange(length(ks)), uc, LengthChecked)
    elseif is_fixed(ks)
        return Index{name}(ks, OneTo(length(ks)), uc, LengthChecked)
    else
        return Index{name}(ks, OneToMRange(length(ks)), uc, LengthChecked)
    end
end

function Index(
    ks::AbstractVector,
    vs::AbstractUnitRange{Int},
    uc::Uniqueness=UnkownUnique,
    lc::AbstractLengthCheck=LengthNotChecked
   )
    return Index{nothing}(ks, vs, uc, lc)
end

function Index{name}(
    ks::AbstractVector{K},
    vs::AbstractUnitRange{Int},
    uc::Uniqueness=UnkownUnique,
    lc::AbstractLengthCheck=LengthNotChecked
   ) where {name,K}
    return Index{name,K}(ks, vs, uc, lc)
end

function Index{name,K}(
    ks::Ks,
    vs::AbstractUnitRange{V},
    uc::Uniqueness=UnkownUnique,
    lc::AbstractLengthCheck=LengthNotChecked
   ) where {name,K,V,Ks<:AbstractVector{K}}
    return Index{name,K,V,Ks}(ks, vs, uc, lc)
end

function Index{name,K,V,Ks}(
    ks::Ks,
    vs::AbstractUnitRange{V},
    uc::Uniqueness=UnkownUnique,
    lc::AbstractLengthCheck=LengthNotChecked
   ) where {name,K,V,Ks<:AbstractVector{K}}
    return Index{name,K,V,Ks,typeof(vs)}(ks, vs, uc, lc)
end

function Index(
    ks::AbstractIndex,
    vs::AbstractUnitRange,
    uc::Uniqueness=UnkownUnique,
    lc::AbstractLengthCheck=LengthNotChecked
   )
    return Index{indnames(ks)}(keys(ks), vs, uc, lc)
end

#function Index{name1,K,V,Ks,Vs}(idx::Index) where {name1,name2,K,V,Ks,Vs}
#    return Index{}()
#end

Index(idx::Index) = Index(keys(idx), values(idx), AllUnique, LengthChecked)

###
### Interface
###
Base.keys(idx::Index) = getfield(idx, :_keys)

Base.values(idx::Index) = getfield(idx, :_values)

function StaticRanges.similar_type(
    idx::Index{name},
    ks_type::Type=keys_type(idx),
    vs_type::Type=values_type(idx)
   ) where {name}
    return Index{name,eltype(ks_type),eltype(vs_type),ks_type,vs_type}
end

function Base.setproperty!(idx::Index{name,K,V,Ks,Vs}, p::Symbol, val) where {name,K,V,Ks,Vs}
    if is_dynamic(idx)
        if p === :keys
            if val isa Ks
                if length(val) == length(idx)
                    return setfield!(idx, :_keys, val)
                else
                    error("`val` must be the same length as the provided index.")
                end
            else
                return setproperty!(idx, p, convert(Ks, val))
            end
        elseif p === :values
            if val isa Vs
                if length(val) == length(idx)
                    return setfield!(idx, :_values, val)
                else
                    error("`val` must be the same length as the provided index.")
                end
            else
                return setproperty!(idx, p, convert(Vs, val))
            end
        else
            error("No property named $p.")
        end
    else
        error("The keys and values of an Index must be mutable to use `setproperty!`, got $(typeof(idx)).")
    end
end

# TODO: is this the best fall back for and AbstractIndex?
AbstractIndex(x) = Index(x)

function Base.show(io::IO, idx::Index)
    if isnothing(dimnames(idx))
        print(io, "Index($(keys(idx)) => $(values(idx)))")
    else
        print(io, "Index{$(dimnames(idx))}($(keys(idx)) => $(values(idx)))")
    end
end

###
### pop
###
function Base.pop!(a::AbstractIndex)
    can_set_last(a) || error("Cannot change size of index of type $(typeof(a)).")
    pop!(keys(a))
    return pop!(values(a))
end

function Base.popfirst!(a::AbstractIndex)
    can_set_first(a) || error("Cannot change size of index of type $(typeof(a)).")
    popfirst!(keys(a))
    return popfirst!(values(a))
end


text:
# functions to colorize numeric arrays or grayscale images with color maps

"""Colorize an array by applying a colormap from `ColorSchemes`

For available colormaps see [this overview](http://juliagraphics.github.io/ColorSchemes.jl/stable/assets/figures/colorschemes.png).
"""
function colorize(A::AbstractArray;
    contrast::Bool=false, # enhance contrast by cutting the 2% on each side, like the QGIS default
    cmap::Symbol=:viridis, # matplotlib default
    nodata=-9999.0,
    vmin=nothing, # costly to filter twice
    vmax=nothing,
    bg=colorant"black") # black, color of nodata

    novalues = (A .== nodata) .| isnan.(A)
    values = broadcast(~, novalues)

    # calculate vmin and vmax
    # taking into account completely empty arrays
    if vmin == nothing || vmax == nothing
        if all(novalues)
            # return black image
            # shortcut because cannot reduce over empty collection
            return fill(bg, size(A))
        else
            # do it here to prevent passing twice over the data
            # if both vmin and vmax are nothing
            vmin_calculated, vmax_calculated = if contrast
                quantile(A[values], [0.02, 0.98])
            else
                extrema(A[values])
            end
        end
    end
    if vmin == nothing
        vmin = vmin_calculated
    end
    if vmax == nothing
        vmax = vmax_calculated
    end

    # get the right colormap from ColorSchemes
    colormap = getfield(ColorSchemes, cmap)
    if (vmin ≈ vmax) || (vmin > vmax)
        return fill(bg, size(A))
    end
    sc_img = (A .- vmin) .* (1.0 / (vmax - vmin))
    sc_img[novalues] = 0.0  # set to background color later, just to avoid get NaN InexactError
    color(x) = get(colormap, x)
    cimg = color.(sc_img)
    cimg[novalues] = bg
    cimg
end

colorize(A::Union{BitMatrix, Matrix{Bool}}) = Gray.(A)


text:
SharedStrings() = SharedStrings(EzXML.ElementNode("sst"))

function SharedStrings(root::EzXML.Node)
    @assert EzXML.nodename(root) == "sst"

    unformatted_strings = Vector{String}()
    for el in EzXML.eachelement(root)
        @assert EzXML.nodename(el) == "si" "Unsupported node $(EzXML.nodename(el)) in sst table."
        push!(unformatted_strings, unformatted_text(el))
    end

    return SharedStrings(root, unformatted_strings)
end

"""
    unformatted_text(el::EzXML.Node) :: String

Helper function to gather unformatted text from Excel data files.
It looks at all childs of `el` for tag name `t` and returns
a join of all the strings found.
"""
function unformatted_text(el::EzXML.Node) :: String

    function gather_strings!(v::Vector{String}, e::EzXML.Node)
        if EzXML.nodename(e) == "t"
            push!(v, EzXML.nodecontent(e))
        end

        for ch in EzXML.eachelement(e)
            gather_strings!(v, ch)
        end

        nothing
    end

    v_string = Vector{String}()
    gather_strings!(v_string, el)
    return join(v_string)
end

"""
    sst_unformatted_string(wb, index) :: String

Looks for a string inside the Shared Strings Table (sst).
`index` starts at 0.
"""
@inline sst_unformatted_string(sst::SharedStrings, index::Int) = sst.unformatted_strings[index+1]
@inline sst_unformatted_string(wb::Workbook, index::Int) :: String = sst_unformatted_string(wb.sst, index)
@inline sst_unformatted_string(xl::XLSXFile, index::Int) :: String = sst_unformatted_string(xl.workbook, index)
@inline sst_unformatted_string(ws::Worksheet, index::Int) :: String = sst_unformatted_string(ws.package, index)
sst_unformatted_string(target::Union{Workbook, XLSXFile, Worksheet}, index_str::String) :: String = sst_unformatted_string(target, parse(Int, index_str))