text:
using RegistryTools: DEFAULT_REGISTRY_URL,
    parse_registry,
    showsafe,
    registration_branch,
    get_registry,
    gitcmd
using LibGit2
using Pkg.TOML
using Pkg.Types: Project

using Test

const TEST_GITCONFIG = Dict(
    "user.name" => "RegistratorTests",
    "user.email" => "ci@juliacomputing.com",
)
const TEST_SIGNATURE = LibGit2.Signature(
    TEST_GITCONFIG["user.name"],
    TEST_GITCONFIG["user.email"],
)

@testset "RegistryTools" begin

@testset "Utilities" begin
    @testset "showsafe" begin
        @test string(showsafe(3)) == "3"
        @test string(showsafe(nothing)) == "nothing"
    end

    @testset "registration_branch" begin
        example = Project(Dict(
            "name" => "Example", "version" => "1.10.2",
            "uuid" => "698ec630-83b2-4a6d-81d4-a10176273030"
        ))
        @test registration_branch(example) == "registrator/example/698ec630/v1.10.2"
    end
end

@testset "RegistryCache" begin
    @testset "get_registry" begin
        mktempdir(@__DIR__) do temp_cache_dir
            # test when registry is not in the cache and not downloaded
            cache = RegistryTools.RegistryCache(temp_cache_dir)
            with(get_registry(DEFAULT_REGISTRY_URL, cache=cache, gitconfig=TEST_GITCONFIG)) do repo
                @test LibGit2.path(repo) == replace(RegistryTools.path(cache, DEFAULT_REGISTRY_URL), '\\'=>'/')
                @test LibGit2.branch(repo) == "master"
                @test !LibGit2.isdirty(repo)
                @test LibGit2.url(LibGit2.lookup_remote(repo, "origin")) == DEFAULT_REGISTRY_URL
            end

            # test when registry is in the cache but not downloaded
            registry_path = RegistryTools.path(cache, DEFAULT_REGISTRY_URL)
            rm(registry_path, recursive=true, force=true)
            @test !ispath(registry_path)
            with(get_registry(DEFAULT_REGISTRY_URL, cache=cache, gitconfig=TEST_GITCONFIG)) do repo
                @test LibGit2.path(repo) == replace(RegistryTools.path(cache, DEFAULT_REGISTRY_URL), '\\'=>'/')
                @test LibGit2.branch(repo) == "master"
                @test !LibGit2.isdirty(repo)
                @test LibGit2.url(LibGit2.lookup_remote(repo, "origin")) == DEFAULT_REGISTRY_URL

                # test when registry is in the cache, downloaded, but mutated
                orig_hash = LibGit2.GitHash()
                LibGit2.branch!(repo, "newbranch", force=true)
                LibGit2.remove!(repo, "Registry.toml")
                LibGit2.commit(
                    repo,
                    "Removing Registry.toml in Registrator tests";
                    author=TEST_SIGNATURE,
                    committer=TEST_SIGNATURE,
                )
                @test LibGit2.GitObject(repo, "HEAD") != LibGit2.GitObject(repo, "master")
                @test ispath(registry_path)
            end
            with(get_registry(DEFAULT_REGISTRY_URL, cache=cache, gitconfig=TEST_GITCONFIG)) do repo
                @test LibGit2.path(repo) == replace(RegistryTools.path(cache, DEFAULT_REGISTRY_URL), '\\'=>'/')
                @test LibGit2.branch(repo) == "master"
                @test !LibGit2.isdirty(repo)
                @test LibGit2.url(LibGit2.lookup_remote(repo, "origin")) == DEFAULT_REGISTRY_URL
            end
        end
    end
end

@testset "RegistryData" begin
    blank = RegistryTools.RegistryData("BlankRegistry", "d4e2f5cd-0f48-4704-9988-f1754e755b45")

    example = Project(Dict(
        "name" => "Example", "uuid" => "7876af07-990d-54b4-ab0e-23690620f79a"
    ))

    @testset "I/O" begin
        registry = """
            name = "General"
            uuid = "23338594-aafe-5451-b93e-139f81909106"
            repo = "https://github.com/JuliaRegistries/General.git"

            description = \"\"\"
            Official general Julia package registry where people can
            register any package they want without too much debate about
            naming and without enforced standards on documentation or
            testing. We nevertheless encourage documentation, testing and
            some amount of consideration when choosing package names.
            \"\"\"

            [packages]
            00701ae9-d1dc-5365-b64a-a3a3ebf5695e = { name = "BioAlignments", path = "B/BioAlignments" }
            00718b61-6157-5045-8849-3d4c4093d022 = { name = "Convertible", path = "C/Convertible" }
            0087ddc6-3964-5e57-817f-9937aefb0357 = { name = "MathOptInterfaceMosek", path = "M/MathOptInterfaceMosek" }
            """

        registry_data = parse_registry(IOBuffer(registry))
        @test registry_data isa RegistryTools.RegistryData
        written_registry = sprint(TOML.print, registry_data)
        written_registry_data = parse_registry(IOBuffer(written_registry))

        @test written_registry_data == registry_data
        @test written_registry == registry
    end

    @testset "Package Operations" begin
        registry_data = copy(blank)

        @test isempty(registry_data.packages)
        @test push!(registry_data, example) == registry_data
        @test length(registry_data.packages) == 1
        @test haskey(registry_data.packages, string(example.uuid))
        @test registry_data.packages[string(example.uuid)]["name"] == "Example"
        @test registry_data.packages[string(example.uuid)]["path"] == "E/Example"
    end
end

@testset "check_version!" begin
    import RegistryTools: ReturnStatus, check_version!, haserror
    import Pkg.Types: Project
    function hascheck(status::ReturnStatus, check)
        return check in (c.id for c in status.triggered_checks)
    end

    for ver in [v"0.0.2", v"0.3.2", v"4.3.2"]
        status = ReturnStatus()
        check_version!(ver, VersionNumber[], status)
        @test hascheck(status, :not_standard_first_version)
        @test hascheck(status, :new_package_label)
        @test length(status.triggered_checks) == 2
        @test !haserror(status)
    end

    for ver in [v"0.0.1", v"0.1.0", v"1.0.0"]
        status = ReturnStatus()
        check_version!(ver, VersionNumber[], status)
        @test hascheck(status, :new_package_label)
        @test length(status.triggered_checks) == 1
        @test !haserror(status)
    end

    versions_list = [v"0.0.5", v"0.1.0", v"0.1.5", v"1.0.0"]
    let    # Less than least existing version
        ver = v"0.0.4"
        status = ReturnStatus()
        check_version!(ver, versions_list, status)
        @test hascheck(status, :version_less_than_all_existing)
        @test length(status.triggered_checks) == 1
        @test !haserror(status)
    end

    let    # Existing version
        ver = v"0.0.5"
        status = ReturnStatus()
        check_version!(ver, versions_list, status)
        @test hascheck(status, :version_exists)
        @test length(status.triggered_checks) == 1
        @test haserror(status)
    end

    # Non-existing version
    for (ver, type) in [(v"0.1.1", "patch"), (v"0.1.6", "patch"), (v"1.0.1", "patch"), (v"1.1.0", "minor")]
        status = ReturnStatus()
        check_version!(ver, versions_list, status)
        @test hascheck(status, Symbol(type, "_release"))
        @test length(status.triggered_checks) == 1
        @test !haserror(status)
    end
    for (ver, type) in [(v"0.0.6", "patch"), (v"0.2.0", "minor"), (v"2.0.0", "major")]
        status = ReturnStatus()
        check_version!(ver, versions_list, status)
        @test hascheck(status, Symbol(type, "_release"))
        @test hascheck(status, :breaking)
        @test length(status.triggered_checks) == 2
        @test !haserror(status)
    end

    # Skip a version
    for (ver, type) in [(v"0.1.2", "patch"), (v"0.1.7", "patch"), (v"1.0.2", "patch"), (v"1.2.0", "minor")]
        status = ReturnStatus()
        check_version!(ver, versions_list, status)
        @test hascheck(status, Symbol(type, "_release"))
        @test hascheck(status, :version_skip)
        @test length(status.triggered_checks) == 2
        @test !haserror(status)
    end
    for (ver, type) in [(v"0.0.7", "patch"), (v"0.3.0", "minor"), (v"3.0.0", "major")]
        status = ReturnStatus()
        check_version!(ver, versions_list, status)
        @test hascheck(status, Symbol(type, "_release"))
        @test hascheck(status, :breaking)
        @test hascheck(status, :version_skip)
        @test length(status.triggered_checks) == 3
        @test !haserror(status)
    end
end

@testset "package_file" begin
    import RegistryTools: check_package!, update_package_file
    mktempdir(@__DIR__) do package_path
        uuid = "698ec630-83b2-4a6d-81d4-a10176273030"
        pkg = Project(Dict("name" => "Example", "uuid" => uuid))
        repo = "https://example.com/example.git"
        status = ReturnStatus()
        @test check_package!("", package_path, status) == ""
        @test haserror(status)
        status = ReturnStatus()
        @test check_package!(repo, package_path, status) == repo
        @test !haserror(status)
        update_package_file(pkg, repo, package_path)
        package_file = joinpath(package_path, "Package.toml")
        @test isfile(package_file)
        @test read(package_file, String) == """
                                            name = "Example"
                                            uuid = "$uuid"
                                            repo = "$repo"
                                            """
        status = ReturnStatus()
        @test check_package!("", package_path, status) == repo
        @test !haserror(status)
        status = ReturnStatus()
        @test check_package!(repo, package_path, status) == repo
        @test !haserror(status)
    end
end

function create_empty_registry(registry_path, registry_name, registry_uuid)
    mkpath(registry_path)
    registry_file = joinpath(registry_path, "Registry.toml")
    registry_data = RegistryTools.RegistryData(registry_name, registry_uuid)
    RegistryTools.write_registry(registry_file, registry_data)
end

@testset "find_package_in_registry" begin
    import RegistryTools: find_package_in_registry
    mktempdir(@__DIR__) do temp_dir
        registry_path = temp_dir
        registry_uuid = "d4e2f5cd-0f48-4704-9988-f1754e755b45"
        create_empty_registry(registry_path, "TestRegistry", registry_uuid)
        registry_file = joinpath(registry_path, "Registry.toml")
        registry_data = parse_registry(registry_file)
        package_uuid = "698ec630-83b2-4a6d-81d4-a10176273030"
        pkg = Project(Dict("name" => "Example",
                           "uuid" => package_uuid))
        status = ReturnStatus()

        find_package_in_registry(pkg, registry_file, registry_path,
                                 registry_data, status)
        @test read(registry_file, String) == """
                                             name = "TestRegistry"
                                             uuid = "$(registry_uuid)"

                                             [packages]
                                             $(package_uuid) = { name = "Example", path = "E/Example" }
                                             """
        @test isdir(joinpath(registry_path, "E", "Example"))
    end
end

@testset "versions_file" begin
    import RegistryTools: get_versions_file, update_versions_file, check_versions!
    mktempdir(@__DIR__) do temp_dir
        pkg = Project(version = v"1.0.0")
        tree_hash = repeat("0", 40)
        status = ReturnStatus()

        filename, data = get_versions_file(temp_dir)
        @test filename == joinpath(temp_dir, "Versions.toml")
        @test data isa Dict
        @test isempty(data)
        check_versions!(pkg, data, status)
        # No previous version registered, this version is fine.
        @test !haserror(status)

        update_versions_file(pkg, filename, data, tree_hash)

        _, data = get_versions_file(temp_dir)
        @test data isa Dict
        @test collect(keys(data)) == ["1.0.0"]
        @test data["1.0.0"] isa Dict
        @test collect(keys(data["1.0.0"])) == ["git-tree-sha1"]
        @test data["1.0.0"]["git-tree-sha1"] == tree_hash

        check_versions!(pkg, data, status)
        # This version was just registered, should be a complaint now.
        @test haserror(status)
    end
end

@testset "deps_file" begin
    import RegistryTools: update_deps_file
    mktempdir(@__DIR__) do temp_dir
        uuid = Base.UUID("8dfed614-e22c-5e08-85e1-65c5234f0b40")
        deps = Dict("Test" => uuid)
        pkg = Project(version = v"1.0.0", deps = deps)
        update_versions_file(pkg, joinpath(temp_dir, "Versions.toml"),
                             Dict{String, Any}(), repeat("0", 40))
        update_deps_file(pkg, temp_dir, VersionNumber[])
        deps_file = joinpath(temp_dir, "Deps.toml")
        @test isfile(deps_file)
        @test read(deps_file, String) == """
                                         [1]
                                         Test = "$uuid"
                                         """
    end
end

@testset "compat_file" begin
    import RegistryTools: update_compat_file
    mktempdir(@__DIR__) do temp_dir
        compat = Dict("julia" => "1.3")
        pkg = Project(version = v"1.0.0", compat = compat)
        update_versions_file(pkg, joinpath(temp_dir, "Versions.toml"),
                             Dict{String, Any}(), repeat("0", 40))
        update_compat_file(pkg, temp_dir, VersionNumber[])
        compat_file = joinpath(temp_dir, "Compat.toml")
        @test isfile(compat_file)
        @test read(compat_file, String) == """
                                           [1]
                                           julia = "1.3.0-1"
                                           """
    end
