text:
module Artifacts

import Base: get, SHA1
import ..depots1, ..depots, ..set_readonly
import ..GitTools
using ..BinaryPlatforms
import ..TOML
import ..Types: parse_toml, write_env_usage
using ..PlatformEngines
using SHA

export create_artifact, artifact_exists, artifact_path, remove_artifact, verify_artifact,
       artifact_meta, artifact_hash, bind_artifact!, unbind_artifact!, download_artifact,
       find_artifacts_toml, ensure_artifact_installed, @artifact_str, archive_artifact

# keep in sync with Base.project_names and Base.manifest_names
const artifact_names = ("JuliaArtifacts.toml", "Artifacts.toml")

const ARTIFACTS_DIR_OVERRIDE = Ref{Union{String,Nothing}}(nothing)
"""
    with_artifacts_directory(f::Function, artifacts_dir::String)

Helper function to allow temporarily changing the artifact installation and search
directory.  When this is set, no other directory will be searched for artifacts, and new
artifacts will be installed within this directory.  Similarly, removing an artifact will
only effect the given artifact directory.  To layer artifact installation locations, use
the typical Julia depot path mechanism.
"""
function with_artifacts_directory(f::Function, artifacts_dir::String)
    try
        ARTIFACTS_DIR_OVERRIDE[] = artifacts_dir
        f()
    finally
        ARTIFACTS_DIR_OVERRIDE[] = nothing
    end
end

"""
    artifacts_dirs(args...)

Return a list of paths joined into all possible artifacts directories, as dictated by the
current set of depot paths and the current artifact directory override via the method
`with_artifacts_dir()`.
"""
function artifacts_dirs(args...)
    if ARTIFACTS_DIR_OVERRIDE[] === nothing
        return [abspath(depot, "artifacts", args...) for depot in depots()]
    else
        # If we've been given an override, use _only_ that directory.
        return [abspath(ARTIFACTS_DIR_OVERRIDE[], args...)]
    end
end

"""
    ARTIFACT_OVERRIDES

Artifact locations can be overridden by writing `Override.toml` files within the artifact
directories of Pkg depots.  For example, in the default depot `~/.julia`, one may create
a `~/.julia/artifacts/Override.toml` file with the following contents:

    78f35e74ff113f02274ce60dab6e92b4546ef806 = "/path/to/replacement"
    c76f8cda85f83a06d17de6c57aabf9e294eb2537 = "fb886e813a4aed4147d5979fcdf27457d20aa35d"

    [d57dbccd-ca19-4d82-b9b8-9d660942965b]
    c_simple = "/path/to/c_simple_dir"
    libfoo = "fb886e813a4aed4147d5979fcdf27457d20aa35d""

This file defines four overrides; two which override specific artifacts identified
through their content hashes, two which override artifacts based on their bound names
within a particular package's UUID.  In both cases, there are two different targets of
the override: overriding to an on-disk location through an absolutet path, and
overriding to another artifact by its content-hash.
"""
const ARTIFACT_OVERRIDES = Ref{Union{Dict,Nothing}}(nothing)
function load_overrides(;force::Bool = false)
    if ARTIFACT_OVERRIDES[] !== nothing && !force
        return ARTIFACT_OVERRIDES[]
    end

    # We organize our artifact location overrides into two camps:
    #  - overrides per UUID with artifact names mapped to a new location
    #  - overrides per hash, mapped to a new location.
    #
    # Overrides per UUID/bound name are intercepted upon Artifacts.toml load, and new
    # entries within the "hash" overrides are generated on-the-fly.  Thus, all redirects
    # mechanisticly happen through the "hash" overrides.
    overrides = Dict(
        # Overrides by UUID
        :UUID => Dict{Base.UUID,Dict{String,Union{String,SHA1}}}(),

        # Overrides by hash
        :hash => Dict{SHA1,Union{String,SHA1}}(),
    )

    for override_file in reverse(artifacts_dirs("Overrides.toml"))
        !isfile(override_file) && continue

        # Load the toml file
        depot_override_dict = parse_toml(override_file)

        function parse_mapping(mapping::String, name::String)
            if !isabspath(mapping) && !isempty(mapping)
                try
                    mapping = Base.SHA1(mapping)
                catch e
                    @error("Invalid override in '$(override_file)': entry '$(name)' must map to an absolute path or SHA1 hash!")
                    rethrow()
                end
            end
            return mapping
        end
        function parse_mapping(mapping::Dict, name::String)
            return Dict(k => parse_mapping(v, name) for (k, v) in mapping)
        end

        for (k, mapping) in depot_override_dict
            # First, parse the mapping. Is it an absolute path, a valid SHA1-hash, or neither?
            try
                mapping = parse_mapping(mapping, k)
            catch
                @error("Invalid override in '$(override_file)': failed to parse entry `$(k)`")
                continue
            end

            # Next, determine if this is a hash override or a UUID/name override
            if isa(mapping, String) || isa(mapping, SHA1)
                # if this mapping is a direct mapping (e.g. a String), store it as a hash override
                hash = try
                    Base.SHA1(hex2bytes(k))
                catch
                    @error("Invalid override in '$(override_file)': Invalid SHA1 hash '$(k)'")
                    continue
                end

                # If this mapping is the empty string, un-override it
                if mapping == ""
                    delete!(overrides[:hash], hash)
                else
                    overrides[:hash][hash] = mapping
                end
            elseif isa(mapping, Dict)
                # Convert `k` into a uuid
                uuid = try
                    Base.UUID(k)
                catch
                    @error("Invalid override in '$(override_file)': Invalid UUID '$(k)'")
                    continue
                end

                # If this mapping is itself a dict, store it as a set of UUID/artifact name overrides
                if !haskey(overrides[:UUID], uuid)
                    overrides[:UUID][uuid] = Dict{String,Union{String,SHA1}}()
                end

                # For each name in the mapping, update appropriately
                for name in keys(mapping)
                    # If the mapping for this name is the empty string, un-override it
                    if mapping[name] == ""
                        delete!(overrides[:UUID][uuid], name)
                    else
                        # Otherwise, store it!
                        overrides[:UUID][uuid][name] = mapping[name]
                    end
                end
            end
        end
    end

    ARTIFACT_OVERRIDES[] = overrides
end

# Helpers to map an override to an actual path
map_override_path(x::String) = x
map_override_path(x::SHA1) = artifact_path(x)
map_override_path(x::Nothing) = nothing

"""
    query_override(hash::SHA1; overrides::Dict = load_overrides())

Query the loaded `<DEPOT>/artifacts/Overrides.toml` settings for artifacts that should be
redirected to a particular path or another content-hash.
"""
function query_override(hash::SHA1; overrides::Dict = load_overrides())
    return map_override_path(get(overrides[:hash], hash, nothing))
end
function query_override(pkg::Base.UUID, artifact_name::String; overrides::Dict = load_overrides())
    if haskey(overrides[:UUID], pkg)
        return map_override_path(get(overrides[:UUID][pkg], artifact_name, nothing))
    end
    return nothing
end

"""
    create_artifact(f::Function)

Creates a new artifact by running `f(artifact_path)`, hashing the result, and moving it
to the artifact store (`~/.julia/artifacts` on a typical installation).  Returns the
identifying tree hash of this artifact.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function create_artifact(f::Function)
    # Ensure the `artifacts` directory exists in our default depot
    artifacts_dir = first(artifacts_dirs())
    mkpath(artifacts_dir)

    # Temporary directory where we'll do our creation business
    temp_dir = mktempdir(artifacts_dir)

    try
        # allow the user to do their work inside the temporary directory
        f(temp_dir)

        # Calculate the tree hash for this temporary directory
        artifact_hash = SHA1(GitTools.tree_hash(temp_dir))

        # If we created a dupe, just let the temp directory get destroyed. It's got the
        # same contents as whatever already exists after all, so it doesn't matter.  Only
        # move its contents if it actually contains new contents.  Note that we explicitly
        # set `honor_overrides=false` here, as we wouldn't want to drop things into the
        # system directory by accidentally creating something with the same content-hash
        # as something that was foolishly overridden.  This should be virtually impossible
        # unless the user has been very unwise, but let's be cautious.
        new_path = artifact_path(artifact_hash; honor_overrides=false)
        if !isdir(new_path)
            # Move this generated directory to its final destination, set it to read-only
            mv(temp_dir, new_path)
            set_readonly(new_path)
        end

        # Give the people what they want
        return artifact_hash
    finally
        # Always attempt to cleanup
        rm(temp_dir; recursive=true, force=true)
    end
end

"""
    artifact_paths(hash::SHA1; honor_overrides::Bool=true)

Return all possible paths for an artifact given the current list of depots as returned
by `Pkg.depots()`.  All, some or none of these paths may exist on disk.
"""
function artifact_paths(hash::SHA1; honor_overrides::Bool=true)
    # First, check to see if we've got an override:
    if honor_overrides
        override = query_override(hash)
        if override !== nothing
            return [override]
        end
    end

    return artifacts_dirs(bytes2hex(hash.bytes))
end

"""
    artifact_path(hash::SHA1; honor_overrides::Bool=true)

Given an artifact (identified by SHA1 git tree hash), return its installation path.  If
the artifact does not exist, returns the location it would be installed to.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function artifact_path(hash::SHA1; honor_overrides::Bool=true)
    # Get all possible paths (rooted in all depots)
    possible_paths = artifact_paths(hash; honor_overrides=honor_overrides)

    # Find the first path that exists and return it
    for p in possible_paths
        if isdir(p)
            return p
        end
    end

    # If none exist, then just return the one that would exist within `depots1()`.
    return first(possible_paths)
end

"""
    artifact_exists(hash::SHA1; honor_overrides::Bool=true)

Returns whether or not the given artifact (identified by its sha1 git tree hash) exists
on-disk.  Note that it is possible that the given artifact exists in multiple locations
(e.g. within multiple depots).

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function artifact_exists(hash::SHA1; honor_overrides::Bool=true)
    return any(isdir.(artifact_paths(hash; honor_overrides=honor_overrides)))
end

"""
    remove_artifact(hash::SHA1; honor_overrides::Bool=false)

Removes the given artifact (identified by its SHA1 git tree hash) from disk.  Note that
if an artifact is installed in multiple depots, it will be removed from all of them.  If
an overridden artifact is requested for removal, it will be silently ignored; this method
will never attempt to remove an overridden artifact.

In general, we recommend that you use `Pkg.gc()` to manage artifact installations and do
not use `remove_artifact()` directly, as it can be difficult to know if an artifact is
being used by another package.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function remove_artifact(hash::SHA1)
    if query_override(hash) !== nothing
        # We never remove overridden artifacts.
        return
    end

    # Get all possible paths (rooted in all depots)
    possible_paths = artifacts_dirs(bytes2hex(hash.bytes))
    for path in possible_paths
        if isdir(path)
            rm(path; recursive=true, force=true)
        end
    end
end

"""
    verify_artifact(hash::SHA1; honor_overrides::Bool=false)

Verifies that the given artifact (identified by its SHA1 git tree hash) is installed on-
disk, and retains its integrity.  If the given artifact is overridden, skips the
verification unless `honor_overrides` is set to `true`.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function verify_artifact(hash::SHA1; honor_overrides::Bool=false)
    # Silently skip overridden artifacts unless we really ask for it
    if !honor_overrides
        if query_override(hash) !== nothing
            return true
        end
    end

    # If it doesn't even exist, then skip out
    if !artifact_exists(hash)
        return false
    end

    # Otherwise actually run the verification
    return hash.bytes == GitTools.tree_hash(artifact_path(hash))
end

"""
    archive_artifact(hash::SHA1, tarball_path::String; honor_overrides::Bool=false)

