text:
import Base.diff

# Because diff is already documented in markdown, we should just use MD documentation
"""
    diff(A::OrgDoc, B::OrgDoc)
    diff(A::Component, B::Component)
    diff(A::Vector{<:Component}, B::Vector{<:Component})
Compare two different Org components or ASTs.

The diff method can detect insertions, deletions, and reshuffling of
immediate decendents of a component.
"""
function diff(A::OrgDoc, B::OrgDoc)
    if A == B
        print(stdout, "Identical\n")
    else
        diff(A.contents, B.contents, [])
    end
    nothing
end

function printbreadcrums(path, prefix=nothing, prefixcolor=:default, postfix='\n')
    # Path is a vector of (component, counter, color) tuples
    # counter can be an integer or a (value, bold?, color) tuple
    if length(path) > 0
        isnothing(prefix) || printstyled(prefix, ' ', color=prefixcolor)
        for (component, counter, color) in path
            (component, counter, color) == first(path) ||
                printstyled(" > ", bold=true)
            printstyled(component; color)
            if counter isa Integer
                printstyled('#', counter, bold=true, color=:light_black)
            elseif counter isa Tuple
                printstyled('#', counter[1], bold=counter[2], color=counter[3])
            end
        end
        print(postfix)
    end
end

function printdiffheader(message, aval, bval)
    printstyled("  × ", bold=true, color=:magenta)
    printstyled(message, '\n', color=:magenta)
    if isnothing(aval)
    elseif aval isa Component
        printstyled("    A (Org representation):\n", color=:red)
        org(stdout, aval, 7)
        print('\n')
    else
        printstyled("    A: ", color=:red)
        print(aval, '\n')
    end
    if isnothing(bval)
    elseif bval isa Component
        printstyled("    B (Org representation):\n", color=:green)
        org(stdout, bval, 7)
        print('\n')
    else
        printstyled("    B: ", color=:green)
        print(bval, '\n')
    end
end

function diff(A::Component, B::Component, path::Vector=[], counter=1)
    pathA = vcat(path, (nameof(typeof(A)), counter, :yellow))
    if A == B
    elseif nameof(typeof(A)) != nameof(typeof(B))
        printbreadcrums(pathA)
        printdiffheader("Component type mismatch", typeof(A), typeof(B))
    else
        props = fieldnames(typeof(A))
        for prop in props
            ap, bp = getproperty(A, prop), getproperty(B, prop)
            if ap != bp
                if (ap isa Component && bp isa Component) ||
                    (ap isa Vector{<:Component} && ap isa Vector{<:Component})
                    diff(ap, bp, pathA, 1)
                elseif ap isa Component || ap isa Vector{<:Component}
                    printbreadcrums(pathA)
                    printdiffheader("$(nameof(typeof(A))) $prop type mismatch", typeof(ap), typeof(bp))
                else
                    printbreadcrums(pathA)
                    printdiffheader("$(nameof(typeof(A))) $prop mismatch", string(ap), string(bp))
                end
            end
        end
    end
    nothing
end

function diff(As::Vector{<:Component}, Bs::Vector{<:Component}, path::Vector=[], counter=0)
    # Detect reorderings
    diffmat = [A == B for A in As, B in Bs]
    reorderings = []
    while maximum(diffmat) > 0
        a, b = Tuple(argmax(diffmat))
        if a != b && a > b
            push!(reorderings, (a, b))
        end
        diffmat[a, :] .= 0
        diffmat[:, b] .= 0
    end
    for (a, b) in reorderings
        printbreadcrums(vcat(path, (nameof(typeof(As[a])), (a, false, :red), :yellow)), "swapped", :blue, "")
        printstyled(" <-> ", color=:blue)
        printbreadcrums([(nameof(typeof(Bs[b])), (b, false, :green), :yellow)])
    end
    reordered = reorderings |> Iterators.flatten |> collect
    # Remaining differences
    if As == Bs
    elseif length(As) == length(Bs)
        for i in setdiff(1:length(As), reordered)
            diff(As[i], Bs[i], path, i)
        end
    else
        minlen = min(length(As), length(Bs))
        for i in setdiff(1:minlen, reordered)
            diff(As[i], Bs[i], path, i)
        end
        printdiffheader("Missing components", "$(length(As)) components", "$(length(Bs)) components")
        if length(As) > length(Bs)
            for i in (minlen+1):length(As)
                printbreadcrums(vcat(path, (nameof(typeof(As[i])), i, :red)), "    !", :red)
            end
        else
            for i in (minlen+1):length(Bs)
                printbreadcrums(vcat(path, (nameof(typeof(Bs[i])), i, :green)), "    !", :green)
            end
        end
    end
    nothing
end


text:
"R_ParseVector wrapped by R_tryCatchError. It catches possible R's `stop` calls
which may cause longjmp in c."
function safe_parseVector(st::Ptr{StrSxp}, status::Ref{Cint}, sf::Ptr{S}=sexp(Const.NilValue)) where S<:Sxp
    protect(st)
    protect(sf)
    fcfunction = @cfunction(
        $((_) -> convert(
            Ptr{UnknownSxp},
            ccall((:R_ParseVector,libR),Ptr{UnknownSxp}, (Ptr{StrSxp},Cint,Ptr{Cint},Ptr{UnknownSxp}), st,-1,status,sf))),
        Ptr{UnknownSxp},
        (Ptr{Cvoid}, ))
    ecfunction = @cfunction(
        $((c, _) -> convert(Ptr{UnknownSxp}, c)),
        Ptr{UnknownSxp},
        (Ptr{UnknownSxp}, Ptr{Cvoid}))
    ret = ccall((:R_tryCatchError, libR), Ptr{UnknownSxp},
          (Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),
          fcfunction.ptr,
          C_NULL,
          ecfunction.ptr,
          C_NULL)
    unprotect(2)
    sexp(ret)
end

"A pure julia wrapper of R_ParseVector"
function parseVector(st::Ptr{StrSxp}, status::Ref{Cint}, sf::Ptr{S}=sexp(Const.NilValue)) where S<:Sxp
    protect(st)
    protect(sf)
    val = ccall((:R_ParseVector,libR),Ptr{UnknownSxp},
                (Ptr{StrSxp},Cint,Ptr{Cint},Ptr{UnknownSxp}),
                st,-1,status,sf)
    unprotect(2)
    sexp(val)
end

"Get the R parser error msg for the previous parsing result."
function getParseErrorMsg()
    unsafe_string(cglobal((:R_ParseErrorMsg, libR), UInt8))
end

"Parse a string as an R expression, returning a Sxp pointer."
function rparse_p(st::Ptr{StrSxp}, sf::Ptr{S}=sexp(Const.NilValue))  where S<:Sxp
    protect(st)
    protect(sf)
    status = Ref{Cint}()
    # use R_tryCatchError to evaluate parseVector as parseVector may longjmp
    result = protect(safe_parseVector(st, status, sf))

    try
        if "error" in rcopy(Array, getclass(result))
            throw(RParseError("Error: " * rcopy(result["message"])))
        else
            if status[] == 0
                throw(RParseError())
            elseif status[] == 2
                throw(RParseIncomplete("Error: " * getParseErrorMsg()))
            elseif status[] == 3
                throw(RParseError("Error: " * getParseErrorMsg()))
            elseif status[] == 4
                throw(RParseEOF())
            end
            sexp(result)
        end
    finally
        unprotect(3)
    end
end
rparse_p(st::AbstractString, sf::Ptr{S}=sexp(Const.NilValue)) where S<:Sxp = rparse_p(sexp(st), sf)
rparse_p(s::Symbol, sf::Ptr{S}=sexp(Const.NilValue)) where S<:Sxp = rparse_p(string(s), sf)

"Parse a string as an R expression, returning an RObject."
rparse(st::AbstractString) = RObject(rparse_p(st))

"""
A pure julia wrapper of R_tryEval.
"""
function tryEval(expr::Ptr{S}, env::Ptr{EnvSxp}=sexp(Const.GlobalEnv), status::Ref{Cint}=Ref{Cint}()) where S<:Sxp
    disable_sigint() do
        protect(expr)
        protect(env)
        val = ccall((:R_tryEval,libR),Ptr{UnknownSxp},(Ptr{S},Ptr{EnvSxp},Ref{Cint}),expr,env,status)
        unprotect(2)
        val
    end
end

"""
Evaluate an R symbol or language object (i.e. a function call) in an R
try/catch block, returning a Sxp pointer.
"""
function reval_p(expr::Ptr{S}, env::Ptr{EnvSxp}=sexp(Const.GlobalEnv)) where S<:Sxp
    status = Ref{Cint}()
    val = tryEval(expr, env, status)
    succeed = status[] == 0
    handle_eval_stdout()
    handle_eval_stderr(as_warning=succeed)
    # always throw an error if status is not zero
    !succeed && throw(REvalError())
    sexp(val)
end

"""
Evaluate an R expression array iteratively. If `throw_error` is `false`,
the error message and warning will be thrown to stderr.
"""
function reval_p(expr::Ptr{ExprSxp}, env::Ptr{EnvSxp})
    local val
    protect(expr)
    protect(env)
    try
        for e in expr
            val = reval_p(e, env)
        end
    finally
        unprotect(2)
    end
    # set .Last.value
    if env == Const.GlobalEnv.p
        set_last_value(val)
    end
    val
end

"""
Evaluate an R symbol or language object (i.e. a function call) in an R
try/catch block, returning an RObject.
"""
reval(r::RObject, env=Const.GlobalEnv) = RObject(reval_p(sexp(r), sexp(env)))
reval(str::T, env=Const.GlobalEnv) where T <: Union{AbstractString, Symbol} =
    RObject(reval_p(rparse_p(str), sexp(env)))

ext: jl


text:
# ## BioModel7d

# ## Model

module biomodel7d #jl

using ReachabilityAnalysis

@taylorize function biomodel7d!(dx, x, p, t)
    dx[1] = -0.4 * x[1] + 5.0 * x[3] * x[4]
    dx[2] = 0.4 * x[1] - x[2]
    dx[3] = x[2] - 5.0 * x[3] * x[4]
    dx[4] = 5.0 * x[5] * x[6] - 5.0 * x[3] * x[4]
    dx[5] = -5.0 * x[5] * x[6] + 5.0 * x[3] * x[4]
    dx[6] = 0.5 * x[7] - 5.0 * x[5] * x[6]
    dx[7] = -0.5 * x[7] + 5.0 * x[5] * x[6]
    return dx
end

