text:
module TestComponents

using Mimi
using Test

import Mimi:
    compdefs, compdef, compkeys, has_comp, first_period,
    last_period, compmodule, compname, compinstance, dim_keys, dim_values

my_model = Model()

# Try running model with no components
@test length(compdefs(my_model)) == 0
@test length(my_model) == 0
@test_throws ErrorException run(my_model)

# Now add several components to the module
@defcomp testcomp1 begin
    var1 = Variable(index=[time])
    par1 = Parameter(index=[time])

    """
    Test docstring.
    """
    function run_timestep(p, v, d, t)
        v.var1[t] = p.par1[t]
    end
end

@defcomp testcomp2 begin
    var1 = Variable(index=[time])
    par1 = Parameter(index=[time])

    function run_timestep(p, v, d, t)
        v.var1[t] = p.par1[t]
    end
end

@defcomp testcomp3 begin
    var1 = Variable(index=[time])
    par1 = Parameter(index=[time])
    cbox = Variable(index=[time, 5])    # anonymous dimension

    function run_timestep(p, v, d, t)
        v.var1[t] = p.par1[t]
    end
end

# Start building up the model
set_dimension!(my_model, :time, 2015:5:2110)
add_comp!(my_model, testcomp1)

# Testing that you cannot add two components of the same name
@test_throws ErrorException add_comp!(my_model, testcomp1)

# Testing to catch adding component twice
@test_throws ErrorException add_comp!(my_model, testcomp1)

# Testing to catch if before or after does not exist
@test_throws ErrorException add_comp!(my_model, testcomp2, before=:testcomp3)

# N.B. Throws ArgumentError in v1.0, but ErrorException in 0.7!
@test_throws ArgumentError add_comp!(my_model, testcomp2, after=:testcomp3)

# Add more components to model
add_comp!(my_model, testcomp2)
add_comp!(my_model, testcomp3)

# Check addition of anonymous dimension
dvals = dim_values(my_model.md, Symbol(5))
dkeys = dim_keys(my_model.md, Symbol(5))
@test dvals == dkeys

comps = collect(compdefs(my_model))

# Test compdefs, compdef, compkeys, etc.
@test comps == collect(compdefs(my_model.md))
@test length(comps) == 3
@test compdef(my_model, :testcomp3).comp_id == comps[3].comp_id
@test_throws KeyError compdef(my_model, :testcomp4) #this component does not exist
@test [compkeys(my_model.md)...] == [:testcomp1, :testcomp2, :testcomp3]
@test has_comp(my_model.md, :testcomp1) == true
@test has_comp(my_model.md, :testcomp4) == false

@test compmodule(testcomp3) == Main.TestComponents
@test compname(testcomp3) == :testcomp3

@test length(my_model) == 3
add_comp!(my_model, testcomp3, :testcomp3_v2)
@test length(my_model) == 4

#------------------------------------------------------------------------------
#   Tests for component run periods when resetting the model's time dimension
#------------------------------------------------------------------------------

@defcomp testcomp1 begin
    var1 = Variable(index=[time])
    par1 = Parameter(index=[time])

    function run_timestep(p, v, d, t)
        v.var1[t] = p.par1[t]
    end
end

# 1. Test resetting the time dimension without explicit first/last values

cd = testcomp1
@test cd.first === nothing   # original component definition's first and last values are unset
@test cd.last === nothing

m = Model()
set_dimension!(m, :time, 2001:2005)
add_comp!(m, testcomp1, :C) # Don't set the first and last values here
cd = compdef(m.md, :C)      # Get the component definition in the model

# These tests are not valid in the composite world...
#@test cd.first === nothing   # First and last values should still be nothing because they were not explicitly set
#@test cd.last === nothing

set_param!(m, :C, :par1, zeros(5))
Mimi.build(m)               # Build the model
ci = compinstance(m, :C) # Get the component instance
@test ci.first == 2001      # The component instance's first and last values should match the model's index
@test ci.last == 2005

set_dimension!(m, :time, 2005:2020) # Reset the time dimension
cd = compdef(m.md, :C)       # Get the component definition in the model

# These tests are not valid in the composite world...
#@test cd.first === nothing   # First and last values should still be nothing
#@test cd.last === nothing

update_param!(m, :par1, zeros(16); update_timesteps=true)
Mimi.build(m)               # Build the model
ci = compinstance(m, :C) # Get the component instance
@test ci.first == 2005      # The component instance's first and last values should match the model's index
@test ci.last == 2020

# 2. Test resetting the time dimension with explicit first/last values

m = Model()
set_dimension!(m, :time, 2000:2100)

add_comp!(m, testcomp1, :C)

cd = compdef(m.md, :C)      # Get the component definition in the model

# first and last are disabled currently
# @test cd.first == 2010      # First and last values are defined in the comp def because they were explicitly given
# @test cd.last == 2090

# Verify that they didn't change
#@test cd.first === nothing
#@test cd.last === nothing

set_dimension!(m, :time, 2010:2090)

set_param!(m, :C, :par1, zeros(81))
Mimi.build(m)               # Build the model
ci = compinstance(m, :C) # Get the component instance
@test ci.first == 2010      # The component instance's first and last values are the same as in the comp def
@test ci.last == 2090

set_dimension!(m, :time, 2000:2200) # Reset the time dimension
cd = compdef(m.md, :C)      # Get the component definition in the model
# @test cd.first == 2010      # First and last values should still be the same
# @test cd.last == 2090

Mimi.build(m)               # Build the model
ci = compinstance(m, :C) # Get the component instance
# @test ci.first == 2010      # The component instance's first and last values are the same as the comp def
# @test ci.last == 2090

end #module


text:
using DataStructures

function sumby_accumulator{T,S}(by::AbstractVector{T}, val::AbstractVector{S})
  a = Accumulator(T,S)
  for (byi, vali) = zip(by, val)
    push!(a, byi, vali)
  end
  return a
end

@time sumby_accumulator(id6, v1) # 100 seconds
#@time sumby_dict(id6, v1)

ext: jl


text:
# this user supplied file creates a dictionary of arguments
# if a key is repeated, the last use of the key is used
# it is a little bit dangerous letting the user run arbitrary code
# as part of the solver
# now that this file is read inside a function, it is better encapsulated
arg_dict = Dict{Any, Any}(
"physics" => "Euler",
"run_type" => 1,
"jac_type" => 3,
"newton_globalize_euler" => true,
"order" => 1,
"real_time" => false,
"variable_type" => :conservative,
"IC_name" => "ICIsentropicVortex",
"numBC" => 4,
"BC1" => [0],
"BC1_name" => "isentropicVortexBC",
"BC2" => [1],
"BC2_name" => "isentropicVortexBC",
"BC3" => [2],
"BC3_name" => "isentropicVortexBC",
"BC4" => [3],
"BC4_name" => "isentropicVortexBC",
"smb_name" => "src/mesh_files/gvortex1np2.smb",
"dmg_name" => ".null",
"use_DG" => true,
"operator_type" => "SBPOmega",
"Flux_name" => "RoeFlux",
"t_max" => 5.0,
"res_tol" => 1e-10,
"itermax" => 20,
"res_abstol" => 1e-10,
"res_reltol" => 1e-9,
"step_tol" => 1e-10,
"itermax" => 30,
"calc_functional" => true,
"functional_error" => true,
"num_functionals" => 1,
"functional_name1" => "drag",
"geom_edges_functional1" => [3],
"analytical_functional_val" => -1/1.4,
"calc_adjoint" => false,
"writeq" => false,
"write_edge_vertnums" => false,
"write_face_vertnums" => false,
"write_qic" => false,
"writeboundary" => false,
"write_res" => false,
"print_cond" => false,
"write_counts" => false,
"write_vis" => false,
"output_freq" => 1,
"solve" => true,
"do_postproc" => true,
"exact_soln_func" => "ICIsentropicVortex"
)


text:
include("utils.jl")
include("RBMBase.jl")
include("SamplingGibbs.jl")
include("SamplingEMF.jl")
include("Scoring.jl")
include("Monitor.jl")
include("Training.jl")
include("nets.jl")
include("charts.jl")
include("io.jl")


text:
## A SortedDict is a wrapper around balancedTree with
## methods similiar to those of Julia container Dict.

type SortedDict{K, D, Ord <: Ordering} <: Associative{K,D}
    bt::BalancedTree23{K,D,Ord}

## Zero-argument constructor, or possibly one argument to specify order.

    function SortedDict(o::Ord=Forward)
        bt1 = BalancedTree23{K,D,Ord}(o)
        new(bt1)
    end

end

## external constructor to take an associative and infer
## argument types

function SortedDict{K, D, Ord <: Ordering}(d::Associative{K,D}, o::Ord=Forward)
    h = SortedDict{K,D,Ord}(o)
    for (k,v) in d
        h[k] = v
    end
    h
end

if VERSION >= v"0.4.0-dev"

    ## More constructors based on those in dict.jl:
    ## Take pairs and infer argument
    ## types.  Note:  this works only for the Forward ordering.

    function SortedDict{K,D}(ps::Pair{K,D}...)
        h = SortedDict{K,D,ForwardOrdering}()
        for p in ps
            h[p.first] = p.second
        end
        h
    end

    ## Take pairs and infer argument
    ## types.  Ordering parameter must be explicit first argument.

    function SortedDict{K,D, Ord <: Ordering}(o::Ord, ps::Pair{K,D}...)
        h = SortedDict{K,D,Ord}(o)
        for p in ps
            h[p.first] = p.second
        end
        h
    end

    ## This one takes an iterable; ordering type is optional.

    SortedDict{Ord <: Ordering}(kv, o::Ord=Forward) =
    sorteddict_with_eltype(kv, eltype(kv), o)

    function sorteddict_with_eltype{K,D,Ord}(kv, ::Type{Tuple{K,D}}, o::Ord)
        h = SortedDict{K,D,Ord}(o)
        for (k,v) in kv
            h[k] = v
        end
        h
    end

else   #if VERSION < v"0.4.0-dev"
    function SortedDict{K,D,Ord <: Ordering}(ks::AbstractArray{K},
                                             vs::AbstractArray{D},
                                             o::Ord = Forward)
        h = SortedDict{K,D,Ord}(o)
        l = length(ks)
        if length(vs) != l
            error("ks and vs arrays in two-array SortedDict constructor must have the same length")
        end
        for i = 1 : l
            h[ks[i]] = vs[i]
        end
        h
    end

    function SortedDict{K,D,Ord <: Ordering}(kv::AbstractArray{(K,D),1},
                                             o::Ord = Forward)
        h = SortedDict{K,D,Ord}(o)
        for pr in kv
            h[pr[1]] = pr[2]
        end
        h
    end
end

typealias SDSemiToken IntSemiToken

typealias SDToken @compat Tuple{SortedDict,IntSemiToken}

## This function implements m[k]; it returns the
## data item associated with key k.

@inline function getindex(m::SortedDict, k_)
    i, exactfound = findkey(m.bt, convert(keytype(m),k_))
    !exactfound && throw(KeyError(k_))
    return m.bt.data[i].d
end

## This function implements m[k]=d; it sets the
## data item associated with key k equal to d.

@inline function setindex!{K, D, Ord <: Ordering}(m::SortedDict{K,D,Ord}, d_, k_)
    insert!(m.bt, convert(K,k_), convert(D,d_), false)
    m
end

## This function looks up a key in the tree;
## if not found, then it returns a marker for the
## end of the tree.

@inline function find(m::SortedDict, k_)
    ll, exactfound = findkey(m.bt, convert(keytype(m),k_))
    IntSemiToken(exactfound? ll : 2)
end

## This function inserts an item into the tree.
## Unlike m[k]=d, it also returns a bool and a token.
## The bool is true if the inserted item is new.
## It is false if there was already an item
## with that key.
## The token points to the newly inserted item.

@inline function insert!{K,D, Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, d_)
    b, i = insert!(m.bt, convert(K,k_), convert(D,d_), false)
    b, IntSemiToken(i)
end

## First and last return the first and last (key,data) pairs
## in the SortedDict.  It is an error to invoke them on an
## empty SortedDict.

@inline function first(m::SortedDict)
    i = beginloc(m.bt)
    i == 2 && throw(BoundsError())
    return m.bt.data[i].k, m.bt.data[i].d
end

@inline function last(m::SortedDict)
    i = endloc(m.bt)
    i == 1 && throw(BoundsError())
    return m.bt.data[i].k, m.bt.data[i].d
end

@inline function in{K,D,Ord <: Ordering}(pr::(@compat Tuple{Any,Any}), m::SortedDict{K,D,Ord})
    i, exactfound = findkey(m.bt,convert(K,pr[1]))
    return exactfound && isequal(m.bt.data[i].d,convert(D,pr[2]))
end

@inline eltype{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}) = @compat Tuple{K,D}
@inline eltype{K,D,Ord <: Ordering}(::Type{SortedDict{K,D,Ord}}) = @compat Tuple{K,D}
@inline orderobject(m::SortedDict) = m.bt.ord
@inline keytype{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}) = K
@inline datatype{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}) = D

@inline function haskey(m::SortedDict, k_)
    i, exactfound = findkey(m.bt,convert(keytype(m),k_))
    exactfound
end

@inline function get{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, default_)
    i, exactfound = findkey(m.bt, convert(K,k_))
   return exactfound? m.bt.data[i].d : convert(D,default_)
end

function get!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, default_)
    k = convert(K,k_)
    i, exactfound = findkey(m.bt, k)
    if exactfound
        return m.bt.data[i].d
    else
        default = convert(D,default_)
        insert!(m.bt,k, default, false)
        return default
    end
end

function getkey{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord}, k_, default_)
    i, exactfound = findkey(m.bt, convert(K,k_))
    exactfound? m.bt.data[i].k : convert(K, default_)
end

## Function delete! deletes an item at a given
## key

@inline function delete!(m::SortedDict, k_)
    i, exactfound = findkey(m.bt,convert(keytype(m),k_))
    !exactfound && throw(KeyError(k_))
    delete!(m.bt, i)
    m