end

@testset "registry updates" begin
    import RegistryTools: RegistryData, ReturnStatus, haserror,
                          write_registry, check_and_update_registry_files,
                          RegBranch, set_metadata!
    import Pkg.Types: read_project
    registry_update_tests =
        [
            # Add one simple package.
            (project_files = ["Example1"],
             status = Symbol[:new_package, :new_package_label],
             regbranch = (error = false, warning = false,
                          kind = "New package", labels = String["new package"]))

            # Increase patch revision.
            (project_files = ["Example1", "Example2"],
             status = Symbol[:new_version, :patch_release],
             regbranch = (error = false, warning = false,
                          kind = "New version",
                          labels = String["patch_release"]))

            # Increase minor revision.
            (project_files = ["Example1", "Example2", "Example3"],
             status = Symbol[:new_version, :minor_release],
             regbranch = (error = false, warning = false,
                          kind = "New version",
                          labels = String["minor_release"]))

            # Increase major revision.
            (project_files = ["Example1", "Example4"],
             status = Symbol[:new_version, :major_release, :breaking],
             regbranch = (error = false, warning = false,
                          kind = "New version",
                          labels = String["major_release", "BREAKING"]))

            # Adding a minor revision after next major revision has
            # been registered.
            (project_files = ["Example1", "Example4", "Example3"],
             status = Symbol[:new_version, :minor_release],
             regbranch = (error = false, warning = false,
                          kind = "New version",
                          labels = String["minor_release"]))

            # Adding a revision that comes before all registered versions.
            (project_files = ["Example2", "Example4", "Example3", "Example1"],
             status = Symbol[:new_version, :version_less_than_all_existing],
             regbranch = (error = true, warning = false,
                          kind = "New version", labels = String[]))

            # Non-standard first version.
            (project_files = ["Example2"],
             status = Symbol[:new_package, :new_package_label,
                             :not_standard_first_version],
             regbranch = (error = false, warning = true,
                          kind = "New package", labels = String["new package"]))

            # Version zero.
            (project_files = ["Example5"],
             status = Symbol[:new_package, :new_package_label, :version_zero,
                             :not_standard_first_version],
             regbranch = (error = true, warning = true,
                          kind = "New package", labels = String["new package"]))

            # Skipped versions.
            (project_files = ["Example1", "Example6"],
             status = Symbol[:new_version, :minor_release, :version_skip],
             regbranch = (error = false, warning = true,
                          kind = "New version",
                          labels = String["minor_release"]))

            # Adding an existing version.
            (project_files = ["Example1", "Example1"],
             status = Symbol[:new_version, :version_exists],
             regbranch = (error = true, warning = false,
                          kind = "New version",
                          labels = String[]))

            # Changing name.
            (project_files = ["Example1", "Example7"],
             status = Symbol[:change_package_name],
             regbranch = (error = true, warning = false,
                          kind = "",
                          labels = String[]))

            # Changing uuid.
            (project_files = ["Example1", "Example8"],
             status = Symbol[:change_package_uuid],
             regbranch = (error = true, warning = false,
                          kind = "",
                          labels = String[]))

            # Incorrect stdlib uuid.
            (project_files = ["Example9"],
             status = Symbol[:new_package, :new_package_label,
                             :wrong_stdlib_uuid],
             regbranch = (error = true, warning = false,
                          kind = "New package",
                          labels = String["new package"]))

            # Dependency on itself.
            (project_files = ["Example10"],
             status = Symbol[:new_package, :new_package_label,
                             :package_self_dep],
             regbranch = (error = true, warning = false,
                          kind = "New package",
                          labels = String["new package"]))

            # Incorrect spelling of an stdlib.
            # (Arguably this should complain differently since the
            # UUID is known.)
            (project_files = ["Example11"],
             status = Symbol[:new_package, :new_package_label,
                             :dependency_not_found],
             regbranch = (error = true, warning = false,
                          kind = "New package",
                          labels = String["new package"]))

            # Missing dependency.
            (project_files = ["Example12"],
             status = Symbol[:new_package, :new_package_label,
                             :dependency_not_found],
             regbranch = (error = true, warning = false,
                          kind = "New package",
                          labels = String["new package"]))

            # Dependency no longer missing.
            (project_files = ["Dep1", "Example12"],
             status = Symbol[:new_package, :new_package_label],
             regbranch = (error = false, warning = false,
                          kind = "New package",
                          labels = String["new package"]))

            # Dependency not missing but incorrectly spelled.
            (project_files = ["Dep1", "Example13"],
             status = Symbol[:new_package, :new_package_label, :name_mismatch],
             regbranch = (error = true, warning = false,
                          kind = "New package",
                          labels = String["new package"]))

            # Compatibility with Julia before 0.7.
            (project_files = ["Example14"],
             status = Symbol[:new_package, :new_package_label,
                             :julia_before_07_in_compat],
             regbranch = (error = true, warning = false,
                          kind = "New package",
                          labels = String["new package"]))

            # Compat for non-dependency. On Julia 1.2 and later this
            # gives an error already in `read_project`, so only run it
            # on Julia 1.1.
            (project_files = ["Example15"],
             skip_for_newer_julia = true,
             status = Symbol[:new_package, :new_package_label,
                             :invalid_compat],
             regbranch = (error = true, warning = false,
                          kind = "New package",
                          labels = String["new package"]))

            # Compat entry for package in deps.
            (project_files = ["Dep1", "Example16"],
             status = Symbol[:new_package, :new_package_label],
             regbranch = (error = false, warning = false,
                          kind = "New package",
                          labels = String["new package"]))

            # Compat entry for package in extras.
            (project_files = ["Dep1", "Example17"],
             status = Symbol[:new_package, :new_package_label],
             regbranch = (error = false, warning = false,
                          kind = "New package",
                          labels = String["new package"]))

            # Change package repo.
            (project_files = ["Example1", "Example2"],
             modify_package_repo = "Example2",
             status = Symbol[:new_version, :patch_release, :change_package_url],
             regbranch = (error = true, warning = false,
                          kind = "New version",
                          labels = String["patch_release"]))

            # Register new package with empty package repo string.
            (project_files = ["Example1"],
             no_package_repo = "Example1",
             status = Symbol[:new_package, :package_url_missing],
             regbranch = (error = true, warning = false,
                          kind = "New package",
                          labels = String[]))

            # Register new version with empty package repo string.
            (project_files = ["Example1", "Example2"],
             no_package_repo = "Example2",
             status = Symbol[:new_version, :patch_release],
             regbranch = (error = false, warning = false,
                          kind = "New version",
                          labels = String["patch_release"]))
        ]

    mktempdir(@__DIR__) do temp_dir
        registry_path = joinpath(temp_dir, "registry")
        projects_path = joinpath(@__DIR__, "project_files")
        registry_deps_paths = String[]
        tree_hash = repeat("0", 40)
        for test_data in registry_update_tests
            if haskey(test_data, :skip_for_newer_julia) && VERSION >= v"1.2"
                continue
            end
            # Clean up from previous iteration.
            isdir(registry_path) && rm(registry_path, recursive = true)
            # Start with an empty registry.
            create_empty_registry(registry_path, "TestRegistry",
                                  "d4e2f5cd-0f48-4704-9988-f1754e755b45")
            local status, regbr
            # Register some packages.
            for project in test_data.project_files
                project_file = joinpath(projects_path, "$(project).toml")
                pkg = read_project(project_file)
                # Create status object.
                status = ReturnStatus()
                regbr = RegBranch(pkg, "")
                package_repo = "http://example.com/$(pkg.name).git"
                if get(test_data, :modify_package_repo, "") == project
                    package_repo = "http://example.org/$(pkg.name).git"
                elseif get(test_data, :no_package_repo, "") == project
                    package_repo = ""
                end
                check_and_update_registry_files(pkg, package_repo, tree_hash,
                                                registry_path,
                                                registry_deps_paths, status)
                haserror(status) && break
            end
            # Test the return status of the last package registration.
            @test sort([check.id for check in status.triggered_checks]) == sort(test_data.status)
            set_metadata!(regbr, status)
            @test haskey(regbr.metadata, "error") == test_data.regbranch.error
            @test haskey(regbr.metadata, "warning") == test_data.regbranch.warning
            @test get(regbr.metadata, "kind", "") == test_data.regbranch.kind
            @test sort(get(regbr.metadata, "labels", String[])) == sort(test_data.regbranch.labels)
        end
    end
end

# Helper function for the next testset.
function create_and_populate_registry(registry_path, registry_name,
                                      registry_uuid, package)
    # Create an empty registry.
    create_empty_registry(registry_path, registry_name, registry_uuid)

    # Add a package.
    projects_path = joinpath(@__DIR__, "project_files")
    project_file = joinpath(projects_path, "$(package).toml")
    pkg = read_project(project_file)
    package_repo = "http://example.com/$(pkg.name).git"
    tree_hash = repeat("0", 40)
    registry_deps_paths = String[]
    status = ReturnStatus()
    check_and_update_registry_files(pkg, package_repo, tree_hash,
                                    registry_path,
                                    registry_deps_paths, status)

    # Turn the registry into a git repository.
    git = gitcmd(registry_path, TEST_GITCONFIG)
    run(`$git init`)
    run(`$git add .`)
    run(`$git commit -m .`)

    return status
end

# Test the `register` function in its entirety as well as having a
# dependency in a different registry.
#
# The strategy is to first create two registries and populate them
# with one package each. Then `register` is called to register a new
# version of one of the packages with a dependency to the other
# package. The registration is pushed to a new branch via a file URL.
@testset "register" begin
    mktempdir(@__DIR__) do temp_dir
        registry1_path = joinpath(temp_dir, "Registry1")
        status = create_and_populate_registry(registry1_path, "Registry1",
                                              "7e1d4fce-5fe6-405e-8bac-078d4138e9a2",
                                              "Example1")
        @test !haserror(status)

        registry2_path = joinpath(temp_dir, "Registry2")
        status = create_and_populate_registry(registry2_path, "Registry2",
                                              "a5a8be26-c942-4674-beee-533a0e81ac1d",
                                              "Dep1")
        @test !haserror(status)

        projects_path = joinpath(@__DIR__, "project_files")
        project_file = joinpath(projects_path, "Example18.toml")
        pkg = read_project(project_file)
        package_repo = "http://example.com/$(pkg.name).git"
        tree_hash = repeat("0", 40)
        registry_repo = "file://$(registry1_path)"
        deps_repo = "file://$(registry2_path)"
        regbr = register(package_repo, pkg, tree_hash, registry = registry_repo,
                         registry_deps = [deps_repo], push = true,
                         gitconfig = TEST_GITCONFIG)
        @test !haskey(regbr.metadata, "error") && !haskey(regbr.metadata, "warning")
        git = gitcmd(registry1_path, TEST_GITCONFIG)
        branches = readlines(`$git branch`)
        @test length(branches) == 2
    end
end

@testset "find_registered_version" begin
    mktempdir(@__DIR__) do temp_dir
        registry_path = temp_dir

        # Create an empty registry.
        create_empty_registry(registry_path, "Registry1",
                              "7e1d4fce-5fe6-405e-8bac-078d4138e9a2")

        # Add a package.
        projects_path = joinpath(@__DIR__, "project_files")
        project_file = joinpath(projects_path, "Example1.toml")
        pkg = read_project(project_file)
        @test find_registered_version(pkg, registry_path) == ""

        package_repo = string("http://example.com/$(pkg.name).git")
        tree_hash = "7dd821daaae58ddf9fee53e00aa1aab33794d130"
        registry_deps_paths = String[]
        status = ReturnStatus()
        check_and_update_registry_files(pkg, package_repo, tree_hash,
                                        registry_path,
                                        registry_deps_paths, status)

        @test find_registered_version(pkg, registry_path) == tree_hash
        project_file = joinpath(projects_path, "Example2.toml")
        pkg = read_project(project_file)
        @test find_registered_version(pkg, registry_path) == ""
    end
end

@testset "Sorted (recursive) TOML.print for Versions.toml file" begin
    pkg = Project(version = v"2.0.0")
    # versions_file = "Versions.toml"
    version_data = Dict{String,Any}(
         "1.0.0" => Dict("yanked"=>true,"git-tree-sha1"=>"b04b6c6bfd3a607aa1b85362b4854ef612137f3e"),
         "3.0.0" => Dict("git-tree-sha1"=>"96429a372b5c4ad63fa9cbff6ba4178a85939705","foo"=>"bar")
    )
    tree_hash = "20cd0a2651eaf28c8a76c8d7fea4f1107f20174b"
    mktemp() do path, io
        RegistryTools.update_versions_file(pkg, path, version_data, tree_hash::AbstractString)
        close(io)
        @test read(path, String) ==
        """
        ["1.0.0"]
        git-tree-sha1 = "b04b6c6bfd3a607aa1b85362b4854ef612137f3e"
        yanked = true

        ["2.0.0"]
        git-tree-sha1 = "20cd0a2651eaf28c8a76c8d7fea4f1107f20174b"

        ["3.0.0"]
        git-tree-sha1 = "96429a372b5c4ad63fa9cbff6ba4178a85939705"
        foo = "bar"
        """
    end
