text:
### A Pluto.jl notebook ###
# v0.12.16

using Markdown
using InteractiveUtils

# ╔═╡ 3e0c149a-3350-11eb-0b17-7bea5d1a28d4
begin
	DAY = 4
	md"# 2020 Day $DAY"
end

# ╔═╡ 18080fe2-339b-11eb-34f0-f393638506fc
md"## Part 1"

# ╔═╡ 4e084a4e-339b-11eb-3ffc-41118215f226
md"*Allow passports without a 'cir' key to pass*"

# ╔═╡ 24e1613c-339b-11eb-2976-59884053d1ee
input = read("../data/day_$DAY.txt", String)

# ╔═╡ 3a0fa07a-3605-11eb-30ef-4598213ce6ec
test_input = """
ecl:gry pid:860033327 eyr:2020 hcl:#fffffd
byr:1937 iyr:2017 cid:147 hgt:183cm

iyr:2013 ecl:amb cid:350 eyr:2023 pid:028048884
hcl:#cfa07d byr:1929

hcl:#ae17e1 iyr:2013
eyr:2024
ecl:brn pid:760753108 byr:1931
hgt:179cm

hcl:#cfa07d eyr:2025 pid:166559648
iyr:2011 ecl:brn hgt:59in
"""

# ╔═╡ 08cd4892-362c-11eb-28d8-81540aa5945e
line2keys(line) = Set([split(entry,":")[1] for entry in split(line)])

# ╔═╡ 34846d78-362f-11eb-3547-4978e97e89d9
fields_raw = """
    byr (Birth Year)
    iyr (Issue Year)
    eyr (Expiration Year)
    hgt (Height)
    hcl (Hair Color)
    ecl (Eye Color)
    pid (Passport ID)
    cid (Country ID)
"""

# ╔═╡ 46ee77ec-362f-11eb-3c13-f384ed78a348
required_fields = Set([string(split(line)[1]) for line in split(strip(fields_raw),"\n")])

# ╔═╡ f08c3ad2-362f-11eb-295d-4ddac22b387f
delete!(required_fields,"cid")

# ╔═╡ 673d8218-366b-11eb-2a4d-f548f9d131ee
has_required_fields(line) = required_fields ⊆ line2keys(line)

# ╔═╡ 8c68eb76-362e-11eb-300f-6b92a90010fa
function part1(input::String)
	sum([has_required_fields(line) for line in split(input,"\n\n")])
end

# ╔═╡ 4b3eb86a-3630-11eb-39c0-6d67c0f84a22
@assert part1(test_input) == 2

# ╔═╡ 53200a48-3630-11eb-05e9-99811ae65d11
part1(input)

# ╔═╡ 5f797e22-339d-11eb-1ff3-dda90f946838
md"## Part 2

Validate the values of fields too.
"

# ╔═╡ 6a522dd4-3641-11eb-1c9e-45074da49c82
line2pair(line) = [split(entry,":") for entry in split(line)]

# ╔═╡ 9ac14414-3641-11eb-33c4-75b9e4357068
function is_valid(key,value)
	if key == "byr"
		return 1920 ≤ parse(Int,value) ≤ 2002
	elseif key == "iyr"
		return 2010 ≤ parse(Int,value) ≤ 2020
	elseif key == "eyr"
		return 2020 ≤ parse(Int,value) ≤ 2030
	elseif key == "hgt"
		if last(value,2)=="cm"
			return return 150 ≤ parse(Int, value[1:end-2]) ≤ 193
		elseif last(value,2)=="in"
			return return 59 ≤ parse(Int, value[1:end-2]) ≤ 76
		else
			return false
		end
	elseif key == "hcl"
		return (
			(value[1] == '#') &
			( Set(value[2:end]) ⊆ Set("0123456789abcdef")) &
			( length(value) == 7)
			)
	elseif key == "ecl"
		return value ∈ split("amb blu brn gry grn hzl oth")
	elseif key == "pid"
		# 9 digit number
		return length(value) == 9
	else
		return true
	end
end

# ╔═╡ 5189e182-3664-11eb-0dcc-f508fe33be4c
@assert is_valid("byr",   "2002") ==true

# ╔═╡ 6947baf8-3664-11eb-0ace-3f6c7f808068
@assert is_valid("byr", "2003") == false

# ╔═╡ 694808aa-3664-11eb-0443-c162e7388740
@assert is_valid("hgt",   "60in") ==true

# ╔═╡ 694e56b0-3664-11eb-0d26-b9e7b641b9f6
@assert is_valid("hgt",   "190cm") ==true

# ╔═╡ 694ee300-3664-11eb-11f7-49462e6aa1d2
@assert is_valid("hgt", "190in") == false

# ╔═╡ 69599412-3664-11eb-26c7-21afa8e93b9c
@assert is_valid("hgt", "190") == false

# ╔═╡ 695a2420-3664-11eb-0978-5be1861a5f49
@assert is_valid("hcl",   "#123abc") ==true

# ╔═╡ 69657fa2-3664-11eb-3056-a98b35ba51e3
@assert is_valid("hcl", "#123abz") == false

# ╔═╡ 696628c6-3664-11eb-2643-fd542cd9eab5
@assert is_valid("hcl", "123abc") == false

# ╔═╡ 6971b038-3664-11eb-2c31-c5245b309553
@assert is_valid("ecl",   "brn") ==true

# ╔═╡ 6972e2b4-3664-11eb-1052-af8fccaeac05
@assert is_valid("ecl", "wat") == false

# ╔═╡ 697f53fa-3664-11eb-3962-030d46f24ba7
@assert is_valid("pid",   "000000001") ==true

# ╔═╡ 69800c82-3664-11eb-2d6e-155e57cdc323
@assert is_valid("pid", "0123456789") == false

# ╔═╡ edc481a6-3664-11eb-1c02-797a3670af3f
[line2pair(line) for line in split(input,"\n\n")]

# ╔═╡ 3da5c12e-3669-11eb-1a79-73e9fad032e8
all_valid(pair) = all([is_valid(t[1], t[2]) for t in pair])

# ╔═╡ 05f62636-366b-11eb-1435-632abff666d1
is_valid_passport(line) = all_valid(line2pair(line)) & has_required_fields(line)

# ╔═╡ 42b157f0-3669-11eb-2abc-f16c5da1a521
function part2(input)
	sum(is_valid_passport.(split(input,"\n\n")))
end

# ╔═╡ dfa0bb9c-3668-11eb-26ae-59f47e74033f
invalid_passports = """
eyr:1972 cid:100
hcl:#18171d ecl:amb hgt:170 pid:186cm iyr:2018 byr:1926

iyr:2019
hcl:#602927 eyr:1967 hgt:170cm
ecl:grn pid:012533040 byr:1946

hcl:dab227 iyr:2012
ecl:brn hgt:182cm pid:021572410 eyr:2020 byr:1992 cid:277

hgt:59cm ecl:zzz
eyr:2038 hcl:74454a iyr:2023
pid:3556412378 byr:2007
"""

# ╔═╡ fd483df0-3668-11eb-09f4-bb26699ab6bf
valid_passports = """
pid:087499704 hgt:74in ecl:grn iyr:2012 eyr:2030 byr:1980
hcl:#623a2f

eyr:2029 ecl:blu cid:129 byr:1989
iyr:2014 pid:896056539 hcl:#a97842 hgt:165cm

hcl:#888785
hgt:164cm byr:2001 iyr:2015 cid:88
pid:545766238 ecl:hzl
eyr:2022
"""

# ╔═╡ d95d9cbc-3665-11eb-131c-3f53f873ce71
@assert part2(invalid_passports) == 0

# ╔═╡ 5b8c1242-3669-11eb-2b51-2d878f04f456
@assert part2(valid_passports) == 3

# ╔═╡ dfff04d0-3669-11eb-37a8-67ede57fa4f6
part2(input)

# ╔═╡ Cell order:
# ╠═3e0c149a-3350-11eb-0b17-7bea5d1a28d4
# ╟─18080fe2-339b-11eb-34f0-f393638506fc
# ╠═4e084a4e-339b-11eb-3ffc-41118215f226
# ╠═24e1613c-339b-11eb-2976-59884053d1ee
# ╠═3a0fa07a-3605-11eb-30ef-4598213ce6ec
# ╠═08cd4892-362c-11eb-28d8-81540aa5945e
# ╠═34846d78-362f-11eb-3547-4978e97e89d9
# ╠═46ee77ec-362f-11eb-3c13-f384ed78a348
# ╠═f08c3ad2-362f-11eb-295d-4ddac22b387f
# ╠═673d8218-366b-11eb-2a4d-f548f9d131ee
# ╠═8c68eb76-362e-11eb-300f-6b92a90010fa
# ╠═4b3eb86a-3630-11eb-39c0-6d67c0f84a22
# ╠═53200a48-3630-11eb-05e9-99811ae65d11
# ╠═5f797e22-339d-11eb-1ff3-dda90f946838
# ╠═6a522dd4-3641-11eb-1c9e-45074da49c82
# ╠═9ac14414-3641-11eb-33c4-75b9e4357068
# ╠═5189e182-3664-11eb-0dcc-f508fe33be4c
# ╠═6947baf8-3664-11eb-0ace-3f6c7f808068
# ╠═694808aa-3664-11eb-0443-c162e7388740
# ╠═694e56b0-3664-11eb-0d26-b9e7b641b9f6
# ╠═694ee300-3664-11eb-11f7-49462e6aa1d2
# ╠═69599412-3664-11eb-26c7-21afa8e93b9c
# ╠═695a2420-3664-11eb-0978-5be1861a5f49
# ╠═69657fa2-3664-11eb-3056-a98b35ba51e3
# ╠═696628c6-3664-11eb-2643-fd542cd9eab5
# ╠═6971b038-3664-11eb-2c31-c5245b309553
# ╠═6972e2b4-3664-11eb-1052-af8fccaeac05
# ╠═697f53fa-3664-11eb-3962-030d46f24ba7
# ╠═69800c82-3664-11eb-2d6e-155e57cdc323
# ╠═edc481a6-3664-11eb-1c02-797a3670af3f
# ╠═3da5c12e-3669-11eb-1a79-73e9fad032e8
# ╠═05f62636-366b-11eb-1435-632abff666d1
# ╠═42b157f0-3669-11eb-2abc-f16c5da1a521
# ╠═dfa0bb9c-3668-11eb-26ae-59f47e74033f
# ╠═fd483df0-3668-11eb-09f4-bb26699ab6bf
# ╠═d95d9cbc-3665-11eb-131c-3f53f873ce71
# ╠═5b8c1242-3669-11eb-2b51-2d878f04f456
# ╠═dfff04d0-3669-11eb-37a8-67ede57fa4f6

max_stars_repo_name: RensDimmendaal/advent-of-code-julia

max_issues_repo_name: RensDimmendaal/advent-of-code-julia

max_forks_repo_name: RensDimmendaal/advent-of-code-julia


text:
module MPIPreferences

using Preferences, Libdl

"""
    MPIPreferences.binary :: String

The currently selected binary. The possible values are

- `"MPICH_jll"`: use the binary provided by [MPICH_jll](https://github.com/JuliaBinaryWrappers/MPICH_jll.jl)
- `"OpenMPI_jll"`: use the binary provided by [OpenMPI_jll](https://github.com/JuliaBinaryWrappers/OpenMPI_jll.jl)
- `"MicrosoftMPI_jll"`: use binary provided by [MicrosoftMPI_jll](https://github.com/JuliaBinaryWrappers/MicrosoftMPI_jll.jl/)
- `"MPItrampoline_jll"`: use the binary provided by [MPItrampoline_jll](https://github.com/JuliaBinaryWrappers/MPItrampoline_jll.jl/)
- `"system"`: use a system-provided binary.

"""
const binary = @load_preference("binary", Sys.iswindows() ? "MicrosoftMPI_jll" : "MPICH_jll")

"""
    MPIPreferences.abi :: String

The ABI of the currently selected binary. Supported values are:

- `"MPICH"`: MPICH-compatible ABI (https://www.mpich.org/abi/)
- `"OpenMPI"`: Open MPI compatible ABI (Open MPI, IBM Spectrum MPI, Fujitsu MPI)
- `"MicrosoftMPI"`: Microsoft MPI
- `"MPItrampoline"`: MPItrampoline
- `"HPE MPT"`: HPE MPT
"""
const abi = if binary == "system"
    @load_preference("abi")
