text:
"""
    mirror_tarball(registry, upstreams::AbstractVector, static_dir::String; kwargs...)

Generating all static contents for registry `registry`.

After building, all static content data are saved to `static_dir`:

* `/registries`: map of registry uuids at this server to their current tree hashes
* `/registry/\$uuid/\$hash`: tarball of registry uuid at the given tree hash
* `/package/\$uuid/\$hash`: tarball of package uuid at the given tree hash
* `/artifact/\$hash`: tarball of an artifact with the given tree hash

Set environment variable `JULIA_NUM_THREADS` to enable multi-threads.

# Keyword parameters

- `http_parameters::Dict{Symbol, Any}`: any parameters that need to pass to `HTTP.get` when fetching resources.
- `packages::AbstractVector{Package}`: manually create and specify a set of packages that needed to be stored. This
   can be used to build only a partial of the complete storage. Check [`read_packages`](@ref read_packages) for how
   to build packages info.
- `show_progress::Bool`: `true` to show an additional progress meter. By default it is `true`.
- `skip_duration::Real`: specify for how long (hours) you want to skip the resources listed in `/failed_resources.txt`
   until the next try of them. By default it is `24` hours.

"""
mirror_tarball(registry::Tuple, args...; kwargs...) = mirror_tarball(RegistryMeta(registry...), args...;kwargs...)

function mirror_tarball(
        registry::RegistryMeta,
        upstreams::AbstractVector,
        static_dir::AbstractString;
        http_parameters::Dict{Symbol, Any} = Dict{Symbol, Any}(),
        packages::Union{AbstractVector, Nothing} = nothing,
        registry_hash::Union{AbstractString, Nothing} = nothing,
        skip_duration = 24,
        show_progress = true
)
    ### Except for the complex error handling strategy, the mirror routine
    ###   1. query upstreams to get the latest hash
    ###   2. fetch registry tarball and save it as `registry/$uuid/$hash`
    ###   3. read pacakges information from registry tarball, download and save all packages (but not
    ###      their artifacts) to `packages/$uuid/$hash`
    ###   4. decompress and extract `*Artifacts.toml` files from each tarball in
    ###      `packages/$uuid/$hash`, and download all artifacts to `artifact/$hash`
    ###   5. update `/registries` so that downstream Pkg client/server knows that
    ###
    ### There will be some caching files:
    ###   * `/failed_resources.txt` records resources that failed to download in the previous mirror
    ###     routine
    ###   * `.cache/package/$uuid/$hash/*Artifacts.toml` caches all `*Artifacts.toml` files so that
    ###     we don't need to waste the time and resources in decompressing and extracting them from
    ###     `/package/$uuid/$hash` -- This is a CPU and IO hotspot

    uuid = registry.uuid
    name = registry.name
    failed_logfile = joinpath(static_dir, "failed_resources.txt")
    blocklist_file = joinpath(static_dir, "blocklist.txt")
    upstream_str = join(upstreams, ", ")

    # Some resources vanishes and can never be downloaded, we skip them in the next 24 hours
    last_try_time = query_last_try_datetime(failed_logfile)
    skipped = now() - last_try_time < Hour(skip_duration)
    skipped_records = skipped ? read_records(failed_logfile) : Set()
    blocked_resourced = read_records(blocklist_file)
    function _download(resource, tarball; throw_warnings=true)
        if resource in skipped_records
            @info "fetching resource has failed in the last $skip_duration hours, skip it" resource date=now()
            return false
        end

        if resource in blocked_resourced
            @info "resource is listed in blocklist, skip it" resource date=now()
            return false
        end

        try
            show_progress && @info "downloading resource..." resource=resource tarball=tarball date=now()
            rst = download_and_verify(upstreams, resource, tarball; http_parameters=http_parameters, throw_warnings=throw_warnings)
            if rst !== false
                return true
            else
                show_progress && @info "skip resource fetching in next $skip_duration hours" resource date=now()
                log_to_failure(resource, failed_logfile)
                return false
            end
        catch err
            throw_warnings && @warn err resource=resource
            rm(tarball; force=true)
            show_progress && @info "skip resource fetching in next $skip_duration hours" resource date=now()
            log_to_failure(resource, failed_logfile)
            return false
        end
    end

    # 1. query latest registry hash
    upstreams = normalize_upstream.(upstreams)
    latest_hash = isnothing(registry_hash) ? query_latest_hash(registry, upstreams) : registry_hash
    if isnothing(latest_hash)
        @error "failed to query the latest registry hash" uuid=uuid hash=latest_hash upstreams=upstreams date=now()
        error("stop mirroring.")
    end

    # 2. fetch registry tarball
    resource = "/registry/$(uuid)/$(latest_hash)"
    tarball = joinpath(static_dir, "registry", uuid, latest_hash)
    rst = _download(resource, tarball)
    if !rst
        @error "failed to fetch registry tarball" uuid=uuid hash=latest_hash upstreams=upstreams date=now()
        error("stop mirroring.")
    end

    # 3. read and download `/package/$uuid/$hash`
    if isnothing(packages)
        packages = mktempdir() do tmpdir
            open(tarball, "r") do io
                Tar.extract(decompress(io), tmpdir)
            end
            # only returns packages that are not stored in static_dir
            read_packages(tmpdir; fetch_full_registry=false, static_dir=static_dir)
        end
    end

    num_versions = mapreduce(x -> length(x.versions), +, packages)
    @info "Start mirrorring" date=now() registry=name uuid=uuid hash=latest_hash num_versions=num_versions upstreams=upstream_str
    skipped && @info "$(length(skipped_records)) previously failed resources are skipped during this build" date=now()

    p = show_progress ? Progress(num_versions; desc="$name: Pulling packages: ") : nothing
    if !isempty(packages)
        function pkg_kernel(pkg)
            for (ver, hash_info) in pkg.versions
                tree_hash = hash_info["git-tree-sha1"]
                resource = "/package/$(pkg.uuid)/$(tree_hash)"
                tarball = joinpath(static_dir, "package", pkg.uuid, tree_hash)
                _download(resource, tarball)
                isnothing(p) || ProgressMeter.next!(p; showvalues = [(:date, now()), (:package, pkg.name), (:version, ver), (:uuid, pkg.uuid), (:hash, tree_hash)])
            end
        end
        if use_threads()
            ThreadPools.@qthreads for pkg in packages
                pkg_kernel(pkg)
            end
        else
            foreach(pkg_kernel, packages)
        end
    end

    # 4. read and download `/artifact/$hash`
    # only returns artifacts that are not stored in static_dir
    @info "querying artifacts info, this may take some time..." date=now()
    artifacts = query_artifacts(static_dir; fetch_full=false)
    @info "found $(length(artifacts)) new artifacts" date=now()
    p = show_progress ? Progress(length(artifacts); desc="$name: Pulling artifacts: ") : nothing
    batch_size = min(length(artifacts), max(20, ceil(Int, length(artifacts)/(5*Threads.nthreads()))))
    if !isempty(artifacts)
        function artifact_kernel(artifact)
            if is_valid(artifact)
                tree_hash = artifact.hash
                resource = "/artifact/$(tree_hash)"
                tarball = joinpath(static_dir, "artifact", tree_hash)

                _download(resource, tarball)
                isnothing(p) || ProgressMeter.next!(p; showvalues = [(:date, now()), (:artifact, tree_hash), (:resource, resource), (:tarball, tarball)])
            end
        end
        if use_threads()
            ThreadPools.@qthreads for artifact in artifacts
                artifact_kernel(artifact)
            end
        else
            foreach(artifact_kernel, artifacts)
        end
    end

    # update /registries after updates
    registries_file = joinpath(static_dir, "registries")
    update_registries(registries_file, uuid, latest_hash)

    # clean up
    # Glob.jl doesn't support ** syntax yet: https://github.com/vtjnash/Glob.jl/issues/19
    foreach(glob(glob"artifact/*.tmp.*", static_dir)) do tarball
        rm(tarball; force=true)
    end
    foreach(glob(glob"*/*/*.tmp.*", static_dir)) do tarball
        # /registry and /package
        rm(tarball; force=true)
    end

    last_try_time = skipped ? last_try_time : now()
    failed_records = read_records(failed_logfile)
    buffer = IOBuffer()
    println(buffer, last_try_time)
    foreach(failed_records) do resource
        if !isfile(joinpath(static_dir, resource[2:end])) # joinpath(pwd(), "/a") == "/a"
            println(buffer, resource)
        end
    end
    open(failed_logfile, "w") do io
        write(io, take!(buffer))
    end
    @info("Mirror completed. There are $(length(failed_records) - length(skipped_records)) new fail-to-fetch resources and will be skipped in the next $(skip_duration) hours.",
        date=now(), registry=name, uuid=uuid, hash=latest_hash, upstreams=upstream_str)

    return latest_hash
end

### helpers
function log_to_failure(resource::AbstractString, logfile)
    try
        if VERSION < v"1.5"
            open(logfile, "a") do io
                println(io, resource)
            end
        else
            open(logfile, "a"; lock=true) do io
                println(io, resource)
            end
        end
    catch err
        @warn err resource=resource
    end
end

function query_last_try_datetime(logfile)
    isfile(logfile) || return DateTime(0)

    line = readline(logfile)
    try
        return DateTime(line)
    catch err
        if err isa ArgumentError
            return DateTime(0)
        end
        rethrow(err)
    end
end

function read_records(logfile)
    isfile(logfile) || return Set()
    records = readlines(logfile)
    isempty(records) && return Set()

    filter(sort(collect(Set(records)))) do resource
        !isnothing(match(resource_re, resource))
    end
end

function update_registries(registries_file, uuid, hash)
    # each line is recorded in format: `/registry/$uuid/$hash`
    get_uuid(line) = split(line, "/")[3]

    if isfile(registries_file)
        registries = Dict(get_uuid(item) => item for item in readlines(registries_file))
    else
        registries = Dict()
    end
    registries[uuid] = "/registry/$uuid/$hash"
    open(registries_file, write = true) do io
        for (uuid, registry_path) in sort!(collect(registries))
            println(io, registry_path)
        end
    end
end

function query_artifacts(static_dir; fetch_full=false)
    cache_root = joinpath(static_dir, ".cache")
    package_re = Regex("package/($uuid_re)/($hash_re)")
    tarball_glob_pattern = ["package", Regex(uuid_re), Regex(hash_re)]
    artifact_glob_pattern = ["package", Regex(uuid_re), Regex(hash_re), r"[Julia]?Artifacts\.toml"]

    # incrementally extract Artifacts.toml
    function artifact_kernel(pkg_tarball)
        pkg_uuid, pkg_hash = match(package_re, replace(pkg_tarball, Base.Filesystem.path_separator=>"/")).captures
        cache_dir = joinpath(cache_root, "package", pkg_uuid, pkg_hash)

        with_cache_dir(cache_dir) do
            # if cache_dir already exists, skip the decompress and extract
            # https://github.com/johnnychen94/StorageMirrorServer.jl/issues/5
            mkpath(cache_dir)
            try
                @info "extract Artifacts.toml from source package" date=now() uuid=pkg_uuid hash=pkg_hash
                open(pkg_tarball, "r") do io
                    Tar.extract(decompress(io), cache_dir) do hdr
                        splitpath(hdr.path)[end] in artifact_names
                    end
                end
            catch err
                @warn "failed to extract package tarball" uuid=pkg_uuid hash=pkg_hash
                rm(pkg_tarball; force=true) # remove it if this tarball is broken
                rm(cache_dir; force=true, recursive=true)
                return false
            end

            # Although not all package has `*Artifacts.toml`, we still keep the empty `cache_dir`
            # for them as a placeholder to notify `with_cache_dir` that it can directly use caches
            # and skip the extraction call
            return true # to notify with_cache_dir that this function call success
        end
    end
    if use_threads()
        ThreadPools.@qbthreads for pkg_tarball in glob(tarball_glob_pattern, static_dir)
            artifact_kernel(pkg_tarball)
        end
    else
        foreach(artifact_kernel, glob(tarball_glob_pattern, static_dir))
    end

    # read all artifacts
    all_artifacts = Artifact[]
    for artifact_toml in glob(artifact_glob_pattern, cache_root)
        for (key, artifact_info) in TOML.parsefile(artifact_toml)
            if artifact_info isa Dict
                # platform independent artifacts, e.g., TestImages
                tree_hash = artifact_info["git-tree-sha1"]
                tarball = joinpath(static_dir, "artifact", tree_hash)
                if fetch_full || !isfile(tarball)
                    push!(all_artifacts, Artifact(tree_hash))
                end
            elseif artifact_info isa AbstractVector
                # platform dependent artifacts, e.g., FFTW_jll
                for x in artifact_info
                    tree_hash = x["git-tree-sha1"]
                    tarball = joinpath(static_dir, "artifact", tree_hash)
                    if fetch_full || !isfile(tarball)
                        push!(all_artifacts, Artifact(tree_hash))
                    end
                end
            else
                @warn "invalid artifact toml" package=pkg version=ver artifact=artifact_info
            end
        end
    end

    # the same artifacts might be recorded in multiple versions
    return Set(all_artifacts)