function model(X0)
    S = @system(x' = biomodel7d!(x), dim:7)
    return IVP(S, X0)
end

end #jl

# ## References

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using MLOPF
using Test

@testset "MLOPF.jl" begin
    # Write your tests here.
end


text:
import Base.convert

# strings
# convert(::Type{String}, a::Array{Uint8}) = bytestring(a)

# ints
convert(::Type{Array{Uint8}}, n::Int) = error("not implemented")
convert(::Type{Int}, a::Type{Array{Uint8}}) = error("not implemented")


text:
using Test
using Quiqbox
using Quiqbox.Molden

@testset "Molden.jl test" begin

mols = [
        ["H", "H"],
        ["H", "F"],
        ["N", "H", "H", "H"]
        ]
molNames = [
            "H2",
            "HF",
            "NH3"
            ]
br = 0.529177210903
# Data from CCCBDB: https://cccbdb.nist.gov
molCoords = [
             [[0.3705,0.0,0.0], [-0.3705,0.0,0.0]],
             [[0.4585,0.0,0.0], [-0.4585,0.0,0.0]],
             [[0.0, 0.0, 0.1111], [0.0, 0.9316, -0.2592],
              [0.8068, -0.4658, -0.2592], [-0.8068, -0.4658, -0.2592]]
             ] ./ br

bfCoords = [molCoords..., GridBox(1, 1.2) |> gridCoords]
bfs = ["STO-3G", "STO-2G"]
bsNames = push!(("-" .*molNames), "-Grid")
HFtypes = [:RHF, :UHF]
dir = @__DIR__
prefix1 = dir*"/"
prefix2 = dir*"/Moldens/"
for (nuc, nucCoords, molName, iMol) in zip(mols, molCoords, molNames, 1:length(mols)),
    (bfCoord, bsName) in zip(bfCoords[iMol:end], bsNames[iMol:end]),
    HFtype in HFtypes,
    bf in bfs

    flag = (bfCoord == nucCoords)
    if flag
        nucConfig = [(bf, i) for i in nuc]
        bs = genBasisFunc.(bfCoord, nucConfig) |> flatten
    else
        bs = genBasisFunc.(bfCoord, bf) |> flatten
        bsName = "-Float"*bsName
    end

    # Number of spin-orbitals must not be smaller than numbers of electrons.
    fVars = try runHF(bs, nuc, nucCoords; HFtype, printInfo=false) catch; continue end

    mol = Molecule(bs, nuc, nucCoords, fVars)
    fn = "Test_"*molName*"_"*bf*bsName*"_"*string(HFtype)
    fd = makeMoldenFile(mol; roundDigits=5, recordUMO=true, fileName=prefix1*fn)
    str1, str2 = replace.(read.((fd, prefix2*fn*".molden"), String),
                          r"[0-9]+\.[0-9]{5}(?![0-9])"=>"X.XXXXX")
    str1, str2 = replace.((str1, str2), "-X.XXXXX"=>" X.XXXXX")
    @test str1 == str2
    rm(fd)
end

end


text:
using DataStructures, DataArrays , DataFrames, StatsFuns, GLM , Distributions, MixedModels, StatsBase, JSON, StatLib, NamedArrays
using DataStructures, DataArrays , DataFrames, StatsFuns, GLM , Distributions, MixedModels, StatsBase, JSON, StatLib, NamedArrays

# ************ CUSTOMIZE *********************
#root="/mnt/resource/analytics/models/Jenny-o"
#root="/mnt/resource/analytics/models/NatChoice"
#root="/mnt/resource/analytics/models/rev"
#root="/mnt/resource/analytics/models/ALL#10"
#root="/mnt/resource/analytics/models/ALL#11"
#root="/mnt/resource/analytics/models/CDW#6"
#root="/mnt/resource/analytics/models/HormelChili#8"
#root="/mnt/resource/analytics/models/Rev#8"
#root="/mnt/resource/analytics/models/Jenny-o"
#root="/mnt/resource/analytics/models/hunts"
#root="/mapr/mapr04p/analytics0001/analytic_users/jmdl/rev"
#root="/mapr/mapr04p/analytics0001/analytic_users/jmdl/NatChoice"
#root="/mapr/mapr04p/analytics0001/analytic_users/jmdl/ALL10"

# ************ END CUSTOMIZE *********************
root = !isdefined(:root) ? pwd() : root

modelsDict = read_modelsDict(root)

function runGlm(root::String, modelname::Symbol, ranef::Symbol=:empty)
    println("Loding data : ")
    dfd = read_dfd(root)
    modelsDict = read_modelsDict(root)
    poolit!(dfd,modelsDict[:factors])
    m=modelsDict[modelname]
    return runGlm(dfd, m, ranef)
end

function runGlm(dfd::DataFrame,m::Dict, ranef::Symbol=:empty)
    if ranef==:empty
        f=genF(m[:y_var],m[:finalvars])
        cols = convert(Array{Symbol},vcat(m[:finalvars], [m[:y_var], m[:logvar]]))
        println("Ok - Running ",m[:modelName]," -  Glm on proc ",myid()," with : ",f)
    else
        f=genF(m[:y_var],setdiff(vcat(m[:finalvars],ranef),[:group]))
        dfd[ranef] = relevel(dfd[ranef], "none")
        cols = convert(Array{Symbol},vcat(m[:finalvars], [m[:y_var], m[:logvar]], ranef ))
        println("Ok - Running ",m[:modelName]," - ",ranef," Glm on proc ",myid()," with : ",f)
    end
    if m[:Buyer_Pos_P1_is1]
        return glm(f, dfd[(dfd[:buyer_pos_p1].==1),cols] , m[:dist], m[:lnk])
    else
        return return glm(f, dfd[cols] , m[:dist], m[:lnk])
    end
end

function runGlmm(root::String, modelname::Symbol, ranef::Symbol=:empty, runθ::Bool=true)
    println("Loding data : Glmm : ",modelname," ~~~ ",ranef)
    dfd = read_dfd(root)   #dfd = readtable(mod_fname,header=true);
    modelsDict = read_modelsDict(root) #modelsDict = readModels(jmod_fname)
    poolit!(dfd,modelsDict[:factors])
    m=modelsDict[modelname]
    if ranef==:empty
        v_out=Dict()
        for r in m[:raneff]
            v_out[r] = runGlmm(dfd, m, [r], runθ)
        end
        return v_out
    else
        return runGlmm(dfd, m, [ranef], runθ)
    end
end

function runGlmm(dfd::DataFrame,m::Dict, ranef::Array{Symbol}=[:empty], runθ::Bool=true)
    ranef = ranef==[:empty] ? m[:raneff]  : ranef
    for r in ranef
        dfd[r] = relevel(dfd[r], "none")
    end
    f=genF(m[:y_var],setdiff(m[:finalvars],[:group]),ranef)
    cols = convert(Array{Symbol},vcat(m[:finalvars], [m[:y_var], m[:logvar]], ranef ))
    println("Ok - Running ",m[:modelName]," - ",ranef," Glmm on proc ",myid()," with : ",f)
    if m[:Buyer_Pos_P1_is1]
        #gmm1 = fit!(glmm(f, dfd[(dfd[:buyer_pos_p1].==1),cols] , m[:dist], m[:lnk]), false, 0)   #gmm1 = fit!(glmm(f, dfd[(dfd[:buyer_pos_p1].==1),cols] , m[:dist], m[:lnk]), nAGQ=0)
        gmm1 = fit!(glmm(f, dfd[(dfd[:buyer_pos_p1].==1),cols] , m[:dist], m[:lnk]), verbose=false, fast=true)
        if runθ if sum(raneffect(gmm1)[:coef]) == 0 profileθ(gmm1) end  end
        return gmm1
    else
        #gmm1 = fit!(glmm(f, dfd[cols], m[:dist], m[:lnk]), false, 0)   #gmm1 = fit!(glmm(f, dfd[cols], m[:dist], m[:lnk]), nAGQ=0)
        gmm1 = fit!(glmm(f, dfd[cols] , m[:dist], m[:lnk]), verbose=false, fast=true)
        if runθ if sum(raneffect(gmm1)[:coef]) == 0 profileθ(gmm1) end end
        return gmm1
    end
end

# *******************************************
function runModels(root::String, modelsDict::Dict)
    dfd = read_dfd(root) #dfd = readtable(mod_fname,header=true);
    poolit!(dfd,modelsDict[:factors])
    res = OrderedDict()
    q = Symbol[]
    cnt=0
    for (k,v) in OrderedDict(m=>modelsDict[m][:raneff] for m in [:ipen,:idolocc,:iocc])
         for r in v push!(q,k); push!(q,r); cnt=cnt+1 end
    end
    x=reshape(q, (2,cnt))
    ml = permutedims(x, [2, 1])
    m=:empty
    for i in 1:length(ml[:,1])
        if m != ml[i,:][1]
            gr=Symbol("glm_"*string(ml[i,:][1]))
            res[gr] = runGlm(dfd,modelsDict[ml[i,:][1]]) #res[gr] = runGlm(mod_fname, jmod_fname, ml[i,:][1])
            m=ml[i,:][1]
        end
        r=Symbol("glmm_"*string(ml[i,:][1])*"_"*string(ml[i,:][2]))
        res[r] = runGlmm(dfd,modelsDict[m], [ml[i,:][2]], true)   #res[r] = runGlmm(mod_fname, jmod_fname, ml[i,:][1], ml[i,:][2])
        covg=Symbol("covglm_"*string(ml[i,:][1])*"_"*string(ml[i,:][2]))
        res[covg] = runGlm(dfd,modelsDict[m], ml[i,:][2])   #res[covg] = runGlm(mod_fname, jmod_fname, ml[i,:][1], ml[i,:][2])
    end
    return res
end
#dx = runModels(root,modelsDict)
dx = runModels(root,modelsDict)

"""
function consolidateResultsOLD(modelsDict::Dict,dx::OrderedDict)
  sdf = DataFrame(parameter=String[], coef=Int64[], stderr=Int64[], zval=Int64[], pval=Int64[], model=Symbol[], ranef=Symbol[], modelType=Symbol[])
    for m in [:iocc, :idolocc, :ipen]
        g = dx[Symbol("glm_"*string(m))]
        dfx = coefDF(g,false)
        dfx[:model] = Symbol(replace(string(m),"i",""))
        dfx[:ranef] = :none
        dfx[:modelType] = :GLM
        sdf = vcat(sdf,dfx)
        grp = sdf[(sdf[:modelType].==:GLM)&(sdf[:model].==Symbol(replace(string(m),"i","")))&(sdf[:parameter].=="group"),:coef][1]
        err = sdf[(sdf[:modelType].==:GLM)&(sdf[:model].==Symbol(replace(string(m),"i","")))&(sdf[:parameter].=="group"),:stderr][1]
        for r in  modelsDict[m][:raneff]
            k="_"*string(m)*"_"*string(r)
            g = dx[Symbol("glmm"*k)]
            dfx = raneffect(g)
            dfx[:zval]=0.0
            dfx[:pval]=0.0
            dfx[:model] = Symbol(replace(string(m),"i",""))
            dfx[:ranef]=r
            dfx[:modelType]=:GLMM
            sdf = vcat(sdf,dfx)
            g = dx[Symbol("covglm"*k)]
            dfx = coefDF(g,false)
            dfx[:model] = Symbol(replace(string(m),"i",""))
            dfx[:ranef] = r
            dfx[:modelType] = :RanCov
            dfx = dfx[findin(dfx[:parameter],filter(x->contains(x,string(r)) ,Array(dfx[:parameter]))),:]
            #dfx[:coef]=dfx[:coef]-grp
            #dfx[:stderr]=sqrt(dfx[:stderr].^2+err^2)
            dfx[:parameter] = map(x-> replace(x,string(r)*": ",""),dfx[:parameter])
            sdf = vcat(sdf,dfx)
        end
    end
    return sdf
end
dx = runModels(root,modelsDict)

function consolidateResults(modelsDict::Dict,dx::OrderedDict)
  sdf = DataFrame(parameter=String[], coef=Int64[], stderr=Int64[], zval=Int64[], pval=Int64[], model=Symbol[], ranef=String[], modelType=String[])
    for m in [:iocc, :idolocc, :ipen]
        g = dx[Symbol("glm_"*string(m))]
        dfx = coefDF(g,false)
        dfx[:model] = replace(string(m),"i","")
        dfx[:ranef] = "none"
        dfx[:modelType] = "GLM"
        sdf = vcat(sdf,dfx)
        grp = sdf[(sdf[:modelType].=="GLM")&(sdf[:model].==replace(string(m),"i",""))&(sdf[:parameter].=="group"),:coef][1]
        err = sdf[(sdf[:modelType].=="GLM")&(sdf[:model].==replace(string(m),"i",""))&(sdf[:parameter].=="group"),:stderr][1]
        for r in  modelsDict[m][:raneff]
            k="_"*string(m)*"_"*string(r)
            g = dx[Symbol("glmm"*k)]
            dfx = raneffect(g)
            dfx[:zval]=0.0
            dfx[:pval]=0.0
            dfx[:model] = replace(string(m),"i","")
            dfx[:ranef]=string(r)
            dfx[:modelType]="GLMM"
            sdf = vcat(sdf,dfx)
            g = dx[Symbol("covglm"*k)]
            dfx = coefDF(g,false)
            dfx[:model] = replace(string(m),"i","")
                dfx[:ranef] = string(r)
            dfx[:modelType] = "RanCov"
            dfx = dfx[findin(dfx[:parameter],filter(x->contains(x,string(r)) ,Array(dfx[:parameter]))),:]
            #dfx[:coef]=dfx[:coef]-grp
            #dfx[:stderr]=sqrt(dfx[:stderr].^2+err^2)
            dfx[:parameter] = map(x-> replace(x,string(r)*": ",""),dfx[:parameter])
            sdf = vcat(sdf,dfx)
        end
    end
    return sdf
end
"""

function consolidateResults(modelsDict::Dict,dx::OrderedDict=OrderedDict())
  sdf = DataFrame(parameter=String[], coef=Int64[], stderr=Int64[], zval=Int64[], pval=Int64[], model=String[], ranef=String[], modelType=String[])
    for m in [:iocc, :idolocc, :ipen]
        g = length(dx) > 0 ? dx[Symbol("glm_"*string(m))] : getw( Symbol("glm_"*string(m)) )
        #g = takew(  Symbol("glm_"*string(m))    )
        dfx = coefDF(g,false)
        dfx[:model] = replace(string(m),"i","")
        dfx[:ranef] = "none"
        dfx[:modelType] = "GLM"
        sdf = vcat(sdf,dfx)
        grp = sdf[(sdf[:modelType].=="GLM")&(sdf[:model].==replace(string(m),"i",""))&(sdf[:parameter].=="group"),:coef][1]
        err = sdf[(sdf[:modelType].=="GLM")&(sdf[:model].==replace(string(m),"i",""))&(sdf[:parameter].=="group"),:stderr][1]
        for r in  modelsDict[m][:raneff]
            k="_"*string(m)*"_"*string(r)
            g = length(dx) > 0 ?  dx[Symbol("glmm"*k)]  : getw(Symbol("glmm"*k))
            dfx = raneffect(g)
            dfx[:zval]=0.0
            dfx[:pval]=0.0
            dfx[:model] = replace(string(m),"i","")
            dfx[:ranef]=string(r)
            dfx[:modelType]="GLMM"
            sdf = vcat(sdf,dfx)
            g = length(dx) > 0 ?  dx[Symbol("covglm"*k)]    : getw(Symbol("covglm"*k))
            dfx = coefDF(g,false)
            dfx[:model] = replace(string(m),"i","")
                dfx[:ranef] = string(r)
            dfx[:modelType] = "RanCov"
            dfx = dfx[findin(dfx[:parameter],filter(x->contains(x,string(r)) ,Array(dfx[:parameter]))),:]
            #dfx[:coef]=dfx[:coef]-grp
            #dfx[:stderr]=sqrt(dfx[:stderr].^2+err^2)
            dfx[:parameter] = map(x-> replace(x,string(r)*": ",""),dfx[:parameter])
            sdf = vcat(sdf,dfx)
        end
    end
    return sdf
end
#dfx = consolidateResults(modelsDict)

#dx = runModels(root,modelsDict)
#dfx = consolidateResults(modelsDict, dx)
#save_dfx(root,dfx)   #writetable(root*"/campaign.csv", dfx)
dfx = consolidateResults(modelsDict, dx)
save_dfx(root,dfx)   #writetable(root*"/campaign.csv", dfx)

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
using .PyPlot
using .PyPlot.PyCall
import .PyPlot: loglog, plot, semilogx, semilogy, figure, fill_between

### overloaded 1D plotting

for plot in (:plot, :loglog, :semilogx, :semilogy)

	# Cℓs
    @eval function ($plot)(ic::InterpolatedCℓs, args...; kwargs...)
		($plot)(ic.ℓ, ic.Cℓ, args...; kwargs...)
	end
	@eval function ($plot)(ic::NamedTuple{<:Any,<:NTuple{<:Any,<:InterpolatedCℓs}}, args...; kwargs...)
		($plot).(values(ic), args...; kwargs...)
	end
	@eval function ($plot)(ic::InterpolatedCℓs{<:Measurement}, args...; kwargs...)
		errorbar(ic.ℓ, Measurements.value.(ic.Cℓ), Measurements.uncertainty.(ic.Cℓ), args...; marker=".", ls="", capsize=2, kwargs...)
		($plot) in [:loglog,:semilogx] && xscale("log")
		($plot) in [:loglog,:semilogy] && yscale("log")
	end

	# Loess-interpolated
	@eval function ($plot)(f::Function, m::Loess.LoessModel, args...; kwargs...)
	    l, = ($plot)(m.xs, f.(m.ys), ".", args...; kwargs...)
        xs′ = vcat(map(1:length(m.xs)-1) do i
		    collect(range(m.xs[i],m.xs[i+1],length=10))[1:end-1]
		end..., [last(m.xs)])
	    ($plot)(xs′, f.(m.(xs′)), args...; c=l.get_color(), kwargs...)
	end
	@eval ($plot)(m::Loess.LoessModel, args...; kwargs...) = ($plot)(identity, m, args...; kwargs...)

	# 1D KDE
	@eval function ($plot)(f::Function, k::GetDistKDE{1}, args...; kwargs...)
	    ($plot)(k.kde.x, f.(k.kde.P), args...; kwargs...)
	end
	@eval ($plot)(k::GetDistKDE{1}, args...; kwargs...) = ($plot)(identity, k, args...; kwargs...)

end

# 2D KDE
function plot(k::GetDistKDE{2}, args...; color=nothing, kwargs...)
	@unpack colors = pyimport("matplotlib")
	args = k.kde.x, k.kde.y, k.kde.P, [k.kde.getContourLevels([0.95,0.68]); Inf]
	if color == nothing
		color = gca()._get_lines.get_next_color()
	end
	contourf(args...; colors=[(colors.to_rgb(color)..., α) for α in (0.4, 0.8)], kwargs...)
	contour(args...;  colors=color, kwargs...)
end

# Cℓ band
function fill_between(ic::InterpolatedCℓs{<:Measurement}, args...; kwargs...)
	fill_between(
		ic.ℓ,
		((@. Measurements.value(ic.Cℓ) - x * Measurements.uncertainty(ic.Cℓ)) for x in (-1,1))...,
		args...; kwargs...
	)
end

### plotting FlatFields

plotsize₀ = 4

pretty_name(s) = pretty_name(Val.(Symbol.(split(string(s),"")))...)
pretty_name(::Val{s}, b::Val) where {s} = "$s "*pretty_name(b)
pretty_name(::Val{:x}) = "Map"
pretty_name(::Val{:l}) = "Fourier"

function _plot(f, ax, k, title, vlim, vscale, cmap; cbar=true, units=:deg, ticklabels=true, axeslabels=false, kwargs...)

	@unpack Nx, Ny, θpix = fieldinfo(f)
	ismap = endswith(string(k), "x")

	# default values
	if title == nothing
		if f isa FlatS0
			title = pretty_name(string(k)[2])
		else
			title = pretty_name(k)
		end
		title *= " ($(Ny)x$(Nx) @ $(θpix)')"
	end
	if vlim == nothing
		vlim = ismap ? :sym : :asym
	end
	if vscale == nothing
		vscale = ismap ? :linear : :log
	end
	if cmap == nothing
		if ismap
			cmap = get_cmap("RdBu_r")
		else
			cmap = get_cmap("viridis")
			cmap.set_bad("lightgray")
		end
	end

	# build array
	if ismap
		arr = Array(f[k])
	else
		arr = abs.(ifftshift(unfold(Array(f[k]), Ny)))
	end
	if vscale == :log
		arr[arr .== 0] .= NaN
	end

	# auto vlim's
	if vlim==:sym
        vmax = quantile(abs.(arr[@. !isnan(arr)][:]),0.999)
        vmin = -vmax
    elseif vlim==:asym
        vmin, vmax = (quantile(arr[@. !isnan(arr)][:],q) for q=(0.001,0.999))
    elseif isa(vlim,Tuple)
        vmin, vmax = vlim
    else
        vmax = vlim
        vmin = -vmax
    end

	# make the plot
	if ismap
		extent = [-Nx,Nx,-Ny,Ny] .* θpix/Dict(:deg=>60,:arcmin=>1)[units]/2
	else
		extent = [-1,1,-1,1] .* fieldinfo(f).nyquist
	end
	norm = vscale == :log ? matplotlib.colors.LogNorm() : nothing
	cax = ax.matshow(
		clamp.(arr, vmin, vmax);
		vmin=vmin, vmax=vmax, extent=extent,
		cmap=cmap, rasterized=true, norm=norm,
		kwargs...
	)

	# annonate
    if cbar
		colorbar(cax,ax=ax,pad=0.01)
	end
    ax.set_title(title, y=1)
    if ticklabels
		if ismap
			@pydef mutable struct MyFmt <: pyimport(:matplotlib).ticker.ScalarFormatter
				__call__(self,v,p=nothing) = py"super"(MyFmt,self).__call__(v,p)*Dict(:deg=>"°",:arcmin=>"′")[units]
			end
			ax.xaxis.set_major_formatter(MyFmt())
			ax.yaxis.set_major_formatter(MyFmt())
			if axeslabels
				ax.set_xlabel("RA")
				ax.set_ylabel("Dec")
			end
		else
	        ax.set_xlabel(raw"$\ell_x$")
	        ax.set_ylabel(raw"$\ell_y$")
	        ax.tick_params(axis="x", rotation=45)
		end
        ax.tick_params(labeltop=false, labelbottom=true)
    else
        ax.tick_params(labeltop=false, labelleft=false)
    end

end

@doc doc"""
    plot(f::Field; kwargs...)
    plot(fs::VecOrMat{\<:Field}; kwarg...)

Plotting fields.
"""
plot(f::Field; kwargs...) = plot([f]; kwargs...)
function plot(D::DiagOp; kwargs...)
	props = _sub_components[findfirst(((k,v),)->diag(D) isa @eval($k), _sub_components)][2]
	plot(
		[diag(D)];
		which = permutedims([Symbol(p) for (p,_) in props if endswith(p,r"[lx]")]),
		kwargs...
	)
end

function plot(
	fs :: AbstractVecOrMat{F};
	plotsize = plotsize₀,
	which = default_which(fs),
	title = nothing,
	vlim = nothing,
	vscale = nothing,
	cmap = nothing,
	return_all = false,
	kwargs...
) where {F<:FlatField}

    (m,n) = size(tuple.(fs, which)[:,:])
    fig,axs = subplots(m, n; figsize=plotsize.*[1.4*n,m], squeeze=false)
    axs = getindex.(Ref(axs), 1:m, (1:n)') # see https://github.com/JuliaPy/PyCall.jl/pull/487#issuecomment-456998345
    _plot.(fs,axs,which,title,vlim,vscale,cmap; kwargs...)

	if return_all
		(fig, axs, which)
	elseif isdefined(Main,:IJulia) && Main.IJulia.inited
		nothing # on IJulia, returning the figure can lead to it getting displayed twice
	else
		fig # returning the figure works on Juno/VSCode/Pluto
	end

end

default_which(::AbstractVecOrMat{<:FlatS0})  = [:Ix]
default_which(::AbstractVecOrMat{<:FlatS2})  = [:Ex :Bx]
default_which(::AbstractVecOrMat{<:FlatS02}) = [:Ix :Ex :Bx]
function default_which(fs::AbstractVecOrMat{<:FlatField})
    try
        ensuresame((default_which([f]) for f in fs)...)
    catch x
        x isa AssertionError ? throw(ArgumentError("Must specify `which` argument by hand for plotting this combination of fields.")) : rethrow()
    end
end

### animations of FlatFields

@doc doc"""
    animate(fields::Vector{\<:Vector{\<:Field}}; interval=50, motionblur=false, kwargs...)

"""
animate(f::AbstractVecOrMat{<:FlatField}; kwargs...) = animate([f]; kwargs...)
animate(annonate::Function, args...; kwargs...) = animate(args...; annonate=annonate, kwargs...)
function animate(fields::AbstractVecOrMat{<:AbstractVecOrMat{<:FlatField}}; fps=25, motionblur=false, annonate=nothing, filename=nothing, kwargs...)
    fig, axs, which = plot(first.(fields); return_all=true, kwargs...)
    motionblur = (motionblur == true) ? [0.1, 0.5, 1, 0.5, 0.1] : (motionblur == false) ? [1] : motionblur

    if (annonate!=nothing); annonate(fig,axs,which); end

    ani = pyimport("matplotlib.animation").FuncAnimation(fig,
        i->begin
            for (f,ax,k) in tuple.(fields,axs,which)
                if length(f)>1
                    img = ax.images[1]
                    img.set_data(sum(x*getindex(f[mod1(i-j+1,length(f))],k) for (j,x) in enumerate(motionblur)) / sum(motionblur))
                end
            end
            first.(getproperty.(axs,:images))[:]
        end,
        1:maximum(length.(fields)[:]),
        interval=1000/fps, blit=true
        )
    close()
    if filename!=nothing
        ani.save(filename,writer="imagemagick",savefig_kwargs=Dict(:facecolor=>fig.get_facecolor()))
        if endswith(filename,".gif")
            run(`convert -layers Optimize $filename $filename`)
        end
    end
    HTML(ani.to_html5_video())
end

### plotting HealpixFields

plot(f::HealpixMap; kwargs...) = hp.mollview(collect(f.arr); cmap="RdBu_r", kwargs...)

### convenience
# for plotting in environments that only show a plot if its the last thing returned

function figure(plotfn::Function, args...; kwargs...)
	figure(args...; kwargs...)
	plotfn()
	gcf()
end

### chains

"""
    plot_kde(samples; [boundary, normalize, smooth_scale_2D], kwargs...)

Plot a Kernel Density Estimate PDF for a set of 1D or 2D samples.
`boundary`, `normalize`, and `smooth_scale_2D` keyword arguments are
passed to to the underlying call to [`kde`](@ref), and all others are
passed to the underlying call to `plot`.

Based on Python [GetDist](https://getdist.readthedocs.io/en/latest/intro.html),
which must be installed.
"""
function plot_kde(samples; boundary=nothing, normalize=nothing, smooth_scale_2D=nothing, kwargs...)
	kde_kwargs = filter(!isnothing∘last, Dict(pairs((;boundary,normalize,smooth_scale_2D))))
    plot(kde(samples; kde_kwargs...); kwargs...)
end


text:
module MedEval3D

include(joinpath("utils","CUDAAtomicUtils.jl"))

include(joinpath( "utils","BitWiseUtils.jl"))
include(joinpath( "structs","BasicStructs.jl"))
include(joinpath( "utils","CUDAGpuUtils.jl"))
include(joinpath( "utils","IterationUtils.jl"))
include(joinpath( "utils","ReductionUtils.jl"))
include(joinpath( "utils","MemoryUtils.jl"))
include(joinpath( "distanceMetrics","Housdorff","verB","MetaDataUtils.jl"))
include(joinpath( "distanceMetrics","Housdorff","verB","ResultListUtils.jl"))

include(joinpath( "overLap","MainOverlap.jl"))
include(joinpath( "PairCounting","RandIndex.jl"))

include(joinpath( "Probabilistic","ProbabilisticMetrics.jl"))
include(joinpath( "InformationTheorhetic","InformationTheorhetic.jl"))
include(joinpath( "volume","VolumeMetric.jl"))

include(joinpath( "distanceMetrics","MeansMahalinobis.jl"))
include(joinpath( "distanceMetrics","Housdorff","verB","HFUtils.jl"))

include(joinpath( "distanceMetrics","Housdorff","verB","PrepareArrtoBool.jl"))
include(joinpath( "distanceMetrics","Housdorff","verB","WorkQueueUtils.jl"))
include(joinpath( "distanceMetrics","Housdorff","verB","ProcessMainDataVerB.jl"))
include(joinpath( "distanceMetrics","Housdorff","verB","ScanForDuplicates.jl"))
include(joinpath( "distanceMetrics","Housdorff","verB","MetadataAnalyzePass.jl"))
include(joinpath( "distanceMetrics","Housdorff","verB","MainLoopKernel.jl"))
include(joinpath( "distanceMetrics","Housdorff","verB","Housdorff.jl"))

include(joinpath( "kernels","TpfpfnKernel.jl"))
include(joinpath( "kernels","InterClassCorrKernel.jl"))
include(joinpath( "MainAbstractions.jl"))

end#MedEval3D


text:
module Pseudospectral

using FastGaussQuadrature, ApproxFun, Parameters, JuMP, SymEngine, Ipopt

include("diff_matrix.jl")
include("problem.jl")

export diff_matrix!, Bounds, UninitializedProblem, PSProblem, build_problem

end # module


text:
@testset "Composite helper functions" begin
    x1 = origin = (0,0)
    x2 = (1,1)
    x3 = (0,4)
    x4 = (0,2)
    binlimits = (0,2)
    gridspacing = 1
    @test AvailablePotentialEnergyFramework.distance(x1,x2,gridspacing) ==  AvailablePotentialEnergyFramework.distance(x2,x1,gridspacing) == sqrt(2)
    @test AvailablePotentialEnergyFramework.isindexindistancebin(binlimits,x2,x1) == true
    @test AvailablePotentialEnergyFramework.isindexindistancebin(binlimits,x3,x1) == false
    @test AvailablePotentialEnergyFramework.isindexindistancebin(binlimits,x4,x1) == true

    array_allindistance_2d = let
        array_allindistance_2d = zeros(100,100)
        for i in 1:100, j in 1:100
            array_allindistance_2d[i,j] = hypot(i,j)
        end
        array_allindistance_2d
    end
    @test isapprox(AvailablePotentialEnergyFramework.averageallindistance((99,100),array_allindistance_2d,origin), 100; rtol = 1)

    array_allindistance_3d = let
        array_allindistance_3d = zeros(100,100,10)
        for i in 1:100, j in 1:100
            array_allindistance_3d[i,j,:] .= hypot(i,j)
        end
        array_allindistance_3d
    end

    #@test isapprox.(AvailablePotentialEnergyFramework.averageallindistance((99,100),array_allindistance_3d,origin), 100; rtol = 1) == fill(true,10)

    @test all(isapprox.(AvailablePotentialEnergyFramework.velocity_topolar(10,0,x4,origin),(0,10),atol=0.1))

    @test all(isapprox.(AvailablePotentialEnergyFramework.velocity_topolar(-10,0,x4,origin),(0,-10),atol=0.1))

    @test all(isapprox.(AvailablePotentialEnergyFramework.velocity_topolar(0,10,x4,origin),(10,0),atol=0.1))

    @test all(isapprox.(AvailablePotentialEnergyFramework.velocity_topolar(0,-10,x4,origin),(-10,0),atol=0.1))

end

@testset "Cyclone detection tools" begin

    psfc = Dataset(joinpath(@__DIR__,"testfiles/test_composite_PSFC_TABS.nc")) do ds
        variable(ds, "PSFC")[:,:,:]
    end
    psfc = cat(psfc,fill(mean(psfc) ,size(psfc)) + 5rand(size(psfc)...),dims=3)

    TABS = Dataset(joinpath(@__DIR__,"testfiles/test_composite_PSFC_TABS.nc")) do ds
        variable(ds, "TABS")[:,:,:,:]
    end
    TABS = cat(TABS,fill(mean(TABS) ,size(TABS)) + 1rand(size(TABS)...),dims=3)
    @testset "Detect centers" begin
        pressure_anomaly = psfc .- mean(psfc,dims=(1,2))
        centerstest =  AvailablePotentialEnergyFramework.findcyclonecenters(AvailablePotentialEnergyFramework.PressureMinima(),pressure_anomaly[:,:,1],-5)

        @test length(centerstest) == 5

        centers_and_labels,cyclones = detect_cyclones(AvailablePotentialEnergyFramework.PressureMinima(),pressure_anomaly[:,:,1],-5,2000)

        @test length(centers_and_labels) == 6
        @test all([cyclones.segment_pixel_count[key] > 1 for key in keys(cyclones.segment_pixel_count)])

        (cyclonecount_2d,addition_2d) = AvailablePotentialEnergyFramework.add_allcyclones(psfc[:,:,1],cyclones,centers_and_labels; maskcyclones = false)
        @test cyclonecount_2d == 3

        (cyclonecount_2d,addition_2d) = AvailablePotentialEnergyFramework.add_allcyclones(psfc[:,:,1],cyclones,centers_and_labels; maskcyclones = true)
        @test cyclonecount_2d == 3

        (cyclonecount_3d,addition_3d) = AvailablePotentialEnergyFramework.add_allcyclones(TABS[:,:,:,1],cyclones,centers_and_labels; maskcyclones = false)
        @test cyclonecount_3d == 3

        (cyclonecount_3d,addition_3d) = AvailablePotentialEnergyFramework.add_allcyclones(TABS[:,:,:,1],cyclones,centers_and_labels; maskcyclones = true)
        @test cyclonecount_3d == 3

        binlimits = 0:2000:300000
        bins = [(binlimits[ind] , binlimits[ind + 1]) for ind in 1:(length(binlimits) - 1) ]

        @test_nowarn [averageallindistance(bin,addition_2d,(128,128),4000) for bin in bins]

        @test_nowarn [averageallindistance(bin,addition_3d,(128,128),4000) for bin in bins]
### Try and add two frames, one with and one without TC
        centers_labels_and_cyclones = [detect_cyclones(AvailablePotentialEnergyFramework.PressureMinima(),pressure_anomaly[:,:,i],-5,2000) for i in 1:2]
        @test 3 == begin
            totalcyclonecount = 0
            for timeindex in 1:2
                centers_and_labels,cyclones = centers_labels_and_cyclones[timeindex]
                if !isnothing(centers_and_labels)
                    count, _ = AvailablePotentialEnergyFramework.add_allcyclones(TABS[:,:,timeindex],cyclones,centers_and_labels;maskcyclones = false)
                    totalcyclonecount += count
                end
            end
            totalcyclonecount
        end
    end
end


text:
module Day8
export Node, parse_input, build_tree, part_1, part_2

struct Node
    metadata::Array{Int, 1}
    children::Array{Node}
end
Node(; metadata=[], children=[]) = Node(metadata, children)
Base.:(==)(x::Node, y::Node) = (x.metadata == y.metadata) & (x.children == y.children)

function parse_input(input::String)::Array{Int, 1}
    return reverse(map(x -> parse(Int, x), split(input)))
end

function build_tree(numbers::Array{Int, 1})::Node

    node = Node()
    n_children = pop!(numbers)
    n_metadata = pop!(numbers)

    for i in 1:n_children
        push!(node.children, build_tree(numbers))
    end

    for i in 1:n_metadata
        push!(node.metadata, pop!(numbers))
    end

    return node

end

function part_1(tree::Node)::Int
    total = sum(tree.metadata)
    for child in tree.children
        total += part_1(child)
    end
    return total
end

function part_2(tree::Node)::Int
    if length(tree.children) == 0
        return sum(tree.metadata)
    end

    value = 0
    for child_index in tree.metadata
        if child_index > length(tree.children)
            continue
        end
        value += part_2(tree.children[child_index])
    end

    return value
end

end


text:
using ParticleInCell

nstep = 250
time = zeros(nstep + 1)
energy = zeros(nstep + 1)
dimx, dimy = 2π, 1.0
nx, ny = 128, 16
nbpart = 100 * nx * ny
dx, dy = dimx / nx, dimy / ny

# particles = landau( dimx, nbpart)

@time vm2d2v(1, time, energy)

nbpart = 10
particles = zeros(7, nbpart)
fields = ones(3, nx, ny)
interpolation!(particles, fields, nbpart, nx, ny, dx, dy)


text:
using DFTK
import DFTK: mpi_mean!
using Test
using Random
using MPI
include("testcases.jl")

@testset "Forces on semiconductor (using total energy)" begin
    function energy_forces(pos)
        Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.psp))
        atoms = [Si => pos]
        model = model_DFT(silicon.lattice, atoms, :lda_xc_teter93)
        basis = PlaneWaveBasis(model; Ecut=5, kgrid=[2, 2, 2], kshift=[0, 0, 0])

        is_converged = DFTK.ScfConvergenceDensity(1e-10)
        scfres = self_consistent_field(basis; is_converged=is_converged)
        scfres.energies.total, compute_forces(scfres), compute_forces_cart(scfres)
    end

    # symmetrical positions, forces should be 0
    pos0  = [(ones(3)) / 8, -ones(3) / 8]
    _, F0 = energy_forces(pos0)
    @test norm(F0) < 1e-4

    pos1 = [([1.01, 1.02, 1.03]) / 8, -ones(3) / 8]  # displace a bit from equilibrium
    disp = rand(3)
    mpi_mean!(disp, MPI.COMM_WORLD)  # must be identical on all processes
    ε = 1e-7
    pos2 = [pos1[1] + ε * disp, pos1[2]]

    E1, F1, Fc1 = energy_forces(pos1)
    E2,  _,  _  = energy_forces(pos2)

    diff_findiff = -(E2 - E1) / ε
    diff_forces = dot(F1[1][1], disp)
    @test abs(diff_findiff - diff_forces) < 1e-4

    # Rough test against QE reference (using PZ functional)
    reference = [[[-0.005622025, -0.00445816, -0.003278985],
                  [ 0.005622025,  0.00445816,  0.003278985]]]
    @test maximum(v -> maximum(abs, v), reference[1] - Fc1[1]) < 1e-5