elseif binary == "MicrosoftMPI_jll"
    "MicrosoftMPI"
elseif binary == "MPICH_jll"
    "MPICH"
elseif binary == "OpenMPI_jll"
    "OpenMPI"
elseif binary == "MPItrampoline_jll"
    "MPItrampoline"
else
    error("Unknown binary: $binary")
end

@static if binary == "system"
    include("system.jl")
end

"""
    MPIPreferences.use_jll_binary(binary; export_prefs=false, force=true)

Switches the underlying MPI implementation to one provided by JLL packages.
Available options are:
- `"MicrosoftMPI_jll"` (Only option and default on Winddows)
- `"MPICH_jll"` (Default on all other platform)
- `"OpenMPI_jll"`
- `"MPItrampoline_jll"`

The `export_prefs` option determines whether the preferences being set should be stored
within `LocalPreferences.toml` or `Project.toml`.
"""
function use_jll_binary(binary = Sys.iswindows() ? "MicrosoftMPI_jll" : "MPICH_jll"; export_prefs=false, force=true)
    binary in ["MicrosoftMPI_jll", "MPICH_jll", "OpenMPI_jll", "MPItrampoline_jll"] ||
        error("Unknown jll: $binary")
    set_preferences!(MPIPreferences,
        "binary" => binary,
        "libmpi" => nothing,
        "abi" => nothing,
        "mpiexec" => nothing;
        export_prefs=export_prefs,
        force=force
    )

    @warn "The underlying MPI implementation has changed. You will need to restart Julia for this change to take effect" binary

    if VERSION <= v"1.6.5" || VERSION == v"1.7.0"
        @warn """
        Due to a bug in Julia (until 1.6.5 and 1.7.1), setting preferences in transitive dependencies
        is broken (https://github.com/JuliaPackaging/Preferences.jl/issues/24). To fix this either update
        your version of Julia, or add MPIPreferences as a direct dependency to your project.
        """
    end

    return nothing
end

"""
    MPIPreferences.use_system_binary(;
        library_names = ["libmpi", "libmpi_ibm", "msmpi", "libmpich", "libmpitrampoline"],
        mpiexec = "mpiexec",
        abi = nothing,
        export_prefs = false,
        force = true)

Switches the underlying MPI implementation to a system provided one.

Options:

- `library_names`: a name or collection of names of the MPI library, passed to
  [`Libdl.find_library`](https://docs.julialang.org/en/v1/stdlib/Libdl/#Base.Libc.Libdl.find_library).
  If the library isn't in the library search path, you can specify the full path
  to the library.

- `mpiexec`: the MPI launcher executable. The default is `mpiexec`, but some
  clusters require using the scheduler launcher interface (e.g. `srun` on Slurm,
  `aprun` on PBS). It is also possible to pass a [`Cmd`
  object](https://docs.julialang.org/en/v1/manual/running-external-programs/#Cmd-Objects)
  to include specific command line options.

- `abi`: the ABI of the MPI library. By default this is determined automatically
  using [`identify_abi`](@ref). See [`abi`](@ref) for currently supported values.

- `export_prefs`: if `true`, the preferences into the `Project.toml` instead of `LocalPreferences.toml`.

- `force`: if `true`, the preferences are set even if they are already set.
"""
function use_system_binary(;
        library_names=["libmpi", "libmpi_ibm", "msmpi", "libmpich", "libmpitrampoline"],
        mpiexec="mpiexec",
        abi=nothing,
        export_prefs=false,
        force=true,
    )
    binary = "system"
    # Set `ZES_ENABLE_SYSMAN` to work around https://github.com/open-mpi/ompi/issues/10142
    libmpi = withenv("ZES_ENABLE_SYSMAN" => "1") do
        find_library(library_names)
    end
    if libmpi == ""
        error("MPI library could not be found")
    end
    if isnothing(abi)
        abi = identify_abi(libmpi)
    end
    if mpiexec isa Cmd
        mpiexec = collect(mpiexec)
    end
    set_preferences!(MPIPreferences,
        "binary" => binary,
        "libmpi" => libmpi,
        "abi" => abi,
        "mpiexec" => mpiexec,
        export_prefs=export_prefs,
        force=force
    )

    @warn "The underlying MPI implementation has changed. You will need to restart Julia for this change to take effect" binary libmpi abi mpiexec

    if VERSION <= v"1.6.5" || VERSION == v"1.7.0"
        @warn """
        Due to a bug in Julia (until 1.6.5 and 1.7.1), setting preferences in transitive dependencies
        is broken (https://github.com/JuliaPackaging/Preferences.jl/issues/24). To fix this either update
        your version of Julia, or add MPIPreferences as a direct dependency to your project.
        """
    end

    return nothing
end

function identify_implementation_version_abi(version_string::AbstractString)
    impl = "unknown"
    version = v"0"

    if startswith(version_string, "MPICH")
        impl = "MPICH"
        # "MPICH Version:\t%s\n" /  "MPICH2 Version:\t%s\n"
        if (m = match(r"^MPICH2? Version:\t(\d+.\d+(?:.\d+)?\w*)\n", version_string)) !== nothing
            version = VersionNumber(m.captures[1])
        end

    elseif startswith(version_string, "Open MPI")
        # Open MPI / Spectrum MPI
        impl = occursin("IBM Spectrum MPI", version_string) ? "IBMSpectrumMPI" : "OpenMPI"

        if (m = match(r"^Open MPI v(\d+.\d+.\d+\w*)", version_string)) !== nothing
            version = VersionNumber(m.captures[1])
        end

    elseif startswith(version_string, "Microsoft MPI")
        impl = "MicrosoftMPI"
        # "Microsoft MPI %u.%u.%u.%u%S"
        # ignore last 2 (build numbers)
        if (m = match(r"^Microsoft MPI (\d+.\d+)", version_string)) !== nothing
            version = VersionNumber(m.captures[1])
        end

    elseif startswith(version_string, "Intel")
        impl = "IntelMPI"

        # "Intel(R) MPI Library 2019 Update 4 for Linux* OS"
        if (m = match(r"^Intel\(R\) MPI Library (\d+)(?: Update (\d+))?", version_string)) !== nothing
            if m.captures[2] === nothing
                version = VersionNumber(m.captures[1])
            else
                version = VersionNumber(m.captures[1]*"."*m.captures[2])
            end
        end

    elseif startswith(version_string, "MVAPICH2")
        impl = "MVAPICH"
        # "MVAPICH2 Version      :\t%s\n")
        if (m = match(r"^MVAPICH2? Version\s*:\t(\S*)\n", version_string)) !== nothing
            version = VersionNumber(m.captures[1])
        end

    elseif occursin("CRAY MPICH", version_string)
        impl = "CrayMPICH"
        # "MPI VERSION    : CRAY MPICH version 7.7.10 (ANL base 3.2)\n"
        if (m = match(r"CRAY MPICH version (\d+.\d+.\d+)", version_string)) !== nothing
            version = VersionNumber(m.captures[1])
        end

    elseif startswith(version_string, "FUJITSU MPI")
        impl = "FujitsuMPI"
        # "FUJITSU MPI Library 4.0.0 (4.0.1fj4.0.0)\0"
        if (m = match(r"^FUJITSU MPI Library (\d+.\d+.\d+)", version_string)) !== nothing
            version = VersionNumber(m.captures[1])
        end

    elseif startswith(version_string, "MPIwrapper")
        impl = "MPIwrapper"
        # MPIwrapper 2.2.2
        if (m = match(r"^MPIwrapper Version:\t(\d+.\d+.\d+\w*)", version_string)) !== nothing
            version = VersionNumber(m.captures[1])
        end

    elseif startswith(version_string, "HPE MPT")
        impl = "HPE MPT"
        # HPE MPT 2.23  08/26/20 02:54:49-root
        if (m = match(r"^HPE MPT (\d+.\d+)", version_string)) !== nothing
            version = VersionNumber(m.captures[1])
        end

    elseif startswith(version_string, "HPE HMPT")
        impl = "HPE HMPT"
        # HPE HMPT 2.23  08/26/20 02:59:48-root
        if (m = match(r"^HPE HMPT (\d+.\d+)", version_string)) !== nothing
            version = VersionNumber(m.captures[1])
        end
    end

    # 3) determine the abi from the implementation + version
    if (impl == "MPICH" && version >= v"3.1" ||
        impl == "IntelMPI" && version > v"2014" ||
        impl == "MVAPICH" && version >= v"2" ||
        impl == "CrayMPICH" && version >= v"7" ||
        # https://www.mpich.org/abi/
        impl == "HPE HMPT")
        abi = "MPICH"
    elseif impl == "OpenMPI" || impl == "IBMSpectrumMPI" || impl == "FujitsuMPI"
        abi = "OpenMPI"
    elseif impl == "MicrosoftMPI"
        abi = "MicrosoftMPI"
    elseif impl == "MPIwrapper"
        abi = "MPItrampoline"
    elseif impl == "HPE MPT"
        abi = "HPE MPT"
    else
        abi = "unknown"
    end

    return (impl, version, abi)
end

"""
    identify_abi(libmpi)

Identify the MPI implementation from the library version string
"""
function identify_abi(libmpi)
    # 1) query MPI_Get_version
    # There is no way to query at runtime what the length of the buffer should be.
    # https://github.com/mpi-forum/mpi-issues/issues/159
    # 8192 is the maximum value of MPI_MAX_LIBRARY_VERSION_STRING across known
    # implementations.
    buf = Array{UInt8}(undef, 8192)
    buflen = Ref{Cint}()

    hndl = @static if Sys.isunix()
        # Again, work around https://github.com/open-mpi/ompi/issues/10142
        withenv("ZES_ENABLE_SYSMAN" => "1") do
            # need to open libmpi with RTLD_GLOBAL flag for Linux, before
            # any ccall cannot use RTLD_DEEPBIND; this leads to segfaults
            # at least on Ubuntu 15.10
            Libdl.dlopen(libmpi, Libdl.RTLD_LAZY | Libdl.RTLD_GLOBAL)
        end
    else
        Libdl.dlopen(libmpi)
    end

    try
        ptr = dlsym(hndl, :MPI_Get_library_version)
        ccall(ptr, Cint, (Ptr{UInt8}, Ref{Cint}), buf, buflen)
    finally
        Libdl.dlclose(hndl)
    end

    @assert buflen[] < 8192
    resize!(buf, buflen[])
    version_string = String(buf)

    # 2) try to identify the MPI implementation
    impl, version, abi = identify_implementation_version_abi(version_string)

    @info "MPI implementation" libmpi version_string impl version abi

    return abi
end

end # module

ext: jl


text:
q2vec(q::UnitQuaternion) = return [q.w, q.x, q.y, q.z]

getStates(mechanism::Mechanism) = return [deepcopy(mechanism.bodies[id].state) for id in 1:length(mechanism.bodies)]


text:
function mixedparsevaltestxk(L,N,j)
X,K,dX,dK = makearrays(L,N)
T = makeT(X,K)

ψ = randn(N...) + im*randn(N...)

n1 = sum(abs2.(ψ))*prod(dX)

ϕ = T.Mxk[j]*ψ

n2 = sum(abs2.(ϕ))*prod(dX)*dK[j]/dX[j]
return n1,n2
end

function mixedparsevaltestkx(L,N,j)
X,K,dX,dK = makearrays(L,N)
T = makeT(X,K)

ϕ = randn(N...) + im*randn(N...)

n1 = sum(abs2.(ϕ))*prod(dK)

ψ = T.Mkx[j]*ϕ

n2 = sum(abs2.(ψ))*prod(dK)*dX[j]/dK[j]
return n1,n2
end

L = (88.)
N = (128)
n1,n2 = mixedparsevaltestxk(L,N,1)
@test n1 ≈ n2

L = (30.,20.)
N = (30,34)
n1,n2 = mixedparsevaltestxk(L,N,1)
@test n1 ≈ n2

n1,n2 = mixedparsevaltestxk(L,N,2)
@test n1 ≈ n2

L = (30.,20.,10)
N = (30,34,24)
n1,n2 = mixedparsevaltestxk(L,N,1)
@test n1 ≈ n2

n1,n2 = mixedparsevaltestxk(L,N,2)
@test n1 ≈ n2

n1,n2 = mixedparsevaltestxk(L,N,3)
@test n1 ≈ n2