end

function is_valid(artifact::Artifact)
    isempty(artifact.hash) && return false

    return true
end

use_threads() = get(ENV, "JULIA_NUM_THREADS", 1) != 1

ext: jl


text:
#!/home/fng514/bin/julia
#$ -S /home/fng514/bin/julia
#$ -cwd
# Christian Theil Have, 2017.
#
using ArgParse
using BGZFStreams

@everywhere function process_line(line1,lineno)
	if (ismatch(r"^#", line1))
		(0, line1)
	else
		(lineno, line1)
	end
end

# FIXME: We may need to do this with gzip readers instead
function main(args)
	s = ArgParseSettings("Example 2 for merge.jl: " *  # description
		"flags, options help, " *
		"required arguments.")

	@add_arg_table s begin
		"input-vcf"
		"--lines"
		arg_type = Int
		default = 1000
		"--prefix"
	end

	parsed_args = parse_args(s) # the result is a Dict{String,Any}
	println("Parsed args:")
	for (key,val) in parsed_args
		println("  $key  =>  $(repr(val))")
	end

	vcf1=BGZFStream(parsed_args["input-vcf"])
	out = BGZFStream(parsed_args["prefix"],"w")
	lines1 = eachline(vcf1)

	state1 = start(lines1)

	keep_positions1 = nothing

	header_lines = []
	content_lines = []

	fileno = 1
	lineno = 0

	while !done(lines1, nothing)
		(current_line1, state1) = next(lines1, state1)
		(lineno, processed_line)  =  process_line(current_line1, lineno+1)

		if lineno == 0
			push!(header_lines, processed_line)
		elseif lineno % parsed_args["lines"] == 0
			push!(content_lines,processed_line)
			out_name = string(parsed_args["prefix"],lpad(fileno, 5, 0),".vcf.gz")
			out = BGZFStream(out_name, "w")
			for line in header_lines
				write(out,line)
			end
			for line in content_lines
				write(out,line)
			end
			close(out)
			content_lines = []
			fileno = fileno + 1
			println("wrote $out_name")
		else
			push!(content_lines,processed_line)
		end
	end

	if length(content_lines) > 0
		out_name = string(parsed_args["prefix"],lpad(fileno, 5, 0),".vcf.gz")
		out = BGZFStream(out_name, "w")
		for line in header_lines
			write(out,line)
		end
		for line in content_lines
			write(out,line)
		end
		close(out)
		println("wrote $out_name")
	end

	close(vcf1)
end

main(ARGS)

max_stars_repo_name: cth/merge-imputations

max_issues_repo_name: cth/merge-imputations

max_forks_repo_name: cth/merge-imputations


text:
module EnsemblesTest

using OnlineAI, FactCheck, Distributions

function xor_data()
  inputs = [0 0; 0 1; 1 0; 1 1]
  targets = float(sum(inputs,2) .== 1)

  # all sets are the same
  inputs = inputs .- mean(inputs,1)
  DataPoints(inputs, targets)
end

# function testxor(; hiddenLayerNodes = [2],
#                    hiddenMapping = SigmoidMapping(),
#                    finalMapping = IdentityMapping(),
#                    params = NetParams(),
#                    solverParams = SolverParams(maxiter=10000))

#   # all xor inputs and results
#   inputs = [0 0; 0 1; 1 0; 1 1]
#   targets = float(sum(inputs,2) .== 1)

#   # all sets are the same
#   inputs = inputs .- mean(inputs,1)
#   data = DataPoints(inputs, targets)
#   sampler = SimpleSampler(data)

#   # hiddenLayerNodes = [2]
#   net = buildRegressionNet(ncols(inputs),
#                            ncols(targets),
#                            hiddenLayerNodes;
#                            hiddenMapping = hiddenMapping,
#                            finalMapping = finalMapping,
#                            params = params)
#   show(net)

#   stats = solve!(net, solverParams, sampler, sampler)

#   # output = Float64[predict(net, d.input)[1] for d in data]
#   output = vec(predict(net, inputs))
#   for (o, d) in zip(output, data)
#     println("Result: input=$(d.x) target=$(d.y) output=$o")
#   end

#   net, output, stats
# end

facts("Ensembles") do

  # minerror = 0.05
  # solverParams = SolverParams(maxiter=10000, minerror=minerror*0.8)

  # net, output, stats = testxor(params=NetParams(μ=0.0, λ=0.0, mloss=L1CostModel()), solverParams=solverParams)
  # @fact net --> anything
  # @fact output --> roughly([0., 1., 1., 0.], atol=0.05)

  ps = ParameterSampler(Dict(zip([:x, :y, :z, :a], [Constant(1), Normal(), Uniform(3,4), VectorSampler(2:4, [0.8,0.1,0.1])])))
  show(ps.d)
  for i in 1:20
    println(sample(ps))
  end
  f(;kwargs...) = println(Dict(kwargs))
  f(;sample(ps)...)

end # facts

end # module


text:
using InvariantCausal
using DelimitedFiles
using Test

include(joinpath(@__DIR__, "test_search.jl"))
include(joinpath(@__DIR__, "test_regression.jl"))

X = readdlm(joinpath(@__DIR__, "X1.dat"))
env = Vector{Int}(X[:,1])
X = X[:,2:end]
S = 1:size(X,2)

@time @testset "causal search" begin
    result = map(i -> causalSearch(X, X[:, i], env, setdiff(S,i) , α=0.01), S)
    @test result[2].S == [5]
    @test result[3].S == [5]
    @test result[7].S == [4] || result[7].S == []
    for i in [1, 4, 5, 6]
        @test result[i].model_reject == true
    end
end

@time @testset "causal search with selection_only" begin
    result = map(i -> causalSearch(X, X[:, i], env, setdiff(S,i), α=0.01, selection_only=true), S)
    @test result[2].S == [5]
    @test result[3].S == [5]
    @test result[7].S == [4] || result[7].S == []
    for i in [1, 4, 5, 6]
        @test result[i].model_reject == true
    end
end

@time @testset "causal search with limited # of true causes" begin
    result = map(i -> causalSearch(X, X[:, i], env, setdiff(S,i), α=0.01, max_num_true_causes=3), S)
    @test result[2].S == [5]
    @test result[3].S == [5]
    @test result[7].S == [4] || result[7].S == []
    for i in [1, 4, 5, 6]
        @test result[i].model_reject == true
    end
end

@time @testset "causal search with subsampling" begin
    result = map(i -> causalSearch(X, X[:, i], env, setdiff(S,i), α=0.01, n_max_for_exact=90), S)
    @test result[2].S == [5]
    @test result[3].S == [5] || result[3].S == []
    @test result[7].S == [4] || result[7].S == []
    for i in [1, 4, 5, 6]
        @test result[i].model_reject == true
    end
end

function generate_setting(setting_configs; n_environments=2)
    n_int = rand(setting_configs["n_int"])
    n_obs = rand(setting_configs["n_obs"])
    p = rand(setting_configs["p"])
    k = rand(setting_configs["k"])
    lb_obs = rand(setting_configs["lb_obs"])
    ub_obs = lb_obs + rand(setting_configs["ub_lb_delta_obs"])
    _var_1 = rand(setting_configs["err_var_min"])
    _var_2 = rand(setting_configs["err_var_min"])
    err_var_min = min(_var_1, _var_2)
    err_var_max = max(_var_1, _var_2)
    noise_multiplier_min = rand(setting_configs["noise_multiplier_min"])
    if rand() < setting_configs["prob_fixed_noise_multiplier"]
        noise_multiplier_max = noise_multiplier_min
    else
        noise_multiplier_max = noise_multiplier_min + rand(setting_configs["noise_multiplier_max_min_delta"])
    end
    _l = rand(setting_configs["lb_int"])
    _u = rand(setting_configs["ub_int"])
    lb_int = min(_l, _u)
    ub_int = max(_l, _u)
    if rand() < setting_configs["prob_int_single"]
        n_int_variables = 1
        frac_int_variables = 1 / p
    else
        frac_int_variables = rand(setting_configs["frac_int_variables"])
        n_int_variables = floor(Int64, p * frac_int_variables)
    end
    prob_coefficient_unchanged = setting_configs["prob_coefficient_unchanged"]
    target = rand(1:p)
    SEMs = Vector{GaussianSEM}(undef, n_environments)
    env = repeat([1], inner=n_obs)
    intervened_variables = Vector{Vector{Int64}}(undef, n_environments)
    intervened_variables[1] = []
    SEMs[1] = random_gaussian_SEM(p, k,
              lb=lb_obs, ub=ub_obs, var_min=err_var_min, var_max=err_var_max)
    true_β = SEMs[1].B[target, setdiff(1:p, target)]
    true_causes = (1:(p-1))[true_β .!= 0]
    for j in 2:n_environments
        SEMs[j], intervened_variables[j] = random_noise_intervened_SEM(SEMs[1], p_intervened=n_int_variables,
                noise_multiplier_min=noise_multiplier_min, noise_multiplier_max=noise_multiplier_max,
                avoid=[target], prob_coeff_unchanged=prob_coefficient_unchanged,
                lb=lb_int, ub=ub_int)
        env = vcat(env, repeat([j], inner=n_int))
    end
    return Dict(
        "n_environments" => n_environments,
        "env" => env,
        "SEMs" => SEMs,
        "intervened_variables" => intervened_variables,
        "target" => target,
        "true_causes" => true_causes,
        "true_β" => true_β,
        "n_int" => n_int,
        "n_obs" => n_obs,
        "p" => p,
        "k" => k,
        "lb_obs" => lb_obs,
        "ub_obs" => ub_obs,
        "err_var_min" => err_var_min,
        "err_var_max" => err_var_max,
        "noise_multiplier_min" => noise_multiplier_min,
        "noise_multiplier_max" => noise_multiplier_max,
        "lb_int" => lb_int,
        "ub_int" => ub_int,
        "prob_coefficient_unchanged" => prob_coefficient_unchanged,
        "n_int_variables" => n_int_variables,
        "frac_int_variables" => frac_int_variables
    )
end

@time @testset "random instances" begin
    setting_configs = Dict(
    "n_int" => [100, 200, 300, 400, 500],
    "n_obs" => [100, 200, 300, 400, 500],
    "p" =>  collect(5:9),
    "k" => collect(1:4),
    "lb_obs" => collect(0.1:0.1:2),
    "ub_lb_delta_obs" => collect(0.1:0.1:1),
    "err_var_min" => collect(0.1:0.1:2),
    "err_var_max" => collect(0.1:0.1:2),
    "noise_multiplier_min" => collect(0.1:0.1:4),
    "noise_multiplier_max_min_delta" => collect(0.1:0.1:2),
    "lb_int" => collect(0.1:0.1:2),
    "ub_int" => collect(0.1:0.1:2),
    "prob_fixed_noise_multiplier" => 1/3,
    "prob_coefficient_unchanged" => 2/3,
    "prob_int_single" => 1/6,
    "frac_int_variables" => 1. ./ collect(1.1:0.1:3)
    )
    for trial in 1:10
        for n_environments in [2, 3, 5]
            setting = generate_setting(setting_configs, n_environments=n_environments)
            env = setting["env"]
            p = setting["p"]
            target = setting["target"]
            true_β = setting["true_β"]
            true_causes = setting["true_causes"]
            intervened_variables = setting["intervened_variables"]
            SEM_obs = setting["SEMs"][1]
            X = simulate(SEM_obs, setting["n_obs"])
            for j in 2:n_environments
                _X = simulate(setting["SEMs"][j], setting["n_int"])
                X = vcat(X, _X)
            end
            y = X[:, target]
            X = X[:, setdiff(1:p, target)]
            result = causalSearch(X, y, env, verbose=false)
            println("S = $(result.S), truth = $true_causes")
            @test issubset(result.S, true_causes)
            for k in result.S
                @test result.confint[k, 1] < true_β[k] < result.confint[k, 2]
            end
        end
    end
end