end

@testset "Forces on metal (using free energy)" begin
    function silicon_energy_forces(pos; smearing=Smearing.FermiDirac())
        Si = ElementPsp(silicon.atnum, psp=load_psp(silicon.psp))
        atoms = [Si => pos]
        model = model_DFT(silicon.lattice, atoms, :lda_xc_teter93;
                          temperature=0.03, smearing=smearing, spin_polarization=:collinear)
        basis = PlaneWaveBasis(model; Ecut=4, kgrid=[4, 1, 2], kshift=[1/2, 0, 0])

        n_bands = 10
        is_converged = DFTK.ScfConvergenceDensity(5e-10)
        scfres = self_consistent_field(basis, n_bands=n_bands,
                                       is_converged=is_converged)
        scfres.energies.total, compute_forces(scfres)
    end

    pos1 = [([1.01, 1.02, 1.03]) / 8, -ones(3) / 8] # displace a bit from equilibrium
    disp = rand(3)
    mpi_mean!(disp, MPI.COMM_WORLD)  # must be identical on all processes
    ε = 1e-6
    pos2 = [pos1[1] + ε * disp, pos1[2]]

    for (tol, smearing) in [(0.003, Smearing.FermiDirac), (5e-5, Smearing.Gaussian)]
        E1, F1 = silicon_energy_forces(pos1; smearing=smearing())
        E2, _ = silicon_energy_forces(pos2; smearing=smearing())

        diff_findiff = -(E2 - E1) / ε
        diff_forces = dot(F1[1][1], disp)

        @test abs(diff_findiff - diff_forces) < tol
    end
