text:
# Top level test script for Stan.jl
using CmdStan, StanMamba, Test, Statistics

println("Running tests for StanMamba-j1.0-v4.0.2:\n")

# Run execution_tests only if cmdstan is installed and CMDSTAN_HOME is set correctly.
execution_tests = [
  "test_bernoulli.jl",
  "test_bernoulli_nochains.jl"
]

if CMDSTAN_HOME != ""
  println("CMDSTAN_HOME set. Try to run tests.")
  @testset "StanMamba.jl" begin

    for my_test in execution_tests
        println("\n\n  * $(my_test) *\n")
        include(my_test)
    end

    println("\n")
  end
else
  println("\n\nCMDSTAN_HOME not set or found.")
  println("Skipping all tests that depend on CmdStan!\n")
end


text:
using Distributions, Documenter, GLM

makedocs(
    format = :html,
    sitename = "GLM",
    modules = [GLM]
)

deploydocs(
    repo   = "github.com/JuliaStats/GLM.jl.git",
    julia  = "0.6",
    target = "build",
    deps   = nothing,
    make   = nothing
)

ext: jl


text:
# discovering binary CUDA dependencies

using CompilerSupportLibraries_jll
using LazyArtifacts
import Libdl

## global state

const __toolkit_origin = Ref{Symbol}()

"""
    toolkit_origin()

Returns the origin of the CUDA toolkit in use (either :artifact, or :local).
"""
toolkit_origin() = @after_init(__toolkit_origin[])

const __toolkit_version = Ref{VersionNumber}()

"""
    toolkit_version()

Returns the version of the CUDA toolkit in use.
"""
toolkit_version() = @after_init(__toolkit_version[])

"""
    toolkit_release()

Returns the CUDA release part of the version as returned by [`version`](@ref).
"""
toolkit_release() = @after_init(VersionNumber(__toolkit_version[].major, __toolkit_version[].minor))

# if using a local toolkit, this contains a list of relevant directories
const __toolkit_dirs = Ref{Vector{String}}()

const __nvdisasm = Ref{String}()
const __compute_sanitizer = Ref{Union{Nothing,String}}()
const __libdevice = Ref{String}()
const __libcudadevrt = Ref{String}()
const __libcupti = Ref{Union{Nothing,String}}()
const __libnvtx = Ref{Union{Nothing,String}}()
const __libcublas = Ref{String}()
const __libcusparse = Ref{String}()
const __libcusolver = Ref{String}()
const __libcufft = Ref{String}()
const __libcurand = Ref{String}()
const __libcusolverMg = Ref{Union{Nothing,String}}(nothing)

nvdisasm() = @after_init(__nvdisasm[])
function compute_sanitizer()
    @after_init begin
        has_compute_sanitizer() ||
            error("This functionality is unavailabe as compute-sanitizer is missing.")
        __compute_sanitizer[]
    end
end
libdevice() = @after_init(__libdevice[])
libcudadevrt() = @after_init(__libcudadevrt[])
function libcupti()
    @after_init begin
        has_cupti() ||
            error("This functionality is unavailable as CUPTI is missing.")
        __libcupti[]
    end
end
function libnvtx()
    @after_init begin
        has_nvtx() ||
            error("This functionality is unavailable as NVTX is missing.")
        __libnvtx[]
    end
end

export has_compute_sanitizer, has_cupti, has_nvtx
has_compute_sanitizer() = @after_init(__compute_sanitizer[]) !== nothing
has_cupti() = @after_init(__libcupti[]) !== nothing
has_nvtx()  = @after_init(__libnvtx[]) !== nothing

libcublas()   = @after_init(__libcublas[])
libcusparse() = @after_init(__libcusparse[])
libcusolver() = @after_init(__libcusolver[])
libcufft()    = @after_init(__libcufft[])
libcurand()   = @after_init(__libcurand[])
function libcusolvermg()
    @after_init begin
        has_cusolvermg() ||
            error("This functionality is unavailabe as CUSOLVERMG is missing.")
        __libcusolverMg[]
    end
end

export has_cusolvermg
has_cusolvermg() = @after_init(__libcusolverMg[]) !== nothing

## discovery

# utilities to look up stuff in the artifact (at known locations)
artifact_binary(artifact_dir, name) = joinpath(artifact_dir, "bin", Sys.iswindows() ? "$name.exe" : name)
artifact_static_library(artifact_dir, name) = joinpath(artifact_dir, "lib", Sys.iswindows() ? "$name.lib" : "lib$name.a")
artifact_file(artifact_dir, path) = joinpath(artifact_dir, path)
function artifact_library(artifact, name, version)
    dir = joinpath(artifact, Sys.iswindows() ? "bin" : "lib")
    all_names = library_versioned_names(name, version)
    for name in all_names
        path = joinpath(dir, name)
        ispath(path) && return path
    end
    error("Could not find $name ($(join(all_names, ", ", " or "))) in $dir")
end

function artifact_cuda_library(artifact, library, toolkit_version)
    version = cuda_library_version(library, toolkit_version)
    name = get(cuda_library_names, library, library)
    artifact_library(artifact, name, version)
end

# CUDA

# workaround @artifact_str eagerness on unsupported platforms by passing a variable
function cuda_artifact(id, cuda::VersionNumber)
    platform = Base.BinaryPlatforms.HostPlatform()
    platform.tags["cuda"] = "$(cuda.major).$(cuda.minor)"
    @artifact_str(id, platform)
end

# NOTE: we don't use autogenerated JLLs, because we have multiple artifacts and need to
#       decide at run time (i.e. not via package dependencies) which one to use.
const cuda_toolkits = [
    (release=v"11.3", version=v"11.3.0",   preferred=true),
    (release=v"11.2", version=v"11.2.2",   preferred=true),
    (release=v"11.1", version=v"11.1.1",   preferred=true),
    (release=v"11.0", version=v"11.0.3",   preferred=true),
    (release=v"10.2", version=v"10.2.89",  preferred=true),
    (release=v"10.1", version=v"10.1.243", preferred=true),
]

function use_artifact_cuda()
    @debug "Trying to use artifacts..."

    # select compatible artifacts
    if haskey(ENV, "JULIA_CUDA_VERSION")
        wanted = VersionNumber(ENV["JULIA_CUDA_VERSION"])
        @debug "Selecting artifacts based on requested $wanted"
        candidate_toolkits = filter(cuda_toolkits) do toolkit
            toolkit.release == wanted || toolkit.version == wanted
        end
        isempty(candidate_toolkits) && @debug "Requested CUDA version $wanted is not provided by any artifact"
    else
        driver_release = release()
        @debug "Selecting artifacts based on driver compatibility $driver_release"
        candidate_toolkits = filter(cuda_toolkits) do toolkit
            toolkit.preferred && toolkit.release <= driver_release
        end
        isempty(candidate_toolkits) && @debug "CUDA driver compatibility $driver_release is not compatible with any artifact"
    end

    # download and install
    artifact = nothing
    for cuda in sort(candidate_toolkits; rev=true, by=toolkit->toolkit.version)
        try
            artifact = (version=cuda.version, dir=cuda_artifact("CUDA", cuda.version))
            break
        catch ex
            @debug "Could not load the CUDA $(cuda.version) artifact" exception=(ex,catch_backtrace())
        end
    end
    if artifact == nothing
        @debug "Could not find a compatible artifact."
        return false
    end

    __toolkit_version[] = artifact.version

    assert_artifact_file(path) =
        isfile(path) ||
            error("""Could not find $(basename(path)) in $(dirname(path))!
                     This is a bug; please file an issue with a verbose directory listing of $(dirname(path)).""")

    __nvdisasm[] = artifact_binary(artifact.dir, "nvdisasm")
    assert_artifact_file(__nvdisasm[])
    __compute_sanitizer[] = artifact_binary(artifact.dir, "compute-sanitizer")

    __libcupti[] = artifact_cuda_library(artifact.dir, "cupti", artifact.version)
    assert_artifact_file(__libcupti[])
    __libnvtx[] = artifact_cuda_library(artifact.dir, "nvtx", artifact.version)
    assert_artifact_file(__libnvtx[])

    __libcudadevrt[] = artifact_static_library(artifact.dir, "cudadevrt")
    assert_artifact_file(__libcudadevrt[])
    __libdevice[] = artifact_file(artifact.dir, joinpath("share", "libdevice", "libdevice.10.bc"))
    assert_artifact_file(__libdevice[])

    # HACK: eagerly load cublasLt, required by cublas (but with the same version), to prevent
    #       a system CUDA installation from messing with our artifacts (JuliaGPU/CUDA.jl#609)
    if artifact.version >= v"10.1"
        version = cuda_library_version("cublas", artifact.version)
        path = artifact_library(artifact.dir, "cublasLt", version)
        Libdl.dlopen(path)
    end

    if artifact.version >= v"10.1"
        lib_list =  ("cublas", "cusparse", "cusolver", "cufft", "curand", "cusolverMg")
    else
        lib_list =  ("cublas", "cusparse", "cusolver", "cufft", "curand")
    end
    for library in lib_list
        handle = getfield(CUDA, Symbol("__lib$library"))

        handle[] = artifact_cuda_library(artifact.dir, library, artifact.version)
        Libdl.dlopen(handle[])
    end

    @debug "Using CUDA $(__toolkit_version[]) from an artifact at $(artifact.dir)"
    __toolkit_origin[] = :artifact
    return true
end

function use_local_cuda()
    @debug "Trying to use local installation..."

    __toolkit_dirs[] = find_toolkit()
    let path = find_cuda_binary("nvdisasm", __toolkit_dirs[])
        if path === nothing
            @debug "Could not find nvdisasm"
            return false
        end
        __nvdisasm[] = path
    end

    __compute_sanitizer[] = find_cuda_binary("compute-sanitizer", __toolkit_dirs[])

    cuda_version = parse_toolkit_version("nvdisasm", __nvdisasm[])
    if cuda_version === nothing
        return false
    end

    for library in ("cublas", "cusparse", "cusolver", "cufft", "curand")
        handle = getfield(CUDA, Symbol("__lib$library"))

        path = find_cuda_library(library, __toolkit_dirs[], cuda_version)
        if path === nothing
            @debug "Could not find $library"
            return false
        end
        handle[] = path
    end

    # CUDA 11.1 Update 1 ships the same `nvdisasm` as 11.1 GA, so look at the version of
    # CUSOLVER (which has a handle-less version getter that does not initialize)
    # to be sure which CUDA we're dealing with (it only matters for CUPTI).
    if cuda_version == v"11.1.0"
        __libcusolver[] = find_cuda_library("cusolver", __toolkit_dirs[], v"11.1.0")
        if __libcusolver[] === nothing
            __libcusolver[] = find_cuda_library("cusolver", __toolkit_dirs[], v"11.1.1")
        end
        if __libcusolver[] === nothing
            @debug "Could not disambiguate CUDA 11.1 from Update 1 due to not finding CUSOLVER"
        else
            # nothing is initialized at this point, so we need to use raw ccalls.
            version = Ref{Cint}()
            @assert 0 == ccall((:cusolverGetVersion, __libcusolver[]), Cint, (Ref{Cint},), version)
            if version[] == 11001
                cuda_version = v"11.1.1"
            elseif version[] != 11000
                @debug "Could not disambiguate CUDA 11.1 from Update 1 with CUSOLVER version $(version[])"
            end
        end
    end
    __toolkit_version[] = cuda_version

    __libcupti[] = find_cuda_library("cupti", __toolkit_dirs[], cuda_version)
    __libnvtx[] = find_cuda_library("nvtx", __toolkit_dirs[], cuda_version)

    let path = find_libcudadevrt(__toolkit_dirs[])
        if path === nothing
            @debug "Could not find libcudadevrt"
            return false
        end
        __libcudadevrt[] = path
    end
    let path = find_libdevice(__toolkit_dirs[])
        if path === nothing
            @debug "Could not find libdevice"
            return false
        end
        __libdevice[] = path
    end

    for (library, required) in (("cublas", true), ("cusparse", true), ("cusolver", true),
                                ("cufft", true), ("curand", true), ("cusolverMg", false))
        handle = getfield(CUDA, Symbol("__lib$library"))
        path = find_cuda_library(library, __toolkit_dirs[], cuda_version)
        if path === nothing
            @debug "Could not find $library"
            if required
                return false
            else
                handle[] = nothing
            end
        else
            handle[] = path
        end
    end
    @debug "Found local CUDA $(cuda_version) at $(join(__toolkit_dirs[], ", "))"
    __toolkit_origin[] = :local
    return true
end

# CUDNN

export has_cudnn

const __libcudnn = Ref{Union{String,Nothing}}()
function find_cudnn()
    try
        version = v"8"
        if toolkit_origin() == :artifact
            artifact_dir = cuda_artifact("CUDNN", toolkit_version())
            path = artifact_library(artifact_dir, "cudnn", version)

            # HACK: eagerly open CUDNN sublibraries to avoid dlopen discoverability issues
            for sublibrary in ("ops_infer", "ops_train",
                            "cnn_infer", "cnn_train",
                            "adv_infer", "adv_train")
                sublibrary_path = artifact_library(artifact_dir, "cudnn_$(sublibrary)", version)
                Libdl.dlopen(sublibrary_path)
            end
            Libdl.dlopen(path)

            @debug "Using CUDNN from an artifact at $(artifact_dir)"
            __libcudnn[] = path
        else
            path = find_library("cudnn", version; locations=__toolkit_dirs[])
            path === nothing && error("Could not find a local CUDNN")

            # HACK: eagerly open CUDNN sublibraries to avoid dlopen discoverability issues
            for sublibrary in ("ops_infer", "ops_train",
                            "cnn_infer", "cnn_train",
                            "adv_infer", "adv_train")
                sublibrary_path = find_library("cudnn_$(sublibrary)", version; locations=__toolkit_dirs[])
                sublibrary_path === nothing && error("Could not find local CUDNN sublibrary $sublibrary")
                Libdl.dlopen(sublibrary_path)
            end
            Libdl.dlopen(path)

            @debug "Using local CUDNN at $(path)"
            __libcudnn[] = path
        end
    catch err
        @debug "Failed to load CUDNN" exception=(err, catch_backtrace())
        @warn "Could not find or load CUDNN; run with JULIA_DEBUG=CUDA for more details."
        __libcudnn[] = nothing
        return
    end

    # __libcudnn is set, so the CUDNN module is functional now

    cudnn_release = VersionNumber(CUDNN.version().major, CUDNN.version().minor)
    if cudnn_release < v"8.0"
        @warn "This version of CUDA.jl only supports CUDNN 8.0 or higher"
    end

    CUDNN.__runtime_init__()