@time @testset "random instances, many variables with screening" begin
    setting_configs = Dict(
    "n_int" => [100, 200, 300, 400, 500],
    "n_obs" => [100, 200, 300, 400, 500],
    "p" =>  collect(20:50),
    "k" => collect(1:3),
    "lb_obs" => collect(0.1:0.1:2),
    "ub_lb_delta_obs" => collect(0.1:0.1:1),
    "err_var_min" => collect(0.1:0.1:2),
    "err_var_max" => collect(0.1:0.1:2),
    "noise_multiplier_min" => collect(0.1:0.1:4),
    "noise_multiplier_max_min_delta" => collect(0.1:0.1:2),
    "lb_int" => collect(0.1:0.1:2),
    "ub_int" => collect(0.1:0.1:2),
    "prob_fixed_noise_multiplier" => 1/3,
    "prob_coefficient_unchanged" => 2/3,
    "prob_int_single" => 1/6,
    "frac_int_variables" => 1. ./ collect(1.1:0.1:3)
    )
    for trial in 1:10
        for n_environments in [2, 3, 5]
            setting = generate_setting(setting_configs, n_environments=n_environments)
            env = setting["env"]
            p = setting["p"]
            target = setting["target"]
            true_β = setting["true_β"]
            true_causes = setting["true_causes"]
            intervened_variables = setting["intervened_variables"]
            SEM_obs = setting["SEMs"][1]
            X = simulate(SEM_obs, setting["n_obs"])
            for j in 2:n_environments
                _X = simulate(setting["SEMs"][j], setting["n_int"])
                X = vcat(X, _X)
            end
            y = X[:, target]
            X = X[:, setdiff(1:p, target)]
            result = causalSearch(X, y, env, verbose=false, p_max=8)
            # performance only guaranteed if lasso selected all the causes
            println("S = $(result.S), truth = $true_causes")
            if issubset(true_causes, result.variables_considered)
                @test issubset(result.S, true_causes)
            end
        end
    end
end

max_issues_repo_name: richardkwo/InvariantCausal

max_forks_repo_name: richardkwo/InvariantCausal


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.
##########################################################################

# P-median problem

pMedianDefaultOptions = Dict([:x_bin => true, :y_bin => false, :solver => "gurobi", :solver_args => [], :solver_kwargs => []])
try Gurobi catch; pMedianDefaultOptions[:solver] = "cbc" end

"""
	solvePMedian(n::Int, c::Array{Float,2}; options::Dict{Symbol,T} = pMedianDefaultOptions, results::Dict = Dict()) where T <: Any
Solves the p-median problem for `n` facilities and costs `c` where `c[i,j]` is the cost of serving demand point `j` from facility location `i`.
The p-median problem is to locate facilities in order to minimise the total cost of serving demand, where each demand point is served by the facility which can serve the point at lowest cost.
Returns a vector indicating which facility locations to use.
`results` will store results such as the objective value and decision variable values.
"""
function solvePMedian(n::Int, c::Array{Float,2}; options::Dict{Symbol,T} = pMedianDefaultOptions, results::Dict = Dict()) where T <: Any
	@assert(n >= 1)
	s, d = size(c) # number of potential facilities, number of demand points

	global pMedianDefaultOptions
	defaultOptions = copy(pMedianDefaultOptions)
	merge!(defaultOptions, options)

	m = model = Model()

	# solver
	solver = options[:solver]
	args = options[:solver_args]
	kwargs = options[:solver_kwargs]
	if solver == "cbc" set_optimizer(model, with_optimizer(Cbc.Optimizer, logLevel=0, args...; kwargs...))
	elseif solver == "glpk" set_optimizer(model, with_optimizer(GLPK.Optimizer, args...; kwargs...))
	elseif solver == "gurobi" @stdout_silent(set_optimizer(model, with_optimizer(Gurobi.Optimizer, OutputFlag=0, args...; kwargs...)))
	end

	# variables
	options[:x_bin] ? @variable(m, x[i=1:s], Bin) : @variable(m, 0 <= x[i=1:s] <= 1) # x[i] = 1 if one facility at location i, 0 otherwise
	options[:y_bin] ? @variable(m, y[i=1:s,j=1:d], Bin) : @variable(m, 0 <= y[i=1:s,j=1:d] <= 1) # y[i,j] = 1 if facility at location i serves point j

	@constraints(model, begin
		(maxFacilityCount, sum(x) <= n)
		(serveDemandFromOpenFacilities[i=1:s,j=1:d], y[i,j] <= x[i]) # can only serve point j from facility location i if it is opened
		(demandServedBySingleFacility[j=1:d], sum(y[:,j]) == 1)
	end)

	@expressions(model, begin
		cost, sum(c.*y)
	end)

	# solve
	@objective(model, Min, cost)
	@stdout_silent optimize!(model)
	@assert(termination_status(model) == MOI.OPTIMAL)

	results[:cost] = JuMP.value.(cost)
	results[:x] = JuMP.value.(x)
	results[:y] = JuMP.value.(y)
	results[:model] = model

	solution = convert(Vector{Int}, round.(results[:x])) # solution[i] = 1 if facility location i should be used, 0 otherwise
	return solution
end


text:
@testset "58.length-of-last-word.jl" begin
    @test length_of_last_word("Hello World") == 5
    @test length_of_last_word("   fly me   to   the moon  ") == 4
    @test length_of_last_word("luffy is still joyboy") == 6
end


text:
"""
    function rand_mate(Sires, Dams, nSib)
---
Ramdomly mate `Sires` and `Dams`, which are name vectors.
Each family will have `nSib` full sibs.
If the bigger to smaller ratio of length of `Sires` and `Dams` is
not integer, some ID may be used more often than others.
"""
function rand_mate(Sires, Dams, nSib)
    np, nm = length(Sires), length(Dams) # number of pa and ma
    nf = max(np, nm)                     # number of full sib families
    nid = nSib * nf
    pa = repeat(shuffle(Sires), inner=nSib, outer=Int(ceil(nf/np)))
    ma = repeat(shuffle(Dams), inner=nSib, outer=Int(ceil(nf/nm)))
    pa[1:nid], ma[1:nid]
end

"""
    function gene_drop(df, r, psnp, osnp, qtl, h²)
---
Given a dataframe `df` of a generation,

- this function drop genes in `prt_snp` specified by `df.sir` and `df.dam`
int `off_snp` specified by `df.id`.
When the genotypes of this genration are created,
the breeding values, and phenotypes are also calculated.
Before phenotypes were calculated, editing is done with editing information `ei`.
"""
function gene_drop(df, r, psnp, osnp, qtl, h², ei=(lc=0, sr=1))
    Threads.@threads for x in eachrow(df)
        gp = view(psnp, :, 2x.sir-1:2x.sir) # pa's genotypes
        gm = view(psnp, :, 2x.dam-1:2x.dam) # ma's genotypes
        op = view(osnp, :, 2x.id-1)         # id's haplotype from pa
        om = view(osnp, :, 2x.id)           # id's haplotype from ma
        ix = sample_alleles(r)
        copy!(op, [gp[i, ix[i]] for i in 1:length(op)])
        ix = sample_alleles(r)
        copy!(om, [gm[i, ix[i]] for i in 1:length(om)])
    end
    fra = 2first(df.id) - 1      # from column in SNP genotypes
    til = 2last(df.id)           # end column
    snp = view(osnp, :, fra:til)
    ei.lc != 0 && begin         # editing right after gene dropping
        tov = ei.lc > 0 ? 1 : 0
        br = rand(Bernoulli(ei.sr), til-fra+1)
        snp[abs(ei.lc), br] .= tov
    end
    df[:, :tbv], df[:, :p7e] = phenotype(snp, qtl, h²)
end

"""
    function toBinary(df, p17h)
---
convert the continuous phenotype (`:p7e`) in `df` to binary,
according to percentage-of-death (`:p17h`).
`:p17h` happened to be the quantile of the phenotypes.
"""
function toBinary(df, p17h)
    m = quantile(df.p7e, p17h)
    for i in 1:nrow(df)
        if df.p7e[i] < m
            df.p7e[i] = 0 # death
        else
            df.p7e[i] = 1 # live
        end
    end
end

"""
    function create_storage(base, par)
---
This is a refactorized function for creating storage.
My function before was too long, making it difficult for debugging.
This function is still very long, but well blockwised.

The storage is created in one-go.
This is to avoid too many times of memory allocations,
which can crash the program.

Finally knew how to use `view`.
Use `view`s for gene-dropping later.

This function:
1. create storage for SNP of production and challenge group
   - base SNP are copied to the production group as generation 0.
2. create pedigree generation 0 → 1 for the production group
   - which include the `base` population.
3. create pedigree for the challenge group.
4. generations are marked
5. sex sampling determined in good advance
6. determine parents of the 1st generation.
"""
function create_storage(base, par, qtl)
    # the pedigree are of 3 parts: the base, 1st generation, 2:end generations
    nSnp = size(base.hap)[1]

    # Generation 0/base: (for simplicity, I copy base SNP into `snp_prd` later)
    n₀    = size(base.hap)[2] ÷ 2   # base population size
    fmr   = par.nDam / par.nSire      # female:male ratio
    n₀Sir = Int(ceil(n₀ / (1 + fmr))) # valid even when nSires > nDams
    n₀Dam = n₀ - n₀Sir
    g₀Sex = shuffle([ones(Int8, n₀Sir); zeros(Int8, n₀Dam)])
    g₀Sir = (1:n₀)[g₀Sex .== 1]
    g₀Dam = (1:n₀)[g₀Sex .== 0]

    # Generation 2:end
    n₂Fam = max(par.nDam, par.nSire) # in generation 2:end
    # - nID in each of generation 2:end
    n₂Prd, n₂Clg = n₂Fam .* [par.nSib - par.nC7e, par.nC7e]

    # Generation 1: generate similar number of offspring in 2+ generation
    n₁Fam = max(n₀Sir, n₀Dam)   # full sib groups from base
    # - sibship sizes in 1st generation
    n₁P8n, n₁C7e = Int.(ceil.([n₂Prd, n₂Clg]./n₁Fam))
    # - total offspring in 1st generation
    n₁Prd, n₁Clg = [n₁P8n, n₁C7e] .* n₁Fam

    # Overall number of offspring
    nPrd, nClg = [n₁Prd + n₀, n₁Clg] + [n₂Prd, n₂Clg] .* (par.nG8n - 1)

    # storage For SNP, after the sizes are determined above
    snp = begin
        prd = BitArray(undef, nSnp, 2nPrd)
        snp₀ = view(prd, :, 1:2n₀)
        copy!(snp₀, base.hap)     # the SNP copy
        clg = BitArray(undef, nSnp, 2nClg)
        dic = Dict(:prd => prd, :clg => clg)
        (; dic...)
    end

    # For pedigree
    ped = begin
        prd = begin
            pa, ma = rand_mate(g₀Sir, g₀Dam, n₁P8n)
            DataFrame(
                g8n = Int16.([zeros(n₀);
                              ones(n₁Prd);
                              repeat(2:par.nG8n, inner = n₂Prd)]),
                id  = Int32.(1:nPrd),
                sex = [g₀Sex; rand(Int8.(0:1), nPrd - n₀)],
                sir = Int32.([zeros(Int32, n₀); pa; zeros(Int32, nPrd-n₀-n₁Prd)]),
                dam = Int32.([zeros(Int32, n₀); ma; zeros(Int32, nPrd-n₀-n₁Prd)]),
                tbv = zeros(Float32, nPrd), # TBV for the production trait
                p7e = zeros(Float32, nPrd), # production phenotypes
                val = zeros(Float32, nPrd)  # selection index
            )
        end
        clg = begin
            pa, ma = rand_mate(g₀Sir, g₀Dam, n₁C7e)
            DataFrame(
                g8n = Int16.([ones(n₁Clg); repeat(2:par.nG8n, inner = n₂Clg)]),
                id  = Int32.(1:nClg),
                sir = Int32.([pa; zeros(Int32, nClg - n₁Clg)]),
                dam = Int32.([ma; zeros(Int32, nClg - n₁Clg)]),
                tbv = zeros(Float32, nClg),
                p7e = zeros(Float32, nClg)
            )
        end
        dic = Dict(:prd => prd, :clg => clg)
        (; dic...)
    end

    prd = @view ped.prd[(ped.prd.g8n .== 0), :] # generation 0
    prd[:, :tbv], prd[:, :p7e] = phenotype(base.hap, qtl[1], par.h²[1])
    prd = @view ped.prd[(ped.prd.g8n .== 1), :]

    clg = @view ped.clg[(ped.clg.g8n .== 1), :]
    gene_drop(prd, base.r, snp.prd, snp.prd, qtl[1], par.h²[1])
    gene_drop(clg, base.r, snp.prd, snp.clg, qtl[2], par.h²[2])
    par.b4y && toBinary(clg, par.p17h)
    ped, snp
