text:
Random.seed!(134)
n = 7
p = 5
x = randn(n)
y = randn(n)
θ = randn(p)

δ  = x .- y
δ1 = norm(δ, 1)
δ2 = norm(δ, 2)^2 / 2
δ3 = norm(δ, 3)^3 / 3
θ1 = norm(θ, 1)
θ2 = norm(θ, 2)^2 / 2
θ3 = norm(θ, 3)^3 / 3
y2 = (sign.(randn(n)) .+ 1) ./ 2

@testset "LPCost basics" begin
    noloss = NoLoss()
    nopenalty = NoPenalty()

    @test noloss(x, y) == 0
    @test nopenalty(θ) == 0

    l1 = L1Loss()
    l2 = L2Loss()
    l3 = LPLoss{3}()
    p1 = L1Penalty()
    p2 = L2Penalty()
    p3 = LPPenalty{3}()

    @test l1(x, y) ≈ δ1
    @test l2(x, y) ≈ δ2
    @test l3(x, y) ≈ δ3
    @test p1(θ) ≈ θ1
    @test p2(θ) ≈ θ2
    @test p3(θ) ≈ θ3

    @test R.getp(p3) == 3
    @test R.lp(θ, 2) ≈ θ2
    @test_throws DomainError R.lp(θ, -1)
end

@testset "Comp LPs" begin
    l1 = L1Loss()
    l2 = L2Loss()
    p1 = L1Penalty()
    p2 = L2Penalty()

    lc = l2 - l1 + NoLoss()
    pc = p1 + p2 - NoPenalty()

    @test (2l2)(x, y) ≈ 2δ2
    @test (l2+l1)(x, y) ≈ δ2 + δ1
    @test (l2+2l1)(x, y) ≈ δ2 + 2δ1
    @test (l2-l1)(x, y) ≈ δ2 - δ1
    @test (l2+2l2-l1)(x, y) ≈ 3δ2-δ1
    @test (lc + lc)(x, y) ≈ 2(δ2 - δ1)
    @test (l2 * 2)(x, y) ≈ 2δ2
    @test 2 * (2 * l2)(x, y) ≈ 4δ2
    @test ((l1 + 2l2) + 2l1)(x, y) ≈ 3δ1 + 2δ2
    @test (l1/3)(x, y) ≈ δ1/3
    @test ((2l1)/3)(x, y) ≈ 2δ1/3
    @test ((l1+l2)/3)(x, y) ≈ (δ1+δ2)/3
    @test ((2l1)+l2)(x, y) ≈ 2δ1+δ2
    @test ((2l1)-3(l1+l2))(x, y) ≈ -δ1 - 3δ2
    @test ((2l1 + 3l2) - (l1 + 2l2))(x, y) ≈ δ1 + δ2

    @test (2p1)(θ) ≈ 2θ1
    @test (p1+p2)(θ) ≈ θ1 + θ2
    @test (p1+2p2)(θ) ≈ θ1 + 2θ2
    @test (p2-p1+3p2)(θ) ≈ θ2 - θ1 + 3θ2
    @test (pc + pc)(θ) ≈ 2(θ1 + θ2)
    @test (p1 * 2)(θ) ≈ 2θ1
    @test 2 * (p1 * 3)(θ) ≈ 6θ1
    @test ((6p1 + 3p2) - 5p1)(θ) ≈ θ1 + 3θ2
    @test (p1/3)(θ) ≈ θ1/3
    @test ((2p1)/4)(θ) ≈ θ1/2
    @test ((2p1 - 5p2)/3)(θ) ≈ (2/3)θ1 - (5/3)θ2
    @test ((2p1) + (p1 + 2p2))(θ) ≈ 3θ1 + 2θ2
    @test ((2p1) - (p1 + p2))(θ) ≈ θ1 - θ2
    @test ((p1 * 2) + 2(p1 - p2))(θ) ≈ 4θ1 - 2θ2
    @test ((4p1 + p2) - (3p2 + p1))(θ) ≈ 3θ1 - 2θ2
end

@testset "Compose" begin
    nl = NoLoss()
    np = NoPenalty()

    @test (nl * 5) isa NoLoss
    @test (nl * 5)(x, y) == 0.0 == (5 * nl)(x, y)
    @test (np * pi) isa NoPenalty
    @test (np * pi)(θ) == 0.0 == (pi * np)(θ)

    @test (np + L1Penalty())(θ) ≈ (L1Penalty() + np)(θ) ≈ θ1

    sl = 2.0 * L1Penalty()
    @test (3sl)(θ) ≈ 6θ1

    @test (nl + L1Loss())(x, y) ≈ (L1Loss() + nl)(x, y) ≈ δ1

    a = 0.5 * L1Loss() + L2Loss() + L1Loss() + 0.3 * L2Loss()
    @test isa(a, CompositeLoss)
    @test length(a.losses) == 2
    @test isa(a.losses[1], ScaledLoss{L2Loss})
    @test a.losses[1].scale == 1.3
    @test isa(a.losses[2], ScaledLoss{L1Loss})
    @test a.losses[2].scale == 1.5

    b = 0.5 * L1Penalty() + L2Penalty() + L1Penalty() + 0.3 * L2Penalty()
    @test isa(b, CompositePenalty)
    @test length(b.penalties) == 2
    @test isa(b.penalties[1], ScaledPenalty{L2Penalty})
    @test b.penalties[1].scale == 1.3
    @test isa(b.penalties[2], ScaledPenalty{L1Penalty})
    @test b.penalties[2].scale == 1.5
end

@testset "Logistic" begin
    ll = LogisticLoss()

    @test sum(log.(1 .+ exp.(-(x .* y)))) ≈ ll(x, y)
end

@testset "Multinomial" begin
    # comparison with sklearn
    θ = [-0.04843, 0.99519, -0.67237, 1.08812, 0.13362, 0.77136]
    X = [ 0.78843 -0.28336;
         -0.75568  0.22546;
         -0.09012  0.68069;
         -0.34437 -0.98773;
          1.09285 -0.37161 ]
    y = [1, 2, 3, 1, 3]
    mnl = MultinomialLoss()
    P = X * reshape(θ, 2, 3)
    @test mnl(P, y) ≈ 4.994135263476211
    @test mnl(y, P) == mnl(P, y)
end


text:
# Check consistency with SPTK

module SPTKTest

using MelGeneralizedCepstrums
using Base.Test

import SPTK
import MelGeneralizedCepstrums: lpc2c!, b2c, frqtr, frqtr!

function test_mcep(order::Int, α::Float64)
    srand(98765)
    x = rand(512)
    mc = SPTK.mcep(x, order, α)
    mĉ = MelGeneralizedCepstrums._mcep(x, order, α)
    @test_approx_eq mc mĉ
end

function test_mgcep(order::Int, α::Float64, γ::Float64)
    srand(98765)
    x = rand(512)
    mgc = SPTK.mgcep(x, order, α, γ; threshold=0.001)
    mgĉ = MelGeneralizedCepstrums._mgcep(x, order, α, γ; criteria=0.001)
    @test_approx_eq mgc mgĉ
end

function test_lpc2c()
    srand(98765)
    a = rand(21)

    c = SPTK.lpc2c(a)
    ĉ = lpc2c(a)
    @test_approx_eq c ĉ
    ĉ = similar(a)
    lpc2c!(ĉ, a)
    @test_approx_eq c ĉ
end

function test_gnorm(γ::Float64)
    srand(98765)
    mc = rand(21)

    g = SPTK.gnorm(mc, γ)
    ĝ = gnorm(mc, γ)
    @test_approx_eq g ĝ
    ĝ = copy(mc)
    gnorm!(ĝ, γ)
    @test_approx_eq g ĝ
end

function test_ignorm(γ::Float64)
    srand(98765)
    mc = rand(21)

    g = SPTK.ignorm(mc, γ)
    ĝ = ignorm(mc, γ)
    @test_approx_eq g ĝ
    ĝ = copy(mc)
    ignorm!(ĝ, γ)
    @test_approx_eq g ĝ
end

function test_mc2b(α::Float64)
    srand(98765)
    mc = rand(21)

    b = SPTK.mc2b(mc, α)
    b̂ = mc2b(mc, α)
    @test_approx_eq b b̂
    b̂ = copy(mc)
    mc2b!(b̂, α)
    @test_approx_eq b b̂
end

function test_b2mc(α::Float64)
    srand(98765)
    b = rand(21)

    mc = SPTK.b2mc(b, α)
    mĉ = b2mc(b, α)
    @test_approx_eq mc mĉ
end

function test_c2ir(len::Int=512)
    srand(98765)
    c = rand(21)

    ir = SPTK.c2ir(c, len)
    ir̂ = c2ir(c, len)
    @test_approx_eq ir ir̂
end

function test_freqt(order::Int, α::Float64)
    srand(98765)
    mc = rand(21)

    m = SPTK.freqt(mc, order, α)
    m̂ = freqt(mc, order, α)
    @test_approx_eq m m̂
end

function test_b2c(order::Int, α::Float64)
    srand(98765)
    mc = rand(21)

    m = SPTK.b2c(mc, order, α)
    m̂2 = b2c(mc, order, α)
    @test_approx_eq m m̂2
end

function test_frqtr(order::Int, α::Float64)
    srand(98765)
    mc = rand(21)

    m = SPTK.frqtr(mc, order, α)
    m̂ = frqtr(mc, order, α)
    @test_approx_eq m m̂
end

function test_gc2gc(order::Int, γ::Float64)
    srand(98765)
    gc = rand(21)

    gc2 = SPTK.gc2gc(gc, 0.0, order, γ)
    gc2̂ = gc2gc(gc, 0.0, order, γ)
    @test_approx_eq gc2 gc2̂
end

function test_mgc2mgc(order::Int, α::Float64, γ::Float64)
    srand(98765)
    mgc = rand(21)

    mgc2 = SPTK.mgc2mgc(mgc, 0.41, 0.0, order, α, γ)
    mgc2̂ = mgc2mgc(mgc, 0.41, 0.0, order, α, γ)
    @test_approx_eq mgc2 mgc2̂
end

function test_mgc2sp(α::Float64, γ::Float64)
    srand(98765)
    mgc = rand(21)

    fftlen = 512
    sp = SPTK.mgc2sp(mgc, α, γ, fftlen)
    sp̂ = mgc2sp(mgc, α, γ, fftlen)
    @test_approx_eq sp sp̂
    @test length(sp̂) == fftlen>>1 + 1
    @test eltype(sp̂) == Complex{eltype(mgc)}
end

for order in 10:2:30
    for α in [-0.544, -0.41, -0.35, 0.0, 0.35, 0.41, 0.544]
        println("mcep: testing with order=$order, α=$α")
        test_mcep(order, α)
    end
end

for order in 25:5:35
    for α in [-0.544, -0.41, -0.35, 0.0, 0.35, 0.41, 0.544]
        for γ in [-1.0, -0.75, -0.5, -0.25, 0.0]
            println("mgcep: testing with order=$order, α=$α, γ=$γ")
            test_mgcep(order, α, γ)
        end
    end
end

for γ in [-1.0, -0.75, -0.5, -0.25, 0.0]
    println("gnorm: testing with γ=$γ")
    test_gnorm(γ)
end

for γ in [-1.0, -0.75, -0.5, -0.25, 0.0]
    println("ignorm: testing with γ=$γ")
    test_ignorm(γ)
end

for α in [-0.544, -0.41, -0.35, 0.0, 0.35, 0.41, 0.544]
    println("mc2b: testing with α=$α")
    test_mc2b(α)
end

for α in [-0.544, -0.41, -0.35, 0.0, 0.35, 0.41, 0.544]
    println("mc2b: testing with α=$α")
    test_b2mc(α)
end

for len in [128, 256, 512, 1024]
    println("c2ir: testing with len=$len")
    test_c2ir(len)
end

for order in 10:2:30
    for α in [-0.544, -0.41, -0.35, 0.0, 0.35, 0.41, 0.544]
        println("freqt: testing with order=$order, α=$α")
        test_freqt(order, α)
    end
end

for order in 10:2:30
    for α in [-0.544, -0.41, -0.35, 0.0, 0.35, 0.41, 0.544]
        println("b2c: testing with order=$order, α=$α")
        test_b2c(order, α)
    end
end

for order in 10:2:30
    for α in [-0.544, -0.41, -0.35, 0.0, 0.35, 0.41, 0.544]
        println("frqtr: testing with order=$order, α=$α")
        test_frqtr(order, α)
    end
end

for order in 15:5:35
    for γ in [-1.0, -0.75, -0.5, -0.25, 0.0]
        println("gc2gc: testing with order=$order, γ=$γ")
        test_gc2gc(order, γ)
    end
end

for order in 15:5:35
    for α in [-0.544, -0.41, -0.35, 0.35, 0.41, 0.544]
        for γ in [-1.0, -0.75, -0.5, -0.25, 0.0]
            println("mgc2mgc: testing with order=$order, α=$α, γ=$γ")
            test_mgc2mgc(order, α, γ)
        end
    end
end

for α in [-0.544, -0.41, -0.35, 0.35, 0.41, 0.544]
    for γ in [-0.75, -0.5, -0.25, 0.0] # TODO(ryuichi) -1.0
        println("mgc2sp: testing with α=$α and γ=$γ")
        test_mgc2sp(α, γ)
    end
end

println("lpc2c: testing")
test_lpc2c()

end # module SPTKTestModule

ext: jl


text:
function _precompile_()
  ccall(:jl_generating_output, Cint, ()) == 1 || return nothing
  for T in (Bool, Int, Float32, Float64)
    for A in (Vector, Matrix)
      precompile(stridedpointer, (A{T},))
      precompile(stridedpointer, (LinearAlgebra.Adjoint{T,A{T}},))
    end
  end
end