end

function has_cudnn()
    @after_init begin
        isassigned(__libcudnn) || find_cudnn()
        __libcudnn[] !== nothing
    end
end

function libcudnn()
    has_cudnn() || error("This functionality is unavailabe as CUDNN is missing.")
    __libcudnn[]
end

# CUTENSOR

export has_cutensor

const __libcutensor = Ref{Union{String,Nothing}}()
function find_cutensor()
    try
        version = Sys.iswindows() ? nothing : v"1"  # cutensor.dll is unversioned on Windows
        if toolkit_origin() == :artifact
            artifact_dir = cuda_artifact("CUTENSOR", toolkit_version())
            path = artifact_library(artifact_dir, "cutensor", version)
            Libdl.dlopen(path)

            @debug "Using CUTENSOR from an artifact at $(artifact_dir)"
            __libcutensor[] = path
        else
            path = find_library("cutensor", version; locations=__toolkit_dirs[])
            if path === nothing
                path = find_library("cutensor"; locations=__toolkit_dirs[])
            end
            path === nothing && error("Could not find a local CUTENSOR")
            Libdl.dlopen(path)

            @debug "Using local CUTENSOR at $(path)"
            __libcutensor[] = path
        end
    catch err
        @debug "Failed to load CUTENSOR" exception=(err, catch_backtrace())
        @warn "Could not find or load CUTENSOR; run with JULIA_DEBUG=CUDA for more details."
        __libcutensor[] = nothing
        return
    end

    # __libcutensor is set, so the CUTENSOR module is functional now

    cutensor_release = VersionNumber(CUTENSOR.version().major, CUTENSOR.version().minor)
    if cutensor_release != v"1.3"
        @warn "This version of CUDA.jl only supports CUTENSOR 1.3"
    end
end

function has_cutensor()
    @after_init begin
        isassigned(__libcutensor) || find_cutensor()
        __libcutensor[] !== nothing
    end
end

function libcutensor()
    has_cutensor() || error("This functionality is unavailabe as CUTENSOR is missing.")
    __libcutensor[]
end

## lazy initialization

function __init_dependencies__()
    found = false

    # CI runs in a well-defined environment, so prefer a local CUDA installation there
    if getenv("CI", false) && !haskey(ENV, "JULIA_CUDA_USE_BINARYBUILDER")
        found = use_local_cuda()
    end

    if !found && getenv("JULIA_CUDA_USE_BINARYBUILDER", true)
        found = use_artifact_cuda()
    end

    # if the user didn't specifically request an artifact version, look for a local installation
    if !found && !haskey(ENV, "JULIA_CUDA_VERSION")
        found = use_local_cuda()
    end

    return found
end


text:
module AirSeaFluxes

using Dierckx
using LinearAlgebra

export simpleflux, bulkformulae, AOGCM1D

include("simpleflux.jl")
include("bulkformulae.jl")
include("holtslag.jl")
include("kpp.jl")
include("AOGCM1D_helpers.jl")
include("AOGCM1D.jl")

end # module


text:
pairs(n::Int64) = ((i, j) for i = 1:n for j = 1:i)
rpairs(n::Int64) = ((i,j) for i in 1:n for j in 1:n) # rectangular option to old pairs
spairs(n::Int64) = ((i, j) for i = 1:n for j = 1:(i-1)) # subdiagonal option to old pairs

triangle(i::Int64) = div(i*(i+1),2)
triangle(i::Int64,j::Int64) = i<j ? triangle(j-1)+i : triangle(i-1)+j

iiterator(n::Int64) = ((i,j,k,l) for (i,j) in pairs(n) for (k,l) in pairs(n) if triangle(i,j) <= triangle(k,l))

iindex(i::Int64,j::Int64,k::Int64,l::Int64) = triangle(triangle(i,j),triangle(k,l))
trace2(A,B) = sum(A.*B)

#=
function all_1e_ints( bfs::BasisSet, atoms::Atoms )
    n = length(bfs)
    S = Array{Float64}(undef,n,n)
    T = Array{Float64}(undef,n,n)
    V = Array{Float64}(undef,n,n)
    for j = 1:n
        for i = j:n
            S[i,j] = overlap(bfs[i], bfs[j])
            S[j,i] = S[i,j]
            T[i,j] = kinetic(bfs[i], bfs[j])
            T[j,i] = T[i,j]
            V[i,j] = nuclear_attraction(bfs[i], bfs[j], atoms)
            V[j,i] = V[i,j]
        end
    end
    return S,T,V
end
=#

function all_1e_ints(bfs::BasisSet, atoms::Atoms)
    n = length(bfs)
    S = Array{Float64}(undef,n,n)
    T = Array{Float64}(undef,n,n)
    V = Array{Float64}(undef,n,n)
    for (i,j) in pairs(n)
        a,b = bfs[i], bfs[j]
        S[i,j] = S[j,i] = overlap(a,b)
        T[i,j] = T[j,i] = kinetic(a,b)
        V[i,j] = V[j,i] = nuclear_attraction(a,b,atoms)
    end
    return S,T,V
end

#=
function all_twoe_ints(bfs,ERI=coulomb)
    n = length(bfs)
    totlen = div(n*(n+1)*(n*n+n+2),8)
    ints2e = zeros(totlen)  # using zeros instead of using Array constructor
    for i = 1:n
    for j = 1:n
        ij = triangle(i-1,j)
        for k = 1:n
        for l = 1:n
            kl = triangle(k-1,l)
            if ij <= kl
                ints2e[iindex(i,j,k,l)] = ERI(bfs[i], bfs[j], bfs[k], bfs[l])
            end
        end
        end
    end
    end
    return ints2e
end
=#

function all_twoe_ints(bfs,ERI=coulomb)
    n = length(bfs)
    totlen = div(n*(n+1)*(n*n+n+2),8)
    ints2e = Array{Float64}(undef,totlen)
    for (i,j,k,l) in iiterator(n)
        ints2e[iindex(i,j,k,l)] = ERI(bfs[i], bfs[j], bfs[k], bfs[l])
    end
    return ints2e
end

#=
function make2JmK(D::Array{Float64,2},Ints::Array{Float64,1})
    n = size(D,1)
    G = Array{Float64}(undef,n,n)
    D1 = reshape(D,n*n)
    temp = Array{Float64}(undef,n*n)
    for j = 1:n
    for i = j:n
        kl = 1
        for k = 1:n
            start = 1
            for l = start:n
                temp[kl] = 2*Ints[iindex(i,j,k,l)]-Ints[iindex(i,k,j,l)]
                kl += 1
            end
        end
        G[i,j] = dot(D1,temp)
        G[j,i] = G[i,j]
    end
    end
    return G
end
=#

function make2JmK(D::Array{Float64,2},Ints::Array{Float64,1})
    n = size(D,1)
    G = Array{Float64}(undef,n,n)
    D1 = reshape(D,n*n)
    temp = Array{Float64}(undef,n*n)
    for (i,j) in pairs(n)
        kl = 1
        for (k,l) in rpairs(n)
            temp[kl] = 2*Ints[iindex(i,j,k,l)]-Ints[iindex(i,k,j,l)]
            kl += 1
        end
        G[i,j] = G[j,i] = dot(D1,temp)
    end
    return G
end

dmat(U::Array{Float64,2},nocc::Int64) = U[:,1:nocc]*U[:,1:nocc]'

function RHF( atoms::Atoms, MaxIter::Int64=20, verbose::Bool=true)
    @printf("Starting RHF calculation:\n")

    Nelectrons = get_Nelectrons(atoms)
    nclosed,nopen = divrem( Int64(Nelectrons), 2 )
    @printf("Nel = %d, Nclosed = %d, nopen = %d\n", Nelectrons, nclosed, nopen)

    @printf("Building basis ...\n")
    bfs = build_basis_new(atoms)

    @printf("1e integrals ...\n")
    S,T,V = all_1e_ints(bfs,atoms)

    @printf("2e integrals ...\n")
    Ints = all_twoe_ints(bfs)

    h = T+V
    E,U = eigen(h,S)
    println("E = ", E)
    Enuke = nuclear_repulsion(atoms)
    Eold = 0.0
    Energy = 0.0

    IS_CONVERGED = false
    dEtot = 1.0
    for iter in 1:MaxIter
        D = dmat(U,nclosed)
        G = make2JmK(D,Ints)
        H = h + G
        E,U = eigen(H,S)
        Eone = trace2(D,h)
        Etwo = trace2(D,H)
        Energy = Enuke + Eone + Etwo
        dEtot = abs(Energy - Eold)
        #println("HF: $iter  $Energy : $Enuke  $Eone  $Etwo $dEtot")
        @printf("HF: %4d %18.10f %18.10f %18.10f: %18.10f %18.10e\n",
                iter, Enuke, Eone, Etwo, Energy, dEtot)
        if isapprox(Energy,Eold)
            IS_CONVERGED = true
            break
        end
        Eold  = Energy
    end

    if !IS_CONVERGED
        println("WARNING: RHF is not converged")
    end

    return Energy,E,U
end


text:
## Config (changes allowed)
latlen = 220          # lattice length
iters = 220           # iterations
start_rule = 113      # starting rule
order2_rule = ">|+"   # the rule of changing the rule
           #      - no changes (regular auomata)
           # >    - shift one to the right
           # >>   - shift two to the right
           # <    - shift one to the left
           # <<   - shift two to the left
           # +    - next rule (increases Wolfram's number)
           # -    - previous rule (decreases Wolfram's number)
           # |    - flip the rule by bits (mirror)
           # !    - perform "not" operation on each bit
           # +>>  - combination (increase rule and move two rules to the right). + operation cometh first
init = "r"        # initial distribution of cells
           # .  - dot
           # r  - random
           # i  - inverse
           # :  - alternating black and white cells
           # .i - combination (inverse dot)
           # a  - get from init_arr
init_arr = []     # initial distribution of cells (leave blank if init is not "a").
                  # watch out, init_arr needs to be of length latlen

## Some interesting configurations
#=
starting with dot:
# 60 >>
# 60 >>|<<
# 60 >|<
# 60 <|>
# 60 >+<
# 101
# 101 >
# 101 >>
# 101 >|
# 101 |
# 105 >>
# 105 >|>
# 105 >|<
random start:
# 40 +
# 40 >+< any even?
# 50 >|<
# 60 >|<
# 105 >>
# 105 >>+
# 110 +<
# 110 <
# 110 +
# 110 -
# 110 --
# 110 >+++<
=#

## Config processing (changes not recommended)
# rule is a vector of integers (1s and 0s)
rule = [Int(i)-48 for i::Char in Base.bin(Unsigned(start_rule), 8, false)]

# a function to change the rule number by one
function next_rule!(rulevec::Vector{Int}, step::Int)::Nothing
    num::Int = parse(Int, join(rulevec), base=2)
    num += step
    if num > 255
        num -= 255
    elseif num < 0
        num = 256 + num
    end
    b::String = Base.bin(Unsigned(num), 8, false)
    for i in 1:length(rulevec)
        rulevec[i] = Int(b[i])-48
    end
    nothing
end

# a function to change the rule with a bitwise "not"
function not!(rulevec::Vector{Int})::Nothing
    for i in 1:length(rulevec)
        rulevec[i] = Int(rulevec[i] == 0)
    end
    nothing
end

# a mapping from characters to actions (order2 rules)
char2rule = Dict([
    '+' => r -> next_rule!(r, 1),
    '-' => r -> next_rule!(r, -1),
    '>' => r -> circshift!(r, copy(r), 1),
    '<' => r -> circshift!(r, copy(r), -1),
    '|' => r -> reverse!(r),
    '!' => r -> not!(r)
])

## Import and setup for visualization
using PyCall
import PyPlot; const plt = PyPlot
rcParams = PyDict(plt.matplotlib."rcParams")
rcParams["image.cmap"] = "binary"

## Data init (first line)
data = [] # comment this line to append the new auomaton graph to the previous one
if occursin("r", init)
    push!(data, [rand((0, 1)) for i::Int in 1:latlen])
elseif occursin(".", init)
    push!(data, [((i==round(latlen/2)) ? 1 : 0) for i::Int in 1:latlen])
elseif occursin("a", init)
    push!(data, init_arr)
elseif occursin(":", init)
    push!(data, [i%2 for i::Int in 1:latlen])
end
if occursin("i", init)
    data = [[abs(i-1) for i::Int in data[1]]]
end

## The magic
past = length(data)-1
for iter::Int in 1:iters
    println(parse(Int, join(rule), base=2), rule)
    new_line::Vector{Int} = []
    for i::Int in 1:length(data[past+iter])
        neigh::String = ""
        for j::Int in i-1:i+1
            aj::Unsigned = j # actual index to look at
            if (j <= 0)
                aj = length(data[past+iter])-1+j
            elseif (j > length(data[past+iter]))
                aj = j - length(data[past+iter])
            end
            neigh *= string(data[past+iter][aj])
        end
        push!(new_line, rule[8-parse(Int, neigh, base=2)])
    end
    push!(data, new_line)
    # change the rule according to the order2_rule
    for c::Char in order2_rule
        char2rule[c](rule)
    end
end

## Visualization
fig, ax = plt.subplots(figsize=(80, 80))
ax.matshow(data)
ax.axis(false)
display(plt.gcf()) # for Juno Plots window
plt.show()         # for terminal usage

max_stars_repo_name: gojakuch/dynamic-rule-cellular-automata

max_issues_repo_name: gojakuch/dynamic-rule-cellular-automata

max_forks_repo_name: gojakuch/dynamic-rule-cellular-automata


text:
# conversions among concrete types
# ================================

# AbstractVector
# --------------

# Because `AbstractVector <: AbstractArray`, and more especially `Base.typename(Vector) == Array`,
# we are dispatching on AbstractArray instead of AbstractVector