end

end


text:
# FASTA Record
# ============

mutable struct Record
    # data and filled range
    data::Vector{UInt8}
    filled::UnitRange{Int}
    # indexes
    identifier::UnitRange{Int}
    description::UnitRange{Int}
    sequence::UnitRange{Int}
end

"""
    FASTA.Record()

Create an unfilled FASTA record.
"""
function Record()
    return Record(UInt8[], 1:0, 1:0, 1:0, 1:0)
end

"""
    FASTA.Record(data::Vector{UInt8})

Create a FASTA record object from `data`.

This function verifies and indexes fields for accessors.

!!! warning
    Note that the ownership of `data` is transferred to a new record object.
    Editing the input data will edit the record, and is not advised after
    construction of the record.
"""
function Record(data::Vector{UInt8})
    record = Record(data, 1:0, 1:0, 1:0, 1:0)
    index!(record)
    return record
end

"""
    FASTA.Record(str::AbstractString)

Create a FASTA record object from `str`.

This function verifies and indexes fields for accessors.
"""
Record(str::AbstractString) = Record(Vector{UInt8}(str))

Base.parse(::Record, str::AbstractString) = Record(str)

"""
    FASTA.Record(identifier, sequence)

Create a FASTA record object from `identifier` and `sequence`.
"""
function Record(identifier::AbstractString, sequence)
    return Record(identifier, nothing, sequence)
end

"""
    FASTA.Record(identifier, description, sequence)

Create a FASTA record object from `identifier`, `description` and `sequence`.
"""
function Record(identifier::AbstractString, description::Union{AbstractString,Nothing}, sequence)
    buf = IOBuffer()
    print(buf, '>', strip(identifier))
    if description != nothing
        print(buf, ' ', description)
    end
    print(buf, '\n')
    print(buf, sequence, '\n')
    return Record(take!(buf))
end

function Base.:(==)(record1::Record, record2::Record)
    if isfilled(record1) == isfilled(record2) == true
        r1 = record1.filled
        r2 = record2.filled
        return length(r1) == length(r2) && memcmp(pointer(record1.data, first(r1)), pointer(record2.data, first(r2)), length(r1)) == 0
    else
        return isfilled(record1) == isfilled(record2) == false
    end
end

function Base.copy(record::Record)
    return Record(
        record.data[record.filled],
        record.filled,
        record.identifier,
        record.description,
        record.sequence)
end

function Base.write(io::IO, record::Record)
    return unsafe_write(io, pointer(record.data, first(record.filled)), length(record.filled))
end

function Base.print(io::IO, record::Record)
    write(io, record)
    return nothing
end

function Base.show(io::IO, record::Record)
    print(io, summary(record), ':')
    if isfilled(record)
        println(io)
        println(io, "   identifier: ", hasidentifier(record) ? identifier(record) : "<missing>")
        println(io, "  description: ", hasdescription(record) ? description(record) : "<missing>")
          print(io, "     sequence: ", hassequence(record) ? truncate(sequence(String, record), 40) : "<missing>")
    else
        print(io, " <not filled>")
    end
end

function truncate(s::String, len::Integer)
    if length(s) > len
        return "$(String(collect(Iterators.take(s, len - 1))))…"
    else
        return s
    end
end

function initialize!(record::Record)
    record.filled = 1:0
    record.identifier = 1:0
    record.description = 1:0
    record.sequence = 1:0
    return record
end

function BioGenerics.isfilled(record::Record)
    return !isempty(record.filled)
end

function memcmp(p1::Ptr, p2::Ptr, n::Integer)
    return ccall(:memcmp, Cint, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t), p1, p2, n)
end

# Accessor functions
# ------------------

"""
    identifier(record::Record)::Union{String, Nothing}

Get the sequence identifier of `record`.

!!! note
    Returns nothing if the record has no identifier.
"""
function identifier(record::Record)::Union{String, Nothing}
    checkfilled(record)
    if !hasidentifier(record)
        return nothing
    end
    return String(record.data[record.identifier])
end

"""
    hasidentifier(record::Record)

Checks whether or not the `record` has an identifier.
"""
function hasidentifier(record)
    return isfilled(record) && !isempty(record.identifier)
end

function BioGenerics.seqname(record::Record)
    return identifier(record)
end

function BioGenerics.hasseqname(record::Record)
    return hasidentifier(record)
end

"""
    description(record::Record)::Union{String, Nothing}

Get the description of `record`.

!!! note
    Returns `nothing` if record has no description.
"""
function description(record::Record)::Union{String, Nothing}
    checkfilled(record)
    if !hasdescription(record)
        return nothing
    end
    return String(record.data[record.description])
end

"""
    hasdescription(record::Record)

Checks whether or not the `record` has a description.
"""
function hasdescription(record)
    return isfilled(record) && !isempty(record.description)
end

"""
    header(record::Record)::Union{String, Nothing}

Returns the stripped header line of `record`, or `nothing` if it was empty.
"""
function header(record::Record)
    hasidentifier(record) || return nothing
    id, desc = record.identifier, record.description
    range = first(id) : max(last(id), last(desc))
    return String(record.data[range])
end

"""
    sequence_iter(T, record::Record)

Yields an iterator of the sequence, with elements of type `T`. `T` is constructed
through `T(Char(x))` for each byte `x`. E.g. `sequence_iter(DNA, record)`.
Mutating the record will corrupt the iterator.
"""
function sequence_iter(::Type{T}, record::Record,
    part::UnitRange{<:Integer}=1:lastindex(record.sequence)) where {T <: BioSymbols.BioSymbol}
    checkfilled(record)
    seqpart = record.sequence[part]
    data = record.data
    return (T(Char(@inbounds (data[i]))) for i in seqpart)
end

"""
    sequence(::Type{S}, record::Record, [part::UnitRange{Int}])::S

Get the sequence of `record`.

`S` can be either a subtype of `BioSequences.BioSequence` or `String`.
If `part` argument is given, it returns the specified part of the sequence.

!!! note
    This method makes a new sequence object every time.
    If you have a sequence already and want to fill it with the sequence
    data contained in a fasta record, you can use `Base.copyto!`.
"""
function sequence(::Type{S}, record::Record, part::UnitRange{Int}=1:lastindex(record.sequence))::S where S <: BioSequences.LongSequence
    checkfilled(record)
    seqpart = record.sequence[part]
    return S(record.data, first(seqpart), last(seqpart))
end

function sequence(::Type{String}, record::Record, part::UnitRange{Int}=1:lastindex(record.sequence))::String
    checkfilled(record)
    return String(record.data[record.sequence[part]])
end

"""
    sequence(record::Record, [part::UnitRange{Int}])

Get the sequence of `record`.

This function infers the sequence type from the data. When it is wrong or
unreliable, use `sequence(::Type{S}, record::Record)`.  If `part` argument is
given, it returns the specified part of the sequence.

!!! note
    This method makes a new sequence object every time.
    If you have a sequence already and want to fill it with the sequence
    data contained in a fasta record, you can use `Base.copyto!`.
"""
function sequence(record::Record, part::UnitRange{Int}=1:lastindex(record.sequence))
    checkfilled(record)
    S = predict_seqtype(record.data, record.sequence)
    return sequence(S, record, part)
end

"""
    hassequence(record::Record)

Checks whether or not a sequence record contains a sequence.
"""
function hassequence(record::Record)
    # zero-length sequence may exist
    return isfilled(record)
end

"Get the length of the fasta record's sequence."
@inline seqlen(record::Record) = last(record.sequence) - first(record.sequence) + 1

function Base.copy!(dest::BioSequences.LongSequence, src::Record)
    resize!(dest, seqlen(src) % UInt)
    copyto!(dest, 1, src, 1, seqlen(src))
end

"""
    Base.copyto!(dest::BioSequences.BioSequence, src::Record)

Copy all of the sequence data from the fasta record `src` to a biological
sequence `dest`. `dest` must have a length greater or equal to the length of
the sequence represented in the fastq record. The first n elements of `dest` are
overwritten, the other elements are left untouched.
"""
function Base.copyto!(dest::BioSequences.LongSequence, src::Record)
    return copyto!(dest, 1, src, 1, seqlen(src))
end

"""
    Base.copyto!(dest::BioSequences.BioSequence, doff, src::Record, soff, N)

Copy an N long block of sequence data from the fasta record `src`, starting at
position `soff`, to the `BioSequence` dest, starting at position `doff`.
"""
function Base.copyto!(dest::BioSequences.LongSequence, doff, src::Record, soff, N)
    checkfilled(src)
    if !hassequence(src)
        missingerror(:sequence)
    end

    # This check is here to prevent boundserror when indexing src.sequence
    iszero(N) && return dest
    return copyto!(dest, doff, src.data, src.sequence[soff], N)
end

function BioGenerics.sequence(record::Record)
    return sequence(record)
end

function BioGenerics.sequence(::Type{S}, record::Record) where S <: BioSequences.LongSequence
    return sequence(S, record)
end

function BioGenerics.hassequence(record::Record)
    return hassequence(record)
end

function checkfilled(record)
    if !isfilled(record)
        throw(ArgumentError("unfilled FASTA record"))
    end
end

# Predict sequence type based on character frequencies in `seq[start:stop]`.
function predict_seqtype(seq::Vector{UInt8}, range)
    # count characters
    a = c = g = t = u = n = alpha = 0
    for i in range
        @inbounds x = seq[i]
        if x == 0x41 || x == 0x61
            a += 1
        elseif x == 0x43 || x == 0x63
            c += 1
        elseif x == 0x47 || x == 0x67
            g += 1
        elseif x == 0x54 || x == 0x74
            t += 1
        elseif x == 0x55 || x == 0x75
            u += 1
        elseif x == 0x4e || x == 0x6e
            n += 1
        end
        if 0x41 ≤ x ≤ 0x5a || 0x61 ≤ x ≤ 0x7a
            alpha += 1
            if alpha ≥ 300 && t + u > 0 && a + c + g + t + u + n == alpha
                # pretty sure that the sequence is either DNA or RNA
                break
            end
        end
    end

    # the threshold (= 0.95) is somewhat arbitrary
    if (a + c + g + t + u + n) / alpha > 0.95
        if t ≥ u
            return BioSequences.LongDNASeq
        else
            return BioSequences.LongRNASeq
        end
    else
        return BioSequences.LongAminoAcidSeq
    end
end

function Base.hash(record::Record, h::UInt)
    return hash(identifier(record), hash(description(record), hash(sequence(record), h)))
end

ext: jl


text:
struct HyperCubic{N, B<:NTuple{N, AbstractBoundary}} <: AbstractLattice
    dims::NTuple{N, Int}
    bcs::B
    translation_vectors::Dict{Int, Vector{Coordinate{N, Int}}}
    function HyperCubic{N, B}(dims::NTuple{N, Int}, bcs::B) where {N, B <: NTuple{N, AbstractBoundary}}
        check_boundaries(bcs)
        new{N, B}(dims, bcs, Dict{Int, Vector{Coordinate{N, Int}}}())
    end
end

==(a::HyperCubic{N, B}, b::HyperCubic{N, B}) where {N, B <: NTuple{N, AbstractBoundary}} = (
    a.dims === b.dims
    && a.bcs === b.bcs
    && a.translation_vectors == b.translation_vectors
)
==(a::HyperCubic, b::HyperCubic) = false

ndims(::HyperCubic{N}) where N = N

################
## Constructors
################

HyperCubic{N}(dims::NTuple{N, Int}, bcs::NTuple{N, AbstractBoundary}) where N = HyperCubic{N, typeof(bcs)}(dims, bcs)
function HyperCubic{N}(dims::NTuple{N, Int}, bc::AbstractBoundary=Periodic()) where N
    bcs = ntuple(x->bc, N)
    HyperCubic{N, typeof(bcs)}(dims, bcs)
end

HyperCubic(dims::NTuple{N, Int}, bcs) where N = HyperCubic{N}(dims, bcs)
HyperCubic(dims::NTuple{N, Int}) where N = HyperCubic{N}(dims)

const Chain = HyperCubic{1}
Chain(dim::Int, bcs) = Chain((dim,), bcs)
Chain(dim::Int) = Chain((dim,))

const Square = HyperCubic{2}
const Cubic = HyperCubic{3}

#############################
## Translation/Basis Vectors
#############################

metric(::HyperCubic{N}) where N = I(N)
ismetricdiag(::HyperCubic) = true

function precompute_translation_vectors!(lattice::HyperCubic{1}, maximum_k::Int)
    d = lattice.translation_vectors
    for i in 1:maximum_k
        d[i] = [Coordinate(i), Coordinate(-i)]
    end
    return d
end

#############
## Neighbors
#############