end

"""
    function sample_alleles(r)
---
This function returns a vector of 1:2 of the same size of `r`,
which is the vector recombination rates.
`r` for the first locus of each chromosome is 0.5,
such that they are independent.
"""
function sample_alleles(r)
    n = length(r)
    t = rand(n)
    v = ones(Int, n)
    k = 1
    @inbounds for i in 1:n
        (t[i] < r[i]) && (k = 3-k)
        v[i] = k
    end
    v
end

"""
    function select_nuclear(df, nSire, nDam)
---
This is to select `nSire`, and `nDam` from `df` on df.val.
DataFrame `df` has 3 columns:
1. ID number
2. sex, 0 for female, 1 for male.
3. val

This is constructed from a evaluation procedure,
where a method can be chosen to calculate `df.val`.
Be careful that no range check is performed.
"""
function select_nuclear(df, nSire, nDam)
    gp = groupby(df, :sex)
    sires = last(sort(gp[(sex=1,)], :val), nSire).id
    dams  = last(sort(gp[(sex=0,)], :val), nDam ).id
    return sort(sires), sort(dams)
end

"""
    function calc_idx(ped, snp, cur, lqtl, par)
---
Calculate selection index SNP and phenotype information up to current generation.
Store the results in `cur[:, :val]`.

- Generations are treated as a fixed effect in the challenge group.
- If `lqtl` is not empty, then its specified loci are fixed as fixed effects
- `twop` are those in the first generation, which were stored in `base`.
"""
function calc_idx(ped, snp, cur, Q, par)
    # About the production trait
    til = last(cur).id
    fra = par.u38y ?
        cur.id[1] - 1 :
        findlast(x -> x == 0, ped.prd.g8n)
    gp = alleles2gt(view(snp.prd, :, 2fra+1:2til))
    mp, sp = snp_blup(gp,
                      ped.prd[:, :p7e][fra+1:til],
                      par.h²[1],
                      dd = par.dd)

    # About the binary trait
    til = findlast(x -> x == first(cur).g8n, ped.clg.g8n)
    if par.u38y
        fra = findfirst(x -> x == first(cur).g8n, ped.clg.g8n)
        F = []                  # no need to fit generation as a fixed effect
    else
        fra = 1
        F = ped.clg.g8n[1:til]
    end
    gc = alleles2gt(view(snp.clg, :, 2fra-1:2til))
    mc, sc = snp_blup(gc,
                      ped.clg[:, :p7e][fra:til],
                      par.h²[2],
                      Q = Q,    # whether to emphsize known QTL
                      F = F,    # whether have more than 1 generations.
                      dd = par.dd)

    # Assign `:val` for selection
    tid = size(gp)[2]
    cgp = view(gp, :, tid-nrow(cur)+1:tid) # genotypes of current generation
    cur[:, :val] = cgp'sp + cgp'sc .* par.w4t
end

"""
    function assign_pama(ped, sires, dams, nsib)
---
Randomly mate `sires` and `dams` according to the sibship sizes `nsib`,
and assign them into the parent columns of `ped`.
"""
function assign_pama(ped, sires, dams, nsib)
    pa, ma = rand_mate(sires, dams, nsib)
    ped[:, :sir] = pa
    ped[:, :dam] = ma
end

"""
    function breeding(ped, snp, qtl, r, par)
---
The breeding program.
The `base`, which is small, is expanded to generation 1, with geno- and phenotypes.
This creates a more general/common base for later breeding with different methods.
"""
function breeding(ped, snp, base, qtl, par)
    rkq = rank_QTL(base.hap, qtl[2]) # QTL for binary trait ranking in ↓ order
    Q = par.fix ? abs.(rkq[1:par.nK3n]) : Int[]
    prd, clg = groupby(ped.prd, :g8n), groupby(ped.clg, :g8n) # shorthands

    # note, there are `nG8n+1` generations, including the base.
    @info "Evaluating generation 1"
    calc_idx(ped, snp, prd[2], Q, par) # evaluation of generation one
    for ig in 2:par.nG8n
        @info "Breeding generation $ig of $(par.nG8n)"
        sires, dams = select_nuclear(prd[ig], par.nSire, par.nDam)
        assign_pama(prd[ig+1], sires, dams, par.nSib - par.nC7e)
        elc = par.edit ? popfirst!(rkq) : 0
        gene_drop(prd[ig+1], base.r, snp.prd, snp.prd, qtl[1], par.h²[1],
                    (lc = elc, sr = par.e19e))
        assign_pama(clg[ig], sires, dams, par.nC7e)
        gene_drop(clg[ig], base.r, snp.prd, snp.clg, qtl[2], par.h²[2])
        par.b4y && toBinary(clg[ig], par.p17h) # convert to binary phenotype
        calc_idx(ped, snp, prd[ig+1], Q, par)
    end
end

"""
    function simple_breeding(ped, snp, qtl, par, op)
---
To select on `TBV`, `phenotypes`, or `EBV` on the production trait.
"""
function simple_breeding(ped, snp, base, qtl, par, op)
    prd, clg = groupby(ped.prd, :g8n), groupby(ped.clg, :g8n) # shorthands
    for ig in 2:par.nG8n
        @info "Breeding generation $ig of $(par.nG8n)"
        cur = prd[ig+1]
        sires, dams = select_nuclear(prd[ig], par.nSire, par.nDam)
        assign_pama(cur, sires, dams, par.nSib - par.nC7e)
        gene_drop(cur, base.r, snp.prd, snp.prd, qtl[1], par.h²[1])
        if op == 1
            cur[:, :val] = cur.tbv
        elseif op == 2
            cur[:, :val] = cur.p7e
        else
            fra = findlast(x -> x == 0, ped.prd.g8n)
            til = last(cur).id
            gp = alleles2gt(view(snp.prd, :, 2fra+1:2til))
            mp, sp = snp_blup(gp, ped.prd[:, :p7e][fra+1:til],
                              par.h²[1], dd = par.dd)
            tid = size(gp)[2]
            cgp = view(gp, :, tid-nrow(cur)+1:tid) # genotypes of current g8n
            cur[:, :val] = cgp'sp
        end
    end
end

"""
    function which_to_edit(snp, qtl, n)
---
Check which SNP in `snp` are segregating.
Return the biggest `n` segregating QTL according to `qtl`.
Ranking is on subsitution effects.
"""
function which_to_edit(snp, qtl, n)
    @info "Under construction"
end

"""
    function qtl_breeding(base, qtl)
---
## Intrduction
This function is to duplicate simulations ofJenko(2015).
It studied on only one trait.
Only QTL were used for SNP genotypes, as if we know all the QTL.
The breeding is then only on these QTL.
No BLUP evaluation was evolved.

## Notes
- To use the `DataFrame` structure, a dummy column `:val` is copy of `:tbv`.
- Summarize `:tbv` in the pedigree dataframe for the genetic improvements.
- n_male ≤ n_females always
- n-females should be even
- `eSires` is a vector indicating which sire to edit
"""
function qtl_breeding(base, par, qtl; eSires=[])
    par.nSire > par.nDam && error("nSire > nDam")
    par.nDam ÷ 2 * 2 == par.nDam || error("nDam not even")
    haskey(par, :nEdit) && println("Debugging")

    snp = begin                 # SNP/QTL of alternative generations
        a = copy(base.hap)
        b = copy(base.hap)
        (a, b)                  # no matter who is first
    end
    cg = 1              # the current generation. 3 - cg → alternative
    ped = begin
        nid = par.nDam * par.nSib # Nid per generation
        tid = par.nG8n * nid      # Nid all together
        g8n = repeat(1:par.nG8n, inner = nid)
        id  = repeat(1:nid, outer = par.nG8n)
        sex = zeros(Int, tid) # needs to be half male/female in every generation
        mdam = max(par.nDam, nid ÷ 2) # to garantee to have enough dams.
        for i in 0:par.nG8n-1
            sex[i*nid+1:i*nid+nid] = shuffle([zeros(Int, mdam); ones(Int, nid - mdam)])
        end
        tmp = DataFrame(g8n = g8n,
                        id = id,
                        sex = sex,
                        sir = zeros(Int, tid),
                        dam = zeros(Int, tid),
                        tbv = zeros(tid),
                        p7e = zeros(tid),
                        val = zeros(tid) # a duplicate of tbv to use codes before
                        )
        groupby(tmp, :g8n)
    end
    # TBV of the first generation
    ped[1].tbv[:] = begin
        gt = alleles2gt(snp[1])
        gt'qtl.effect
    end
    print("Breeding of generation: ")
    for i in 2:par.nG8n
        print(" $i")
        ped[i-1].val[:] = ped[i-1].tbv
        sires, dams = select_nuclear(ped[i-1], par.nSire, par.nDam)
        assign_pama(ped[i], sires, dams, par.nSib)
        gene_drop(ped[i], base.r, snp[cg], snp[3-cg], qtl, par.h²)
        cg = 3 - cg             # change to the alternative generation
    end
    println()
    ped, snp[cg]                # for base of the next simulation
end


text:
module TestT3FF

include("other/FEMMShellT3DSGMTModule.jl")

end # module


text:
using DFTK
using MPI
using JSON
using JLD2
using Printf
using Plots
using Statistics
using LinearAlgebra
using Dates
using Logging
include("quantum_espresso.jl")

function build_magnetic_moments(atoms::AbstractArray; magmoms...)
    magmoms = Dict(magmoms)
    map(atoms) do (element, positions)
        element => fill(magmoms[element.symbol], length(positions))
    end
end

function attach_pseudos(atoms::AbstractArray; pseudomap...)
    pseudomap = Dict(pseudomap)

    map(atoms) do (element, position)
        pspfile = get(pseudomap, element.symbol, nothing)
        ElementPsp(element.symbol, psp=load_psp(pspfile)) => position
    end
end

function standard_cases()
    cases = Dict{String, Any}(
        "adaptive"    => (damping=DFTK.AdaptiveDamping(), ),  # Parameter-free adaptive damping
        "adaptive_05" => (damping=DFTK.AdaptiveDamping(0.5), ),
    )
    for α in 1:9
        cases["standard_0$α"] = (damping=DFTK.FixedDamping(α / 10), )
    end
    cases["standard_10"] = (damping=DFTK.FixedDamping(1.0), )
    cases
end

function random_logfile(prefix)
    logfile = "$(prefix)_$(@sprintf "%05d" rand(UInt16)).log"
    while isfile(logfile)
        logfile = "$(prefix)_$(@sprintf "%05d" rand(UInt16)).log"
    end
    logfile
end

function run_step_on_cases(payload, cases)
    spay = string(payload)
    header_printed = false

    anymissing = false
    for dir in sort(collect(keys(cases)))
        config = cases[dir]

        MPI.Barrier(MPI.COMM_WORLD)
        if isfile(joinpath(config.directory, string(payload) * ".done"))
            continue
        end

        for i in 1:10  # To avoid infinite loop
            MPI.Barrier(MPI.COMM_WORLD)
            res = nothing
            logfile = nothing
            if mpi_master()
                if !header_printed
                    println("#####" * "#"^length(spay) * "#####")
                    println("#--  $(string(payload))  --#")
                    println("#####" * "#"^length(spay) * "#####")
                    header_printed = true
                end

                logfile = random_logfile(joinpath(config.directory, string(payload)))
                println(basename(dir), "  ->  ", basename(dir) * "/" * basename(logfile))
                res = open(logfile, "w") do fp
                    redirect_stdout(fp) do
                        with_logger(SimpleLogger(stdout, Logging.Debug)) do
                            payload(config)
                        end
                    end
                end
                println("    ... $res\n")
            else
                res = payload(config)
            end

            if res == :exhausted
                if mpi_master()
                    rm(logfile)
                    anymissing || touch(joinpath(config.directory, string(payload) * ".done"))
                end
                break
            elseif res == :missing
                anymissing = true
                mpi_master() && rm(logfile)
            end
        end
    end

    MPI.Barrier(MPI.COMM_WORLD)
    nothing
end

function has_results(fileprefix)
    resultsfile = fileprefix * ".scfres.json"
    isfile(resultsfile) && open(JSON.parse, resultsfile, "r")["done"]
end

function basis_from_config(config)
    model = model_PBE(config.lattice, config.atoms;
                      magnetic_moments=config.magnetic_moments, config.kwargs_model...)
    fft_size = compute_fft_size(model, config.Ecut, supersampling=config.supersampling)
    if !haskey(config, :kgrid)
        kgrid = kgrid_size_from_minimal_spacing(config.lattice, 2π * 0.022)
    else
        kgrid=config.kgrid
    end
    PlaneWaveBasis(model, config.Ecut, kgrid=kgrid, kshift=[0, 0, 0], fft_size=fft_size)