text:
function gamess_to_julia_key(gamess_key)
  gamess_to_julia_mapping = Dict([ 1 => 1,
                                   2 => 2,
                                   3 => 3,
                                   4 => 4,
                                   5 => 5,
                                   6 => 6,
                                   7 => 7,
                                   8 => 8,
                                   9 => 9,
                                   10 => 10,
                                   11 => 13,
                                   12 => 15,
                                   13 => 11,
                                   14 => 12,
                                   15 => 14,
                                   16 => 16,
                                   17 => 17,
                                   18 => 18,
                                   19 => 19,
                                   20 => 20
                                 ])

  ijkl_match = eachmatch(r"[0-9]+", gamess_key)
  ijkl = collect(ijkl_match)
  i = parse(Int,ijkl[1].match)
  j = parse(Int,ijkl[2].match)
  k = parse(Int,ijkl[3].match)
  l = parse(Int,ijkl[4].match)

  i_new = gamess_to_julia_mapping[i]
  j_new = gamess_to_julia_mapping[j]
  k_new = gamess_to_julia_mapping[k]
  l_new = gamess_to_julia_mapping[l]

  if i_new <= j_new
    tmp = i_new
    i_new = j_new
    j_new = tmp
  end

  if k_new <= l_new
    tmp = k_new
    k_new = l_new
    l_new = tmp
  end

  ij = i_new > j_new ? i_new*(i_new-1)/2 + j_new : j_new*(j_new-1)/2 + i_new
  kl = k_new > l_new ? k_new*(k_new-1)/2 + l_new : l_new*(l_new-1)/2 + k_new

  if ij <= kl
    tmp = i_new
    i_new = k_new
    k_new = tmp

    tmp = j_new
    j_new = l_new
    l_new = tmp
  end

  return "$i_new,$j_new,$k_new,$l_new"
  #return "$i,$j,$k,$l"
end

function julia_to_gamess_key(julia_key)
  julia_to_gamess_mapping = Dict([ 1 => 1,
                                   2 => 2,
                                   3 => 3,
                                   4 => 4,
                                   5 => 5,
                                   6 => 6,
                                   7 => 7,
                                   8 => 8,
                                   9 => 9,
                                   10 => 10,
                                   13 => 11,
                                   15 => 12,
                                   11 => 13,
                                   12 => 14,
                                   14 => 15,
                                   16 => 16,
                                   17 => 17,
                                   18 => 18,
                                   19 => 19,
                                   20 => 20
                                 ])

  ijkl_match = eachmatch(r"[0-9]+", julia_key)
  ijkl = collect(ijkl_match)
  i = parse(Int,ijkl[1].match)
  j = parse(Int,ijkl[2].match)
  k = parse(Int,ijkl[3].match)
  l = parse(Int,ijkl[4].match)

  i_new = julia_to_gamess_mapping[i]
  j_new = julia_to_gamess_mapping[j]
  k_new = julia_to_gamess_mapping[k]
  l_new = julia_to_gamess_mapping[l]

  if i_new <= j_new
    tmp = i_new
    i_new = j_new
    j_new = tmp
  end

  if k_new <= l_new
    tmp = k_new
    k_new = l_new
    l_new = tmp
  end

  ij = i_new > j_new ? i_new*(i_new-1)/2 + j_new : j_new*(j_new-1)/2 + i_new
  kl = k_new > l_new ? k_new*(k_new-1)/2 + l_new : l_new*(l_new-1)/2 + k_new

  if ij <= kl
    tmp = i_new
    i_new = k_new
    k_new = tmp

    tmp = j_new
    j_new = l_new
    l_new = tmp
  end

  return "$i_new,$j_new,$k_new,$l_new"
  #return "$i,$j,$k,$l"
end

function get_two_electron_integrals(input_file_string::String)
  #== read in input file ==#
  input_file::Array{String,1} = []
  open(input_file_string) do file
    input_file = readlines(file)
  end

  #== extract integrals from input file ==#
  starting::Int64 = 1
  ending::Int64 = 1

  integrals = Vector{String}([])
  for line in input_file
    if occursin("ERI(",line)
      push!(integrals,line)
    end
  end

  #== set up integral dictionary ==#
  integral_dict = Dict([])
  num_duplicates = 0
  for integral in integrals
    #==get ijkl numbers ==#
    ijkl_match = eachmatch(r"[0-9]+", integral)
    ijkl = collect(ijkl_match)
    i = ijkl[1].match
    j = ijkl[2].match
    k = ijkl[3].match
    l = ijkl[4].match

    #== get actual integral ==#
    integral_match = match(r"0","0")
    if typeof(match(r"(-?[\d]{1,}\.[\d]{1,}[e()E][+()-][\d]{1,})",integral)) != Nothing
      integral_match = match(r"(-?[\d]{1,}\.[\d]{1,}[e()E][+()-][\d]{1,})",integral)
    elseif typeof(match(r"(-\.[\d]{1,})", integral)) != Nothing
      integral_match = match(r"(-\.[\d]{1,})",integral)
    else
      integral_match = match(r"(-?[\d]\.[\d]{1,})",integral)
    end

    integral_value = parse(Float64,integral_match.match)
    if haskey(integral_dict,"$i,$j,$k,$l")
      println("Duplicate Julia key at: "*"$i,$j,$k,$l")
      num_duplicates += 1
    else
      integral_dict["$i,$j,$k,$l"] = integral_value
    end
  end

  println("Number of duplicate keys: ", num_duplicates)

  return integral_dict
end

julia_dict = get_two_electron_integrals(ARGS[1])
gamess_dict = get_two_electron_integrals(ARGS[2])

println(length(keys(julia_dict)))
println(length(keys(gamess_dict)))

for gamess_key in keys(gamess_dict)
  julia_key = gamess_to_julia_key(gamess_key)
  has_julia_key = haskey(julia_dict,julia_key)
  if !has_julia_key
    println("Missing Julia key:", julia_key)
  else
    same = gamess_dict[gamess_key] < 1.0E-10 && julia_dict[julia_key] < 1.0E-10 ? true :
      isapprox(gamess_dict[gamess_key], julia_dict[julia_key], atol = 1e-9)
    if !same
      println("Not same!: $gamess_key => $julia_key: ", gamess_dict[gamess_key], ", ", julia_dict[julia_key])
    end
  end
end

for julia_key in keys(julia_dict)
  gamess_key = julia_to_gamess_key(julia_key)
  has_gamess_key = haskey(gamess_dict,gamess_key)
  if !has_gamess_key
    println("Missing GAMESS key:", gamess_key)
  else
    same = gamess_dict[gamess_key] < 1.0E-10 && julia_dict[julia_key] < 1.0E-10 ? true :
      isapprox(gamess_dict[gamess_key], julia_dict[julia_key], atol = 1e-9)
    if !same
      println("Not same!: $julia_key => $gamess_key: ", julia_dict[julia_key], ", ", gamess_dict[gamess_key])
    end
  end
end


text:
# Lift a value and a diff instance to a wrapped Diffed instance.
Δ(v, d) = Diffed(v, d)
Δ(v, d::Type) = Diffed(v, d())

# A bunch of nice passes which clean up the IR after tracing.
# Other cleaning passes can be found in transforms.jl.
# This function _does not_ remove dead code (because this interferes with inference semantics).
partial_cleanup!(ir) = ir |> inline_consts! |> partials! |> ssa! |> prune! |> IRTools.renumber

# Convenience - run trace with DiffInterpreter primitives.
function _propagate(f, Dfs, args)
    ir = lower_to_ir(f, args...)
    typed_args = map(zip(IRTools.arguments(ir)[2 : end], Dfs)) do (a, t)
        a => t
    end
    env = Dict{Any, Any}(typed_args)
    env[IRTools.var(1)] = f
    tr = infer!(DiffInterpreter(), env, ir)
    tr
end

# Lift a runtime diff type to the simple indicator Change/NoChange types.
function create_flip_diff(a::Type{Diffed{K, DV}}) where {K, DV}
    DV != NoChange && return Change
    NoChange
end


text:
function AnywhereLevelTwo()
    x = trace(:x, Normal(0.0, 1.0))
    trace(:loop, AnywhereInLoop)
    return x
end

function AnywhereLevelOne()
    x = trace(:x, Normal(0.0, 1.0))
    y = trace(:y, AnywhereLevelTwo)
    return x
end

function AnywhereTopLevel()
    x = trace(:x, Normal(0.0, 1.0))
    y = trace(:y, AnywhereLevelOne)
    return x
end

function AnywhereInLoop()
    x = trace(:q, Normal(0.0, 1.0))
    for i in 1:50
        z = trace(:q => i, Normal(0.0, 1.0))
    end
end

@testset "Constrained targets" begin

    @testset "Constrained by address" begin
        sel = target([(:x, ) => 5.0,
                      (:y, :q, :z) => 10.0])
        @test sel == sel
        @test haskey(sel, (:x, ))
        @test haskey(sel, :x)
        @test haskey(sel, (:y, :q, :z))
        @test sel[:y, :q, :z] == 10.0

        sel2 = target([(:z, :m) => 15.0])
        merge!(sel, sel2)
        @test sel == target([(:x, ) => 5.0,
                             (:m, :q, :z) => 10.0,
                             (:z, :m) => 15.0])
        arr = array(sel, Float64)
        @test arr == [5.0, 10.0, 15.0]
        back = target(sel, arr)
        @test back == sel
    end

    @testset "Unconstrained by address" begin
        sel = target([(:x, ), (:x, :q, :z)])
        @test sel == sel
        @test haskey(sel, (:x, ))
        @test haskey(sel, :x)
        @test haskey(sel, (:x, :q, :z))

        sel2 = target([(:z, :m)])
        merge!(sel, sel2)
        @test sel == target([(:x, ), (:x, :q, :z), (:z, :m)])
    end

    #@testset "Filtering" begin
    #    observations = target([(:x, ) => 5.0,
    #                           (:z, :x) => 5.0,
    #                           (:z, :z, :y) => 5.0])
    #    filtered = filter(x -> x == :y, x -> true, observations)
    #    @test has_top(filtered, (:z, :z, :y))
    #    filtered = filter(x -> x == :x, x -> true, observations)
    #    @test !has_top(filtered, (:z, :z, :y))
    #    @test has_top(filtered, :x)
    #    @test has_top(filtered, (:z, :x))
    #end
end


text:
"""
    Multiply(input=inputDimension, output=outputDimension, winit=xavier, atype=KnetLayers.arrtype)

Creates a matrix multiplication layer based on `inputDimension` and `outputDimension`.
    (m::Multiply) = m.w * x

By default parameters initialized with xavier, you may change it with `winit` argument.

# Keywords
* `input=inputDimension`: input dimension
* `output=outputDimension`: output dimension
* `winit=xavier`: weight initialization distribution
* `atype=KnetLayers.arrtype` : array type for parameters.
   Default value is KnetArray{Float32} if you have gpu device. Otherwise it is Array{Float32}
"""
mutable struct Multiply <: Layer
    weight
end

Multiply(;input::Int, output::Int, winit=xavier, atype=arrtype) =  Multiply(param(output, input; init=winit, atype=atype))

(m::Multiply)(x::Array{<:Integer}) = m.weight[:,x] # Lookup (EmbedLayer)

# TODO: Find a faster (or compound) way for tensor-product
function (m::Multiply)(x; keepsize=true)
    if ndims(x) > 2
        s = size(x)
        y = m.weight * reshape(x, s[1], prod(s[2:end]))
        return (keepsize ? reshape(y, size(y, 1), s[2:end]...) : y)
    else
        return m.weight * x
    end
end

"""
    Embed(input=inputSize, output=embedSize, winit=xavier, atype=KnetLayers.arrtype)
Creates an embedding layer according to given `inputSize` and `embedSize` where `inputSize` is your number of unique items you want to embed, and `embedSize` is the size of output vectors.
By default parameters initialized with xavier, you yam change it with `winit` argument.

    (m::Embed)(x::Array{T}) where T<:Integer
    (m::Embed)(x; keepsize=true)

Embed objects are callable with an input which is either and integer array
(one hot encoding) or an N-dimensional matrix. For N-dimensional matrix,
`size(x,1)==inputSize`

# Keywords

* `input=inputDimension`: input dimension
* `output=embeddingDimension`: output dimension
* `winit=xavier`: weight initialization distribution
* `atype=KnetLayers.arrtype` : array type for parameters.
   Default value is KnetArray{Float32} if you have gpu device. Otherwise it is Array{Float32}

"""
Embed = Multiply

# TODO: find a better documentation style e.g put input : and output etc.
# Input: Type of its input is an ::Array{T} where T<: Integer

"""
Gemm
"""

mutable struct Gemm
    alpha
    beta
    transA
    transB
end

# inputs: matrices to be multiplied, handle C somehow (maybe zeros() with the same size as A*B)
function (g::Gemm)(A, B, C)
    if g.transA != 0; A = transpose(A); end
    if g.transB != 0; B = transpose(B); end
    g.alpha * A * B + (g.beta * C)
end

"""
    Linear(input=inputSize, output=outputSize, winit=xavier, binit=zeros, atype=KnetLayers.arrtype)

Creates and linear layer according to given `inputSize` and `outputSize`.

# Keywords
* `input=inputSize`   input dimension
* `output=outputSize` output dimension
* `winit=xavier`: weight initialization distribution
* `bias=zeros`: bias initialization distribution
* `atype=KnetLayers.arrtype` : array type for parameters.
   Default value is KnetArray{Float32} if you have gpu device. Otherwise it is Array{Float32}
"""
mutable struct Linear <: Layer
    mult::Multiply
    bias
end

function Linear(;input::Int, output::Int, winit=xavier, binit=zeros, atype=arrtype)
    Linear(Multiply(input=input, output=output, winit=winit, atype=atype),param(output, init=binit, atype=atype))
end

(m::Linear)(x) = m.mult(x) .+ m.bias