# in 1D, Helical and Periodic BCs are equivalent
for BC in (:Periodic, :Helical), k in (:1, :K)
    @eval function neighbors(lattice::Chain{NTuple{1,$BC}}, site::Coordinate{1, Int}, ::Val{$k}) where K
        l = lattice.dims[1]
        c = site.coordinates[1]
        return [
            Coordinate(apply_boundary(Periodic(), l, c - $k)),
            Coordinate(apply_boundary(Periodic(), l, c + $k)),
        ]
    end
end

function neighbors(lattice::Chain{NTuple{1,Open}}, site::Coordinate{1, Int}, ::Val{K}) where K
    l = lattice.dims[1]
    sp, sm = site + Coordinate(K), site - Coordinate(K)
    if sp.coordinates[1] > l
        if sm.coordinates[1] < 1
            return []
        else
            return [sm]
        end
    else
        if sm.coordinates[1] < 1
            return [sp]
        else
            return [sm, sp]
        end
    end
end

# function neighbors(lattice::HyperCubic{N, NTuple{N, Periodic}}, site::Coordinate{N, Int}, ::Val{K}) where {N,K}
#     return [apply_boundary_conditions(lattice, n) for n in _neighbors(lattice, site, K)]
# end

# function neighbors(lattice::HyperCubic{N, NTuple{N, Helical}}, site::Int, ::Val{1}) where {N}
#     shifts = cumprod(lattice.dims)
#     linear_size = shifts[end]
#     shifts = [1; shifts[1:end-1]]

#     ns = Vector{Int}(undef, 2*N)

#     for i in 1:N
#         ns[i] = mod(site + shifts[i], 1:linear_size)
#         ns[i + N] = mod(site - shifts[i], 1:linear_size)
#     end
#     return ns
# end
# neighbors(lattice::HyperCubic{N, NTuple{N, Helical}}, site::Coordinate{N, Int}, k) where N = neighbors(
#     lattice,
#     sum(map(*, [1; cumprod(lattice.dims[1:end-1])], site.coordinates .- 1)),
#     k
# )


text:
using LinearAlgebra

using ModelingToolkit: isdifferenceeq, has_continuous_events, generate_rootfinding_callback, generate_difference_cb, merge_cb

const MAX_INLINE_NLSOLVE_SIZE = 8

function torn_system_with_nlsolve_jacobian_sparsity(state, var_eq_matching, var_sccs, nlsolve_scc_idxs, eqs_idxs, states_idxs)
    fullvars = state.fullvars
    graph = state.structure.graph

    # The sparsity pattern of `nlsolve(f, u, p)` w.r.t `p` is difficult to
    # determine in general. Consider the "simplest" case, a linear system. We
    # have
    #                   A u = p.
    # Clearly, the sparsity of `u` depends on the sparsity of both `p` and `A`
    # in a non-trivial way. However, in the generic case, `u` is dense even when
    # `A` and `p` are sparse. For instance
    #
    # ```julia
    # julia> using Random, SparseArrays
    #
    # julia> A = sprand(MersenneTwister(1234), 100, 100, 0.1);
    #
    # julia> p = sprand(MersenneTwister(12345), 100, 0.05);
    #
    # julia> count(x->abs(x) < 1e-5, A \ Vector(p))
    # 0
    # ```
    #
    # Let 𝑇 be the set of tearing variables and 𝑉 be the set of all *states* in
    # the residual equations. In the following code, we are going to assume the
    # connection between 𝑇 (the `u` in from above) and 𝑉 ∖ 𝑇 (the `p` in from
    # above) has full incidence.
    #
    # Note that as we are reducing algebraic equations numerically, it could be
    # the case that a later partition (a BLT block) contains tearing variables
    # from other partitions.
    #
    # We know that partitions are BLT ordered. Hence, the tearing variables in
    # each partition is unique, and all states in a partition must be
    # either differential variables or algebraic tearing variables that are
    # from previous partitions. Hence, we can build the dependency chain as we
    # traverse the partitions.

    var_rename = ones(Int64, ndsts(graph))
    nlsolve_vars = Int[]
    for i in nlsolve_scc_idxs, c in var_sccs[i]
        append!(nlsolve_vars, c)
        for v in c
            var_rename[v] = 0
        end
    end
    masked_cumsum!(var_rename)

    dig = DiCMOBiGraph{true}(graph, var_eq_matching)

    fused_var_deps = map(1:ndsts(graph)) do v
        BitSet(v′ for v′ in neighborhood(dig, v, Inf; dir=:in) if var_rename[v′] != 0)
    end

    for scc in var_sccs[nlsolve_scc_idxs]
        if length(scc) >= 2
            deps = fused_var_deps[scc[1]]
            for c in 2:length(scc)
                union!(deps, fused_var_deps[c])
                fused_var_deps[c] = deps
            end
        end
    end

    var2idx = Dict{Int,Int}(v => i for (i, v) in enumerate(states_idxs))
    eqs2idx = Dict{Int,Int}(v => i for (i, v) in enumerate(eqs_idxs))
    nlsolve_vars_set = BitSet(nlsolve_vars)

    I = Int[]; J = Int[]
    s = state.structure
    for ieq in 𝑠vertices(graph)
        nieq = get(eqs2idx, ieq, 0)
        nieq == 0 && continue
        for ivar in 𝑠neighbors(graph, ieq)
            isdervar(s, ivar) && continue
            if var_rename[ivar] != 0
                push!(I, nieq)
                push!(J, var2idx[ivar])
            else
                for dvar in fused_var_deps[ivar]
                    isdervar(s, dvar) && continue
                    niv = get(var2idx, dvar, 0)
                    niv == 0 && continue
                    push!(I, nieq)
                    push!(J, niv)
                end
            end
        end
    end
    sparse(I, J, true)
end

function gen_nlsolve!(is_not_prepended_assignment, eqs, vars, u0map::AbstractDict, assignments, (deps, invdeps), var2assignment; checkbounds=true)
    isempty(vars) && throw(ArgumentError("vars may not be empty"))
    length(eqs) == length(vars) || throw(ArgumentError("vars must be of the same length as the number of equations to find the roots of"))
    rhss = map(x->x.rhs, eqs)
    # We use `vars` instead of `graph` to capture parameters, too.
    paramset = ModelingToolkit.vars(r for r in rhss)

    # Compute necessary assignments for the nlsolve expr
    init_assignments = [var2assignment[p] for p in paramset if haskey(var2assignment, p)]
    if isempty(init_assignments)
        needed_assignments_idxs = Int[]
        needed_assignments = similar(assignments, 0)
    else
        tmp = [init_assignments]
        # `deps[init_assignments]` gives the dependency of `init_assignments`
        while true
            next_assignments = reduce(vcat, deps[init_assignments])
            isempty(next_assignments) && break
            init_assignments = next_assignments
            push!(tmp, init_assignments)
        end
        needed_assignments_idxs = reduce(vcat, unique(reverse(tmp)))
        needed_assignments = assignments[needed_assignments_idxs]
    end

    # Compute `params`. They are like enclosed variables
    rhsvars = [ModelingToolkit.vars(r.rhs) for r in needed_assignments]
    vars_set = Set(vars)
    outer_set = BitSet()
    inner_set = BitSet()
    for (i, vs) in enumerate(rhsvars)
        j = needed_assignments_idxs[i]
        if isdisjoint(vars_set, vs)
            push!(outer_set, j)
        else
            push!(inner_set, j)
        end
    end
    init_refine = BitSet()
    for i in inner_set
        union!(init_refine, invdeps[i])
    end
    intersect!(init_refine, outer_set)
    setdiff!(outer_set, init_refine)
    union!(inner_set, init_refine)

    next_refine = BitSet()
    while true
        for i in init_refine
            id = invdeps[i]
            isempty(id) && break
            union!(next_refine, id)
        end
        intersect!(next_refine, outer_set)
        isempty(next_refine) && break
        setdiff!(outer_set, next_refine)
        union!(inner_set, next_refine)

        init_refine, next_refine = next_refine, init_refine
        empty!(next_refine)
    end
    global2local = Dict(j=>i for (i, j) in enumerate(needed_assignments_idxs))
    inner_idxs = [global2local[i] for i in collect(inner_set)]
    outer_idxs = [global2local[i] for i in collect(outer_set)]
    extravars = reduce(union!, rhsvars[inner_idxs], init=Set())
    union!(paramset, extravars)
    setdiff!(paramset, vars)
    setdiff!(paramset, [needed_assignments[i].lhs for i in inner_idxs])
    union!(paramset, [needed_assignments[i].lhs for i in outer_idxs])
    params = collect(paramset)

    # splatting to tighten the type
    u0 = []
    for v in vars
        v in keys(u0map) || (push!(u0, 1e-3); continue)
        u = substitute(v, u0map)
        for i in 1:length(u0map)
            u = substitute(u, u0map)
            u isa Number && (push!(u0, u); break)
        end
        u isa Number || error("$v doesn't have a default.")
    end
    u0 = [u0...]
    # specialize on the scalar case
    isscalar = length(u0) == 1
    u0 = isscalar ? u0[1] : SVector(u0...)

    fname = gensym("fun")
    # f is the function to find roots on
    f = Func(
        [
         DestructuredArgs(vars, inbounds=!checkbounds)
         DestructuredArgs(params, inbounds=!checkbounds)
        ],
        [],
        Let(
            needed_assignments[inner_idxs],
            isscalar ? rhss[1] : MakeArray(rhss, SVector),
            false
           )
    ) |> SymbolicUtils.Code.toexpr

    # solver call contains code to call the root-finding solver on the function f
    solver_call = LiteralExpr(quote
            $numerical_nlsolve(
                               $fname,
                               # initial guess
                               $u0,
                               # "captured variables"
                               ($(params...),)
                              )
        end)

    preassignments = []
    for i in outer_idxs
        ii = needed_assignments_idxs[i]
        is_not_prepended_assignment[ii] || continue
        is_not_prepended_assignment[ii] = false
        push!(preassignments, assignments[ii])
    end

    nlsolve_expr = Assignment[
                    preassignments
                    fname ← @RuntimeGeneratedFunction(f)
                    DestructuredArgs(vars, inbounds=!checkbounds) ← solver_call
                   ]

    nlsolve_expr
end

function build_torn_function(
        sys;
        expression=false,
        jacobian_sparsity=true,
        checkbounds=false,
        max_inlining_size=nothing,
        kw...
    )

    max_inlining_size = something(max_inlining_size, MAX_INLINE_NLSOLVE_SIZE)
    rhss = []
    eqs = equations(sys)
    eqs_idxs = Int[]
    for (i, eq) in enumerate(eqs)
        isdiffeq(eq) || continue
        push!(eqs_idxs, i)
        push!(rhss, eq.rhs)
    end

    state = get_or_construct_tearing_state(sys)
    fullvars = state.fullvars
    var_eq_matching, var_sccs = algebraic_variables_scc(state)
    condensed_graph = MatchedCondensationGraph(
        DiCMOBiGraph{true}(complete(state.structure.graph), complete(var_eq_matching)), var_sccs)
    toporder = topological_sort_by_dfs(condensed_graph)
    var_sccs = var_sccs[toporder]

    states_idxs = collect(diffvars_range(state.structure))
    mass_matrix_diag = ones(length(states_idxs))

    assignments, deps, sol_states = tearing_assignments(sys)
    invdeps = map(_->BitSet(), deps)
    for (i, d) in enumerate(deps)
        for a in d
            push!(invdeps[a], i)
        end
    end
    var2assignment = Dict{Any,Int}(eq.lhs => i for (i, eq) in enumerate(assignments))
    is_not_prepended_assignment = trues(length(assignments))

    torn_expr = Assignment[]

    defs = defaults(sys)
    nlsolve_scc_idxs = Int[]

    needs_extending = false
    @views for (i, scc) in enumerate(var_sccs)
        torn_vars_idxs = Int[var for var in scc if var_eq_matching[var] !== unassigned]
        torn_eqs_idxs = [var_eq_matching[var] for var in torn_vars_idxs]
        isempty(torn_eqs_idxs) && continue
        if length(torn_eqs_idxs) <= max_inlining_size
            nlsolve_expr = gen_nlsolve!(is_not_prepended_assignment, eqs[torn_eqs_idxs], fullvars[torn_vars_idxs], defs, assignments, (deps, invdeps), var2assignment, checkbounds=checkbounds)
            append!(torn_expr, nlsolve_expr)
            push!(nlsolve_scc_idxs, i)
        else
            needs_extending = true
            append!(eqs_idxs, torn_eqs_idxs)
            append!(rhss, map(x->x.rhs, eqs[torn_eqs_idxs]))
            append!(states_idxs, torn_vars_idxs)
            append!(mass_matrix_diag, zeros(length(torn_eqs_idxs)))
        end
    end

    mass_matrix = needs_extending ? Diagonal(mass_matrix_diag) : I

    out = Sym{Any}(gensym("out"))
    funbody = SetArray(
        !checkbounds,
        out,
        rhss
    )

    states = fullvars[states_idxs]
    syms = map(Symbol, states_idxs)

    pre = get_postprocess_fbody(sys)

    expr = SymbolicUtils.Code.toexpr(
        Func(
             [
              out
              DestructuredArgs(states, inbounds=!checkbounds)
              DestructuredArgs(parameters(sys), inbounds=!checkbounds)
              independent_variables(sys)
             ],
             [],
             pre(Let(
                     [torn_expr; assignments[is_not_prepended_assignment]],
                     funbody,
                     false
                ))
            ),
        sol_states
    )
    if expression
        expr, states
    else
        observedfun = let state = state, dict=Dict(), assignments=assignments, deps=(deps, invdeps), sol_states=sol_states, var2assignment=var2assignment
            function generated_observed(obsvar, u, p, t)
                obs = get!(dict, value(obsvar)) do
                    build_observed_function(state, obsvar, var_eq_matching, var_sccs,
                                            assignments, deps, sol_states, var2assignment,
                                            checkbounds=checkbounds,
                                           )
                end
                obs(u, p, t)
            end
        end

        ODEFunction{true}(
                          @RuntimeGeneratedFunction(expr),
                          sparsity = jacobian_sparsity ? torn_system_with_nlsolve_jacobian_sparsity(state, var_eq_matching, var_sccs, nlsolve_scc_idxs, eqs_idxs, states_idxs) : nothing,
                          syms = syms,
                          observed = observedfun,
                          mass_matrix = mass_matrix,
                         ), states
    end