end

function scf(config)
    basis = basis_from_config(config)

    scfkey = nothing
    if mpi_master()
        for scf in sort(collect(keys(config.cases)))
            lockfile = joinpath(config.directory, scf * ".running")
            fileprefix = joinpath(config.directory, scf)
            (isfile(lockfile) || has_results(fileprefix))  && continue
            write(lockfile, read(`hostname`, String))
            scfkey = scf
            break
        end
    end
    scfkey = MPI.bcast(scfkey, 0, MPI.COMM_WORLD)

    scfkey === nothing && return :exhausted
    fileprefix  = joinpath(config.directory, scfkey)
    lockfile    = joinpath(config.directory, scfkey * ".running")
    resultsfile = joinpath(config.directory, scfkey * ".scfres.json")
    guessfile   = joinpath(config.directory, "init.jld2")

    !isfile(guessfile) && return :missing
    ρ, ψ0, V = load_guess(basis, guessfile)

    try
        if mpi_master()
            println("#\n#-- SCF $(basename(config.directory)) -- case=$scfkey\n#")
            println("    keywords: $((; config.cases[scfkey]..., config.kwargs_scf...))")
        end

        MPI.Barrier(MPI.COMM_WORLD)
        mpi_master() && DFTK.reset_timer!(DFTK.timer)
        if haskey(config.cases[scfkey], :type) && config.cases[scfkey].type == :qe
            caseargs = Dict(k => v for (k, v) in pairs(config.cases[scfkey]) if k != :type)
            run_qe(basis, fileprefix, resultsfile;
                   magnetic_moments=config.magnetic_moments, caseargs...)
        else
            run_dftk(basis, resultsfile; ρ=ρ, ψ0=ψ0, V=V,
                     mixing=config.mixing, config.cases[scfkey]..., config.kwargs_scf...)
        end
        MPI.Barrier(MPI.COMM_WORLD)
        mpi_master() && println(DFTK.timer)
    catch
        mpi_master() && rm(joinpath(config.directory, "scf.done"), force=true)
        rethrow()
    finally
        mpi_master() && rm(lockfile)
    end
    mpi_master() && plot_convergence_summary(config)

    :done
end

function dump_guess(config)
    basis     = basis_from_config(config)
    guessfile = joinpath(config.directory, "init.jld2")
    isfile(guessfile) && return :exhausted

    ρ = guess_density(basis, config.magnetic_moments)
    _, ham = energy_hamiltonian(basis, nothing, nothing; ρ=ρ)

    V = nothing
    isnothing(config.potential) || (V = config.potential(basis))
    if !isnothing(V)
        ham = DFTK.hamiltonian_with_total_potential(ham, V)
    end

    n_bands = (haskey(config, :n_bands) ? (n_bands = config.n_bands)
                                        : DFTK.default_n_bands(basis.model))
    res_ψ0 = DFTK.next_density(ham; n_bands=DFTK.default_n_bands(basis.model),
                               n_ep_extra=3, miniter=2, tol=0.0001)
    diagiter = DFTK.mpi_mean(mean(res_ψ0.diagonalization.iterations), basis.comm_kpts)
    mpi_master() && println("guess:    diag = $diagiter  diagtol = 0.0001")
    isnothing(V) && (V = DFTK.total_local_potential(ham))

    # Gather ψ array on master
    ψ = DFTK.gather_kpts(res_ψ0.ψ, basis)
    if mpi_master()
        JLD2.jldopen(guessfile, "w") do jld
            jld["ρ"] = ρ
            jld["V"] = V
            jld["ψ"] = ψ
        end
    end

    :done
end

function load_guess(basis, guessfile)
    JLD2.jldopen(guessfile, "r") do jld
        ψ = jld["ψ"][basis.krange_thisproc]
        jld["ρ"], ψ, jld["V"]
    end
end

function run_dftk(basis, resultsfile;
                  ρ, ψ0, V, mixing=SimpleMixing(),
                  damping, tol=1e-11, type=:potmix, diagtol_max=1e-3, kwargs...)

    n_bands = (haskey(config, :n_bands) ? (n_bands = config.n_bands)
                                        : DFTK.default_n_bands(basis.model))
    if mpi_master()
        println()
        println("hostname        = $(read(`hostname`, String))")
        println("started on      = $(Dates.now())")
        println("julia threads   = $(Threads.nthreads())")
        println("BLAS threads    = $(BLAS.get_num_threads())")
        println("MPI procs       = $(mpi_nprocs())")
        println()
        println("temperature     = $(basis.model.temperature)")
        println("smearing        = $(basis.model.smearing)")
        println("lattice         = $(round.(basis.model.lattice, sigdigits=4))")
        println("Ecut            = $(basis.Ecut)")
        println("fft_size        = $(basis.fft_size)")
        println("kgrid           = $(basis.kgrid)")
        println("kshift          = $(basis.kshift)")
        println("n_irreducible_k = $(sum(length, basis.krange_allprocs))")
        println("n_bands         = $(n_bands)")
        println("n_electrons     = $(basis.model.n_electrons)")
        println("type            = $type")
        println("mixing          = $mixing")
        println("damping         = $damping")
        flush(stdout)
    end

    # Collection of data to dump later
    save_dict = Dict(
        "temperature"  => basis.model.temperature,
        "smearing"     => string(basis.model.smearing),
        "lattice"      => basis.model.lattice,
        "Ecut"         => basis.Ecut,
        "fft_size"     => basis.fft_size,
        "kgrid"        => basis.kgrid,
        "n_kirred"     => sum(length, basis.krange_allprocs),
        "n_bands"      => n_bands,
        "n_electr"     => basis.model.n_electrons,
        "mixing"       => string(mixing),
        "damping"      => string(damping),
        "type"         => string(type),
        "kwargs"       => string(kwargs),
        "done"         => false,
        #
        "energies"     => Float64[],
        "residuals"    => Float64[],
        "ndiag"        => Int[],
        "α"            => Float64[],
    )

    # Empty existing results
    if mpi_master()
        open(fp -> JSON.print(fp, save_dict), resultsfile, "w")
    end

    function ExtractCallback()
        function callback(info)
            if info.stage != :finalize
                push!(save_dict["energies"],  info.energies.total)
                push!(save_dict["ndiag"],     length(info.diagonalization))
                push!(save_dict["α"], isnan(info.α) ? DFTK.trial_damping(damping) : info.α)

                if type == :potmix
                    # Note: This misses a factor sqrt(info.basis.dvol)), which is taken into
                    # account in the analyse.jl script
                    push!(save_dict["residuals"], norm(info.Vout - info.Vin))
                elseif type == :densmix
                    # Note: This misses a factor sqrt(info.basis.dvol)), which is taken into
                    # account in the analyse.jl script
                    push!(save_dict["residuals"], norm(info.ρout - info.ρin))
                end
            end

            if mpi_master()
                open(fp -> JSON.print(fp, save_dict), resultsfile, "w")
            end
        end
    end

    determine_diagtol = DFTK.ScfDiagtol(; ratio_ρdiff=0.01, diagtol_max, diagtol_min=1e-12)
    callback = ExtractCallback() ∘ DFTK.ScfDefaultCallback()

    if type == :potmix
        run_scf = ((damping isa DFTK.AdaptiveDamping) ? DFTK.scf_potential_mixing_adaptive
                                                      : DFTK.scf_potential_mixing)
        scfres = run_scf(basis; tol=tol, ψ=ψ0, determine_diagtol=determine_diagtol,
                         mixing=mixing, callback=callback, damping=damping, ρ=ρ, V=V,
                         diag_miniter=2, kwargs...)
    elseif type == :densmix
        @assert damping isa Number
        scfres = self_consistent_field(basis; tol=tol, ψ=ψ0,
                                       determine_diagtol=determine_diagtol,
                                       mixing=mixing, callback=callback,
                                       damping=damping, ρ=ρ, kwargs...)
    else
        error("Unsupported SCF type: $type")
    end

    if mpi_master()
        println()
        println(scfres.energies)
        flush(stdout)
        save_dict["done"] = true
        open(fp -> JSON.print(fp, save_dict), resultsfile, "w")
    end

    scfres
end

function run_qe(basis, qe_directory, resultsfile; magnetic_moments,
                mixing_beta=0.7, kwargs...)
    if !mpi_master()
        MPI.Barrier(MPI.COMM_WORLD)
        return nothing
    end
    # This only runs on MPI master

    println()
    println("hostname        = $(read(`hostname`, String))")
    println("started on      = $(Dates.now())")
    println()
    flush(stdout)

    model = basis.model
    outputfile = joinpath(qe_directory, "espresso.pwo")
    if !isfile(outputfile)
        qeres = scf_quantum_espresso(model.lattice, model.atoms;
                                     Ecut=basis.Ecut,
                                     fileprefix=joinpath(qe_directory, "espresso"),
                                     temperature=model.temperature,
                                     mixing_beta=mixing_beta,
                                     kgrid=basis.kgrid,
                                     magnetic_moments=magnetic_moments,
                                     kwargs...)
    @assert qeres.outputfile == outputfile
    end

    # Copy output file to stdout
    println.(readlines(outputfile))
    flush(stdout)

    # Extract results
    save_dict = Dict(
        "done"         => true,
        "energies"     => Float64[],
        "residuals"    => Float64[],
        "ndiag"        => Int[],
        "α"            => Float64[],
    )
    for iteration in parse_quantum_espresso_pwo(outputfile)
        push!(save_dict["energies"],  iteration.total_energy)
        push!(save_dict["residuals"], iteration.accuracy)
        push!(save_dict["ndiag"],     1)
        push!(save_dict["α"],         mixing_beta)
    end
    open(fp -> JSON.print(fp, save_dict), resultsfile, "w")

    MPI.Barrier(MPI.COMM_WORLD)
    nothing
end

function collect_minimal_total_energy(config::NamedTuple; kwargs...)
    collect_minimal_total_energy(config.directory, config.cases; kwargs...)
end
function collect_minimal_total_energy(directory::AbstractString, cases::AbstractDict;
                                      kwargs...)
    collect_minimal_total_energy(directory, collect(keys(cases)); kwargs...)
end
function collect_minimal_total_energy(directory::AbstractString, cases::AbstractVector;
                                      tol_resid=0.05, tol_energy=1e-4)
    mintotals = Float64[]
    for scfkey in sort(cases)
        startswith(scfkey, "espresso_") && continue
        resultfile = joinpath(directory, scfkey * ".scfres.json")
        isfile(resultfile) || continue

        residuals = open(JSON.parse, resultfile, "r")["residuals"]
        energies  = open(JSON.parse, resultfile, "r")["energies"]
        energies === nothing && continue
        isempty(energies) && continue

        minidx = argmin(energies[max(end-4, 1):end]) + max(length(energies) - 4, 1) - 1
        residuals[minidx] > tol_resid && continue

        min_matching = findfirst(e -> abs(e - energies[minidx]) < tol_energy, mintotals)
        if min_matching === nothing
            push!(mintotals, energies[minidx])
        else
            mintotals[min_matching] = min(mintotals[min_matching], energies[minidx])
        end
    end

    mintotals
end

function parse_scfkey_label(scfkey)
    if scfkey == "adaptive"
        ("adaptive", 0.2)
    elseif startswith(scfkey, "adaptive_") && all(isnumeric, scfkey[10:end])
        digits = length(scfkey[10:end]) - 1
        alpha  = parse(Int, scfkey[10:end]) / 10^digits
        ("adaptive", alpha)
    elseif startswith(scfkey, "adaptive_")
        ("adaptive", scfkey[10:end])
    elseif startswith(scfkey, "standard_") && all(isnumeric, scfkey[10:end])
        digits = length(scfkey[10:end]) - 1
        alpha  = parse(Int, scfkey[10:end]) / 10^digits
        ("fixed", alpha)
    elseif startswith(scfkey, "espresso_") && all(isnumeric, scfkey[10:end])
        digits = length(scfkey[10:end]) - 1
        alpha  = parse(Int, scfkey[10:end]) / 10^digits
        ("espresso", alpha)
    elseif startswith(scfkey, "highprec_") && all(isnumeric, scfkey[10:end])
        digits = length(scfkey[10:end]) - 1
        alpha  = parse(Int, scfkey[10:end]) / 10^digits
        ("highprec", alpha)
    else
        (scfkey, nothing)
    end
end