end

@testset "Forces on spin-polarised case" begin
    function oxygen_energy_forces(pos)
        O = ElementPsp(o2molecule.atnum, psp=load_psp("hgh/pbe/O-q6.hgh"))
        magnetic_moments = [O => [1.0, 1.0]]
        model = model_PBE(diagm([7.0, 7.0, 7.0]), [O => pos],
                          temperature=0.02, smearing=smearing=Smearing.Gaussian(),
                          magnetic_moments=magnetic_moments)
        basis = PlaneWaveBasis(model; Ecut=4, kgrid=[1, 1, 1])

        scfres = self_consistent_field(basis;
                                       is_converged=DFTK.ScfConvergenceDensity(1e-7),
                                       ρ=guess_density(basis, magnetic_moments),
                                       damping=0.7, mixing=SimpleMixing())
        scfres.energies.total, compute_forces(scfres)
    end

    pos1 = [[0, 0, 0.1155], [0.01, -2e-3, -0.2]]
    disp = rand(3)
    mpi_mean!(disp, MPI.COMM_WORLD)  # must be identical on all processes
    ε = 1e-6
    pos2 = [pos1[1] + ε * disp, pos1[2]]

    E1, F1 = oxygen_energy_forces(pos1)
    E2, _  = oxygen_energy_forces(pos2)

    diff_findiff = -(E2 - E1) / ε
    diff_forces  = dot(F1[1][1], disp)

    @test abs(diff_findiff - diff_forces) < 1e-4
end


text:
"""
    cudnnActivationForward(x; mode, nanOpt, coef, alpha)
    cudnnActivationForward(x, d::cudnnActivationDescriptor; alpha)
    cudnnActivationForward!(y, x; mode, nanOpt, coef, alpha, beta)
    cudnnActivationForward!(y, x, d::cudnnActivationDescriptor; alpha, beta)

Return the result of the specified elementwise activation operation applied to `x`.
Optionally `y` holds the result and `d` specifies the operation. `y` should be similar to
`x` if specified. Keyword arguments `alpha=1, beta=0` can be used for scaling, i.e. `y .=
alpha * op.(x) .+ beta * y`.  The following keyword arguments specify the operation if `d` is
not given:

* `mode = CUDNN_ACTIVATION_RELU`: Options are `SIGMOID`, `RELU`, `TANH`, `CLIPPED_RELU`, `ELU`, `IDENTITY`
* `nanOpt = CUDNN_NOT_PROPAGATE_NAN`: NaN propagation policy, the other option is `CUDNN_PROPAGATE_NAN`
* `coef=1`: When the activation mode is set to `CUDNN_ACTIVATION_CLIPPED_RELU`, this input specifies the clipping threshold; and when the activation mode is set to `CUDNN_ACTIVATION_ELU`, this input specifies the `alpha` parameter.
"""
cudnnActivationForward, cudnnActivationForward!