Base.convert(T::Type{<:AbstractArray}, x::Identity) = convert(T, [x.value])
Base.convert(T::Type{<:AbstractArray}, x::Const) = convert(T, [])
# ContextManager is executed if someone asks for a Vector from ContextManager
Base.convert(T::Type{<:AbstractArray}, x::ContextManager) = convert(T, [x(identity)])

# Const
# -----

function Base.convert(::Type{Const}, x::Vector)
  @assert isempty(x) "can only convert empty Vector to Nothing, got $(x)"
  Const(nothing)  # == Option()
end

# Identity
# --------

# ContextManager is executed
Base.convert(::Type{<:Identity}, x::ContextManager) = Identity(run(x))

# ContextManager
# --------------

Base.convert(::Type{<:ContextManager}, x::Identity) = @ContextManager cont -> cont(x.value)

# Note that we do not support convert from Const to ContextManager.
# ContextManager is similar to Identity a Container which always contains a single element,
# however, unlike Identity, we cannot simply pass-through Const and preserve the empty container by this trick.
# To map over ContextManager, you NEED to create a new ContextManager, if you don't want to run the contextmanager
# immediately. Hence, also flatmap needs to create another ContextManager, i.e. single element container.

# Vector is neither supported, as it rather corresponds to Continuable which can have several calls to `cont`
# while ContextManager should only have one call to `cont`.

# Option, Either, Try
# -------------------

Base.convert(::Type{Either{L, R}}, x::Identity) where {L, R} = Identity(Base.convert(R, x.value))
Base.convert(::Type{Either{L, R}}, x::Identity{R}) where {L, R} = x  # nothing to convert
Base.convert(::Type{Either{L, R}}, x::Const) where {L, R} = Const(Base.convert(L, x.value))
Base.convert(::Type{Either{L, R}}, x::Const{L}) where {L, R} = x  # nothing to convert

Base.convert(::Type{Either{<:Any, R}}, x::Identity{R}) where {R} = x  # nothing to convert
Base.convert(::Type{Either{<:Any, R}}, x::Identity) where {R} = Identity(Base.convert(R, x.value))
Base.convert(::Type{Either{L, <:Any}}, x::Const{L}) where {L} = x  # nothing to convert
Base.convert(::Type{Either{L, <:Any}}, x::Const) where {L} = Const(Base.convert(L, x.value))

Base.convert(::Type{Either}, x::Identity) = x
Base.convert(::Type{Either}, x::Const) = x