Archive an artifact into a tarball stored at `tarball_path`, returns the SHA256 of the
resultant tarball as a hexidecimal string. Throws an error if the artifact does not
exist.  If the artifact is overridden, throws an error unless `honor_overrides` is set.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function archive_artifact(hash::SHA1, tarball_path::String; honor_overrides::Bool=false)
    if !honor_overrides
        if query_override(hash) !== nothing
            error("Will not archive an overridden artifact unless `honor_overrides` is set!")
        end
    end

    if !artifact_exists(hash)
        error("Unable to archive artifact $(bytes2hex(hash.bytes)): does not exist!")
    end

    probe_platform_engines!()

    # Package it up
    package(artifact_path(hash), tarball_path)

    # Calculate its sha256 and return that
    return open(tarball_path, "r") do io
        return bytes2hex(sha256(io))
    end
end

"""
    unpack_platform(entry::Dict, name::String, artifacts_toml::String)

Given an `entry` for the artifact named `name`, located within the file `artifacts_toml`,
returns the `Platform` object that this entry specifies.  Returns `nothing` on error.
"""
function unpack_platform(entry::Dict, name::String, artifacts_toml::String)
    if !haskey(entry, "os")
        @error("Invalid artifacts file at '$(artifacts_toml)': platform-specific artifact entry '$name' missing 'os' key")
        return nothing
    end

    if !haskey(entry, "arch")
        @error("Invalid artifacts file at '$(artifacts_toml)': platform-specific artifact entrty '$name' missing 'arch' key")
        return nothing
    end

    # Helpers to pull out `Symbol`s and `VersionNumber`s while preserving `nothing`.
    nosym(x::Nothing) = nothing
    nosym(x) = Symbol(lowercase(x))
    nover(x::Nothing) = nothing
    nover(x) = VersionNumber(x)

    # First, extract OS; we need to build a mapping here
    os_map = Dict(
        "windows" => Windows,
        "macos" => MacOS,
        "freebsd" => FreeBSD,
        "linux" => Linux,
    )
    P = get(os_map, lowercase(entry["os"]), UnknownPlatform)

    # Next, architecture, libc, libgfortran version and cxxabi (if given)
    arch = nosym(get(entry, "arch", nothing))
    libc = nosym(get(entry, "libc", nothing))
    libgfortran_version = nover(get(entry, "libgfortran_version", nothing))
    libstdcxx_version = nover(get(entry, "libstdcxx_version", nothing))
    cxxstring_abi = nosym(get(entry, "cxxstring_abi", nothing))

    # Construct the actual Platform object
    return P(arch;
        libc=libc,
        compiler_abi=CompilerABI(
            libgfortran_version=libgfortran_version,
            libstdcxx_version=libstdcxx_version,
            cxxstring_abi=cxxstring_abi
        ),
    )
end

function pack_platform!(meta::Dict, p::Platform)
    os_map = Dict(
        Windows => "windows",
        MacOS => "macos",
        FreeBSD => "freebsd",
        Linux => "linux",
    )
    meta["os"] = os_map[typeof(p)]
    meta["arch"] = string(arch(p))
    if libc(p) !== nothing
        meta["libc"] = string(libc(p))
    end
    if libgfortran_version(p) !== nothing
        meta["libgfortran_version"] = string(libgfortran_version(p))
    end
    if libstdcxx_version(p) !== nothing
        meta["libstdcxx_version"] = string(libstdcxx_version(p))
    end
    if cxxstring_abi(p) !== nothing
        meta["cxxstring_abi"] = string(cxxstring_abi(p))
    end
end

"""
    load_artifacts_toml(artifacts_toml::String;
                        pkg_uuid::Union{UUID,Nothing}=nothing)

Loads an `(Julia)Artifacts.toml` file from disk.  If `pkg_uuid` is set to the `UUID` of the
owning package, UUID/name overrides stored in a depot `Overrides.toml` will be resolved.
"""
function load_artifacts_toml(artifacts_toml::String;
                             pkg_uuid::Union{Base.UUID,Nothing} = nothing)
    artifact_dict = parse_toml(artifacts_toml)

    # Process overrides for this `pkg_uuid`
    process_overrides(artifact_dict, pkg_uuid)
    return artifact_dict
end

"""
    process_overrides(artifact_dict::Dict, pkg_uuid::Base.UUID)

When loading an `Artifacts.toml` file, we must check `Override.toml` files to see if any
of the artifacts within it have been overridden by UUID.  If they have, we honor the
overrides by inspecting the hashes of the targeted artifacts, then overriding them to
point to the given override, punting the actual redirection off to the hash-based
override system.  This does not modify the `artifact_dict` object, it merely dynamically
adds more hash-based overrides as `Artifacts.toml` files that are overridden are loaded.
"""
function process_overrides(artifact_dict::Dict, pkg_uuid::Base.UUID)
    # Insert just-in-time hash overrides by looking up the names of anything we need to
    # override for this UUID, and inserting new overrides for those hashes.
    overrides = load_overrides()
    if haskey(overrides[:UUID], pkg_uuid)
        pkg_overrides = overrides[:UUID][pkg_uuid]

        for name in keys(artifact_dict)
            # Skip names that we're not overriding
            if !haskey(pkg_overrides, name)
                continue
            end

            # If we've got a platform-specific friend, override all hashes:
            if isa(artifact_dict[name], Array)
                for entry in artifact_dict[name]
                    hash = SHA1(entry["git-tree-sha1"])
                    overrides[:hash][hash] = overrides[:UUID][pkg_uuid][name]
                end
            elseif isa(artifact_dict[name], Dict)
                hash = SHA1(artifact_dict[name]["git-tree-sha1"])
                overrides[:hash][hash] = overrides[:UUID][pkg_uuid][name]
            end
        end
    end
    return artifact_dict
end

# If someone tries to call process_overrides() with `nothing`, do exactly that
process_overrides(artifact_dict::Dict, pkg_uuid::Nothing) = nothing

"""
    artifact_meta(name::String, artifacts_toml::String;
                  platform::Platform = platform_key_abi(),
                  pkg_uuid::Union{Base.UUID,Nothing}=nothing)

Get metadata about a given artifact (identified by name) stored within the given
`(Julia)Artifacts.toml` file.  If the artifact is platform-specific, use `platform` to choose the
most appropriate mapping.  If none is found, return `nothing`.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function artifact_meta(name::String, artifacts_toml::String;
                       platform::Platform = platform_key_abi(),
                       pkg_uuid::Union{Base.UUID,Nothing}=nothing)
    if !isfile(artifacts_toml)
        return nothing
    end

    # Parse the toml of the artifacts_toml file
    artifact_dict = load_artifacts_toml(artifacts_toml; pkg_uuid=pkg_uuid)
    return artifact_meta(name, artifact_dict, artifacts_toml; platform=platform)
end

function artifact_meta(name::String, artifact_dict::Dict, artifacts_toml::String;
                       platform::Platform = platform_key_abi())
    if !haskey(artifact_dict, name)
        return nothing
    end
    meta = artifact_dict[name]

    # If it's an array, find the entry that best matches our current platform
    if isa(meta, Array)
        dl_dict = Dict(unpack_platform(x, name, artifacts_toml) => x for x in meta)
        meta = select_platform(dl_dict, platform)
    # If it's NOT a dict, complain
    elseif !isa(meta, Dict)
        @error("Invalid artifacts file at $(artifacts_toml): artifact '$name' malformed, must be array or dict!")
        return nothing
    end

    # This is such a no-no, we are going to call it out right here, right now.
    if meta !== nothing && !haskey(meta, "git-tree-sha1")
        @error("Invalid artifacts file at $(artifacts_toml): artifact '$name' contains no `git-tree-sha1`!")
        return nothing
    end

    # Return the full meta-dict.
    return meta
end

"""
    artifact_hash(name::String, artifacts_toml::String; platform::Platform = platform_key_abi())

Thin wrapper around `artifact_meta()` to return the hash of the specified, platform-
collapsed artifact.  Returns `nothing` if no mapping can be found.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function artifact_hash(name::String, artifacts_toml::String;
                       platform::Platform = platform_key_abi(),
                       pkg_uuid::Union{Base.UUID,Nothing}=nothing)
    meta = artifact_meta(name, artifacts_toml; platform=platform)
    if meta === nothing
        return nothing
    end

    return SHA1(meta["git-tree-sha1"])
end

"""
    bind_artifact!(artifacts_toml::String, name::String, hash::SHA1;
                   platform::Union{Platform,Nothing} = nothing,
                   download_info::Union{Vector{Tuple},Nothing} = nothing,
                   lazy::Bool = false,
                   force::Bool = false)

Writes a mapping of `name` -> `hash` within the given `(Julia)Artifacts.toml` file. If
`platform` is not `nothing`, this artifact is marked as platform-specific, and will be
a multi-mapping.  It is valid to bind multiple artifacts with the same name, but
different `platform`s and `hash`'es within the same `artifacts_toml`.  If `force` is set
to `true`, this will overwrite a pre-existant mapping, otherwise an error is raised.

`download_info` is an optional tuple that contains a vector of URLs and a hash.  These
URLs will be listed as possible locations where this artifact can be obtained.  If `lazy`
is set to `true`, even if download information is available, this artifact will not be
downloaded until it is accessed via the `artifact"name"` syntax, or
`ensure_artifact_installed()` is called upon it.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function bind_artifact!(artifacts_toml::String, name::String, hash::SHA1;
                        platform::Union{Platform,Nothing} = nothing,
                        download_info::Union{Vector{<:Tuple},Nothing} = nothing,
                        lazy::Bool = false,
                        force::Bool = false)
    # First, check to see if this artifact is already bound:
    if isfile(artifacts_toml)
        artifact_dict = parse_toml(artifacts_toml)

        if !force && haskey(artifact_dict, name)
            meta = artifact_dict[name]
            if !isa(meta, Array)
                error("Mapping for '$name' within $(artifacts_toml) already exists!")
            elseif any((unpack_platform(x, name, artifacts_toml) for x in meta) .== Ref(platform))
                error("Mapping for '$name'/$(triplet(platform)) within $(artifacts_toml) already exists!")
            end
        end
    else
        artifact_dict = Dict()
    end

    # Otherwise, the new piece of data we're going to write out is this dict:
    meta = Dict{String,Any}(
        "git-tree-sha1" => bytes2hex(hash.bytes),
    )

    # If we're set to be lazy, then lazy we shall be
    if lazy
        meta["lazy"] = true
    end

    # Integrate download info, if it is given.  We represent the download info as a
    # vector of dicts, each with its own `url` and `sha256`, since different tarballs can
    # expand to the same tree hash.
    if download_info !== nothing
        meta["download"] = [
            Dict("url" => dl[1],
                 "sha256" => dl[2],
            ) for dl in download_info
        ]
    end

    if platform === nothing
        artifact_dict[name] = meta
    else
        # Add platform-specific keys to our `meta` dict
        pack_platform!(meta, platform)

        # Insert this entry into the list of artifacts
        if !haskey(artifact_dict, name)
            artifact_dict[name] = [meta]
        else
            # Delete any entries that contain identical platforms
            artifact_dict[name] = filter(
                x -> unpack_platform(x, name, artifacts_toml) != platform,
                artifact_dict[name]
            )
            push!(artifact_dict[name], meta)
        end
    end

    # Spit it out onto disk
    open(artifacts_toml, "w") do io
        TOML.print(io, artifact_dict, sorted=true)
    end

    # Mark that we have used this Artifact.toml
    write_env_usage(artifacts_toml, "artifact_usage.toml")
    return
end

"""
    unbind_artifact!(artifacts_toml::String, name::String; platform = nothing)