# Public methods
cudnnActivationForward(x; o...)     = cudnnActivationForwardWithDefaults(x; o...)
cudnnActivationForward!(y, x; o...) = cudnnActivationForwardWithDefaults(x; y, o...)
cudnnActivationForward(x, d::cudnnActivationDescriptor; o...)     = cudnnActivationForwardWithDefaults(x; activationDesc=d, o...)
cudnnActivationForward!(y, x, d::cudnnActivationDescriptor; o...) = cudnnActivationForwardWithDefaults(x; y, activationDesc=d, o...)

# Private method
function cudnnActivationForwardWithDefaults(
    x;
    y = similar(x),
    mode::cudnnActivationMode_t = CUDNN_ACTIVATION_RELU,
    nanOpt::cudnnNanPropagation_t = CUDNN_NOT_PROPAGATE_NAN,
    coef::Real=1,
    activationDesc::cudnnActivationDescriptor = cudnnActivationDescriptor(mode, nanOpt, Cdouble(coef)),
    alpha::Real=1,
    beta::Real=0,
    xDesc::cudnnTensorDescriptor = cudnnTensorDescriptor(x),
    yDesc::cudnnTensorDescriptor = xDesc,
)
    T = eltype(x)
    alpha, beta = scalingParameter(T,alpha), scalingParameter(T,beta)
    cudnnActivationForwardAD(x; activationDesc, alpha, xDesc, beta, yDesc, y)
end

# AD method:
function cudnnActivationForwardAD(x; activationDesc, alpha, xDesc, beta, yDesc, y)
    cudnnActivationForward(handle(), activationDesc, alpha, xDesc, x, beta, yDesc, y)
    return y
end

# Deprecated:
function cudnnActivationForward(x::DenseCuArray{T,N}, y::DenseCuArray{T,N}; o...) where {T,N}
    @warn "`cudnnActivationForward(x,y)` is deprecated, please use one of the methods in `@doc cudnnActivationForward`." maxlog=1
    cudnnActivationForward!(y, x; o...)
end


text:
import DomainSets: boundary

"""
    abstract type AbstractSubGrid{T,N} <: AbstractGrid{T,N} end

A subgrid of an underlying grid.
"""
abstract type AbstractSubGrid{T,N} <: AbstractGrid{T,N} end

"""
    supergrid(g::AbstractSubGrid)

The underlying grid of the subgrid.
"""
supergrid(g::AbstractSubGrid) = g.supergrid

include("indexsubgrid.jl")
include("maskedsubgrid.jl")
include("product.jl")
include("boundary.jl")

subgrid(grid::AbstractGrid, domain::Domain) = MaskedGrid(grid, domain)

function subgrid(grid::AbstractEquispacedGrid, domain::AbstractInterval)
    a = infimum(domain)
    b = supremum(domain)
    h = step(grid)
    idx_a = convert(Int, ceil( (a-grid[1])/step(grid))+1 )
    idx_b = convert(Int, floor( (b-grid[1])/step(grid))+1 )
    idx_a = max(idx_a, 1)
    idx_b = min(idx_b, length(grid))
    IndexSubGrid(grid, idx_a:idx_b, domain)
end

function subgrid(grid::ScatteredGrid, domain::Domain)
    mask = in.(grid, Ref(domain))
    points = grid.points[mask]
    ScatteredGrid(points, domain)
end

subgrid(grid::ProductGrid, domain::ProductDomain) =
    ProductGrid(map(subgrid, elements(grid), elements(domain))...)

# supergrid(grid::ProductGrid) = grid
# supergrid(grid::ProductGrid{Tuple{Vararg{G}}}) where G<:Union{MaskedGrid,IndexSubGrid} = ProductGrid(map(supergrid, elements(grid))...)


text:
m = AnSchorfheide()
DSGE.default_settings!(m)

@testset "Check default settings" begin
    # I/O File locations
    @test get_setting(m, :saveroot) == normpath(joinpath(dirname(@__FILE__), "..","save"))
    @test get_setting(m, :dataroot) == normpath(joinpath(dirname(@__FILE__), "..","save","input_data"))

    # Data settings for released and conditional data. Default behavior is to set vintage
    # of data to today's date.
    @test get_setting(m, :data_vintage) == Dates.format(now(), DSGE.DSGE_DATE_FORMAT)
    @test get_setting(m, :data_id) == 3
    @test get_setting(m, :cond_vintage) == Dates.format(now(), DSGE.DSGE_DATE_FORMAT)
    @test get_setting(m, :cond_full_names) ==  [:obs_gdp, :obs_corepce, :obs_spread,
                                                            :obs_nominalrate, :obs_longrate]
    @test get_setting(m, :cond_semi_names) == [:obs_spread, :obs_nominalrate, :obs_longrate]
    @test get_setting(m, :use_population_forecast) == false
    @test get_setting(m, :population_mnemonic).value == :CNP16OV__FRED
    @test get_setting(m, :hpfilter_population) == true
    @test get_setting(m, :rate_expectations_source) == :ois

    # Dates
    @test get_setting(m, :date_presample_start) == DSGE.quartertodate("1959-Q3")
    @test get_setting(m, :date_mainsample_start) == DSGE.quartertodate("1960-Q1")
    @test get_setting(m, :date_zlb_start) == DSGE.quartertodate("2008-Q4")
    @test get_setting(m, :date_forecast_start) == Dates.lastdayofquarter(Dates.today())
    @test get_setting(m, :date_conditional_end) == Dates.lastdayofquarter(Dates.today())

    # Anticipated shocks
    @test get_setting(m, :n_mon_anticipated_shocks) == 0
    @test get_setting(m, :n_mon_anticipated_shocks_padding) == 20

    # General computation
    @test get_setting(m, :use_parallel_workers) == true

    # Optimization
    @test get_setting(m, :reoptimize) == true
    @test get_setting(m, :optimization_method) == :csminwel
    @test get_setting(m, :optimization_iterations) == 100
    @test get_setting(m, :optimization_step_size) == .01
    @test get_setting(m, :simulated_annealing_temperature) == Optim.log_temperature
    @test get_setting(m, :simulated_annealing_block_proportion) == .3
    @test get_setting(m, :optimization_ftol) == 1e-10
    @test get_setting(m, :optimization_xtol) == 1e-10
    @test get_setting(m, :optimization_gtol) == 1e-10
    @test get_setting(m, :combined_optimizer_max_cycles) == 4
    @test get_setting(m, :optimization_attempts) == 4

    # Estimation
    @test get_setting(m, :sampling_method) == :MH
    @test get_setting(m, :calculate_hessian) == true
    @test get_setting(m, :n_hessian_test_params) == typemax(Int)
    @test get_setting(m, :use_chand_recursion) == false

    # Metropolis-Hastings
    @test get_setting(m, :n_mh_simulations) == 5000
    @test get_setting(m, :n_mh_blocks) == 22
    @test get_setting(m, :n_mh_param_blocks) == 1
    @test get_setting(m, :mh_adaptive_accept) == false
    @test get_setting(m, :mh_target_accept) == .25
    @test get_setting(m, :mh_c) == 0.5
    @test get_setting(m, :mh_α) == 1.0
    @test get_setting(m, :n_mh_burn) == 2
    @test get_setting(m, :mh_thin) == 5
    @test get_setting(m, :mh_cc) == 0.09
    @test get_setting(m, :mh_cc0) == 0.01

    # Forecast
    @test get_setting(m, :forecast_block_size) == 5000
    @test get_setting(m, :forecast_start_block) == 1
    @test get_setting(m, :forecast_input_file_overrides) == Dict{Symbol,String}()
    @test get_setting(m, :forecast_jstep) == 5
    @test get_setting(m, :forecast_uncertainty_override).hasvalue == false
    @test get_setting(m, :forecast_smoother) == :durbin_koopman
    @test get_setting(m, :forecast_horizons) == 60
    @test get_setting(m, :forecast_tdist_shocks) == false
    @test get_setting(m, :forecast_tdist_df_val) == 15
    @test get_setting(m, :forecast_zlb_value) == 0.13 / 4
    @test get_setting(m, :shockdec_startdate).hasvalue == false
    @test get_setting(m, :shockdec_enddate).hasvalue == false
    @test get_setting(m, :impulse_response_horizons) == 40
    @test get_setting(m, :compute_shockdec_bands) == false

    # Sequential Monte Carlo
    @test get_setting(m, :n_particles) == 2000
    @test get_setting(m, :n_Φ) == 300
    @test get_setting(m, :λ) == 2.1
    @test get_setting(m, :n_smc_blocks) == 1
    @test get_setting(m, :step_size_smc) == .5
    @test get_setting(m, :n_mh_steps_smc) == 1
    @test get_setting(m, :target_accept) == 0.25
    @test get_setting(m, :resampler_smc) == :multinomial
    @test get_setting(m, :mixture_proportion) == 1.

    # Adaptive ϕ Schedule
    @test get_setting(m, :use_fixed_schedule) == true
    @test get_setting(m, :adaptive_tempering_target_smc) == 0.97
    @test get_setting(m, :resampling_threshold) == 0.5

    # Temporary setting to save different output files
    @test get_setting(m, :smc_iteration) == 1
    @test get_setting(m, :previous_data_vintage) == Dates.format(now(), DSGE.DSGE_DATE_FORMAT)

    # Alternative policy
    @test !haskey(DSGE.get_settings(m), :alternative_policy)
end

m = AnSchorfheide(testing = true)
@testset "Check default test settings" begin
    # I/O
    @test get_setting(m, :dataroot) == normpath(joinpath(dirname(@__FILE__), "..", "test", "reference", "input_data"))

    # General
    @test get_setting(m, :data_vintage) == "REF"
    @test get_setting(m, :cond_vintage) == "REF"
    @test get_setting(m, :use_parallel_workers) == false
    @test get_setting(m, :n_hessian_test_params) == 3

    # Metropolis-Hastings
    @test get_setting(m, :n_mh_simulations) == 100
    @test get_setting(m, :n_mh_blocks) == 1
    @test get_setting(m, :n_mh_burn) == 0
    @test get_setting(m, :mh_thin) == 1

    # Forecast
    @test get_setting(m, :date_forecast_start) == DSGE.quartertodate("2015-Q4")
    @test get_setting(m, :forecast_horizons) == 2
    @test get_setting(m, :forecast_jstep) == 1
    @test get_setting(m, :impulse_response_horizons) == 2
end


text:
# test the second GLS implementation

"""
  Compute GLS in an inefficent way and compare to the efficient implementation
  in the code.  This is not a test function, but is called by test functions.
"""
function getGLS_test(mesh, sbp, eqn, opts)
#    println("----- Entered get GLS_test -----")
  Tmsh = eltype(mesh.dxidx)
  Tsol = eltype(eqn.q)
  Tres = eltype(eqn.res)

  for i=1:1  # DEBUGGING
    alpha_x = eqn.params.alpha_x
    alpha_y = eqn.params.alpha_y
    # do only element 1
    Dxi = diagm(1./sbp.w)*sbp.Q[:, :, 1]
    Deta = diagm(1./sbp.w)*sbp.Q[:, :, 2]

    dxidx_true = mesh.dxidx[:, :, 1, i]*mesh.jac[1,i]
    jac = mesh.jac[:, i]
    Dx = Dxi*dxidx_true[1, 1] + Deta*dxidx_true[2, 1]
    Dy = Dxi*dxidx_true[1, 2] + Deta*dxidx_true[2, 2]
    q = [eqn.q[1,1,i], eqn.q[1,2,i], eqn.q[1,3,i]]
    tau = zeros(Tres, 3,3)
    tau[1,1] = AdvectionEquationMod.getTau(alpha_x, alpha_y, jac[1], mesh.min_node_dist)
    tau[2,2] = AdvectionEquationMod.getTau(alpha_x, alpha_y, jac[2], mesh.min_node_dist)
    tau[3,3] = AdvectionEquationMod.getTau(alpha_x, alpha_y, jac[3], mesh.min_node_dist)

    tmp = alpha_x*Dx*q + alpha_y*Dy*q
    tmp2 = tau*diagm(sbp.w)*diagm(1./jac)*tmp
    gls_test = -(alpha_x*Dx.'*tmp2 + alpha_y*Dy.'*tmp2)

    fill!(eqn.res, 0.0)
    AdvectionEquationMod.applyGLS2(mesh, sbp, eqn, opts, eqn.src_func)

    gls_code = reshape(eqn.res[:, :, i], 3)
    @test isapprox( gls_code, gls_test) atol=1e-14
  end

#    println("----- finished get GLS_test -----")

  return nothing

end  # end function