"""
    Dense(input=inputSize, output=outputSize, activation=ReLU(), winit=xavier, binit=zeros, atype=KnetLayers.arrtype)

Creates and deense layer according to given `input=inputSize` and `output=outputSize`.
If activation is `nothing`, it acts like a `Linear` Layer.
# Keywords
* `input=inputSize`   input dimension
* `output=outputSize` output dimension
* `winit=xaiver`: weight initialization distribution
* `bias=zeros`:   bias initialization distribution
* `activation=ReLU()`  activation function(it does not broadcast) or an  activation layer
* `atype=KnetLayers.arrtype` : array type for parameters.
   Default value is KnetArray{Float32} if you have gpu device. Otherwise it is Array{Float32}
"""
mutable struct Dense <: Layer
    linear::Linear
    activation
end

function Dense(;input::Int, output::Int, activation=ReLU(), winit=xavier, binit=zeros, atype=arrtype)
    Dense(Linear(input=input, output=output, winit=winit, binit=binit, atype=atype), activation)
end

(m::Dense)(x) = m.activation===nothing ? m.linear(x) : m.activation(m.linear(x))

#TO-DO: Remove after the issue is resolved:
#https://github.com/denizyuret/Knet.jl/issues/418
"""
    BatchNorm(channels:Int;options...)
    (m::BatchNorm)(x;training=false) #forward run
# Options
* `momentum=0.1`: A real number between 0 and 1 to be used as the scale of
 last mean and variance. The existing running mean or variance is multiplied by
 (1-momentum).
* `mean=nothing': The running mean.
* `var=nothing`: The running variance.
* `meaninit=zeros`: The function used for initialize the running mean. Should either be `nothing` or
of the form `(eltype, dims...)->data`. `zeros` is a good option.
* `varinit=ones`: The function used for initialize the run
* `dataType=eltype(KnetLayers.arrtype)` : element type ∈ {Float32,Float64} for parameters. Default value is `eltype(KnetLayers.arrtype)`
* `usegpu=KnetLayers.arrtype <: KnetArray` :
# Keywords
* `training`=nothing: When training is true, the mean and variance of x are used and moments
 argument is modified if it is provided. When training is false, mean and variance
 stored in the moments argument are used. Default value is true when at least one
 of x and params is AutoGrad.Value, false otherwise.
"""
mutable struct BatchNorm <: Layer
    params
    moments::Knet.BNMoments
end

function BatchNorm(channels::Int; usegpu = arrtype <: KnetArray, dataType=eltype(arrtype), o...)
    w = bnparams(dataType,channels)
    m = bnmoments(;o...)
    BatchNorm(usegpu ? Param(KnetArray(w)) : Param(w),m)
end
(m::BatchNorm)(x;o...) = batchnorm(x,m.moments,m.params;o...)

#3D -> 2D
struct Flatten <: Layer
end

function (f::Flatten)(x)
    s = size(x)
    # batch mode (..., 1) at the end
    if s[end] == 1
        #reshape(x, :, size(x)[end-1], 1)
        reshape(x, :, size(x)[end])

    # no batch
    else
        reshape(x, :, size(x)[end])
    end
end

"""
Adds any number of inputs together.
Doesn't check size or type.
"""
struct Add; end
(a::Add)(args...) = +(args...)

"""
Mul
Performs element-wise binary multiplication (with Numpy-style broadcasting support).
"""
struct Mul; end
(m::Mul)(A, B) = A .* B

"""
Div
Performs element-wise binary division (with Numpy-style broadcasting support).
"""
struct Div; end
(d::Div)(A, B) = A ./ B

"""
Sub
Performs element-wise binary subtraction (with Numpy-style broadcasting support).
"""
struct Sub; end
(s::Sub)(A, B) = A .- B

struct Exp; end
(e::Exp)(A) = A.^2

struct Transpose
    perm
end
(t::Transpose)(A) = permutedims(A, t.perm)

struct Reshape
    shape
end
(r::Reshape)(A) = reshape(A, r.shape)

struct Cast
    to
end
(c::Cast)(A) = convert(c.to, A)

# Not working now!
struct Slice
    starts
    ends
    axes
    steps
end

function (s::Slice)(A)
    A[s.axes][s.starts:s.steps:s.ends]
end

struct MatMul
end
(m::MatMul)(A,B) = A * B


text:
module Grid

"""
Data type to represent 2D grid.
# Fields
* `x` : horizontal values
* `z` : vertical values
* `nx` : number of samples in horizontal direction
* `nz` : number of samples in vertical direction
* `δx` : sampling interval in horizontal direction
* `δz` : sampling interval in vertical direction
"""
mutable struct M2D
	x::Array{Float64,1}
	z::Array{Float64,1}
	nx::Int64
	nz::Int64
	δx::Float64
	δz::Float64
	"adding conditions that are to be false while construction"
	M2D(x, z, nx, nz, δx, δz) =
		any([
       		  δx < 0.0, length(x) != nx,
       		  δz < 0.0, length(z) != nz
		  ]) ?
		error("error in M2D construction") : new(x, z, nx, nz, δx, δz)
end

"Logical operation for `M2D`"
function Base.isequal(grid1::M2D, grid2::M2D)
	return all([(isequal(getfield(grid1, name),getfield(grid2, name))) for name in fieldnames(M2D)])
end

"""
Construct 2-D grid based on number of samples.

# Arguments

* `xmin::Float64` : first value of second dimension
* `xmax::Float64` : last value of second dimension
* `zmin::Float64` : first value of first dimension
* `zmax::Float64` : firs value of first dimension
* `nx::Int64` : size of second dimension
* `nz::Int64` : size of first dimension

# Return

* a `M2D` grid
"""
function M2D(xmin::Float64, xmax::Float64,
	zmin::Float64, zmax::Float64,
	nx::Int64, nz::Int64,
	)
	x = Array(range(xmin,stop=xmax,length=nx));
	z = Array(range(zmin,stop=zmax,length=nz));
	return M2D(x, z, nx, nz, x[2]-x[1], z[2]-z[1])
end

"""
Construct 2-D grid based on sampling intervals.

# Arguments

* `xmin::Float64` : first value of second dimension
* `xmax::Float64` : last value of second dimension
* `zmin::Float64` : first value of first dimension
* `zmax::Float64` : firs value of first dimension
* `δx::Float64` : second sampling interval
* `δz::Float64` : first sampling interval

# Return

* a `M2D` grid
"""
function M2D(xmin::Float64, xmax::Float64,
	zmin::Float64, zmax::Float64,
	δx::Float64, δz::Float64,
	)
	x = [xx for xx in xmin:δx:xmax]
	z = [zz for zz in zmin:δz:zmax]
	nx, nz = size(x,1), size(z,1);
	return M2D(x, z, nx, nz,  x[2]-x[1], z[2]-z[1])
end

"""
Resample a 2-D grid.

# Arguments

* `grid::M2D` : input grid this is to be resampled
* `δx::Float64` : new second sampling interval
* `δz::Float64` : new first sampling interval

# Return

* a `M2D` resampled grid
"""
M2D_resamp(grid::M2D, δx::Float64, δz::Float64) = M2D(grid.x[1], grid.x[end],
					grid.z[1], grid.z[end], δx, δz)

"""
Extend M2D by on its PML grid points on all sides.
"""
function M2D_pml_pad_trun(mgrid::M2D; flag::Int64=1, npml=0)

	if(isequal(flag,1))
		xmin = mgrid.x[1] - npml*mgrid.δx
		xmax = mgrid.x[end] + npml*mgrid.δx
		zmin = mgrid.z[1] - npml*mgrid.δz
		zmax = mgrid.z[end] + npml*mgrid.δz
		return M2D(xmin,xmax,zmin,zmax, mgrid.nx+2*npml,mgrid.nz+2*npml)
	elseif(isequal(flag,-1))
		xmin = mgrid.x[1] + npml*mgrid.δx
		xmax = mgrid.x[end] - npml*mgrid.δx
		zmin = mgrid.z[1] + npml*mgrid.δz
		zmax = mgrid.z[end] - npml*mgrid.δz
		return M2D(xmin,xmax,zmin,zmax, mgrid.nx-2*npml,mgrid.nz-2*npml)
	else
		error("invalid flag")
	end

end

"""
Return the X and Z positions of the boundary of mgrid
* `attrib::Symbol` :
  * `=:inner`
  * `=:outer`
* `onlycount::Bool=false`
"""
function M2D_boundary(mgrid::M2D, nlayer::Int64, attrib::Symbol; onlycount::Bool=false)
	if(attrib == :inner)
		x = mgrid.x; z = mgrid.z;
	elseif(attrib == :outer)
		# extending x and z and then choosing inner layers
		x = vcat(
	   		[mgrid.x[1]-(ilayer)*mgrid.δx for ilayer=1:nlayer:-1],
	   		mgrid.x,
			[mgrid.x[end]+(ilayer)*mgrid.δx for ilayer=1:nlayer],
			)
		z = vcat(
	   		[mgrid.z[1]-(ilayer)*mgrid.δz for ilayer=1:nlayer:-1],
	   		mgrid.z,
			[mgrid.z[end]+(ilayer)*mgrid.δz for ilayer=1:nlayer],
			)
	end
	nx = length(x); nz = length(z);
	bx = vcat(
		  x, x[end]*ones(nz-1),	  x[end-1:-1:1], x[1]*ones(nz-2))
	bz = vcat(
		  z[1]*ones(nx), z[2:end],  z[end]*ones(nx-1), z[end-1:-1:2])

	if(nlayer > 1)
		bx = vcat(bx,  x[2:end-1], x[end-1]*ones(nz-3),
		  x[end-2:-1:2], x[2]*ones(nz-4))

		bz = vcat(bz,  z[2]*ones(nx-2), z[3:end-1],
		  z[end-1]*ones(nx-3), z[end-2:-1:3])
	end
	if(nlayer > 2)
		bx = vcat(bx,  x[3:end-2], x[end-2]*ones(nz-5),
		  x[end-3:-1:3], x[3]*ones(nz-6)
		  )
		bz = vcat(bz,  z[3]*ones(nx-4), z[4:end-2],
		  z[end-2]*ones(nx-5), z[end-3:-1:4]
		  )
	end
	isequal(length(bz), length(bx)) ? nothing : error("unequal dimensions")
	if(onlycount)
		return length(bz)
	else
		return bz, bx, length(bz)
	end
end

"""
Data type to represent 1D grid.
# Fields
* `x` : values
* `nx` : number of samples
* `δx` : sampling interval
"""
mutable struct M1D
	x::Array{Float64}
	nx::Int64
	δx::Float64
	"adding conditions that are to be false while construction"
	M1D(x, nx, δx) =
		any([δx < 0.0, length(x) != nx]) ?
			error("error in M1D construction") : new(x, nx, δx)
end

"Logical operation for `M1D`"
function Base.isequal(grid1::M1D, grid2::M1D)
	return all([(isequal(getfield(grid1, name),getfield(grid2, name))) for name in fieldnames(M1D)])
end

"""
Construct 1-D grid based on number of samples.

# Arguments

* `xbeg::Float64` : first value
* `xend::Float64` : last value
* `nx::Int64` : number of samples
"""
function M1D(xbeg::Float64, xend::Float64, nx::Int64)
	x = Array(range(xbeg, stop=xend, length=nx))
	δx = length(x)==1 ? 0. : x[2]-x[1]
	return M1D(x, nx, δx)
end

"""
Construct 1-D grid based on sampling interval.

# Arguments

* `xbeg::Float64` : first value
* `xend::Float64` : last value
* `δx::Float64` : number of samples
"""
function M1D(xbeg::Float64, xend::Float64, δx::Float64)
	x = [tt for tt in xbeg:δx:xend]
	δx = length(x)==1 ? 0. : x[2]-x[1]
	return M1D(x, size(x,1), δx)
end

"""
Grid with both positive and negative samples for a given lag.
Construction makes sure that the number of output grid samples is odd.

* xlag if scalar use same for both +ve and -ve lags, otherwise give vector

# Output
* Grid with both +ve and-ve lags
* number of +ve and -ve lags
"""
function M1D_lag(xlag, δx::Float64)
	if(length(xlag)==1)
		xlag=[xlag[1],xlag[1]]
	end
	nplag=round(Int,xlag[1]/δx)
	nnlag=round(Int,xlag[2]/δx)
	x=[0.0]
	x1=[(it-1)*δx for it in 2:nplag]
	if(x1≠[])
		x=vcat(x,x1)
	end
	nplag=length(x1)
	x2=[(it-1)*δx for it in 2:nnlag]
	if(x2≠[])
		x=vcat(-1.0*reverse(x2,dims=1),x)
	end
	nnlag=length(x2)
	if(x==[])
		δx=0.
		x=zeros(1)
	end
	if(xlag[1]==xlag[2]≠0)
		isodd(length(x)) ? nothing : error("error in creating lag grid")
	end
	return M1D(x, size(x,1), δx), [nplag, nnlag] # lags are one less because of zero lag
end

"""
1-D grid with a different sampling interval
* Not yet implemented for fft grids
"""
function M1D_resamp(grid::M1D, δx::Float64)
	return M1D(grid.x[1], grid.x[end], δx)
end

"""
1-D grid which is has a different size
"""
function M1D_truncate(grid::M1D, xbeg::Float64, xend::Float64)
	ix1 = indmin((grid.x - xbeg).^2.0);
	ix2 = indmin((grid.x - xend).^2.0);
	ixmin = minimum([ix1, ix2]);
	ixmax = maximum([ix1, ix2]);
	x = grid.x[ixmin:ixmax];
	return M1D(x, size(x,1), x[2]-x[1])
end

"grid after FFT"
M1D_fft(grid::M1D) = M1D_fft(grid.nx, inv(grid.nx*grid.δx))
M1D_rfft(grid::M1D) = M1D_rfft(grid.nx, inv(grid.nx*grid.δx))

"""
Frequency grid after rfft
"""
function M1D_rfft(nx::Int64, δ::Float64)
	vec = zeros(div(nx,2)+1);
	# zero lag
	vec[1] = 0.0;

	# +ve
	for i = 1: div(nx,2)
		vec[1+i] = δ * i
	end
	return M1D(vec, length(vec), δ)
end