end

@inline function pop!(m::SortedDict, k_)
    i, exactfound = findkey(m.bt,convert(keytype(m),k_))
    !exactfound && throw(KeyError(k_))
    d = m.bt.data[i].d
    delete!(m.bt, i)
    d
end

## Check if two SortedDicts are equal in the sense of containing
## the same (K,D) pairs.  This sense of equality does not mean
## that semitokens valid for one are also valid for the other.

function isequal(m1::SortedDict, m2::SortedDict)
    ord = orderobject(m1)
    if !isequal(ord, orderobject(m2)) || !isequal(eltype(m1), eltype(m2))
        error("Cannot use isequal for two SortedDicts unless their element types and ordering objects are equal")
    end
    p1 = startof(m1)
    p2 = startof(m2)
    while true
        if p1 == pastendsemitoken(m1)
            return p2 == pastendsemitoken(m2)
        end
        if p2 == pastendsemitoken(m2)
            return false
        end
        k1,d1 = deref((m1,p1))
        k2,d2 = deref((m2,p2))
        if !eq(ord,k1,k2) || !isequal(d1,d2)
            return false
        end
        p1 = advance((m1,p1))
        p2 = advance((m2,p2))
    end
end

function mergetwo!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord},
                                        m2::Associative{K,D})
    for (k,v) in m2
        m[convert(K,k)] = convert(D,v)
    end
end

function packcopy{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord})
    w = SortedDict(Dict{K,D}(),orderobject(m))
    mergetwo!(w,m)
    w
end

function packdeepcopy{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord})
    w = SortedDict(Dict{K,D}(),orderobject(m))
    for (k,v) in m
        newk = deepcopy(k)
        newv = deepcopy(v)
        w[newk] = newv
    end
    w
end

function merge!{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord},
                                     others::Associative{K,D}...)
    for o in others
        mergetwo!(m,o)
    end
end

function merge{K,D,Ord <: Ordering}(m::SortedDict{K,D,Ord},
                                    others::Associative{K,D}...)
    result = packcopy(m)
    merge!(result, others...)
    result
end

similar{K,D,Ord<:Ordering}(m::SortedDict{K,D,Ord}) =
SortedDict{K,D,Ord}(orderobject(m))


text:
@testset "TransitionalMCMC" begin

    disable_logging(Logging.Info)

    @testset "1D" begin
        Random.seed!(123456)
        lb, ub  = -15, 15

        μ = [0 5]
        σ = [1 0.2]
        w = [0.7 0.3]

        fT(x) = logpdf(Uniform(lb, ub), x[1])
        sample_fT(Nsamples) = rand(Uniform(lb, ub), Nsamples, 1)

        log_fD_T(x) = log(w[1] * pdf(Normal(μ[1], σ[1]), x[1]) + w[2] * pdf(Normal(μ[2], σ[2]), x[1]))

        Nsamples = 1000
        X, _ = tmcmc(log_fD_T, fT, sample_fT, Nsamples)

        m = sum(w .* μ)
        s = sqrt(sum(w .* (σ.^2 + μ.^2 .- m^2)))

        h0 = ExactOneSampleKSTest(vec(X), Normal(m, s))

        @test pvalue(h0) < 1e-4
    end

    @testset "2D" begin
        Random.seed!(1234)

        lb, ub  = -15, 15

        w = [0.6, 0.4]
        μ = [0 5; 0 5]

        Σ1 = [1 -0.5; -0.5 1]
        Σ2 = [1 0.5; 0.5 1]

        # Compute mean and cov of resulting gaussian mixture
        m = sum(w .* μ, dims=2) |> vec
        C = w[1] * Σ1 + w[2] * Σ2
        C += w[1] * (μ[:, 1] - m) * (μ[:, 1] - m)'
        C += w[2] * (μ[:, 2] - m) * (μ[:, 2] - m)'

        Nsamples = 10000
        Y = rand(MvNormal(m, C), Nsamples)'

        fT(x) = logpdf(Uniform(lb, ub), x[1]) .+ logpdf(Uniform(lb, ub), x[2])
        sample_fT(Nsamples) = rand(Uniform(lb, ub), Nsamples, 2)

        log_fD_T(x) = log.(w[1] * pdf(MvNormal(μ[:, 1], Σ1), x) + w[2] * pdf(MvNormal(μ[:, 2], Σ2), x))

        X, _ = tmcmc(log_fD_T, fT, sample_fT, Nsamples)

        h0 = OneSampleHotellingT2Test(X, m)
        @test pvalue(h0) < 0.05

        h0 = BartlettTest(X, Y)
        @test pvalue(h0) < 0.05
    end

    @testset "Himmelblau" begin
        Random.seed!(123456)
        lb, ub  = -5, 5

        # Prior Density and sampler
        logprior(x) = logpdf(Uniform(lb, ub), x[1]) + logpdf(Uniform(lb, ub), x[2])
        priorRnd(Nsamples) = rand(Uniform(lb, ub), Nsamples, 2)

        # Log Likelihood
        logLik(x) = -1 * ((x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7)^2)

        samps, acc = tmcmc(logLik, logprior, priorRnd, 20000)

        μ = mean(samps, dims=1) |> vec
        σ = std(samps, dims=1) |> vec
        corrs = cor(samps)

        # Reference values obtained using 1e7 samples
        @test [0.7, 0.2] < μ < [1, 0.4] # [0.832681 0.286942]
        @test [2.8, 2.1] < σ < [3.4, 2.8] # [3.16236 2.45603]
        @test corrs ≈ [1.0 0; 0 1.0] atol = 0.1 # independent
    end

    @testset "Himmelblau parallel" begin
        addprocs(2; exeflags="--project")
        @everywhere begin

            using TransitionalMCMC, Distributions, Random

            Random.seed!(123456)

            # Prior Bounds
            lb, ub  = -5, 5

            # Prior Density and sampler
            logprior(x) = logpdf(Uniform(lb, ub), x[1]) + logpdf(Uniform(lb, ub), x[2])
            priorRnd(Nsamples) = rand(Uniform(lb, ub), Nsamples, 2)

            # Log Likelihood
            logLik(x) = -1 * ((x[1]^2 + x[2] - 11)^2 + (x[1] + x[2]^2 - 7).^2)

        end

        Nsamples = 20000

        samps, acc = tmcmc(logLik, logprior, priorRnd, Nsamples)

        μ = mean(samps, dims=1) |> vec
        σ = std(samps, dims=1) |> vec
        corrs = cor(samps)

        # Reference values obtained using 1e7 samples
        @test [0.7, 0.2] < μ < [1, 0.4] # [0.832681 0.286942]
        @test [2.8, 2.1] < σ < [3.4, 2.8] # [3.16236 2.45603]
        @test corrs ≈ [1.0 0; 0 1.0] atol = 0.1 # independent

        rmprocs(workers())
    end

    Random.seed!()
end


text:
# generate the spatial kernel for an n-dimensional wave equation solver
include("generate_loops.jl")

# This function generates code to solve the n-dimensional
# constant coefficient wave equation in first order form, for dimensions 1
# to maxdim.  It uses 4th order accuate finite difference stencils.
# The function perf_hdindexing5 loops over all dimensions and calls the
# function kernel5 to evaluate the stencils at each point.  Dispatch on the
# array dimension ties the n different solvers together.
function make_stencil(fname, maxdim)
# generate kernels and loops for arrays of dimension 2 to maxdim + 1 and write
# to file named fname

  stencil = ["1/12", "-2/3", "0", "2/3", "-1/12"]
  neq = 2  # number of equations

  npts = length(stencil)
  nghost = div(npts, 2)

  str = generate_kernel(maxdim, stencil, neq)
  str *= "\n"
  str *= generate_loops(maxdim, npts)

  f = open(fname, "w")
  println(f, str)
  close(f)

end

function generate_kernel(maxdim, stencil, neq)
# maxdim is maximum dimension
# stencil is the stencil (vector of strings that are the coefficients)
# neq is the number of equations
# prefix is prepended to the file name
  npts = length(stencil)

  str = ""
  for i=1:maxdim
    func_sig = getKernelSignature(i, npts)
    str *= func_sig*"\n"

    func_prologue = getKernelPrologue(i)
    func_prologue = indentString(2, func_prologue)
    str *= func_prologue*"\n"

    func_body = generate_body(i, stencil, neq)
    func_body = indentString(2, func_body)
    str *= func_body*"\n"

    func_end = getKernelEpilogue()
    str *= func_end*"\n"

  end

  return str
end

function generate_body(maxdim, stencil, neq)
# writes to file kernel_dim_npts
# neq is the number of equations in the system

  npts = length(stencil)
  order = npts - 1
  second_line_indent = " "^8

  str = ""

  for eq=1:neq  # loop over equations in the system
    for dim=1:maxdim
      var_name = string("u", dim, dim, "_", eq)
      stencil_dim = getStencil(maxdim, stencil, dim, eq, second_line_indent)
      str *= string(var_name, " = ", stencil_dim)*"\n"

    end

    # sum the results into the receiving arrays
    kernel_assembly = getKernelAssembly(maxdim, eq, neq - eq + 1)
    str *= kernel_assembly*"\n"
  end

  return str

end

function getStencil(maxdim::Integer, stencil, dim::Integer, src_eq::Integer,
                    second_line_indent::String="")

# Generates a string that multiplies the given stencil coefficients by
#   the corresponding values from an array called u_i.
# The stencil coefficients are list from lowest index to highest index
#   maxdim is the number of indices the array u_i has
# The indices are named d1 through dn.
# dim is the dimension of u_i upon which to apply the array.
# Each term of the stencil is put on its own line.  second_line_indent
# is the string prepended to all lines except the first to indent them

# the src_eq an dest_eq are the equation indices for the right hand side
#   and left hand side arrays, and is the final index

# the stencil array is expected to be an array of strings
# the strings are enclosed in parenthesis before multiplication
# great pains are taken to ensure alignment of all terms of the stencil

# the returned string is of the form:
#    udd = stencil calculation...

  npts = length(stencil)
  stencil_startoffset = -div(npts - 1, 2)
  stencil_endoffset = div(npts - 1, 2)
  stencil_offsets = stencil_startoffset:stencil_endoffset

  # figure out maximum size of a stencil coefficient
  max_coeff_len = 0
  for i=1:npts
    len_i = length(stencil[i])
    if len_i > max_coeff_len
      max_coeff_len = len_i
    end
  end

  str = ""

  # index names are d1 through dmaxdim
  for pt=1:npts

    if pt == 1
      indentlevel = ""
    else
      indentlevel = second_line_indent
    end

    stencil_pt = stencil[pt]
    stencil_pt_len = length(stencil_pt)
    stencil_offset = stencil_offsets[pt]

    padding = max_coeff_len - stencil_pt_len
    coeff_str = string( " "^padding, "(", stencil_pt, ")")

    str_pt = indentlevel*coeff_str*"*u_i[ "

    # insert all indices before dim
    for i=1:(dim-1)
      idx_i = "d$i"
      str_pt = str_pt*idx_i*", "
    end

    # insert index dim and modifier

    # figure out if this should be a plus or minus
    stencil_offset_inner = abs(stencil_offset)
    if stencil_offset < 0
      mod_str = "d$dim - $stencil_offset_inner, "
    elseif stencil_offset > 0
      mod_str = "d$dim + $stencil_offset_inner, "
    else  # stencil_offset == 0
      mod_str = "d$dim    , "
    end

    # apply the index dim and modifier
    str_pt *= mod_str

    # insert all indices after dim
    for i=(dim+1):maxdim
      idx_i = "d$i"
      str_pt *= idx_i*", "
    end

    # remove last comma and space
    str_pt = str_pt[1:end-2]

    # close bracket
    str_pt = str_pt*", $src_eq"*"]"

    # figure out whether or not to add a plus sign for another term
    if pt != npts
      mod_str = " +\n"
    else
      mod_str = "\n"
    end

    str *= str_pt*mod_str

  end  # end loop over npts

  return str

end

function getKernelAssembly(maxdim::Integer, src_eq::Integer, dest_eq::Integer)
# sums the terms from each dimension into the receiving array
# src_eq tells which equation the terms come from, dest_eq tells which
# equation they are summed into

  second_line_indent = " "^(12 + 4*maxdim)
  # write receiving location
  str = "u_ip1[ "
  for i=1:maxdim
    str *= "d$i, "
  end

  str *= string(dest_eq, "] = ")
  # remove final comma and space, add close brack
#  str = str[1:end-2]
#  str *= "] = "

  # multiply stencil in each direction by grid spacing square
  for i=1:maxdim
    delta_name = string("delta_", i, "2")
    stencil_fac_name = string("u", i, i, "_", src_eq)
    str *= delta_name*"*"*stencil_fac_name*" + "

    if ((i % 3) == 0) && i != maxdim
      str *= "\n"*second_line_indent
    end

  end

  # remove last " + "
  str = str[1:end-3]
  str *= "\n"

  return str
end

function getKernelSignature(dim::Integer, npts::Integer)
# get the function signature

  kernel_name = string("kernel", npts)
  array_typetag = string("::AbstractArray{T,", dim+1, "}")

  str = string("function ", kernel_name, "(idx, ")
  str *= string("u_i", array_typetag, ", u_ip1", array_typetag, ") where T\n")

  return str
end

function getKernelPrologue(maxdim::Integer)
# unpack all the variables

  str = ""

  for i=1:maxdim
    varname = "d$i"
    str *= string(varname, " = idx[", i, "]\n")
  end

  str *= "\n"

  for i=1:maxdim
    varname = string("delta_", i, "2")
#    fieldname = string("deltax_invs2[", i, "]")
    str *= string(varname, " = ", "0.5", "\n")
  end

  str *= "\n"

  return str
end

function getKernelEpilogue()
  # get the return and end statement

  str = "  return nothing\n"
  str *= "end"

  return str
end