function plot_convergence_summary(config; iteration=Inf, tol=1e-10, damping=true,
                                  outputfile=joinpath(config.directory, "summary.svg"),
                                  scfkeys=nothing, mintotals=nothing, kwargs...)
    function deduce_label(scfkey)
        method, alpha = parse_scfkey_label(scfkey)
        method * (isnothing(alpha) ? "" : " (α = $alpha)")
    end

    !mpi_master() && error("Plotting only on MPI master supported.")
    if isnothing(mintotals)
        mintotals = collect_minimal_total_energy(config)
    end
    isempty(mintotals) && return plot()  # Return on no data

    if scfkeys === nothing
        scfkeys = collect(keys(config.cases))
    end

    errors   = Dict{String, Vector{Float64}}()
    cumdiags = Dict{String, Vector{Int}}()
    dampings = Dict{String, Vector{Float64}}()
    whichmin = Dict{String, Int}()
    for scfkey in sort(scfkeys)
        method, _ = parse_scfkey_label(scfkey)
        resultfile = joinpath(config.directory, scfkey * ".scfres.json")
        isfile(resultfile) || continue

        energies  = open(JSON.parse, resultfile, "r")["energies"]
        ndiags    = open(JSON.parse, resultfile, "r")["ndiag"]
        damps     = open(JSON.parse, resultfile, "r")["α"]

        energies === nothing && continue
        if method == "espresso"
            values = open(JSON.parse, resultfile, "r")["residuals"]
            mintrial  = minimum(energies[max(end-4, 1):end])
            bestmatch = sortperm(abs.(mintrial .- mintotals))[1]
        else
            isempty(energies) && continue
            # Find best matching minimum
            mintrial  = minimum(energies[max(end-4, 1):end])
            bestmatch = sortperm(abs.(mintrial .- mintotals))[1]
            values = abs.(energies .- mintotals[bestmatch])
            idx = findfirst(values .< 1e-16)
            if !isnothing(idx)
                values = values[1:idx-1]
                ndiags = ndiags[1:idx-1]
                damps  = damps[1:idx-1]
            end
        end
        isempty(values) && continue

        whichmin[scfkey] = bestmatch
        errors[scfkey]   = values
        dampings[scfkey] = damps
        cumdiags[scfkey] = cumsum(ndiags)
    end

    # If no data return
    isempty(errors) && return plot()

    any_adaptive = false
    ymax = 0
    maxiter = 0
    colors = Dict{Any, Int}()
    p = plot(; yaxis=:log, legend=:topright, title=basename(config.directory), kwargs...)
    for (i, scfkey) in enumerate(sort(collect(keys(errors))))
        label = deduce_label(scfkey)
        maxiter = max(maxiter, length(errors[scfkey]))
        ymax = max(ymax, maximum(errors[scfkey]))

        if length(mintotals) > 1
            label *= " (min $(whichmin[scfkey]))"
        end

        method, alpha = parse_scfkey_label(scfkey)
        if isnothing(alpha)
            color = i
        elseif alpha in keys(colors)
            color = colors[alpha]
        else
            colors[alpha] = i
            color = colors[alpha]
        end

        mark = :x
        if method == "adaptive"
            mark = :o
            any_adaptive = true
        elseif method == "espresso"
            mark = :rtriangle
        elseif method == "highprec"
            mark = :ltriangle
        end
        plot!(p, cumdiags[scfkey], errors[scfkey];
              label=label, color=color, mark=mark, lw=1.3, markersize=4)
    end
    ylims!(p, (tol, 2 * ymax))
    ylabel!(p, "Total energy absolute error")

    if damping && any_adaptive
        q = twinx()
        for scfkey in sort(collect(keys(errors)))
            method, alpha = parse_scfkey_label(scfkey)
            method == "adaptive" || continue
            (alpha in keys(colors)) || continue
            plot!(q, cumdiags[scfkey], abs.(dampings[scfkey]), mark=:+,
                  color=colors[alpha], label="", linestyle=:dot)
        end
        ylims!(q, (0.0, 1.0))
        ylabel!(q, "Damping parameter")
    end
    xlims!(p, (0, min(iteration, maxiter)))
    xlabel!(p, "Iteration / Hamiltonian diagonalization")

    if isnothing(outputfile)
        p
    else
        savefig(p, outputfile)
    end
end

max_stars_repo_name: mfherbst/supporting-adaptive-damping

max_issues_repo_name: mfherbst/supporting-adaptive-damping

max_forks_repo_name: mfherbst/supporting-adaptive-damping


text:
#CORRELATION AND SIMULATION TOOLS FOR MC HAMMER
#by Eric Torkia, April 17th 2019

# using Distributions
# using StatsBase
# using DataFrames
# using LinearAlgebra
# using Random

function cormat(ArrayName, RankOrder)
cor_mat = []
M_Size = size(ArrayName,2)
if RankOrder == 1
      for i=1:M_Size
      cor_vector = []
            for i2=1:M_Size
             cor_i = cor(tiedrank(ArrayName[i]),tiedrank(ArrayName[i2]));
             push!(cor_vector,cor_i)
             #print(cov_i)
            end
       push!(cor_mat,cor_vector)
      end
else
      for i=1:M_Size
      cor_vector = []
            for i2=1:M_Size
             cor_i = cor(ArrayName[i], ArrayName[i2]);
             push!(cor_vector,cor_i)
             #print(cov_i)
      end
      push!(cor_mat,cor_vector)
      end
end

return hcat(cor_mat...)
end

function covmat(ArrayName)
M_Size = size(ArrayName,2)
cov_mat = []
      for i=1:M_Size
            for i2=1:M_Size
             cov_i = cov(ArrayName[i],ArrayName[i2]);
             push!(cov_mat,cov_i)
             #print(cov_i)
            end
      end
return reshape(cov_mat, M_Size, M_Size)
end

#corvar function correlates simulation inputs unsing the Iman Connover Method
function corvar(ArrayName, n_trials, correl_matrix)

#Define how many columns of ISNs are required
array_dims = length(ArrayName)

#calc cholesky transform to transfer correlation to ISNs
P = cholesky(correl_matrix)

#Normal() Returns Standard Normals (ISNs)
R = rand(Normal(),n_trials,array_dims)
ISN_Matrix = R*P.U
ISN_Matrix_DF = DataFrame(ISN_Matrix)

#apply ranks to create independant correlation rankings matrix
ISN_Ranked = []
for i = 1:array_dims
      temp_ranks = ordinalrank(ISN_Matrix_DF[i])
      push!(ISN_Ranked, temp_ranks)
end
ISN_Ranked_DF = DataFrame(ISN_Ranked)
#Results_Ranked = []
#println(ISN_Ranked_DF)

#Reindex the array of samples using the ISN_Ranks. Sort(Array)[OrderingVector]
final_array=[]
for i = 1:array_dims
      sorted_array = sort(ArrayName[i])[ISN_Ranked[i]]
      push!(final_array, sorted_array)
end
Final_DF=DataFrame(final_array)
end

function GetCertainty(ArrayName, x, AboveBelow)
      if AboveBelow == 1
            certainty = count(i ->(i>=x), ArrayName)/length(ArrayName)
      else
            certainty = count(i ->(i<=x), ArrayName)/length(ArrayName)
      end
      return certainty
end

#Shell /Dos Command wrapper
cmd(x) = run(`cmd /C $x`)

function P_Buffet(ArrayName, Increment=0.1)
names = map(i -> "P" * string(i), RoundDown(collect(0:Increment:1)*100,0))
P_Vals = collect(0:Increment:1)
#convert(Array{Float64}, ArrayName)
fractiles = quantile(collect(Float64, ArrayName), P_Vals)
return hcat(names, fractiles)
end

#convert(Array{Float64}, ArrayName)
# convert(Array{Int64}, collect(0:Increment:1)*100)

#Testing Script
#--------------------------
#Must setup as vectors using commas in the []

#Final_Ranks = []
# for i = 1:array_dims
#       temp_ranks = ordinalrank(Final_DF[i])
#       push!(Final_Ranks, temp_ranks)
# end
# Final_Ranks_DF = DataFrame(Final_Ranks)

#Useful Variables
#Phi = (sqrt(5)+1)/2
# n_trials = 1000
# correl_matrix = [1 0 0 0 0 0; 0 1 0 0 0 0; 0 0 1 0 0 0;0 0 0 1 0.75 -0.7; 0 0 0 0.75 1 -0.95; 0 0 0 -0.7 -0.95 1 ]
#
# #Test Data
# test_logNormal = [rand(LogNormal(0, 0.5),n_trials), rand(Normal(3,2),n_trials), rand(Gamma(1, 0.5),n_trials), rand(LogNormal(0, 0.5),n_trials), rand(Normal(3,2),n_trials), rand(Gamma(1, 0.5),n_trials)]


text:
"Base type for all parameters"
abstract Parameters

export linearize_parameters!

# --- Parameter types

"Layer parameters - contains the values, the gradient and optimization specific information"
type ArrayParameters{D<:Device, F<:Float, N} <: Parameters
  # The parameters
  values::denseTensor(D,F,N)

  # Gradient, or null if the gradient
  gradient::Nullable{denseTensor(D,F,N)}

  # FIXME: remove? should not be part of operators?
  optimization_state::Nullable{Any} # Used by an optimization method to store any parameter related information

  ArrayParameters(values, gradient) = new(values, gradient, Nullable())

  function ArrayParameters(d::D, ::Type{F}, dims::Int64...)
    @assert length(dims) == N
    new(array(d, F, dims...), Nullable(array(d, F, dims...)), Nullable())
  end
end

typealias VectorParameters{D<:Device, F<:Float} ArrayParameters{D, F, 1}
typealias MatrixParameters{D<:Device, F<:Float} ArrayParameters{D, F, 2}

matrixParameters(D::TypeVar, F::TypeVar) = ArrayParameters{D,F,2}
matrixParameters{D<:Device, F<:Float}(::Type{D}, ::Type{F}) = ArrayParameters{D,F,2}

vectorParameters(D::TypeVar, F::TypeVar) = ArrayParameters{D,F,1}
vectorParameters{D<:Device, F<:Float}(::Type{D}, ::Type{F}) = ArrayParameters{D,F,1}

arrayParameters{D<:Device, F<:Float}(d::D, ::Type{F}, dims::Int64...) =
  ArrayParameters{D,F,length(dims)}(d, F, dims...)

matrixParameters{D<:Device, F<:Float}(d::D, ::Type{F}, dims::Int64...) =
  ArrayParameters{D,F,2}(d, F, dims...)

vectorParameters{D<:Device, F<:Float}(d::D, ::Type{F}, dims::Int64...) =
  ArrayParameters{D,F,1}(d, F, dims...)

length(m::ArrayParameters) = length(m.values)
size(m::ArrayParameters) = size(m.values)

# Caches the parameter fields for each type
const parametersMap = Dict{Type, Array{Symbol}}()

function _parameters(layer::Layer)
  t = typeof(layer)
  return get!(parametersMap, t) do
    p = Any[]
    for field in fieldnames(layer)
        if fieldtype(t, field) <: Parameters
            push!(p, field)
        end
    end
    p
  end
end

"""Returns the parameters of a given layer"""
function parameters(m::Layer)
    function _it()
        for field in _parameters(m)
            produce(m.(field))
        end
    end
    Task(_it)
end

# --- Methods for parameters

init!(p::ArrayParameters) = randn!(p.values)
init_gradient!(p::ArrayParameters) = fill!(get(p.gradient), 0.)

"""Linearize parameters and their gradients

- Parameters with null gradient are ignored
- All parameters should have the same storage and type

"""
function linearize_parameters!{F<:Float, D<:Device}(d::D, ::Type{F}, m::Layer, copy::Bool = false)
  @assert(!copy, "not implemented")

  # Compute the total size needed
  totalsize::Int = 0
  for p in parameters(m)
    if typeof(p) <: ArrayParameters && !isnull(p.gradient)
      @assert typeof(p) <: ArrayParameters{D,F}
      totalsize += length(p)
    end
  end

  # Allocate the array
  values = array(d, F, totalsize)
  gradient = array(d, F, totalsize)

  offset::Int = 0
  for p in parameters(m)
    if typeof(p) <: ArrayParameters && !isnull(p.gradient)
      psize = length(p.values)
      r = offset:(offset+psize-1)

      p.values = reshape_view(view(values, r), size(p.values))
      p.gradient = Nullable(reshape_view(view(gradient, r), size(get(p.gradient))))

      offset += psize
    end
  end

  return values, gradient
end


text:
using DataFrames
using Dates
using Base