L = (30.,20.,10,50.)
N = (30,34,24,22)
n1,n2 = mixedparsevaltestxk(L,N,1)
@test n1 ≈ n2

n1,n2 = mixedparsevaltestxk(L,N,2)
@test n1 ≈ n2

n1,n2 = mixedparsevaltestxk(L,N,3)
@test n1 ≈ n2

n1,n2 = mixedparsevaltestxk(L,N,4)
@test n1 ≈ n2

L = (88.)
N = (128)
n1,n2 = mixedparsevaltestkx(L,N,1)
@test n1 ≈ n2

L = (30.,20.)
N = (30,34)
n1,n2 = mixedparsevaltestkx(L,N,1)
@test n1 ≈ n2

n1,n2 = mixedparsevaltestkx(L,N,2)
@test n1 ≈ n2

L = (30.,20.,10)
N = (30,34,24)
n1,n2 = mixedparsevaltestkx(L,N,1)
@test n1 ≈ n2

n1,n2 = mixedparsevaltestkx(L,N,2)
@test n1 ≈ n2

n1,n2 = mixedparsevaltestkx(L,N,3)
@test n1 ≈ n2

L = (30.,20.,10,50.)
N = (30,34,24,22)
n1,n2 = mixedparsevaltestkx(L,N,1)
@test n1 ≈ n2

n1,n2 = mixedparsevaltestkx(L,N,2)
@test n1 ≈ n2

n1,n2 = mixedparsevaltestkx(L,N,3)
@test n1 ≈ n2

n1,n2 = mixedparsevaltestkx(L,N,4)
@test n1 ≈ n2


text:
module AccurateArithmetic

# Adapted from AccurateArithmetic.jl

export two_sum, apply_offset

@inline function two_sum(a, b)
    hi = a + b
    a1 = hi - b
    b1 = hi - a1
    lo = (a - a1) + (b - b1)
    return hi, lo
end

@inline function two_hilo_sum(a, b)
    hi = a + b
    lo = b - (hi - a)
    return hi, lo
end

function four_sum(a, b, c, d)
    t0, t1 = two_sum(a,  b)
    t2, t3 = two_sum(c,  d)
    hi, t4 = two_sum(t0, t2)
    t5, lo = two_sum(t1, t3)
    hm, ml = two_sum(t4, t5)
    ml, lo = two_hilo_sum(ml, lo)
    hm, ml = two_hilo_sum(hm, ml)
    hi, hm = two_hilo_sum(hi,hm)
    return hi, hm, ml, lo
end

function handle_infinity(fraction)
    second = ifelse(fraction < 0, typemin(Int64), typemax(Int64))
    return second, fraction, zero(fraction)
end

function apply_offset(s1::Int64, f1, e1, s2::Int64, f2, e2)
    isfinite(f1 + f2) || return handle_infinity(f1 + f2)

    sum, residual, _ = four_sum(f1, f2, e1, e2)
    int_seconds = floor(Int64, sum)
    second = s1 + s2 + int_seconds
    fraction = sum - int_seconds
    return second, fraction, residual
end

function apply_offset(s1::Int64, f1, e1, offset)
    isfinite(f1 + offset) || return handle_infinity(f1 + offset)

    s2 = floor(Int64, offset)
    f2 = offset - s2
    return apply_offset(s1, f1, e1, s2, f2, 0.0)
end

end

max_stars_repo_name: JuliaAstrodynamics/AstronomicalTime

max_issues_repo_name: JuliaAstrodynamics/AstronomicalTime

max_forks_repo_name: JuliaAstrodynamics/AstronomicalTime


text:
# One dimensional wave equation
# Test systems of conservation laws

@testset "1D Sytems Algorithms" begin

using ConservationLawsDiffEq
using OrdinaryDiffEq
CL = ConservationLawsDiffEq

CFL = 0.45
Tend = 1.0
cc = 1.0

Jf(u::AbstractVector) = [0.0 cc;cc 0.0]
f(u::AbstractVector) = [0.0 cc;cc 0.0]*u
f0(x) = [sin(4*π*x), 0.0]

Nflux(ϕl, ϕr) = 0.5*(f(ϕl)+f(ϕr))
exact_sol(x::Float64, t::Float64) = hcat(0.5*(sin.(4*π*(-t+x))+sin.(4*π*(t+x))),
0.5*(sin.(4*π*(-t+x))-sin.(4*π*(t+x))))

mesh = Uniform1DFVMesh(20, [-1.0, 1.0])

function get_problem(alg::CL.AbstractFVAlgorithm, use_threads, mesh, Tend,CFL)
  # Now get a explicit semidiscretization (discrete in space) du_h(t)/dt = f_h(u_h(t))
  f_h = getSemiDiscretization(f,alg,mesh,[Periodic()]; Df = Jf, use_threads = use_threads,numvars = 2)

  #Compute discrete initial data
  u0 = getInitialState(mesh,f0,use_threads = true)

  #Setup ODE problem for a time interval = [0.0,1.0]
  ode_prob = ODEProblem(f_h,u0,(0.0,Tend))

  #Setup callback in order to fix CFL constant value
  cb = getCFLCallback(f_h, CFL)

  #Estimate an initial dt
  dt = update_dt!(u0, f_h, CFL)
  return ode_prob, cb, dt
end

function test_scheme(alg, maxerror, threaded_mode=true)
  println("Testing ",split("$alg","(")[1],": ")
  prob,cb,dt = get_problem(alg, false, mesh, Tend,CFL)
  @time sol = solve(prob,SSPRK22(); dt = dt, callback = cb, save_everystep = false)
  @test get_L1_error(exact_sol, fv_solution(sol,mesh)) < maxerror
  if threaded_mode
    println("in threaded mode: ")
    prob,cb,dt = get_problem(alg, true, mesh, Tend,CFL)
    @time sol1 = solve(prob,SSPRK22(); dt = dt, callback = cb, save_everystep = false)
    @test get_L1_error(exact_sol, fv_solution(sol1,mesh)) ≈ get_L1_error(exact_sol, fv_solution(sol,mesh))
  end
end

# Test Schemes error
test_scheme(FVSKTScheme(), 1.16)
test_scheme(FVTecnoScheme(Nflux;order=2), 1.152, false)
test_scheme(FVTecnoScheme(Nflux;order=3), 1.14, false)
test_scheme(FVTecnoScheme(Nflux;order=5), 0.92, false)
test_scheme(LaxFriedrichsScheme(), 1.16)
test_scheme(LocalLaxFriedrichsScheme(), 1.16, false)
test_scheme(GlobalLaxFriedrichsScheme(), 1.16)
test_scheme(LaxWendroffScheme(), 37.0)    #way off
test_scheme(LaxWendroff2sScheme(), 1.152)
test_scheme(FVCompWENOScheme(), 0.94)
test_scheme(FVCompMWENOScheme(), 0.68)
test_scheme(FVSpecMWENOScheme(), 0.68, false)
test_scheme(FVCUScheme(), 1.16)
test_scheme(FVDRCUScheme(), 1.17)
test_scheme(FVDRCU5Scheme(), 0.48)

end


text:
module XRFast

export nnls, subX, lasso, lasso_sparse, approximate_ksvd, dictionary_learning, dictionary_learning_1, LLS, INV_LLS, Savitsky_Golay, SNIP, initial_dictionary, approx_nn_dictionary_learning, full_nn_dictionary_learning

using Lasso
using LinearAlgebra
using DSP
using ProgressMeter
using SparseArrays
using NonNegLeastSquares
using Statistics

#subX function randomly selects 5% of spectra from original data cube and saves as a new array
#Spectral dimensions in y (dim 1), pixels in x (dim 2)
function subX(a)
    channels = size(a,2)
    pixels = size(a,1)
    percent = trunc(Int64, (size(a,1))*0.05) #select 5% of pixels to form
    subX =zeros(0) #empty array
    for i in 1:percent #loop to create array of randomly selected pixels. Stored in an array called subX
        n = size(a,1)
        idx = rand(1:n)
        append!(subX,a[idx,:])
    end
    return reshape(subX,channels,percent)
end

function initial_dictionary(a;
    Atoms::Int = default_Atoms)

    default_Atoms = 40

    channels = size(a,2)
    pixels = size(a,1)
    Prc = Atoms/pixels

    percent = trunc(Int64, (size(a,1))*Prc)
    subX =zeros(0) #empty array
    for i in 1:percent
        n = size(a,1)
        idx = rand(1:n)
        append!(subX,a[idx,:])
    end
    return reshape(subX,channels,percent)
end

function lasso_sparse(subx, D, Threshold)
    #p = size(subx,2) # numb of pixels in image
    w = zeros(0)
    for i in 1:size(subx,2)
        q = subx[:,i]
        W = fit(LassoPath, D, q;intercept=false, cd_tol=1e-5)
        J = Matrix(coef(W))
        append!(w,J[:,size(J,2)])
    end
    W = reshape(w,size(D,2),size(subx,2))
    T = [broadcast(abs, x) > Threshold ? x : 0 for x in W] #Threshold parameter to make sparse, necessary?
    T =[val*1.0 for val in T]
    return T
end

function lasso(subx, D, cd_tolerance)

    #p = size(subx,2) # numb of pixels in image
    w = zeros(0)
    for i in 1:size(subx,2)
        q = subx[:,i]
        W = fit(LassoPath, D, q;intercept=false, cd_tol=cd_tolerance)
        J = Matrix(coef(W))
        append!(w,J[:,size(J,2)])
    end
    W = reshape(w,size(D,2),size(subx,2))
    return W
end