function indentString(indent::Integer, str::String)
# indent all lines of a string a certain number of spaces
  indent_str = " "^indent
  newline_indent = "\n"*indent_str
  str = indent_str*str
  str = replace(str, "\n" => newline_indent)

  return str
end


text:
export Biss_SecA_ls
function Biss_SecA_ls(h :: LineModel, stop_ls :: LS_Stopping;
                     f_meta = LS_Function_Meta(),
                     h₀ = obj(h, 0.0), g₀ = grad(h, 0.0),
                     φ_dφ :: Function = (x, y) -> phi_dphi(x, y),
                     γ :: Float64 = 0.8, verboseLS :: Bool = false,
                     kwargs...)

  state = stop_ls.current_state

  t = 1.0
  update!(state, x = t)
  φt, dφt = φ_dφ(h, state; kwargs...)
  tqnp = 0.0; dφtqnp = h₀ * (1.0 + 0.01)
  ht = obj(h, t); gt = grad(h, t)
  OK = update_and_start!(stop_ls, ht = ht, gt = gt, h0 = h₀, g0 = g₀)

  (ta, φta, dφta, tb, φtb, dφtb) = find_interval_ls(h, stop_ls,
                                                    verboseLS = verboseLS;
                                                    kwargs...)
  verboseLS && println("ta = $ta tb = $tb")

  t = ta
  tp = tb
  tqnp = tb

  iter = 0

  dφt = dφta; dφa = dφta; dφb = dφtb;

  verboseLS && @printf(" iter     tqnp        t         dφtqnp     dφt \n")
  verboseLS && @printf(" %4d %7.2e  %7.2e  %7.2e  %7.2e\n",
                        iter, tqnp, t, dφa, dφb)

  while !OK             # admissible: satisfies Armijo & Wolfe,
                        # tired: exceeds maximum number of iterations
    s = t - tqnp; y = dφt - dφtqnp;
    Γ = 3.0 * (dφt + dφtqnp) * s - 6.0 * (φt - dφtqnp)
    if y * s + Γ < eps(Float64) .* (s ^ 2)
      yt = y + Γ / s
    else
      yt = y
    end
    dN = -dφt * s / yt

    if ((tp - t) * dN > 0.0) && (dN / (tp - t) < γ)
      tplus = t + dN
      verboseLS && println("N")
    else
      tplus = (t + tp) / 2
      verboseLS && println("B")
    end
    update!(state, x = tplus)
    φplus, dφplus = φ_dφ(h, state; kwargs...)

    if t > tp
      if dφplus < 0.0
        tp = t
        tqnp = t
      else
        tqnp = t
      end
    else
      if dφplus > 0.0
        tp = t
        tqnp = t
      else
        tqnp = t
      end
    end

    t = tplus
    dφtqnp = dφt
    dφt = dφplus

    iter += 1
    OK = update_and_stop!(stop_ls, ht = obj(h, t), gt = grad(h, t))

    verboseLS && @printf(" %4d %7.2e  %7.2e  %7.2e  %7.2e\n",
                          iter, tqnp, t, dφtqnp, dφt)
  end

  return (stop_ls.current_state, stop_ls.meta.optimal)
end


text:
using TypeClasses
using DataTypesBasic
using Test

# MonoidAlternative
# =================

@test neutral(Option) == Option(nothing)

@test Option(3) ⊘ Option(4) == Option(3)  # take the first non-nothing
@test Option(nothing) ⊘ Option(4) == Option(4)  # take the first non-nothing
@test Option(nothing) ⊘ Option(4) == Option(4)  # take the first non-nothing
@test Option("hi") ⊕ Option("ho") == Option("hiho")

# FunctorApplicativeMonad
# =======================

@test map(Option(3)) do x
  x*x
end == Option(9)

@test map(Option(nothing)) do x
  x*x
end == Option(nothing)

@test mapn(Option(3), Option("hi")) do x, y
  "$x, $y"
end == Option("3, hi")

@test mapn(Option(3), Option(nothing)) do x, y
  "$x, $y"
end == Option()

ho = @syntax_flatmap begin
  a = Option(3)
  b = Option("hi")
  @pure "$a, $b"
end
@test ho == Option("3, hi")

ho = @syntax_flatmap begin
  a = Option(3)
  b = Option(nothing)
  @pure "$a, $b"
end
@test ho == Option(nothing)

@test pure(Option, 4) == Option(4)

# FlipTypes
# =========

@test flip_types(Identity([1,2,3])) == Identity.([1,2,3])
@test_throws MethodError flip_types(Const(nothing))


text:
#=
Copyright (c) 2015, Intel Corporation
All rights reserved.

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:
- Redistributions of source code must retain the above copyright notice,
  this list of conditions and the following disclaimer.
- Redistributions in binary form must reproduce the above copyright notice,
  this list of conditions and the following disclaimer in the documentation
  and/or other materials provided with the distribution.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
THE POSSIBILITY OF SUCH DAMAGE.
=#

module Comprehension

export @comprehend
import CompilerTools
import ..API

"""
Translate an ast whose head is :comprehension into equivalent code that uses cartesianarray call.
"""
function comprehension_to_cartesianarray(ast)
  assert(ast.head == :comprehension)
  body = ast.args[1]
  if isa(body, Expr) && body.head == :generator
    body = ast.args[1].args[1]
    args = ast.args[1].args[2:end]
    ndim = length(args)
  else
    body = ast.args[1]
    args = ast.args[2:end]
    ndim = length(args)
  end
  params = Array(Symbol, ndim)
  indices = Array(Symbol, ndim)
  ranges = Array(Any, ndim)
  headers = Array(Any, ndim)
  for i = 1:ndim
    r = args[i]
    assert(r.head == :(=))
    indices[i] = r.args[1]
    params[i] = gensym(string(indices[i]))
    ranges[i] = r.args[2]
    #headers[i] = Expr(:(=), indices[i], Expr(:call, GlobalRef(Base, :getindex), ranges[i], params[i]))
    #first(v) + (i-1)*step(v)
    first = GlobalRef(Base, :first)
    step = GlobalRef(Base, :step)
    headers[i] = :($(indices[i]) = $first($(ranges[i])) + ($(params[i]) - 1) * $step($(ranges[i])))
    # if range is 1:N, generate simple assignment instead of formula
    # this helps copy propagation for nested comprehensions (motivated by K-Means)
    if isa(ranges[i], Expr) && isStart1UnitRange(ranges[i])
        headers[i] = :($(indices[i]) = $(params[i]))
    end
  end
  args = Expr(:tuple, params...)
  dims = [ Expr(:call, GlobalRef(Base, :length), r) for r in ranges ]
  dimsExpr = Expr(:tuple, [ Expr(:call, GlobalRef(Base, :length), r) for r in ranges ]...)
  tmpret = gensym("tmp")
  tmpinits = [ :($idx = 1) for idx in params ]
  typetest = :(local $tmpret; if 1<0 let $(tmpinits...); $(headers...); $tmpret=$body end end)
  if VERSION < v"0.5.0-dev+5381"
      ast = Expr(:call, GlobalRef(API, :cartesianarray),
                    :($args -> let $(headers...); $body end), Expr(:call, GlobalRef(Core, :apply_type), GlobalRef(Base, :Tuple), Expr(:static_typeof, tmpret)), :($dimsExpr))
      ast = Expr(:block, typetest, ast)
  else
      ast = Expr(:call,  GlobalRef(API, :cartesianarray),
                    :($args -> let $(headers...); $body end), dims...)
  end
  ast
end

function isStart1UnitRange(node::Expr)
    if node.head==:(:) && length(node.args)==2 && node.args[1]==1
        return true
    end
    return false
end

"""
This function is a AstWalker callback.
"""
function process_node(node, state, top_level_number, is_top_level, read)
  if !isa(node,Expr)
    return CompilerTools.AstWalker.ASTWALK_RECURSE
  end
  if node.head == :typed_comprehension
    typ = node.args[1]
    # Transform into untyped :comprehension because type will be passed as a
    # parameter the same way we handle untyped :comprehensions
    node.head = :comprehension
    node.args = node.args[2:end]
  end
  return (node.head == :comprehension) ? comprehension_to_cartesianarray(node) : CompilerTools.AstWalker.ASTWALK_RECURSE
end

"""
Translate all comprehension in an AST into equivalent code that uses cartesianarray call.
"""
macro comprehend(ast)
  CompilerTools.AstWalker.AstWalk(ast, process_node, nothing)
  Core.eval(current_module(), ast)
end

end


text:
using Givre
using Test
using Base.GC: gc

# ENV["JULIA_DEBUG"] = "all"

ENV["DISPLAY"] = ":0"

function main()
    Givre.reset_timer!(Givre.to)
    app = Application()
    run(app)
    print(Givre.to)
end

main()


text:
module trunc_cyl_shell_examples
using FinEtools

function trunc_cyl_shell_mesh()
    h = 0.05 * phun("M");
    l = 10 * h;
    Rmed = h / 0.2;
    psi   = 0;    # Cylinder
    nh = 12; nl  = 40; nc = 120;
    nh = 6; nl  = 20; nc = 60;
    # nh = 3; nl  = 8; nc = 30;
    tolerance = h / nh / 100;

    t0 = time()
    MR = DeforModelRed3D
    fens, fes  = H8block(h, l, 2.0 * pi, nh, nl, nc)
    # Shape into a cylinder
    R = zeros(3, 3)
    for i = 1:count(fens)
        x, y, z = fens.xyz[i,:];
        rotmat3!(R, [0, z, 0])
        Q = [cos(psi * pi / 180) sin(psi * pi / 180) 0;
        -sin(psi * pi / 180) cos(psi * pi / 180) 0;
        0 0 1]
        fens.xyz[i,:] = reshape([x + Rmed - h / 2, y - l / 2, 0], 1, 3) * Q * R;
    end
    println("  before merging  = $(count(fens))")
    # File =  "unit_cube_modes.vtk"
    # vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)

    candidates = selectnode(fens, box = boundingbox([Rmed - h -Inf 0.0; Rmed + h +Inf 0.0]), inflate = tolerance)
    fens, fes = mergenodes(fens, fes,  tolerance, candidates);

    # fens,fes = mergenodes(fens, fes,  tolerance);

    println("  after merging  = $(count(fens))")

    println("Mesh generation ($(time() - t0) sec)")

    File =  "trunc_cyl_shell_mesh.vtk"
    vtkexportmesh(File, fens, fes)
    @async run(`"paraview.exe" $File`)

end # trunc_cyl_shell_mesh

function trunc_cyl_shell_mesh_2()
    h = 0.05*phun("M");
    l = 10*h;
    Rmed = h/0.2;
    psi   = 0;    # Cylinder
    nh = 12; nl  = 40; nc = 120;
    nh = 6; nl  = 20; nc = 60;
    # nh = 3; nl  = 8; nc = 30;
    tolerance = h/nh/100;

    t0 = time()
    MR = DeforModelRed3D
    fens,fes  = H8block(h,l,2.0*pi,nh,nl,nc)
    # Shape into a cylinder
    R = zeros(3, 3)
    for i = 1:count(fens)
        x, y, z = fens.xyz[i,:];
        rotmat3!(R, [0, z, 0])
        Q = [cos(psi*pi/180) sin(psi*pi/180) 0;
        -sin(psi*pi/180) cos(psi*pi/180) 0;
        0 0 1]
        fens.xyz[i,:] = reshape([x+Rmed-h/2, y-l/2, 0], 1, 3)*Q*R;
    end
    println("  before merging  = $(count(fens))")
    # File =  "unit_cube_modes.vtk"
    # vtkexportmesh(File, fens, fes)
    # @async run(`"paraview.exe" $File`)

    candidates = selectnode(fens, plane = [0.0 0.0 1.0 0.0], inflate = tolerance)
    fens,fes = mergenodes(fens, fes,  tolerance, candidates);

    # fens,fes = mergenodes(fens, fes,  tolerance);

    println("  after merging  = $(count(fens))")

    println("Mesh generation ($(time() - t0) sec)")

    File =  "trunc_cyl_shell_mesh_2.vtk"
    vtkexportmesh(File, fens, fes)
    @async run(`"paraview.exe" $File`)

end # trunc_cyl_shell_mesh_2

function allrun()
    println("#####################################################")
    println("# trunc_cyl_shell_mesh ")
    trunc_cyl_shell_mesh()
    println("#####################################################")
    println("# trunc_cyl_shell_mesh_2 ")
    trunc_cyl_shell_mesh_2()
    return true
end # function allrun

end # module trunc_cyl_shell_examples


text:
using MicMods
using Documenter

DocMeta.setdocmeta!(MicMods, :DocTestSetup, :(using MicMods); recursive=true)

makedocs(;
    modules=[MicMods],
    authors="Thomas Wutzler <progtw@arcor.de> and contributors",
    repo="https://github.com/bgctw/MicMods.jl/blob/{commit}{path}#{line}",
    sitename="MicMods.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://bgctw.github.io/MicMods.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/bgctw/MicMods.jl",
    devbranch="main",
)


text:
### A Pluto.jl notebook ###
# v0.14.1

using Markdown
using InteractiveUtils

# ╔═╡ 5eebac68-891b-11eb-37ba-c1d9481c6134
begin
    using Pkg; Pkg.activate("MLJ_env", shared=true)
	using Test
	using Random
	using Distributions
	using PlutoUI
	using Printf
	using JSON
	using Plots

	push!(LOAD_PATH, "./src")
	using YaCounter
	using YaDeepLearning
	using YaDistances: cosine_similarity
	using .TensorDT: Tensor
end

# ╔═╡ 87e5e2ec-8915-11eb-362b-6bf13a36b8e4
md"""
## NLP, part 2 - Word Vectors

ref. from book **"Data Science from Scratch"**, Chap 21

$(html"
© Pascal, April 2021
")
"""

# ╔═╡ 83e01f80-884a-11eb-1ef7-e9cb5ab3ed20
PlutoUI.TableOfContents(indent=true, depth=4, aside=true)