end

"""
    find_solve_sequence(sccs, vars)

given a set of `vars`, find the groups of equations we need to solve for
to obtain the solution to `vars`
"""
function find_solve_sequence(sccs, vars)
    subset = filter(i -> !isdisjoint(sccs[i], vars), 1:length(sccs))
    isempty(subset) && return []
    vars′ = mapreduce(i->sccs[i], union, subset)
    if vars′ == vars
        return subset
    else
        return find_solve_sequence(sccs, vars′)
    end
end

function build_observed_function(
        state, ts, var_eq_matching, var_sccs,
        assignments,
        deps,
        sol_states,
        var2assignment;
        expression=false,
        output_type=Array,
        checkbounds=true,
    )

    is_not_prepended_assignment = trues(length(assignments))
    if (isscalar = !(ts isa AbstractVector))
        ts = [ts]
    end
    ts = Symbolics.scalarize.(value.(ts))

    vars = Set()
    sys = state.sys
    foreach(Base.Fix1(vars!, vars), ts)
    ivs = independent_variables(sys)
    dep_vars = collect(setdiff(vars, ivs))

    fullvars = state.fullvars
    s = state.structure
    graph = s.graph
    diffvars = map(i->fullvars[i], diffvars_range(s))
    algvars = map(i->fullvars[i], algvars_range(s))

    required_algvars = Set(intersect(algvars, vars))
    obs = observed(sys)
    observed_idx = Dict(map(x->x.lhs, obs) .=> 1:length(obs))
    # FIXME: this is a rather rough estimate of dependencies.
    maxidx = 0
    sts = Set(states(sys))
    for (i, s) in enumerate(dep_vars)
        idx = get(observed_idx, s, nothing)
        if idx === nothing
            if !(s in sts)
                throw(ArgumentError("$s is either an observed nor a state variable."))
            end
            continue
        end
        idx > maxidx && (maxidx = idx)
    end
    vs = Set()
    for idx in 1:maxidx
        vars!(vs, obs[idx].rhs)
        union!(required_algvars, intersect(algvars, vs))
        empty!(vs)
    end

    varidxs = findall(x->x in required_algvars, fullvars)
    subset = find_solve_sequence(var_sccs, varidxs)
    if !isempty(subset)
        eqs = equations(sys)

        nested_torn_vars_idxs = []
        for iscc in subset
            torn_vars_idxs = Int[var for var in var_sccs[iscc] if var_eq_matching[var] !== unassigned]
            isempty(torn_vars_idxs) || push!(nested_torn_vars_idxs, torn_vars_idxs)
        end
        torn_eqs = [[eqs[var_eq_matching[i]] for i in idxs] for idxs in nested_torn_vars_idxs]
        torn_vars = [fullvars[idxs] for idxs in nested_torn_vars_idxs]
        u0map = defaults(sys)
        assignments = copy(assignments)
        solves = map(zip(torn_eqs, torn_vars)) do (eqs, vars)
            gen_nlsolve!(is_not_prepended_assignment, eqs, vars,
                         u0map, assignments, deps, var2assignment; checkbounds=checkbounds)
        end
    else
        solves = []
    end

    subs = []
    for sym in vars
        eqidx = get(observed_idx, sym, nothing)
        eqidx === nothing && continue
        push!(subs, sym ← obs[eqidx].rhs)
    end
    pre = get_postprocess_fbody(sys)

    ex = Code.toexpr(Func(
        [
         DestructuredArgs(diffvars, inbounds=!checkbounds)
         DestructuredArgs(parameters(sys), inbounds=!checkbounds)
         independent_variables(sys)
        ],
        [],
        pre(Let(
            [
             collect(Iterators.flatten(solves))
             map(eq -> eq.lhs←eq.rhs, obs[1:maxidx])
             subs
             assignments[is_not_prepended_assignment]
            ],
            isscalar ? ts[1] : MakeArray(ts, output_type),
            false
           ))
       ), sol_states)

    expression ? ex : @RuntimeGeneratedFunction(ex)
end

struct ODAEProblem{iip}
end

ODAEProblem(args...; kw...) = ODAEProblem{true}(args...; kw...)

"""
    ODAEProblem{iip}(sys, u0map, tspan, parammap = DiffEqBase.NullParameters(); kw...)

This constructor acts similar to the one for [`ODEProblem`](@ref) with the following changes:
`ODESystem`s can sometimes be further reduced if `structural_simplify` has
already been applied to them. This is done this constructor.
In these cases, the constructor uses the knowledge of the strongly connected
components calculated during the process of simplification as the basis for
building pre-simplified nonlinear systems in the implicit solving.
In summary: these problems are structurally modified, but could be
more efficient and more stable. Note, the returned object is still of type [`ODEProblem`](@ref).
"""
function ODAEProblem{iip}(
                          sys,
                          u0map,
                          tspan,
                          parammap=DiffEqBase.NullParameters();
                          callback = nothing,
                          kwargs...
                         ) where {iip}
    fun, dvs = build_torn_function(sys; kwargs...)
    ps = parameters(sys)
    defs = defaults(sys)

    u0 = ModelingToolkit.varmap_to_vars(u0map, dvs; defaults=defs)
    p = ModelingToolkit.varmap_to_vars(parammap, ps; defaults=defs)

    has_difference = any(isdifferenceeq, equations(sys))
    if has_continuous_events(sys)
        event_cb = generate_rootfinding_callback(sys; kwargs...)
    else
        event_cb = nothing
    end
    difference_cb = has_difference ? generate_difference_cb(sys; kwargs...) : nothing
    cb = merge_cb(event_cb, difference_cb)
    cb = merge_cb(cb, callback)

    if cb === nothing
        ODEProblem{iip}(fun, u0, tspan, p; kwargs...)
    else
        ODEProblem{iip}(fun, u0, tspan, p; callback=cb, kwargs...)
    end
end


text:
module TestOutliernicer

using StatsBase: sample, mean
using Test
using TSML

function test_artificialdata()
  Random.seed!(123)
  mdates = DateTime(2017,1,1):Dates.Hour(1):DateTime(2017,6,1)
  mvals = rand(1:1000,length(mdates))
  # create some outliers
  soutliers = rand([500:10000;-10000:500],div(length(mdates),10))
  soutndx = sample(1:length(mdates),length(soutliers))
  mvals[soutndx] = soutliers
  df = DataFrame(Date=mdates,Value=mvals)
  outnicer = Outliernicer(Dict(:dateinterval => Dates.Hour(1)))
  fit!(outnicer,df)
  resdf = transform!(outnicer,df)
  @test round(mean(resdf.Value),digits=2) == 510.32
end
@testset "Outliernicer: using artificial data" begin
  test_artificialdata()
end

function test_basicoutlier()
  Random.seed!(123)
  fname = joinpath(dirname(pathof(TSML)),"../data/testdata.csv")
  csvfilter = CSVDateValReader(Dict(:filename=>fname,:dateformat=>"dd/mm/yyyy HH:MM"))
  valgator = DateValgator(Dict(:dateinterval=>Dates.Hour(1)))
  valnner = DateValNNer(Dict(:dateinterval=>Dates.Hour(1)))
  stfier = Statifier(Dict(:processmissing=>true))
  mono = Monotonicer(Dict())
  outliernicer = Outliernicer(Dict(:dateinterval=>Dates.Hour(1)))

  mpipeline1 = @pipeline csvfilter |> valgator |> mono |> valnner |> outliernicer |> stfier
  respipe1 = fit_transform!(mpipeline1)
  @test round(sum(respipe1[1,3:20])) == -213862.0

  mpipeline2 = @pipeline csvfilter |> valgator |> mono |> outliernicer |> stfier
  respipe2 = fit_transform!(mpipeline2)
  @test round(sum(respipe2[1,3:20])) == -219595.0

  mpipeline3 = @pipeline csvfilter |> valgator |> valnner |> mono |> outliernicer |> stfier
  respipe3 = fit_transform!(mpipeline3)
  @test round(sum(respipe3[1,3:20])) == -213862.0
end
@testset "Outliernicer: readcsv |> valgator |> valnner |> mono |> outliernicer |> stfier" begin
  test_basicoutlier()
end

function test_typesoutliernicer()
  regularfile = joinpath(dirname(pathof(TSML)),"../data/typedetection/regular.csv")
  monofile = joinpath(dirname(pathof(TSML)),"../data/typedetection/monotonic.csv")
  dailymonofile = joinpath(dirname(pathof(TSML)),"../data/typedetection/dailymonotonic.csv")
  regularfilecsv = CSVDateValReader(Dict(:filename=>regularfile,:dateformat=>"dd/mm/yyyy HH:MM"))
  monofilecsv = CSVDateValReader(Dict(:filename=>monofile,:dateformat=>"dd/mm/yyyy HH:MM"))
  dailymonofilecsv = CSVDateValReader(Dict(:filename=>dailymonofile,:dateformat=>"dd/mm/yyyy HH:MM"))

  valgator = DateValgator(Dict(:dateinterval=>Dates.Hour(1)))
  valnner = DateValNNer(Dict(:dateinterval=>Dates.Hour(1)))
  stfier = Statifier(Dict(:processmissing=>true))
  mono = Monotonicer(Dict())
  outliernicer = Outliernicer(Dict(:dateinterval=>Dates.Hour(1)))

  regpipeline = @pipeline regularfilecsv |> valgator |> valnner |> mono |> outliernicer |> stfier
  fit!(regpipeline)
  regulardf=transform!(regpipeline)
  @test round(sum(regulardf[1,3:20])) == -61184.0

  monopipeline = @pipeline monofilecsv |> valgator |> valnner |> mono |> outliernicer |> stfier
  fit!(monopipeline)
  monodf=transform!(monopipeline)
  @test round(sum(monodf[1,3:20])) == -890049.0

  dailymonopipeline = @pipeline dailymonofilecsv |> valgator |> valnner |> mono |> outliernicer |> stfier
  fit!(dailymonopipeline)
  dailymonodf=transform!(dailymonopipeline)

  @test round(sum(dailymonodf[1,3:20])) == -294446.0

end
@testset "Outliernicer: monotonic type and outlier detections" begin
  test_typesoutliernicer()
end

end


text:
function interpolate(s)
    isescape = false
    i1 = 1
    i2 = 1
    res = Expr(:string)
    while i2 <= ncodeunits(s)
        c = s[i2]
        if c == '\\'
            isescape = !isescape
            i2 = nextind(s, i2)
        elseif c == '$' && !isescape
            i1 != i2 && push!(res.args, s[i1:prevind(s, i2)])
            i2 = nextind(s, i2)
            data, i1 = Meta.parse(s, i2; greedy = false)
            i2 = i1
            push!(res.args, data)
        else
            isescape = false
            i2 = nextind(s, i2)
        end
    end

    i1 != i2 && push!(res.args, s[i1:prevind(s, i2)])
    return res
end

macro interpolate(s)
    esc(esc(interpolate(s)))
end


text:
module ArrayInterface

using Requires
using LinearAlgebra
using SparseArrays

function ismutable end

"""
    ismutable(x::DataType)

Query whether a type is mutable or not, see
https://github.com/JuliaDiffEq/RecursiveArrayTools.jl/issues/19.
"""
Base.@pure ismutable(x::DataType) = x.mutable
ismutable(x) = ismutable(typeof(x))

ismutable(::Type{<:Array}) = true
ismutable(::Type{<:Number}) = false

# Piracy
function Base.setindex(x::AbstractArray,v,i...)
  _x = copy(x)
  _x[i...] = v
  _x
end

function Base.setindex(x::AbstractVector,v,i::Int)
  n = length(x)
  x .* (i .!== 1:n) .+ v .* (i .== 1:n)