function M1D_fft(nx::Int64, δ::Float64)
	vec = zeros(nx);
	# zero lag
	vec[1] = 0.0;

	if(isodd(nx))
		# +ve
		for i = 1: div(nx,2)
			vec[1+i] = δ * i
		end
		# -ve
		for i = 1: div(nx,2)
			vec[nx-i+1] = -δ * i
		end
	else
		# +ve
		for i = 1: div(nx,2)
			vec[1+i] = δ * i
		end
		# -ve one less the number of +ve lags
		for i = 1: div(nx,2)-1
			vec[nx-i+1] = -δ * i
		end
	end
	return M1D(vec, nx, δ)
end

"""
Return grid after cross-correlation.
The positive and negative lags are approximately given by lags.
"""
function M1D_xcorr(tgrid; lags=[1.,1.].*abs(tgrid.x[end]-tgrid.x[1]))

	plags=round(Int,  lags[1]*inv(tgrid.δx))
	nlags=round(Int,  lags[2]*inv(tgrid.δx))

	vec=vcat(-1.0*collect(nlags:-1:0),collect(1:plags)).*tgrid.δx
	return M1D(vec, length(vec), tgrid.δx)
end

end # module


text:
export ExchangeRule

"""
    ExchangeRule(graph)

Transition rule for Metropolis-Hastings sampling where at every step a random
couple of sites i,j is switched.

Couples of sites are generated from the graph `graph` (or operator), where all
sites that are connected (or coupled by a 2-body term) are considered for switches.
"""
struct ExchangeRule{T<:Vector} <: MCMCRule
    distances::T
end

function ExchangeRule(H::AbsLinearOperator)
    N_sites = nsites(basis(H))

    all_sites = sites(H)
    if !(eltype(all_sites) <: Vector)
        all_sites = [all_sites]
    end

    couplings = Vector{Tuple{Int,Int}}()
    for coupling=all_sites
        length(coupling) == 1 && continue
        if length(coupling) > 2
            @warn "Can't exchange between 3-site couplings. This coupling is ignored"
            continue
        end
        push!(couplings, (coupling[1], coupling[2]))
    end

    return ExchangeRule{typeof(couplings)}(couplings)
end

function propose_step!(σp::Union{AState,ADoubleState}, s::MetropolisSampler{<:ExchangeRule},
                       net::Union{MatrixNet,KetNet}, c, rc)
    dist = c.rule.distances

    couple_i = rand(c.rng, 1:length(dist))
    i, j = dist[couple_i]
    if σp isa AState
        @inbounds tmp = σp[i]
        @inbounds σp[i] = σp[j]
        @inbounds σp[j] = tmp
    else
        throw("not implemented")
    end
end

function propose_step!(σp_b::Union{AStateBatch,ADoubleStateBatch}, s::MetropolisSampler{<:ExchangeRule},
                       net::Union{MatrixNet,KetNet}, c, rc)
    dist = s.rule.distances

    for σp=states(σp_b)

        couple_i = rand(c.rng, 1:length(dist))
        i, j = dist[couple_i]
        if σp isa AState
            @inbounds tmp = σp[i]
            @inbounds σp[i] = σp[j]
            @inbounds σp[j] = tmp
        else
            throw("not implemented")
        end
    end
end


text:
# ## Turbulence Closures
# In `turbulence.jl` we specify turbulence closures. Currently,
# pointwise models of the eddy viscosity/eddy diffusivity type are
# supported for turbulent shear and tracer diffusivity. Methods currently supported
# are:\
# [`ConstantViscosityWithDivergence`](@ref constant-viscosity)\
# [`SmagorinskyLilly`](@ref smagorinsky-lilly)\
# [`Vreman`](@ref vreman)\
# [`AnisoMinDiss`](@ref aniso-min-diss)\

#md # !!! note
#md #     Usage: This is a quick-ref guide to using turbulence models as a subcomponent
#md #     of `AtmosModel` \
#md #     $\nu$ is the kinematic viscosity, $C_smag$ is the Smagorinsky Model coefficient,
#md #     `turbulence=ConstantViscosityWithDivergence(ν)`\
#md #     `turbulence=SmagorinskyLilly(C_smag)`\
#md #     `turbulence=Vreman(C_smag)`\
#md #     `turbulence=AnisoMinDiss(C_poincare)`

using DocStringExtensions
using CLIMAParameters.Atmos.SubgridScale: inv_Pr_turb
export ConstantViscosityWithDivergence, SmagorinskyLilly, Vreman, AnisoMinDiss
export turbulence_tensors

# ### Abstract Type
# We define a `TurbulenceClosure` abstract type and
# default functions for the generic turbulence closure
# which will be overloaded with model specific functions.
abstract type TurbulenceClosure end

vars_state_conservative(::TurbulenceClosure, FT) = @vars()
vars_state_auxiliary(::TurbulenceClosure, FT) = @vars()

"""
    atmos_init_aux!
Initialise auxiliary variables for turbulence models.
Overload for specific type of turbulence closure.
"""
function atmos_init_aux!(
    ::TurbulenceClosure,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
) end