"""
    function datatype(str::String)
This function evaluates an input string as a DataType if it is defined.
Otherwise, it will return false.
See: (thread on discourse.julialang.org)[https://discourse.julialang.org/t/parse-string-to-datatype/7118/9]
"""
function datatype(str::String)
    type = :($(Symbol(str)))
    return isdefined(SLiDE, type) ? eval(type) : nothing
end

"""
    Base.broadcastable(x::InvertedIndex{T}) where {T<:Any}
!!!!
"""
Base.broadcastable(x::InvertedIndex{T}) where {T<:Any} = [x];

"""
    Base.split(x::Missing)
    Base.split(x::Number)
Extends `split` to ignore missing fields.
"""
Base.split(str::Missing) = str
Base.split(str::Missing, splitter::Any) = str

"""
    Base.strip(x::Missing)
    Base.strip(x::Number)
Extends `strip` to ignore missing fields and numbers.
"""
# Base.strip(x::String) = replace(x, r"^\s*\"*|\"*\s*$" => "")
Base.strip(x::Missing) = x
Base.strip(x::Number) = x

"""
    Base.lowercase(x::Symbol)
Extends `lowercase` to handle other data types.
"""
Base.lowercase(x::Int) = x
Base.lowercase(x::Symbol) = Symbol(lowercase(string(x)))
Base.lowercase(x::Missing) = missing

"""
    Base.titlecase(x::Symbol)
Extends `titlecase` to handle other data types.
"""
Base.titlecase(x::Int) = x
Base.titlecase(x::Symbol) = Symbol(titlecase(string(x)))
Base.titlecase(x::Missing) = missing

"""
    Base.uppercase(x::Symbol)
Extends `uppercase` to handle other data types.
"""
Base.uppercase(x::Int) = x
Base.uppercase(x::Symbol) = Symbol(uppercase(string(x)))
Base.uppercase(x::Missing) = missing

"""
    Base.uppercasefirst(x::Symbol)
Extends `uppercasefirst` to handle other data types.
"""
Base.uppercasefirst(x::Int) = x
Base.uppercasefirst(x::Symbol) = Symbol(uppercasefirst(string(x)))
Base.uppercasefirst(x::Missing) = missing

"""
    Base.occursin(x::Symbol, y::Symbol)
    Base.occursin(x::String, y::Symbol)
Extends `occursin` to work for symbols. Potentially helpful for DataFrame columns.
"""
Base.occursin(x::Symbol, y::Symbol) = occursin(string(x), y)
Base.occursin(x::String, y::Symbol) = occursin(x, string(y))

"""
"""
function dropvalue!(df::DataFrame, x::Float64)
    cols = find_oftype(df, typeof(x));
    if isnan(x); [filter!(row -> .!isnan.(row[col]), df) for col in cols]
    else;        [filter!(row -> row[col] .!== x, df) for col in cols]
    end
    return df
end

dropvalue(df::DataFrame, x::Float64) = dropvalue!(copy(df), x)

"""
    dropzero!(df::DataFrame)
Returns a DataFrame without zero values in columns of type AbstractFloat.
"""
dropzero!(df::DataFrame) = dropvalue!(df, 0.0)
dropzero(df::DataFrame) = dropzero!(copy(df))

dropnan!(df::DataFrame) = dropvalue!(df, NaN)
dropnan(df::DataFrame) = dropnan!(copy(df))

"""
    convert_type(::Type{T}, x::Any)
    convert_type(::Dict{Any,Any}, df::DataFrame, value_col::Symbol; kwargs...)
Converts `x` into the specified `Type{xT}`.

Consider extending [convert](https://docs.julialang.org/en/v1/base/base/#Base.convert)
function (!!!!)

# Arguments
- `::Type{T}`: target DataType.
- `x<:Any`: value to convert.

# Keyword Arguments
Options available when converting a DataFrame into a dictionary of keys pointing to a value:
- `drop_cols = []`: Columns not to include in
    the keys. By default, no columns are dropped.
- `value_col::Symbol = :end`: If converting

# Returns
Data in specified type
"""
convert_type(::Type{T}, x::Any) where T<:AbstractString = string(x)
convert_type(::Type{T}, x::Date) where T<:Integer = Dates.year(x)

convert_type(::Type{Map}, x::Group) = Map(x.file, [x.from], x.to, [x.input], x.output, :inner)

convert_type(::Type{T}, x::AbstractString) where T<:AbstractString = string(strip(x))

function convert_type(::Type{T}, x::AbstractString) where T<:Integer
    return convert_type(T, convert_type(Float64, x))
end

function convert_type(::Type{T}, x::AbstractString) where T<:Real
    return parse(T, reduce(replace, ["," => "", "\"" => ""], init = x))
end

convert_type(::Type{T}, x::Symbol) where T<:Real = convert_type(T, convert_type(String, x))

convert_type(::Type{DataFrame}, lst::Array{Dict{Any,Any},1}) = [DataFrame.(lst)...;]

function convert_type(::Type{DataFrame}, arr::JuMP.Containers.DenseAxisArray; cols = [])
    cols = ensurearray(cols)

    val = JuMP.value.(arr.data)
    ind = permute(arr.axes);
    val = collect(Iterators.flatten(val));

    df = hcat(DataFrame(ensuretuple.(ind)), DataFrame([val], [:value]))
    return edit_with(df, Rename.(propertynames(df)[1:length(cols)], cols))
end

function convert_type(::Type{Dict}, df::DataFrame; drop_cols = [], value_col::Symbol = :Float)
    # Find and save the column containing values and that/those containing keys.
    # If no value column indicator is specified, find the first DataFrame column of floats.
    value_col == :Float && (value_col = find_oftype(df, AbstractFloat)[1])
    key_cols = setdiff(propertynames(df), convert_type.(Symbol, ensurearray(drop_cols)), [value_col])
    ONEKEY = length(key_cols) == 1

    d = Dict((ONEKEY ? row[key_cols[1]] : (row[key_cols]...,)) => row[value_col]
        for row in eachrow(df))
    return d
end

convert_type(::Type{DataType}, x::AbstractString) = datatype(x)

convert_type(::Type{Array{T}}, x::Any) where T<:Any = convert_type.(T, x)
convert_type(::Type{Array{T,1}}, x::Any) where T<:Any = convert_type.(T, x)
convert_type(::Type{Array}, d::Dict) = [collect(values(d))...;]

convert_type(::Type{T}, x::Missing) where T<:Real = x;
convert_type(::Type{T}, x::Missing) where T<:AbstractString = x
convert_type(::Type{Any}, x::AbstractString) = "missing" == lowercase(x) ? missing : x
convert_type(::Type{Any}, x::Any) = x

convert_type(::Type{T}, x::Any) where T = T(x)

convert_type(::Type{Bool}, x::AbstractString) = lowercase(x) == "true" ? true : false

# [@printf("%-8s %s\n", T, fieldpropertynames(T)[T.types .== Any]) for T in subtypes(Edit) if Any in T.types]

"""
    isarray(x::Any)
Returns true/false if the the DataType or object is an array.
"""
isarray(::Type{Array{T,1}}) where T <: Any = true
isarray(x::Array{T,1}) where T <: Any = true
isarray(::Any) = false

"""
    ensurearray(x::Any)
"""
ensurearray(x::Array{T,1}) where T <: Any = x
ensurearray(x::Tuple{Vararg{Any}}) = collect(x)
ensurearray(x::UnitRange) = collect(x)
ensurearray(x::Any) = [x]

"""
    ensuretuple(x::Any)
Returns `x` in a tuple.
"""
ensuretuple(x::Tuple{Vararg{Any}}) = x
ensuretuple(x::Any) = tuple(x)

"""
"""
istype(df::DataFrame, T::DataType) = broadcast(<:, eltype.(eachcol(dropmissing(df))), T)
# eltype.(eachcol(df))

function hasnames(df::DataFrame, cols::Array{Symbol,1})
    col_in = setdiff(propertynames(df),cols)
    col_out = setdiff(cols, propertynames(df))
    return (length(col_in) == 0) && (length(col_out) == 0)
end

function ensurenames!(df::DataFrame, cols::Array{Symbol,1})
    size(df,2) !== length(cols) && @error("Can only ensure column names of the data frame length")
    cols_in = setdiff(propertynames(df),cols)
    cols_out = setdiff(cols, propertynames(df))

    [rename!(df, col_in => col_out) for (col_in,col_out) in zip(cols_in,cols_out)]
    return df[:,cols]
end
ensurenames(df::DataFrame, cols::Array{Symbol,1}) = ensurenames!(copy(df), cols)

"""
    find_oftype(df::DataFrame, T::DataType)
    find_oftype(df::Dict, T::DataType)
Returns DataFrame column names of the specified type.
"""
find_oftype(df::DataFrame, T::DataType) = propertynames(df)[istype(df, T)]
find_oftype(df::DataFrame, T::InvertedIndex{DataType}) = propertynames(df)[.!istype(df, T.skip)]

function find_oftype(d::Dict, T::DataType)
    return Dict(k => v for (k,v) in d if any(broadcast(<:, typeof.(ensurearray(v)), T)))
end

"""
    permute(x::Any)
This function finds all possible permutations of the input arrays.

# Arguments
- `x::Tuple` or `x::NamedTuple{}` or `x::Array`: list of arrays to permute.
    If `x` is a NamedTuple, its values will be permuted.

# Returns
- `x::Array{Tuple,1}`: list of all possible permutations of the input values.
    If `x` does not contain at least one array, there will be nothing to permute and the function will return `x`.
"""
function permute(df::DataFrame)
    cols = propertynames(df)
    df = sort(DataFrame(ensuretuple.(permute(unique.(eachcol(df))))))
    df = edit_with(df, Rename.(propertynames(df), cols))
    return df
end

function permute(x::Tuple)
    xperm = if length(x) == 1
        sort(unique(x[1]))
    else
        [collect(Base.Iterators.product(sort.(unique.(ensurearray.(x)))...))...;]
    end
    return xperm
end

function permute(x::NamedTuple)
    cols = keys(x)
    xperm = eachcol(sort(DataFrame(Tuple.(ensurearray.(permute(ensurearray.(values(x))))))))
    return NamedTuple{Tuple(cols,)}(xperm,)
end

function permute(x::Array)
    xperm = if any(isarray.(x))
        permute(Tuple(x))
    elseif length(unique(length.(x))) .== 1 && all(length.(x) .> 1)
        permute(unique.(eachcol(DataFrame(x))))
    else
        sort(unique(x))
    end
    return xperm
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
struct PrefIntensitySpec{T} <: SpatQSimSpec
    realization::Int
    pref_power::T
    prep_file::String

    function PrefIntensitySpec(realization::Int,
                                pref_power::T,
                                prep_file::String = "prep.h5") where T
        new{T}(realization, pref_power, prep_file)
    end
end

sim_value(spec::PrefIntensitySpec) = spec.pref_power

function comm_targeting(spec::PrefIntensitySpec, prep::SpatQSimPrep)
    comm_q = comm_catchability(spec, prep)
    DynamicPreferentialTargeting(init_pop(prep).P,
                                 p -> (comm_q.catchability .* p) .^ sim_value(spec))
end


text:
using Oceananigans.Solvers: solve!, MatrixIterativeSolver, spai_preconditioner
using Oceananigans.Fields: interior_copy
using Oceananigans.Operators: volume, Δyᶠᶜᵃ, Δyᶜᶠᵃ, Δxᶠᶜᵃ, Δxᶜᶠᵃ, Δyᵃᶜᵃ, Δxᶜᵃᵃ, Δzᵃᵃᶠ, Δzᵃᵃᶜ, ∇²ᶜᶜᶜ
using Oceananigans.Architectures: arch_array
using KernelAbstractions: @kernel, @index
using Statistics, LinearAlgebra, SparseArrays

function calc_∇²!(∇²ϕ, ϕ, arch, grid)
    fill_halo_regions!(ϕ, arch)
    event = launch!(arch, grid, :xyz, ∇²!, ∇²ϕ, grid, ϕ)
    wait(event)
    fill_halo_regions!(∇²ϕ, arch)
    return nothing
end

function identity_operator!(b, x)
    parent(b) .= parent(x)
    return nothing
end

function run_identity_operator_test(arch, grid)

    N = size(grid)
    M = prod(N)

    b = arch_array(arch, zeros(M))

    A = zeros(N...)
    C =  ones(N...)
    D = arch_array(arch, zeros(N...))

    solver = MatrixIterativeSolver((A, A, A, C, D), grid = grid)

    fill!(b, rand())

    initial_guess = solution = CenterField(arch, grid)
    set!(initial_guess, (x, y, z) -> rand())

    solve!(initial_guess, solver, b, 1.0)

    b = reshape(b, size(grid)...)
    @test norm(interior(solution) .- b) .< solver.tolerance