end

function Base.setindex(x::AbstractMatrix,v,i::Int,j::Int)
  n,m = size(x)
  x .* (i .!== 1:n) .* (j  .!== i:m)' .+ v .* (i .== 1:n) .* (j  .== i:m)'
end

"""
    can_setindex(x::DataType)

Query whether a type can use `setindex!`
"""
can_setindex(x) = true
can_setindex(x::AbstractArray) = can_setindex(typeof(x))

"""
    fast_scalar_indexing(x)

Query whether an array type has fast scalar indexing
"""
fast_scalar_indexing(x) = true
fast_scalar_indexing(x::AbstractArray) = fast_scalar_indexing(typeof(x))

"""
    allowed_getindex(x,i...)

A scalar getindex which is always allowed
"""
allowed_getindex(x,i...) = x[i...]

"""
    allowed_setindex!(x,v,i...)

A scalar setindex! which is always allowed
"""
allowed_setindex!(x,v,i...) = setindex!(x,v,i...)

"""
    isstructured(x::DataType)

Query whether a type is a representation of a structured matrix
"""
isstructured(x) = false
isstructured(x::AbstractArray) = isstructured(typeof(x))
isstructured(::Symmetric) = true
isstructured(::Hermitian) = true
isstructured(::UpperTriangular) = true
isstructured(::LowerTriangular) = true
isstructured(::Tridiagonal) = true
isstructured(::SymTridiagonal) = true
isstructured(::Bidiagonal) = true
isstructured(::Diagonal) = true

"""
    has_sparsestruct(x::AbstractArray)

determine whether `findstructralnz` accepts the parameter `x`
"""
has_sparsestruct(x) = false
has_sparsestruct(x::AbstractArray) = has_sparsestruct(typeof(x))
has_sparsestruct(x::Type{<:AbstractArray}) = false
has_sparsestruct(x::Type{<:SparseMatrixCSC}) = true
has_sparsestruct(x::Type{<:Diagonal}) = true
has_sparsestruct(x::Type{<:Bidiagonal}) = true
has_sparsestruct(x::Type{<:Tridiagonal}) = true
has_sparsestruct(x::Type{<:SymTridiagonal}) = true

"""
    findstructralnz(x::AbstractArray)

Return: (I,J) #indexable objects
Find sparsity pattern of special matrices, the same as the first two elements of findnz(::SparseMatrixCSC)
"""
function findstructralnz(x::Diagonal)
  n=size(x,1)
  (1:n,1:n)
end

abstract type MatrixIndex end

struct BidiagonalIndex <: MatrixIndex
  count::Int
  isup::Bool
end

struct TridiagonalIndex <: MatrixIndex
  count::Int#count==nsize+nsize-1+nsize-1
  nsize::Int
  isrow::Bool
end

struct BandedMatrixIndex <: MatrixIndex
  count::Int
  rowsize::Int
  colsize::Int
  bandinds::Array{Int,1}
  bandsizes::Array{Int,1}
  isrow::Bool
end

function _bandsize(bandind,rowsize,colsize)
  -(rowsize-1) <= bandind <= colsize-1 || throw(ErrorException("Invalid Bandind"))
  if (bandind*(colsize-rowsize)>0) & (abs(bandind)<=abs(colsize-rowsize))
    return min(rowsize,colsize)
  elseif bandind*(colsize-rowsize)<=0
    return min(rowsize,colsize)-abs(bandind)
  else
    return min(rowsize,colsize)-abs(bandind)+abs(colsize-rowsize)
  end
end

function BandedMatrixIndex(rowsize,colsize,lowerbandwidth,upperbandwidth,isrow)
  upperbandwidth>-lowerbandwidth || throw(ErrorException("Invalid Bandwidths"))
  bandinds=upperbandwidth:-1:-lowerbandwidth
  bandsizes=[_bandsize(band,rowsize,colsize) for band in bandinds]
  BandedMatrixIndex(sum(bandsizes),rowsize,colsize,bandinds,bandsizes,isrow)
end

struct BlockBandedMatrixIndex <: MatrixIndex
  count::Int
  refinds::Array{Int,1}
  refcoords::Array{Int,1}#storing col or row inds at ref points
  isrow::Bool
end

function BlockBandedMatrixIndex(nrowblock,ncolblock,rowsizes,colsizes,l,u)
  blockrowind=BandedMatrixIndex(nrowblock,ncolblock,l,u,true)
  blockcolind=BandedMatrixIndex(nrowblock,ncolblock,l,u,false)
  sortedinds=sort([(blockrowind[i],blockcolind[i]) for i in 1:length(blockrowind)],by=x->x[1])
  sort!(sortedinds,by=x->x[2],alg=InsertionSort)#stable sort keeps the second index in order
  refinds=Array{Int,1}()
  refrowcoords=Array{Int,1}()
  refcolcoords=Array{Int,1}()
  rowheights=cumsum(pushfirst!(copy(rowsizes),1))
  blockheight=0
  blockrow=1
  blockcol=1
  currenti=1
  lastrowind=sortedinds[1][1]-1
  lastcolind=sortedinds[1][2]
  for ind in sortedinds
    rowind,colind=ind
    if colind==lastcolind
      if rowind>lastrowind
        blockheight+=rowsizes[rowind]
      end
    else
      for j in blockcol:blockcol+colsizes[lastcolind]-1
        push!(refinds,currenti)
        push!(refrowcoords,blockrow)
        push!(refcolcoords,j)
        currenti+=blockheight
      end
      blockcol+=colsizes[lastcolind]
      blockrow=rowheights[rowind]
      blockheight=rowsizes[rowind]
    end
    lastcolind=colind
    lastrowind=rowind
  end
  for j in blockcol:blockcol+colsizes[lastcolind]-1
    push!(refinds,currenti)
    push!(refrowcoords,blockrow)
    push!(refcolcoords,j)
    currenti+=blockheight
  end
  push!(refinds,currenti)#guard
  push!(refrowcoords,-1)
  push!(refcolcoords,-1)
  rowindobj=BlockBandedMatrixIndex(currenti-1,refinds,refrowcoords,true)
  colindobj=BlockBandedMatrixIndex(currenti-1,refinds,refcolcoords,false)
  rowindobj,colindobj
end

struct BandedBlockBandedMatrixIndex <: MatrixIndex
  count::Int
  refinds::Array{Int,1}
  refcoords::Array{Int,1}#storing col or row inds at ref points
  reflocalinds::Array{BandedMatrixIndex,1}
  isrow::Bool
end

function BandedBlockBandedMatrixIndex(nrowblock,ncolblock,rowsizes,colsizes,l,u,lambda,mu)
  blockrowind=BandedMatrixIndex(nrowblock,ncolblock,l,u,true)
  blockcolind=BandedMatrixIndex(nrowblock,ncolblock,l,u,false)
  sortedinds=sort([(blockrowind[i],blockcolind[i]) for i in 1:length(blockrowind)],by=x->x[1])
  sort!(sortedinds,by=x->x[2],alg=InsertionSort)#stable sort keeps the second index in order
  rowheights=cumsum(pushfirst!(copy(rowsizes),1))
  colwidths=cumsum(pushfirst!(copy(colsizes),1))
  currenti=1
  refinds=Array{Int,1}()
  refrowcoords=Array{Int,1}()
  refcolcoords=Array{Int,1}()
  reflocalrowinds=Array{BandedMatrixIndex,1}()
  reflocalcolinds=Array{BandedMatrixIndex,1}()
  for ind in sortedinds
    rowind,colind=ind
    localrowind=BandedMatrixIndex(rowsizes[rowind],colsizes[colind],lambda,mu,true)
    localcolind=BandedMatrixIndex(rowsizes[rowind],colsizes[colind],lambda,mu,false)
    push!(refinds,currenti)
    push!(refrowcoords,rowheights[rowind])
    push!(refcolcoords,colwidths[colind])
    push!(reflocalrowinds,localrowind)
    push!(reflocalcolinds,localcolind)
    currenti+=localrowind.count
  end
  push!(refinds,currenti)
  push!(refrowcoords,-1)
  push!(refcolcoords,-1)
  rowindobj=BandedBlockBandedMatrixIndex(currenti-1,refinds,refrowcoords,reflocalrowinds,true)
  colindobj=BandedBlockBandedMatrixIndex(currenti-1,refinds,refcolcoords,reflocalcolinds,false)
  rowindobj,colindobj
end

Base.firstindex(ind::MatrixIndex)=1
Base.lastindex(ind::MatrixIndex)=ind.count
Base.length(ind::MatrixIndex)=ind.count
function Base.getindex(ind::BidiagonalIndex,i::Int)
  1 <= i <= ind.count || throw(BoundsError(ind, i))
  if ind.isup
    ii=i+1
  else
    ii=i+1+1
  end
  convert(Int,floor(ii/2))
end

function Base.getindex(ind::TridiagonalIndex,i::Int)
  1 <= i <= ind.count || throw(BoundsError(ind, i))
  offsetu= ind.isrow ? 0 : 1
  offsetl= ind.isrow ? 1 : 0
  if 1 <= i <= ind.nsize
    return i
  elseif ind.nsize < i <= ind.nsize+ind.nsize-1
    return i-ind.nsize+offsetu
  else
    return i-(ind.nsize+ind.nsize-1)+offsetl
  end
end

function Base.getindex(ind::BandedMatrixIndex,i::Int)
  1 <= i <= ind.count || throw(BoundsError(ind, i))
  _i=i
  p=1
  while _i-ind.bandsizes[p]>0
    _i-=ind.bandsizes[p]
    p+=1
  end
  bandind=ind.bandinds[p]
  startfromone=!xor(ind.isrow,(bandind>0))
  if startfromone
    return _i
  else
    return _i+abs(bandind)
  end
end

function Base.getindex(ind::BlockBandedMatrixIndex,i::Int)
  1 <= i <= ind.count || throw(BoundsError(ind, i))
  p=1
  while i-ind.refinds[p]>=0
    p+=1
  end
  p-=1
  _i=i-ind.refinds[p]
  if ind.isrow
    return ind.refcoords[p]+_i
  else
    return ind.refcoords[p]
  end
end

function Base.getindex(ind::BandedBlockBandedMatrixIndex,i::Int)
  1 <= i <= ind.count || throw(BoundsError(ind, i))
  p=1
  while i-ind.refinds[p]>=0
    p+=1
  end
  p-=1
  _i=i-ind.refinds[p]+1
  ind.reflocalinds[p][_i]+ind.refcoords[p]-1
end

function findstructralnz(x::Bidiagonal)
  n=size(x,1)
  isup= x.uplo=='U' ? true : false
  rowind=BidiagonalIndex(n+n-1,isup)
  colind=BidiagonalIndex(n+n-1,!isup)
  (rowind,colind)
end

function findstructralnz(x::Union{Tridiagonal,SymTridiagonal})
  n=size(x,1)
  rowind=TridiagonalIndex(n+n-1+n-1,n,true)
  colind=TridiagonalIndex(n+n-1+n-1,n,false)
  (rowind,colind)
end

function findstructralnz(x::SparseMatrixCSC)
  rowind,colind,_=findnz(x)
  (rowind,colind)
end

abstract type ColoringAlgorithm end

"""
    fast_matrix_colors(A)

    Query whether a matrix has a fast algorithm for getting the structural
    colors of the matrix.
"""
fast_matrix_colors(A) = false
fast_matrix_colors(A::AbstractArray) = fast_matrix_colors(typeof(A))
fast_matrix_colors(A::Type{<:Union{Diagonal,Bidiagonal,Tridiagonal,SymTridiagonal}}) = true

"""
    matrix_colors(A::Union{Array,UpperTriangular,LowerTriangular})

    The color vector for dense matrix and triangular matrix is simply
    `[1,2,3,...,size(A,2)]`
"""
function matrix_colors(A::Union{Array,UpperTriangular,LowerTriangular})
    eachindex(1:size(A,2)) # Vector size matches number of rows
end

function _cycle(repetend,len)
    repeat(repetend,div(len,length(repetend))+1)[1:len]
end

function matrix_colors(A::Diagonal)
    fill(1,size(A,2))
end

function matrix_colors(A::Bidiagonal)
    _cycle(1:2,size(A,2))
end

function matrix_colors(A::Union{Tridiagonal,SymTridiagonal})
    _cycle(1:3,size(A,2))
end