text:
@testset "Gaussian Quadrature" begin
    # NOTE: It is assumed p = 1:5 is sufficient for testing
    @testset "RefLine" begin
        for T ∈ Floats
            for p = 1:5
                weights, points = gauss_quadrature(Val(:legendre), RefLine(), Val(p), T)
                @test abs(sum(weights) - 1) < 10*eps(T)
                # Test the integration of a polynomial of degree 2p-1, which should
                # be exact for p integration points.
                # 1
                # ∫xᵖdx = 1/(p+1) ≈ ∑ wᵢxᵢᵖ
                # 0
                exact = 1/big(p+1)
                approx = zero(T)
                for i = 1:p
                    approx += weights[i]*(points[i][1]^p)
                end
                rel_err = T((approx - exact)/exact)
                @test rel_err < 10*eps(T)
            end
        end
    end

    @testset "RefSquare" begin
        for T ∈ Floats
            for p = 1:5
                weights, points = gauss_quadrature(Val(:legendre), RefSquare(), Val(p), T)
                @test abs(sum(weights) - 1) < 10*eps(T)
                #
                # q = 2p-1
                #
                # 11                                   5 + (-1)^q + 2q
                # ∫∫(r^q + s^q + (1 - r - s)^q ds dr = --------------
                # 00                                     (q+1)(q+2)
                q = 2p-1
                exact = big((5 + (-1)^q + 2q))/big((q+1)*(q+2))
                approx = zero(T)
                for i = 1:p^2
                    r, s = points[i]
                    approx += weights[i]*(r^q + s^q + (1 - r - s)^q)
                end
                rel_err = T((approx - exact)/exact)
                @test rel_err < 10*eps(T)
            end
        end
    end

    @testset "RefCube" begin
        for T ∈ Floats
            for p = 1:5
                weights, points = gauss_quadrature(Val(:legendre), RefCube(), Val(p), T)
                @test abs(sum(weights) - 1) < 10*eps(T)
                #
                # q = 2p-1
                #
                # 111
                # ∫∫∫(r^q + s^q + t^q) dt ds dr = 3/(q+1)
                # 000
                q = 2p-1
                exact = 3/big(q + 1)
                approx = zero(T)
                for i = 1:p^3
                    r, s, t = points[i]
                    approx += weights[i]*(r^q + s^q + t^q)
                end
                rel_err = T((approx - exact)/exact)
                @test rel_err < 10*eps(T)
            end
        end
    end

    @testset  "RefTriangle" begin
        for T ∈ Floats
            for p = 1:5
                weights, points = gauss_quadrature(Val(:legendre), RefTriangle(), Val(p), T)
                @test abs(sum(weights) - 1//2) < 10*eps(T)
                # This quadrature is only exact to polynomial degree p, not 2p-1!
                #
                # 1 1-r
                # ∫  ∫ (r^p + s^p + (1 - r - s)^p) ds dr = 3/((p+1)(p+2))
                # 0  0
                #
                exact = 3/big((p+1)*(p+2))
                approx = zero(T)
                for i = 1:length(weights)
                    r, s = points[i]
                    approx += weights[i]*(r^p + s^p + (1 - r - s)^p)
                end
                rel_err = T((approx - exact)/exact)
                @test rel_err < 10*eps(T)
            end
        end
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
include("property.jl");
include("geometry.jl");
include("model.jl");

using Plots;
pyplot();
begin
    local n_particle = 8192;
    local grid_size = 128;
    local rho = 1;
    local nu = 0.3;
    local E = 400;
    local dx::Float32 = 1 / grid_size;
    local dt = 2f-4;
    local g = -9.8f0;

    local m = Material(E, nu, rho);
    local grid = Grid2D(grid_size, grid_size, dx);
    local p = Property2D(m);
    local model = Model(grid, g);
    for i  = 1 : n_particle
        local x = Vector{Float32}([rand() * 0.4 + 0.2, rand() * 0.4 + 0.2]);
        local v = Vector{Float32}([0, -1])
        local particle = Particle2D(x, v, (dx * 0.5f0)^2, p);
        addParticle(model, particle);
    end
    for i = 1 : 5000
        solve(model, dt);
        println(i);
        if i % 50 == 0
            local j::Int32 = floor(i / 50);
            local r = getPositions(model);
            local fig = plot(r[:, 1], r[:, 2], seriestype=:scatter, aspect_ratio=:equal, xlims = (0, 1), ylims = (0, 1), size = (800, 800));
            #local r = getGridPositions(model);
            #local v = getGridVelocities(model);
            #fig = quiver(r[:, 1], r[:, 2], quiver=(v[:, 1] * 0.05, v[:, 2] * 0.05), size = (600 * 4, 600 * 4), aspect_ratio=:equal, xlims = (0, 1), ylims = (0, 1));
            savefig(fig, "image/$j.png");
        end
    end
end


text:
module ImportAll

macro importall(moduleIn, verbose=false)
    quote
        try
            $(esc(moduleIn))
        catch
            import $(moduleIn)
        end
        for name in names($(esc(moduleIn)))
            @eval import .$(moduleIn): $(Expr(:$, :name))
        end
    end
end

export @importall

end


text:
using DiffEqBase, ClimaTimeSteppers, LinearAlgebra, StaticArrays

if !@isdefined(ArrayType)
    ArrayType = Array
end

const_prob = ODEProblem{true}(
    (du,u,p,t,α=true,β=false) -> (du .= α .* p .+ β .* du),
    [0.0],(0.0,1.0),2.0)
const_prob_fe = ODEProblem(
        ForwardEulerODEFunction((un,u,p,t,dt) -> (un .= u .+ dt.* p)),
        [0.0],(0.0,1.0),2.0)
function const_sol(u0,p,t)
    u0 .+ p*t
end

"""
Single variable linear ODE

```math
\\frac{du}{dt} = αu
```
with initial condition ``u_0=\\frac{1}{2}``, parameter ``α=1.01``, and solution
```math
u(t) = u_0 e^{αt}
```

This is an in-place variant of the one from DiffEqProblemLibrary.jl.
"""
linear_prob = IncrementingODEProblem{true}(
    (du,u,p,t,α=true,β=false) -> (du .= α .* p .* u .+ β .* du),
    [1/2],(0.0,1.0),1.01)
linear_prob_fe = ODEProblem(
    ForwardEulerODEFunction((un,u,p,t,dt) -> (un .= u .+ dt .* p .* u)),
    [1.0],(0.0,1.0),-0.2)

linear_prob_wfactt = ODEProblem(
        ODEFunction(
          (du,u,p,t,α=true,β=false) -> (du .= α .* p .* u .+ β .* du);
          jac_prototype=zeros(1,1),
          Wfact_t = (W,u,p,γ,t) -> (W[1,1]=1/γ-p),
        ),
        [1/2],(0.0,1.0),-0.2)

split_linear_prob_wfact_split = ODEProblem(
    SplitFunction(
        ODEFunction(
            (du,u,p,t) -> (du .= real(p) .* u);
            jac_prototype=zeros(ComplexF64,1,1),
            Wfact = (W,u,p,γ,t) -> (W[1,1]=1-γ*real(p)),
        ),
        (ux, u, p, t, dt) -> (ux .= ux .+ dt .* imag(p) * im .* u),
    ),
    [1/2 + 0.0*im],(0.0,1.0),-0.2+0.1*im)

# DiffEqProblemLibrary.jl uses the `analytic` argument to ODEFunction to store the exact solution
# IncrementingODEFunction doesn't have that
function linear_sol(u0,p,t)
    u0 .* exp(p*t)
end

"""
Two variable linear ODE

```math
\\frac{du_1}{dt} = α u_2 \\quad \\frac{du_2}{dt} = -α u_1
```
with initial condition ``u_0=[0,1]``, parameter ``α=2``, and solution
```math
u(t) = [cos(αt) sin(αt); -sin(αt) cos(αt) ] u_0
```
"""
sincos_prob = IncrementingODEProblem{true}(
    (du,u,p,t,α=true,β=false) -> (du[1] = α*p*u[2]+β*du[1]; du[2] = -α*p*u[1]+β*du[2]),
    [0.0,1.0], (0.0,1.0), 2.0)
sincos_prob_fe = ODEProblem(
    ForwardEulerODEFunction((un,u,p,t,dt) -> (un[1] = u[1] + dt*p*u[2]; un[2] = u[2]-dt*p*u[1])),
    [0.0,1.0], (0.0,1.0), 2.0)

function sincos_sol(u0,p,t)
    s,c = sincos(p*t)
    [c s; -s c]*u0
end

"""
IMEX problem with autonomous linear part

```math
\\frac{du}{dt} = u + cos(t)/α
```
with initial condition ``u_0 = 1/2``, parameter ``α=4``, and solution
```math
u(t) = \frac{e^t + sin(t) - cos(t)}{2 α} + u_0 e^t
```
"""
imex_autonomous_prob = SplitODEProblem(
    (du, u, p, t, α=true, β=false) -> (du .= α .* u        .+ β .* du),
    (du, u, p, t, α=true, β=false) -> (du .= α .* cos(t)/p .+ β .* du),
    ArrayType([0.5]), (0.0,1.0), 4.0)

function linsolve_direct(::Type{Val{:init}}, f, u0; kwargs...)
    function _linsolve!(x, A, b, update_matrix = false; kwargs...)
        x .= A \ b
    end
end

imex_autonomous_prob_jac = ODEProblem(
        ODEFunction(
            (du, u, p, t, α=true, β=false) -> (du .= α .* (u .+ cos(t)/p) .+ β .* du),
            jac_prototype = zeros(1,1),
            Wfact = (W,u,p,γ,t) -> W[1,1] = 1-γ,
        ),
        ArrayType([0.5]), (0.0,1.0), 4.0)

function imex_autonomous_sol(u0,p,t)
    (exp(t)  + sin(t) - cos(t)) / 2p .+ exp(t) .* u0
end

imex_nonautonomous_prob = SplitODEProblem(
    (du, u, p, t, α=true, β=false) -> (du .= α .* cos(t) * u .+ β .* du),
    (du, u, p, t, α=true, β=false) -> (du .= α .* cos(t)/p   .+ β .* du),
    ArrayType([0.5]), (0.0,2.0), 4.0)

function imex_nonautonomous_sol(u0,p,t)
    (exp(sin(t)) .* (1 .+ p .* u0) .- 1) ./ p
end

"""
Test problem (4.2) from RobertsSarsharSandu2018arxiv
@article{RobertsSarsharSandu2018arxiv,
    title={Coupled Multirate Infinitesimal GARK Schemes for Stiff Systems with
            Multiple Time Scales},
    author={Roberts, Steven and Sarshar, Arash and Sandu, Adrian},
    journal={arXiv preprint arXiv:1812.00808},
    year={2019}
}

Note: The actual rates are all over the place with this test and passing largely
        depends on final dt size
"""
kpr_param = (
    ω = 100,
    λf = -10,
    λs = -1,
    ξ = 0.1,
    α = 1,
)

function kpr_rhs(Q,param,t)
    yf, ys = Q
    ω, λf, λs, ξ, α = param

    ηfs = ((1 - ξ) / α) * (λf - λs)
    ηsf = -ξ * α * (λf - λs)
    Ω = @SMatrix [
        λf ηfs
        ηsf λs
    ]

    g = @SVector [
        (-3 + yf^2 - cos(ω * t)) / 2yf,
        (-2 + ys^2 - cos(t)) / 2ys,
    ]
    h = @SVector [
        ω * sin(ω * t) / 2yf,
        sin(t) / 2ys,
    ]
    return Ω * g - h
end

function kpr_fast!(dQ, Q, param, t, α=1, β=0)
    dQf,_ = kpr_rhs(Q,param,t)
    dQ[1] = α*dQf + β*dQ[1]
    dQ[2] = β*dQ[2]
end

function kpr_slow!(dQ, Q, param, t, α=1, β=0)
    _,dQs = kpr_rhs(Q,param,t)
    dQ[1] = β*dQ[1]
    dQ[2] = α*dQs + β*dQ[2]
end

function kpr_sol(u0,param,t)
    #@assert u0 == [sqrt(4) sqrt(3)]
    ω, λf, λs, ξ, α = param
    return [
        sqrt(3 + cos(ω * t))
        sqrt(2 + cos(t))
    ]
end

kpr_multirate_prob = SplitODEProblem(
    IncrementingODEFunction{true}(kpr_fast!),
    IncrementingODEFunction{true}(kpr_slow!),
    [sqrt(4), sqrt(3)], (0.0, 5π/2), kpr_param,
)

kpr_singlerate_prob = IncrementingODEProblem{true}(
    (dQ,Q,param,t,α=1,β=0) -> (dQ .= α .* kpr_rhs(Q,param,t) .+ β .* dQ),
    [sqrt(4), sqrt(3)], (0.0, 5π/2), kpr_param,
)


text:
# Deprecations in OIFITS module.

import Base: @deprecate

# Deprecated in v0.2
@deprecate readtable(ff::FITSFile) read_table(ff::FITSFile)


text:
"""
    visit(operation, obj; print::Bool=false)

Scan `obj` and all of its "sub-objects" (e.g., functions if `obj::Module`,
methods if `obj::Function`, etc.) recursively.
`operation(x)` should return `true` if `visit` should descend
into "sub-objects" of `x`.

If `print` is `true`, each visited object is printed to standard output.

# Example

To collect all MethodInstances of a function,

```jldoctest; setup=:(using MethodAnalysis), filter=r"[0-9][0-9]"
julia> mis = Core.MethodInstance[];

julia> visit(findfirst) do x
           if isa(x, Core.MethodInstance)
               push!(mis, x)
               return false
           end
           true
       end

julia> length(mis)
34
```

The exact number of MethodInstances will depend on what code you've run in your Julia session.
"""
visit(@nospecialize(operation), @nospecialize(obj); print::Bool=false) =
    _visit(operation, obj, IdSet{Any}(), print)

"""
    visit(operation; print::Bool=false)

Scan all loaded modules with `operation`.
See [`visit(operation, obj)`](@ref) for further detail.

# Example

Collect all loaded modules, even if they are internal.

```jldoctest; setup=:(using MethodAnalysis), filter=r"[0-9][0-9][0-9]-element"
julia> mods = Module[];

julia> visit() do x
           if isa(x, Module)
               push!(mods, x)
               return true
           end
           false
       end
```

julia> mods
113-element Array{Module,1}:
 Random
 Random.DSFMT
[...]
"""
function visit(@nospecialize(operation); print::Bool=false)
    visited = IdSet{Any}()
    for mod in Base.loaded_modules_array()
        _visit(operation, mod, visited, print)
    end
    return nothing
end

# These are non-keyword functions due to https://github.com/JuliaLang/julia/issues/34516

function _visit(@nospecialize(operation), mod::Module, visited::IdSet{Any}, print::Bool)
    mod ∈ visited && return nothing
    push!(visited, mod)
    print && println("Module ", mod)
    if operation(mod)
        for nm in names(mod; all=true)
            if isdefined(mod, nm)
                obj = getfield(mod, nm)
                _visit(operation, obj, visited, print)
            end
        end
    end
    return nothing
end

function _visit(@nospecialize(operation), @nospecialize(f::Callable), visited::IdSet{Any}, print::Bool)
    f ∈ visited && return nothing
    f === Vararg && return nothing   # methods(Varargs) errors due to Type{Vararg}
    push!(visited, f)
    print && println("Callable ", f)
    if operation(f)
        ml = methods(f)
        _visit(operation, ml.mt, visited, print)
        for m in ml.ms
            _visit(operation, m, visited, print)
        end
    end
    return nothing
end

function _visit(@nospecialize(operation), ml::Base.MethodList, visited::IdSet{Any}, print::Bool)
    ml ∈ visited && return nothing
    push!(visited, ml)
    print && println("MethodList ", ml)
    _visit(operation, ml.mt, visited, print)
    for m in ml.ms
        _visit(operation, m, visited, print)
    end
    return nothing
end

function _visit(@nospecialize(operation), mt::MethodTable, visited::IdSet{Any}, print::Bool)
    mt ∈ visited && return nothing
    push!(visited, mt)
    print && println("MethodTable ", mt)
    operation(mt)
    return nothing
end

function _visit(@nospecialize(operation), m::Method, visited::IdSet{Any}, print::Bool)
    m ∈ visited && return nothing
    push!(visited, m)
    print && println("Method ", m)
    if operation(m)
        for fn in (:specializations, :invokes)
            if isdefined(m, fn)
                spec = getfield(m, fn)
                if spec === nothing
                elseif isa(spec, Core.TypeMapEntry) || isa(spec, Core.TypeMapLevel)
                    Base.visit(spec) do mi
                        _visit(operation, mi, visited, print)
                    end
                elseif isa(spec, Core.SimpleVector)
                    _visit(operation, spec, visited, print)
                else
                    error("unhandled type ", typeof(spec), ": ", spec)
                end
            end
        end
    end
    return nothing
end

function _visit(@nospecialize(operation), sv::SimpleVector, visited::IdSet{Any}, print::Bool)
    for i = 1:length(sv)
        if isassigned(sv, i)
            _visit(operation, sv[i], visited, print)
        end
    end
    return nothing
end

function _visit(@nospecialize(operation), mi::MethodInstance, visited::IdSet{Any}, print::Bool)
    mi ∈ visited && return nothing
    push!(visited, mi)
    print && println(mi)
    if operation(mi)
        if isdefined(mi, :cache)
            _visit(operation, mi.cache, visited, print)
        end
    end
    return nothing
end

if isdefined(Core, :CodeInstance)
    function _visit(@nospecialize(operation), ci::Core.CodeInstance, visited::IdSet{Any}, print::Bool)
        ci ∈ visited && return nothing
        push!(visited, ci)
        print && println(ci)
        if operation(ci)
            if isdefined(ci, :next)
                _visit(operation, ci.next, visited, print)
            end
        end
        return nothing
    end
end

_visit(@nospecialize(operation), @nospecialize(x), visited::IdSet{Any}, print::Bool) = nothing

"""
    visit_backedges(operation, obj)

Visit the backedges of `obj` and apply `operation` to each.
`operation` may need to be able to handle two call forms, `operation(mi)` and
`operation(sig=>mi)`, where `mi` is a `MethodInstance` and `sig` is a `Tuple`-type.
The latter arises from `MethodTable` backedges and can be ignored if `obj` does not
contain `MethodTable`s.

`operation(edge)` should return `true` if the backedges of `edge` should in turn be visited,
`false` otherwise.
However, no `MethodInstance` will be visited more than once.

The set of visited objects includes `obj` itself.
For example, `visit_backedges(operation, f::Function)` will visit all methods of `f`,
and this in turn will visit all MethodInstances of these methods.
"""
function visit_backedges(@nospecialize(operation), obj)
    visited = IdSet{MethodInstance}()
    visit_backedges(operation, obj, visited)
end

function visit_backedges(@nospecialize(operation), obj, visited::IdSet{MethodInstance})
    function opwrapper(@nospecialize(x))
        if isa(x, MethodInstance)
            _visit_backedges(operation, x, visited)
            return false
        elseif isa(x, MethodTable)
            mt = x::MethodTable
            if isdefined(mt, :backedges)
                sigmis = mt.backedges::Vector{Any}
                for i = 1:2:length(sigmis)
                    sig, mi = sigmis[i], sigmis[i+1]
                    _visit_backedges(operation, Pair{Any,MethodInstance}(sig, mi), visited)
                end
            end
            return false
        end
        return true
    end

    visit(opwrapper, obj)
end

function _visit_backedges(@nospecialize(operation), mi::MethodInstance, visited)
    mi ∈ visited && return nothing
    push!(visited, mi)
    if operation(mi) && isdefined(mi, :backedges)
        for edge in mi.backedges
            _visit_backedges(operation, edge, visited)
        end
    end
    return nothing
end

function _visit_backedges(@nospecialize(operation), misig::Pair{Any,MethodInstance}, visited)
    if operation(misig)
        _visit_backedges(operation, misig.second, visited)
    end
    return nothing
end


text:
@inline function createSegments(coords::Vector{Position})
    segments = []
    if length(coords) === 2
        push!(segments, Linestring([coords[1], coords[2]]))
        return segments
    end

    for i in 1:length(coords)-1
        line = LineString([coords[i], coords[i + 1]])
        push!(segments, line)
    end
    return segments
end

function linesegmentFeature(geojson::Union{LineString, Polygon}, result::Vector{LineString})
    coords = geojson.coordinates

    segments = createSegments(coords)

    for s in segments
        push!(result, s)
    end

    return result
end

"""
    linesegment(geojson::LineString)

Create a 2-vertex LineString segments from a LineString.

# Examples
```jldoctest
julia> line = LineString([[0, 1], [2, 3], [3, 3]])
LineString(Array{Float64,1}[[0.0, 1.0], [2.0, 3.0], [3.0, 3.0]])

julia> linesegment(line)
2-element Array{LineString,1}:
 LineString(Array{Float64,1}[[0.0, 1.0], [2.0, 3.0]])
 LineString(Array{Float64,1}[[2.0, 3.0], [3.0, 3.0]])
```
"""
function linesegment(geojson::LineString)
    result::Vector{LineString} = []
    linesegmentFeature(geojson, result)

    return result
end

function to360(a::Real)
    b = a % 360
    b < 0 && (b += 360)

    return b
end

"""
    linearc(center::Point, radius::Real, bearing1::Real, bearing2::Real, steps::Real=64., units::String="kilometers")

Create a circular arc, of a circle of the given radius and center point, between bearing1 and bearing2;
0 bearing is North of center point, positive clockwise.

# Examples
```jldoctest
julia> point = Point([35, 45])
Point([35.0, 45.0])

julia> linearc(point, 5, 0, 10, 5)
LineString(Array{Float64,1}[[35.0, 45.045], [35.0111, 45.0443]])

julia> linearc(point, 5, 0, 270, 5)
LineString(Array{Float64,1}[[35.0, 45.045], [35.0605, 45.0139], [35.0374, 44.9636], [34.9626, 44.9636], [34.9364, 45.0]])
```
"""
function linearc(center::Point, radius::Real, bearing1::Real, bearing2::Real, steps::Real=64., units::String="kilometers")
    angle1 = to360(bearing1)
    angle2 = to360(bearing2)

    angle1 === angle2 && return LineString(circle(center=center, radius=radius, steps=steps, units=units).coordinates[1])

    startdeg = angle1
    enddeg = (angle1 < angle2) ?  angle2 :  angle2 + 360
    α = startdeg

    coords = []
    i = 0

    while α < enddeg
        push!(coords, destination(center.coordinates, radius, α).coordinates)
        i += 1
        α = startdeg + i * 360 / steps
    end

    α > enddeg && push!(coords, destination(center.coordinates, radius, enddeg).coordinates)

    return LineString(coords)
end


text:
using Pkg
Pkg.activate(pwd())

# # More complex networks
# ## Preparing data

using MLDatasets
using MLDatasets: MNIST

T = Float32
X_train, y_train = MNIST.traindata(T)
X_test, y_test = MNIST.testdata(T)

# ### Exercise:
# Plot the first 15 images of the digit 0 from the training set.
#
# **Hints:**:
# - The `ImageInspector` package written earlier provides the function
# `imageplot(X_train, inds; nrows=3)`, where `inds` are the desired indices.
# - To find the correct indices, use the function `findall`.
#
# ---
# ### Solution:

# ---
#
# ### Exercise:
# Write function `reshape_data`, which reshapes `X_train` and `X_test` into the correct
# size required by Flux.
#
# **Hint**: The function should work only on inputs with the correct size. This can be
# achieved by specifying the correct input type `X::AbstractArray{<:Real, 3}`.
#
# ---
# ### Solution:

# ---

using Flux
using Flux: onehotbatch, onecold

function load_data(dataset; T=Float32, onehot=false, classes=0:9)
    X_train, y_train = dataset.traindata(T)
    X_test, y_test = dataset.testdata(T)

    X_train = reshape_data(X_train)
    X_test = reshape_data(X_test)

    if onehot
        y_train = onehotbatch(y_train, classes)
        y_test = onehotbatch(y_test, classes)
    end

    return X_train, y_train, X_test, y_test
end

#+

X_train, y_train, X_test, y_test = load_data(MLDatasets.MNIST; T=T, onehot=true)

# ### Exercise:
# Try to load the CIFAR10 dataset via the `load_data` function and fix the error in one
# line of code.
#
# **Hint**: Use ` dataset = MLDatasets.CIFAR10`.
#
# ---
# ### Solution:

# ---
#
# ## Training and storing the network
# ### Exercise:
# Use the help of the function `DataLoader` to split the dataset into minibatches.
#
# **Hint**: It needs to be imported from Flux via `using Flux.Data: DataLoader`.
#
# ---
# ### Solution:

# ---

using Base.Iterators: partition
using Random

batches = map(partition(randperm(size(y_train, 2)), batchsize)) do inds
    return (X_train[:, :, :, inds], y_train[:, inds])
end

#+

batches2 = [(X_train[:, :, :, inds], y_train[:, inds]) for inds in partition(randperm(size(y_train, 2)), batchsize)]

#+

using Random
using BSON

Random.seed!(666)
m = Chain(
    Conv((2,2), 1=>16, relu),
    MaxPool((2,2)),
    Conv((2,2), 16=>8, relu),
    MaxPool((2,2)),
    flatten,
    Dense(288, size(y_train,1)),
    softmax,
)

#+

using Flux: crossentropy

L(X, y) = crossentropy(m(X), y)

#+

function train_model!(m, L, X, y;
        opt = Descent(0.1),
        batchsize = 128,
        n_epochs = 10,
        file_name = "")

    batches = DataLoader((X, y); batchsize, shuffle = true)

    for _ in 1:n_epochs
        Flux.train!(L, params(m), batches, opt)
    end

    !isempty(file_name) && BSON.bson(file_name, m=m)

    return
end

# ### Exercise:
# Train the model for one epoch and save it to `MNIST_simple.bson`. Print the accuracy on
# the testing set.
#
# ---
# ### Solution:

# ---
#
# ### Exercise:
# Write a function `train_or_load!(file_name, m, args...; ???)` checking whether the file
# `file_name` exists.
# - If it exists, it loads it and then copies its parameters into `m` using the function
# `Flux.loadparams!`.
# - If it does not exist, it trains it using `train_model!`.
# In both cases, the model `m` should be modified inside the `train_or_load!` function. Pay
# special attention to the optional arguments `???`.
#
# Use this function to load the model from `data/mnist.bson` and evaluate the performance
# at the testing set.
#
# ---
# ### Solution:


text:
@testset "Blades" begin
    @test blade.index == 4
    @test eltype(blade) == Int
    @test eltype(Blade(1.0, 4)) == Float64
    @test Blade(1.0, 4) == blade
    @test Blade(1.0, 5) ≉ blade
    @test Blade(1.0, 4) ≈ blade
    @test Blade(1.0, 5) ≉ blade
    @test is_homogeneous(blade)
    @test scalar(1.0) == Blade(1.0, 1)
end


text:
export SpUri, SpId, SpCategoryId, SpUserId, SpUrl
"""
All web API arguments are strings, but types
`SpUri`, `SpId`, `CategoryId`, `SpUserId`, `SpUrl`
aids in picking default values. [format](
https://developer.spotify.com/documentation/web-api/#spotify-uris-and-ids)

PARAMETER  | DESCRIPTION | VALUE
:----------:| :---------- |:------------:
SpUri       |The resource identifier that you can enter, for example, in the |
            | Spotify Desktop client’s search box to locate an artist,|
            | album, or track. To find a Spotify URI simply right-|
            | click (on Windows) or Ctrl-Click (on a Mac) on |spotify:track:
            | the artist’s or album’s or track’s name.| 6rqhFgbbKwnb9MLmUQDhG6|
SpId        |The base-62 identifier that you can find at the end of the |
            | Spotify URI (see above) for an artist, track, album, |
            | playlist, etc. Unlike a Spotify URI, a Spotify ID does|
            | not clearly identify the type of resource; that information is |
            | provided elsewhere in the call.| 6rqhFgbbKwnb9MLmUQDhG6
SpCategoryId|The unique string identifying the Spotify category.| party|
SpUserId    |The unique string identifying the Spotify user that you can |
            | find at the end of the Spotify URI for the user. The ID |
            | of the current user can be obtained via the Web API endpoint.| wizzler
SpUrl       |An HTML link that opens a track, album, app, playlist or other |
            | Spotify resource in a Spotify client (which  client |
            | is determined by the user’s device and |  http://open.spotify.com/
            | account settings at play.spotify.com. |   track/6rqhFgbbKwnb9MLmUQDhG6

"""
SpUri, SpId, SpCategoryId, SpUserId, SpUrl

mutable struct SpUri
    s::String
    SpUri(s) = isuri(s) ? new(s) : error("must be 'spotify:track:<base 62 string>")
end
SpUri() = SpUri("spotify:track:6rqhFgbbKwnb9MLmUQDhG6")
mutable struct SpId
    s::String
    SpId(s) =  isid(s) ? new(s) : error("must be <base 62 string>")
end
SpId() = SpId("6rqhFgbbKwnb9MLmUQDhG6")
mutable struct SpCategoryId
    s::String
end
SpCategoryId() = SpCategoryId("party")
mutable struct SpUserId
    s::String
end
SpUserId() = get_user_name() == "" ? throw(".ini file has no user name") : SpUserId(get_user_name())

mutable struct SpUrl
    s::String
    SpUrl(s) = isurl(s) ? new(s) : error("must be an url")
end
SpUrl() =  SpUrl("http://open.spotify.com/track/6rqhFgbbKwnb9MLmUQDhG6")
isurl(s) = !isnothing(match(r"(http(s)?:\/\/.)?(www\.)?[-a-zA-Z0-9@:%._\+~#=]{2,256}\.[a-z]{2,6}\b([-a-zA-Z0-9@:%_\+.~#?&//=]*)", s))
isid(s) = !isnothing(match(r"\b[a-zA-Z0-9]{22}", s))
isuri(s) = count(==( ':'), s) == 2 && length(s) - findlast(':', s) == 22
# These types should work with string interpolation using $.
# So we need to alter their undecorated representations.
# For argument inspection, we don't want to hide them not
# being ordinary strings, so we add color. The color will
# not be used in contexts which can't display color (i.e.
# when used as an url)
show(io::IO, s::SpUri) =        printstyled(io, s.s; color = :blue)
show(io::IO, s::SpId) =         printstyled(io, s.s; color = 176)
show(io::IO, s::SpCategoryId) = printstyled(io, s.s; color = 172)
show(io::IO, s::SpUrl) =        printstyled(io, s.s; color = :blue, bold = true)
show(io::IO, s::SpUserId) =     printstyled(io, s.s; color = :green)
#
show(io::IO,  ::MIME"text/plain", s::SpUri) =       printstyled(io, '"' * s.s * '"'; color = :blue)
show(io::IO, ::MIME"text/plain", s::SpId) =         printstyled(io, '"' * s.s * '"'; color = 176)
show(io::IO, ::MIME"text/plain", s::SpCategoryId) = printstyled(io, '"' * s.s * '"'; color = 172)
show(io::IO, ::MIME"text/plain", s::SpUrl) =        printstyled(io, '"' * s.s * '"'; color = :blue, bold = true)
show(io::IO, ::MIME"text/plain", s::SpUserId) =     printstyled(io, '"' * s.s * '"'; color = :green)

show(io::IO, v::Vector{SpUri}) =        show_vector(io, v, "", "")
show(io::IO, v::Vector{SpId}) =         show_vector(io, v, "", "")
show(io::IO, v::Vector{SpCategoryId}) = show_vector(io, v, "", "")
show(io::IO, v::Vector{SpUrl}) =         show_vector(io, v, "", "")
show(io::IO, v::Vector{SpUserId}) =         show_vector(io, v, "", "")

typeinfo_implicit(::Type{SpUri}) = true
typeinfo_implicit(::Type{SpId}) = true
typeinfo_implicit(::Type{SpCategoryId}) = true
typeinfo_implicit(::Type{SpUrl}) = true
typeinfo_implicit(::Type{SpUserId}) = true


text:
##### Buoyancy models

@inline buoyancy(param_set, α_0::FT, α::FT) where {FT} =
    FT(grav(param_set)) * (α - α_0) / α_0

"""
    compute_buoyancy!

Define buoyancy field

 - `aux[:buoy, k, i]`

 for all `k` and all `i`
"""
function compute_buoyancy! end

function compute_buoyancy!(grid, q, aux, params)
    gm, en, ud, sd, al = allcombinations(q)
    @unpack param_set = params
    @inbounds for i in (ud..., en)
        @inbounds for k in over_elems_real(grid)
            q_tot = q[:q_tot, k, i]
            θ_liq = q[:θ_liq, k, i]
            p_0 = aux[:p_0, k]
            ts = PhaseEquil_pθq(param_set, p_0, θ_liq, q_tot)
            α_i = specific_volume(ts)
            aux[:buoy, k, i] = buoyancy(param_set, aux[:α_0, k], α_i)
        end
    end

    # Filter buoyancy
    @inbounds for i in ud
        @inbounds for k in over_elems_real(grid)
            weight = aux[:HVSD_a, k, i]
            aux[:buoy, k, i] =
                weight * aux[:buoy, k, i] + (1 - weight) * aux[:buoy, k, en]
        end
    end

    # Subtract grid mean buoyancy
    @inbounds for k in over_elems_real(grid)
        aux[:buoy, k, gm] = sum([q[:a, k, i] * aux[:buoy, k, i] for i in sd])
        @inbounds for i in sd
            aux[:buoy, k, i] -= aux[:buoy, k, gm]
        end
    end
end

function compute_tke_buoy!(
    grid::Grid{FT},
    q,
    aux,
    aux_O2,
    cv,
    params,
) where {FT}
    gm, en, ud, sd, al = allcombinations(q)
    @unpack param_set = params

    _R_v::FT = R_v(param_set)
    _R_d::FT = R_d(param_set)
    _grav::FT = grav(param_set)

    # Note that source terms at the first interior point are not really used because that is where tke boundary condition is
    # enforced (according to MO similarity). Thus here I am being sloppy about lowest grid point
    @inbounds for k in over_elems_real(grid)
        q_tot_dry = aux[:q_tot_dry, k]
        θ_dry = aux[:θ_dry, k]
        t_cloudy = aux[:t_cloudy, k]
        q_vap_cloudy = aux[:q_vap_cloudy, k]
        q_tot_cloudy = aux[:q_tot_cloudy, k]
        θ_cloudy = aux[:θ_cloudy, k]
        p_0 = aux[:p_0, k]

        lh = latent_heat_vapor(param_set, t_cloudy)
        cpm = cp_m(param_set, PhasePartition(q_tot_cloudy))
        grad_θ_liq = ∇_dn(q[:θ_liq, Cut(k), en], grid)
        grad_q_tot = ∇_dn(q[:q_tot, Cut(k), en], grid)

        prefactor = _R_d * TD.exner_given_pressure(param_set, p_0) / p_0
        ε_vi = _R_v / _R_d

        d_alpha_θ_liq_dry = prefactor * (1 + (ε_vi - 1) * q_tot_dry)
        d_alpha_q_tot_dry = prefactor * θ_dry * (ε_vi - 1)
        CF_env = aux[:CF, k]

        if CF_env > 0
            d_alpha_θ_liq_cloudy = (
                prefactor * (
                    1 + ε_vi * (1 + lh / _R_v / t_cloudy) * q_vap_cloudy -
                    q_tot_cloudy
                ) / (
                    1 +
                    lh * lh / cpm / _R_v / t_cloudy / t_cloudy * q_vap_cloudy
                )
            )
            d_alpha_q_tot_cloudy =
                (lh / cpm / t_cloudy * d_alpha_θ_liq_cloudy - prefactor) *
                θ_cloudy
        else
            d_alpha_θ_liq_cloudy = 0
            d_alpha_q_tot_cloudy = 0
        end

        d_alpha_θ_liq_total =
            (CF_env * d_alpha_θ_liq_cloudy + (1 - CF_env) * d_alpha_θ_liq_dry)
        d_alpha_q_tot_total =
            (CF_env * d_alpha_q_tot_cloudy + (1 - CF_env) * d_alpha_q_tot_dry)

        K_h_k = aux[:K_h, k, gm]
        term_1 = -K_h_k * grad_θ_liq * d_alpha_θ_liq_total
        term_2 = -K_h_k * grad_q_tot * d_alpha_q_tot_total

        # TODO - check
        aux_O2[cv][:buoy, k] =
            _grav / aux[:α_0, k] *
            q[:a, k, en] *
            aux[:ρ_0, k] *
            (term_1 + term_2)
    end
end


text:
using Test
using Printf: @printf, @sprintf
using StaticArrays: SVector
using QuantumLattices.Essentials.DegreesOfFreedom
using QuantumLattices.Essentials.QuantumOperators: ID, Operator, Operators, sequence, id, LaTeX, latexformat
using QuantumLattices.Essentials.Spatials: AbstractPID, PID, CPID, Point, Bond, Bonds, Lattice, pidtype, rcoord, icoord
using QuantumLattices.Essentials: kind, update!, reset!
using QuantumLattices.Interfaces: rank, expand, ⊗, ⊕
using QuantumLattices.Prerequisites: Float, decimaltostr
using QuantumLattices.Prerequisites.Traits: parameternames, isparameterbound, contentnames, getcontent

import QuantumLattices.Essentials.QuantumOperators: latexname, script
import QuantumLattices.Essentials.DegreesOfFreedom: ishermitian, statistics, couplingcenters, abbr
import QuantumLattices.Prerequisites.VectorSpaces: shape

struct DID{N<:Union{Int, Symbol}} <: SimpleIID
    nambu::N
end
@inline Base.show(io::IO, did::DID) = @printf io "DID(%s)" did.nambu
@inline Base.adjoint(sl::DID{Int}) = DID(3-sl.nambu)
@inline statistics(::Type{<:DID}) = :f
@inline DID(did::DID, ::CompositeInternal) = did

struct DFock <: SimpleInternal{DID{Int}}
    nnambu::Int
end
@inline shape(f::DFock) = (1:f.nnambu,)
@inline DID(i::CartesianIndex, vs::DFock) = DID(i.I...)
@inline CartesianIndex(did::DID{Int}, vs::DFock) = CartesianIndex(did.nambu)
@inline shape(iidspace::IIDSpace{DID{Symbol}, DFock}) = (1:iidspace.internal.nnambu,)
@inline shape(iidspace::IIDSpace{DID{Int}, DFock}) = (iidspace.iid.nambu:iidspace.iid.nambu,)
@inline shape(iidspace::IIDSpace{DID{Symbol}, DFock}) = (1:iidspace.internal.nnambu,)
@inline shape(iidspace::IIDSpace{DID{Int}, DFock}) = (iidspace.iid.nambu:iidspace.iid.nambu,)

@inline script(::Val{:site}, index::Index{<:AbstractPID, <:DID}; kwargs...) = index.pid.site
@inline script(::Val{:nambu}, index::Index{<:AbstractPID, <:DID}; kwargs...) = index.iid.nambu==2 ? "\\dagger" : ""

@inline latexname(::Type{<:OID{<:Index{<:AbstractPID, <:DID}}}) = Symbol("OID{Index{AbstractPID, DID}}")
latexformat(OID{<:Index{<:AbstractPID, <:DID}}, LaTeX{(:nambu,), (:site,)}('d'))

const DCoupling{V, I<:ID{DID}, C<:Subscripts} = Coupling{V, I, C}
@inline Base.repr(tc::(Coupling{V, <:ID{DID}} where V)) = @sprintf "%s ph(%s)" decimaltostr(tc.value) join(tc.cid.nambus, ", ")
@inline couplingcenters(::(Coupling{V, <:ID{DID}} where V), ::Bond, ::Val) = (1, 2)
@inline DCoupling(value, nambus::Tuple{Vararg{Int}}) = Coupling(value, ID(DID, nambus), Subscripts((nambu=Subscript(nambus),)))
@inline DCoupling(value, nambus::Subscript) = Coupling(value, ID(DID, Tuple(nambus)), Subscripts((nambu=nambus,)))

@inline abbr(::Type{<:Term{:Mu}}) = :mu
@inline abbr(::Type{<:Term{:Hp}}) = :hp
@inline ishermitian(::Type{<:Term{:Mu}}) = true
@inline ishermitian(::Type{<:Term{:Hp}}) = false

@testset "IID" begin
    did = DID(1)
    @test statistics(did) == statistics(typeof(did)) == :f

    did₁, did₂ = DID(1), DID(2)
    ciid = CompositeIID(did₁, did₂)
    @test length(ciid) == length(typeof(ciid)) == 2
    @test rank(ciid) == rank(typeof(ciid)) == 2
    @test iidtype(ciid, 1) == iidtype(typeof(ciid), 1) == DID{Int}
    @test iidtype(ciid, 2) == iidtype(typeof(ciid), 2) == DID{Int}
    @test ciid[1]==did₁ && ciid[2]==did₂
    @test ciid.nambus == (1, 2)
    @test string(ciid) == "DID(1) ⊗ DID(2)"
    @test did₁⊗did₂ == ciid
    @test did₁⊗ciid == CompositeIID(did₁, did₁, did₂)
    @test ciid⊗did₁ == CompositeIID(did₁, did₂, did₁)
    @test ciid⊗ciid == CompositeIID(did₁, did₂, did₁, did₂)
end

@testset "SimpleInternal" begin
    it = DFock(2)
    @test it|>eltype == DID{Int}
    @test it|>typeof|>eltype == DID{Int}
    @test it == deepcopy(it)
    @test isequal(it, deepcopy(it))
    @test it|>string == "DFock(nnambu=2)"
    @test it|>collect == [DID(1), DID(2)]
    @test statistics(it) == statistics(typeof(it)) == :f
    @test match(DID(1), it) && match(DID, DFock)
    @test filter(DID(1), it) == filter(DID, it) == it
    @test filter(DID(1), DFock) == filter(DID, DFock) == DFock
end

@testset "CompositeInternal" begin
    it₁, it₂ = DFock(2), DFock(3)

    ci = CompositeInternal{:⊕}(it₁, it₂)
    @test eltype(ci) == eltype(typeof(ci)) == DID{Int}
    @test rank(ci) == rank(typeof(ci)) == 2
    @test string(ci) == "DFock(nnambu=2) ⊕ DFock(nnambu=3)"
    for i = 1:2
        @test ci[i] == it₁[i]
    end
    for i = 1:3
        @test ci[2+i] == it₂[i]
    end
    @test it₁⊕it₂ == ci
    @test it₁⊕ci == CompositeInternal{:⊕}(it₁, it₁, it₂)
    @test ci⊕it₁ == CompositeInternal{:⊕}(it₁, it₂, it₁)
    @test ci⊕ci == CompositeInternal{:⊕}(it₁, it₂, it₁, it₂)
    @test filter(DID(1), ci) == filter(DID, ci) == ci
    @test filter(DID(1), typeof(ci)) == filter(DID, typeof(ci)) == typeof(ci)

    ci = CompositeInternal{:⊗}(it₁, it₂)
    @test eltype(ci) == eltype(typeof(ci)) == CompositeIID{Tuple{DID{Int}, DID{Int}}}
    @test rank(ci) == rank(typeof(ci)) == 2
    @test string(ci) == "DFock(nnambu=2) ⊗ DFock(nnambu=3)"
    count = 1
    for i = 1:3
        for j = 1:2
            @test ci[count] == CompositeIID(it₁[j], it₂[i])
            count += 1
        end
    end
    @test it₁⊗it₂ == ci
    @test it₁⊗ci == CompositeInternal{:⊗}(it₁, it₁, it₂)
    @test ci⊗it₁ == CompositeInternal{:⊗}(it₁, it₂, it₁)
    @test ci⊗ci == CompositeInternal{:⊗}(it₁, it₂, it₁, it₂)
    @test filter(DID(1), ci) == filter(DID, ci) == ci
    @test filter(DID(1), typeof(ci)) == filter(DID, typeof(ci)) == typeof(ci)
end

@testset "Index" begin
    @test parameternames(Index) == (:pid, :iid)
    @test isparameterbound(Index, Val(:pid), PID) == false
    @test isparameterbound(Index, Val(:iid), DID) == true

    index = Index(CPID('S', 4), DID(1))
    @test index|>pidtype == CPID{Char}
    @test index|>typeof|>pidtype == CPID{Char}
    @test index|>iidtype == DID{Int}
    @test index|>typeof|>iidtype == DID{Int}
    @test index|>adjoint == Index(CPID('S', 4), DID(2))
    @test statistics(index) == statistics(typeof(index)) == :f
    @test ishermitian(ID(index', index)) == true
    @test ishermitian(ID(index, index)) == false
end

@testset "OID" begin
    @test contentnames(CompositeOID) == (:index,)
    @test parameternames(CompositeOID) == (:index,)
    @test isparameterbound(CompositeOID, Val(:index), Index) == true
    @test isparameterbound(CompositeOID, Val(:index), Index{PID, DID{Int}}) == false

    @test contentnames(OID) == (:index, :rcoord, :icoord)
    @test parameternames(OID) == (:index, :coord)
    @test isparameterbound(OID, Val(:index), Index) == true
    @test isparameterbound(OID, Val(:index), Index{PID, DID{Int}}) == false
    @test isparameterbound(OID, Val(:coord), SVector) == true
    @test isparameterbound(OID, Val(:coord), SVector{2, Float}) == false

    oid = OID(Index(PID(1), DID(1)), [0.0, -0.0], [0.0, 0.0])
    @test indextype(oid) == indextype(typeof(oid)) == Index{PID, DID{Int}}
    @test statistics(oid) == statistics(typeof(oid)) == :f
    @test oid' == OID(Index(PID(1), DID(2)), rcoord=SVector(0.0, 0.0), icoord=SVector(0.0, 0.0))
    @test hash(oid, UInt(1)) == hash(OID(Index(PID(1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 1.0)), UInt(1))
    @test propertynames(ID(oid)) == (:indexes, :rcoords, :icoords)
    @test fieldnames(OID) == (:index, :rcoord, :icoord)
    @test string(oid) == "OID(Index(PID(1), DID(1)), [0.0, 0.0], [0.0, 0.0])"
    @test ID(oid', oid)' == ID(oid', oid)
    @test ishermitian(ID(oid', oid)) == true
    @test ishermitian(ID(oid, oid)) == false
    @test oidtype(DFock, Point{2, PID, Float}, Val(:default)) == OID{Index{PID, DID{Int}}, SVector{2, Float}}
end

@testset "Operator" begin
    opt = Operator(1.0,
        OID(Index(PID(1), DID(2)), SVector(0.5, 0.5), SVector(1.0, 1.0)),
        OID(Index(PID(1), DID(1)), SVector(0.0, 0.5), SVector(0.0, 1.0))
        )
    @test rcoord(opt) == SVector(0.5, 0.0)
    @test icoord(opt) == SVector(1.0, 0.0)

    opt = Operator(1.0, OID(Index(PID(1), DID(2)), SVector(0.5, 0.0), SVector(1.0, 0.0)))
    @test rcoord(opt) == SVector(0.5, 0.0)
    @test icoord(opt) == SVector(1.0, 0.0)
end

@testset "Hilbert" begin
    map = pid->DFock((pid.site-1)%2+1)
    hilbert = Hilbert(map, [CPID(1, 1), CPID(1, 2)])
    @test hilbert == Hilbert{DFock}(map, [CPID(1, 1), CPID(1, 2)])
    reset!(hilbert, (CPID(2, 1), CPID(2, 2)))

    hilbert = Hilbert(pid=>DFock(2) for pid in [PID(1), PID(2)])
    @test hilbert[PID(1)] == hilbert[PID(2)] == DFock(2)

    hilbert = Hilbert(PID(1)=>DFock(2), PID(2)=>DFock(3))
    @test hilbert[PID(1)]==DFock(2) && hilbert[PID(2)]==DFock(3)
end

@testset "IIDSpace" begin
    DID₁, DID₂, it = DID(2), DID(:σ), DFock(2)
    iidspace = IIDSpace(DID₁⊗DID₂, it⊗it)
    @test eltype(iidspace) == eltype(typeof(iidspace)) == CompositeIID{Tuple{DID{Int}, DID{Int}}}
    @test kind(iidspace) == kind(typeof(iidspace)) == :info
    @test length(iidspace) == 2
    @test iidspace[1]==DID(2)⊗DID(1) && iidspace[2]==DID(2)⊗DID(2)
    @test expand((DID₁, DID₂), (it, it)) == iidspace
    @test collect(iidspace) == [DID(2)⊗DID(1), DID(2)⊗DID(2)]
end

@testset "Metric" begin
    valtype(Metric, AbstractOID) = AbstractOID

    index = Index(CPID('S', 4), DID(1))
    oid = OID(index, SVector(0.5, 0.0), SVector(1.0, 0.0))

    m = OIDToTuple((:scope, :site, :nambu))
    @test m == OIDToTuple(:scope, :site, :nambu)
    @test isequal(m, OIDToTuple(:scope, :site, :nambu))
    @test keys(m) == keys(typeof(m)) == (:scope, :site, :nambu)
    @test filter(≠(:nambu), m) == OIDToTuple(:scope, :site)
    @test OIDToTuple(Index{PID, DID{Int}}) == OIDToTuple(:site, :nambu)
    @test OIDToTuple(OID{Index{CPID{Int}, DID{Int}}}) == OIDToTuple(:scope, :site, :nambu)
    @test OIDToTuple(Hilbert{DFock, CPID{Int}}) == OIDToTuple(:scope, :site, :nambu)
    @test m(index) == ('S', 4, 1) == m(oid)
    @test filter(≠(:scope), m)(index) == (4, 1)
    @test filter(≠(:nambu), m)(index) == ('S', 4)
    @test filter(∉((:site, :nambu)), m)(index) == ('S',)
end

@testset "Table" begin
    @test contentnames(Table) == (:by, :contents)

    by = filter(≠(:nambu), OIDToTuple(Index{PID, DID{Int}}))

    table = Table([Index(PID(1), DID(1)), Index(PID(1), DID(2))], by)
    @test empty(table) == Table{Index{PID, DID{Int}}}(by)
    @test table[Index(PID(1), DID(1))]==1 && table[Index(PID(1), DID(2))]==1

    hilbert = Hilbert(pid=>DFock(2) for pid in [PID(1), PID(2)])
    inds1 = (Index(PID(1), iid) for iid in DFock(2))|>collect
    inds2 = (Index(PID(2), iid) for iid in DFock(2))|>collect
    @test Table(hilbert, by) == Table([inds1; inds2], by)
    @test Table(hilbert, by) == union(Table(inds1, by), Table(inds2, by))

    opt = Operator(1.0im,
        OID(Index(PID(1), DID(2)), SVector(0.0, 0.0), SVector(1.0, 0.0)),
        OID(Index(PID(1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0))
        )
    @test sequence(opt, table) == (1, 1)
    @test haskey(table, opt.id) == (true, true)

    table = Table(hilbert)
    @test table == Table([inds1; inds2])
    @test reset!(empty(table), [inds1; inds2]) == table
    @test reset!(empty(table), hilbert) == table
end

@testset "Subscript" begin
    sub = Subscript(4)
    @test contentnames(typeof(sub)) == (:contents, :rep, :constraint)
    @test getcontent(sub, :contents) == sub.pattern
    @test sub==deepcopy(sub) && isequal(sub, deepcopy(sub))
    @test string(sub) == "[* * * *]"
    @test rank(sub) == rank(typeof(sub)) == 4
    @test match(sub, (2, 2, 2, 2))

    sub = Subscript((1, 2, 3, 4))
    @test string(sub) == "[1 2 3 4]"
    @test sub == subscript"[1 2 3 4]" == subscript"[1 2 3 4]; false"
    @test rank(sub) == 4
    @test match(sub, (1, 3, 3, 4))

    sub = Subscript((1, 2, 2, 1), true)
    @test string(sub) == "[1 2 2 1]"
    @test sub == subscript"[1 2 2 1]; true"
    @test rank(sub) == 4
    @test match(sub, (1, 2, 2, 1)) && !match(sub, (1, 1, 2, 1))

    sub = subscript"[x₁ x₂ x₁ x₂]"
    @test string(sub) == "[x₁ x₂ x₁ x₂]"
    @test rank(sub) == 4
    @test match(sub, (1, 2, 1, 2)) && !match(sub, (1, 2, 2, 1))

    sub = subscript"[x₁ 4 4 x₂](x₁ < x₂)"
    @test string(sub) == "[x₁ 4 4 x₂](x₁ < x₂)"
    @test rank(sub) == rank(typeof(sub)) == 4
    @test match(sub, (1, 4, 4, 2)) && !match(sub, (2, 4, 4, 1))
end

@testset "Subscripts" begin
    subscripts = Subscripts((nambu=subscript"[x y](x > y)",), (nambu=subscript"[x x y y](x < y)",))
    @test subscripts.rep == (1:2=>("[x y](x > y)",), 3:6=>("[x x y y](x < y)",))
    @test hash(subscripts) == hash(subscripts.rep)
    @test string(subscripts) == "nambu[x y](x > y) × nambu[x x y y](x < y)"
    @test repr(subscripts, 1:2, :nambu) == "[x y](x > y)×[x x y y](x < y)"
    @test repr(subscripts, 1, :nambu) == "[x y](x > y)"
    @test repr(subscripts, 2, :nambu) == "[x x y y](x < y)"

    @test rank(subscripts) == rank(typeof(subscripts)) == 6
    @test rank(subscripts, 1) == rank(typeof(subscripts), 1) == 2
    @test rank(subscripts, 2) == rank(typeof(subscripts), 2) == 4
    @test match(subscripts, (DID(2), DID(1), DID(2), DID(2), DID(3), DID(3)))
    @test match(subscripts, DID(2)⊗DID(1)⊗DID(2)⊗DID(2)⊗DID(3)⊗DID(3))
    @test !match(subscripts, DID(1)⊗DID(2)⊗DID(2)⊗DID(2)⊗DID(3)⊗DID(3))
    @test !match(subscripts, DID(2)⊗DID(1)⊗DID(3)⊗DID(3)⊗DID(2)⊗DID(2))
    @test subscripts == Subscripts((nambu=subscript"[x y](x > y)",))*Subscripts((nambu=subscript"[x x y y](x < y)",))
end

@testset "couplings" begin
    @test parameternames(Coupling) == (:value, :cid, :subscripts)
    @test isparameterbound(Coupling, :cid, Tuple{DID{Int}}) == false
    @test isparameterbound(Coupling, :cid, ID{DID{Int}}) == true
    @test isparameterbound(Coupling, :subscripts, Subscripts{Tuple{NamedTuple{(:nambu,), Tuple{Subscript{Tuple{Int}}}}}}) == true
    @test isparameterbound(Coupling, :subscripts, Subscripts) == true
    @test isparameterbound(Coupling, :subscripts, Subscripts{Tuple{}, Tuple{}}) == false

    tc = DCoupling(2.0, (2,))
    @test id(tc) == ID(CompositeIID(tc.cid), tc.subscripts)
    @test rank(tc) == rank(typeof(tc)) == 1
    @test tc == Coupling(2.0, id(tc))
    @test ID{SimpleIID}(tc) == tc.cid
    @test Subscripts(tc) == tc.subscripts

    point = Point(CPID(1, 1), (0.0, 0.0), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>DFock(2))
    @test couplingcenters(tc, point, Val(:Mu)) == (1,)
    @test couplingpoints(tc, point, Val(:Mu)) == (point,)
    @test couplinginternals(tc, point, hilbert, Val(:Mu)) == (DFock(2),)

    tc₁ = DCoupling(1.5, (1, 2))
    tc₂ = DCoupling(2.0, subscript"[a b](a < b)")
    @test tc₁*tc₂ == Coupling(3.0, ID(DID(1), DID(2), DID(:a), DID(:b)), Subscripts((nambu=Subscript((1, 2)),), (nambu=subscript"[a b](a < b)",)))

    ex = expand(tc₁, point, hilbert, Val(:info))
    @test eltype(ex) == eltype(typeof(ex)) == Operator{Float64, NTuple{2, OID{Index{CPID{Int}, DID{Int}}, SVector{2, Float64}}}}
    @test collect(ex) == [
        Operator(1.5, ID(
            OID(Index(CPID(1, 1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0)),
            OID(Index(CPID(1, 1), DID(2)), SVector(0.0, 0.0), SVector(0.0, 0.0))
            ))
        ]

    ex = expand(tc₂, point, hilbert, Val(:info))
    @test eltype(ex) == eltype(typeof(ex)) == Operator{Float64, NTuple{2, OID{Index{CPID{Int}, DID{Int}}, SVector{2, Float64}}}}
    @test collect(ex) == [
        Operator(2.0, ID(
            OID(Index(CPID(1, 1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0)),
            OID(Index(CPID(1, 1), DID(2)), SVector(0.0, 0.0), SVector(0.0, 0.0))
            ))
        ]

    @test @couplings(DCoupling(1.0, (1, 1))) == Couplings(DCoupling(1.0, (1, 1)))
    @test @couplings(DCoupling(1.0, (1,))+DCoupling(1.0, (2,))) == DCoupling(1.0, (1,))+DCoupling(1.0, (2,))
end

@testset "TermFunction" begin
    ta = TermAmplitude()
    @test ta() == 1

    ta = TermAmplitude(x->x+3.0)
    @test ta(1.0) == 4.0

    tcs = DCoupling(1.0, (1, 1)) + DCoupling(2.0, (2, 2))
    termcouplings = TermCouplings(tcs)
    @test termcouplings == deepcopy(TermCouplings(tcs))
    @test isequal(termcouplings, deepcopy(TermCouplings(tcs)))
    @test valtype(termcouplings) == valtype(typeof(termcouplings)) == typeof(tcs)
    @test termcouplings() == tcs

    fx = i -> i%2==1 ? (DCoupling(1.0, (1, 1))+DCoupling(1.0, (2, 2))) : (DCoupling(1.0, (2, 1))+DCoupling(1.0, (1, 2)))
    termcouplings = TermCouplings(fx)
    @test termcouplings == TermCouplings{typejoin(typeof(fx(1)), typeof(fx(2)))}(fx)
    @test valtype(termcouplings) == valtype(typeof(termcouplings)) == typejoin(typeof(fx(1)), typeof(fx(2)))
    @test termcouplings(1) == fx(1)
    @test termcouplings(2) == fx(2)

    @test ismodulatable(TermModulate{Val{true}}) == ismodulatable(TermModulate{<:Function}) == true
    @test ismodulatable(TermModulate{Val{false}}) == false
    termmodulate = TermModulate(:t)
    @test termmodulate(t=1) == 1
    @test isnothing(termmodulate(mu=1))
    @test isnothing(termmodulate())
    @test ismodulatable(termmodulate) == true

    termmodulate = TermModulate(:t, t->t*2.0)
    @test ismodulatable(termmodulate) == true
    @test termmodulate(2) == 4
end

@testset "Term" begin
    point = Point(CPID(1, 1), (0.0, 0.0), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>DFock(2))
    term = Term{:Mu}(:mu, 1.5, 0, couplings=@couplings(DCoupling(1.0, (2, 1))), amplitude=(bond->3), modulate=false)
    @test term|>kind == term|>typeof|>kind == :Mu
    @test term|>id == term|>typeof|>id == :mu
    @test term|>valtype == term|>typeof|>valtype == Float
    @test term|>rank == term|>typeof|>rank == 2
    @test term|>abbr == term|>typeof|>abbr == :mu
    @test term|>ismodulatable == term|>typeof|>ismodulatable == false
    @test term|>ishermitian == term|>typeof|>ishermitian == true
    @test term == deepcopy(term)
    @test isequal(term, deepcopy(term))
    @test string(term) == "Mu{2}(id=mu, value=1.5, bondkind=0, factor=1.0)"
    @test repr(term, point, hilbert) == "mu: 4.5 ph(2, 1)"
    @test +term == term
    @test -term == term*(-1) == replace(term, factor=-term.factor)
    @test 2*term == term*2 == replace(term, factor=term.factor*2)
    @test term/2 == term*(1/2) == replace(term, factor=term.factor/2)

    p1 = Point(PID(1), (0.0, 0.0), (0.0, 0.0))
    p2 = Point(PID(2), (0.0, 0.0), (0.0, 0.0))
    hilbert = Hilbert(pid=>DFock(2) for pid in [PID(1), PID(2)])
    term = Term{:Mu}(:mu, 1.5, 0,
        couplings=bond->bond.pid.site%2==0 ? Couplings(DCoupling(1.0, (2, 2))) : Couplings(DCoupling(1.0, (1, 1))),
        amplitude=bond->3,
        modulate=true
    )
    @test term|>ismodulatable == term|>typeof|>ismodulatable == true
    @test repr(term, p1, hilbert) == "mu: 4.5 ph(1, 1)"
    @test repr(term, p2, hilbert) == "mu: 4.5 ph(2, 2)"
    @test one(term) == replace(term, value=1.0)
    @test zero(term) == replace(term, value=0.0)
    @test term.modulate(mu=4.0) == 4.0
    @test isnothing(term.modulate(t=1.0))
    @test update!(term, mu=4.25) == replace(term, value=4.25)
    @test term.value == 4.25
end

@testset "expand" begin
    point = Point(PID(1), (0.0, 0.0), (0.0, 0.0))
    hilbert = Hilbert(point.pid=>DFock(2))
    term = Term{:Mu}(:mu, 1.5, 0, couplings=@couplings(DCoupling(1.0, (2, 1))), amplitude=bond->3.0, modulate=true)
    operators = Operators(
        Operator(+2.25,
            OID(Index(PID(1), DID(2)), SVector(0.0, 0.0), SVector(0.0, 0.0)),
            OID(Index(PID(1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0))
            )
        )
    @test expand(term, point, hilbert, half=true) == operators
    @test expand(term, point, hilbert, half=false) == operators*2

    bond = Bond(1, Point(PID(2), (1.5, 1.5), (1.0, 1.0)), Point(PID(1), (0.5, 0.5), (0.0, 0.0)))
    hilbert = Hilbert(pid=>DFock(2) for pid in [PID(1), PID(2)])
    term = Term{:Hp}(:t, 1.5, 1, couplings=@couplings(DCoupling(1.0, (2, 1))), amplitude=bond->3.0, modulate=true)
    operators = Operators(
        Operator(4.5,
            OID(Index(PID(1), DID(2)), SVector(0.5, 0.5), SVector(0.0, 0.0)),
            OID(Index(PID(2), DID(1)), SVector(1.5, 1.5), SVector(1.0, 1.0))
            )
        )
    @test expand(term, bond, hilbert, half=true) == operators
    @test expand(term, bond, hilbert, half=false) == operators+operators'

    lattice = Lattice(:Tuanzi, [Point(CPID(1, 1), (0.0, 0.0), (0.0, 0.0))], vectors=[[1.0, 0.0]], neighbors=1)
    bonds = Bonds(lattice)
    hilbert = Hilbert(pid=>DFock(2) for pid in lattice.pids)
    term = Term{:Mu}(:mu, 1.5, 0, couplings=@couplings(DCoupling(1.0, (2, 1))), amplitude=bond->3.0, modulate=true)
    operators = Operators(
        Operator(+2.25,
            OID(Index(CPID(1, 1), DID(2)), SVector(0.0, 0.0), SVector(0.0, 0.0)),
            OID(Index(CPID(1, 1), DID(1)), SVector(0.0, 0.0), SVector(0.0, 0.0))
            )
        )
    @test expand(term, bonds, hilbert, half=true) == operators
    @test expand(term, bonds, hilbert, half=false) == operators*2
end

@testset "script" begin
    oid = OID(Index(CPID('d', 1), DID(2)), SVector(0.0, 0.0), SVector(1.0, 0.0))
    latex = LaTeX{(:nambu,), (:site,)}('c', vectors=(SVector(1.0, 0.0), SVector(0.0, 1.0)))
    @test script(Val(:rcoord), oid) == "[0.0, 0.0]"
    @test script(Val(:icoord), oid) == "[1.0, 0.0]"
    @test script(Val(:integralicoord), oid; vectors=get(latex.options, :vectors, nothing)) == "[1, 0]"
    @test script(Val(:site), oid) == 1
    @test script(Val(:nambu), oid) == "\\dagger"
end


text:
include("header.jl")
using Printf

# x* = f(1, · · · , 1)
# f(x*) = 0

rosenbrock(x) = sum(abs2.(1 .- x[1:end-1]) .+ 100 .* abs2.(x[2:end] .- abs2.(x[1:end-1])))

function rosenmulti(x)
    v = value(x)
    if isbitstype(eltype(v))
        rosenbrock(x)
    elseif isa(v, AbstractDict)
        rosenbrock(x[:a]) + rosenbrock(x[:b])
    else
        rosenbrock(x[1]) + rosenbrock(x[2])
    end
end

rosengrad = gradloss(rosenmulti)
Random.seed!(123456789) # TODO: tests sensitive to random seed ???

function rosenopt(w, params; verbose=false, ftol = 1e-3, xtol = 1e-10, maxiter = 12000)
    i = 1
    prev = 0
    current = 1
    t0 = time()
    while i <= maxiter && abs(current - prev) > xtol && current > ftol
	prev = current
	g, current = rosengrad(w)
	update!(w, g, params)
	i += 1
    end
    t1 = time()
    if verbose
        @printf("%s: f=%f iter=%-5d time=%.2f type=%s opt=%s\n",
                (current <= ftol ? "PASS" : "FAIL"),
                current, i-1, t1-t0, typeof(w), typeof(params))
    end
    return current <= ftol
end

@testset "update!" begin

    dims = 6
    w = randn(dims)
    # CPU Tests
    @test rosenopt(copy(w),SGD(lr=0.0005))
    @test rosenopt(copy(w),Momentum(lr=0.00025, gamma=0.95))
    @test rosenopt(copy(w),Nesterov(lr=0.00025, gamma=0.95))
    @test rosenopt(copy(w),Adam(lr=0.005, beta1=0.9, beta2=0.95, eps=1e-8))
    @test rosenopt(copy(w),Adagrad(lr=0.35, eps=1e-6))
    @test rosenopt(copy(w),Adadelta(lr=0.01, rho=0.5, eps=1e-6))
    @test rosenopt(copy(w),Rmsprop(lr=0.0005, rho=0.9, eps=1e-6))
    adam() = Adam(lr=0.005, beta1=0.9, beta2=0.95, eps=1e-8)
    v = 10*randn(dims)
    @test rosenopt((copy(w),copy(v)), (adam(),adam()))
    @test rosenopt(Any[copy(w),copy(v)], [adam(),adam()])
    @test rosenopt(Dict(:a=>copy(w),:b=>copy(v)), Dict(:a=>adam(),:b=>adam()))
    if gpu() >= 0
        Knet.gc()
        w = KnetArray(w) #GPU Tests
        @test rosenopt(copy(w),SGD(lr=0.0005))
        @test rosenopt(copy(w),Momentum(lr=0.00025, gamma=0.95))
        @test rosenopt(copy(w),Nesterov(lr=0.00025, gamma=0.95))
        @test rosenopt(copy(w),Adam(lr=0.005, beta1=0.9, beta2=0.95, eps=1e-8))
        @test rosenopt(copy(w),Adagrad(lr=0.35, eps=1e-6))
        @test rosenopt(copy(w),Adadelta(lr=0.01, rho=0.5, eps=1e-6))
        @test rosenopt(copy(w),Rmsprop(lr=0.0005, rho=0.9, eps=1e-6))
        v = KnetArray(v)
        @test rosenopt((copy(w),copy(v)), (adam(),adam()))
        @test rosenopt(Any[copy(w),copy(v)], [adam(),adam()])
        @test rosenopt(Dict(:a=>copy(w),:b=>copy(v)), Dict(:a=>adam(),:b=>adam()))
    end
end

# Issue 305: Using Knet with StaticArrays
# PR 306: optimizer for AbstractArrays
mutable struct A305 <: AbstractArray{Float32,1} end

@testset "optimizers" begin
    opt = optimizers(A305(), Adam)
    @test typeof(opt) == Adam
    opt = optimizers([A305(), A305()], Adam)
    @test typeof.(opt) == [Adam, Adam]
end

nothing


text:
export  excelSubpopAgeReport

function excelSubpopAgeReport( mpSim::MPsim, timeGrid::Vector{T}, ageRes::Real,
    ageType::Symbol, subpops::Subpopulation...;
    filename::AbstractString = "compReport", overwrite::Bool = true,
    timeFactor::Real = 12.0 ) where T <: Real

    filename, overwrite = setupFile( filename, overwrite )
    tStart = now()
    ageReport = subpopulationAgeReport( mpSim, timeGrid, ageRes, ageType,
        subpops... )
    reportGenerationTime = ( now() - tStart ).value / 1000.0

    if isempty( ageReport ) || all( subpop -> isempty( ageReport[ subpop ] ),
        collect( keys( ageReport ) ) )
        return
    end  # if isempty( ageReport ) || ...

    tStart = now()

    XLSX.openxlsx( filename, mode = overwrite ? "w" : "rw" ) do xf
        dateStr = Dates.format( now(), "yyyymmdd HHMMSS" )
        ws = xf[ 1 ]
        nReports = 1

        # Prep the new sheet.
        if overwrite
            XLSX.rename!( ws, string( "Report 1  ", dateStr ) )
        else
            nReports = length( filter( name -> startswith( name, "Report" ),
                XLSX.sheetnames( xf ) ) ) + 1
            ws = XLSX.addsheet!( xf, string( "Report ", nReports, "  ",
                dateStr ) )
        end  # if overwrite

        # Sheet header
        ws[ "A1" ] = "Simulation length"
        ws[ "B1" ] = mpSim.simLength / timeFactor
        ws[ "C1" ] = "years"
        ws[ "A2" ] = "Data generation time"
        ws[ "B2" ] = reportGenerationTime
        ws[ "C2" ] = "seconds"
        ws[ "A3" ] = "Excel generation time"
        ws[ "C3" ] = "seconds"
        ws[ "A4" ] = "Subpopulations"

        subpopNames = getfield.( vcat( subpops... ), :name )
        ii = 3

        for subpop in filter( subpop -> haskey( ageReport, subpop ),
            subpopNames )
            report = ageReport[ subpop ]

            if isempty( report )
                continue
            end  # if isempty( report )

            ii += 1
            ws[ ii, 2 ] = subpop

            rws = XLSX.addsheet!( xf,
                string( subpop, " ages  Report ", nReports ) )
            nrows, ncols = size( report )
            jj = 0

            for name in names( report )
                jj += 1
                val = tryparse( Float64, string( name ) )
                rws[ 1, jj ] = val isa Float64 ? val / timeFactor :
                    string( name )

                for kk in 1:nrows
                    rws[ kk + 1, jj ] = report[ kk, name ] /
                        ( val isa Float64 ? 1.0 : timeFactor )
                end  # for kk in 1:nrows
            end
        end  # for subpop in filter( ... )

        excelGenerationTime = ( now() - tStart ).value / 1000.0
        ws[ "B3" ] = excelGenerationTime
    end  # XLSX.openxlsx( ... ) do xf

end  # excelSubpopAgeReport( mpSim, timeGrid, ageRes, ageType, nodes, filename,
     #   overwrite, timeFactor )

excelSubpopAgeReport( mpSim::MPsim, timeRes::Real, ageRes::Real,
    ageType::Symbol, subpops::Subpopulation...;
    filename::AbstractString = "compReport", overwrite::Bool = true,
    timeFactor::Real = 12.0 ) =
    excelSubpopAgeReport( mpSim, MP.generateTimeGrid( mpSim, timeRes ), ageRes,
    ageType, subpops...; filename = filename, overwrite = overwrite,
    timeFactor = timeFactor )


text:
using Test
using ComplexMixtures, PDBTools
const CM = ComplexMixtures

@testset "PDB" begin

  #
  # Tests with trajectory given in a PDB file
  #

  dir="./data/PDB"
  atoms = readPDB("$dir/trajectory.pdb","model 1")
  options = Options(stride=1,seed=321,StableRNG=true,nthreads=1,silent=true)

  # Example 1: protein-tmao

  # save(R,"$dir/protein_tmao.json")
  R_save = load("$dir/protein_tmao.json")
  protein = Selection(select(atoms,"protein"),nmols=1)
  tmao = Selection(select(atoms,"resname TMAO"),natomspermol=14)
  traj = Trajectory("$dir/trajectory.pdb",protein,tmao,format="PDBTraj")
  R = mddf(traj,options)
  @test isapprox(R,R_save,debug=true)

  # Example 3: tmao-tmao

  # save(R,"$dir/tmao_tmao.json")
  R_save = load("$dir/tmao_tmao.json")
  tmao = Selection(select(atoms,"resname TMAO"),natomspermol=14)
  traj = Trajectory("$dir/trajectory.pdb",tmao,format="PDBTraj")
  R = mddf(traj,options)
  @test isapprox(R,R_save,debug=true)

end


text:
using DataDrivenControl
using Test

@testset "DataDrivenControl.jl" begin
    include("utils/utils.jl")
    include("irl/irl.jl")
end


text:
module SampleMockedTests

include("utils.jl")
using .Utils: @defmock

@defmock test
@defmock test_2

module TestNested
using ..Utils: @defmock
@defmock test
@defmock test_2
end  # module TestNested

module TestNested2
module TestNested3
using ...Utils: @defmock
@defmock test
@defmock test_2
end  # module TestNested3
end  # module TestNested2

module TestShouldTest
using ..Utils: @defmock, @getbool
should_test_module() = @getbool("should_test_module")
@defmock test
end  # module TestShouldTest

module NotATestModule
test() = error("this should not be executed")
end  # module NotATestModule

include("test_module_context.jl")

end  # module SampleMockedTests


text:
"""
	readtsf(filein;unitsonly)
Load tsf file and store it in DataFrame

**Input**
* filein: full name of the input file
* unitsonly: switch to output either data or units only (by default false=>data only)
* channelname: "measurement" (default)= set dataframe columns to measurement name, "site_measurement" = set columns to site_measurement

**Output**
* dataframe containing columns + timevector in DateTime format
* data units: if input `unitsonly` paramterer set to `true`

**Example**
```
data = readtsf("../test/input/tsf_data.tsf");
units = readtsf("../test/input/tsf_data.tsf",unitsonly=true);
```
"""
function readtsf(filein::String;unitsonly::Bool=false,channelname="measurement")
	# Set default/intial values
	undetval = 9999.999;
	countinfo::Int64 = 0; # count data lines
	count_header::Int64 = 0;
	channels = Vector{String}();
	units = Vector{String}();

	# Set aux Function for channel and unit line extraction
	function extract_channels(fid::IOStream,row::String,out::Vector{String})
		row = readline(fid);count_header += 1;
		while !occursin("[",row)
			if length(row) > 0
				push!(out,row);
			end
			row = readline(fid);count_header += 1;
			if eof(fid)
				row = "[";# =>will stop while loop
			end
		end
		return out,row
	end

	# Open and read the file
	open(filein) do fid
		row = readline(fid);count_header += 1;
	    while !eof(fid)
			# Get important header info
			if isempty(row)
				row = "empty row not to be parsed"
			elseif occursin("[UNDETVAL]",row)
	            undetval = Base.parse(Float64,row[11:end]);
				row = readline(fid);count_header += 1;
			elseif occursin("[COUNTINFO]",row)
				countinfo = Base.parse(Int,row[12:end]);
				row = readline(fid);count_header += 1;
			elseif occursin("[CHANNELS]",row)
				# Now get channel names
				channels,row = extract_channels(fid,row,channels);
			elseif occursin("[UNITS]",row)
				# Now get channel units
				units,row = extract_channels(fid,row,units);
			elseif occursin("[DATA]",row)
				# do not read using readline function (low performance)
				break;
			else
				row = readline(fid);count_header += 1;
			end
	    end
	end # close reading line by line
	# Correct/get channel names and units
	channels = correct_channels(channels,":",channelname=channelname);
	units = correct_channels(units,":");
	if unitsonly
		return units;
	else
		# Read all data (will be sorted afterwards)
		dataall = readdlm(filein,skipstart=count_header,header=false);
		# Convert time to datetime
		data = DataFrame(datetime=mat2time(dataall[:,1:6]))
		# Fill columns
		for i = 7:size(dataall,2)
			channame = i-6 <= length(channels) ? channels[i-6] : "measurement"*string(i-6)
			data[Symbol(channame)] = dataall[:,i];
			# Remove Undetval
			data[Symbol(channame)][findall(x->x.==undetval,data[Symbol(channame)])] .= NaN;
		end
		return data
	end
end # readtsf

"""
	writetsf(data,fileout,units,decimal,comment)
Write DataFrame to tsf format

**Input**
* data: dataframe containing columns + timevector in DateTime format
* fileout: full name of the output file
* units: optional output units for all columns
* decimal: optional output precision for either all or individual columns. Maximum precision is %.10g!
* comment: optional output comment

**Example**
```
data = DataFrame(temp=[10.,11.,12.,14.],grav=[9.8123,9.9,NaN,9.7],
       datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,1),
           DateTime(2010,1,1,2),DateTime(2010,1,1,4)]);
writetsf(data,"../test/output/tsf_data.tsf",units=["degC","nm/s^2"],
			comment=["first line","second line"],decimal=[1,3]);
```
"""
function writetsf(data::DataFrame,fileout::String;
					units=[],decimal=[4],comment=[]);
	# find datetime column + use only numberic values for output
	channels,timei = findchannels(data);
    # Write header
	open(fileout,"w") do fid
    	@printf(fid,"[TSF-file] v01.0\n\n");
		flagval = "9999.999";
    	@printf(fid,"[UNDETVAL] %s\n\n",flagval);
	    @printf(fid,"[TIMEFORMAT] DATETIME\n\n");
	    # Compute time resolution and write to header
		increment = round(Dates.value(data[timei][2]-data[timei][1])./1000);
    	@printf(fid,"[INCREMENT] %6.0f\n\n",increment);
     	# Add channel names and units
    	@printf(fid,"[CHANNELS]\n");
	    for i in channels
			@printf(fid,"Location:Instrument:%s\n",i);
        end
        @printf(fid,"\n[UNITS]\n");
		for i = 1:length(channels)
			if !isempty(units)
				@printf(fid,"%s\n",units[i]);
			else
				@printf(fid,"units?\n");
			end
		end
		# Add final comment if on input
		@printf(fid,"\n[COMMENT]\n");
	    if !isempty(comment)
            for i in comment
                @printf(fid,"%s\n",i)
            end
	    end
		# Number of data point + start Data
		nlines = size(data,1);
    	@printf(fid,"\n\n[COUNTINFO] %i\n\n",nlines);
    	@printf(fid,"[DATA]\n");
		# Round data to output precision
		dataout = round2output(data,decimal,channels)
		for i = 1:nlines
			# Write date + time (do not use Dates.format == very slow)
			@printf(fid,"%04i %02i %02i %02i %02i %02i ",
						Dates.year(data[timei][i]),
						Dates.month(data[timei][i]),
						Dates.day(data[timei][i]),
						Dates.hour(data[timei][i]),
						Dates.minute(data[timei][i]),
						Dates.second(data[timei][i]));
			# add data
			for j in channels
				if ismissing(dataout[j][i]) || isnan(dataout[j][i])
					@printf(fid," %s",flagval);
				else
					@printf(fid," %.10g",dataout[j][i]);
				end
			end
			@printf(fid,"\n");
		end
	end # open file
end

"""
	mat2time(mat)
Auxiliary function to convert time in matrix to datetime format
"""
function mat2time(mat::Matrix{Float64})
	if eltype(mat) != Int
		temp = round.(Int64,mat)
	end
	if size(mat,2) == 3
		return DateTime.(mat[:,1],mat[:,2],mat[:,3])
	else
		return DateTime.(mat[:,1],mat[:,2],mat[:,3],
						 mat[:,4],mat[:,5],mat[:,6])
	end
end

"""
	correct_channels(in_string)
Auxiliary function to get channel names/units from input vector of strings
"""
function correct_channels(in_text::Vector{String},sp::String;
							channelname="measurement")
	out = Vector{String}();
	c = 2; # count channels with identical name
	for i in in_text
		temp = split(i,sp);
		# remove empty spaces and other symbols
		name = replace(temp[end],r" |\r|\""=>"")
		if channelname == "site_measurement"
			site = replace(temp[1],r" |\r|\""=>"")*"_";
			chan_number = ""
		else
			site = "";
			chan_number = string(c);
		end
		# check if the name already exists
		if any(occursin.(name,out)) || channelname == "site_measurement"
			push!(out,site*name*chan_number);c += 1;
		else
			push!(out,name)
		end
	end
	return out;
end

"""
	findchannels(data)
Auxiliary function to find columns with channel names and datetime
"""
function findchannels(data::DataFrame)
	# find datetime column + use only numberic values for output
	channels = eltype(names(data))==String ? Vector{String}() : Vector{Symbol}();
	timei = 1;
	for i in names(data)
		if eltype(data[i]) == DateTime  || eltype(data[i]) == Union{DateTime, Missings.Missing}
			timei = i;
		elseif eltype(data[i]) == Int || eltype(data[i]) == Float64 ||
			   eltype(data[i]) == Union{Int, Missings.Missing} ||
			   eltype(data[i]) == Union{Float64, Missings.Missing}
			push!(channels,i);
		end
	end
	return channels,timei;
end

"""
	round2output(data,decimal,channels)
Auxiliary function to round data to output precision
"""
function round2output(data::DataFrame,decimal,channels)
	dataout = deepcopy(data);
	for (i,v) in enumerate(channels)
		if length(channels) == length(decimal)
			dataout[v] = round.(data[v].*10^decimal[i])./10^decimal[i];
		elseif length(decimal) == 1
			dataout[v] = round.(data[v].*10^decimal[1])./10^decimal[1];
		end
	end
	return dataout;
end