Unbind the given `name` from an `(Julia)Artifacts.toml` file.
Silently fails if no such binding exists within the file.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function unbind_artifact!(artifacts_toml::String, name::String;
                         platform::Union{Platform,Nothing} = nothing)
    artifact_dict = parse_toml(artifacts_toml)
    if !haskey(artifact_dict, name)
        return
    end

    if platform === nothing
        delete!(artifact_dict, name)
    else
        artifact_dict[name] = filter(
            x -> unpack_platform(x, name, artifacts_toml) != platform,
            artifact_dict[name]
        )
    end

    open(artifacts_toml, "w") do io
        TOML.print(io, artifact_dict, sorted=true)
    end
    return
end

"""
    download_artifact(tree_hash::SHA1, tarball_url::String, tarball_hash::String;
                      verbose::Bool = false)

Download/install an artifact into the artifact store.  Returns `true` on success.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function download_artifact(tree_hash::SHA1, tarball_url::String, tarball_hash::String;
                           verbose::Bool = false)
    if artifact_exists(tree_hash)
        return true
    end

    # Ensure that we're ready to download things
    probe_platform_engines!()

    if Sys.iswindows()
        # The destination directory we're hoping to fill:
        dest_dir = artifact_path(tree_hash; honor_overrides=false)

        # On Windows, we have some issues around stat() and chmod() that make properly
        # determining the git tree hash problematic; for this reason, we use the "unsafe"
        # artifact unpacking method, which does not properly verify unpacked git tree
        # hash.  This will be fixed in a future Julia release which will properly interrogate
        # the filesystem ACLs for executable permissions, which git tree hashes care about.
        try
            download_verify_unpack(tarball_url, tarball_hash, dest_dir, ignore_existence=true, verbose=verbose)
        catch e
            # Clean that destination directory out if something went wrong
            rm(dest_dir; force=true, recursive=true)

            if isa(e, InterruptException)
                rethrow(e)
            end
            return false
        end
    else
        # We download by using `create_artifact()`.  We do this because the download may
        # be corrupted or even malicious; we don't want to clobber someone else's artifact
        # by trusting the tree hash that has been given to us; we will instead download it
        # to a temporary directory, calculate the true tree hash, then move it to the proper
        # location only after knowing what it is, and if something goes wrong in the process,
        # everything should be cleaned up.  Luckily, that is precisely what our
        # `create_artifact()` wrapper does, so we use that here.
        calc_hash = try
            create_artifact() do dir
                download_verify_unpack(tarball_url, tarball_hash, dir, ignore_existence=true, verbose=verbose)
            end
        catch e
            if isa(e, InterruptException)
                rethrow(e)
            end
            # If something went wrong during download, return false
            return false
        end

        # Did we get what we expected?  If not, freak out.
        if calc_hash.bytes != tree_hash.bytes
            msg  = "Tree Hash Mismatch!\n"
            msg *= "  Expected git-tree-sha1:   $(bytes2hex(tree_hash.bytes))\n"
            msg *= "  Calculated git-tree-sha1: $(bytes2hex(calc_hash.bytes))"
            @error(msg)
            return false
        end
    end

    return true
end

"""
    find_artifacts_toml(path::String)

Given the path to a `.jl` file, (such as the one returned by `__source__.file` in a macro
context), find the `(Julia)Artifacts.toml` that is contained within the containing project (if it
exists), otherwise return `nothing`.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function find_artifacts_toml(path::String)
    if !isdir(path)
        path = dirname(path)
    end

    # Run until we hit the root directory.
    while dirname(path) != path
        for f in artifact_names
            artifacts_toml_path = joinpath(path, f)
            if isfile(artifacts_toml_path)
                return abspath(artifacts_toml_path)
            end
        end

        # Does a `(Julia)Project.toml` file exist here, in the absence of an Artifacts.toml?
        # If so, stop the search as we've probably hit the top-level of this package,
        # and we don't want to escape out into the larger filesystem.
        for f in Base.project_names
            if isfile(joinpath(path, f))
                return nothing
            end
        end

        # Move up a directory
        path = dirname(path)
    end

    # We never found anything, just return `nothing`
    return nothing
end

"""
    ensure_artifact_installed(name::String, artifacts_toml::String;
                              platform::Platform = platform_key_abi(),
                              pkg_uuid::Union{Base.UUID,Nothing}=nothing)

Ensures an artifact is installed, downloading it via the download information stored in
`artifacts_toml` if necessary.  Throws an error if unable to install.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function ensure_artifact_installed(name::String, artifacts_toml::String;
                                   platform::Platform = platform_key_abi(),
                                   pkg_uuid::Union{Base.UUID,Nothing}=nothing,
                                   verbose::Bool = false)
    meta = artifact_meta(name, artifacts_toml; pkg_uuid=pkg_uuid, platform=platform)
    if meta === nothing
        error("Cannot locate artifact '$(name)' in '$(artifacts_toml)'")
    end

    return ensure_artifact_installed(name, meta, artifacts_toml; platform=platform, verbose=verbose)
end

function ensure_artifact_installed(name::String, meta::Dict, artifacts_toml::String;
                                   platform::Platform = platform_key_abi(),
                                   verbose::Bool = false)
    hash = SHA1(meta["git-tree-sha1"])

    if !artifact_exists(hash)
        # If this artifact does not exist on-disk already, ensure it has download
        # information, then download it!
        if !haskey(meta, "download")
            error("Cannot automatically install '$(name)'; no download section in '$(artifacts_toml)'")
        end

        # Attempt to download from all sources
        for entry in meta["download"]
            url = entry["url"]
            tarball_hash = entry["sha256"]
            if download_artifact(hash, url, tarball_hash; verbose=verbose)
                return artifact_path(hash)
            end
        end
        error("Unable to automatically install '$(name)' from '$(artifacts_toml)'")
    else
        return artifact_path(hash)
    end
end

"""
    ensure_all_artifacts_installed(artifacts_toml::String;
                                   platform = platform_key_abi(),
                                   pkg_uuid = nothing,
                                   include_lazy = false,
                                   verbose = false)

Installs all non-lazy artifacts from a given `(Julia)Artifacts.toml` file. `package_uuid` must
be provided to properly support overrides from `Overrides.toml` entries in depots.

If `include_lazy` is set to `true`, then lazy packages will be installed as well.

!!! compat "Julia 1.3"
    This function requires at least Julia 1.3.
"""
function ensure_all_artifacts_installed(artifacts_toml::String;
                                        platform::Platform = platform_key_abi(),
                                        pkg_uuid::Union{Nothing,Base.UUID} = nothing,
                                        include_lazy::Bool = false,
                                        verbose::Bool = false)
    if !isfile(artifacts_toml)
        return
    end
    artifact_dict = load_artifacts_toml(artifacts_toml; pkg_uuid=pkg_uuid)

    for name in keys(artifact_dict)
        # Get the metadata about this name for the requested platform
        meta = artifact_meta(name, artifact_dict, artifacts_toml; platform=platform)

        # If there are no instances of this name for the desired platform, skip it
        meta === nothing && continue

        # If this mapping doesn't have a `download` stanza or is lazy, skip it
        if !haskey(meta, "download") || (get(meta, "lazy", false) && !include_lazy)
            continue
        end

        # Otherwise, let's try and install it!
        ensure_artifact_installed(name, meta, artifacts_toml; platform=platform, verbose=verbose)
    end
end

"""
    extract_all_hashes(artifacts_toml::String;
                       platform = platform_key_abi(),
                       pkg_uuid = nothing,
                       include_lazy = false)

Extract all hashes from a given `(Julia)Artifacts.toml` file. `package_uuid` must
be provided to properly support overrides from `Overrides.toml` entries in depots.

If `include_lazy` is set to `true`, then lazy packages will be installed as well.
"""
function extract_all_hashes(artifacts_toml::String;
                            platform::Platform = platform_key_abi(),
                            pkg_uuid::Union{Nothing,Base.UUID} = nothing,
                            include_lazy::Bool = false)
    hashes = Base.SHA1[]
    if !isfile(artifacts_toml)
        return hashes
    end

    artifact_dict = load_artifacts_toml(artifacts_toml; pkg_uuid=pkg_uuid)

    for name in keys(artifact_dict)
        # Get the metadata about this name for the requested platform
        meta = artifact_meta(name, artifact_dict, artifacts_toml; platform=platform)

        # If there are no instances of this name for the desired platform, skip it
        meta === nothing && continue

        # If it's a lazy one and we aren't including lazy ones, skip
        if get(meta, "lazy", false) && !include_lazy
            continue
        end

        # Otherwise, add it to the list!
        push!(hashes, Base.SHA1(meta["git-tree-sha1"]))
    end

    return hashes
end

function do_artifact_str(name, artifact_dict, artifacts_toml, __module__)
    local pkg_uuid = nothing
    if haskey(Base.module_keys, __module__)
        # Process overrides for this UUID, if we know what it is
        process_overrides(artifact_dict, Base.module_keys[__module__].uuid)
    end

    # Get platform once to avoid extra work
    platform = platform_key_abi()

    # Get the metadata about this name for the requested platform
    meta = artifact_meta(name, artifact_dict, artifacts_toml; platform=platform)

    if meta === nothing
        error("Cannot locate artifact '$(name)' in '$(artifacts_toml)'")
    end

    # This is the resultant value at the end of all things
    return ensure_artifact_installed(name, meta, artifacts_toml; platform=platform)
end

"""
    macro artifact_str(name)

Macro that is used to automatically ensure an artifact is installed, and return its
location on-disk.  Automatically looks the artifact up by name in the project's
`(Julia)Artifacts.toml` file.  Throws an error on inability to install the requested artifact.

!!! compat "Julia 1.3"
    This macro requires at least Julia 1.3.