# ╔═╡ 853ab02c-884a-11eb-01bf-1712dc177ffb
html"""

"""

# ╔═╡ 8cf46618-9109-11eb-2094-d99ba49a5b89
begin
  const D_SI = Dict{String, Integer}
  const D_IS = Dict{Integer, String}
  const S_N = Union{String, Nothing}
  const I_N = Union{Integer, Nothing}
  const T_N = Union{Tensor, Nothing}
  const F = Float64
end

# ╔═╡ 83d71b10-90d3-11eb-1c48-1ddf8be8bb57
begin
  import .AbstractLayers: AbstractLayer, forward, backward, parms, ∇parms
  using .Initializations: init_rand_normal

  struct Embedding <: AbstractLayer
    # @add_embedding_fields
    num_embedding::Integer
    embedding_dim::Integer
    embeddings::Tensor
    ∇::Tensor
    last_input_id::Vector{I_N}
    _type::Symbol
    store::Dict # {Any. Any}

    function Embedding(num_embedding::Integer, embedding_dim::Integer)
      shape = (embedding_dim, )
      embeddings = [init_rand_normal(shape) for _ ∈ 1:num_embedding]
      ∇ = [zeros(eltype(embeddings[1]), shape[1]) for _ ∈ 1:num_embedding]
      #
      new(num_embedding, embedding_dim, embeddings, ∇, [nothing], :embedding,
        Dict())
    end
  end

  function forward(self::Embedding, vinput_id::Vector{T})::Tensor where T <: Integer
    input_id = vinput_id[1]           ## only 1 element
    self.store[:input_id] = input_id  ## requires for backprop.
    self.embeddings[input_id]
  end

  function backward(self::Embedding, ∇p::Tensor)
    ## Zero out the gradient corresponding to the last input.
    ## This is cheaper than creating a new all-zero tensor each time.
    if !isnothing(self.last_input_id[1])
      zero_row = zeros(eltype(self.embeddings[1]), self.embedding_dim)
      self.∇[self.last_input_id[1]] = zero_row
    end
    self.last_input_id[1] = self.store[:input_id]
    self.∇[self.store[:input_id]] = vec(∇p)
  end

  parms(self::Embedding) = [self.embeddings]
  ∇parms(self::Embedding) = [self.∇]

  ## Extra
  idims(self::Embedding) = self.num_embedding
  odims(self::Embedding) = self.embedding_dim
end

## ================================

# ╔═╡ 5f5e6d02-892a-11eb-32d8-dfc1a7f259fa
md"""
### Vocabulary
"""

# ╔═╡ c1cb6562-90c8-11eb-2095-cbf63bb8b59d
md"""
We will start wit h a toy dataset to learn some word vectors.
"""

# ╔═╡ 073f6442-892e-11eb-0fc2-bbeb1d12619e
begin
  const colors = ["red", "green", "blue", "yellow", "black", ""]
  const nouns = ["bed", "car", "boat", "cat"]
  const verbs = ["is", "was", "seems"]
  const adverbs = ["very", "quite", "extremely", ""]
  const adjectives = ["slow", "fast", "soft", "hard"]

  const Num_Sentences = 50
  Random.seed!(70);
end

# ╔═╡ 2ec41c48-8a95-11eb-13b6-61f4a07a31db
function make_sentences()::String
  λ = coll -> coll[rand(1:length(coll))]
  [
    "The",
    λ(colors),
    λ(nouns),
    λ(verbs),
    λ(adverbs),
    λ(adjectives),
    "."
  ] |> a -> join(a, " ")
end

# ╔═╡ 045023be-8aac-11eb-29c1-3f4da066ff9b
const Sentences = [make_sentences() for _ in 1:Num_Sentences];

# ╔═╡ 04ab0ce8-90d3-11eb-3567-7f8f35e74e6a
md"""
We will want to one-hot-encode our words, which means we will need to convert them to IDs. We will introduce a Vocabulary struct to keep track of this mapping:
"""

# ╔═╡ b1c0a416-8aba-11eb-140a-9f3c4eb93346
begin
  const JSON_DATA_DIR = "./vocab_data"

  struct Vocabulary
    word2ix::D_SI
    ix2word::D_IS

    Vocabulary() = new(D_SI(), D_IS())

    Vocabulary(w2ix::D_SI, ix2w::D_IS) = new(w2ix, ix2w)

    function Vocabulary(words::Vector{String})
      @assert length(words) > 0
      new(_add(words)...)
    end
  end

  ## Public API

  function get_id(self::Vocabulary, w::String)::I_N
    !haskey(self.word2ix, w) && (return nothing)
    self.word2ix[w]
  end

  function get_word(self::Vocabulary, ix::Integer)::S_N
    !(1 ≤ ix ≤ length(self.word2ix)) && (return nothing)
    self.ix2word[ix]
  end

  size(self::Vocabulary) = length(self.word2ix)

  function one_hot_encode(self::Vocabulary, w::String)::T_N
    !haskey(self.word2ix, w) && (return nothing)
    w_ix = self.word2ix[w]
    [ix == w_ix ? 1. : 0. for ix ∈ 1:length(self.word2ix)]
    # FIXME: use a BitArray?
  end

  function add(self::Vocabulary, word::String)
    if !haskey(self.word2ix, word)
      w_ix = length(self.word2ix) + 1
      self.word2ix[word] = w_ix
      self.ix2word[w_ix] = word
    end
  end

  function save_vocab(self::Vocabulary, filename::String)
    open(string(JSON_DATA_DIR, "/", filename), "w") do f
        JSON.print(f, self.word2ix)
    end
  end

  function load_vocab(filename::String)
    word2ix = open(string(JSON_DATA_DIR, "/", filename), "r") do f
          JSON.parse(f)
    end |> h -> D_SI(h)
    ix2word = D_IS(ix => w for (w, ix) ∈ word2ix)
    Vocabulary(word2ix, ix2word)
  end

  ## Internal

  function _add(words::Vector{String})
    w2ix = D_SI()
    ix2w = D_IS()
    for w ∈ words
      haskey(w2ix, w) && continue
      w_ix = length(w2ix) + 1
      w2ix[w] = w_ix
      ix2w[w_ix] = w
    end
    (w2ix, ix2w)
  end
end

# ╔═╡ 65ec8554-8953-11eb-3973-b56039754312
begin
  const vlst = ["B", "A", "G", "E", "A", "M", "P", "M"]
  vocab = Vocabulary(vlst)

  @test size(vocab) == (length ∘ unique)(vlst)
  @test get_id(vocab, "A") == 2
  @test one_hot_encode(vocab, "B") == [1., 0., 0., 0., 0., 0.]
  @test isnothing(get_id(vocab, "z"))
  @test get_word(vocab, 6) == "P"
  @test isnothing(get_word(vocab, 10))

  add(vocab, "Y")
  @test size(vocab) == (length ∘ unique)(vlst) + 1
  @test one_hot_encode(vocab, "B") == [1., 0., 0., 0., 0., 0., 0.]
  @test one_hot_encode(vocab, "Y") == [0., 0., 0., 0., 0., 0., 1.]
  @test isnothing(one_hot_encode(vocab, "z"))
end

# ╔═╡ 0c8055b8-89b2-11eb-3046-331119a0dc9b
begin
  # use previous example to test save / load JSON
  save_vocab(vocab, "vocab.json")

  l_vocab = load_vocab("vocab.json")

  @test size(l_vocab) == (length ∘ unique)(vlst) + 1
  @test one_hot_encode(l_vocab, "B") == [1., 0., 0., 0., 0., 0., 0.]
  @test get_id(l_vocab, "A") == 2
  @test one_hot_encode(l_vocab, "Y") == [0., 0., 0., 0., 0., 0., 1.]
  @test isnothing(one_hot_encode(l_vocab, "z"))
end

# ╔═╡ ce9e2400-91a9-11eb-1e99-6d75ec2a3078
html"""

  back to TOC

"""

## ================================

# ╔═╡ cdf7edf6-91a9-11eb-1d6c-a5288d3fd950
md"""
### Embedding Layer
"""

# ╔═╡ d728d526-8956-11eb-3c22-e788d025e8b4
md"""
We will be using a word vector model called *skip-gram) that takes as input a word and generates probabilities for what words are likely to be seen near it.

We will feed this model training pairs (word, nearby_word) and try to minimize the SoftmaxCrossEntropy loss.

Let us design our neural network to achieve this. At its heart will be an embedding layer that takes as input a word ID and returns a word vector. Under the covers we can just use a lookup table for this (yes, an embedding is a lookup table).

Then we will pass the word vector to a Linear layer with the same number of outputs as we have words in our vocabulary.
As before, we will use softmax to convert these outputs to probabilities over nearby words. As we use gradient descent to train the model, we will be updating the vectors in the lookup table. Once we have finished the training, that lookup table gives us our word vectors.

Let us create that embedding layer. In practice we might want to embed things other than words, so we will construct a more general Embedding laye (later we will write a TextEmbedding layer that ss specifically for word vectors). In its constructor we will provide the number and dimension of our embedding vectors, so it can create the embeddings (which will be standard random normals, initially):
"""

# ╔═╡ c228148c-90e0-11eb-1007-43a7de1d22ee
# ## define a macro which build a macro to define fieldnames that can be added
# ## across several structs

# macro def(name, definition)
#   quote
#       macro $(esc(name))()
#           esc($(Expr(:quote, definition)))
#       end
#   end
# end

# ╔═╡ e6634800-90e0-11eb-2d3d-c58c08564720
# @def add_embedding_fields begin
#   num_embedding::Integer
#   embedding_dim::Integer
#   embeddings::Tensor
#   ∇::Tensor
#   last_input_id::I_N
#   _type::Symbol
#   store::Dict # {Any. Any}
# end

# ╔═╡ ec95a8d0-90d6-11eb-3713-9febeeceb15c
md"""
Remarks:
  - In our case we will only be embedding one word at a time.
  - For the backward pass we will get a gradient corresponding to the chosen embedding vector, and we will need to construct the corresponding gradient for self.embeddings, which is zero for every embedding other than the chosen one.
"""

# ╔═╡ f3fa0216-90da-11eb-157a-e3941af7117b
md"""
Now let us create a specific struct for `TextEmbdding`
"""

# ╔═╡ 17768334-8a98-11eb-1172-3d5079435dcf
begin
  struct TextEmbedding <: AbstractLayer
    # @add_embedding_fields
    vocab::Vocabulary
    embedding::Embedding
    _type::Symbol

    function TextEmbedding(vocab::Vocabulary, embedding_dim::Integer)
      embedding = Embedding(size(vocab), embedding_dim)
      new(vocab, embedding, :embedding)
    end
  end

  ## delegate
  for (fn, type_) ∈ [
      (:forward, Vector{Int64}),
      (:backward, Tensor),
      (:parms, nothing),
      (:∇parms, nothing),
      (:idims, nothing),
      (:odims, nothing)]

    @eval begin
      if isnothing($(type_))
        ($fn)(self::TextEmbedding) = ($fn)(self.embedding)
      else
        ($fn)(self::TextEmbedding, p::$(type_)) = ($fn)(self.embedding, p)
      end
    end
  end

  ## Specific TextEmbedding
  function Base.getindex(self::TextEmbedding, word::String)::Tensor
    """
    to allow us to call something like: embedding["black"]
    """
    w_id = get_id(self.vocab, word)
    isnothing(w_id) && (return nothing)
    self.embedding.embeddings[w_id]
  end

  get_embeddings(self::TextEmbedding) = self.embedding.embeddings

  function closest(self::TextEmbedding, word::String; n=5)::Vector{Tuple{F, String}}
    """
    returns the n closest words based on cosine similarity
    """
    vect = self[word]  ## getindex

    scores = [
      (cosine_similarity(vect, get_embeddings(self)[ix]), o_word)
      for (o_word, ix) ∈ self.vocab.word2ix
    ]
    sort(scores, rev=true)[1:n]
  end
end

# ╔═╡ e27e6ed6-90e7-11eb-17a7-79104767ea17
md"""
Our embedding layer is outputting vectors which can be fed into a linear(dense) layer.

But first we need to tokenize the sentences.
"""

# ╔═╡ 1ed2ac2a-91aa-11eb-2b74-7d00f8d91a53
html"""

  back to TOC

"""

## ================================

# ╔═╡ 00982dfc-91aa-11eb-2e9d-f5831922a61d
md"""
### Run it
"""

# ╔═╡ f3d7e0be-90f3-11eb-00c7-036c1168da05
tokenized_sentences =  split.(Sentences) |>
  vvs -> [lowercase.(vs) for vs ∈ vvs] |>
  vvs -> filter.(w -> occursin(r"\A(?:[a-z]+|\.)\z", w), vvs)

# ╔═╡ e5532350-90e9-11eb-2e60-591759fb2f28
# Create a vocabulary (that is, a mapping word -> word_id) based on our text.
the_vocab = tokenized_sentences |>
  vvs -> reduce((rv, vs) -> rv = [rv..., vs...], vvs; init=[]) |>
  Vocabulary

# ╔═╡ ba364aec-90f5-11eb-3093-71206176896b
begin
  using .Layers

  Random.seed!(42)

  const EMBEDDING_DIM = 5
  embedding_layer = TextEmbedding(the_vocab, EMBEDDING_DIM)

  the_model = Layers.Sequential([
      ## 1. Given a word (vector of word_ids), look up its embedding.
      embedding_layer,

      ## 2. Use a linear layer to compute scores for "nearby" words
      Layers.Linear(EMBEDDING_DIM, size(the_vocab))
  ]);
end

# ╔═╡ d6b82be6-9112-11eb-110c-036fa8880c8a
begin
  using YaWorkingData: pca, transform

  comps = pca(embedding_layer.embedding.embeddings, 2)  # first 2 components
  transfo = transform(embedding_layer.embedding.embeddings, comps);

  coll = collect(zip(transfo...));

  v = []
  for (word, ix) ∈ the_vocab.word2ix
    (x, y) = transfo[ix]
    push!(v,
      (x, y, Plots.text(word, :red, :left, 8)))
  end