"""
  Test GLS against reference implemntation, also tau for some cases
  where it is known exactly.
"""
function test_GLS2_term()
  @testset "----- Testing GLS2 -----" begin

    fname = "input_vals_GLS2.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    eqn.params.alpha_x = 2.0
    eqn.params.alpha_y = 1.0
    dxidx_true = mesh.dxidx[:, :, 1, 2]*mesh.jac[1,2]
    println("dxidx_true = ", dxidx_true)
    tau_code = AdvectionEquationMod.getTau(2.0, 1.0, dxidx_true, 2)
    @test isapprox( tau_code, 1/sqrt(8)) atol=1e-14

    fill!(eqn.q, 1.0)
    fill!(eqn.res, 0.0)
    AdvectionEquationMod.applyGLS2(mesh, sbp, eqn, opts, eqn.src_func)
    vals_code = reshape(eqn.res[:, :, 1], 3)
    @test isapprox( vals_code, zeros(Float64, 3)) atol=1e-14

    # set q to something interesting
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        x = mesh.coords[1, j, i]
        y = mesh.coords[2, j, i]
        eqn.q[1, j, i] = sin(x)*cos(y)
      end
    end

    # set source term
    eqn.src_func = AdvectionEquationMod.SRCDict["SRC0"]
    getGLS_test(mesh, sbp, eqn, opts)

    eqn.q[1, :, 1] = [1., 0, 0]

    fill!(eqn.res, 0.0)
    AdvectionEquationMod.applyGLS2(mesh, sbp, eqn, opts, eqn.src_func)
    res_tmp = reshape(eqn.res[1, :, 1], 3)

    # now do complex step
    println("----- Doing Complex Step -----")
    fname = "input_vals_GLS2.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    arg_dict["run_type"] = 5  # complex step run
    arg_dict["jac_method"] = 2  # complex step run
    f = open("input_vals_GLS2c.jl", "w")
    println(f, arg_dict)
    close(f)
    fname = "input_vals_GLS2c.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    # set q to something interesting
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        x = mesh.coords[1, j, i]
        y = mesh.coords[2, j, i]
        eqn.q[1, j, i] = sin(x)*cos(y)
      end
    end

  end  # end facts block

  return nothing
end  # end function

#test_GLS2_term()
add_func1!(AdvectionTests, test_GLS2_term, [TAG_SHORTTEST])

"""
  Finite difference and complex step test for GLS term
"""
function test_GLS2_jac()
  # finite difference checks
  @testset "----- Performing GLS2 Finite Difference Checks -----" begin

    fname = "input_vals_GLS2.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    eqn.src_func = AdvectionEquationMod.SRCDict["SRC0"]
    eqn.params.alpha_x = 2.0
    eqn.params.alpha_y = 1.0
    # set q to something interesting
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        x = mesh.coords[1, j, i]
        y = mesh.coords[2, j, i]
        eqn.q[1, j, i] = sin(x)*cos(y)
      end
    end

    fill!(eqn.res, 0.0)

    jac_fd = zeros(Float64, 3, 3)

    eps_fd = 1e-7
    # do initial calculation
    fill!(eqn.res, 0.0)
    AdvectionEquationMod.applyGLS2(mesh, sbp, eqn, opts, eqn.src_func)
    res_0 = copy(reshape(eqn.res[1, :, 1], 3))
    for i=1:3
      eqn.q[1, i, 1] += eps_fd
      fill!(eqn.res, 0.0)
      getGLS_test(mesh, sbp, eqn, opts)
      res_i = reshape(eqn.res[1, :, 1], 3)
      for j=1:3
        jac_fd[j, i] = (res_i[j] - res_0[j])/eps_fd
      end

      # undo perturbation
      eqn.q[1, i, 1] -= eps_fd
    end

    eqn.q[1, :, 1] = [1., 0, 0]

    fill!(eqn.res, 0.0)
    AdvectionEquationMod.applyGLS2(mesh, sbp, eqn, opts, eqn.src_func)
    res_tmp = reshape(eqn.res[1, :, 1], 3)

    # now do complex step
    println("----- Doing Complex Step -----")
    fname = "input_vals_GLS2.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)
    arg_dict["run_type"] = 5  # complex step run
    arg_dict["jac_method"] = 2
    f = open("input_vals_GLS2c.jl", "w")
    println(f, arg_dict)
    close(f)
    fname = "input_vals_GLS2c.jl"
    mesh, sbp, eqn, opts = solvePDE(fname)

    # set q to something interesting
    for i=1:mesh.numEl
      for j=1:mesh.numNodesPerElement
        x = mesh.coords[1, j, i]
        y = mesh.coords[2, j, i]
        eqn.q[1, j, i] = sin(x)*cos(y)
      end
    end

    eqn.params.alpha_x = 2.0
    eqn.params.alpha_y = 1.0

    # now we have complex equation object
    jac_c = zeros(Float64, 3,3)
    eps_complex = complex(0, 1e-20)
    for i=1:3
      eqn.q[1, i, 1] += eps_complex
      fill!(eqn.res, 0.0)
      getGLS_test(mesh, sbp, eqn, opts)  # do checks on complex version
  #    AdvectionEquationMod.applyGLS2(mesh, sbp, eqn, opts)
      res_i = reshape(eqn.res[1, :, 1], 3)
      for j=1:3
        jac_c[j, i] = imag(res_i[j])/abs(eps_complex)
      end

      eqn.q[1, i, 1] -= eps_complex
    end

    @test isapprox( jac_c, jac_fd) atol=1e-6

  end  # end facts block

  return nothing
end  # end function

#test_GLS2_jac()
add_func1!(AdvectionTests, test_GLS2_jac, [TAG_SHORTTEST])


text:
module Helpme

export @helpme

features(q::String) = unique(flatten(quote;$(parse(q));end))
features(q)         = features(repr(q))
flatten(s::Symbol)  = {s}
function flatten(e::Expr)
	flat = {}
	push!(flat, string(e.head))
	for arg in e.args
		append!(flat, isa(arg, Expr)? flatten(arg) : [arg])
	end

	push!(flat, e.typ)
	return flat
end

type Example
	suggestion::Symbol
	efeatures::Array
	sfeatures::Array
	weights::Dict{Any,Int}
	Example(s::Symbol) = new(s, {}, {}, Dict())
	function Example(e::Exception, f::String, s::Symbol)
		weights = Dict{Any,Int}()
		efeatures = features(e)
		sfeatures = features(f)
		for feature in efeatures
			weights[feature] = 1
		end

		for feature in sfeatures
			weights[feature] = 1
		end

		new(s, efeatures, sfeatures, weights)
	end
end

################################################################################

const database = Example[]
const suggbase = Dict{Symbol, String}()
const testbase = (Exception, String, Symbol)[]
macro example(suggestion, ex)
	exstr = string(ex)
	quote
		try
			$(esc(ex))
			warn("An example did not raise an error: "*$suggestion)
		catch e
			push!(database, Example(e, $exstr, $suggestion))
		end
	end
end

macro testcase(suggestion, ex)
	exstr = string(ex)
	quote
		try
			$ex
			warn("A testcase did not raise an error: "*$suggestion)
		catch e
			push!(testbase, (e, $exstr, $suggestion))
		end
	end
end

include("suggestions.jl")
include("examples.jl")
include("testcases.jl")

################################################################################

function distance2(example, efeatures, sfeatures, erepr, s)
	dist = 0
	for feature in example.efeatures
		if !(feature in efeatures)
			dist += example.weights[feature]
		end
	end

	for feature in example.sfeatures
		if !(feature in sfeatures)
			dist += example.weights[feature]
		end
	end

	return dist
end

function search(e::Exception, s::String)
	examples = Example[]
	dists = Float64[]
	max_dist = 0
	efeatures = features(e)
	sfeatures = features(s)
	erepr = repr(e)
	for example in database
		dist = distance2(example, efeatures, sfeatures, erepr, s)
		if length(examples) < 3
			push!(examples, example)
			push!(dists, dist)
			max_dist = max(dist, max_dist)
		elseif dist < max_dist
			for i in 1:3
				if dist < dists[i]
					example, examples[i] = examples[i], example
					dist, dists[i] = dists[i], dist
				end

				max_dist = maximum(dists)
			end
		end
	end

	for i in [1, 2, 1]
		if dists[i+1] < dists[i]
			examples[i+1], examples[i] = examples[i], examples[i+1]
			dists[i+1], dists[i] = dists[i], dists[i+1]
		end
	end

	return examples
end

################################################################################

function reroll(example, adj)
	for (feature, weight) in example.weights
		if rand() < 0.5
			example.weights[feature] += adj
		end
	end
end

function train()
	info("Training database")
	clean = false
	while !clean
		clean = true
		for (error, fnstr, suggestion) in testbase
			results = search(error, fnstr)
			if !(suggestion in map(x->x.suggestion, results))
				clean = false
				for example in database
					if suggestion == example.suggestion
						reroll(example, -1)
					end
				end

				for result in results
					reroll(result, 1)
				end
			else
				# nothing
			end
		end
	end
end

function save()
	info("Saving database")
	open(joinpath(Pkg.dir(), "Helpme", "src", string("database.",VERSION,".jl")), "w") do f
		for (i, example) in enumerate(database)
			for (j, feature) in enumerate(example.efeatures)
				w = database[i].weights[database[i].efeatures[j]]
				println(f, "database[$i].weights[database[$i].efeatures[$j]]=$w")
			end

			for (j, feature) in enumerate(example.sfeatures)
				w = database[i].weights[database[i].sfeatures[j]]
				println(f, "database[$i].weights[database[$i].sfeatures[$j]]=$w")
			end
		end
	end
end

if isfile(joinpath(Pkg.dir(), "Helpme", "src", string("database.",VERSION,".jl")))
	include(string("database.",VERSION,".jl"))
else
	train()
	save()
end

################################################################################

macro helpme(ex)
	s = string(ex)
	quote
		try
			$(esc(ex))
		catch e
			results = search(e, $s)
			deletes = Int[]
			messages = Dict()
			for i in 1:length(results)
				if haskey(messages, results[i].suggestion)
					push!(deletes, i)
				else
					messages[results[i].suggestion] = true
				end
			end

			for del in reverse(deletes)
				deleteat!(results, del)
			end

			if length(results) < 3
				push!(results, Example(:DEFAULTMSG))
			end

			for example in results
				info(suggbase[example.suggestion])
			end

			rethrow(e)
		end
	end
end

end


text:
# Shallow ice approximation (SIA) implicit solver for Greenland (steady state)

using JLD, Plots, Printf, LinearAlgebra

# enable plotting & saving by default
if !@isdefined do_visu; do_visu = true end
if !@isdefined do_save; do_save = true end

# finite difference stencil operation support functions
@views av(A)    = 0.25*(A[1:end-1,1:end-1].+A[2:end,1:end-1].+A[1:end-1,2:end].+A[2:end,2:end]) # average
@views av_xa(A) = 0.5.*(A[1:end-1,:].+A[2:end,:]) # average x-dir
@views av_ya(A) = 0.5.*(A[:,1:end-1].+A[:,2:end]) # average y-dir
@views inn(A)   = A[2:end-1,2:end-1] # inner points

@views function iceflow(dx, dy, Zbed, Hice, Mask, grad_b, z_ELA, b_max)
    println("Initialising ice flow model ... ")
    # physics
    s2y      = 3600*24*365.25  # seconds to years
    rho_i    = 910.0           # ice density
    g        = 9.81            # gravity acceleration
    npow     = 3.0             # Glen's power law exponent
    a0       = 1.5e-24         # Glen's law enhancement term
    # numerics
    @assert (dx>0 && dy>0) "dx and dy need to be positive"
    nx, ny   = size(Zbed,1), size(Zbed,2) # numerical grid resolution
    @assert (nx, ny) == size(Zbed) == size(Hice) == size(Mask) "Sizes don't match"
    itMax    = 1e5             # number of iteration (max)
    nout     = 200             # error check frequency
    tolnl    = 1e-6            # nonlinear tolerance
    epsi     = 1e-4            # small number
    damp     = 0.85            # convergence accelerator (this is a tuning parameter, dependent on e.g. grid resolution)
    dtausc   = 1.0/3.0         # iterative dtau scaling
    # derived physics
    a        = 2.0*a0/(npow+2)*(rho_i*g)^npow*s2y
    # derived numerics
    cfl      = max(dx^2,dy^2)/4.1
    # array initialization
    Err      = zeros(nx  , ny  )
    dSdx     = zeros(nx-1, ny  )
    dSdy     = zeros(nx  , ny-1)
    gradS    = zeros(nx-1, ny-1)
    D        = zeros(nx-1, ny-1)
    qHx      = zeros(nx-1, ny-2)
    qHy      = zeros(nx-2, ny-1)
    dtau     = zeros(nx-2, ny-2)
    ResH     = zeros(nx-2, ny-2)
    dHdt     = zeros(nx-2, ny-2)
    Vx       = zeros(nx-1, ny-1)
    Vy       = zeros(nx-1, ny-1)
    M        = zeros(nx  , ny  )
    B        = zeros(nx  , ny  )
    H        = zeros(nx  , ny  )
    S        = zeros(nx  , ny  )
    # initial conditions
    B       .= Zbed
    H       .= Hice
    S       .= B .+ H
    # iteration loop
    println(" starting iteration loop:")
    iter = 1; err = 2*tolnl
    while err>tolnl && iter<itMax
        Err   .= H
        # mass balance
        M     .= min.(grad_b.*(S .- z_ELA), b_max)
        # compute diffusivity
        dSdx  .= diff(S, dims=1)/dx
        dSdy  .= diff(S, dims=2)/dy
        gradS .= sqrt.(av_ya(dSdx).^2 .+ av_xa(dSdy).^2)
        D     .= a*av(H).^(npow+2) .* gradS.^(npow-1)
        # compute flux
        qHx   .= .-av_ya(D).*diff(S[:,2:end-1], dims=1)/dx
        qHy   .= .-av_xa(D).*diff(S[2:end-1,:], dims=2)/dy
        # update ice thickness
        dtau  .= dtausc*min.(10.0, cfl./(epsi .+ av(D)))
        ResH  .= .-(diff(qHx, dims=1)/dx .+ diff(qHy, dims=2)/dy) .+ inn(M)
        dHdt  .= dHdt.*damp .+ ResH
        H[2:end-1,2:end-1] .= max.(0.0, inn(H) .+ dtau.*dHdt)
        # apply mask (a very poor-man's calving law)
        H[Mask.==0] .= 0.0
        # update surface
        S     .= B .+ H
        # error check
        if mod(iter, nout)==0
            Err .= Err .- H
            err = norm(Err)/length(Err)
            @printf(" iter = %d, error = %1.2e \n", iter, err)
            if isnan(err)
                error("""NaNs encountered.  Try a combination of:
                             decreasing `damp` and/or `dtausc`, more smoothing steps""")
            end
        end
        iter += 1
    end
    # compute velocities
    Vx .= -D./(av(H) .+ epsi).*av_ya(dSdx)
    Vy .= -D./(av(H) .+ epsi).*av_xa(dSdy)
    # return as GeoArrays
    return  as_geoarray(H,  Zbed, name=:thickness),
            as_geoarray(S,  Zbed, name=:surface),
            as_geoarray(M,  Zbed, name=:smb),
            as_geoarray(Vx, Zbed, name=:vel_x, staggerd=true),
            as_geoarray(Vy, Zbed, name=:vel_y, staggerd=true)