"""
macro artifact_str(name)
    # Load Artifacts.toml at compile time, so that we don't have to use `__source__.file`
    # at runtime, which gets stale if the `.ji` file is relocated.
    local artifacts_toml = find_artifacts_toml(string(__source__.file))
    if artifacts_toml === nothing
        error(string(
            "Cannot locate '(Julia)Artifacts.toml' file when attempting to use artifact '",
            name,
            "' in '",
            __module__,
            "'",
        ))
    end

    local artifact_dict = load_artifacts_toml(artifacts_toml)
    return quote
        # Invalidate .ji file if Artifacts.toml file changes
        Base.include_dependency($(artifacts_toml))

        # Use invokelatest() to introduce a compiler barrier, preventing many backedges from being added
        # and slowing down not only compile time, but also `.ji` load time.  This is critical here, as
        # artifact"" is used in other modules, so we don't want to be spreading backedges around everywhere.
        Base.invokelatest(do_artifact_str, $name, $(artifact_dict), $(artifacts_toml), $__module__)
    end
end

end # module Artifacts


text:
using DifferentiableStateSpaceModels, Symbolics, LinearAlgebra, Test
using DifferentiableStateSpaceModels.Examples
using DifferentiableStateSpaceModels: order_vector_by_symbols, fill_array_by_symbol_dispatch, all_fields_equal

# # # Use while testing internals
# m = @include_example_module(Examples.rbc_observables)
# # Basic Steady State
# p_f = (ρ=0.2, δ=0.02, σ=0.01, Ω_1=0.01)
# p_d = (α=0.5, β=0.95)
# c = SolverCache(m, Val(2), p_d)
# sol = generate_perturbation(m, p_d, p_f, Val(2); cache = c) # manually passing in f
# ex = DifferentiableStateSpaceModels.exfiltrated

@testset "Second Order Construction" begin
    m = @include_example_module(Examples.rbc_observables)

    # bookkeeping tests
    @test m.n_y == 2
    @test m.n_x == 2
    @test m.n_p == 6
    @test m.n_ϵ == 1
    @test m.n_z == 2
    @test m.η == reshape([0; -1], 2, m.n_ϵ)
    @test m.Q == [1.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0]

    # function tests (steady state)
    # Basic Steady State
    p_f = (ρ=0.2, δ=0.02, σ=0.01, Ω_1=0.01)
    p_d = (α=0.5, β=0.95)
    sol = generate_perturbation(m, p_d, p_f, Val(2))
    @inferred generate_perturbation(m, p_d, p_f, Val(2))
    @test sol.y ≈ [5.936252888048733, 6.884057971014498]
    @test sol.x ≈ [47.39025414828825, 0.0]
    @test sol.retcode == :Success

    # Call all variables differentiated
    sol = generate_perturbation(m, merge(p_d, p_f), nothing, Val(2))
    @test sol.y ≈ [5.936252888048733, 6.884057971014498]
    @test sol.x ≈ [47.39025414828825, 0.0]

    # With a prebuilt cache
    c = SolverCache(m, Val(2), p_d)
    sol = generate_perturbation(m, p_d, p_f, Val(2); cache=c)
    @test sol.y ≈ [5.936252888048733, 6.884057971014498]
    @test sol.x ≈ [47.39025414828825, 0.0]
    @inferred generate_perturbation(m, p_d, p_f, Val(2); cache=c)
end

@testset "Second Order Function Evaluation" begin
    m = @include_example_module(Examples.rbc_observables)
    p_f = (ρ=0.2, δ=0.02, σ=0.01, Ω_1=0.01)
    p_d = (α=0.5, β=0.95)
    p_d_symbols = collect(Symbol.(keys(p_d)))  #The order of derivatives in p_d
    c = SolverCache(m, Val(2), p_d)

    # Create parameter vector in the same ordering the internal algorithms would
    p = order_vector_by_symbols(merge(p_d, p_f), m.mod.p_symbols)

    y = zeros(m.n_y)
    x = zeros(m.n_x)

    m.mod.ȳ!(y, p)
    m.mod.x̄!(x, p)
    @test y ≈ [5.936252888048733, 6.884057971014498]
    @test x ≈ [47.39025414828825, 0.0]

    m.mod.H_yp!(c.H_yp, y, x, p)
    @test c.H_yp ≈ [0.028377570562199098 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]

    m.mod.H_y!(c.H_y, y, x, p)
    @test c.H_y ≈ [-0.0283775705621991 0.0; 1.0 -1.0; 0.0 1.0; 0.0 0.0]

    m.mod.H_xp!(c.H_xp, y, x, p)
    @test c.H_xp ≈ [0.00012263591151906127 -0.011623494029190608
                    1.0 0.0
                    0.0 0.0
                    0.0 1.0]

    m.mod.H_x!(c.H_x, y, x, p)
    @test c.H_x ≈ [0.0 0.0
                   -0.98 0.0
                   -0.07263157894736837 -6.884057971014498
                   0.0 -0.2]

    m.mod.Ψ!(c.Ψ, y, x, p)
    @test c.Ψ[1] ≈
          [-0.009560768753410337 0.0 0.0 0.0 -2.0658808482697935e-5 0.0019580523687917364 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.009560768753410338 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           -2.0658808482697935e-5 0.0 0.0 0.0 -3.881681383327978e-6 0.00012263591151906127 0.0 0.0
           0.0019580523687917364 0.0 0.0 0.0 0.00012263591151906127 -0.011623494029190608 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
    @test c.Ψ[2] ≈ zeros(8, 8)
    @test c.Ψ[3] ≈ [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0007663134567721225 -0.07263157894736837
           0.0 0.0 0.0 0.0 0.0 0.0 -0.07263157894736837 -6.884057971014498]
    @test c.Ψ[4] ≈ zeros(8, 8)

    m.mod.Γ!(c.Γ, p)
    @test c.Γ ≈ [0.01]

    # The derivative ones dispatch by the derivative symbol
    fill_array_by_symbol_dispatch(m.mod.H_x_p!, c.H_x_p, p_d_symbols, y, x, p)
    @test c.H_x_p ≈ [[0.0 0.0
            0.0 0.0
            -0.4255060477077458 -26.561563542978472
            0.0 0.0], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]]
    fill_array_by_symbol_dispatch(m.mod.H_yp_p!, c.H_yp_p, p_d_symbols, y, x, p)
    @test c.H_yp_p ≈ [[0.011471086498795562 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0],
           [0.029871126907577997 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]]

    fill_array_by_symbol_dispatch(m.mod.H_y_p!, c.H_y_p, p_d_symbols, y, x, p)
    @test c.H_y_p ≈
          [[0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]]

    fill_array_by_symbol_dispatch(m.mod.H_xp_p!, c.H_xp_p, p_d_symbols, y, x, p)
    @test c.H_xp_p ≈ [[0.000473180436623283 -0.06809527035753198
            0.0 0.0
            0.0 0.0
            0.0 0.0], [0.00012909043317795924 -0.01223525687283222
                       0.0 0.0
                       0.0 0.0
                       0.0 0.0]]

    fill_array_by_symbol_dispatch(m.mod.Γ_p!, c.Γ_p, p_d_symbols, p)

    @test c.Γ_p ≈ [[0.0], [0.0]]

    fill_array_by_symbol_dispatch(m.mod.H_p!, c.H_p, p_d_symbols, y, x, p)

    @test c.H_p ≈ [[-0.06809527035753199, 0.0, -26.561563542978472, 0.0],
           [-0.1773225633743801, 0.0, 0.0, 0.0]]
    fill_array_by_symbol_dispatch(m.mod.Ω_p!, c.Ω_p, p_d_symbols, p)
    @test c.Ω_p ≈ [[0.0, 0.0], [0.0, 0.0]]

    # Second order checks!
    m.mod.Ψ_yp!(c.Ψ_yp, y, x, p)
    @test c.Ψ_yp[1] ≈
          [[0.0048317190660755365 0.0 0.0 0.0 6.960218465183541e-6 -0.0006596930439853133 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            6.960218465183541e-6 0.0 0.0 0.0 6.53894208439611e-7 -2.0658808482697952e-5 0.0 0.0
            -0.0006596930439853133 0.0 0.0 0.0 -2.0658808482697952e-5 0.0019580523687917372 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], zeros(8, 8), zeros(8, 8), zeros(8, 8)]
    @test c.Ψ_yp[2] ≈ [zeros(8, 8), zeros(8, 8), zeros(8, 8), zeros(8, 8)]

    m.mod.Ψ_y!(c.Ψ_y, y, x, p)
    @test c.Ψ_y[1] ≈ [[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 -0.0048317190660755365 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], zeros(8, 8), zeros(8, 8), zeros(8, 8)]

    @test c.Ψ_y[2] ≈ [zeros(8, 8), zeros(8, 8), zeros(8, 8), zeros(8, 8)]

    m.mod.Ψ_xp!(c.Ψ_xp, y, x, p)
    @test c.Ψ_xp[1] ≈
          [[6.9602184651835404e-6 0.0 0.0 0.0 6.53894208439611e-7 -2.0658808482697952e-5 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            6.53894208439611e-7 0.0 0.0 0.0 2.0477213369556226e-7 -3.881681383327981e-6 0.0 0.0
            -2.0658808482697952e-5 0.0 0.0 0.0 -3.881681383327981e-6 0.00012263591151906138 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], zeros(8, 8), zeros(8, 8), zeros(8, 8)]
    @test c.Ψ_xp[2] ≈
          [[-0.0006596930439853132 0.0 0.0 0.0 -2.0658808482697952e-5 0.0019580523687917372 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            -2.0658808482697952e-5 0.0 0.0 0.0 -3.881681383327981e-6 0.00012263591151906138 0.0 0.0
            0.0019580523687917372 0.0 0.0 0.0 0.00012263591151906138 -0.011623494029190612 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], zeros(8, 8), zeros(8, 8), zeros(8, 8)]
    m.mod.Ψ_x!(c.Ψ_x, y, x, p)
    @test c.Ψ_x[1] ≈ [zeros(8, 8), zeros(8, 8),
           [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 -2.4255412970806024e-5 0.000766313456772123
            0.0 0.0 0.0 0.0 0.0 0.0 0.000766313456772123 -0.07263157894736838], zeros(8, 8)]
    @test c.Ψ_x[2] ≈ [zeros(8, 8), zeros(8, 8),
           [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.000766313456772123 -0.07263157894736838
            0.0 0.0 0.0 0.0 0.0 0.0 -0.07263157894736838 -6.884057971014497], zeros(8, 8)]

    fill_array_by_symbol_dispatch(m.mod.Ψ_p!, c.Ψ_p, p_d_symbols, y, x, p)
    @test c.Ψ_p[1] ≈
          [[-0.0038647566790457792 0.0 0.0 0.0 -7.971028956261657e-5 0.011471086498795566 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            -7.971028956261657e-5 0.0 0.0 0.0 -1.238935629209052e-5 0.00047318043662328334 0.0 0.0
            0.011471086498795566 0.0 0.0 0.0 0.00047318043662328334 -0.068095270357532 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], zeros(8, 8),
           [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.002956756561550659 -0.42550604770774586
            0.0 0.0 0.0 0.0 0.0 0.0 -0.42550604770774586 -26.56156354297847], zeros(8, 8)]

    @test c.Ψ_p[2] ≈
          [[-0.010063967108852991 0.0 0.0 0.0 -2.1746114192313637e-5 0.0020611077566228815 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            -2.1746114192313637e-5 0.0 0.0 0.0 -4.085980403503138e-6 0.00012909043317795935 0.0 0.0
            0.0020611077566228815 0.0 0.0 0.0 0.00012909043317795935 -0.012235256872832223 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], zeros(8, 8), zeros(8, 8), zeros(8, 8)]
end

@testset "Evaluation Second Order into cache" begin
    m = @include_example_module(Examples.rbc_observables)
    p_f = (ρ=0.2, δ=0.02, σ=0.01, Ω_1=0.01)
    p_d = (α=0.5, β=0.95)
    p_d_symbols = collect(Symbol.(keys(p_d)))  #The order of derivatives in p_d
    c = SolverCache(m, Val(2), p_d)
    sol = generate_perturbation(m, p_d, p_f, Val(2); cache=c)
    # Create parameter vector in the same ordering the internal algorithms would

    @test c.y ≈ [5.936252888048733, 6.884057971014498]
    @test c.x ≈ [47.39025414828825, 0.0]
    @test c.H_yp ≈ [0.028377570562199098 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]
    @test c.H_y ≈ [-0.0283775705621991 0.0; 1.0 -1.0; 0.0 1.0; 0.0 0.0]
    @test c.H_xp ≈ [0.00012263591151906127 -0.011623494029190608
                    1.0 0.0
                    0.0 0.0
                    0.0 1.0]
    @test c.H_x ≈ [0.0 0.0
                   -0.98 0.0
                   -0.07263157894736837 -6.884057971014498
                   0.0 -0.2]
    @test c.Ψ[1] ≈
          [-0.009560768753410337 0.0 0.0 0.0 -2.0658808482697935e-5 0.0019580523687917364 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.009560768753410338 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           -2.0658808482697935e-5 0.0 0.0 0.0 -3.881681383327978e-6 0.00012263591151906127 0.0 0.0
           0.0019580523687917364 0.0 0.0 0.0 0.00012263591151906127 -0.011623494029190608 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0]
    @test c.Ψ[2] ≈ zeros(8, 8)
    @test c.Ψ[3] ≈ [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
           0.0 0.0 0.0 0.0 0.0 0.0 0.0007663134567721225 -0.07263157894736837
           0.0 0.0 0.0 0.0 0.0 0.0 -0.07263157894736837 -6.884057971014498]
    @test c.Ψ[4] ≈ zeros(8, 8)
    @test c.Γ ≈ [0.01]
    @test c.Ω ≈ [0.01, 0.01]

    # Second order checks!
    @test c.Ψ_yp[1] ≈
          [[0.0048317190660755365 0.0 0.0 0.0 6.960218465183541e-6 -0.0006596930439853133 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            6.960218465183541e-6 0.0 0.0 0.0 6.53894208439611e-7 -2.0658808482697952e-5 0.0 0.0
            -0.0006596930439853133 0.0 0.0 0.0 -2.0658808482697952e-5 0.0019580523687917372 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], zeros(8, 8), zeros(8, 8), zeros(8, 8)]
    @test c.Ψ_yp[2] ≈ [zeros(8, 8), zeros(8, 8), zeros(8, 8), zeros(8, 8)]

    @test c.Ψ_y[1] ≈ [[0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 -0.0048317190660755365 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], zeros(8, 8), zeros(8, 8), zeros(8, 8)]

    @test c.Ψ_y[2] ≈ [zeros(8, 8), zeros(8, 8), zeros(8, 8), zeros(8, 8)]

    @test c.Ψ_xp[1] ≈
          [[6.9602184651835404e-6 0.0 0.0 0.0 6.53894208439611e-7 -2.0658808482697952e-5 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            6.53894208439611e-7 0.0 0.0 0.0 2.0477213369556226e-7 -3.881681383327981e-6 0.0 0.0
            -2.0658808482697952e-5 0.0 0.0 0.0 -3.881681383327981e-6 0.00012263591151906138 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], zeros(8, 8), zeros(8, 8), zeros(8, 8)]
    @test c.Ψ_xp[2] ≈
          [[-0.0006596930439853132 0.0 0.0 0.0 -2.0658808482697952e-5 0.0019580523687917372 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            -2.0658808482697952e-5 0.0 0.0 0.0 -3.881681383327981e-6 0.00012263591151906138 0.0 0.0
            0.0019580523687917372 0.0 0.0 0.0 0.00012263591151906138 -0.011623494029190612 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], zeros(8, 8), zeros(8, 8), zeros(8, 8)]
    @test c.Ψ_x[1] ≈ [zeros(8, 8), zeros(8, 8),
           [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 -2.4255412970806024e-5 0.000766313456772123
            0.0 0.0 0.0 0.0 0.0 0.0 0.000766313456772123 -0.07263157894736838], zeros(8, 8)]
    @test c.Ψ_x[2] ≈ [zeros(8, 8), zeros(8, 8),
           [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.000766313456772123 -0.07263157894736838
            0.0 0.0 0.0 0.0 0.0 0.0 -0.07263157894736838 -6.884057971014497], zeros(8, 8)]

    # solution tests
    @test c.y ≈ [5.936252888048733, 6.884057971014498]
    @test c.x ≈ [47.39025414828824, 0.0]
    @test c.g_x ≈ [0.0957964300241661 0.6746869652586178
           0.07263157894736878 6.884057971014507]
    @test c.h_x ≈ [0.9568351489232028 6.209371005755889; -1.5076865909646354e-18 0.2]
    @test c.Σ ≈ [1e-4]
    @test c.Ω ≈ [0.01, 0.01]
    @test c.η == reshape([0; -1], 2, m.n_ϵ)
    @test c.B ≈ [0.0; -0.01]
    @test c.Q ≈ [1.0 0.0 0.0 0.0; 0.0 0.0 1.0 0.0]
    @test c.C_1 ≈ [0.0957964300241661 0.6746869652585828; 1.0 0.0]
    @test Array(c.V) ≈ [0.07005411173180227 0.00015997603451513485
           0.00015997603451513485 0.00010416666666666667]

    # Check the solution type matches these all
    fields_to_compare = (:y, :x, :g_x, :B, :Q, :η, :Γ, :g_σσ, :g_xx, :C_0, :C_1, :C_2)
    @test all_fields_equal(c, sol, fields_to_compare)
    # Some special cases
    @test 0.5 * c.h_σσ ≈ sol.A_0
    @test c.h_x ≈ sol.A_1
    @test 0.5 * c.h_xx ≈ sol.A_2
    @test c.Ω ≈ sol.D.σ
end

@testset "Evaluate 2nd Order Derivatives into cache" begin
    m = @include_example_module(Examples.rbc_observables)
    p_f = (ρ=0.2, δ=0.02, σ=0.01, Ω_1=0.01)
    p_d = (α=0.5, β=0.95)
    p_d_symbols = collect(Symbol.(keys(p_d)))  #The order of derivatives in p_d
    c = SolverCache(m, Val(2), p_d)
    sol = generate_perturbation(m, p_d, p_f, Val(2); cache=c)
    generate_perturbation_derivatives!(m, p_d, p_f, c)  # Solves and fills the cache
    @inferred generate_perturbation(m, p_d, p_f, Val(2); cache=c)
    @inferred generate_perturbation_derivatives!(m, p_d, p_f, c)

    @test c.H_x_p ≈ [[0.0 0.0
            0.0 0.0
            -0.4255060477077458 -26.561563542978472
            0.0 0.0], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]]
    @test c.H_yp_p ≈ [[0.011471086498795562 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0],
           [0.029871126907577997 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]]
    @test c.H_y_p ≈
          [[0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0], [0.0 0.0; 0.0 0.0; 0.0 0.0; 0.0 0.0]]
    @test c.H_xp_p ≈ [[0.000473180436623283 -0.06809527035753198
            0.0 0.0
            0.0 0.0
            0.0 0.0], [0.00012909043317795924 -0.01223525687283222
                       0.0 0.0
                       0.0 0.0
                       0.0 0.0]]
    @test c.Γ_p ≈ [[0.0], [0.0]]
    @test c.H_p ≈ [[-0.06809527035753199, 0.0, -26.561563542978472, 0.0],
           [-0.1773225633743801, 0.0, 0.0, 0.0]]
    @test c.Ω_p ≈ [[0.0, 0.0], [0.0, 0.0]]
    @test c.Ψ_p[1] ≈
          [[-0.0038647566790457792 0.0 0.0 0.0 -7.971028956261657e-5 0.011471086498795566 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            -7.971028956261657e-5 0.0 0.0 0.0 -1.238935629209052e-5 0.00047318043662328334 0.0 0.0
            0.011471086498795566 0.0 0.0 0.0 0.00047318043662328334 -0.068095270357532 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], zeros(8, 8),
           [0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.002956756561550659 -0.42550604770774586
            0.0 0.0 0.0 0.0 0.0 0.0 -0.42550604770774586 -26.56156354297847], zeros(8, 8)]
    @test c.Ψ_p[2] ≈
          [[-0.010063967108852991 0.0 0.0 0.0 -2.1746114192313637e-5 0.0020611077566228815 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            -2.1746114192313637e-5 0.0 0.0 0.0 -4.085980403503138e-6 0.00012909043317795935 0.0 0.0
            0.0020611077566228815 0.0 0.0 0.0 0.00012909043317795935 -0.012235256872832223 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
            0.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0], zeros(8, 8), zeros(8, 8), zeros(8, 8)]

    # solution
    @test hcat(c.y_p...) ≈ [
        55.78596896689701 76.10141579073955
        66.89124302798608 105.01995379122064
    ]
    @test hcat(c.x_p...) ≈ [555.2637030544529 1445.9269000240533; 0.0 0.0]
    @test c.g_x_p ≈ [
        [
            -0.12465264193058262 5.596211904442805
            -1.2823781479976832e-15 66.89124302798608
        ],
        [
            -1.6946742377792863 -0.8343618226192915
            -1.1080332409972313 105.01995379122064
        ],
    ]
    @test c.h_x_p ≈ [
        [0.12465264193058134 61.29503112354326; 0.0 0.0],
        [0.586640996782055 105.85431561383992; 0.0 0.0],
    ]
    @test c.Σ_p ≈ [[0.0], [0.0]]
    @test hcat(c.Ω_p...) ≈ [0.0 0.0; 0.0 0.0]
    @test c.B_p ≈ [[0.0; 0.0], [0.0; 0.0]]
    @test c.C_1_p ≈ [
        [-0.12465264193057919 5.596211904442171; 0.0 0.0],
        [-1.6946742377792825 -0.8343618226202246; 0.0 0.0],
    ]
    @test c.V_p ≈ [
        [1.584528257999749 0.0015841155991973127; 0.0015841155991973127 0.0],
        [3.336643488330957 0.002750404724942799; 0.002750404724942799 0.0],
    ]
end

@testset "Dense RBC 2nd Order, sigma derivatives" begin
        m = @include_example_module(Examples.rbc)
        p_f = (ρ=0.2, δ=0.02)
        p_d = (α=0.5, β=0.95, σ=0.01)
        c = SolverCache(m, Val(2), p_d)
        sol = generate_perturbation(m, p_d, p_f, Val(2); cache=c)
        generate_perturbation_derivatives!(m, p_d, p_f, c)  # Solves and fills the cache
        @inferred generate_perturbation(m, p_d, p_f, Val(2); cache=c)
        @inferred generate_perturbation_derivatives!(m, p_d, p_f, c)

        @test c.g_σσ_p ≈
        [0.001363945590429837 0.0035331253439556264 0.03129961925086118; 0.0 0.0 0.0]
        @test c.h_σσ_p ≈
        [-0.001363945590429837 -0.0035331253439556264 -0.03129961925086118; 0.0 0.0 0.0]
end

@testset "RBC second order" begin
        m = @include_example_module(Examples.rbc)
        p_f = (ρ=0.2, δ=0.02, σ=0.01)
        p_d = (α=0.5, β=0.95)
        c = SolverCache(m, Val(2), p_d)
        sol = generate_perturbation(m, p_d, p_f, Val(2); cache=c)
        generate_perturbation_derivatives!(m, p_d, p_f, c)  # Solves and fills the cache
        @inferred generate_perturbation(m, p_d, p_f, Val(2); cache=c)
        @inferred generate_perturbation_derivatives!(m, p_d, p_f, c)

    # sol tests
    @test c.y ≈ [5.936252888048733, 6.884057971014498]
    @test c.x ≈ [47.39025414828824, 0.0]
    @test hcat(c.y_p...) ≈ [
        55.78596896689701 76.10141579073955
        66.89124302798608 105.01995379122064
    ]
    @test hcat(c.x_p...) ≈ [555.2637030544529 1445.9269000240533; 0.0 0.0]
    @test c.g_x ≈ [
        0.0957964300241661 0.6746869652586178
        0.07263157894736878 6.884057971014507
    ]
    @test c.h_x ≈ [0.9568351489232028 6.209371005755889; -1.5076865909646354e-18 0.2]
    @test c.g_x_p ≈ [
        [
            -0.12465264193058262 5.596211904442805
            -1.2823781479976832e-15 66.89124302798608
        ],
        [
            -1.6946742377792863 -0.8343618226192915
            -1.1080332409972313 105.01995379122064
        ],
    ]
    @test c.h_x_p ≈ [
        [0.12465264193058134 61.29503112354326; 0.0 0.0],
        [0.586640996782055 105.85431561383992; 0.0 0.0],
    ]
    @test c.Σ ≈ [1e-4]
    @test c.Σ_p ≈ [[0.0], [0.0]]

    @test c.g_xx[:, :, 1] ≈ [
        -0.000371083339499955 0.005130472630563616
        -0.0007663134567721218 0.07263157894736846
    ]
    @test c.g_xx[:, :, 2] ≈
          [0.005130472630563628 0.6265410073784347; 0.07263157894736846 6.8840579710144985]
    @test c.h_xx[:, :, 1] ≈ [-0.0003952301172721667 0.06750110631680481; 0.0 0.0]
    @test c.h_xx[:, :, 2] ≈ [0.06750110631680481 6.257516963636062; 0.0 0.0]
    @test c.g_σσ ≈ [0.0001564980962543059, 0.0]
    @test c.h_σσ ≈ [-0.0001564980962543059, 0.0]

    @test c.g_xx_p[1] ≈ cat(
        [
            0.005919879353027914 0.0028179768532923975
            0.010511393863734047 1.0255059778850425e-15
        ],
        [
            0.0028179768532926456 5.321141188794126
            1.0255059778850425e-15 66.89124302798595
        ],
        dims = 3,
    )
    @test c.g_xx_p[2] ≈ cat(
        [
            0.017520602529482888 -0.1729419353982087
            0.03507155408568066 -1.1080332409972282
        ],
        [
            -0.17294193539820824 -0.804951434933394
            -1.1080332409972282 105.01995379122047
        ],
        dims = 3,
    )
    @test c.h_xx_p[1] ≈ cat(
        [0.0045915145107061316 -0.0028179768532913723; 0.0 0.0],
        [-0.0028179768532916203 61.570101839191814; 0.0 0.0],
        dims = 3,
    )
    @test c.h_xx_p[2] ≈ cat(
        [0.017550951556197774 -0.9350913055990194; 0.0 0.0],
        [-0.9350913055990198 105.82490522615385; 0.0 0.0],
        dims = 3,
    )
    @test c.g_σσ_p ≈ [0.001363945590429837 0.0035331253439556264; 0.0 0.0]
    @test c.h_σσ_p ≈ [-0.001363945590429837 -0.0035331253439556264; 0.0 0.0]

    @test c.Ω === nothing
    @test c.Ω_p === nothing
    @test c.η == reshape([0; -1], 2, m.n_ϵ)
    @test c.Q == I
    @test c.y ≈ sol.y
    @test c.x ≈ sol.x
    @test c.g_x ≈ sol.g_x
    @test c.h_x ≈ sol.A_1
    @test c.B ≈ sol.B
    @test c.Ω === sol.D  # should be nothing
    @test c.Q ≈ sol.Q
    @test c.η ≈ sol.η
    @test c.g_σσ ≈ sol.g_σσ
    @test c.h_σσ ≈ sol.A_0 * 2
    @test c.g_xx ≈ sol.g_xx
    @test c.h_xx ≈ sol.A_2 * 2
    @test sol.retcode == :Success
end

@testset "Second Order Pullback inference" begin
    m = @include_example_module(Examples.rbc_observables)
    p_f = (ρ=0.2, δ=0.02, σ=0.01, Ω_1=0.01)
    p_d = (α=0.5, β=0.95)

    c = SolverCache(m, Val(2), p_d)
    sol = generate_perturbation(m, p_d, p_f, Val(2); cache = c)

    _, pb = Zygote.pullback(generate_perturbation, m, p_d, p_f, Val(2))
    @inferred pb(sol)
end

ext: jl


text:
using Statistics
using OffsetArrays
using JLD2
using Plots
using LaTeXStrings
pgfplotsx()

const L = 20
MCS = 230_000

const Δϕ = 10 #°
# const T = 2
const ξ = 20.0
const n₀ = 1.5
const nₑ = 1.7

# const P_func(β) = (3 * cosd(β)^2 -1)/2
const P₂ = OffsetArray(map(β ->  (3 * cosd(β)^2 -1)/2, -180:180), -180:180)
Iₙ = collect(2:L+1)
Iₚ = collect(0:L-1)
Iₙ[20] = 1
Iₚ[1] = L

const randoms_Δϕ = Int16.(cat(-Δϕ÷2:-1, 1:Δϕ÷2, dims=1))

n_eff(φ::Number, n₀ = n₀, nₑ = nₑ)::Float64 = n₀ * nₑ / √(n₀^2*cosd(φ)^2 + nₑ^2*sind(φ)^2)

function nlc(E::Number, MCS::Integer=230_000, L::Integer=20, skipfirst::Integer=30_000, probe_every::Integer=100, P₂=P₂, Iₙ=Iₙ, Iₚ=Iₚ, ξ=ξ, randoms_Δϕ=randoms_Δϕ)::Float64
    #ordered initial conditions - all ϕ=0°
    ϕ = zeros(Int16,L,L)
    n_eff_sum = 0.0

    for k ∈ 1:MCS
        for j ∈ 1:L, i ∈ 2:L-1
            #metropolis algorithm
            @inbounds ϕ_new::Int16 = ϕ[i,j] + rand(randoms_Δϕ)
            ϕ_new > 90  && (ϕ_new -= 180) # a && b  => if a then run b
            ϕ_new < -90 && (ϕ_new += 180)

            @inbounds U_old::Float64 = -ξ * (P₂[ϕ[i,j] - ϕ[Iₙ[i],j]] + P₂[ϕ[i,j] - ϕ[Iₚ[i],j]] + P₂[ϕ[i,j] - ϕ[i,Iₙ[j]]] + P₂[ϕ[i,j] - ϕ[i,Iₚ[j]]]) - E^2 * P₂[90-ϕ[i,j]]

            @inbounds U_new::Float64 = -ξ * (P₂[ϕ_new - ϕ[Iₙ[i],j]] + P₂[ϕ_new - ϕ[Iₚ[i],j]] + P₂[ϕ_new - ϕ[i,Iₙ[j]]] + P₂[ϕ_new - ϕ[i,Iₚ[j]]]) - E^2 * P₂[90-ϕ_new]

            ΔU::Float64 = U_new - U_old
            if ΔU < 0 || rand() ≤ exp(-ΔU)
                @inbounds ϕ[i,j] = ϕ_new
            end
        end
        if k > skipfirst && k%probe_every == 0
            n_eff_sum += sum(n_eff, ϕ)/L^2
        end
    end
    return n_eff_sum / ((MCS-skipfirst)÷probe_every) #return mean <n_eff>
    #return ϕ
end

function nlc_cos2_profile(E::Number, MCS::Integer=230_000, L::Integer=20, skipfirst::Integer=30_000, probe_every::Integer=100, P₂=P₂, Iₙ=Iₙ, Iₚ=Iₚ, ξ=ξ, randoms_Δϕ=randoms_Δϕ)::Array{Float64}
    #ordered initial conditions - all ϕ=0°
    ϕ = zeros(Int16,L,L)
    cos2_profile = zeros(L)

    for k ∈ 1:MCS
        for j ∈ 1:L, i ∈ 2:L-1
            #metropolis algorithm
            @inbounds ϕ_new::Int16 = ϕ[i,j] + rand(randoms_Δϕ)
            ϕ_new > 90  && (ϕ_new -= 180) # a && b  => if a then run b
            ϕ_new < -90 && (ϕ_new += 180)

            @inbounds U_old::Float64 = -ξ * (P₂[ϕ[i,j] - ϕ[Iₙ[i],j]] + P₂[ϕ[i,j] - ϕ[Iₚ[i],j]] + P₂[ϕ[i,j] - ϕ[i,Iₙ[j]]] + P₂[ϕ[i,j] - ϕ[i,Iₚ[j]]]) - E^2 * P₂[90-ϕ[i,j]]

            @inbounds U_new::Float64 = -ξ * (P₂[ϕ_new - ϕ[Iₙ[i],j]] + P₂[ϕ_new - ϕ[Iₚ[i],j]] + P₂[ϕ_new - ϕ[i,Iₙ[j]]] + P₂[ϕ_new - ϕ[i,Iₚ[j]]]) - E^2 * P₂[90-ϕ_new]

            ΔU::Float64 = U_new - U_old
            if ΔU < 0 || rand() ≤ exp(-ΔU)
                @inbounds ϕ[i,j] = ϕ_new
            end
        end
        if k > skipfirst && k%probe_every == 0
            cos2_profile += mean(x -> cosd(x)^2, ϕ, dims=2)
        end
    end
    return cos2_profile ./ ((MCS-skipfirst)÷probe_every)
end

function run_nlc(Es)
    result = Dict{Float64,Float64}()
    Threads.@threads for E ∈ Es
        result[E] = @time nlc(E)
    end
    return result
end

Es = 0:0.1:3.5

result_ret = @time run_nlc(Es)
@save "result_002.jld2" result_ret
@load "result_002.jld2" result_ret

# serialize("result_002.jls", result)

scatter(result_ret, markershape=:auto, markersize=2, xlabel="\$E^*\$", ylabel="\$<n_{eff}>\$",
        title="Effective refracting index as a function of reduced external electric field - ordered initial conditions",
             titlefontsize=8, markerstrokealpha=0.0, label="L=20",  legend=true)
vline!([0.73], alpha=0.7, label=L"E_F^*")
lens!([0.65,0.8],[1.68,1.7], minorticks=5, minorgrid=true, legend=false, inset = (1, bbox(0.6, 0.1, 0.3, 0.5)))
savefig("neff_E.tex")

#@save "result.jld2" result_ret
#@load "result.jld2"

using Profile
using StatProfilerHTML

# Profile.clear()
# result_ret = @@profilehtml run_nlc(Es)

cos2_profile = nlc_cos2_profile(1.2*0.73)
scatter(cos2_profile, grid=:both, xticks=1:20, xlabel="z", ylabel="<\$\\cos^2(\\phi)\$>", legend=:none,
        title="Ensemble averaged particles' orientation profile")
savefig("cos2_profile.tex")


text:
push!(ARGS, "../../input_files/dynamic/basin_switching/8GPa/switch2.dat")
include("../../Basin.jl")


text:
import CSV

devices = Dict{Symbol, DeviceModel}(
    :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
    :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
)
branches = Dict{Symbol, DeviceModel}(
    :L => DeviceModel(Line, StaticLine),
    :T => DeviceModel(Transformer2W, StaticTransformer),
    :TT => DeviceModel(TapTransformer, StaticTransformer),
)
services = Dict{Symbol, ServiceModel}()

function get_deserialized(op_problem::OperationsProblem; kwargs...)
    orig = pwd()
    path = mktempdir()
    cd(path)

    try
        filename = "test_op_problem.bin"
        PSI.serialize(op_problem, filename)
        return OperationsProblem(filename; kwargs...)
    finally
        cd(orig)
    end
end

@testset "Solving ED with CopperPlate" begin
    template = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
    parameters_value = [true, false]
    c_sys5 = build_system("c_sys5")
    c_sys14 = build_system("c_sys14")
    systems = [c_sys5, c_sys14]
    test_results = Dict{System, Float64}(c_sys5 => 240000.0, c_sys14 => 142000.0)
    @info "Test solve ED with CopperPlatePowerModel network"
    for sys in systems, p in parameters_value
        @testset "ED CopperPlatePowerModel model use_parameters = $(p)" begin
            ED = OperationsProblem(
                TestOpProblem,
                template,
                sys;
                optimizer = OSQP_optimizer,
                use_parameters = p,
            )
            psi_checksolve_test(ED, [MOI.OPTIMAL], test_results[sys], 10000)

            # Serialize, deserialize, rebuild, re-run.
            ED2 = get_deserialized(ED; optimizer = OSQP_optimizer)
            psi_checksolve_test(ED2, [MOI.OPTIMAL], test_results[sys], 10000)
        end
    end
    c_sys5_re = build_system("c_sys5_re")
    ED = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5_re;
        optimizer = GLPK_optimizer,
        slack_variables = true,
    )
    psi_checksolve_test(ED, [MOI.OPTIMAL], 240000.0, 10000)
end

@testset "Solving ED with PTDF Models" begin
    template = OperationsProblemTemplate(StandardPTDFModel, devices, branches, services)
    parameters_value = [true, false]
    c_sys5 = build_system("c_sys5")
    c_sys14 = build_system("c_sys14")
    c_sys14_dc = build_system("c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    PTDF_ref = Dict{UUIDs.UUID, PTDF}(
        IS.get_uuid(c_sys5) => build_PTDF5(),
        IS.get_uuid(c_sys14) => build_PTDF14(),
        IS.get_uuid(c_sys14_dc) => build_PTDF14_dc(),
    )
    test_results = Dict{UUIDs.UUID, Float64}(
        IS.get_uuid(c_sys5) => 340000.0,
        IS.get_uuid(c_sys14) => 142000.0,
        IS.get_uuid(c_sys14_dc) => 142000.0,
    )

    @info "Test solve ED with StandardPTDFModel network"
    for sys in systems, p in parameters_value
        @testset "ED StandardPTDFModel model use_parameters = $(p)" begin
            ED = OperationsProblem(
                TestOpProblem,
                template,
                sys;
                optimizer = OSQP_optimizer,
                use_parameters = p,
                PTDF = PTDF_ref[IS.get_uuid(sys)],
            )
            psi_checksolve_test(ED, [MOI.OPTIMAL], test_results[IS.get_uuid(sys)], 10000)
        end
    end
end

@testset "Solving ED With PowerModels with loss-less convex models" begin
    c_sys5 = build_system("c_sys5")
    c_sys14 = build_system("c_sys14")
    c_sys14_dc = build_system("c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    parameters_value = [true, false]
    networks = [DCPPowerModel, NFAPowerModel]
    test_results = Dict{System, Float64}(
        c_sys5 => 330000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )

    for net in networks, p in parameters_value, sys in systems
        @info("Test solve ED with $(net) network")
        @testset "ED model $(net) and use_parameters = $(p)" begin
            template = OperationsProblemTemplate(net, devices, branches, services)
            ED = OperationsProblem(
                TestOpProblem,
                template,
                sys;
                optimizer = ipopt_optimizer,
                use_parameters = p,
            )
            #The tolerance range here is large because NFA has a much lower objective value
            psi_checksolve_test(
                ED,
                [MOI.OPTIMAL, MOI.LOCALLY_SOLVED],
                test_results[sys],
                35000,
            )
        end
    end

end

@testset "Solving ED With PowerModels with linear convex models" begin
    c_sys5 = build_system("c_sys5")
    c_sys14 = build_system("c_sys14")
    c_sys14_dc = build_system("c_sys14_dc")
    systems = [c_sys5, c_sys14]
    parameters_value = [true, false]
    networks = [DCPLLPowerModel, LPACCPowerModel]
    test_results = IdDict{System, Float64}(
        c_sys5 => 340000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )

    for net in networks, p in parameters_value, sys in systems
        @info("Test solve ED with $(net) network")
        @testset "ED model $(net) and use_parameters = $(p)" begin
            template = OperationsProblemTemplate(net, devices, branches, services)
            ED = OperationsProblem(
                TestOpProblem,
                template,
                sys;
                optimizer = ipopt_optimizer,
                use_parameters = p,
            )
            #The tolerance range here is large because NFA has a much lower objective value
            psi_checksolve_test(
                ED,
                [MOI.OPTIMAL, MOI.LOCALLY_SOLVED],
                test_results[sys],
                10000,
            )

        end
    end

end

@testset "Operation Model Constructors with Slacks" begin

    networks = [StandardPTDFModel, DCPPowerModel, ACPPowerModel]

    thermal_gens = [ThermalDispatch]

    c_sys5_re = build_system("c_sys5_re")
    systems = [c_sys5_re]
    for net in networks, thermal in thermal_gens, system in systems
        devices = Dict{Symbol, DeviceModel}(
            :Generators => DeviceModel(ThermalStandard, thermal),
            :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
            :RE => DeviceModel(RenewableDispatch, FixedOutput),
        )
        branches = Dict{Symbol, DeviceModel}(:L => DeviceModel(Line, StaticLine))
        template = OperationsProblemTemplate(net, devices, branches, services)
        op_problem = OperationsProblem(
            TestOpProblem,
            template,
            system;
            slack_variables = true,
            optimizer = ipopt_optimizer,
            PTDF = build_PTDF5(),
        )
        res = solve!(op_problem)
        @test termination_status(op_problem.psi_container.JuMPmodel) in
              [MOI.OPTIMAL, MOI.LOCALLY_SOLVED]
    end
end
#=
@testset "Solving ED With PowerModels with convex SOC and QC models" begin
    systems = [c_sys5, c_sys14]
    parameters_value = [true, false]
    networks = [SOCWRPowerModel,
                 QCRMPowerModel,
                 QCLSPowerModel,]
    test_results = Dict{System, Float64}(c_sys5 => 320000.0,
                                             c_sys14 => 142000.0)
    for  net in networks, p in parameters_value, sys in systems
        @info("Test solve ED with $(net) network")
        @testset "ED model $(net) and use_parameters = $(p)" begin
        template = OperationsProblemTemplate(net, devices, branches, services);
        ED = OperationsProblem(TestOpProblem, template, sys; optimizer = ipopt_optimizer, use_parameters = p);
        #The tolerance range here is large because Relaxations have a lower objective value
        psi_checksolve_test(ED, [MOI.OPTIMAL, MOI.LOCALLY_SOLVED], test_results[sys], 25000)
        end
    end
end
=#

@testset "Solving ED With PowerModels Non-Convex Networks" begin
    c_sys5 = build_system("c_sys5")
    c_sys14 = build_system("c_sys14")
    c_sys14_dc = build_system("c_sys14_dc")
    systems = [c_sys5, c_sys14, c_sys14_dc]
    parameters_value = [true, false]
    networks = [
        ACPPowerModel,
        #ACRPowerModel,
        ACTPowerModel,
    ]
    test_results = Dict{System, Float64}(
        c_sys5 => 340000.0,
        c_sys14 => 142000.0,
        c_sys14_dc => 142000.0,
    )

    for net in networks, p in parameters_value, sys in systems
        @info("Test solve ED with $(net) network")
        @testset "ED model $(net) and use_parameters = $(p)" begin
            template = OperationsProblemTemplate(net, devices, branches, services)
            ED = OperationsProblem(
                TestOpProblem,
                template,
                sys;
                optimizer = ipopt_optimizer,
                use_parameters = p,
            )
            psi_checksolve_test(
                ED,
                [MOI.OPTIMAL, MOI.LOCALLY_SOLVED],
                test_results[sys],
                10000,
            )
        end
    end

end

@testset "Solving UC Linear Networks" begin
    devices = Dict{Symbol, DeviceModel}(
        :Generators => DeviceModel(ThermalStandard, ThermalStandardUnitCommitment),
        :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
    )
    c_sys5 = build_system("c_sys5")
    c_sys5_dc = build_system("c_sys5_dc")
    parameters_value = [true, false]
    systems = [c_sys5, c_sys5_dc]
    networks = [DCPPowerModel, NFAPowerModel, StandardPTDFModel, CopperPlatePowerModel]
    PTDF_ref = Dict{UUIDs.UUID, PTDF}(
        IS.get_uuid(c_sys5) => build_PTDF5(),
        IS.get_uuid(c_sys5_dc) => build_PTDF5_dc(),
    )

    for net in networks, p in parameters_value, sys in systems
        @info("Test solve UC with $(net) network")
        @testset "UC model $(net) and use_parameters = $(p)" begin
            template = OperationsProblemTemplate(net, devices, branches, services)
            UC = OperationsProblem(
                TestOpProblem,
                template,
                sys;
                optimizer = GLPK_optimizer,
                use_parameters = p,
                PTDF = PTDF_ref[IS.get_uuid(sys)],
            )
            psi_checksolve_test(UC, [MOI.OPTIMAL, MOI.LOCALLY_SOLVED], 340000, 100000)
        end
    end
end
################################################################

function test_write_functions(file_path, op_problem, res)
    @testset "Test write optimizer problem" begin
        path = mkdir(joinpath(file_path, "op_problem"))
        file = joinpath(path, "op_problem.json")
        PSI.export_op_model(op_problem, file)
        PSI.write_data(op_problem, path)
        list = sort!(collect(readdir(path)))
        @test ["P__ThermalStandard.feather", "op_problem.json"] == list
    end

    @testset "Test write_data functions" begin
        PSI.write_data(get_variables(res), mkdir(joinpath(file_path, "one")))
        readdir(joinpath(file_path, "one"))
        for (k, v) in get_variables(res)
            @test isfile(joinpath(file_path, "one", "$k.feather"))
        end

        PSI.write_data(
            get_variables(res),
            res.time_stamp,
            mkdir(joinpath(file_path, "two"));
            file_type = CSV,
        )
        for (k, v) in get_variables(res)
            @test isfile(joinpath(file_path, "two/$k.csv"))
        end

        PSI.write_data(
            get_variables(res),
            res.time_stamp,
            mkdir(joinpath(file_path, "three")),
        )
        for (k, v) in get_variables(res)
            @test isfile(joinpath(file_path, "three", "$k.feather"))
        end

        var_name = PSI.variable_name(PSI.ACTIVE_POWER, PSY.ThermalStandard)
        PSI.write_data(
            get_variables(res)[var_name],
            mkdir(joinpath(file_path, "four")),
            string(var_name),
        )
        @test isfile(joinpath(file_path, "four", "$(var_name).feather"))

        #testing if directory is a file
        PSI.write_data(
            get_variables(res)[var_name],
            joinpath(file_path, "four", "$(var_name).feather"),
            string(var_name),
        )
        @test isfile(joinpath(file_path, "four", "$(var_name).feather"))

        PSI.write_optimizer_log(get_optimizer_log(res), mkdir(joinpath(file_path, "five")))
        @test isfile(joinpath(file_path, "five", "optimizer_log.json"))

        PSI.write_to_CSV(res, mkdir(joinpath(file_path, "six")))
        @test !isempty(joinpath(file_path, "six", "results"))
    end

    @testset "Test write result functions" begin
        new_path = joinpath(file_path, "seven")
        IS.write_results(res, mkdir(new_path))
        @test !isempty(new_path)
    end

    @testset "Test parameter values" begin
        c_sys5_re = build_system("c_sys5_re")
        system = op_problem.sys
        params =
            PSI.get_parameter_array(op_problem.psi_container.parameters[:P__get_maxactivepower__PowerLoad])
        params = PSI.axis_array_to_dataframe(params)
        devices = collect(PSY.get_components(PSY.PowerLoad, c_sys5_re))
        multiplier = [PSY.get_activepower(devices[1])]
        for d in 2:length(devices)
            multiplier = hcat(multiplier, PSY.get_activepower(devices[d]))
        end
        extracted = -multiplier .* params
        @test IS.get_base_power(res) .* extracted == res.parameter_values[:P_PowerLoad]
    end

    @testset "Set optimizer at solve call" begin
        c_sys5 = build_system("c_sys5")
        devices = Dict{Symbol, DeviceModel}(
            :Generators => DeviceModel(ThermalStandard, ThermalStandardUnitCommitment),
            :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
        )
        template = OperationsProblemTemplate(DCPPowerModel, devices, branches, services)
        UC = OperationsProblem(TestOpProblem, template, c_sys5;)
        set_services_template!(
            UC,
            Dict(:Reserve => ServiceModel(VariableReserve{ReserveUp}, RangeReserve)),
        )
        res = solve!(UC; optimizer = GLPK_optimizer)
        @test isapprox(get_total_cost(res)[:OBJECTIVE_FUNCTION], 340000.0; atol = 100000.0)
    end

    @testset "Test get_variable function" begin
        variable = PSI.get_result_variable(res, :P, ThermalStandard)
        @test isa(variable, DataFrames.DataFrame)
    end
end

@testset "Miscellaneous OperationsProblem" begin
    duals = [:CopperPlateBalance]
    template = OperationsProblemTemplate(CopperPlatePowerModel, devices, branches, services)
    c_sys5_re = build_system("c_sys5_re")
    op_problem = OperationsProblem(
        TestOpProblem,
        template,
        c_sys5_re;
        optimizer = OSQP_optimizer,
        use_parameters = true,
        constraint_duals = duals,
    )
    res = solve!(op_problem)
    @testset "Test print methods" begin
        list = [template, op_problem, op_problem.psi_container, res, services]
        _test_plain_print_methods(list)
        list = [services]
        _test_html_print_methods(list)
    end

    @testset "test constraint duals in the operations problem" begin
        name = PSI.constraint_name("CopperPlateBalance")
        for i in 1:ncol(get_time_stamp(res))
            dual = JuMP.dual(op_problem.psi_container.constraints[name][i])
            @test isapprox(dual, get_duals(res)[name][i, 1])
        end
        dual_results = get_dual_values(op_problem.psi_container, duals)
        @test dual_results == res.dual_values
    end

    @testset "test get variable function" begin
        @test_throws IS.ConflictingInputsError PSI.get_variable(res, :fake)
        @test res.variable_values[:P__ThermalStandard] ==
              PSI.get_variable(res, :P__ThermalStandard)
    end

    @testset "Test writing functions" begin
        path = joinpath(pwd(), "test_writing")
        try
            !isdir(path) && mkdir(path)
            test_write_functions(path, op_problem, res)
        finally
            @info("removing test files")
            rm(path, recursive = true)
        end
    end
end


text:
using LinearAlgebra, Optim

export tracenorm, tracedist, frobdist, isunitary, ishermitian_tol,
        isspecial, offdiag, jacobi_simultaneous_diag

ϵ = sqrt(eps(Float64))

"""
Compute the trace norm of a matrix ``A``
"""
function tracenorm(A::Matrix{T}) where T<:Number
    #Compute the trace norm of a matrix A
    if ishermitian(A)
        return sum(abs.(eigvals(A)))
    else
        return sum(svd(A).S)
    end
end

"""
Compute the trace distance between matrices ``A``, ``B``.
"""
function tracedist(A::Matrix{T}, B::Matrix{T}) where T<:Number
    #Compute trace distance between matrices A, B
    return 0.5*tracenorm(A-B)
end

"""
Compute the Frobenius distance between matrices ``A``, ``B``.
"""
function frobdist(A::Matrix{T}, B::Matrix{T}) where T<:Number
    return norm(A-B, p=2)
end

"""
Check if a matrix is unitary up to a tolerance ``\epsilon``.
"""
function isunitary(A::Matrix{T}; tol=ϵ) where T <: Number
    return tracenorm(A'A - I) < tol
end

"""
Check if a matrix is Hermitian up to a tolerance ``\epsilon``.
"""
function ishermitian_tol(A::Matrix{T}; tol=ϵ) where T <: Number
    return tracedist(A, Matrix(A')) < tol
end

"""
Check if a matrix is special up to a tolerance ``\epsilon``.
"""
function isspecial(A::Matrix{T}; tol=ϵ) where T <: Number
    return isapprox(det(A), 1.0, rtol=tol)
end

"""
Sum of the off-diagonal elements of a matrix.
"""
function offdiag(A::Matrix{T}) where T <: Number
    return sum(tril(A,-1).^2) + sum(triu(A,1).^2)
end

"""
Function which implements the Jacobi algorithm for simultaneous (possibly approxiamate)
diagonalization of two real, symmetric matrices A, B.