"""
    atmos_nodal_update_auxiliary_state!
Update auxiliary variables for turbulence models.
Overload for specific turbulence closure type
"""
function atmos_nodal_update_auxiliary_state!(
    ::TurbulenceClosure,
    ::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end

"""
    compute_gradient_argument!
Assign pre-gradient-transform variables specific to turbulence models.
"""
function compute_gradient_argument!(
    ::TurbulenceClosure,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end
"""
    compute_gradient_flux!(::TurbulenceClosure, _...)
Post-gradient-transformed variables specific to turbulence models.
"""
function compute_gradient_flux!(
    ::TurbulenceClosure,
    ::Orientation,
    diffusive,
    ∇transform,
    state,
    aux,
    t,
) end

function turbulence_tensors end

turbulence_tensors(atmos::AtmosModel, args...) =
    turbulence_tensors(atmos.turbulence, atmos, args...)

"""
    ν, D_t, τ = turbulence_tensors(
                    ::TurbulenceClosure,
                    orientation::Orientation,
                    param_set::AbstractParameterSet,
                    state::Vars,
                    diffusive::Vars,
                    aux::Vars,
                    t::Real
                )
Compute the kinematic viscosity (`ν`), the diffusivity (`D_t`) and SGS momentum flux tensor (`τ`)
for a given turbulence closure. Each closure overloads this method with the appropriate calculations
for the returned quantities.

# Arguments

- `::TurbulenceClosure` = Struct identifier for turbulence closure model
- `orientation` = `AtmosModel.orientation`
- `param_set` = `AtmosModel.param_set`
- `state` = Array of prognostic (state) variables. See `vars_state_conservative` in `AtmosModel`
- `diffusive` = Array of diffusive variables
- `aux` = Array of auxiliary variables
- `t` = time
"""
turbulence_tensors(m::TurbulenceClosure, atmos::AtmosModel, args...) =
    turbulence_tensors(m, atmos.orientation, atmos.param_set, args...)

# We also provide generic math functions for use within the turbulence closures,
# commonly used quantities such as the [principal tensor invariants](@ref tensor-invariants), handling of
# [symmetric tensors](@ref symmetric-tensors) and [tensor norms](@ref tensor-norms)are addressed.

# ### [Pricipal Invariants](@id tensor-invariants)
# ```math
# \textit{I}_{1} = \mathrm{tr(X)} \\
# \textit{I}_{2} = (\mathrm{tr(X)}^2 - \mathrm{tr(X^2)}) / 2 \\
# \textit{I}_{3} = \mathrm{det(X)} \\
# ```
"""
    principal_invariants(X)

Calculates principal invariants of a tensor `X`. Returns 3 element tuple containing the invariants.
"""
function principal_invariants(X)
    first = tr(X)
    second = (first^2 - tr(X .^ 2)) / 2
    third = det(X)
    return (first, second, third)
end

# ### [Symmetrize](@id symmetric-tensors)
# ```math
# \frac{\mathrm{X} + \mathrm{X}^{T}}{2} \\
# ```
"""
    symmetrize(X)

Given a (3,3) second rank tensor X, compute `(X + X')/2`, returning a
symmetric `SHermitianCompact` object.
"""
function symmetrize(X::StaticArray{Tuple{3, 3}})
    SHermitianCompact(SVector(
        X[1, 1],
        (X[2, 1] + X[1, 2]) / 2,
        (X[3, 1] + X[1, 3]) / 2,
        X[2, 2],
        (X[3, 2] + X[2, 3]) / 2,
        X[3, 3],
    ))
end

# ### [2-Norm](@id tensor-norms)
# Given a tensor X, return the tensor dot product
# ```math
# \sum_{i,j} S_{ij}^2
# ```
"""
    norm2(X)
Given a tensor `X`, computes its tensor dot product.
"""
function norm2(X::SMatrix{3, 3, FT}) where {FT}
    abs2(X[1, 1]) +
    abs2(X[2, 1]) +
    abs2(X[3, 1]) +
    abs2(X[1, 2]) +
    abs2(X[2, 2]) +
    abs2(X[3, 2]) +
    abs2(X[1, 3]) +
    abs2(X[2, 3]) +
    abs2(X[3, 3])
end
function norm2(X::SHermitianCompact{3, FT, 6}) where {FT}
    abs2(X[1, 1]) +
    2 * abs2(X[2, 1]) +
    2 * abs2(X[3, 1]) +
    abs2(X[2, 2]) +
    2 * abs2(X[3, 2]) +
    abs2(X[3, 3])
end

# ### [Strain-rate Magnitude](@id strain-rate-magnitude)
# By definition, the strain-rate magnitude, as defined in
# standard turbulence modelling is computed such that
# ```math
# |\mathrm{S}| = \sqrt{2 \sum_{i,j} \mathrm{S}_{ij}^2}
# ```
# where
# ```math
# \vec{S}(\vec{u}) = \frac{1}{2}  \left(\nabla\vec{u} +  \left( \nabla\vec{u} \right)^T \right)
# ```
# \mathrm{S} is the rate-of-strain tensor. (Symmetric component of the velocity gradient). Note that the
# skew symmetric component (rate-of-rotation) is not currently computed.
"""
    strain_rate_magnitude(S)
Given the rate-of-strain tensor `S`, computes its magnitude.
"""
function strain_rate_magnitude(S::SHermitianCompact{3, FT, 6}) where {FT}
    return sqrt(2 * norm2(S))
end

# ### [Constant Viscosity Model](@id constant-viscosity)
# `ConstantViscosityWithDivergence` requires a user to specify the constant viscosity (kinematic)
# and appropriately computes the turbulent stress tensor based on this term. Diffusivity can be
# computed using the turbulent Prandtl number for the appropriate problem regime.
# ```math
# \tau = - 2 \nu \mathrm{S}
# ```
"""
    ConstantViscosityWithDivergence <: TurbulenceClosure

Turbulence with constant dynamic viscosity (`ρν`).
Divergence terms are included in the momentum flux tensor.

# Fields

$(DocStringExtensions.FIELDS)
"""
struct ConstantViscosityWithDivergence{FT} <: TurbulenceClosure
    "Dynamic Viscosity [kg/m/s]"
    ρν::FT
end

vars_state_gradient(::ConstantViscosityWithDivergence, FT) = @vars()
vars_state_gradient_flux(::ConstantViscosityWithDivergence, FT) =
    @vars(S::SHermitianCompact{3, FT, 6})

function compute_gradient_flux!(
    ::ConstantViscosityWithDivergence,
    ::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

    diffusive.turbulence.S = symmetrize(∇transform.u)
end

function turbulence_tensors(
    m::ConstantViscosityWithDivergence,
    orientation::Orientation,
    param_set::AbstractParameterSet,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)

    FT = eltype(state)
    _inv_Pr_turb::FT = inv_Pr_turb(param_set)
    S = diffusive.turbulence.S
    ν = m.ρν / state.ρ
    D_t = ν * _inv_Pr_turb
    τ = (-2 * ν) * S + (2 * ν / 3) * tr(S) * I
    return ν, D_t, τ
end

# ### [Smagorinsky-Lilly](@id smagorinsky-lilly)
# The Smagorinsky turbulence model, with Lilly's correction to
# stratified atmospheric flows, is included in CLIMA.
# The input parameter to this model is the Smagorinsky coefficient.
# For atmospheric flows, the coefficient `C_smag` typically takes values between
# 0.15 and 0.23. Flow dependent `C_smag` are currently not supported (e.g. Germano's
# extension). The Smagorinsky-Lilly model does not contain explicit filtered terms.
# #### Equations
# ```math
# \nu = (C_{s} \mathrm{f}_{b} \Delta)^2 \sqrt{|\mathrm{S}|}
# ```
# with the stratification correction term
# ```math
# f_{b} =
#    \begin{cases}
#    1 & \mathrm{Ri} \leq 0 ,\\
#    \max(0, 1 - \mathrm{Ri} / \mathrm{Pr}_{t})^{1/4} & \mathrm{Ri} > 0 .
#    \end{cases}
# ```
# ```math
# \mathrm{Ri} =  \frac{N^2}{{|S|}^2}
# ```
# ```math
# N = \left( \frac{g}{\theta_v} \frac{\partial \theta_v}{\partial z}\right)^{1/2}
# ```
# Here, $\mathrm{Ri}$ and $\mathrm{Pr}_{t}$ are the Richardson and
# turbulent Prandtl numbers respectively.  $\Delta$ is the mixing length in the
# relevant coordinate direction. We use the DG metric terms to determine the
# local effective resolution (see `src/Mesh/Geometry.jl`), and modify the vertical lengthscale by the
# stratification correction factor $\mathrm{f}_{b}$ so that $\Delta_{vert} = \Delta z f_b$.
"""
    SmagorinskyLilly <: TurbulenceClosure

# Fields

$(DocStringExtensions.FIELDS)

# Smagorinsky Model Reference
    article{doi:10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2,
      author = {Smagorinksy, J.},
      title = {General circulation experiments with the primitive equations},
      journal = {Monthly Weather Review},
      volume = {91},
      number = {3},
      pages = {99-164},
      year = {1963},
      doi = {10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2},
      URL = {https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2},
      eprint = {https://doi.org/10.1175/1520-0493(1963)091<0099:GCEWTP>2.3.CO;2}
      }

# Lilly Model Reference
    article{doi:10.1111/j.2153-3490.1962.tb00128.x,
      author = {LILLY, D. K.},
      title = {On the numerical simulation of buoyant convection},
      journal = {Tellus},
      volume = {14},
      number = {2},
      pages = {148-172},
      doi = {10.1111/j.2153-3490.1962.tb00128.x},
      url = {https://onlinelibrary.wiley.com/doi/abs/10.1111/j.2153-3490.1962.tb00128.x},
      eprint = {https://onlinelibrary.wiley.com/doi/pdf/10.1111/j.2153-3490.1962.tb00128.x},
      year = {1962}
      }

# Brunt-Väisälä Frequency Reference
    Brunt-Vaisala frequency N² defined as in equation (1b) in
      Durran, D.R. and J.B. Klemp, 1982:
      On the Effects of Moisture on the Brunt-Väisälä Frequency.
      J. Atmos. Sci., 39, 2152–2158,
      https://doi.org/10.1175/1520-0469(1982)039<2152:OTEOMO>2.0.CO;2

"""
struct SmagorinskyLilly{FT} <: TurbulenceClosure
    "Smagorinsky Coefficient [dimensionless]"
    C_smag::FT
end

vars_state_auxiliary(::SmagorinskyLilly, FT) = @vars(Δ::FT)
vars_state_gradient(::SmagorinskyLilly, FT) = @vars(θ_v::FT)
vars_state_gradient_flux(::SmagorinskyLilly, FT) =
    @vars(S::SHermitianCompact{3, FT, 6}, N²::FT)

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

function compute_gradient_argument!(
    m::SmagorinskyLilly,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.turbulence.θ_v = aux.moisture.θ_v
end

function compute_gradient_flux!(
    ::SmagorinskyLilly,
    orientation::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)

    diffusive.turbulence.S = symmetrize(∇transform.u)
    ∇Φ = ∇gravitational_potential(orientation, aux)
    diffusive.turbulence.N² =
        dot(∇transform.turbulence.θ_v, ∇Φ) / aux.moisture.θ_v
end

function turbulence_tensors(
    m::SmagorinskyLilly,
    orientation::Orientation,
    param_set::AbstractParameterSet,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)

    FT = eltype(state)
    _inv_Pr_turb::FT = inv_Pr_turb(param_set)
    S = diffusive.turbulence.S
    normS = strain_rate_magnitude(S)
    k̂ = vertical_unit_vector(orientation, param_set, aux)

    # squared buoyancy correction
    Richardson = diffusive.turbulence.N² / (normS^2 + eps(normS))
    f_b² = sqrt(clamp(FT(1) - Richardson * _inv_Pr_turb, FT(0), FT(1)))
    ν₀ = normS * (m.C_smag * aux.turbulence.Δ)^2 + FT(1e-5)
    ν = SVector{3, FT}(ν₀, ν₀, ν₀)
    ν_v = k̂ .* dot(ν, k̂)
    ν_h = ν₀ .- ν_v
    ν = SDiagonal(ν_h + ν_v .* f_b²)
    D_t = diag(ν) * _inv_Pr_turb
    τ = -2 * ν * S
    return ν, D_t, τ
end

# ### [Vreman Model](@id vreman)
# Vreman's turbulence model for anisotropic flows, which provides a
# less dissipative solution (specifically in the near-wall and transitional regions)
# than the Smagorinsky-Lilly method. This model
# relies of first derivatives of the velocity vector (i.e., the gradient tensor).
# By design, the Vreman model handles transitional as well as fully turbulent flows adequately.
# The input parameter to this model is the Smagorinsky coefficient - the coefficient is modified
# within the model functions to account for differences in model construction.
# #### Equations
# ```math
# \nu_{t} = 2.5 C_{s}^2 \sqrt{\frac{B_{\beta}}{u_{i,j}u_{i,j}}},
# ```
# where ($i,j, m = (1,2,3)$)
# ```math
# \begin{align}
# B_{\beta} &= \beta_{11}\beta_{22} + \beta_{11}\beta_{33} + \beta_{22}\beta_{33} - (\beta_{13}^2 + \beta_{12}^2 + \beta_{23}^2) \\
# \beta_{ij} &= \Delta_{m}^2 u_{i, m} u_{j, m} \\
# u_{i,j} &= \frac{\partial u_{i}}{\partial x_{j}}.
# \end{align}
# ```
"""
    Vreman{FT} <: TurbulenceClosure

Filter width Δ is the local grid resolution calculated from the mesh metric tensor. A Smagorinsky coefficient
is specified and used to compute the equivalent Vreman coefficient.

1) ν_e = √(Bᵦ/(αᵢⱼαᵢⱼ)) where αᵢⱼ = ∂uⱼ∂uᵢ with uᵢ the resolved scale velocity component.
2) βij = Δ²αₘᵢαₘⱼ
3) Bᵦ = β₁₁β₂₂ + β₂₂β₃₃ + β₁₁β₃₃ - β₁₂² - β₁₃² - β₂₃²
βᵢⱼ is symmetric, positive-definite.
If Δᵢ = Δ, then β = Δ²αᵀα

# Fields

$(DocStringExtensions.FIELDS)

# Reference

    @article{Vreman2004,
      title={An eddy-viscosity subgrid-scale model for turbulent shear flow: Algebraic theory and applications},
      author={Vreman, AW},
      journal={Physics of fluids},
      volume={16},
      number={10},
      pages={3670--3681},
      year={2004},
      publisher={AIP}
    }
"""
struct Vreman{FT} <: TurbulenceClosure
    "Smagorinsky Coefficient [dimensionless]"
    C_smag::FT
end
vars_state_auxiliary(::Vreman, FT) = @vars(Δ::FT)
vars_state_gradient(::Vreman, FT) = @vars(θ_v::FT)
vars_state_gradient_flux(::Vreman, FT) = @vars(∇u::SMatrix{3, 3, FT, 9}, N²::FT)

function atmos_init_aux!(::Vreman, ::AtmosModel, aux::Vars, geom::LocalGeometry)
    aux.turbulence.Δ = lengthscale(geom)
end
function compute_gradient_argument!(
    m::Vreman,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.turbulence.θ_v = aux.moisture.θ_v
end
function compute_gradient_flux!(
    ::Vreman,
    orientation::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    diffusive.turbulence.∇u = ∇transform.u
    ∇Φ = ∇gravitational_potential(orientation, aux)
    diffusive.turbulence.N² =
        dot(∇transform.turbulence.θ_v, ∇Φ) / aux.moisture.θ_v
end

function turbulence_tensors(
    m::Vreman,
    orientation::Orientation,
    param_set::AbstractParameterSet,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    _inv_Pr_turb::FT = inv_Pr_turb(param_set)
    α = diffusive.turbulence.∇u
    S = symmetrize(α)
    k̂ = vertical_unit_vector(orientation, param_set, aux)

    normS = strain_rate_magnitude(S)
    Richardson = diffusive.turbulence.N² / (normS^2 + eps(normS))
    f_b² = sqrt(clamp(1 - Richardson * _inv_Pr_turb, 0, 1))

    β = f_b² * (aux.turbulence.Δ)^2 * (α' * α)
    Bβ = principal_invariants(β)[2]

    ν₀ = m.C_smag^2 * FT(2.5) * sqrt(abs(Bβ / (norm2(α) + eps(FT))))

    ν = SVector{3, FT}(ν₀, ν₀, ν₀)
    ν_v = k̂ .* dot(ν, k̂)
    ν_h = ν₀ .- ν_v
    ν = SDiagonal(ν_h + ν_v .* f_b²)
    D_t = diag(ν) * _inv_Pr_turb
    τ = -2 * ν * S
    return ν, D_t, τ
end

# ### [Anisotropic Minimum Dissipation](@id aniso-min-diss)
# This method is based Vreugdenhil and Taylor's minimum-dissipation eddy-viscosity model.
# The principles of the Rayleigh quotient minimizer are applied to the energy dissipation terms in the
# conservation equations, resulting in a maximum dissipation bound, and a model for
# eddy viscosity and eddy diffusivity.
# ```math
# \nu_e = (\mathrm{C}\delta)^2  \mathrm{max}\left[0, - \frac{\hat{\partial}_k \hat{u}_{i} \hat{\partial}_k \hat{u}_{j} \mathrm{\hat{S}}_{ij}}{\hat{\partial}_p \hat{u}_{q}} \right]
# ```
"""
    AnisoMinDiss{FT} <: TurbulenceClosure

Filter width Δ is the local grid resolution calculated from the mesh metric tensor. A Poincare coefficient
is specified and used to compute the equivalent AnisoMinDiss coefficient (computed as the solution to the
eigenvalue problem for the Laplacian operator).

# Fields
$(DocStringExtensions.FIELDS)

# Reference

    @article{
        doi:10.1063/1.5037039,
        author = {Vreugdenhil,Catherine A.  and Taylor,John R. },
        title = {Large-eddy simulations of stratified plane Couette flow using the anisotropic minimum-dissipation model},
        journal = {Physics of Fluids},
        volume = {30},
        number = {8},
        pages = {085104},
        year = {2018},
        doi = {10.1063/1.5037039},
        URL = {https://doi.org/10.1063/1.5037039}
    }

"""
struct AnisoMinDiss{FT} <: TurbulenceClosure
    C_poincare::FT
end
vars_state_auxiliary(::AnisoMinDiss, FT) = @vars(Δ::FT)
vars_state_gradient(::AnisoMinDiss, FT) = @vars(θ_v::FT)
vars_state_gradient_flux(::AnisoMinDiss, FT) =
    @vars(∇u::SMatrix{3, 3, FT, 9}, N²::FT)
function atmos_init_aux!(
    ::AnisoMinDiss,
    ::AtmosModel,
    aux::Vars,
    geom::LocalGeometry,
)
    aux.turbulence.Δ = lengthscale(geom)
end
function compute_gradient_argument!(
    m::AnisoMinDiss,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    transform.turbulence.θ_v = aux.moisture.θ_v
end
function compute_gradient_flux!(
    ::AnisoMinDiss,
    orientation::Orientation,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ∇Φ = ∇gravitational_potential(orientation, aux)
    diffusive.turbulence.∇u = ∇transform.u
    diffusive.turbulence.N² =
        dot(∇transform.turbulence.θ_v, ∇Φ) / aux.moisture.θ_v
end
function turbulence_tensors(
    m::AnisoMinDiss,
    orientation::Orientation,
    param_set::AbstractParameterSet,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
)
    FT = eltype(state)
    k̂ = vertical_unit_vector(orientation, param_set, aux)
    _inv_Pr_turb::FT = inv_Pr_turb(param_set)

    ∇u = diffusive.turbulence.∇u
    S = symmetrize(∇u)
    normS = strain_rate_magnitude(S)

    δ = aux.turbulence.Δ
    Richardson = diffusive.turbulence.N² / (normS^2 + eps(normS))
    f_b² = sqrt(clamp(1 - Richardson * _inv_Pr_turb, 0, 1))

    δ_vec = SVector(δ, δ, δ)
    δ_m = δ_vec ./ transpose(δ_vec)
    ∇û = ∇u .* δ_m
    Ŝ = symmetrize(∇û)
    ν₀ =
        (m.C_poincare .* δ_vec) .^ 2 * max(
            FT(1e-5),
            -dot(transpose(∇û) * (∇û), Ŝ) / (dot(∇û, ∇û) .+ eps(normS)),
        )

    ν_v = k̂ .* dot(ν₀, k̂)
    ν_h = ν₀ .- ν_v
    ν = SDiagonal(ν_h + ν_v .* f_b²)
    D_t = diag(ν) * _inv_Pr_turb
    τ = -2 * ν * S
    return ν, D_t, τ
end


text:
using Test

@testset "./src/file_io.jl" begin

using BellScenario

@testset "pretty_print_txt()" begin
    test_dir = "./test/integration/files/"
    filename = test_dir * "pretty_print_txt_test"
    filename_txt = filename * ".txt"

    if isfile(filename_txt)
        rm(filename_txt)
    end

    bell_games = BellGame.([[1 0 0;1 0 0;0 0 0],[1 0 0;0 1 0;0 0 1]], [1,2])
    match_string = "(1)  1 ≥ \n\t[1, 0, 0]\n\t[1, 0, 0]\n\t[0, 0, 0]\n\n(2)  2 ≥ \n\t[1, 0, 0]\n\t[0, 1, 0]\n\t[0, 0, 1]\n\n"

    @testset "no .txt extension in `filename`" begin
        @test !isfile(filename_txt)

        pretty_print_txt(bell_games, filename)

        @test isfile(filename_txt)

        file_string = read(filename_txt, String)

        @test file_string == match_string

        rm(filename_txt)
    end

    @testset ".txt extension in `filename`" begin
        @test !isfile(filename_txt)

        pretty_print_txt(bell_games, filename_txt)

        @test isfile(filename_txt)

        file_string = read(filename_txt, String)

        @test file_string == match_string

        rm(filename_txt)
    end
end

end


text:
function hello()
    print("Hello, World!")
end

max_stars_repo_name: Rigante-pl/Exercism-Julia

max_issues_repo_name: Rigante-pl/Exercism-Julia

max_forks_repo_name: Rigante-pl/Exercism-Julia


text:
# create recursion coefficient of logistic density numerically
# and compare against analytic solution

using PolyChaos, Test
import LinearAlgebra: norm

function myquad(N::Int64,pos::Bool)
    a,b = rm_laguerre(N+1)
    n,w = golubwelsch(a,b)
    w = w ./ ((1 .+ exp.(-n)).^2)
    pos ? (n,w) : (-n,w)
end

ε, tol = 1e-10, 1e-6
N = 1:3:15
quads = [ n->myquad(n,true); n->myquad(n,false) ]

@testset "Logistic density" begin
    for n in N, disc in [stieltjes, lanczos]
        a,b = mcdiscretization(n,quads;discretization=disc,Nmax=300,gaussquad=true,ε=ε)
        aa,bb = rm_logistic(n)
        @test isapprox(norm(b-bb,Inf),0.;atol=tol)
    end
end

function hrhermite(t)
    exp(-t^2)
end

AB=[[0. 3.];[3. 6.];[6. 9.];[9. Inf]];
quads = [ n->quadgp(hrhermite,AB[i,1],AB[i,2],n;quadrature=fejer) for i=1:size(AB,1)]
α, β = mcdiscretization(40,quads,Nmax=400,gaussquad=false,discretization=lanczos,ε=1e-7)


text:
# Computes K_v(k) and L_v(k) vectors from Agol & Luger (2018), along
# the derivatives with respect to k using recursion.
include("transit_structure.jl")

function dKLv_raise_dk!(t::Transit_Struct{T})  where {T <: Real}
kc = t.kc; k=t.k; k2 = t.k2; kc2 = t.kc2; kck=t.kck; kap=t.kap; Eofk=t.Eofk; Em1mKdm=t.Em1mKdm
# Compute K_v, L_v for 0 <= v <= v_max = l_max+2
# Define k:
# Iterate upwards in v:
v = t.v_max
# Compute I_v via upward iteration on v:
if k2 < 1
# First, compute value for v=0:
  t.Kv[1] = kap
# Try something else:
# Next, iterate upwards in v:
  f0 = kck
# Loop over v, computing K_v and L_v from higher v:
  @inbounds for v=1:t.v_max
    t.Kv[v+1]=((v-0.5)*t.Kv[v]+f0)/v
    f0 *= kc2
  end
  if t.grad
    # Now compute compute derivatives:
    t.dKvdk[1] = 2/kc  # TBD
    @inbounds for v=1:t.v_max
      t.dKvdk[v+1] = k2*t.dKvdk[v] # TB
    end
  end
else # k^2 >= 1
  # Compute v=0
  t.Kv[1] = pi
  @inbounds for v=1:t.v_max
    t.Kv[v+1]=t.Kv[v]*(1-1/(2v))
  end
  if t.grad
    # Derivatives of I_v are zero:
    fill!(t.dKvdk,zero(T))
  end
end
# Need to compute J_v for v=0, 1:
v= 0
if k2 < 1
  # Use elliptic integrals (already computed for s_2/S_1):
  if k2 > 0
    t.Lv[v+1]=2/(3k2*k)*(k2*(3k2-2)*Em1mKdm+k2*Eofk)
    t.Lv[v+2]= 2/(15k2*k)*(k2*(3k2+1)*Eofk-2k2*(1-3k2)*Em1mKdm)
    if t.grad
      t.dLvdk[v+1] = 2/k2*(-Eofk+2*Em1mKdm) # TBD
      t.dLvdk[v+2] = -3*t.Lv[v+2]/k+k2*t.dLvdk[v+1] # TBD
    end
  else
    t.Lv[v+1]= 0.0
    t.Lv[v+2]= 0.0
    if t.grad
      t.dLvdk[v+1] = 0.0
      t.dLvdk[v+2] = 0.0
    end
  end
else # k^2 >=1
  k2inv = inv(k2)
  t.Lv[v+1]=2/3*((3-2*k2inv)*Eofk+k2inv*Em1mKdm)
  t.Lv[v+2]=2*((3+k2inv)*Em1mKdm+(6-2*k2inv)*Eofk)/15
  if t.grad
    t.dLvdk[v+1] = 2/(k2*k)*(2*Eofk-Em1mKdm) # TBD
    t.dLvdk[v+2] = -3*t.Lv[v+2]/k+k2*t.dLvdk[v+1] # TBD
  end
end
@inbounds for v=2:t.v_max
#  t.Lv[v+1] = (2*(2v+(1-v)*k2)*t.Lv[v]+(k2-1)*(2v-3)*t.Lv[v-1])/(2v+3)
  t.Lv[v+1] = k2*(2*(2v*inv(k2)+(1-v))*t.Lv[v]+(1-inv(k2))*(2v-3)*t.Lv[v-1])/(2v+3)
  if t.grad
    t.dLvdk[v+1] = -3*t.Lv[v+1]/k+k2*t.dLvdk[v] # TBD
  end
end
return
end

max_stars_repo_name: tagordon/limbdark

max_issues_repo_name: tagordon/limbdark

max_forks_repo_name: tagordon/limbdark


text:
"""
    import_data(
        path::String
        )

Import species collections `Z` and compositions `P`.

# Arguments:
- `path::String` : location of data
- `p::Float64` : optional partition
"""
function import_data(
            path::String,
            p::Float64=1.0
            )
    P = readdlm(path,',',Float64,'\n')
    Z = (x-> x>0 ? 1 : 0 ).(P)
    reps = unique([ Set(findall([ zz == z for z in eachcol(Z) ])) for zz in eachcol(Z) ])
    inx = [ rand([r...]) for r in reps]
    P = P[:,inx]
    Z = Z[:,inx]
    P = hcat([ normalize(p,1) for p in eachcol(P)]...)
    Z = hcat([ normalize(z,1) for z in eachcol(Z)]...)
    if p<1.0
        (Z,P),_ = splitobs(shuffleobs((Z,P)),at=p)
    end

    return Z, P
end

"""
    split_data(
        Z::Array{Float64},
        P::Array{Float64},
        p::Float64,
        q::Float64
    )

Split data into train, validation and test sets.

# Arguments:
- `Z::Array{Float64}` : collections
- `P::Array{Float64}` : compositions
- `p::Float64` : percentage for training
- `q::Float64` : partition into training (`q`%) and validation (`1-q`%) sets
"""
function split_data(
                Z::Array{Float64},
                P::Array{Float64},
                p::Float64=1.0,
                q::Float64=0.8
                )
    (ztrn,ptrn),(ZT,PT) = splitobs(shuffleobs((Z,P)),at = q)
    (ztrn,ptrn),_ = splitobs(shuffleobs((ztrn,ptrn)),at = p - 1e-5)
    (zval,pval),(ztst,ptst) = splitobs(shuffleobs((ZT,PT)),at = 1/2)
    return ztrn,ptrn,zval,pval,ztst,ptst
end


text:
using Documenter
using Bukdu
using Bukdu.HTML5.Form

makedocs(
    build = joinpath(@__DIR__, "local" in ARGS ? "build_local" : "build"),
    modules = [Bukdu],
    clean = false,
    format = Documenter.HTML(
        prettyurls = !("local" in ARGS),
        assets = ["assets/custom.css"],
    ),
    sitename = "Bukdu.jl 🌌",
    authors = "WooKyoung Noh",
    pages = Any[
        "Home" => "index.md",
        "server" => "server.md",           # Bukdu.start Bukdu.stop
        "controllers" => "controllers.md", # ApplicationController Conn
                                           # Assoc
        "renders" => "renders.md",         # Render
                                           # render
        "routes" => "routes.md",           # routes pipeline
                                           # get post delete patch put resources
                                           # Routing.reset!
        "plugs" => "plugs.md",             # Plug plug
        "Actions" => "Actions.md",         # Actions
        "HTML5.Form" => "HTML5/Form.md",   # change
                                           # form_for text_input submit
        "Changeset" => "changeset.md",     # Changeset
        "CLI" => "CLI.md",                 # CLI
        "System" => "System.md",           # System
        "Router" => "Router.md",           # Router
        "Utils" => "Utils.md",             # Utils
    ],
)


text:
if !isdefined(:number)
    number() = rand(0:100)
end

facts("Call counting") do
    a = Signal(0)
    b = Signal(0)

    c = map(+, a, b)
    d = merge(a, b)
    e = map(+, a, map(x->2x, a)) # Both depend on a
    f = map(+, a, b, c, e)

    count = s -> foldp((x, y) -> x+1, 0, s)

    ca = count(a)
    cb = count(b)
    cc = count(c)
    cd = count(d)
    ce = count(e)
    cf = count(f)

    for i in 1:100
        push!(a, number())
        step()
        push!(b, number())
        step()

        @fact ca.value --> i
        @fact cb.value --> i
        @fact cc.value --> 2i
        @fact cd.value --> 2i
        @fact ce.value --> i
        @fact cf.value --> 2i
    end

end


text:
#----------------------------------------------------------------------------------------------------------------------------------------
# Implements solution to the diffusion equation in the steady-state and time-domain through a turbid parallelepiped [1]

# [1] Alwin Kienle, "Light diffusion through a turbid parallelepiped," J. Opt. Soc. Am. A 22, 1883-1888 (2005)
#----------------------------------------------------------------------------------------------------------------------------------------

#--------------------------------------
# Steady-State Fluence
#--------------------------------------
"""
    fluence_DA_paralpip_CW(μa, μsp; n_ext = 1.0, n_med = 1.0, rd = [4.0, 5.0, 0.0], rs = [5.0, 5.0], L = [10.0, 10.0, 10.0], xs = 10)

Compute the steady-state fluence in a parallelepiped [lx, ly, lz].

# Arguments
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Optional arguments
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `rd`: target location within medium [x,y,z] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `rs`: the location of the source [x,y] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `L`: the dimensions [lx, ly, lz] of the parallelepiped
- `xs`: the number of sources to compute in the series

# Examples
julia> `fluence_DA_paralpip_CW(0.1, 10.0, n_ext = 1.0, n_med = 1.0, rd = [24.0, 25.0, 0.0], rs = [25.0,25.0], L = [50.0,50.0,50.0], xs = 20)`
"""
function fluence_DA_paralpip_CW(μa, μsp; n_ext = 1.0, n_med = 1.0, rd = [4.0, 5.0, 0.0], rs = [5.0, 5.0], L = [10.0, 10.0, 10.0], xs = 10)
    params = DiffusionKernelParams(μsp, n_med, n_ext)
    μeff = sqrt(3 * μa * μsp)

    return _kernel_fluence_DA_paralpip_CW(μeff, params.D, params.zb, rd, rs, L, params.z0, xs)
end

function _kernel_fluence_DA_paralpip_CW(μeff, D, zb, rd, rs, L, z0, xs)
    x = rd[1]
    y = rd[2]
    z = rd[3]
    xu = rs[1]
    yu = rs[2]
    lx = L[1]
    ly = L[2]
    lz = L[3]

    ϕ = zero(eltype(D))

    for l in -xs:xs, m in -xs:xs, n in -xs:xs
    	ϕ += _sources_DA_paralpip_CW(l, m, n, x, y, z, lx, ly, lz, xu, yu, z0, zb, μeff)
    end

    return ϕ / (4 * π * D)
end
function _sources_DA_paralpip_CW(l, m, n, x, y, z, lx, ly, lz, xu, yu, z0, zb, μeff)

    x1l = 2 * l * lx + 4 * l * zb + xu
    x2l = 2 * l * lx + (4 * l - 2) * zb - xu
    y1m = 2 * m * ly + 4 * m * zb + yu
    y2m = 2 * m * ly + (4 * m - 2) * zb - yu
    z1n = 2 * n * lz + 4 * n * zb + z0
    z2n = 2 * n * lz + (4 * n - 2) * zb - z0

    r1 = sqrt((x - x1l)^2 + (y - y1m)^2 + (z - z1n)^2)
    r2 = sqrt((x - x1l)^2 + (y - y1m)^2 + (z - z2n)^2)
    r3 = sqrt((x - x1l)^2 + (y - y2m)^2 + (z - z1n)^2)
    r4 = sqrt((x - x1l)^2 + (y - y2m)^2 + (z - z2n)^2)
    r5 = sqrt((x - x2l)^2 + (y - y1m)^2 + (z - z1n)^2)
    r6 = sqrt((x - x2l)^2 + (y - y1m)^2 + (z - z2n)^2)
    r7 = sqrt((x - x2l)^2 + (y - y2m)^2 + (z - z1n)^2)
    r8 = sqrt((x - x2l)^2 + (y - y2m)^2 + (z - z2n)^2)

    ϕ  = exp(-μeff * r1) / r1 - exp(-μeff * r2) / r2 - exp(-μeff * r3) / r3 + exp(-μeff * r4) / r4
    ϕ += -exp(-μeff * r5) / r5 + exp(-μeff * r6) / r6 + exp(-μeff * r7) / r7 - exp(-μeff * r8) / r8

    return ϕ
end

#--------------------------------------
# Time-Domain Fluence
#--------------------------------------
"""
    fluence_DA_paralpip_TD(t, μa, μsp; n_ext = 1.0, n_med = 1.0, rd = [4.0, 5.0, 0.0], rs = [5.0, 5.0], L = [10.0, 10.0, 10.0], xs = 10)

Compute the time-domain fluence in a parallelepiped [lx, ly, lz].

# Arguments
- `t`: the time vector (ns).
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Optional arguments
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `rd`: target location within medium [x,y,z] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `rs`: the location of the source [x,y] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `L`: the dimensions [lx, ly, lz] of the parallelepiped
- `xs`: the number of sources to compute in the series

# Examples
```
julia> fluence_DA_paralpip_TD(0.5, 0.1, 10.0, n_ext = 1.0, n_med = 1.0, rd = [24.0, 25.0, 0.0], rs = [25.0,25.0], L = [50.0,50.0,50.0], xs = 20)
```
"""
function fluence_DA_paralpip_TD(t, μa, μsp; n_ext = 1.0, n_med = 1.0, rd = [4.0, 5.0, 0.0], rs = [5.0, 5.0], L = [10.0, 10.0, 10.0], xs = 10)
    params = DiffusionKernelParams(μsp, n_med, n_ext)
    return map(t -> _kernel_fluence_DA_paralpip_TD(params.D, params.ν, t, μa, params.zb, rd, rs, L, params.z0, xs), t)
end
function _kernel_fluence_DA_paralpip_TD(D, ν, t, μa, zb, rd, rs, L, z0, xs)
    x = rd[1]
    y = rd[2]
    z = rd[3]
    xu = rs[1]
    yu = rs[2]
    lx = L[1]
    ly = L[2]
    lz = L[3]

    tmp = 4 * D * ν * t
    ϕ1 = ν * exp(-μa * ν * t) / ((π * tmp)^(3/2))
    ϕ2 = zero(eltype(ϕ1))
    ϕ3 = zero(eltype(ϕ1))
    ϕ4 = zero(eltype(ϕ1))

    for m in -xs:xs
    	ϕ2, ϕ3, ϕ4 = _sources_DA_paralpip_TD!(ϕ2, ϕ3, ϕ4, m, zb, lx, ly, lz, xu, yu, z0, x, y, z, tmp)
    end

    return ϕ1 * ϕ2 * ϕ3 * ϕ4
end
@inline function _sources_DA_paralpip_TD!(ϕ2, ϕ3, ϕ4, m, zb, lx, ly, lz, xu, yu, z0, x, y, z, tmp)
    tmp1 = 4 * m * zb
    tmp2 = (4 * m - 2) * zb
    tmp3 = 2 * m

    x1l = tmp3 * lx + tmp1 + xu
    x2l = tmp3 * lx + tmp2 - xu
    y1m = tmp3 * ly + tmp1 + yu
    y2m = tmp3 * ly + tmp2 - yu
    z1n = tmp3 * lz + tmp1 + z0
    z2n = tmp3 * lz + tmp2 - z0

    ϕ2 += exp(-(x - x1l)^2 / tmp) - exp(-(x - x2l)^2 / tmp)
    ϕ3 += exp(-(y - y1m)^2 / tmp) - exp(-(y - y2m)^2 / tmp)
    ϕ4 += exp(-(z - z1n)^2 / tmp) - exp(-(z - z2n)^2 / tmp)

    return ϕ2, ϕ3, ϕ4
end

### TD Reflectance ###
"""
    flux_DA_paralpip_TD(t, μa, μsp; n_ext = 1.0, n_med = 1.0, rd = [4.0, 5.0, 0.0], rs = [5.0, 5.0], L = [10.0, 10.0, 10.0], xs = -10:10)

Compute the time-domain flux from a parallelepiped [lx, ly, lz] applying Fick's law to Eqn. 6 Kienle 05.

# Arguments
- `t`: the time vector (ns).
- `μa`: absorption coefficient (cm⁻¹)
- `μsp`: reduced scattering coefficient (cm⁻¹)

# Optional arguments
- `n_ext`: the boundary's index of refraction (air or detector)
- `n_med`: the sample medium's index of refraction
- `rd`: target location within medium [x,y,z] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `rs`: the location of the source [x,y] with origin x,y = 0 at corner of parallelepiped; z assumed to be 0
- `L`: the dimensions [lx, ly, lz] of the parallelepiped
- `xs`: the number of sources to compute in the series

# Examples
julia> `flux_DA_paralpip_TD(0.1:0.1:2.0, 0.1, 10.0)`
"""
function flux_DA_paralpip_TD(t, μa, μsp; n_ext = 1.0, n_med = 1.0, rd = [4.0, 5.0, 0.0], rs = [5.0, 5.0], L = [10.0, 10.0, 10.0], xs = 10)
    params = DiffusionKernelParams(μsp, n_med, n_ext)
    return map(t -> _flux_DA_paralpip_TD(t, μa, params.D, params.ν, params.z0, params.zb, rd, rs, L, xs), t)
end

function _flux_DA_paralpip_TD(t, μa, D, ν, z0, zb, rd, rs, L, xs)
    x = rd[1]
    y = rd[2]
    xu = rs[1]
    yu = rs[2]
    lx = L[1]
    ly = L[2]
    lz = L[3]

    if rd[3] == zero(eltype(μa))
        f = _refl_DA_paralpip_TD
    elseif rd[3] == L[3]
        f =  _trans_DA_paralpip_TD
    else
        @assert z == zero(eltype(z)) || z == s
    end

    Rt1 = exp(-μa * ν * t) / (2 * (4 * π * D * ν)^(3/2) * t^(5/2))
    Rt2 = zero(eltype(Rt1))
    Rt3 = zero(eltype(Rt1))
    Rt4 = zero(eltype(Rt1))
    tmp = 4 * D * ν * t

    for m in -xs:xs
        x1l = 2 * m * lx + 4 * m * zb + xu
        x2l = 2 * m * lx + (4 * m - 2) * zb - xu
        y1m = 2 * m * ly + 4 * m * zb + yu
        y2m = 2 * m * ly + (4 * m - 2) * zb - yu
        z1n = 2 * m * lz + 4 * m * zb + z0
        z2n = 2 * m * lz + (4 * m - 2) * zb - z0

        Rt2 += exp(-(x - x1l)^2 / tmp) - exp(-(x - x2l)^2 / tmp)
        Rt3 += exp(-(y - y1m)^2 / tmp) - exp(-(y - y2m)^2 / tmp)
        Rt4 += f(z1n, tmp, z2n, lz)
    end

    Rt1 *= Rt2 * Rt3 * Rt4
	return Rt1
end

_refl_DA_paralpip_TD(z1n, tmp, z2n, lz) = z1n * exp(-z1n^2 / tmp) - z2n * exp(-z2n^2 / tmp)
_trans_DA_paralpip_TD(z1n, tmp, z2n, lz) = (lz - z1n) * exp(-(lz - z1n)^2 / tmp) - (lz - z2n) * exp(-(lz - z2n)^2 / tmp)

max_issues_repo_name: heltonmc/LightPropagation


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

const ZERO = false

sum_ladjs(ladjs::Real) = ladjs
sum_ladjs(ladjs::AbstractVector{<:Real}) = sum(ladjs)
sum_ladjs(ladjs::AbstractMatrix{<:Real}) = vec(sum(ladjs, dims = 1))'

function similar_zeros(x::AbstractArray{<:Number}, sz::Dims)
    r = similar(x, sz)
    fill!(r, zero(eltype(x)))
end

function rrule(::typeof(similar_zeros), x::AbstractArray{<:Number}, sz::Dims)
    _similar_zeros_pullback(ΔΩ) = (NoTangent(), ZeroTangent(), NoTangent())
    return similar_zeros(x, sz), _similar_zeros_pullback
end

function similar_fill(value::T, x::AbstractArray{<:U}, sz::Dims) where {T<:Number,U<:Number}
    R = promote_type(T, U)
    v = convert(R, value)
    r = similar(x, R, sz)
    fill!(r, value)
end

function rrule(::typeof(similar_fill), value::T, x::AbstractArray{<:U}, sz::Dims) where {T<:Number,U<:Number}
    _similar_zeros_pullback(ΔΩ) = (NoTangent(), sum(unthunk(ΔΩ)), ZeroTangent(), NoTangent())
    return similar_zeros(x, sz), _similar_zeros_pullback
end


text:
module TestPDA
    using Base.Test
    using DataArrays

    p = @pdata [9, 9, 8, NA, 1, 1]
    pcopy = copy(p)
    @assert levels(p) == [1, 8, 9]
    @assert levels(setlevels(p, ["a", "b", "c"])) == ["a", "b", "c"]
    @assert dropna(setlevels(p, (@data ["a", "b", NA]))) == ["b", "a", "a"]
    @assert dropna(setlevels(p, (@data ["a", "b", "a"]))) == ["a", "a", "b", "a", "a"]
    @assert levels(setlevels(p, (@data ["a", "b", "a"]))) == ["a", "b"]
    @assert levels(setlevels(p, [1 => 111])) == [111, 8, 9]
    @assert levels(setlevels(p, [1 => 111, 8 => NA])) == [111, 9]
    @assert levels(PooledDataArray(p, [9, 8, 1])) == [9, 8, 1]
    @assert levels(PooledDataArray(p, [9, 8])) == [9, 8]
    @assert dropna(PooledDataArray(p, [9, 8])) == [9, 9, 8]
    @assert levels(PooledDataArray(p, levels(p)[[3,2,1]])) == [9,8,1]
    v = [1:6]
    @assert isequal(p, reorder(p))
    # @assert levels(reorder(p, v)) == [9,8,1]
    @assert isequal(p, pcopy)

    @assert levels(setlevels!(copy(p), [10,80,90])) == [10, 80, 90]
    @assert levels(setlevels!(copy(p), [1,8,1])) == [1, 8]
    @assert levels(setlevels!(copy(p), (@data [1, 8, NA]))) == [1, 8]
    @assert levels(setlevels!(copy(p), [1,8,9, 10])) == [1, 8, 9, 10]
    @assert levels(setlevels!(copy(p), [1 => 111])) == [111, 8, 9]
    @assert levels(setlevels!(copy(p), [1 => 111, 8 => NA])) == [111, 9]

    pp = PooledDataArray(Any[])
    @assert length(pp) == 0
    @assert length(levels(pp)) == 0

    # Test explicitly setting refs type
    testarray = [1, 1, 2, 2, 0, 0, 3, 3]
    testdata = @data [1, 1, 2, 2, 0, 0, 3, 3]
    for t in Any[testarray, testdata]
        for R in [Uint8, Uint16, Uint32, Uint64]
            @assert eltype(PooledDataArray(t, R).refs) == R
            @assert eltype(PooledDataArray(t, [1,2,3], R).refs) == R
            @assert eltype(PooledDataArray(t, [1,2,3], t .== 0, R).refs) == R
        end
    end

    pcopy = copy(p)
    @assert levels(append!(pcopy, @pdata [4, NA, 6, 5])) == [1, 8, 9, 4, 5, 6]

    x = PooledDataArray([9, 9, 8])
    y = PooledDataArray([1, 9, 3, 2, 2])
    @assert append!(x, y) == [9, 9, 8, 1, 9, 3, 2, 2]

    x = PooledDataArray([9, 9, 8])
    y = [1, 9, 3, 2, 2]
    @assert append!(x, y) == [9, 9, 8, 1, 9, 3, 2, 2]

    # convert methods
    for from in (@pdata(ones(5, 5)), @data(ones(5, 5)), ones(5, 5))
        for (to, totype) in ((PooledDataArray{Float32,Uint16,2}, PooledDataArray{Float32,Uint16,2}),
                             (PooledDataArray{Float32,Uint32,2}, PooledDataArray{Float32,Uint32,2}),
                             (PooledDataArray{Float64,Uint16,2}, PooledDataArray{Float64,Uint16,2}),
                             (PooledDataArray{Float64,Uint32,2}, PooledDataArray{Float64,Uint32,2}),
                             (PooledDataArray{Float32,Uint16}, PooledDataArray{Float32,Uint16,2}),
                             (PooledDataArray{Float32,Uint32}, PooledDataArray{Float32,Uint32,2}),
                             (PooledDataArray{Float64,Uint16}, PooledDataArray{Float64,Uint16,2}),
                             (PooledDataArray{Float64,Uint32}, PooledDataArray{Float64,Uint32,2}),
                             (PooledDataArray{Float32}, PooledDataArray{Float32,Uint32,2}),
                             (PooledDataArray{Float64}, PooledDataArray{Float64,Uint32,2}),
                             (PooledDataArray, PooledDataArray{Float64,Uint32,2}))
            rettype = typeof(convert(to, from))
            if rettype != totype
                error("convert($to, ::$(typeof(from))) returned $rettype (expected $totype)")
            end
        end
    end

    da = convert(DataArray, @pdata(ones(5, 5)))
    @assert isequal(da, @data(ones(5, 5)))
    @assert typeof(da) == DataArray{Float64,2}
    da = convert(DataArray{Float32}, @pdata(ones(5, 5)))
    @assert isequal(da, @data(ones(Float32, 5, 5)))
    @assert typeof(da) == DataArray{Float32,2}
    da = convert(DataArray{Float32,2}, @pdata(ones(5, 5)))
    @assert isequal(da, @data(ones(Float32, 5, 5)))
    @assert typeof(da) == DataArray{Float32,2}
end


text:
macro assert0(expression)
	:(@assert all($expression .<= 1e-12))
end

macro unpack_dictionary(expression)
	lines = []
	dict = eval(expression)
	for (key, value) in dict
		:(local $key = $value)
	end
end

function fill_dict(;kwargs...)
	return Dict(kwargs)
end

function fill_dict!(dict; kwargs...)
	merge!(dict, Dict(kwargs))
end

function array_transpose(array)
	index = [i for i in 1:ndims(array)]
	index[1], index[2] = 2, 1
	return permutedims(array, index)
end

function rotate_along_dimension(rotator, array, n)
	index = [i for i in 1:ndims(array)]
	index[1], index[n] = n, 1
	temp = permutedims(array, index)
	S1 = size(temp)
	A = zeros()
	for i in CartesianRange(S1[2:end])
		@show i
		A[:,i] = rotator*temp[:,i]
	end
	return A
end

function share_within_dimension(X, i)
	return X ./ sum(X,i)
end

function eigen_share(A)
	# the solution to x = Ax with x'1 = 1
	D, V = eig(A)
	i = find(abs.(D - 1.0) .< 1e-6)[1]
	return V[:,i] ./ sum(V[:,i])
end

function non_random_variable(y, t)
	# array coersion is going to take care of rest
	B = y[:,:,:,t]
	return cat(ndims(B)+1, B)
end

function remove_shock!(parameters, symbol)
	for t=1:length(parameters[:A_njs])
		shock_to_remove = exp.(parameters[symbol][t])
		parameters[:A_njs][t] = parameters[:A_njs][t] ./ shock_to_remove
	end
end


text:
using Distributions, Statistics, LinearAlgebra

import Base: size, length, hcat, +, -, fill!

import Statistics: mean, var, std

export EnsembleState!, EnsembleState, deviation!, deviation, cut, fill!, std

"""
    EnsembleState

An structure  for Ensemble.

Fields:

- 'N' : Number of ensemble members"

- 'S' : Array of the different members
"""
mutable struct EnsembleState{N,TS}
    " Array of the different ensemble members"
    S::Array{TS,1}
end

# Return an ensemble of P members where each member is an
# array of type TS and dimension NS

function EnsembleState(N::Int, u)
    TS = typeof(u)
    return EnsembleState{N,TS}([zero(u) for i = 1:N])
end

function EnsembleState(States::Array{TS,1}) where {TS}
    N = length(States)
    return EnsembleState{N, TS}(States)
end

size(ENS::EnsembleState{N, TS})  where {N, TS} = N, size(ENS.S[1])

length(ENS::EnsembleState{N, TS})  where {N, TS} = N

"""
    Return the mean of all the ensemble members
"""
mean(ENS::EnsembleState{N, TS})  where {N, TS} = mean(ENS.S)
std(ENS::EnsembleState{N, TS})  where {N, TS} = std(ENS.S)

function deviation!(ENS::EnsembleState{N, TS})  where {N, TS}
    S̄ = deepcopy(mean(ENS))
     for s in ENS.S
         s .-= S̄
     end
     return ENS
 end

function deviation(ENSfluc::EnsembleState{N, TS}, ENS::EnsembleState{N, TS})  where {N, TS}
 S̄ = deepcopy(mean(ENS))
 ENSfluc.S .= deepcopy(ENS.S)
  for s in ENSfluc.S
      s .-= S̄
  end
  return ENSfluc
end

function deviation(tabfluc::Array{TS,2}, ENS::EnsembleState{N, TS})  where {N, TS}
 S̄ = deepcopy(mean(ENS))
 tabfluc .= hcat(deepcopy(ENS))
 tabfluc .-= S̄
   return tabfluc
end

"""
    Extend definition of hcat to EnsembleState
"""
function hcat(ENS::EnsembleState{N, TS})  where {N, TS}
   return hcat(ENS.S...)
end

"""
    Cut an array along the different columns and create an EnsembleState
    variable with these columns
"""
function cut(A::AbstractMatrix{TR}) where {TR}
    # Get size of A = (length of state vector, number of ensemble members)
        NS, N = size(A)
        B = deepcopy(A)
    # Allocate space
    ENS = EnsembleState(N, zeros(NS))
    for i  = 1:N
        ENS.S[i] = B[:,i]
    end
    return ENS
end

function fill!(ENS::EnsembleState{N, TS}, A::TS)  where {N, TS}
    B = deepcopy(A)
    for s in ENS.S
        s .= B
    end
    return ENS
end

"""
    Define addition of two EnsembleState
"""
function (+)(A::EnsembleState{N, TS}, B::EnsembleState{N, TS}) where {N, TS}
    C = deepcopy(A)
    C.S .+= B.S
    return C
end

"""
    Define addition of an Array and an EnsembleState
"""
function (+)(A::EnsembleState{N, TS}, B::TS) where {N, TS}
    C = deepcopy(A)
    for s in C.S
        s .+=B
    end
    return C
end

"""
    Define substraction of two EnsembleState
"""
function (-)(A::EnsembleState{N, TS}, B::EnsembleState{N, TS}) where {N, TS}
    C = deepcopy(A)
    C.S .-= B.S
    return C
end

"""
    Define substraction of an Array from an EnsembleState
"""
function (-)(A::EnsembleState{N, TS}, B::TS) where {N, TS}
    C = deepcopy(A)
    for s in C.S
        s .-=B
    end
    return C
end

max_stars_repo_name: mleprovost/EnKF

max_issues_repo_name: mleprovost/EnKF

max_forks_repo_name: mleprovost/EnKF


text:
module ParallelKMeans

using StatsBase
import MLJModelInterface
import Base.Threads: @spawn
import Distances

const MMI = MLJModelInterface

include("seeding.jl")
include("kmeans.jl")
include("lloyd.jl")
include("hamerly.jl")
include("elkan.jl")
include("mlj_interface.jl")

export kmeans
export Lloyd, Hamerly, Elkan

end # module


text:
include("../src/Embedding.jl")

TimeSteps = 100
EmbDims = 128
Vocab = Integer(1e4)

Embed = Embedding(Vocab, EmbDims)
X = map(a->Int32.(a), [[3,7,3,5], [8,6,4,5,2,1], [10,7,5], [11,3]])
@time output = Embed(X)


text:
module CellMaps

using Test
using CellwiseValues
using CellwiseValues.CellValues: _test_iter_cell_value

export CellMap
export IterCellMap
export IndexCellMap
import CellwiseValues: evaluate
export test_iter_cell_map
export test_index_cell_map
export test_index_cell_map_with_index_arg

import Base: iterate
import Base: length
import Base: size
import Base: getindex

"""
Abstract object that traverses a set of cells and at every cell returns a
`Map{S,M,T,N}`
"""
const IterCellMap{S,M,T,N,R<:Map{S,M,T,N}} = IterCellValue{R}

"""
Abstract array indexed by cells that returns a `Map{S,M,T,N}`
"""
const IndexCellMap{S,M,T,N,C,R<:Map{S,M,T,N}} = IndexCellValue{R,C}

"""
Abstract object that for a given cell index returns a `Map{S,M,T,N}`
"""
const CellMap{S,M,T,N} = Union{IterCellMap{S,M,T,N},IndexCellMap{S,M,T,N}}

# evaluation

"""
Return the cellwise maps of a `CellMap` on a cellwise set of points
"""
function evaluate(cm::CellMap{S,M,T,N},ca::CellArray{S,M}) where {S,M,T,N}
  IterCellMapValue(cm,ca)
end

function evaluate(cm::IndexCellMap{S,M,T,N},ca::IndexCellArray{S,M}) where {S,M,T,N}
  IndexCellMapValue(cm,ca)
end

struct IterCellMapValue{T,N,V,A} <: IterCellArray{T,N,CachedArray{T,N,Array{T,N}}}
  cm::V
  ca::A
end

function IterCellMapValue(cm::CellMap{S,M,T,N},ca::CellArray{S,M}) where {S,M,T,N}
  V = typeof(cm)
  A = typeof(ca)
  IterCellMapValue{T,N,V,A}(cm,ca)
end

length(v::IterCellMapValue) = length(v.ca)

@inline function iterate(v::IterCellMapValue{T,N}) where {T,N}
  cmnext = iterate(v.cm)
  canext = iterate(v.ca)
  r = CachedArray(T,N)
  _iterate(cmnext,canext,r)
end

@inline function iterate(v::IterCellMapValue,state)
  cmstate, castate, r = state
  cmnext = iterate(v.cm,cmstate)
  canext = iterate(v.ca,castate)
  _iterate(cmnext,canext,r)
end

@inline function _iterate(cmnext,canext,r)
  if cmnext === nothing; return nothing end
  if canext === nothing; return nothing end
  m, cmstate = cmnext
  a, castate = canext
  s = return_size(m,size(a))
  setsize!(r,s)
  evaluate!(m,a,r)
  state = (cmstate,castate,r)
  (r, state)
end

struct IndexCellMapValue{T,N,V,A} <: IndexCellArray{T,N,CachedArray{T,N,Array{T,N}},1}
  cm::V
  ca::A
  r::CachedArray{T,N,Array{T,N}}
end

function IndexCellMapValue(cm::IndexCellMap{S,M,T,N},ca::IndexCellArray{S,M}) where {S,M,T,N}
  V = typeof(cm)
  A = typeof(ca)
  r = CachedArray(T,N)
  IndexCellMapValue{T,N,V,A}(cm,ca,r)
end

length(v::IndexCellMapValue) = length(v.ca)

size(v::IndexCellMapValue) = (length(v),)

function getindex(v::IndexCellMapValue,i::Integer)
  m = v.cm[i]
  a = v.ca[i]
  s = return_size(m,size(a))
  setsize!(v.r,s)
  evaluate!(m,a,v.r)
  v.r
end

# Testers

function test_iter_cell_map(
  m::CellMap{S,M,T,N},
  a::CellArray{S,M},
  b::AbstractArray{<:AbstractArray{T,N}}) where {S,M,T,N}

  @test length(m) == length(b)
  @test length(m) == length(a)

  c = evaluate(m,a)
  _test_iter_cell_value(c,b)

  for (mi,ai,bi) in zip(m,a,b)
    @assert isa(mi,Map{S,M,T,N})
    @assert evaluate(mi,ai) ≈ bi
    @assert typeof(mi) == eltype(m)
  end

end

function test_index_cell_map(
  m::IndexCellMap{S,M,T,N},
  a::CellArray{S,M},
  b::AbstractArray{<:AbstractArray{T,N}}) where {S,M,T,N}

  @test size(m) == size(b)

  @test IndexStyle(m) == IndexStyle(b)

  test_iter_cell_map(m,a,b)

  for (i,ai) in enumerate(a)
    mi = m[i]
    bi = b[i]
    @assert isa(mi,Map{S,M,T,N})
    @assert evaluate(mi,ai) ≈ bi
    @assert typeof(mi) == eltype(m)
  end

end

function test_index_cell_map_with_index_arg(
  m::IndexCellMap{S,M,T,N},
  a::IndexCellArray{S,M},
  b::AbstractArray{<:AbstractArray{T,N}}) where {S,M,T,N}

  test_index_cell_map(m,a,b)

  c = evaluate(m,a)
  @test isa(c,IndexCellArray)

  @test length(c) == length(m)
  @test length(c) == length(a)

  for i in 1:length(c)
    mi = m[i]
    bi = b[i]
    ci = c[i]
    @assert bi ≈ ci
  end

  for i in 1:length(m)
    mi = m[i]
    ai = a[i]
    bi = b[i]
    @assert isa(mi,Map{S,M,T,N})
    @assert evaluate(mi,ai) ≈ bi
    @assert typeof(mi) == eltype(m)
  end

end

end # module CellMaps


text:
dispersion_equation(ω::AbstractFloat, source::AbstractSource, material::Material; kws...) = dispersion_equation(ω, source.medium, material.species, setupsymmetry(source,material); kws...)

function dispersion_equation(ω::T, medium::PhysicalMedium{T,Dim}, species::Species{T,Dim}, symmetry::AbstractSetupSymmetry = PlanarSymmetry{Dim}();
        basis_order = 3 * Int(round(maximum(outer_radius.(species)) * ω / abs(medium.c) )) + 1,
        tol::T = 1e-4, low_tol::T = max(1e-4, tol), kws...
    ) where {T<:Number, Dim}

    # low_tol: a tolerance used for a first pass with time_limit

    MM = eigensystem(ω, medium, species, symmetry; basis_order=basis_order, kws... )

    # the constraint uses keff_vec[2] < -low_tol to better specify solutions where imag(k_effs)~0 and imag(k_effs)<0
    constraint(k_eff::Complex{T}) = (imag(k_eff) < -low_tol) ? (-one(T) + exp(-T(100.0) * imag(k_eff))) : zero(T)

    function detMM(k_eff::Complex{T})
        constraint(k_eff) + abs(det(MM(k_eff)))
    end

    return detMM
end

function dispersion_complex(ω::T, medium::PhysicalMedium{T,Dim}, species::Species{T,Dim}, symmetry::AbstractSetupSymmetry = PlanarSymmetry{Dim}();
        basis_order = 3 * Int(round(maximum(outer_radius.(species)) * ω / abs(medium.c) )) + 1,
        tol::T = 1e-4, kws...
    ) where {T<:Number, Dim}

    MM = eigensystem(ω, medium, species, symmetry; basis_order=basis_order, kws... )
    detMM(k_eff::Complex{T})::Complex{T} = det(MM(k_eff))

    return detMM
end


text:
using MurmurHash3: mmhash128_a

function Base.hash(x::ShortString, h::UInt)
    h += Base.memhash_seed
    last(mmhash128_a(sizeof(x), bswap(x.size_content), h%UInt32)) + h
end


text:
export multiple_trajectories

function one_trajectory()::Simulation
    s = simulation_init();
    simulation_constructor_x_v!(s);
    simulation_constructor_nn!(s);
    simulation_constructor_energy(s);
    simulation_constructor_force(s);
    propagate_init!(s);
    simulate!(s);
    return s
end

function multiple_trajectory()
    filepath_parent = "W:\\ALL\\Theory Group\\iesh\\data\\"
    current_time = Dates.now()
    str_time = string(Dates.format(current_time, "yyyy-dd-mm_HH-MM-SS"))
    filepath_run = filepath_parent*curr_vers*"\\"*str_time *"\\"
    mkpath(filepath_run)
    Threads.@threads for traj in 1:numtraj
        println(traj)
        s = one_trajectory()
        hdf5filename = "traj_"*string(traj)*".h5"
        hdf5filepath = filepath_run*hdf5filename
        fid = h5open(hdf5filepath, "w")
        create_dataset(fid, "x_no", s.storage_xno)
        create_dataset(fid, "v", s.storage_vno)
        create_dataset(fid, "adsorbate orbital population", s.storage_aop)
        create_dataset(fid, "orbital population", s.storage_op)
        create_dataset(fid, "electronic excitation energies", s.storage_e)
        create_dataset(fid, "inst temp", s.storage_temp)
        create_dataset(fid, "max hopping probability", s.storage_phop)
        create_dataset(fid, "kinetic energy", s.storage_K)
        create_dataset(fid, "potential energy", s.storage_P)
        create_dataset(fid, "hopping states", s.storage_state)
        create_dataset(fid, "deltaK_hop NO", s.storage_deltaKNO)
        create_dataset(fid, "deltaK_hop Au", s.storage_deltaKAu)
        create_dataset(fid, "hoptimes", s.storage_hoptimes)
        create_dataset(fid, "nf", s.nf)
        create_dataset(fid, "attnum", s.attnum)
        create_dataset(fid, "exnum", s.exnum)
        create_dataset(fid, "temperature", tsurf)
        create_dataset(fid, "e_trans_no",  e_trans_i)
        create_dataset(fid, "e_vib_no", e_vib_i)
        create_dataset(fid, "e_rot_no", e_rot_i)
        create_dataset(fid, "incident angle", phi_inc)
        create_dataset(fid, "trajectories", numtraj)
        create_dataset(fid, "time stepsize", dt)
        create_dataset(fid, "time steps", tsteps)
        create_dataset(fid, "number electrons", Ne)
        create_dataset(fid, "Number of states", Ms)
        close(fid)
    end
end