end

# ╔═╡ bb601544-90f3-11eb-1f6a-9758edfa0c59
md"""
From which we can create the training data:
"""

# ╔═╡ bb3fa48a-90f3-11eb-184e-7b04752115c1
begin
  # using .TensorDT: Tensor

  inputs = Vector{Int}()
  targets = Vector{Tensor}()

  for sentence ∈ tokenized_sentences, (ix, word) ∈ enumerate(sentence)
    for jx ∈ ix-2:ix+2
      jx == 0 && continue
      if 0 ≤ jx < length(sentence)
        nearby_word = sentence[jx]
        push!(inputs, get_id(the_vocab, word))
        push!(targets, one_hot_encode(the_vocab, nearby_word))
      end
    end
  end
end

# ╔═╡ ba1677d2-90f5-11eb-1779-05df1ee68db6
begin
  using .Loss
  using .Optimizers

  the_loss_fn = Loss.SoftmaxXEntropy()
  the_opt_fn = Optimizers.GD(0.01)

  for epoch ∈ 1:500
    epoch_loss = 0.
    for (x, y) ∈ zip(inputs, targets)
      ŷ = forward(the_model, [x])
      epoch_loss += Loss.loss(the_loss_fn, ŷ, y)
      ∇p = Loss.∇loss(the_loss_fn, ŷ, y)
      Layers.backward(the_model, ∇p)
      Optimizers.astep(the_opt_fn, the_model)
    end
    if epoch % 50 == 0
      @printf("%3d => %7.5f\n", epoch, epoch_loss) ## Print the loss
      println(closest(embedding_layer, "black"))   ## and also a few nearest words
      println(closest(embedding_layer, "slow"))    ## so we can see what's being
      println(closest(embedding_layer, "car"))     ## learned
    end
  end
end

# ╔═╡ bb248f9c-90f3-11eb-3766-e9c5cfae3920
length(inputs), length(targets)

# ╔═╡ bb0b45c8-90f3-11eb-2087-dd660f4305d4
md"""
And now we can build our model:
"""

# ╔═╡ da063a92-9110-11eb-1613-f51d8744fd95
begin
  pairs = [
    (cosine_similarity(embedding_layer[w1], embedding_layer[w2]), w1, w2)
    for (w1,) ∈ the_vocab.word2ix for (w2,) ∈ the_vocab.word2ix if w1 < w2
  ]
  sort!(pairs, rev=true)
  with_terminal() do
    for p ∈ pairs[1:7]
      println(p)
    end
  end
end

# ╔═╡ 78e07c8c-9133-11eb-2721-cdb8c99043b5
md"""
Let us extract the first two principal components and plot them:
"""

# ╔═╡ 8153d3ac-9113-11eb-1e9c-5dabd3012046
begin
  scatter(F[coll[1]...], F[coll[2]...], color=:lightgreen, marker=:circle,
    markersize=2, size=(600, 500), legend=false)
  annotate!(v...)
end

# ╔═╡ c484c36e-9133-11eb-22ef-dfdcfa0d75ec
md"""
This shows that similar words tends to cluster together (almost...)
"""

# ╔═╡ Cell order:
# ╟─87e5e2ec-8915-11eb-362b-6bf13a36b8e4
# ╟─83e01f80-884a-11eb-1ef7-e9cb5ab3ed20
# ╟─853ab02c-884a-11eb-01bf-1712dc177ffb
# ╠═5eebac68-891b-11eb-37ba-c1d9481c6134
# ╠═8cf46618-9109-11eb-2094-d99ba49a5b89
# ╟─5f5e6d02-892a-11eb-32d8-dfc1a7f259fa
# ╟─c1cb6562-90c8-11eb-2095-cbf63bb8b59d
# ╠═073f6442-892e-11eb-0fc2-bbeb1d12619e
# ╠═2ec41c48-8a95-11eb-13b6-61f4a07a31db
# ╠═045023be-8aac-11eb-29c1-3f4da066ff9b
# ╟─04ab0ce8-90d3-11eb-3567-7f8f35e74e6a
# ╠═b1c0a416-8aba-11eb-140a-9f3c4eb93346
# ╠═65ec8554-8953-11eb-3973-b56039754312
# ╠═0c8055b8-89b2-11eb-3046-331119a0dc9b
# ╟─ce9e2400-91a9-11eb-1e99-6d75ec2a3078
# ╟─cdf7edf6-91a9-11eb-1d6c-a5288d3fd950
# ╟─d728d526-8956-11eb-3c22-e788d025e8b4
# ╠═c228148c-90e0-11eb-1007-43a7de1d22ee
# ╠═e6634800-90e0-11eb-2d3d-c58c08564720
# ╠═83d71b10-90d3-11eb-1c48-1ddf8be8bb57
# ╟─ec95a8d0-90d6-11eb-3713-9febeeceb15c
# ╟─f3fa0216-90da-11eb-157a-e3941af7117b
# ╠═17768334-8a98-11eb-1172-3d5079435dcf
# ╟─e27e6ed6-90e7-11eb-17a7-79104767ea17
# ╟─1ed2ac2a-91aa-11eb-2b74-7d00f8d91a53
# ╟─00982dfc-91aa-11eb-2e9d-f5831922a61d
# ╠═f3d7e0be-90f3-11eb-00c7-036c1168da05
# ╠═e5532350-90e9-11eb-2e60-591759fb2f28
# ╟─bb601544-90f3-11eb-1f6a-9758edfa0c59
# ╠═bb3fa48a-90f3-11eb-184e-7b04752115c1
# ╠═bb248f9c-90f3-11eb-3766-e9c5cfae3920
# ╟─bb0b45c8-90f3-11eb-2087-dd660f4305d4
# ╠═ba364aec-90f5-11eb-3093-71206176896b
# ╠═ba1677d2-90f5-11eb-1779-05df1ee68db6
# ╠═da063a92-9110-11eb-1613-f51d8744fd95
# ╟─78e07c8c-9133-11eb-2721-cdb8c99043b5
# ╠═d6b82be6-9112-11eb-110c-036fa8880c8a
# ╠═8153d3ac-9113-11eb-1e9c-5dabd3012046
# ╟─c484c36e-9133-11eb-22ef-dfdcfa0d75ec

max_stars_repo_path: Data Science From Scratch/21-NLP-2.jl

max_stars_repo_name: pascal-p/julia-notebooks

max_issues_repo_path: Data Science From Scratch/21-NLP-2.jl

max_issues_repo_name: pascal-p/julia-notebooks

max_forks_repo_path: Data Science From Scratch/21-NLP-2.jl

max_forks_repo_name: pascal-p/julia-notebooks


text:
#! /usr/bin/julia

# Rosetta Code, Non-decimal radices/Output

hi = 50
println("Primes \u2264 ", hi, " written in common bases.")
println("     bin   oct  dec  hex")
for i in primes(hi)
    println(@sprintf("%8s %4s   %2s   %2s", bin(i), oct(i), dec(i), hex(i)))
end

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Deldir"
version = v"0.1.15"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://cran.r-project.org/src/contrib/Archive/deldir/deldir_0.1-15.tar.gz", "571011c208829f47ecd9f92a19fd94a8eb59de5a2645ab8c62e73926ade30710")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/deldir/src

for f in *.f; do
    ${FC} -fPIC -O2 -pipe -g -c "${f}" -o "$(basename "${f}" .f).o"
done

mkdir -p "${libdir}"
${CC} -shared -o ${libdir}/libdeldir.${dlext} *.o

install_license /usr/share/licenses/GPL3
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libdeldir", :libdeldir)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# Chain Potts Stein Distribution
#
# Represents a chain Potts distribution. Components are marginally +1 or -1
# with equal probability. The probability of (x_1,...x_p) is proportional to
#   exp(-beta * (x_1 x_2 + ... + x_{p-1} x_p)).
# The normalizing constant is 2*(2 cosh beta)^(p-1). beta is an affinity
# parameter; if beta > 0, then adjacent nodes tend to have opposite spins,
# and the opposite is true for beta < 0.
#
# We encode -1 as 0.

using StatsBase: wsample

type SteinPotts <: SteinDistribution
    # Affinity parameter
    beta::Float64
    # Dimensionality (number of nodes in the chain)
    p::Int64
    # Stein factors
    c1::Float64
    c2::Float64
    c3::Float64
end

# Constructor with (c1,c2,c3) = (1,1,1)
SteinPotts(beta::Float64, p::Int64) = SteinPotts(
    beta,
    p,
    1.0,
    1.0,
    1.0
);

# Draw n independent samples from distribution
function rand(d::SteinPotts, n::Int64)
    p = d.p;
    beta = d.beta;
    x = zeros(n,p);
    x[:,1] = rand(0:1,n);
    for i=1:n
        for j=2:p
            if x[i,j-1] == 1
                x[i,j] = wsample([0,1],[exp(beta),exp(-beta)])
            else
                x[i,j] = wsample([0,1],[exp(-beta),exp(beta)])
            end
        end
    end
    x
end

# Joint distribution evaluated at a single point
function jointdistribution(d::SteinPotts, y::Array{Float64})
    y = 2*y-1;     # Convert to +1 and -1
    exp(-d.beta*dot(y[1:d.p-1],y[2:d.p]))
end

# Conditional distribution of jth coordinate of y, given all other coordinates,
# represented as a probability vector
function conddistribution(d::SteinPotts, y::Array{Float64,1}, j::Int64)
    conddist = [jointdistribution(d, [y[1:(j-1)], yj, y[(j+1):end]])
                for yj in 0:1]
    conddist / sum(conddist)
end

function supportlowerbound(d::SteinPotts, j::Int64)
    0
end

function supportupperbound(d::SteinPotts, j::Int64)
    1
end

function numdimensions(d::SteinPotts)
    d.p
end


text:
using QuantumOpticsBase
using StaticArrays
using LinearAlgebra

# using BenchmarkTools
# using Dates

include("./functions.jl")
using .functions

# piecewise constant propagator
function full_propagator(c_vec, H0)
    exp((H0 + c_vec[1] .* sx .+ c_vec[2] .* sy) .* (-Δt*π*im))
end

# note the definitions that follow are mostly const, this means that they are
# considered constant and should have no performance impact, globals will
# considerably slow performance.