Algorithm is from:

V. Kuleshov, A. Chaganty, P. Liang. "Simultaneous diagonalization: the asymmetric,
low-rank, and noisy settings" arXiv: 1051.06318 (2015)

J.-F. Cardoso, A. Soulomiac. "Jacobi angles for simultaneous diagonalization"
SIAM J. Matrix Anal. Appl. 17(1), 161-164. (1996)

"""
function jacobi_simultaneous_diag(A::Matrix{T}, B::Matrix{T};
        tol=sqrt(eps(T)), maxiter=1000) where T<:Real

    @assert ishermitian_tol(A) && ishermitian_tol(B) "A, B must be symmetric!"

    (m,mm) = size(A)
    done = false

    U = Float64.(Matrix(I, (m,mm)))
    Λa = A
    Λb = B

    #Cost function for minimization of "off-diagonal-ness"
    #TODO: I was not able to succesfully implement the exact solution for the minimizer from the
    #Cardoso paper...

    cost(X,Y) = offdiag(X)+offdiag(Y)

    start_cost = cost(A,B)

    function minims_lsq(X, Y, i, j)
        function f(x)
            g = LinearAlgebra.Givens(i,j, cos(x), sin(x))
            return cost(g'*X*g, g'*Y*g)
        end
        xopt = Optim.minimizer(optimize(f, -π/4, π/4)) #are these the correct limits?
        return (cos(xopt), sin(xopt))
    end

    target = tol*(norm(A,2)+norm(B,2))

    n = 1
    while ((cost(Λa,Λb) > tol) && (n < maxiter))
        for j=1:m-1
            for k=j+1:m
                (c, s) = minims_lsq(Λa, Λb, j, k)
                Γ = LinearAlgebra.Givens(j, k, c, s)
                U = U*Γ
                Λa = Γ'*Λa*Γ
                Λb = Γ'*Λb*Γ
            end
        end
    end

    final_cost = cost(Λa,Λb)

    success = true

    if ((n >= maxiter) || (final_cost > tol))
        @warn "Was unable to find a simultaneous diagonalization: start
                cost function: $(start_cost), final: $(final_cost)."
        success = false
    end

    return (U, Diagonal(Λa), Diagonal(Λb), success)
end


text:
# Hello world example, similar to the Boost.Python hello world

using CxxWrap
using Base.Test

# Wrap the functions defined in C++
wrap_modules(CxxWrap._l_inheritance)

using CppInheritance

b = B()
c = C()
d = D()

@test message(b) == "B"
@test message(c) == "C"
@test message(d) == "D"

@test take_ref(b) == "B"
@test take_ref(c) == "C"
@test take_ref(d) == "D"

# factory function returning an abstract type A
@test message(create_abstract()) == "B"

@test dynamic_message_c(c) == "C"

# shared ptr variants
b_ptr = shared_b()
c_ptr = shared_c()
d_ptr = shared_d()

@test shared_ptr_message(b_ptr) == "B"
@test shared_ptr_message(c_ptr) == "C"
@test shared_ptr_message(d_ptr) == "D"

@test message(b_ptr) == "B"
@test message(c_ptr) == "C"
@test message(d_ptr) == "D"

@test weak_ptr_message_b(b_ptr) == "B"
@test weak_ptr_message_a(b_ptr) == "B"
@test weak_ptr_message_a(c_ptr) == "C"
@test weak_ptr_message_a(d_ptr) == "D"

a = VirtualSolver.E()
VirtualSolver.solve(a)

b = VirtualSolver.F(safe_cfunction(x -> 2x, Float64, (Float64,)))
VirtualSolver.solve(b)


text:
@testset "Distances" begin
  p = P2(0,1)
  l = Line(P2(0,0), P2(1,0))
  @test evaluate(Euclidean(), p, l) == T(1)
  @test evaluate(Euclidean(), l, p) == T(1)

  p1, p2 = P2(1,0), P2(0,1)
  @test evaluate(Chebyshev(), p1, p2) == T(1)
end


text:
module LossFunctionsTests

using LossFunctions
using CategoricalArrays
using DualNumbers
using Statistics
using Random
using Test

tests = [
    "tst_loss.jl",
    "tst_api.jl",
    "tst_properties.jl"
]

perf = [
]

# for deterministic testing
Random.seed!(1234)

distance_losses = [
    L2DistLoss(), LPDistLoss(2.0), L1DistLoss(), LPDistLoss(1.0),
    LPDistLoss(0.5), LPDistLoss(1.5), LPDistLoss(3),
    LogitDistLoss(), L1EpsilonInsLoss(0.5), EpsilonInsLoss(1.5),
    L2EpsilonInsLoss(0.5), L2EpsilonInsLoss(1.5),
    PeriodicLoss(1), PeriodicLoss(1.5), HuberLoss(1),
    HuberLoss(1.5), QuantileLoss(.2), QuantileLoss(.5),
    QuantileLoss(.8), LogCoshLoss()
]

margin_losses = [
    LogitMarginLoss(), L1HingeLoss(), L2HingeLoss(),
    PerceptronLoss(), SmoothedL1HingeLoss(.5),
    SmoothedL1HingeLoss(1), SmoothedL1HingeLoss(2),
    ModifiedHuberLoss(), ZeroOneLoss(), L2MarginLoss(),
    ExpLoss(), SigmoidLoss(), DWDMarginLoss(.5),
    DWDMarginLoss(1), DWDMarginLoss(2)
]

other_losses = [
    MisclassLoss(), PoissonLoss(), CrossEntropyLoss()
]

for t in tests
    @testset "$t" begin
        include(t)
    end
end

for p in perf
    @testset "$p" begin
        include(p)
    end
end

end # module


text:
function integrate(f, a, b, n)
	h = (b - a) / n
	x = zeros(1, 3)
	x[1] = f(a) + f(b)
	for i ∈ 1:(n - 1)
		X = a + i*h
		if iseven(i)
			x[3] += f(X)
		else
			x[2] += f(X)
		end
	end
	return h*(x[1] + 2*x[3] + 4*x[2])/3
end

f(x) = x
println(integrate(f, 0, 3, 100))