function __init__()

  @require StaticArrays="90137ffa-7385-5640-81b9-e52037218182" begin
    ismutable(::Type{<:StaticArrays.StaticArray}) = false
    can_setindex(::Type{<:StaticArrays.StaticArray}) = false
    ismutable(::Type{<:StaticArrays.MArray}) = true
  end

  @require LabelledArrays="2ee39098-c373-598a-b85f-a56591580800" begin
    ismutable(::Type{<:LabelledArrays.LArray{T,N,Syms}}) where {T,N,Syms} = ismutable(T)
    can_setindex(::Type{<:LabelledArrays.SLArray}) = false
  end

  @require Tracker="9f7883ad-71c0-57eb-9f7f-b5c9e6d3789c" begin
    ismutable(::Type{<:Tracker.TrackedArray}) = false
    can_setindex(::Type{<:Tracker.TrackedArray}) = false
    fast_scalar_indexing(::Type{<:Tracker.TrackedArray}) = false
  end

  @require CuArrays="3a865a2d-5b23-5a0f-bc46-62713ec82fae" begin
    fast_scalar_indexing(::Type{<:CuArrays.CuArray}) = false
    @inline allowed_getindex(x::CuArrays.CuArray,i...) = CuArrays._getindex(x,i...)
    @inline allowed_setindex!(x::CuArrays.CuArray,v,i...) = CuArrays._setindex!(x,v,i...)

    function Base.setindex(x::CuArrays.CuArray,v,i::Int)
      _x = copy(x)
      allowed_setindex!(_x,v,i)
      _x
    end
  end

  @require BandedMatrices="aae01518-5342-5314-be14-df237901396f" begin
    function findstructralnz(x::BandedMatrices.BandedMatrix)
      l,u=BandedMatrices.bandwidths(x)
      rowsize,colsize=size(x)
      rowind=BandedMatrixIndex(rowsize,colsize,l,u,true)
      colind=BandedMatrixIndex(rowsize,colsize,l,u,false)
      (rowind,colind)
    end

    has_sparsestruct(::Type{<:BandedMatrices.BandedMatrix}) = true
    is_structured(::Type{<:BandedMatrices.BandedMatrix}) = true
    fast_matrix_colors(::Type{<:BandedMatrices.BandedMatrix}) = true

    function matrix_colors(A::BandedMatrices.BandedMatrix)
        l,u=BandedMatrices.bandwidths(A)
        width=u+l+1
        _cycle(1:width,size(A,2))
    end

  end

  @require BlockBandedMatrices="ffab5731-97b5-5995-9138-79e8c1846df0" begin
    function findstructralnz(x::BlockBandedMatrices.BlockBandedMatrix)
      l,u=BlockBandedMatrices.blockbandwidths(x)
      nrowblock=BlockBandedMatrices.nblocks(x,1)
      ncolblock=BlockBandedMatrices.nblocks(x,2)
      rowsizes=[BlockBandedMatrices.blocksize(x,(i,1))[1] for i in 1:nrowblock]
      colsizes=[BlockBandedMatrices.blocksize(x,(1,i))[2] for i in 1:ncolblock]
      BlockBandedMatrixIndex(nrowblock,ncolblock,rowsizes,colsizes,l,u)
    end

    function findstructralnz(x::BlockBandedMatrices.BandedBlockBandedMatrix)
      l,u=BlockBandedMatrices.blockbandwidths(x)
      lambda,mu=BlockBandedMatrices.subblockbandwidths(x)
      nrowblock=BlockBandedMatrices.nblocks(x,1)
      ncolblock=BlockBandedMatrices.nblocks(x,2)
      rowsizes=[BlockBandedMatrices.blocksize(x,(i,1))[1] for i in 1:nrowblock]
      colsizes=[BlockBandedMatrices.blocksize(x,(1,i))[2] for i in 1:ncolblock]
      BandedBlockBandedMatrixIndex(nrowblock,ncolblock,rowsizes,colsizes,l,u,lambda,mu)
    end

    has_sparsestruct(::Type{<:BlockBandedMatrices.BlockBandedMatrix}) = true
    has_sparsestruct(::Type{<:BlockBandedMatrices.BandedBlockBandedMatrix}) = true
    is_structured(::Type{<:BlockBandedMatrices.BlockBandedMatrix}) = true
    is_structured(::Type{<:BlockBandedMatrices.BandedBlockBandedMatrix}) = true
    fast_matrix_colors(::Type{<:BlockBandedMatrices.BlockBandedMatrix}) = true
    fast_matrix_colors(::Type{<:BlockBandedMatrices.BandedBlockBandedMatrix}) = true

    function matrix_colors(A::BlockBandedMatrices.BlockBandedMatrix)
        l,u=BlockBandedMatrices.blockbandwidths(A)
        blockwidth=l+u+1
        nblock=BlockBandedMatrices.nblocks(A,2)
        cols=[BlockBandedMatrices.blocksize(A,(1,i))[2] for i in 1:nblock]
        blockcolors=_cycle(1:blockwidth,nblock)
        #the reserved number of colors of a block is the maximum length of columns of blocks with the same block color
        ncolors=[maximum(cols[i:blockwidth:nblock]) for i in 1:blockwidth]
        endinds=cumsum(ncolors)
        startinds=[endinds[i]-ncolors[i]+1 for i in 1:blockwidth]
        colors=[(startinds[blockcolors[i]]:endinds[blockcolors[i]])[1:cols[i]] for i in 1:nblock]
        vcat(colors...)
    end

    function matrix_colors(A::BlockBandedMatrices.BandedBlockBandedMatrix)
        l,u=BlockBandedMatrices.blockbandwidths(A)
        lambda,mu=BlockBandedMatrices.subblockbandwidths(A)
        blockwidth=l+u+1
        subblockwidth=lambda+mu+1
        nblock=BlockBandedMatrices.nblocks(A,2)
        cols=[BlockBandedMatrices.blocksize(A,(1,i))[2] for i in 1:nblock]
        blockcolors=_cycle(1:blockwidth,nblock)
        #the reserved number of colors of a block is the min of subblockwidth and the largest length of columns of blocks with the same block color
        ncolors=[min(subblockwidth,maximum(cols[i:blockwidth:nblock])) for i in 1:min(blockwidth,nblock)]
        endinds=cumsum(ncolors)
        startinds=[endinds[i]-ncolors[i]+1 for i in 1:min(blockwidth,nblock)]
        colors=[_cycle(startinds[blockcolors[i]]:endinds[blockcolors[i]],cols[i]) for i in 1:nblock]
        vcat(colors...)
    end

  end
end

end


text:
import MathOptInterface
using Test

const MOI = MathOptInterface
const MOIU = MOI.Utilities
const MPS = MOI.FileFormats.MPS
const MPS_TEST_FILE = "test.mps"

@test sprint(show, MPS.Model()) == "A Mathematical Programming System (MPS) model"

function test_model_equality(model_string, variables, constraints)
    model = MPS.Model()
    MOIU.loadfromstring!(model, model_string)
    MOI.write_to_file(model, MPS_TEST_FILE)
    model_2 = MPS.Model()
    MOI.read_from_file(model_2, MPS_TEST_FILE)
    MOIU.test_models_equal(model, model_2, variables, constraints)
end

@testset "Errors" begin
    failing_models_dir = joinpath(@__DIR__, "failing_models")

    @testset "Quadratic Objective" begin
        model = MPS.Model()
        @test_throws(
            MOI.UnsupportedAttribute,
            MOIU.loadfromstring!(model, """
               variables: x
               minobjective: 1.0*x*x
               """)
        )
    end

    @testset "Non-empty model" begin
        model = MPS.Model()
        @test MOI.is_empty(model)
        MOI.add_variable(model)
        @test !MOI.is_empty(model)
        MOI.empty!(model)
        @test MOI.is_empty(model)
        MOI.add_variable(model)
        @test_throws Exception MOI.read_from_file(
            model, joinpath(failing_models_dir, "bad_name.mps"))
    end

    @testset "$(filename)" for filename in filter(
        f -> endswith(f, ".mps"), readdir(failing_models_dir))
        @test_throws Exception MOI.read_from_file(MPS.Model(),
            joinpath(failing_models_dir, filename))
    end
end

@testset "ROWS - empty ROW name" begin
    model = MPS.Model()
    x = MOI.add_variable(model)
    MOI.add_constraint(model,
        MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
        MOI.LessThan(1.0)
    )
    @test_throws Exception sprint(MPS.write_rows, model)
end

@testset "SOS" begin
    model = MPS.Model()
    x = MOI.add_variables(model, 3)
    names = Dict{MOI.VariableIndex, String}()
    for i in 1:3
        MOI.set(model, MOI.VariableName(), x[i], "x$(i)")
        names[x[i]] = "x$(i)"
    end
    MOI.add_constraint(model,
        MOI.VectorOfVariables(x),
        MOI.SOS1([1.5, 2.5, 3.5])
    )
    MOI.add_constraint(model,
        MOI.VectorOfVariables(x),
        MOI.SOS2([1.25, 2.25, 3.25])
    )
    @test sprint(MPS.write_sos, model, names) ==
        "SOS\n" *
        " S1 SOS1\n" *
        "    x1        1.5\n" *
        "    x2        2.5\n" *
        "    x3        3.5\n" *
        " S2 SOS2\n" *
        "    x1        1.25\n" *
        "    x2        2.25\n" *
        "    x3        3.25\n"
end

@testset "Maximization problems" begin
    model = MPS.Model()
    x = MOI.add_variable(model)
    MOI.set(model, MOI.VariableName(), x, "x")
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.set(model, MOI.ObjectiveFunction{MOI.SingleVariable}(),
        MOI.SingleVariable(x))
    @test sprint(MPS.write_columns, model, ["x"], Dict(x => "x")) ==
        "COLUMNS\n    x         OBJ       -1\n"
end

@testset "stacked_data" begin
    model = MPS.Model()
    MOI.read_from_file(model, joinpath(@__DIR__, "stacked_data.mps"))
    MOI.set(model, MOI.ConstraintName(), MOI.get(model,
            MOI.ListOfConstraintIndices{MOI.SingleVariable, MOI.Integer}())[1],
        "con5")
    MOI.set(model, MOI.ConstraintName(), MOI.get(model,
            MOI.ListOfConstraintIndices{MOI.SingleVariable, MOI.Interval{Float64}}())[1],
        "con6")
    MOI.set(model, MOI.ConstraintName(), MOI.get(model,
            MOI.ListOfConstraintIndices{MOI.SingleVariable, MOI.ZeroOne}())[1],
        "con7")
    model_2 = MPS.Model()
    MOIU.loadfromstring!(model_2, """
    variables: x, y, z
    minobjective: x + y + z
    con1: 1.0 * x in Interval(1.0, 5.0)
    con2: 1.0 * x in Interval(2.0, 6.0)
    con3: 1.0 * x in Interval(3.0, 7.0)
    con4: 2.0 * x in Interval(4.0, 8.0)
    con5: y in Integer()
    con6: y in Interval(1.0, 4.0)
    con7: z in ZeroOne()
    """)
    MOI.set(model_2, MOI.Name(), "stacked_data")
    MOIU.test_models_equal(model, model_2, ["x", "y", "z"],
        ["con1", "con2", "con3", "con4", "con5", "con6", "con7"])
end

@testset "free_integer" begin
    model = MPS.Model()
    MOI.read_from_file(model, joinpath(@__DIR__, "free_integer.mps"))
    MOI.set(model, MOI.ConstraintName(), MOI.get(model,
            MOI.ListOfConstraintIndices{MOI.SingleVariable, MOI.Integer}())[1],
        "con2")
    model_2 = MPS.Model()
    MOIU.loadfromstring!(model_2, """
    variables: x
    minobjective: x
    con1: 1.0 * x >= 1.0
    con2: x in Integer()
    """)
    MOIU.test_models_equal(model, model_2, ["x"], ["con1", "con2"])
end