end

@kernel function _multiply_by_volume!(r, grid)
    i, j, k = @index(Global, NTuple)
    r[i, j, k] *= volume(i, j, k, grid, Center(), Center(), Center())
end

function compute_poisson_weights(grid)
    N = size(grid)
    Ax = zeros(N...)
    Ay = zeros(N...)
    Az = zeros(N...)
    C  = zeros(N...)
    D  = arch_array(grid.architecture, zeros(N...))
    for i =1:grid.Nx, j = 1:grid.Ny, k = 1:grid.Nz
        Ax[i, j, k] = Δzᵃᵃᶜ(i, j, k, grid) * Δyᶠᶜᵃ(i, j, k, grid) / Δxᶠᶜᵃ(i, j, k, grid)
        Ay[i, j, k] = Δzᵃᵃᶜ(i, j, k, grid) * Δxᶜᶠᵃ(i, j, k, grid) / Δyᶜᶠᵃ(i, j, k, grid)
        Az[i, j, k] = Δxᶜᵃᵃ(i, j, k, grid) * Δyᵃᶜᵃ(i, j, k, grid) / Δzᵃᵃᶠ(i, j, k, grid)
    end
    return (Ax, Ay, Az, C, D)
end

function poisson_rhs!(r, grid)
    event = launch!(grid.architecture, grid, :xyz, _multiply_by_volume!, r, grid)
    wait(event)
    return
end

function run_poisson_equation_test(arch, grid)
    # Solve ∇²ϕ = r
    ϕ_truth = Field(Center, Center, Center, arch, grid)

    # Initialize zero-mean "truth" solution with random numbers
    set!(ϕ_truth, (x, y, z) -> rand())
    parent(ϕ_truth) .-= mean(ϕ_truth)
    fill_halo_regions!(ϕ_truth, arch)

    # Calculate Laplacian of "truth"
    ∇²ϕ = Field(Center, Center, Center, arch, grid)
    calc_∇²!(∇²ϕ, ϕ_truth, arch, grid)

    rhs = deepcopy(∇²ϕ)
    poisson_rhs!(rhs, grid)
    rhs = interior_copy(rhs)[:]
    weights = compute_poisson_weights(grid)
    solver  = MatrixIterativeSolver(weights, grid = grid)

    # Solve Poisson equation
    ϕ_solution = CenterField(arch, grid)

    solve!(ϕ_solution, solver, rhs, 1.0)

    # Diagnose Laplacian of solution
    ∇²ϕ_solution = CenterField(arch, grid)
    calc_∇²!(∇²ϕ_solution, ϕ_solution, arch, grid)

    parent(ϕ_solution) .-= mean(ϕ_solution)

    CUDA.@allowscalar begin
        @test all(interior(∇²ϕ_solution) .≈ interior(∇²ϕ))
        @test all(interior(ϕ_solution)   .≈ interior(ϕ_truth))
    end

    return nothing
end

@testset "MatrixIterativeSolver" begin
    topologies = [(Periodic, Periodic, Flat), (Bounded, Bounded, Flat), (Periodic, Bounded, Flat), (Bounded, Periodic, Flat)]
    for arch in archs, topo in topologies
        @info "Testing 2D MatrixIterativeSolver [$(typeof(arch)) $topo]..."

        grid = RectilinearGrid(arch, size=(4, 8), extent=(1, 3), topology = topo)
        run_identity_operator_test(arch, grid)
        run_poisson_equation_test(arch, grid)
    end
    topologies = [(Periodic, Periodic, Periodic), (Bounded, Bounded, Periodic), (Periodic, Bounded, Periodic), (Bounded, Periodic, Bounded)]
    for arch in archs, topo in topologies
        @info "Testing 3D MatrixIterativeSolver [$(typeof(arch)) $topo]..."

        grid = RectilinearGrid(arch, size=(4, 8, 6), extent=(1, 3, 4), topology = topo)
        run_identity_operator_test(arch, grid)
        run_poisson_equation_test(arch, grid)
    end
    stretch_coord = [0, 1.5, 3, 7, 8.5, 10]
    for arch in archs
        grids = [RectilinearGrid(arch, size=(5, 5, 5), x = stretch_coord, y = (0, 10), z = (0, 10), topology = (Periodic, Periodic, Periodic)),
                 RectilinearGrid(arch, size=(5, 5, 5), x = (0, 10), y = stretch_coord, z = (0, 10), topology = (Periodic, Periodic, Periodic)),
                 RectilinearGrid(arch, size=(5, 5, 5), x = (0, 10), y = (0, 10), z = stretch_coord, topology = (Periodic, Periodic, Periodic))]

        for grid in grids
            @info "Testing stretched grid MatrixIterativeSolver [$(typeof(arch))]..."
            run_poisson_equation_test(arch, grid)
        end
    end

    @info "Testing Sparse Approximate Inverse Preconditioner"

    A   = sprand(100, 100, 0.1)
    A   = A + A' + 1I
    A⁻¹ = sparse(inv(Array(A)))
    M   = spai_preconditioner(A, ε = 0.0, nzrel = size(A, 1))

    @test all(M .≈ A⁻¹)

end


text:
# === PRIVATE METHODS BELOW HERE ====================================================================================== #
# === PRIVATE METHODS ABOVE HERE ====================================================================================== #

# === PUBLIC METHODS BELOW HERE ======================================================================================= #
function compute_price_return(price_vector::Array{Float64,1})::Array{Float64,1}

    try

        # initialize -
        number_of_items = length(price_vector)
        price_return_array = Array{Float64,1}(undef, (number_of_items - 1))

        # main loop -
        for time_index = 1:(number_of_items - 1)

            # grab the prices -
            today_price = price_vector[time_index]
            tomorrow_price = price_vector[time_index + 1]

            # calculate the return -
            r = log(tomorrow_price / today_price)

            # cache -
            price_return_array[time_index] = r
        end

        # return -
        return price_return_array
    catch error
        rethrow(error)
    end
end

function compute_price_return(price_array::Array{Float64,2})::Array{Float64,2}

    try

        # initialize -
        (number_of_steps, number_of_cols) = size(price_array)
        price_return_array = Array{Float64,2}(undef, (number_of_steps - 1), number_of_cols)

        # main loop -
        for col_index = 1:number_of_cols

            # grab a col -
            data_vector = price_array[:,col_index]

            # compute the return for this sequence -
            local_return_array = compute_price_return(data_vector)

            # package -
            for time_index = 1:(number_of_steps - 1)
                price_return_array[time_index,col_index] = local_return_array[time_index]
            end
        end

        # return -
        return price_return_array
    catch error
        rethrow(error)
    end
end

function compute_price_return(simulation_result_array::Array{DataFrame,1})::Array{Float64,2}

    try

        # how many sample tables do we have?
        number_of_samples = length(simulation_result_array)

        # grab a table to get the dimensions -
        tmp_table = simulation_result_array[1]
        (number_of_steps, _) = size(tmp_table)

        # initialize storage -
        price_return_array = Array{Float64,2}(undef, (number_of_steps - 1), number_of_samples)

        # process each sample -
        for sample_index = 1:number_of_samples

            # grab the price data -
            P = simulation_result_array[sample_index][!,:price]

            # compute the return -
            r_array = compute_price_return(P) # this has number_of_steps - 1 elements

            # copy the data into the storage -
            for step_index = 1:(number_of_steps - 1)
                price_return_array[step_index, sample_index] = r_array[step_index]
            end
        end

        # return -
        return price_return_array
    catch error
        rethrow(error)
    end
end

function compute_autocorrelation_array(price_return_vector::Array{Float64,1},
    lag_array::Array{Int64,1})::Array{Float64,1}

    try

        # return -
        return autocor(price_return_vector, lag_array)
    catch error
        rethrow(error)
    end
end

function compute_autocorrelation_array(price_return_array::Array{Float64,2},
    lag_array::Array{Int64,1})::Array{Float64,2}

    try

        # initialize -
        (number_of_steps, number_of_samples) = size(price_return_array)
        number_of_lags = length(lag_array)
        autocor_array = Array{Float64,2}(undef, number_of_lags, number_of_samples)

        # process each sample -
        for sample_index = 1:number_of_samples

            # data -
            data_col = price_return_array[:, sample_index]
            local_autocor_array = compute_autocorrelation_array(data_col, lag_array)

            # package -
            for lag_index = 1:number_of_lags
                autocor_array[lag_index, sample_index] = local_autocor_array[lag_index]
            end
        end

        # return -
        return autocor_array
    catch error
        rethrow(error)
    end
end
# === PUBLIC METHODS ABOVE HERE ======================================================================================= #


text:
@testset "NoiseWrapper" begin

using StochasticDiffEq,  DiffEqBase, DiffEqNoiseProcess, Test

f1 = (u,p,t) -> 1.01u
g1 = (u,p,t) -> 1.01u
dt = 1//2^(4)
prob1 = SDEProblem(f1,g1,1.0,(0.0,1.0))
sol1 = solve(prob1,EM(),dt=dt,save_noise = true)

W2 = NoiseWrapper(sol1.W)

prob1 = SDEProblem(f1,g1,1.0,(0.0,1.0),noise=W2)

sol2 = solve(prob1,EM(),dt=dt)

@test sol1.u ≈ sol2.u

W3 = NoiseWrapper(sol1.W)
prob2 = SDEProblem(f1,g1,1.0,(0.0,1.0),noise=W3)

dt = 1//2^(5)
sol3 = solve(prob2,EM(),dt=dt)

f1 = (du,u,p,t) -> du.=1.01u
g1 = (du,u,p,t) -> du.= 1.01u
dt = 1//2^(4)
prob1 = SDEProblem(f1,g1,ones(4),(0.0,1.0))
sol1 = solve(prob1,EM(),dt=dt,save_noise = true)

W2 = NoiseWrapper(sol1.W)

prob1 = SDEProblem(f1,g1,ones(4),(0.0,1.0),noise=W2)

sol2 = solve(prob1,EM(),dt=dt)

for i in 1:length(sol1)
  @test sol1.u[i] ≈ sol2.u[i]
end

W3 = NoiseWrapper(sol1.W)
prob2 = SDEProblem(f1,g1,ones(4),(0.0,1.0),noise=W3)

dt = 1//2^(5)
sol3 = solve(prob2,EM(),dt=dt)

prob = SDEProblem(f1,g1,ones(2),(0.0,1.0))
sol4 = solve(prob,SRI(),abstol=1e-6,save_noise = true)

W2 = NoiseWrapper(sol4.W)
prob2 = SDEProblem(f1,g1,ones(2),(0.0,1.0),noise=W2)
sol5 = solve(prob2,SRIW1(),abstol=1e-8)

@test sol4.t != sol5.t
@test ≈(sol5[end],sol4[end],atol=1e-1)

end


text:
"""
Provides virtual interface for the black-box system under test (SUT).
"""
module BlackBox

include("GrayBox.jl")

export initialize!,
       evaluate!,
       distance,
       isevent,
       isterminal

"""
    initialize!(sim::GrayBox.Simulation)

Reset state to its initial state.
"""
function initialize!(sim::GrayBox.Simulation)::Nothing end

"""
    evaluate!(sim::GrayBox.Simulation)::Tuple{logprob::Real, miss_distance::Real, isevent::Bool}
    evaluate!(sim::GrayBox.Simulation, sample::GrayBox.EnvironmentSample)::Tuple{logprob::Real, miss_distance::Real, isevent::Bool}

Evaluate the SUT given some input seed and current state, returns `logprob`, `miss_distance`, and `isevent` indication.
If the `sample` version is implemented, then `ASTSampleAction` will be used instead of `ASTSeedAction`.
"""
function evaluate!(sim::GrayBox.Simulation)::Tuple{Real, Real, Bool} end
function evaluate!(sim::GrayBox.Simulation, sample::GrayBox.EnvironmentSample)::Tuple{Real, Real, Bool} end

"""
    distance(sim::GrayBox.Simulation)::Real

Return how close to an event a terminal state was (i.e. some measure of "miss distance" to the event of interest).
"""
function distance(sim::GrayBox.Simulation)::Real end

"""
    isevent(sim::GrayBox.Simulation)::Bool

Return a boolean indicating if the SUT reached a failure event of interest.
"""
function isevent(sim::GrayBox.Simulation)::Bool end

"""
    isterminal(sim::GrayBox.Simulation)::Bool

Return an indication that the simulation is in a terminal state.
"""
function isterminal(sim::GrayBox.Simulation)::Bool end

end # module BlackBox