function approximate_ksvd(Y::AbstractMatrix, D::AbstractMatrix, W::AbstractMatrix)
    R = Y - D*W
    for k in 1:size(D,2)
        I = findall( x->(x > 0), W[k,:])
        Ri = R[:,I] + D[:,k]*W[k,I']
        dk = Ri * W[k,I]
        dk = dk/sqrt(dk'*dk)  # normalize
        D[:,k] = dk
        W[k,I] = dk'*Ri;
        R[:,I] = Ri - D[:,k]*W[k,I']
    end
    return D
end

function nn_approximate_ksvd(Y::AbstractMatrix, D::AbstractMatrix, W::AbstractMatrix)
    R = Y - D*W
    for k in 1:size(D,2)
        I = findall( x->(x > 0), W[k,:])
	    #R[:,I] = [x > 0 ? x : 0 for x in R[:,I]]
        Ri = R[:,I] + D[:,k]*W[k,I']
        dk = Ri * W[k,I]
        dk = dk/sqrt(dk'*dk)  # normalize
	    @. dk[isnan(dk)] = 0
        D[:,k] = dk
        D[:,k] = [x > 0 ? x : 0 for x in D[:,k]]
	    W[k,I] = D[:,k]'*Ri;
        R[:,I] = Ri - D[:,k]*W[k,I']
    end
    return D
end

function nn_full_ksvd(Y::AbstractMatrix, D::AbstractMatrix, W::AbstractMatrix)
    R = Y - D*W
    for k in 1:size(D,2)
        I = findall( x->(x > 0), W[k,:])
        Ri = R[:,I] + D[:,k]*W[k,I']
        U,S,V = svd(Ri)
        #@. U[isnan(U)] = 0
        D[:,k] = U[:,1]
        D[:,k] = [x > 0 ? x : 0 for x in D[:,k]]
        W[k,I] = V*S
        R[:,I] = Ri - D[:,k]*W[k,I']
    end
    return D
end

function nnls(subx, D, Lambda)
    w = zeros(0)
    for i in 1:size(subx,2)
        q = subx[:,i]
        J = nonneg_lsq(D,subx[:,i];alg=:nnls)
        J = [x > quantile!(J[:,1], Lambda) ? x : 0 for x in J]
        append!(w,J[:,size(J,2)])
    end
    W = reshape(w,size(D,2),size(subx,2))
    return W
end

function dictionary_learning(D::AbstractMatrix, L::AbstractMatrix;
        noIt_Sub::Int = default_It_Sub,
        noIt_KSVD::Int = default_It_KSVD,
        Threshold::Float64 = default_Thresh)

    default_It_Sub = 10
    default_It_KSVD = 1
    default_Thresh = 0.25

    p = Progress(noIt_Sub, 1, "Optimizing Dictionary...")
    for i in 1:noIt_Sub
        y = subX(L)
        for j in 1:noIt_KSVD
            W = lasso_sparse(y,D,Threshold)
            D = approximate_ksvd(y, D, W)
        end
    next!(p)
    end
    return(D)
end

function dictionary_learning_1(D::AbstractMatrix, L::AbstractMatrix;
        noIt_Sub::Int = default_It_Sub,
        noIt_KSVD::Int = default_It_KSVD,
        cd_tolerance = default_It_KSVD)

    default_It_Sub = 10
    default_It_KSVD = 1
    default_cd_tolerance = 1e-6

    p = Progress(noIt_Sub, 1, "Optimizing Dictionary...")
    for i in 1:noIt_Sub
        y = subX(L)
        for j in 1:noIt_KSVD
            W = lasso(y,D,cd_tolerance)
            D = approximate_ksvd(y, D, W)
        end
    next!(p)
    end
    return(D)
end

function approx_nn_dictionary_learning(D::AbstractMatrix, L::AbstractMatrix;
        noIt_Sub::Int = default_It_Sub,
        noIt_KSVD::Int = 1,
        Lambda::Float64= default_Th)

    default_It_Sub = 10
    default_It_KSVD = 1
    default_Th = 0.75

    p = Progress(noIt_Sub, 1, "Optimizing Dictionary...")
    for i in 1:noIt_Sub
        y = subX(L)
        for j in 1:noIt_KSVD
            W = nnls(y,D,Lambda)
            D = nn_approximate_ksvd(y, D, W)
        end
    next!(p)
    end
    return(D)
end

function full_nn_dictionary_learning(D::AbstractMatrix, L::AbstractMatrix;
        noIt_Sub::Int = default_It_Sub,
        noIt_KSVD::Int = = 1,
        Lambda::Float64= default_Th)

    default_It_Sub = 10
    default_It_KSVD = 1
    default_Th = 0.75

    p = Progress(noIt_Sub, 1, "Optimizing Dictionary...")
    for i in 1:noIt_Sub
        y = subX(L)
        for j in 1:noIt_KSVD
            W = nnls(y,D,Lambda)
            D = nn_full_ksvd(y, D, W)
        end
    next!(p)
    end
    return(D)
end

function LLS(a)
    channels = size(a,2)
    pixels = size(a,1)
    LLS =zeros(0) #empty array
    for i in 1:size(a,1)
        T = [x > 0 ? x : 0 for x in a[i,:]]
        dr = [log(log(sqrt(val+1)+1)+1) for val in T]
        append!(LLS,dr)
    end
    return reshape(LLS,channels,pixels)
end

function INV_LLS(a)
    channels = size(a,2)
    pixels = size(a,1)
    INV_LLS =zeros(0) #empty array
    for i in 1:size(a,1)
        df = [(exp(exp(val) - 1) - 1)^2 - 1 for val in a[i,:]]
        append!(INV_LLS,df)
    end
    return reshape(INV_LLS,channels,pixels)
end

function Savitsky_Golay(d::AbstractMatrix;
        windowSize::Int = default_windowSize,
        polyOrder::Int = default_polyOrder,
        deriv::Int = default_deriv)
    pixels = size(d,1) #number of pixels in image
    default_windowSize = 33 #must be odd
    default_polyOrder = 2 #must be less than the window size
    default_deriv = 0
    for i in 1:pixels
        x = d[i,:]
        isodd(windowSize) || throw("Window size must be an odd integer.")
        polyOrder < windowSize || throw("Polynomial order must me less than window size.")
        halfWindow = Int( ceil((windowSize-1)/2) )
        S = zeros.(windowSize, polyOrder+1)
        for ct = 0:polyOrder
        S[:,ct+1] = (-halfWindow:halfWindow).^(ct)
        end
        G = S * pinv(S' * S)
        filterCoeffs = G[:, deriv+1] * factorial(deriv)
        paddedX = [x[1]*ones(halfWindow); x; x[end]*ones(halfWindow)]
        y = conv(filterCoeffs[end:-1:1], paddedX)
        xj = y[2*halfWindow+1:end-2*halfWindow]
        d[i,:] = [yal > 0 ? yal : 0 for yal in xj]
    end
    return(d)
end

function SNIP(d::AbstractMatrix)
    M = 36 # distance of M
    iter = 15; # number of iterations

    Mb = 14 # distance of M
    iterb = 5

    Mc = 7 # distance of M
    iterc = 4

    Md = 3 # distance of M
    iterd = 2

    for j in 1:iter, i in M+1:4096-M
            M1 = d[:,i-M]
            M2 = d[:,i+M]
            J = [(x+y)/2 for (x,y) in zip(M1, M2)]
            d[:,i] = [x < y ? x : y for (x, y) in zip(J, d[:,i])]
    end
    for j in 1:iterb, i in Mb+1:4096-Mb
            M1 = d[:,i-Mb]
            M2 = d[:,i+Mb]
            J = [(x+y)/2 for (x,y) in zip(M1, M2)]
            d[:,i] = [x < y ? x : y for (x, y) in zip(J, d[:,i])]
    end
    for j in 1:iterc, i in Mc+1:4096-Mc
            M1 = d[:,i-Mc]
            M2 = d[:,i+Mc]
            J = [(x+y)/2 for (x,y) in zip(M1, M2)]
            d[:,i] = [x < y ? x : y for (x, y) in zip(J, d[:,i])]
    end
    for j in 1:iterd, i in Md+1:4096-Md
            M1 = d[:,i-Md]
            M2 = d[:,i+Md]
            J = [(x+y)/2 for (x,y) in zip(M1, M2)]
            d[:,i] = [x < y ? x : y for (x, y) in zip(J, d[:,i])]
    end
    return(d)
end

end # module

max_stars_repo_name: NU-ACCESS/XRFast

max_issues_repo_name: NU-ACCESS/XRFast

max_forks_repo_name: NU-ACCESS/XRFast


text:
# class constructors
Quotient(X::Integer,::Type{T}) where T<:Integer = T / T(X)
Polynomial(::Type{Q}) where {P,Q<:Quotient{P}} = P

# convenience type constructor
# enable `Z / m` for anonymous quotient class constructor
function /(::Type{R}, m) where R<:Ring

    ideal = pseudo_ideal(R, m)
    p, r = characteristic(R), deg(ideal)
    b = order(basetype(R))
    new_class(Quotient{R,typeof(ideal),sintern(m),(p,r)}, ideal)
end

# Constructors
category_trait(Z::Type{<:QuotientRing}) = isprimemod(Z) ? FieldTrait : CommutativeRingTrait
basetype(::Type{<:Quotient{T}}) where T = T
depth(::Type{<:Quotient{T}}) where T = depth(T) + 1

function Quotient{R,I,X,Id}(a::R) where {I,X,R<:Ring,Id}
    m = modulus(Quotient{R,I,X,Id})
    v = rem(a, m)
    Quotient{R,I,X,Id}(v, NOCHECK)
end

# convert argument to given R
Quotient{R,I,X,Id}(v::Quotient{R,I,X,Id}) where {I,X,R<:Ring,Id} = Quotient{R,I,X,Id}(v.val)
Quotient{R,I,X,Id}(v) where {I,X,R<:Ring,Id} = Quotient{R,I,X,Id}(R(v))

# promotion and conversion
_promote_rule(::Type{<:Quotient}, ::Type{<:Quotient}) = Base.Bottom
_promote_rule(::Type{Quotient{R,I,X,Id}}, ::Type{S}) where {I,X,R,S<:Ring,Id} = Quotient{promote_type(R,S),I,X,Id}
promote_rule(::Type{Quotient{R,I,X,Id}}, ::Type{S}) where {I,X,R,S<:Integer,Id} = Quotient{R,I,X,Id}

(::Type{Q})(a::S) where {S,R,Q<:Quotient{R}} = Q(R(a))

Base.isless(p::T, q::T) where T<:Quotient = isless(p.val, q.val)

"""
    iszerodiv(q::R) where R<:Ring

Return true iff q is not a zero divisor in `R`
"""
iszerodiv(q::Q) where Q<:QuotientRing = !isone(gcd(value(q), modulus(Q)))
iszerodiv(q::R) where R<:Union{Ring,Number} = iszero(q)

## Arithmetic

==(a::T, b::T) where T<:Quotient =  a.val == b.val
==(a::Quotient, b::Quotient) =  false
+(a::T, b::T) where T<:Quotient =  T(a.val + b.val)
-(a::T, b::T) where T<:Quotient =  T(a.val - b.val)
*(a::T, b::T) where T<:Quotient =  T(a.val * b.val)
*(a::Integer, b::T) where T<:Quotient =  T(a * b.val)
*(a::T, b::Integer) where T<:Quotient =  T(a.val * b)
-(a::T) where T<:Quotient =  T(-a.val)
inv(a::T) where T<:Quotient = T(invert(a.val, modulus(T)), NOCHECK)

isunit(a::T) where T<:Quotient = isunit(a.val) || isinvertible(modulus(T), a.val)
iszero(x::Quotient) = iszero(x.val)
isone(x::Quotient) = isone(x.val)
zero(::Type{Q}) where {S,Q<:Quotient{S}} = Q(zero(S), NOCHECK)
one(::Type{Q}) where {S,Q<:Quotient{S}} = Q(one(S), NOCHECK)
value(a::QuotientRing) = a.val
characteristic(::Type{Quotient{R,I,X,Id}}) where {R,I,X,Id} = Id[1]
dimension(::Type{Quotient{R,I,X,Id}}) where {R,I,X,Id} = Id[2]
function order(::Type{Quotient{R,I,X,Id}}) where {R,I,X,Id}
    r = Id[2]
    b = order(basetype(R))
    iszero(r * b) ? 0 : uptype(intpower(b, r), Int)
end

# induced homomorphism - invalid if Q = R/I and I not in kernel(F)
function (h::Hom{F,R,S})(a::Q) where {F,R,S,Q<:Quotient{<:R}}
    iszero(h.f(modulus(Q))) || throw(DomainError((F,R), "ideal not in kernel of homomorphism"))
    h.f(a.val)
end

# note:
# the real work is in the functions `Ideal`, `rem`, `invert`, `isinvertible` which
# have all been delegated to Ideal

## Help functions

# return the ideal associated uniquely with this quotient ring
modulus(t::Type{<:Quotient{R}}) where R = gettypevar(t).modulus

# standard functions
==(a::Quotient{S,I,X},b::Quotient{T,I,X}) where {I,X,S,T} = a.val == b.val
hash(a::Quotient, h::UInt) = hash(a.val, hash(modulus(a), h))

function Base.show(io::IO, a::Quotient)
    v = a.val
    m = modulus(a)
    if m isa UnivariatePolynomial && deg(m) == 2 &&
         iszero(m.coeff[2]) && isone(m.coeff[3])

        x = string(varnames(m)[1])
        y = string('\u23b7', -m.coeff[1])
        vs = replace(sprint(show, v), x => y)
        print(io, vs)
    else
        print(io, v, " mod(", m, ")")
    end
end

//(a::G, b::G) where G<:QuotientRing = a / b
div(a::G, b::G) where G<:QuotientRing = a / b
rem(a::G, b::G) where G<:QuotientRing = zero(G)
gcd(a::G, b::G) where G<:QuotientRing = one(G)
gcdx(a::G, b::G) where G<:QuotientRing = one(G), zero(G), zero(G)
pgcd(a::G, b::G) where G<:QuotientRing = gcd(a, b)
pgcdx(a::G, b::G) where G<:QuotientRing = gcdx(a, b)


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# RUN: julia --startup-file=no %s | opt -enable-new-pm=0 -load libjulia-codegen%shlibext -AllocOpt -S - | FileCheck %s
# RUN: julia --startup-file=no %s | opt -enable-new-pm=1 --load-pass-plugin=libjulia-codegen%shlibext -passes='function(AllocOpt)' -S - | FileCheck %s

isz = sizeof(UInt) == 8 ? "i64" : "i32"

println("""
@tag = external addrspace(10) global {}
""")

# Test that the gc_preserve intrinsics are deleted directly.

# CHECK-LABEL: @preserve_branches
# CHECK: call {}*** @julia.ptls_states()
# CHECK: L1:
# CHECK-NOT: @llvm.julia.gc_preserve_begin
# CHECK-NEXT: @external_function()
# CHECK-NEXT: br i1 %b2, label %L2, label %L3

# CHECK: L2:
# CHECK: @external_function()
# CHECK-NEXT: br label %L3

# CHECK: L3:
println("""
define void @preserve_branches(i8* %fptr, i1 %b, i1 %b2) {
  %pgcstack = call {}*** @julia.get_pgcstack()
  %ptls = call {}*** @julia.ptls_states()
  %ptls_i8 = bitcast {}*** %ptls to i8*
  br i1 %b, label %L1, label %L3

L1:
  %v = call noalias {} addrspace(10)* @julia.gc_alloc_obj(i8* %ptls_i8, $isz 8, {} addrspace(10)* @tag)
  %tok = call token (...) @llvm.julia.gc_preserve_begin({} addrspace(10)* %v)
  call void @external_function()
  br i1 %b2, label %L2, label %L3

L2:
  call void @external_function()
  br label %L3

L3:
  ret void
}
""")
# CHECK-LABEL: }{{$}}

# CHECK-LABEL: @preserve_branches2
# CHECK: call {}*** @julia.ptls_states()
# CHECK: L1:
# CHECK-NEXT: @llvm.julia.gc_preserve_begin{{.*}}{} addrspace(10)* %v2
# CHECK-NEXT: @external_function()
# CHECK-NEXT: br i1 %b2, label %L2, label %L3

# CHECK: L2:
# CHECK: @external_function()
# CHECK-NEXT: br label %L3

# CHECK: L3:
println("""
define void @preserve_branches2(i8* %fptr, i1 %b, i1 %b2) {
  %pgcstack = call {}*** @julia.get_pgcstack()
  %ptls = call {}*** @julia.ptls_states()
  %ptls_i8 = bitcast {}*** %ptls to i8*
  %v2 = call {} addrspace(10)* @external_function2()
  br i1 %b, label %L1, label %L3

L1:
  %v = call noalias {} addrspace(10)* @julia.gc_alloc_obj(i8* %ptls_i8, $isz 8, {} addrspace(10)* @tag)
  %tok = call token (...) @llvm.julia.gc_preserve_begin({} addrspace(10)* %v, {} addrspace(10)* %v2)
  call void @external_function()
  br i1 %b2, label %L2, label %L3

L2:
  call void @external_function()
  br label %L3

L3:
  ret void
}
""")
# CHECK-LABEL: }{{$}}

# CHECK-LABEL: @legal_int_types
# CHECK: alloca [12 x i8]
# CHECK-NOT: alloca i96
# CHECK: ret void
println("""
define void @legal_int_types() {
  %pgcstack = call {}*** @julia.get_pgcstack()
  %ptls = call {}*** @julia.ptls_states()
  %ptls_i8 = bitcast {}*** %ptls to i8*
  %var1 = call {} addrspace(10)* @julia.gc_alloc_obj(i8* %ptls_i8, $isz 12, {} addrspace(10)* @tag)
  %var2 = addrspacecast {} addrspace(10)* %var1 to {} addrspace(11)*
  %var3 = call {}* @julia.pointer_from_objref({} addrspace(11)* %var2)
  ret void
}
""")
# CHECK-LABEL: }{{$}}

println("""
declare void @external_function()
declare {} addrspace(10)* @external_function2()
declare {}*** @julia.ptls_states()
declare {}*** @julia.get_pgcstack()
declare noalias {} addrspace(10)* @julia.gc_alloc_obj(i8*, $isz, {} addrspace(10)*)
declare {}* @julia.pointer_from_objref({} addrspace(11)*)
declare void @llvm.memcpy.p11i8.p0i8.i64(i8 addrspace(11)* nocapture writeonly, i8* nocapture readonly, i64, i32, i1)
declare token @llvm.julia.gc_preserve_begin(...)
declare void @llvm.julia.gc_preserve_end(token)
""")

# CHECK-LABEL: @memref_collision
# CHECK: call {}*** @julia.ptls_states()
# CHECK-NOT: store {}
# CHECK: store i
# CHECK-NOT: store {}
# CHECK: L1:
# CHECK: load {}
# CHECK: L2:
# CHECK: load i
println("""
define void @memref_collision($isz %x) {
  %pgcstack = call {}*** @julia.get_pgcstack()
  %ptls = call {}*** @julia.ptls_states()
  %ptls_i8 = bitcast {}*** %ptls to i8*
  %v = call noalias {} addrspace(10)* @julia.gc_alloc_obj(i8* %ptls_i8, $isz 8, {} addrspace(10)* @tag)
  %v_p = bitcast {} addrspace(10)* %v to $isz addrspace(10)*
  store $isz %x, $isz addrspace(10)* %v_p
  br i1 0, label %L1, label %L2

L1:
  %v1 = bitcast {} addrspace(10)* %v to {} addrspace(10)* addrspace(10)*
  %v1_x = load {} addrspace(10)*, {} addrspace(10)* addrspace(10)* %v1
  ret void

L2:
  %v2 = bitcast {} addrspace(10)* %v to $isz addrspace(10)*
  %v2_x = load i64, i64 addrspace(10)* %v2
  ret void
}
""")
# CHECK-LABEL: }{{$}}

max_stars_repo_name: JuliaLabs/julia

max_issues_repo_name: JuliaLabs/julia

max_forks_repo_name: JuliaLabs/julia


text:
import Knet.Ops20: batchnorm4, batchnorm4_back
using Knet.Ops20: _wsize, _lazy_init!
using CUDA: CuArray
using Knet.KnetArrays: KnetArray
using AutoGrad: AutoGrad, recording

const BN_MODE_SPATIAL = 1
const BN_MODE_ACTIVATION = 0
const CUDNN_BN_MIN_EPS = 1e-5

batchnorm4(g::KnetArray{T}, b::KnetArray{T}, x::KnetArray{T}; o...) where T = _batchnorm4(T,g,b,x;o...)
batchnorm4(g::CuArray{T}, b::CuArray{T}, x::CuArray{T}; o...) where T = _batchnorm4(T,g,b,x;o...)

# Only spatial mode is supported
# TODO: support per-activation mode
function _batchnorm4(T, g, b, x;
                     training=AutoGrad.recording(),
                     cache=nothing,
                     moments=nothing,
                     eps=1e-5,
                     alpha=1, beta=0,
                     handle = CUDNN.handle(),
                     cache_verbose=false, #reporting cache uses
                     o...)
    y = similar(x)
    weight_size = _wsize(y)
    # TODO: implement other bn mode
    bnmode = BN_MODE_SPATIAL
    # epsilon fix
    if eps < CUDNN_BN_MIN_EPS
        eps = CUDNN_BN_MIN_EPS
        @warn("eps $eps is too small for cudnn, so it is set to $CUDNN_BN_MIN_EPS", maxlog=1)
    end
    # moments
    if moments !== nothing
        _lazy_init!(moments, x)
        running_mean = moments.mean
        running_var = moments.var
        momentum = moments.momentum
    else
        running_mean, running_var = CU_NULL, CU_NULL
        momentum = .1
    end
    # The training mode
    if training
        # Cache the mean and ivar for later
        if cache !== nothing
            mean = similar(x,weight_size)
            ivar = similar(x,weight_size)
        else
            mean = CU_NULL
            ivar = CU_NULL
        end
        bnmode = CUDNN.cudnnBatchNormMode_t(bnmode)
        CUDNN.cudnnBatchNormalizationForwardTraining(handle, bnmode, Ref(T(alpha)), Ref(T(beta)), TD(x), x, TD(y), y, TD(g), g, b, momentum, running_mean, running_var, eps, mean, ivar)

        # Cache the resulting mean and inverse variance
        if cache != nothing
            cache_verbose && info("mean and ivar data saved to cache")
            cache.mean = mean
            cache.ivar = ivar
        end
    else
        @assert (moments!==nothing) "You must provide moments for the test mode!"
        bnmode = CUDNN.cudnnBatchNormMode_t(bnmode)
        CUDNN.cudnnBatchNormalizationForwardInference(handle, bnmode, Ref(T(alpha)), Ref(T(beta)), TD(x), x, TD(y), y, TD(g), g, b, running_mean, running_var, eps)
    end
    return y
end

function batchnorm4(x::Union{KnetArray,CuArray};o...)
    # Dummy buffers
    #  (cudnn doesn't support bn w/o affine
    #  although it is used in many applications)
    g = oftype(x, ones(_wsize(x)...))
    b = oftype(x, zeros(_wsize(x)...))
    return batchnorm4(g, b, x; o...)
end

batchnorm4_back(g::Union{KnetArray{T},Nothing}, x::KnetArray{T}, dy::KnetArray{T}; o...) where T = _batchnorm4_back(T,g,x,dy;o...)
batchnorm4_back(g::Union{CuArray{T},Nothing}, x::CuArray{T}, dy::CuArray{T}; o...) where T = _batchnorm4_back(T,g,x,dy;o...)

function _batchnorm4_back(T, g, x, dy;
                          training=AutoGrad.recording(),
                          cache=nothing,
                          moments=nothing,
                          grad_cache_disabled=false,
                          eps=1e-5, alpha=1, beta=0,
                          dalpha=1, dbeta=0,
                          handle = CUDNN.handle(),
                          cache_verbose=false,
                          o...)
    if training
        dx = similar(x)
        weight_size = _wsize(dy)
        if g==nothing; g=oftype(x, ones(weight_size)); end
        dg = similar(x,weight_size)
        db = similar(x,weight_size)
        # TODO: support other modes
        bnmode = BN_MODE_SPATIAL
        if cache !== nothing # (Assume cache still exists)
            mean, ivar = cache.mean, cache.ivar
            cache_verbose && info("mean and ivar are fetched from the cache")
        else
            mean, ivar = CU_NULL, CU_NULL
        end
        bnmode = CUDNN.cudnnBatchNormMode_t(bnmode)
        CUDNN.cudnnBatchNormalizationBackward(handle, bnmode, Ref(T(alpha)), Ref(T(beta)), Ref(T(dalpha)), Ref(T(dbeta)), TD(x), x, TD(dy), dy, TD(dx), dx, TD(g), g, dg, db, eps, mean, ivar)

    else
        # At test mode, g .*( x ./ sqrt(var) - mean ./ sqrt(var)) .+ beta
        # is performed;
        # so the derivative dx = dy .* g. / sqrt(var + eps) since mean and var
        # are constants
        # Note: moments must exist since otherwise forward pass fails
        ivar = 1 ./ sqrt.(moments.var .+ eps)
        dx = (g !== nothing) ? (dy .* g .* ivar) : (dy .* ivar)
        if g !== nothing
            dg = sum(dy .* (x .- moments.mean) .* ivar, dims=_reddims(dy))
            db = sum(dy, dims=_reddims(dy))
        else
            dg, db = nothing, nothing
        end
    end
    return dg, db, dx
end


text:
# module TripletSampling

using Plots
using Distributions
using TripletEmbeddings

import RecipesBase: plot

function Embedding(d::Int, n::Int, classes::Int; σ::Real=1.0, bias::Vector{Float64}=5*ones(d))
    d == length(bias) || throw(DimensionMismatch("d and length(bias) must be equal."))

    R = [[cos(θ) -sin(θ); sin(θ) cos(θ)] for θ in rand(Uniform(0, 2π), classes)]

    X = [broadcast(+, σ * randn(d, n), (c-1) * R[c] * bias) for c in 1:classes]
    return [TripletEmbeddings.Embedding(X[c]) for c in 1:classes]
end

function Embedding(d::Int, n::Vector{Int}; σ::Real=1.0, bias::Vector{Float64}=5*ones(d))
    d == length(bias) || throw(DimensionMismatch("d and length(bias) must be equal."))
    classes = length(n)

    R = [[cos(θ) -sin(θ); sin(θ) cos(θ)] for θ in rand(Uniform(0, 2π), classes)]
    X = [broadcast(+, σ * randn(d, n[c]), (c-1) * R[c] * bias) for c in 1:classes]
    return [TripletEmbeddings.Embedding(X[c]) for c in 1:classes]
end

# function core(n::Vector{Int})
#     classes == length(n) || throw(DimensionMismatch("classes != length(n)"))

#     for c in 1:classes
#         c_k = mod(c, classes) + 1
#         n_ij = n[c]   # i is always in the same class as j
#         n_k = n[c_k] # k never is in the same class as i or j
#         for k = 1:n_k, j = 1:n_ij, i = 1:n_ij
#             if i != j
#                 t_i = c > 1 ? sum(n[1:c-1]) + i : i
#                 t_j = c > 1 ? sum(n[1:c-1]) + j : j
#                 t_k = c_k > 1 ? sum(n[1:c]) + k : k
#                 @show (c, c_k, t_i, t_j, t_k, n)
#             end
#         end
#     end
# end

function TripletEmbeddings.Triplets(X::Vector{TripletEmbeddings.Embedding{Float64}})
    n = nitems.(X)
    classes = length(X)
    all(y -> y == ndims(X[1]), ndims.(X)) || throw(DimensionMismatch("All embeddings must have the same dimension"))

    triplets = Vector{Tuple{Int,Int,Int}}(undef, sum(n) * binomial(sum(n) - 1, 2))
    counter = 0

    for c in 1:classes
        c_k = mod(c, classes) + 1
        n_ij = n[c]   # i is always in the same class as j
        n_k = n[c_k] # k never is in the same class as i or j
        for k = 1:n_k, j = 1:n_ij, i = 1:n_ij
            if i != j
                counter += 1
                t_i = c > 1 ? sum(n[1:c-1]) + i : i
                t_j = c > 1 ? sum(n[1:c-1]) + j : j
                t_k = c_k > 1 ? sum(n[1:c]) + k : k
                @inbounds triplets[counter] = (t_i, t_j, t_k)
            end
        end
    end

    return Triplets(triplets[1:counter])
end

function RecipesBase.plot(X::Vector{TripletEmbeddings.Embedding{Float64}}; kwargs...)
    plot()
    foreach(x -> scatter!(x[1,:], x[2,:]; kwargs...), X)
end

function fit(loss::TripletEmbeddings.AbstractLoss, triplets::Triplets, d::Int)
    d > 0 || throw(ArgumentError("dimension d must be > 0"))

    X = TripletEmbeddings.Embedding(d, maximum(getindex.(triplets, [1 2 3])))

    violations = fit!(loss, triplets, X)

    return X
end

# end

max_stars_repo_name: kmundnic/nonuniform-triplet-sampling

max_issues_repo_name: kmundnic/nonuniform-triplet-sampling

max_forks_repo_name: kmundnic/nonuniform-triplet-sampling


text:
"""
$(TYPEDSIGNATURES)

Create Grid from TetGen data.

See [`default_options`](@ref) for available `kwargs`.

"""
function ExtendableGrids.simplexgrid(::Type{TetGenType}, TetGen, input; kwargs...)

    opts=blendoptions!(default_options();kwargs...)

    flags=makeflags(opts,:tetgen)

    if opts[:verbose]
        @show flags
    end

    if !isnothing(opts[:unsuitable])
        TetGen.tetunsuitable(opts[:unsuitable])
    end

    tetout=TetGen.tetrahedralize(input,flags)

    pointlist=tetout.pointlist

    tetrahedronlist=tetout.tetrahedronlist

    if size(tetout.tetrahedronattributelist,2)==0
        cellregions=ones(Int32,size(tetrahedronlist,2))
    else
        cellregions=Vector{Int32}(vec(tetout.tetrahedronattributelist))
    end

    segmentlist=tetout.trifacelist

    segmentmarkerlist=tetout.trifacemarkerlist

    ExtendableGrids.simplexgrid(pointlist,tetrahedronlist,cellregions,segmentlist,segmentmarkerlist)
end

"""
$(TYPEDSIGNATURES)
Create a RawTetGenIO structure
from a number of input arrays.
The 2D input arrays are transposed if necessary and converted to
the proper data types for TetGen.

This conversion is not performed if the data types are those
indicated in the defaults and the leading dimension of 2D arrays
corresponds to the space dimension.

"""
function tetgenio(TetGen;points=Array{Cdouble,2}(undef,0,0),
                  bfaces=Array{Cint,2}(undef,0,0),
                  bfaceregions=Array{Cint,1}(undef,0),
                  regionpoints=Array{Cdouble,2}(undef,0,0),
                  regionnumbers=Array{Cint,1}(undef,0),
                  regionvolumes=Array{Cdouble,1}(undef,0)
                  )
    @assert ndims(points)==2
    if size(points,2)==3
        points=transpose(points)
    end
    if typeof(points)!=Array{Cdouble,2}
        points=Array{Cdouble,2}(points)
    end
    @assert(size(points,2)>2)

    # if  ndims(bfaces)==2
    #     if size(bfaces,2)==2
    #         bfaces=transpose(bfaces)
    #     end
    #     if typeof(bfaces)!=Array{Cint,2}
    #         bfaces=Array{Cint,2}(bfaces)
    #     end
    # end
    # @assert(size(bfaces,2)>0)

    @assert ndims(bfaceregions)==1
    if  ndims(bfaces)==2
        @assert size(bfaceregions,1)==size(bfaces,2)
    else
        @assert size(bfaceregions,1)==size(bfaces,1)
    end
    if typeof(bfaceregions)!=Array{Cint,1}
        bfaceregions=Array{Cint,1}(bfaceregions)
    end

    @assert ndims(regionpoints)==2
    if size(regionpoints,1)!=3
        regionpoints=transpose(regionpoints)
    end
    if typeof(regionpoints)!=Array{Cdouble,2}
        regionpoints=Array{Cdouble,2}(regionpoints)
    end
    # @assert(size(regionpoints,2)>0)

    @assert ndims(regionnumbers)==1
    @assert ndims(regionvolumes)==1
    @assert size(regionnumbers,1)==size(regionpoints,2)
    @assert size(regionvolumes,1)==size(regionpoints,2)

    nholes=0
    nregions=0
    for i=1:length(regionnumbers)
        if regionnumbers[i]==0
            nholes+=1
        else
            nregions+=1
        end
    end

    regionlist=Array{Cdouble,2}(undef,5,nregions)
    holelist=Array{Cdouble,2}(undef,3,nholes)

    ihole=1
    iregion=1
    for i=1:length(regionnumbers)
        if regionnumbers[i]==0
            holelist[1,ihole]=regionpoints[1,i]
            holelist[2,ihole]=regionpoints[2,i]
            holelist[3,ihole]=regionpoints[3,i]
            ihole+=1
        else
            regionlist[1,iregion]=regionpoints[1,i]
            regionlist[2,iregion]=regionpoints[2,i]
            regionlist[3,iregion]=regionpoints[3,i]
            regionlist[4,iregion]=regionnumbers[i]
            regionlist[5,iregion]=regionvolumes[i]
            iregion+=1
        end
    end
    tio=TetGen.RawTetGenIO{Float64}()
    tio.pointlist=points
    if size(bfaces,2)>0
        TetGen.facetlist!(tio,bfaces)
    end
    if size(bfaceregions,1)>0
        tio.facetmarkerlist=bfaceregions
    end
    if size(regionlist,2)>0
        tio.regionlist=regionlist
    end
    if size(holelist,2)>0
        tio.holelist=holelist
    end
    tio
end

"""
$(TYPEDSIGNATURES)

Create tetgen input from the current state of the builder.
"""
function tetgenio(this::SimplexGridBuilder)
    dim_space(this)=3 || throw(error("dimension !=2 not implemented"))

   tetgenio(this.Generator,points=this.pointlist.points,
            bfaces=this.facets,
            bfaceregions=this.facetregions,
            regionpoints=this.regionpoints,
            regionnumbers=this.regionnumbers,
            regionvolumes=this.regionvolumes)
end


text:
using Sog
using Test
using Random

@testset "Sog.jl" begin
    @test isequal(sog(), nothing)
    @test isequal(sog([]), nothing)
    @test isequal(sog(['o']), Bool[1])
    @test isequal(sog(["one"]), Bool[1])
    @test isequal(sog(["one", "two", "one", "one", "two"]), Bool[1, 1, 1, 0, 1])
    Random.seed!(127)
    @test isequal(sog(rand(1:10, 20)), Bool[1, 1, 1, 1, 1, 1, 0, 1, 1, 0, 1, 1, 0, 1, 1, 1, 1, 1, 0, 1])
    @test isequal(sog([missing, missing, missing]), Bool[1, 0, 0])
    @test isequal(sog([missing, 1, missing, missing, 1, 1]), Bool[1, 1, 1, 0, 1, 0])
    @test isequal(sog(["a", "a", "a", 'a']), Bool[1, 0, 0, 1])
    @test isequal(sog([1, 1.0, 1.00000]), Bool[1, 0, 0])
    @test isequal(sog([Inf, Inf32, Inf64]), Bool[1, 0, 0])

    #For Vectors
    @test isequal(sog([[], []]),nothing)
    @test isequal(sog([[1, 1, 1, 2, 2], [1.5, 1.5, 20.0, 3.0, 3.0]]), Bool[1, 0, 1, 1, 0])
    @test isequal(sog([["abc", "ab", "ab"],["ab", "ab", "ab"]]), Bool[1, 1, 0])

    #xxx = [round.(rand(10^6)), rand(1:10, 10^6), [randstring('a':'e', 2) for _ in 1:10^6]]
    #@time sog(xxx)
end


text:
# Testing file for creating useful circuit decompositions

import QuantumCircuitOpt as QCOpt
using JuMP
using CPLEX
using Gurobi

# using Cbc
include("optimizers.jl")
include("2qubit_gates.jl")
include("3qubit_gates.jl")

#-------------------------------#
#      User-defined inputs      #
#-------------------------------#
function input_circuit()
         # [(depth, gate)]
    return [(1, "CNot_2_1"),
            (2, "S_1"),
            (3, "H_2"),
            (4, "S_2")
            ]
end

function target_gate(gate::Int)
    if gate == 1
        return QCOpt.kron_single_qubit_gate(2, QCOpt.U3Gate(0,0,π/4), "q1")
    elseif gate == 2
        return QCOpt.kron_single_qubit_gate(2, QCOpt.RXGate(π/4), "q1") * QCOpt.kron_single_qubit_gate(2, QCOpt.RYGate(π/4), "q2") * QCOpt.kron_single_qubit_gate(2, QCOpt.RZGate(π/4), "q1")
    elseif gate == 3
        return QCOpt.iSwapGate()
    elseif gate == 4
        return QCOpt.kron_single_qubit_gate(3, QCOpt.RXGate(π/4), "q3")
    elseif gate == 5
        return Array{Complex{Float64},2}([1 0 0 0; 0 1 0 0; 0 0 1 0; 0 0 0 -1])
    end
end

params = Dict{String, Any}(
"num_qubits" => 3,
"maximum_depth" => 7,
"elementary_gates" => ["RY_3", "CNot_1_2", "CNot_2_3", "CNot_1_3", "Identity"],
# "elementary_gates" => ["CV_1_2", "CV_2_3", "CV_1_3", "CVdagger_1_2", "CVdagger_2_3", "CVdagger_1_3", "CNot_1_2", "CNot_3_2", "CNot_2_3", "CNot_1_3", "Identity"],
"target_gate" => target_gate(5),
"objective" => "minimize_depth",

"RY_discretization" => -π:π/4:π,

# "input_circuit" => toffoli_circuit(),
)

#------------------------------#
#      Optimization model      #
#------------------------------#
qcopt_optimizer = get_gurobi()

# data = QCOpt.get_data(params)
result_qc = QCOpt.run_QCModel(params, qcopt_optimizer)


text:
# TODO:
# - Allow passing the system path if the simulation wasn't serialized
# - Handle PER-UNIT conversion of variables according to type
# - Enconde Variable/Parameter/Dual from other inputs to avoid passing Symbol

const ResultsByTime = SortedDict{Dates.DateTime, DataFrames.DataFrame}
const FieldResultsByTime = Dict{Symbol, ResultsByTime}

"""Holds the results of a simulation problem for plotting or exporting"""
mutable struct SimulationProblemResults <: PSIResults
    problem::String
    base_power::Float64
    execution_path::String
    results_output_folder::String
    existing_timestamps::StepRange{Dates.DateTime, Dates.Millisecond}
    results_timestamps::Vector{Dates.DateTime}
    system::Union{Nothing, PSY.System}
    system_uuid::Base.UUID
    resolution::Dates.TimePeriod
    forecast_horizon::Int
    end_of_interval_step::Int
    variable_values::FieldResultsByTime
    dual_values::FieldResultsByTime
    parameter_values::FieldResultsByTime
end

function SimulationProblemResults(
    store::SimulationStore,
    problem_name::AbstractString,
    problem_params::SimulationStoreProblemParams,
    sim_params::SimulationStoreParams,
    path;
    results_output_path = nothing,
)
    if results_output_path === nothing
        results_output_path = joinpath(path, "results")
    end

    time_steps = range(
        sim_params.initial_time,
        length = problem_params.num_executions * sim_params.num_steps,
        step = problem_params.interval,
    )
    name = Symbol(problem_name)
    variables = list_fields(store, name, STORE_CONTAINER_VARIABLES)
    parameters = list_fields(store, name, STORE_CONTAINER_PARAMETERS)
    duals = list_fields(store, name, STORE_CONTAINER_DUALS)

    return SimulationProblemResults(
        problem_name,
        problem_params.base_power,
        path,
        results_output_path,
        time_steps,
        Vector{Dates.DateTime}(),
        nothing,
        problem_params.system_uuid,
        get_resolution(problem_params),
        get_horizon(problem_params),
        get_end_of_interval_step(problem_params),
        _fill_result_value_container(variables),
        _fill_result_value_container(duals),
        _fill_result_value_container(parameters),
    )
end

function Base.empty!(res::SimulationProblemResults)
    foreach(empty!, _get_dicts(res))
    empty!(res.results_timestamps)
    return
end

Base.isempty(res::SimulationProblemResults) = all(isempty, _get_dicts(res))

# This returns the number of timestamps stored in all containers.
Base.length(res::SimulationProblemResults) = mapreduce(length, +, _get_dicts(res))

get_problem_name(res::SimulationProblemResults) = res.problem
get_system(res::SimulationProblemResults) = res.system
get_resolution(res::SimulationProblemResults) = res.resolution
get_forecast_horizon(res::SimulationProblemResults) = res.forecast_horizon
get_end_of_interval_step(res::SimulationProblemResults) = res.end_of_interval_step
get_execution_path(res::SimulationProblemResults) = res.execution_path
get_existing_variables(res::SimulationProblemResults) = collect(keys(res.variable_values))
get_existing_duals(res::SimulationProblemResults) = collect(keys(res.dual_values))
get_existing_parameters(res::SimulationProblemResults) = collect(keys(res.parameter_values))
get_existing_timestamps(res::SimulationProblemResults) = res.existing_timestamps
get_model_base_power(res::SimulationProblemResults) = res.base_power
IS.get_timestamp(result::SimulationProblemResults) = result.results_timestamps

get_interval(res::SimulationProblemResults) = res.existing_timestamps.step
IS.get_variables(result::SimulationProblemResults) = result.variable_values
get_duals(result::SimulationProblemResults) = result.dual_values
IS.get_parameters(result::SimulationProblemResults) = result.parameter_values
IS.get_base_power(result::SimulationProblemResults) = result.base_power

"""
Return the system used for the problem. If the system hasn't already been deserialized or
set with [`set_system!`](@ref) then deserialize and store it.
"""
function get_system!(results::SimulationProblemResults)
    file = joinpath(
        results.execution_path,
        "problems",
        make_system_filename(results.system_uuid),
    )
    results.system = PSY.System(file, time_series_read_only = true)
    return results.system
end

"""
Set the system in the results instance.

Throws InvalidValue if the system UUID is incorrect.
"""
function set_system!(results::SimulationProblemResults, system::PSY.System)
    sys_uuid = IS.get_uuid(system)
    if sys_uuid != results.system_uuid
        throw(
            IS.InvalidValue(
                "System mismatch. $sys_uuid does not match the stored value of $results.system_uuid",
            ),
        )
    end

    results.system = system
end

_get_containers(x::SimulationProblemResults) =
    (x.variable_values, x.parameter_values, x.dual_values)
_get_dicts(res::SimulationProblemResults) =
    (y for x in _get_containers(res) for y in values(x))

function _get_store_value(
    res::SimulationProblemResults,
    field::Symbol,
    names::Vector{Symbol},
    timestamps,
    ::Nothing,
)
    simulation_store_path = joinpath(get_execution_path(res), "data_store")
    return h5_store_open(simulation_store_path, "r") do store
        _get_store_value(res, field, names, timestamps, store)
    end
end

function _get_store_value(
    res::SimulationProblemResults,
    field::Symbol,
    names::Vector{Symbol},
    timestamps,
    store::SimulationStore,
)
    results = Dict{Symbol, SortedDict{Dates.DateTime, DataFrames.DataFrame}}()
    problem_name = Symbol(get_problem_name(res))
    problem_interval = get_interval(res)
    resolution = get_resolution(res)
    horizon = get_forecast_horizon(res)
    for name in names
        _results = SortedDict{Dates.DateTime, DataFrames.DataFrame}()
        for ts in timestamps
            out = read_result(DataFrames.DataFrame, store, problem_name, field, name, ts)
            time_col = range(ts, length = horizon, step = resolution)
            DataFrames.insertcols!(out, 1, :DateTime => time_col)
            _results[ts] = out
        end
        results[name] = _results
    end

    return results
end

function _validate_names(existing_names::Vector{Symbol}, names::Vector{Symbol})
    for name in names
        if name ∉ existing_names
            @error("$name is not stored", sort(existing_names))
            throw(IS.InvalidValue("$name is not stored"))
        end
    end
    nothing
end

function _process_timestamps(
    res::SimulationProblemResults,
    initial_time::Union{Nothing, Dates.DateTime},
    count::Union{Int, Nothing},
)
    if initial_time === nothing
        initial_time = first(get_existing_timestamps(res))
    end
    existing_timestamps = get_existing_timestamps(res)

    if initial_time ∉ existing_timestamps
        invalid_timestamps = [initial_time]
    else
        if count === nothing
            requested_range = [v for v in existing_timestamps if v >= initial_time]
        else
            requested_range =
                collect(range(initial_time, length = count, step = get_interval(res)))
        end
        invalid_timestamps = [v for v in requested_range if v ∉ existing_timestamps]
    end
    if !isempty(invalid_timestamps)
        @error "Timestamps $(invalid_timestamps) not stored" get_existing_timestamps(res)
        throw(IS.InvalidValue("Timestamps not stored"))
    end
    return requested_range
end

function _read_variables(
    res::SimulationProblemResults,
    names::Vector{Symbol},
    timestamps,
    store,
)
    isempty(names) && return FieldResultsByTime()
    existing_names = get_existing_variables(res)
    _validate_names(existing_names, names)
    same_time_stamps = isempty(setdiff(res.results_timestamps, timestamps))
    names_with_values = [k for (k, v) in res.variable_values if !isempty(v)]
    same_names = isempty([n for n in names if n ∉ names_with_values])
    if same_time_stamps && same_names
        @info "reading variables from SimulationsResults"
        vals = filter(p -> (p.first ∈ names), res.variable_values)
    else
        @info "reading variables from data store"
        vals = _get_store_value(res, STORE_CONTAINER_VARIABLES, names, timestamps, store)
    end
    return vals
end

"""
    Returns the values for the requested variable names. Accepts a vector of names for the
    return of the values. If the time stamps and names are loaded using the [load_results!](@ref)
    function it will read from memory.

    # Accepted Key Words
    - `names::Vector{Symbol}` : names of desired results
    - `initial_time::Dates.DateTime` : initial of the requested results
    - `count::Int`: Number of results
    - `store::SimulationStore`: a store that has been opened for reading
"""
function read_variables(
    res::SimulationProblemResults;
    names::Union{Vector{Symbol}, Nothing} = nothing,
    initial_time::Union{Nothing, Dates.DateTime} = nothing,
    count::Union{Int, Nothing} = nothing,
    store = nothing,
)
    names = isnothing(names) ? collect(keys(res.variable_values)) : names
    timestamps = _process_timestamps(res, initial_time, count)
    values = _read_variables(res, names, timestamps, store)
    return values
end

function _read_duals(
    res::SimulationProblemResults,
    names::Vector{Symbol},
    timestamps,
    store,
)
    isempty(names) && return FieldResultsByTime()
    existing_names = get_existing_duals(res)
    _validate_names(existing_names, names)
    same_time_stamps = isempty(setdiff(res.results_timestamps, timestamps))
    names_with_values = [k for (k, v) in res.dual_values if !isempty(v)]
    same_names = isempty([n for n in names if n ∉ names_with_values])
    if same_time_stamps && same_names
        @debug "reading duals from SimulationsResults"
        vals = filter(p -> (p.first ∈ names), res.dual_values)
    else
        @debug "reading duals from data store"
        vals = _get_store_value(res, STORE_CONTAINER_DUALS, names, timestamps, store)
    end
    return vals
end

"""
    Returns the values for the requested dual names. It must match the duals requested in the simulation problem definition.
    It keeps requests when performing multiple retrievals. Accepts a vector of names for the return of the values

    # Accepted Key Words
    - `names::Vector{Symbol}` : names of desired results
    - `initial_time::Dates.DateTime` : initial of the requested results
    - `count::Int`: Number of results
    - `store::SimulationStore`: a store that has been opened for reading
"""
function read_duals(
    res::SimulationProblemResults;
    names::Union{Vector{Symbol}, Nothing} = nothing,
    initial_time::Union{Nothing, Dates.DateTime} = nothing,
    count::Union{Int, Nothing} = nothing,
    store = nothing,
)
    names = isnothing(names) ? collect(keys(res.dual_values)) : names
    timestamps = _process_timestamps(res, initial_time, count)
    values = _read_duals(res, names, timestamps, store)
    return values
end

function _read_parameters(
    res::SimulationProblemResults,
    names::Vector{Symbol},
    timestamps,
    store,
)
    isempty(names) && return FieldResultsByTime()
    existing_names = get_existing_parameters(res)
    _validate_names(existing_names, names)
    same_time_stamps = isempty(setdiff(res.results_timestamps, timestamps))
    names_with_values = [k for (k, v) in res.parameter_values if !isempty(v)]
    same_names = isempty([n for n in names if n ∉ names_with_values])
    if same_time_stamps && same_names
        @info "reading parameters from SimulationsResults"
        vals = filter(p -> (p.first ∈ names), res.parameter_values)
    else
        @info "reading parameters from data store"
        vals = _get_store_value(res, STORE_CONTAINER_PARAMETERS, names, timestamps, store)
    end
    return vals
end

"""
    Returns the values for the parameters used in the simulation. It keeps requests when performing multiple retrievals. Accepts a vector of names for the return of the values

    # Accepted Key Words
    - `names::Vector{Symbol}` : names of desired results
    - `initial_time::Dates.DateTime` : initial time of the requested results
    - `count::Int`: Number of results
    - `store::SimulationStore`: a store that has been opened for reading
"""
function read_parameters(
    res::SimulationProblemResults;
    names::Union{Vector{Symbol}, Nothing} = nothing,
    initial_time::Union{Nothing, Dates.DateTime} = nothing,
    count::Union{Int, Nothing} = nothing,
    store = nothing,
)
    names = isnothing(names) ? collect(keys(res.parameter_values)) : names
    timestamps = _process_timestamps(res, initial_time, count)
    values = _read_parameters(res, names, timestamps, store)
    return values
end

"""
    Returns the values for the requested variable name. It keeps requests when performing multiple retrievals. Accepts a variable name to return the result.

    # Accepted Key Words
    - `initial_time::Dates.DateTime` : initial of the requested results
    - `count::Int`: Number of results
    - `store::SimulationStore`: a store that has been opened for reading
"""
function read_variable(res::SimulationProblemResults, name::Symbol; kwargs...)
    return read_variables(res; names = [name], kwargs...)[name]
end

"""
    Returns the values for the requested dual name. It keeps requests when performing multiple retrievals. Accepts a dual name to return the result.
    # Accepted Key Words
    - `initial_time::Dates.DateTime` : initial of the requested results
    - `count::Int`: Number of results
    - `store::SimulationStore`: a store that has been opened for reading
"""
function read_dual(res::SimulationProblemResults, name::Symbol; kwargs...)
    return read_duals(res; names = [name], kwargs...)[name]
end

"""
    Returns the values for the requested parameter name. It keeps requests when performing multiple retrievals. Accepts a parameter name to return the result.
    # Accepted Key Words
    - `initial_time::Dates.DateTime` : initial of the requested results
    - `count::Int`: Number of results
"""
function read_parameter(res::SimulationProblemResults, name::Symbol; kwargs...)
    return read_parameters(res; names = [name], kwargs...)[name]
end

"""
Returns the optimizer stats for the problem as a DataFrame.

# Accepted keywords
- `store::SimulationStore`: a store that has been opened for reading
"""
function read_optimizer_stats(res::SimulationProblemResults; store = nothing)
    return _read_optimizer_stats(res, store)
end

function _read_optimizer_stats(res::SimulationProblemResults, ::Nothing)
    h5_store_open(joinpath(get_execution_path(res), "data_store"), "r") do store
        _read_optimizer_stats(res, store)
    end
end

function _read_optimizer_stats(res::SimulationProblemResults, store::SimulationStore)
    return read_problem_optimizer_stats(store, Symbol(res.problem))
end

struct RealizedMeta
    initial_time::Dates.DateTime
    resolution::Dates.TimePeriod
    count::Int
    start_offset::Int
    end_offset::Int
    interval_len::Int
    end_of_interval_step::Int
end

function RealizedMeta(
    res::SimulationProblemResults;
    initial_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Int, Nothing} = nothing,
)
    existing_timestamps = get_existing_timestamps(res)
    interval = existing_timestamps.step
    resolution = get_resolution(res)
    interval_len = Int(interval / resolution)
    end_of_interval_step = get_end_of_interval_step(res)
    realized_timestamps =
        get_realized_timestamps(res, initial_time = initial_time, len = len)

    result_initial_time = existing_timestamps[findlast(
        x -> x .<= first(realized_timestamps),
        existing_timestamps,
    )]
    result_end_time = existing_timestamps[findlast(
        x -> x .<= last(realized_timestamps),
        existing_timestamps,
    )]

    count = length(result_initial_time:interval:result_end_time)

    start_offset = length(result_initial_time:resolution:first(realized_timestamps))
    end_offset = length(
        (last(realized_timestamps) + resolution):resolution:(result_end_time + interval - resolution),
    )

    return RealizedMeta(
        result_initial_time,
        resolution,
        count,
        start_offset,
        end_offset,
        interval_len,
        end_of_interval_step,
    )
end

function get_realized_timestamps(
    res::SimulationProblemResults;
    initial_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Int, Nothing} = nothing,
)
    existing_timestamps = get_existing_timestamps(res)
    interval = existing_timestamps.step
    resolution = get_resolution(res)
    horizon = get_forecast_horizon(res)
    initial_time = isnothing(initial_time) ? first(existing_timestamps) : initial_time
    end_time =
        isnothing(len) ? last(existing_timestamps) + interval - resolution :
        initial_time + (len - 1) * resolution

    requested_range = initial_time:resolution:end_time
    available_range = first(existing_timestamps):resolution:(last(
        existing_timestamps,
    ) + (horizon - 1) * resolution)
    invalid_timestamps = setdiff(requested_range, available_range)

    if !isempty(invalid_timestamps)
        msg = "Requested time does not match available results"
        @error msg
        throw(IS.InvalidValue(msg))
    end

    return requested_range
end

function get_realization(
    result_values::Dict{Symbol, SortedDict{Dates.DateTime, DataFrames.DataFrame}},
    meta::RealizedMeta,
    timestamps,
)
    realized_values = Dict{Symbol, DataFrames.DataFrame}()
    for (key, result_value) in result_values
        results_concat = Dict{Symbol, Vector{Float64}}()
        for (step, (t, df)) in enumerate(result_value)
            first_id = step > 1 ? 1 : meta.start_offset
            last_id =
                step == meta.count ? meta.interval_len - meta.end_offset : meta.interval_len
            result_length = length(first_id:last_id)
            for colname in propertynames(df)
                colname == :DateTime && continue
                if meta.end_of_interval_step == 1 # indicates RH
                    col = ones(result_length) .* df[!, colname][1] # realization is first period setpoint
                else
                    col = df[!, colname][first_id:last_id]
                end
                if !haskey(results_concat, colname)
                    results_concat[colname] = col
                else
                    results_concat[colname] = vcat(results_concat[colname], col)
                end
            end
        end
        realized_values[key] = DataFrames.DataFrame(results_concat, copycols = false)
        DataFrames.insertcols!(realized_values[key], 1, :DateTime => timestamps)
    end
    return realized_values
end

"""
    Returns the final values for the requested variable names for each time step for a problem.
    Accepts a vector of names for the return of the values. If the time stamps and names are
    loaded using the [load_results!](@ref) function it will read from memory.

    # Accepted Key Words
    - `names::Vector{Symbol}` : names of desired results
    - `initial_time::Dates.DateTime` : initial time of the requested results
    - `len::Int`: length of results
"""
function read_realized_variables(
    res::SimulationProblemResults;
    names::Union{Vector{Symbol}, Nothing} = nothing,
    initial_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Int, Nothing} = nothing,
)
    names = isnothing(names) ? collect(keys(res.variable_values)) : names
    meta = RealizedMeta(res, initial_time = initial_time, len = len)
    result_values = read_variables(
        res,
        names = names,
        initial_time = meta.initial_time,
        count = meta.count,
    )
    timestamps = get_realized_timestamps(res, initial_time = initial_time, len = len)
    return get_realization(result_values, meta, timestamps)
end

"""
    Returns the final values for the requested parameter names for each time step for a problem.
    Accepts a vector of names for the return of the values. If the time stamps and names are
    loaded using the [load_results!](@ref) function it will read from memory.

    # Accepted Key Words
    - `names::Vector{Symbol}` : names of desired results
    - `initial_time::Dates.DateTime` : initial time of the requested results
    - `len::Int`: length of results
"""
function read_realized_parameters(
    res::SimulationProblemResults;
    names::Union{Vector{Symbol}, Nothing} = nothing,
    initial_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Int, Nothing} = nothing,
)
    names = isnothing(names) ? collect(keys(res.parameter_values)) : names
    meta = RealizedMeta(res, initial_time = initial_time, len = len)
    result_values = read_parameters(
        res,
        names = names,
        initial_time = meta.initial_time,
        count = meta.count,
    )
    timestamps = get_realized_timestamps(res, initial_time = initial_time, len = len)
    return get_realization(result_values, meta, timestamps)
end

"""
    Returns the final values for the requested dual names for each time step for a problem.
    Accepts a vector of names for the return of the values. If the time stamps and names are
    loaded using the [load_results!](@ref) function it will read from memory.

    # Accepted Key Words
    - `names::Vector{Symbol}` : names of desired results
    - `initial_time::Dates.DateTime` : initial time of the requested results
    - `len::Int`: length of results
"""
function read_realized_duals( # TODO: Should this be get_realized_duals_values?
    res::SimulationProblemResults;
    names::Union{Vector{Symbol}, Nothing} = nothing,
    initial_time::Union{Nothing, Dates.DateTime} = nothing,
    len::Union{Int, Nothing} = nothing,
)
    names = isnothing(names) ? collect(keys(res.dual_values)) : names
    meta = RealizedMeta(res, initial_time = initial_time, len = len)
    result_values =
        read_duals(res, names = names, initial_time = meta.initial_time, count = meta.count)
    timestamps = get_realized_timestamps(res, initial_time = initial_time, len = len)
    return get_realization(result_values, meta, timestamps)
end

"""
    Loads the simulation results into memory for repeated reads. Running this function twice
    overwrites the previously loaded results. This is useful when loading results from remote
    locations over network connections

    # Required Key Words
    - `initial_time::Dates.DateTime` : initial of the requested results
    - `count::Int`: Number of results
    # Accepted Key Words
    - `variables::Vector{Symbol}`: Variables names to load into results
    - `duals::Vector{Symbol}`: Dual names to load into results
    - `parameters::Vector{Symbol}`: Parameter names to load into results
"""
function load_results!(
    res::SimulationProblemResults,
    count::Int;
    initial_time::Union{Dates.DateTime, Nothing} = nothing,
    variables::Vector{Symbol} = Symbol[],
    duals::Vector{Symbol} = Symbol[],
    parameters::Vector{Symbol} = Symbol[],
)
    initial_time =
        isnothing(initial_time) ? first(get_existing_timestamps(res)) : initial_time

    res.results_timestamps = _process_timestamps(res, initial_time, count)

    simulation_store_path = joinpath(res.execution_path, "data_store")
    h5_store_open(simulation_store_path, "r") do store
        merge!(
            res.variable_values,
            _read_variables(res, variables, res.results_timestamps, store),
        )
        merge!(res.dual_values, _read_duals(res, duals, res.results_timestamps, store))
        merge!(
            res.parameter_values,
            _read_parameters(res, parameters, res.results_timestamps, store),
        )
    end

    return nothing
end

#= NEEDS RE-IMPLEMENTATION
""" Exports the results in the SimulationProblemResults object to  CSV files"""
function write_to_CSV(res::SimulationProblemResults; kwargs...)
    folder_path = res.results_output_folder
    if !isdir(folder_path)
        throw(IS.ConflictingInputsError("Specified path is not valid. Set up results folder."))
    end
    variables_export = Dict()
    for (k, v) in IS.get_variables(res)
        start = decode_symbol(k)[1]
        if start !== "ON" || start !== "START" || start !== "STOP"
            variables_export[k] = get_model_base_power(res) .* v
        else
            variables_export[k] = v
        end
    end
    parameters_export = Dict()
    for (p, v) in IS.get_parameters(res)
        parameters_export[p] = get_model_base_power(res) .* v
    end
    write_data(variables_export, res.time_stamp, folder_path; file_type = CSV, kwargs...)
    write_optimizer_stats(IS.get_total_cost(res), folder_path)
    write_data(IS.get_timestamp(res), folder_path, "time_stamp"; file_type = CSV, kwargs...)
    write_data(get_duals(res), folder_path; file_type = CSV, kwargs...)
    write_data(parameters_export, folder_path; file_type = CSV, kwargs...)
    files = readdir(folder_path)
    compute_file_hash(folder_path, files)
    @info("Files written to $folder_path folder.")
    return
end
=#


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiStorageV1beta1CSIDriverSpec <: SwaggerModel
    attachRequired::Any # spec type: Union{ Nothing, Bool } # spec name: attachRequired
    podInfoOnMount::Any # spec type: Union{ Nothing, Bool } # spec name: podInfoOnMount
    volumeLifecycleModes::Any # spec type: Union{ Nothing, Vector{String} } # spec name: volumeLifecycleModes

    function IoK8sApiStorageV1beta1CSIDriverSpec(;attachRequired=nothing, podInfoOnMount=nothing, volumeLifecycleModes=nothing)
        o = new()
        validate_property(IoK8sApiStorageV1beta1CSIDriverSpec, Symbol("attachRequired"), attachRequired)
        setfield!(o, Symbol("attachRequired"), attachRequired)
        validate_property(IoK8sApiStorageV1beta1CSIDriverSpec, Symbol("podInfoOnMount"), podInfoOnMount)
        setfield!(o, Symbol("podInfoOnMount"), podInfoOnMount)
        validate_property(IoK8sApiStorageV1beta1CSIDriverSpec, Symbol("volumeLifecycleModes"), volumeLifecycleModes)
        setfield!(o, Symbol("volumeLifecycleModes"), volumeLifecycleModes)
        o
    end
end # type IoK8sApiStorageV1beta1CSIDriverSpec

const _property_map_IoK8sApiStorageV1beta1CSIDriverSpec = Dict{Symbol,Symbol}(Symbol("attachRequired")=>Symbol("attachRequired"), Symbol("podInfoOnMount")=>Symbol("podInfoOnMount"), Symbol("volumeLifecycleModes")=>Symbol("volumeLifecycleModes"))
const _property_types_IoK8sApiStorageV1beta1CSIDriverSpec = Dict{Symbol,String}(Symbol("attachRequired")=>"Bool", Symbol("podInfoOnMount")=>"Bool", Symbol("volumeLifecycleModes")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiStorageV1beta1CSIDriverSpec }) = collect(keys(_property_map_IoK8sApiStorageV1beta1CSIDriverSpec))
Swagger.property_type(::Type{ IoK8sApiStorageV1beta1CSIDriverSpec }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiStorageV1beta1CSIDriverSpec[name]))}
Swagger.field_name(::Type{ IoK8sApiStorageV1beta1CSIDriverSpec }, property_name::Symbol) =  _property_map_IoK8sApiStorageV1beta1CSIDriverSpec[property_name]

function check_required(o::IoK8sApiStorageV1beta1CSIDriverSpec)
    true
end

function validate_property(::Type{ IoK8sApiStorageV1beta1CSIDriverSpec }, name::Symbol, val)
end