end
# ------------------------------------------------------------------------------
include(joinpath(@__DIR__, "helpers.jl"))

# load the data
print("Loading the data ... ")
Zbed, Hice, Mask, dx, dy, xc, yc = load_data(; nx=96) # nx=96,160 are included in the repo
                                                      # other numbers will trigger a 2GB download
println("done.")

# apply some smoothing
print("Applying some smoothing ... ")
for is=1:2 # two smoothing steps
    smooth!(Zbed)
    smooth!(Hice)
end
println("done.")

# calculate mass balance coefficients for given spatial grid
grad_b, z_ELA, b_max = mass_balance_constants(xc, yc)

# run the SIA flow model
H, S, M, Vx, Vy = iceflow(dx, dy, Zbed, Hice, Mask, grad_b, z_ELA, b_max)

# visualization and save
nx, ny = size(H)
if do_visu
    !ispath("../output") && mkdir("../output")

    H_v = copy(H); H_v[Mask.==0].=NaN
    Hice_v = copy(Hice); Hice_v[Mask.==0].=NaN
    S_v = copy(S); S_v[Mask.==0].=NaN
    M_v = copy(M); M_v[Mask.==0].=NaN
    V_v = sqrt.(Vx.^2 .+ Vy.^2)

    # outputs
    fontsize  = 7
    opts = (aspect_ratio=1, yaxis=font(fontsize, "Courier"), xaxis=font(fontsize, "Courier"),
            ticks=nothing, framestyle=:box, titlefontsize=fontsize, titlefont="Courier", colorbar_title="",
            xlabel="", ylabel="", xlims=(dims(H_v)[1][1],dims(H_v)[1][end]), ylims=(dims(H_v)[2][end],dims(H_v)[2][1]) )
    p1 = heatmap(S_v; c=:davos, title="Surface elev. [m]", opts...)
    p2 = heatmap(H_v; c=:davos, title="Ice thickness [m]", opts...)
    p3 = heatmap(log10.(V_v); clims=(0.1, 2.0), title="log10(vel) [m/yr]", opts...)
    p4 = heatmap(M_v; c=:devon, title="Mass Bal. rate [m/yr]", opts...)
    p = plot(p1, p2, p3, p4, size=(400,400), dpi=200) #background_color=:transparent, foreground_color=:white
    ## uncomment if you want a pop-up plot pane showing:
    # display(p)
    savefig("../output/iceflow_out1.png")

    # error
    H_diff = Hice_v.-H_v
    fontsize = 7
    p1 = heatmap(Hice_v; c=:davos, title="H data [m]", opts...)
    p2 = heatmap(H_v; c=:davos, title="H model [m]", opts...)
    clim = max(abs.(extrema(H_diff[.!isnan.(H_diff)]))...)
    p3 = heatmap(H_diff; title="H (data-model) [m]",  clims=(-clim,clim), seriescolor=:balance, opts...)
    p = plot(p1, p2, p3, layout=(1, 3), size=(500,160), dpi=200) #background_color=:transparent, foreground_color=:white
    ## uncomment if you want a pop-up plot pane showing:
    # display(p)
    savefig("../output/iceflow_out2.png")
end

if do_save
    save("../output/iceflow_$(nx)x$(ny).jld", "Hice", Hice,
                                              "Mask", Mask,
                                              "H"   , H,
                                              "S"   , S,
                                              "M"   , M,
                                              "Vx"  , Vx,
                                              "Vy"  , Vy,
                                              "xc", xc, "yc", yc)
end

println("... done.")


text:
# key words in forward flight
V #forward velocity
  #thus the rotor blade loading and motion are periodic with a fundamental
  #frequency equal to the rotor speed Ω
A = π*R^2 #the area of a circle with 1/2 diameter
          #also: Cdi = CL^2/(π*A*R)
vi #the induced velocity of the helicopter rotor in high-speed forward flight
vi = T/(2*ρ*A*V)
AR = R/c = (N/π)*σ #aspect ratio should be high
U = sqrt(uT^2+uP^2) #
Φ = atan(uP/UT)
α = θ - Φ #the aerodynamic angle-of-attack of the blade element
L = 1/2*ρ*U^2*c*cl
D = 1/2*ρ*U^2*c*cd
Fz = L*cos(Φ)-D*sin(Φ)
Fx = L*sin(Φ)+D*cos(Φ)
dT = N*Fz*dr
dQ = N*Fx*r*dr
dP = Ω*dQ = N*Fx*Ω*r*dr

max_stars_repo_path: src/Forwa Flight Wake.jl

max_stars_repo_name: dingzhw/RotorTrim

max_issues_repo_path: src/Forwa Flight Wake.jl

max_issues_repo_name: helicopterliu/RotorTrim

max_forks_repo_path: src/Forwa Flight Wake.jl

max_forks_repo_name: helicopterliu/RotorTrim


text:
@testset "lotka volterra" begin
    @system LotkaVolterra(Controller) begin
        t(context.clock.time) ~ track(u"yr")
        N(N, P, b, a): prey_population => b*N - a*N*P ~ accumulate(init=N0)
        P(N, P, c, a, m): predator_population => c*a*N*P - m*P ~ accumulate(init=P0)
        N0: prey_initial_population ~ preserve(parameter)
        P0: predator_initial_population ~ preserve(parameter)
        b: prey_birth_rate ~ preserve(u"yr^-1", parameter)
        a: predation_rate ~ preserve(u"yr^-1", parameter)
        c: predator_reproduction_rate ~ preserve(parameter)
        m: predator_mortality_rate ~ preserve(u"yr^-1", parameter)
    end
    config = @config (
        :Clock => (;
            step = 1u"d",
        ),
        :LotkaVolterra => (;
            b = 0.6,
            a = 0.02,
            c = 0.5,
            m = 0.5,
            N0 = 20,
            P0 = 30,
        ),
    )
    stop = 20u"yr"
    r = simulate(LotkaVolterra; config, stop)
    @test r.t[end] >= stop
    visualize(r, :t, [:N, :P], names=["prey", "predator"]) |> display
end


text:
"""
    SCHash

Separate chaining hash table, using dynamic arrays.
"""
mutable struct SCHash{K, V}
    table::Vector{Vector{Pair{K, V}}}
    nkeys::Int # number of items
    tablesize::Int # table size
    function SCHash{K, V}(m::Int = 10) where {K, V}
        table = Vector{Vector{Pair{K, V}}}(undef, m)
        for i in 1:m
            table[i] = []
        end
        new(table, 0, m)
    end
end

"""
    SCLHash

Separate chaining hash table, using linked list.
"""
mutable struct SCLHash{K, V}
    table::Vector{DLinkedList{Pair{K, V}}}
    nkeys::Int
    tablesize::Int
    function SCLHash{K, V}(m::Int = 10) where {K, V}
        table = Vector{DLinkedList{Pair{K, V}}}(undef, m)
        for i in 1:m
            table[i] = DLinkedList{Pair{K, V}}()
        end
        new(table, 0, m)
    end
end

"""
    hash(map, key)

Hash function for key, implement for allowed key types, for now Numbers and Strings.
"""
function hash(map::Union{SCHash, SCLHash}, key::K) where K <: Number
    return compress(Int(ceil(abs(key))), map.tablesize)
end

function hash(map::Union{SCHash, SCLHash}, key::K) where K <: AbstractString
    h = 0
    for c in key
        h  = 37h + UInt8(c)
    end
    return compress(h, map.tablesize)
end

"""
    compress(h, tablesize)

Compress the hash, quite simple.
"""
function compress(h::Int, tablesize::Int)
    return h % tablesize + 1
end

"""
    get(map, key)

Get value from key.
"""
function get(map::SCHash, key::K) where K
    if map.nkeys == 0 return nothing end

    bucket = map.table[hash(map, key)]
    for (i, v) in enumerate(bucket)
        if key == v[1]
            return v[2]
        end
    end

    return nothing
end

function get(map::SCLHash, key::K) where K
    if map.nkeys == 0 return nothing end

    bucket = map.table[hash(map, key)]
    n = bucket.head
    while n != n.next
        if key == n.value[1]
            return n.value[2]
        end
    end

    return nothing
end

"""
    put(map, key, value)

Add value to hash table with key.
"""
function put!(map::SCHash, key::K, value::V) where {K, V}
    bucket = map.table[hash(map, key)]
    for (i, v) in enumerate(bucket)
        if key == v[1]
            bucket[i] = Pair(key, value)
            return nothing
        end
    end
    push!(bucket, Pair(key, value))
    map.nkeys += 1

    #if map.nkeys > 8 * map.tablesize resize(map) end

    return nothing
end

function put!(map::SCLHash, key::K, value::V) where {K, V}
    bucket = map.table[hash(map, key)]
    n = bucket.head
    while n != n.next
        if key == n.value[1]
            n.value = Pair(key, value)
            return nothing
        end
        n = n.next
    end
    addlast!(bucket, Pair(key, value))
    map.nkeys += 1
    map.tablesize = length(map.table)

    return nothing
end

"""
    delete!(map, key)

Delete kay value pair from map.
"""
function delete!(map::SCHash, key::K)  where K
    bucket = map.table[hash(map, key)]
    for (i, v) in enumerate(bucket)
        if key == v[1]
            bucket[i][2] = value
            return nothing
        end
    end

    return nothing
end

function delete!(map::SCLHash, key::K)  where K
    h = hash(map, key)
    bucket = map.table[h]
    n = bucket.head
    while n != n.next
        if key == n.value[1]
            n.prev.next = n.next
            n.next.prev = n.prev
            return nothing
        end
    end

    return nothing
end

"""
    resize(map, factor)

Resize the map. Resizing should be to a close prime (not implemented).
"""
function resize(map::SCHash{K, V}, factor) where {K,V}
    oldtable = map.table
    map.table = Vector{Vector{Pair{K, V}}}(undef, factor * map.tablesize)
    for bucket in oldtable
        for i in bucket
            put!(map, i)
        end
    end
end

function resize(map::SCLHash{K, V}, factor) where {K,V}
    oldtable = map.table
    map.table = Vector{DLinkedList{Pair{K, V}}}(undef, factor * map.tablesize)
    for bucket in oldtable
        n = bucket.head
        while n != n.next
            put!(map, n.value)
            n = n.next
        end
    end
end

"""
    OAHash

Open addressing hash table.
"""
struct OAHash
    table
end


text:
#=
  http://www.cs.ubc.ca/~murphyk/Bayes/bnintro.html

  Cf ~/problog/sprinkler.pl
     ~/blog/sprinkler.blog
     ~/psi/sprinkler.psi
     ~/webppl/sprinkler.wppl

=#

using Turing, StatsPlots, DataFrames
include("jl_utils.jl")

@model function sprinkler()
    cloudy ~ flip(0.5)

    sprinkler ~ cloudy ? flip(0.1) : flip(0.5)

    rain ~ cloudy ? flip(0.8) : flip(0.2)

    wet_grass ~ (sprinkler==false && rain==false) ? flip(0.0) :
                (sprinkler==true && rain==false) ? flip(0.9) :
                (sprinkler==false && rain==true) ? flip(0.9) :
                (sprinkler==true && rain==true) ? flip(0.99) : flip(0.0)

    true ~ Dirac(wet_grass == true)
    # true ~ Dirac (wet_grass == false);

    return cloudy, sprinkler, rain, wet_grass
end

model = sprinkler()

num_chains = 4

# chains = sample(model, Prior(), MCMCThreads(), 10_000, num_chains)

# chains = sample(model, MH(), 10_000)
# chains = sample(model, MH(), MCMCThreads(), 40_000, num_chains)

# chains = sample(model, PG(15), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, PG(20), 1_000)

chains = sample(model, SMC(1000), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, SMC(1000), 40_000)

# chains = sample(model, IS(), MCMCThreads(), 10_000, num_chains)
# chains = sample(model, IS(), 10_000)

# chains = sample(model, NUTS(1000,0.65), 1_000)
# chains = sample(model, HMC(0.1,5), 1_000)
# chains = sample(model, Gibbs(MH(:gender),NUTS(1000,0.65,:height)), 1_000)
# chains = sample(model, Gibbs(MH(:gender),NUTS(10,0.65,:height)), 1_000)
# chains = sample(model, Gibbs(MH(:gender),HMC(0.1,5,:height)), 1_000)
# chains = sample(model, Gibbs(PG(10,:gender),HMC(0.1,5,:height)), 1_000)
# chains = sample(model, Gibbs(MH(:gender),NUTS(1_000,0.65,:height)), 1_000)