# construct a spin basis
b = SpinBasis(1//2)
# gather the Pauli matrices
const sx = SMatrix{2,2}(DenseOperator(sigmax(b)).data)
const sy = SMatrix{2,2}(DenseOperator(sigmay(b)).data)
const sz = SMatrix{2,2}(DenseOperator(sigmaz(b)).data)
const i2 = SMatrix{2,2}(DenseOperator(identity(b)).data)

# define two spin states
const ρ0 = @SVector [1 + 0.0im, 0.0+0.0im]
# target state
const σ0 = @SVector [0.0 + 0.0im, 1.0 + 0.0im]

# for use with the ensemble
const N_ensemble = 3
# define the weight array
const weights = @SArray ones(N_ensemble)
const weights = weights/sum(weights)

# define a list of drift Hamiltonians
const H_drift = collect(range(-100, 100, length = N_ensemble)) .* [sz]
# define a list of target states
const σ0_ens = repeat([ρ0], N_ensemble)
# specify some states that we want to flip
σ0_ens[Int(floor(N_ensemble/2))] = σ0
σ0_ens[Int(floor(N_ensemble/4))] = σ0
σ0_ens[Int(floor(N_ensemble/4*3))] = σ0

const N = 5 # number of time slices
const K = 2 # number of controls

const T = 2 # total duration of the sequence in us
const Δt = T/N # time step

# initial random controls
c_mat = (rand(K, N) .* 2 .- 1)*10

# functional that maps control vector x into fidelity, here averaged over the ensemble
function J4(x)
    err = 0
    x = complex.(real.(x))
    for i = 1:N_ensemble
        U = stack_props(x, H_drift[i])
        err += (1 - abs2(σ0_ens[i]' * U * ρ0)) * weights[i]
    end
    err
end

# function that computes the propagator across many time slices
# ends up being faster than using mapreduce
function stack_props(c_mat,H0)
    init = I(2)
    for i in 1:size(c_mat)[2]
        init = init * full_propagator(c_mat[:,i], H0)
    end
    init
end

J4(c_mat)

using Zygote
Zygote.gradient(J4, c_mat)

# Zygote.refresh() # useful to have laying around

# optimise
out = GRAPE(J4, c_mat*0.01, K, N)

J4(c_mat)

# make some quick plots to check the output looks okay
using Plots

x = out.minimizer #c_mat * 0.1
res = 501
dist = 100
z = []
z2 = []
for i = -res:res
    i = i/res * dist
    H0 = -i * sz
    #U = mapreduce(c->full_propagator(c, H0), *, eachcol(x))
    U = stack_props(x, H0)
    state = U * ρ0
    append!(z, real(state' * sz * state))
end

x = range(-dist, dist, length = length(z))
plot(x, z, xlabel = "Detuning (MHz)", ylabel = "z projection", label = "State")
x2 = range(-dist, dist, length = N_ensemble)
z45 = [real(i[1]) for i in σ0_ens]
plot!(x2, z45.*2 .-1, label = "Ideal state")


text:
@testset "reg-mdl" begin
    m = ARDRegressor()
    @test m isa MB.Deterministic
    @test m.n_iter > 0
    # set
    m = ARDRegressor(n_iter=100)
    @test m.n_iter == 100
    # clean method
    @test @test_logs (:warn, r"Constraint") ARDRegressor(n_iter=-5).n_iter > 0

    # traits
    @test MB.load_path(m) == "MLJScikitLearnInterface.ARDRegressor"
    @test MB.package_name(m) == "ScikitLearn"
end

@testset "clf-mdl" begin
    m = LogisticClassifier()
    @test m isa MB.Probabilistic
    @test m.penalty == "l2"
    # set
    m = LogisticClassifier(fit_intercept=false)
    @test !m.fit_intercept
    # clean
    @test @test_logs (:warn, r"Constraint") LogisticClassifier(
            penalty="foo").penalty == "l2"

    # traits
    @test MB.load_path(m) == "MLJScikitLearnInterface.LogisticClassifier"
    @test MB.package_license(m) == "BSD"
end


text:
@testset "LinShrink: target F with LW (ref⭒) " begin
    lw = LinearShrinkage(ConstantCorrelation())
    testTransposition(lw, X)
    testUncorrelated(lw)
    testTranslation(lw, X)
    testDims(lw, X)
    for X̂ ∈ test_matrices
        ref_results = matlab_ledoitwolf_covcor(X̂)
        lwfixed = LinearShrinkage(ConstantCorrelation(), ref_results["shrinkage"])
        @test cov(lw, X̂) ≈ ref_results["lwcov"]
        @test cov(lwfixed, X̂) ≈ ref_results["lwcov"]
    end
end

@testset "LinShrink: target D with SS (ref⭒) " begin
    p = ifelse(endswith(pwd(), "CovarianceEstimation"), "test", "")
    ## R Script used to compare:
    # require(corpcor)
    # tm1 = read.table("20x100.csv")
    # tm2 = read.table("100x20.csv")
    # tm3 = read.table("50x50.csv")
    # c1 = cov.shrink(tm1, lambda.var=0.0)
    # c2 = cov.shrink(tm2, lambda.var=0.0)
    # c3 = cov.shrink(tm3, lambda.var=0.0)

    ss = LinearShrinkage(target=DiagonalUnequalVariance(),
                                  shrinkage=:ss; corrected=true)
    test_mat1 = readdlm(joinpath(p, "test_matrices/20x100.csv"))
    ref_cov1  = readdlm(joinpath(p, "test_matrices/20x100_corpcor.csv"))
    test_mat2 = readdlm(joinpath(p, "test_matrices/100x20.csv"))
    ref_cov2  = readdlm(joinpath(p, "test_matrices/100x20_corpcor.csv"))
    test_mat3 = readdlm(joinpath(p, "test_matrices/50x50.csv"))
    ref_cov3  = readdlm(joinpath(p, "test_matrices/50x50_corpcor.csv"))
    @test cov(ss, test_mat1) ≈ ref_cov1
    @test cov(ss, test_mat2) ≈ ref_cov2
    @test cov(ss, test_mat3) ≈ ref_cov3
end

@testset "LinShrink: target ABCDE with LW    " begin
    # TARGET A
    lwa = LinearShrinkage(DiagonalUnitVariance())
    for X̂ ∈ test_matrices
        n, p = size(X̂)
        S = cov(SimpleCovariance(), X̂)
        Xtmp = centercols(X̂)
        shrinkage  = sum_var_sij(Xtmp, S, n)
        shrinkage /= sum((S-Diagonal(S)).^2) + sum((diag(S).-1).^2)
        shrinkage = clamp(shrinkage, 0.0, 1.0)
        @test cov(lwa, X̂) ≈ (1.0-shrinkage) * S + shrinkage * I
        lwafixed = LinearShrinkage(DiagonalUnitVariance(), shrinkage)
        @test cov(lwafixed, X̂) ≈ (1.0 - shrinkage) * S + shrinkage * I
    end
    # TARGET B
    lwb = LinearShrinkage(DiagonalCommonVariance())
    for X̂ ∈ test_matrices
        n, p = size(X̂)
        S = cov(SimpleCovariance(), X̂)
        Xtmp = centercols(X̂)
        v = tr(S)/p
        F = v * I
        shrinkage  = sum_var_sij(Xtmp, S, n)
        shrinkage /= sum((S-Diagonal(S)).^2) + sum((diag(S).-v).^2)
        shrinkage = clamp(shrinkage, 0.0, 1.0)
        @test cov(lwb, X̂) ≈ (1.0-shrinkage) * S + shrinkage * F
        lwbfixed = LinearShrinkage(DiagonalCommonVariance(), shrinkage)
        @test cov(lwbfixed, X̂) ≈ (1.0 - shrinkage) * S + shrinkage * F
    end
    # TARGET C
    lwc = LinearShrinkage(CommonCovariance())
    for X̂ ∈ test_matrices
        n, p = size(X̂)
        S = cov(SimpleCovariance(), X̂)
        Xtmp = centercols(X̂)
        v = tr(S)/p
        c = sum(S-Diagonal(S))/(p*(p-1))
        F = v * I + c * (ones(p, p) - I)
        shrinkage  = sum_var_sij(Xtmp, S, n)
        shrinkage /= sum(((S-Diagonal(S)) - c*(ones(p, p)-I)).^2) + sum((diag(S) .- v).^2)
        shrinkage = clamp(shrinkage, 0.0, 1.0)
        @test cov(lwc, X̂) ≈ (1.0-shrinkage) * S + shrinkage * F
        lwcfixed = LinearShrinkage(CommonCovariance(), shrinkage)
        @test cov(lwcfixed, X̂) ≈ (1.0-shrinkage) * S + shrinkage * F
    end
    # TARGET D
    lwd = LinearShrinkage(DiagonalUnequalVariance())
    for X̂ ∈ test_matrices
        n, p = size(X̂)
        S = cov(SimpleCovariance(), X̂)
        Xtmp = centercols(X̂)
        F = Diagonal(S)
        shrinkage  = sum_var_sij(Xtmp, S, n, false; with_diag=false)
        shrinkage /= sum((S-Diagonal(S)).^2)
        shrinkage = clamp(shrinkage, 0.0, 1.0)
        @test cov(lwd, X̂) ≈ (1.0-shrinkage) * S + shrinkage * F
        lwdfixed = LinearShrinkage(DiagonalUnequalVariance(), shrinkage)
        @test cov(lwdfixed, X̂) ≈ (1.0-shrinkage) * S + shrinkage * F
    end
    # TARGET E
    lwe = LinearShrinkage(PerfectPositiveCorrelation())
    for X̂ ∈ test_matrices
        n, p = size(X̂)
        S = cov(SimpleCovariance(), X̂)
        Xtmp = centercols(X̂)
        d = diag(S)
        F = sqrt.(d*d')
        shrinkage  = sum_var_sij(Xtmp, S, n; with_diag=false)
        shrinkage -= CE.sum_fij(Xtmp, S, n, n)
        shrinkage /= sum((S - F).^2)
        shrinkage = clamp(shrinkage, 0.0, 1.0)
        @test cov(lwe, X̂) ≈ (1.0-shrinkage) * S + shrinkage * F
        lwefixed = LinearShrinkage(PerfectPositiveCorrelation(), shrinkage)
        @test cov(lwefixed, X̂) ≈ (1.0-shrinkage) * S + shrinkage * F
    end
end

@testset "LinShrink: target B with RBLW+OAS  " begin
    rblw = LinearShrinkage(DiagonalCommonVariance(), :rblw)
    testTransposition(rblw, X)
    testUncorrelated(rblw)
    testTranslation(rblw, X)
    testDims(rblw, X)

    oas = LinearShrinkage(DiagonalCommonVariance(), :oas)
    testTransposition(oas, X)
    testUncorrelated(oas)
    testTranslation(oas, X)
    testDims(oas, X)

    for X̂ ∈ test_matrices
        Ŝ_rblw = cov(rblw, X̂)
        Ŝ_oas  = cov(oas, X̂)

        X̂ = centercols(X̂)
        n, p = size(X̂)
        Ŝ    = cov(SimpleCovariance(), X̂)

        F_ref = tr(Ŝ)/p * I
        # https://arxiv.org/pdf/0907.4698.pdf eq 17
        λ_rblw_ref = ((n-2)/n*tr(Ŝ^2)+tr(Ŝ)^2)/((n+2)*(tr(Ŝ^2)-tr(Ŝ)^2/p))
        λ_rblw_ref = clamp(λ_rblw_ref, 0.0, 1.0)
        # https://arxiv.org/pdf/0907.4698.pdf eq 23
        λ_oas_ref = ((1-2/p)*tr(Ŝ^2)+tr(Ŝ)^2)/((n+1-2/p)*(tr(Ŝ^2)-tr(Ŝ)^2/p))
        λ_oas_ref = clamp(λ_oas_ref, 0.0, 1.0)

        @test cov(rblw, X̂) ≈ CE.linshrink(F_ref, Ŝ, λ_rblw_ref)
        @test cov(oas, X̂) ≈ CE.linshrink(F_ref, Ŝ, λ_oas_ref)

        rblw_fixed = LinearShrinkage(DiagonalCommonVariance(), λ_rblw_ref)
        oas_fixed = LinearShrinkage(DiagonalCommonVariance(), λ_oas_ref)
        @test cov(rblw_fixed, X̂) ≈ CE.linshrink(F_ref, Ŝ, λ_rblw_ref)
        @test cov(oas_fixed, X̂) ≈ CE.linshrink(F_ref, Ŝ, λ_oas_ref)
    end
end

@testset "LinShrink: all targets, std (SS)   " begin
    for target ∈ [
            DiagonalUnitVariance(),
            DiagonalCommonVariance(),
            DiagonalUnequalVariance(),
            CommonCovariance(),
            PerfectPositiveCorrelation(),
            ConstantCorrelation()
            ]
        #=
        :ss assumes that the shrinkage should be the same if the data
        matrix is standardised. This is what is tested.
        =#
        for c ∈ [false, true], X̂ ∈ test_matrices
            Xcs = X̂
            Xcs = centercols(Xcs)
            n = size(Xcs, 1)
            for i ∈ 1:size(Xcs, 2)
                Xcs[:, i] ./= sqrt(var(Xcs[:, i], corrected=c))
            end
            LSEss = LinearShrinkage(target, :ss, corrected=c)
            LSElw = LinearShrinkage(target, :lw, corrected=c)
            @test cov(LSEss, Xcs) ≈ cov(LSElw, Xcs)
            # weights are not currently exported but it seems to be a good
            # idea to ensure proper errors
            fw1 = FrequencyWeights(rand(1:10, size(Xcs, 1)))
            @test_throws ErrorException cov(LSEss, Xcs, fw1)
            @test_throws ErrorException cov(LSEss, Xcs, fw1; dims=2)
            @test_throws ErrorException cov(LSEss, Xcs, fw1; mean=nothing)
        end
    end
end

@testset "LinShrink: all targets, mean arg   " begin
    for target ∈ [
            DiagonalUnitVariance(),
            DiagonalCommonVariance(),
            DiagonalUnequalVariance(),
            CommonCovariance(),
            PerfectPositiveCorrelation(),
            ConstantCorrelation()
            ]
        for c ∈ [true, false], dim ∈ [1, 2], X̂ ∈ test_matrices[1:3]
            Xcs = X̂
            Xcs = centercols(Xcs)
            n = size(Xcs, dim)
            LSEss = LinearShrinkage(target, :ss, corrected=c)
            LSElw = LinearShrinkage(target, :lw, corrected=c)
            m = mean(cov(LSEss, Xcs, dims=dim, mean=nothing))
            @test cov(LSEss, Xcs, dims=dim, mean=mean(Xcs, dims=dim)) ≈ cov(LSEss, Xcs, dims=dim, mean=nothing)
            @test cov(LSElw, Xcs, dims=dim, mean=mean(Xcs, dims=dim)) ≈ cov(LSElw, Xcs, dims=dim, mean=nothing)
            if dim == 2
                meanvec = vec(mean(Xcs, dims=dim))
                @test cov(LSEss, Xcs, dims=dim, mean=meanvec) ≈ cov(LSEss, Xcs, dims=dim, mean=nothing)
                @test cov(LSElw, Xcs, dims=dim, mean=meanvec) ≈ cov(LSElw, Xcs, dims=dim, mean=nothing)
            end
        end
    end
end


text:
# We wish to extend operations to identically named methods dispatched
# on `Machine`s. For example, we have from the model API
#
# `predict(model::M, fitresult, X) where M<:Supervised`
#
# but want also want to define
#
# 1. `predict(machine::Machine, X)` where `X` is concrete data
#
# and we would like the syntactic sugar (for `X` a node):
#
# 2. `predict(machine::Machine, X::Node) = node(predict, machine, X)`
#
# Finally, for a `model` that is `ProbabilisticComposite`,
# `DetermisiticComposite`, or `UnsupervisedComposite`, we want
#
# 3. `predict(model, fitresult, X) = fitresult.predict(X)`
#
# which makes sense because `fitresult` in those cases is a named
# tuple keyed on supported operations and with nodes as values.

## TODO: need to add checks on the arguments of
## predict(::Machine, ) and transform(::Machine, )

const OPERATIONS = (:predict, :predict_mean, :predict_mode, :predict_median, :predict_joint,
                    :transform, :inverse_transform)

for operation in OPERATIONS

    if operation != :inverse_transform

        ex = quote
            # 0. operations on machs, given empty data:
            function $(operation)(mach::Machine; rows=:)
                # Base.depwarn("`$($operation)(mach)` and "*
                #              "`$($operation)(mach, rows=...)` are "*
                #              "deprecated. Data or nodes "*
                #              "should be explictly specified, "*
                #              "as in `$($operation)(mach, X)`. ",
                #              Base.Core.Typeof($operation).name.mt.name)
                if isempty(mach.args) # deserialized machine with no data
                    throw(ArgumentError("Calling $($operation) on a "*
                                        "deserialized machine with no data "*
                                        "bound to it. "))
                end
                return ($operation)(mach, mach.args[1](rows=rows))
            end
        end
        eval(ex)

    end
end

for operation in (:inverse_transform,)
    ex = quote
        # 0. operations on machines, given empty data:
        $operation(mach::Machine; rows=:) =
            throw(ArgumentError("`$($operation)(mach)` and "*
                                "`$($operation)(mach, rows=...)` is "*
                                "not supported. Data or nodes "*
                                "must be explictly specified, "*
                                "as in `$($operation)(mach, X)`. "))
    end
    eval(ex)
end

_symbol(f) = Base.Core.Typeof(f).name.mt.name

for operation in OPERATIONS

    ex = quote
        # 1. operations on machines, given *concrete* data:
        function $operation(mach::Machine, Xraw)
            if mach.state > 0
                return $(operation)(mach.model, mach.fitresult,
                                    Xraw)
            else
                error("$mach has not been trained.")
            end
        end

        function $operation(mach::Machine{<:Static}, Xraw, Xraw_more...)
            isdefined(mach, :fitresult) || (mach.fitresult = nothing)
            return $(operation)(mach.model, mach.fitresult,
                                    Xraw, Xraw_more...)
        end

        # 2. operations on machines, given *dynamic* data (nodes):
        $operation(mach::Machine, X::AbstractNode) =
            node($(operation), mach, X)

        $operation(mach::Machine{<:Static}, X::AbstractNode, Xmore::AbstractNode...) =
            node($(operation), mach, X, Xmore...)
    end
    eval(ex)
end

## SURROGATE AND COMPOSITE MODELS

for operation in [:predict, :predict_joint, :transform, :inverse_transform]
    ex = quote
        $operation(model::Union{Composite,Surrogate}, fitresult,X) =
            fitresult.$operation(X)
    end
    eval(ex)
end

for (operation, fallback) in [(:predict_mode, :mode), (:predict_mean, :mean), (:predict_median, :median)]
    ex = quote
        function $(operation)(m::Union{ProbabilisticComposite,ProbabilisticSurrogate},
                                    fitresult,
                                    Xnew)
        if haskey(fitresult, $(QuoteNode(operation)))
                return fitresult.$(operation)(Xnew)
            end
            return $(fallback).(predict(m, fitresult, Xnew))
        end
    end
    eval(ex)
end


text:
# https://juliamolsim.github.io/DFTK.jl/stable/guide/tutorial/

using ASEPotential
using DFTK
using Unitful
using UnitfulAtomic

lattice = austrip(5.431u"Å") / 2 * [[0 1 1.]; [1 0 1.]; [1 1 0.]]
Si = ElementCoulomb(:Si)
dftk_atoms = [Si => [ones(3)/8, -ones(3)/8]]
atoms = ase_atoms(lattice, dftk_atoms)
parameters = DFTKCalculatorParameters(
    ecut=7u"hartree",
    kpts=[4,4,4],
    scftol=1e-8
)

@time println(get_potential_energy(atoms, parameters))
@time println(get_forces(atoms, parameters))


text:
#=
# AST of a Julia function definition

In this example we are going to plot an abstract syntax tree of a Julia function using the
Bucheim Layout from [`NetworkLayout.jl`](https://github.com/JuliaGraphs/NetworkLayout.jl).
=#
using CairoMakie
CairoMakie.activate!(type="png") #hide
CairoMakie.inline!(true) #hide
set_theme!(resolution=(800, 600)) #hide
using LightGraphs
using GraphMakie
using NetworkLayout
using CairoMakie

#=
The following code, which walks the AST and creats a `SimpleDiGraph` was taken and slightly
modified from [`TreeView.jl`](https://github.com/JuliaTeX/TreeView.jl). Thanks!
=#
function walk_tree(ex; show_call=true)
    g = SimpleDiGraph()
    labels = Any[]
    walk_tree!(g, labels, ex, show_call)
    return (g, labels)
end

function walk_tree!(g, labels, ex, show_call)
    add_vertex!(g)
    top_vertex = vertices(g)[end]

    where_start = 1  # which argument to start with

    if !(show_call) && ex.head == :call
        f = ex.args[1]   # the function name
        push!(labels, f)
        where_start = 2   # drop "call" from tree
    else
        push!(labels, ex.head)
    end

    for i in where_start:length(ex.args)
        if isa(ex.args[i], Expr)
            child = walk_tree!(g, labels, ex.args[i], show_call)
            add_edge!(g, top_vertex, child)
        elseif !isa(ex.args[i], LineNumberNode)
            add_vertex!(g)
            n = vertices(g)[end]
            add_edge!(g, top_vertex, n)
            push!(labels, ex.args[i])
        end
    end

    return top_vertex
end
nothing #hide

#=
The expression we want to look at is the recursive definition of the Fibonacci sequence.
=#
expr = quote
    function fib(n)
        if n > 1
            return fib(n-1) + fib(n-2)
        else
            return n
        end
    end
end

g, labels = walk_tree(expr, show_call=true)
nlabels_align = [(:left, :bottom) for v in vertices(g)]
fig, ax, p = graphplot(g; layout=Buchheim(),
                       nlabels=repr.(labels),
                       nlabels_distance=5,
                       nlabels_align)
hidedecorations!(ax); hidespines!(ax)
fig #hide

# This does not look nice yet! Lets tweak the `align` parameter of the nodes labels...
for v in vertices(g)
    if isempty(inneighbors(g, v)) # root
        nlabels_align[v] = (:center,:bottom)
    elseif isempty(outneighbors(g, v)) #leaf
        nlabels_align[v] = (:center,:top)
    else
        self = p[:node_positions][][v]
        parent = p[:node_positions][][inneighbors(g, v)[1]]
        if self[1] < parent[1] # left branch
            nlabels_align[v] = (:right,:bottom)
        end
    end
end
p.nlabels_align = nlabels_align
fig # hide


text:
module ShaderManager #GLShaderManager

export ShaderProperty
export ShaderProgram
export Shader

using CoreExtended
using ModernGL
using FileManager

using JLGEngine.GraphicsManager
using ..GLLists
using ..GLDebugControl
using ..GLExtendedFunctions

""" TODO """
type ShaderProperty <: IGraphicsShaderProperty
	enabled::Bool
	linked::Bool

	name::String
	typname::String
	category::String

	typ::DataType
	elements::Int32

	categoryid::Int32
	typid::Int32
	location::Int32 #GLuint
	offset::Int32

	arrayStride::Int32
	matrixStride::Int32
	topSize::Int32
	topStride::Int32

	size::UInt32
	code::UInt32

	index::Int32
	binding::Int32
	id::UInt32

	data::Any
	update::Function

	properties::Dict{String, ShaderProperty}

	ShaderProperty() = new(true,false,"","","",Void,0,-1,-1,-1,0,-1,-1,-1,-1,0,0,-1,-1,0,nothing,()->nothing,Dict{String,ShaderProperty}())
end

""" TODO """
type Shader <: IGraphicsShader
	enabled ::Bool
	created ::Bool

	typ			::Symbol
	file 		::FileSourcePart
	cache		::String

	id      ::GLuint # GL
	typid   ::GLuint # GL

	Shader(typ::Symbol, source::FileSource) = Shader(typ, FileSourcePart(source))
  Shader(typ::Symbol, file::FileSourcePart)	= new(true,false,typ,file,"",0,0)
end

""" TODO """
type ShaderProgram <: IGraphicsShaderProgram
	id			::GLuint

	enabled	::Bool
	created ::Bool
	linked	::Bool
	bound		::Bool

	#source	::FileSource
	shaders ::Dict{Symbol,Shader}
	properties ::Dict{String, ShaderProperty}

	#ShaderProgram() = new(0,true,false,false,false,FileSource(),Dict(),Dict())
  ShaderProgram() = new(0,true,false,false,false,Dict(),Dict())
end

#listShaderProgramListener = createSortedDict(Dict{String,Function})
#listenOnShaderProgramUpdate(id::Symbol, t::Tuple) = listenOnShaderProgramUpdate(id, Dict{String,Function}(t))
#listenOnShaderProgramUpdate(id::Symbol, d::Dict{String,Function}) = listShaderProgramListener[id]=d
#notlistenOnShaderPropertyUpdate(id::Symbol) = (listShaderProgramListener[id]=nothing)

""" TODO """
link(this::Shader, p::AbstractGraphicsShaderProgram) = (this.program = p; true)

""" TODO """
function compile(this::Shader)
	println("compile: ", this.typ);
	f = this.file
	glGetShaderSource(this.id)
	cache = f.cache #source.cache #FileManager.getCache(f.source)

	if cache == ""
		warn("OpenGL Error! Compile shader '$(f.source.path)': Failed loading content!")
		return false
	end

	#cache = cache[f.range]
	#println(f.source.path, " ", f.range, ":(",length(cache), ")\n--------------------------\n", cache,"\n--------------------------\n")

	glCompile(this.id,cache)
	r = validate(:SHADER, this.id, GL_COMPILE_STATUS)
	if !r
		warn("OpenGL Error! Compile shader '$(f.source.path)': $(getInfoLog(:SHADER, this.id))")
		glCompile(this.id,this.cache)
		r = validate(:SHADER, this.id, GL_COMPILE_STATUS)
		if !r
			warn("OpenGL Error! Compile shader cache '$(f.source.path)': $(getInfoLog(:SHADER, this.id))")
		end
	else
		this.cache = cache
	end
	r
end

""" TODO """
function reload(this::AbstractGraphicsShaderProgram, source::FileSource)
	if this == nothing return end

	clear(this)

	if printError("Before readShadersFromSource") end

	println("[ readShadersFromSource ]")
	readShaders(this,source)
	println("---------------------------------------")
	println("[ setUp ]")
	setUp(this)
	println("---------------------------------------")

	#programID = this.id
	#getProgramInfo(this)
	#getAttributesInfo(this)
	#getUniformsInfo(this)

	if isValid(this)
		#glUseProgram(this.id)
		println("[ findActiveRessources ]")
		findActiveRessources(this)
		println("---------------------------------------")
		println("[ findShaderProperties ]")
		findShaderProperties(this)
		println("---------------------------------------")
	end
	link(this)
	#useLinkedShaderProgram() # switch back to linked shader
end

""" TODO """
createShaderProgram() = ShaderProgram()

""" TODO """
createShader(program::ShaderProgram, typ::Symbol, file::FileSourcePart) = (this=Shader(typ, file); program.shaders[typ]=this; this)

#getSource(this::AbstractGraphicsShaderProgram) = this.source
#setPath(this::AbstractGraphicsShaderProgram, path::String) = (this.source.path = path)

""" TODO """
clear(this::AbstractGraphicsShaderProgram) = this.shaders = Dict()

""" TODO """
hasShader(this::AbstractGraphicsShaderProgram, id::Symbol) = haskey(this.shaders,id)

linked = nothing

""" TODO """
unlink() = link(nothing)
#isValid(this::ShaderProgram) = this.created && this.bound

""" TODO """
isValid(this::AbstractGraphicsShaderProgram) = this != nothing && isa(this,ShaderProgram) && this.created && this.bound

""" TODO """
isLinked(this::AbstractGraphicsShaderProgram) =	isValid(this) && this == linked

""" TODO """
function create(this::Shader)
	if this.created return true end
	this.created = true
	println("create: ", this.typ);
	this.typid = LIST_SHADER[this.typ][:TYPE]
	this.id = glCreateShader(this.typid)
	hasNoError()
end

""" TODO """
function delete(this::Shader)
	if !this.created return true end
	println("delete: ", this.typ);
	glDeleteShader(this.id)
	this.created = false
	hasNoError()
end

""" TODO """
attach(this::Shader, programID::GLuint) = (println("attach: ", this.typ); glAttachShader(programID, this.id); hasNoError())

""" TODO """
detach(this::Shader, programID::GLuint) = (println("detach: ", this.typ); glDetachShader(programID, this.id); hasNoError())

""" TODO """
function link(this::AbstractGraphicsShaderProgram)
	global linked
	if linked == this return end
	linked = this
	result = isValid(this)
	id = 0
	if result	id = this.id end
	glUseProgram(id)
	if result	result=hasNoError() end
	result
end

""" TODO """
function create(this::ShaderProgram)
	if this.created return true end
	this.created = true
	this.id = glCreateProgram()
	hasNoError()
end

""" TODO """
function delete(this::ShaderProgram)
	if !this.created return true end
	this.created = false
	unlink()
	glDeleteProgram(this.id)
	hasNoError()
end

""" TODO """
function bind(this::ShaderProgram)
	glLinkProgram(this.id)
  #result = validate(:PROGRAM, this.id, GL_LINK_STATUS)
	#if !result	warn("OpenGL Error! Link ShaderProgram $(this.id): $(getInfoLog(:PROGRAM, this.id))")	end
	#result
	#this.bound=result
	#result
	this.bound=hasNoError()
end

""" TODO """
function setUp(this::ShaderProgram)
	result = (
	(reset(this) ? true : true) &&
	(create(this) && detach(this) && compile(this) && attach(this) && bind(this) && clearList(this))
	) ? true : (reset(this) ? false : false)
	#if !result warn("OpenGL Error! SetUp ShaderProgram $(this.id): $(getInfoLog(:PROGRAM, this.id))") end
	result
end

""" TODO """
function invokeList(list::Dict, f::Function, args...)
	if length(list) <= 0 return false end;
	for (k,e) in list
		if !f(e,args...) return false end
	end
	true
end

""" TODO """
reset() = unlink()
#(create(this) && compileAll(this) && link(this, true) && cleanAll(this))

""" TODO """
reset(this::ShaderProgram) = (unlink(); clearList(this) && delete(this))

""" TODO """
compile(this::ShaderProgram) = invokeList(this.shaders, (s)->(create(s) && compile(s)))

""" TODO """
detach(this::ShaderProgram) = (shader_ids = glGetAttachedShaders(this.id); foreach(glDetachShader, shader_ids); true)

""" TODO """
attach(this::ShaderProgram) = invokeList(this.shaders, attach, this.id)

""" TODO """
clearList(this::ShaderProgram) = invokeList(this.shaders, delete)

""" TODO """
clear(this::ShaderProgram) = (this.properties = Dict{String, ShaderProperty}())

""" TODO """
function render(this::ShaderProgram, buffer, data)
	if !isLinked(this) return end
	if !buffer.linked return end

	mode = data.mode
	if hasShader(this,:TESS_CONTROL) ||  hasShader(this,:TESS_EVALUATION)
		mode=GL_PATCHES
	end

	if buffer.typ == GL_ELEMENT_ARRAY_BUFFER
		glDrawElements(mode, data.elementsCount, data.typ, C_NULL)
	else
		glDrawArrays(mode, 0, data.elementsCount) #first::GLint
	end
end

include("GLSLParser.jl")

end #GLShaderManager

max_stars_repo_name: Gilga/JGE

max_issues_repo_name: Gilga/JGE

max_forks_repo_name: Gilga/JGE


text:
using Base.Test

include("difference-of-squares.jl")

@testset "Square the sum of the numbers up to the given number" begin
    @test square_of_sum(5) == 225
    @test square_of_sum(10) == 3025
    @test square_of_sum(100) == 25502500
end

@testset "Sum the squares of the numbers up to the given number" begin
    @test sum_of_squares(5) == 55
    @test sum_of_squares(10) == 385
    @test sum_of_squares(100) == 338350
end

@testset "Subtract sum of squares from square of sums" begin
    @test difference(0) == 0
    @test difference(5) == 170
    @test difference(10) == 2640
    @test difference(100) == 25164150
end


text:
abstract type AbstractBinaryTree end
"""
array-based binary tree
"""
abstract type AbstractBinaryTree_arr <: AbstractBinaryTree end
mutable struct BinaryTree{T} <: AbstractBinaryTree_arr where {T}
    n::Int
    root::Int
    par::Vector{Int}
    lch::Vector{Int}
    rch::Vector{Int}
    val::Vector{T}
end
function BinaryTree{T}(size::Int, rootval::T) where {T}
    x = BinaryTree(
        1,
        1,
        zeros(T, size),
        zeros(T, size),
        zeros(T, size),
        Vector{T}(undef, size),
    )
    x.val[1] = rootval
    return x
end

function Base.insert!(tree::BinaryTree, from::Int, value, place::Bool)
    tree.n += 1
    tree.par[tree.n] = from
    if place
        tree.rch[from] = tree.n
    else
        tree.lch[from] = tree.n
    end
    return tree.val[tree.n] = value
end

function Base.insert!(tree::BinaryTree, from::Int, value)
    if tree.lch[from] == 0
        insert!(tree, from, value, false)
    elseif tree.rch[from] == 0
        insert!(tree, from, value, true)
    else
        return false
    end
    return true
end

# for AbstractBinaryTree_arr
function isleaf(
    tree::Tree,
    from::Int,
)::Bool where {Tree<:AbstractBinaryTree_arr}
    return tree.lch[from] == 0 && tree.rch[from] == 0
end

function ch(
    tree::Tree,
    from::Int,
    select::Bool,
) where {Tree<:AbstractBinaryTree_arr}
    if select
        return tree.rch[from]
    else
        return tree.lch[from]
    end
end

function set_ch(
    tree::Tree,
    from::Int,
    select::Bool,
    value::T,
) where {Tree<:AbstractBinaryTree_arr,T}
    if select
        tree.rch[from] = value
    else
        tree.lch[from] = value
    end
end

function height(
    tree::Tree,
    from::Int = tree.root,
) where {Tree<:AbstractBinaryTree_arr}
    count = 0
    if tree.lch[from] != 0
        count = height(tree, tree.lch[from]) + 1
    end
    if tree.rch[from] != 0
        count = max(height(tree, tree.rch[from]) + 1, count)
    end
    return count
end

function depth(tree::Tree, from::Int) where {Tree<:AbstractBinaryTree_arr}
    count = 0
    while from != tree.root
        count += 1
        from = tree.par[from]
    end
    return count
end

function left(
    tree::Tree,
    from::Int = tree.root,
) where {Tree<:AbstractBinaryTree_arr}
    if tree.lch[from] == 0
        return from
    else
        return left(tree, tree.lch[from])
    end
end

function right(
    tree::Tree,
    from::Int = tree.root,
) where {Tree<:AbstractBinaryTree_arr}
    if tree.rch[from] == 0
        return from
    else
        return right(tree, tree.rch[from])
    end
end

max_issues_repo_name: AugustoCL/Julia

max_forks_repo_name: AugustoCL/Julia


text:
"""
Genetic algorithm is a search algorithm, which a computational
metaphor of the natural evolution. Genetic algorithm is useful to solve
a wide range of problems that are hard to solve using brute-force.
"""
module GeneticAlgorithm

export crossoverAtIndex, pickBestIndividual, selectIndividual
export runGA
export mutateAtIndex
export gaLog

using Random

"""
    createIndividual()

Create an individual
"""
function createIndividual(createGene, numberOfGenes)
    return [createGene() for i in 1:numberOfGenes]
end

"""
    createPopulation()

Create population
"""
function createPopulation(createGene, numberOfGenes)
    return [createIndividual(createGene, numberOfGenes) for i in 1:10]
end

"""
    crossover(ind1, ind2)

Genetic operation: Crossover
"""
function crossover(ind1, ind2)
    size = length(ind1)
    randomIndex = floor(Int64, rand() * size)
    return crossoverAtIndex(ind1, ind2, randomIndex)
end

"""
    crossover(ind1, ind2, index)

Perform a crossover at a particular index
"""
function crossoverAtIndex(ind1, ind2, index)
    return vcat(ind1[1:index], ind2[index + 1:length(ind2)])
end

"""
    mutate(ind)

Perform a mutation operation on the provided individual
"""
function mutate(createGene, ind)
    size = length(ind)
    randomIndex = floor(Int64, rand() * size + 1)
    return mutateAtIndex(createGene, ind, randomIndex)
end

"""
    mutateAtIndex(ind, index)

Perform a mutation operation on the provided individual, at a provided index
"""
function mutateAtIndex(createGene, ind, index)
    t = copy(ind)
    t[index] = createGene()
    return t
end

"""
    bestFitnessOf(fitness, population)

Return the maximum fitness for a given population of individuals
"""
function bestFitnessOf(fitness, population)
    return maximum(map(fitness, population))
end

"""
    pickBestIndividual(population)

Pick best individual of a given population.
"""
function pickBestIndividual(fitness, population)
    indexMaxFitness = argmin(map(fitness, population))
    bestIndividual = population[indexMaxFitness]
    return bestIndividual
end

"""
    selectIndividual(population; k=5)

Selection using the tournament of size k. Assume that the smallest fitness is the best.
"""
function selectIndividual(fitness, population; k=5)
    populationSize = length(population)
    selectedIndividuals = []
    for tmp in 1:k
        anInd = population[floor(Int64, rand()*populationSize) + 1]
        push!(selectedIndividuals, anInd)
    end
    return pickBestIndividual(fitness, selectedIndividuals)
end

function gaLog(aString::String, shouldLog=true)
    if(shouldLog)
        printstyled(aString, color=:blue)
    end
end

"""
    run(;maxNumberOfIterations=10, probMutation=0.2, seed=42)

This is the main function of the genetic algorithm.

# Examples

```julia
julia> runGA(maxNumberOfIterations=40)
```
"""
function runGA(fitness, createGene, numberOfGenes; maxNumberOfIterations=10, probMutation=0.2, seed=42, logging=true)
    Random.seed!(seed)
    population = createPopulation(createGene, numberOfGenes)
    numberOfIndividuals = length(population)
    fitnesses = []
    for iteration in 1:maxNumberOfIterations
        gaLog("Begining of iteration = $(iteration)/$(maxNumberOfIterations) ... ", logging)
        newPopulation = []
        for it in 1:numberOfIndividuals
            ind1 = selectIndividual(fitness, population)
            ind2 = selectIndividual(fitness, population)
            newIndividual = crossover(ind1, ind2)
            if (rand() <= probMutation)
                newIndividual = mutate(createGene, newIndividual)
            end
            push!(newPopulation, newIndividual)
        end
        population = newPopulation
        bestFitness = bestFitnessOf(fitness, population)
        push!(fitnesses, bestFitness)
        gaLog("end (best fitness = $(bestFitness))\n", logging)
    end
    return pickBestIndividual(fitness, population), fitnesses
end

"""
    plusTwo(x)

Sum the numeric "2" to whatever it receives as input

A more detailed explanation can go here, although I guess it is not needed in this case

# Arguments
* `x`: The amount to which we want to add 2

# Notes
* Notes can go here

# Examples
```julia
julia> five = plusTwo(3)
5
```
"""
plusTwo(x) = return x + 2

end # module


text:
## Demonstrates use of BAMF on a single region with fixed Number of emitters

# !this example may be broken

using Revise
using Plots
println(pwd())
#include("../src/RJMCMC.jl")
include("../src/BayesMultiFit.jl")
using .RJMCMC

## create a dataset
n=Int32(3)
sz=Int32(16)
sigma=1.3f0
x=Vector{Float32}([sz/2,sz/2-3,sz/2+3])+randn(Float32,3)
y=Vector{Float32}([sz/2,sz/2+3,sz/2+3])+randn(Float32,3)
photons=Vector{Float32}([777,500,800])
bg=1f-4
datastate=BAMF.StateFlatBg(n,x,y,photons,bg)

##  create a model then make noisy
dataroi=BAMF.ArrayDD(sz)
BAMF.genmodel_2Dgauss!(datastate,sz,sigma,dataroi.data)
noisyroi=BAMF.poissrnd(dataroi)

## create a BAMF-type RJMCMC structure
xystd=Float32(sigma/10)
istd=Float32(10)
myRJ=BAMF.RJStructDD(sz,sigma,xystd,istd,noisyroi)

## setup the RJMCMC.jl model
#Jumptypes are: move, bg, add, remove, split, merge
njumptypes=6
jumpprobability=[1,0,0,0,0,0] #Move only
jumpprobability=jumpprobability/sum(jumpprobability)

# create an RJMCMC structure with all model info
iterations=5000
burnin=100
acceptfuns=[BAMF.accept_move] #array of functions
propfuns=[BAMF.propose_move] #array of functions
myRJMCMC=RJMCMC.RJMCMCStruct(burnin,iterations,njumptypes,jumpprobability,propfuns,acceptfuns)

#create an intial state
#state1=BAMF.calcintialstate(myRJ)
n=Int32(3)
sz=Int32(16)
sigma=1.3f0
x=Vector{Float32}([sz/2,sz/2,sz/2])+randn(Float32,3)
y=Vector{Float32}([sz/2,sz/2,sz/2])+randn(Float32,3)
photons=Vector{Float32}([777,500,800])
bg=1f-4
state1=BAMF.StateFlatBg(n,x,y,photons,bg)

## run chain
@time mychain=RJMCMC.buildchain(myRJMCMC,myRJ,state1)

## Display
zoom=Int32(20)
jts=mychain.jumptypes
plotly()
BAMF.histogram2D(mychain.states,sz,zoom,noisyroi,datastate)


text:
@testset "`request_agenda_items`" begin
    meeting_link = "http://somervillecityma.iqm2.com/Citizens/Detail_Meeting.aspx?ID=3163"
    items = SomervilleCouncilParser.request_agenda_items(meeting_link)
    @test nrow(items) == 4

    @test count(items.is_heading) == 3
    @test names(items) == ["item", "content", "is_heading"]

    @test request_agenda_items(meeting_link) == request_agenda_items(3163)
    @test request_agenda_items("3163") == request_agenda_items(3163)

    # support additional special cases
    meeting_link2 = "http://somervillecityma.iqm2.com/Citizens/Detail_Meeting.aspx?ID=3421"
    items = SomervilleCouncilParser.request_agenda_items(meeting_link2)
    @test nrow(items) == 0
end

@testset "`get_agenda_items`" begin
    meeting_link = "http://somervillecityma.iqm2.com/Citizens/Detail_Meeting.aspx?ID=3163"

    @test get_agenda_items(meeting_link) == get_agenda_items(3163)
    @test get_agenda_items("3163") == get_agenda_items(3163)

    # Test caching
    cache_dir = mktempdir()
    @test length(readdir(cache_dir)) == 0
    items = get_agenda_items(meeting_link; cache_dir=nothing)
    @test length(readdir(cache_dir)) == 0
    items = get_agenda_items(meeting_link; cache_dir)
    @test length(readdir(cache_dir)) == 1

    # Test that current caching hasn't changed
    # If this test _fails_, need to (a) bump `agenda_version` and (b) resave current
    # version of cached asset:
    # get_agenda_items(meeting_link; cache_dir=TEST_ASSETS) # Uncomment to save new version
    test_agenda_path = SomervilleCouncilParser.agenda_cache_path(TEST_ASSETS, 3163)
    @test isfile(test_agenda_path)
    previous_items = DataFrame(Legolas.read(test_agenda_path))
    @test previous_items == items

    # When cached item exists, load it!
    # ...to test, override table saved at cache path and make sure _that's_ what loads
    junk = DataFrame([Agenda(; item=3, content="foo")])
    p = only(readdir(joinpath(cache_dir, "v$(agenda_version)"); join=true))
    Legolas.write(p, junk, Legolas.Schema("agenda@1"))
    cached_items = get_agenda_items(meeting_link; cache_dir)
    @test cached_items == junk
end

@testset "`agenda_contains`" begin
    meeting_link = "http://somervillecityma.iqm2.com/Citizens/Detail_Meeting.aspx?ID=3163"
    agenda = get_agenda_items(meeting_link)
    @test nrow(filter_agenda(agenda, ["Director", "hawk"])) == 1

    @test isempty(filter_agenda(agenda, ["hawk"]))

    @test isempty(filter_agenda(agenda, ["dIrEcToR"]; case_invariant=false))
    @test nrow(filter_agenda(agenda, ["dIrEcToR"]; case_invariant=true)) == 1
end

@testset "`search_agendas_for_content`" begin
    cache_dir = mktempdir()
    x = search_agendas_for_content("9/1/2021", "10/1/2021", ["dpw"]; cache_dir)
    @test nrow(x.items) == 4
    @test nrow(x.meetings) == 13
    @test all(isnothing.(x.meetings.failed_parsing))
    @test isnothing(display_items_by_meeting(x.items))
end