@testset "Round trips" begin
    @testset "min objective" begin
        test_model_equality("""
            variables: x
            minobjective: x
        """, ["x"], String[])
    end
    @testset "default RHS >=" begin
        test_model_equality("""
            variables: x
            minobjective: x
            c1: 2.0 * x >= 0.0
        """, ["x"], ["c1"])
    end
    @testset "default RHS <=" begin
        test_model_equality("""
            variables: x
            minobjective: x
            c1: 2.0 * x <= 0.0
        """, ["x"], ["c1"])
    end
    @testset "default RHS ==" begin
        test_model_equality("""
            variables: x
            minobjective: x
            c1: 2.0 * x == 0.0
        """, ["x"], ["c1"])
    end
    @testset "min scalaraffine" begin
        test_model_equality("""
            variables: x
            minobjective: 1.2x
        """, ["x"], String[])
    end

    @testset "ScalarAffine-in-GreaterThan" begin
        test_model_equality("""
            variables: x
            minobjective: 1.2x
            c1: 1.1 * x >= 2.0
        """, ["x"], ["c1"])
    end
    @testset "ScalarAffine-in-LessThan" begin
        test_model_equality("""
            variables: x
            minobjective: 1.2x
            c1: 1.1 * x <= 2.0
        """, ["x"], ["c1"])
    end
    @testset "ScalarAffine-in-EqualTo" begin
        test_model_equality("""
            variables: x
            minobjective: 1.2x
            c1: 1.1 * x == 2.0
        """, ["x"], ["c1"])
    end
    @testset "ScalarAffine-in-Interval" begin
        test_model_equality("""
            variables: x
            minobjective: 1.2x
            c1: 1.1 * x in Interval(1.0, 2.0)
        """, ["x"], ["c1"])
    end
    @testset "MARKER INT" begin
        model = MPS.Model()
        MOIU.loadfromstring!(model, """
            variables: x, y, z
            minobjective: x + y + z
            c1: x in Integer()
            c2: 2 * x + -1.0 * z <= 1.0
            c3: z in ZeroOne()
            c4: x >= 1.0
        """)
        MOI.write_to_file(model, MPS_TEST_FILE)
        model_2 = MPS.Model()
        MOI.read_from_file(model_2, MPS_TEST_FILE)
        for (set_type, constraint_name) in [(MOI.Integer, "c1"),
                                            (MOI.ZeroOne, "c3"),
                                            (MOI.GreaterThan{Float64}, "c4")]
            MOI.set(model_2, MOI.ConstraintName(), MOI.get(model_2,
                MOI.ListOfConstraintIndices{MOI.SingleVariable, set_type}())[1],
                constraint_name)
        end
        MOIU.test_models_equal(
            model, model_2, ["x", "y", "z"], ["c1", "c2", "c3", "c4"])
    end
    @testset "Zero variable bounds" begin
        model = MPS.Model()
        MOIU.loadfromstring!(model, """
            variables: x, y, z
            minobjective: x + y + z
            c1: x >= 0.0
            c2: y <= 0.0
        """)
        MOI.write_to_file(model, MPS_TEST_FILE)
        model_2 = MPS.Model()
        MOI.read_from_file(model_2, MPS_TEST_FILE)
        for (set_type, constraint_name) in [(MOI.GreaterThan{Float64}, "c1"),
                                            (MOI.LessThan{Float64}, "c2")]
            MOI.set(model_2, MOI.ConstraintName(), MOI.get(model_2,
                MOI.ListOfConstraintIndices{MOI.SingleVariable, set_type}())[1],
                constraint_name)
        end
        MOIU.test_models_equal(
            model, model_2, ["x", "y", "z"], ["c1", "c2"])
    end
    @testset "Non-zero variable bounds" begin
        model = MPS.Model()
        MOIU.loadfromstring!(model, """
            variables: w, x, y, z
            minobjective: w + x + y + z
            c1: x == 1.0
            c2: y >= 2.0
            c3: z <= 3.0
            c4: w in Interval(4.0, 5.0)
        """)
        MOI.write_to_file(model, MPS_TEST_FILE)
        model_2 = MPS.Model()
        MOI.read_from_file(model_2, MPS_TEST_FILE)
        for (set_type, constraint_name) in [(MOI.EqualTo{Float64}, "c1"),
                                            (MOI.GreaterThan{Float64}, "c2"),
                                            (MOI.LessThan{Float64}, "c3"),
                                            (MOI.Interval{Float64}, "c4")]
            MOI.set(model_2, MOI.ConstraintName(), MOI.get(model_2,
                MOI.ListOfConstraintIndices{MOI.SingleVariable, set_type}())[1],
                constraint_name)
        end
        MOIU.test_models_equal(
            model, model_2, ["w", "x", "y", "z"], ["c1", "c2", "c3", "c4"])
    end
    @testset "Multiple variable bounds" begin
        model = MPS.Model()
        MOIU.loadfromstring!(model, """
            variables: a_really_long_name
            minobjective: a_really_long_name
            c1: a_really_long_name >= 1.0
            c2: a_really_long_name <= 2.0
        """)
        MOI.write_to_file(model, MPS_TEST_FILE)
        @test read(MPS_TEST_FILE, String) ==
            "NAME          \n" *
            "ROWS\n" *
            " N  OBJ\n" *
            "COLUMNS\n" *
            "    a_really_long_name OBJ       1\n" *
            "RHS\n" *
            "RANGES\n" *
            "BOUNDS\n" *
            " LO bounds    a_really_long_name 1\n" *
            " UP bounds    a_really_long_name 2\n" *
            "ENDATA\n"
    end
    @testset "Un-used variable" begin
        # In this test, `x` will not be written to the file since it does not
        # appear in the objective or in the constriants.
        model = MPS.Model()
        MOIU.loadfromstring!(model, """
            variables: x, y
            minobjective: y
            c1: 2.0 * y >= 1.0
            c2: x >= 0.0
        """)
        MOI.write_to_file(model, MPS_TEST_FILE)
        @test MOI.get(model, MOI.NumberOfVariables()) == 2
        model2 = MPS.Model()
        MOI.read_from_file(model2, MPS_TEST_FILE)
        @test MOI.get(model2, MOI.NumberOfVariables()) == 1
    end
    @testset "Names with spaces" begin
        model = MPS.Model()
        x = MOI.add_variable(model)
        MOI.set(model, MOI.VariableName(), x, "x[1, 2]")
        c = MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction(
                [MOI.ScalarAffineTerm(1.0, x)],
                0.0
            ),
            MOI.EqualTo(1.0)
        )
        MOI.set(model, MOI.ConstraintName(), c, "c c")
        @test sprint(write, model) ==
            "NAME          \n" *
            "ROWS\n" *
            " N  OBJ\n" *
            " E  c_c\n" *
            "COLUMNS\n" *
            "    x[1,_2]   c_c       1\n" *
            "RHS\n" *
            "    rhs       c_c       1\n" *
            "RANGES\n" *
            "BOUNDS\n" *
            " FR bounds    x[1,_2]\n" *
            "ENDATA\n"
    end
end

# Clean up
sleep(1.0)  # Allow time for unlink to happen.
rm(MPS_TEST_FILE, force = true)


text:
# Autogenerated wrapper script for gitoxide_jll for aarch64-apple-darwin
export ein, gix

JLLWrappers.@generate_wrapper_header("gitoxide")
JLLWrappers.@declare_executable_product(ein)
JLLWrappers.@declare_executable_product(gix)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        ein,
        "bin/ein",
    )

    JLLWrappers.@init_executable_product(
        gix,
        "bin/gix",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
export Lens, set, get, modify
export @lens
export set, get, modify
export setproperties

import Base: get
using Base: setindex, getproperty

"""
    Lens

A `Lens` allows to access or replace deeply nested parts of complicated objects.

# Example
```jldoctest
julia> using Setfield

julia> struct T;a;b; end

julia> obj = T("AA", "BB")
T("AA", "BB")

julia> lens = @lens _.a
(@lens _.a)

julia> get(obj, lens)
"AA"

julia> set(obj, lens, 2)
T(2, "BB")

julia> obj
T("AA", "BB")

julia> modify(lowercase, obj, lens)
T("aa", "BB")
```

# Interface
Concrete subtypes of `Lens` have to implement
* `set(obj, lens, val)`
* `get(obj, lens)`

These must be pure functions, that satisfy the three lens laws:

```jldoctest; output = false, setup = :(using Setfield; obj = (a="A", b="B"); lens = @lens _.a; val = 2; val1 = 10; val2 = 20)
@assert get(set(obj, lens, val), lens) == val
        # You get what you set.
@assert set(obj, lens, get(obj, lens)) == obj
        # Setting what was already there changes nothing.
@assert set(set(obj, lens, val1), lens, val2) == set(obj, lens, val2)
        # The last set wins.

# output

```

See also [`@lens`](@ref), [`set`](@ref), [`get`](@ref), [`modify`](@ref).
"""
abstract type Lens end

"""
    modify(f, obj, l::Lens)

Replace a deeply nested part `x` of `obj` by `f(x)`. See also [`Lens`](@ref).
"""
function modify end

"""
    get(obj, l::Lens)

Access a deeply nested part of `obj`. See also [`Lens`](@ref).
"""
function get end

"""
    set(obj, l::Lens, val)

Replace a deeply nested part of `obj` by `val`. See also [`Lens`](@ref).
"""
function set end

@inline function modify(f, obj, l::Lens)
    old_val = get(obj, l)
    new_val = f(old_val)
    set(obj, l, new_val)
end

struct IdentityLens <: Lens end
get(obj, ::IdentityLens) = obj
set(obj, ::IdentityLens, val) = val

struct PropertyLens{fieldname} <: Lens end

function get(obj, l::PropertyLens{field}) where {field}
    getproperty(obj, field)
end

@generated function set(obj, l::PropertyLens{field}, val) where {field}
    Expr(:block,
         Expr(:meta, :inline),
        :(setproperties(obj, ($field=val,)))
       )
end

@generated constructor_of(::Type{T}) where T =
    getfield(parentmodule(T), nameof(T))

function assert_hasfields(T, fnames)
    for fname in fnames
        if !(fname in fieldnames(T))
            msg = "$T has no field $fname"
            throw(ArgumentError(msg))
        end
    end
end

"""
    setproperties(obj, patch)

Return a copy of `obj` with attributes updates accoring to `patch`.

# Examples
```jldoctest
julia> using Setfield

julia> struct S;a;b;c; end

julia> s = S(1,2,3)
S(1, 2, 3)

julia> setproperties(s, (a=10,c=4))
S(10, 2, 4)

julia> setproperties((a=1,c=2,b=3), (a=10,c=4))
(a = 10, c = 4, b = 3)
```
"""
function setproperties end

@generated function setproperties(obj, patch)
    assert_hasfields(obj, fieldnames(patch))
    args = map(fieldnames(obj)) do fn
        if fn in fieldnames(patch)
            :(patch.$fn)
        else
            :(obj.$fn)
        end
    end
    Expr(:block,
        Expr(:meta, :inline),
        Expr(:call,:(constructor_of($obj)), args...)
    )
end

@generated function setproperties(obj::NamedTuple, patch)
    # this function is only generated to force the following check
    # at compile time
    assert_hasfields(obj, fieldnames(patch))
    Expr(:block,
        Expr(:meta, :inline),
        :(merge(obj, patch))
    )
end

struct ComposedLens{LO, LI} <: Lens
    outer::LO
    inner::LI
end

compose() = IdentityLens()
compose(l::Lens) = l
compose(::IdentityLens, ::IdentityLens) = IdentityLens()
compose(::IdentityLens, l::Lens) = l
compose(l::Lens, ::IdentityLens) = l
compose(outer::Lens, inner::Lens) = ComposedLens(outer, inner)
function compose(l1::Lens, ls::Lens...)
    # We can build _.a.b.c as (_.a.b).c or _.a.(b.c)
    # The compiler prefers (_.a.b).c
    compose(l1, compose(ls...))
end

"""
    lens₁ ∘ lens₂
    compose([lens₁, [lens₂, [lens₃, ...]]])

Compose lenses `lens₁`, `lens₂`, ..., `lensₙ` to access nested objects.

# Example
```jldoctest
julia> using Setfield

julia> obj = (a = (b = (c = 1,),),);

julia> la = @lens _.a
       lb = @lens _.b
       lc = @lens _.c
       lens = la ∘ lb ∘ lc
(@lens _.a.b.c)

julia> get(obj, lens)
1
```
"""
Base.:∘(l1::Lens, l2::Lens) = compose(l1, l2)

function get(obj, l::ComposedLens)
    inner_obj = get(obj, l.outer)
    get(inner_obj, l.inner)
end

function set(obj,l::ComposedLens, val)
    inner_obj = get(obj, l.outer)
    inner_val = set(inner_obj, l.inner, val)
    set(obj, l.outer, inner_val)
end

struct IndexLens{I <: Tuple} <: Lens
    indices::I
end

Base.@propagate_inbounds function get(obj, l::IndexLens)
    getindex(obj, l.indices...)
end
Base.@propagate_inbounds function set(obj, l::IndexLens, val)
    setindex(obj, val, l.indices...)
end

"""
    ConstIndexLens{I}

Lens with index stored in type parameter.  This is useful for type-stable
[`get`](@ref) and [`set`](@ref) operations on tuples and named tuples.

This lens can be constructed by, e.g., `@lens _[\$1]`.  Complex expression
must be wrapped with `\$(...)` like `@lens _[\$(length(xs))]`.

# Examples
```jldoctest
julia> using Setfield

julia> get((1, 2.0), @lens _[\$1])
1

julia> Base.promote_op(get, typeof.(((1, 2.0), @lens _[\$1]))...)
Int64

julia> Base.promote_op(get, typeof.(((1, 2.0), @lens _[1]))...) !== Int
true
```
"""
struct ConstIndexLens{I} <: Lens end

Base.@propagate_inbounds get(obj, ::ConstIndexLens{I}) where I = obj[I...]

Base.@propagate_inbounds set(obj, ::ConstIndexLens{I}, val) where I =
    setindex(obj, val, I...)

@generated function set(obj::Union{Tuple, NamedTuple},
                        ::ConstIndexLens{I},
                        val) where I
    if length(I) == 1
        n, = I
        args = map(1:length(obj.types)) do i
            i == n ? :val : :(obj[$i])
        end
        quote
            $(Expr(:meta, :inline))
            ($(args...),)
        end
    else
        quote
            throw(ArgumentError($(string(
                "A `Tuple` and `NamedTuple` can only be indexed with one ",
                "integer.  Given: $I"))))
        end
    end
end

Base.@deprecate get(lens::Lens, obj)       get(obj, lens)
Base.@deprecate set(lens::Lens, obj, val)  set(obj, lens, val)
Base.@deprecate modify(f, lens::Lens, obj) modify(f, obj, lens)