display(chains)
# display(plot(chains))

gen = generated_quantities(model, chains)
show_var_dist_pct(gen, 120)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
# This file was generated, do not modify it. # hide
figure(figsize=(8,6))
plot(curves.parameter_values, curves.measurements)
ylabel("Root Mean Squared error", fontsize=16)
xlabel("Number of trees", fontsize=16)
xticks([10, 250, 500, 750, 1000], fontsize=14)
yticks(fontsize=14)

savefig(joinpath(@OUTPUT, "A-ensembles-2-curves.svg")) # hide


text:
module Kaleidoscope

using Unicode
import LLVM

include("lexer.jl")
include("ast.jl")
include("scope.jl")
include("codegen.jl")
include("run.jl")

end # module


text:
"""
ASGD Weight-Dropped LSTM
"""

using Flux
import Flux: gate, tanh, σ, Tracker, params, gpu, cpu, _testmode!, rand!, _dropout_kernel

cd(@__DIR__)
include("utils.jl")

gpu!(entity) = nothing
cpu!(entity) = nothing
reset_masks!(entity) = nothing
reset_probability!(entity) = nothing

# Generates Mask
function drop_mask(x, p)
    y = similar(x, size(x))
    rand!(y)
    y .= _dropout_kernel.(y, p, 1 - p)
    return y
end

drop_mask(shape::Tuple, p; type = Float32) = (mask = rand(type, shape...);mask .= _dropout_kernel.(mask, p, 1 - p))

#################### Weight-Dropped LSTM Cell ######################
mutable struct WeightDroppedLSTMCell{A, V, M}
    Wi::A
    Wh::A
    b::V
    h::V
    c::V
    p::Float64
    maskWi::M
    maskWh::M
    active::Bool
end

function WeightDroppedLSTMCell(in::Integer, out::Integer, probability::Float64=0.0;
    init = Flux.glorot_uniform)
    cell = WeightDroppedLSTMCell(
        param(init(out*4, in)),
        param(init(out*4, out)),
        param(init(out*4)),
        param(zeros(Float32, out)),
        param(zeros(Float32, out)),
        probability,
        drop_mask((out*4, in), probability),
        drop_mask((out*4, out), probability),
        true
    )
    cell.b.data[gate(out, 2)] .= 1
    return cell
end

function (m::WeightDroppedLSTMCell)((h, c), x)
    b, o = m.b, size(h, 1)
    Wi = m.active ? m.Wi .* m.maskWi : m.Wi
    Wh = m.active ? m.Wh .* m.maskWh : m.Wh
    g = Wi*x .+ Wh*h .+ b
    input = σ.(gate(g, o, 1))
    forget = σ.(gate(g, o, 2))
    cell = tanh.(gate(g, o, 3))
    output = σ.(gate(g, o, 4))
    c = forget .* c .+ input .* cell
    h′ = output .* tanh.(c)
    return (h′, c), h′
end

Flux.@treelike WeightDroppedLSTMCell

_testmode!(m::WeightDroppedLSTMCell, test) = (m.active = !test)

# Weight-Dropped LSTM [stateful]
function WeightDroppedLSTM(a...; kw...)
    cell = WeightDroppedLSTMCell(a...;kw...)
    hidden = (cell.h, cell.c)
    return Flux.Recur(cell, hidden, hidden)
end

function reset_masks!(wd::T) where T <: Flux.Recur{<:WeightDroppedLSTMCell}
    wd.cell.maskWi = drop_mask(wd.cell.Wi, wd.cell.p)
    wd.cell.maskWh = drop_mask(wd.cell.Wh, wd.cell.p)
    return
end
####################################################################

################## ASGD Weight-Dropped LSTM Layer ##################
mutable struct AWD_LSTM
    layer::Flux.Recur
    T::Integer
    accum
end

AWD_LSTM(in::Integer, out::Integer, probability::Float64=0.0; kw...) = AWD_LSTM(WeightDroppedLSTM(in, out, probability; kw...), -1, [])

Flux.@treelike AWD_LSTM

(m::AWD_LSTM)(in) = m.layer(in)

set_trigger!(t, m) = nothing
set_trigger!(trigger_point::Integer, m::AWD_LSTM) = m.T = trigger_point;

function gpu!(m::AWD_LSTM)
    m.layer = gpu(m.layer)
    isempty(m.accum) || (m.accum = gpu(m.accum))
    return
end

function cpu!(m::AWD_LSTM)
    m.layer = cpu(m.layer)
    isempty(m.accum) || (m.accum = cpu(m.accum))
    return
end

reset_masks!(awd::AWD_LSTM) = reset_masks!(awd.layer)

# Averaged Stochastic Gradient Descent Step
asgd_step!(i, l) = nothing

function asgd_step!(iter::Integer, layer::AWD_LSTM)
    if iter >= layer.T
        p = get_trainable_params([layer])
        avg_fact = 1/max(iter - layer.T + 1, 1)
        if avg_fact != 1
            layer.accum = layer.accum .+ Tracker.data.(p)
            for (ps, accum) in zip(p, layer.accum)
                Tracker.data(ps) .= avg_fact*accum
            end
        else
            layer.accum = deepcopy(Tracker.data.(p))   # Accumulator for ASGD
        end
    end
    return
end
####################################################################

"""
Variational Dropout
"""

########################## Varitional DropOut ######################
mutable struct VarDrop{F}
    p::F
    mask
    reset::Bool
    active::Bool
    VarDrop(probability::Float64=0.0) = new{typeof(probability)}(probability, Array{Float32, 2}(UndefInitializer(), 0, 0), true, true)
end

function (vd::VarDrop)(inp)
    vd.active || return inp
    if vd.reset
        vd.mask = drop_mask(inp, vd.p)
        vd.reset = false
    end
    return inp .* vd.mask
end

_testmode!(vd::VarDrop, test) = (vd.active = !test)

reset_masks!(vd::VarDrop) = (vd.reset = true;)

function gpu!(vd::VarDrop)
    vd.mask = gpu(vd.mask);
    return
end

function cpu!(vd::VarDrop)
    vd.mask = cpu(vd.mask);
    return
end
####################################################################

"""
Embeddings with varitional dropout
"""

################# Varitional Dropped Embeddings ####################
mutable struct DroppedEmbeddings{F}
    emb::TrackedArray
    p::F
    mask
    active::Bool
end

function DroppedEmbeddings(in::Integer, embed_size::Integer, probability::Float64=0.0;
    init = Flux.glorot_uniform)
        de = DroppedEmbeddings{typeof(probability)}(
            param(init(in, embed_size)),
            probability,
            drop_mask((in, embed_size), probability),
            true
        )
    return de
end

function (de::DroppedEmbeddings)(in::AbstractArray, tying::Bool=false)
    dropped = de.active ? de.emb .* de.mask : de.emb
    return tying ? dropped * in : transpose(dropped[in, :])
end

Flux.@treelike DroppedEmbeddings

_testmode!(de::DroppedEmbeddings, test) = (de.active = !test)

function gpu!(de::DroppedEmbeddings)
    de.emb = gpu(de.emb)
    de.mask = gpu(de.mask)
    return
end

function cpu!(de::DroppedEmbeddings)
    de.emb = cpu(de.emb)
    de.mask = cpu(de.mask)
    return
end

function reset_masks!(de::DroppedEmbeddings)
    de.mask = drop_mask(de.emb, de.p)
    return
end
####################################################################

"""
Concat-Pooled linear layer
"""

################# Concat Pooling Dense layer #######################
mutable struct PooledDense{F, S, T}
    W::S
    b::T
    σ::F
end

PooledDense(W, b) = PooledDense(W, b, identity)

function PooledDense(hidden_sz::Integer, out::Integer, σ = identity;
             initW = Flux.glorot_uniform, initb = (dims...) -> zeros(Float32, dims...))
return PooledDense(param(initW(out, hidden_sz*3)), param(initb(out)), σ)
end

Flux.@treelike PooledDense

function (a::PooledDense)(in)
    W, b, σ = a.W, a.b, a.σ
    in = cat(in..., dims=3)
    maxpool = maximum(in, dims=3)[:, :, 1]
    meanpool = (sum(in, dims=3)/size(in, 3))[:, :, 1]
    hc = cat(in[:, :, 1], maxpool, meanpool, dims=1)
    σ.(W*hc .+ b)
end

####################################################################

# Get the trainable params in the given layers
function get_trainable_params(layers)
    p = []
    function get_awd_params(awd::AWD_LSTM)
        return [awd.layer.cell.Wi,
        awd.layer.cell.Wh,
        awd.layer.cell.b]
    end
    for layer in layers
        layer isa AWD_LSTM && (append!(p, get_awd_params(layer)); continue)
        push!(p, layer)
    end
    return params(p...)
end


text:
if !@isdefined(BodePlots)
include("BodePlots.jl")
end

module DSPFilters

using DSP
using InspectDR
import InspectDR: Plot2D
using NumericIO
using FFTW
using Colors

#Colour constants
#-------------------------------------------------------------------------------
const black = RGB24(0, 0, 0); const white = RGB24(1, 1, 1)
const red = RGB24(1, 0, 0); const green = RGB24(0, 1, 0); const blue = RGB24(0, 0, 1)

#Convenience aliases
#-------------------------------------------------------------------------------
mag(x) = abs(x); phase(x) = rad2deg(angle(x))
SI(x) = formatted(x, :SI, ndigits=3)
const j = im

#==Equations
===============================================================================#

#==Main functions
===============================================================================#
function newplot()
	w = 500 #WANTCONST
	mplot = InspectDR.Multiplot(title="Filter Response")
	#Mag/phase plots look better with wider aspect ratio:
#	mplot.layout[:halloc_plot] = w; mplot.layout[:valloc_plot] = w/2.3 #Not needed with multi-y strips
	mplot.layout[:ncolumns] = 1

	#Initialize as Bode plot:
	plot_bode = add(mplot, Main.BodePlots.new(InspectDR.Plot))
	strip_mag, strip_phase = plot_bode.strips
		strip_mag.yext_full = InspectDR.PExtents1D(-90, 10)
		strip_phase.yext_full = InspectDR.PExtents1D(-200, 200)
		plot_bode.annotation.title = ""

	#Use linear frequency instead of actual Bode plot:
	plot_bode.xscale = InspectDR.AxisScale(:lin)

	return mplot
end

#Modify previously-generated plot:
#-------------------------------------------------------------------------------
#TODO: Look for off by 1 errors (esp. for FFT vector lengths, etc).
#(Frequencies: in Hz)
function update(mplot::InspectDR.Multiplot, fmax, filttypeid::Symbol, filtimplid::Symbol, order, fl, fh, passrip, stoprip)
	npts = 500 #WANTCONST
	ldata = line(color=blue, width=3, style=:solid) #WANTCONST
	lspec = line(color=red, width=1, style=:solid) #WANTCONST: Spectrum
	plot_bode = mplot.subplots[1] #WANTCONST
	rlow = 1/100; rhigh = 1-rlow #WANTCONST
	limfreqr(ratio) = clamp(ratio, rlow, rhigh)

	if fh < fl #Re-order fl/fh:
		(fl, fh) = (fh, fl)
	end

	#Generate filter type object:
	filttype = nothing
	if :highpass == filttypeid
		filttype = Highpass(limfreqr(fh/fmax))
	elseif :bandstop == filttypeid
		filttype = Bandstop(limfreqr(fl/fmax), limfreqr(fh/fmax))
	elseif :bandpass == filttypeid
		filttype = Bandpass(limfreqr(fl/fmax), limfreqr(fh/fmax))
	else #:lowpass
		filttype = Lowpass(limfreqr(fl/fmax))
	end

	#Generate filter implementation object
	filtimpl = nothing
	if :elliptic == filtimplid
		filtimpl = Elliptic(order, passrip, stoprip)
	elseif :chebyshevpass == filtimplid
		filtimpl = Chebyshev1(order, passrip)
	elseif :chebyshevstop == filtimplid
		filtimpl = Chebyshev2(order, stoprip)
	else
		filtimpl = Butterworth(order)
	end

	_filt = digitalfilter(filttype, filtimpl)

	#Generate mag/phase plot of discrete-time system:
	#TODO: find proper name for frequency variable.
	fz = collect(range(0, stop=pi, length=npts))
	tf_graph = freqz(_filt, fz)
	f = fz*(fmax/pi)
	Main.BodePlots.update(plot_bode, f, tf_graph)
	plot_bode.xext = InspectDR.PExtents1D() #Reset
	strip_mag, strip_phase = plot_bode.strips
		strip_mag.yext = InspectDR.PExtents1D() #Reset
		strip_phase.yext = InspectDR.PExtents1D() #Reset

	#Overlay spectrum using random data:
	xin = rand(100000)
	xout = filt(_filt, xin)
	spec = rfft(xout)/sqrt(length(xin)/2)
	fz = collect(range(0, stop=pi, length=length(spec)))
	f = fz*(fmax/pi)
	wfrm = add(plot_bode, f, spec, strip=1)
	wfrm.line = lspec

	#Hack to plot spectrum *before* (under) frequency graph:
	plot_bode.data = reverse(plot_bode.data)

	return mplot
end

end #module